The Museum of HP Calculators


Weierstrass Elliptic Functions for the HP-41

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