HP Forums
WP34s contributed programs. Where to send? - Printable Version

+- HP Forums (https://archived.hpcalc.org/museumforum)
+-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html)
+--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html)
+--- Thread: WP34s contributed programs. Where to send? (/thread-244016.html)



WP34s contributed programs. Where to send? - Eduardo Duenez - 05-21-2013

I wrote a WP34s program to evaluate elliptic integrals of the second kind (i.e., to find the perimeter of ellipses) using the recent beautiful recursive formulas by Semjon Adlaj:

http://www.ams.org/notices/201208/rtx120801094p.pdf

Is there any repository or official place to submit code contributions?

Eduardo


Re: WP34s contributed programs. Where to send? - Thomas Klemm - 05-21-2013

Software Wanted!


Re: WP34s contributed programs. Where to send? - Namir - 05-22-2013

Many years ago I opted to get my ow web site where I can promptly post and update articles and listings for various calculators and mathematical software. If you feel that you have much to add, then having your own web site, under your direct control, is a good approach.

Eddie Shore is another "member" of this web site who also has his own web site where he posts listings for various calculators.

Namir

Edited: 22 May 2013, 1:08 a.m.


Re: WP34s contributed programs. Where to send? - Paul Dale - 05-22-2013

Send the code to Marcus, Walter or myself and we'll add it to the library over at the sourceforge site.


- Pauli


Re: WP34s contributed programs. Where to send? - Eduardo Duenez - 05-22-2013

Yes, I have seen the announcement on the forum, but I don't see any WP34s software listed there...

For the time being I've sent the code to Pauli. Perhaps it'll make it to the Sourceforge repository eventually.

Eduardo


Re: WP34s contributed programs. Where to send? - Eduardo Duenez - 05-22-2013

Pauli,

I sent you the program via "contact user" at the forum, but "less than" signs were removed from the message by the forum software, breaking the source code...

Any alternatives?

Eduardo


Re: WP34s contributed programs. Where to send? - Paul Dale - 05-22-2013

I fixed things.

The programs are on the sourceforge site (elliptic.wp34s and modified-AGM.wp34s) in the library directory.

- Pauli


Re: WP34s contributed programs. Where to send? - Eduardo Duenez - 05-23-2013

Oh, thanks!

Eduardo


Re: WP34s contributed programs. Where to send? - Eduardo Duenez - 05-23-2013

Quote:
I fixed things.

The programs are on the sourceforge site (elliptic.wp34s and modified-AGM.wp34s) in the library directory.

- Pauli


I found one minor error, not in the code, but in the comment documenting the functionality of EI2:

Quote:
int_0^1 sqrt((1-x*t^2)/(1-t^2))dt

should read

Quote:
int_0^1 sqrt((1-x^2*t^2)/(1-t^2))dt




Re: WP34s contributed programs. Where to send? - Reth - 05-23-2013

I wish someone developed a new solver program that gives results like any HP calculator with solver.


Re: WP34s contributed programs. Where to send? - Thomas Klemm - 05-23-2013

Links for the lazy:




Re: WP34s contributed programs. Where to send? - Walter B - 05-23-2013

Reminds me on an old fairy tale starter: "At a time when wishing had helped ..."

The WP 34S software is open source. Everybody is free to write a better solver routine, e.g. you.

d:-/


Re: WP34s contributed programs. Where to send? - Thomas Klemm - 05-23-2013

This might be a good starting point:
Personal Calculator Has Key to Solve Any Equation f(x) = 0, by William M. Kahan

And then maybe debugging an emulator of either HP-34C or HP-15C. Or maybe the HP-48 series would be easier?

Kind regards

Thomas




Re: WP34s contributed programs. Where to send? - Paul Dale - 05-23-2013

I don't have time at the moment. So it isn't going to be me.

However, the source is available, the solver is nothing more than a keystroke program so go ahead and write your own and submit it or convince someone else to.


_ Pauli


Re: WP34s contributed programs. Where to send? - Thomas Klemm - 05-24-2013

Arithmetic-Geometric Mean

This is a Python-script:

eps = 1e-14
def agm(y):
x = 1.0
while math.fabs(x - y) > eps:
x, y = (x + y)/2, math.sqrt(x * y)
return x

And that's the program for the HP-42S:

00 { 29-Byte Prgm }
01 LBL "AGM"
02 4
03 X<>Y
04 1
05 LBL 00
06 RCL+ ST Y
07 X<>Y
08 RCL* ST L
09 SQRT
10 X<>Y
11 2
12 /
13 DSE ST Z
14 GTO 00

