htswanalysis/MACS/lib/gsl/gsl-1.11/integration/qk61.c

   1 /* integration/qk61.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <gsl/gsl_integration.h>
  22
  23 /* Gauss quadrature weights and kronrod quadrature abscissae and
  24    weights as evaluated with 80 decimal digit arithmetic by
  25    L. W. Fullerton, Bell Labs, Nov. 1981. */
  26
  27 static const double xgk[31] =   /* abscissae of the 61-point kronrod rule */
  28 {
  29   0.999484410050490637571325895705811,
  30   0.996893484074649540271630050918695,
  31   0.991630996870404594858628366109486,
  32   0.983668123279747209970032581605663,
  33   0.973116322501126268374693868423707,
  34   0.960021864968307512216871025581798,
  35   0.944374444748559979415831324037439,
  36   0.926200047429274325879324277080474,
  37   0.905573307699907798546522558925958,
  38   0.882560535792052681543116462530226,
  39   0.857205233546061098958658510658944,
  40   0.829565762382768397442898119732502,
  41   0.799727835821839083013668942322683,
  42   0.767777432104826194917977340974503,
  43   0.733790062453226804726171131369528,
  44   0.697850494793315796932292388026640,
  45   0.660061064126626961370053668149271,
  46   0.620526182989242861140477556431189,
  47   0.579345235826361691756024932172540,
  48   0.536624148142019899264169793311073,
  49   0.492480467861778574993693061207709,
  50   0.447033769538089176780609900322854,
  51   0.400401254830394392535476211542661,
  52   0.352704725530878113471037207089374,
  53   0.304073202273625077372677107199257,
  54   0.254636926167889846439805129817805,
  55   0.204525116682309891438957671002025,
  56   0.153869913608583546963794672743256,
  57   0.102806937966737030147096751318001,
  58   0.051471842555317695833025213166723,
  59   0.000000000000000000000000000000000
  60 };
  61
  62 /* xgk[1], xgk[3], ... abscissae of the 30-point gauss rule.
  63    xgk[0], xgk[2], ... abscissae to optimally extend the 30-point gauss rule */
  64
  65 static const double wg[15] =    /* weights of the 30-point gauss rule */
  66 {
  67   0.007968192496166605615465883474674,
  68   0.018466468311090959142302131912047,
  69   0.028784707883323369349719179611292,
  70   0.038799192569627049596801936446348,
  71   0.048402672830594052902938140422808,
  72   0.057493156217619066481721689402056,
  73   0.065974229882180495128128515115962,
  74   0.073755974737705206268243850022191,
  75   0.080755895229420215354694938460530,
  76   0.086899787201082979802387530715126,
  77   0.092122522237786128717632707087619,
  78   0.096368737174644259639468626351810,
  79   0.099593420586795267062780282103569,
  80   0.101762389748405504596428952168554,
  81   0.102852652893558840341285636705415
  82 };
  83
  84 static const double wgk[31] =   /* weights of the 61-point kronrod rule */
  85 {
  86   0.001389013698677007624551591226760,
  87   0.003890461127099884051267201844516,
  88   0.006630703915931292173319826369750,
  89   0.009273279659517763428441146892024,
  90   0.011823015253496341742232898853251,
  91   0.014369729507045804812451432443580,
  92   0.016920889189053272627572289420322,
  93   0.019414141193942381173408951050128,
  94   0.021828035821609192297167485738339,
  95   0.024191162078080601365686370725232,
  96   0.026509954882333101610601709335075,
  97   0.028754048765041292843978785354334,
  98   0.030907257562387762472884252943092,
  99   0.032981447057483726031814191016854,
 100   0.034979338028060024137499670731468,
 101   0.036882364651821229223911065617136,
 102   0.038678945624727592950348651532281,
 103   0.040374538951535959111995279752468,
 104   0.041969810215164246147147541285970,
 105   0.043452539701356069316831728117073,
 106   0.044814800133162663192355551616723,
 107   0.046059238271006988116271735559374,
 108   0.047185546569299153945261478181099,
 109   0.048185861757087129140779492298305,
 110   0.049055434555029778887528165367238,
 111   0.049795683427074206357811569379942,
 112   0.050405921402782346840893085653585,
 113   0.050881795898749606492297473049805,
 114   0.051221547849258772170656282604944,
 115   0.051426128537459025933862879215781,
 116   0.051494729429451567558340433647099
 117 };
 118
 119 void
 120 gsl_integration_qk61 (const gsl_function * f, double a, double b,
 121                       double *result, double *abserr,
 122                       double *resabs, double *resasc)
 123 {
 124   double fv1[31], fv2[31];
 125   gsl_integration_qk (31, xgk, wg, wgk, fv1, fv2, f, a, b, result, abserr, resabs, resasc);
 126 }