uclibc-ng/libm/double/bernum.c

/* 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 );
}


}