SPHINTEG

# Integrals over a Sphere & its Surface for the HP-41

Overview

1°)  Integrals over a Sphere

1-a)  Program#1
1-b)  Program#2

2°)  Integrals on the Surface of a Sphere

2-a)  Program#1
2-b)  Program#2

3°)  Integrals over the Unit Hypersphere

3-a)  Program#1 ( 5th degree of accuracy )
3-b)  Program#2 ( 7th degree of accuracy )

4°)  Integrals on the Surface of the Unit Hypersphere

4-a)  Program#1 ( 5th degree of accuracy )
4-b)  Program#2 ( 7th degree of accuracy )

-We could employ usual formulae - like Gaussian integration - in several directions
-But specific formulae for spheres and hyperspheres require much less evaluations of the function.

1°)  Integrals over a Sphere

a)  Program#1

-We want to estimate   I = §§§x^2+y^2+z^2<=R^2    f(x,y,z)  dx dy dz

Data Registers:           •  R00 = function name                  ( Register R00 is to be initialized before executing "I3SPH" )

R01 = Integral   R02 = R
Flags: /
Subroutine:  A program that takes x , y , z in registers  X , Y , Z  respectively and returns  f(x,y,z) in X-register.

 01  LBL "I3SPH"  02  STO 02  03  CLST  04  XEQ IND 00  05  4  06  *  07  STO 01  08  CLST  09  RCL 02  10  XEQ IND 00 11  ST+ 01  12  CLST  13  RCL 02  14  CHS  15  XEQ IND 00  16  ST+ 01  17  CLST  18  RCL 02  19  X<>Y  20  XEQ IND 00 21  ST+ 01  22  CLST  23  RCL 02  24  CHS  25  X<>Y  26  XEQ IND 00  27  ST+ 01  28  RCL 02  29  0  30  ENTER 31  XEQ IND 00  32  ST+ 01  33  RCL 02  34  CHS  35  0  36  ENTER  37  XEQ IND 00  38  ST+ 01  39  RCL 02  40  3 41  Y^X  42  .4  43  *  44  PI  45  *  46  3  47  /  48  ST* 01  49  RCL 01         50  END

( 75 bytes / SIZE 003 )

 STACK INPUT OUTPUT X R Integral

where R = radius of the sphere

Example:     I = §§§x^2+y^2+z^2<=R^2    Ln ( 16 + x + y2 + z3 )  dx dy dz     with  R = 2

-This short routine calculates  Ln ( 16 + x + y2 + z3 )

 01  LBL "T" 02  X<>Y 03  X^2 04  + 05  X<>Y 06  3 07  Y^X 08  + 09  16 10  + 11  LN 12  RTN

"T"  ASTO 00

2   XEQ "I3SPH"  >>>>  I = 93.3891

b)  Program#2

-The same integral  I = §§§x^2+y^2+z^2<=R^2    f(x,y,z)  dx dy dz  is evaluated with a better accuracy

Data Registers:           •  R00 = function name                  ( Register R00 is to be initialized before executing "I3SPH" )

R01 = Integral   R02 = R   R03-R04: temp
Flags: /
Subroutine:  A program that takes x , y , z in registers  X , Y , Z  respectively and returns  f(x,y,z) in X-register.

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

( 317 bytes / SIZE 005 )

 STACK INPUT OUTPUT X R Integral

where R = radius of the sphere

Example:    I = §§§x^2+y^2+z^2<=R^2    Ln ( 16 + x + y2 + z3 )  dx dy dz     with  R = 2

With the same subroutine,

"T"  ASTO 00

2   XEQ "I3SPH"  >>>>  I = 94.2545     ( instead of  93.3891  with the 1st program )

-A substantial improvement since the exact result - rounded to 4 decimals - is  94.2541

2°)  Integrals on the Surface of a Sphere

a)  Program#1

-Now we want to estimate   I = §§x^2+y^2+z^2=R^2    f(x,y,z)  ds    where ds is the element of surface on a sphere

Data Registers:           •  R00 = function name                             ( Register R00 is to be initialized before executing "I2SPH" )

R01 = Integral   R02 = R   R03-R04: temp
Flags: /
Subroutine:  A program that takes x , y , z in registers  X , Y , Z  respectively and returns  f(x,y,z) in X-register.

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

