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 P_{2}(x) = (3x^{2}  1) / 2. Even in ALL mode it came back within a few seconds with a result of 8,16 E13, i.e. approx. zero for a 12digit 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 3E14 and 6E15 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 (3E16, 4E17, 1E16, ...) 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 010, 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:
x[approx]? Y
JMP int_done
RCL r_SkAbs
RCL[times] r_deltau
/
ABS
Num 1
ULP
x<? Y
JMP int_next_size
/* else exit */
/* Two matches in a row is goodness. */
int_done::
[cmplx]RCL cr_limits
STO L
[cmplx]x[<>] Z
TOP?
MSG MSG_INTEGRATE
RTN
You should check the way ULP works in XROM code  it should return 1E15 in SP and 1E33 in DP mode.
In ALL 03 mode, the integral of the originally posted function is approximated this way:
0,0634765625
2,0751953125 E4
9,15527343744 E6
8,021902 E17
... and after this the iteration terminates with a final (standard precision) result of 7,81753878289 E17.
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:
5,68490384437 E19
5,73024595021 E19
5,72957175023 E19
5,72957800022 E19
5,72957795128 E19
5,72957795131 E19
This matches the true result of 10^{20} * 180/Pi.
Might this be a feasible (and correct) solution?
Dieter
Addendum:
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.
