# hp41programs

Quternionic Functions

# Quaternionic Special Functions for the HP-41

Overview

1°)  Gamma Function

a)  Asymptotic Expansion
b)  Continued Fraction
c)  Reciprocal of the Gamma Function

2°)  Digamma Function
3°)  Error Function
4°)  Exponential Integrals

a)  Ei(q)
b)  En(q)  n # 1 , 2 , 3 , .........

5°)  Sine & Hyperbolic Sine Integrals
6°)  Cosine & Hyperbolic Cosine Integrals
7°)  Bessel Functions - Real order

a)  Bessel Functions of the 1st kind
b)  Bessel Functions of the 2nd kind - non-integer order
c)  Bessel Functions of the 2nd kind - integer order

8°)  Kummer's Function
9°)  Parabolic Cylinder Functions

a)  Via Kummer's Function
b)  Asymptotic Expansion

10°)  Struve Functions
11°)  Hypergeometric Function
12°)  Associated Legendre Functions of the 1st kind

a)  Program#1 - Functions of Type 2 / Type 3
b)  Integer Order and Degree - Functions of Type 2 / Type 3

13°)  Associated Legendre Functions of the 2nd kind

a)  Program#1 - Function of Type 2 - | q | < 1
b)  Program#2 - Function of Type 3 - | q | > 1

14°)  Anger & Weber Functions
15°)  Regular Coulomb Wave Function

-When ascending series are employed to compute these functions, accurate results are obtained for "small" arguments.
-For "large" quaternions, the precision is much more uncertain!

-The variables are quaternions but the indexes and orders are restricted to real values.
-When the function involves one index, it's to be stored in register R00
-When several indexes and/or orders are involved, they are to be stored in R09 R10 , ....
-The content of these registers may be altered during the calculations, but they are restored at the end.

1°)  Gamma Function

a)  Asymptotic Expansion

-The following asymptotic series is used:

Gam(q) ~  exp { (q-1/2) Ln q + Ln (2.PI)1/2  - q + ( 1/12.q ) [ 1 - ( 1/30.q2 ) + 1/( 105.q4 ) - 1/( 140.q6 ) + 1/( 99.q8 ) ] }    for  Re(q) > 5

-The relation   Gam(q+1) = q Gam(q)  is used recursively if Re(q) < 5 until   Re(q+1+.......+1)  > 5

Data Registers:   R00 unused     R01 to R16: temp
Flags: /
Subroutines:   "Q*Q"  "1/Q"  "LNQ"  "E^Q"  ( cf  "Quaternions for the HP-41" )

 01  LBL "GAMQ"   02  STO 01            03  RDN   04  STO 02   05  STO 10   06  RDN   07  STO 03   08  STO 11   09  X<>Y   10  STO 04   11  STO 12   12  CLX   13  STO 06   14  STO 07   15  STO 08   16  SIGN   17  STO 05   18  LBL 01   19  XEQ "Q*Q"   20  STO 05   21  RDN   22  STO 06   23  RDN   24  STO 07   25  X<>Y   26  STO 08   27  1   28  ST+ 01   29  RCL 01   30  5   31  X>Y?   32  GTO 01   33  RCL 05   34  STO 13   35  RCL 06   36  STO 14   37  RCL 07 38  STO 15    39  RCL 08            40  STO 16   41  RCL 01   42  STO 05   43  STO 09   44  RCL 02   45  STO 06   46  RCL 03   47  STO 07   48  RCL 04   49  STO 08   50  XEQ "Q*Q"   51  XEQ "1/Q"   52  STO 01   53  STO 05   54  RDN   55  STO 02   56  STO 06   57  RDN   58  STO 03   59  STO 07   60  X<>Y   61  STO 04   62  99   63  ST/ 05   64  ST/ 06   65  ST/ 07   66  /   67  STO 08   68  140   69  1/X   70  ST- 05   71  XEQ "Q*Q"   72  STO 05   73  RDN   74  STO 06 75  RDN   76  STO 07            77  X<>Y   78  STO 08   79  105   80  1/X   81  ST+ 05   82  XEQ "Q*Q"    83  STO 05   84  RDN   85  STO 06   86  RDN   87  STO 07   88  X<>Y   89  STO 08   90  30   91  1/X   92  ST- 05   93  XEQ "Q*Q"   94  STO 05   95  RDN   96  STO 06   97  RDN   98  STO 07   99  X<>Y 100  STO 08 101  1 102  ST+ 05 103  RCL 12 104  RCL 11 105  RCL 10 106  RCL 09 107  XEQ "1/Q" 108  STO 01 109  CLX 110  12 111  ST/ 01 112  ST/ T 113  ST/ Z 114  / 115  STO 02          116  RDN 117  STO 03 118  X<>Y 119  STO 04 120  XEQ "Q*Q"  121  X<> 09 122  STO 01 123  ST- 09 124  RDN 125  X<> 10 126  STO 02 127  ST- 10 128  RDN 129  X<> 11 130  STO 03 131  ST- 11 132  RDN 133  X<> 12 134  STO 04 135  ST- 12 136  RDN 137  XEQ "LNQ" 138  STO 05 139  RDN 140  STO 06 141  RDN 142  STO 07 143  X<>Y 144  STO 08 145  PI 146  .5 147  ST- 01 148  / 149  SQRT 150  LN 151  ST+ 09 152  XEQ "Q*Q"  153  ST+ 09 154  RDN 155  ST+ 10 156  RDN 157  ST+ 11 158  CLX 159  RCL12          160  + 161  RCL 11 162  RCL 10 163  RCL 09 164  XEQ "E^Q" 165  STO 05 166  RDN 167  STO 06 168  RDN 169  STO 07 170  X<>Y 171  STO 08 172  RCL 16 173  RCL 15 174  RCL 14 175  RCL 13 176  XEQ "1/Q" 177  STO 01 178  RDN 179  STO 02 180  RDN 181  STO 03 182  X<>Y 183  STO 04 184  XEQ "Q*Q" 185  END

