Weierstrass

# Weierstrass Elliptic Functions for the HP-41

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

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

 01  LBL "WEFZ"   02  R-P   03  X^2   04  STO 03   05  RDN   06  ST+ X   07  STO 04   08  CLX   09  20   10  /   11  STO 07   12  CLX   13  28   14  /   15  STO 08   16  RCL 03   17  *   18  RCL 04   19  X<>Y   20  P-R   21  RCL 07   22  + 23  STO 01   24  X<>Y   25  STO 02   26  7   27  STO 00   28  LBL 01   29  2   30  STO 06   31  LBL 02   32  RCL 00   33  .1   34  %   35  7   36  +   37  STO 05   38  CLX   39  LBL 03   40  RCL IND Y   41  RCL IND 05   42  *   43  +   44  DSE Y 45  ISG 05   46  GTO 03   47  3   48  *   49  RCL 00   50  ENTER^   51  ST+ Y   52  DSE X   53  5   54  ST- Z   55  -   56  *   57  /   58  ISG 05   59  CLX   60  STO IND 05   61  RCL 04   62  RCL 03   63  RCL 00   64  5   65  -   66  ST* Z 67  Y^X   68  RCL IND 05   69  *   70  P-R   71  RCL 01   72  +   73  STO 01   74  LASTX   75  -   76  X<>Y   77  RCL 02   78  +   79  STO 02   80  LASTX   81  -   82  R-P   83  ISG 00   84  CLX   85  X#0?   86  GTO 01   87  RCL 02   88  RCL 01 89  R-P   90  RCL 03   91  *   92  X<>Y   93  RCL 04          94  +   95  X<>Y   96  P-R   97  RCL 04   98  CHS   99  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 if  flag F01 is clear ,  R10 = Omega'2  if flag F01 is set  ,   R12: temp

Flags: F01 - F02 - F05

Subroutines:

"CEI"  ( Complete Ellptic Integrals ) see "Jacobian Elliptic Functions for the HP-41"
"JEF" ( Jacobian Elliptic Functions )  see "Jacobian Elliptic Functions for the HP-41"
"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.

 01  LBL "WEF2"   02  STO 08   03  RDN   04  CHS   05  0   06  X<> Z   07  CHS   08  4   09  RDN   10  XEQ "P3"   11  FS? 01   12  GTO 02   13  RCL Z   14  STO 11   15  ST- Z   16  -   17  STO 12   18  XEQ 01   19  ST/ Y   20  X^2   21  STO 01   22  /   23  * 24  RCL 01   25  1/X   26  GTO 03   27  LBL 01   28  CF 05   29  /   30  STO 10   31  1   32  X<>Y   33  X=Y?   34  SF 05   35  FC? 05   36  XEQ "CEI"     37  X<>Y   38  FS?C 05   39   E90   40  STO 09   41  1   42  STO Y   43  RCL 10   44  -   45  X=Y?   46  SF 05 47  FC? 05   48  XEQ "CEI"   49  X<>Y   50  FS?C 05   51   E90   52  X<> 10   53  RCL 12   54  SQRT   55  ST/ 09   56  ST/ 10   57  RCL 08   58  FS? 01   59  ST+ X   60  *   61  XEQ "JEF"     62  RTN   63  LBL 02   64  STO 11   65  X^2   66  ST+ X   67  RDN   68  R-P   69  X^2 70  R^   71  +   72  SQRT   73  STO 12   74  RCL 11   75  3   76  *   77  X<>Y   78  /   79  2   80  -   81  CHS   82  4   83  XEQ 01          84  SF 01   85  X<>Y   86  STO 01   87  1   88  -   89  STO 02   90  X^2   91  /   92  * 93  ST+ X   94  1   95  RCL 01   96  +   97  RCL 02   98  CHS   99  / 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:

[1] Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4