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

   1 /* cheb/eval.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_chebyshev.h>
  25
  26 /* For efficiency there are separate implementations of each of these
  27    functions */
  28
  29 double
  30 gsl_cheb_eval (const gsl_cheb_series * cs, const double x)
  31 {
  32   size_t i;
  33   double d1 = 0.0;
  34   double d2 = 0.0;
  35
  36   double y = (2.0 * x - cs->a - cs->b) / (cs->b - cs->a);
  37   double y2 = 2.0 * y;
  38
  39   for (i = cs->order; i >= 1; i--)
  40     {
  41       double temp = d1;
  42       d1 = y2 * d1 - d2 + cs->c[i];
  43       d2 = temp;
  44     }
  45
  46   return y * d1 - d2 + 0.5 * cs->c[0];
  47 }
  48
  49 double
  50 gsl_cheb_eval_n (const gsl_cheb_series * cs, const size_t n, const double x)
  51 {
  52   size_t i;
  53   double d1 = 0.0;
  54   double d2 = 0.0;
  55
  56   size_t eval_order = GSL_MIN (n, cs->order);
  57
  58   double y = (2.0 * x - cs->a - cs->b) / (cs->b - cs->a);
  59   double y2 = 2.0 * y;
  60
  61   for (i = eval_order; i >= 1; i--)
  62     {
  63       double temp = d1;
  64       d1 = y2 * d1 - d2 + cs->c[i];
  65       d2 = temp;
  66     }
  67
  68   return y * d1 - d2 + 0.5 * cs->c[0];
  69 }
  70
  71
  72 int
  73 gsl_cheb_eval_err (const gsl_cheb_series * cs, const double x,
  74                    double *result, double *abserr)
  75 {
  76   size_t i;
  77   double d1 = 0.0;
  78   double d2 = 0.0;
  79
  80   double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
  81   double y2 = 2.0 * y;
  82
  83   double absc = 0.0;
  84
  85   for (i = cs->order; i >= 1; i--)
  86     {
  87       double temp = d1;
  88       d1 = y2 * d1 - d2 + cs->c[i];
  89       d2 = temp;
  90     }
  91
  92   *result = y * d1 - d2 + 0.5 * cs->c[0];
  93
  94   /* Estimate cumulative numerical error */
  95
  96   for (i = 0; i <= cs->order; i++)
  97     {
  98       absc += fabs(cs->c[i]);
  99     }
 100
 101   /* Combine truncation error and numerical error */
 102
 103   *abserr = fabs (cs->c[cs->order]) + absc * GSL_DBL_EPSILON;
 104
 105   return GSL_SUCCESS;
 106 }
 107
 108 int
 109 gsl_cheb_eval_n_err (const gsl_cheb_series * cs,
 110                      const size_t n, const double x,
 111                      double *result, double *abserr)
 112 {
 113   size_t i;
 114   double d1 = 0.0;
 115   double d2 = 0.0;
 116
 117   double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
 118   double y2 = 2.0 * y;
 119
 120   double absc = 0.0;
 121
 122   size_t eval_order = GSL_MIN (n, cs->order);
 123
 124   for (i = eval_order; i >= 1; i--)
 125     {
 126       double temp = d1;
 127       d1 = y2 * d1 - d2 + cs->c[i];
 128       d2 = temp;
 129     }
 130
 131   *result = y * d1 - d2 + 0.5 * cs->c[0];
 132
 133   /* Estimate cumulative numerical error */
 134
 135   for (i = 0; i <= eval_order; i++)
 136     {
 137       absc += fabs(cs->c[i]);
 138     }
 139
 140   /* Combine truncation error and numerical error */
 141
 142   *abserr = fabs (cs->c[eval_order]) + absc * GSL_DBL_EPSILON;
 143
 144   return GSL_SUCCESS;
 145 }
 146
 147 int
 148 gsl_cheb_eval_mode_e (const gsl_cheb_series * cs,
 149                       const double x, gsl_mode_t mode,
 150                       double *result, double *abserr)
 151 {
 152   size_t i;
 153   double d1 = 0.0;
 154   double d2 = 0.0;
 155
 156   double y = (2. * x - cs->a - cs->b) / (cs->b - cs->a);
 157   double y2 = 2.0 * y;
 158
 159   double absc = 0.0;
 160
 161   size_t eval_order;
 162
 163   if (GSL_MODE_PREC (mode) == GSL_PREC_DOUBLE)
 164     eval_order = cs->order;
 165   else
 166     eval_order = cs->order_sp;
 167
 168   for (i = eval_order; i >= 1; i--)
 169     {
 170       double temp = d1;
 171       d1 = y2 * d1 - d2 + cs->c[i];
 172       d2 = temp;
 173     }
 174
 175   *result = y * d1 - d2 + 0.5 * cs->c[0];
 176
 177   /* Estimate cumulative numerical error */
 178
 179   for (i = 0; i <= eval_order; i++)
 180     {
 181       absc += fabs(cs->c[i]);
 182     }
 183
 184   /* Combine truncation error and numerical error */
 185
 186   *abserr = fabs (cs->c[eval_order]) + absc * GSL_DBL_EPSILON;
 187
 188   return GSL_SUCCESS;
 189 }
 190
 191 double
 192 gsl_cheb_eval_mode (const gsl_cheb_series * cs,
 193                     const double x, gsl_mode_t mode)
 194 {
 195   double result, abserr;
 196   int status = gsl_cheb_eval_mode_e (cs, x, mode, &result, &abserr);
 197
 198   if (status != GSL_SUCCESS)
 199     {
 200       GSL_ERROR_VAL("gsl_cheb_eval_mode", status, result);
 201     };
 202
 203   return result;
 204 }
 205
 206