The Museum of HP Calculators

# Numerical Integration for the HP-41

This program is Copyright © 2005 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) Continuous Data

a)  Gauss-Legendre 3-point Formula
b)  Gauss-Legendre 16-point ( & 48-point ) Formula
c)  Gauss-Chebyshev Formula
d)  Filon's Formulae
e)  An Example of Curvilinear Integral

2°) Discrete Data

a)  Equally Spaced Abscissas

a-1) Simpson's Rule
a-2) Other Newton-Cotes Formulae

b)  Unequally Spaced Abscissas

b-1) Trapezoidal Rule  ( new )
b-2) Connecting Parabolic Segments
b-3) Connecting Cubic Segments ( new )
b-4) Natural Cubic Spline ( new )
b-5) Integration of the Lagrange Polynomial ( new )

c)  Double Integrals ( Equally Spaced Arguments only )
d)  Triple Integrals ( Equally Spaced Arguments only )

"Triple Integrals for the HP-41"
and "Multiple Integration for the HP-41"

-Integrals are denoted §      §ab f(x).dx

1°) Continuous Data

a) Gauss-Legendre 3-Point  Formula

-Gauss-Legendre formulas are defined by  §-11 f(x).dx  ~  w1.x1 + w2.x2 + ....... + wn.xn
where wi and xi  are choosen so that the formula produces perfect accuracy when  f(x) is a polynomial of degree < 2n
-The  xi  are the zeros of the Legendre polynomials and they can't be easily computed for large n-values.

-The 3-point Gauss-Legendre formula is   §-11 f(x).dx  ~  (5/9).f(-(3/5)1/2) + (8/9).f(0) + (5/9).f((3/5)1/2)
-Then, a linear change of variable transforms  [-1;1]  into any ( finite ) interval [a;b]

-Actually, the interval [a;b] is divided into  n sub-intervals  [ a+k.(b-a)/n ; a+(k+1).(b-a)/n ]    k = 0 , 1 , ...... , n-1
and, at least theoretically, the results tend to the exact integral as n tends to infinity
( provided f and its derivative have no "strong" singularity ).

Data Registers:       •  R00 = Function name                             ( Registers R00 thru R03 are to be initialized before executing "GL3" )
•  R01 = a
•  R02 = b
•  R03 = n = number of subintervals over which the formula is applied.          R04 =  §ab f(x).dx        R05 thru R08: temp
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.

01  LBL "GL3"
02  RCL 02
03  RCL 01
04  STO 07
05  -
06  RCL 03
07  STO 05
08  /
09  STO 06
10  2
11  /
12  ST+ 07
13  .6
14  SQRT
15  *
16  STO 08
17  CLX
18  STO 04
19  LBL 01
20  RCL 07
21  RCL 08
22  -
23  XEQ IND 00
24  ST+ 04
25  RCL 07
26  XEQ IND 00
27  1.6
28  *
29  ST+ 04
30  RCL 07
31  RCL 08
32  +
33  XEQ IND 00
34  ST+ 04
35  RCL 06
36  ST+ 07
37  DSE 05
38  GTO 01
39  RCL 04
40  *
41  3.6
42  /
43  STO 04
44  END

( 67 bytes / SIZE 009 )

Example:     Evaluate  §13  e-x^2 dx

01  LBL "T"   ( any global Label , maximum 6 characters )
02  X^2
03  CHS
04  E^X
05  RTN

"T"  ASTO 00
1    STO 01
3    STO 02       if  n = 2    2  STO 03  XEQ "GL3"  >>>>  0.139390854
n = 4    4  STO 03       R/S          >>>>  0.139383255
n = 8    8  STO 03       R/S          >>>>  0.139383216      ( execution time = 17 seconds )

-If a or b is infinite, make a change of variable ( like x = Tan u ) to reduce [ a ; b ] to a finite interval.

b) Gauss-Legendre 16-Point Formula

-The 16 coefficients wi and xi  are to be stored in registers R11 to R26 as listed below.

Data Registers:       •  R00 = Function name            ( Registers R00 thru R03 and R11 thru R26  are to be initialized before executing "GL16" )
•  R01 = a
•  R02 = b
•  R03 = n = number of subintervals over which the formula is applied.    R04 = Integral         R05 thru R10: temp

•  R11 = 0.02715245941         •  R19 = 0.1495959888
•  R12 = 0.9894009350           •  R20 = 0.6178762444
•  R13 = 0.06225352394         •  R21 = 0.1691565194
•  R14 = 0.9445750231           •  R22 = 0.4580167777
•  R15 = 0.09515851168         •  R23 = 0.1826034150
•  R16 = 0.8656312024           •  R24 = 0.2816035508
•  R17 = 0.1246289713           •  R25 = 0.1894506105
•  R18 = 0.7554044084           •  R26 = 0.09501250984
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.

01  LBL "GL16"
02  CLX
03  STO 04
04  RCL 02
05  RCL 01
06  STO 07
07  -
08  RCL 03
09  STO 05
10  ST+ X
11  /
12  STO 06
13  LBL 01
14  ST+ 07
15  26.01
16  STO 08
17  LBL 02
18  RCL 07
19  RCL IND 08
20  RCL 06
21  *
22  STO 09
23  -
24  XEQ IND 00
25  STO 10
26  RCL 07
27  RCL 09
28  +
29  XEQ IND 00
30  RCL 10
31  +
32  DSE 08
33  RCL IND 08
34  *
35  ST+ 04
36  DSE 08
37  GTO 02
38  RCL 06
39  ST+ 07
40  DSE 05
41  GTO 01
42  ST* 04
43  RCL 04
44  END

( 71 bytes / SIZE 027 )

Example:     Evaluate  §03  e-x^4 dx

01  LBL "T"   ( any global Label , maximum 6 characters )
02  X^2
03  X^2
04  CHS
05  E^X
06  RTN

"T"  ASTO 00
0    STO 01
3    STO 02       if  n = 1    1  STO 03  XEQ "GL16"  >>>>  0.906402825
n = 2    2  STO 03        R/S           >>>>  0.906402476    ( in 28 seconds )

Remark:

-If you want to use the 48-point formula,
replace line 01 by  LBL "GL48"
----------- 15 by  58.01            ( instead of 26.01 )
and store these 48 coefficients in registers R11 to R58

R11 = 0.003153346052           R23 = 0.02742650971             R35 = 0.04761665849          R47 = 0.06070443917
R12 = 0.9987710073               R24 = 0.9058791367               R36 = 0.6778723796            R48 = 0.3487558863
R13 = 0.007327553901           R25 = 0.03116722783             R37 = 0.05035903555          R49 = 0.06203942316
R14 = 0.9935301723               R26 = 0.8765720203               R38 = 0.6288673968            R50 = 0.2873624874
R15 = 0.01147723458             R27 = 0.03477722256             R39 = 0.05289018949          R51 = 0.06311419229
R16 = 0.9841245837               R28 = 0.8435882616               R40 = 0.5772247261            R52 = 0.2247637904
R17 = 0.01557931572             R29 = 0.03824135107             R41 = 0.05519950370          R53 = 0.06392423858
R18 = 0.9705915925               R30 = 0.8070662040               R42 = 0.5231609747            R54 = 0.1612223561
R19 = 0.01961616046             R31 = 0.04154508294             R43 = 0.05727729210          R55 = 0.06446616444
R20 = 0.9529877032               R32 = 0.7671590325               R44 = 0.4669029048            R56 = 0.09700469921
R21 = 0.02357076084             R33 = 0.04467456086             R45 = 0.05911483970          R57 = 0.06473769681
R22 = 0.9313866907               R34 = 0.7240341309               R46 = 0.4086864820            R58 = 0.03238017096

c) Gauss-Chebyshev Formula

-This program calculates the singular integral:    §ab  ((x-a)(x-b))-1/2  f(x).dx  ~  (pi/n) [ f(x1) + f(x2) + ..... + f(xn) ]
where   xi = (a+b)/2 + ((b-a)/2).cos((2i-1)pi/(2n))
-Unlike "GL3" and "GL16", "GC" determines the coefficients wi and xi directly since they are much easier to obtain.

Data Registers:       •  R00 = Function name        ( Registers R00 thru R03 are to be initialized before executing "GC" )
•  R01 = a
•  R02 = b
•  R03 = n          R04 = Integral         R05 thru R07: temp
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.

01  LBL "GC"
02  RCL 03
03  STO 05
04  RCL 02
05  RCL 01
06  -
07  2
08  /
09  STO 06
10  CLX
11  STO 04
12  LBL 01
13  1                       CLX  SIGN  is faster
14  ASIN
15  RCL 05
16  ST+ X
17  DSE X              or  LASTX   -
18  RCL 03
19  /
20  *
21  COS
22  RCL 06
23  ST* Y
24  +
25  RCL 01
26  +
27  XEQ IND 00
28  ST+ 04
29  DSE 05
30  GTO 01
31  PI
32  RCL 03
33  /
34  ST* 04
35  RCL 04
36  END

( 51 bytes / SIZE 007 )

Example:     Compute  §13  ex(-x2+4.x-3)-1/2 dx  =  §13  ex((3-x)(x-1))-1/2 dx

01  LBL "T"   ( any global Label , maximum 6 characters )
02  E^X
03  RTN

"T"  ASTO 00
1    STO 01
3    STO 02       if  n = 2    2  STO 03  XEQ "GC"  >>>>  29.262
n = 4    4  STO 03       R/S         >>>>  29.389695
n = 8    8  STO 03       R/S         >>>>  29.38969917
n = 16  16  STO 03      R/S         >>>>  29.38969918    ( in 25 seconds )

Note:   This program works in all angular mode,
but "GC" and the subroutine must be executed in the same angular mode.

d) Filon's Formulae

-This method computes integrals of the form:    §ab  f(x).cos(k.x).dx   and   §ab  f(x).sin(k.x).dx

§ab  f(x).cos(k.x).dx  ~  h.[ A(k.h) (  f(x2n).sin(k.x2n)  -  f(x0).sin(k.x0) ) + B(k.h) C2n + C(k.h) C2n-1 ]                 a = x0 ;  b = x2n ;  h = (b-a)/(2n)
§ab  f(x).sin(k.x).dx   ~  h.[ A(k.h) ( -f(x2n).cos(k.x2n) + f(x0).cos(k.x0) ) + B(k.h) S2n + C(k.h) S2n-1 ]

where   C2n = f(x0).cos(k.x0) + f(x2).cos(k.x2) + ......... + f(x2n).cos(k.x2n)  -  (1/2) ( f(x2n).cos(k.x2n) + f(x0).cos(k.x0) )
S2n = f(x0).sin(k.x0) + f(x2).sin(k.x2) + ......... + f(x2n).sin(k.x2n)  -  (1/2) ( f(x2n).sin(k.x2n) + f(x0).sin(k.x0) )

C2n-1 = f(x1).cos(k.x1) + f(x3).cos(k.x3) + ......... + f(x2n-1).cos(k.x2n-1)
S2n-1 = f(x1).sin(k.x1) + f(x3).sin(k.x3) + ......... + f(x2n-1).sin(k.x2n-1)

