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
|
/// special functions coded for MathPiper by Serge Winitzki
/////////////////////////////////////////////////
/// Riemann's Zeta function
/////////////////////////////////////////////////
/// See: Bateman, Erdelyi: <i>Higher Transcendental Functions<i>, vol. 1;
/// P. Borwein, <i>An efficient algorithm for Riemann Zeta function<i> (1995).
/// Numerical computation of Zeta function using Borwein's "third" algorithm
/// The value of $n$ must be large enough to ensure required precision
/// Also $s$ must satisfy $Re(s)+n+1 > 0$
Internal'ZetaNum(_s, n_IsInteger) <-- [
Local(result, j, sign);
If (InVerboseMode(), Echo({"Internal'ZetaNum: Borwein's method, precision ", BuiltinPrecisionGet(), ", n = ", n}));
result := 0;
sign := 1; // flipping sign
For(j:=0, j<=2*n-1, j++)
[ // this is suboptimal b/c we can compute the coefficients a lot faster in this same loop, but ok for now
result := N(result + sign*Internal'ZetaNumCoeffEj(j,n)/(1+j)^s );
sign := -sign;
];
N(result/(2^n)/(1-2^(1-s)));
];
/// direct method -- only good for large s
Internal'ZetaNum1(s, limit) := [
Local(i, sum);
If (InVerboseMode(), Echo({"Internal'ZetaNum: direct method (sum), precision ", BuiltinPrecisionGet(), ", N = ", limit}));
sum := 0;
limit := Ceil(N(limit));
For(i:=2, i<=limit, i++) sum := sum+N(1/PowerN(i, s));
// sum := sum + ( N( 1/PowerN(limit, s-1)) + N(1/PowerN(limit+1, s-1)) )/2/(s-1); // these extra terms don't seem to help much
sum+1; // add small terms together and then add 1
];
/// direct method -- using infinite product. For internal math, Internal'ZetaNum2 is faster for Bernoulli numbers > 250 or so.
Internal'ZetaNum2(s, limit) :=
[
Local(i, prod);
If (InVerboseMode(), Echo({"Internal'ZetaNum: direct method (product), precision ", BuiltinPrecisionGet(), ", N = ", limit}));
prod := N( (1-1/PowerN(2, s))*(1-1/PowerN(3,s)) );
limit := Ceil(N(limit));
For(i:=5, i<=limit, i:= NextPrime(i))
prod := prod*N(1-1/PowerN(i, s));
1/prod;
];
/// Compute coefficients e[j] (see Borwein -- excluding (-1)^j )
Internal'ZetaNumCoeffEj(j,n) := [
Local(k);
2^n-If(j<n,
0,
Sum(k,0,j-n,Bin(n,k)) // this is suboptimal but ok for now
);
];
/// fast numerical calculation of Zeta(3) using a special series
Zeta3() :=
[
Local(result, old'result, k, term);
N([
For(
[
k:=1;
result := 1;
old'result := -1;
term := 1;
],
old'result!=result,
k++
)
[
old'result := result;
term := -term * k^2 / ((2*k+1)*(2*k));
result := result + term/(k+1)^2;
];
result := 5/4*result;
], BuiltinPrecisionGet()+IntLog(BuiltinPrecisionGet(),10)+1);
result;
];
/// User interface for Internal'ZetaNum(s,n)
10 # Internal'ZetaNum(_s) _ (N(s)=0) <-- -0.5;
10 # Internal'ZetaNum(_s) _ (N(s)=1) <-- Infinity;
20 # Internal'ZetaNum(_s) <-- [
Local(n, prec, result);
prec := BuiltinPrecisionGet();
If( // use identity if s<1/2 to replace with larger s. Also must be sn!=0 or else we get infinity * zero
N(Re(s)) < 0.5,
// call ourselves with a different argument
[
If(InVerboseMode(), Echo({"Internal'ZetaNum: using s->1-s identity, s=", s, ", precision ", prec}));
result := 2*Exp(Internal'LnGammaNum(1-s)-(1-s)*Ln(2*Internal'Pi()))*Sin(Internal'Pi()*s/2) * Internal'ZetaNum(1-s);
],
// choose between methods
If (N(Re(s)) > N(1+(prec*Ln(10))/(Ln(prec)+0.1), 6),
[ // use direct summation
n:= N(10^(prec/(s-1)), 6)+2; // 2 guard terms
BuiltinPrecisionSet(prec+2); // 2 guard digits
result := Internal'ZetaNum1(s, n);
],
[ // use Internal'ZetaNum(s, n)
n := Ceil( N( prec*Ln(10)/Ln(8) + 2, 6 ) ); // add 2 digits just in case
BuiltinPrecisionSet(prec+2); // 2 guard digits
result := Internal'ZetaNum(s, n);
]
)
);
BuiltinPrecisionSet(prec);
result;
];
|
