1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
|
/* 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 );
}
}
|