( 245 bytes / SIZE 005 )

 STACK INPUT OUTPUT X R Integral

where R = radius of the sphere

Example:    I = §§x^2+y^2+z^2=R^2    Ln ( 16 + x + y2 + z3 )  dx dy dz     with  R = 2

With the same subroutine,

"T"  ASTO 00

2   XEQ "I2SPH"  >>>>  I = 142.1916

b)  Program#2

-The same integral is now estimated with a better precision

Data Registers:           •  R00 = function name                                         ( Register R00 is to be initialized before executing "I2SPH" )

R01 = Integral   R02 = R   R03-R04: temp
Flags: /
Subroutine:  A program that takes x , y , z in registers  X , Y , Z  respectively and returns  f(x,y,z) in X-register.

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

( 430 bytes / SIZE 006 )

 STACK INPUT OUTPUT X R Integral

where R = radius of the sphere

Example:    I = §§x^2+y^2+z^2=R^2    Ln ( 16 + x + y2 + z3 )  dx dy dz     with  R = 2

With the same subroutine,

"T"  ASTO 00

2   XEQ "I2SPH"  >>>>  I = 142.2271    ( instead of  142.1916  with the 1st program )

-The exact value is  I = 142.2311    rounded to 4 decimals.

3°)  Integrals over the Unit HyperSphere

a)  Program#1 ( 5th degree of accuracy )

"IHS5" evaluates   I = §Sn   f(x1,x2,.....,xn)  dx1 dx2 ...... dxn   where  (Sn):  x12 + x22 + ........ + xn2 <= 1

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

R01 = x1   R02 = x2   .....................   Rnn = xn
Flags: /
Subroutine:  A program LBL "T"  that takes x1 , x2 , ............ ,  xn  in registers  R01 R02 ...... Rnn   respectively
and returns  f(x1,x2,.....,xn) in X-register.

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

( 234 bytes / SIZE 006 )

 STACK INPUT OUTPUT X / I

where  I  = §Sn   f(x1,x2,.....,xn)  dx1 dx2 ...... dxn   with   (Sn):  x12 + x22 + ........ + xn2 <= 1

Example1:       I = §§§§S4  ( x2 + y2 + z2 + u2 )  dx dy dz du

 01  LBL "T"  02  RCL 01  03  X^2  04  RCL 02  05  X^2  06  +  07  RCL 03  08  X^2  09  +  10  RCL 04  11  X^2  12  +  13  RTN

4  STO 00   XEQ "IHS5"   >>>>   I = 3.289868134              ---Execution time = 27s---

-The exact result is PI^2 / 3 , so all the decimals are correct ( rounded to 9 D )
-The formula gives perfect accuracy if f is a polynomial of degree < 6

Example2:       I = §§§§§§S6  Ln (  PI2 + x + y2 + z3 + u4 + v5 + w6 )  dx dy dz du dv dw

 01  LBL "T"  02  RCL 01  03  RCL 02  04  X^2  05  +  06  RCL 03  07  3  08  Y^X  09  +  10  RCL 04  11  X^2  12  X^2  13  +  14  RCL 05  15  5  16  Y^X  17  +  18  RCL 06  19  6  20  Y^X  21  +  22  PI  23  X^2  24  +  25  LN  26  RTN

6  STO 00   XEQ "IHS5"   >>>>   I = 11.9174              ---Execution time = 2mn15s---

Notes:

-The formula is exact if f is a polynomial of degree < 6

-The subroutine must begin by LBL "T" otherwise change lines 23-28-88-93-107
-The function f is evaluated ( 2 n2 + 1 ) times

-Registers R01 ..... Rnn  are cleared by this program

-Lines 121 to 142 may be replaced by XEQ "HS" ( cf "Hyperspheres for the HP-41' )
-They return 0 if n > 187, so don't use this program in this case.
-They actually calculate PIn/2/(n/2)!  so the program may be simplified if you have an M-Code Gamma function.

b)  Program#2 ( 7th degree of accuracy )

"IHS7" also evaluates   I = §Sn   f(x1,x2,.....,xn)  dx1 dx2 ...... dxn   where  (Sn):  x12 + x22 + ........ + xn2 <= 1  but with a better precision

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

