htswanalysis/MACS/lib/gsl/gsl-1.11/interpolation/poly.c

   1 /* interpolation/interp_poly.c
   2  *
   3  * Copyright (C) 2001 DAN, HO-JIN
   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 <math.h>
  22 #include <gsl/gsl_errno.h>
  23 #include <gsl/gsl_poly.h>
  24 #include <gsl/gsl_interp.h>
  25
  26 typedef struct
  27 {
  28   double *d;
  29   double *coeff;
  30   double *work;
  31 }
  32 polynomial_state_t;
  33
  34 static void *
  35 polynomial_alloc (size_t size)
  36 {
  37   polynomial_state_t *state =
  38     (polynomial_state_t *) malloc (sizeof (polynomial_state_t));
  39
  40   if (state == 0)
  41     {
  42       GSL_ERROR_NULL ("failed to allocate space for polynomial state",
  43                       GSL_ENOMEM);
  44     }
  45
  46   state->d = (double *) malloc (sizeof (double) * size);
  47
  48   if (state->d == 0)
  49     {
  50       free (state);
  51       GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
  52     }
  53
  54   state->coeff = (double *) malloc (sizeof (double) * size);
  55
  56   if (state->coeff == 0)
  57     {
  58       free (state->d);
  59       free (state);
  60       GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
  61     }
  62
  63   state->work = (double *) malloc (sizeof (double) * size);
  64
  65   if (state->work == 0)
  66     {
  67       free (state->coeff);
  68       free (state->d);
  69       free (state);
  70       GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
  71     }
  72
  73   return state;
  74 }
  75
  76 static int
  77 polynomial_init (void *vstate,
  78                  const double xa[], const double ya[], size_t size)
  79 {
  80   polynomial_state_t *state = (polynomial_state_t *) vstate;
  81
  82   int status = gsl_poly_dd_init (state->d, xa, ya, size);
  83
  84   return status;
  85 }
  86
  87 static int
  88 polynomial_eval (const void *vstate,
  89                  const double xa[], const double ya[], size_t size, double x,
  90                  gsl_interp_accel * acc, double *y)
  91 {
  92   const polynomial_state_t *state = (const polynomial_state_t *) vstate;
  93
  94   *y = gsl_poly_dd_eval (state->d, xa, size, x);
  95
  96   return GSL_SUCCESS;
  97 }
  98
  99
 100 static int
 101 polynomial_deriv (const void *vstate,
 102                   const double xa[], const double ya[], size_t size, double x,
 103                   gsl_interp_accel * acc, double *y)
 104 {
 105   const polynomial_state_t *state = (const polynomial_state_t *) vstate;
 106
 107   gsl_poly_dd_taylor (state->coeff, x, state->d, xa, size, state->work);
 108
 109   *y = state->coeff[1];
 110
 111   return GSL_SUCCESS;
 112 }
 113
 114 static int
 115 polynomial_deriv2 (const void *vstate,
 116                    const double xa[], const double ya[], size_t size,
 117                    double x, gsl_interp_accel * acc, double *y)
 118 {
 119   const polynomial_state_t *state = (const polynomial_state_t *) vstate;
 120
 121   gsl_poly_dd_taylor (state->coeff, x, state->d, xa, size, state->work);
 122
 123   *y = 2.0 * state->coeff[2];
 124
 125   return GSL_SUCCESS;
 126 }
 127
 128 static int
 129 polynomial_integ (const void *vstate, const double xa[], const double ya[],
 130                   size_t size, gsl_interp_accel * acc, double a, double b,
 131                   double *result)
 132 {
 133   const polynomial_state_t *state = (const polynomial_state_t *) vstate;
 134   size_t i;
 135   double sum;
 136
 137   gsl_poly_dd_taylor (state->coeff, 0.0, state->d, xa, size, state->work);
 138
 139   sum = state->coeff[0] * (b - a);
 140
 141   for (i = 1; i < size; i++)
 142     {
 143       sum += state->coeff[i] * (pow (b, i + 1) - pow (a, i + 1)) / (i + 1.0);
 144     }
 145
 146   *result = sum;
 147
 148   return GSL_SUCCESS;
 149 }
 150
 151 static void
 152 polynomial_free (void *vstate)
 153 {
 154   polynomial_state_t *state = (polynomial_state_t *) vstate;
 155   free (state->d);
 156   free (state->coeff);
 157   free (state->work);
 158   free (state);
 159 }
 160
 161 static const gsl_interp_type polynomial_type = {
 162   "polynomial",
 163   3,
 164   &polynomial_alloc,
 165   &polynomial_init,
 166   &polynomial_eval,
 167   &polynomial_deriv,
 168   &polynomial_deriv2,
 169   &polynomial_integ,
 170   &polynomial_free,
 171 };
 172
 173 const gsl_interp_type *gsl_interp_polynomial = &polynomial_type;