and     A(µ) = 1/µ + (sin(2µ))/(2µ2) - 2(sin2µ)/µ3
B(µ) =  2(1+cos2µ)/µ2 - 2(sin(2µ))/µ3
C(µ) =  4(sinµ)/µ3 - 4(cosµ)/µ2

-If  k = 0 , Simpson's rule is used

Data Registers:       •  R00 = Function name        ( Registers R00 thru R03 are to be initialized before executing "FIL" )
•  R01 = a
•  R02 = b         R06 = k             R07 to R12: temp
•  R03 = n
When the program stops:   R04 =   §ab  f(x).cos(k.x).dx  = X-register
R05 =   §ab  f(x).sin(k.x).dx   = Y-register
Flags:  /
Subroutine:   A program which computes  f(x)  assuming x is in X-register upon entry.

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

( 167 bytes / SIZE 013 )

 STACK INPUT OUTPUTS Y / §ab  f(x).sin(k.x).dx X k §ab  f(x).cos(k.x).dx

Example:   Evaluate   I = §16  Ln(x).cos(10x).dx    and   J = §16  Ln(x).sin(10x).dx

01  LBL "T"
02  LN
03  RTN

"T"  ASTO 00
1    STO 01
6    STO 02

with n = 8   ;    8   STO 03        10      R/S     >>>>    I =  -0.047890755
X<>Y   J =   0.175512930

with n = 16  ;  16  STO 03        10      R/S     >>>>    I =  -0.047429223
X<>Y   J =   0.174731804

with n = 32  ;  32  STO 03        10      R/S     >>>>    I =  -0.047453034
X<>Y   J =   0.174714501

with n = 64  ;  64  STO 03        10      R/S     >>>>    I =  -0.047454443     ( exact results are  I = -0.047454534
X<>Y   J =   0.174713854                           and  J =  0.174713817   to 9 decimals )

