WP34s contributed programs. Where to send?



#30

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


#31

Software Wanted!


#32

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

#33

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.

#34

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


- Pauli


#35

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


#36

I fixed things.

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

- Pauli


#37

Oh, thanks!

Eduardo

#38

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

#39

Links for the lazy:

#40

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.


#41

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.


#42

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


#43

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.

#44

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


#45

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:-)


#46

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


#47

Is there still room for complex Zeta? Thanks!

Gerson.


#48

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


#49

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

#50

Why not using the 8-level stack?

Cheers

Thomas

#51

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


#52

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.


#53

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


#54

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

#55

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


#56

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:-/


#57

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

#58

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


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP Prime: do not send programs and shuts down by pressing Apps Davi Ribeiro de Oliveira 1 446 11-12-2013, 11:05 AM
Last Post: Joseph Ec
  Hp PRIME - how to send a list to the connectivity Kit giancarlo 1 482 11-10-2013, 11:50 AM
Last Post: Tim Wessman
  What does the Send menu key in the Apps window do? John Colvin 3 532 10-23-2013, 10:56 AM
Last Post: Tim Wessman
  How to save WP34s programs to a file? Harald 9 1,096 04-26-2013, 01:44 PM
Last Post: Walter B
  [WP34S] WP34S firmware on the AT91SAM7L-STK dev kit? jerome ibanes 1 427 10-04-2012, 04:59 PM
Last Post: Paul Dale
  [wp34s] Incomplete Gamma on the wp34s Les Wright 18 1,669 12-06-2011, 11:07 AM
Last Post: Namir
  [wp34s] Romberg Integration on the wp34s Les Wright 2 538 12-04-2011, 10:49 PM
Last Post: Les Wright
  pleas help me to send "txt" file to 48sx by usb irda from the computer samiroch 9 815 05-28-2010, 07:55 PM
Last Post: Vieira, Luiz C. (Brazil)
  HP User Contributed Programs Bill (Smithville, NJ) 0 263 03-15-2008, 07:06 PM
Last Post: Bill (Smithville, NJ)
  Plea for HP-41 User Contributed Programs Bill (Smithville, NJ) 11 926 03-02-2008, 03:26 PM
Last Post: Bill (Smithville, NJ)

Forum Jump: