# hp41programs

Lagrange Lagrange Polynomials for the HP-41

Overview

1°) One Dimension Problem

a) Univariate Interpolation
b) Derivative of the Lagrange Polynomial
c) Integration of the Lagrange Polynomial
d) 7 Points with Equidistant Abscissas: f(x) = 0 & f'(x) = 0

2°) Bivariate Interpolation

a) non-synthetic program
b) synthetic program

Latest Update:   §1°)d)

1°) One Dimension Problem

a) Univariate Interpolation

-Let be a set of n data points:   M1 ( x1 ; y1 )  ,   M2 ( x2 ; y2 )  ,  ..........  ,  Mn ( xn ; yn )
There is a unique polynomial L of degree < n  such that:     L(xi) = yi       i = 1 ; 2 ; ..... ; n
-In other words, L is the collocation polynomial for these arguments.
( xi values can be unequally-spaced.)

-The following program calculates  L(x) = y1 L1(x) + y2 L2(x) + ..... + yn Ln(x)

where    Li(x) = ( x - x1 ) ...... ( x - xi-1 ) ( x - xi+1 ) ...... ( x - xn ) / [ ( xi - x1 ) ...... ( xi - xi-1 ) ( xi - xi+1 ) ( xi - xn ) ]
are the Lagrange's multiplier functions.

-This formula is also applied to an extrapolation to the limit in example2.

Data Registers:

-You have to store the 2n numbers   x1 , y1 , x2 , y2 , ........ , xn , yn  into contiguous registers,  starting at any register  Rbb:

 x1 x2 ........... xn y1 y2 ........... yn

into

 Rbb Rbb+2 ............. Ree-1 Rbb+1 Rbb+3 ............. Ree

( ee = bb + 2n -1 )

Flags:  /
Subroutines:  /

-Thanks to synthetic programming,  "LAGR" only uses the registers containing the  xi and yi numbers,
but you can replace status registers M N O P Q by the standard R00 R01 R02 R03 R04 ( in this case, add  CLX  STO 03  after line 06 and of course, choose bb > 04 ).

 01  LBL "LAGR" 02  CLA 03  STO M          04  X<>Y 05  STO O 06  STO Q 07  LBL 01 08  RCL Q 09  STO N 10  SIGN 11  LBL 02 12  RCL M          13  RCL IND O 14  RCL IND N 15  ST- Z 16  - 17  X#0? 18  GTO 03 19  SIGN 20  STO Y 21  LBL 03          22  / 23  * 24  ISG N 25  ISG N           26  GTO 02 27  ISG O 28  RCL IND O 29  * 30  ST+ P 31  ISG O 32  GTO 01 33  RCL Q  34  RCL M          35  SIGN 36  X<> P 37  CLA 38  END

( 69 bytes / SIZE eee+1 )

 STACK INPUTS OUTPUTS Y bbb.eee bbb.eee X x L(x) L / x

Example1:     Evaluate  y(3) and y(5) for the collocation polynomial that takes the values prescribed below:

 xi 0 1 2 4 7 yi 3 2 4 6 5

-For instance, you store these 10 numbers ( 0 3 1 2 2 4 4 6 7 5 ) into R01 thru R10, then

1.010  ENTER^
3     XEQ "LAGR"  produces:  y(3) = 5.847619048    ( in 14 seconds )

RDN    5   R/S     >>>>     y(5) = 4.523809520

Example2:

-Lagrange's polynomial can also be used for an extrapolation to the limit:
Suppose we are using the trapezoidal rule to evaluate the integral:       I  =  §02  ex dx
The interval  [0;2]  is divided into n subintervals and

n = 1  yields   I1 = 8.389056101
n = 2  yields   I2 = 6.912809880
n = 4  yields   I4 = 6.521610110
n = 8  yields   I8 = 6.422297820

We would like to know the result with n = infinity !  i-e  1/n = 0.
The error of the trapezoidal rule is nearly proportional to  1/n2  , therefore, we can use Lagrange interpolation formula with

x1 = 1   x2 = 1/4   x3 = 1/16   x4 = 1/64
y1 = I1   y2 = I2      y3 = I4       y4 =  I8      and  evaluate   L (0)

If these 8 numbers are stored in R11 thru R18 for instance,

11.018  ENTER^
0       XEQ "LAGR"  produces   6.389056387

Exact result is of course e2 - 1 = 6.389056099 to ten places
and the error of our Lagrange approximation is only  3. 10-7  although I8 was only exact to 1 decimal!

N.B:

- If you are using Simpson's rule instead of the trapezoidal rule, store  xn = 1/n4  : it's a 4th order method ( at least far from any singularity )
-  n is usually doubled at every step, but it's not necessary here.
-  The same idea can be applied to any other problem of this kind.

b) Derivative of the Lagrange Polynomial

-"DLAGR" takes x in X-register and returns  L'(x)  in X-register.

Data Registers:

-The same 2n numbers   x1 , y1 , x2 , y2 , ........ , xn , yn  are to be stored into contiguous registers,  starting at any register  Rbb:

 x1 x2 ........... xn y1 y2 ........... yn

into

 Rbb Rbb+2 ............. Ree-1 Rbb+1 Rbb+3 ............. Ree

( ee = bb + 2n -1 )

Flags:  /
Subroutines:  /

 01  LBL "DLAGR"  02  CLA  03  STO M             04  X<>Y  05  STO Q  06  STO a  07  LBL 01  08  RCL a  09  STO O  10  CLX  11  LBL 02  12  RCL O  13  RCL Q  14  X=Y?  15  GTO 02 16  X<> Z  17  RCL a  18  STO N             19  SIGN  20  LBL 03  21  RCL N  22  RCL O  23  X=Y?  24  GTO 03  25  CLX  26  RCL Q  27  X=Y?  28  GTO 03  29  RDN  30  CLX 31  RCL M             32  RCL IND N  33  -  34  *  35  ENTER^  36  ENTER^  37  LBL 03  38  R^  39  R^  40  ISG N  41  ISG N  42  GTO 03  43  +  44  STO Z  45  LBL 02 46  X<> Z  47  ISG O  48  ISG O  49  GTO 02            50  RCL a  51  X<>Y  52  LBL 04  53  RCL IND Q  54  RCL IND Z  55  -  56  X=0?  57  SIGN  58  /  59  ISG Y  60  ISG Y 61  GTO 04  62  ISG Q  63  RCL IND Q  64  *  65  ST+ P  66  ISG Q  67  GTO 01            68  RCL a  69  RCL M  70  SIGN  71  X<> P  72  CLA  73  END

( 123 bytes / SIZE eee+1 )

 STACK INPUTS OUTPUTS Y bbb.eee bbb.eee X x L' (x) L / x

Example:

Evaluate  y(PI)  for the collocation polynomial that takes the values prescribed below:

 xi 0 1 2 4 7 yi 3 2 4 6 5

-If, for instance, you store these 10 numbers ( 0 3 1 2 2 4 4 6 7 5 ) into R01 thru R10, then

1.010  ENTER^
PI    XEQ "DLAGR"  >>>>  L' (PI) = 0.857353871                                  ---Execution time = 67s---

Notes:

-Synthetic register a is used, so this program cannot be called as more than a first level subroutine.
-Registers M N O P Q a  may of course be replaced by standard registers like R00 thru R05
-In this case, choose  bbb > 005.

c) Integration of the Lagrange Polynomial

-This program is also listed in "Numerical Integration for the HP-41"

n data points are given and stored in registers  R01 thru R2n

 x1 x2 ........... xn y1 y2 ........... yn

-The program below calculates its coefficients c0 , c1 , ........... , cn-1  when L(x) is written:

L(x) = c0 + c1 ( x-x1 ) + c2 ( x-x1 )2 + .................. + cn-1 ( x-x1 )n-1

-First, x1 is subtracted from all other xi 's
-We have obviously  c0 = y1  and  the other yi 's  are replaced by  ( yi - y1 ) / ( xi - x1 )
-Then, we find the value of the new Lagrange polynomial at zero, thus producing c1 which is stored in register R04
-The process is repeated until  cn-1 is computed and stored in register R2n
-Finally, lines 42 to 57 evaluate  §x1xn  L(x).dx

Data Registers:       •  R00 = n = number of points                                ( Registers R00 thru R2n  are to be initialized before executing "LGI" )

•  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
•  R02 = y1     •  R04 = y2      .......   •  R2n   = yn