-The complete Filon's formulas involve the 3rd derivative of the function f
but this term is omitted here ( f'''(x) is difficult to evaluate on an HP-41 ... )

e) An Example of Curvilinear Integral

-The following program computes    §(C)  f(x,y).ds    where  (C) is the cicumference of the circle:  x2 + y2 = R2     (  ds = ( dx2+dy2)1/2 )
-The function f is evaluated 2n times.

Data Registers:        •  R00:  function name       ( Registers R00 thru R02 are to be initialized before executing "CINT" )
•  R01 = R
•  R02 = n
when the program stops:  R03 = §(C)  f(x,y).dx         R04-R05: temp

Flags: /
Subroutine:  A program which calculates  f(x,y)  assuming x is in X-register and y is in Y-register upon entry.

01  LBL "CINT"
02  RCL 02
03  ST+ X
04  STO 04
05  CLX
06  STO 03
07  SIGN
08  CHS
09  ACOS
10  RCL 02
11  /
12  STO 05
13  LBL 01
14  RCL 04
15  RCL 05
16  *
17  RCL 01
18  P-R
19  XEQ IND 00
20  ST+ 03
21  DSE 04
22  GTO 01
23  RCL 03
24  PI
25  *
26  RCL 01
27  *
28  RCL 02
29  /
30  STO 03
31  END

( 45 bytes / SIZE 006 )

Example:   Evaluate    §(C)  Ln(3+x.y).ds    where  (C):  x2 + y2 = 1

01  LBL "T"
02  *
03  3
04  +
05  LN
06  RTN

"T"  ASTO 00    R = 1  whence  1  STO 01

-with  n =  4    4   STO 02    XEQ "CINT"  >>>>  6.858533883
-with  n =  8    8   STO 02          R/S           >>>>  6.858689700
-with  n = 16  16  STO 02          R/S           >>>>  6.858689706

2°) Discrete Data

a) Equally Spaced Abscissas

a-1) Simpson's Rule

-We want to evaluate  §x1xn  f(x).dx  but we only know   f(x1) , f(x2) , ....... , f(xn)  where the  xi  are equally spaced     xi+1 - xi = h    for i = 1,2,...,n-1

-If n is odd, we use the formula:    §x1x3  f(x).dx  ~  h[f(x1)+4f(x2)+f(x3)]/3    over the intervals  [x1;x3] ,  [x3;x5] , ..... ,  [xn-2;xn]
-If n is even,  we use  §x1x4  f(x).dx  ~  3h[f(x1)+3f(x2)+3f(x3)+f(x4)]/8  and the same formula as above over  [x4;x6] , ..... ,  [xn-2;xn]

-These formulas produce perfect accuracy if f is a polynomial of degree < 4

-We assume  n > 2

Data Registers:       •  R00 = h                                                              ( Registers R00 thru Rnn  are to be initialized before executing "SIMP" )
•  R01 = f(x1)  •  R02 = f(x2)  .......  •  Rnn = f(xn)
Flags:  /
Subroutines:  /

-Synthetic register M may of course be replaced by any unused data register, for instance R99 ( provided n < 99 )

01  LBL "SIMP"
02  STO M
03  4
04  SQRT
05  RCL IND Y
06  CHS
07  LBL 01
08  RCL IND M
09  LASTX
10  X<> T
11  *
12  ST+ Y
13  RDN
14  DSE M
15  GTO 01
16  X<>Y
17  4
18  -
19  X=0?
20  GTO 02
21  CLX                  Lines 21 to 36 are only useful when n is even
22  RCL 02
23  11
24  *
25  RCL 01
26  15
27  *
28  -
29  RCL 03
30  5
31  *
32  -
33  RCL 04
34  +
35  8
36  /
37  LBL 02
38  +
39  RCL 01
40  -
41  RCL 00
42  *
43  3
44  /
45  END

( 64 bytes / SIZE n+1 )

 STACK INPUT OUTPUT X n §x1xn  f(x).dx

where n = the number of points ,  n > 2

Example:   Evaluate  §0pi/2  f(x).dx  and  §05pi/12  f(x).dx  from the following values

 x 0 pi/12 pi/6 pi/4 pi/3 5pi/12 pi/2 f(x) 0 0.2588190 0.5 0.7071068 0.8660254 0.9659258 1

0  STO 01   0.2588190  STO 02  ................  1  STO 07

PI  12  /   STO 00

7   XEQ "SIMP"  >>>>  1.0000263   ( exact result = §0pi/2  sin(x).dx = 1 )
6         R/S          >>>>   0.7412102   ( exact result = §05pi/12  sin(x).dx = 0.7411810 )

a-2) Other Newton-Cotes Formulae

-We assume  n = 6k+1  ( where k is an integer )   The  xi  are still  equally spacedxi+1 - xi = h    for i = 1,2,...,n-1

Formula:   §x1x7  f(x).dx  ~  h[41f(x1)+216f(x2)+27f(x3)+272f(x4)+27f(x5)+216f(x6)+41f(x7)]/140

Data Registers:       •  R00 = h                                                                ( Registers R00 thru Rnn  are to be initialized before executing "NC6" )
•  R01 = f(x1)  •  R02 = f(x2)  .......  •  Rnn = f(xn)
Flags:  /
Subroutines:  /

01  LBL "NC6"
02  RCL IND X
03  41
04  *
05  DSE Y
06  LBL 01
07  RCL IND Y
08  216
09  *
10  +
11  DSE Y
12  RCL IND Y
13  27
14  *
15  +
16  DSE Y
17  RCL IND Y
18  272
19  *
20  +
21  DSE Y
22  RCL IND Y
23  27
24  *
25  +
26  DSE Y
27  RCL IND Y
28  216
29  *
30  +
31  DSE Y
32  RCL IND Y
33  82
34  *
35  +
36  DSE Y
37  GTO 01
38  RCL 01
39  41
40  *
41  -
42  RCL 00
43  *
44  140
45  /
46  END

( 82 bytes / SIZE nnn+1 )

 STACK INPUT OUTPUT X n §x1xn  f(x).dx

where n = the number of points ;  n = 7 , 13 , 19 , .....

Example:   Evaluate  §0pi/2  f(x).dx   from the following values

 x 0 pi/12 pi/6 pi/4 pi/3 5pi/12 pi/2 f(x) 0 0.2588190 0.5 0.7071068 0.8660254 0.9659258 1

0  STO 01   0.2588190  STO 02  ................  1  STO 07

PI  12  /   STO 00
7   XEQ "NC6"  >>>>  1.0000000   ( the result is correct to 7 decimals! )

-Newton-Cotes 10 point formula is also interesting for sets of (9k+1) equally spaced arguments ( if you want to write a "NC9" program ):

§x1x10  f(x).dx  ~  9h[2857(f(x1)+f(x10))+15741(f(x2)+f(x9))+1080(f(x3)+f(x8))+19344(f(x4)+f(x7))+5778(f(x5)+f(x6))]/89600

-Higher degree formulas are theoretically better, but they involve large coefficients of alternate signs which may lead to poor accuracy...

b) Unequally Spaced Abscissas

b-1) Trapezoidal Rule

-We assume a function f is only known by n data points

 x1 x2 ........... xn f(x1) f(x2) ........... f(xn)

-The trapezoidal rule is the easiest formula to evaluate  §x1xn  f(x).dx
-We simply add the areas of n trapezoids:  Sumi=1,2,...,n-1   (xi+1-xi).[ f(xi+1)+f(xi) ]/2

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

•  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
•  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
Flags:  /
Subroutines:  /

01  LBL "TRAP"
02  RCL 00
03  DSE X
04  ST+ X
05  STO M
06  CLX
07  LBL 01
08  RCL M
09  2
10  +
11  X<>Y
12  RCL IND Y
13  RCL IND M
14  +
15  DSE Z
16  DSE M
17  RCL IND Z
18  RCL IND M
19  -
20  *
21  +
22  DSE M
23  GTO 01
24  2
25  /
26  END

( 47 bytes / SIZE 2n+1 )

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

where n = the number of points.

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

 x 1 2.4 4 5.2 7 8 f(x) 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   ( f(x) )    respectively

-There are 6 points so,     6   STO 00    XEQ "TRAP"    >>>>    §18  f(x).dx  ~  29.2

b-2) Connecting Parabolic Segments

-More accurate results may be obtained if we use several connected parabolic segments to compute  §x1xn  f(x).dx
-However, if n is even,  §x1x2  f(x).dx  is calculated by a polynomial of degree 3
-Thus, "ITG" yields the same results as "SIMP"  when the  xi's are equally spaced.

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

•  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
•  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
Flags:  /
Subroutines:  /

-Synthetic registers M N O P may be replaced by any unused data registers, for instance R96  R97  R98  R99 if  n < 48
-In this case, replace line 02 ( CLA ) by  CLX  STO 98

01  LBL "ITG"
02  CLA
03  RCL 00
04  2
05  -
06  ST+ X
07  STO M
08  4
09  MOD
10  X#0?                        Lines 12 to 86 are only useful when n is even. A slightly less accurate result may be obtained by replacing lines 12 to 86 by:
11  GTO 01
12  RCL 01                     RCL 06
13  RCL 03                     RCL 05
14  +                                RCL 01
15  2                                -
16  /                                 RCL 03
17  STO O                      RCL 05
18  RCL 05                     -
19  -                                *
20  STO N                      /
21  RCL 07                     RCL 03
22  RCL 01                     RCL 01
23  -                                -
24  /                                STO O
25  RCL 07                     *
26  RCL 03                     RCL 04
27  -                                RCL 05
28  /                                RCL 03
29  RCL 08                     -
30  *                                /
31  RCL O                      +
32  RCL 07                     RCL 02
33  -                                RCL 05
34  ST* N                       RCL 01
35  RCL 05                     -
36  RCL 01                     /
37  -                                -
38  /                                RCL O
39  RCL 05                     *
40  RCL 03                     RCL 02
41  -                                RCL 04
42  /                                 +
43  RCL 06                      3
44  *                                *
45  -                                 +
46  RCL 05                      ST* O
47  RCL 07
48  -
49  /
50  RCL 03
51  RCL 01
52  -
53  X^2
54  *
55  STO O
56  RCL 02
57  RCL 01
58  RCL 05
59  -
60  /
61  RCL 01
62  RCL 07
63  -
64  /
65  RCL 04
66  RCL 03
67  RCL 05
68  -
69  /
70  RCL 03
71  RCL 07
72  -
73  /
74  +
75  RCL N
76  *
77  ST+ X
78  RCL 02
79  +
80  RCL 04
81  +
82  ST+ O
83  RCL 03
84  RCL 01
85  -
86  ST* O
87  LBL 01
88  DSE M
89  4
90  RCL M
91  ST+ Y
92  2
93  +
94  RCL IND M
95  RCL IND Y
96  -
97  STO P                   ( synthetic )
98  RCL IND Z
99  LASTX
100  -
101  STO N
102  /
103  PI
104  INT
105  ST+ Z
106  DSE X
107  +
108  RCL IND Y
109  *
110  RCL P
111  RCL N
112  -
113  ST* X
114  LASTX
115  /
116  RCL P
117  /
118  DSE Z
119  DSE Z
120  RCL IND Z
121  *
122  -
123  RCL N
124  RCL P
125  ST- N
126  /
127  2
128  ST- T
129  ST- M
130  +
131  RCL IND Z
132  *
133  +
134  RCL N
135  *
136  ST+ O
137  DSE M
138  GTO 01
139  X<> O
140  6
141  /
142  CLA
143  END

( 192 bytes / SIZE 2n+1 )

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

where n = the number of points.  ( n > 2 )

Example:     Calculate  §17  f(x).dx  and  §18  f(x).dx   from the following values

 x 1 2.4 4 5.2 7 8 f(x) 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   ( f(x) )    respectively

-Five points:   5   STO 00   XEQ "ITG"  >>>>  §17  f(x).dx  ~  26.4226
-Six points:     6   STO 00        R/S         >>>>  §18  f(x).dx  ~  30.5339

-With the modification listed on the right, it gives  §18  f(x).dx  ~  30.6457

Remark:   We can also use Lagrange interpolation to obtain equally spaced arguments ( cf for instance "Lagrange Interpolation Formula for the HP-41" )

-With this set of 6 data points, it yields:

 x 1 2 3 4 5 6 7 8 f(x) 1 2.857 5.3453 6 5.2069 4.3568 4 2

and "SIMP" now produces   §18  f(x).dx  ~  29.6997

-In this example, we can even use the Newton-Cotes 8 point formula

§x1x8  f(x).dx  ~  7h[751f(x1)+3577f(x2)+1323f(x3)+2989f(x4)+2989f(x5)+1323f(x6)+3577f(x7)+751f(x8)]/17280
§x1x8  f(x).dx  ~  29.6179

which is probably the best evaluation of this integral without further information about the function.

-For larger sets of data points, one may likewise use Lagrange interpolation to transform unequally spaced abscissas
into (6k+1) or (9k+1) equally spaced abcsissas, and then apply "NC6" or "NC9" ( cf 2°) a-2) above ).
-Take for instance successive groups of  7 or 10 points with "LAGR"
but too many points could lead to worse rather than better accuracy!

b-3) Connecting Cubic Segments  ( X-Functions module required )

-The following program uses several segments of cubic curves to evaluate  §x1xn  f(x).dx
so that the method produces a perfect accuracy if the function is a polynomial of degree < 4

-Registers R00 to R2n are temporarily modified during the calculations, but their contents are finally restored.

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

•  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn
•  R02 = f(x1)  •  R04 = f(x2)  .......  •  R2n   = f(xn)
Flags:  /
Subroutines:  /

01  LBL "ITG3"
02  CLA
03  RCL 00
04  RCL 00
05   E3
06  ST/ 00
07  SIGN
08  -
09  3
10  MOD
11  -
12  ST+ 00
13  LBL 01
14  ISG 00
15  FS? 30
16  GTO 02
17  RCL 01
18  RCL 03
19  +
20  2
21  /
22  STO O
23  RCL 05
24  -
25  STO N
26  RCL 07
27  RCL 01
28  -
29  /
30  RCL 07
31  RCL 03
32  -
33  /
34  RCL 08
35  *
36  RCL O
37  RCL 07
38  -
39  ST* N
40  RCL 05
41  RCL 01
42  -
43  /
44  RCL 05
45  RCL 03
46  -
47  /
48  RCL 06
49  *
50  -
51  RCL 05
52  RCL 07
53  -
54  /
55  RCL 03
56  RCL 01
57  -
58  X^2
59  *
60  STO O
61  RCL 02
62  RCL 01
63  RCL 05
64  -
65  /
66  RCL 01
67  RCL 07
68  -
69  /
70  RCL 04
71  RCL 03
72  RCL 05
73  -
74  /
75  RCL 03
76  RCL 07
77  -
78  /
79  +
80  RCL N
81  *
82  ST+ X
83  RCL 02
84  +
85  RCL 04
86  +
87  RCL O
88  +
89  RCL 03
90  RCL 01
91  -
92  *
93  ST+ M
94  XEQ 00
95  GTO 01
96  LBL 00
97  RCL 00
98  FRC
99   E-3
100  ST- Y
101  *
102  ST+ X
103  3
104  LASTX
105  +
106  +
107  REGSWAP
108  RTN
109  LBL 02
110  RCL 00
111  2
112  -
113  3
114  /
115  INT
116  ST- 00
117  LBL 03
118  RCL 04
119  RCL 01
120  RCL 03
121  -
122  STO N
123  LASTX
124  RCL 07
125  -
126  *
127  /
128  RCL 06
129  RCL 01
130  RCL 05
131  -
132  STO O
133  /
134  RCL 05
135  RCL 07
136  -
137  /
138  -
139  RCL 01
140  RCL 07
141  -
142  *
143  RCL 03
144  RCL 05
145  -
146  /
147  RCL 02
148  RCL N
149  /
150  RCL O
151  /
152  -
153  RCL 08
154  RCL 03
155  RCL 07
156  -
157  /
158  RCL 05
159  RCL 07
160  -
161  /
162  -
163  2
164  /
165  RCL O
166  RCL 04
167  *
168  RCL N
169  /
170  RCL 03
171  RCL 05
172  -
173  /
174  RCL 03
175  RCL 07
176  -
177  /
178  -
179  RCL N
180  RCL 06
181  *
182  RCL O
183  /
184  RCL 03
185  RCL 05
186  -
187  /
188  RCL 05
189  RCL 07
190  -
191  /
192  +
193  RCL 01
194  RCL 07
195  -
196  *
197  RCL N
198  1/X
199  RCL O
200  1/X
201  +
202  RCL 02
203  *
204  +
205  RCL 03
206  RCL 07
207  -
208  1/X
209  RCL 05
210  RCL 07
211  -
212  1/X
213  +
214  RCL 08
215  *
216  +
217  RCL 01
218  RCL 07
219  -
220  *
221  RCL 02
222  RCL 08
223  +
224  3
225  *
226  -
227  RCL 01
228  RCL 07
229  -
230  *
231  ST+ M
232  XEQ 00
233  XEQ 00
234  XEQ 00
235  ISG 00
236  GTO 03                 a three-byte GTO
237  REGSWAP
238  RCL 00
239  INT
240  DSE X
241  STO 00
242  X<> M
243  6
244  /
245  CLA
246  END

( 302 bytes / SIZE 2n+1 )

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

where n = the number of points.  ( n > 3 )

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

 x 1 2.4 4 5.2 7 8 f(x) 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   ( f(x) )    respectively

-There are 6 points,     6   STO 00    XEQ "ITG3"    >>>>    §18  f(x).dx  ~  30.2135   ( in 14 seconds )

b-4) Natural Cubic Spline

-In the cubic spline integration, we also connect arcs of cubic polynomials,
but in such a way that the first derivative y' and the second derivative y" are continuous over the whole interval  [ x1 , xn ]
-The spline is called "natural" if we set  y"1 = y"n = 0

-This leads to a tridiagonal linear system of  (n-2)  equations in  (n-2)  unknowns.

(1/6).( xk - xk-1 ).y"k-1 + (1/3).( xk+1 - xk-1 ).y"k + (1/6).( xk+1 - xk ).y"k+1 =  ( yk+1 - yk )/( xk+1 - xk ) - ( yk - yk-1 )/( xk - xk-1 )

-Then    §xkxk+1  y.dx = ( xk+1 - xk ).( yk+1 + yk )/2 - ( xk+1 - xk )3.( y"k+1 + y"k )/24

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

•  R01 = x1     •  R03 = x2      .......  •  R2n-1 = xn             R2n+1 thru R6n-7: temp
•  R02 = y1     •  R04 = y2      .......   •  R2n   = yn

>>>>   If n > 2 , when the program stops:  R5n-8 = y"1 = 0 ,  R5n-7 = y"2 , ............... ,  R6n-10 = y"n-1 , R6n-9 = y"n = 0

Flags:  /
Subroutine:   none if  n < 4
otherwise, "3DLS"  "Tridiagonal Linear Systems"  ( cf "Linear and non-Linear Systems for the HP-41" )

01  LBL "NCSI"
02  RCL 00
03  2
04  X#Y?
05  GTO 00
06  RCL 02                         Lines 06 to 14 simply apply the trapezoidal rule if n = 2
07  RCL 04
08  +
09  RCL 03
10  RCL 01
11  -
12  *
13  2
14  /
15  RTN
16  LBL 00
17  3
18  STO M
19  RCL 00
20  ST+ X
21  1
22  +
23  RCL 00
24  5
25  *
26  7
27  -
28  STO O
29  DSE X
30   E3
31  /
32  +
33  STO N
34  RCL 05
35  RCL 01
36  -
37  ST+ X
38  STO IND N
39  RCL 06
40  RCL 04
41  -
42  RCL 05
43  RCL 03
44  -
45  /
46  RCL 04
47  RCL 02
48  -
49  RCL 03
50  RCL 01
51  -
52  /
53  -
54  6
55  *
56  STO IND O
57  LBL 01
58  ISG N
59  FS? 30
60  GTO 02
61  2
62  RCL M
63  ST+ Y
64  4
65  +
66  RCL IND Y
67  RCL IND M
68  -
69  STO P                      ( synthetic )
70  STO IND N
71  ISG N
72  STO IND N
73  CLX
74  RCL IND Y
75  LASTX
76  -
77  STO Q
78  ST+ X
79  ISG N
80  STO IND N
81  ISG M
82  CLX
83  CLX
84  SIGN
85  -
86  RCL IND M
87  RCL IND Y
88  -
89  RCL P
90  ST- Q
91  /
92  ISG O
93  CLX
94  STO IND O
95  X<> Q
96  RCL Y
97  2
98  +
99  RDN
100  RCL IND T
101  RCL IND Z
102  -
103  X<>Y
104  /
105  ST+ IND O
106  6
107  ST* IND O
108  SIGN
109  ST+ M
110  GTO 01
111  LBL 02
112  RCL 00
113  3
114  X#Y?
115  GTO 03
116  CLX                         If n = 3 the "system" has only one equation which is solved directly.
117  STO 09
118  RCL 07
119  ST/ 08
120  7.009
121  GTO 04
122  LBL 03
123  SIGN
124  ST+ O
125  CLX
126  STO IND O
127  RCL M
128  4
129  +
130  STO Y
131  RCL 00
132  LASTX
133  *
134  +
135  11
136  -
137   E3
138  /
139  +
140  XEQ "3DLS"
141  LASTX                    R00 is disturbed by "3DLS"
142  2                              but its content is now restored
143  +                              R00 = n
144  STO 00
145  CLX
146  .999
147  -
148  LBL 04
149  CLA
150  STO N
151  CLX
152  STO IND N
153  RCL 00
154  ST+ X
155   E3
156  /
157  3
158  +
159  STO O
160  LBL 05
161  RCL O
162  2
163  -
164  RCL IND O
165  RCL IND Y
166  -
167  RCL IND N
168  ISG N
169  RCL IND N
170  +
171  RCL Y
172  X^2
173  *
174  ISG Z
175  RCL IND Z
176  ISG O
177  RCL IND O
178  +
179  12
180  *
181  -
182  *
183  ST- M
184  ISG O
185  GTO 05
186  RCL M
187  24
188  /
189  CLA
190  END

( 283 bytes / SIZE 6n-8 )

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

where n = the number of points.

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

 x 1 2.4 4 5.2 7 8 f(x) 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   ( f(x) )    respectively

-SIZE 028
-There are 6 points,     6   STO 00    XEQ "NCSI"    >>>>    §18  f(x).dx  =  29.99938860   ( in 21 seconds )

-And we have  R22 = y"1 = 0 ,  R23 = y"2 = -0.237729622 , R24 = y"3 = -2.456728203 ,
R25 = y"4 = 1.365037775 ,  R26 = y"5 = -1.986381189 ,  R27 = y"6 = 0

b-5) Integration of the Lagrange Polynomial

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

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

-The Lagrange polynomial  is the unique polynomial L(x) of degree < n  such that:     L(xi) = yi       i = 1 , 2 , ..... , n
-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 "Lagrange's Interpolation Formula for the HP-41" )

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 )
( this value was already found more laboriously in paragraph 2°) b-2) above )

-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

c) Double Integrals  ( Equally Spaced Arguments )

-We want to evaluate  §x1xn  §y1ym  f(x,y).dx.dy    assuming  n & m are odd integers  ( n > 1 , m > 1 )
-The following program uses Simpson's rule in the direction of  x-axis and y-axis

-  xi  must be  equally spaced:  xi+1 - xi = h     for i = 1,2,...,n-1
-  yj  must be  equally spaced:  yj+1 - yj = k    for j = 1,2,...,m-1

Data Registers:       •  R00   = h.k                                              ( Registers R00 thru Rn.m  are to be initialized before executing "INT2" )

•  R01   = f(x1,y1)  •  Rn+1 = f(x1,y2)  .......  •  Rnm-n+1 = f(x1,ym)
•  R02   = f(x2,y1)  •  Rn+2 = f(x2,y2)  .......  •  Rnm-n+2 = f(x2,ym)
.......................................................................................
•  Rnn   = f(xn,y1)  •  R2n = f(xn,y2)    .......  •  Rnm = f(xn,ym)
Flags:  /
Subroutines:  /

-Synthetic registers M N O may be replaced by any unused data registers, for instance  R97  R98  R99 if  n.m < 97

01  LBL "INT2"
02  STO N
03  X<>Y
04  STO M
05  *
06  STO O
07  CLX
08  LBL 01
09  RCL O
10  RCL N
11  MOD
12  2
13  X>Y?
14  CLX
15  MOD
16  LASTX
17  X<>Y
18  -
19  ST+ X
20  X<=0?
21  SIGN
22  ABS
23  RCL O
24  1
25  -
26  RCL N
27  /
28  INT
29  1
30  +
31  RCL M
32  MOD
33  2
34  X>Y?
35  CLX
36  MOD
37  LASTX
38  X<>Y
39  -
40  ST+ X
41  X<=0?
42  SIGN
43  ABS
44  *
45  RCL IND O
46  *
47  +
48  DSE O
49  GTO 01
50  RCL 00
51  *
52  9
53  /
54  CLA
55  END

( 77 bytes / SIZE n.m+1 )

 STACK INPUTS OUTPUTS Y m / X n §x1xn §y1ym f(x,y).dx.dy

Example:     Calculate   §26 §15  f(x,y).dx.dy   from the following f(x,y) values

 x\y 1 2 3 4 5 2 3 4 7 6 3 4 1 2 4 5 3 6 4 1 3 4 6

3  4  7  6  3                                   R01  R04  R07  R10  R13
We store     1  2  4  5  3        into registers        R02  R05  R08  R11  R14      respectively
4  1  3  4  6                                   R03  R06  R09  R12  R15

-Here, we have  h = 2 and  k = 1   whence  1*2 = 2  STO 00

n = 3  and  m = 5    whence    5  ENTER^
3  XEQ "INT2"  >>>>   §26 §15  f(x,y).dx.dy  ~  56.8889

-Perfect accuracy is achieved if f is a polynomial and  deg(f) < 4

d) Triple Integrals  ( Equally Spaced Arguments )

-Now we evaluate  §x1xn  §y1ym  §z1zp  f(x,y,z).dx.dy.dz    assuming  n , m , p are odd integers  ( n > 1 , m > 1 , p > 1 )
-Simpson's rule in the direction of  x- , y- and z-axis is used.

-  xi  must be  equally spaced:  xi+1 - xi = h         for i = 1,2,...,n-1
-  yj  must be  equally spaced:  yj+1 - yj = h'       for j = 1,2,...,m-1
-  zk  must be  equally spaced:  zk+1 - yk = h''    for j = 1,2,...,p-1

Data Registers:       •  R00   = h.h'.h''                                                ( Registers R00 thru Rn.m.p  are to be initialized before executing "INT3" )

•  R01   = f(x1,y1,z1)  •  Rn+1 = f(x1,y2,z1)  .......  •  Rnm-n+1 = f(x1,ym,z1)
•  R02   = f(x2,y1,z1)  •  Rn+2 = f(x2,y2,z1)  .......  •  Rnm-n+2 = f(x2,ym,z1)
......................................................................................................
•  Rn     = f(xn,y1,z1)  •  R2n = f(xn,y2,z1)    .......  •  Rnm = f(xn,ym,z1)

•  Rnm+1   = f(x1,y1,z2)  •  Rnm+n+1 = f(x1,y2,z2)  .......  •  R2nm-n+1 = f(x1,ym,z2)
•  Rnm+2   = f(x2,y1,z2)  •  Rnm+n+2 = f(x2,y2,z2)  .......  •  R2nm-n+2 = f(x2,ym,z2)
......................................................................................................
•  Rnm+n   = f(xn,y1,z2)  •  Rnm+2n = f(xn,y2,z2)    .......  •  R2nm = f(xn,ym,z2)

.............................................................................................................
.............................................................................................................
.............................................................................................................

•  Rnm(p-1)+1   = f(x1,y1,zp)  •  Rnm(p-1)+n+1 = f(x1,y2,zp)  .......  •  Rnmp-n+1 = f(x1,ym,zp)
•  Rnm(p-1)+2   = f(x2,y1,zp)  •  Rnm(p-1)+n+2 = f(x2,y2,zp)  .......  •  Rnmp-n+2 = f(x2,ym,zp)
......................................................................................................
•  Rnm(p-1)+n   = f(xn,y1,zp)  •  Rnm(p-1)+2n = f(xn,y2,zp)    .......  •  Rnmp = f(xn,ym,zp)
Flags:  /
Subroutines:  /

-Synthetic registers M N O P Q may be replaced by any unused data registers, for instance  R95 R95 R97  R98  R99 if  n.m.p < 95

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

( 124 bytes / SIZE nmp+1 )

 STACK INPUTS OUTPUTS Z p / Y m / X n §x1xn§y1ym§z1zp f(x,y,z)dx.dy.dz

Example:     Evaluate   §13 §15 §17  f(x,y,z).dx.dy.dz   from the following values

f(1,1,1) =  4      f(1,3,1) =  6     f(1,5,1) =  8                                                     R01   R04   R07
f(2,1,1) =  7      f(2,3,1) =  9     f(2,5,1) = 11           to be stored in registers      R02   R05   R08     respectively
f(3,1,1) = 10     f(3,3,1) = 12    f(3,5,1) = 14                                                    R03   R06   R09

f(1,1,4) =  64       f(1,3,4) =  96      f(1,5,4) = 128                                                    R10   R13   R16
f(2,1,4) = 112      f(2,3,4) = 144     f(2,5,4) = 176           to be stored in registers      R11   R14   R17     respectively
f(3,1,4) = 160      f(3,3,4) = 192     f(3,5,4) = 224                                                    R12   R15   R18

f(1,1,7) =  196      f(1,3,7) = 294     f(1,5,7) = 392                                                    R19   R22   R25
f(2,1,7) =  343      f(2,3,7) = 441     f(2,5,7) = 539           to be stored in registers      R20   R23   R26     respectively
f(3,1,7) =  490      f(3,3,7) = 588     f(3,5,7) = 686                                                    R21   R24   R27

h.h'.h'' = 1*2*3 = 6     6  STO 00

n = m = p = 3   whence   3  ENTER^   ENTER^   XEQ "INT3"   >>>>   §13 §15 §17  f(x,y,z).dx.dy.dz  =  8208

-These values were actually computed from  f(x,y,z) =  (3x+y).z2   and the result is perfectly accurate since f is a polynomial and deg(f) < 4.

References:

[1]  Abramowitz and Stegun , "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  F.Scheid - "Numerical Analysis" - Mc Graw Hill - ISBN  07-055197-9
[3]  Press , Vetterling , Teukolsky , Flannery , "Numerical Recipes in C++" - Cambridge University Press - ISBN  0-521-75033-4