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

   1 /* specfunc/mathieu_charv.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 <stdio.h>
  25 #include <math.h>
  26 #include <gsl/gsl_math.h>
  27 #include <gsl/gsl_eigen.h>
  28 #include <gsl/gsl_errno.h>
  29 #include <gsl/gsl_sf_mathieu.h>
  30
  31
  32 /* prototypes */
  33 static double solve_cubic(double c2, double c1, double c0);
  34
  35
  36 static double ceer(int order, double qq, double aa, int nterms)
  37 {
  38
  39   double term, term1;
  40   int ii, n1;
  41
  42
  43   if (order == 0)
  44       term = 0.0;
  45   else
  46   {
  47       term = 2.0*qq*qq/aa;
  48
  49       if (order != 2)
  50       {
  51           n1 = order/2 - 1;
  52
  53           for (ii=0; ii<n1; ii++)
  54               term = qq*qq/(aa - 4.0*(ii+1)*(ii+1) - term);
  55       }
  56   }
  57
  58   term += order*order;
  59
  60   term1 = 0.0;
  61
  62   for (ii=0; ii<nterms; ii++)
  63       term1 = qq*qq/
  64         (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
  65
  66   if (order == 0)
  67       term1 *= 2.0;
  68
  69   return (term + term1 - aa);
  70 }
  71
  72
  73 static double ceor(int order, double qq, double aa, int nterms)
  74 {
  75   double term, term1;
  76   int ii, n1;
  77
  78   term = qq;
  79   n1 = (int)((float)order/2.0 - 0.5);
  80
  81   for (ii=0; ii<n1; ii++)
  82       term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);
  83   term += order*order;
  84
  85   term1 = 0.0;
  86   for (ii=0; ii<nterms; ii++)
  87       term1 = qq*qq/
  88         (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
  89
  90   return (term + term1 - aa);
  91 }
  92
  93
  94 static double seer(int order, double qq, double aa, int nterms)
  95 {
  96   double term, term1;
  97   int ii, n1;
  98
  99   term = 0.0;
 100   n1 = order/2 - 1;
 101
 102   for (ii=0; ii<n1; ii++)
 103       term = qq*qq/(aa - 4.0*(ii + 1)*(ii + 1) - term);
 104   term += order*order;
 105
 106   term1 = 0.0;
 107   for (ii=0; ii<nterms; ii++)
 108       term1 = qq*qq/
 109         (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
 110
 111   return (term + term1 - aa);
 112 }
 113
 114
 115 static double seor(int order, double qq, double aa, int nterms)
 116 {
 117   double term, term1;
 118   int ii, n1;
 119
 120
 121   term = -1.0*qq;
 122   n1 = (int)((float)order/2.0 - 0.5);
 123   for (ii=0; ii<n1; ii++)
 124       term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);
 125   term += order*order;
 126
 127   term1 = 0.0;
 128   for (ii=0; ii<nterms; ii++)
 129       term1 = qq*qq/
 130         (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
 131
 132   return (term + term1 - aa);
 133 }
 134
 135
 136 /*----------------------------------------------------------------------------
 137  * Asymptotic and approximation routines for the characteristic value.
 138  *
 139  * Adapted from F.A. Alhargan's paper,
 140  * "Algorithms for the Computation of All Mathieu Functions of Integer
 141  * Orders," ACM Transactions on Mathematical Software, Vol. 26, No. 3,
 142  * September 2000, pp. 390-407.
 143  *--------------------------------------------------------------------------*/
 144 static double asymptotic(int order, double qq)
 145 {
 146   double asymp;
 147   double nn, n2, n4, n6;
 148   double hh, ah, ah2, ah3, ah4, ah5;
 149
 150
 151   /* Set up temporary variables to simplify the readability. */
 152   nn = 2*order + 1;
 153   n2 = nn*nn;
 154   n4 = n2*n2;
 155   n6 = n4*n2;
 156
 157   hh = 2*sqrt(qq);
 158   ah = 16*hh;
 159   ah2 = ah*ah;
 160   ah3 = ah2*ah;
 161   ah4 = ah3*ah;
 162   ah5 = ah4*ah;
 163
 164   /* Equation 38, p. 397 of Alhargan's paper. */
 165   asymp = -2*qq + nn*hh - 0.125*(n2 + 1);
 166   asymp -= 0.25*nn*(                          n2 +     3)/ah;
 167   asymp -= 0.25*   (             5*n4 +    34*n2 +     9)/ah2;
 168   asymp -= 0.25*nn*(            33*n4 +   410*n2 +   405)/ah3;
 169   asymp -=         ( 63*n6 +  1260*n4 +  2943*n2 +   486)/ah4;
 170   asymp -=      nn*(527*n6 + 15617*n4 + 69001*n2 + 41607)/ah5;
 171
 172   return asymp;
 173 }
 174
 175
 176 /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
 177 static double solve_cubic(double c2, double c1, double c0)
 178 {
 179   double qq, rr, ww, ss, tt;
 180
 181
 182   qq = (3*c1 - c2*c2)/9;
 183   rr = (9*c2*c1 - 27*c0 - 2*c2*c2*c2)/54;
 184   ww = qq*qq*qq + rr*rr;
 185
 186   if (ww >= 0)
 187   {
 188       double t1 = rr + sqrt(ww);
 189       ss = fabs(t1)/t1*pow(fabs(t1), 1/3.);
 190       t1 = rr - sqrt(ww);
 191       tt = fabs(t1)/t1*pow(fabs(t1), 1/3.);
 192   }
 193   else
 194   {
 195       double theta = acos(rr/sqrt(-qq*qq*qq));
 196       ss = 2*sqrt(-qq)*cos((theta + 4*M_PI)/3.);
 197       tt = 0.0;
 198   }
 199
 200   return (ss + tt - c2/3);
 201 }
 202
 203
 204 /* Compute an initial approximation for the characteristic value. */
 205 static double approx_c(int order, double qq)
 206 {
 207   double approx;
 208   double c0, c1, c2;
 209
 210
 211   if (order < 0)
 212   {
 213     GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);
 214   }
 215
 216   switch (order)
 217   {
 218       case 0:
 219           if (qq <= 4)
 220               return (2 - sqrt(4 + 2*qq*qq)); /* Eqn. 31 */
 221           else
 222               return asymptotic(order, qq);
 223           break;
 224
 225       case 1:
 226           if (qq <= 4)
 227               return (5 + 0.5*(qq - sqrt(5*qq*qq - 16*qq + 64))); /* Eqn. 32 */
 228           else
 229               return asymptotic(order, qq);
 230           break;
 231
 232       case 2:
 233           if (qq <= 3)
 234           {
 235               c2 = -8.0;  /* Eqn. 33 */
 236               c1 = -48 - 3*qq*qq;
 237               c0 = 20*qq*qq;
 238           }
 239           else
 240               return asymptotic(order, qq);
 241           break;
 242
 243       case 3:
 244           if (qq <= 6.25)
 245           {
 246               c2 = -qq - 8;  /* Eqn. 34 */
 247               c1 = 16*qq - 128 - 2*qq*qq;
 248               c0 = qq*qq*(qq + 8);
 249           }
 250           else
 251               return asymptotic(order, qq);
 252           break;
 253
 254       default:
 255           if (order < 70)
 256           {
 257               if (1.7*order > 2*sqrt(qq))
 258               {
 259                   /* Eqn. 30 */
 260                   double n2 = (double)(order*order);
 261                   double n22 = (double)((n2 - 1)*(n2 - 1));
 262                   double q2 = qq*qq;
 263                   double q4 = q2*q2;
 264                   approx = n2 + 0.5*q2/(n2 - 1);
 265                   approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));
 266                   approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/
 267                       (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));
 268                   if (1.4*order < 2*sqrt(qq))
 269                   {
 270                       approx += asymptotic(order, qq);
 271                       approx *= 0.5;
 272                   }
 273               }
 274               else
 275                   approx = asymptotic(order, qq);
 276
 277               return approx;
 278           }
 279           else
 280               return order*order;
 281   }
 282
 283   /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
 284   approx = solve_cubic(c2, c1, c0);
 285
 286   if ( approx < 0 && sqrt(qq) > 0.1*order )
 287       return asymptotic(order-1, qq);
 288   else
 289       return (order*order + fabs(approx));
 290 }
 291
 292
 293 static double approx_s(int order, double qq)
 294 {
 295   double approx;
 296   double c0, c1, c2;
 297
 298
 299   if (order < 1)
 300   {
 301     GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);
 302   }
 303
 304   switch (order)
 305   {
 306       case 1:
 307           if (qq <= 4)
 308               return (5 - 0.5*(qq + sqrt(5*qq*qq + 16*qq + 64))); /* Eqn. 35 */
 309           else
 310               return asymptotic(order-1, qq);
 311           break;
 312
 313       case 2:
 314           if (qq <= 5)
 315               return (10 - sqrt(36 + qq*qq)); /* Eqn. 36 */
 316           else
 317               return asymptotic(order-1, qq);
 318           break;
 319
 320       case 3:
 321           if (qq <= 6.25)
 322           {
 323               c2 = qq - 8; /* Eqn. 37 */
 324               c1 = -128 - 16*qq - 2*qq*qq;
 325               c0 = qq*qq*(8 - qq);
 326           }
 327           else
 328               return asymptotic(order-1, qq);
 329           break;
 330
 331       default:
 332           if (order < 70)
 333           {
 334               if (1.7*order > 2*sqrt(qq))
 335               {
 336                   /* Eqn. 30 */
 337                   double n2 = (double)(order*order);
 338                   double n22 = (double)((n2 - 1)*(n2 - 1));
 339                   double q2 = qq*qq;
 340                   double q4 = q2*q2;
 341                   approx = n2 + 0.5*q2/(n2 - 1);
 342                   approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));
 343                   approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/
 344                       (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));
 345                   if (1.4*order < 2*sqrt(qq))
 346                   {
 347                       approx += asymptotic(order-1, qq);
 348                       approx *= 0.5;
 349                   }
 350               }
 351               else
 352                   approx = asymptotic(order-1, qq);
 353
 354               return approx;
 355           }
 356           else
 357               return order*order;
 358   }
 359
 360   /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
 361   approx = solve_cubic(c2, c1, c0);
 362
 363   if ( approx < 0 && sqrt(qq) > 0.1*order )
 364       return asymptotic(order-1, qq);
 365   else
 366       return (order*order + fabs(approx));
 367 }
 368
 369
 370 int gsl_sf_mathieu_a(int order, double qq, gsl_sf_result *result)
 371 {
 372   int even_odd, nterms = 50, ii, counter = 0, maxcount = 200;
 373   double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;
 374
 375
 376   even_odd = 0;
 377   if (order % 2 != 0)
 378       even_odd = 1;
 379
 380   /* If the argument is 0, then the coefficient is simply the square of
 381      the order. */
 382   if (qq == 0)
 383   {
 384       result->val = order*order;
 385       result->err = 0.0;
 386       return GSL_SUCCESS;
 387   }
 388
 389   /* Use symmetry characteristics of the functions to handle cases with
 390      negative order and/or argument q.  See Abramowitz & Stegun, 20.8.3. */
 391   if (order < 0)
 392       order *= -1;
 393   if (qq < 0.0)
 394   {
 395       if (even_odd == 0)
 396           return gsl_sf_mathieu_a(order, -qq, result);
 397       else
 398           return gsl_sf_mathieu_b(order, -qq, result);
 399   }
 400
 401   /* Compute an initial approximation for the characteristic value. */
 402   aa = approx_c(order, qq);
 403
 404   /* Save the original approximation for later comparison. */
 405   aa_orig = aa;
 406
 407   /* Loop as long as the final value is not near the approximate value
 408      (with a max limit to avoid potential infinite loop). */
 409   while (counter < maxcount)
 410   {
 411       a1 = aa + 0.001;
 412       ii = 0;
 413       if (even_odd == 0)
 414           fa1 = ceer(order, qq, a1, nterms);
 415       else
 416           fa1 = ceor(order, qq, a1, nterms);
 417
 418       for (;;)
 419       {
 420           if (even_odd == 0)
 421               fa = ceer(order, qq, aa, nterms);
 422           else
 423               fa = ceor(order, qq, aa, nterms);
 424
 425           a2 = a1;
 426           a1 = aa;
 427
 428           if (fa == fa1)
 429           {
 430               result->err = GSL_DBL_EPSILON;
 431               break;
 432           }
 433           aa -= (aa - a2)/(fa - fa1)*fa;
 434           dela = fabs(aa - a2);
 435           if (dela < GSL_DBL_EPSILON)
 436           {
 437               result->err = GSL_DBL_EPSILON;
 438               break;
 439           }
 440           if (ii > 20)
 441           {
 442               result->err = dela;
 443               break;
 444           }
 445           fa1 = fa;
 446           ii++;
 447       }
 448
 449       /* If the solution found is not near the original approximation,
 450          tweak the approximate value, and try again. */
 451       if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)))
 452       {
 453           counter++;
 454           if (counter == maxcount)
 455           {
 456               result->err = fabs(aa - aa_orig);
 457               break;
 458           }
 459           if (aa > aa_orig)
 460               aa = aa_orig - da*counter;
 461           else
 462               aa = aa_orig + da*counter;
 463
 464           continue;
 465       }
 466       else
 467           break;
 468   }
 469
 470   result->val = aa;
 471
 472   /* If we went through the maximum number of retries and still didn't
 473      find the solution, let us know. */
 474   if (counter == maxcount)
 475   {
 476       GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);
 477   }
 478
 479   return GSL_SUCCESS;
 480 }
 481
 482
 483 int gsl_sf_mathieu_b(int order, double qq, gsl_sf_result *result)
 484 {
 485   int even_odd, nterms = 50, ii, counter = 0, maxcount = 200;
 486   double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;
 487
 488
 489   even_odd = 0;
 490   if (order % 2 != 0)
 491       even_odd = 1;
 492
 493   /* The order cannot be 0. */
 494   if (order == 0)
 495   {
 496       GSL_ERROR("Characteristic value undefined for order 0", GSL_EFAILED);
 497   }
 498
 499   /* If the argument is 0, then the coefficient is simply the square of
 500      the order. */
 501   if (qq == 0)
 502   {
 503       result->val = order*order;
 504       result->err = 0.0;
 505       return GSL_SUCCESS;
 506   }
 507
 508   /* Use symmetry characteristics of the functions to handle cases with
 509      negative order and/or argument q.  See Abramowitz & Stegun, 20.8.3. */
 510   if (order < 0)
 511       order *= -1;
 512   if (qq < 0.0)
 513   {
 514       if (even_odd == 0)
 515           return gsl_sf_mathieu_b(order, -qq, result);
 516       else
 517           return gsl_sf_mathieu_a(order, -qq, result);
 518   }
 519
 520   /* Compute an initial approximation for the characteristic value. */
 521   aa = approx_s(order, qq);
 522
 523   /* Save the original approximation for later comparison. */
 524   aa_orig = aa;
 525
 526   /* Loop as long as the final value is not near the approximate value
 527      (with a max limit to avoid potential infinite loop). */
 528   while (counter < maxcount)
 529   {
 530       a1 = aa + 0.001;
 531       ii = 0;
 532       if (even_odd == 0)
 533           fa1 = seer(order, qq, a1, nterms);
 534       else
 535           fa1 = seor(order, qq, a1, nterms);
 536
 537       for (;;)
 538       {
 539           if (even_odd == 0)
 540               fa = seer(order, qq, aa, nterms);
 541           else
 542               fa = seor(order, qq, aa, nterms);
 543
 544           a2 = a1;
 545           a1 = aa;
 546
 547           if (fa == fa1)
 548           {
 549               result->err = GSL_DBL_EPSILON;
 550               break;
 551           }
 552           aa -= (aa - a2)/(fa - fa1)*fa;
 553           dela = fabs(aa - a2);
 554           if (dela < 1e-18)
 555           {
 556               result->err = GSL_DBL_EPSILON;
 557               break;
 558           }
 559           if (ii > 20)
 560           {
 561               result->err = dela;
 562               break;
 563           }
 564           fa1 = fa;
 565           ii++;
 566       }
 567
 568       /* If the solution found is not near the original approximation,
 569          tweak the approximate value, and try again. */
 570       if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)))
 571       {
 572           counter++;
 573           if (counter == maxcount)
 574           {
 575               result->err = fabs(aa - aa_orig);
 576               break;
 577           }
 578           if (aa > aa_orig)
 579               aa = aa_orig - da*counter;
 580           else
 581               aa = aa_orig + da*counter;
 582
 583           continue;
 584       }
 585       else
 586           break;
 587   }
 588
 589   result->val = aa;
 590
 591   /* If we went through the maximum number of retries and still didn't
 592      find the solution, let us know. */
 593   if (counter == maxcount)
 594   {
 595       GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);
 596   }
 597
 598   return GSL_SUCCESS;
 599 }
 600
 601
 602 /* Eigenvalue solutions for characteristic values below. */
 603
 604
 605 /*  figi.c converted from EISPACK Fortran FIGI.F.
 606  *
 607  *   given a nonsymmetric tridiagonal matrix such that the products
 608  *    of corresponding pairs of off-diagonal elements are all
 609  *    non-negative, this subroutine reduces it to a symmetric
 610  *    tridiagonal matrix with the same eigenvalues.  if, further,
 611  *    a zero product only occurs when both factors are zero,
 612  *    the reduced matrix is similar to the original matrix.
 613  *
 614  *    on input
 615  *
 616  *       n is the order of the matrix.
 617  *
 618  *       t contains the input matrix.  its subdiagonal is
 619  *         stored in the last n-1 positions of the first column,
 620  *         its diagonal in the n positions of the second column,
 621  *         and its superdiagonal in the first n-1 positions of
 622  *         the third column.  t(1,1) and t(n,3) are arbitrary.
 623  *
 624  *    on output
 625  *
 626  *       t is unaltered.
 627  *
 628  *       d contains the diagonal elements of the symmetric matrix.
 629  *
 630  *       e contains the subdiagonal elements of the symmetric
 631  *         matrix in its last n-1 positions.  e(1) is not set.
 632  *
 633  *       e2 contains the squares of the corresponding elements of e.
 634  *         e2 may coincide with e if the squares are not needed.
 635  *
 636  *       ierr is set to
 637  *         zero       for normal return,
 638  *         n+i        if t(i,1)*t(i-1,3) is negative,
 639  *         -(3*n+i)   if t(i,1)*t(i-1,3) is zero with one factor
 640  *                    non-zero.  in this case, the eigenvectors of
 641  *                    the symmetric matrix are not simply related
 642  *                    to those of  t  and should not be sought.
 643  *
 644  *    questions and comments should be directed to burton s. garbow,
 645  *    mathematics and computer science div, argonne national laboratory
 646  *
 647  *    this version dated august 1983.
 648  */
 649 static int figi(int nn, double *tt, double *dd, double *ee,
 650                 double *e2)
 651 {
 652   int ii;
 653
 654   for (ii=0; ii<nn; ii++)
 655   {
 656       if (ii != 0)
 657       {
 658           e2[ii] = tt[3*ii]*tt[3*(ii-1)+2];
 659
 660           if (e2[ii] < 0.0)
 661           {
 662               /* set error -- product of some pair of off-diagonal
 663                  elements is negative */
 664               return (nn + ii);
 665           }
 666
 667           if (e2[ii] == 0.0 && (tt[3*ii] != 0.0 || tt[3*(ii-1)+2] != 0.0))
 668           {
 669               /* set error -- product of some pair of off-diagonal
 670                  elements is zero with one member non-zero */
 671               return (-1*(3*nn + ii));
 672           }
 673
 674           ee[ii] = sqrt(e2[ii]);
 675       }
 676
 677       dd[ii] = tt[3*ii+1];
 678   }
 679
 680   return 0;
 681 }
 682
 683
 684 int gsl_sf_mathieu_a_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])
 685 {
 686   unsigned int even_order = work->even_order, odd_order = work->odd_order,
 687       extra_values = work->extra_values, ii, jj;
 688   int status;
 689   double *tt = work->tt, *dd = work->dd, *ee = work->ee, *e2 = work->e2,
 690          *zz = work->zz, *aa = work->aa;
 691   gsl_matrix_view mat, evec;
 692   gsl_vector_view eval;
 693   gsl_eigen_symmv_workspace *wmat = work->wmat;
 694
 695   if (order_max > work->size || order_max <= order_min || order_min < 0)
 696     {
 697       GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);
 698     }
 699
 700   /* Convert the nonsymmetric tridiagonal matrix to a symmetric tridiagonal
 701      form. */
 702
 703   tt[0] = 0.0;
 704   tt[1] = 0.0;
 705   tt[2] = qq;
 706   for (ii=1; ii<even_order-1; ii++)
 707   {
 708       tt[3*ii] = qq;
 709       tt[3*ii+1] = 4*ii*ii;
 710       tt[3*ii+2] = qq;
 711   }
 712   tt[3*even_order-3] = qq;
 713   tt[3*even_order-2] = 4*(even_order - 1)*(even_order - 1);
 714   tt[3*even_order-1] = 0.0;
 715
 716   tt[3] *= 2;
 717
 718   status = figi((signed int)even_order, tt, dd, ee, e2);
 719
 720   if (status)
 721     {
 722       GSL_ERROR("Internal error in tridiagonal Mathieu matrix", GSL_EFAILED);
 723     }
 724
 725   /* Fill the period \pi matrix. */
 726   for (ii=0; ii<even_order*even_order; ii++)
 727       zz[ii] = 0.0;
 728
 729   zz[0] = dd[0];
 730   zz[1] = ee[1];
 731   for (ii=1; ii<even_order-1; ii++)
 732   {
 733       zz[ii*even_order+ii-1] = ee[ii];
 734       zz[ii*even_order+ii] = dd[ii];
 735       zz[ii*even_order+ii+1] = ee[ii+1];
 736   }
 737   zz[even_order*(even_order-1)+even_order-2] = ee[even_order-1];
 738   zz[even_order*even_order-1] = dd[even_order-1];
 739
 740   /* Compute (and sort) the eigenvalues of the matrix. */
 741   mat = gsl_matrix_view_array(zz, even_order, even_order);
 742   eval = gsl_vector_subvector(work->eval, 0, even_order);
 743   evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);
 744   gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
 745   gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
 746
 747   for (ii=0; ii<even_order-extra_values; ii++)
 748       aa[2*ii] = gsl_vector_get(&eval.vector, ii);
 749
 750   /* Fill the period 2\pi matrix. */
 751   for (ii=0; ii<odd_order*odd_order; ii++)
 752       zz[ii] = 0.0;
 753   for (ii=0; ii<odd_order; ii++)
 754       for (jj=0; jj<odd_order; jj++)
 755       {
 756           if (ii == jj)
 757               zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);
 758           else if (ii == jj + 1 || ii + 1 == jj)
 759               zz[ii*odd_order+jj] = qq;
 760       }
 761   zz[0] += qq;
 762
 763   /* Compute (and sort) the eigenvalues of the matrix. */
 764   mat = gsl_matrix_view_array(zz, odd_order, odd_order);
 765   eval = gsl_vector_subvector(work->eval, 0, odd_order);
 766   evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);
 767   gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
 768   gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
 769
 770   for (ii=0; ii<odd_order-extra_values; ii++)
 771       aa[2*ii+1] = gsl_vector_get(&eval.vector, ii);
 772
 773   for (ii = order_min ; ii <= order_max ; ii++)
 774     {
 775       result_array[ii - order_min] = aa[ii];
 776     }
 777
 778   return GSL_SUCCESS;
 779 }
 780
 781
 782 int gsl_sf_mathieu_b_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])
 783 {
 784   unsigned int even_order = work->even_order-1, odd_order = work->odd_order,
 785       extra_values = work->extra_values, ii, jj;
 786   double *zz = work->zz, *bb = work->bb;
 787   gsl_matrix_view mat, evec;
 788   gsl_vector_view eval;
 789   gsl_eigen_symmv_workspace *wmat = work->wmat;
 790
 791   if (order_max > work->size || order_max <= order_min || order_min < 0)
 792     {
 793       GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);
 794     }
 795
 796   /* Fill the period \pi matrix. */
 797   for (ii=0; ii<even_order*even_order; ii++)
 798       zz[ii] = 0.0;
 799   for (ii=0; ii<even_order; ii++)
 800       for (jj=0; jj<even_order; jj++)
 801       {
 802           if (ii == jj)
 803               zz[ii*even_order+jj] = 4*(ii + 1)*(ii + 1);
 804           else if (ii == jj + 1 || ii + 1 == jj)
 805               zz[ii*even_order+jj] = qq;
 806       }
 807
 808   /* Compute (and sort) the eigenvalues of the matrix. */
 809   mat = gsl_matrix_view_array(zz, even_order, even_order);
 810   eval = gsl_vector_subvector(work->eval, 0, even_order);
 811   evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);
 812   gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
 813   gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
 814
 815   bb[0] = 0.0;
 816   for (ii=0; ii<even_order-extra_values; ii++)
 817       bb[2*(ii+1)] = gsl_vector_get(&eval.vector, ii);
 818
 819   /* Fill the period 2\pi matrix. */
 820   for (ii=0; ii<odd_order*odd_order; ii++)
 821       zz[ii] = 0.0;
 822   for (ii=0; ii<odd_order; ii++)
 823       for (jj=0; jj<odd_order; jj++)
 824       {
 825           if (ii == jj)
 826               zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);
 827           else if (ii == jj + 1 || ii + 1 == jj)
 828               zz[ii*odd_order+jj] = qq;
 829       }
 830
 831   zz[0] -= qq;
 832
 833   /* Compute (and sort) the eigenvalues of the matrix. */
 834   mat = gsl_matrix_view_array(zz, odd_order, odd_order);
 835   eval = gsl_vector_subvector(work->eval, 0, odd_order);
 836   evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);
 837   gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
 838   gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
 839
 840   for (ii=0; ii<odd_order-extra_values; ii++)
 841       bb[2*ii+1] = gsl_vector_get(&eval.vector, ii);
 842
 843   for (ii = order_min ; ii <= order_max ; ii++)
 844     {
 845       result_array[ii - order_min] = bb[ii];
 846     }
 847
 848   return GSL_SUCCESS;
 849 }