Overview
1°) Via the Arithmetic-Geometric Mean
2°) Using Gauss' Transformations
a) Program#1
b) M-Code Routine ( 0 <= m < 1
)
c) Program#2
3°) Jacobian Elliptic Function Sn ( x | m ) 0 <= m < 1
a) Focal Program
b) M-Code Routine
1°) Via the Arithmetic-Geometric Mean
-There are 12 Jacobian elliptic functions.They can be defined like this:
Let m be a number that verifies: 0 < = m
< = 1 ( m is called the parameter )
if u = §0v
( 1 - m sin2 t ) -1/2 dt ( 0 <
= v < = 90° )
the angle v = am u
is called the amplitude
and the 3 Jacobian elliptic functions sn , cn ,
dn are defined by:
sn ( u | m ) = sin v , cn ( u | m) = cos v , dn ( u | m ) = ( 1 - m sin2 v )1/2
-The 9 other Jacobian elliptic functions can be obtained with the relations:
cd = cn / dn ; dc = dn / cn ;
ns = 1 / sn
sd = sn / dn ; nc = 1 / cn
; ds = dn / sn
nd = 1 / dn ; sc = sn /
cn ; cs = cn / sn
The Arithmetic-Geometric
Mean:
-The following program uses the process of the Arithmetic-Geometric Mean (AGM):
Starting with a0 = 1 ; b0 = (1-m)1/2
; c0 = m1/2
we determine a1
, b1 , c1 ; ....... ; an , bn
, cn with ak= (ak-1+
bk-1)/2 ; bk= (ak-1.bk-1)1/2
; ck = (ak-1- bk-1)/2
and we stop when cn = 0 (
in fact, when an-1= an )
-This program calculates only the 3 copolar trio: sn ( u | m )
; cn ( u | m ) ; dn ( u | m ).
-From these, all the 9 other elliptic functions can be determined by
the relations mentioned above.
-When the AGM process stops, the HP-41 calculates vn
= 2nanu*180/pi ( in degrees )
and then, vn-1, ... , v0 from
the recurrence relation: sin (2vk-1-vk)
= ( ck/ak ) sin vk
then, sn ( u | m ) = sin v0 ;
cn ( u | m ) = cos v0 and dn ( u | m
) = ( 1 - m sin2 v0 )1/2
-If m > 1, the program uses the relations:
sn ( u | m ) = sn ( u m1/2
| 1/m ) / m1/2
cn ( u | m ) = dn ( u m1/2
| 1/m )
dn ( u | m ) = cn ( u m1/2
| 1/m )
-If m = 1: sn ( u | m ) = tanh u ; cn ( u | m ) = dn ( u | m ) = sech u
-If m < 0:
sn ( u | m ) = sd ( u (1-m)1/2
| -m/(1-m) ) / (1-m)1/2
cn ( u | m ) = cd ( u (1-m)1/2
| -m/(1-m) )
dn ( u | m ) = nd ( u (1-m)1/2
| -m/(1-m) )
Data Registers: R00 = 1 , 2 ,
... , n and then n-1, ... , 1 , 0 ;
R01 = c1/a1 , ..... , Rnn = cn/an
Flags: F01 & F02
Subroutines: /
01 LBL "JEF"
02 DEG 03 CF 01 04 CF 02 05 X<>Y 06 X#0? 07 X>0? 08 GTO 02 09 CHS 10 STO Z 11 1 12 + 13 ST/ Z 14 SQRT 15 * 16 X<>Y 17 SF 01 18 LBL 02 19 1 20 X>Y? 21 GTO 04 22 X#Y? 23 GTO 03 24 X<> Z 25 ENTER^ 26 E^X 27 LASTX 28 CHS |
29 E^X
30 + 31 2 32 X<>Y 33 / 34 X<>Y 35 ST+ X 36 E^X-1 37 RCL X 38 2 39 + 40 / 41 GTO 06 42 LBL 03 43 RDN 44 STO Z 45 SQRT 46 * 47 X<>Y 48 1/X 49 1 50 SF 02 51 LBL 04 52 STO 00 53 X<> Z 54 STO M 55 RDN 56 STO N |
57 -
58 SQRT 59 SIGN 60 LASTX 61 LBL 01 62 RCL Y 63 RCL Y 64 - 65 STO IND 00 66 X<> T 67 RCL Y 68 + 69 ST/ IND 00 70 2 71 / 72 R^ 73 X=Y? 74 GTO 02 75 RCL Z 76 * 77 SQRT 78 ISG 00 79 CLX 80 GTO 01 81 LBL 02 82 2 83 RCL 00 84 Y^X |
85 *
86 RCL M 87 * 88 R-D 89 LBL 10 90 ENTER^ 91 SIN 92 RCL IND 00 93 * 94 ASIN 95 + 96 2 97 / 98 DSE 00 99 GTO 10 100 COS 101 LAST X 102 SIN 103 ENTER^ 104 X^2 105 RCL N 106 CHS 107 * 108 1 109 + 110 SQRT 111 STO T 112 X<> M |
113 SIGN
114 X<> N 115 RDN 116 FS?C 01 117 GTO 05 118 FC?C 02 119 GTO 06 120 X<>Y 121 X<> T 122 SQRT 123 * 124 GTO 06 125 LBL 05 126 R^ 127 R^ 128 ST/ T 129 ST/ Z 130 1/X 131 RDN 132 SIGN 133 ST- L 134 X<> L 135 CHS 136 SQRT 137 * 138 LBL 06 139 CLA 140 END |
( 193 bytes / SIZE 008 )
STACK | INPUTS | OUTPUTS |
Z | / | dn ( u | m ) |
Y | m | cn ( u | m ) |
X | u | sn ( u | m ) |
If 0 < = m < 1 , m is saved in T and u is saved in L
Examples:
1- Evaluate sn ( 0.7 | 0.3 ) cn ( 0.7 | 0.3 ) dn ( 0.7 | 0.3 )
0.3 ENTER^
0.7 XEQ "JEF"
gives sn ( 0.7 | 0.3 )
= 0.632304777
--Execution time = 11s---
RDN cn ( 0.7 | 0.3 ) = 0.774719736
RDN dn ( 0.7 | 0.3 ) = 0.938113640
2-Likewise,
sn ( 0.7 | 1 ) = 0.604367777
sn ( 0.7 | 2 ) = 0.564297008
sn ( 0.7 | -3 ) = 0.759113422
cn ( 0.7 | 1 ) = 0.796705460
cn ( 0.7 | 2 ) = 0.825571855
cn ( 0.7 | -3 ) = 0.650958381
dn ( 0.7 | 1 ) = 0.796705460
dn ( 0.7 | 2 ) = 0.602609139
dn ( 0.7 | -3 ) =1.651895747
Notes:
-It seems that n is never greater than 7.
-Therefore, synthetic registers M and N can be replaced by R08 and
R09
-If m < -9999999999 the program can give wrong results!
2°) Using Gauss' Transformations
a) Program#1
-If 0 <= m < 1 , this variant employs the following formulae:
-With m' = 1-m , let be µ = [ ( 1-sqrt(m') / ( 1+sqrt(m') ]2 and v = u / ( 1+sqrt(µ) ] , we have:
sn ( u | m ) = [ ( 1 + sqrt(µ) ) sn ( v | µ
) ] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]
cn ( u | m ) = [ cn ( v | µ ) dn ( v | µ )
] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]
dn ( u | m ) = [ 1 - sqrt(µ) sn2 ( v
| µ ) ] / [ 1 + sqrt(µ) sn2 ( v | µ ) ]
-These formulas are applied recursively until µ is small enough to use
sn ( v | 0 ) = Sin v
cn ( v | 0 ) = Cos v
dn ( v | 0 ) = 1
Data Registers: R00 = 1 , 2 ,
... , n and then n-1, ... , 1 , 0 ;
R01 = sqrt(µ1) , ..... , Rnn = sqrt(µn)
Flags: F06-F07
Subroutines: /
-Line 41 should be a three-byte GTO 04. It may also be replaced by
RTN thus saving 2 bytes.
01 LBL "JEF"
02 DEG 03 CF 06 04 CF 07 05 X<>Y 06 X#0? 07 X>0? 08 GTO 00 09 CHS 10 STO Z 11 1 12 + 13 ST/ Z 14 SQRT 15 * 16 X<>Y 17 SF 06 18 LBL 00 19 1 20 X>Y? 21 GTO 01 22 X#Y? 23 GTO 00 24 X<> Z 25 ENTER^ 26 E^X 27 LASTX 28 CHS |
29 E^X
30 + 31 2 32 X<>Y 33 / 34 X<>Y 35 ST+ X 36 E^X-1 37 RCL X 38 2 39 + 40 / 41 GTO 04 42 LBL 00 43 RDN 44 STO Z 45 SQRT 46 * 47 X<>Y 48 1/X 49 ENTER^ 50 SF 07 51 LBL 01 52 CLX 53 STO 00 54 RDN 55 STO M 56 X<>Y |
57 LBL 02
58 1 59 RCL Z 60 - 61 SQRT 62 1 63 X<>Y 64 - 65 X=0? 66 GTO 00 67 1 68 ST+ 00 69 LASTX 70 + 71 / 72 STO IND 00 73 X^2 74 X<>Y 75 1 76 LASTX 77 + 78 / 79 GTO 02 80 LBL 00 81 SIGN 82 X<>Y 83 R-D 84 1 |
85 P-R
86 X<>Y 87 LBL 03 88 STO N 89 RDN 90 * 91 R^ 92 X^2 93 RCL IND 00 94 * 95 ENTER^ 96 ENTER^ 97 SIGN 98 ST- Z 99 + 100 ST/ N 101 ST/ Z 102 / 103 CHS 104 X<>Y 105 RCL IND 00 106 1 107 + 108 RCL N 109 * 110 DSE 00 111 GTO 03 112 FS?C 06 |
113 GTO 00
114 FC?C 07 115 GTO 04 116 RDN 117 X<>Y 118 R^ 119 RCL M 120 SQRT 121 * 122 GTO 04 123 LBL 00 124 X<> Z 125 ST/ Y 126 ST/ Z 127 1/X 128 X<> Z 129 1 130 ST- M 131 X<> M 132 CHS 133 SQRT 134 * 135 LBL 04 136 CLA 137 END |
( 187 bytes / SIZE 008 )
STACK | INPUTS | OUTPUTS |
Z | / | dn ( u | m ) |
Y | m | cn ( u | m ) |
X | u | sn ( u | m ) |
Examples:
1- Evaluate sn ( 0.7 | 0.3 ) cn ( 0.7 | 0.3 ) dn ( 0.7 | 0.3 )
0.3 ENTER^
0.7 XEQ "JEF"
gives sn ( 0.7 | 0.3 )
= 0.632304776
--Execution time = 7s---
RDN cn ( 0.7 | 0.3 ) = 0.774719736
RDN dn ( 0.7 | 0.3 ) = 0.938113640
2-Likewise,
sn ( 0.7 | 1 ) = 0.604367777
sn ( 0.7 | 2 ) = 0.564297007
sn ( 0.7 | -3 ) = 0.759113421
cn ( 0.7 | 1 ) = 0.796705460
cn ( 0.7 | 2 ) = 0.825571855
cn ( 0.7 | -3 ) = 0.650958382
dn ( 0.7 | 1 ) = 0.796705460
dn ( 0.7 | 2 ) = 0.602609138
dn ( 0.7 | -3 ) =1.651895746
Notes:
-It also seems that n is never greater than 7.
-Therefore, synthetic registers M and N can be replaced by R08 and
R09
-Storing 1 in synthetic register O and replacing 1 by RCL O in
the loops would reduce execution time.
-If m < -9999999999 the program can give wrong results!
-If m = 0.9999999999 , execution time = 14 seconds.
b) M-Code Routine (
0 <= m < 1 )
-Lines 57 to 111 are the slow part of the program above.
-The following routine "AJF" does the same job much faster.
086 "F"
00A "J"
001 "A"
378 C=c
03C RCR 3
106 A=C S&X
1BC RCR 11
130 LDI S&X
1F8 1F8h=504d
Correct value for HP-41 CX , HP-41CV or HP-41C with a quadmemory
module
306 ?A<C S&X
381 ?NCGO
00A NEXIST
The existence of the required registers is checked
0A6 A<>C S&X
228 P=C
0B8 C=Y
2FE ?C#0 MS
0B5 ?C GO
0A3 DATAERROR
DATA ERROR is displayed if m < 0
028 T=C
10E A=C ALL
04E C
35C =
050 1
36E ?A#C ALL
0B5 ?NC GO
0A2 DATAERROR
DATA ERROR is displayed if m = 1
068 Z=C
0F8 C=X
128 L=C
2A0 SETDEC
LOOP1
0B8 C=Y
2BE C=-C
10E A=C ALL
078 C=Z
01D C=
060 A+C
305 C=
If m > 1 , the HP-41 will display DATA ERROR
060 sqrt(AB)
when it tries to calculate sqrt(1-m)
1E8 O=C
2BE C=-C
10E A=C ALL
078 C=Z
01D C=
060 A+C
2EE ?C#0 ALL
133 JNC+38d
GOTO LBL 1
070 N=C ALL
1F8 C=O
10E A=C ALL
078 C=Z
01D C=
060 A+C
10E A=C ALL
0B0 C=N ALL
0AE A<>C ALL
261 C=
060 A/C
1E8 O=C
001 C=
060 AB+1
10E A=C ALL
0F8 C=X
0AE A<>C ALL
261 C=
060 A/C
0E8 X=C
1F8 C=O
10E A=C ALL
135 C=
060 A*C
0A8 Y=C
1F8 C=O
10E A=C ALL
238 C=P
260 SETHEX
226 C=C+1 S&X
228 P=C
270 RAMSLCT
0AE A<>C ALL
2F0 WRITDATA
046 C=0 S&X
270 RAMSLCT
25B JNC-53d
GOTO LOOP1
0F8 C=X
LBL 1
231 performs
064 R-D
0A8 Y=C
10E A=C
078 C=Z
0E8 X=C
070 C&N
351
084 =
1D5
078 P-R(C,A)
0A8 Y=C
0B0 C=N
0E8 X=C
238 C=P
106 A=C S&X
LOOP2
03C RCR 3
366 ?A#C S&X
3A0 ?NC RTN
The routine stops here
078 C=Z
10E A=C ALL
0B8 C=Y
135 C=
060 A*C
0A8 Y=C
0F8 C=X
10E A=C ALL
135 C=
060 A*C
238 C=P
270 RAMSLCT
038 READATA
070 N=C ALL
046 C=0 S&X
270 RAMSLCT
0B0 C=N ALL
13D C=
060 AB*C
2BE C=-C
1E8 O=C
001 C=
060 AB+1
1A8 N=C
04E C
35C =
050 1
268 Q=C
10E A=C ALL
1F8 C=O
01D C=
060 A+C
1B8 C=N
269 C=
060 AB/C
068 Z=C
0B8 C=Y
10E A=C ALL
1B8 C=N
261 C=
060 A/C
0A8 Y=C
278 C=Q
10E A=C ALL
0B0 C=N ALL
01D C=
060 A+C
0F8 C=X
13D C=
060 AB*C
1B8 C=N
269 C=
060 AB/C
0E8 X=C
260 SETHEX
238 C=P
266 C=C-1 S&X
228 P=C
2A0 SETDEC
20B JNC-63d
GOTO LOOP2
( 163 words / SIZE 008 )
STACK | INPUTS | OUTPUTS |
T | / | m |
Z | / | dn ( u | m ) |
Y | m | cn ( u | m ) |
X | u | sn ( u | m ) |
L | / | u |
with 0 <= m < 1
Example:
-Evaluate sn ( 0.7 | 0.3 ) cn ( 0.7 | 0.3 ) dn ( 0.7 | 0.3 )
0.3 ENTER^
0.7 XEQ "AJF"
gives sn ( 0.7 | 0.3 )
= 0.632304777
--Execution time = 3.5s---
RDN cn ( 0.7 | 0.3 ) = 0.774719736
RDN dn ( 0.7 | 0.3 ) = 0.938113640
Notes:
-If m = 0.9999999999 , execution time = 5.5 seconds instead
of 14 seconds.
-The alpha register is not cleared and synthetic register Q
is used.
-There is no check for alpha data, underflow or underflow.
c) Program#2
-The program listed in paragraph 2-a) may now be simplified:
Data Registers: R00 is unused
, R01 = sqrt(µ1) , ..... ,
Rnn = sqrt(µn)
Flags: F06-F07
Subroutines: /
01 LBL "JEF"
02 DEG 03 CF 06 04 CF 07 05 X<>Y 06 X#0? 07 X>0? 08 GTO 00 09 CHS 10 STO Z 11 1 12 + 13 ST/ Z 14 SQRT 15 * 16 X<>Y 17 SF 06 |
18 LBL 00
19 1 20 X>Y? 21 GTO 01 22 X#Y? 23 GTO 00 24 X<> Z 25 ENTER^ 26 E^X 27 LASTX 28 CHS 29 E^X 30 + 31 2 32 X<>Y 33 / 34 X<>Y |
35 ST+ X
36 E^X-1 37 RCL X 38 2 39 + 40 / 41 GTO 02 42 LBL 00 43 RDN 44 STO Z 45 SQRT 46 * 47 X<>Y 48 1/X 49 ENTER^ 50 SF 07 51 LBL 01 |
52 RDN
53 STO M 54 X<>Y 55 AJF 56 FS?C 06 57 GTO 00 58 FC?C 07 59 GTO 04 60 RDN 61 X<>Y 62 R^ 63 RCL M 64 SQRT 65 * 66 GTO 02 67 LBL 00 68 X<> Z |
69 ST/ Y
70 ST/ Z 71 1/X 72 X<> Z 73 1 74 ST- M 75 X<> M 76 CHS 77 SQRT 78 * 79 LBL 02 80 CLA 81 END |
( 117 bytes / SIZE 008 )
STACK | INPUTS | OUTPUTS |
Z | / | dn ( u | m ) |
Y | m | cn ( u | m ) |
X | u | sn ( u | m ) |
Examples:
1- Evaluate sn ( 0.7 | 0.3 ) cn ( 0.7 | 0.3 ) dn ( 0.7 | 0.3 )
0.3 ENTER^
0.7 XEQ "JEF"
gives sn ( 0.7 | 0.3 )
= 0.632304777
--Execution time = 4s---
RDN cn ( 0.7 | 0.3 ) = 0.774719736
RDN dn ( 0.7 | 0.3 ) = 0.938113640
2-Likewise,
sn ( 0.7 | 1 ) = 0.604367777
sn ( 0.7 | 2 ) = 0.564297007
sn ( 0.7 | -3 ) = 0.759113420
cn ( 0.7 | 1 ) = 0.796705460
cn ( 0.7 | 2 ) = 0.825571855
cn ( 0.7 | -3 ) = 0.650958382
dn ( 0.7 | 1 ) = 0.796705460
dn ( 0.7 | 2 ) = 0.602609138
dn ( 0.7 | -3 ) =1.651895745
Notes:
-Since R00 is unused by the M-code routine AJF, you can replace synthetic register M by R00 in the listing above.
-If m < -9999999999 the program can give wrong results !
-If m = 0.9999999999 , execution time = 6 seconds instead of
14 seconds.
-If 0 <= m < 1 , m is saved in T and u is saved in L-register.
3°) Jacobian Elliptic Function Sn (
x | m ) 0 <= m < 1
a) Focal Program
-If you only have to compute sn ( x | m ) for 0 <= m
< 1 - for instance to solve the Ellipsoidal Wave differential equation
-
the following program is faster than "JEF"
-Gauss' transformation is used
Data Registers: R00 = 1 , 2 ,
... , n and then n-1, ... , 1 , 0 ;
R01 = sqrt(µ1) , ..... , Rnn = sqrt(µn)
Flags: /
Subroutines: /
01 LBL "SNX"
02 DEG 03 0 04 STO 00 05 RDN 06 LBL 01 07 1 08 RCL Z 09 - 10 SQRT |
11 1
12 X<>Y 13 - 14 X=0? 15 GTO 00 16 1 17 ST+ 00 18 LASTX 19 + 20 / |
21 STO IND 00
22 X^2 23 X<>Y 24 1 25 LASTX 26 + 27 / 28 GTO 01 29 LBL 00 30 X<>Y |
31 R-D
32 SIN 33 LBL 02 34 ENTER^ 35 X^2 36 RCL IND 00 37 * 38 1 39 + 40 X<>Y |
41 RCL IND 00
42 LASTX 43 + 44 * 45 X<>Y 46 / 47 DSE 00 48 GTO 02 49 END |
( 66 bytes / SIZE 008 )
STACK | INPUTS | OUTPUTS |
Y | m | / |
X | u | sn ( u | m ) |
Example:
Evaluate sn ( 0.7 | 0.8 )
0.8 ENTER^
0.7 XEQ "SNX" >>>>
sn ( 0.7 | 0.8 ) = 0.612365484
---Execution time = 7s---
Notes:
-There will be an infinte loop if m = 1
-If m = 0.9999999999 , execution time = 11 seconds instead of
14 seconds.
b) M-Code Routine
-Likewise, the M-Code routine may be simplified to return sn( u | m
) only:
098 "X"
00E "N"
013 "S"
378 C=c
03C RCR 3
106 A=C S&X
1BC RCR 11
130 LDI S&X
1F8 1F8h=504d
Correct value for HP-41 CX , HP-41CV or HP-41C with a quadmemory
module
306 ?A<C S&X
381 ?NCGO
00A NEXIST
The existence of the required registers is checked
0A6 A<>C S&X
228 P=C
0B8 C=Y
2FE ?C#0 MS
0B5 ?C GO
0A3 DATAERROR
DATA ERROR is displayed if m < 0
168 M=C
10E A=C ALL
04E C
35C =
050 1
36E ?A#C ALL
0B5 ?NC GO
0A2 DATAERROR
DATA ERROR is displayed if m = 1
1A8 N=C
0F8 C=X
128 L=C
2A0 SETDEC
LOOP1
0B8 C=Y
2BE C=-C
10E A=C ALL
1B8 C=N
01D C=
060 A+C
305 C=
If m > 1 , the HP-41 will display DATA ERROR
060 sqrt(AB)
when it tries to calculate sqrt(1-m)
1E8 O=C
2BE C=-C
10E A=C ALL
1B8 C=N
01D C=
060 A+C
2EE ?C#0 ALL
133 JNC+38d
GOTO LBL 1
070 N=C ALL
1F8 C=O
10E A=C ALL
1B8 C=N
01D C=
060 A+C
10E A=C ALL
0B0 C=N ALL
0AE A<>C ALL
261 C=
060 A/C
1E8 O=C
001 C=
060 AB+1
10E A=C ALL
0F8 C=X
0AE A<>C ALL
261 C=
060 A/C
0E8 X=C
1F8 C=O
10E A=C ALL
135 C=
060 A*C
0A8 Y=C
1F8 C=O
10E A=C ALL
238 C=P
260 SETHEX
226 C=C+1 S&X
228 P=C
270 RAMSLCT
0AE A<>C ALL
2F0 WRITDATA
046 C=0 S&X
270 RAMSLCT
25B JNC-53d
GOTO LOOP1
130 LDI S&X
LBL 1
090 090h
these instructions
358 ST=C XP
set the required
144 CLRF 6
CPU flags
284 CLRF 7
so that Sin(C)
0F8 C=X
is executed in radians
070 N=C ALL
without changing the current angular mode
229 C=
048 Sin(C)
0E8 X=C
238 C=P
106 A=C S&X
LOOP2
03C RCR 3
366 ?A#C S&X
11B JNC+35d
GTO LBL2
0F8 C=X
10E A=C ALL
135 C=
060 A*C
238 C=P
270 RAMSLCT
038 READATA
070 N=C ALL
046 C=0 S&X
270 RAMSLCT
0B0 C=N ALL
13D C=
060 AB*C
001 C=
060 AB+1
0A8 Y=C
1B8 C=N
10E A=C ALL
0B0 C=N ALL
01D C=
060 A+C
0F8 C=X
13D C=
060 AB*C
0B8 C=N
269 C=
060 AB/C
0E8 X=C
260 SETHEX
238 C=P
266 C=C-1 S&X
228 P=C
2A0 SETDEC
2DB JNC-37d
GTO LOOP2
178 C=M
LBL2
0A8 Y=C
345 ?NCGO
042 CLA
( 136 words / SIZE 008 )
STACK | INPUTS | OUTPUTS |
T | T | T |
Z | Z | Z |
Y | m | m |
X | u | sn ( u | m ) |
L | L | u |
with 0 <= m < 1
Example:
-Evaluate sn ( 0.7 | 0.3 )
0.3 ENTER^
0.7 XEQ "SNX"
gives sn ( 0.7 | 0.3 )
= 0.632304777
--Execution time = 2.2s---
Notes:
-If m = 0.8 , execution time = 2.5 seconds.
-If m = 0.9999999999 , execution time = 4 seconds.
-Thus, SNX is significantly faster than AJF.
-Moreover, registers Z & T are saved.
-The alpha register is cleared and synthetic register Q is used.
-There is no check for alpha data, underflow or underflow.
References:
[1] Abramowitz and Stegun - "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4
[2] http://functions.wolfram.com/