zetac.c
Go to the documentation of this file.
1 /* zetac.c
2  *
3  * Riemann zeta function
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double x, y, zetac();
10  *
11  * y = zetac( x );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  *
18  *
19  * inf.
20  * - -x
21  * zetac(x) = > k , x > 1,
22  * -
23  * k=2
24  *
25  * is related to the Riemann zeta function by
26  *
27  * Riemann zeta(x) = zetac(x) + 1.
28  *
29  * Extension of the function definition for x < 1 is implemented.
30  * Zero is returned for x > log2(INFINITY).
31  *
32  * ACCURACY:
33  *
34  * Tabulated values have full machine accuracy.
35  *
36  * Relative error:
37  * arithmetic domain # trials peak rms
38  * IEEE 1,50 10000 9.8e-16 1.3e-16
39  *
40  *
41  */
42 
43 /*
44  * Cephes Math Library Release 2.1: January, 1989
45  * Copyright 1984, 1987, 1989 by Stephen L. Moshier
46  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
47  */
48 
49 #include "mconf.h"
50 #include "lanczos.h"
51 
52 /* Riemann zeta(x) - 1
53  * for integer arguments between 0 and 30.
54  */
55 static const double azetac[] = {
56  -1.50000000000000000000E0,
57  0.0, /* Not used; zetac(1.0) is infinity. */
58  6.44934066848226436472E-1,
59  2.02056903159594285400E-1,
60  8.23232337111381915160E-2,
61  3.69277551433699263314E-2,
62  1.73430619844491397145E-2,
63  8.34927738192282683980E-3,
64  4.07735619794433937869E-3,
65  2.00839282608221441785E-3,
66  9.94575127818085337146E-4,
67  4.94188604119464558702E-4,
68  2.46086553308048298638E-4,
69  1.22713347578489146752E-4,
70  6.12481350587048292585E-5,
71  3.05882363070204935517E-5,
72  1.52822594086518717326E-5,
73  7.63719763789976227360E-6,
74  3.81729326499983985646E-6,
75  1.90821271655393892566E-6,
76  9.53962033872796113152E-7,
77  4.76932986787806463117E-7,
78  2.38450502727732990004E-7,
79  1.19219925965311073068E-7,
80  5.96081890512594796124E-8,
81  2.98035035146522801861E-8,
82  1.49015548283650412347E-8,
83  7.45071178983542949198E-9,
84  3.72533402478845705482E-9,
85  1.86265972351304900640E-9,
86  9.31327432419668182872E-10
87 };
88 
89 /* 2**x (1 - 1/x) (zeta(x) - 1) = P(1/x)/Q(1/x), 1 <= x <= 10 */
90 static double P[9] = {
91  5.85746514569725319540E11,
92  2.57534127756102572888E11,
93  4.87781159567948256438E10,
94  5.15399538023885770696E9,
95  3.41646073514754094281E8,
96  1.60837006880656492731E7,
97  5.92785467342109522998E5,
98  1.51129169964938823117E4,
99  2.01822444485997955865E2,
100 };
101 
102 static double Q[8] = {
103  /* 1.00000000000000000000E0, */
104  3.90497676373371157516E11,
105  5.22858235368272161797E10,
106  5.64451517271280543351E9,
107  3.39006746015350418834E8,
108  1.79410371500126453702E7,
109  5.66666825131384797029E5,
110  1.60382976810944131506E4,
111  1.96436237223387314144E2,
112 };
113 
114 /* log(zeta(x) - 1 - 2**-x), 10 <= x <= 50 */
115 static double A[11] = {
116  8.70728567484590192539E6,
117  1.76506865670346462757E8,
118  2.60889506707483264896E10,
119  5.29806374009894791647E11,
120  2.26888156119238241487E13,
121  3.31884402932705083599E14,
122  5.13778997975868230192E15,
123  -1.98123688133907171455E15,
124  -9.92763810039983572356E16,
125  7.82905376180870586444E16,
126  9.26786275768927717187E16,
127 };
128 
129 static double B[10] = {
130  /* 1.00000000000000000000E0, */
131  -7.92625410563741062861E6,
132  -1.60529969932920229676E8,
133  -2.37669260975543221788E10,
134  -4.80319584350455169857E11,
135  -2.07820961754173320170E13,
136  -2.96075404507272223680E14,
137  -4.86299103694609136686E15,
138  5.34589509675789930199E15,
139  5.71464111092297631292E16,
140  -1.79915597658676556828E16,
141 };
142 
143 /* (1-x) (zeta(x) - 1), 0 <= x <= 1 */
144 static double R[6] = {
145  -3.28717474506562731748E-1,
146  1.55162528742623950834E1,
147  -2.48762831680821954401E2,
148  1.01050368053237678329E3,
149  1.26726061410235149405E4,
150  -1.11578094770515181334E5,
151 };
152 
153 static double S[5] = {
154  /* 1.00000000000000000000E0, */
155  1.95107674914060531512E1,
156  3.17710311750646984099E2,
157  3.03835500874445748734E3,
158  2.03665876435770579345E4,
159  7.43853965136767874343E4,
160 };
161 
162 static double TAYLOR0[10] = {
163  -1.0000000009110164892,
164  -1.0000000057646759799,
165  -9.9999983138417361078e-1,
166  -1.0000013011460139596,
167  -1.000001940896320456,
168  -9.9987929950057116496e-1,
169  -1.000785194477042408,
170  -1.0031782279542924256,
171  -9.1893853320467274178e-1,
172  -1.5,
173 };
174 
175 #define MAXL2 127
176 #define SQRT_2_PI 0.79788456080286535587989
177 
178 extern double MACHEP;
179 
180 static double zeta_reflection(double);
181 static double zetac_smallneg(double);
182 static double zetac_positive(double);
183 
184 
185 /*
186  * Riemann zeta function, minus one
187  */
188 double zetac(double x)
189 {
190  if (isnan(x)) {
191  return x;
192  }
193  else if (x == -INFINITY) {
194  return NAN;
195  }
196  else if (x < 0.0 && x > -0.01) {
197  return zetac_smallneg(x);
198  }
199  else if (x < 0.0) {
200  return zeta_reflection(-x) - 1;
201  }
202  else {
203  return zetac_positive(x);
204  }
205 }
206 
207 
208 /*
209  * Riemann zeta function
210  */
211 double riemann_zeta(double x)
212 {
213  if (isnan(x)) {
214  return x;
215  }
216  else if (x == -INFINITY) {
217  return NAN;
218  }
219  else if (x < 0.0 && x > -0.01) {
220  return 1 + zetac_smallneg(x);
221  }
222  else if (x < 0.0) {
223  return zeta_reflection(-x);
224  }
225  else {
226  return 1 + zetac_positive(x);
227  }
228 }
229 
230 
231 /*
232  * Compute zetac for positive arguments
233  */
234 static inline double zetac_positive(double x)
235 {
236  int i;
237  double a, b, s, w;
238 
239  if (x == 1.0) {
240  return INFINITY;
241  }
242 
243  if (x >= MAXL2) {
244  /* because first term is 2**-x */
245  return 0.0;
246  }
247 
248  /* Tabulated values for integer argument */
249  w = floor(x);
250  if (w == x) {
251  i = x;
252  if (i < 31) {
253 #ifdef UNK
254  return (azetac[i]);
255 #else
256  return (*(double *) &azetac[4 * i]);
257 #endif
258  }
259  }
260 
261  if (x < 1.0) {
262  w = 1.0 - x;
263  a = polevl(x, R, 5) / (w * p1evl(x, S, 5));
264  return a;
265  }
266 
267  if (x <= 10.0) {
268  b = pow(2.0, x) * (x - 1.0);
269  w = 1.0 / x;
270  s = (x * polevl(w, P, 8)) / (b * p1evl(w, Q, 8));
271  return s;
272  }
273 
274  if (x <= 50.0) {
275  b = pow(2.0, -x);
276  w = polevl(x, A, 10) / p1evl(x, B, 10);
277  w = exp(w) + b;
278  return w;
279  }
280 
281  /* Basic sum of inverse powers */
282  s = 0.0;
283  a = 1.0;
284  do {
285  a += 2.0;
286  b = pow(a, -x);
287  s += b;
288  }
289  while (b / s > MACHEP);
290 
291  b = pow(2.0, -x);
292  s = (s + b) / (1.0 - b);
293  return s;
294 }
295 
296 
297 /*
298  * Compute zetac for small negative x. We can't use the reflection
299  * formula because to double precision 1 - x = 1 and zetac(1) = inf.
300  */
301 static inline double zetac_smallneg(double x)
302 {
303  return polevl(x, TAYLOR0, 9);
304 }
305 
306 
307 /*
308  * Compute zetac using the reflection formula (see DLMF 25.4.2) plus
309  * the Lanczos approximation for Gamma to avoid overflow.
310  */
311 static inline double zeta_reflection(double x)
312 {
313  double base, large_term, small_term, hx, x_shift;
314 
315  hx = x / 2;
316  if (hx == floor(hx)) {
317  /* Hit a zero of the sine factor */
318  return 0;
319  }
320 
321  /* Reduce the argument to sine */
322  x_shift = fmod(x, 4);
323  small_term = -SQRT_2_PI * sin(0.5 * M_PI * x_shift);
324  small_term *= lanczos_sum_expg_scaled(x + 1) * zeta(x + 1, 1);
325 
326  /* Group large terms together to prevent overflow */
327  base = (x + lanczos_g + 0.5) / (2 * M_PI * M_E);
328  large_term = pow(base, x + 0.5);
329  if (isfinite(large_term)) {
330  return large_term * small_term;
331  }
332  /*
333  * We overflowed, but we might be able to stave off overflow by
334  * factoring in the small term earlier. To do this we compute
335  *
336  * (sqrt(large_term) * small_term) * sqrt(large_term)
337  *
338  * Since we only call this method for negative x bounded away from
339  * zero, the small term can only be as small sine on that region;
340  * i.e. about machine epsilon. This means that if the above still
341  * overflows, then there was truly no avoiding it.
342  */
343  large_term = pow(base, 0.5 * x + 0.25);
344  return (large_term * small_term) * large_term;
345 }
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
base
Annotation indicating that a class derives from another given type.
Definition: attr.h:64
M_E
#define M_E
Definition: mconf.h:112
s
RealScalar s
Definition: level1_cplx_impl.h:126
ceres::sin
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
P
static double P[9]
Definition: zetac.c:90
MAXL2
#define MAXL2
Definition: zetac.c:175
SQRT_2_PI
#define SQRT_2_PI
Definition: zetac.c:176
b
Scalar * b
Definition: benchVecAdd.cpp:17
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
B
Definition: test_numpy_dtypes.cpp:299
R
static double R[6]
Definition: zetac.c:144
zetac_positive
static double zetac_positive(double)
Definition: zetac.c:234
isnan
#define isnan(X)
Definition: main.h:93
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
zetac
double zetac(double x)
Definition: zetac.c:188
TAYLOR0
static double TAYLOR0[10]
Definition: zetac.c:162
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
A
Definition: test_numpy_dtypes.cpp:298
azetac
static const double azetac[]
Definition: zetac.c:55
Eigen::bfloat16_impl::fmod
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 fmod(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:567
riemann_zeta
double riemann_zeta(double x)
Definition: zetac.c:211
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
isfinite
#define isfinite(X)
Definition: main.h:95
lanczos.h
lanczos_sum_expg_scaled
double lanczos_sum_expg_scaled(double x)
Definition: lanczos.c:25
Eigen::ArrayBase::pow
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Definition: GlobalFunctions.h:143
Eigen::Quaternion
The quaternion class used to represent 3D orientations and rotations.
Definition: ForwardDeclarations.h:293
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
S
static double S[5]
Definition: zetac.c:153
zeta
double zeta(double x, double q)
Definition: zeta.c:89
zeta_reflection
static double zeta_reflection(double)
Definition: zetac.c:311
zetac_smallneg
static double zetac_smallneg(double)
Definition: zetac.c:301
M_PI
#define M_PI
Definition: mconf.h:117
MACHEP
double MACHEP
Definition: const.c:54
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
floor
const EIGEN_DEVICE_FUNC FloorReturnType floor() const
Definition: ArrayCwiseUnaryOps.h:481
lanczos_g
static const double lanczos_g
Definition: lanczos.h:131


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:43:14