Modified Arithmetic-Geometric Mean

I've made some modifications to the formula presented in the paper. Be aware that x and y correspond to x - z and y - z.

This is a Python-script:

eps = 1e-14
def magm(y):
x, z = 1.0, 0.0
while math.fabs(x - y) > eps:
w = math.sqrt(x * y)
x, y, z = (x + y)/2 + w, 2 * w, z - w
return y + z

And that's the program for the HP-42S:

00 { 37-Byte Prgm }
01 LBL "MAGM"
02 0
03 X<>Y
04 1
05 LBL 00
06 RCL+ ST Y
07 X<>Y
08 RCL* ST L
09 SQRT
10 STO- ST Z
11 STO+ ST X
12 STO+ ST Y
13 X<>Y
14 2
15 /
16 X#Y?
17 GTO 00
18 RDN
19 +
20 END

Perimeter of Ellipse


00 { 37-Byte Prgm }
01 LBL "ELIPER"
02 STO 00
03 /
04 STO 01
05 XEQ "AGM"
06 X<> 01
07 X^2
08 XEQ "MAGM"
09 RCL/ 01
10 PI
11 *
12 STO+ ST X
13 RCL* 00
14 END

Using the example (a = 3, b = 2) that Gerson used in
this other thread I got:

2 ENTER 3
XEQ "ELIPER"

15.8654395897

Kind regards

Thomas

Edited: 24 May 2013, 7:40 p.m.


Re: WP34s contributed programs. Where to send? - Gerson W. Barbosa - 05-24-2013

Very nice! Accurate to at least 11 digits and less than 2 seconds on my HP-42S. This runs "ellipses" around the approximation I've used :-)

Cheers,

Gerson.


Re: WP34s contributed programs. Where to send? - Thomas Klemm - 05-24-2013

This is Hugh Steers' C-program written in Python:

def agm(b):
a, t, s = 1, 1, (1 + b*b)/2
while a != b:
c = (a - b)/2
s -= c*c*t
t *= 2
a, b = (a + b)/2, math.sqrt(a*b)
return s/a

def ellipse(a, b):
return 2 * math.pi * a * agm(b / a)

And that's the program for the HP-42S:

00 { 71-Byte Prgm }      
01 LBL "Ellipse" 21 /
02 STO 0 22 RCL 02
03 / 23 RCL* ST Z
04 1 24 SQRT
05 X<>Y 25 RDN
06 ENTER 26 RDN
07 X^2 27 RCL- ST T
08 RCL+ ST Z 28 X^2
09 2 29 RCL* 01
10 STO 03 30 -
11 STO+ 03 31 X<> ST Z
12 1/X 32 DSE 03
13 STO 01 33 GTO 00
14 * 34 RDN
15 X<> ST Z 35 /
16 LBL 00 36 RCL* 00
17 STO 02 37 PI
18 RCL + ST Y 38 *
19 2 39 STO+ ST X
20 STO* 01 40 END

Using your example (a = 3, b = 2) I got:

2 ENTER 3
XEQ "Ellipse"

15.8654395892

Kind regards

Thomas




Re: WP34s contributed programs. Where to send? - Gerson W. Barbosa - 05-25-2013

Thanks for converting Hugh's program!
For a=1 and b=0.1 it returns

4.06397418824
The previous program returns
4.0639741721
The exact result is
4.06397418010089574.. 
There's loss of accuracy as b gets smaller, but the results are consistently greater than 4, unlike the previous program which will eventually give results below 4. IMO, this should go to the Software Library.

Cheers,

Gerson.


Re: WP34s contributed programs. Where to send? - Eduardo Duenez - 05-29-2013

Quote:

Arithmetic-Geometric Mean

This is a Python-script:

eps = 1e-14
def agm(y):
x = 1.0
while math.fabs(x - y) > eps:
x, y = (x + y)/2, math.sqrt(x * y)
return x

Note that the WP34s modified-AGM algorithm I wrote, although longer, is very robust. It keeps track of the last two iterations. As soon as one iteration fails to improve over the last one, the last one is returned. If at any point xn = yn, it also stops. The algorithm is guaranteed to stop, and to return the best possible answer limited only by the calculator's working precision. (I.e., no need to set up any small acceptable tolerance eps beforehand.) The algorithm is oblivious to the WP34s being in standard or double precision. The price paid is the use of 7 local registers (though only 2 stack levels).

Eduardo


Re: WP34s contributed programs. Where to send? - Walter B - 05-29-2013

