|Re: WP34s integration trapped in infinite loop|
Message #8 Posted by Dieter on 13 Oct 2013, 12:10 p.m.,
in response to message #7 by Paul Dale
I tried the originally posted function on my trusted 35s with P2(x) = (3x2 - 1) / 2. Even in ALL mode it came back within a few seconds with a result of 8,16 E-13, i.e. approx. zero for a 12-digit machine like this one. HP's integration algorithm uses a modified Romberg method, so there must be a way to overcome the limitations of the current 34s implementation.
The crucial point here is the "approx. equal" test in the current 34s Romberg code. In FIX n mode, this simply tests if the two last approximations agree if rounded to n decimals. This also works for values close to zero - for instance 3E-14 and 6E-15 are both rounded to zero, the two values agree, et voilà: the iteration terminates.
On the other hand, ALL, SCI and ENG will round the two last approximations to a certain number of significant digits. So if the iteration oscillates through several values very close to zero (3E-16, 4E-17, 1E-16, ...) they will almost never agree in some or even all significant digits, and the iteration will not terminate (at least not before 16000+ function calls).
I examined the current 34s integration code, transferred it to the emulator (with a few adjustments, e.g. it uses the regular registers 0-10, and 11ff for the Romberg matrix), and finally I changed the exit condition in the following way:
After the call of the user function and the summation via STO+ r_Sk, also the absolute value of this increment is accumulated in a new register r_SkAbs, which of course has been set to zero at the beginning. This way, SkAbs times deltaU represents a kind of average function value over the interval. If now the final x[approx]y test fails, but the last approximation is very small compared to this average, the iteration can terminate:
/* else exit */
/* Two matches in a row is goodness. */
You should check the way ULP works in XROM code - it should return 1E-15 in SP and 1E-33 in DP mode.
In ALL 03 mode, the integral of the originally posted function is approximated this way:
... and after this the iteration terminates with a final (standard precision) result of 7,81753878289 E-17.
On the other hand, an integral that really is that small is evaluated correctly without premature break. For instance, the integral of sin(x) * 10-20 for x = 0..90 degrees (!) returns the following approximation sequence:
This matches the true result of 10-20 * 180/Pi.
Might this be a feasible (and correct) solution?
The x[approx] test compares two rounded values. It seems that, just like a manual ROUND command, this way at most 12 digits get compared. So the result will usually not carry 16 valid digits in SP mode, let alone 34 digits in DP mode: everything beyond 12 digits is not guaranteed. IMHO this is one more reason to get rid of the x[approx] test in the 34s integration routine.
Edited: 13 Oct 2013, 6:10 p.m.