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

   1 /* integration/qmomo.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 initialise (double * ri, double * rj, double * rg, double * rh,
  27             double alpha, double beta);
  28
  29 gsl_integration_qaws_table *
  30 gsl_integration_qaws_table_alloc (double alpha, double beta, int mu, int nu)
  31 {
  32   gsl_integration_qaws_table * t;
  33
  34   if (alpha < -1.0)
  35     {
  36       GSL_ERROR_VAL ("alpha must be greater than -1.0", GSL_EINVAL, 0);
  37     }
  38
  39   if (beta < -1.0)
  40     {
  41       GSL_ERROR_VAL ("beta must be greater than -1.0", GSL_EINVAL, 0);
  42     }
  43
  44   if (mu != 0 && mu != 1)
  45     {
  46       GSL_ERROR_VAL ("mu must be 0 or 1", GSL_EINVAL, 0);
  47     }
  48
  49   if (nu != 0 && nu != 1)
  50     {
  51       GSL_ERROR_VAL ("nu must be 0 or 1", GSL_EINVAL, 0);
  52     }
  53
  54   t = (gsl_integration_qaws_table *)
  55     malloc(sizeof(gsl_integration_qaws_table));
  56
  57   if (t == 0)
  58     {
  59       GSL_ERROR_VAL ("failed to allocate space for qaws_table struct",
  60                         GSL_ENOMEM, 0);
  61     }
  62
  63   t->alpha = alpha;
  64   t->beta = beta;
  65   t->mu = mu;
  66   t->nu = nu;
  67
  68   initialise (t->ri, t->rj, t->rg, t->rh, alpha, beta);
  69
  70   return t;
  71 }
  72
  73
  74 int
  75 gsl_integration_qaws_table_set (gsl_integration_qaws_table * t,
  76                                 double alpha, double beta, int mu, int nu)
  77 {
  78   if (alpha < -1.0)
  79     {
  80       GSL_ERROR ("alpha must be greater than -1.0", GSL_EINVAL);
  81     }
  82
  83   if (beta < -1.0)
  84     {
  85       GSL_ERROR ("beta must be greater than -1.0", GSL_EINVAL);
  86     }
  87
  88   if (mu != 0 && mu != 1)
  89     {
  90       GSL_ERROR ("mu must be 0 or 1", GSL_EINVAL);
  91     }
  92
  93   if (nu != 0 && nu != 1)
  94     {
  95       GSL_ERROR ("nu must be 0 or 1", GSL_EINVAL);
  96     }
  97
  98   t->alpha = alpha;
  99   t->beta = beta;
 100   t->mu = mu;
 101   t->nu = nu;
 102
 103   initialise (t->ri, t->rj, t->rg, t->rh, alpha, beta);
 104
 105   return GSL_SUCCESS;
 106 }
 107
 108
 109 void
 110 gsl_integration_qaws_table_free (gsl_integration_qaws_table * t)
 111 {
 112   free (t);
 113 }
 114
 115 static void
 116 initialise (double * ri, double * rj, double * rg, double * rh,
 117             double alpha, double beta)
 118 {
 119   const double alpha_p1 = alpha + 1.0;
 120   const double beta_p1 = beta + 1.0;
 121
 122   const double alpha_p2 = alpha + 2.0;
 123   const double beta_p2 = beta + 2.0;
 124
 125   const double r_alpha = pow (2.0, alpha_p1);
 126   const double r_beta = pow (2.0, beta_p1);
 127
 128   size_t i;
 129
 130   double an, anm1;
 131
 132   ri[0] = r_alpha / alpha_p1;
 133   ri[1] = ri[0] * alpha / alpha_p2;
 134
 135   an = 2.0;
 136   anm1 = 1.0;
 137
 138   for (i = 2; i < 25; i++)
 139     {
 140       ri[i] = -(r_alpha + an * (an - alpha_p2) * ri[i - 1])
 141         / (anm1 * (an + alpha_p1));
 142       anm1 = an;
 143       an = an + 1.0;
 144     }
 145
 146   rj[0] = r_beta / beta_p1;
 147   rj[1] = rj[0] * beta / beta_p2;
 148
 149   an = 2.0;
 150   anm1 = 1.0;
 151
 152   for (i = 2; i < 25; i++)
 153     {
 154       rj[i] = -(r_beta + an * (an - beta_p2) * rj[i - 1])
 155         / (anm1 * (an + beta_p1));
 156       anm1 = an;
 157       an = an + 1.0;
 158     }
 159
 160   rg[0] = -ri[0] / alpha_p1;
 161   rg[1] = -rg[0] - 2.0 * r_alpha / (alpha_p2 * alpha_p2);
 162
 163   an = 2.0;
 164   anm1 = 1.0;
 165
 166   for (i = 2; i < 25; i++)
 167     {
 168       rg[i] = -(an * (an - alpha_p2) * rg[i - 1] - an * ri[i - 1]
 169                 + anm1 * ri[i]) / (anm1 * (an + alpha_p1));
 170       anm1 = an;
 171       an = an + 1.0;
 172     }
 173
 174   rh[0] = -rj[0] / beta_p1;
 175   rh[1] = -rh[0] - 2.0 * r_beta / (beta_p2 * beta_p2);
 176
 177   an = 2.0;
 178   anm1 = 1.0;
 179
 180   for (i = 2; i < 25; i++)
 181     {
 182       rh[i] = -(an * (an - beta_p2) * rh[i - 1] - an * rj[i - 1]
 183                 + anm1 * rj[i]) / (anm1 * (an + beta_p1));
 184       anm1 = an;
 185       an = an + 1.0;
 186     }
 187
 188   for (i = 1; i < 25; i += 2)
 189     {
 190       rj[i] *= -1;
 191       rh[i] *= -1;
 192     }
 193 }