diff options
Diffstat (limited to 'libm/double/bernum.c')
-rw-r--r-- | libm/double/bernum.c | 74 |
1 files changed, 0 insertions, 74 deletions
diff --git a/libm/double/bernum.c b/libm/double/bernum.c deleted file mode 100644 index e401ff5df..000000000 --- a/libm/double/bernum.c +++ /dev/null @@ -1,74 +0,0 @@ -/* This program computes the Bernoulli numbers. - * See radd.c for rational arithmetic. - */ - -typedef struct{ - double n; - double d; - }fract; - -#define PD 44 -fract x[PD+1] = {0.0}; -fract p[PD+1] = {0.0}; -#include <math.h> -#ifdef ANSIPROT -extern double fabs ( double ); -extern double log10 ( double ); -#else -double fabs(), log10(); -#endif -extern double MACHEP; - -main() -{ -int nx, np, nu; -int i, j, k, n, sign; -fract r, s, t; - - -for(i=0; i<=PD; i++ ) - { - x[i].n = 0.0; - x[i].d = 1.0; - p[i].n = 0.0; - p[i].d = 1.0; - } -p[0].n = 1.0; -p[0].d = 1.0; -p[1].n = 1.0; -p[1].d = 1.0; -np = 1; -x[0].n = 1.0; -x[0].d = 1.0; - -for( n=1; n<PD-2; n++ ) -{ - -/* Create line of Pascal's triangle */ -/* multiply p = u * p */ -for( k=0; k<=np; k++ ) - { - radd( &p[np-k+1], &p[np-k], &p[np-k+1] ); - } -np += 1; - -/* B0 + nC1 B1 + ... + nCn-1 Bn-1 = 0 */ -s.n = 0.0; -s.d = 1.0; - -for( i=0; i<n; i++ ) - { - rmul( &p[i], &x[i], &t ); - radd( &s, &t, &s ); - } - - -rdiv( &p[n], &s, &x[n] ); /* x[n] = -s/p[n] */ -x[n].n = -x[n].n; -nx += 1; -printf( "%2d %.15e / %.15e\n", n, x[n].n, x[n].d ); -} - - -} - |