htswanalysis/MACS/lib/gsl/gsl-1.11/roots/test_funcs.c

   1 /* roots/test_funcs.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Reid Priedhorsky, 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 #include <config.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_roots.h>
  25
  26 #include "test.h"
  27
  28 gsl_function create_function (double (*f)(double, void *))
  29 {
  30   gsl_function F ;
  31   F.function = f;
  32   F.params = 0;
  33   return F ;
  34 }
  35
  36 gsl_function_fdf create_fdf (double (*f)(double, void *),
  37                              double (*df)(double, void *),
  38                              void (*fdf)(double, void *, double *, double *))
  39 {
  40   gsl_function_fdf FDF ;
  41   FDF.f = f ;
  42   FDF.df = df ;
  43   FDF.fdf = fdf ;
  44   FDF.params = 0 ;
  45   return FDF ;
  46 }
  47
  48 /* f(x) = x^{20} - 1 */
  49 /* f'(x) = 20x^{19} */
  50 /* zero at x = 1 or -1 */
  51
  52 double
  53 func1 (double x, void *p)
  54 {
  55   return pow (x, 20.0) - 1;
  56 }
  57
  58 double
  59 func1_df (double x, void * p)
  60 {
  61   return 20.0 * pow (x, 19.0);
  62 }
  63
  64 void
  65 func1_fdf (double x, void * p, double *y, double *yprime)
  66 {
  67   *y = func1 (x, p);
  68   *yprime = 20.0 * pow (x, 19.0);
  69 }
  70
  71 /* f(x) = sqrt(abs(x))*sgn(x) */
  72 /* f'(x) = 1 / sqrt(abs(x) */
  73 /* zero at x = 0 */
  74 double
  75 func2 (double x, void * p)
  76 {
  77   double delta;
  78
  79   if (x > 0)
  80     delta = 1.0;
  81   else if (x < 0)
  82     delta = -1.0;
  83   else
  84     delta = 0.0;
  85
  86   return sqrt (fabs (x)) * delta;
  87 }
  88
  89 double
  90 func2_df (double x, void * p)
  91 {
  92   return 1 / sqrt (fabs (x));
  93 }
  94
  95 void
  96 func2_fdf (double x, void * p, double *y, double *yprime)
  97 {
  98   *y = func2 (x, p);
  99   *yprime = 1 / sqrt (fabs (x));
 100 }
 101
 102
 103 /* f(x) = x^2 - 1e-8 */
 104 /* f'(x) = 2x */
 105 /* zero at x = sqrt(1e-8) or -sqrt(1e-8) */
 106 double
 107 func3 (double x, void * p)
 108 {
 109   return pow (x, 2.0) - 1e-8;
 110 }
 111
 112 double
 113 func3_df (double x, void * p)
 114 {
 115   return 2 * x;
 116 }
 117
 118 void
 119 func3_fdf (double x, void * p, double *y, double *yprime)
 120 {
 121   *y = func3 (x, p);
 122   *yprime = 2 * x;
 123 }
 124
 125 /* f(x) = x exp(-x) */
 126 /* f'(x) = exp(-x) - x exp(-x) */
 127 /* zero at x = 0 */
 128 double
 129 func4 (double x, void * p)
 130 {
 131   return x * exp (-x);
 132 }
 133
 134 double
 135 func4_df (double x, void * p)
 136 {
 137   return exp (-x) - x * exp (-x);
 138 }
 139
 140 void
 141 func4_fdf (double x, void * p, double *y, double *yprime)
 142 {
 143   *y = func4 (x, p);
 144   *yprime = exp (-x) - x * exp (-x);
 145 }
 146
 147 /* f(x) = 1/(1+exp(x)) */
 148 /* f'(x) = -exp(x) / (1 + exp(x))^2 */
 149 /* no roots! */
 150 double
 151 func5 (double x, void * p)
 152 {
 153   return 1 / (1 + exp (x));
 154 }
 155
 156 double
 157 func5_df (double x, void * p)
 158 {
 159   return -exp (x) / pow (1 + exp (x), 2.0);
 160 }
 161
 162 void
 163 func5_fdf (double x, void * p, double *y, double *yprime)
 164 {
 165   *y = func5 (x, p);
 166   *yprime = -exp (x) / pow (1 + exp (x), 2.0);
 167 }
 168
 169 /* f(x) = (x - 1)^7 */
 170 /* f'(x) = 7 * (x - 1)^6 */
 171 /* zero at x = 1 */
 172 double
 173 func6 (double x, void * p)
 174 {
 175   return pow (x - 1, 7.0);
 176 }
 177
 178 double
 179 func6_df (double x, void * p)
 180 {
 181   return 7.0 * pow (x - 1, 6.0);
 182 }
 183
 184 void
 185 func6_fdf (double x, void * p, double *y, double *yprime)
 186 {
 187   *y = func6 (x, p);
 188   *yprime = 7.0 * pow (x - 1, 6.0);
 189 }
 190
 191 /* sin(x) packaged up nicely. */
 192 double
 193 sin_f (double x, void * p)
 194 {
 195   return sin (x);
 196 }
 197
 198 double
 199 sin_df (double x, void * p)
 200 {
 201   return cos (x);
 202 }
 203
 204 void
 205 sin_fdf (double x, void * p, double *y, double *yprime)
 206 {
 207   *y = sin (x);
 208   *yprime = cos (x);
 209 }
 210
 211 /* cos(x) packaged up nicely. */
 212 double
 213 cos_f (double x, void * p)
 214 {
 215   return cos (x);
 216 }
 217
 218 double
 219 cos_df (double x, void * p)
 220 {
 221   return -sin (x);
 222 }
 223
 224 void
 225 cos_fdf (double x, void * p, double *y, double *yprime)
 226 {
 227   *y = cos (x);
 228   *yprime = -sin (x);
 229 }