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

   1 /* integration/qcheb.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 <stdlib.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_integration.h>
  24
  25
  26 /* This function computes the 12-th order and 24-th order Chebyshev
  27    approximations to f(x) on [a,b] */
  28
  29 void
  30 gsl_integration_qcheb (gsl_function * f, double a, double b, double *cheb12, double *cheb24)
  31 {
  32   size_t i;
  33   double fval[25], v[12];
  34
  35   /* These are the values of cos(pi*k/24) for k=1..11 needed for the
  36      Chebyshev expansion of f(x) */
  37
  38   const double x[11] = { 0.9914448613738104,
  39                          0.9659258262890683,
  40                          0.9238795325112868,
  41                          0.8660254037844386,
  42                          0.7933533402912352,
  43                          0.7071067811865475,
  44                          0.6087614290087206,
  45                          0.5000000000000000,
  46                          0.3826834323650898,
  47                          0.2588190451025208,
  48                          0.1305261922200516 };
  49
  50   const double center = 0.5 * (b + a);
  51   const double half_length =  0.5 * (b - a);
  52
  53   fval[0] = 0.5 * GSL_FN_EVAL (f, b);
  54   fval[12] = GSL_FN_EVAL (f, center);
  55   fval[24] = 0.5 * GSL_FN_EVAL (f, a);
  56
  57   for (i = 1; i < 12; i++)
  58     {
  59       const size_t j = 24 - i;
  60       const double u = half_length * x[i-1];
  61       fval[i] = GSL_FN_EVAL(f, center + u);
  62       fval[j] = GSL_FN_EVAL(f, center - u);
  63     }
  64
  65   for (i = 0; i < 12; i++)
  66     {
  67       const size_t j = 24 - i;
  68       v[i] = fval[i] - fval[j];
  69       fval[i] = fval[i] + fval[j];
  70     }
  71
  72   {
  73     const double alam1 = v[0] - v[8];
  74     const double alam2 = x[5] * (v[2] - v[6] - v[10]);
  75
  76     cheb12[3] = alam1 + alam2;
  77     cheb12[9] = alam1 - alam2;
  78   }
  79
  80   {
  81     const double alam1 = v[1] - v[7] - v[9];
  82     const double alam2 = v[3] - v[5] - v[11];
  83     {
  84       const double alam = x[2] * alam1 + x[8] * alam2;
  85
  86       cheb24[3] = cheb12[3] + alam;
  87       cheb24[21] = cheb12[3] - alam;
  88     }
  89
  90     {
  91       const double alam = x[8] * alam1 - x[2] * alam2;
  92       cheb24[9] = cheb12[9] + alam;
  93       cheb24[15] = cheb12[9] - alam;
  94     }
  95   }
  96
  97   {
  98     const double part1 = x[3] * v[4];
  99     const double part2 = x[7] * v[8];
 100     const double part3 = x[5] * v[6];
 101
 102     {
 103       const double alam1 = v[0] + part1 + part2;
 104       const double alam2 = x[1] * v[2] + part3 + x[9] * v[10];
 105
 106       cheb12[1] = alam1 + alam2;
 107       cheb12[11] = alam1 - alam2;
 108     }
 109
 110     {
 111       const double alam1 = v[0] - part1 + part2;
 112       const double alam2 = x[9] * v[2] - part3 + x[1] * v[10];
 113       cheb12[5] = alam1 + alam2;
 114       cheb12[7] = alam1 - alam2;
 115     }
 116   }
 117
 118   {
 119     const double alam = (x[0] * v[1] + x[2] * v[3] + x[4] * v[5]
 120                 + x[6] * v[7] + x[8] * v[9] + x[10] * v[11]);
 121     cheb24[1] = cheb12[1] + alam;
 122     cheb24[23] = cheb12[1] - alam;
 123   }
 124
 125   {
 126     const double alam = (x[10] * v[1] - x[8] * v[3] + x[6] * v[5]
 127                 - x[4] * v[7] + x[2] * v[9] - x[0] * v[11]);
 128     cheb24[11] = cheb12[11] + alam;
 129     cheb24[13] = cheb12[11] - alam;
 130   }
 131
 132   {
 133     const double alam = (x[4] * v[1] - x[8] * v[3] - x[0] * v[5]
 134                 - x[10] * v[7] + x[2] * v[9] + x[6] * v[11]);
 135     cheb24[5] = cheb12[5] + alam;
 136     cheb24[19] = cheb12[5] - alam;
 137   }
 138
 139   {
 140     const double alam = (x[6] * v[1] - x[2] * v[3] - x[10] * v[5]
 141                 + x[0] * v[7] - x[8] * v[9] - x[4] * v[11]);
 142     cheb24[7] = cheb12[7] + alam;
 143     cheb24[17] = cheb12[7] - alam;
 144   }
 145
 146   for (i = 0; i < 6; i++)
 147     {
 148       const size_t j = 12 - i;
 149       v[i] = fval[i] - fval[j];
 150       fval[i] = fval[i] + fval[j];
 151     }
 152
 153   {
 154     const double alam1 = v[0] + x[7] * v[4];
 155     const double alam2 = x[3] * v[2];
 156
 157     cheb12[2] = alam1 + alam2;
 158     cheb12[10] = alam1 - alam2;
 159   }
 160
 161   cheb12[6] = v[0] - v[4];
 162
 163   {
 164     const double alam = x[1] * v[1] + x[5] * v[3] + x[9] * v[5];
 165     cheb24[2] = cheb12[2] + alam;
 166     cheb24[22] = cheb12[2] - alam;
 167   }
 168
 169   {
 170     const double alam = x[5] * (v[1] - v[3] - v[5]);
 171     cheb24[6] = cheb12[6] + alam;
 172     cheb24[18] = cheb12[6] - alam;
 173   }
 174
 175   {
 176     const double alam = x[9] * v[1] - x[5] * v[3] + x[1] * v[5];
 177     cheb24[10] = cheb12[10] + alam;
 178     cheb24[14] = cheb12[10] - alam;
 179   }
 180
 181   for (i = 0; i < 3; i++)
 182     {
 183       const size_t j = 6 - i;
 184       v[i] = fval[i] - fval[j];
 185       fval[i] = fval[i] + fval[j];
 186     }
 187
 188   cheb12[4] = v[0] + x[7] * v[2];
 189   cheb12[8] = fval[0] - x[7] * fval[2];
 190
 191   {
 192     const double alam = x[3] * v[1];
 193     cheb24[4] = cheb12[4] + alam;
 194     cheb24[20] = cheb12[4] - alam;
 195   }
 196
 197   {
 198     const double alam = x[7] * fval[1] - fval[3];
 199     cheb24[8] = cheb12[8] + alam;
 200     cheb24[16] = cheb12[8] - alam;
 201   }
 202
 203   cheb12[0] = fval[0] + fval[2];
 204
 205   {
 206     const double alam = fval[1] + fval[3];
 207     cheb24[0] = cheb12[0] + alam;
 208     cheb24[24] = cheb12[0] - alam;
 209   }
 210
 211   cheb12[12] = v[0] - v[2];
 212   cheb24[12] = cheb12[12];
 213
 214   for (i = 1; i < 12; i++)
 215     {
 216       cheb12[i] *= 1.0 / 6.0;
 217     }
 218
 219   cheb12[0] *= 1.0 / 12.0;
 220   cheb12[12] *= 1.0 / 12.0;
 221
 222   for (i = 1; i < 24; i++)
 223     {
 224       cheb24[i] *= 1.0 / 12.0;
 225     }
 226
 227   cheb24[0] *= 1.0 / 24.0;
 228   cheb24[24] *= 1.0 / 24.0;
 229 }