htswanalysis/MACS/lib/gsl/gsl-1.11/diff/diff.c

   1 /* diff/diff.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 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 <gsl/gsl_math.h>
  23 #include <gsl/gsl_errno.h>
  24 #include <gsl/gsl_diff.h>
  25
  26 int
  27 gsl_diff_backward (const gsl_function * f,
  28                    double x, double *result, double *abserr)
  29 {
  30   /* Construct a divided difference table with a fairly large step
  31      size to get a very rough estimate of f''.  Use this to estimate
  32      the step size which will minimize the error in calculating f'. */
  33
  34   int i, k;
  35   double h = GSL_SQRT_DBL_EPSILON;
  36   double a[3], d[3], a2;
  37
  38   /* Algorithm based on description on pg. 204 of Conte and de Boor
  39      (CdB) - coefficients of Newton form of polynomial of degree 2. */
  40
  41   for (i = 0; i < 3; i++)
  42     {
  43       a[i] = x + (i - 2.0) * h;
  44       d[i] = GSL_FN_EVAL (f, a[i]);
  45     }
  46
  47   for (k = 1; k < 4; k++)
  48     {
  49       for (i = 0; i < 3 - k; i++)
  50         {
  51           d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
  52         }
  53     }
  54
  55   /* Adapt procedure described on pg. 282 of CdB to find best value of
  56      step size. */
  57
  58   a2 = fabs (d[0] + d[1] + d[2]);
  59
  60   if (a2 < 100.0 * GSL_SQRT_DBL_EPSILON)
  61     {
  62       a2 = 100.0 * GSL_SQRT_DBL_EPSILON;
  63     }
  64
  65   h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));
  66
  67   if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
  68     {
  69       h = 100.0 * GSL_SQRT_DBL_EPSILON;
  70     }
  71
  72   *result = (GSL_FN_EVAL (f, x) - GSL_FN_EVAL (f, x - h)) / h;
  73   *abserr = fabs (10.0 * a2 * h);
  74
  75   return GSL_SUCCESS;
  76 }
  77
  78 int
  79 gsl_diff_forward (const gsl_function * f,
  80                   double x, double *result, double *abserr)
  81 {
  82   /* Construct a divided difference table with a fairly large step
  83      size to get a very rough estimate of f''.  Use this to estimate
  84      the step size which will minimize the error in calculating f'. */
  85
  86   int i, k;
  87   double h = GSL_SQRT_DBL_EPSILON;
  88   double a[3], d[3], a2;
  89
  90   /* Algorithm based on description on pg. 204 of Conte and de Boor
  91      (CdB) - coefficients of Newton form of polynomial of degree 2. */
  92
  93   for (i = 0; i < 3; i++)
  94     {
  95       a[i] = x + i * h;
  96       d[i] = GSL_FN_EVAL (f, a[i]);
  97     }
  98
  99   for (k = 1; k < 4; k++)
 100     {
 101       for (i = 0; i < 3 - k; i++)
 102         {
 103           d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
 104         }
 105     }
 106
 107   /* Adapt procedure described on pg. 282 of CdB to find best value of
 108      step size. */
 109
 110   a2 = fabs (d[0] + d[1] + d[2]);
 111
 112   if (a2 < 100.0 * GSL_SQRT_DBL_EPSILON)
 113     {
 114       a2 = 100.0 * GSL_SQRT_DBL_EPSILON;
 115     }
 116
 117   h = sqrt (GSL_SQRT_DBL_EPSILON / (2.0 * a2));
 118
 119   if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
 120     {
 121       h = 100.0 * GSL_SQRT_DBL_EPSILON;
 122     }
 123
 124   *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x)) / h;
 125   *abserr = fabs (10.0 * a2 * h);
 126
 127   return GSL_SUCCESS;
 128 }
 129
 130 int
 131 gsl_diff_central (const gsl_function * f,
 132                   double x, double *result, double *abserr)
 133 {
 134   /* Construct a divided difference table with a fairly large step
 135      size to get a very rough estimate of f'''.  Use this to estimate
 136      the step size which will minimize the error in calculating f'. */
 137
 138   int i, k;
 139   double h = GSL_SQRT_DBL_EPSILON;
 140   double a[4], d[4], a3;
 141
 142   /* Algorithm based on description on pg. 204 of Conte and de Boor
 143      (CdB) - coefficients of Newton form of polynomial of degree 3. */
 144
 145   for (i = 0; i < 4; i++)
 146     {
 147       a[i] = x + (i - 2.0) * h;
 148       d[i] = GSL_FN_EVAL (f, a[i]);
 149     }
 150
 151   for (k = 1; k < 5; k++)
 152     {
 153       for (i = 0; i < 4 - k; i++)
 154         {
 155           d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i]);
 156         }
 157     }
 158
 159   /* Adapt procedure described on pg. 282 of CdB to find best
 160      value of step size. */
 161
 162   a3 = fabs (d[0] + d[1] + d[2] + d[3]);
 163
 164   if (a3 < 100.0 * GSL_SQRT_DBL_EPSILON)
 165     {
 166       a3 = 100.0 * GSL_SQRT_DBL_EPSILON;
 167     }
 168
 169   h = pow (GSL_SQRT_DBL_EPSILON / (2.0 * a3), 1.0 / 3.0);
 170
 171   if (h > 100.0 * GSL_SQRT_DBL_EPSILON)
 172     {
 173       h = 100.0 * GSL_SQRT_DBL_EPSILON;
 174     }
 175
 176   *result = (GSL_FN_EVAL (f, x + h) - GSL_FN_EVAL (f, x - h)) / (2.0 * h);
 177   *abserr = fabs (100.0 * a3 * h * h);
 178
 179   return GSL_SUCCESS;
 180 }