R01 = x1   R02 = x2   .....................   Rnn = xn   Rn+1: temp
Flags: /
Subroutine:  A program LBL "T"  that takes x1 , x2 , ............ ,  xn  in registers  R01 R02 ...... Rnn   respectively
and returns  f(x1,x2,.....,xn) in X-register.

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

( 443 bytes / SIZE n+2 )

 STACK INPUT OUTPUT X / I

where  I  = §Sn   f(x1,x2,.....,xn)  dx1 dx2 ...... dxn   with   (Sn):  x12 + x22 + ........ + xn2 <= 1

Example1:       I = §§§§§§S6   Ln ( PI + x2 + y2 + z2 + u2 + v2 + w2 )  dx dy dz du dv dw

 01  LBL "T"  02  RCL 01  03  X^2  04  RCL 02  05  X^2  06  +  07  RCL 03  08  X^2  09  +  10  RCL 04  11  X^2  12  +  13  RCL 05  14  X^2  15  +  16  RCL 06  17  X^2  18  +  19  PI  20  +  21  LN  22  RTN

6  STO 00   XEQ "IHS7"   >>>>   I = 7.015497950              ---Execution time = 5m47s---

-The exact result is I = 7.015376313...

Example2:       I = §§§§§§S6  Ln (  PI2 + x + y2 + z3 + u4 + v5 + w6 )  dx dy dz du dv dw

 01  LBL "T"  02  RCL 01  03  RCL 02  04  X^2  05  +  06  RCL 03  07  3  08  Y^X  09  +  10  RCL 04  11  X^2  12  X^2  13  +  14  RCL 05  15  5  16  Y^X  17  +  18  RCL 06  19  6  20  Y^X  21  +  22  PI  23  X^2  24  +  25  LN  26  RTN

6  STO 00   XEQ "IHS7"   >>>>   I = 11.91901135              ---Execution time = 7mn47s---

Notes:

-The formula is exact if f is a polynomial of degree < 8
-The subroutine must begin by LBL "T"
-Registers R01 ..... Rnn  are cleared by this program

>>>  n must be at least 3

-The function f is evaluated ( 4 n3 - 6 n2 + 14 n + 3 ) / 3  times

 n 3 4 5 6 7 8 9 10 eval 33 73 141 245 393 593 853 1181

4°)  Integrals on the Surface of the Unit HyperSphere

a)  Program#1 ( 5th degree of accuracy )

"ISHS5" evaluates   I = §dSn   f(x1,x2,.....,xn)  ds                  where  (dSn):  x12 + x22 + ........ + xn2 = 1

-After solving the equations (1) to (4) given in reference [3], we find the following weights: w

for the 2n  points              ( x , 0 , ........... , 0 )       x2 = 1        w = A ( 4 - n ) / ( 2 n2 + 4 n )
for the 2n(2n-1) points      ( y , y , 0 , ..... , 0 )       y2 = 1/2     w =  A / ( n2 + 2 n )

where  A =  Surface of the hypersphere

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

R01 = x1   R02 = x2   .....................   Rnn = xn
Flags: /
Subroutine:  A program LBL "T"  that takes x1 , x2 , ............ ,  xn  in registers  R01 R02 ...... Rnn   respectively
and returns  f(x1,x2,.....,xn) in X-register.

 01  LBL "ISHS5"  02  RCL 00  03   E3  04  /  05  ISG X  06  CLRGX  07  .5  08  SQRT  09  STO 01  10  STO 02  11  CLA  12  SIGN  13  +  14  STO N  15  1.001  16  -  17  STO M  18  LBL 01  19  XEQ "T"  20  ST+ O  21  RCL IND M   22  CHS  23  STO IND M 24  X<0?  25  GTO 01  26  RCL IND N  27  CHS  28  STO IND N  29  X<0?  30  GTO 01  31  CLX  32  X<> IND N  33  ISG N  34  X=0?  35  GTO 00  36  STO IND N  37  GTO 01  38  LBL 00  39  CLX  40  X<> IND M   41  ISG M  42  X=0?  43  GTO 00  44  STO IND M  45  RCL M  46  INT 47  RCL N  48  FRC  49  +  50  ISG X  51  STO N  52  X<>Y  53  STO IND Y  54  GTO 01  55  LBL 00  56  1  57  STO 01  58  CLX  59  X<> N  60  FRC  61  ISG X  62  STO M  63  LBL 02  64  XEQ "T"  65  ST+ N  66  RCL IND M   67  CHS  68  STO IND M  69  X<0? 70  GTO 02  71  CLX  72  X<> IND M  73  ISG M  74  X=0?  75  GTO 00  76  STO IND M   77  GTO 02  78  LBL 00  79  X<> N  80  4  81  RCL 00  82  -  83  *  84  2  85  /  86  ST+ O  87  RCL 00  88  RCL 00  89  RCL 00  90  2  91  ST/ Z  92  MOD 93  1  94  +  95  X<>Y  96  INT  97  STO Z  98  LBL 03  99  CLX 100  PI 101  ST* Y 102  X<> L 103  ST/ Y 104  SIGN 105  ST- L 106  DSE Z 107  GTO 03  108  X<> O 109  * 110  2 111  RCL 00        112  + 113  / 114  CLA 115  END

