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

   1 /* deriv/test.c
   2  *
   3  * Copyright (C) 2000 David Morrison
   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_deriv.h>
  26 #include <gsl/gsl_errno.h>
  27 #include <gsl/gsl_test.h>
  28 #include <gsl/gsl_ieee_utils.h>
  29
  30 double
  31 f1 (double x, void *params)
  32 {
  33   return exp (x);
  34 }
  35
  36 double
  37 df1 (double x, void *params)
  38 {
  39   return exp (x);
  40 }
  41
  42 double
  43 f2 (double x, void *params)
  44 {
  45   if (x >= 0.0)
  46     {
  47       return x * sqrt (x);
  48     }
  49   else
  50     {
  51       return 0.0;
  52     }
  53 }
  54
  55 double
  56 df2 (double x, void *params)
  57 {
  58   if (x >= 0.0)
  59     {
  60       return 1.5 * sqrt (x);
  61     }
  62   else
  63     {
  64       return 0.0;
  65     }
  66 }
  67
  68 double
  69 f3 (double x, void *params)
  70 {
  71   if (x != 0.0)
  72     {
  73       return sin (1 / x);
  74     }
  75   else
  76     {
  77       return 0.0;
  78     }
  79 }
  80
  81 double
  82 df3 (double x, void *params)
  83 {
  84   if (x != 0.0)
  85     {
  86       return -cos (1 / x) / (x * x);
  87     }
  88   else
  89     {
  90       return 0.0;
  91     }
  92 }
  93
  94 double
  95 f4 (double x, void *params)
  96 {
  97   return exp (-x * x);
  98 }
  99
 100 double
 101 df4 (double x, void *params)
 102 {
 103   return -2.0 * x * exp (-x * x);
 104 }
 105
 106 double
 107 f5 (double x, void *params)
 108 {
 109   return x * x;
 110 }
 111
 112 double
 113 df5 (double x, void *params)
 114 {
 115   return 2.0 * x;
 116 }
 117
 118 double
 119 f6 (double x, void *params)
 120 {
 121   return 1.0 / x;
 122 }
 123
 124 double
 125 df6 (double x, void *params)
 126 {
 127   return -1.0 / (x * x);
 128 }
 129
 130 typedef int (deriv_fn) (const gsl_function * f, double x, double h, double * res, double *abserr);
 131
 132 void
 133 test (deriv_fn * deriv, gsl_function * f, gsl_function * df, double x,
 134       const char * desc)
 135 {
 136   double result, abserr;
 137   double expected = GSL_FN_EVAL (df, x);
 138   (*deriv) (f, x, 1e-4, &result, &abserr);
 139
 140   gsl_test_abs (result, expected, GSL_MIN(1e-4,fabs(expected)) + GSL_DBL_EPSILON, desc);
 141
 142   if (abserr < fabs(result-expected))
 143     {
 144       gsl_test_factor (abserr, fabs(result-expected), 2, "%s error estimate", desc);
 145     }
 146   else if (result == expected || expected == 0.0)
 147     {
 148       gsl_test_abs (abserr, 0.0, 1e-6, "%s abserr", desc);
 149     }
 150   else
 151     {
 152       double d = fabs(result - expected);
 153       gsl_test_abs (abserr, fabs(result-expected), 1e6*d, "%s abserr", desc);
 154     }
 155 }
 156
 157 int
 158 main ()
 159 {
 160   gsl_function F1, DF1, F2, DF2, F3, DF3, F4, DF4, F5, DF5, F6, DF6;
 161
 162   gsl_ieee_env_setup ();
 163
 164   F1.function = &f1;
 165   DF1.function = &df1;
 166
 167   F2.function = &f2;
 168   DF2.function = &df2;
 169
 170   F3.function = &f3;
 171   DF3.function = &df3;
 172
 173   F4.function = &f4;
 174   DF4.function = &df4;
 175
 176   F5.function = &f5;
 177   DF5.function = &df5;
 178
 179   F6.function = &f6;
 180   DF6.function = &df6;
 181
 182   test (&gsl_deriv_central, &F1, &DF1, 1.0, "exp(x), x=1, central deriv");
 183   test (&gsl_deriv_forward, &F1, &DF1, 1.0, "exp(x), x=1, forward deriv");
 184   test (&gsl_deriv_backward, &F1, &DF1, 1.0, "exp(x), x=1, backward deriv");
 185
 186   test (&gsl_deriv_central, &F2, &DF2, 0.1, "x^(3/2), x=0.1, central deriv");
 187   test (&gsl_deriv_forward, &F2, &DF2, 0.1, "x^(3/2), x=0.1, forward deriv");
 188   test (&gsl_deriv_backward, &F2, &DF2, 0.1, "x^(3/2), x=0.1, backward deriv");
 189
 190   test (&gsl_deriv_central, &F3, &DF3, 0.45, "sin(1/x), x=0.45, central deriv");
 191   test (&gsl_deriv_forward, &F3, &DF3, 0.45, "sin(1/x), x=0.45, forward deriv");
 192   test (&gsl_deriv_backward, &F3, &DF3, 0.45, "sin(1/x), x=0.45, backward deriv");
 193
 194   test (&gsl_deriv_central, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, central deriv");
 195   test (&gsl_deriv_forward, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, forward deriv");
 196   test (&gsl_deriv_backward, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, backward deriv");
 197
 198   test (&gsl_deriv_central, &F5, &DF5, 0.0, "x^2, x=0, central deriv");
 199   test (&gsl_deriv_forward, &F5, &DF5, 0.0, "x^2, x=0, forward deriv");
 200   test (&gsl_deriv_backward, &F5, &DF5, 0.0, "x^2, x=0, backward deriv");
 201
 202   test (&gsl_deriv_central, &F6, &DF6, 10.0, "1/x, x=10, central deriv");
 203   test (&gsl_deriv_forward, &F6, &DF6, 10.0, "1/x, x=10, forward deriv");
 204   test (&gsl_deriv_backward, &F6, &DF6, 10.0, "1/x, x=10, backward deriv");
 205
 206   exit (gsl_test_summary ());
 207 }
 208
 209