Long Live the 33C!!!



#2

One of my recent acquistions is a 33C, and I positively love it--it 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 ST-7 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 46-step 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


#3

Hi, Les --

Based on lack of LBL, DSE, and ISG, it sounds like the very limited parogramming paradigm developed for the HP-55 and carried over to the HP-25, HP-10C, and HP-12C, as well. (That's a shame...)
I'm surprised that the HP-33C has GSB, because at least some of the others don't.

I don't know of any shortcuts. The scientific Spice-series model that is most worth having, of course, is the HP-34C. DSE, ISG, LBL, 70-210 steps, and the first with SOLVE and INTEG. A 1982-83 version with better construction techniques is more desirable.

-- KS


#4

Hi Karl,

You wrote:

Quote:
Based on lack of LBL, DSE, and ISG, it sounds like the very limited parogramming paradigm developed for the HP-55 and carried over to the HP-25, HP-10C, and HP-12C, 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 CMOS-model was a way to get a reliable PRogrammable for less.

Regards, Walter

#5

If you only need 2 working registers (R0 and R1) and one counter (R2) you can use [Sigma-][X=0] as a DSE-like instruction and LASTx will do what you need.

#6

Without further ado, here is a 2-point Gaussian quadrature routine for the 33C, which is basically Valentin Abillo's excellent 3-point 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 three-point Gaussian quadrature and the 5-point Newton-Cotes 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 three-point Gaussian quadrature requires 12. It is easy to see that if we divide up the area of integration into 6 subintervals for a 3-point 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 2-point 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. 3-point Newton-Cotes 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 impressive--Simpson's rule computes three times on the first subinterval, and twice on every subinterval thereafter, whereas repeated application of 2-point 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 digits--7 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.


#7

Nice Christmas Eve work, Les!

And you have left enough room to key in Prof. Kahan's weird example in Hewlett-Packard Journal (Aug-1980):

39 ENTER
40 SQRT
41 LST x
42 1
43 -
44 /
45 x<>y
46 LN
47 1/x
48 -
49 RTN

The HP-34C shows the correct answer in 23 minutes (in FIX 7): 0.03649 +/- 0.00000007 (On the 48GX in STD mode: 0.0364899739831 +/- 3.649E-13). On my HP-33C 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.


#8

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 3-point 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 3-point routine--in 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 mine--that 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

#9

Merry Christmas, Gerson (and to all!) --

Quote:
And you have left enough room to key in Prof. Kahan's weird example in Hewlett-Packard Journal (Aug-1980):

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
x-1 ln(x)

which has a nearly-vertical "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 HP-34C will return a value near 0.03662 using FIX 5 as the setting for input-function 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 HP-34C 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 HP-15C, HP-41 Advantage Pac, HP-32S, and HP-32SII will return results identical to those of the HP-34C, although at widely differing speeds of computation. The HP-42S, the HP-48 series and their derivatives, as well as the HP-33S, 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


#10

Hi Karl,

Quote:

The integral from the HP Journal article is of the function

sqrt(x) 1
f(x) = ------- - ----- between 0 and 1
x-1 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:

             2w2          w
f(w) = ---------- - ----- between 0 and 1
(w-1)(w+1) ln(w)

According to the article, the HP-34C 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
| w2 - 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 TI-51III, when just introduced to Calculus: with 32 steps only the built-in 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.


#11

Hello again, Gerson --

Quote:
            2w2          w
f(w) = ---------- - ----- between 0 and 1
(w-1)(w+1) ln(w)

According to the article, the HP-34C 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.


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 w2 = x, with 2w.dw = dx.

Strictly speaking, the recommendation in the H-P Journal article is to replace (w2 - 1) with (w - 1)(w + 1) after substituting w2 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.

#12

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 3-point routine on my 11C.

Using 50 subintervals (150 evaluations of the function), I get .036489941. Takes a generous 5 minutes and 40 seconds.

Les


#13

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.5E-12, not 3.6E-13.

Les


#14

Sorry, Les, for the typo. That's the actual 48GX estimate:

3.64899739883E-12

It's 30 minutes past midnight here but still in time to wish you a Merry Christmas in Canada :-)

Gerson.

#15

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 sped-up 15C... in 3 minutes and 26 seconds (no optimized routine)! That's what sped-up 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 speed-up factor is 1.996, 3 minutes and 26 second is exactly what we would have expected.


Edited: 24 Dec 2006, 10:32 p.m.

#16

Quote:
According to the article, the HP-34C 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 round-off 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 (w-1)(w+1) form as recommended solves the problem.

Happy New Year!

Les


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP Prime: Long integers (continued) Helge Gabert 2 786 11-07-2013, 11:24 AM
Last Post: Helge Gabert
  HP Prime: Pass "Long" Integers to a Program Helge Gabert 6 1,330 11-03-2013, 01:12 PM
Last Post: Helge Gabert
  HP Prime polynomial long division bluesun08 13 1,931 10-30-2013, 03:29 AM
Last Post: parisse
  Live from HHC2013 Namir 4 824 09-22-2013, 12:19 AM
Last Post: Jim Horn
  Battery holder for my 33C Palmer O. Hanson, Jr. 1 596 07-08-2013, 04:03 PM
Last Post: Gerson W. Barbosa
  A very long HP-17BII equation Gerson W. Barbosa 22 3,025 04-19-2013, 12:37 AM
Last Post: Gerson W. Barbosa
  A long WP-34S night Siegfried (Austria) 10 1,627 04-16-2013, 02:11 AM
Last Post: Siegfried (Austria)
  HP-25 left on for a long, long, while Matt Agajanian 12 1,768 04-10-2013, 11:33 PM
Last Post: Steve Leibson
  RPL long vs. short names peacecalc 5 1,055 10-30-2012, 01:25 PM
Last Post: peacecalc
  HP41 Long term storage Steve Hunt 6 1,041 09-27-2012, 11:08 PM
Last Post: db (martinez, ca.)

Forum Jump: