▼
Posts: 56
Threads: 10
Joined: May 2006
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
▼
Posts: 735
Threads: 34
Joined: May 2007
▼
Posts: 56
Threads: 10
Joined: May 2006
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
Send the code to Marcus, Walter or myself and we'll add it to the library over at the sourceforge site.
 Pauli
▼
Posts: 56
Threads: 10
Joined: May 2006
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
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
I fixed things.
The programs are on the sourceforge site (elliptic.wp34s and modifiedAGM.wp34s) in the library directory.
 Pauli
▼
Posts: 56
Threads: 10
Joined: May 2006
Posts: 56
Threads: 10
Joined: May 2006
Quote:
I fixed things.
The programs are on the sourceforge site (elliptic.wp34s and modifiedAGM.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((1x*t^2)/(1t^2))dt
should read
Quote:
int_0^1 sqrt((1x^2*t^2)/(1t^2))dt
Posts: 735
Threads: 34
Joined: May 2007
Posts: 735
Threads: 34
Joined: May 2007
ArithmeticGeometric Mean
This is a Pythonscript:
eps = 1e14
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 HP42S:
00 { 29Byte 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 ArithmeticGeometric 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 Pythonscript:
eps = 1e14
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 HP42S:
00 { 37Byte 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 { 37Byte 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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Very nice! Accurate to at least 11 digits and less than 2 seconds on my HP42S. This runs "ellipses" around the approximation I've used :)
Cheers,
Gerson.
▼
Posts: 735
Threads: 34
Joined: May 2007
This is Hugh Steers' Cprogram 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 HP42S:
00 { 71Byte 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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 56
Threads: 10
Joined: May 2006
Quote:
ArithmeticGeometric Mean
This is a Pythonscript:
eps = 1e14
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 modifiedAGM 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
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
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:)
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
Is there still room for complex Zeta? Thanks!
Gerson.
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
▼
Posts: 735
Threads: 34
Joined: May 2007
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 zetafunctions is the remarkable ability of the Riemann zetafunction and other, similar, functions, such as the Dirichlet Lfunctions, to approximate arbitrary nonvanishing 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
Posts: 735
Threads: 34
Joined: May 2007
Why not using the 8level stack?
Cheers
Thomas
Posts: 735
Threads: 34
Joined: May 2007
Did you notice that my HP42S 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 WP34S we could use x~? Y.
Quote:
x~? will be true if the rounded values of x and a are equal (see ROUND).
Cheers
Thomas
▼
Posts: 56
Threads: 10
Joined: May 2006
Quote:
But I admit that I was a little sloppy with "MAGM": it might be possible that it doesn't terminate. But with the WP34S 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 modifiedAGM, 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.
▼
Posts: 735
Threads: 34
Joined: May 2007
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
▼
Posts: 56
Threads: 10
Joined: May 2006
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
Posts: 556
Threads: 9
Joined: Jul 2007
I wish someone developed a new solver program that gives results like any HP calculator with solver.
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
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:/
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