Quote:
The price paid is the use of 7 local registers (though only 2 stack levels).
Little is as cheap as local registers. Each and every program may use up to 144 of them (as long as memory is sufficient).

d:-)


Re: WP34s contributed programs. Where to send? - Paul Dale - 05-29-2013

In XROM code, even the memory part is free :)

Not that I'm planning on including the elliptic integral code as a built in function.


- Pauli


Re: WP34s contributed programs. Where to send? - Gerson W. Barbosa - 05-29-2013

Is there still room for complex Zeta? Thanks!

Gerson.


Complex zeta function - Paul Dale - 05-29-2013

There possibly is room, the problem is finding an algorithm that converges in reasonable time. I don't know of any.

For real zeta, I'm using Borwein's algorithm that converges rapidly & was state of the art when I implement this on the 34S.

Unfortunately, for the complex number a + i b, this algorithm requires a number of steps proportional to |b|. I originally coded this in C using a fixed number of steps. The results were great for complex numbers near the real line but went wrong as they got further away. Adapting the code to use a variable number of terms would have been possible & the number of terms required can be calculated before starting the computation. However, we'd then have had a function whose run time was open ended. Calculating zeta(3 + 1e23 i) would have taken more than a while -- Wolfram alpha fails to calculate this for me. Try out the HP 41 program to see how long it takes :-)

The might be some hope using the Lerch transcendent for which some series are known but I've not investigated their convergence characteristics. Implementing this function would buy us a pile of others for free so it would be my preferred route.


- Pauli


Re: WP34s contributed programs. Where to send? - Thomas Klemm - 05-29-2013

Why not using the 8-level stack?

Cheers

Thomas


Re: WP34s contributed programs. Where to send? - Thomas Klemm - 05-29-2013

Did you notice that my HP-42S program "AGM" uses a fixed amount of 4 loops? Very robust and guaranteed to stop. But I admit that I was a little sloppy with "MAGM": it might be possible that it doesn't terminate. But with the WP-34S we could use x~? Y.

Quote:
x~? will be true if the rounded values of x and a are equal (see ROUND).

Cheers

Thomas


Re: Complex zeta function - Thomas Klemm - 05-29-2013

I could imagine that the problem to calculate zeta for complex values is related to the Zeta function universality:

Quote:
In mathematics, the universality of zeta-functions is the remarkable ability of the Riemann zeta-function and other, similar, functions, such as the Dirichlet L-functions, to approximate arbitrary non-vanishing holomorphic functions arbitrarily well.

It isn't premised that the holomorphic function converges fast on U. Therefore I assume zeta converges there at about the same speed.

Kind regards

Thomas


Re: WP34s contributed programs. Where to send? - Eduardo Duenez - 05-29-2013

Quote:
But I admit that I was a little sloppy with "MAGM": it might be possible that it doesn't terminate. But with the WP-34S we could use x~? Y.

Thomas


While one would (optimistically) expect the sequences xn and yn to pass the test x~? Y after a few iterations, an algorithm relying on that test is not guaranteed to terminate. My "trial runs" debugging the WP34s code showed that errors of at least 1 ULP between xn and yn persist in some cases. Also, x~? Y would invalidate any extended precision digits (further left of those displayed in the calculator's current setting). I like (expect, actually) my calculators to compute expressions like exp(1) to the internal full precision of the calculator even while working, say, in SCI 3 mode. Of course there are reasons to perform expensive calculations such as SOLVE and INTEGRATE to lesser precision, but that is not the case here.

In my implementation of the modified-AGM, I tried to achieve maximum precision, whether the digits be visible or not. Since the algorithm converges very quickly (quadratically), I think aiming for maximum precision is the right thing to do. The only relative drawbacks are sightly longer code and a few more local registers used.

Eduardo

Edited: 29 May 2013, 9:27 p.m.


Re: WP34s contributed programs. Where to send? - Thomas Klemm - 05-30-2013

You might still be able to reduce the size of the program using the simpler formulas after the transformation:

Quote:
Be aware that x and y correspond to x - z and y - z.
        w = math.sqrt(x * y)
x, y, z = (x + y)/2 + w, 2 * w, z - w

Cheers

Thomas


Re: WP34s contributed programs. Where to send? - Eduardo Duenez - 05-30-2013

Quote:
You might still be able to reduce the size of the program using the simpler formulas after the transformation:


Cheers

Thomas


Yes, this change of variables may shorten the program. I might implement some changes in the WP34s code later. (My original goal in implementing the algorithm was maximal accuracy, not code optimization. I wasn't too worried about speed due to the inherent efficiency of the algorithm.)

Thanks,
Eduardo