The Museum of HP Calculators

HP Forum Archive 21

 WP34s integration trapped in infinite loopMessage #1 Posted by Bernd Grubert on 10 Oct 2013, 4:32 a.m. Hello, Sorry if that has been adressed before. I tried to integrate the 2nd Lengendre polynom over -1 ... 1, which should give a result close to zero. With "ALL 00" the integration seems to be trapped in an infinite loop. The problem doesn't occur with "FIX 11". Here is the function, that I'm integrating: ```0001 ****LBL'LE2' 0002 **LBL A 0003 2 0004 x[<->] Y 0005 P[sub-n] 0006 END ``` The integration is done the following way: ```1 ENTER[^] +/- [integral]'LE2' ``` I'm using Version 3458 with printer support. BTW: The WP 34S is a really fine machine. Thanks to all who made it possible! Regards Bernd

 Re: WP34s integration trapped in infinite loopMessage #2 Posted by Paul Dale on 10 Oct 2013, 5:37 a.m.,in response to message #1 by Bernd Grubert It isn't in an infinite loop. The Romberg integration is having difficulties and it will call your user subroutine 16383 times before returning the answer 6.612851843479552E-16. The exact integral here is (x3 - x) / 2 + C, when taken over the interval specified gives zero. The integration routine has difficulties with integrals that evaluate to zero when in ALL mode. The problem being related to the approximately equal function which rounds to the current display mode -- with zero, the rounding doesn't do much of anything except in FIX mode. Interestingly, version 2 firmware, which uses a fixed point Gauss-Kronrod quadrature, will get the correct result very quickly. This quadrature is known exact for polynomials of such low degree. - Pauli Edited: 10 Oct 2013, 5:38 a.m.

 Re: WP34s integration trapped in infinite loopMessage #3 Posted by Paul Dale on 10 Oct 2013, 5:40 a.m.,in response to message #1 by Bernd Grubert A work around would be to break the interval of integration somewhere that doesn't produce a zero result. E.g. integrate -1 .. 0.5 and 0.5 .. 1 and add the results together. Unfortunately, -1 .. 0 and 0 .. 1 both run into the same problem. - Pauli

 Re: WP34s integration trapped in infinite loopMessage #4 Posted by Bernd Grubert on 10 Oct 2013, 4:55 p.m.,in response to message #3 by Paul Dale Thanks for the answer. I understand that that 16383 is not infinite. But it is quite close - at least with the real calc running on battries :) At first I thought that this problem has to do with the ALL -setting, but I noticed that this problem occurs also with SCI 10. So maybe for me the best solution is to set the accuracy not too high when calculating integrals. Regards Bernd

 Re: WP34s integration trapped in infinite loopMessage #5 Posted by Dieter on 10 Oct 2013, 5:19 p.m.,in response to message #4 by Bernd Grubert Assuming the current implementation, the problem will occur with any display setting except FIX, i.e. in ALL, SCI and ENG. Try FIX 11 and the result (approx. 1,07 E-16) appears after a few seconds. Dieter

 Re: WP34s integration trapped in infinite loopMessage #6 Posted by fhub on 10 Oct 2013, 5:38 p.m.,in response to message #5 by Dieter Quote: Assuming the current implementation, the problem will occur with any display setting except FIX, i.e. in ALL, SCI and ENG. Well, then the current implementation is not very clever - needing >16000 steps for a Romberg integration is definitely unacceptable. Franz

 Re: WP34s integration trapped in infinite loopMessage #7 Posted by Paul Dale on 10 Oct 2013, 6:18 p.m.,in response to message #6 by fhub Anyone is most welcome to submit a replacement integration routine. The source code is in trunk/xrom/integrate.wp34s and is no more than a standard keystroke program. I've implemented two different integrators so far and don't really want to do a third. - Pauli

 Re: WP34s integration trapped in infinite loopMessage #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: ``` x[approx]? Y JMP int_done RCL r_SkAbs RCL[times] r_deltau / ABS Num 1 ULP x] Z TOP? MSG MSG_INTEGRATE RTN ``` 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: ``` 0,0634765625 -2,0751953125 E-4 -9,15527343744 E-6 -8,021902 E-17 ``` ... 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: ``` 5,68490384437 E-19 5,73024595021 E-19 5,72957175023 E-19 5,72957800022 E-19 5,72957795128 E-19 5,72957795131 E-19 ``` 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.

 Re: WP34s integration trapped in infinite loopMessage #9 Posted by walter b on 13 Oct 2013, 7:13 p.m.,in response to message #8 by Dieter Hallo Dieter, Thanks for investigating. I leave the case to Pauli. d:-)

 Re: WP34s integration trapped in infinite loopMessage #10 Posted by Paul Dale on 13 Oct 2013, 8:01 p.m.,in response to message #8 by Dieter I've add these changes but can't build new images at the moment. Assuming I got everything correct. - Pauli

 Re: WP34s integration trapped in infinite loopMessage #11 Posted by Dieter on 14 Oct 2013, 8:13 a.m.,in response to message #10 by Paul Dale There still is the question of how the ULP function works in SP resp. DP mode for XROM code. Since the x[approx] test disregards any digits beyond the 12th anyway, I would suggest replacing the 1 ULP sequence by a simple 1 SDR 15, i.e. a constant 1E-15. A better solution would get rid of the x[approx] test completely. Dieter

 Re: WP34s integration trapped in infinite loopMessage #12 Posted by Paul Dale on 14 Oct 2013, 8:46 a.m.,in response to message #11 by Dieter ULP should work the same in xrom code as elsewhere. The only thing xrom changes when activated is to force a four level stack. The xIN command does switch to double precision and an eight level stack, but the integration code doesn't call this function. The integration code uses the approx test so that it has some kind of compatibility with HP's earlier integration code -- the screen settings changed the result. We could lose this compatibility and do something different I guess. The approx test only fails when the integral is convergent to zero and the code you've already supplied seems to handle this case. - Pauli

 Re: WP34s integration trapped in infinite loopMessage #13 Posted by Marcus von Cube, Germany on 14 Oct 2013, 11:40 a.m.,in response to message #12 by Paul Dale I've uploaded a new build (3460) to SF with Pauli's latest version of the integrator for testing.

 Re: WP34s integration trapped in infinite loopMessage #14 Posted by Dieter on 14 Oct 2013, 2:45 p.m.,in response to message #13 by Marcus von Cube, Germany Say hello to v. 3461 - see below. ;-) Dieter

 Re: WP34s integration trapped in infinite loopMessage #15 Posted by Marcus von Cube, Germany on 14 Oct 2013, 3:32 p.m.,in response to message #14 by Dieter Quote: Say hello to v. 3461 If Pauli changes the code and I'll do the build, the final revision will be 3462.

 Re: WP34s integration trapped in infinite loopMessage #16 Posted by Dieter on 14 Oct 2013, 2:43 p.m.,in response to message #12 by Paul Dale Pauli, Quote: ULP should work the same in xrom code as elsewhere. I thought XROM routines always run in DP mode, so that ULP(1) will always return 1E-33. But anyway - here is a bugfix. There still is an error that returns the wrong values on the stack if the integral is close to zero. Also, since the x[approx] test ignores everything beyond the first 12 digits anyway, I now prefer this version: ``` x[approx]? Y JMP int_done ENTER RCL/ r_SkAbs RCL/ r_deltau ABS SDL 15 x>1? JMP int_next_size DROP /* 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 ``` If you want to save 15 powers of 10 in the result, here is an alternative code: ``` x[approx]? Y JMP int_done ENTER RCL/ r_SkAbs RCL/ r_deltau ABS Num 1 SDR 15 x] Z TOP? MSG MSG_INTEGRATE RTN ``` Time for a new build. :-) Dieter Edited: 14 Oct 2013, 2:45 p.m.

 Re: WP34s integration trapped in infinite loopMessage #17 Posted by Dieter on 14 Oct 2013, 3:03 p.m.,in response to message #16 by Dieter An old German proverb says: "Das einzig Beständige ist der Wandel". So here's a more conservative version that uses an additional (initially cleared) flag that is set for a final iteration. In other words: if the algorithm has detected that the integral is very close to zero, it runs one final iteration - just to be sure. This way the last two close-to-zero-results are returned in X and Y. ``` FC?C flag_final /* "final" flag set */ x[approx]? Y /* OR x~y ? */ JMP int_done RCL/ r_SkAbs RCL/ r_deltau ABS Num 1 SDR 15 x>? Y /* integral close to zero? */ SF flag_final /* then do one final iteration */ JMP int_next_size    /* Two matches in a row is goodness. */   int_done:: [cmplx]RCL cr_limits STO L [cmplx]x[<->] Z TOP? MSG MSG_INTEGRATE RTN ``` Dieter Edited: 14 Oct 2013, 3:09 p.m.

 Re: WP34s integration trapped in infinite loopMessage #18 Posted by Paul Dale on 14 Oct 2013, 5:35 p.m.,in response to message #17 by Dieter I've implemented this version. - Pauli

 Re: WP34s integration trapped in infinite loopMessage #19 Posted by Marcus von Cube, Germany on 15 Oct 2013, 1:28 p.m.,in response to message #18 by Paul Dale A new revision (3463) is built and online on SF.

 Re: WP34s integration trapped in infinite loopMessage #20 Posted by David Maier on 18 Oct 2013, 10:47 a.m.,in response to message #17 by Dieter Just thought I would point out that, according to the printed version of the WP 34S Owner's Manual, integration returns "an upper limit of the uncertainty in Y". There's also a footnote that mentions the 34C and 15C Owner's Handbooks. The 34C and 15C do return a calculated uncertainty in Y, and their handbooks include a discussion of how this is calculated.

 Re: WP34s integration trapped in infinite loopMessage #21 Posted by Dieter on 18 Oct 2013, 3:05 p.m.,in response to message #20 by David Maier That's right - the 34C and 15C (and newer machines like the 35s, for that matter) return an estimate for the uncertainty in Y. However, the 34s does not. Neither before nor after the last modification of the integration routine. Actually the exit code remained unchanged. On exit, the stack contains these values: ``` T lower limit a Z upper limit b Y previous approximation X result (final approximation) L upper limit b ``` Yes, the 34s manual contains various references to other HP manuals, and there are cases where the 34s behaves differently compared to earlier HPs. This is one example, and it should be corrected in the manual. However, my PDF manual does not seem to contain the quoted "upper limit of the uncertainty in Y". Maybe someone knows a good way to estimate the integral's uncertainty for the Romberg algorithm used here. Dieter

 Re: WP34s integration trapped in infinite loopMessage #22 Posted by Paul Dale on 18 Oct 2013, 8:28 p.m.,in response to message #21 by Dieter STO- Y would give an uncertainty estimate of sorts. - Pauli

 Re: WP34s integration trapped in infinite loopMessage #23 Posted by Paul Dale on 14 Oct 2013, 5:27 p.m.,in response to message #16 by Dieter Quote:I thought XROM routines always run in DP mode, so that ULP(1) will always return 1E-33. It is the xIN that starts most XROM commands that switches to DP mode. Before that it stays in whatever mode the device is currently in. Solve, integrate, product, sum and the derivatives don't call xIN anywhere so they stay single precision and don't gain the benefits of automatic argument processing and Last X handling. - Pauli

 Re: WP34s integration trapped in infinite loopMessage #24 Posted by Dieter on 14 Oct 2013, 2:50 p.m.,in response to message #8 by Dieter I wrote: Quote: ... and after this the iteration terminates with a final (standard precision) result of 7,81753878289 E-17. This is caused by an error that is corrected in the bugfix I posted below. With this correction, the last approximation -8,021902 E-17 is returned in X (and the second last -9,155 E-6 in Y). Dieter

 Re: WP34s integration trapped in infinite loopMessage #25 Posted by Miguel Toro on 16 Oct 2013, 9:31 a.m.,in response to message #24 by Dieter Hi, This is what I like of this project and why I continue to enjoy working (and playing, of course) with my trusted WP34s : the constant and really fast improvement of an already marvelous little machine. Thank you Paul, Walter, Marcus and Dieter Regards, Miguel

 Re: WP34s integration trapped in infinite loopMessage #26 Posted by Dieter on 17 Oct 2013, 8:50 a.m.,in response to message #25 by Miguel Toro Thank you very much, Miguel. We should now discuss whether the new version includes a valid and realiable exit condition. In the updated algorithm this is done with an additional approximation for the integral of |f(x)|: ``` Integral |f(x)| ~= 1,5 * r_ba4 * r_SkAbs * r_deltau ``` The basic idea now is simple: compare this integral of |f(x)| with the integral of f(x). If the latter is, say, 15 magnitudes smaller than the former, the algorithm assumes that the true value of the integral is zero, so it exits after one final iteration. Since we compare just magnitudes, the factor 1,5 in the formula is omitted in the program. Also, it does not divide by r_ba4, which may cause problems if the interval is very small or very large (|b-a| << 1 or >> 1). So I think this should be added: ``` ... FC?C flag_final /* "final" flag set */ x[approx]? Y /* OR x~y ? */ JMP int_done RCL/ r_SkAbs RCL/ r_deltau RCL/ r_ba4 ABS ... ``` Finally, there still is one case where the program may exit with an error. This would happen if the sum SkAbs of all sampled |(f(x)| is zero, i.e. f(x) is zero for each and every sampled x-value. This results in a division by zero. On the other hand, I think that in this case the two initial approximations of the integral would both be zero as well so that the algorithm exits with a result of zero. At least I tried this with f(x) = 0 and it worked. ;-) Any objections or suggestions? Dieter

Go back to the main exhibit hall