libm/double/simpsn.c

   1 /*                                                      simpsn.c        */
   2 /* simpsn.c
   3  * Numerical integration of function tabulated
   4  * at equally spaced arguments
   5  */
   6
   7 /* Coefficients for Cote integration formulas */
   8
   9 /* Note: these numbers were computed using 40-decimal precision. */
  10
  11 #define NCOTE 8
  12
  13 /* 6th order formula */
  14 /*
  15 static double simcon[] =
  16 {
  17   4.88095238095238095E-2,
  18   2.57142857142857142857E-1,
  19   3.2142857142857142857E-2,
  20   3.2380952380952380952E-1,
  21 };
  22 */
  23
  24 /* 8th order formula */
  25 static double simcon[] =
  26 {
  27   3.488536155202821869E-2,
  28   2.076895943562610229E-1,
  29  -3.27336860670194003527E-2,
  30   3.7022927689594356261E-1,
  31  -1.6014109347442680776E-1,
  32 };
  33
  34 /* 10th order formula */
  35 /*
  36 static double simcon[] =
  37 {
  38   2.68341483619261397039E-2,
  39   1.77535941424830313719E-1,
  40  -8.1043570626903960237E-2,
  41   4.5494628827962161295E-1,
  42  -4.3515512265512265512E-1,
  43   7.1376463043129709796E-1,
  44 };
  45 */
  46 \f
  47 /*                                                      simpsn.c 2      */
  48 /* 20th order formula */
  49 /*
  50 static double simcon[] =
  51 {
  52   1.182527324903160319E-2,
  53   1.14137717644606974987E-1,
  54  -2.36478370511426964E-1,
  55   1.20618689348187566E+0,
  56  -3.7710317267153304677E+0,
  57   1.03367982199398011435E+1,
  58  -2.270881584397951229796E+1,
  59   4.1828057422193554603E+1,
  60  -6.4075279490154004651555E+1,
  61   8.279728347247285172085E+1,
  62  -9.0005367135242894657916E+1,
  63 };
  64 */
  65 \f
  66 /*                                                      simpsn.c 3      */
  67 double simpsn( f, delta )
  68 double f[];     /* tabulated function */
  69 double delta;   /* spacing of arguments */
  70 {
  71 extern double simcon[];
  72 double ans;
  73 int i;
  74
  75
  76 ans = simcon[NCOTE/2] * f[NCOTE/2];
  77 for( i=0; i < NCOTE/2; i++ )
  78         ans += simcon[i] * ( f[i] + f[NCOTE-i] );
  79
  80 return( ans * delta * NCOTE );
  81 }