1 /*
2 Copyright (C) 2012 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "arb.h"
13
14 /* With parameter n, the error is bounded by 3/(3+sqrt(8))^n */
15 #define ERROR_A 1.5849625007211561815 /* log2(3) */
16 #define ERROR_B 2.5431066063272239453 /* log2(3+sqrt(8)) */
17
18 typedef struct
19 {
20 arb_t A;
21 arb_t B;
22 arb_t C;
23 arb_t D;
24 arb_t Q1;
25 arb_t Q2;
26 arb_t Q3;
27 }
28 zeta_bsplit_state;
29
30 typedef zeta_bsplit_state zeta_bsplit_t[1];
31
32 static __inline__ void
zeta_bsplit_init(zeta_bsplit_t S)33 zeta_bsplit_init(zeta_bsplit_t S)
34 {
35 arb_init(S->A);
36 arb_init(S->B);
37 arb_init(S->C);
38 arb_init(S->D);
39 arb_init(S->Q1);
40 arb_init(S->Q2);
41 arb_init(S->Q3);
42 }
43
44 static __inline__ void
zeta_bsplit_clear(zeta_bsplit_t S)45 zeta_bsplit_clear(zeta_bsplit_t S)
46 {
47 arb_clear(S->A);
48 arb_clear(S->B);
49 arb_clear(S->C);
50 arb_clear(S->D);
51 arb_clear(S->Q1);
52 arb_clear(S->Q2);
53 arb_clear(S->Q3);
54 }
55
56
57 static __inline__ void
zeta_coeff_k(zeta_bsplit_t S,slong k,slong n,slong s)58 zeta_coeff_k(zeta_bsplit_t S, slong k, slong n, slong s)
59 {
60 arb_set_si(S->D, 2 * (n + k));
61 arb_mul_si(S->D, S->D, n - k, ARF_PREC_EXACT);
62 arb_set_si(S->Q1, k + 1);
63 arb_mul_si(S->Q1, S->Q1, 2*k + 1, ARF_PREC_EXACT);
64
65 if (k == 0)
66 {
67 arb_zero(S->A);
68 arb_one(S->Q2);
69 }
70 else
71 {
72 arb_set_si(S->A, k % 2 ? 1 : -1);
73 arb_mul(S->A, S->A, S->Q1, ARF_PREC_EXACT);
74 arb_ui_pow_ui(S->Q2, k, s, ARF_PREC_EXACT);
75 }
76
77 arb_mul(S->Q3, S->Q1, S->Q2, ARF_PREC_EXACT);
78 arb_zero(S->B);
79 arb_set(S->C, S->Q1);
80 }
81
82 static void
zeta_bsplit(zeta_bsplit_t L,slong a,slong b,slong n,slong s,int cont,slong bits)83 zeta_bsplit(zeta_bsplit_t L, slong a, slong b,
84 slong n, slong s, int cont, slong bits)
85 {
86 if (a + 1 == b)
87 {
88 zeta_coeff_k(L, a, n, s);
89 }
90 else
91 {
92 zeta_bsplit_t R;
93
94 slong m = (a + b) / 2;
95
96 zeta_bsplit(L, m, b, n, s, 1, bits);
97
98 zeta_bsplit_init(R);
99 zeta_bsplit(R, a, m, n, s, 1, bits);
100
101 arb_mul(L->B, L->B, R->D, bits);
102 arb_addmul(L->B, L->A, R->C, bits);
103
104 arb_mul(L->B, L->B, R->Q2, bits);
105 arb_addmul(L->B, R->B, L->Q3, bits);
106
107 arb_mul(L->A, L->A, R->Q3, bits);
108 arb_addmul(L->A, R->A, L->Q3, bits);
109
110 arb_mul(L->C, L->C, R->D, bits);
111 arb_addmul(L->C, R->C, L->Q1, bits);
112
113 if (cont)
114 {
115 arb_mul(L->D, L->D, R->D, bits);
116 arb_mul(L->Q2, L->Q2, R->Q2, bits);
117 }
118
119 arb_mul(L->Q1, L->Q1, R->Q1, bits);
120 arb_mul(L->Q3, L->Q3, R->Q3, bits);
121
122 zeta_bsplit_clear(R);
123 }
124 }
125
126 /* The error for eta(s) is bounded by 3/(3+sqrt(8))^n */
127 void
mag_borwein_error(mag_t err,slong n)128 mag_borwein_error(mag_t err, slong n)
129 {
130 /* upper bound for 1/(3+sqrt(8)) */
131 mag_set_ui_2exp_si(err, 736899889, -32);
132
133 mag_pow_ui(err, err, n);
134 mag_mul_ui(err, err, 3);
135 }
136
137 void
arb_zeta_ui_borwein_bsplit(arb_t x,ulong s,slong prec)138 arb_zeta_ui_borwein_bsplit(arb_t x, ulong s, slong prec)
139 {
140 zeta_bsplit_t sum;
141 mag_t err;
142 slong wp, n;
143
144 /* zeta(0) = -1/2 */
145 if (s == 0)
146 {
147 arb_set_si(x, -1);
148 arb_mul_2exp_si(x, x, -1);
149 return;
150 }
151
152 if (s == 1)
153 {
154 flint_printf("zeta_ui_borwein_bsplit: zeta(1)");
155 flint_abort();
156 }
157
158 n = prec / ERROR_B + 2;
159 wp = prec + 30;
160
161 zeta_bsplit_init(sum);
162 zeta_bsplit(sum, 0, n + 1, n, s, 0, wp);
163
164 /* A/Q3 - B/Q3 / (C/Q1) = (A*C - B*Q1) / (Q3*C) */
165 arb_mul(sum->A, sum->A, sum->C, wp);
166 arb_mul(sum->B, sum->B, sum->Q1, wp);
167 arb_sub(sum->A, sum->A, sum->B, wp);
168 arb_mul(sum->Q3, sum->Q3, sum->C, wp);
169 arb_div(sum->C, sum->A, sum->Q3, wp);
170
171 mag_init(err);
172 mag_borwein_error(err, n);
173 mag_add(arb_radref(sum->C), arb_radref(sum->C), err);
174 mag_clear(err);
175
176 /* convert from eta(s) to zeta(s) */
177 arb_div_2expm1_ui(x, sum->C, s - 1, wp);
178 arb_mul_2exp_si(x, x, s - 1);
179
180 zeta_bsplit_clear(sum);
181 }
182
183