Posts: 1,368
Threads: 212
Joined: Dec 2006
One of my recent acquistions is a 33C, and I positively love itit is probably one of the best purchases of the several I have recently made.
I am really growing more interested in compact, powerful programs, but the lack of DSE and ISG features on 33C makes writing with finite loops a challenge. Let's say my counter is in register 7. My 33C equivalent of a DSE instruction is 1 ST7 RCL7 x<>0?so three extra steps. Also, in my particular program of interest now I have to restore the contents of X register for the further calculations to work, so that takes an extra step.
In a 49 step calculator this is pricey indeed! Anyone out there have an ingenious and more economical work around? Fortunately (though it is confusing to keep track) GTO and GSB have direct line addressing so that saves a step here and there since LBL doesn't exist to take up a line.
Eager for direction. FWIW, I am borrowing heavily from Valentin's ingenious 46step Gaussian quadrature routine for the 11C. His is an iteration of the 3 point quadrature, whereas mine uses only 2 points (iterated over subintervals) to save 10 steps so that there are some program lines left to store the function to be integrated. When I iron it out I will share it. I think it does a little better than Simpson's, though of course not as well as Valentin's original routine.
Best,
Les
Posts: 1,792
Threads: 62
Joined: Jan 2005
Hi, Les 
Based on lack of LBL, DSE, and ISG, it sounds like the very limited parogramming paradigm developed for the HP55 and carried over to the HP25, HP10C, and HP12C, as well. (That's a shame...)
I'm surprised that the HP33C has GSB, because at least some of the others don't.
I don't know of any shortcuts. The scientific Spiceseries model that is most worth having, of course, is the HP34C. DSE, ISG, LBL, 70210 steps, and the first with SOLVE and INTEG. A 198283 version with better construction techniques is more desirable.
 KS
Posts: 1,477
Threads: 71
Joined: Jan 2005
If you only need 2 working registers (R0 and R1) and one counter (R2) you can use [Sigma][X=0] as a DSElike instruction and LASTx will do what you need.
Posts: 4,587
Threads: 105
Joined: Jul 2005
Hi Karl,
You wrote:
Quote:
Based on lack of LBL, DSE, and ISG, it sounds like the very limited parogramming paradigm developed for the HP55 and carried over to the HP25, HP10C, and HP12C, as well. (That's a shame...)
Sorry, I agree for the latest models only. Please look into THE museum to judge the early calculators fairly. AFAIR, not many people had the $$$ to buy an HP65 in the seventies. So, a smaller CMOSmodel was a way to get a reliable PRogrammable for less.
Regards, Walter
Posts: 1,368
Threads: 212
Joined: Dec 2006
Without further ado, here is a 2point Gaussian quadrature routine for the 33C, which is basically Valentin Abillo's excellent 3point routine with a few steps saved since one of the function computations per subinterval is removed to make room for more steps from step 39 on for the function to be integrated:
01 STO 1
02 
03 RCL 7
04 /
05 STO 0
06 2
07 /
08 ST+ 1
09 3
10 SQRT
11 1/X
12 *
13 STO 3
14 CLX
15 STO 2
16 RCL 1
17 RCL 3
18 +
19 GSB 39
20 ST 2
21 RCL 1
22 RCL 3
23 
24 GSB 39
25 ST 2
26 RCL 0
27 ST+ 1
28 1
29 ST 7
30 RCL 7
31 X!=0?
32 GTO 16
33 RCL 0
34 RCL 2
35 *
36 2
37 /
38 RTN
To use, put the function to be integrated from line 39 on, making sure the result is returned to the X register, and end that segment with RTN. Then store the number of desired subintervals in R7, and the limits of integration on the stack as lower ENTER upper. Then GSB 01 and voila!
This is not nearly as effective as Valentin's routine. On a given interval of say width h, both threepoint Gaussian quadrature and the 5point NewtonCotes formula (Bode's Rule) have roughly the same theoretical error, though of course the Gaussian version saves two function evaluations on the first subinterval and one on every subinterval thereafter (since with Bode's Rule the connecting endpoints need be computed only once). So if you divide up the area of integration into 4 subintervals, Bode's Rule requires 5 + 3*4 = 17 function evaluations, whereas threepoint Gaussian quadrature requires 12. It is easy to see that if we divide up the area of integration into 6 subintervals for a 3point Gaussian integration, the accuracy should be much better than Bode's Rule at the cost of only one extra evaluation of the function. I would guess that it is also superior Simpson's rule computed over nine subintervals (19 function evalutions), but since the subinterval size is changing this is a little harder for me to demonstrate to myself.
With 2point Gaussian quadrature, the theoretical error on an interval of width h is of the order h^5/4320 (times some value of the fourth derivative of the function being integrated), vs. 3point NewtonCotes quadrature (Simpson's Rule), which is of order h^5/2880. So, at least in theory, perceptibly better. However, the savings in function evaluation isn't as impressiveSimpson's rule computes three times on the first subinterval, and twice on every subinterval thereafter, whereas repeated application of 2point Gaussian quadrature is simply twice per subinterval. So there is a saving of only one function evalution. However, I would guess, on balance, this routine is preferable to Simpson's rule as you get more accuracy with one fewer function evaluations, and one need not compute the function at the endpoints.
As a quick comparison, I computed the integral of sin(x)/x from 0 to 2 using the HP41 Math Pac's INTG routine (Simpson's Rule) with n = 12 (13 function evaluations), Valentin's original routine with 4 sub intervals (12 function evaluations), and the above routine with 6 subintervals (12 function evaluations). This is what I get:
Simpson's: 1.605413995
Abillo's: 1.605412978
Les's: 1.605412298
The actual 10 digit answer is 1.605412977, so Valentin's routine gives a full 10 digits within 1 ULP. Both mine and Simpson's get it to 7 digits within 1 ULP after rounding, but mine is a wee bit closer when you take it out to 8 digits7 ULP vs. 10 ULP. This is not earthshattering, but it is gratifying since it got a little closer with one less evaluation of the function.
I doubt that this routine has any practical implications given the superior integrators in HP calculators starting from the 34C and right up to the otherwise much maligned 33S and of course delightful 42S. But I do like the challenge of the small program, and I positively love the 33C, despite its modest powers. It was a great pleasure for me to adapt Valentin's excellent work to a more modest end, and I hope that the pleasure can be shared.
Les
Edited: 24 Dec 2006, 4:03 a.m.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Nice Christmas Eve work, Les!
And you have left enough room to key in Prof. Kahan's weird example in HewlettPackard Journal (Aug1980):
39 ENTER
40 SQRT
41 LST x
42 1
43 
44 /
45 x<>y
46 LN
47 1/x
48 
49 RTN
The HP34C shows the correct answer in 23 minutes (in FIX 7): 0.03649 +/ 0.00000007 (On the 48GX in STD mode: 0.0364899739831 +/ 3.649E13). On my HP33C your program gives 0.036492416 in 23 minutes (N=300 !). But that's a weird integral. Anyway, for N=5 the answer 0.036493448 is pretty good!
Best regards,
Gerson.
Posts: 1,368
Threads: 212
Joined: Dec 2006
Hey, none of these basic "carve up the area of integration into equal strips" approaches is any match for the adaptive quadrature routines built into the 34C, 42S, 41 Advantage Pac, 33S, 15C, etc.... Also, the good performance with 5 subintervals is a bit of a fluke, since I find for higher values it doesn't do much better and sometimes a little worse. BTW, on my 11C Valentin's original 3point routine fares about as well, and doesn't seem to have the clear advantage it does with the well behaved Sine integral.
In Mathematica I differentiated Prof Kahan's expression four times and plotted the fourth derivative. For x less than about 0.05 the 4th derivative, which is a strong contributor to the error term in both 2 point Gaussian quadrature and Simpson's Rule, just whizzes off to negative infinity, so the error in that subinterval will figure prominently even if one carves the area up into smaller strips. Likewise with Valentin's 3point routinein this case, the plot of the sixth derivative does the same thing, which explains why his routine doesn't seem to confer much advantage over minethat pesky first interval with the sharp peak and the wild high derivatives defies any hope that a polynomial can well estimate the function close to the lower limit of integration.
Long live adaptive quadrature routines!!!
Les
Posts: 1,792
Threads: 62
Joined: Jan 2005
Merry Christmas, Gerson (and to all!) 
Quote:
And you have left enough room to key in Prof. Kahan's weird example in HewlettPackard Journal (Aug1980):
39 ENTER
40 SQRT
41 LST x
42 1
43 
44 /
45 x<>y
46 LN
47 1/x
48 
49 RTN
The integral from the HP Journal article is of the function
sqrt(x) 1
f(x) =    between 0 and 1
x1 ln(x)
which has a nearlyvertical "lover's leap" trace near x = 0, although the value of the function is finite, as is the integral.
The correct answer is near 0.03649, but the article points out that the HP34C will return a value near 0.03662 using FIX 5 as the setting for inputfunction uncertainty. This result is not within the maximum estimated error at this setting (about 0.000005) of the correct answer for the integral. The integration algorithm fails to capture the squirrely behavior of the function near x = 0 at fhe FIX 5 uncertainty setting.
However, using FIX 6, the HP34C will give a correct answer within the maximum estimated error (about 0.0000005).
The Journal article goes on to show that speed and accuracy can be improved by transforming the integrand into a function that has a smoother trace.
The HP15C, HP41 Advantage Pac, HP32S, and HP32SII will return results identical to those of the HP34C, although at widely differing speeds of computation. The HP42S, the HP48 series and their derivatives, as well as the HP33S, will give different results because they use a "proportional" method of specifying the uncertainty.
For more details, see my article (#556) in the Articles Archive.
 KS
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hi Karl,
Quote:
The integral from the HP Journal article is of the function
sqrt(x) 1
f(x) =    between 0 and 1
x1 ln(x)
Thanks for taking to time to typing the function in algebraic notation. I ought to have done that as it's not so apparent from the listing.
Quote:
The Journal article goes on to show that speed and accuracy can be improved by transforming the integrand into a function that has a smoother trace
Namely:
2w^{2} w
f(w) =    between 0 and 1
(w1)(w+1) ln(w)
According to the article, the HP34C computes this as 0.03649 +/ 0.000005 in 100 seconds at FIX 5. It is also recommended not to replace (w  1)(w + 1) by (w^{2}  1) in order to avoid loss of half of the significant digits and a spike when w approaches 1. However I used the simplification below (not optimized routine) and obtained the same result in 90 seconds at FIX 5:
 
 2w 1 
f(w) = w *      between 0 and 1
 w^{2}  1 ln(w) 
 
Unfortunately it appears these don't fit in eleven steps so they cannot be tried in Les's program. Anyway the situation is better than in my first numerical integration program on the TI51III, when just introduced to Calculus: with 32 steps only the builtin functions could be integrated, which was really useless :)
Quote:
Merry Christmas, Gerson (and to all!) 
I wish you a Merry Christmas to you and your dear ones!
Gerson.
Edited: 24 Dec 2006, 7:06 p.m.
Posts: 1,792
Threads: 62
Joined: Jan 2005
Hello again, Gerson 
Quote:
2w^{2} w
f(w) =    between 0 and 1
(w1)(w+1) ln(w)
According to the article, the HP34C computes this as 0.03649 +/ 0.000005 in 100 seconds at FIX 5. It is also recommended not to replace (w  1)(w + 1) by (w^{2}  1) in order to avoid loss of half of the significant digits and a spike when w approaches 1.
Thanks  I was too lazy to look up the article on my hard drive or printouts. (BTW, the article is available on the MoHPC CD/DVD set.)
The substitution is w^{2} = x, with 2w.dw = dx.
Strictly speaking, the recommendation in the HP Journal article is to replace (w^{2}  1) with (w  1)(w + 1) after substituting w^{2} in place of x, in order to reduce roundoff errors near w = 0.
A few years ago, I benchmarked a number of HP calc's to see how they handled this calcuation  results and computational time. I never posted the incomplete results. Since then, I acquired a few more units that include numerical integration. The thing I discovered, though, is that the differing methods of specifying uncertainty accounts almost entirely for the differences in results among HP units, and plays a role for the differences in timings. It seems as though there have been few changes to the mathematical algorithms. So, my benchmark comparisons have little meaningful value.
 KS
Warm Christmas wishes in kind for you and your family.
Edited: 1 Jan 2007, 1:24 a.m.
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
Unfortunately it appears these don't fit in eleven steps so they cannot be tried in Les's program.
No, but I was able to check out Valentin's original 3point routine on my 11C.
Using 50 subintervals (150 evaluations of the function), I get .036489941. Takes a generous 5 minutes and 40 seconds.
Les
Posts: 1,368
Threads: 212
Joined: Dec 2006
P.S. I would question the optimistic error bounds of the 48GX STD mode determination. In both Maple and Mathematica, I get the following as the first 15 digits of the integral: 0.0364899739785765. The 48GX result begins to disagree in about the 10th or 11th digit, which seems to put the actual error outside of the 48GX's estimate. Actually, it is more like 4.5E12, not 3.6E13.
Les
Posts: 2,761
Threads: 100
Joined: Jul 2005
Sorry, Les, for the typo. That's the actual 48GX estimate:
3.64899739883E12
It's 30 minutes past midnight here but still in time to wish you a Merry Christmas in Canada :)
Gerson.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Quote:
Using 50 subintervals (150 evaluations of the function), I get .036489941. Takes a generous 5 minutes and 40 seconds.
That's exactly what I get on my spedup 15C... in 3 minutes and 26 seconds (no optimized routine)! That's what spedup units are for :)
Gerson.

...and 6 minutes and 45 seconds on a normal 15C. Considering the 11C is faster than the 15C (1.14x in average) and the speedup factor is 1.996, 3 minutes and 26 second is exactly what we would have expected.
Edited: 24 Dec 2006, 10:32 p.m.
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
According to the article, the HP34C computes this as 0.03649 +/ 0.000005 in 100 seconds at FIX 5. It is also recommended not to replace (w  1)(w + 1) by (w2  1) in order to avoid loss of half of the significant digits and a spike when w approaches 1. However I used the simplification below (not optimized routine) and obtained the same result in 90 seconds at FIX 5:
 
 2w 1 
f(w) = w *      between 0 and 1
 w^2  1 ln(w) 
 
Hi Gerson
I have been working with this version in my explorations of the PPC ROM IG routine on my HP41CV. This factorization does not avoid the roundoff error when I try to compute the integral for SCI 6, which of course gives 7 significant digits (FIX 5 yields only 4). The PPC routine permits printing out of intermediate results and after a point the thing actually diverges. Changing the denominator in question to the (w1)(w+1) form as recommended solves the problem.
Happy New Year!
Les
