w32tex
About: TeX Live provides a comprehensive TeX system including all the major TeX-related programs, macro packages, and fonts that are free software. Windows sources.
  Fossies Dox: w32tex-src.tar.xz  ("unofficial" and yet experimental doxygen-generated source code documentation)  

zeta.c
Go to the documentation of this file.
1 /* mpfr_zeta -- compute the Riemann Zeta function
2 
3 Copyright 2003-2020 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5 
6 This file is part of the GNU MPFR Library.
7 
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17 
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 
23 #include <float.h> /* for DBL_MAX */
24 
25 #define MPFR_NEED_LONGLONG_H
26 #include "mpfr-impl.h"
27 
28 /*
29  Parameters:
30  s - the input floating-point number
31  n, p - parameters from the algorithm
32  tc - an array of p floating-point numbers tc[1]..tc[p]
33  Output:
34  b is the result, i.e.
35  sum(tc[i]*product((s+2j)*(s+2j-1)/n^2,j=1..i-1), i=1..p)*s*n^(-s-1)
36 */
37 static void
39 {
40  mpfr_t s1, d, u;
41  unsigned long n2;
42  int l, t;
44 
45  if (p == 0)
46  {
47  MPFR_SET_ZERO (b);
48  MPFR_SET_POS (b);
49  return;
50  }
51 
52  n2 = n * n;
54 
55  /* t equals 2p-2, 2p-3, ... ; s1 equals s+t */
56  t = 2 * p - 2;
57  mpfr_set (d, tc[p], MPFR_RNDN);
58  for (l = 1; l < p; l++)
59  {
60  mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l) */
61  mpfr_mul (d, d, s1, MPFR_RNDN);
62  t = t - 1;
63  mpfr_add_ui (s1, s, t, MPFR_RNDN); /* s + (2p-2l-1) */
64  mpfr_mul (d, d, s1, MPFR_RNDN);
65  t = t - 1;
66  mpfr_div_ui (d, d, n2, MPFR_RNDN);
67  mpfr_add (d, d, tc[p-l], MPFR_RNDN);
68  /* since s is positive and the tc[i] have alternate signs,
69  the following is unlikely */
70  if (MPFR_UNLIKELY (mpfr_cmpabs (d, tc[p-l]) > 0))
71  mpfr_set (d, tc[p-l], MPFR_RNDN);
72  }
73  mpfr_mul (d, d, s, MPFR_RNDN);
75  mpfr_neg (s1, s1, MPFR_RNDN);
76  mpfr_ui_pow (u, n, s1, MPFR_RNDN);
77  mpfr_mul (b, d, u, MPFR_RNDN);
78 
80 }
81 
82 /* Input: p - an integer
83  Output: fills tc[1..p], tc[i] = bernoulli(2i)/(2i)!
84  tc[1]=1/12, tc[2]=-1/720, tc[3]=1/30240, ...
85  Assumes all the tc[i] have the same precision.
86 
87  Uses the recurrence (4.60) from the book "Modern Computer Arithmetic"
88  by Brent and Zimmermann for C_k = bernoulli(2k)/(2k)!:
89  sum(C_k/(2k+1-2j)!/4^(k-j), j=0..k) = 1/(2k)!/4^k
90  If we put together the terms involving C_0 and C_1 we get:
91  sum(D_k/(2k+1-2j)!/4^(k-j), j=1..k) = 0
92  with D_1 = C_0/4/(2k+1)/(2k)+C_1-1/(2k)/4=(k-1)/(12k+6),
93  and D_k = C_k for k >= 2.
94 
95  FIXME: we have C_k = (-1)^(k-1) 2/(2pi)^(2k) * zeta(2k),
96  see for example formula (4.65) from the above book,
97  thus since |zeta(2k)-1| < 2^(1-2k) for k >= 2, we have:
98  |C_k - E_k| < E_k * 2^(1-2k) for k >= 2 and E_k := (-1)^(k-1) 2/(2pi)^(2k).
99  Then if 2k-1 >= prec we can evaluate E_k instead, which only requires one
100  multiplication per term, instead of O(k) small divisions.
101 */
102 static void
103 mpfr_zeta_c (int p, mpfr_t *tc)
104 {
105  if (p > 0)
106  {
107  mpfr_t d;
108  int k, l;
109  mpfr_prec_t prec = MPFR_PREC (tc[1]);
110 
111  mpfr_init2 (d, prec);
112  mpfr_div_ui (tc[1], __gmpfr_one, 12, MPFR_RNDN);
113  for (k = 2; k <= p; k++)
114  {
115  mpfr_set_ui (d, k-1, MPFR_RNDN);
116  mpfr_div_ui (d, d, 12*k+6, MPFR_RNDN);
117  for (l=2; l < k; l++)
118  {
119  mpfr_div_ui (d, d, 4*(2*k-2*l+3)*(2*k-2*l+2), MPFR_RNDN);
120  mpfr_add (d, d, tc[l], MPFR_RNDN);
121  }
122  mpfr_div_ui (tc[k], d, 24, MPFR_RNDN);
123  MPFR_CHANGE_SIGN (tc[k]);
124  }
125  mpfr_clear (d);
126  }
127 }
128 
129 /* Input: s - a floating-point number
130  n - an integer
131  Output: sum - a floating-point number approximating sum(1/i^s, i=1..n-1) */
132 static void
134 {
135  mpfr_t u, s1;
136  int i;
138 
139  MPFR_GROUP_INIT_2 (group, MPFR_PREC (sum), u, s1);
140 
141  mpfr_neg (s1, s, MPFR_RNDN);
142  mpfr_ui_pow (u, n, s1, MPFR_RNDN);
143  mpfr_div_2ui (u, u, 1, MPFR_RNDN);
144  mpfr_set (sum, u, MPFR_RNDN);
145  for (i=n-1; i>1; i--)
146  {
147  mpfr_ui_pow (u, i, s1, MPFR_RNDN);
148  mpfr_add (sum, sum, u, MPFR_RNDN);
149  }
150  mpfr_add (sum, sum, __gmpfr_one, MPFR_RNDN);
151 
153 }
154 
155 /* Input: s - a floating-point number >= 1/2.
156  rnd_mode - a rounding mode.
157  Assumes s is neither NaN nor Infinite.
158  Output: z - Zeta(s) rounded to the precision of z with direction rnd_mode
159 */
160 static int
162 {
163  mpfr_t b, c, z_pre, f, s1;
164  double beta, sd, dnep;
165  mpfr_t *tc1;
166  mpfr_prec_t precz, precs, d, dint;
167  int p, n, l, add;
168  int inex;
171 
172  MPFR_ASSERTD (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0);
173 
174  precz = MPFR_PREC (z);
175  precs = MPFR_PREC (s);
176 
177  /* Zeta(x) = 1+1/2^x+1/3^x+1/4^x+1/5^x+O(1/6^x)
178  so with 2^(EXP(x)-1) <= x < 2^EXP(x)
179  So for x > 2^3, k^x > k^8, so 2/k^x < 2/k^8
180  Zeta(x) = 1 + 1/2^x*(1+(2/3)^x+(2/4)^x+...)
181  = 1 + 1/2^x*(1+sum((2/k)^x,k=3..infinity))
182  <= 1 + 1/2^x*(1+sum((2/k)^8,k=3..infinity))
183  And sum((2/k)^8,k=3..infinity) = -257+128*Pi^8/4725 ~= 0.0438035
184  So Zeta(x) <= 1 + 1/2^x*2 for x >= 8
185  The error is < 2^(-x+1) <= 2^(-2^(EXP(x)-1)+1) */
186  if (MPFR_GET_EXP (s) > 3)
187  {
188  mpfr_exp_t err;
189  err = MPFR_GET_EXP (s) - 1;
190  if (err > (mpfr_exp_t) (sizeof (mpfr_exp_t)*CHAR_BIT-2))
191  err = MPFR_EMAX_MAX;
192  else
193  err = ((mpfr_exp_t)1) << err;
194  err = 1 - (-err+1); /* GET_EXP(one) - (-err+1) = err :) */
196  rnd_mode, {});
197  }
198 
199  d = precz + MPFR_INT_CEIL_LOG2(precz) + 10;
200 
201  /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */
202  dint = (mpfr_uexp_t) MPFR_GET_EXP (s);
203  mpfr_init2 (s1, MAX (precs, dint));
204  inex = mpfr_sub (s1, s, __gmpfr_one, MPFR_RNDN);
205  MPFR_ASSERTD (inex == 0);
206 
207  /* case s=1 should have already been handled */
209 
210  MPFR_GROUP_INIT_4 (group, MPFR_PREC_MIN, b, c, z_pre, f);
211 
212  MPFR_ZIV_INIT (loop, d);
213  for (;;)
214  {
215  /* Principal loop: we compute, in z_pre,
216  an approximation of Zeta(s), that we send to can_round */
217  if (MPFR_GET_EXP (s1) <= -(mpfr_exp_t) ((mpfr_prec_t) (d-3)/2))
218  /* Branch 1: when s-1 is very small, one
219  uses the approximation Zeta(s)=1/(s-1)+gamma,
220  where gamma is Euler's constant */
221  {
222  dint = MAX (d + 3, precs);
223  /* branch 1, with internal precision dint */
224  MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
225  mpfr_div (z_pre, __gmpfr_one, s1, MPFR_RNDN);
227  mpfr_add (z_pre, z_pre, f, MPFR_RNDN);
228  }
229  else /* Branch 2 */
230  {
231  size_t size;
232 
233  /* branch 2 */
234  /* Computation of parameters n, p and working precision */
235  dnep = (double) d * LOG2;
236  sd = mpfr_get_d (s, MPFR_RNDN);
237  /* beta = dnep + 0.61 + sd * log (6.2832 / sd);
238  but a larger value is OK */
239 #define LOG6dot2832 1.83787940484160805532
240  beta = dnep + 0.61 + sd * (LOG6dot2832 - LOG2 *
241  __gmpfr_floor_log2 (sd));
242  if (beta <= 0.0)
243  {
244  p = 0;
245  /* n = 1 + (int) (exp ((dnep - LOG2) / sd)); */
246  n = 1 + (int) __gmpfr_ceil_exp2 ((d - 1.0) / sd);
247  }
248  else
249  {
250  p = 1 + (int) beta / 2;
251  n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
252  }
253  /* add = 4 + floor(1.5 * log(d) / log (2)).
254  We should have add >= 10, which is always fulfilled since
255  d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */
256  add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
257  MPFR_ASSERTD(add >= 10);
258  dint = d + add;
259  if (dint < precs)
260  dint = precs;
261 
262  /* internal precision is dint */
263 
264  size = (p + 1) * sizeof(mpfr_t);
265  tc1 = (mpfr_t*) mpfr_allocate_func (size);
266  for (l=1; l<=p; l++)
267  mpfr_init2 (tc1[l], dint);
268  MPFR_GROUP_REPREC_4 (group, dint, b, c, z_pre, f);
269 
270  /* precision of z is precz */
271 
272  /* Computation of the coefficients c_k */
273  mpfr_zeta_c (p, tc1);
274  /* Computation of the 3 parts of the function Zeta. */
275  mpfr_zeta_part_a (z_pre, s, n);
276  mpfr_zeta_part_b (b, s, n, p, tc1);
277  /* s1 = s-1 is already computed above */
279  mpfr_ui_pow (f, n, s1, MPFR_RNDN);
280  mpfr_div (c, c, f, MPFR_RNDN);
281  mpfr_add (z_pre, z_pre, c, MPFR_RNDN);
282  mpfr_add (z_pre, z_pre, b, MPFR_RNDN);
283  for (l=1; l<=p; l++)
284  mpfr_clear (tc1[l]);
285  mpfr_free_func (tc1, size);
286  /* End branch 2 */
287  }
288 
289  if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, d-3, precz, rnd_mode)))
290  break;
291  MPFR_ZIV_NEXT (loop, d);
292  }
294 
295  inex = mpfr_set (z, z_pre, rnd_mode);
296 
298  mpfr_clear (s1);
299 
300  return inex;
301 }
302 
303 /* return add = 1 + floor(log(c^3*(13+m1))/log(2))
304  where c = (1+eps)*(1+eps*max(8,m1)),
305  m1 = 1 + max(1/eps,2*sd)*(1+eps),
306  eps = 2^(-precz-14)
307  sd = abs(s-1)
308  */
309 static long
311 {
312  mpfr_t t, u, m1;
313  long add;
314 
315  mpfr_inits2 (64, t, u, m1, (mpfr_ptr) 0);
316  if (mpfr_cmp_ui (s, 1) >= 0)
317  mpfr_sub_ui (t, s, 1, MPFR_RNDU);
318  else
319  mpfr_ui_sub (t, 1, s, MPFR_RNDU);
320  /* now t = sd = abs(s-1), rounded up */
321  mpfr_set_ui_2exp (u, 1, - precz - 14, MPFR_RNDU);
322  /* u = eps */
323  /* since 1/eps = 2^(precz+14), if EXP(sd) >= precz+14, then
324  sd >= 1/2*2^(precz+14) thus 2*sd >= 2^(precz+14) >= 1/eps */
325  if (mpfr_get_exp (t) >= precz + 14)
326  mpfr_mul_2ui (t, t, 1, MPFR_RNDU);
327  else
328  mpfr_set_ui_2exp (t, 1, precz + 14, MPFR_RNDU);
329  /* now t = max(1/eps,2*sd) */
330  mpfr_add_ui (u, u, 1, MPFR_RNDU); /* u = 1+eps, rounded up */
331  mpfr_mul (t, t, u, MPFR_RNDU); /* t = max(1/eps,2*sd)*(1+eps) */
332  mpfr_add_ui (m1, t, 1, MPFR_RNDU);
333  if (mpfr_get_exp (m1) <= 3)
334  mpfr_set_ui (t, 8, MPFR_RNDU);
335  else
336  mpfr_set (t, m1, MPFR_RNDU);
337  /* now t = max(8,m1) */
338  mpfr_div_2ui (t, t, precz + 14, MPFR_RNDU); /* eps*max(8,m1) */
339  mpfr_add_ui (t, t, 1, MPFR_RNDU); /* 1+eps*max(8,m1) */
340  mpfr_mul (t, t, u, MPFR_RNDU); /* t = c */
341  mpfr_add_ui (u, m1, 13, MPFR_RNDU); /* 13+m1 */
342  mpfr_mul (u, u, t, MPFR_RNDU); /* c*(13+m1) */
343  mpfr_sqr (t, t, MPFR_RNDU); /* c^2 */
344  mpfr_mul (u, u, t, MPFR_RNDU); /* c^3*(13+m1) */
345  add = mpfr_get_exp (u);
346  mpfr_clears (t, u, m1, (mpfr_ptr) 0);
347  return add;
348 }
349 
350 /* return in z a lower bound (for rnd = RNDD) or upper bound (for rnd = RNDU)
351  of |zeta(s)|/2, using:
352  log(|zeta(s)|/2) = (s-1)*log(2*Pi) + lngamma(1-s)
353  + log(|sin(Pi*s/2)| * zeta(1-s)).
354  Assumes s < 1/2 and s1 = 1-s exactly, thus s1 > 1/2.
355  y and p are temporary variables.
356  At input, p is Pi rounded down.
357  The comments in the code are for rnd = RNDD. */
358 static void
361 {
362  mpz_t sint;
363 
365 
366  /* Since log is increasing, we want lower bounds on |sin(Pi*s/2)| and
367  zeta(1-s). */
368  mpz_init (sint);
369  mpfr_get_z (sint, s, MPFR_RNDD); /* sint = floor(s) */
370  /* We first compute a lower bound of |sin(Pi*s/2)|, which is a periodic
371  function of period 2. Thus:
372  if 2k < s < 2k+1, then |sin(Pi*s/2)| is increasing;
373  if 2k-1 < s < 2k, then |sin(Pi*s/2)| is decreasing.
374  These cases are distinguished by testing bit 0 of floor(s) as if
375  represented in two's complement (or equivalently, as an unsigned
376  integer mod 2):
377  0: sint = 0 mod 2, thus 2k < s < 2k+1 and |sin(Pi*s/2)| is increasing;
378  1: sint = 1 mod 2, thus 2k-1 < s < 2k and |sin(Pi*s/2)| is decreasing.
379  Let's recall that the comments are for rnd = RNDD. */
380  if (mpz_tstbit (sint, 0) == 0) /* |sin(Pi*s/2)| is increasing: round down
381  Pi*s to get a lower bound. */
382  {
383  mpfr_mul (y, p, s, rnd);
384  if (rnd == MPFR_RNDD)
385  mpfr_nextabove (p); /* we will need p rounded above afterwards */
386  }
387  else /* |sin(Pi*s/2)| is decreasing: round up Pi*s to get a lower bound. */
388  {
389  if (rnd == MPFR_RNDD)
390  mpfr_nextabove (p);
392  }
393  mpfr_div_2ui (y, y, 1, MPFR_RNDN); /* exact, rounding mode doesn't matter */
394  /* The rounding direction of sin depends on its sign. We have:
395  if -4k-2 < s < -4k, then -2k-1 < s/2 < -2k, thus sin(Pi*s/2) < 0;
396  if -4k < s < -4k+2, then -2k < s/2 < -2k+1, thus sin(Pi*s/2) > 0.
397  These cases are distinguished by testing bit 1 of floor(s) as if
398  represented in two's complement (or equivalently, as an unsigned
399  integer mod 4):
400  0: sint = {0,1} mod 4, thus -2k < s/2 < -2k+1 and sin(Pi*s/2) > 0;
401  1: sint = {2,3} mod 4, thus -2k-1 < s/2 < -2k and sin(Pi*s/2) < 0.
402  Let's recall that the comments are for rnd = RNDD. */
403  if (mpz_tstbit (sint, 1) == 0) /* -2k < s/2 < -2k+1; sin(Pi*s/2) > 0 */
404  {
405  /* Round sin down to get a lower bound of |sin(Pi*s/2)|. */
406  mpfr_sin (y, y, rnd);
407  }
408  else /* -2k-1 < s/2 < -2k; sin(Pi*s/2) < 0 */
409  {
410  /* Round sin up to get a lower bound of |sin(Pi*s/2)|. */
412  mpfr_abs (y, y, MPFR_RNDN); /* exact, rounding mode doesn't matter */
413  }
414  mpz_clear (sint);
415  /* now y <= |sin(Pi*s/2)| when rnd=RNDD, y >= |sin(Pi*s/2)| when rnd=RNDU */
416  mpfr_zeta_pos (z, s1, rnd); /* zeta(1-s) */
417  mpfr_mul (z, z, y, rnd);
418  /* now z <= |sin(Pi*s/2)|*zeta(1-s) */
419  mpfr_log (z, z, rnd);
420  /* now z <= log(|sin(Pi*s/2)|*zeta(1-s)) */
421  mpfr_lngamma (y, s1, rnd);
422  mpfr_add (z, z, y, rnd);
423  /* z <= lngamma(1-s) + log(|sin(Pi*s/2)|*zeta(1-s)) */
424  /* since s-1 < 0, we want to round log(2*pi) upwards */
428  mpfr_sub (z, z, y, rnd);
429  mpfr_exp (z, z, rnd);
430  if (rnd == MPFR_RNDD)
431  mpfr_nextbelow (p); /* restore original p */
432 }
433 
434 int
436 {
437  mpfr_t z_pre, s1, y, p;
438  long add;
439  mpfr_prec_t precz, prec1, precs, precs1;
440  int inex;
443  MPFR_SAVE_EXPO_DECL (expo);
444 
445  MPFR_LOG_FUNC (
446  ("s[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (s), mpfr_log_prec, s, rnd_mode),
447  ("z[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (z), mpfr_log_prec, z, inex));
448 
449  /* Zero, Nan or Inf ? */
451  {
452  if (MPFR_IS_NAN (s))
453  {
454  MPFR_SET_NAN (z);
455  MPFR_RET_NAN;
456  }
457  else if (MPFR_IS_INF (s))
458  {
459  if (MPFR_IS_POS (s))
460  return mpfr_set_ui (z, 1, MPFR_RNDN); /* Zeta(+Inf) = 1 */
461  MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
462  MPFR_RET_NAN;
463  }
464  else /* s iz zero */
465  {
467  return mpfr_set_si_2exp (z, -1, -1, rnd_mode);
468  }
469  }
470 
471  /* s is neither Nan, nor Inf, nor Zero */
472 
473  /* check tiny s: we have zeta(s) = -1/2 - 1/2 log(2 Pi) s + ... around s=0,
474  and for |s| <= 2^(-4), we have |zeta(s) + 1/2| <= |s|.
475  EXP(s) + 1 < -PREC(z) is a sufficient condition to be able to round
476  correctly, for any PREC(z) >= 1 (see algorithms.tex for details). */
477  if (MPFR_GET_EXP (s) + 1 < - (mpfr_exp_t) MPFR_PREC(z))
478  {
479  int signs = MPFR_SIGN(s);
480 
481  MPFR_SAVE_EXPO_MARK (expo);
482  mpfr_set_si_2exp (z, -1, -1, rnd_mode); /* -1/2 */
483  if (rnd_mode == MPFR_RNDA)
484  rnd_mode = MPFR_RNDD; /* the result is around -1/2, thus negative */
485  if ((rnd_mode == MPFR_RNDU || rnd_mode == MPFR_RNDZ) && signs < 0)
486  {
487  mpfr_nextabove (z); /* z = -1/2 + epsilon */
488  inex = 1;
489  }
490  else if (rnd_mode == MPFR_RNDD && signs > 0)
491  {
492  mpfr_nextbelow (z); /* z = -1/2 - epsilon */
493  inex = -1;
494  }
495  else
496  {
497  if (rnd_mode == MPFR_RNDU) /* s > 0: z = -1/2 */
498  inex = 1;
499  else if (rnd_mode == MPFR_RNDD)
500  inex = -1; /* s < 0: z = -1/2 */
501  else /* (MPFR_RNDZ and s > 0) or MPFR_RNDN: z = -1/2 */
502  inex = (signs > 0) ? 1 : -1;
503  }
504  MPFR_SAVE_EXPO_FREE (expo);
505  return mpfr_check_range (z, inex, rnd_mode);
506  }
507 
508  /* Check for case s= -2n */
509  if (MPFR_IS_NEG (s))
510  {
511  mpfr_t tmp;
512  tmp[0] = *s;
513  MPFR_EXP (tmp) = MPFR_GET_EXP (s) - 1;
514  if (mpfr_integer_p (tmp))
515  {
516  MPFR_SET_ZERO (z);
517  MPFR_SET_POS (z);
518  MPFR_RET (0);
519  }
520  }
521 
522  /* Check for case s=1 before changing the exponent range */
523  if (mpfr_cmp (s, __gmpfr_one) == 0)
524  {
525  MPFR_SET_INF (z);
526  MPFR_SET_POS (z);
527  MPFR_SET_DIVBY0 ();
528  MPFR_RET (0);
529  }
530 
531  MPFR_SAVE_EXPO_MARK (expo);
532 
533  /* Compute Zeta */
534  if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
535  inex = mpfr_zeta_pos (z, s, rnd_mode);
536  else /* use reflection formula
537  zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
538  {
539  int overflow = 0;
540 
541  precz = MPFR_PREC (z);
542  precs = MPFR_PREC (s);
543 
544  /* Precision precs1 needed to represent 1 - s, and s + 2,
545  without any truncation */
546  precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
547  /* Precision prec1 is the precision on elementary computations;
548  it ensures a final precision prec1 - add for zeta(s) */
549  add = compute_add (s, precz);
550  prec1 = precz + add;
551  /* FIXME: To avoid that the working precision (prec1) depends on the
552  input precision, one would need to take into account the error made
553  when s1 is not exactly 1-s when computing zeta(s1) and gamma(s1)
554  below, and also in the case y=Inf (i.e. when gamma(s1) overflows).
555  Make sure that underflows do not occur in intermediate computations.
556  Due to the limited precision, they are probably not possible
557  in practice; add some MPFR_ASSERTN's to be sure that problems
558  do not remain undetected? */
559  prec1 = MAX (prec1, precs1) + 10;
560 
561  MPFR_GROUP_INIT_4 (group, prec1, z_pre, s1, y, p);
562  MPFR_ZIV_INIT (loop, prec1);
563  for (;;)
564  {
565  mpfr_exp_t ey;
566  mpfr_t z_up;
567 
568  mpfr_const_pi (p, MPFR_RNDD); /* p is Pi */
569 
570  mpfr_sub (s1, __gmpfr_one, s, MPFR_RNDN); /* s1 = 1-s */
571  mpfr_gamma (y, s1, MPFR_RNDN); /* gamma(1-s) */
572  if (MPFR_IS_INF (y)) /* zeta(s) < 0 for -4k-2 < s < -4k,
573  zeta(s) > 0 for -4k < s < -4k+2 */
574  {
575  /* FIXME: An overflow in gamma(s1) does not imply that
576  zeta(s) will overflow. A solution:
577  1. Compute
578  log(|zeta(s)|/2) = (s-1)*log(2*pi) + lngamma(1-s)
579  + log(abs(sin(Pi*s/2)) * zeta(1-s))
580  (possibly sharing computations with the normal case)
581  with a rather good accuracy (see (2)).
582  Memorize the sign of sin(...) for the final sign.
583  2. Take the exponential, ~= |zeta(s)|/2. If there is an
584  overflow, then this means an overflow on the final result
585  (due to the multiplication by 2, which has not been done
586  yet).
587  3. Ziv test.
588  4. Correct the sign from the sign of sin(...).
589  5. Round then multiply by 2. Here, an overflow in either
590  operation means a real overflow. */
591  mpfr_reflection_overflow (z_pre, s1, s, y, p, MPFR_RNDD);
592  /* z_pre is a lower bound of |zeta(s)|/2, thus if it overflows,
593  or has exponent emax, then |zeta(s)| overflows too. */
594  if (MPFR_IS_INF (z_pre) || MPFR_GET_EXP(z_pre) == __gmpfr_emax)
595  { /* determine the sign of overflow */
596  mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */
597  mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */
598  overflow = (mpfr_cmp_si_2exp (s1, -1, -1) > 0) ? -1 : 1;
599  break;
600  }
601  else /* EXP(z_pre) < __gmpfr_emax */
602  {
603  int ok = 0;
604  mpfr_t z_down;
605  mpfr_init2 (z_up, mpfr_get_prec (z_pre));
607  /* if the lower approximation z_pre does not overflow, but
608  z_up does, we need more precision */
609  if (MPFR_IS_INF (z_up) || MPFR_GET_EXP(z_up) == __gmpfr_emax)
610  goto next_loop;
611  /* check if z_pre and z_up round to the same number */
612  mpfr_init2 (z_down, precz);
613  mpfr_set (z_down, z_pre, rnd_mode);
614  /* Note: it might be that EXP(z_down) = emax here, in that
615  case we will have overflow below when we multiply by 2 */
616  mpfr_prec_round (z_up, precz, rnd_mode);
617  ok = mpfr_cmp (z_down, z_up) == 0;
618  mpfr_clear (z_up);
619  mpfr_clear (z_down);
620  if (ok)
621  {
622  /* get correct sign and multiply by 2 */
623  mpfr_div_2ui (s1, s, 2, MPFR_RNDN); /* s/4, exact */
624  mpfr_frac (s1, s1, MPFR_RNDN); /* exact, -1 < s1 < 0 */
625  if (mpfr_cmp_si_2exp (s1, -1, -1) > 0)
626  mpfr_neg (z_pre, z_pre, rnd_mode);
627  mpfr_mul_2ui (z_pre, z_pre, 1, rnd_mode);
628  break;
629  }
630  else
631  goto next_loop;
632  }
633  }
634  mpfr_zeta_pos (z_pre, s1, MPFR_RNDN); /* zeta(1-s) */
635  mpfr_mul (z_pre, z_pre, y, MPFR_RNDN); /* gamma(1-s)*zeta(1-s) */
636 
637  /* multiply z_pre by 2^s*Pi^(s-1) where p=Pi, s1=1-s */
638  mpfr_mul_2ui (y, p, 1, MPFR_RNDN); /* 2*Pi */
639  mpfr_neg (s1, s1, MPFR_RNDN); /* s-1 */
640  mpfr_pow (y, y, s1, MPFR_RNDN); /* (2*Pi)^(s-1) */
641  mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
642  mpfr_mul_2ui (z_pre, z_pre, 1, MPFR_RNDN);
643 
644  /* multiply z_pre by sin(Pi*s/2) */
645  mpfr_mul (y, s, p, MPFR_RNDN);
646  mpfr_div_2ui (p, y, 1, MPFR_RNDN); /* p = s*Pi/2 */
647  /* FIXME: sinpi will be available, we should replace the mpfr_sin
648  call below by mpfr_sinpi(s/2), where s/2 will be exact.
649  Can mpfr_sin underflow? Moreover, the code below should be
650  improved so that the "if" condition becomes unlikely, e.g.
651  by taking a slightly larger working precision. */
652  mpfr_sin (y, p, MPFR_RNDN); /* y = sin(Pi*s/2) */
653  ey = MPFR_GET_EXP (y);
654  if (ey < 0) /* take account of cancellation in sin(p) */
655  {
656  mpfr_t t;
657 
658  MPFR_ASSERTN (- ey < MPFR_PREC_MAX - prec1);
659  mpfr_init2 (t, prec1 - ey);
661  mpfr_mul (t, s, t, MPFR_RNDN);
662  mpfr_div_2ui (t, t, 1, MPFR_RNDN);
663  mpfr_sin (y, t, MPFR_RNDN);
664  mpfr_clear (t);
665  }
666  mpfr_mul (z_pre, z_pre, y, MPFR_RNDN);
667 
668  if (MPFR_LIKELY (MPFR_CAN_ROUND (z_pre, prec1 - add, precz,
669  rnd_mode)))
670  break;
671 
672  next_loop:
673  MPFR_ZIV_NEXT (loop, prec1);
674  MPFR_GROUP_REPREC_4 (group, prec1, z_pre, s1, y, p);
675  }
677  if (overflow != 0)
678  {
679  inex = mpfr_overflow (z, rnd_mode, overflow);
681  }
682  else
683  inex = mpfr_set (z, z_pre, rnd_mode);
685  }
686 
687  MPFR_SAVE_EXPO_FREE (expo);
688  return mpfr_check_range (z, inex, rnd_mode);
689 }
#define CHAR_BIT
Definition: ChkTeX.h:91
#define mpfr_exp
#define n
Definition: t4ht.c:1290
#define b
Definition: jpegint.h:372
const mpfr_t __gmpfr_one
Definition: constant.c:26
int mpfr_div_2ui(mpfr_ptr y, mpfr_srcptr x, unsigned long n, mpfr_rnd_t rnd_mode)
Definition: div_2ui.c:26
int z
Definition: dviconv.c:26
mpfr_exp_t __gmpfr_emax
Definition: exceptions.c:27
int mpfr_overflow(mpfr_ptr x, mpfr_rnd_t rnd_mode, int sign)
Definition: exceptions.c:406
int mpfr_frac(mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
Definition: frac.c:31
int mpfr_gamma(mpfr_ptr gamma, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
Definition: gamma.c:113
mpz_t * f
Definition: gen-fib.c:34
#define s
Definition: afcover.h:80
#define t
Definition: afcover.h:96
int mpfr_get_z(mpz_ptr z, mpfr_srcptr f, mpfr_rnd_t rnd)
Definition: get_z.c:27
#define mpz_init
Definition: gmp.h:851
#define mpz_clear
Definition: gmp.h:679
#define mpz_tstbit
Definition: gmp.h:1100
__mpz_struct mpz_t[1]
Definition: gmp.h:164
#define c(n)
Definition: gpos-common.c:150
#define d(n)
Definition: gpos-common.c:151
void mpfr_inits2(unsigned va_alist)
Definition: inits2.c:45
int mpfr_integer_p(mpfr_srcptr x)
Definition: isinteger.c:27
#define MAX(a, b)
Definition: jpegint.h:267
small capitals from c petite p scientific f u
Definition: afcover.h:88
small capitals from c petite p
Definition: afcover.h:72
small capitals from c petite p scientific i
Definition: afcover.h:80
kerning y
Definition: ttdriver.c:212
#define beta
Definition: gd_nnquant.c:55
void overflow(const char *)
Definition: cwebboot.c:1385
#define loop
Definition: tie.c:8
int mpfr_lngamma(mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
Definition: lngamma.c:721
int mpfr_log(mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
Definition: log.c:42
void mpfr_free_func(void *ptr, size_t size)
Definition: mpfr-gmp.c:330
void * mpfr_allocate_func(size_t alloc_size)
Definition: mpfr-gmp.c:308
#define MPFR_IS_INF(x)
Definition: mpfr-impl.h:1082
#define MPFR_SET_NAN(x)
Definition: mpfr-impl.h:1081
#define LOG2
Definition: mpfr-impl.h:594
#define MPFR_UNLIKELY(x)
Definition: mpfr-impl.h:1490
#define MPFR_LIKELY(x)
Definition: mpfr-impl.h:1489
#define MPFR_IS_NAN(x)
Definition: mpfr-impl.h:1080
#define MPFR_INVERT_RND(rnd)
Definition: mpfr-impl.h:1224
#define MPFR_INT_CEIL_LOG2(x)
Definition: mpfr-impl.h:1558
#define MPFR_ZIV_INIT(_x, _p)
Definition: mpfr-impl.h:2053
#define MPFR_IS_NEG(x)
Definition: mpfr-impl.h:1140
long __gmpfr_floor_log2(double)
Definition: ufloor_log2.c:27
#define MPFR_EMAX_MAX
Definition: mpfr-impl.h:1034
#define MPFR_SAVE_EXPO_UPDATE_FLAGS(x, flags)
Definition: mpfr-impl.h:1740
#define MPFR_SAVE_EXPO_MARK(x)
Definition: mpfr-impl.h:1730
#define MPFR_FAST_COMPUTE_IF_SMALL_INPUT(y, v, err1, err2, dir, rnd, extra)
Definition: mpfr-impl.h:1981
#define MPFR_PREC(x)
Definition: mpfr-impl.h:958
#define MPFR_SET_POS(x)
Definition: mpfr-impl.h:1143
#define MPFR_GROUP_INIT_2(g, prec, x, y)
Definition: mpfr-impl.h:2296
#define MPFR_GET_EXP(x)
Definition: mpfr-impl.h:1058
#define MPFR_IS_POS(x)
Definition: mpfr-impl.h:1141
#define MPFR_GROUP_REPREC_4(g, prec, x, y, z, t)
Definition: mpfr-impl.h:2348
#define mpfr_const_euler(_d, _r)
Definition: mpfr-impl.h:1386
#define MPFR_SET_ZERO(x)
Definition: mpfr-impl.h:1085
#define MPFR_SAVE_EXPO_FREE(x)
Definition: mpfr-impl.h:1736
#define MPFR_IS_ZERO(x)
Definition: mpfr-impl.h:1084
double __gmpfr_ceil_exp2(double)
Definition: uceil_exp2.c:29
#define MPFR_GROUP_DECL(g)
Definition: mpfr-impl.h:2260
#define MPFR_GROUP_INIT_4(g, prec, x, y, z, t)
Definition: mpfr-impl.h:2303
#define MPFR_ZIV_NEXT(_x, _p)
Definition: mpfr-impl.h:2054
#define MPFR_RET_NAN
Definition: mpfr-impl.h:1184
#define MPFR_RET(I)
Definition: mpfr-impl.h:1182
#define mpfr_check_range(x, t, r)
Definition: mpfr-impl.h:1744
#define MPFR_CAN_ROUND(b, err, prec, rnd)
Definition: mpfr-impl.h:1934
#define MPFR_SET_INF(x)
Definition: mpfr-impl.h:1083
#define mpfr_const_pi(_d, _r)
Definition: mpfr-impl.h:1384
#define MPFR_SAVE_EXPO_DECL(x)
Definition: mpfr-impl.h:1729
#define MPFR_ASSERTN(expr)
Definition: mpfr-impl.h:495
#define MPFR_ZIV_DECL(_x)
Definition: mpfr-impl.h:2052
#define MPFR_ZIV_FREE(x)
Definition: mpfr-impl.h:2055
#define MPFR_LOG_FUNC(x, y)
Definition: mpfr-impl.h:2227
#define MPFR_ASSERTD(expr)
Definition: mpfr-impl.h:516
#define MPFR_EXP(x)
Definition: mpfr-impl.h:959
#define MPFR_CHANGE_SIGN(x)
Definition: mpfr-impl.h:1146
#define MPFR_IS_SINGULAR(x)
Definition: mpfr-impl.h:1100
#define MPFR_GROUP_CLEAR(g)
Definition: mpfr-impl.h:2261
#define MPFR_GROUP_INIT_3(g, prec, x, y, z)
Definition: mpfr-impl.h:2299
#define MPFR_SET_DIVBY0()
Definition: mpfr-impl.h:405
int mpfr_add(mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
Definition: add.c:26
int mpfr_add_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode)
Definition: add_ui.c:27
void mpfr_clear(mpfr_ptr m)
Definition: clear.c:26
void mpfr_clears(unsigned va_alist)
Definition: clears.c:40
int mpfr_cmp_si_2exp(mpfr_srcptr b, long int i, mpfr_exp_t f)
Definition: cmp_si.c:34
int mpfr_cmpabs(mpfr_srcptr b, mpfr_srcptr c)
Definition: cmpabs.c:29
int mpfr_div(mpfr_ptr q, mpfr_srcptr u, mpfr_srcptr v, mpfr_rnd_t rnd_mode)
Definition: div.c:748
int mpfr_div_ui(mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode)
Definition: div_ui.c:33
double mpfr_get_d(mpfr_srcptr src, mpfr_rnd_t rnd_mode)
Definition: get_d.c:36
void mpfr_init2(mpfr_ptr x, mpfr_prec_t p)
Definition: init2.c:26
int mpfr_ui_pow(mpfr_ptr, unsigned long, mpfr_srcptr, mpfr_rnd_t)
int mpfr_pow(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)
Definition: pow.c:380
int mpfr_neg(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)
Definition: neg.c:26
unsigned long mpfr_uexp_t
Definition: mpfr.h:196
int mpfr_sqr(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)
Definition: sqr.c:508
int mpfr_sub_ui(mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t)
long mpfr_prec_t
Definition: mpfr.h:168
void mpfr_nextbelow(mpfr_ptr)
Definition: next.c:121
#define mpfr_set(a, b, r)
Definition: mpfr.h:864
__mpfr_struct mpfr_t[1]
Definition: mpfr.h:250
int mpfr_ui_sub(mpfr_ptr, unsigned long, mpfr_srcptr, mpfr_rnd_t)
int mpfr_set_si_2exp(mpfr_ptr, long, mpfr_exp_t, mpfr_rnd_t)
Definition: set_si_2exp.c:28
#define MPFR_FLAGS_OVERFLOW
Definition: mpfr.h:78
#define mpfr_get_prec(_x)
Definition: mpfr.h:848
#define MPFR_PREC_MAX
Definition: mpfr.h:181
int mpfr_set_ui(mpfr_ptr, unsigned long, mpfr_rnd_t)
Definition: set_ui.c:27
#define mpfr_get_exp(_x)
Definition: mpfr.h:849
int mpfr_mul_2ui(mpfr_ptr, mpfr_srcptr, unsigned long, mpfr_rnd_t)
#define MPFR_SIGN(x)
Definition: mpfr.h:260
#define mpfr_abs(a, b, r)
Definition: mpfr.h:865
#define mpfr_cmp_ui(b, i)
Definition: mpfr.h:862
int mpfr_set_ui_2exp(mpfr_ptr, unsigned long, mpfr_exp_t, mpfr_rnd_t)
Definition: set_ui_2exp.c:28
mpfr_rnd_t
Definition: mpfr.h:102
@ MPFR_RNDA
Definition: mpfr.h:107
@ MPFR_RNDN
Definition: mpfr.h:103
@ MPFR_RNDU
Definition: mpfr.h:105
@ MPFR_RNDD
Definition: mpfr.h:106
@ MPFR_RNDZ
Definition: mpfr.h:104
int mpfr_sin(mpfr_ptr, mpfr_srcptr, mpfr_rnd_t)
Definition: sin.c:37
int mpfr_prec_round(mpfr_ptr, mpfr_prec_t, mpfr_rnd_t)
Definition: round_prec.c:53
#define MPFR_PREC_MIN
Definition: mpfr.h:180
long mpfr_exp_t
Definition: mpfr.h:195
#define mpfr_cmp(b, c)
Definition: mpfr.h:869
void mpfr_nextabove(mpfr_ptr)
Definition: next.c:107
int mpfr_mul(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)
Definition: mul.c:732
int mpfr_sub(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)
Definition: sub.c:26
@ err
Definition: mtxline.h:24
int k
Definition: otp-parser.c:70
set set set set set set set set set set set set set set set set set set set set *set set set macro pixldst op &r &cond WK op &r &cond WK op &r &cond WK else op &m &cond &ia op &r &cond WK else op &m &cond &ia elseif elseif else error unsupported base if elseif elseif else error unsupported unaligned pixldst unaligned endm macro pixst base base else pixldst base endif endm macro PF base if bpp PF set rept prefetch_distance PF set OFFSET endr endif endm macro preload_leading_step2 base if bpp ifc DST PF PF else if bpp lsl PF PF lsl PF PF lsl PF PF PF else PF lsl PF add
static int size
Definition: ppmlabel.c:24
Definition: sd.h:82
Definition: dvips.h:235
s1
Definition: t4ht.c:1059
return() int(((double) *(font_tbl[cur_fnt].wtbl+(int)(*(font_tbl[cur_fnt].char_wi+(int)(ch - font_tbl[cur_fnt].char_f)% 256)))/(double)(1L<< 20)) *(double) font_tbl[cur_fnt].scale)
#define rnd(x)
Definition: xim.h:106
static void mpfr_zeta_part_b(mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
Definition: zeta.c:38
#define LOG6dot2832
static void mpfr_zeta_c(int p, mpfr_t *tc)
Definition: zeta.c:103
static void mpfr_reflection_overflow(mpfr_t z, mpfr_t s1, const mpfr_t s, mpfr_t y, mpfr_t p, mpfr_rnd_t rnd)
Definition: zeta.c:359
int mpfr_zeta(mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode)
Definition: zeta.c:435
static void mpfr_zeta_part_a(mpfr_t sum, mpfr_srcptr s, int n)
Definition: zeta.c:133
static long compute_add(mpfr_srcptr s, mpfr_prec_t precz)
Definition: zeta.c:310
static int mpfr_zeta_pos(mpfr_t z, mpfr_srcptr s, mpfr_rnd_t rnd_mode)
Definition: zeta.c:161