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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
/// Simple implementation of the recurrence relation: create an array of Bernoulli numbers
// special cases: n=0 or n=1
10 # Internal'BernoulliArray(n_IsInteger)_(n=0 Or n=1) <-- [
Local(B);
B:=ArrayCreate(n+1,0);
B[1] := 1;
If(n=1, B[2] := -1/2);
B;
];
/// Assume n>=2
20 # Internal'BernoulliArray(n_IsInteger) <-- [
Local(B, i, k, k2, bin);
If (InVerboseMode(), Echo({"Internal'BernoulliArray: using direct recursion, n = ", n}));
B:=ArrayCreate(n+1, 0); // array of B[k], k=1,2,... where B[1] is the 0th Bernoulli number
// it would be better not to store the odd elements but let's optimize this later
// we could also maintain a global cache of Bernoulli numbers computed so far, but it won't really speed up things at large n
// all odd elements after B[2] are zero
B[1] := 1;
B[2] := -1/2;
B[3] := 1/6;
For(i:=4, i<=n, i := i+2) // compute and store B[i]
[ // maintain binomial coefficient
bin := 1; // Bin(i+1,0)
// do not sum over odd elements that are zero anyway - cuts time in half
B[i+1] := 1/2-1/(i+1)*(1 + Sum(k, 1, i/2-1,
[
bin := bin * (i+3-2*k) * (i+2-2*k)/ (2*k-1) / (2*k);
B[2*k+1]*bin; // *Bin(i+1, 2*k)
]
) );
];
B;
];
/// Find the fractional part of Bernoulli number with even index >=2
/// return negative if the sign of the Bernoulli number is negative
BernoulliFracPart(n_IsEven)_(n>=2) <-- [
Local(p, sum);
// always 2 and 3
sum := 1/2+1/3;
// check whether n+1 and n/2+1 are prime
If(IsPrime(n+1), sum := sum+1/(n+1));
If(IsPrime(n/2+1), sum := sum+1/(n/2+1));
// sum over all primes p such that n / p-1 is integer
// enough to check up to n/3 now
For(p:=5, p<=n/3+1, p:=NextPrime(p))
If(Mod(n, p-1)=0, sum := sum + 1/p);
// for negative Bernoulli numbers, let's change sign
// Mod(n/2, 2) is 0 for negative Bernoulli numbers and 1 for positive ones
Div(Numer(sum), Denom(sum)) - sum
+ Mod(n/2,2); // we'll return a negative number if the Bernoulli itself is negative -- slightly against our definitions in the manual
//+ 1; // this would be exactly like the manual says
];
/// Find one Bernoulli number for large index
/// compute Riemann's zeta function and combine with the fractional part
Bernoulli1(n_IsEven)_(n>=2) <-- [
Local(B, prec);
prec := BuiltinPrecisionGet();
// estimate the size of B[n] using Stirling formula
// and compute Ln(B[n])/Ln(10) to find the number of digits
BuiltinPrecisionSet(10);
BuiltinPrecisionSet(
Ceil(N((1/2*Ln(8*Pi*n)-n+n*Ln(n/2/Pi))/Ln(10)))+3 // 3 guard digits
);
If (InVerboseMode(), Echo({"Bernoulli: using zeta funcion, precision ", BuiltinPrecisionSet(), ", n = ", n}));
B := Floor(N( // compute integer part of B
If( // use different methods to compute Zeta function
n>250, // threshold is roughly right for internal math
Internal'ZetaNum2(n, n/17+1), // with this method, a single Bernoulli number n is computed in O(n*M(P)) operations where P = O(n*Ln(n)) is the required precision
// Brent's method requires n^2*P+n*M(P)
// simple array method requires
Internal'ZetaNum1(n, n/17+1) // this gives O(n*Ln(n)*M(P))
)
*N(2*n! /(2*Pi)^n)))
// 2*Pi*e is approx. 17, add 1 to guard precision
* (2*Mod(n/2,2)-1) // sign of B
+ BernoulliFracPart(n); // this already has the right sign
BuiltinPrecisionSet(prec); // restore old precision
B;
];
/// Bernoulli numbers; algorithm from: R. P. Brent, "A FORTRAN multiple-precision arithmetic package", ACM TOMS vol. 4, no. 1, p. 57 (1978).
/// this may be good for floating-point (not exact) evaluation of B[n] at large n
/// but is not good at all for exact evaluation! (too slow)
/// Brent claims that the usual recurrence is numerically unstable
/// but we can't check this because MathPiper internal math is fixed-point and Brent's algorithm needs real floating point (C[k] are very small and then multiplied by (2*k)! )
Internal'BernoulliArray1(n_IsEven) _ (n>=2) <--
[
Local(C, f, k, j, denom, sum);
C := ArrayCreate(n+1, 0);
f := ArrayCreate(n/2, 0);
C[1] := 1;
C[2] := -1/2;
C[3] := 1/12; // C[2*k+1] = B[2*k]/(2*k)!
f[1] := 2; // f[k] = (2k)!
For(k:=2, k<=n/2, k++) // we could start with k=1 but it would be awkward to compute f[] recursively
[
// compute f[k]
f[k] := f[k-1] * (2*k)*(2*k-1);
// compute C[k]
C[2*k+1] := 1/(1-4^(-k))/2*(
[
denom := 4; // = 4^1
sum := 0;
For(j:=1, j<=k-1, j++)
[
sum := sum + C[2*(k-j)+1]/denom/f[j]; // + C[k-j]/(2*j)! /4^j
denom := denom * 4;
];
(2*k-1)/denom/f[k] - sum;
]
);
// Echo({n, k, denom, C[k]});
];
// multiply C's with factorials to get B's
For(k:=1, k<=n/2, k++)
C[2*k+1] := C[2*k+1] * f[k];
// return array object
C;
];
|
