The Museum of HP Calculators
The Museum of HP Calculators
This program is Copyright © 2006 by Jean-Marc Baillard and is used here by permission.
This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.
Overview
1°) A Laurent Series ( real arguments )
2°) Duplication formula ( real arguments
)
3°) A Laurent Series ( complex arguments
)
4°) Duplication formula ( complex arguments
)
5°) Weierstrass & Jacobian Elliptic
Functions - Half Periods
1°) A Laurent Series ( real arguments )
-The Weierstrass Elliptic Function P(x;g2;g3)
satisfies the differential equation: (P')2
= 4.P3 - g2.P - g3
-This program calculates P(x;g2;g3)
by a Laurent series:
P(x;g2;g3) = x-2 + c2.x2 + c3.x4 + ...... + cn.x2n-2 + ....
where c2 = g2/20 ; c3 = g3/28 and cn = 3 ( c2. cn-2 + c3. cn-3 + ....... + cn-2. c2 ) / (( 2n+1 )( n-3 )) ( n > 3 )
-The successive cn are computed and stored into registers R05 R06 ......
Data Registers: R00 to R04 are used for temporary
data storage and R05 = c2 ; R06 = c3 ; R07
= c4 ; ......
Flags: /
Subroutines: /
01 LBL "WEF"
02 X^2
03 STO 02
04 CLX
05 20
06 /
07 STO 05
08 X<>Y
09 28
10 /
11 STO 06
12 RCL 02
13 *
14 +
15 STO 01
16 5
17 STO 00
18 LBL 01
19 2
20 STO 04
21 LBL 02
22 RCL 00
23 .1
24 %
25 5
26 +
27 STO 03
28 CLX
29 LBL 03
30 RCL IND Y
31 RCL IND 03
32 *
33 +
34 DSE Y
35 ISG 03
36 GTO 03
37 3
38 *
39 RCL 00
40 ENTER^
41 ST+ Y
42 DSE Y
43 4
44 -
45 *
46 /
47 ISG 03
48 CLX
49 STO IND 03
50 RCL 02
51 RCL 00
52 3
53 -
54 Y^X
55 *
56 RCL 01
57 +
58 STO 01
59 ISG 00
60 CLX
61 LASTX
62 X#Y?
63 GTO 01
64 DSE 04
65 GTO 02
66 RCL 02
67 ST* Y
68 1/X
69 +
70 END
( 95 bytes / SIZE ? )
| STACK | INPUTS | OUTPUTS |
| Z | g3 | / |
| Y | g2 | / |
| X | x | P(x;g2;g3) |
Example: Calculate P(x;g2;g3) for x = 0.6 g2 = 0.9 g3 = 1.4
1.4 ENTER^
0.9 ENTER^
0.6 XEQ "WEF" >>>> 2.800500784
( in 32 seconds )
-The SIZE must be ( at least ) 016 in this example.
-If you get the message "NON EXISTENT" ( line 49 ) , increase the SIZE
and R/S
-The Laurent series may diverge if x is too great. In this case, one
may use the duplication formula below.
2°) Duplication formula ( real arguments )
-We have P(2x) = -2 P(x) + ( 6 P2(x)
- g2/2 )2 / ( 4 ( 4 P3(x) - g2
P(x) - g3 ) )
-This formula is used recursively ( n times ) to obtain P(2nx)
if we know P(x) ( n = 1 ; 2 ; 3 ;
..... )
Data Registers: • R00 = n ( n
> 0 )
Register R00 is to be initialized before executing "DUPW"
R01: temp R02 = g2 R03 = g3
Flags: /
Subroutine: /
01 LBL "DUPW"
02 RDN
03 STO 02
04 X<>Y
05 STO 03
06 R^
07 LBL 01
08 STO 01
09 X^2
10 6
11 *
12 RCL 02
13 2
14 /
15 -
16 X^2
17 RCL 01
18 ST+ X
19 X^2
20 RCL 02
21 -
22 RCL 01
23 *
24 RCL 03
25 -
26 4
27 *
28 /
29 RCL 01
30 ST+ X
31 -
32 DSE 00
33 GTO 01
34 END
( 47 bytes / SIZE 004 )
| STACK | INPUTS | OUTPUTS |
| Z | g3 | / |
| Y | g2 | / |
| X | P(x;g2;g3) | P(2nx;g2;g3) |
Example: Calculate P(x;g2;g3) for x = 4.8 g2 = 0.9 g3 = 1.4
-The Laurent series diverges for these arguments but 4.8 = 23* 0.6 and we have already found that P(0.6;0.9;1.4) = 2.800500784 , thus:
3 STO 00
1.4 ENTER^
0.9 ENTER^
2.800500784 XEQ "DUPW" >>>> P(4.8;0.9;1.4)
= 1.954704160 ( in 3seconds )
3°) A Laurent Series ( complex arguments )
-The same Laurent series is now applied to z = x + i.y ( but g2 and g3 real )
Data Registers: R00 to R06 are used for temporary
data storage and R07 = c2 ; R08 = c3 ; R09
= c4 ; ......
Flags: /
Subroutines: /
001 LBL "WEFZ"
002 R-P
003 X^2
004 STO 03
005 RDN
006 ST+ X
007 STO 04
008 CLX
009 20
010 /
011 STO 07
012 CLX
013 28
014 /
015 STO 08
016 RCL 03
017 *
018 RCL 04
019 X<>Y
020 P-R
021 RCL 07
022 +
023 STO 01
024 X<>Y
025 STO 02
026 7
027 STO 00
028 LBL 01
029 2
030 STO 06
031 LBL 02
032 RCL 00
033 .1
034 %
035 7
036 +
037 STO 05
038 CLX
039 LBL 03
040 RCL IND Y
041 RCL IND 05
042 *
043 +
044 DSE Y
045 ISG 05
046 GTO 03
047 3
048 *
049 RCL 00
050 ENTER^
051 ST+ Y
052 DSE X
053 5
054 ST- Z
055 -
056 *
057 /
058 ISG 05
059 CLX
060 STO IND 05
061 RCL 04
062 RCL 03
063 RCL 00
064 5
065 -
066 ST* Z
067 Y^X
068 RCL IND 05
069 *
070 P-R
071 RCL 01
072 +
073 STO 01
074 LASTX
075 -
076 X<>Y
077 RCL 02
078 +
079 STO 02
080 LASTX
081 -
082 R-P
083 ISG 00
084 CLX
085 X#0?
086 GTO 01
087 RCL 02
088 RCL 01
089 R-P
090 RCL 03
091 *
092 X<>Y
093 RCL 04
094 +
095 X<>Y
096 P-R
097 RCL 04
098 CHS
099 RCL 03
100 1/X
101 P-R
102 RDN
103 ST+ Z
104 X<> T
105 +
106 END
( 135 bytes / SIZE? )
| STACK | INPUTS | OUTPUTS |
| T | g3 | / |
| Z | g2 | / |
| Y | y | B |
| X | x | A |
with P(x+iy;g2;g3) = A + i.B
Example: Calculate P(z;g2;g3) for z = 0.6 + 0.4 i g2 = 0.9 g3 = 1.4
1.4 ENTER^
0.9 ENTER^
0.4 ENTER^
0.6 XEQ "WEFZ" >>>> 0.7390436447
X<>Y -1.744031391 ( in 48seconds
; SIZE 018 ( at least ) )
Whence P(0.6+0.4i ;0.9;1.4)
= 0.7390436447 -1.744031391 i
4°) Duplication formula ( complex arguments )
-The same duplication formula holds for complex arguments.
Data Registers: • R00 = n ( n
> 0 )
Register R00 is to be initialized before executing "DUPWZ"
R01: temp R02 = g2 R03 = g3
R04 to R06: temp
Flags: /
Subroutines: /
01 LBL "DUPWZ"
02 R^
03 STO 03
04 R^
05 STO 02
06 R^
07 R^
08 LBL 01
09 STO 01
10 X<>Y
11 STO 04
12 X<>Y
13 R-P
14 3
15 ST* Z
16 Y^X
17 4
18 *
19 P-R
20 RCL 04
21 RCL 02
22 *
23 ST- Z
24 CLX
25 RCL 01
26 LASTX
27 *
28 -
29 RCL 03
30 -
31 R-P
32 STO 05
33 X<>Y
34 STO 06
35 RCL 04
36 RCL 01
37 R-P
38 2
39 ST* Z
40 Y^X
41 6
42 *
43 P-R
44 RCL 02
45 2
46 /
47 -
48 R-P
49 2
50 ST/ Y
51 ST* Z
52 Y^X
53 RCL 06
54 ST- Z
55 CLX
56 RCL 05
57 /
58 P-R
59 RCL 04
60 ST+ X
61 ST- Z
62 CLX
63 RCL 01
64 ST+ X
65 -
66 DSE 00
67 GTO 01
68 END
( 89 bytes / SIZE 007 )
| STACK | INPUTS | OUTPUTS |
| T | g3 | / |
| Z | g2 | / |
| Y | B | B' |
| X | A | A' |
with P(z;g2;g3)
= A + i.B and P(2nz;g2;g3)
= A' + i.B'
Example: Deduce P(4.8+3.2
i ;0.9;1.4) from P(0.6+0.4i ;0.9;1.4)
= 0.7390436447 -1.744031391 i
3 STO 00
1.4 ENTER^
0.9 ENTER^
0.4 ENTER^
0.6 XEQ "DUPWZ" >>>> 0.152165025
X<>Y -1.254979593 ( in 22 seconds
)
whence P(4.8+3.2 i ;0.9;1.4)
= 0.152165025 - 1.254979593 i
5°) Weierstrass & Jacobian Elliptic Functions
- Half Periods
-The Weierstrass Elliptic Functions may also be computed by mean of
the Jacobian Elliptic Functions.
-This program also gives the first derivative and the half-periods
of P(x;g2;g3)
-We have P(x+2k*Omega+2ik'*Omega') = P(x)
for all x ( except the poles ) if k
and k' are integers & Omega and i.Omega' are the half-periods.
Formulae: Let r1
; r2 ; r3 the 3 roots of
the polynomial 4x3- g2.x -g3
and m the parameter of the complete elliptic integrals
of the 1st kind: K(m)
IF THE 3 ROOTS ARE REAL ( R1 > R2 > R3 ) IF ONLY ONE ROOT IS REAL ( SAY R1 )
m = ( r2 - r3 )/( r1 - r3
)
m = 1/2 - 3 r1/(4H)
Omega = K(m)/( r1 - r3 )1/2
Omega2 = K(m)/H1/2
Omega' = K(1-m)/( r1 - r3 )1/2
Omega'2 = K(1-m)/H1/2
P(x;g2;g3)
= r3 + ( r1 - r3 )/sn2( y |
m )
P(x;g2;g3) = r1 + H ( 1 + cn(y'|m)
) / ( 1 - cn(y'|m) )
P'(x;g2;g3)
= -2 ( r1 - r3 )3/2 cn (y|m) dn(y|m) /
sn3(y|m)
P'(x;g2;g3) = -4 H3/2 sn(y'|m)
dn(y'|m) / ( 1 - cn(y'|m) )2
with y = ( r1 - r3 )1/2.x
with y' = 2.x.H1/2 and
H = ( 2.r12 + r2.r3 )1/2
Omega2 = Omega + i. Omega'
Omega'2 = Omega' + i. Omega
Data Registers: R00 to R07 are used for temporary
data storage by the Jacobian Elliptic Function and Complete Elliptic Integral
programs
and when the program stops: R08 = x
; R09 = Omega if flag F01 is clear R09 = Omega2
if flag F01 is set R11: temp
R10 = Omega2 ----------------- R10 = Omega'2 ---------------
R12: temp
Flags: F01 - F02 - F05
Subroutines: "CEI" ( Complete Ellptic
Integrals ) see "Jacobian Elliptic Functions for the HP-41"
"JEF" ( Jacobian Elliptic Functions ) ---------------------------------------------
"P3" ( Cubic Equations ) see "Polynomials for the HP-41" or another
program which returns the larger real solution of a cubic in X-register.
Note: If you have used registers R08 & R09 instead
of synthetic registers M & N in the "JEF" program,
replace
R08 and R09 by R13 and R14 in the following listing.
001 LBL "WEF2"
002 STO 08
003 RDN
004 CHS
005 0
006 X<> Z
007 CHS
008 4
009 RDN
010 XEQ "P3"
011 FS? 01
012 GTO 02
013 RCL Z
014 STO 11
015 ST- Z
016 -
017 STO 12
018 XEQ 01
019 ST/ Y
020 X^2
021 STO 01
022 /
023 *
024 RCL 01
025 1/X
026 GTO 03
027 LBL 01
028 CF 05
029 /
030 STO 10
031 1
032 X<>Y
033 X=Y?
034 SF 05
035 FC? 05
036 XEQ "CEI"
037 X<>Y
038 FS?C 05
039 E90
040 STO 09
041 1
042 STO Y
043 RCL 10
044 -
045 X=Y?
046 SF 05
047 FC? 05
048 XEQ "CEI"
049 X<>Y
050 FS?C 05
051 E90
052 X<> 10
053 RCL 12
054 SQRT
055 ST/ 09
056 ST/ 10
057 RCL 08
058 FS? 01
059 ST+ X
060 *
061 XEQ "JEF"
062 RTN
063 LBL 02
064 STO 11
065 X^2
066 ST+ X
067 RDN
068 R-P
069 X^2
070 R^
071 +
072 SQRT
073 STO 12
074 RCL 11
075 3
076 *
077 X<>Y
078 /
079 2
080 -
081 CHS
082 4
083 XEQ 01
084 SF 01
085 X<>Y
086 STO 01
087 1
088 -
089 STO 02
090 X^2
091 /
092 *
093 ST+ X
094 1
095 RCL 01
096 +
097 RCL 02
098 CHS
099 /
100 LBL 03
101 RCL 12
102 SQRT
103 ST+ X
104 CHS
105 ST* Z
106 X<> L
107 ST* Z
108 *
109 RCL 11
110 +
111 END
( 170 bytes / SIZE 013 )
| STACK | INPUTS | OUTPUTS |
| Z | g3 | / |
| Y | g2 | P'(x;g2;g3) |
| X | x | P(x;g2;g3) |
Example1: Calculate P(x;g2;g3) & P'(x;g2;g3) for x = 2 g2 = 4 g3 = 1
1 ENTER^
4 ENTER^
2 XEQ "WEF2" >>>> P(2;4;1)
= 4.950267724 X<>Y
P'(2;4;1)
= 21.55057197 ( in 28 seconds )
-We have R09 = 1.225694692 & R10 = 1.496729323 ( Omega & Omega' because F01 is clear )
Therefore the primitive half-periods
are: 1.225694692 & 1.496729323
i
Example2: Calculate P(x;g2;g3) & P'(x;g2;g3) for x = 1 g2 = 2 g3 = 3
3 ENTER^
2 ENTER^
1 XEQ "WEF2" >>>> P(1;2;3)
= 1.214433709 X<>Y
P'(1;2;3) = -1.317406193 ( in 27
seconds )
-We have R09 = 1.197220889 & R10 = 2.350281226 ( Omega2 & Omega'2 because F01 is set )
Whence:
Omega = 0.598610445 - 1.175140613 i &
Omega' = 0.598610445 + 1.175140613 i
Notes:
-I've called Omega' & Omega'2 what Abramowitz and Stegun call
Omega'/i & (Omega'2)/i
-This program doesn't work if g2 = g3 =
0 but in this case, P(x) = 1/x2
-Lines 28-32 -33-34-35-38-39-42-45-46-47-50-51 are useful only
if 2 roots of the polynomial 4x3- g2.x -g3
are equal:
in this case, one of the half-periods is infinite and a number
near E90 is stored in register R09 or R10.
Without these lines, there would be a DATA ERROR at line 38
of the "CEI" program ( K(1) = infinite )
For instance P(1;48;-64) = 2.181597704 ;
P'(1;48;-64)
= -0.903006185 ; R09 = 4.082 1089 = infinite =
Omega ; R10 = 0.641274915 = Omega' ( CF01 )
P(1;12;8)
= 2.079381536 ; P'(1;12;8) = 1.735214810 ; R09 = 0.906899682
= Omega ; R10 = 5.7735 1089 = infinite = Omega'
( CF01 )
-Numerous other relations with Jacobian Elliptic Functions and Theta
Functions may be found in the reference below
Reference:
Abramowitz and Stegun , "Handbook of Mathematical Functions" - Dover
Publications - ISBN 0-486-61272-4
Go back to the HP-41 software library
Go back to the general software library
Go
back to the main exhibit hall