htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/mathieu_coeff.c

   1 /* specfunc/mathieu_coeff.c
   2  *
   3  * Copyright (C) 2002 Lowell Johnson
   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., 675 Mass Ave, Cambridge, MA 02139, USA.
  18  */
  19
  20 /* Author:  L. Johnson */
  21
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <math.h>
  25 #include <gsl/gsl_sf_mathieu.h>
  26
  27
  28 /*****************************************************************************
  29  * backward_recurse
  30  *
  31  * Purpose:
  32  ****************************************************************************/
  33 static void backward_recurse_c(double aa, double qq, double xx, double *ff,
  34                                double *gx, int even_odd, int ni)
  35 {
  36   int ii, nn;
  37   double g1;
  38
  39
  40   g1 = *gx;
  41   ff[ni] = xx;
  42
  43   if (even_odd == 0)
  44   {
  45       for (ii=0; ii<ni; ii++)
  46       {
  47           nn = GSL_SF_MATHIEU_COEFF - ii - 1;
  48           ff[ni-ii-1] = -1.0/((4*nn*nn - aa)/qq + ff[ni-ii]);
  49       }
  50       if (ni == GSL_SF_MATHIEU_COEFF - 1)
  51           ff[0] *= 2.0;
  52   }
  53   else
  54   {
  55       for (ii=0; ii<ni; ii++)
  56       {
  57           nn = GSL_SF_MATHIEU_COEFF - ii - 1;
  58           ff[ni-ii-1] = -1.0/(((2*nn + 1)*(2*nn + 1) - aa)/qq + ff[ni-ii]);
  59       }
  60   }
  61
  62   *gx = ff[0] - g1;
  63 }
  64
  65
  66 static void backward_recurse_s(double aa, double qq, double xx, double *ff,
  67                                double *gx, int even_odd, int ni)
  68 {
  69   int ii, nn;
  70   double g1;
  71
  72
  73   g1 = *gx;
  74   ff[ni] = xx;
  75
  76   if (even_odd == 0)
  77   {
  78       for (ii=0; ii<ni; ii++)
  79       {
  80           nn = GSL_SF_MATHIEU_COEFF - ii - 1;
  81           ff[ni-ii-1] = -1.0/((4*(nn + 1)*(nn + 1) - aa)/qq + ff[ni-ii]);
  82       }
  83   }
  84   else
  85   {
  86       for (ii=0; ii<ni; ii++)
  87       {
  88           nn = GSL_SF_MATHIEU_COEFF - ii - 1;
  89           ff[ni-ii-1] = -1.0/(((2*nn + 1)*(2*nn + 1) - aa)/qq + ff[ni-ii]);
  90       }
  91   }
  92
  93   *gx = ff[0] - g1;
  94 }
  95
  96
  97 int gsl_sf_mathieu_a_coeff(int order, double qq, double aa, double coeff[])
  98 {
  99   int ni, nn, ii, even_odd;
 100   double eps, g1, g2, x1, x2, e1, e2, de, xh, sum, ratio,
 101          ff[GSL_SF_MATHIEU_COEFF];
 102
 103
 104   eps = 1e-14;
 105   coeff[0] = 1.0;
 106
 107   even_odd = 0;
 108   if (order % 2 != 0)
 109       even_odd = 1;
 110
 111   /* If the coefficient array is not large enough to hold all necessary
 112      coefficients, error out. */
 113   if (order > GSL_SF_MATHIEU_COEFF)
 114       return GSL_FAILURE;
 115
 116   /* Handle the trivial case where q = 0. */
 117   if (qq == 0.0)
 118   {
 119       for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
 120           coeff[ii] = 0.0;
 121
 122       coeff[order/2] = 1.0;
 123
 124       return GSL_SUCCESS;
 125   }
 126
 127   if (order < 5)
 128   {
 129       nn = 0;
 130       sum = 0.0;
 131       if (even_odd == 0)
 132           ratio = aa/qq;
 133       else
 134           ratio = (aa - 1 - qq)/qq;
 135   }
 136   else
 137   {
 138       if (even_odd == 0)
 139       {
 140           coeff[1] = aa/qq;
 141           coeff[2] = (aa - 4)/qq*coeff[1] - 2;
 142           sum = coeff[0] + coeff[1] + coeff[2];
 143           for (ii=3; ii<order/2+1; ii++)
 144           {
 145               coeff[ii] = (aa - 4*(ii - 1)*(ii - 1))/qq*coeff[ii-1] -
 146                                                                   coeff[ii-2];
 147               sum += coeff[ii];
 148           }
 149       }
 150       else
 151       {
 152           coeff[1] = (aa - 1)/qq - 1;
 153           sum = coeff[0] + coeff[1];
 154           for (ii=2; ii<order/2+1; ii++)
 155           {
 156               coeff[ii] = (aa - (2*ii - 1)*(2*ii - 1))/qq*coeff[ii-1] -
 157                                                                   coeff[ii-2];
 158               sum += coeff[ii];
 159           }
 160       }
 161
 162       nn = ii - 1;
 163
 164       ratio = coeff[nn]/coeff[nn-1];
 165   }
 166
 167   ni = GSL_SF_MATHIEU_COEFF - nn - 1;
 168
 169   /* Compute first two points to start root-finding. */
 170   if (even_odd == 0)
 171       x1 = -qq/(4.0*GSL_SF_MATHIEU_COEFF*GSL_SF_MATHIEU_COEFF);
 172   else
 173       x1 = -qq/((2.0*GSL_SF_MATHIEU_COEFF + 1.0)*(2.0*GSL_SF_MATHIEU_COEFF + 1.0));
 174   g1 = ratio;
 175   backward_recurse_c(aa, qq, x1, ff, &g1, even_odd, ni);
 176   x2 = g1;
 177   g2 = ratio;
 178   backward_recurse_c(aa, qq, x2, ff, &g2, even_odd, ni);
 179
 180   /* Find the root. */
 181   while (1)
 182   {
 183       /* Compute the relative error. */
 184       e1 = g1 - x1;
 185       e2 = g2 - x2;
 186       de = e1 - e2;
 187
 188       /* If we are close enough to the root, break... */
 189       if (fabs(de) < eps)
 190           break;
 191
 192       /* Otherwise, determine the next guess and try again. */
 193       xh = (e1*x2 - e2*x1)/de;
 194       x1 = x2;
 195       g1 = g2;
 196       x2 = xh;
 197       g2 = ratio;
 198       backward_recurse_c(aa, qq, x2, ff, &g2, even_odd, ni);
 199   }
 200
 201   /* Compute the rest of the coefficients. */
 202   sum += coeff[nn];
 203   for (ii=nn+1; ii<GSL_SF_MATHIEU_COEFF; ii++)
 204   {
 205       coeff[ii] = ff[ii-nn-1]*coeff[ii-1];
 206       sum += coeff[ii];
 207
 208       /* If the coefficients are getting really small, set the remainder
 209          to zero. */
 210       if (fabs(coeff[ii]) < 1e-20)
 211       {
 212           for (; ii<GSL_SF_MATHIEU_COEFF;)
 213               coeff[ii++] = 0.0;
 214       }
 215   }
 216
 217   /* Normalize the coefficients. */
 218   for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
 219       coeff[ii] /= sum;
 220
 221   return GSL_SUCCESS;
 222 }
 223
 224
 225 int gsl_sf_mathieu_b_coeff(int order, double qq, double aa, double coeff[])
 226 {
 227   int ni, nn, ii, even_odd;
 228   double eps, g1, g2, x1, x2, e1, e2, de, xh, sum, ratio,
 229          ff[GSL_SF_MATHIEU_COEFF];
 230
 231
 232   eps = 1e-10;
 233   coeff[0] = 1.0;
 234
 235   even_odd = 0;
 236   if (order % 2 != 0)
 237       even_odd = 1;
 238
 239   /* If the coefficient array is not large enough to hold all necessary
 240      coefficients, error out. */
 241   if (order > GSL_SF_MATHIEU_COEFF)
 242       return GSL_FAILURE;
 243
 244   /* Handle the trivial case where q = 0. */
 245   if (qq == 0.0)
 246   {
 247       for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
 248           coeff[ii] = 0.0;
 249
 250       coeff[(order-1)/2] = 1.0;
 251
 252       return GSL_SUCCESS;
 253   }
 254
 255   if (order < 5)
 256   {
 257       nn = 0;
 258       sum = 0.0;
 259       if (even_odd == 0)
 260           ratio = (aa - 4)/qq;
 261       else
 262           ratio = (aa - 1 - qq)/qq;
 263   }
 264   else
 265   {
 266       if (even_odd == 0)
 267       {
 268           coeff[1] = (aa - 4)/qq;
 269           sum = 2*coeff[0] + 4*coeff[1];
 270           for (ii=2; ii<order/2; ii++)
 271           {
 272               coeff[ii] = (aa - 4*ii*ii)/qq*coeff[ii-1] - coeff[ii-2];
 273               sum += 2*(ii + 1)*coeff[ii];
 274           }
 275       }
 276       else
 277       {
 278           coeff[1] = (aa - 1)/qq + 1;
 279           sum = coeff[0] + 3*coeff[1];
 280           for (ii=2; ii<order/2+1; ii++)
 281           {
 282               coeff[ii] = (aa - (2*ii - 1)*(2*ii - 1))/qq*coeff[ii-1] -
 283                                                                   coeff[ii-2];
 284               sum += (2*(ii + 1) - 1)*coeff[ii];
 285           }
 286       }
 287
 288       nn = ii - 1;
 289
 290       ratio = coeff[nn]/coeff[nn-1];
 291   }
 292
 293   ni = GSL_SF_MATHIEU_COEFF - nn - 1;
 294
 295   /* Compute first two points to start root-finding. */
 296   if (even_odd == 0)
 297       x1 = -qq/(4.0*(GSL_SF_MATHIEU_COEFF + 1.0)*(GSL_SF_MATHIEU_COEFF + 1.0));
 298   else
 299       x1 = -qq/((2.0*GSL_SF_MATHIEU_COEFF + 1.0)*(2.0*GSL_SF_MATHIEU_COEFF + 1.0));
 300   g1 = ratio;
 301   backward_recurse_s(aa, qq, x1, ff, &g1, even_odd, ni);
 302   x2 = g1;
 303   g2 = ratio;
 304   backward_recurse_s(aa, qq, x2, ff, &g2, even_odd, ni);
 305
 306   /* Find the root. */
 307   while (1)
 308   {
 309       /* Compute the relative error. */
 310       e1 = g1 - x1;
 311       e2 = g2 - x2;
 312       de = e1 - e2;
 313
 314       /* If we are close enough to the root, break... */
 315       if (fabs(de) < eps)
 316           break;
 317
 318       /* Otherwise, determine the next guess and try again. */
 319       xh = (e1*x2 - e2*x1)/de;
 320       x1 = x2;
 321       g1 = g2;
 322       x2 = xh;
 323       g2 = ratio;
 324       backward_recurse_s(aa, qq, x2, ff, &g2, even_odd, ni);
 325   }
 326
 327   /* Compute the rest of the coefficients. */
 328   sum += 2*(nn + 1)*coeff[nn];
 329   for (ii=nn+1; ii<GSL_SF_MATHIEU_COEFF; ii++)
 330   {
 331       coeff[ii] = ff[ii-nn-1]*coeff[ii-1];
 332       sum += 2*(ii + 1)*coeff[ii];
 333
 334       /* If the coefficients are getting really small, set the remainder
 335          to zero. */
 336       if (fabs(coeff[ii]) < 1e-20)
 337       {
 338           for (; ii<GSL_SF_MATHIEU_COEFF;)
 339               coeff[ii++] = 0.0;
 340       }
 341   }
 342
 343   /* Normalize the coefficients. */
 344   for (ii=0; ii<GSL_SF_MATHIEU_COEFF; ii++)
 345       coeff[ii] /= sum;
 346
 347   return GSL_SUCCESS;
 348 }