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

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. *

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

- Pauli

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

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

I fixed things.

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

- Pauli

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

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

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

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

## 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. *

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.

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

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.

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

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

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

Is there still room for complex Zeta? Thanks!

Gerson.

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

Why not using the 8-level stack?

Cheers

Thomas

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

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

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. *

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

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