Statistics
| Revision:

root / branches / threaded2 / examples / threaded / mforest / qromb.c @ 1099

History | View | Annotate | Download (696 Bytes)

1
#include <math.h>
2
#define EPS 1.0e-6
3
#define JMAX 20
4
#define JMAXP (JMAX+1)
5
#define K 5
6

    
7
float qromb(float (*func)(float), float a, float b)
8
{
9
        void polint(float xa[], float ya[], int n, float x, float *y, float *dy);
10
        float trapzd(float (*func)(float), float a, float b, int n);
11
        void nrerror(char error_text[]);
12
        float ss,dss;
13
        float s[JMAXP],h[JMAXP+1];
14
        int j;
15

    
16
        h[1]=1.0;
17
        for (j=1;j<=JMAX;j++) {
18
                s[j]=trapzd(func,a,b,j);
19
                if (j >= K) {
20
                        polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss);
21
                        if (fabs(dss) <= EPS*fabs(ss)) return ss;
22
                }
23
                h[j+1]=0.25*h[j];
24
        }
25
        nrerror("Too many steps in routine qromb");
26
        return 0.0;
27
}
28
#undef EPS
29
#undef JMAX
30
#undef JMAXP
31
#undef K