>>>>    when the program stops,  R02 = c0 = y1 ,  R04 = c1 , ............... ,  R2n = cn-1

Flags:  /
Subroutine:   "LAGR"  ( cf paragrph 1-a) )

 01  LBL "LGI" 02  RCL 00 03  ST+ X 04  STO Y 05  RCL 01          06  LBL 00 07  DSE Z 08  ST- IND Z  09  DSE Z 10  GTO 00 11  CLX 12   E3 13  / 14  ISG X 15  STO 01 16  RCL 02 17  GTO 02          18  LBL 01 19  2 20  ST+ 01 21  RCL 01 22  0 23  XEQ "LAGR" 24  LBL 02 25  ISG Y 26  STO IND Y 27  ISG Y 28  FS? 30 29  GTO 04          30  LBL 03 31  RCL IND Y  32  ISG Z 33  X<>Y 34  ST- IND Z 35  X<>Y 36  ST/ IND Z 37  RDN 38  ISG Y 39  GTO 03          40  GTO 01 41  LBL 04 42  RCL 00 43  ST+ X 44  0 45  LBL 05 46  RCL IND Y  47  RCL 00 48  / 49  + 50  RCL IND 01  51  * 52  2 53  ST- Z 54  RDN 55  DSE 00 56  GTO 05          57  END

( 98 bytes / SIZE 2n+1 )

 STACK INPUT OUTPUT X / §x1xn  L(x).dx

Example:     Calculate  §18  L(x).dx   from the following values

 x 1 2.4 4 5.2 7 8 y 1 4 6 5 4 2

-Store these 12 numbers into registers

R01  R03  R05  R07  R09  R11   ( x )
R02  R04  R06  R08  R10  R12    ( y )    respectively

-There are 6 points:     6   STO 00      XEQ "LGI"     >>>>    §18  L(x).dx  =   29.61789480  ( in 43 seconds )

-We have also the coefficients of the Lagrange polynomial ( as a sum of powers of ( x-x1 ) )  in registers  R02  R04  .............  R12
-Here, x1 = 1 and we find:

L(x) =  1 - 0.362103178 ( x-1 ) + 3.623795356 ( x-1 )2 - 1.661873944 ( x-1 )3 + 0.272598127 ( x-1 )4 - 0.015381483 ( x-1 )5

Remarks:

-For large n-values, the coefficients of the Lagrange polynomial are often great in magnitude
and of alternate signs, thus leading to inaccurate results.
-Moreover, the coefficients themselves are not obtained very accurately if n is too great.

-For instance, with the following data:

 x 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 y 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40

"LGI"  yields  797.9971774  wheras the exact result is 798

Execution time t:

 n t 4 15s 6 43s 10 2mn51s 20 21mn13s

d) 7 Points with Equidistant Abscissas:  f(x) = 0 & f'(x) = 0

-The following program uses the Lagrange 7-point interpolation formula.
-We assume that the xi's are equally spaced. ( xi+1 - xi = h  ,   i  = -3 , -2 , -1 , 0 , 1 , 2 )

-So,   f(x0+p.h) - f(x0) = A.p + B.p2 + C.p3 + D.p4 + E.p5 + F.p6   with:

A = (1/720) [ 540 ( y1 - y-1 ) - 108 ( y2 - y-2 )  + 12 ( y3 - y-3 ) ]
B = (1/720) [ -980 y0 + 540 ( y1 + y-1 ) - 54 ( y2 + y-2 )  + 4 ( y3 + y-3 ) ]
C = (1/720) [ -195 ( y1 - y-1 ) + 120 ( y2 - y-2 )  - 15 ( y3 - y-3 ) ]
D = (1/720) [ 280 y0 - 195 ( y1 + y-1 ) + 60 ( y2 + y-2 )  - 5 ( y3 + y-3 ) ]                         yk = f(xk)
E = (1/720) [ 15 ( y1 - y-1 ) - 12 ( y2 - y-2 )  + 3 ( y3 - y-3 ) ]
F = (1/720) [ -20 y0 + 15 ( y1 + y-1 ) - 6 ( y2 + y-2 )  + ( y3 + y-3 ) ]

Data Registers:              R00 = x0            R16 = h             ( Registers R01 thru R07 are to be initialized before executing "SLVD7" )

