|
 |
On 03/30/2018 08:21 AM, William F Pokorny wrote:
>...
Sturm / polysove() investigation. Coefficients wrongly seen as zeroes.
In this chapter we venture up into the Solve_Polynomial() function which
calls polysove(), solve_quadratic(), solve_cubic() and solve_quartic().
Consider this code added back in 1994 and part of POV-Ray v3.0.
/*
* Determine the "real" order of the polynomial, i.e.
* eliminate small leading coefficients.
*/
i = 0;
while ((i < n) && (fabs(c0[i]) < SMALL_ENOUGH))
{
i++;
}
n -= i;
c = &c0[i];
My last commit deletes it! Having so large a SMALL_ENOUGH (1e-10) and
any fixed value given equation order varies as does distance along the
ray where the equations not normalized is a bad idea(1).
Perhaps reason enough that when any even order polynomial with no roots
is reduced to an odd order polynomial it turns into a polynomial with
roots. All odd order polynomials have at least one root.
In the attached image we have two examples - neither using sturm /
polysolve - where the code above goes bad.
In the top row we've got the shipped primativ.pov scene and specifically
a poly shape which is not explicitly bounded. It creates 4th order
equations so should use solve_quartic(). However the code up top
sometimes strips the leading coefficient and in this case we go from a
quartic which has no roots to an odd third order cubic which has one.
This is the vertical bar artifact in the upper left.
In the bottom row we have the shipped subsurface.pov scene. The
differences on the bottom right (at 4x actual difference values) happen
again where the code at top drops a single leading coefficient. In this
case for pixels and rays where I've examine actual roots, subsurface
related rays go from 4th order equations with >=2 roots with one of
those roots inside the blob interval being kept, to a 3rd order equation
with a single root outside the active blob interval. In other words we
are missing subsurface roots/intersections where we see differences.
Notes:
- The subsurface effect here is visible with no difference multiplier
and many times larger than the returned root threshold differences
recently discussed.
- Because the ray-surface equations are based upon scene geometry the
'zeroes-not-zeroes-dropped' issue here is often scene-object or axis
aligned in some fashion. I'm reminded of a rounded box canyon subsurface
scene Robert McGregor posted where I noticed one box face having a
strange look. Someone joked it was perhaps a door to Narnia or similar
after nobody had an idea what caused the effect. I'm thinking perhaps
the code above - or the still present order reduction by ratio code -
played a part...
NOTE! This update corrupts some previous scenes while fixing some
previous bad ones. Some planned but not yet defined or implemented
equation tuning on ray DX,DY,DZ values is to come.
Updates at:
https://github.com/wfpokorny/povray/tree/fix/polynomialsolverAccuracy
Performance and specific commit details below.
Bill P.
(1) - Yes, this strongly suggests order reduction based upon the largest
two coefficients is also a bad idea. However leaving that feature in
place for now as some primitives require this type of reduction to
function at all. Perhaps in these cases push the reductions out of the
solvers and into the shape code.
Performance info:
------------ /usr/bin/time povray -j -wt1 -fn -p -d -c lemonSturm.pov
0) master 30.22user 0.04system 0:30.89elapsed
16) 21387e2 14.88user 0.02system 0:15.51elapsed -50.76%
17) 80860cd --- NA --- (1e-15 constant calc fix)
(PRECISE_FLOAT)
18) aa6a0a6 --- Bad Result(a) --- (float)
14.95user 0.02system 0:15.56elapsed +0.47% (double)
25.12user 0.02system 0:25.73elapsed +68.82% (long double)
282.97user 0.13system 4:44.13elapsed +1801.68% (__float128)
19) 4e16623 14.94user 0.01system 0:15.56elapsed ---- -50.56%
(a) - Not exhaustively evaluated, but some shapes render OK at 'float'
though those which did were slower (+150% ballpark). Unsure why so slow,
but don't plan to dig further. Sturm / polysolve is for those looking to
be more accurate, not less.
17) Correcting two double 1e-15 constant calculations.
In recent commits started setting two 1e-15 values with the calculation:
= 1/std::numeric_limits<DBL>::digits10;
when I meant to code:
= 1/pow(10,std::numeric_limits<DBL>::digits10);
Both forms bad practice in not being specific as to the types for the
explicit values. The first returns integer 0, later cast to double and I
got away with it in all but my full set of test cases. Corrected:
= (DBL)1.0/pow((DBL)10.0,std::numeric_limits<DBL>::digits10);
18) Implementing suggested PRECISE_FLOAT macro mechanism for polysolve().
Primarily enabling 'long double' support for the polysolve() / 'sturm'
solver given 'long double' is part of the C+11 standard and gets most
users 80 bits or more as opposed to 64. Further, the additional bits are
almost always hardware backed making for reasonable performance.
With GNU g++ 64 bit environments '__float128' is a valid setting with no
additional library required - though it's very slow. Other extended
accuracy floats such as 'double double' enabled, via the PRECISE_
mechanisms, but they involve additional libraries and settings. See code.
Note! With this update the polysolve sturm chain is no longer pruned
with a default 1e-10 value, but rather a value appropriate for the
floating point type. Means better accuracy for some previously bad root
results and less accuracy for others. Coming updates will better tune
the chain equation pruning.
19) Removing leading, near zero coefficient reduction in Solve_Polynomial().
Leading coefficients of < SMALL_ENOUGH (1e-10) were being dropped
reducing the order of the incoming polynomial when the near zero values
had meaning with respect to roots. Led to both false roots and missed roots.
Post a reply to this message
Attachments:
Download 'notzeroesdropped.png' (326 KB)
Preview of image 'notzeroesdropped.png'
|
|
