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

   1 /* cheb/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   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_errno.h>
  24 #include <gsl/gsl_test.h>
  25 #include <gsl/gsl_ieee_utils.h>
  26 #include <gsl/gsl_chebyshev.h>
  27
  28 double f_T0 (double x, void * p) {
  29   p = 0;
  30   return 1.0;
  31 }
  32
  33 double f_T1 (double x, void * p) {
  34   p = 0;
  35   return x;
  36 }
  37
  38 double f_T2 (double x, void * p) {
  39   p = 0;
  40   return 2*x*x - 1;
  41 }
  42
  43 double f_sin (double x, void * p) {
  44   p = 0;
  45   return sin(x);
  46 }
  47
  48 int
  49 main(void)
  50 {
  51   double tol = 100.0 * GSL_DBL_EPSILON;
  52   double ftol = 20.0;
  53   double x;
  54   size_t i;
  55
  56   gsl_cheb_series * cs = gsl_cheb_alloc(40);
  57   gsl_cheb_series * csd = gsl_cheb_alloc(40);
  58   gsl_cheb_series * csi = gsl_cheb_alloc(40);
  59
  60   gsl_function F_sin, F_T0, F_T1, F_T2;
  61
  62   F_sin.function = f_sin;
  63   F_sin.params = 0;
  64
  65   F_T0.function = f_T0;
  66   F_T0.params = 0;
  67
  68   F_T1.function = f_T1;
  69   F_T1.params = 0;
  70
  71   F_T2.function = f_T2;
  72   F_T2.params = 0;
  73
  74   gsl_ieee_env_setup();
  75
  76   gsl_cheb_init(cs, &F_T0, -1.0, 1.0);
  77
  78   for (i = 0; i<cs->order; i++)
  79     {
  80       double c_exp = (i == 0) ? 2.0 : 0.0;
  81       gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_0(x)", i);
  82     }
  83
  84   gsl_cheb_init(cs, &F_T1, -1.0, 1.0);
  85
  86   for (i = 0; i<cs->order; i++)
  87     {
  88       double c_exp = (i == 1) ? 1.0 : 0.0;
  89       gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_1(x)", i);
  90     }
  91
  92   gsl_cheb_init(cs, &F_T2, -1.0, 1.0);
  93
  94   for (i = 0; i<cs->order; i++)
  95     {
  96       double c_exp = (i == 2) ? 1.0 : 0.0;
  97       gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_2(x)", i);
  98     }
  99
 100   gsl_cheb_init(cs, &F_sin, -M_PI, M_PI);
 101
 102   gsl_test_abs (cs->c[0], 0.0, tol, "c[0] for F_sin(x)");
 103   gsl_test_abs (cs->c[1], 5.69230686359506e-01, tol, "c[1] for F_sin(x)");
 104   gsl_test_abs (cs->c[2], 0.0, tol, "c[2] for F_sin(x)");
 105   gsl_test_abs (cs->c[3], -6.66916672405979e-01, tol, "c[3] for F_sin(x)");
 106   gsl_test_abs (cs->c[4], 0.0, tol, "c[4] for F_sin(x)");
 107   gsl_test_abs (cs->c[5], 1.04282368734237e-01, tol, "c[5] for F_sin(x)");
 108
 109   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 110     double r = gsl_cheb_eval(cs, x);
 111     gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval, sin(%.3g)", x);
 112   }
 113
 114   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 115     double r, e;
 116     gsl_cheb_eval_err(cs, x, &r, &e);
 117     gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval_err, sin(%.3g)", x);
 118     gsl_test_factor(fabs(r-sin(x)) + GSL_DBL_EPSILON, e, ftol,
 119                     "gsl_cheb_eval_err, error sin(%.3g)", x);
 120   }
 121
 122   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 123     double r = gsl_cheb_eval_n(cs, 25, x);
 124     gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval_n, sin(%.3g)", x);
 125   }
 126
 127   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 128     double r, e;
 129     gsl_cheb_eval_n_err(cs, 25, x, &r, &e);
 130     gsl_test_abs(r, sin(x), 100.0 * tol, "gsl_cheb_eval_n_err, deriv sin(%.3g)", x);
 131     gsl_test_factor(fabs(r-sin(x)) + GSL_DBL_EPSILON, e, ftol,
 132                     "gsl_cheb_eval_n_err, error sin(%.3g)", x);
 133   }
 134
 135   /* Test derivative */
 136
 137   gsl_cheb_calc_deriv(csd, cs);
 138
 139   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 140     double r = gsl_cheb_eval(csd, x);
 141     gsl_test_abs(r, cos(x), 1600 * tol, "gsl_cheb_eval, deriv sin(%.3g)", x);
 142   }
 143
 144 #ifdef TEST_DERIVATIVE_ERR
 145   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 146     double r, e;
 147     gsl_cheb_eval_err(csd, x, &r, &e);
 148     gsl_test_abs(r, cos(x), tol, "gsl_cheb_eval_err, deriv sin(%.3g)", x);
 149     gsl_test_factor(fabs(r-cos(x)) + GSL_DBL_EPSILON, e, ftol,
 150                     "gsl_cheb_eval_err, deriv error sin(%.3g)", x);
 151   }
 152 #endif
 153
 154   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 155     double r = gsl_cheb_eval_n(csd, 25, x);
 156     gsl_test_abs(r, cos(x), 1600 * tol, "gsl_cheb_eval_n, deriv sin(%.3g)", x);
 157   }
 158
 159 #ifdef TEST_DERIVATIVE_ERR
 160   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 161     double r, e;
 162     gsl_cheb_eval_n_err(csd, 25, x, &r, &e);
 163     gsl_test_abs(r, cos(x), 100.0 * tol, "gsl_cheb_eval_n_err, deriv sin(%.3g)", x);
 164     gsl_test_factor(fabs(r-cos(x)) + GSL_DBL_EPSILON, e, ftol,
 165                     "gsl_cheb_eval_n_err, deriv error sin(%.3g)", x);
 166   }
 167 #endif
 168
 169   /* Test integral */
 170
 171   gsl_cheb_calc_integ(csi, cs);
 172
 173   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 174     double r = gsl_cheb_eval(csi, x);
 175     gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval, integ sin(%.3g)", x);
 176   }
 177
 178 #ifdef TEST_INTEGRAL_ERR
 179   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 180     double r, e;
 181     gsl_cheb_eval_err(csi, x, &r, &e);
 182     gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval_err, integ sin(%.3g)", x);
 183     gsl_test_factor(fabs(r-(-1-cos(x))) + GSL_DBL_EPSILON, e, ftol,
 184                     "gsl_cheb_eval_err, integ error sin(%.3g)", x);
 185   }
 186 #endif
 187
 188   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 189     double r = gsl_cheb_eval_n(csi, 25, x);
 190     gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval_n, integ sin(%.3g)", x);
 191   }
 192
 193 #ifdef TEST_INTEGRAL_ERR
 194   for(x=-M_PI; x<M_PI; x += M_PI/100.0) {
 195     double r, e;
 196     gsl_cheb_eval_n_err(csi, 25, x, &r, &e);
 197     gsl_test_abs(r, -(1+cos(x)), 100.0 * tol, "gsl_cheb_eval_n_err, integ sin(%.3g)", x);
 198     gsl_test_factor(fabs(r-(-1-cos(x))) + GSL_DBL_EPSILON, e, ftol,
 199                     "gsl_cheb_eval_n_err, integ error sin(%.3g)", x);
 200   }
 201 #endif
 202
 203   gsl_cheb_free(csi);
 204   gsl_cheb_free(csd);
 205   gsl_cheb_free(cs);
 206
 207   exit (gsl_test_summary());
 208 }