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

   1 /* integration/qmomof.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_integration.h>
  23 #include <gsl/gsl_errno.h>
  24
  25 static void
  26 compute_moments (double par, double * cheb);
  27
  28 static int
  29 dgtsl (size_t n, double *c, double *d, double *e, double *b);
  30
  31 gsl_integration_qawo_table *
  32 gsl_integration_qawo_table_alloc (double omega, double L,
  33                                   enum gsl_integration_qawo_enum sine,
  34                                   size_t n)
  35 {
  36   gsl_integration_qawo_table *t;
  37   double * chebmo;
  38
  39   if (n == 0)
  40     {
  41       GSL_ERROR_VAL ("table length n must be positive integer",
  42                         GSL_EDOM, 0);
  43     }
  44
  45   t = (gsl_integration_qawo_table *)
  46     malloc (sizeof (gsl_integration_qawo_table));
  47
  48   if (t == 0)
  49     {
  50       GSL_ERROR_VAL ("failed to allocate space for qawo_table struct",
  51                         GSL_ENOMEM, 0);
  52     }
  53
  54   chebmo = (double *)  malloc (25 * n * sizeof (double));
  55
  56   if (chebmo == 0)
  57     {
  58       free (t);
  59       GSL_ERROR_VAL ("failed to allocate space for chebmo block",
  60                         GSL_ENOMEM, 0);
  61     }
  62
  63   t->n = n;
  64   t->sine = sine;
  65   t->omega = omega;
  66   t->L = L;
  67   t->par = 0.5 * omega * L;
  68   t->chebmo = chebmo;
  69
  70   /* precompute the moments */
  71
  72   {
  73     size_t i;
  74     double scale = 1.0;
  75
  76     for (i = 0 ; i < t->n; i++)
  77       {
  78         compute_moments (t->par * scale, t->chebmo + 25*i);
  79         scale *= 0.5;
  80       }
  81   }
  82
  83   return t;
  84 }
  85
  86 int
  87 gsl_integration_qawo_table_set (gsl_integration_qawo_table * t,
  88                                     double omega, double L,
  89                                     enum gsl_integration_qawo_enum sine)
  90 {
  91   t->omega = omega;
  92   t->sine = sine;
  93   t->L = L;
  94   t->par = 0.5 * omega * L;
  95
  96   /* recompute the moments */
  97
  98   {
  99     size_t i;
 100     double scale = 1.0;
 101
 102     for (i = 0 ; i < t->n; i++)
 103       {
 104         compute_moments (t->par * scale, t->chebmo + 25*i);
 105         scale *= 0.5;
 106       }
 107   }
 108
 109   return GSL_SUCCESS;
 110 }
 111
 112
 113 int
 114 gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * t,
 115                                            double L)
 116 {
 117   /* return immediately if the length is the same as the old length */
 118
 119   if (L == t->L)
 120     return GSL_SUCCESS;
 121
 122   /* otherwise reset the table and compute the new parameters */
 123
 124   t->L = L;
 125   t->par = 0.5 * t->omega * L;
 126
 127   /* recompute the moments */
 128
 129   {
 130     size_t i;
 131     double scale = 1.0;
 132
 133     for (i = 0 ; i < t->n; i++)
 134       {
 135         compute_moments (t->par * scale, t->chebmo + 25*i);
 136         scale *= 0.5;
 137       }
 138   }
 139
 140   return GSL_SUCCESS;
 141 }
 142
 143
 144 void
 145 gsl_integration_qawo_table_free (gsl_integration_qawo_table * t)
 146 {
 147   free (t->chebmo);
 148   free (t);
 149 }
 150
 151 static void
 152 compute_moments (double par, double *chebmo)
 153 {
 154   double v[28], d[25], d1[25], d2[25];
 155
 156   const size_t noeq = 25;
 157
 158   const double par2 = par * par;
 159   const double par4 = par2 * par2;
 160   const double par22 = par2 + 2.0;
 161
 162   const double sinpar = sin (par);
 163   const double cospar = cos (par);
 164
 165   size_t i;
 166
 167   /* compute the chebyschev moments with respect to cosine */
 168
 169   double ac = 8 * cospar;
 170   double as = 24 * par * sinpar;
 171
 172   v[0] = 2 * sinpar / par;
 173   v[1] = (8 * cospar + (2 * par2 - 8) * sinpar / par) / par2;
 174   v[2] = (32 * (par2 - 12) * cospar
 175           + (2 * ((par2 - 80) * par2 + 192) * sinpar) / par) / par4;
 176
 177   if (fabs (par) <= 24)
 178     {
 179       /* compute the moments as the solution of a boundary value
 180          problem using the asyptotic expansion as an endpoint */
 181
 182       double an2, ass, asap;
 183       double an = 6;
 184       size_t k;
 185
 186       for (k = 0; k < noeq - 1; k++)
 187         {
 188           an2 = an * an;
 189           d[k] = -2 * (an2 - 4) * (par22 - 2 * an2);
 190           d2[k] = (an - 1) * (an - 2) * par2;
 191           d1[k + 1] = (an + 3) * (an + 4) * par2;
 192           v[k + 3] = as - (an2 - 4) * ac;
 193           an = an + 2.0;
 194         }
 195
 196       an2 = an * an;
 197
 198       d[noeq - 1] = -2 * (an2 - 4) * (par22 - 2 * an2);
 199       v[noeq + 2] = as - (an2 - 4) * ac;
 200       v[3] = v[3] - 56 * par2 * v[2];
 201
 202       ass = par * sinpar;
 203       asap = (((((210 * par2 - 1) * cospar - (105 * par2 - 63) * ass) / an2
 204                 - (1 - 15 * par2) * cospar + 15 * ass) / an2
 205                - cospar + 3 * ass) / an2
 206               - cospar) / an2;
 207       v[noeq + 2] = v[noeq + 2] - 2 * asap * par2 * (an - 1) * (an - 2);
 208
 209       dgtsl (noeq, d1, d, d2, v + 3);
 210
 211     }
 212   else
 213     {
 214       /* compute the moments by forward recursion */
 215       size_t k;
 216       double an = 4;
 217
 218       for (k = 3; k < 13; k++)
 219         {
 220           double an2 = an * an;
 221           v[k] = ((an2 - 4) * (2 * (par22 - 2 * an2) * v[k - 1] - ac)
 222                   + as - par2 * (an + 1) * (an + 2) * v[k - 2])
 223             / (par2 * (an - 1) * (an - 2));
 224           an = an + 2.0;
 225         }
 226     }
 227
 228
 229   for (i = 0; i < 13; i++)
 230     {
 231       chebmo[2 * i] = v[i];
 232     }
 233
 234   /* compute the chebyschev moments with respect to sine */
 235
 236   v[0] = 2 * (sinpar - par * cospar) / par2;
 237   v[1] = (18 - 48 / par2) * sinpar / par2 + (-2 + 48 / par2) * cospar / par;
 238
 239   ac = -24 * par * cospar;
 240   as = -8 * sinpar;
 241
 242   if (fabs (par) <= 24)
 243     {
 244       /* compute the moments as the solution of a boundary value
 245          problem using the asyptotic expansion as an endpoint */
 246
 247       size_t k;
 248       double an2, ass, asap;
 249       double an = 5;
 250
 251       for (k = 0; k < noeq - 1; k++)
 252         {
 253           an2 = an * an;
 254           d[k] = -2 * (an2 - 4) * (par22 - 2 * an2);
 255           d2[k] = (an - 1) * (an - 2) * par2;
 256           d1[k + 1] = (an + 3) * (an + 4) * par2;
 257           v[k + 2] = ac + (an2 - 4) * as;
 258           an = an + 2.0;
 259         }
 260
 261       an2 = an * an;
 262
 263       d[noeq - 1] = -2 * (an2 - 4) * (par22 - 2 * an2);
 264       v[noeq + 1] = ac + (an2 - 4) * as;
 265       v[2] = v[2] - 42 * par2 * v[1];
 266
 267       ass = par * cospar;
 268       asap = (((((105 * par2 - 63) * ass - (210 * par2 - 1) * sinpar) / an2
 269                 + (15 * par2 - 1) * sinpar
 270                 - 15 * ass) / an2 - sinpar - 3 * ass) / an2 - sinpar) / an2;
 271       v[noeq + 1] = v[noeq + 1] - 2 * asap * par2 * (an - 1) * (an - 2);
 272
 273       dgtsl (noeq, d1, d, d2, v + 2);
 274
 275     }
 276   else
 277     {
 278       /* compute the moments by forward recursion */
 279       size_t k;
 280       double an = 3;
 281       for (k = 2; k < 12; k++)
 282         {
 283           double an2 = an * an;
 284           v[k] = ((an2 - 4) * (2 * (par22 - 2 * an2) * v[k - 1] + as)
 285                   + ac - par2 * (an + 1) * (an + 2) * v[k - 2])
 286             / (par2 * (an - 1) * (an - 2));
 287           an = an + 2.0;
 288         }
 289     }
 290
 291   for (i = 0; i < 12; i++)
 292     {
 293       chebmo[2 * i + 1] = v[i];
 294     }
 295
 296 }
 297
 298 static int
 299 dgtsl (size_t n, double *c, double *d, double *e, double *b)
 300 {
 301   /* solves a tridiagonal matrix A x = b
 302
 303      c[1 .. n - 1]   subdiagonal of the matrix A
 304      d[0 .. n - 1]   diagonal of the matrix A
 305      e[0 .. n - 2]   superdiagonal of the matrix A
 306
 307      b[0 .. n - 1]   right hand side, replaced by the solution vector x */
 308
 309   size_t k;
 310
 311   c[0] = d[0];
 312
 313   if (n == 0)
 314     {
 315       return GSL_SUCCESS;
 316     }
 317
 318   if (n == 1)
 319     {
 320       b[0] = b[0] / d[0] ;
 321       return GSL_SUCCESS;
 322     }
 323
 324   d[0] = e[0];
 325   e[0] = 0;
 326   e[n - 1] = 0;
 327
 328   for (k = 0; k < n - 1; k++)
 329     {
 330       size_t k1 = k + 1;
 331
 332       if (fabs (c[k1]) >= fabs (c[k]))
 333         {
 334           {
 335             double t = c[k1];
 336             c[k1] = c[k];
 337             c[k] = t;
 338           };
 339           {
 340             double t = d[k1];
 341             d[k1] = d[k];
 342             d[k] = t;
 343           };
 344           {
 345             double t = e[k1];
 346             e[k1] = e[k];
 347             e[k] = t;
 348           };
 349           {
 350             double t = b[k1];
 351             b[k1] = b[k];
 352             b[k] = t;
 353           };
 354         }
 355
 356       if (c[k] == 0)
 357         {
 358           return GSL_FAILURE ;
 359         }
 360
 361       {
 362         double t = -c[k1] / c[k];
 363
 364         c[k1] = d[k1] + t * d[k];
 365         d[k1] = e[k1] + t * e[k];
 366         e[k1] = 0;
 367         b[k1] = b[k1] + t * b[k];
 368       }
 369
 370     }
 371
 372   if (c[n - 1] == 0)
 373     {
 374       return GSL_FAILURE;
 375     }
 376
 377
 378   b[n - 1] = b[n - 1] / c[n - 1];
 379
 380   b[n - 2] = (b[n - 2] - d[n - 2] * b[n - 1]) / c[n - 2];
 381
 382   for (k = n ; k > 2; k--)
 383     {
 384       size_t kb = k - 3;
 385       b[kb] = (b[kb] - d[kb] * b[kb + 1] - e[kb] * b[kb + 2]) / c[kb];
 386     }
 387
 388   return GSL_SUCCESS;
 389 }