On 03/30/2018 08:21 AM, William F Pokorny wrote:
....

Sturm / polysove() investigation. Those Difficult Coefficients.

When we have 4th order equations we are today running this piece of code:

   /* Test for difficult coeffs. */

   if (difficult_coeffs(4, c))
   {
       sturm = true;
   }

in the Solve_Polynomial wrapper for sturm/polysolve() and three specific 
solvers including the one for 4th order equations in solve_quartic(). 
The difficult_coeffs() function scans the coefficients looking for large 
differences in values. If found, user end up running sturm for some rays 
no matter what was specified in the scene SDL. The code existed in the 
original POV-Ray 1.0, but in a slightly modified form and used with an 
older version of solve_quartic().

I've had the feeling for a while the code was probably doing nothing 
useful for us today. I was almost completely right. Close enough it's 
going into the bit bucket with this set of commits.

The code was set up to crash anytime difficult_coeffs() returned true. 
Scanned hundreds of scenes. Turned up helix.pov and a handful of scenes 
with blobs. Ran those hard coded to use only solve_quartic() and not a 
single difference in result. Changed the code above to return 
immediately with no roots on difficult_coeffs() returning true. The 
helix scene was unchanged as the 'rays' were not part of the resultant 
shape.

The blobs were a different story. Many rays are being re-directed to 
sturm / polysolve for no benefit given the current solve_quartic code. 
To oblivion with that code you all say!

I thought the same, but remembered we have known blob accuracy issues in 
https://github.com/POV-Ray/povray/issues/187, maybe half a dozen related 
problem scenes in my test space. Plus to support the subsurface feature 
Chrisoph added this bit of code in blob.cpp to avoid visible seams:

   DBL depthTolerance = (ray.IsSubsurfaceRay()? 0 : DEPTH_TOLERANCE);

I've been running with DEPTH_TOLERANCE at 1e-4 for a long time in my 
personal version of POV-Ray which fixes all the blob accuracy issues 
except for the subsurface one.

Given difficult_coeffs() is almost exclusively active with blobs, 
updated blob.cpp to run with accuracy addressing the issues in the 
previous paragraph. With this update in place, we do occasionally see 
cases where the bump into to sturm / polysolve method is necessary. 
However, the automatic difficult_coeffs() method doesn't currently catch 
all of the fringe case rays.

It could be adjusted so as to push more solve_quartic rays into sturm / 
polysolve, but it's aleady the case most difficult-coefficient rays are 
unnecessarily run in sturm / polysolve when solve_quartic would work 
just fine. The right thing to do is let users decide when to use 'sturm' 
and this set of commit dumps the difficult_coeffs() code.

Examples.

The attached image is a scene originally created to explore default 
texturing related to negative blobs, but it demonstrates nicely the 
issues described above.

In the top row we have on the left the image rendered with the blob 
accuracy unmodified and all the difficult_coeffs() code running as it 
does today. In the middle image the 'if difficult_coeffs()' test has 
been hard coded to return no roots. On the right the difference image 
shows all the rays running in the sturm / polysolve method which all run 
fine - and faster given the blob set up - with solve_quartic().

In the second row we are running with the new more accurate blob.cpp. 
Further, the blobs in the scene have been changed to have a threshold of 
0.001 instead of 0.1 - values this low are not common because you don't 
get much blobbing.

The first column of the second row was rendered with a version of 
POV-Ray compiled to always force the use of solve_quartic(). The middle 
column restores the 'if difficult_coeffs()' test and internal swap to 
sturm / polysolve. It does help as shown in column 3. However column 2 
shows the already too aggressive use of sturm / polysolve with blobs is 
not enough to automatically 'fix' all the solve_quartic issues. Running 
with sturm is clean in all cases.

Aside: The potentials as glows image I did sometime back had what I 
thought were media related speckles. I applied very aggressive and 
expensive AA to hide them as I usually do with media. Turns out using 
sturm would have been the better choice. I wanted no blobbing for the 
blob cylinders and used a threshold of 1e-5. It was a case of outrunning 
the effectiveness of the 'if difficult_coeffs() patch' with existing 
blob accuracy, but I didn't realize it at the time.

Updates at:

https://github.com/wfpokorny/povray/tree/fix/polynomialsolverAccuracy

Performance and specific commit details below.

Bill P.

Performance info:
------------ /usr/bin/time povray -j -wt1 -fn -p -d -c lemonSturm.pov

0)  master    30.22user 0.04system 0:30.89elapsed

10) 0f6afcc   14.77user 0.02system 0:15.36elapsed  -51.13%
11) 0f9509b   14.77user 0.01system 0:15.36elapsed
12) c590c0e   15.08user 0.01system 0:15.71elapsed   +2.10%
13) 4ab17df   14.88user 0.02system 0:15.51elapsed   -1.33%  -50.76%
14) 78e664c   --- NA ---
15) 5d8cedc   --- NA ---
16) 21387e2   --- NA ---

---- povray subsurface.pov -p +wt2 -j -d -c -fn
master  sturm on       540.89user 0.06system 4:31.15elapsed
21387e2 sturm on       347.29user 0.07system 2:54.42elapsed  -35.79%

master  sturm off      286.99user 0.06system 2:24.28elapsed
21387e2 sturm off      278.54user 0.05system 2:20.12elapsed   -2.94%

If blob container needs sturm.  278.54 -> 347.29  +24.68%


11) Changing polysolve visible_roots and numchanges zero test to <=0.

Depending upon value used for reducing and trimming the sturm chain 
equations in modp, the sphere_sweep b_spline functions can return 
negative root counts when there are no roots.

12) Adding sanity test of upper bound Cauchy method result.

In practice have occasionally found roots slightly above upper bound. 
Perhaps related to how sturm chain is tweaked during pruning. Sort it later.

13) Removing pointless operations in polysolve and Solve_Polynomial.

visible_roots() passed a couple integer pointers which it set prior to 
return but the values were never used in polysolve itself.

When fourth order equations seen a test call to difficult_coeffs() was 
made which if true would set sturm true - even when sturm was already true.

14) Updating blob accuracy so better synchronized with other shapes.

Fix for GitHub issue #187

Returned minimum intersection depths now set to MIN_ISECT_DEPTH.

Internal determine_influences() intersect_<sub-element> functions using 
 >0.0 as was already the subsurface feature. Anything more opens up gaps 
or jumps where the sub-element density influences are added too late or 
dropped too earlier. The change does create ray equations slighly more 
difficult to solve which was likely the reason for the original, 
largish, v1.0 1e-2 value. The sturm / polysolve solver is sometimes 
necessary over solve_quartic().

15) Removing difficult_coeffs() and related code in Solve_Polynomial().

Code not completely effective where sturm / polysolve needed and very 
costly in most scenes with 4th order equations due use of sturm over 
solve_quartic() where solver_quartic works fine and is much faster.

16) Reversing change to sor.cpp as unneeded with updated polysolve().

Not verified, but suspect bump made in commit 3f2e1a9 to get sign change 
was previously needed where regula_falsa() handed single root problem 
with 0 to MAX_DISTANCE range. In any case, the patch is no longer needed 
with updated polysolve() code.

Post a reply to this message

Attachments:
Download 'storydifficultcoeff.png' (110 KB)

Preview of image 'storydifficultcoeff.png'