•  R01 = f(x-3)    •  R02 = f(x-2)    •  R03 = f(x-1)    •  R04 = f(x0)    •  R05 = f(x1)    •  R04 = f(x2)    •  R05 = f(x3)                R08 to R15:  temp

Flag:   F01                       CF 01  to solve f(x) = 0
SF 01  to solve f'(x) = 0   ( extrema )
Subroutines:  /

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

( 232 bytes / SIZE 017 )

 STACK INPUTS OUTPUTS Y h f(x) X x0 x

If CF 01   x is a solution of f(x) = 0
If SF 01   x is a solution of f'(x) = 0   ( extremum )

Example1:     f(x) = 0   from these 7 values

 x 1 3 5 7 9 11 13 y 16 10 6 1 -3 -8 -20

16  STO 01   10  STO 02   6  STO 03   1  STO 04  -3  STO 05  -8  STO 06  -20  STO 07

CF 01           h = 2  &  x0 = 7   so:

2  ENTER^
7  XEQ "SLVD7"   >>>>  x = 7.446084649                     ---Execution time = 20s---

Example2:     Extremum from these 7 values

 x 1 3 5 7 9 11 13 y 16 10 6 4 5 8 20

16  STO 01   10  STO 02   6  STO 03   4  STO 04   5  STO 05   8  STO 06   20  STO 07

SF 01           h = 2  &  x0 = 7   so:

2  ENTER^
7  XEQ "SLVD7"   >>>>      x  =  7.235889441                     ---Execution time = 16s---
X<>Y  f(x) =  3.977546065

-So  f(x)  has a minimum  3.977546065  for x = 7.235889441

Notes:

-The method may fail , for example if  A = 0  when CF 01

-It's also possible to have an infinite loop with  X#Y?  ( line 159 )
-It would be better to replace line 159 by the M-Code routine  X#Y??  ( cf "A few M-Code routines for the HP41" )

-Similar programs with 5-point Lagrange polynomials are listed in "Linear & Non-Linear systems for the HP41" ( paragraph 2°)a4) )
and "Extrema for the HP41" ( paragraph 1°)e) )

2°) Bivariate Interpolation

a) Non-synthetic program

-Now we have a grid of  nxm  values:   f(xi,yj)   with   i = 1,2,......,n  &  j = 1,2,.....,m   and we want to estimate f(x,y)

-"LAGR2" uses the Lagrange Polynomial in the directions of x-axis and y-axis.
-You have to store the data in registers Rbb to Ree as shown below ( bb > 10 ):

Data Registers:           •  R00 = bbb.eeenn                             ( Registers R00 and Rbb thru Ree are to be initialized before executing "LAGR2"  )

R01 = x    R03 = f(x,y)
R02 = y    R04 thru R10: temp                                             We must choose bb > 10

x  \  y |        •  Rbb+n = y1                •  Rbb+2n+1 = y2       ..............................    •  Ree-n = ym
---------------------------------------------------------------------------------------------------------------------
•  Rbb         = x1 |   •  Rbb+n+1 = f(x1,y1)    •  Rbb+2n+2 = f(x1,y2)    ..............................     •  Ree-n+1 = f(x1,ym)
•  Rbb+1     = x2 |   •  Rbb+n+2 = f(x2,y1)    •  Rbb+2n+3 = f(x2,y2)    ..............................     •  Ree-n+2 = f(x2,ym)
....................    |     .............................................................................................................................................
|
•  Rbb+n-1 = xn  |   •  Rbb+2n  = f(xn,y1)     •  Rbb+3n+1 = f(xn,y2)   ...............................     •  Ree  = f(xn,ym)

Flags: /
Subroutines: /

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

( 121 bytes / SIZE? )

 STACK INPUTS OUTPUTS Y y y X x f(x,y) L / x

Example:      A function f(x,y) is only known by the following values ( meaning for example:  f(1,4) = 3   f(4,6) = 9  ...etc... )

>>>   if, for instance, you choose bb = 11 ,  store these  19 numbers  into

x \ y |   2    3    4    6                                                          |   R14   R18   R22   R26
--------------------------                                              ---------------------------------
1   |   4    3    3    5                                                  R11  |  R15   R19   R23   R27                             respectively
2   |   3    1    2    6                                                  R12  |  R16   R20   R24   R28
4   |   1    0    4    9                                                  R13  |  R17   R21   R25   R29

