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

   1 /* roots/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Reid Priedhorsky, 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 <gsl/gsl_math.h>
  23 #include <gsl/gsl_test.h>
  24 #include <gsl/gsl_roots.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_ieee_utils.h>
  27
  28 #include "roots.h"
  29 #include "test.h"
  30
  31 /* stopping parameters */
  32 const double EPSREL = (10 * GSL_DBL_EPSILON);
  33 const double EPSABS = (10 * GSL_DBL_EPSILON);
  34 const unsigned int MAX_ITERATIONS = 150;
  35
  36 void my_error_handler (const char *reason, const char *file,
  37                        int line, int err);
  38
  39 #define WITHIN_TOL(a, b, epsrel, epsabs) \
  40  ((fabs((a) - (b)) < (epsrel) * GSL_MIN(fabs(a), fabs(b)) + (epsabs)))
  41
  42 int
  43 main (void)
  44 {
  45   gsl_function F_sin, F_cos, F_func1, F_func2, F_func3, F_func4,
  46     F_func5, F_func6;
  47
  48   gsl_function_fdf FDF_sin, FDF_cos, FDF_func1, FDF_func2, FDF_func3, FDF_func4,
  49     FDF_func5, FDF_func6;
  50
  51   const gsl_root_fsolver_type * fsolver[4] ;
  52   const gsl_root_fdfsolver_type * fdfsolver[4] ;
  53
  54   const gsl_root_fsolver_type ** T;
  55   const gsl_root_fdfsolver_type ** S;
  56
  57   gsl_ieee_env_setup();
  58
  59   fsolver[0] = gsl_root_fsolver_bisection;
  60   fsolver[1] = gsl_root_fsolver_brent;
  61   fsolver[2] = gsl_root_fsolver_falsepos;
  62   fsolver[3] = 0;
  63
  64   fdfsolver[0] = gsl_root_fdfsolver_newton;
  65   fdfsolver[1] = gsl_root_fdfsolver_secant;
  66   fdfsolver[2] = gsl_root_fdfsolver_steffenson;
  67   fdfsolver[3] = 0;
  68
  69   F_sin = create_function (sin_f) ;
  70   F_cos = create_function (cos_f) ;
  71   F_func1 = create_function (func1) ;
  72   F_func2 = create_function (func2) ;
  73   F_func3 = create_function (func3) ;
  74   F_func4 = create_function (func4) ;
  75   F_func5 = create_function (func5) ;
  76   F_func6 = create_function (func6) ;
  77
  78   FDF_sin = create_fdf (sin_f, sin_df, sin_fdf) ;
  79   FDF_cos = create_fdf (cos_f, cos_df, cos_fdf) ;
  80   FDF_func1 = create_fdf (func1, func1_df, func1_fdf) ;
  81   FDF_func2 = create_fdf (func2, func2_df, func2_fdf) ;
  82   FDF_func3 = create_fdf (func3, func3_df, func3_fdf) ;
  83   FDF_func4 = create_fdf (func4, func4_df, func4_fdf) ;
  84   FDF_func5 = create_fdf (func5, func5_df, func5_fdf) ;
  85   FDF_func6 = create_fdf (func6, func6_df, func6_fdf) ;
  86
  87   gsl_set_error_handler (&my_error_handler);
  88
  89   for (T = fsolver ; *T != 0 ; T++)
  90     {
  91       test_f (*T, "sin(x) [3, 4]", &F_sin, 3.0, 4.0, M_PI);
  92       test_f (*T, "sin(x) [-4, -3]", &F_sin, -4.0, -3.0, -M_PI);
  93       test_f (*T, "sin(x) [-1/3, 1]", &F_sin, -1.0 / 3.0, 1.0, 0.0);
  94       test_f (*T, "cos(x) [0, 3]", &F_cos, 0.0, 3.0, M_PI / 2.0);
  95       test_f (*T, "cos(x) [-3, 0]", &F_cos, -3.0, 0.0, -M_PI / 2.0);
  96       test_f (*T, "x^20 - 1 [0.1, 2]", &F_func1, 0.1, 2.0, 1.0);
  97       test_f (*T, "sqrt(|x|)*sgn(x)", &F_func2, -1.0 / 3.0, 1.0, 0.0);
  98       test_f (*T, "x^2 - 1e-8 [0, 1]", &F_func3, 0.0, 1.0, sqrt (1e-8));
  99       test_f (*T, "x exp(-x) [-1/3, 2]", &F_func4, -1.0 / 3.0, 2.0, 0.0);
 100       test_f (*T, "(x - 1)^7 [0.9995, 1.0002]", &F_func6, 0.9995, 1.0002, 1.0);
 101
 102       test_f_e (*T, "invalid range check [4, 0]", &F_sin, 4.0, 0.0, M_PI);
 103       test_f_e (*T, "invalid range check [1, 1]", &F_sin, 1.0, 1.0, M_PI);
 104       test_f_e (*T, "invalid range check [0.1, 0.2]", &F_sin, 0.1, 0.2, M_PI);
 105     }
 106
 107   for (S = fdfsolver ; *S != 0 ; S++)
 108     {
 109       test_fdf (*S,"sin(x) {3.4}", &FDF_sin, 3.4, M_PI);
 110       test_fdf (*S,"sin(x) {-3.3}", &FDF_sin, -3.3, -M_PI);
 111       test_fdf (*S,"sin(x) {0.5}", &FDF_sin, 0.5, 0.0);
 112       test_fdf (*S,"cos(x) {0.6}", &FDF_cos, 0.6, M_PI / 2.0);
 113       test_fdf (*S,"cos(x) {-2.5}", &FDF_cos, -2.5, -M_PI / 2.0);
 114       test_fdf (*S,"x^{20} - 1 {0.9}", &FDF_func1, 0.9, 1.0);
 115       test_fdf (*S,"x^{20} - 1 {1.1}", &FDF_func1, 1.1, 1.0);
 116       test_fdf (*S,"sqrt(|x|)*sgn(x) {1.001}", &FDF_func2, 0.001, 0.0);
 117       test_fdf (*S,"x^2 - 1e-8 {1}", &FDF_func3, 1.0, sqrt (1e-8));
 118       test_fdf (*S,"x exp(-x) {-2}", &FDF_func4, -2.0, 0.0);
 119       test_fdf_e (*S,"max iterations x -> +Inf, x exp(-x) {2}", &FDF_func4, 2.0, 0.0);
 120       test_fdf_e (*S,"max iterations x -> -Inf, 1/(1 + exp(-x)) {0}", &FDF_func5, 0.0, 0.0);
 121     }
 122
 123   test_fdf (gsl_root_fdfsolver_steffenson,
 124             "(x - 1)^7 {0.9}", &FDF_func6, 0.9, 1.0);
 125
 126   /* now summarize the results */
 127
 128   exit (gsl_test_summary ());
 129 }
 130
 131
 132 /* Using gsl_root_bisection, find the root of the function pointed to by f,
 133    using the interval [lower_bound, upper_bound]. Check if f succeeded and
 134    that it was accurate enough. */
 135
 136 void
 137 test_f (const gsl_root_fsolver_type * T, const char * description, gsl_function *f,
 138         double lower_bound, double upper_bound, double correct_root)
 139 {
 140   int status;
 141   size_t iterations = 0;
 142   double r, a, b;
 143   double x_lower, x_upper;
 144   gsl_root_fsolver * s;
 145
 146   x_lower = lower_bound;
 147   x_upper = upper_bound;
 148
 149   s = gsl_root_fsolver_alloc(T);
 150   gsl_root_fsolver_set(s, f, x_lower, x_upper) ;
 151
 152   do
 153     {
 154       iterations++ ;
 155
 156       gsl_root_fsolver_iterate (s);
 157
 158       r = gsl_root_fsolver_root(s);
 159
 160       a = gsl_root_fsolver_x_lower(s);
 161       b = gsl_root_fsolver_x_upper(s);
 162
 163       if (a > b)
 164         gsl_test (GSL_FAILURE, "interval is invalid (%g,%g)", a, b);
 165
 166       if (r < a || r > b)
 167         gsl_test (GSL_FAILURE, "r lies outside interval %g (%g,%g)", r, a, b);
 168
 169       status = gsl_root_test_interval (a,b, EPSABS, EPSREL);
 170     }
 171   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
 172
 173
 174   gsl_test (status, "%s, %s (%g obs vs %g expected) ",
 175             gsl_root_fsolver_name(s), description,
 176             gsl_root_fsolver_root(s), correct_root);
 177
 178   if (iterations == MAX_ITERATIONS)
 179     {
 180       gsl_test (GSL_FAILURE, "exceeded maximum number of iterations");
 181     }
 182
 183   /* check the validity of the returned result */
 184
 185   if (!WITHIN_TOL (r, correct_root, EPSREL, EPSABS))
 186     {
 187       gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)",
 188                 r, correct_root);
 189
 190     }
 191
 192   gsl_root_fsolver_free(s);
 193 }
 194
 195 void
 196 test_f_e (const gsl_root_fsolver_type * T,
 197           const char * description, gsl_function *f,
 198           double lower_bound, double upper_bound, double correct_root)
 199 {
 200   int status;
 201   size_t iterations = 0;
 202   double x_lower, x_upper;
 203   gsl_root_fsolver * s;
 204
 205   x_lower = lower_bound;
 206   x_upper = upper_bound;
 207
 208   s = gsl_root_fsolver_alloc(T);
 209   status = gsl_root_fsolver_set(s, f, x_lower, x_upper) ;
 210
 211   gsl_test (status != GSL_EINVAL, "%s (set), %s", T->name, description);
 212
 213   if (status == GSL_EINVAL)
 214     {
 215       gsl_root_fsolver_free(s);
 216       return ;
 217     }
 218
 219   do
 220     {
 221       iterations++ ;
 222       gsl_root_fsolver_iterate (s);
 223       x_lower = gsl_root_fsolver_x_lower(s);
 224       x_upper = gsl_root_fsolver_x_lower(s);
 225       status = gsl_root_test_interval (x_lower, x_upper,
 226                                       EPSABS, EPSREL);
 227     }
 228   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
 229
 230   gsl_test (!status, "%s, %s", gsl_root_fsolver_name(s), description,
 231             gsl_root_fsolver_root(s) - correct_root);
 232
 233   gsl_root_fsolver_free(s);
 234 }
 235
 236 void
 237 test_fdf (const gsl_root_fdfsolver_type * T, const char * description,
 238         gsl_function_fdf *fdf, double root, double correct_root)
 239 {
 240   int status;
 241   size_t iterations = 0;
 242   double prev = 0 ;
 243
 244   gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc(T);
 245   gsl_root_fdfsolver_set (s, fdf, root) ;
 246
 247   do
 248     {
 249       iterations++ ;
 250       prev = gsl_root_fdfsolver_root(s);
 251       gsl_root_fdfsolver_iterate (s);
 252       status = gsl_root_test_delta(gsl_root_fdfsolver_root(s), prev,
 253                                    EPSABS, EPSREL);
 254     }
 255   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
 256
 257   gsl_test (status, "%s, %s (%g obs vs %g expected) ",
 258             gsl_root_fdfsolver_name(s), description,
 259             gsl_root_fdfsolver_root(s), correct_root);
 260
 261   if (iterations == MAX_ITERATIONS)
 262     {
 263       gsl_test (GSL_FAILURE, "exceeded maximum number of iterations");
 264     }
 265
 266   /* check the validity of the returned result */
 267
 268   if (!WITHIN_TOL (gsl_root_fdfsolver_root(s), correct_root,
 269                    EPSREL, EPSABS))
 270     {
 271       gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)",
 272                 gsl_root_fdfsolver_root(s), correct_root);
 273
 274     }
 275   gsl_root_fdfsolver_free(s);
 276 }
 277
 278 void
 279 test_fdf_e (const gsl_root_fdfsolver_type * T,
 280             const char * description, gsl_function_fdf *fdf,
 281             double root, double correct_root)
 282 {
 283   int status;
 284   size_t iterations = 0;
 285   double prev = 0 ;
 286
 287   gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc(T);
 288   status = gsl_root_fdfsolver_set (s, fdf, root) ;
 289
 290   gsl_test (status, "%s (set), %s", T->name, description);
 291
 292   do
 293     {
 294       iterations++ ;
 295       prev = gsl_root_fdfsolver_root(s);
 296       gsl_root_fdfsolver_iterate (s);
 297       status = gsl_root_test_delta(gsl_root_fdfsolver_root(s), prev,
 298                                    EPSABS, EPSREL);
 299     }
 300   while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);
 301
 302   gsl_test (!status, "%s, %s", gsl_root_fdfsolver_name(s),
 303             description, gsl_root_fdfsolver_root(s) - correct_root);
 304   gsl_root_fdfsolver_free(s);
 305 }
 306
 307 void
 308 my_error_handler (const char *reason, const char *file, int line, int err)
 309 {
 310   if (0)
 311     printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
 312 }
 313
 314
 315