( 184 bytes / SIZE n+1 )

 STACK INPUT OUTPUT X / I

where  I  = §dSn   f(x1,x2,.....,xn)  ds              with   (dSn):  x12 + x22 + ........ + xn2 = 1

Example:    Calculate  I = §§§dS4   Ln ( Pi + x2 y + z t ) ds

 01  LBL "T"  02  RCL 01  03  X^2  04  RCL 02  05  *  06  RCL 03  07  RCL 04  08  *  09  +  10  PI  11  +  12  LN  13  RTN

4  STO 00

XEQ "ISHS5"  >>>>   I = 22.53289243                           ---Execution time = 35s---

Notes:

-The formula is exact if f is a polynomial of degree < 6
-The subroutine must begin by LBL "T"
-Registers R01 thru Rnn are cleared.
-There are 2 n2  evaluations of the function

b)  Program#2 ( 7th degree of accuracy )

"ISHS7" also evaluates   I = §dSn   f(x1,x2,.....,xn)  ds          where  (dSn):  x12 + x22 + ........ + xn2 = 1             but with a better precision.

-It generalizes the program listed in paragraph 2°)a)

-After solving the equations (1) to (7)  given in reference [3], we find the following weights: w

for the 2n  points                     ( x , 0 , ............... , 0 )       x2 = 1        w = A ( n2 - 9 n + 38 ) / 4 / ( n3 + 6 n2 + 8 n )
for the 2n(2n-1) points            ( y , y , 0 , ........... , 0 )       y2 = 1/2     w = A ( - 2 n + 10 ) / ( n3 + 6 n2 + 8 n )
for the 4n(n-1)(n-2)/3  points  ( z , z , z , 0 , ....... , 0 )       z2 = 1/3      w = 27 A / ( n3 + 6 n2 + 8 n )

where  A =  Surface of the hypersphere

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

R01 = x1   R02 = x2   .....................   Rnn = xn   Rn+1: temp
Flags: /
Subroutine:  A program LBL "T"  that takes x1 , x2 , ............ ,  xn  in registers  R01 R02 ...... Rnn   respectively
and returns  f(x1,x2,.....,xn) in X-register.

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

( 346 bytes / SIZE n+2 )

 STACK INPUT OUTPUT X / I

where  I  = §dSn   f(x1,x2,.....,xn)  ds              with   (dSn):  x12 + x22 + ........ + xn2 = 1

Example:    Calculate  I = §§§dS4   Ln ( Pi + x2 y + z t ) ds

 01  LBL "T"  02  RCL 01  03  X^2  04  RCL 02  05  *  06  RCL 03  07  RCL 04  08  *  09  +  10  PI  11  +  12  LN  13  RTN

4  STO 00

XEQ "ISHS7"  >>>>   I = 22.53840629                           ---Execution time = 80s---

Notes:

-The formula is exact if f is a polynomial of degree < 8
-The subroutine must begin by LBL "T"
-Registers R01 thru Rnn are cleared.
-There are ( 4 n3 - 6 n2 + 8 n ) / 3   evaluations of the function

 n 3 4 5 6 7 8 9 10 eval 26 64 130 232 378 576 834 1160

References:

[1]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[2]  Christopher A. Feuchter - "Numerical Integration Over a Sphere"
[3]  Preston C. Hammer and Arthur H. Stroud - "Numerical Evaluation of Multiple Integrals II"
[4]  L. N. Dobrodeev - A cubature formula of the seventh degree of accuracy for a hypersphere and a hypercube.