-The control number of this table is   11.02903  so   11.02903  STO 00
( note that nn = 03 is the number of x-values which may be different from m = the number of y-values ( here m = 4 ) )

-Then, to compute, say  f(3,5)

5  ENTER^
3  XEQ "LAGR2"  >>>>   f(3,5) ~  5.8333   ( in 29 seconds )

-Likewise, to obtain  f(1.6,2.7)

2.7  ENTER^
1.6     R/S        >>>>   f(1.6,2.7) ~  1.9038

Note:   Execution time is approximately proportional to nxm

b) Synthetic program

-Synthetic programming allows to store the data into registers R01 thru Rn.m+n+m
-However, synthetic register a is used, so this program cannot be called as more than a first level subroutine.

Data Registers:           •  R00 = n.mmm = n + m/103         ( Registers R00 thru Rnm+m+n are to be initialized before executing "BVI"  )

x  \  y |       •  Rnn+1 = y1             •  R2n+2 = y2       ..............................    •  Rnm+m = ym
---------------------------------------------------------------------------------------------------------------------
•  R01  = x1 |   •  Rnn+2 = f(x1,y1)    •  R2n+3 = f(x1,y2)    ..............................     •  Rnm+m+1 = f(x1,ym)
•  R02  = x2 |   •  Rnn+3 = f(x2,y1)    •  R2n+4 = f(x2,y2)    ..............................     •  Rnm+m+2 = f(x2,ym)
...........    |      .............................................................................................................................................
|
•  Rnn = xn  |   •  R2n+1  = f(xn,y1)   •  R3n+2 = f(xn,y2)   ..............................      •  Rnm+m+n  = f(xn,ym)

Flags: /
Subroutines: /

 01  LBL "BVI"  02  CLA 03  STO M           04  X<>Y 05  STO N 06  RCL 00 07  INT 08  STO P 09  LBL 01 10  RCL 00 11  INT 12  LASTX 13  FRC 14  ENTER^ 15  SIGN 16  10^X 17  CHS 18  10^X 19  SQRT 20  + 21  ST* Y 22  + 23  ISG X 24  STO Q           25  CLX 26  LBL 02 27  RCL Q 28  FRC 29  ISG X 30  STO a 31  SIGN 32  LBL 03 33  RCL N 34  RCL IND a 35  - 36  RCL IND Q 37  ST- L 38  X<> L 39  CHS 40  X#0? 41  GTO 03 42  SIGN 43  STO Y 44  LBL 03 45  / 46  * 47  ISG a 48  GTO 03 49  RCL P           50  RCL Q 51  + 52  RDN 53  RCL IND T 54  * 55  + 56  ISG Q           57  GTO 02 58  RCL 00 59  INT 60  STO Q 61  X<>Y 62  LBL 04 63  RCL M 64  RCL IND P 65  RCL IND Q 66  ST- Z 67  - 68  X#0? 69  GTO 04 70  SIGN 71  STO Y 72  LBL 04 73  / 74  * 75  DSE Q           76  GTO 04 77  ST+ O 78  DSE P 79  GTO 01 80  RCL M 81  SIGN 82  RCL N 83  RCL O 84  CLA 85  END

( 131 bytes / SIZE (n+1)(m+1) )

 STACK INPUTS OUTPUTS Y y y X x f(x,y) L / x

Example:      With the same data as above:

x \ y |   2    3    4    6                                          |   R04   R08   R12   R16
--------------------------                              --------------------------------
1   |   4    3    3    5                =                R01  |  R05   R09   R13   R17            respectively
2   |   3    1    2    6                                  R02  |  R06   R10   R14   R18
4   |   1    0    4    9                                  R03  |  R07   R11   R15   R19

-Here,  n = 3 and  m = 4  so   3.004  STO 00

-Then, to compute  f(3,5)

5  ENTER^
3  XEQ "BVI"  >>>>   f(3,5) ~  5.8333   ( in 31 seconds )

-And to obtain  f(1.6,2.7)

2.7  ENTER^
1.6     R/S        >>>>   f(1.6,2.7) ~  1.9038

References:

[1]   Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]   Press, Vetterling, Teukolsky, Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4