htswanalysis/MACS/lib/gsl/gsl-1.11/poly/gsl_poly.h

   1 /* poly/gsl_poly.h
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007 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 #ifndef __GSL_POLY_H__
  21 #define __GSL_POLY_H__
  22
  23 #include <stdlib.h>
  24 #include <gsl/gsl_complex.h>
  25
  26 #undef __BEGIN_DECLS
  27 #undef __END_DECLS
  28 #ifdef __cplusplus
  29 # define __BEGIN_DECLS extern "C" {
  30 # define __END_DECLS }
  31 #else
  32 # define __BEGIN_DECLS /* empty */
  33 # define __END_DECLS /* empty */
  34 #endif
  35
  36 __BEGIN_DECLS
  37
  38
  39 /* Evaluate polynomial
  40  *
  41  * c[0] + c[1] x + c[2] x^2 + ... + c[len-1] x^(len-1)
  42  *
  43  * exceptions: none
  44  */
  45
  46 /* real polynomial, real x */
  47 double gsl_poly_eval(const double c[], const int len, const double x);
  48
  49 /* real polynomial, complex x */
  50 gsl_complex gsl_poly_complex_eval (const double c [], const int len, const gsl_complex z);
  51
  52 /* complex polynomial, complex x */
  53 gsl_complex gsl_complex_poly_complex_eval (const gsl_complex c [], const int len, const gsl_complex z);
  54
  55
  56 #ifdef HAVE_INLINE
  57 extern inline
  58 double
  59 gsl_poly_eval(const double c[], const int len, const double x)
  60 {
  61   int i;
  62   double ans = c[len-1];
  63   for(i=len-1; i>0; i--) ans = c[i-1] + x * ans;
  64   return ans;
  65 }
  66
  67 extern inline
  68 gsl_complex
  69 gsl_poly_complex_eval(const double c[], const int len, const gsl_complex z)
  70 {
  71   int i;
  72   gsl_complex ans;
  73   GSL_SET_COMPLEX (&ans, c[len-1], 0.0);
  74   for(i=len-1; i>0; i--) {
  75     /* The following three lines are equivalent to
  76        ans = gsl_complex_add_real (gsl_complex_mul (z, ans), c[i-1]);
  77        but faster */
  78     double tmp = c[i-1] + GSL_REAL (z) * GSL_REAL (ans) - GSL_IMAG (z) * GSL_IMAG (ans);
  79     GSL_SET_IMAG (&ans, GSL_IMAG (z) * GSL_REAL (ans) + GSL_REAL (z) * GSL_IMAG (ans));
  80     GSL_SET_REAL (&ans, tmp);
  81   }
  82   return ans;
  83 }
  84
  85 extern inline
  86 gsl_complex
  87 gsl_complex_poly_complex_eval(const gsl_complex c[], const int len, const gsl_complex z)
  88 {
  89   int i;
  90   gsl_complex ans = c[len-1];
  91   for(i=len-1; i>0; i--) {
  92     /* The following three lines are equivalent to
  93        ans = gsl_complex_add (c[i-1], gsl_complex_mul (x, ans));
  94        but faster */
  95     double tmp = GSL_REAL (c[i-1]) + GSL_REAL (z) * GSL_REAL (ans) - GSL_IMAG (z) * GSL_IMAG (ans);
  96     GSL_SET_IMAG (&ans, GSL_IMAG (c[i-1]) + GSL_IMAG (z) * GSL_REAL (ans) + GSL_REAL (z) * GSL_IMAG (ans));
  97     GSL_SET_REAL (&ans, tmp);
  98   }
  99   return ans;
 100 }
 101 #endif /* HAVE_INLINE */
 102
 103 /* Work with divided-difference polynomials, Abramowitz & Stegun 25.2.26 */
 104
 105 int
 106 gsl_poly_dd_init (double dd[], const double x[], const double y[],
 107                   size_t size);
 108
 109 double
 110 gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x);
 111
 112 #ifdef HAVE_INLINE
 113 extern inline
 114 double
 115 gsl_poly_dd_eval(const double dd[], const double xa[], const size_t size, const double x)
 116 {
 117   size_t i;
 118   double y = dd[size - 1];
 119   for (i = size - 1; i--;) y = dd[i] + (x - xa[i]) * y;
 120   return y;
 121 }
 122 #endif /* HAVE_INLINE */
 123
 124
 125 int
 126 gsl_poly_dd_taylor (double c[], double xp,
 127                     const double dd[], const double x[], size_t size,
 128                     double w[]);
 129
 130 /* Solve for real or complex roots of the standard quadratic equation,
 131  * returning the number of real roots.
 132  *
 133  * Roots are returned ordered.
 134  */
 135 int gsl_poly_solve_quadratic (double a, double b, double c,
 136                               double * x0, double * x1);
 137
 138 int
 139 gsl_poly_complex_solve_quadratic (double a, double b, double c,
 140                                   gsl_complex * z0, gsl_complex * z1);
 141
 142
 143 /* Solve for real roots of the cubic equation
 144  * x^3 + a x^2 + b x + c = 0, returning the
 145  * number of real roots.
 146  *
 147  * Roots are returned ordered.
 148  */
 149 int gsl_poly_solve_cubic (double a, double b, double c,
 150                           double * x0, double * x1, double * x2);
 151
 152 int
 153 gsl_poly_complex_solve_cubic (double a, double b, double c,
 154                               gsl_complex * z0, gsl_complex * z1,
 155                               gsl_complex * z2);
 156
 157
 158 /* Solve for the complex roots of a general real polynomial */
 159
 160 typedef struct
 161 {
 162   size_t nc ;
 163   double * matrix ;
 164 }
 165 gsl_poly_complex_workspace ;
 166
 167 gsl_poly_complex_workspace * gsl_poly_complex_workspace_alloc (size_t n);
 168 void gsl_poly_complex_workspace_free (gsl_poly_complex_workspace * w);
 169
 170 int
 171 gsl_poly_complex_solve (const double * a, size_t n,
 172                         gsl_poly_complex_workspace * w,
 173                         gsl_complex_packed_ptr z);
 174
 175 __END_DECLS
 176
 177 #endif /* __GSL_POLY_H__ */