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
75
76
77
78
79
80
81
82
83
84
85
|
/* This program computes the Bernoulli numbers.
* See radd.c for rational arithmetic.
*/
typedef struct{
double n;
double d;
}fract;
#define PD 30
/*
fract x[PD+1] = {0.0};
fract p[PD+1] = {0.0};
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double md_fabs ( double );
extern double md_log10 ( double );
#else
double md_fabs(), md_log10();
#endif
extern double MACHEP;
void bernum_wrap(num, den)
double num[PD-2], den[PD-2];
{
int nx, np;
int i, k, n;
fract s, t;
extern void radd ( fract *, fract *, fract *);
extern void rsub ( fract *, fract *, fract *);
extern void rmul ( fract *, fract *, fract *);
extern void rdiv ( fract *, fract *, fract *);
fract x[PD+1], p[PD+1];
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 );
num[n] = x[n].n;
den[n] = x[n].d;
}
}
|
