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

   1 /* integration/test.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 <stdio.h>
  23 #include <math.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_integration.h>
  26 #include <gsl/gsl_errno.h>
  27 #include <gsl/gsl_test.h>
  28 #include <gsl/gsl_ieee_utils.h>
  29
  30 #include "tests.h"
  31
  32 gsl_function make_function (double (* f) (double, void *), double * p);
  33
  34 gsl_function make_function (double (* f) (double, void *), double * p)
  35 {
  36   gsl_function f_new;
  37
  38   f_new.function = f ;
  39   f_new.params = p ;
  40
  41   return f_new;
  42 }
  43
  44 struct counter_params {
  45   gsl_function * f;
  46   int neval;
  47 } ;
  48
  49 double counter (double x, void * params);
  50 gsl_function make_counter (gsl_function * f, struct counter_params * p);
  51
  52 double
  53 counter (double x, void * params)
  54 {
  55   struct counter_params * p = (struct counter_params *) params;
  56   p->neval++ ; /* increment counter */
  57   return GSL_FN_EVAL(p->f, x);
  58 }
  59
  60 gsl_function make_counter (gsl_function * f, struct counter_params * p)
  61 {
  62   gsl_function f_new;
  63
  64   p->f = f;
  65   p->neval = 0 ;
  66
  67   f_new.function = &counter ;
  68   f_new.params = p ;
  69
  70   return f_new;
  71 }
  72
  73 void my_error_handler (const char *reason, const char *file,
  74                        int line, int err);
  75
  76 int main (void)
  77 {
  78   gsl_ieee_env_setup ();
  79   gsl_set_error_handler (&my_error_handler);
  80
  81   /* Test the basic Gauss-Kronrod rules with a smooth positive function. */
  82
  83   {
  84     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
  85     double exp_result = 7.716049357767090777E-02;
  86     double exp_abserr = 2.990224871000550874E-06;
  87     double exp_resabs = 7.716049357767090777E-02;
  88     double exp_resasc = 4.434273814139995384E-02;
  89
  90     double alpha = 2.6 ;
  91     gsl_function f = make_function(&f1, &alpha) ;
  92
  93     gsl_integration_qk15 (&f, 0.0, 1.0,
  94                                   &result, &abserr, &resabs, &resasc) ;
  95     gsl_test_rel(result,exp_result,1e-15,"qk15(f1) smooth result") ;
  96     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) smooth abserr") ;
  97     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) smooth resabs") ;
  98     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) smooth resasc") ;
  99
 100     gsl_integration_qk15 (&f, 1.0, 0.0,
 101                                   &result, &abserr, &resabs, &resasc) ;
 102
 103     gsl_test_rel(result,-exp_result,1e-15,"qk15(f1) reverse result") ;
 104     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) reverse abserr") ;
 105     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) reverse resabs") ;
 106     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) reverse resasc") ;
 107   }
 108
 109   {
 110     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 111     double exp_result = 7.716049379303084599E-02;
 112     double exp_abserr = 9.424302194248481445E-08;
 113     double exp_resabs = 7.716049379303084599E-02;
 114     double exp_resasc = 4.434311425038358484E-02;
 115
 116     double alpha = 2.6 ;
 117     gsl_function f = make_function(&f1, &alpha);
 118
 119     gsl_integration_qk21 (&f, 0.0, 1.0,
 120                                   &result, &abserr, &resabs, &resasc) ;
 121     gsl_test_rel(result,exp_result,1e-15,"qk21(f1) smooth result") ;
 122     gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) smooth abserr") ;
 123     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) smooth resabs") ;
 124     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) smooth resasc") ;
 125
 126     gsl_integration_qk21 (&f, 1.0, 0.0,
 127                                   &result, &abserr, &resabs, &resasc) ;
 128     gsl_test_rel(result,-exp_result,1e-15,"qk21(f1) reverse result") ;
 129     gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) reverse abserr") ;
 130     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) reverse resabs") ;
 131     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) reverse resasc") ;
 132   }
 133
 134   {
 135     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 136     double exp_result = 7.716049382494900855E-02;
 137     double exp_abserr = 1.713503193600029893E-09;
 138     double exp_resabs = 7.716049382494900855E-02;
 139     double exp_resasc = 4.427995051868838933E-02;
 140
 141     double alpha = 2.6 ;
 142     gsl_function f = make_function(&f1, &alpha);
 143
 144     gsl_integration_qk31 (&f, 0.0, 1.0,
 145                                   &result, &abserr, &resabs, &resasc) ;
 146     gsl_test_rel(result,exp_result,1e-15,"qk31(f1) smooth result") ;
 147     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) smooth abserr") ;
 148     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) smooth resabs") ;
 149     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) smooth resasc") ;
 150
 151     gsl_integration_qk31 (&f, 1.0, 0.0,
 152                                   &result, &abserr, &resabs, &resasc) ;
 153     gsl_test_rel(result,-exp_result,1e-15,"qk31(f1) reverse result") ;
 154     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) reverse abserr") ;
 155     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) reverse resabs") ;
 156     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) reverse resasc") ;
 157   }
 158
 159   {
 160     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 161     double exp_result = 7.716049382681375302E-02;
 162     double exp_abserr = 9.576386660975511224E-11;
 163     double exp_resabs = 7.716049382681375302E-02;
 164     double exp_resasc = 4.421521169637691873E-02;
 165
 166     double alpha = 2.6 ;
 167     gsl_function f = make_function(&f1, &alpha);
 168
 169     gsl_integration_qk41 (&f, 0.0, 1.0,
 170                                   &result, &abserr, &resabs, &resasc) ;
 171     gsl_test_rel(result,exp_result,1e-15,"qk41(f1) smooth result") ;
 172     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) smooth abserr") ;
 173     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) smooth resabs") ;
 174     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) smooth resasc") ;
 175
 176     gsl_integration_qk41 (&f, 1.0, 0.0,
 177                                   &result, &abserr, &resabs, &resasc) ;
 178     gsl_test_rel(result,-exp_result,1e-15,"qk41(f1) reverse result") ;
 179     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) reverse abserr") ;
 180     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) reverse resabs") ;
 181     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) reverse resasc") ;
 182   }
 183
 184   {
 185     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 186     double exp_result = 7.716049382708510540E-02;
 187     double exp_abserr = 1.002079980317363772E-11;
 188     double exp_resabs = 7.716049382708510540E-02;
 189     double exp_resasc = 4.416474291216854892E-02;
 190
 191     double alpha = 2.6 ;
 192     gsl_function f = make_function(&f1, &alpha);
 193
 194     gsl_integration_qk51 (&f, 0.0, 1.0,
 195                                   &result, &abserr, &resabs, &resasc) ;
 196     gsl_test_rel(result,exp_result,1e-15,"qk51(f1) smooth result") ;
 197     gsl_test_rel(abserr,exp_abserr,1e-5,"qk51(f1) smooth abserr") ;
 198     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) smooth resabs") ;
 199     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) smooth resasc") ;
 200
 201     gsl_integration_qk51 (&f, 1.0, 0.0,
 202                                   &result, &abserr, &resabs, &resasc) ;
 203     gsl_test_rel(result,-exp_result,1e-15,"qk51(f1) reverse result") ;
 204     gsl_test_rel(abserr,exp_abserr,1e-5,"qk51(f1) reverse abserr") ;
 205     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) reverse resabs") ;
 206     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) reverse resasc") ;
 207   }
 208
 209   {
 210     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 211     double exp_result = 7.716049382713800753E-02;
 212     double exp_abserr = 1.566060362296155616E-12;
 213     double exp_resabs = 7.716049382713800753E-02;
 214     double exp_resasc = 4.419287685934316506E-02;
 215
 216     double alpha = 2.6 ;
 217     gsl_function f = make_function(&f1, &alpha);
 218
 219     gsl_integration_qk61 (&f, 0.0, 1.0,
 220                                   &result, &abserr, &resabs, &resasc) ;
 221     gsl_test_rel(result,exp_result,1e-15,"qk61(f1) smooth result") ;
 222     gsl_test_rel(abserr,exp_abserr,1e-5,"qk61(f1) smooth abserr") ;
 223     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) smooth resabs") ;
 224     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) smooth resasc") ;
 225
 226     gsl_integration_qk61 (&f, 1.0, 0.0,
 227                                   &result, &abserr, &resabs, &resasc) ;
 228     gsl_test_rel(result,-exp_result,1e-15,"qk61(f1) reverse result") ;
 229     gsl_test_rel(abserr,exp_abserr,1e-5,"qk61(f1) reverse abserr") ;
 230     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) reverse resabs") ;
 231     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) reverse resasc") ;
 232   }
 233
 234   /* Now test the basic rules with a positive function that has a
 235      singularity. This should give large values of abserr which would
 236      find discrepancies in the abserr calculation. */
 237
 238   {
 239     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 240     double exp_result = 1.555688196612745777E+01;
 241     double exp_abserr = 2.350164577239293706E+01;
 242     double exp_resabs = 1.555688196612745777E+01;
 243     double exp_resasc = 2.350164577239293706E+01;
 244
 245     double alpha = -0.9 ;
 246     gsl_function f = make_function(&f1, &alpha);
 247
 248     gsl_integration_qk15 (&f, 0.0, 1.0,
 249                                   &result, &abserr, &resabs, &resasc) ;
 250     gsl_test_rel(result,exp_result,1e-15,"qk15(f1) singular result") ;
 251     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) singular abserr") ;
 252     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) singular resabs") ;
 253     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) singular resasc") ;
 254
 255     gsl_integration_qk15 (&f, 1.0, 0.0,
 256                                   &result, &abserr, &resabs, &resasc) ;
 257     gsl_test_rel(result,-exp_result,1e-15,"qk15(f1) reverse result") ;
 258     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) reverse abserr") ;
 259     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) reverse resabs") ;
 260     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) reverse resasc") ;
 261   }
 262
 263   {
 264     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 265     double exp_result = 1.799045317938126232E+01;
 266     double exp_abserr = 2.782360287710622515E+01;
 267     double exp_resabs = 1.799045317938126232E+01;
 268     double exp_resasc = 2.782360287710622515E+01;
 269
 270     double alpha = -0.9 ;
 271     gsl_function f = make_function(&f1, &alpha);
 272
 273     gsl_integration_qk21 (&f, 0.0, 1.0,
 274                                   &result, &abserr, &resabs, &resasc) ;
 275     gsl_test_rel(result,exp_result,1e-15,"qk21(f1) singular result") ;
 276     gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) singular abserr") ;
 277     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) singular resabs") ;
 278     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) singular resasc") ;
 279
 280     gsl_integration_qk21 (&f, 1.0, 0.0,
 281                                   &result, &abserr, &resabs, &resasc) ;
 282     gsl_test_rel(result,-exp_result,1e-15,"qk21(f1) reverse result") ;
 283     gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) reverse abserr") ;
 284     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) reverse resabs") ;
 285     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) reverse resasc") ;
 286   }
 287
 288   {
 289     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 290     double exp_result = 2.081873305159121657E+01;
 291     double exp_abserr = 3.296500137482590276E+01;
 292     double exp_resabs = 2.081873305159121301E+01;
 293     double exp_resasc = 3.296500137482590276E+01;
 294
 295     double alpha = -0.9 ;
 296     gsl_function f = make_function(&f1, &alpha);
 297
 298     gsl_integration_qk31 (&f, 0.0, 1.0,
 299                                   &result, &abserr, &resabs, &resasc) ;
 300     gsl_test_rel(result,exp_result,1e-15,"qk31(f1) singular result") ;
 301     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) singular abserr") ;
 302     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) singular resabs") ;
 303     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) singular resasc") ;
 304
 305     gsl_integration_qk31 (&f, 1.0, 0.0,
 306                                   &result, &abserr, &resabs, &resasc) ;
 307     gsl_test_rel(result,-exp_result,1e-15,"qk31(f1) reverse result") ;
 308     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) reverse abserr") ;
 309     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) reverse resabs") ;
 310     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) reverse resasc") ;
 311   }
 312
 313   {
 314     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 315     double exp_result = 2.288677623903126701E+01;
 316     double exp_abserr = 3.671538820274916048E+01;
 317     double exp_resabs = 2.288677623903126701E+01;
 318     double exp_resasc = 3.671538820274916048E+01;
 319
 320     double alpha = -0.9 ;
 321     gsl_function f = make_function(&f1, &alpha);
 322
 323     gsl_integration_qk41 (&f, 0.0, 1.0,
 324                                   &result, &abserr, &resabs, &resasc) ;
 325     gsl_test_rel(result,exp_result,1e-15,"qk41(f1) singular result") ;
 326     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) singular abserr") ;
 327     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) singular resabs") ;
 328     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) singular resasc") ;
 329
 330     gsl_integration_qk41 (&f, 1.0, 0.0,
 331                                   &result, &abserr, &resabs, &resasc) ;
 332     gsl_test_rel(result,-exp_result,1e-15,"qk41(f1) reverse result") ;
 333     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) reverse abserr") ;
 334     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) reverse resabs") ;
 335     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) reverse resasc") ;
 336   }
 337
 338   {
 339     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 340     double exp_result = 2.449953612016972215E+01;
 341     double exp_abserr = 3.967771249391228849E+01;
 342     double exp_resabs = 2.449953612016972215E+01;
 343     double exp_resasc = 3.967771249391228849E+01;
 344
 345     double alpha = -0.9 ;
 346     gsl_function f = make_function(&f1, &alpha);
 347
 348     gsl_integration_qk51 (&f, 0.0, 1.0,
 349                                   &result, &abserr, &resabs, &resasc) ;
 350     gsl_test_rel(result,exp_result,1e-15,"qk51(f1) singular result") ;
 351     gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f1) singular abserr") ;
 352     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) singular resabs") ;
 353     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) singular resasc") ;
 354
 355     gsl_integration_qk51 (&f, 1.0, 0.0,
 356                                   &result, &abserr, &resabs, &resasc) ;
 357     gsl_test_rel(result,-exp_result,1e-15,"qk51(f1) reverse result") ;
 358     gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f1) reverse abserr") ;
 359     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) reverse resabs") ;
 360     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) reverse resasc") ;
 361   }
 362
 363   {
 364     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 365     double exp_result = 2.583030240976628988E+01;
 366     double exp_abserr = 4.213750493076978643E+01;
 367     double exp_resabs = 2.583030240976628988E+01;
 368     double exp_resasc = 4.213750493076978643E+01;
 369
 370     double alpha = -0.9 ;
 371     gsl_function f = make_function(&f1, &alpha);
 372
 373     gsl_integration_qk61 (&f, 0.0, 1.0,
 374                                   &result, &abserr, &resabs, &resasc) ;
 375     gsl_test_rel(result,exp_result,1e-15,"qk61(f1) singular result") ;
 376     gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f1) singular abserr") ;
 377     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) singular resabs") ;
 378     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) singular resasc") ;
 379
 380     gsl_integration_qk61 (&f, 1.0, 0.0,
 381                                   &result, &abserr, &resabs, &resasc) ;
 382     gsl_test_rel(result,-exp_result,1e-15,"qk61(f1) reverse result") ;
 383     gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f1) reverse abserr") ;
 384     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) reverse resabs") ;
 385     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) reverse resasc") ;
 386   }
 387
 388   /* Test the basic Gauss-Kronrod rules with a smooth oscillating
 389      function, over an unsymmetric range. This should find any
 390      discrepancies in the abscissae. */
 391
 392   {
 393     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 394     double exp_result =-7.238969575483799046E-01;
 395     double exp_abserr = 8.760080200939757174E-06;
 396     double exp_resabs = 1.165564172429140788E+00;
 397     double exp_resasc = 9.334560307787327371E-01;
 398
 399     double alpha = 1.3 ;
 400     gsl_function f = make_function(&f3, &alpha);
 401
 402     gsl_integration_qk15 (&f, 0.3, 2.71,
 403                                   &result, &abserr, &resabs, &resasc) ;
 404     gsl_test_rel(result,exp_result,1e-15,"qk15(f3) oscill result") ;
 405     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f3) oscill abserr") ;
 406     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f3) oscill resabs") ;
 407     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f3) oscill resasc") ;
 408
 409     gsl_integration_qk15 (&f, 2.71, 0.3,
 410                                   &result, &abserr, &resabs, &resasc) ;
 411     gsl_test_rel(result,-exp_result,1e-15,"qk15(f3) reverse result") ;
 412     gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f3) reverse abserr") ;
 413     gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f3) reverse resabs") ;
 414     gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f3) reverse resasc") ;
 415   }
 416
 417   {
 418     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 419     double exp_result =-7.238969575482959717E-01;
 420     double exp_abserr = 7.999213141433641888E-11;
 421     double exp_resabs = 1.150829032708484023E+00;
 422     double exp_resasc = 9.297591249133687619E-01;
 423
 424     double alpha = 1.3 ;
 425     gsl_function f = make_function(&f3, &alpha);
 426
 427     gsl_integration_qk21 (&f, 0.3, 2.71,
 428                                   &result, &abserr, &resabs, &resasc) ;
 429     gsl_test_rel(result,exp_result,1e-15,"qk21(f3) oscill result") ;
 430     gsl_test_rel(abserr,exp_abserr,1e-5,"qk21(f3) oscill abserr") ;
 431     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f3) oscill resabs") ;
 432     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f3) oscill resasc") ;
 433
 434     gsl_integration_qk21 (&f, 2.71, 0.3,
 435                                   &result, &abserr, &resabs, &resasc) ;
 436     gsl_test_rel(result,-exp_result,1e-15,"qk21(f3) reverse result") ;
 437     gsl_test_rel(abserr,exp_abserr,1e-5,"qk21(f3) reverse abserr") ;
 438     gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f3) reverse resabs") ;
 439     gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f3) reverse resasc") ;
 440   }
 441
 442   {
 443     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 444     double exp_result =-7.238969575482959717E-01;
 445     double exp_abserr = 1.285805464427459261E-14;
 446     double exp_resabs = 1.158150602093290571E+00;
 447     double exp_resasc = 9.277828092501518853E-01;
 448
 449     double alpha = 1.3 ;
 450     gsl_function f = make_function(&f3, &alpha);
 451
 452     gsl_integration_qk31 (&f, 0.3, 2.71,
 453                                   &result, &abserr, &resabs, &resasc) ;
 454     gsl_test_rel(result,exp_result,1e-15,"qk31(f3) oscill result") ;
 455     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f3) oscill abserr") ;
 456     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f3) oscill resabs") ;
 457     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f3) oscill resasc") ;
 458
 459     gsl_integration_qk31 (&f, 2.71, 0.3,
 460                                   &result, &abserr, &resabs, &resasc) ;
 461     gsl_test_rel(result,-exp_result,1e-15,"qk31(f3) reverse result") ;
 462     gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f3) reverse abserr") ;
 463     gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f3) reverse resabs") ;
 464     gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f3) reverse resasc") ;
 465   }
 466
 467   {
 468     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 469     double exp_result =-7.238969575482959717E-01;
 470     double exp_abserr = 1.286535726271015626E-14;
 471     double exp_resabs = 1.158808363486595328E+00;
 472     double exp_resasc = 9.264382258645686985E-01;
 473
 474     double alpha = 1.3 ;
 475     gsl_function f = make_function(&f3, &alpha);
 476
 477     gsl_integration_qk41 (&f, 0.3, 2.71,
 478                                   &result, &abserr, &resabs, &resasc) ;
 479     gsl_test_rel(result,exp_result,1e-15,"qk41(f3) oscill result") ;
 480     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f3) oscill abserr") ;
 481     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f3) oscill resabs") ;
 482     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f3) oscill resasc") ;
 483
 484     gsl_integration_qk41 (&f, 2.71, 0.3,
 485                                   &result, &abserr, &resabs, &resasc) ;
 486     gsl_test_rel(result,-exp_result,1e-15,"qk41(f3) reverse result") ;
 487     gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f3) reverse abserr") ;
 488     gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f3) reverse resabs") ;
 489     gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f3) reverse resasc") ;
 490   }
 491
 492   {
 493     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 494     double exp_result =-7.238969575482961938E-01;
 495     double exp_abserr = 1.285290995039385778E-14;
 496     double exp_resabs = 1.157687209264406381E+00;
 497     double exp_resasc = 9.264666884071264263E-01;
 498
 499     double alpha = 1.3 ;
 500     gsl_function f = make_function(&f3, &alpha);
 501
 502     gsl_integration_qk51 (&f, 0.3, 2.71,
 503                                   &result, &abserr, &resabs, &resasc) ;
 504     gsl_test_rel(result,exp_result,1e-15,"qk51(f3) oscill result") ;
 505     gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f3) oscill abserr") ;
 506     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f3) oscill resabs") ;
 507     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f3) oscill resasc") ;
 508
 509     gsl_integration_qk51 (&f, 2.71, 0.3,
 510                                   &result, &abserr, &resabs, &resasc) ;
 511     gsl_test_rel(result,-exp_result,1e-15,"qk51(f3) reverse result") ;
 512     gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f3) reverse abserr") ;
 513     gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f3) reverse resabs") ;
 514     gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f3) reverse resasc") ;
 515   }
 516
 517   {
 518     double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
 519     double exp_result =-7.238969575482959717E-01;
 520     double exp_abserr = 1.286438572027470736E-14;
 521     double exp_resabs = 1.158720854723590099E+00;
 522     double exp_resasc = 9.270469641771273972E-01;
 523
 524     double alpha = 1.3 ;
 525     gsl_function f = make_function(&f3, &alpha);
 526
 527     gsl_integration_qk61 (&f, 0.3, 2.71,
 528                                   &result, &abserr, &resabs, &resasc) ;
 529     gsl_test_rel(result,exp_result,1e-15,"qk61(f3) oscill result") ;
 530     gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f3) oscill abserr") ;
 531     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f3) oscill resabs") ;
 532     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f3) oscill resasc") ;
 533
 534     gsl_integration_qk61 (&f, 2.71, 0.3,
 535                                   &result, &abserr, &resabs, &resasc) ;
 536     gsl_test_rel(result,-exp_result,1e-15,"qk61(f3) reverse result") ;
 537     gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f3) reverse abserr") ;
 538     gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f3) reverse resabs") ;
 539     gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f3) reverse resasc") ;
 540   }
 541
 542   /* Test the non-adaptive gaussian integrator QNG */
 543
 544   {
 545     int status = 0; size_t neval = 0 ;
 546     double result = 0, abserr = 0 ;
 547     double exp_result = 7.716049379303083211E-02;
 548     double exp_abserr = 9.424302199601294244E-08;
 549     int exp_neval  =  21;
 550     int exp_ier    =   0;
 551
 552     double alpha = 2.6 ;
 553     gsl_function f = make_function(&f1, &alpha);
 554
 555     status = gsl_integration_qng (&f, 0.0, 1.0, 1e-1, 0.0,
 556                                   &result, &abserr, &neval) ;
 557     gsl_test_rel(result,exp_result,1e-15,"qng(f1) smooth result") ;
 558     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) smooth abserr") ;
 559     gsl_test_int((int)neval,exp_neval,"qng(f1) smooth neval") ;
 560     gsl_test_int(status,exp_ier,"qng(f1) smooth status") ;
 561
 562     status = gsl_integration_qng (&f, 1.0, 0.0, 1e-1, 0.0,
 563                                   &result, &abserr, &neval) ;
 564     gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse result") ;
 565     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) reverse abserr") ;
 566     gsl_test_int((int)neval,exp_neval,"qng(f1) reverse neval") ;
 567     gsl_test_int(status,exp_ier,"qng(f1) reverse status") ;
 568   }
 569
 570   {
 571     int status = 0; size_t neval = 0 ;
 572     double result = 0, abserr = 0 ;
 573
 574     double exp_result = 7.716049382706505200E-02;
 575     double exp_abserr = 2.666893044866214501E-12;
 576     int exp_neval  =  43;
 577     int exp_ier    =   0;
 578
 579     double alpha = 2.6 ;
 580     gsl_function f = make_function(&f1, &alpha);
 581
 582     status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-9,
 583                                   &result, &abserr, &neval) ;
 584     gsl_test_rel(result,exp_result,1e-15,"qng(f1) smooth 43pt result") ;
 585     gsl_test_rel(abserr,exp_abserr,1e-5,"qng(f1) smooth 43pt abserr") ;
 586     gsl_test_int((int)neval,exp_neval,"qng(f1) smooth 43pt neval") ;
 587     gsl_test_int(status,exp_ier,"qng(f1) smooth 43pt status") ;
 588
 589     status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-9,
 590                                   &result, &abserr, &neval) ;
 591     gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse 43pt result") ;
 592     gsl_test_rel(abserr,exp_abserr,1e-5,"qng(f1) reverse 43pt abserr") ;
 593     gsl_test_int((int)neval,exp_neval,"qng(f1) reverse 43pt neval") ;
 594     gsl_test_int(status,exp_ier,"qng(f1) reverse 43pt status") ;
 595   }
 596
 597   {
 598     int status; size_t neval = 0 ;
 599     double result = 0, abserr = 0 ;
 600     double exp_result =-7.238969575482961938E-01;
 601     double exp_abserr = 1.277676889520056369E-14;
 602     int exp_neval  =  43;
 603     int exp_ier    =   0;
 604
 605     double alpha = 1.3 ;
 606     gsl_function f = make_function(&f3, &alpha);
 607
 608     status = gsl_integration_qng (&f, 0.3, 2.71, 0.0, 1e-12,
 609                                   &result, &abserr, &neval) ;
 610     gsl_test_rel(result,exp_result,1e-15,"qnq(f3) oscill result") ;
 611     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f3) oscill abserr") ;
 612     gsl_test_int((int)neval,exp_neval,"qng(f3) oscill neval") ;
 613     gsl_test_int(status,exp_ier,"qng(f3) oscill status") ;
 614
 615     status = gsl_integration_qng (&f, 2.71, 0.3, 0.0, 1e-12,
 616                                   &result, &abserr, &neval) ;
 617     gsl_test_rel(result,-exp_result,1e-15,"qnq(f3) reverse result") ;
 618     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f3) reverse abserr") ;
 619     gsl_test_int((int)neval,exp_neval,"qng(f3) reverse neval") ;
 620     gsl_test_int(status,exp_ier,"qng(f3) reverse status") ;
 621   }
 622
 623   {
 624     int status = 0; size_t neval = 0 ;
 625     double result = 0, abserr = 0 ;
 626
 627     double exp_result = 7.716049382716029525E-02;
 628     double exp_abserr = 8.566535680046930668E-16;
 629     int exp_neval  =  87;
 630     int exp_ier    =   0;
 631
 632     double alpha = 2.6 ;
 633     gsl_function f = make_function(&f1, &alpha);
 634
 635     status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-13,
 636                                   &result, &abserr, &neval) ;
 637     gsl_test_rel(result,exp_result,1e-15,"qng(f1) 87pt smooth result") ;
 638     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) 87pt smooth abserr") ;
 639     gsl_test_int((int)neval,exp_neval,"qng(f1) 87pt smooth neval") ;
 640     gsl_test_int(status,exp_ier,"qng(f1) 87pt smooth status") ;
 641
 642     status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-13,
 643                                   &result, &abserr, &neval) ;
 644     gsl_test_rel(result,-exp_result,1e-15,"qng(f1) 87pt reverse result") ;
 645     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) 87pt reverse abserr") ;
 646     gsl_test_int((int)neval,exp_neval,"qng(f1) 87pt reverse neval") ;
 647     gsl_test_int(status,exp_ier,"qng(f1) 87pt reverse status") ;
 648   }
 649
 650   {
 651     int status = 0; size_t neval = 0 ;
 652     double result = 0, abserr = 0 ;
 653
 654     double exp_result = 3.222948711817264211E+01;
 655     double exp_abserr = 2.782360287710622870E+01;
 656     int exp_neval  =  87;
 657     int exp_ier    =  GSL_ETOL;
 658
 659     double alpha = -0.9 ;
 660     gsl_function f = make_function(&f1, &alpha);
 661
 662     status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-3,
 663                                   &result, &abserr, &neval) ;
 664     gsl_test_rel(result,exp_result,1e-15,"qng(f1) sing beyond 87pt result");
 665     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) sing beyond 87pt abserr");
 666     gsl_test_int((int)neval,exp_neval,"qng(f1) sing beyond 87pt neval") ;
 667     gsl_test_int(status,exp_ier,"qng(f1) sing beyond 87pt status") ;
 668
 669     status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-3,
 670                                   &result, &abserr, &neval) ;
 671     gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse beyond 87pt result");
 672     gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) rev beyond 87pt abserr");
 673     gsl_test_int((int)neval,exp_neval,"qng(f1) rev beyond 87pt neval") ;
 674     gsl_test_int(status,exp_ier,"qng(f1) rev beyond 87pt status") ;
 675   }
 676
 677   /* Test the adaptive integrator QAG */
 678
 679   {
 680     int status = 0, i; struct counter_params p;
 681     double result = 0, abserr=0;
 682
 683     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
 684
 685     double exp_result = 7.716049382715854665E-02 ;
 686     double exp_abserr = 6.679384885865053037E-12 ;
 687     int exp_neval  =     165;
 688     int exp_ier    =       0;
 689     int exp_last   =       6;
 690
 691     double a[6] = { 0, 0.5, 0.25, 0.125, 0.0625, 0.03125 } ;
 692     double b[6] = { 0.03125, 1, 0.5, 0.25, 0.125, 0.0625 } ;
 693     double r[6] = { 3.966769831709074375E-06, 5.491842501998222409E-02,
 694                     1.909827770934243926E-02, 2.776531175604360531E-03,
 695                     3.280661030752063693E-04, 3.522704932261797744E-05 } ;
 696     double e[6] = { 6.678528276336181873E-12, 6.097169993333454062E-16,
 697                     2.120334764359736934E-16, 3.082568839745514608E-17,
 698                     3.642265412331439511E-18, 3.910988124757650942E-19 } ;
 699     int order[6] = { 1, 2, 3, 4, 5, 6 } ;
 700
 701     double alpha = 2.6 ;
 702     gsl_function f = make_function(&f1, &alpha) ;
 703
 704     gsl_function fc = make_counter(&f, &p) ;
 705
 706     status = gsl_integration_qag (&fc, 0.0, 1.0, 0.0, 1e-10, w->limit,
 707                                   GSL_INTEG_GAUSS15, w,
 708                                   &result, &abserr) ;
 709
 710     gsl_test_rel(result,exp_result,1e-15,"qag(f1) smooth result") ;
 711     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1) smooth abserr") ;
 712     gsl_test_int((int)(p.neval),exp_neval,"qag(f1) smooth neval") ;
 713     gsl_test_int((int)(w->size),exp_last,"qag(f1) smooth last") ;
 714     gsl_test_int(status,exp_ier,"qag(f1) smooth status") ;
 715
 716     for (i = 0; i < 6 ; i++)
 717         gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f1) smooth alist") ;
 718
 719     for (i = 0; i < 6 ; i++)
 720         gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f1) smooth blist") ;
 721
 722     for (i = 0; i < 6 ; i++)
 723         gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f1) smooth rlist") ;
 724
 725     for (i = 0; i < 6 ; i++)
 726         gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f1) smooth elist") ;
 727
 728     for (i = 0; i < 6 ; i++)
 729         gsl_test_int((int)w->order[i],order[i]-1,"qag(f1) smooth order") ;
 730
 731     p.neval = 0;
 732
 733     status = gsl_integration_qag (&fc, 1.0, 0.0, 0.0, 1e-10, w->limit,
 734                                   GSL_INTEG_GAUSS15, w,
 735                                   &result, &abserr) ;
 736
 737     gsl_test_rel(result,-exp_result,1e-15,"qag(f1) reverse result") ;
 738     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1) reverse abserr") ;
 739     gsl_test_int((int)(p.neval),exp_neval,"qag(f1) reverse neval") ;
 740     gsl_test_int((int)(w->size),exp_last,"qag(f1) reverse last") ;
 741     gsl_test_int(status,exp_ier,"qag(f1) reverse status") ;
 742
 743     gsl_integration_workspace_free (w) ;
 744
 745   }
 746
 747   /* Test the same function using an absolute error bound and the
 748      21-point rule */
 749
 750   {
 751     int status = 0, i; struct counter_params p;
 752     double result = 0, abserr=0;
 753
 754     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
 755
 756     double exp_result = 7.716049382716050342E-02 ;
 757     double exp_abserr = 2.227969521869139532E-15 ;
 758     int exp_neval  =     315;
 759     int exp_ier    =       0;
 760     int exp_last   =       8;
 761
 762     double a[8] = { 0, 0.5, 0.25, 0.125, 0.0625, 0.03125, 0.015625,
 763                     0.0078125 } ;
 764     double b[8] = { 0.0078125, 1, 0.5, 0.25, 0.125, 0.0625, 0.03125,
 765                     0.015625 } ;
 766     double r[8] = { 3.696942726831556522E-08, 5.491842501998223103E-02,
 767                     1.909827770934243579E-02, 2.776531175604360097E-03,
 768                     3.280661030752062609E-04, 3.522704932261797744E-05,
 769                     3.579060884684503576E-06, 3.507395216921808047E-07 } ;
 770     double e[8] = { 1.371316364034059572E-15, 6.097169993333454062E-16,
 771                     2.120334764359736441E-16, 3.082568839745514608E-17,
 772                     3.642265412331439511E-18, 3.910988124757650460E-19,
 773                     3.973555800712018091E-20, 3.893990926286736620E-21 } ;
 774     int order[8] = { 1, 2, 3, 4, 5, 6, 7, 8 } ;
 775
 776     double alpha = 2.6 ;
 777     gsl_function f = make_function(&f1, &alpha);
 778
 779     gsl_function fc = make_counter(&f, &p) ;
 780
 781     status = gsl_integration_qag (&fc, 0.0, 1.0, 1e-14, 0.0, w->limit,
 782                                   GSL_INTEG_GAUSS21, w,
 783                                   &result, &abserr) ;
 784
 785     gsl_test_rel(result,exp_result,1e-15,"qag(f1,21pt) smooth result") ;
 786     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1,21pt) smooth abserr") ;
 787     gsl_test_int((int)(p.neval),exp_neval,"qag(f1,21pt) smooth neval") ;
 788     gsl_test_int((int)(w->size),exp_last,"qag(f1,21pt) smooth last") ;
 789     gsl_test_int(status,exp_ier,"qag(f1,21pt) smooth status") ;
 790
 791     for (i = 0; i < 8 ; i++)
 792         gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f1,21pt) smooth alist") ;
 793
 794     for (i = 0; i < 8 ; i++)
 795         gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f1,21pt) smooth blist") ;
 796
 797     for (i = 0; i < 8 ; i++)
 798         gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f1,21pt) smooth rlist") ;
 799
 800     for (i = 0; i < 8 ; i++)
 801         gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f1,21pt) smooth elist") ;
 802
 803     for (i = 0; i < 8 ; i++)
 804         gsl_test_int((int)w->order[i],order[i]-1,"qag(f1,21pt) smooth order");
 805
 806
 807     p.neval = 0;
 808     status = gsl_integration_qag (&fc, 1.0, 0.0, 1e-14, 0.0, w->limit,
 809                                   GSL_INTEG_GAUSS21, w,
 810                                   &result, &abserr) ;
 811
 812     gsl_test_rel(result,-exp_result,1e-15,"qag(f1,21pt) reverse result") ;
 813     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1,21pt) reverse abserr") ;
 814     gsl_test_int((int)(p.neval),exp_neval,"qag(f1,21pt) reverse neval") ;
 815     gsl_test_int((int)(w->size),exp_last,"qag(f1,21pt) reverse last") ;
 816     gsl_test_int(status,exp_ier,"qag(f1,21pt) reverse status") ;
 817
 818     gsl_integration_workspace_free (w) ;
 819
 820   }
 821
 822   /* Adaptive integration of an oscillatory function which terminates because
 823      of roundoff error, uses the 31-pt rule */
 824
 825   {
 826     int status = 0; struct counter_params p;
 827     double result = 0, abserr=0;
 828
 829     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
 830
 831     double exp_result = -7.238969575482959717E-01;
 832     double exp_abserr =  1.285805464427459261E-14;
 833     int exp_neval   =     31;
 834     int exp_ier     =     GSL_EROUND;
 835     int exp_last    =     1;
 836
 837     double alpha = 1.3 ;
 838     gsl_function f = make_function(&f3, &alpha);
 839
 840     gsl_function fc = make_counter(&f, &p) ;
 841
 842     status = gsl_integration_qag (&fc, 0.3, 2.71, 1e-14, 0.0, w->limit,
 843                                   GSL_INTEG_GAUSS31, w,
 844                                   &result, &abserr) ;
 845
 846     gsl_test_rel(result,exp_result,1e-15,"qag(f3,31pt) oscill result");
 847     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f3,31pt) oscill abserr");
 848     gsl_test_int((int)(p.neval),exp_neval,"qag(f3,31pt) oscill neval") ;
 849     gsl_test_int((int)(w->size),exp_last,"qag(f3,31pt) oscill last") ;
 850     gsl_test_int(status,exp_ier,"qag(f3,31pt) oscill status") ;
 851
 852     p.neval = 0;
 853     status = gsl_integration_qag (&fc, 2.71, 0.3, 1e-14, 0.0, w->limit,
 854                                   GSL_INTEG_GAUSS31, w,
 855                                   &result, &abserr) ;
 856
 857     gsl_test_rel(result,-exp_result,1e-15,"qag(f3,31pt) reverse result");
 858     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f3,31pt) reverse abserr");
 859     gsl_test_int((int)(p.neval),exp_neval,"qag(f3,31pt) reverse neval") ;
 860     gsl_test_int((int)(w->size),exp_last,"qag(f3,31pt) reverse last") ;
 861     gsl_test_int(status,exp_ier,"qag(f3,31pt) reverse status") ;
 862
 863     gsl_integration_workspace_free (w) ;
 864
 865   }
 866
 867   /* Check the singularity detection (singularity at x=-0.1 in this example) */
 868
 869   {
 870     int status = 0; struct counter_params p;
 871     double result = 0, abserr=0;
 872
 873     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
 874
 875     int exp_neval  =     5151;
 876     int exp_ier    =     GSL_ESING;
 877     int exp_last   =     51;
 878
 879     double alpha = 2.0 ;
 880     gsl_function f = make_function(&f16, &alpha);
 881
 882     gsl_function fc = make_counter(&f, &p) ;
 883
 884     status = gsl_integration_qag (&fc, -1.0, 1.0, 1e-14, 0.0, w->limit,
 885                                   GSL_INTEG_GAUSS51, w,
 886                                   &result, &abserr) ;
 887
 888     gsl_test_int((int)(p.neval),exp_neval,"qag(f16,51pt) sing neval") ;
 889     gsl_test_int((int)(w->size),exp_last,"qag(f16,51pt) sing last") ;
 890     gsl_test_int(status,exp_ier,"qag(f16,51pt) sing status") ;
 891
 892     p.neval = 0;
 893     status = gsl_integration_qag (&fc, 1.0, -1.0, 1e-14, 0.0, w->limit,
 894                                   GSL_INTEG_GAUSS51, w,
 895                                   &result, &abserr) ;
 896
 897     gsl_test_int((int)(p.neval),exp_neval,"qag(f16,51pt) rev neval") ;
 898     gsl_test_int((int)(w->size),exp_last,"qag(f16,51pt) rev last") ;
 899     gsl_test_int(status,exp_ier,"qag(f16,51pt) rev status") ;
 900
 901     gsl_integration_workspace_free (w) ;
 902
 903   }
 904
 905   /* Check for hitting the iteration limit */
 906
 907   {
 908     int status = 0, i; struct counter_params p;
 909     double result = 0, abserr=0;
 910
 911     gsl_integration_workspace * w = gsl_integration_workspace_alloc (3) ;
 912
 913     double exp_result =  9.565151449233894709 ;
 914     double exp_abserr =  1.570369823891028460E+01;
 915     int exp_neval  =     305;
 916     int exp_ier    =     GSL_EMAXITER;
 917     int exp_last   =     3;
 918
 919     double a[3] = { -5.000000000000000000E-01,
 920                     0.000000000000000000,
 921                     -1.000000000000000000 } ;
 922
 923     double b[3] = { 0.000000000000000000,
 924                     1.000000000000000000,
 925                     -5.000000000000000000E-01 } ;
 926
 927     double r[3] = { 9.460353469435913709,
 928                     9.090909090909091161E-02,
 929                     1.388888888888888812E-02 } ;
 930
 931     double e[3] = { 1.570369823891028460E+01,
 932                     1.009293658750142399E-15,
 933                     1.541976423090495140E-16 } ;
 934
 935     int order[3] = { 1, 2, 3 } ;
 936
 937     double alpha = 1.0 ;
 938     gsl_function f = make_function(&f16, &alpha);
 939     gsl_function fc = make_counter(&f, &p) ;
 940
 941     status = gsl_integration_qag (&fc, -1.0, 1.0, 1e-14, 0.0, w->limit,
 942                                   GSL_INTEG_GAUSS61, w,
 943                                   &result, &abserr) ;
 944
 945     gsl_test_rel(result,exp_result,1e-15,"qag(f16,61pt) limit result") ;
 946     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f16,61pt) limit abserr") ;
 947     gsl_test_int((int)(p.neval),exp_neval,"qag(f16,61pt) limit neval") ;
 948     gsl_test_int((int)(w->size),exp_last,"qag(f16,61pt) limit last") ;
 949     gsl_test_int(status,exp_ier,"qag(f16,61pt) limit status") ;
 950
 951     for (i = 0; i < 3 ; i++)
 952         gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f16,61pt) limit alist") ;
 953
 954     for (i = 0; i < 3 ; i++)
 955         gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f16,61pt) limit blist") ;
 956
 957     for (i = 0; i < 3 ; i++)
 958         gsl_test_rel(w->rlist[i],r[i],1e-15,"qag(f16,61pt) limit rlist") ;
 959
 960     for (i = 0; i < 3 ; i++)
 961         gsl_test_rel(w->elist[i],e[i],1e-6,"qag(f16,61pt) limit elist") ;
 962
 963     for (i = 0; i < 3 ; i++)
 964         gsl_test_int((int)w->order[i],order[i]-1,"qag(f16,61pt) limit order");
 965
 966     p.neval = 0;
 967     status = gsl_integration_qag (&fc, 1.0, -1.0, 1e-14, 0.0, w->limit,
 968                                   GSL_INTEG_GAUSS61, w,
 969                                   &result, &abserr) ;
 970
 971     gsl_test_rel(result,-exp_result,1e-15,"qag(f16,61pt) reverse result") ;
 972     gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f16,61pt) reverse abserr") ;
 973     gsl_test_int((int)(p.neval),exp_neval,"qag(f16,61pt) reverse neval") ;
 974     gsl_test_int((int)(w->size),exp_last,"qag(f16,61pt) reverse last") ;
 975     gsl_test_int(status,exp_ier,"qag(f16,61pt) reverse status") ;
 976
 977     gsl_integration_workspace_free (w) ;
 978
 979   }
 980
 981   /* Test the adaptive integrator with extrapolation QAGS */
 982
 983   {
 984     int status = 0, i; struct counter_params p;
 985     double result = 0, abserr=0;
 986
 987     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
 988
 989     double exp_result = 7.716049382715789440E-02 ;
 990     double exp_abserr = 2.216394961010438404E-12 ;
 991     int exp_neval  =     189;
 992     int exp_ier    =       0;
 993     int exp_last   =       5;
 994
 995     double a[5] = { 0, 0.5, 0.25, 0.125, 0.0625 } ;
 996     double b[5] = { 0.0625, 1, 0.5, 0.25, 0.125 } ;
 997     double r[5] = { 3.919381915366914693E-05,
 998                     5.491842501998223103E-02,
 999                     1.909827770934243579E-02,
1000                     2.776531175604360097E-03,
1001                     3.280661030752062609E-04 } ;
1002     double e[5] = { 2.215538742580964735E-12,
1003                     6.097169993333454062E-16,
1004                     2.120334764359736441E-16,
1005                     3.082568839745514608E-17,
1006                     3.642265412331439511E-18 } ;
1007     int order[5] = { 1, 2, 3, 4, 5 } ;
1008
1009     double alpha = 2.6 ;
1010     gsl_function f = make_function(&f1, &alpha);
1011     gsl_function fc = make_counter(&f, &p) ;
1012
1013     status = gsl_integration_qags (&fc, 0.0, 1.0, 0.0, 1e-10, w->limit,
1014                                    w,
1015                                    &result, &abserr) ;
1016
1017     gsl_test_rel(result,exp_result,1e-15,"qags(f1) smooth result") ;
1018     gsl_test_rel(abserr,exp_abserr,1e-6,"qags(f1) smooth abserr") ;
1019     gsl_test_int((int)(p.neval),exp_neval,"qags(f1) smooth neval") ;
1020     gsl_test_int((int)(w->size),exp_last,"qags(f1) smooth last") ;
1021     gsl_test_int(status,exp_ier,"qags(f1) smooth status") ;
1022
1023     for (i = 0; i < 5 ; i++)
1024         gsl_test_rel(w->alist[i],a[i],1e-15,"qags(f1) smooth alist") ;
1025
1026     for (i = 0; i < 5 ; i++)
1027         gsl_test_rel(w->blist[i],b[i],1e-15,"qags(f1) smooth blist") ;
1028
1029     for (i = 0; i < 5 ; i++)
1030         gsl_test_rel(w->rlist[i],r[i],1e-15,"qags(f1) smooth rlist") ;
1031
1032     for (i = 0; i < 5 ; i++)
1033         gsl_test_rel(w->elist[i],e[i],1e-6,"qags(f1) smooth elist") ;
1034
1035     for (i = 0; i < 5 ; i++)
1036         gsl_test_int((int)w->order[i],order[i]-1,"qags(f1) smooth order") ;
1037
1038     p.neval = 0;
1039     status = gsl_integration_qags (&fc, 1.0, 0.0, 0.0, 1e-10, w->limit,
1040                                    w,
1041                                    &result, &abserr) ;
1042
1043     gsl_test_rel(result,-exp_result,1e-15,"qags(f1) reverse result") ;
1044     gsl_test_rel(abserr,exp_abserr,1e-6,"qags(f1) reverse abserr") ;
1045     gsl_test_int((int)(p.neval),exp_neval,"qags(f1) reverse neval") ;
1046     gsl_test_int((int)(w->size),exp_last,"qags(f1) reverse last") ;
1047     gsl_test_int(status,exp_ier,"qags(f1) reverse status") ;
1048
1049     gsl_integration_workspace_free (w) ;
1050
1051   }
1052
1053   /* Test f11 using an absolute error bound */
1054
1055   {
1056     int status = 0, i; struct counter_params p;
1057     double result = 0, abserr=0;
1058
1059     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1060
1061     /* All results are for GSL_IEEE_MODE=double-precision */
1062
1063     double exp_result = -5.908755278982136588E+03 ;
1064     double exp_abserr = 1.299646281053874554E-10 ;
1065     int exp_neval  =     357;
1066     int exp_ier    =       0;
1067     int exp_last   =       9;
1068
1069     double a[9] = { 1.000000000000000000E+00,
1070                     5.005000000000000000E+02,
1071                     2.507500000000000000E+02,
1072                     1.258750000000000000E+02,
1073                     6.343750000000000000E+01,
1074                     3.221875000000000000E+01,
1075                     1.660937500000000000E+01,
1076                     8.804687500000000000E+00,
1077                     4.902343750000000000E+00 } ;
1078     double b[9] = { 4.902343750000000000E+00,
1079                     1.000000000000000000E+03,
1080                     5.005000000000000000E+02,
1081                     2.507500000000000000E+02,
1082                     1.258750000000000000E+02,
1083                     6.343750000000000000E+01,
1084                     3.221875000000000000E+01,
1085                     1.660937500000000000E+01,
1086                     8.804687500000000000E+00 } ;
1087     double r[9] = { -3.890977835520834649E+00,
1088                     -3.297343675805121620E+03,
1089                     -1.475904154146372775E+03,
1090                     -6.517404019686431411E+02,
1091                     -2.829354222635842007E+02,
1092                     -1.201692001973227519E+02,
1093                     -4.959999906099650246E+01,
1094                     -1.971441499411640308E+01,
1095                     -7.457032710459004399E+00 } ;
1096     double e[9] = { 6.448276035006137169E-11,
1097                     3.660786868980994028E-11,
1098                     1.638582774073219226E-11,
1099                     7.235772003440423011E-12,
1100                     3.141214202790722909E-12,
1101                     1.334146129098576244E-12,
1102                     5.506706097890446534E-13,
1103                     2.188739744348345039E-13,
1104                     8.278969410534525339E-14 } ;
1105     int order[9] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 } ;
1106
1107     double alpha = 2.0 ;
1108     gsl_function f = make_function(&f11, &alpha);
1109     gsl_function fc = make_counter(&f, &p) ;
1110
1111     status = gsl_integration_qags (&fc, 1.0, 1000.0, 1e-7, 0.0, w->limit,
1112                                    w,
1113                                    &result, &abserr) ;
1114
1115     gsl_test_rel(result,exp_result,1e-15,"qags(f11) smooth result") ;
1116     gsl_test_rel(abserr,exp_abserr,1e-3,"qags(f11) smooth abserr") ;
1117     gsl_test_int((int)(p.neval),exp_neval,"qags(f11) smooth neval") ;
1118     gsl_test_int((int)(w->size),exp_last,"qags(f11) smooth last") ;
1119     gsl_test_int(status,exp_ier,"qags(f11) smooth status") ;
1120
1121     for (i = 0; i < 9 ; i++)
1122         gsl_test_rel(w->alist[i],a[i],1e-15,"qags(f11) smooth alist") ;
1123
1124     for (i = 0; i < 9 ; i++)
1125         gsl_test_rel(w->blist[i],b[i],1e-15,"qags(f11) smooth blist") ;
1126
1127     for (i = 0; i < 9 ; i++)
1128         gsl_test_rel(w->rlist[i],r[i],1e-15,"qags(f11) smooth rlist") ;
1129
1130     for (i = 0; i < 9 ; i++)
1131         gsl_test_rel(w->elist[i],e[i],1e-5,"qags(f11) smooth elist") ;
1132
1133     for (i = 0; i < 9 ; i++)
1134         gsl_test_int((int)w->order[i],order[i]-1,"qags(f11) smooth order");
1135
1136     p.neval = 0;
1137     status = gsl_integration_qags (&fc, 1000.0, 1.0, 1e-7, 0.0, w->limit,
1138                                    w,
1139                                    &result, &abserr) ;
1140
1141     gsl_test_rel(result,-exp_result,1e-15,"qags(f11) reverse result") ;
1142     gsl_test_rel(abserr,exp_abserr,1e-3,"qags(f11) reverse abserr") ;
1143     gsl_test_int((int)(p.neval),exp_neval,"qags(f11) reverse neval") ;
1144     gsl_test_int((int)(w->size),exp_last,"qags(f11) reverse last") ;
1145     gsl_test_int(status,exp_ier,"qags(f11) reverse status") ;
1146
1147     gsl_integration_workspace_free (w) ;
1148
1149   }
1150
1151   /* Test infinite range integral f455 using a relative error bound */
1152
1153   {
1154     int status = 0, i; struct counter_params p;
1155     double result = 0, abserr=0;
1156
1157     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1158
1159     /* All results are for GSL_IEEE_MODE=double-precision */
1160
1161     double exp_result = -3.616892186127022568E-01 ;
1162     double exp_abserr = 3.016716913328831851E-06;
1163     int exp_neval  =      285;
1164     int exp_ier    =        0;
1165     int exp_last   =       10;
1166
1167     double a[10] = { 9.687500000000000000E-01,
1168                      0.000000000000000000E+00,
1169                      5.000000000000000000E-01,
1170                      2.500000000000000000E-01,
1171                      7.500000000000000000E-01,
1172                      1.250000000000000000E-01,
1173                      8.750000000000000000E-01,
1174                      6.250000000000000000E-02,
1175                      9.375000000000000000E-01,
1176                      3.125000000000000000E-02 } ;
1177     double b[10] = { 1.000000000000000000E+00,
1178                      3.125000000000000000E-02,
1179                      7.500000000000000000E-01,
1180                      5.000000000000000000E-01,
1181                      8.750000000000000000E-01,
1182                      2.500000000000000000E-01,
1183                      9.375000000000000000E-01,
1184                      1.250000000000000000E-01,
1185                      9.687500000000000000E-01,
1186                      6.250000000000000000E-02 } ;
1187     double r[10] = { -1.390003415539725340E-01,
1188                      1.429785306003466313E-03,
1189                      -1.229943369113085765E-02,
1190                      2.995321156568048898E-03,
1191                      -4.980050133751051655E-02,
1192                      2.785385934678596704E-03,
1193                      -8.653752279614615461E-02,
1194                      1.736218164975512294E-03,
1195                      -8.398745675010892142E-02,
1196                      1.041689192004495576E-03 } ;
1197     double e[10] = { 2.395037249893453013E-02,
1198                      2.161214992172538524E-04,
1199                      5.720644840858777846E-14,
1200                      3.325474514168701167E-17,
1201                      3.147380432198176412E-14,
1202                      3.092399597147240624E-17,
1203                      9.607595030230581153E-16,
1204                      1.927589382528252344E-17,
1205                      9.324480826368044019E-16,
1206                      1.156507325466566521E-17 } ;
1207     int order[10] = { 1, 2, 3, 5, 7, 9, 4, 6, 8, 10 } ;
1208
1209     gsl_function f = make_function(&f455, 0);
1210     gsl_function fc = make_counter(&f, &p) ;
1211
1212     status = gsl_integration_qagiu (&fc, 0.0, 0.0, 1.0e-3, w->limit,
1213                                     w,
1214                                     &result, &abserr) ;
1215
1216     gsl_test_rel(result,exp_result,1e-14,"qagiu(f455) smooth result") ;
1217     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(f455) smooth abserr") ;
1218     gsl_test_int((int)(p.neval),exp_neval,"qagiu(f455) smooth neval") ;
1219     gsl_test_int((int)(w->size),exp_last,"qagiu(f455) smooth last") ;
1220     gsl_test_int(status,exp_ier,"qagiu(f455) smooth status") ;
1221
1222     for (i = 0; i < 10 ; i++)
1223         gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(f455) smooth alist") ;
1224
1225     for (i = 0; i < 10 ; i++)
1226         gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(f455) smooth blist") ;
1227
1228     for (i = 0; i < 10 ; i++)
1229         gsl_test_rel(w->rlist[i],r[i],1e-15,"qagiu(f455) smooth rlist") ;
1230
1231     for (i = 0; i < 10 ; i++)
1232         gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(f455) smooth elist") ;
1233
1234     for (i = 0; i < 10 ; i++)
1235         gsl_test_int((int)w->order[i],order[i]-1,"qagiu(f455) smooth order");
1236
1237     gsl_integration_workspace_free (w) ;
1238
1239   }
1240
1241   /* Test infinite range integral f15 using a relative error bound */
1242
1243   {
1244     int status = 0, i; struct counter_params p;
1245     double result = 0, abserr=0;
1246
1247     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1248
1249     /* All results are for GSL_IEEE_MODE=double-precision */
1250
1251     double exp_result = 6.553600000000024738E+04;
1252     double exp_abserr = 7.121667111456009280E-04;
1253     int exp_neval  =      285;
1254     int exp_ier    =        0;
1255     int exp_last   =       10;
1256
1257     double a[10] = { 0.000000000000000000E+00,
1258                      5.000000000000000000E-01,
1259                      2.500000000000000000E-01,
1260                      1.250000000000000000E-01,
1261                      6.250000000000000000E-02,
1262                      3.125000000000000000E-02,
1263                      1.562500000000000000E-02,
1264                      7.812500000000000000E-03,
1265                      3.906250000000000000E-03,
1266                      1.953125000000000000E-03 } ;
1267     double b[10] = { 1.953125000000000000E-03,
1268                      1.000000000000000000E+00,
1269                      5.000000000000000000E-01,
1270                      2.500000000000000000E-01,
1271                      1.250000000000000000E-01,
1272                      6.250000000000000000E-02,
1273                      3.125000000000000000E-02,
1274                      1.562500000000000000E-02,
1275                      7.812500000000000000E-03,
1276                      3.906250000000000000E-03 } ;
1277     double r[10] = { 1.099297665754340292E+00,
1278                      3.256176475185617591E-01,
1279                      8.064694554185326325E+00,
1280                      8.873128656118993263E+01,
1281                      6.977679035845269482E+02,
1282                      4.096981198511257389E+03,
1283                      1.574317583220441520E+04,
1284                      2.899418134793237914E+04,
1285                      1.498314766425578091E+04,
1286                      9.225251570832365360E+02 } ;
1287     double e[10] = { 7.101865971621337814E-04,
1288                      1.912660677170175771E-08,
1289                      9.167763417119923333E-08,
1290                      3.769501719163865578E-07,
1291                      6.973493131275552509E-07,
1292                      1.205653952340679711E-07,
1293                      1.380003928453846583E-07,
1294                      1.934652413547325474E-07,
1295                      3.408933028357320364E-07,
1296                      2.132473175465897029E-09 } ;
1297     int order[10] = { 1, 5, 4, 9, 8, 7, 6, 3, 2, 10 } ;
1298
1299     double alpha = 5.0;
1300
1301     gsl_function f = make_function(&f15, &alpha);
1302     gsl_function fc = make_counter(&f, &p) ;
1303
1304     status = gsl_integration_qagiu (&fc, 0.0, 0.0, 1.0e-7, w->limit,
1305                                     w,
1306                                     &result, &abserr) ;
1307
1308     gsl_test_rel(result,exp_result,1e-14,"qagiu(f15) smooth result") ;
1309     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(f15) smooth abserr") ;
1310     gsl_test_int((int)(p.neval),exp_neval,"qagiu(f15) smooth neval") ;
1311     gsl_test_int((int)(w->size),exp_last,"qagiu(f15) smooth last") ;
1312     gsl_test_int(status,exp_ier,"qagiu(f15) smooth status") ;
1313
1314     for (i = 0; i < 10 ; i++)
1315         gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(f15) smooth alist") ;
1316
1317     for (i = 0; i < 10 ; i++)
1318         gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(f15) smooth blist") ;
1319
1320     for (i = 0; i < 10 ; i++)
1321         gsl_test_rel(w->rlist[i],r[i],1e-15,"qagiu(f15) smooth rlist") ;
1322
1323     for (i = 0; i < 10 ; i++)
1324         gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(f15) smooth elist") ;
1325
1326     for (i = 0; i < 10 ; i++)
1327         gsl_test_int((int)w->order[i],order[i]-1,"qagiu(f15) smooth order");
1328
1329     gsl_integration_workspace_free (w) ;
1330
1331   }
1332
1333   /* Test infinite range integral f16 using an absolute error bound */
1334
1335   {
1336     int status = 0, i; struct counter_params p;
1337     double result = 0, abserr=0;
1338
1339     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1340
1341     /* All results are for GSL_IEEE_MODE=double-precision */
1342
1343     double exp_result = 1.000000000006713292E-04;
1344     double exp_abserr = 3.084062020905636316E-09;
1345     int exp_neval  =      165;
1346     int exp_ier    =        0;
1347     int exp_last   =        6;
1348
1349     double a[6] = { 0.000000000000000000E+00,
1350                     5.000000000000000000E-01,
1351                     2.500000000000000000E-01,
1352                     1.250000000000000000E-01,
1353                     6.250000000000000000E-02,
1354                     3.125000000000000000E-02 } ;
1355     double b[6] = { 3.125000000000000000E-02,
1356                     1.000000000000000000E+00,
1357                     5.000000000000000000E-01,
1358                     2.500000000000000000E-01,
1359                     1.250000000000000000E-01,
1360                     6.250000000000000000E-02 } ;
1361     double r[6] = { 7.633587786326674618E-05,
1362                     9.900990099009899620E-07,
1363                     1.922522349322310737E-06,
1364                     3.629434715543053753E-06,
1365                     6.501422186103209199E-06,
1366                     1.062064387653501389E-05 } ;
1367     double e[6] = { 3.084061858351569051E-09,
1368                     3.112064814755089674E-17,
1369                     4.543453652226561245E-17,
1370                     4.908618166361344548E-17,
1371                     3.014338672269481784E-17,
1372                     6.795996738013555461E-18 } ;
1373     int order[6] = { 1, 4, 3, 2, 5, 6 } ;
1374
1375     double alpha = 1.0;
1376
1377     gsl_function f = make_function(&f16, &alpha);
1378     gsl_function fc = make_counter(&f, &p) ;
1379
1380     status = gsl_integration_qagiu (&fc, 99.9, 1.0e-7, 0.0, w->limit,
1381                                     w,
1382                                     &result, &abserr) ;
1383
1384     gsl_test_rel(result,exp_result,1e-14,"qagiu(f16) smooth result") ;
1385     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(f16) smooth abserr") ;
1386     gsl_test_int((int)(p.neval),exp_neval,"qagiu(f16) smooth neval") ;
1387     gsl_test_int((int)(w->size),exp_last,"qagiu(f16) smooth last") ;
1388     gsl_test_int(status,exp_ier,"qagiu(f16) smooth status") ;
1389
1390     for (i = 0; i < 6 ; i++)
1391         gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(f16) smooth alist") ;
1392
1393     for (i = 0; i < 6 ; i++)
1394         gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(f16) smooth blist") ;
1395
1396     for (i = 0; i < 6 ; i++)
1397         gsl_test_rel(w->rlist[i],r[i],1e-15,"qagiu(f16) smooth rlist") ;
1398
1399     for (i = 0; i < 6 ; i++)
1400         gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(f16) smooth elist") ;
1401
1402     for (i = 0; i < 6 ; i++)
1403         gsl_test_int((int)w->order[i],order[i]-1,"qagiu(f16) smooth order");
1404
1405     gsl_integration_workspace_free (w) ;
1406
1407   }
1408
1409   /* Test infinite range integral myfn1 using an absolute error bound */
1410
1411   {
1412     int status = 0, i; struct counter_params p;
1413     double result = 0, abserr=0;
1414
1415     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1416
1417     /* All results are for GSL_IEEE_MODE=double-precision */
1418
1419     double exp_result = 2.275875794468747770E+00;
1420     double exp_abserr = 7.436490118267390744E-09;
1421     int exp_neval  =      270;
1422     int exp_ier    =        0;
1423     int exp_last   =        5;
1424
1425     double a[5] = { 1.250000000000000000E-01,
1426                     5.000000000000000000E-01,
1427                     2.500000000000000000E-01,
1428                     0.000000000000000000E+00,
1429                     3.750000000000000000E-01 } ;
1430     double b[5] = { 2.500000000000000000E-01,
1431                     1.000000000000000000E+00,
1432                     3.750000000000000000E-01,
1433                     1.250000000000000000E-01,
1434                     5.000000000000000000E-01 } ;
1435     double r[5] = { 4.639317228058405717E-04,
1436                     1.691664195356748834E+00,
1437                     1.146307471900291086E-01,
1438                     4.379392477350953574E-20,
1439                     4.691169201991640669E-01 } ;
1440     double e[5] = { 3.169263960393051137E-09,
1441                     4.265988974874425043E-09,
1442                     1.231954072964969637E-12,
1443                     8.360902986775307673E-20,
1444                     5.208244060463541433E-15 } ;
1445     int order[5] = { 2, 1, 3, 5, 4 } ;
1446
1447     gsl_function f = make_function(&myfn1, 0);
1448     gsl_function fc = make_counter(&f, &p) ;
1449
1450     status = gsl_integration_qagi (&fc, 1.0e-7, 0.0, w->limit,
1451                                    w,
1452                                    &result, &abserr) ;
1453
1454     gsl_test_rel(result,exp_result,1e-14,"qagiu(myfn1) smooth result") ;
1455     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(myfn1) smooth abserr") ;
1456     gsl_test_int((int)(p.neval),exp_neval,"qagiu(myfn1) smooth neval") ;
1457     gsl_test_int((int)(w->size),exp_last,"qagiu(myfn1) smooth last") ;
1458     gsl_test_int(status,exp_ier,"qagiu(myfn1) smooth status") ;
1459
1460     for (i = 0; i < 5 ; i++)
1461         gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(myfn1) smooth alist") ;
1462
1463     for (i = 0; i < 5 ; i++)
1464         gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(myfn1) smooth blist") ;
1465
1466     for (i = 0; i < 5 ; i++)
1467         gsl_test_rel(w->rlist[i],r[i],1e-14,"qagiu(myfn1) smooth rlist") ;
1468
1469     for (i = 0; i < 5 ; i++)
1470         gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(myfn1) smooth elist") ;
1471
1472     for (i = 0; i < 5 ; i++)
1473         gsl_test_int((int)w->order[i],order[i]-1,"qagiu(myfn1) smooth order");
1474
1475     gsl_integration_workspace_free (w) ;
1476
1477   }
1478
1479   /* Test infinite range integral myfn1 using an absolute error bound */
1480
1481   {
1482     int status = 0, i; struct counter_params p;
1483     double result = 0, abserr=0;
1484
1485     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1486
1487     /* All results are for GSL_IEEE_MODE=double-precision */
1488
1489     double exp_result = 2.718281828459044647E+00;
1490     double exp_abserr = 1.588185109253204805E-10;
1491     int exp_neval  =      135;
1492     int exp_ier    =        0;
1493     int exp_last   =        5;
1494
1495     double a[5] = { 0.000000000000000000E+00,
1496                     5.000000000000000000E-01,
1497                     2.500000000000000000E-01,
1498                     1.250000000000000000E-01,
1499                     6.250000000000000000E-02 } ;
1500     double b[5] = { 6.250000000000000000E-02,
1501                     1.000000000000000000E+00,
1502                     5.000000000000000000E-01,
1503                     2.500000000000000000E-01,
1504                     1.250000000000000000E-01 } ;
1505     double r[5] = { 8.315287189746029816E-07,
1506                     1.718281828459045091E+00,
1507                     8.646647167633871867E-01,
1508                     1.328565310599463256E-01,
1509                     2.477920647947255521E-03 } ;
1510     double e[5] = { 1.533437090413525935E-10,
1511                     4.117868247943567505E-12,
1512                     7.802455785301941044E-13,
1513                     5.395586026138397182E-13,
1514                     3.713312434866150125E-14 } ;
1515     int order[5] = { 1, 2, 3, 4, 5 } ;
1516
1517     double alpha = 1.0 ;
1518     gsl_function f = make_function(&myfn2, &alpha);
1519     gsl_function fc = make_counter(&f, &p) ;
1520
1521     status = gsl_integration_qagil (&fc, 1.0, 1.0e-7, 0.0, w->limit,
1522                                     w,
1523                                     &result, &abserr) ;
1524
1525     gsl_test_rel(result,exp_result,1e-14,"qagiu(myfn2) smooth result") ;
1526     gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(myfn2) smooth abserr") ;
1527     gsl_test_int((int)(p.neval),exp_neval,"qagiu(myfn2) smooth neval") ;
1528     gsl_test_int((int)(w->size),exp_last,"qagiu(myfn2) smooth last") ;
1529     gsl_test_int(status,exp_ier,"qagiu(myfn2) smooth status") ;
1530
1531     for (i = 0; i < 5 ; i++)
1532         gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(myfn2) smooth alist") ;
1533
1534     for (i = 0; i < 5 ; i++)
1535         gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(myfn2) smooth blist") ;
1536
1537     for (i = 0; i < 5 ; i++)
1538         gsl_test_rel(w->rlist[i],r[i],1e-14,"qagiu(myfn2) smooth rlist") ;
1539
1540     for (i = 0; i < 5 ; i++)
1541         gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(myfn2) smooth elist") ;
1542
1543     for (i = 0; i < 5 ; i++)
1544         gsl_test_int((int)w->order[i],order[i]-1,"qagiu(myfn2) smooth order");
1545
1546     gsl_integration_workspace_free (w) ;
1547
1548   }
1549
1550   /* Test integral f454 with integrable singular points */
1551
1552   {
1553     int status = 0, i; struct counter_params p;
1554     double result = 0, abserr=0;
1555
1556     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1557
1558     /* All results are for GSL_IEEE_MODE=double-precision */
1559
1560     double exp_result = 5.274080611672716401E+01;
1561     double exp_abserr = 1.755703848687062418E-04;
1562     int exp_neval  =        777;
1563     int exp_ier    =          0;
1564     int exp_last   =         20;
1565
1566     double a[20] = { 9.687500000000000000E-01,
1567                     1.401269388548935790E+00,
1568                     1.414213562373095145E+00,
1569                     1.000000000000000000E+00,
1570                     0.000000000000000000E+00,
1571                     2.207106781186547462E+00,
1572                     1.810660171779821415E+00,
1573                     1.207106781186547462E+00,
1574                     5.000000000000000000E-01,
1575                     1.103553390593273731E+00,
1576                     1.612436867076458391E+00,
1577                     1.310660171779821415E+00,
1578                     7.500000000000000000E-01,
1579                     1.051776695296636976E+00,
1580                     1.513325214724776657E+00,
1581                     1.362436867076458391E+00,
1582                     8.750000000000000000E-01,
1583                     1.463769388548935790E+00,
1584                     1.388325214724776657E+00,
1585                     9.375000000000000000E-01} ;
1586     double b[20] = { 1.000000000000000000E+00,
1587                      1.414213562373095145E+00,
1588                      1.463769388548935790E+00,
1589                      1.051776695296636976E+00,
1590                      5.000000000000000000E-01,
1591                      3.000000000000000000E+00,
1592                      2.207106781186547462E+00,
1593                      1.310660171779821415E+00,
1594                      7.500000000000000000E-01,
1595                      1.207106781186547462E+00,
1596                      1.810660171779821415E+00,
1597                      1.362436867076458391E+00,
1598                      8.750000000000000000E-01,
1599                      1.103553390593273731E+00,
1600                      1.612436867076458391E+00,
1601                      1.388325214724776657E+00,
1602                      9.375000000000000000E-01,
1603                      1.513325214724776657E+00,
1604                      1.401269388548935790E+00,
1605                      9.687500000000000000E-01} ;
1606     double r[20] = { -1.125078814079027711E-01,
1607                      -1.565132123531515207E-01,
1608                      -4.225328513207429193E-01,
1609                      -1.830392049835374568E-01,
1610                      6.575875041899758092E-03,
1611                      4.873920540843067783E+01,
1612                      6.032891565603589079E+00,
1613                      -2.991531901645863023E-01,
1614                      -7.326282608704996063E-03,
1615                      -2.431894410706912923E-01,
1616                      5.911661670635662835E-01,
1617                      -2.236786562536174916E-01,
1618                      -5.647871991778510847E-02,
1619                      -1.305470403178642658E-01,
1620                      -1.721363984401322045E-01,
1621                      -1.589345454585119055E-01,
1622                      -7.406626263352669715E-02,
1623                      -2.208730668000830344E-01,
1624                      -1.048692749517999567E-01,
1625                      -6.302287584527696551E-02} ;
1626     double e[20] = { 2.506431410088378817E-02,
1627                      2.730454695485963826E-02,
1628                      1.017446081816190118E-01,
1629                      3.252808038935910834E-02,
1630                      7.300687878575027348E-17,
1631                      5.411138804637469780E-13,
1632                      6.697855121200013106E-14,
1633                      3.321267596107916554E-15,
1634                      1.417509685426979386E-16,
1635                      2.699945168224041491E-15,
1636                      6.573952690524728748E-15,
1637                      2.483331942899818875E-15,
1638                      6.270397525408045936E-16,
1639                      1.449363299575615261E-15,
1640                      1.911097929242846383E-15,
1641                      1.764527917763735212E-15,
1642                      8.223007012367522077E-16,
1643                      2.452183642810224359E-15,
1644                      1.164282836272345215E-15,
1645                      6.996944784151910810E-16} ;
1646     int order[20] = { 3, 4, 2, 1, 6, 7, 11, 8, 10, 12, 18,
1647                      15, 16, 14, 19, 17, 20, 13, 9, 5 } ;
1648
1649     gsl_function f = make_function(&f454, 0);
1650     gsl_function fc = make_counter(&f, &p) ;
1651
1652     double pts[4] ;
1653
1654     pts[0] = 0.0;
1655     pts[1] = 1.0;
1656     pts[2] = sqrt(2.0);
1657     pts[3] = 3.0;
1658
1659     status = gsl_integration_qagp (&fc, pts, 4,
1660                                    0.0, 1.0e-3, w->limit,
1661                                    w,
1662                                    &result, &abserr) ;
1663
1664     gsl_test_rel(result,exp_result,1e-14,"qagp(f454) singular result") ;
1665     gsl_test_rel(abserr,exp_abserr,1e-5,"qagp(f454) singular abserr") ;
1666     gsl_test_int((int)(p.neval),exp_neval,"qagp(f454) singular neval") ;
1667     gsl_test_int((int)(w->size),exp_last,"qagp(f454) singular last") ;
1668     gsl_test_int(status,exp_ier,"qagp(f454) singular status") ;
1669
1670     for (i = 0; i < 20 ; i++)
1671         gsl_test_rel(w->alist[i],a[i],1e-15,"qagp(f454) singular alist") ;
1672
1673     for (i = 0; i < 20 ; i++)
1674         gsl_test_rel(w->blist[i],b[i],1e-15,"qagp(f454) singular blist") ;
1675
1676     for (i = 0; i < 20 ; i++)
1677         gsl_test_rel(w->rlist[i],r[i],1e-14,"qagp(f454) singular rlist") ;
1678
1679     for (i = 0; i < 20 ; i++)
1680         gsl_test_rel(w->elist[i],e[i],1e-4,"qagp(f454) singular elist") ;
1681
1682     for (i = 0; i < 20 ; i++)
1683         gsl_test_int((int)w->order[i],order[i]-1,"qagp(f454) singular order");
1684
1685     gsl_integration_workspace_free (w) ;
1686
1687   }
1688
1689
1690   /* Test cauchy integration using a relative error bound */
1691
1692   {
1693     int status = 0, i; struct counter_params p;
1694     double result = 0, abserr=0;
1695
1696     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1697
1698     /* All results are for GSL_IEEE_MODE=double-precision */
1699
1700     double exp_result = -8.994400695837000137E-02;
1701     double exp_abserr =  1.185290176227023727E-06;
1702     int exp_neval  =      215;
1703     int exp_ier    =        0;
1704     int exp_last   =        6;
1705
1706     double a[6] = { -1.000000000000000000E+00,
1707                     2.500000000000000000E+00,
1708                     1.250000000000000000E+00,
1709                     6.250000000000000000E-01,
1710                     -5.000000000000000000E-01,
1711                     -7.500000000000000000E-01} ;
1712     double b[6] = { -7.500000000000000000E-01,
1713                     5.000000000000000000E+00,
1714                     2.500000000000000000E+00,
1715                     1.250000000000000000E+00,
1716                     6.250000000000000000E-01,
1717                     -5.000000000000000000E-01} ;
1718     double r[6] = { -1.234231128040012976E-01,
1719                     3.579970394639702888E-03,
1720                     2.249831615049339983E-02,
1721                     7.214232992127905808E-02,
1722                     2.079093855884046535E-02,
1723                     -8.553244917962132821E-02} ;
1724     double e[6] = { 1.172832717970022565E-06,
1725                     9.018232896137375412E-13,
1726                     1.815172652101790755E-12,
1727                     1.006998195150956048E-13,
1728                     1.245463873006391609E-08,
1729                     1.833082948207153514E-15 } ;
1730     int order[6] = { 1, 5, 3, 2, 4, 6 } ;
1731
1732     double alpha = 1.0 ;
1733     gsl_function f = make_function(&f459, &alpha);
1734     gsl_function fc = make_counter(&f, &p) ;
1735
1736     status = gsl_integration_qawc (&fc, -1.0, 5.0, 0.0, 0.0, 1.0e-3, w->limit,
1737                                    w,
1738                                    &result, &abserr) ;
1739
1740     gsl_test_rel(result,exp_result,1e-14,"qawc(f459) result") ;
1741     gsl_test_rel(abserr,exp_abserr,1e-6,"qawc(f459) abserr") ;
1742     gsl_test_int((int)(p.neval),exp_neval,"qawc(f459) neval") ;
1743     gsl_test_int((int)(w->size),exp_last,"qawc(f459) last") ;
1744     gsl_test_int(status,exp_ier,"qawc(f459) status") ;
1745
1746     for (i = 0; i < 6 ; i++)
1747         gsl_test_rel(w->alist[i],a[i],1e-15,"qawc(f459) alist") ;
1748
1749     for (i = 0; i < 6 ; i++)
1750         gsl_test_rel(w->blist[i],b[i],1e-15,"qawc(f459) blist") ;
1751
1752     for (i = 0; i < 6 ; i++)
1753         gsl_test_rel(w->rlist[i],r[i],1e-14,"qawc(f459) rlist") ;
1754
1755     for (i = 0; i < 6 ; i++)
1756         gsl_test_rel(w->elist[i],e[i],1e-5,"qawc(f459) elist") ;
1757
1758     for (i = 0; i < 6 ; i++)
1759         gsl_test_int((int)w->order[i],order[i]-1,"qawc(f459) order");
1760
1761     p.neval = 0;
1762     status = gsl_integration_qawc (&fc, 5.0, -1.0, 0.0, 0.0, 1.0e-3, w->limit,
1763                                    w,
1764                                    &result, &abserr) ;
1765
1766     gsl_test_rel(result,-exp_result,1e-14,"qawc(f459) rev result") ;
1767     gsl_test_rel(abserr,exp_abserr,1e-6,"qawc(f459) rev abserr") ;
1768     gsl_test_int((int)(p.neval),exp_neval,"qawc(f459) rev neval") ;
1769     gsl_test_int((int)(w->size),exp_last,"qawc(f459) rev last") ;
1770     gsl_test_int(status,exp_ier,"qawc(f459) rev status") ;
1771
1772     gsl_integration_workspace_free (w) ;
1773
1774   }
1775
1776   /* Test QAWS singular integration using a relative error bound */
1777
1778   {
1779     int status = 0, i; struct counter_params p;
1780     double result = 0, abserr=0;
1781
1782     gsl_integration_qaws_table * t
1783       = gsl_integration_qaws_table_alloc (0.0, 0.0, 1, 0);
1784
1785     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1786
1787     /* All results are for GSL_IEEE_MODE=double-precision */
1788
1789     double exp_result = -1.892751853489401670E-01;
1790     double exp_abserr = 1.129133712015747658E-08;
1791     int exp_neval  =      280;
1792     int exp_ier    =        0;
1793     int exp_last   =        8;
1794
1795     double a[8] = { 0.000000000000000000E+00,
1796                     5.000000000000000000E-01,
1797                     2.500000000000000000E-01,
1798                     1.250000000000000000E-01,
1799                     6.250000000000000000E-02,
1800                     3.125000000000000000E-02,
1801                     1.562500000000000000E-02,
1802                     7.812500000000000000E-03} ;
1803     double b[8] = { 7.812500000000000000E-03,
1804                     1.000000000000000000E+00,
1805                     5.000000000000000000E-01,
1806                     2.500000000000000000E-01,
1807                     1.250000000000000000E-01,
1808                     6.250000000000000000E-02,
1809                     3.125000000000000000E-02,
1810                     1.562500000000000000E-02} ;
1811     double r[8] = { -4.126317299834445824E-05,
1812                     -1.076283950172247789E-01,
1813                     -6.240573216173390947E-02,
1814                     -1.456169844189576269E-02,
1815                     -3.408925115926728436E-03,
1816                     -8.914083918175634211E-04,
1817                     -2.574191402137795482E-04,
1818                     -8.034390712936630608E-05} ;
1819     double e[8] = { 1.129099387465713953E-08,
1820                     3.423394967694403596E-13,
1821                     6.928428071454762659E-16,
1822                     1.616673288784094320E-16,
1823                     3.784667152924835070E-17,
1824                     9.896621209399419425E-18,
1825                     2.857926564445496100E-18,
1826                     8.919965558336773736E-19} ;
1827     int order[8] = { 1, 2, 3, 4, 5, 6, 7, 8 } ;
1828
1829     double alpha = 1.0 ;
1830     gsl_function f = make_function(&f458, &alpha);
1831     gsl_function fc = make_counter(&f, &p) ;
1832
1833     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
1834                                    w,
1835                                    &result, &abserr) ;
1836
1837     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) ln(x-a) result") ;
1838     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) ln(x-a) abserr") ;
1839     gsl_test_int((int)(p.neval),exp_neval,"qaws(f458) ln(x-a) neval") ;
1840     gsl_test_int((int)(w->size),exp_last,"qaws(f458) ln(x-a) last") ;
1841     gsl_test_int(status,exp_ier,"qaws(f458) ln(x-a) status") ;
1842
1843     for (i = 0; i < 6 ; i++)
1844         gsl_test_rel(w->alist[i],a[i],1e-15,"qaws(f458) ln(x-a) alist") ;
1845
1846     for (i = 0; i < 6 ; i++)
1847         gsl_test_rel(w->blist[i],b[i],1e-15,"qaws(f458) ln(x-a) blist") ;
1848
1849     for (i = 0; i < 6 ; i++)
1850         gsl_test_rel(w->rlist[i],r[i],1e-14,"qaws(f458) ln(x-a) rlist") ;
1851
1852     for (i = 0; i < 6 ; i++)
1853         gsl_test_rel(w->elist[i],e[i],1e-4,"qaws(f458) ln(x-a) elist") ;
1854
1855     for (i = 0; i < 6 ; i++)
1856         gsl_test_int((int)w->order[i],order[i]-1,"qaws(f458) ln(x-a) order");
1857
1858     /* Test without logs */
1859
1860     gsl_integration_qaws_table_set (t, -0.5, -0.3, 0, 0);
1861
1862     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
1863                                    w, &result, &abserr) ;
1864
1865     exp_result = 9.896686656601706433E-01;
1866     exp_abserr = 5.888032513201251628E-08;
1867
1868     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) AB result") ;
1869     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) AB abserr") ;
1870
1871     /* Test with ln(x - a) */
1872
1873     gsl_integration_qaws_table_set (t, -0.5, -0.3, 1, 0);
1874
1875     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
1876                                    w, &result, &abserr) ;
1877
1878     exp_result = -3.636679470586539620E-01;
1879     exp_abserr = 2.851348775257054093E-08;
1880
1881     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) AB ln(x-a) result") ;
1882     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) AB ln(x-a) abserr") ;
1883
1884     /* Test with ln(b - x) */
1885
1886     gsl_integration_qaws_table_set (t, -0.5, -0.3, 0, 1);
1887
1888     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
1889                                    w, &result, &abserr) ;
1890
1891     exp_result = -1.911489253363409802E+00;
1892     exp_abserr = 9.854016753016499034E-09;
1893
1894     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) AB ln(b-x) result") ;
1895     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) AB ln(b-x) abserr") ;
1896
1897     /* Test with ln(x - a) ln(b - x) */
1898
1899     gsl_integration_qaws_table_set (t, -0.5, -0.3, 1, 1);
1900
1901     status = gsl_integration_qaws (&fc, 0.0, 1.0, t, 0.0, 1.0e-7, w->limit,
1902                                    w, &result, &abserr) ;
1903
1904     exp_result = 3.159922862811048172E-01;
1905     exp_abserr = 2.336183482198144595E-08;
1906
1907     gsl_test_rel(result,exp_result,1e-14,"qaws(f458) AB ln(x-a)ln(b-x) result") ;
1908     gsl_test_rel(abserr,exp_abserr,1e-6,"qaws(f458) AB ln(x-a)ln(b-x) abserr") ;
1909
1910     gsl_integration_workspace_free (w) ;
1911     gsl_integration_qaws_table_free (t) ;
1912
1913   }
1914
1915
1916   /* Test oscillatory integration using a relative error bound */
1917
1918   {
1919     int status = 0, i; struct counter_params p;
1920     double result = 0, abserr=0;
1921
1922     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
1923     gsl_integration_qawo_table * wo
1924       = gsl_integration_qawo_table_alloc (10.0 * M_PI, 1.0,
1925                                               GSL_INTEG_SINE, 1000) ;
1926
1927     /* All results are for GSL_IEEE_MODE=double-precision */
1928
1929     double exp_result = -1.281368483991674190E-01;
1930     double exp_abserr =  6.875028324415666248E-12;
1931     int exp_neval  =      305;
1932     int exp_ier    =        0;
1933     int exp_last   =        9;
1934
1935     double a[9] = { 0.000000000000000000E+00,
1936                     5.000000000000000000E-01,
1937                     2.500000000000000000E-01,
1938                     1.250000000000000000E-01,
1939                     6.250000000000000000E-02,
1940                     3.125000000000000000E-02,
1941                     1.562500000000000000E-02,
1942                     7.812500000000000000E-03,
1943                     3.906250000000000000E-03 } ;
1944     double b[9] = { 3.906250000000000000E-03,
1945                     1.000000000000000000E+00,
1946                     5.000000000000000000E-01,
1947                     2.500000000000000000E-01,
1948                     1.250000000000000000E-01,
1949                     6.250000000000000000E-02,
1950                     3.125000000000000000E-02,
1951                     1.562500000000000000E-02,
1952                     7.812500000000000000E-03 } ;
1953     double r[9] = { -1.447193692377651136E-03,
1954                     2.190541162282139478E-02,
1955                     -2.587726479625663753E-02,
1956                     5.483209176363500886E-02,
1957                     -3.081695575172510582E-02,
1958                     -9.178321994387816929E-02,
1959                     -3.886716016498160953E-02,
1960                     -1.242306301902117854E-02,
1961                     -3.659495117871544145E-03} ;
1962     double e[9] = { 8.326506625798146465E-07,
1963                     1.302638552580516100E-13,
1964                     7.259224351945759794E-15,
1965                     1.249770395036711102E-14,
1966                     7.832180081562836579E-16,
1967                     1.018998440559284116E-15,
1968                     4.315121611695628020E-16,
1969                     1.379237060008662177E-16,
1970                     4.062855738364339357E-17 } ;
1971     int order[9] = { 1, 2, 4, 3, 6, 5, 7, 8, 9 } ;
1972
1973     double alpha = 1.0 ;
1974     gsl_function f = make_function(&f456, &alpha);
1975     gsl_function fc = make_counter(&f, &p) ;
1976
1977     status = gsl_integration_qawo (&fc, 0.0, 0.0, 1e-7, w->limit,
1978                                    w, wo, &result, &abserr) ;
1979
1980     gsl_test_rel(result,exp_result,1e-14,"qawo(f456) result") ;
1981     gsl_test_rel(abserr,exp_abserr,1e-3,"qawo(f456) abserr") ;
1982     gsl_test_int((int)(p.neval),exp_neval,"qawo(f456) neval") ;
1983     gsl_test_int((int)(w->size),exp_last,"qawo(f456) last") ;
1984     gsl_test_int(status,exp_ier,"qawo(f456) status") ;
1985
1986     for (i = 0; i < 9 ; i++)
1987         gsl_test_rel(w->alist[i],a[i],1e-15,"qawo(f456) alist") ;
1988
1989     for (i = 0; i < 9 ; i++)
1990         gsl_test_rel(w->blist[i],b[i],1e-15,"qawo(f456) blist") ;
1991
1992     for (i = 0; i < 9 ; i++)
1993         gsl_test_rel(w->rlist[i],r[i],1e-14,"qawo(f456) rlist") ;
1994
1995     for (i = 0; i < 9 ; i++)
1996         gsl_test_rel(w->elist[i],e[i],1e-3,"qawo(f456) elist") ;
1997
1998     for (i = 0; i < 9 ; i++)
1999         gsl_test_int((int)w->order[i],order[i]-1,"qawo(f456) order");
2000
2001
2002     /* In reverse, flip limit and sign of length */
2003
2004     gsl_integration_qawo_table_set_length (wo, -1.0);
2005
2006     p.neval = 0;
2007     status = gsl_integration_qawo (&fc, 1.0, 0.0, 1e-7, w->limit,
2008                                    w, wo, &result, &abserr) ;
2009
2010     gsl_test_rel(result,-exp_result,1e-14,"qawo(f456) rev result") ;
2011     gsl_test_rel(abserr,exp_abserr,1e-3,"qawo(f456) rev abserr") ;
2012     gsl_test_int((int)(p.neval),exp_neval,"qawo(f456) rev neval") ;
2013     gsl_test_int((int)(w->size),exp_last,"qawo(f456) rev last") ;
2014     gsl_test_int(status,exp_ier,"qawo(f456) rev status") ;
2015
2016
2017     gsl_integration_qawo_table_free (wo) ;
2018     gsl_integration_workspace_free (w) ;
2019
2020   }
2021
2022   /* Test fourier integration using an absolute error bound */
2023
2024   {
2025     int status = 0, i; struct counter_params p;
2026     double result = 0, abserr=0;
2027
2028     gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;
2029     gsl_integration_workspace * wc = gsl_integration_workspace_alloc (1000) ;
2030     gsl_integration_qawo_table * wo
2031       = gsl_integration_qawo_table_alloc (M_PI / 2.0, 1.0,
2032                                               GSL_INTEG_COSINE, 1000) ;
2033
2034     /* All results are for GSL_IEEE_MODE=double-precision */
2035
2036     double exp_result = 9.999999999279802765E-01;
2037     double exp_abserr = 1.556289974669056164E-08;
2038     int exp_neval  =      590;
2039     int exp_ier    =        0;
2040     int exp_last   =       12;
2041
2042     double r[12] = { 1.013283128125232802E+00,
2043                     -1.810857954748607349E-02,
2044                     7.466754034900931897E-03,
2045                     -4.360312526786496237E-03,
2046                     2.950184068216192904E-03,
2047                     -2.168238443073697373E-03,
2048                     1.680910783140869081E-03,
2049                     -1.352797860944863345E-03,
2050                     1.119354921991485901E-03,
2051                     -9.462367583691360827E-04,
2052                     8.136341270731781887E-04,
2053                     -7.093931338504278145E-04 } ;
2054     double e[12] = { 1.224798040766472695E-12,
2055                     1.396565155187268456E-13,
2056                     1.053844511655910310E-16,
2057                     6.505213034913026604E-19,
2058                     7.155734338404329264E-18,
2059                     1.105886215935214523E-17,
2060                     9.757819552369539906E-18,
2061                     5.854691731421723944E-18,
2062                     4.553649124439220312E-18,
2063                     7.643625316022806260E-18,
2064                     2.439454888092388058E-17,
2065                     2.130457268934021451E-17 } ;
2066
2067     double alpha = 1.0 ;
2068     gsl_function f = make_function(&f457, &alpha);
2069     gsl_function fc = make_counter(&f, &p) ;
2070
2071     status = gsl_integration_qawf (&fc, 0.0, 1e-7, w->limit,
2072                                    w, wc, wo, &result, &abserr) ;
2073
2074     gsl_test_rel(result,exp_result,1e-14,"qawf(f457) result") ;
2075     gsl_test_rel(abserr,exp_abserr,1e-3,"qawf(f457) abserr") ;
2076     gsl_test_int((int)(p.neval),exp_neval,"qawf(f457) neval") ;
2077     gsl_test_int((int)(w->size),exp_last,"qawf(f457) last") ;
2078     gsl_test_int(status,exp_ier,"qawf(f457) status") ;
2079
2080     for (i = 0; i < 9 ; i++)
2081         gsl_test_rel(w->rlist[i],r[i],1e-12,"qawf(f457) rlist") ;
2082
2083     /* We can only get within two orders of magnitude on the error
2084        here, which is very sensitive to the floating point precision */
2085
2086     for (i = 0; i < 9 ; i++)
2087         gsl_test_rel(w->elist[i],e[i],50.0,"qawf(f457) elist") ;
2088
2089
2090     gsl_integration_qawo_table_free (wo) ;
2091     gsl_integration_workspace_free (wc) ;
2092     gsl_integration_workspace_free (w) ;
2093
2094   }
2095
2096   exit (gsl_test_summary());
2097 }
2098
2099 void
2100 my_error_handler (const char *reason, const char *file, int line, int err)
2101 {
2102   if (0) printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err) ;
2103 }
2104
2105