( 281 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Gam (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 4 + 3 i + 2 j + k

1   ENTER^
2   ENTER^
3   ENTER^
4   XEQ "GAMQ"   >>>>     0.541968821                     ---Execution time = 29s---
RDN    -0.740334194
RDN    -0.493556130
RDN    -0.246778065

Whence,   Gam ( 4 + 3 i + 2 j + k ) =  0.541968821 - 0.740334194 i - 0.493556130 j - 0.246778065 k

b)  Continued Fraction

-Another method is   Gam(q) =  exp [ (q-1/2) Ln q + Ln (2.PI)1/2  - q + ( 1/12 )/( q + ( 1/30 )/( q + ( 53/210)/( q + (195/371)/( q + ...  )))) ]

Data Registers:   R00 unused     R01 to R16: temp
Flags: /
Subroutines:   "Q*Q"  "1/Q"  "LNQ"  "E^Q"  ( cf  "Quaternions for the HP-41" )

 01  LBL "GAMQ"   02  STO 01   03  RDN   04  STO 02   05  RDN   06  STO 03   07  X<>Y   08  STO 04   09  CLX   10  STO 06   11  STO 07   12  STO 08   13  SIGN   14  STO 05   15  LBL 01   16  XEQ "Q*Q"   17  STO 05   18  RDN   19  STO 06   20  RDN   21  STO 07   22  X<>Y   23  STO 08   24  1   25  ST+ 01   26  RCL 01   27  5   28  X>Y?   29  GTO 01   30  RCL 04   31  RCL 03   32  RCL 02   33  RCL 01   34  XEQ "1/Q"   35  STO 09 36  RDN   37  STO 10            38  CLX   39  195   40  371   41  /   42  ST* 09   43  ST* 10   44  ST* Z   45  *   46  RCL 04   47  ST+ Z   48  CLX   49  RCL 03   50  +   51  RCL 02   52  RCL 10   53  +   54  RCL 01   55  ST+ 09   56  X<> 09   57  XEQ "1/Q"   58  STO 09   59  RDN   60  STO 10   61  CLX   62  53   63  210   64  /   65  ST* 09   66  ST* 10   67  ST* Z   68  *   69  RCL 04   70  ST+ Z 71  CLX   72  RCL 03    73  +   74  RCL 02            75  RCL 10   76  +   77  RCL 01   78  ST+ 09   79  CLX   80  30   81  ST* 09   82  ST* T   83  ST* Z   84  *   85  RCL 09   86  XEQ "1/Q"   87  STO 09   88  CLX   89  RCL 04   90  ST+ T   91  CLX   92  RCL 03   93  ST+ Z   94  CLX   95  RCL 02   96  +   97  RCL 01   98  ST+ 09   99  CLX 100  12 101  ST* 09 102  ST* T 103  ST* Z 104  * 105  RCL 09 106  XEQ "1/Q" 107  STO 09  108  RDN 109  STO 10          110  RDN 111  STO 11 112  X<>Y 113  STO 12 114  RCL 04 115  ST- 12 116  RCL 03 117  ST- 11 118  RCL 02 119  ST- 10 120  RCL 01 121  ST- 09 122  XEQ "LNQ" 123  X<> 05 124  STO 13 125  RDN 126  X<> 06 127  STO 14 128  RDN 129  X<> 07 130  STO 15 131  X<>Y 132  X<> 08 133  STO 16 134  PI 135  .5 136  ST- 01 137  / 138  SQRT 139  LN 140  ST+ 09 141  XEQ "Q*Q" 142  ST+ 09 143  CLX 144  RCL 12          145  ST+ T 146  CLX 147  RCL 11 148  ST+ Z 149  CLX 150  RCL 10 151  + 152  RCL 09 153  XEQ "E^Q" 154  STO 05 155  RDN 156  STO 06 157  RDN 158  STO 07 159  X<>Y 160  STO 08 161  RCL 16 162  RCL 15 163  RCL 14 164  RCL 13 165  XEQ "1/Q" 166  STO 01 167  RDN 168  STO 02 169  RDN 170  STO 03 171  X<>Y 172  STO 04 173  XEQ "Q*Q" 174  END

( 272 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Gam (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 4 + 3 i + 2 j + k

1   ENTER^
2   ENTER^
3   ENTER^
4   XEQ "GAMQ"   >>>>     0.541968821                     ---Execution time = 22s---
RDN    -0.740334194
RDN    -0.493556130
RDN    -0.246778065

Gam ( 4 + 3 i + 2 j + k ) =  0.541968821 - 0.740334194 i - 0.493556130 j - 0.246778065 k

Note:

-The second version is slightly faster than the first one
-The accuracy is the same.

c)  Reciprocal of the Gamma Function

-The same continued fraction is used to compute 1/Gam(q)
-Lines 02 to 100 are identical.

Data Registers:   R00 unused     R01 to R16: temp
Flags: /
Subroutines:   "Q*Q"  "1/Q"  "LNQ"  "E^Q"  ( cf  "Quaternions for the HP-41" )

 01  LBL "1/GMQ"   02  STO 01             03  RDN   04  STO 02   05  RDN   06  STO 03   07  X<>Y   08  STO 04   09  CLX   10  STO 06   11  STO 07   12  STO 08   13  SIGN   14  STO 05   15  LBL 01   16  XEQ "Q*Q"   17  STO 05   18  RDN   19  STO 06   20  RDN   21  STO 07   22  X<>Y   23  STO 08   24  1   25  ST+ 01   26  RCL 01   27  5   28  X>Y?   29  GTO 01   30  RCL 04   31  RCL 03   32  RCL 02   33  RCL 01   34  XEQ "1/Q"   35  STO 09 36  RDN   37  STO 10             38  CLX   39  195   40  371   41  /   42  ST* 09   43  ST* 10   44  ST* Z   45  *   46  RCL 04   47  ST+ Z   48  CLX   49  RCL 03   50  +   51  RCL 02   52  RCL 10   53  +   54  RCL 01   55  ST+ 09   56  X<> 09   57  XEQ "1/Q"   58  STO 09   59  RDN   60  STO 10   61  CLX   62  53   63  210   64  /   65  ST* 09   66  ST* 10   67  ST* Z   68  *   69  RCL 04   70  ST+ Z 71  CLX   72  RCL 03             73  +   74  RCL 02   75  RCL 10   76  +   77  RCL 01   78  ST+ 09   79  CLX   80  30   81  ST* 09   82  ST* T   83  ST* Z   84  *   85  RCL 09   86  XEQ "1/Q"   87  STO 09   88  CLX   89  RCL 04   90  ST+ T   91  CLX   92  RCL 03   93  ST+ Z   94  CLX   95  RCL 02   96  +   97  RCL 01   98  ST+ 09   99  CLX 100  12 101  CHS 102  ST* 09 103  ST* T 104  ST* Z 105  * 106  RCL 09 107  XEQ "1/Q" 108  STO 09           109  RDN 110  STO 10 111  RDN 112  STO 11 113  X<>Y 114  STO 12 115  RCL 04 116  ST+ 12 117  RCL 03 118  ST+ 11 119  RCL 02 120  ST+ 10 121  RCL 01 122  ST+ 09 123  XEQ "LNQ" 124  X<> 05 125  STO 13 126  RDN 127  X<> 06 128  STO 14 129  RDN 130  X<> 07 131  STO 15 132  X<>Y 133  X<> 08 134  STO 16 135  PI 136  .5 137  ST- 01 138  / 139  SQRT 140  LN 141  ST- 09 142  XEQ "Q*Q"  143  ST- 09 144  RDN 145  ST- 10 146  RDN 147  ST- 11 148  X<>Y 149  ST- 12 150  RCL 12           151  RCL 11 152  RCL 10 153  RCL 09 154  XEQ "E^Q" 155  STO 05 156  RDN 157  STO 06 158  RDN 159  STO 07 160  X<>Y 161  STO 08 162  RCL 13 163  STO 01 164  RCL 14 165  STO 02 166  RCL 15 167  STO 03 168  RCL 16 169  STO 04 170  XEQ "Q*Q" 171  END

( 270 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      1/Gam (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 4 + 3 i + 2 j + k

1   ENTER^
2   ENTER^
3   ENTER^
4   XEQ "1/GMQ"   >>>>     0.472789344                     ---Execution time = 21s---
RDN     0.645834418
RDN     0.430556279
RDN     0.215278139

1/Gam ( 4 + 3 i + 2 j + k ) =  0.472789344 + 0.645834418 i + 0.430556279 j + 0.215278139 k

2°)  Digamma Function

Formula:             Psi(q) ~ Ln q - 1/(2q) -1/(12q2) + 1/(120q4) - 1/(252q6) + 1/(240q8)    is used if  Re(q) > 8

Psi(q+1) = Psi(q) + 1/q    is used recursively until  Re(q+1+....+1) > 8

Data Registers:   R00 unused     R01 to R16: temp   When the program stops,  Psi(q) is in registers R01 thru R04
Flags: /
Subroutines:   "Q*Q"  "1/Q"  "LNQ"  ( cf  "Quaternions for the HP-41" )

 01  LBL "PSIQ"   02  STO 01          03  RDN   04  STO 02   05  STO 06   06  STO 10   07  RDN   08  STO 03   09  STO 07   10  STO 11   11  X<>Y   12  STO 04   13  STO 08   14  STO 12   15  CLX   16  STO 13   17  STO 14   18  STO 15   19  STO 16   20  LBL 01   21  RCL 04   22  RCL 03   23  RCL 02   24  RCL 01   25  XEQ "1/Q"   26  ST- 13   27  RDN   28  ST- 14 29  RDN   30  ST- 15   31  X<>Y   32  ST- 16   33  1   34  ST+ 01   35  RCL 01   36  8   37  X>Y?   38  GTO 01   39  RCL 01   40  STO 05   41  STO 09          42  XEQ "Q*Q"   43  XEQ "1/Q"   44  STO 01   45  STO 05   46  RDN   47  STO 02   48  STO 06   49  RDN   50  STO 03   51  STO 07   52  X<>Y   53  STO 04   54  STO 08   55  20   56  ST/ 05 57  ST/ 06   58  ST/ 07   59  ST/ 08   60  21   61  1/X   62  ST- 05   63  XEQ "Q*Q"   64  STO 05   65  RDN   66  STO 06          67  RDN   68  STO 07   69  X<>Y   70  STO 08   71  .1   72  ST+ 05   73  XEQ "Q*Q"   74  STO 05   75  RDN   76  STO 06   77  RDN   78  STO 07   79  X<>Y   80  STO 08   81  1   82  ST- 05   83  XEQ "Q*Q"   84  STO 01 85  RDN   86  STO 02   87  RDN   88  STO 03   89  X<>Y   90  STO 04   91  6   92  ST/ 01   93  ST/ 02   94  ST/ 03   95  ST/ 04   96  RCL 12          97  RCL 11   98  RCL 10   99  RCL 09 100  XEQ "1/Q"  101  ST- 01 102  RDN 103  ST- 02 104  RDN 105  ST- 03 106  X<>Y 107  ST- 04 108  2 109  ST/ 01 110  ST/ 02 111  ST/ 03 112  ST/ 04 113  RCL 12 114  RCL 11 115  RCL 10 116  RCL 09 117  XEQ "LNQ" 118  ST+ 01 119  CLX 120  RCL 14 121  + 122  ST+ 02 123  CLX 124  RCL 15        125  + 126  ST+ 03 127  CLX 128  RCL 16 129  + 130  ST+ 04 131  RCL 13 132  ST+ 01 133  RCL 04 134  RCL 03 135  RCL 02 136  RCL 01 137  END

( 213 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Psi (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 2 i + 3 j + 4 k

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "PSIQ"   >>>>     1.686531556                     ---Execution time = 26s---
RDN     0.548896352
RDN     0.823344528
RDN     1.097792703

Psi ( 1 + 2 i + 3 j + 4 k ) =  1.686531556 + 0.548896352 i + 0.823344528 j + 1.097792703 k

3°)  Error Function

-This routine uses the formula:

erf q  = (2/pi1/2)  SUMn=0,1,2,.....  (-1)n q2n+1 / (n! (2n+1))

Data Registers:  R00 thru R12: temp

When the program stops, the result is in R13 to R16

Flags:  /
Subroutine:      "Q*Q"  ( cf "Quaternions for the HP-41" )

 01  LBL "ERFQ"   02  STO 01          03  STO 05   04  STO 13   05  RDN   06  STO 02   07  STO 06   08  STO 14   09  RDN   10  STO 03   11  STO 07   12  STO 15   13  X<>Y   14  STO 04   15  STO 08   16  STO 16   17  XEQ "Q*Q"   18  STO 09   19  RDN   20  STO 10   21  RDN 22  STO 11   23  X<>Y   24  STO 12          25  CLX   26  STO 00   27  LBL 01   28  RCL 09   29  STO 05   30  RCL 10   31  STO 06   32  RCL 11    33  STO 07   34  RCL 12   35  STO 08   36  RCL 00   37  1   38  +   39  STO 00   40  ST+ X   41  STO Y   42  LASTX 43  ST- Z   44  +   45  RCL 00          46  *   47  /   48  CHS   49  ST* 05   50  ST* 06   51  ST* 07   52  ST* 08   53  XEQ "Q*Q"   54  STO 01   55  RDN   56  STO 02   57  RCL 14   58  +   59  STO 14   60  LASTX   61  -   62  ABS   63  X<>Y 64  STO 03   65  RCL 15   66  +   67  STO 15          68  LASTX   69  -   70  ABS   71  +   72  X<>Y   73  STO 04   74  RCL 16    75  +   76  STO 16   77  LASTX   78  -   79  ABS   80  +   81  RCL 01   82  RCL 13   83  +   84  STO 13 85  LASTX   86  -   87  ABS   88  +   89  X#0?   90  GTO 01          91  2   92  PI   93  SQRT   94  /   95  ST* 13   96  ST* 14   97  ST* 15   98  ST* 16   99  RCL 16 100  RCL 15  101  RCL 14 102  RCL 13 103  END

( 136 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Erf (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 1.1 i + 1.2 j + 1.3 k

1.3   ENTER^
1.2   ENTER^
1.1   ENTER^
1     XEQ "ERFQ"   >>>>     -2.355991403                     ---Execution time = 104s---
RDN      -3.358261140
RDN      -3.663557610
RDN      -3.968854075

Whence,   Erf ( 1 + 1.1 i + 1.2 j + 1.3 k ) =  -2.355991403 - 3.358261140 i - 3.663557610 j - 3.968854075 k

4°)  Exponential Integrals

a) Ei(q)

-The following program employs

Ei(q)  = C + ln(q) + Sumn=1,2,.....  qn/(n.n!)                where  C = 0.5772156649... = Euler's constant.

Data Registers:  R00 thru R12: temp

When the program stops, the result is in R13 to R16

Flags:  /
Subroutines:     "LNQ"  "Q*Q"  ( cf "Quaternions for the HP-41" )

 01  LBL "EIQ" 02  STO 01 03  STO 09         04  STO 13 05  RDN 06  STO 02 07  STO 10 08  STO 14 09  RDN 10  STO 03 11  STO 11 12  STO 15 13  RDN 14  STO 04 15  STO 12 16  STO 16 17  RDN 18  XEQ "LNQ"  19  ST+ 13 20  RDN 21  ST+ 14 22  RDN 23 ST+ 15 24  X<>Y 25  ST+ 16 26  .5772156649 27  ST+ 13 28  1 29  STO 00         30  LBL 01 31  RCL 09 32  STO 05 33  RCL 10 34  STO 06 35  RCL 11 36  STO 07 37  RCL 12 38  STO 08 39  RCL 00 40  RCL 00 41  1 42  + 43  STO 00         44  X^2 45  / 46  ST* 05 47  ST* 06 48  ST* 07 49  ST* 08 50  XEQ "Q*Q"  51  STO 01 52  RDN 53  STO 02 54  RCL 14 55  + 56  STO 14 57  LASTX 58  - 59  ABS 60  X<>Y 61  STO 03 62  RCL 15         63  + 64  STO 15 65  LASTX 66  - 67  ABS 68  + 69  X<>Y 70  STO 04  71  RCL 16 72  + 73  STO 16 74  LASTX 75  - 76  ABS 77  + 78  RCL 01 79  RCL 13 80  + 81  STO 13         82  LASTX 83  - 84  ABS 85  + 86  X#0? 87  GTO 01 88  RCL 16  89  RCL 15 90  RCL 14 91  RCL 13 92  END

( 132 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Ei (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k

Example:              q = 1 + 2 i + 3 j + 4 k

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "EIQ"   >>>>     -0.381065061                     ---Execution time = 100s---
RDN      1.052055477
RDN      1.578083214
RDN      2.104110954

Whence,   Ei ( 1 + 2 i + 3 j + 4 k ) =  -0.381065061 + 1.052055477 i + 1.578083214 j + 2.104110954 k

a) En(q)  n # 1 , 2 , 3 , ............

Formula:     En(q) = qn-1 Gam(1-n) - [1/(1-n)] M( 1-n , 2-n ; -q )      where    M = Kummer's function

Data Registers:           •  R00 = n                                  ( Register R00 is to be initialized before executing "ENQ" )

R01 thru R15: temp
When the program stops, the result is in R11 to R14
Flags:  /
Subroutines:     "Q^R"                      ( cf "Quaternions for the HP-41" )
"1/G+" or "GAM+"  ( cf "Gamma Function for the HP-41" )
"KUMQ"                 ( cf paragraph 8 below )

 01  LBL "ENQ"  02  CHS 03  STO 01           04  RDN 05  CHS 06  STO 02 07  RDN 08  CHS 09  STO 03 10  X<>Y 11  CHS 12  STO 04 13  1 14  STO 10 15  RCL 00 16  STO 15 17  - 18  STO 09 19  ST+ 10 20  RCL 04 21  RCL 03 22  RCL 02 23  RCL 01 24  XEQ "KUMQ" 25  RCL 09           26  CHS 27  STO 00           28  ST/ 11 29  ST/ 12 30  ST/ 13 31  ST/ 14 32  RCL 09 33  XEQ "1/G+"  34  STO 09 35  RCL 04 36  CHS 37  RCL 03 38  CHS 39  RCL 02 40  CHS 41  RCL 01           42  CHS 43  XEQ "Q^R"  44  X<> 09 45  ST/ 09 46  ST/ T 47  ST/ Z 48  / 49  ST+ 12 50  RDN 51  ST+ 13 52  X<>Y 53  ST+ 14 54  RCL 09           55  ST+ 11 56  RCL 15  57  STO 00 58  RCL 14 59  RCL 13 60  RCL 12 61  RCL 11 62  END

( 97 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      En(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                 with       R00 = n # 1 , 2 , 3 , ...............

Example:         n = PI     q = 1 + i/2 + j/3 + k/4

PI   STO 00

4   1/X
3   1/X
2   1/X
1   XEQ "ENQ"   >>>>    0.066444330                         ---Execution time = 54s---
RDN   -0.059916131
RDN   -0.039944088
RDN   -0.029958064

EPI( 1 + i/2 + j/3 + k/4 ) = 0.066444330 - 0.059916131 i - 0.039944088 j - 0.029958064 k

5°)  Sine & Hyperbolic Sine Integrals

-The following series are used

Si(q)  = Sumn=0,1,2,..... (-1)n q2n+1/((2n+1).(2n+1)!)
Shi(q) = Sumn=0,1,2,.....  q2n+1/((2n+1).(2n+1)!)

Data Registers:  R00 thru R12: temp

When the program stops, the result is in R13 to R16

Flag:   F01           CF 01  to calculate Si(q)
SF 01  to calculate  Shi(q)

Subroutine:     "Q*Q"  ( cf "Quaternions for the HP-41" )

 01  LBL "SIQ" 02  STO 01 03  STO 05 04  STO 13 05  RDN 06  STO 02 07  STO 06 08  STO 14 09  RDN 10  STO 03 11  STO 07 12  STO 15 13  X<>Y 14  STO 04 15  STO 08 16  STO 16 17  XEQ "Q*Q" 18  STO 09       19  RDN 20  STO 10 21  RDN 22  STO 11 23  X<>Y 24  STO 12 25  CLX 26  STO 00 27  LBL 01 28  RCL 09 29  STO 05 30  RCL 10  31  STO 06 32  RCL 11 33  STO 07 34  RCL 12 35  STO 08 36  RCL 00       37  1 38  + 39  STO 00 40  ST+ X 41  ENTER^ 42  ENTER^ 43  ENTER^ 44  SIGN 45  ST- T 46  + 47  X^2 48  * 49  / 50  FC? 01 51  CHS 52  ST* 05 53  ST* 06 54  ST* 07 55  ST* 08 56  XEQ "Q*Q" 57  STO 01  58  RDN 59  STO 02 60  RCL 14 61  + 62  STO 14 63  LASTX 64  - 65  ABS 66  X<>Y 67  STO 03 68  RCL 15 69  + 70  STO 15 71  LASTX 72  - 73  ABS 74  + 75  X<>Y 76  STO 04      77  RCL 16       78  + 79  STO 16 80  LASTX 81  - 82  ABS 83  + 84  RCL 01 85  RCL 13 86  + 87  STO 13 88  LASTX 89  - 90  ABS 91  + 92  X#0? 93  GTO 01 94  RCL 16  95  RCL 15       96  RCL 14 97  RCL 13 98  END

( 126 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Si (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    CF 01
and      Shi (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    SF 01

Example:              q = 1 + 2 i + 3 j + 4 k

•   CF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "CIQ"   >>>>    17.98954627                     ---Execution time = 51s---
RDN      7.075096293
RDN    10.61264444
RDN    14.15019259

So,   Si ( 1 + 2 i + 3 j + 4 k ) =  17.98954627 + 7.075096293 i + 10.61264444 j + 14.15019259 k

•   SF 01

4   ENTER^
3   ENTER^
2   ENTER^
1      R/S             >>>>    -0.160779306                     ---Execution time = 55s---
RDN      0.522153864
RDN      0.783230794
RDN      1.044307726

And   Shi ( 1 + 2 i + 3 j + 4 k ) =  -0.160779306 + 0.522153864 i + 0.783230794 j + 1.044307726 k

6°)  Cosine & Hyperbolic Cosine Integrals

Ci(q) & Chi(q) are computed by

Ci(q)  = C + ln(q) + Sumn=1,2,..... (-1)n q2n/(2n.(2n)!)
Chi(q)= C + ln(q) + Sumn=1,2,.....  q2n/(2n.(2n)!)                 where  C = 0.5772156649... =  Euler's constant.

Data Registers:  R00 thru R12: temp

When the program stops, the result is in R13 to R16

Flag:   F01           CF 01  to calculate Ci(q)
SF 01  to calculate  Chi(q)

Subroutines:     "LNQ"  "Q*Q"  ( cf "Quaternions for the HP-41" )

 01  LBL "CIQ"   02  STO 01   03  STO 05   04  RDN   05  STO 02   06  STO 06   07  RDN   08  STO 03   09  STO 07   10  RDN   11  STO 04   12  STO 08   13  RDN   14  XEQ "LNQ"   15  STO 13   16  RDN   17  STO 14   18  RDN   19  STO 15   20  X<>Y   21  STO 16           22  .5772156649   23  ST+ 13   24  XEQ "Q*Q"   25  STO 01 26  STO 09   27  RDN   28  STO 02   29  STO 10   30  RDN   31  STO 03   32  STO 11   33  X<>Y   34  STO 04   35  STO 12    36  4   37  FC? 01   38  CHS   39  ST/ 01   40  ST/ 02   41  ST/ 03   42  ST/ 04   43  RCL 01           44  ST+ 13   45  RCL 02   46  ST+ 14   47  RCL 03   48  ST+ 15   49  RCL 04   50  ST+ 16 51  1   52  STO 00   53  LBL 01   54  RCL 09   55  STO 05   56  RCL 10   57  STO 06   58  RCL 11   59  STO 07   60  RCL 12   61  STO 08    62  RCL 00   63  ENTER^   64  ST+ X   65  RCL 00   66  1   67  ST+ Z   68  +   69  STO 00           70  X^2   71  ST+ X   72  *   73  /   74  FC? 01   75  CHS 76  ST* 05   77  ST* 06   78  ST* 07   79  ST* 08   80  XEQ "Q*Q"   81  STO 01   82  RDN   83  STO 02   84  RCL 14   85  +   86  STO 14    87  LASTX   88  -   89  ABS   90  X<>Y   91  STO 03   92  RCL 15   93  +   94  STO 15           95  LASTX   96  -   97  ABS   98  +   99  X<>Y 100  STO 04 101  RCL 16 102  + 103  STO 16 104  LASTX 105  - 106  ABS 107  + 108  RCL 01 109  RCL 13  110  + 111  STO 13 112  LASTX 113  - 114  ABS 115  + 116  X#0? 117  GTO 01         118  RCL 16 119  RCL 15 120  RCL 14 121  RCL 13 122  END

( 175 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Ci (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    CF 01
and      Chi (x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    SF 01

Example:              q = 1 + 2 i + 3 j + 4 k

•   CF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "CIQ"   >>>>    19.04999179                     ---Execution time = 56s---
RDN     -6.098016773
RDN     -9.147025157
RDN   -12.19603355

So,   Ci ( 1 + 2 i + 3 j + 4 k ) =  19.04999179 - 6.098016773 i - 9.147025157 j - 12.19603355 k

•   SF 01

4   ENTER^
3   ENTER^
2   ENTER^
1      R/S             >>>>     -0.220285753                     ---Execution time = 56s---
RDN      0.529901613
RDN      0.794852419
RDN      1.059803225

So,   Chi ( 1 + 2 i + 3 j + 4 k ) =  -0.220285753 + 0.529901613 i + 0.794852419 j + 1.059803225 k

7°)  Bessel Functions - Real order

a) Bessel Functions of the 1st kind

-The following routine returns Jn(q) if F01 is clear or In(q) if F01 is set,
where n is a real number and q = x + y i + z j + t k  is a quaternion.

-The formulas are the same as those for real or complex arguments.

Data Registers:      •  R00 = n                                  ( Register R00 is to be initialized before executing "JNQ" )

R01 thru R08: temp
R09 thru R12 = Partial sums and eventually Jn(q) or In(q)
R13 thru R16 = (q/2)2       R17 = index

Flag:  F01                   CF 01  to compute  Jn(q)
SF 01  to compute   In(q)
Subroutines:    "1/G+"  or   "GAM+"    ( cf "Gamma Function for the HP-41" )
"Q*Q"  "Q^R"              ( cf "Quaternions for the HP-41" )

-Lines 46 to 54 may be replaced by

XEQ "GAM+"
ST/ 01
ST/ 02
ST/ 03
ST/ 04
ST/ 09
ST/ 10
ST/ 11
ST/ 12

 01  LBL "JNQ"   02  STO 01   03  STO 05   04  CLX   05  2   06  ST/ 01   07  ST/ 05   08  ST/ T   09  ST/ Z   10  /   11  STO 02   12  STO 06   13  RDN   14  STO 03   15  STO 07   16  X<>Y   17  STO 04   18  STO 08          19  XEQ "Q*Q"   20  STO 13   21  RDN   22  STO 14   23  RDN   24  STO 15 25  X<>Y   26  STO 16   27  RCL 04   28  RCL 03   29  RCL 02   30  RCL 01   31  XEQ "Q^R"   32  STO 01   33  STO 09   34  RDN   35  STO 02   36  STO 10   37  RDN   38  STO 03   39  STO 11   40  X<>Y   41  STO 04          42  STO 12   43  RCL 00   44  1   45  +   46  XEQ "1/G+"   47  ST* 01    48  ST* 02 49  ST* 03    50  ST* 04    51  ST* 09    52  ST* 10    53  ST* 11    54  ST* 12    55  CLX   56  STO 17   57  LBL 01   58  ISG 17   59  CLX   60  RCL 13    61  STO 05   62  RCL 14   63  STO 06          64  RCL 15   65  STO 07   66  RCL 16   67  STO 08   68  RCL 00   69  RCL 17   70  ST+ Y   71  FC? 01   72  CHS 73  *   74  ST/ 05   75  ST/ 06   76  ST/ 07   77  ST/ 08   78  XEQ "Q*Q"   79  STO 01   80  RDN   81  STO 02   82  RCL 10   83  +   84  STO 10          85  LASTX   86  -   87  ABS   88  X<>Y   89  STO 03   90  RCL 11   91  +   92  STO 11   93  LASTX   94  -   95  ABS   96  + 97  X<>Y   98  STO 04   99  RCL 12 100  + 101  STO 12  102  LASTX 103  - 104  ABS 105  + 106  RCL 01 107  RCL 09        108  + 109  STO 09 110  LASTX 111  - 112  ABS 113  + 114  X#0? 115  GTO 01 116  RCL 12 117  RCL 11 118  RCL 10 119  RCL 09 120  END

( 169 bytes / SIZE 018 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Jn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    CF 01            with       R00 = n
and        In(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    SF 01             with       R00 = n

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k                  PI   STO 00

•   CF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "JNQ"   >>>>   -11.22987514 = R09                         ---Execution time = 52s---
RDN    -3.570801501 = R10
RDN    -5.356202254 = R11
RDN    -7.141603004 = R12

So,   JPI( 1 + 2 i + 3 j + 4 k ) = -11.22987514 - 3.570801501 i - 5.356202254 j - 7.141603004 k

•   SF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "JNQ"   >>>>     0.315197913 = R09                         ---Execution time = 56s---
RDN    -0.129988660 = R10
RDN    -0.194982989 = R11
RDN    -0.259977319 = R12

Whence,   IPI( 1 + 2 i + 3 j + 4 k ) = 0.315197913 - 0.129988660 i - 0.194982989 j - 0.259977319 k

Note:

-If the order is a complex number or a quaternion, several definitions could be used,
but they would not lead to the same results because q^q' has 2 definitions and the multiplication is not commutative.

b) Bessel Functions of the 2nd kind - Non-Integer Order

-This program returns Yn(q) if F01 is clear or Kn(q) if F01 is set,
where n is a real number and q = x + y i + z j + t k  is a quaternion.

-The formulas are the same as those for real or complex arguments.

Data Registers:      •  R00 = n                                       ( Register R00 is to be initialized before executing "YNQ" )

R01 thru R17: temp
R18 thru R21 = q
R22 thru R25 = Yn(q) or Kn(q)

Flag:  F01                   CF 01  to compute  Yn(q)
SF 01  to compute   Kn(q)
Subroutine:    "JNQ"    ( cf § 7-a) )

 01  LBL "YNQ" 02  STO 18       03  RDN 04  STO 19 05  RDN 06  STO 20 07  RDN 08  STO 21 09  RDN 10  XEQ "JNQ" 11  STO 22 12  RDN 13  STO 23 14  RDN 15  STO 24 16  X<>Y 17  STO 25 18  FS? 01 19  GTO 00       20  RCL 00  21  1 22  CHS 23  ACOS 24  * 25  COS 26  ST* 22 27  ST* 23 28  ST* 24 29  ST* 25 30  LBL 00 31  RCL 00 32  CHS 33  STO 00       34  RCL 21 35  RCL 20 36  RCL 19 37  RCL 18 38  XEQ "JNQ" 39  ST- 22 40  RDN 41  ST- 23 42  RDN 43  ST- 24 44  X<>Y 45  ST- 25 46  PI 47  2 48  / 49  CHS 50  FC? 01 51  ST/ X 52  RCL 00       53  CHS 54  STO 00 55  1 56  CHS 57  ACOS 58  * 59  SIN 60  / 61  ST* 22 62  ST* 23 63  ST* 24 64  ST* 25 65  RCL 25  66  RCL 24        67  RCL 23 68  RCL 22 69  END

( 117 bytes / SIZE 026 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Yn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    CF 01            with       R00 = n
and        Kn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    SF 01            with       R00 = n

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k                  PI   STO 00

•   CF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "YNQ"   >>>>    9.617564044 = R22                         ---Execution time = 120s---
RDN    -4.171358834 = R23
RDN    -6.257038260 = R24
RDN    -8.342717673 = R25

So,   YPI( 1 + 2 i + 3 j + 4 k ) = 9.617564044 - 4.171358834 i - 6.257038260 j - 8.342717673 k

•   SF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "YNQ"   >>>      0.203067070 = R22                         ---Execution time = 126s---
RDN    -0.055030894 = R23
RDN    -0.082546338 = R24
RDN    -0.110061786 = R25

Whence,   KPI( 1 + 2 i + 3 j + 4 k ) = 0.203067070 - 0.055030894 i - 0.082546338 j - 0.110061786 k

Note:

-This program does not work if n is an integer.

c) Bessel Functions of the 2nd kind - Integer Order

-"YNQ1" returns Yn(q) if F01 is clear or Kn(q) if F01 is set,
where n is a non-negative integer and q = x + y i + z j + t k  is a quaternion.

-The formulas are the same as those for real arguments.

Data Registers:      •  R00 = n                                        ( Register R00 is to be initialized before executing "YNQ1" )

R01 thru R27: temp

R13 thru R16 = Yn(q) or Kn(q)

Flag:  F01                   CF 01  to compute  Yn(q)
SF 01  to compute   Kn(q)

Subroutines:    "JNQ"                           ( cf § 7-a) )
"LNQ"  "Q^R"  "Q*Q"  ( cf "Quaternions for the HP-41" )

-Line 141 =  2 x Euler's Constant

 01  LBL "YNQ1"   02  STO 18           03  RDN   04  STO 19   05  RDN   06  STO 20   07  RDN   08  STO 21   09  RDN   10  XEQ "JNQ"   11  STO 05   12  RDN   13  STO 06   14  RDN   15  STO 07   16  X<>Y   17  STO 08   18  2   19  ST/ 18   20  ST/ 19   21  ST/ 20   22  ST/ 21   23  RCL 21   24  RCL 20   25  RCL 19   26  RCL 18   27  XEQ "LNQ"   28  STO 01   29  RDN   30  STO 02   31  RDN   32  STO 03   33  X<>Y   34  STO 04   35  1   36  FS? 01   37  CHS   38  STO 27   39  RCL 00   40  STO 17   41  STO 22   42  Y^X   43  ST+ X   44  ST* 01   45  ST* 02   46  ST* 03   47  ST* 04   48  XEQ "Q*Q"   49  STO 13   50  RDN 51  STO 14   52  RDN   53  STO 15           54  X<>Y   55  STO 16   56  RCL 18   57  STO 01   58  STO 05   59  RCL 19   60  STO 02   61  STO 06   62  RCL 20   63  STO 03   64  STO 07   65  RCL 21   66  STO 04   67  STO 08   68  XEQ "Q*Q"   69  STO 23   70  RDN   71  STO 24   72  RDN   73  STO 25   74  X<>Y   75  STO 26   76  CLX   77  STO 05   78  STO 06   79  STO 07   80  STO 08   81  LBL 01   82  RCL 17   83  RCL 22   84  -   85  FACT   86  RCL 22   87  1   88  -   89  X<0?   90  CLST   91  STO 00   92  FACT   93  /   94  RCL 27   95  RCL 00   96  Y^X   97  *   98  STO 01   99  RCL 26 100  RCL 25 101  RCL 24 102  RCL 23         103  XEQ "Q^R" 104  X<> 01 105  ST* 01 106  ST* Z 107  ST* T 108  * 109  ST+ 06 110  RDN 111  ST+ 07 112  X<>Y 113  ST+ 08 114  RCL 01 115  ST+ 05 116  DSE 22 117  GTO 01 118  RCL 17 119  CHS 120  STO 00 121  RCL 21 122  RCL 20 123  RCL 19 124  RCL 18 125  XEQ "Q^R" 126  STO 01 127  RDN 128  STO 02 129  RDN 130  STO 03 131  X<>Y 132  STO 04 133  XEQ "Q*Q" 134  ST- 13 135  RDN 136  ST- 14 137  RDN 138  ST- 15 139  X<>Y 140  ST- 16 141  1.15443133 142  X<> 17 143  STO 00 144  ABS 145  LBL 02 146  X#0? 147  1/X 148  ST- 17 149  X<> L 150  DSE X 151  GTO 02 152  RCL 23         153  RCL 27 154  ST* 18 155  ST* 19 156  ST* 20 157  ST* 21 158  CHS 159  * 160  STO 05 161  RCL 24 162  LASTX 163  * 164  STO 06 165  RCL 25 166  LASTX 167  * 168  STO 07 169  RCL 26 170  LASTX 171  * 172  STO 08 173  CLX 174  STO 02 175  STO 03 176  STO 04 177  STO 10 178  STO 11 179  STO 12 180  SIGN 181  STO 22 182  RCL 17 183  RCL 00 184  FACT 185  1/X 186  STO 01 187  * 188  STO 09 189  LBL 03 190  XEQ "Q*Q" 191  STO 01 192  RDN 193  STO 02 194  RDN 195  STO 03 196  X<>Y 197  STO 04 198  RCL 22 199  STO Y 200  RCL 00 201  + 202  * 203  ST/ 01 204  ST/ 02 205  ST/ 03 206  ST/ 04 207  RCL 17         208  LASTX 209  1/X 210  - 211  RCL 22 212  1/X 213  - 214  STO 17 215  ISG 22 216  CLX 217  RCL 01 218  * 219  RCL 09 220  + 221  STO 09 222  LASTX 223  - 224  ABS 225  RCL 02 226  RCL 17 227  * 228  RCL 10 229  + 230  STO 10 231  LASTX 232  - 233  ABS 234  + 235  RCL 03 236  RCL 17 237  * 238  RCL 11 239  + 240  STO 11 241  LASTX 242  - 243  ABS 244  + 245  RCL 04 246  RCL 17 247  * 248  RCL 12 249  + 250  STO 12 251  LASTX 252  - 253  ABS 254  + 255  X#0? 256  GTO 03 257  RCL 09 258  STO 05 259  RCL 10         260  STO 06 261  RCL 11 262  STO 07 263  RCL 12 264  STO 08 265  RCL 21 266  RCL 20 267  RCL 19 268  RCL 18 269  XEQ "Q^R" 270  STO 01 271  RDN 272  STO 02 273  RDN 274  STO 03 275  X<>Y 276  STO 04 277  XEQ "Q*Q" 278  ST+ 13 279  RDN 280  ST+ 14 281  RDN 282  ST+ 15 283  X<>Y 284  ST+ 16 285  PI 286  FS? 01 287  -2 288  ST/ 13 289  ST/ 14 290  ST/ 15 291  ST/ 16 292  RCL 16 293  RCL 15 294  RCL 14 295  RCL 13 296  END

( 456 bytes / SIZE 028 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Yn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    CF 01            with       R00 = n   is a non-negative integer
and        Kn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    SF 01            with       R00 = n   is a non-negative integer

Example:         n = 3     q = 1 + 2 i + 3 j + 4 k                  3   STO 00

•   CF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "YNQ"   >>>>    7.690027238 = R13                         ---Execution time = 141s---
RDN    -5.224812676 = R14
RDN    -7.837219010 = R15
RDN  -10.44962535   = R16

So,   Y3( 1 + 2 i + 3 j + 4 k ) = 7.690027238 - 5.224812676 i - 7.837219010 j - 10.44962535 k

•   SF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "YNQ"   >>>      0.208767799 = R13                         ---Execution time = 152s---
RDN    -0.047983196 = R14
RDN    -0.071974795 = R15
RDN    -0.095966395 = R16

Whence,   K3( 1 + 2 i + 3 j + 4 k ) = 0.208767799 - 0.047983196 i - 0.071974795 j - 0.095966395 k

-Likewise, we find:

Y0( 1 + 2 i + 3 j + 4 k ) = 29.92977303 + 8.767829494 i + 13.15174424 j + 17.53565898 k    ( in 132 seconds )
K0( 1 + 2 i + 3 j + 4 k ) = 0.190884051 + 0.016289913 i + 0.024434869 j + 0.032579825 k    ( in 147 seconds )

Note:

-If  n < 0 , we can use  Y-n(q)  =  (-1)n  Yn(q)

8°)  Kummer's Function

Formula:        M(a;b;q) = 1 +  (a)1/(b)1. q1/1! + ............. +  (a)n/(b)n . qn/n! + ..........  where (a)n = a(a+1)(a+2) ...... (a+n-1)

where   a & b  are real numbers

Data Registers:    •  R09 = a              ( Registers R09 & R10 are to be initialized before executing "KUMQ" )
•  R10 = b

The result is in registers  R11 thru R14 at the end  ,  q  is in R01 thru R04
Flags: /
Subroutine:  "Q*Q"  ( cf "Quaternions for the HP-41" )

 01  LBL "KUMQ" 02  STO 01          03  RDN 04  STO 02 05  RDN 06  STO 03 07  X<>Y 08  STO 04 09  CLX 10  STO 00 11  STO 06 12  STO 07 13  STO 08 14  STO 12 15  STO 13 16  STO 14 17  SIGN 18  STO 05 19  STO 11 20  LBL 01          21  XEQ "Q*Q"  22  STO 05 23  RDN 24  STO 06 25  RDN 26  STO 07 27  X<>Y 28  STO 08 29  RCL 09 30  RCL 10 31  RCL 00 32  ST+ Y 33  ST+ Z 34  1 35  + 36  STO 00          37  * 38  / 39  ST* 05 40  ST* 06 41  ST* 07 42  ST* 08 43  RCL 05 44  RCL 11 45  + 46  STO 11 47  LASTX 48  - 49  ABS 50  RCL 06 51  RCL 12          52  + 53  STO 12  54  LASTX 55  - 56  ABS 57  + 58  RCL 07 59  RCL 13 60  + 61  STO 13 62  LASTX 63  - 64  ABS 65  + 66  RCL 08 67  RCL 14          68  + 69  STO 14  70  LASTX 71  - 72  ABS 73  + 74  X#0? 75  GTO 01 76  RCL 14 77  RCL 13 78  RCL 12 79  RCL 11 80  END

( 100 bytes / SIZE 015 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where     M ( a , b ,  x+y.i+z.j+t.k )  = x'+y'.i+z'.j+t'.k     provided    R09 = a  &  R10 = b

Example:              q = 1 + 2 i + 3 j + 4 k   ,   a = sqrt(2)  ,  b = PI

2   SQRT   STO 09     PI   STO 10

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "KUMQ"   >>>>   -0.533892203                     ---Execution time =99s---
RDN      0.068805302
RDN      0.103207953
RDN      0.137610604

Whence,   M ( sqrt(2) , PI ; 1 + 2 i + 3 j + 4 k ) =  -0.533892203 + 0.068805302 i + 0.103207953 j + 0.137610604 k

9°)  Parabolic Cylinder Functions

a)  Via Kummer's Function

Formula:         Dn(q) = 2n/2 Pi1/2 exp(-q2/4) [ 1/Gam((1-n)/2) M( -n/2 , 1/2 , q2/2 ) - 21/2 ( q / Gam(-n/2) ) M [ (1-n)/2 , 3/2 , q2/2 ]

Data Registers:           •  R00 = n              ( Register R00 is to be initialized before executing "DNQ" )

R01 thru R23: temp
Flags: /
Subroutines:    "KUMQ"
"Q*Q"  "E^Q"  ( cf "Quaternions for the HP-41" )
"1/G" or "1/G+"  ( cf "Gamma Function for the HP-41" )

 01  LBL "DNQ"    02  STO 01   03  STO 05             04  STO 15   05  RDN   06  STO 02   07  STO 06   08  STO 16   09  RDN   10  STO 03   11  STO 07   12  STO 17   13  X<>Y   14  STO 04   15  STO 08   16  STO 18   17  RCL 00   18  STO 19   19  CHS   20  .5   21  ST* 05   22  ST* 06   23  ST* 07   24  ST* 08   25  STO 10   26  * 27  STO 09   28  XEQ "Q*Q"   29  XEQ "KUMQ"   30  STO 20             31  RDN   32  STO 21   33  RDN   34  STO 22   35  X<>Y   36  STO 23   37  RCL 10   38  ST+ 09   39  RCL 09   40  XEQ "1/G"   41  ST* 20   42  ST* 21   43  ST* 22   44  ST* 23   45  ISG 10   46  RCL 04   47  RCL 03   48  RCL 02   49  RCL 01   50  XEQ "KUMQ"   51  STO 05   52  RDN 53  STO 06   54  RDN   55  STO 07   56  X<>Y   57  STO 08   58  RCL 19             59  STO 00   60  2   61  /   62  CHS   63  XEQ "1/G"    64  2   65  SQRT   66  CHS   67  *   68  ST* 05   69  ST* 06   70  ST* 07   71  ST* 08   72  RCL 01   73  X<> 15   74  STO 01   75  RCL 02   76  X<> 16   77  STO 02   78  RCL 03 79  X<>17   80  STO 03   81  RCL 04   82  X<> 18   83  STO 04   84  XEQ "Q*Q"    85  STO 05             86  CLX   87  RCL 21   88  +   89  STO 06   90  CLX   91  RCL 22   92  +   93  STO 07   94  CLX   95  RCL 23   96  +   97  STO 08   98  RCL 20   99  ST+ 05 100  2 101  CHS 102  ST/ 15 103  ST/ 16 104  ST/ 17 105  ST/ 18 106  RCL 18 107  RCL 17 108  RCL 16 109  RCL 15 110  XEQ "E^Q"  111  STO 01           112  RDN 113  STO 02 114  RDN 115  STO 03 116  X<>Y 117  STO 04 118  2 119  RCL 00 120  Y^X 121  PI 122  * 123  SQRT 124  ST* 01 125  ST* 02 126  ST* 03 127  ST* 04 128  XEQ "Q*Q" 129  END

( 215 bytes / SIZE 024 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Dn(x+y.i+z.j+t.k) = x'+y'.i+z'.j+t'.k                  with       R00 = n

Example:         n = 0.4      q = 1 + 1.1 i + 1.2 j + 1.3 k

0.4   STO 00

1.3   ENTER^
1.2   ENTER^
1.1   ENTER^
1     XEQ "DNQ"   >>>>    2.606105105                         ---Execution time = 151s---
RDN    -0.969116887
RDN    -1.057218423
RDN    -1.145319957

D0.4( 1 + 1.1 i + 1.2 j + 1.3 k ) = 2.606105105 - 0.969116887 i - 1.057218423 j - 1.145319957 k

b)  Asymptotic Expansion  ( q Large )

Formula:    Dn(q)  ~  qn exp(-q2/4) [ 1 - n(n-1) / ( 2 q2 ) + n(n-1)(n-2)(n-3) / ( 2 ( 4 q4 ) ) - ....... ]

Data Registers:           •  R00 = n              ( Register R00 is to be initialized before executing "AEDNQ" )

R01 thru R21: temp
Flags: /
Subroutines:    "Q*Q"   "1/Q"   "E^Q"   "Q^R"   ( cf "Quaternions for the HP-41" )

 01  LBL "AEDNQ"   02  STO 01              03  STO 05   04  STO 13   05  RDN   06  STO 02   07  STO 06   08  STO 14   09  RDN   10  STO 03   11  STO 07   12  STO 15   13  X<>Y   14  STO 04   15  STO 08   16  STO 16   17  XEQ "Q*Q"   18  STO 17   19  RDN   20  STO 18   21  RDN   22  STO 19   23  RDN   24  STO 20   25  RDN   26  XEQ "1/Q"   27  STO 01   28  RDN   29  STO 02   30  RDN   31  STO 03 32  X<>Y   33  STO 04              34  CLX   35  STO 06   36  STO 07   37  STO 08   38  STO 10   39  STO 11   40  STO 12   41  STO 21   42  SIGN   43  STO 05   44  STO 09   45  LBL 01   46  XEQ "Q*Q"    47  STO 05   48  RDN   49  STO 06   50  RDN   51  STO 07   52  X<>Y   53  STO 08   54  2   55  RCL 00   56  ST+ Y   57  1   58  ST+ 21   59  +   60  RCL 21   61  ST+ X   62  ST- Y 63  ST- Z   64  /   65  *   66  CHS   67  ST* 05   68  ST* 06   69  ST* 07   70  ST* 08   71  RCL 05              72  RCL 09   73  +   74  STO 09   75  LASTX   76  -   77  ABS   78  RCL 06    79  RCL 10   80  +   81  STO 10   82  LASTX   83  -   84  ABS   85  +   86  RCL 07   87  RCL 11   88  +   89  STO 11   90  LASTX   91  -   92  ABS   93  + 94  RCL 08   95  RCL 12   96  +   97  STO 12              98  LASTX   99  - 100  ABS 101  + 102  X#0? 103  GTO 01 104  4 105  CHS 106  ST/ 17 107  ST/ 18 108  ST/ 19 109  ST/ 20 110  RCL 20 111  RCL 19 112  RCL 18 113  RCL 17 114  XEQ "E^Q"  115  STO 05 116  RDN 117  STO 06 118  RDN 119  STO 07 120  X<>Y 121  STO 08 122  RCL 16 123  RCL 15 124  RCL 14 125  RCL 13            126  XEQ "Q^R" 127  STO 01 128  RDN 129  STO 02 130  RDN 131  STO 03 132  X<>Y 133  STO 04 134  XEQ "Q*Q"  135  STO 01 136  RDN 137  STO 02 138  RDN 139  STO 03 140  X<>Y 141  STO 04 142  RCL 09 143  STO 05 144  RCL 10 145  STO 06 146  RCL 11 147  STO 07 148  RCL 12 149  STO 08 150  XEQ "Q*Q" 151  END

( 215 bytes / SIZE 022 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Dn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                  with       R00 = n

Example:         n = PI       q = 1 + 2 i + 3 j + 4 k

PI   STO 00

4    ENTER^
3    ENTER^
2    ENTER^
1    XEQ "AEDNQ"   >>>>   -33138.87293                     ---Execution time =50s---
RDN      93318.15137
RDN      139977.2270
RDN      186636.3025

DPI( 1 + 2 i + 3 j + 4 k ) =  -33138.87293 + 93318.15137 i + 139977.2270 j + 186636.3025 k

Notes:

-As usual, infinte loops will occur if q is too small.
-Add  RND  after line 101 to get a less accurate result in this case.
-The number of exact digits will be controlled by the display settings.

-This program may also be used if q is relatively small when n is a positive integer.

10°)  Struve Functions

Formulae:

Hn(q) =  (q/2)n+1 SUM k>=0  (-1)k(q/2)2k/ ( Gam(k+3/2).Gam(k+n+3/2) )
Ln(q) =  (q/2)n+1 SUM k>=0  (q/2)2k/ ( Gam(k+3/2).Gam(k+n+3/2) )

Data Registers:           •  R00 = n              ( Register R00 is to be initialized before executing "HLNQ" )

R01 thru R13: temp        At the end, the result is in registers  R09 thru R12

Flag:                   CF 01  to calculate Hn(q)
SF 01  to calculate  Ln(q)

Subroutines:    "Q*Q"  "Q^R"    ( cf "Quaternions for the HP-41" )
"1/G" or "1/G+"  ( cf "Gamma Function for the HP-41" )

 01  LBL "HLNQ"   02  STO 01           03  STO 05   04  CLX   05  2   06  ST/ 01   07  ST/ 05   08  ST/ T   09  ST/ Z   10  /   11  STO 02   12  STO 06   13  RDN   14  STO 03   15  STO 07   16  X<>Y   17  STO 04   18  STO 08   19  XEQ "Q*Q"   20  STO 05   21  RDN   22  STO 06   23  RDN   24  STO 07   25  X<>Y 26  STO 08   27  RCL 00   28  STO 13           29  1   30  ST+ 00   31  RCL 04   32  RCL 03   33  RCL 02   34  RCL 01   35  XEQ "Q^R"    36  STO 01   37  STO 09   38  RDN   39  STO 02   40  STO 10   41  RDN   42  STO 03   43  STO 11   44  X<>Y   45  STO 04   46  STO 12   47  RCL 13   48  STO 00   49  1.5   50  + 51  XEQ "1/G"   52  ST+ X   53  PI   54  SQRT   55  /   56  ST* 01   57  ST* 02   58  ST* 03   59  ST* 04   60  ST* 09   61  ST* 10   62  ST* 11   63  ST* 12   64  .5   65  STO 13           66  LBL 01   67  XEQ "Q*Q"    68  STO 01   69  RDN   70  STO 02   71  RDN   72  STO 03   73  X<>Y   74  STO 04   75  RCL 00 76  RCL 13   77  1   78  +   79  STO 13   80  ST+ Y   81  *   82  FC? 01   83  CHS   84  ST/ 01   85  ST/ 02   86  ST/ 03   87  ST/ 04   88  RCL 01           89  RCL 09   90  +   91  STO 09    92  LASTX   93  -   94  ABS   95  RCL 02   96  RCL 10   97  +   98  STO 10   99  LASTX 100  - 101  ABS 102  + 103  RCL 03 104  RCL 11 105  + 106  STO 11 107  LASTX 108  - 109  ABS 110  + 111  RCL 04 112  RCL 12         113  + 114  STO 12 115  LASTX  116  - 117  ABS 118  + 119  X#0? 120  GTO 01 121  RCL 12 122  RCL 11 123  RCL 10 124  RCL 09 125  END

( 174 bytes / SIZE 014 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Hn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    CF 01            with       R00 = n   &  n+3/2 # 0 , -1 , -2 , ..........
and        Ln(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k       if    SF 01             with       R00 = n   &  n+3/2 # 0 , -1 , -2 , ..........

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k                  PI   STO 00

•   CF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "HLNQ"   >>>>    8.546633222 = R09                         ---Execution time = 48s---
RDN    -4.085079931 = R10
RDN    -6.127619899 = R11
RDN    -8.170159859 = R12

So,   Hpi( 1 + 2 i + 3 j + 4 k ) = 8.546633222 - 4.085079931 i - 6.127619899 j - 8.170159859 k

•   SF 01

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "HLNQ"   >>>      1.864582784 = R09                         ---Execution time = 51s---
RDN    -0.116220199 = R10
RDN    -0.174330298 = R11
RDN    -0.232440396 = R12

Whence,   Lpi( 1 + 2 i + 3 j + 4 k ) = 1.864582784 - 0.116220199 i - 0.174330298 j - 0.232440396 k

11°)  Hypergeometric Function

-This program computes  2F1(a,b;c;q) = 1 +  (a)1(b)1/(c)1. q1/1! + ............. +  (a)n(b)n/(c)n . qn/n! + ..........

where  a , b , c  are real numbers  and  q is a quaternion.         [  (a)n = a(a+1)(a+2) ...... (a+n-1)   is Pochhammer's symbol  ]

-If , however, F does not exist,  the regularized hypergeometric function F tilde is returned ( cf "Hypergeometric Functions for the HP-41" )
-In this case, flag F00 is set.

Data Registers:                                       ( Registers R09-R10-R11 are to be initialized before executing "HGFQ" )

•  R09 = a
•  R10 = b
•  R11 = c

The result is in registers  R12 thru R15 at the end  ,  q  is in R01 thru R04
Flag:   F00
Subroutines:  "Q*Q"  "Q^R"  ( cf "Quaternions for the HP-41" )

 01  LBL "HGFQ"   02  CF 00   03  STO 01           04  RDN   05  STO 02   06  RDN   07  STO 03   08  X<>Y   09  STO 04   10  CLX   11  STO 00   12  STO 06   13  STO 07   14  STO 08   15  SIGN   16  STO 05   17  RCL 11   18  ENTER^   19  FRC   20  X=0?   21  XY?   32  GTO 02 33  LBL 00   34  RCL 10   35  ENTER^   36  FRC   37  X=0?   38  XY?   43  GTO 02   44  LBL 00   45  SF 00   46  1   47  RCL 11   48  -   49  STO 00   50  RCL 04   51  RCL 03   52  RCL 02   53  RCL 01   54  XEQ "Q^R"   55  STO 05   56  RDN   57  STO 06   58  RDN   59  STO 07   60  X<>Y   61  STO 08   62  RCL 00   63  1   64  LBL 01 65  RCL 09   66  1   67  -   68  R^   69  +   70  *   71  RCL 10           72  R^   73  +   74  1   75  -   76  *   77  R^   78  /   79  DSE Y   80  GTO 01   81  ST* 05   82  ST* 06   83  ST* 07   84  ST* 08   85  LBL 02   86  RCL 05   87  STO 12   88  RCL 06   89  STO 13   90  RCL 07   91  STO 14   92  RCL 08   93  STO 15   94  LBL 03   95  XEQ "Q*Q"   96  STO 05 97  RDN   98  STO 06    99  RDN 100  STO 07 101  X<>Y 102  STO 08 103  RCL 09         104  RCL 10 105  RCL 11 106  RCL 00 107  ST+ Y 108  ST+ Z 109  ST+ T 110  ISG X 111  ""  112  STO 00 113  * 114  / 115  * 116  ST* 05 117  ST* 06 118  ST* 07 119  ST* 08 120  RCL 05 121  RCL 12 122  + 123  STO 12 124  LASTX 125  - 126  ABS 127  RCL 06 128  RCL 13 129  + 130  STO 13  131  LASTX 132  - 133  ABS 134  + 135  RCL 07 136  RCL 14         137  + 138  STO 14 139  LASTX 140  - 141  ABS 142  + 143  RCL 08 144  RCL 15 145  + 146  STO 15 147  LASTX 148  - 149  ABS 150  + 151  X#0? 152  GTO 03 153  RCL 15 154  RCL 14 155  RCL 13 156  RCL 12 157  END

( 196 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where     x'+y'.i+z'.j+t'.k = 2F1 ( a , b ; c ;  x+y.i+z.j+t.k )                       if   CF 00                with             R09 = a  ,  R10 = b  ,  R11 = c
or            2F~1 ( a , b ; c ;  x+y.i+z.j+t.k )                      if   SF 00

Example1:              q = 0.1 + 0.2 i + 0.3 j + 0.4 k           a = 1.1   b = 1.2   c = 1.3

1.1   STO 09    1.2   STO 10   1.3   STO 11

0.4   ENTER^
0.3   ENTER^
0.2   ENTER^
0.1   XEQ "HGFQ"   >>>>    0.814236592 = R12                         ---Execution time = 2mn25s---
RDN     0.184421233 = R13
RDN     0.276631850 = R14
RDN     0.368842467 = R15                             Flag F00 is clear

2F1 ( 1.1 , 1.2 ; 1.3 ;  0.1 + 0.2 i + 0.3 j + 0.4 k ) = 0.814236592 + 0.184421233 i + 0.276631850 j + 0.368842467 k

Example2:              q = 0.1 + 0.2 i + 0.3 j + 0.4 k           a = 1   b = 2   c = -4

1   STO 09    2   STO 10   4  CHS   STO 11

0.4   ENTER^
0.3   ENTER^
0.2   ENTER^
0.1   XEQ "HGFQ"   >>>>  -7.152496406 = R12                         ---Execution time = 4mn35s---
RDN   -9.061272386 = R13
RDN   -13.59190854 = R14
RDN   -18.12254485 = R15                             Flag F00 is set
~
2F1 ( 1 , 2 ; -4 ;  0.1 + 0.2 i + 0.3 j + 0.4 k ) = -7.152496406 - 9.061272386 i - 13.59190854 j - 18.12254485 k

Note:

-In general, the hypergeometric series is convergent if  | q | < 1

12°)  Associated Legendre Functions of the 1st kind

a)  Program#1 - Functions of Type 2 / Type 3

Formulae:

Type 2      Pnm(q) = [ (q+1)/(1-q) ]m/2  2F1(-n , n+1 ; 1-m ; (1-q)/2 ) / Gamma(1-m)        (  q # 1 )

Type 3      Pnm(q) = [ (q+1)/(q-1) ]m/2  2F1(-n , n+1 ; 1-m ; (1-q)/2 ) / Gamma(1-m)        (  q # 1 )

Data Registers:    •  R09 = m              ( Registers R09 &  R10 are to be initialized before executing "ALF1Q" )
•  R10 = n

Flags:   F03 & F00 ( which is used by "HGFQ" )

Clear F03 to compute the function of type 2
Set   F03 to compute the function of type 3

Subroutines:

"HGFQ"
"Q*Q"  "Q^R"  "1/Q"  ( cf "Quaternions for the HP-41" )
"1/G+"  or  "GAM+"   ( cf  "Gamma Function" )

 01  LBL "ALF1Q" 02  STO 16          03  RDN 04  STO 17 05  RDN 06  STO 18 07  X<>Y 08  STO 19 09  RCL 10 10  STO 20 11  CHS 12  X<> 09 13  STO 21 14  1 15  ST+ 10 16  X<>Y 17  - 18  STO 11 19  RCL 16 20  1 21  - 22  SIGN 23  RCL 19          24  RCL 18 25  RCL 17 26  2 27  CHS 28  ST/ L 29  ST/ Y 30  ST/ Z 31  ST/ T 32  X<> L 33  XEQ "HGFQ" 34  FS?C 00 35  GTO 00 36  RCL 11 37  XEQ "1/G+" 38  ST* 12 39  ST* 13 40  ST* 14 41  ST* 15 42  LBL 00 43  RCL 16          44  STO 05 45  1 46  ST+ 05 47  - 48  FC? 03 49  CHS 50  RCL 19 51  STO 08 52  FC? 03 53  CHS 54  RCL 18 55  STO 07 56  FC? 03 57  CHS 58  RCL 17 59  STO 06 60  FC? 03 61  CHS 62  R^ 63  XEQ "1/Q"  64  STO 01 65  RDN 66  STO 02          67  RDN 68  STO 03 69  X<>Y 70  STO 04 71  RCL 20 72  STO 10 73  RCL 21 74  STO 09 75  2 76  / 77  STO 00 78  XEQ "Q*Q" 79  XEQ "Q^R" 80  STO 01 81  RDN 82  STO 02 83  RDN 84  STO 03          85  X<>Y 86  STO 04 87  RCL 12 88  STO 05 89  RCL 13 90  STO 06 91  RCL 14 92  STO 07 93  RCL 15 94  STO 08 95  XEQ "Q*Q"  96  END

( 166 bytes / SIZE 022 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Pmn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                  with       R09 = m   &   R10 = n

Function of type 2 if CF 03
Function of type 3 if SF 03

Example1:      m = sqrt(2)    n = sqrt(3)       q = 0.4 + 0.5 i + 0.6 j + 0.7 k

2  SQRT   STO 09    3  SQRT   STO 10

•  CF 03   Legendre Function of type 2

0.7   ENTER^
0.6   ENTER^
0.5   ENTER^
0.4   XEQ "ALF1Q"   >>>>    -0.866899243         ---Execution time = 3m11s---
RDN    -1.808960473
RDN    -2.170752568
RDN    -2.532544661

Psqrt(2)sqrt(3) ( 0.4 + 0.5 i + 0.6 j + 0.7 k ) = -0.866899243 - 1.808960473 i - 2.170752568 j - 2.532544661 k

•  SF 03   Legendre Function of type 3

0.7   ENTER^
0.6   ENTER^
0.5   ENTER^
0.4   XEQ "ALF1Q"   >>>>    -2.494183066         ---Execution time = 3m11s---
RDN      1.424529609
RDN      1.709435533
RDN      1.994341459

Psqrt(2)sqrt(3) ( 0.4 + 0.5 i + 0.6 j + 0.7 k ) = -2.494183066 + 1.424529609 i + 1.709435533 j + 1.994341459 k

Example2:      m = 1    n = 3       q = 0.4 + 0.5 i + 0.6 j + 0.7 k

1   STO 09    3   STO 10

•  CF 03   Legendre Function of type 2

0.7   ENTER^
0.6   ENTER^
0.5   ENTER^
0.4   XEQ "ALF1Q"   >>>>     10.31801144            ---Execution time = 30s---
RDN    -5.472130005
RDN    -6.566555996
RDN    -7.660982026

P13 ( 0.4 + 0.5 i + 0.6 j + 0.7 k ) = 10.31801144 - 5.472130005 i - 6.566555996 j - 7.660982026 k

•  SF 03   Legendre Function of type 3

0.7   ENTER^
0.6   ENTER^
0.5   ENTER^
0.4   XEQ "ALF1Q"   >>>>   -11.47843674           ---Execution time = 30s---
RDN    -4.918918944
RDN    -5.902702745
RDN    -6.886486535

P13 ( 0.4 + 0.5 i + 0.6 j + 0.7 k ) = -11.47843674 - 4.918918944 i - 5.902702745 j - 6.886486535 k

Note:

-In general, the series is convergent if  | 1 - q | < 2

b)  Integer Order and Degree - Functions of Type 2 / Type 3

Formulae:      (n-m) Pnm(q) = q (2n-1) Pn-1m(q) - (n+m-1) Pn-2m(q)

Type 2                Pmm(q) = (-1)m (2m-1)!!  ( 1-q2 )m/2                    where (2m-1)!! = (2m-1)(2m-3)(2m-5).......5.3.1
Type 3                Pmm(q) =   (2m-1)!!  ( q2-1 )m/2

Data Registers:    •  R09 = m              ( Registers R09 &  R10 are to be initialized before executing "PMNQ" )
•  R10 = n

When the program stops, the result is in R05 to 508

Flag:      F03

Clear F03 to compute the function of type 2
Set   F03 to compute the function of type 3

Subroutines:  "Q*Q"  "Q^R"   ( cf "Quaternions for the HP-41" )

 01  LBL "PMNQ"   02  STO 01            03  STO 05   04  RDN   05  STO 02   06  STO 06   07  RDN   08  STO 03   09  STO 07   10  X<>Y   11  STO 04   12  STO 08   13  RCL 09   14  2   15  /   16  STO 00   17  XEQ "Q*Q"   18  SIGN   19  ST* X   20  ST- L   21  X<> L   22  FC? 03   23  CHS   24  R^   25  FC? 03   26  CHS   27  R^ 28  FC? 03   29  CHS   30  R^   31  FC? 03   32  CHS   33  R^   34  XEQ "Q^R"    35  STO 05            36  RDN   37  STO 06   38  RDN   39  STO 07   40  X<>Y   41  STO 08   42  RCL 09   43  ST+ X   44  1   45  ST+ Y   46  X<>Y   47  LBL 01   48  2   49  -   50  X>0?   51  ST* Y   52  X>0?   53  GTO 01   54  CLX 55  SIGN   56  FC? 03   57  CHS   58  RCL 09            59  Y^X   60  *   61  ST* 05   62  ST* 06   63  ST* 07   64  ST* 08   65  CLX   66  STO 11   67  STO 12   68  STO 13   69  STO 14   70  RCL 10   71   E3   72  /   73  RCL 09   74  +   75  STO 00   76  LBL 02   77  ISG 00   78  FS? 30   79  GTO 00   80  1   81  RCL 09 82  -   83  RCL 00            84  INT   85  -   86  ST* 11   87  ST* 12   88  ST* 13   89  ST* 14   90  LASTX   91  ST+ X   92  1   93  -   94  STO 15   95  XEQ "Q*Q"   96  X<> 15   97  ST* 15   98  ST* T   99  ST* Z 100  * 101  RCL 06 102  X<> 12 103  + 104  STO 06 105  X<> 07 106  X<> 13 107  + 108  STO 07 109  X<> 08 110  X<> 14 111  + 112  STO 08 113  RCL 05          114  X<> 11 115  RCL 15 116  + 117  STO 05 118  RCL 00  119  INT 120  RCL 09 121  - 122  ST/ 05 123  ST/ 06 124  ST/ 07 125  ST/ 08 126  GTO 02 127  LBL 00 128  RCL 08 129  RCL 07 130  RCL 06 131  RCL 05 132  END

( 193 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Pmn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                  with       R09 = m   &  R10 = n         m & n integers   0 <= m <= n

Function of type 2 if CF 03
Function of type 3 if SF 03

Example:      m = 3    n = 7       q = 1 + i/2 + j/3 + k/4

3   STO 09    7   STO 10

•  CF 03   Legendre Function of type 2

4   1/X
3   1/X
2   1/X
1   XEQ "PMNQ"   >>>>   -12188.53940            ---Execution time = 23s---
RDN    -8670.127578
RDN    -5780.085053
RDN    -4335.063783

P37 ( 1 + i/2 + j/3 + k/4 ) = -12188.53940 - 8670.127578 i - 5780.085053 j - 4335.063783 k

•  SF 03   Legendre Function of type 3

4   1/X
3   1/X
2   1/X
1   XEQ "PMNQ"   >>>>    11285.97688            ---Execution time = 23s---
RDN    -9363.495320
RDN    -6242.330210
RDN    -4681.747663

P37 ( 1 + i/2 + j/3 + k/4 ) = 11285.97688 - 9363.495320 i - 6242.330210 j - 4681.747663 k

Notes:

-This program doesn't check if m & n are integers that satisfy 0 <= m <= n
Pmn-1(q) is in registers R11 to R14.

13°)  Associated Legendre Functions of the 2nd kind

a)  Program#1 - Function of Type 2 - | q | < 1

Formula:      Qnm(q) = 2m pi1/2  (1-q2)-m/2 [ -Gam((1+m+n)/2)/(2.Gam((2-m+n)/2)) . sin pi(m+n)/2 .  2F1(-n/2-m/2 ; 1/2+n/2-m/2 ; 1/2 ; q2 )
+ q Gam((2+n+m)/2) / Gam((1+n-m)/2) . cos pi(m+n)/2 . 2F1((1-m-n)/2 ; (2+n-m)/2 ; 3/2 ; q2 )  ]

Data Registers:        •  R09 = m         •  R10 = n         ( Registers R09 &  R10 are to be initialized before executing "ALF2Q2" )

Flags: /
Subroutines:

"HGFQ"
"Q*Q"  "Q^R"  "1/Q"  ( cf "Quaternions for the HP-41" )
"1/G+"  or  "GAM+"   ( cf  "Gamma Function" )

 01  LBL "ALF2Q2"   02  STO 16              03  RDN   04  STO 17   05  RDN   06  STO 18   07  RDN   08  STO 19   09  RDN   10  STO 00   11  CLX   12  2   13  X<> 00   14  XEQ "Q^R"   15  STO 21   16  RDN   17  STO 22   18  RDN   19  STO 23   20  X<>Y   21  STO 24   22  RCL 10   23  STO 25   24  RCL 09   25  STO 20   26  ST- 10   27  +   28  STO 30   29  CHS   30  STO 09   31  .5   32  ST* 09   33  ST* 10   34  ST+ 10 35  STO 11   36  RCL 24   37  RCL 23   38  RCL 22              39  RCL 21   40  XEQ "HGFQ"    41  STO 26   42  RDN   43  STO 27   44  RDN   45  STO 28   46  X<>Y   47  STO 29   48  1   49  ASIN   50  RCL 30   51  *   52  STO 31   53  SIN   54  STO 33   55  RCL 11   56  RCL 09   57  -   58  XEQ "1/G+"   59  ST/ 33   60  RCL 10   61  STO 32   62  RCL 11   63  +   64  STO 10   65  XEQ "1/G+"   66  RCL 33   67  *   68  2 69  /   70  ST* 26   71  ST* 27   72  ST* 28   73  ST* 29   74  RCL 11              75  ST+ 09   76  RCL 24   77  RCL 23   78  RCL 22   79  RCL 21   80  ISG 11   81  XEQ "HGFQ"    82  STO 05   83  RDN   84  STO 06   85  RDN   86  STO 07   87  X<>Y   88  STO 08   89  RCL 16   90  STO 01   91  RCL 17   92  STO 02   93  RCL 18   94  STO 03   95  RCL 19   96  STO 04   97  RCL 30   98  2   99  ST+ Y 100  / 101  XEQ "1/G+" 102  STO 33 103  RCL 32 104  XEQ "1/G+" 105  RCL 31 106  COS 107  * 108  RCL 33            109  / 110  ST* 01 111  ST* 02 112  ST* 03 113  ST* 04 114  XEQ"Q*Q"  115  STO 05 116  CLX 117  RCL 27 118  - 119  STO 06 120  CLX 121  RCL 28 122  - 123  STO 07 124  CLX 125  RCL 29 126  - 127  STO 08 128  RCL 26 129  ST- 05 130  RCL 20 131  STO 09 132  2 133  / 134  CHS 135  STO 00 136  RCL 25 137  STO 10 138  1 139  ST- 21 140  RCL 24            141  CHS 142  RCL 23 143  CHS 144  RCL 22 145  CHS 146  RCL 21 147  CHS 148  XEQ "Q^R"  149  STO 01 150  RDN 151  STO02 152  RDN 153  STO 03 154  X<>Y 155  STO 04 156  PI 157  SQRT 158  2 159  RCL 09 160  Y^X 161  * 162  ST* 01 163  ST* 02 164  ST* 03 165  ST* 04 166  XEQ"Q*Q" 167  END

( 295 bytes / SIZE 034 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Qnm(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                  with       R09 = m  ,  R10 = n

Examples:       m = 0.4    n = 1.3       q = 0.1 + 0.2 i + 0.3 j + 0.4 k

0.4   STO 09
1.3   STO 10

0.4    ENTER^
0.3    ENTER^
0.2    ENTER^
0.1    XEQ "ALF2Q2"  >>>>   -0.944581517               ---Execution time = 2mn28s---
RDN   -0.368022565
RDN   -0.552033847
RDN   -0.736045129

Q1.30.4( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = -0.944581517 - 0.368022565 i - 0.552033847 j - 0.736045129 k

-Likewise, "ALF2Q2"  returns in 86 seconds:

Q43( 0.1 + 0.2 i + 0.3 j + 0.4 k ) = 131.9022422 - 13.96653106 i - 20.94979658 j - 27.93306214 k

b)  Program#2 - Function of Type 3 - | q | > 1

Formulae:
~
Qnm(q) =  exp( i m.PI ) 2 -n-1 sqrt(PI) Gam(m+n+1) q -m-n-1 (q2-1)m/2  2F1( (2+m+n)/2 , (1+m+n)/2 ; n+3/2 ; 1/q2 )

-Left-multiplying by exp( i m.PI ) produces a different result from right-multiplying by exp( i m.PI )
-In the following program, the lfte-multiplication has been choosen, but it is only conventional.

Data Registers:          •  R09 = m          •  R10 = n        ( Registers R09 &  R10 are to be initialized before executing "ALF2Q3" )

Flags:     F00  is used by "HGFQ"

Subroutines:

"HGFQ"
"Q*Q"  "Q^R"  "1/Q"  ( cf "Quaternions for the HP-41" )
"1/G+"  or  "GAM+"   ( cf  "Gamma Function" )

 01  LBL "ALF2Q3"   02  STO 16   03  RDN   04  STO 17   05  RDN   06  STO 18   07  X<>Y   08  STO 19   09  RCL 09   10  STO 20   11  RCL 10   12  STO 21   13  +   14  2   15  ST+ Y   16  /   17  STO 09   18  LASTX   19  1/X   20  -   21  STO 10   22  LASTX   23  1   24  +   25  RCL 21   26  +   27  STO 11   28  2 29  STO 00   30  RCL 19   31  RCL 18              32  RCL 17   33  RCL 16   34  XEQ "Q^R"   35  STO 22   36  RDN   37  STO 23   38  RDN   39  STO 24   40  RDN   41  STO 25   42  RDN   43  XEQ "1/Q"   44  XEQ "HGFQ"    45  FS?C 00   46  GTO 00   47  RCL 11   48  XEQ "1/G+"   49  ST* 12   50  ST* 13   51  ST* 14   52  ST* 15   53  LBL 00   54  RCL 20   55  2   56  / 57  STO 00   58  1   59  ST- 22   60  RCL 25              61  RCL 24   62  RCL 23   63  RCL 22   64  XEQ "Q^R"   65  STO 01   66  RDN   67  STO 02   68  RDN   69  STO 03   70  X<>Y   71  STO 04   72  RCL 12   73  STO 05   74  RCL 13   75  STO 06   76  RCL 14   77  STO 07   78  RCL 15   79  STO 08   80  XEQ "Q*Q"    81  STO 05   82  RDN   83  STO 06   84  RDN 85  STO 07   86  X<>Y   87  STO 08   88  RCL 10              89  ST+ X   90  CHS   91  STO 00   92  RCL 19   93  RCL 18   94  RCL 17   95  RCL 16   96  XEQ "Q^R"   97  STO 01   98  RDN   99  STO 02 100  RDN 101  STO 03 102  X<>Y 103  STO 04 104  XEQ "Q*Q"  105  STO 05 106  RDN 107  STO 06 108  RDN 109  STO 07 110  X<>Y 111  STO 08 112  RCL 00 113  CHS 114  XEQ "1/G+" 115  PI 116  SQRT 117  X<>Y 118  / 119  2 120  ST/ Y 121  RCL 21            122  STO 10 123  Y^X 124  / 125  1 126  CHS 127  ACOS 128  RCL 20 129  STO 09 130  * 131  X<>Y 132  P-R 133  STO 01 134  X<>Y 135  STO 02 136  CLX 137  STO 03 138  STO 04 139  XEQ "Q*Q"  140  END

( 229 bytes / SIZE 026 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Qnm(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k                  with       R09 = m  ,  R10 = n

Example:     m = sqrt(2)    n = - sqrt(3)       q = 1 + 2 i + 3 j + 4 k

2    SQRT   STO 09
3    SQRT   STO 10

4    ENTER^
3    ENTER^
2    ENTER^
1    XEQ "ALF2Q3"  >>>>   -0.478904370             ---Execuyion time = 61s---
RDN     1.598255777
RDN   -1.355956855
RDN     1.755128357

Qsqrt(2)-sqrt(3) ( 1 + 2 i + 3 j + 4 k ) = -0.478904370 + 1.598255777 i - 1.355956855 j + 1.755128357 k

14°)  Anger & Weber Functions

Formulae:

Jn(q) = + (q/2) sin ( 90°n ) 1F2( 1 ; (3-n)/2 , (3+n)/2 ; -q2/4 ) / Gam((3-n)/2) / Gam ((3+n)/2)            Anger's functions
+ cos ( 90°n ) 1F2( 1 ; (2-n)/2 , (2+n)/2 ; -q2/4 ) / Gam((2-n)/2) / Gam ((2+n)/2)

En(q) = - (q/2) cos ( 90°n )  1F2( 1 ; (3-n)/2 , (3+n)/2 ; -q2/4 ) / Gam((3-n)/2) / Gam ((3+n)/2)          Weber's functions
+ sin ( 90°n ) 1F2( 1 ; (2-n)/2 , (2+n)/2 ; -q2/4 ) / Gam((2-n)/2) / Gam ((2+n)/2)

Data Registers:         •  R00 = n                                           ( Registers R00 is to be initialized before executing "ANWEBQ" )

When the program stops,

R01 to R04 = Jn(q)         Anger Function
R05 to R08 = En(q)        Weber Function

Flag:  F04
Subroutines:

"Q*Q"  "Q^R"  "1/Q"  ( cf "Quaternions for the HP-41" )
"1/G+"  or  "GAM+"   ( cf  "Gamma Function" )
"HGFQ" ( cf §11 )  modified as follows:

>>>  Replace lines 114-115 ( /  * )  by  FC? 04  /   FS? 04  *  FC? 04  *  FS? 04  /
and  add  FS? 04   GTO 02  after line 16

-Thus,  "HGFQ" will calculate  1F2  if  SF 04  and   2F1  if  CF 04

-Lines 58-59 may be replaced by   RCL 21   STO 25   RCL 22   STO 26   RCL 23   STO 27   RCL 24   STO 28

 01  LBL "ANWEBQ"   02  STO 01                03  STO 05   04  STO 16   05  RDN   06  STO 02   07  STO 06   08  STO 17   09  RDN   10  STO 03   11  STO 07   12  STO 18   13  X<>Y   14  STO 04   15  STO 08   16  STO 19   17  4   18  CHS   19  ST/ 01   20  ST/ 02   21  ST/ 03   22  ST/ 04   23  2   24  RCL 00   25  STO 20   26  -   27  2   28  /   29  STO 10   30  RCL 00   31  LASTX   32  ST+ Y   33  /   34  STO 11 35  1   36  STO 09                37  XEQ "Q*Q"   38  SF 04   39  XEQ "HGFQ"    40  STO 21   41  RDN   42  STO 22   43  RDN   44  STO 23   45  X<>Y   46  STO 24   47  RCL 10   48  XEQ "1/G+"   49  STO 00   50  RCL 11   51  XEQ "1/G+"   52  RCL 00   53  *   54  ST* 21   55  ST* 22   56  ST* 23   57  ST* 24   58  21.025004   59  REGMOVE   60  2   61  1/X   62  ST+ 10   63  ST+ 11   64  RCL 04   65  RCL 03   66  RCL 02   67  RCL 01   68  XEQ "HGFQ" 69  RCL 10   70  XEQ "1/G+"   71  STO 00                72  RCL 11   73  XEQ "1/G+"    74  RCL 00   75  *   76  ST* 12   77  ST* 13   78  ST* 14   79  ST* 15   80  RCL 12   81  STO 05   82  RCL 13   83  STO 06   84  RCL 14   85  STO 07   86  RCL 15   87  STO 08   88  RCL 16   89  STO 01   90  RCL 17   91  STO 02   92  RCL 18   93  STO 03   94  RCL 19   95  STO 04   96  XEQ "Q*Q"   97  STO 12   98  STO 16   99  RDN 100  STO 13 101  STO 17 102  RDN 103  STO 14 104  STO 18              105  X<>Y 106  STO 15 107  STO 19  108  RCL 20 109  STO 00 110  1 111  ASIN 112  * 113  1 114  P-R 115  ST* 21 116  ST* 22 117  ST* 23 118  ST* 24 119  2 120  / 121  ST* 16 122  ST* 17 123  ST* 18 124  ST* 19 125  X<>Y 126  ST* 25 127  ST* 26 128  ST* 27 129  ST* 28 130  LASTX 131  / 132  ST* 12 133  ST* 13 134  ST* 14 135  ST* 15 136  RCL 25 137  RCL 16 138  - 139  STO 05              140  RCL 26 141  RCL 17  142  - 143  STO 06 144  RCL 27 145  RCL 18 146  - 147  STO 07 148  RCL 28 149  RCL 19 150  - 151  STO 08 152  RCL 12 153  RCL 21 154  + 155  STO 01 156  RCL 15 157  RCL 24 158  + 159  STO 04 160  RCL 14 161  RCL 23 162  + 163  STO 03 164  RCL 13 165  RCL 22 166  + 167  STO 02 168  RCL 01 169  CF 04 170  END

( 291 bytes / SIZE 029 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

where      Jn(x+y.i+z.j+t.k)  = x'+y'.i+z'.j+t'.k   = R01 thru R04            with       R00 = n  # +2 , -2 ,  +3 , -3 ,..........
and       En(x+y.i+z.j+t.k)  = x"+y".i+z".j+t".k  = R05 thru R08

Example:         n = PI     q = 1 + 2 i + 3 j + 4 k

PI  STO 00

4   ENTER^
3   ENTER^
2   ENTER^
1   XEQ "ANWEBQ"   >>>>  -11.24234896 = R01                         ---Execution time = 123s---
RDN    -3.564968434 = R02
RDN    -5.347452650 = R03
RDN    -7.129936871 = R04

JPi( 1 + 2 i + 3 j + 4 k ) = -11.24234896 - 3.564968434 i - 5.347452650 j - 7.129936871 k

and      EPi( 1 + 2 i + 3 j + 4 k ) = -9.566817362 + 4.177204891 i + 6.265807334 j + 8.354409784 k             in  R05 to R08

Notes:

-If n is an integer,  Jn(x) = Jn(x) = Bessel function of the first kind.
-However, this routine does not work if n is an integer, except if n = -1 , 0 , +1

-But it will work if we use the regularized hypergeometric function 1F2 tilde.
( cf "Quaternionic Special Functions (II) )

15°)  Regular Coulomb Wave Function

-This program employs the formulae given in "Coulomb Wave Functions for the HP-41" , paragraph 1°) a) & b)

Data Registers:           •  R09 = L          •  R10 = n            ( Registers R09 & R10 are to be initialized before executing "RCWFQ" )

R00 to R08 & R11 to R17: temp
R14 thru R17 contain the partial sums of the series
Flag:  F10 is cleared
Subroutines:  "Q*Q" & "Q^R"  ( cf "Quaternions for the HP-41" )

 01  LBL "RCWFQ"   02  STO 05               03  RDN   04  STO 06   05  RDN   06  STO 07   07  X<>Y   08  STO 08   09  CLX   10  STO 02   11  STO 03   12  STO 04   13  STO 11   14  STO 15   15  STO 16   16  STO 17   17  SIGN   18  STO 01   19  STO 12   20  STO 13   21  STO 14   22  RCL 09   23  STO 00   24  X=0?   25  GTO 00   26  SIGN   27  LBL 01   28  RCL 00   29  X^2   30  RCL 10   31  X^2   32  +   33  *   34  DSE 00   35  GTO 01 36  LBL 00   37  X=0?   38  SIGN   39  RCL 10               40  PI   41  *   42  STO 00   43  ENTER^   44  X=0?   45  ISG Y   46  LBL 00   47  E^X-1   48  LASTX   49  CHS   50  E^X-1   51  -   52  2   53  /   54  X=0?   55  SIGN   56  /   57  *   58  SQRT   59  2   60  ST/ 00   61  RCL 09   62  Y^X   63  *   64  RCL 00   65  E^X   66  /   67  RCL 09   68  ST+ X   69  1   70  + 71  FACT   72  /   73  STO 00               74  LBL 02   75  SF 10   76  LBL 03   77  RCL 10   78  ST+ X   79  RCL 12   80  ST* Y   81  X<> 11   82  -   83  RCL 13   84  ST/ Y   85  RCL 09   86  ST+ X   87  +   88  1   89  ST+ 13   90  +   91  /   92  STO 12   93  XEQ "Q*Q"   94  STO 01   95  RDN   96  STO 02   97  RCL 12   98  *   99  RCL 15 100  + 101  STO 15 102  LASTX 103  - 104  ABS 105  X<>Y 106  STO 03 107  RCL 12             108  * 109  RCL 16 110  + 111  STO 16 112  LASTX 113  - 114  ABS 115  + 116  X<>Y 117  STO 04 118  RCL 12 119  * 120  RCL 17 121  + 122  STO 17 123  LASTX 124  - 125  ABS 126  + 127  RCL 01 128  RCL 12 129  * 130  RCL 14 131  + 132  STO 14 133  LASTX 134  - 135  ABS 136  + 137  X#0? 138  GTO 02 139  FS?C 10 140  GTO 03 141  RCL 00             142  ST* 14 143  ST* 15 144  ST* 16 145  ST* 17 146  RCL 09 147  1 148  + 149  STO 00 150  RCL 08 151  RCL 07 152  RCL 06 153  RCL 05 154  XEQ "Q^R" 155  STO 01 156  RDN 157  STO 02 158  RDN 159  STO 03 160  X<>Y 161  STO 04 162  RCL 14 163  STO 05 164  RCL 15 165  STO 06 166  RCL 16 167  STO 07 168  RCL 17 169  STO 08 170  XEQ "Q*Q" 171  END

( 221 bytes / SIZE 018 )

 STACK INPUTS OUTPUTS T t t' Z z z' Y y y' X x x'

with      FL( h ; x+y.i+z.j+t.k )  = x'+y'.i+z'.j+t'.k
where   L is a non-negative integer and  h is a real number.

Example:         L = 2  ,   h = 0.75  ;   q = 1 + 1.2 i + 1.4 j + 1.6 k

2     STO 09
0.75   STO 10

1.6     ENTER^
1.4     ENTER^
1.2     ENTER^
1      XEQ "RCWFQ"   >>>>    -0.478767859                                  ---Execution time = 53s---
RDN    -0.179694772
RDN    -0.209643900
RDN    -0.239593029

-So,   F2( 0.75 ;  1 + 1.2 i + 1.4 j + 1.6 k ) = -0.478767859 - 0.179694772 i - 0.209643900 j - 0.239593029 k

Note:

-The execution time has been measured with an M-code routine Q*Q
-The focal program "Q*Q" is of course slower.