/*
Programmazione II - Algoritmo di Ridders
<![CDATA[
Ridders algorithm implementation (from Numerical Recipes)
By the way: don't emulate this coding style!
*/
#include <stdio.h>
#include <math.h>
#define CON 1.4 /* Stepsize id decreased by CON at each iteration */
#define CON2 (CON*CON)
#define BIG 1.0e50
#define NTAB 10 /* Sets maximum size of tableau */
#define SAFE 2.0 /* Return when error is SAFE worse than the best so far */
double
dfridr(double (*func)(double), double x, double h, double *err)
/* returns the derivative of a function func at a point x by Ridders'
method of polynomial extrapolation. The value h is input as an estimated
initial stepsize; it need not be small, but rather should be an
increment in x over which func changes substantially.
An estimate of the error in the derivative is returned as err */
{
int i,j;
double errt, fac, hh, a[NTAB+1][NTAB+1], ans;
double test1,test2;
if (h == 0)
{
fprintf(stderr,"Third argument to function dfridr must be nonzero.");
return(0);
}
hh = h;
/* Note: this routine assumes a to be a one-based matrix */
a[1][1]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh);
*err = BIG;
for (i=2;i<=NTAB;i++)
{
/* Successive columns in the Neville tableau will go to smaller
stepsizes and higher orders of extrapolation. */
hh /= CON;
/* Try new, smaller stepsize */
a[1][i]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh);
fac=CON2;
for(j=2;j<=i;j++)
{
/* Compute extrapolations of various orders, requiring no
new function evaluations */
a[j][i] = (a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
fac = CON2*fac;
test1 = fabs(a[j][i]-a[j-1][i]);
test2 = fabs(a[j][i]-a[j-1][i-1]);
if (test1>test2) errt = test1;
else errt = test2;
/* The error strategy is to compare each new extrapolation to one
order lower, both at the present stepsize ans the previous one. */
if (errt <= *err)
{
/* If error is decreased, save the improved answer */
*err=errt;
ans=a[j][i];
}
}
/* If higher order is worse by a significant factor SAFE, then quit early */
if (fabs(a[i][i]-a[i-1][i-1]) >= SAFE*(*err)) break;
}
return (ans);
}
]]></plaintext></span>