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

   1 /* integration/tests.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 <math.h>
  22 #include <gsl/gsl_math.h>
  23
  24 #include "tests.h"
  25
  26 /* These are the test functions from table 4.1 of the QUADPACK book */
  27
  28 /* f1(x) = x^alpha * log(1/x) */
  29 /* integ(f1,x,0,1) = 1/(alpha + 1)^2 */
  30
  31 double f1 (double x, void * params) {
  32   double alpha = *(double *) params ;
  33   return pow(x,alpha) * log(1/x) ;
  34 }
  35
  36 /* f2(x) = 4^-alpha / ((x-pi/4)^2 + 16^-alpha) */
  37 /* integ(f2,x,0,1) = arctan((4-pi)4^(alpha-1)) + arctan(pi 4^(alpha-1)) */
  38
  39 double f2 (double x, void * params) {
  40   double alpha = *(double *) params ;
  41   return pow(4.0,-alpha) / (pow((x-M_PI/4.0),2.0) + pow(16.0,-alpha)) ;
  42 }
  43
  44 /* f3(x) = cos(2^alpha * sin(x)) */
  45 /* integ(f3,x,0,pi) = pi J_0(2^alpha) */
  46
  47 double f3 (double x, void * params) {
  48   double alpha = *(double *) params ;
  49   return cos(pow(2.0,alpha) * sin(x)) ;
  50 }
  51
  52 /* Functions 4, 5 and 6 are duplicates of functions  1, 2 and 3 */
  53 /* ....                                                         */
  54
  55 /* f7(x) = |x - 1/3|^alpha */
  56 /* integ(f7,x,0,1) = ((2/3)^(alpha+1) + (1/3)^(alpha+1))/(alpha + 1) */
  57
  58 double f7 (double x, void * params) {
  59   double alpha = *(double *) params ;
  60   return pow(fabs(x - (1.0/3.0)),alpha) ;
  61 }
  62
  63 /* f8(x) = |x - pi/4|^alpha */
  64 /* integ(f8,x,0,1) =
  65    ((1 - pi/4)^(alpha+1) + (pi/4)^(alpha+1))/(alpha + 1) */
  66
  67 double f8 (double x, void * params) {
  68   double alpha = *(double *) params ;
  69   return pow(fabs(x - (M_PI/4.0)),alpha) ;
  70 }
  71
  72 /* f9(x) = sqrt(1 - x^2) / (x + 1 + 2^-alpha) */
  73 /* integ(f9,x,-1,1) = pi/sqrt((1+2^-alpha)^2-1) */
  74
  75 double f9 (double x, void * params) {
  76   double alpha = *(double *) params ;
  77   return 1 / ((x + 1 + pow(2.0,-alpha)) * sqrt(1-x*x)) ;
  78 }
  79
  80 /* f10(x) = sin(x)^(alpha - 1) */
  81 /* integ(f10,x,0,pi/2) = 2^(alpha-2) ((Gamma(alpha/2))^2)/Gamma(alpha) */
  82
  83 double f10 (double x, void * params) {
  84   double alpha = *(double *) params ;
  85   return pow(sin(x), alpha-1) ;
  86 }
  87
  88 /* f11(x) = log(1/x)^(alpha - 1) */
  89 /* integ(f11,x,0,1) = Gamma(alpha) */
  90
  91 double f11 (double x, void * params) {
  92   double alpha = *(double *) params ;
  93   return pow(log(1/x), alpha-1) ;
  94 }
  95
  96 /* f12(x) = exp(20*(x-1)) * sin(2^alpha * x) */
  97 /* integ(f12,x,0,1) =
  98    (20 sin(2^alpha) - 2^alpha cos(2^alpha) + 2^alpha exp(-20))
  99    /(400 + 4^alpha) */
 100
 101 double f12 (double x, void * params) {
 102   double alpha = *(double *) params ;
 103   return exp(20*(x-1)) * sin(pow(2.0,alpha) * x) ;
 104 }
 105
 106 /* f13(x) = cos(2^alpha * x)/sqrt(x(1 - x)) */
 107 /* integ(f13,x,0,1) = pi cos(2^(alpha-1)) J_0(2^(alpha-1))  */
 108
 109 double f13 (double x, void * params) {
 110   double alpha = *(double *) params ;
 111   return cos(pow(2.0,alpha)*x)/sqrt(x*(1-x)) ;
 112 }
 113
 114 double f14 (double x, void * params) {
 115   double alpha = *(double *) params ;
 116   return exp(-pow(2.0,-alpha)*x)*cos(x)/sqrt(x) ;
 117 }
 118
 119 double f15 (double x, void * params) {
 120   double alpha = *(double *) params ;
 121   return x*x * exp(-pow(2.0,-alpha)*x) ;
 122 }
 123
 124 double f16 (double x, void * params) {
 125   double alpha = *(double *) params ;
 126   if (x==0 && alpha == 1) return 1 ;  /* make the function continuous in x */
 127   if (x==0 && alpha > 1) return 0 ;   /* avoid problems with pow(0,1) */
 128   return pow(x,alpha-1)/pow((1+10*x),2.0) ;
 129 }
 130
 131 double f17 (double x, void * params) {
 132   double alpha = *(double *) params ;
 133   return pow(2.0,-alpha)/(((x-1)*(x-1)+pow(4.0,-alpha))*(x-2)) ;
 134 }
 135
 136 /* f454(x) = x^3 log|(x^2-1)(x^2-2)| */
 137 /* integ(f454,x,0,inf) = 61 log(2) + (77/4) log(7) - 27 */
 138
 139 double f454 (double x, void * params) {
 140   double x2 = x * x;
 141   double x3 = x * x2;
 142   params = 0 ;
 143   return x3 * log(fabs((x2 - 1.0) * (x2 - 2.0))) ;
 144 }
 145
 146 /* f455(x) = log(x)/(1+100*x^2) */
 147 /* integ(f455,x,0,inf) = -log(10)/20 */
 148
 149 double f455 (double x, void * params) {
 150   params = 0 ;
 151   return log(x) / (1.0 + 100.0 * x * x) ;
 152 }
 153
 154 /* f456(x) = log(x) */
 155 /* integ(f456*sin(10 pi x),x,0,1) = -(gamma + log(10pi) - Ci(10pi))/(10pi) */
 156
 157 double f456 (double x, void * params) {
 158   params = 0 ;
 159   if (x == 0.0)
 160     {
 161       return 0;
 162     }
 163   return log(x) ;
 164 }
 165
 166 /* f457(x) = 1/sqrt(x) */
 167 /* integ(f457*cos(pi x / 2),x,0,+inf) = 1 */
 168
 169 double f457 (double x, void * params) {
 170   params = 0 ;
 171   if (x == 0.0)
 172     {
 173       return 0;
 174     }
 175   return 1/sqrt(x) ;
 176 }
 177
 178 /* f458(x) = 1/(1 + log(x)^2)^2 */
 179 /* integ(log(x) f458(x),x,0,1) = (Ci(1) sin(1) + (pi/2 - Si(1)) cos(1))/pi
 180                                = -0.1892752 */
 181
 182 double f458 (double x, void * params) {
 183   params = 0 ;
 184
 185   if (x == 0.0)
 186     {
 187       return 0;
 188     }
 189   else
 190     {
 191       double u = log(x);
 192       double v = 1 + u * u;
 193
 194       return 1.0 / (v * v) ;
 195     }
 196 }
 197
 198 /* f459(x) = 1/(5 x^3 + 6) */
 199 /* integ(f459/(x-0),x,-1,5) = log(125/631)/18 */
 200
 201 double f459 (double x, void * params) {
 202   params = 0 ;
 203   return 1.0 / (5.0 * x * x * x + 6.0) ;
 204 }
 205
 206 /* myfn1(x) = exp(-x - x^2) */
 207 /* integ(myfn1,x,-inf,inf) = sqrt(pi) exp(-1/4) */
 208
 209 double myfn1 (double x, void * params) {
 210   params = 0;
 211   return exp(-x - x*x) ;
 212 }
 213
 214 /* myfn2(x) = exp(alpha*x) */
 215 /* integ(myfn2,x,-inf,b) = exp(alpha*b)/alpha */
 216
 217 double myfn2 (double x, void * params) {
 218   double alpha = *(double *) params ;
 219   return exp(alpha*x) ;
 220 }