1 /* vim:set ts=4 sw=2 sts=2 et: */
2
3 /* This file was imported from a private scientific library
4 * based on GSL coined Home Scientific Libray (HSL) by its author
5 * Jerome Benoit; this very material is itself inspired from the
6 * material written by G. Jungan and distributed by GSL.
7 * Ultimately, some modifications were done in order to render the
8 * imported material independent from the rest of GSL.
9 */
10
11 /* `hsl/specfunc/hzeta.c' C source file
12 // HSL - Home Scientific Library
13 // Copyright (C) 2017-2018 Jerome Benoit
14 //
15 // HSL is free software; you can redistribute it and/or
16 // modify it under the terms of the GNU General Public License
17 // as published by the Free Software Foundation; either version 2
18 // of the License, or (at your option) any later version.
19 //
20 // This program is distributed in the hope that it will be useful,
21 // but WITHOUT ANY WARRANTY; without even the implied warranty of
22 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 // GNU General Public License for more details.
24 //
25 // You should have received a copy of the GNU General Public License
26 // along with this program; if not, write to the Free Software
27 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
28 */
29
30 /*
31 // The material in this file is mainly inspired by the material written by
32 // G. Jungan and distributed under GPLv2 by the GNU Scientific Library (GSL)
33 // ( https://www.gnu.org/software/gsl/ [specfunc/zeta.c]), itself inspired by
34 // the material written by Moshier and distributed in the Cephes Mathematical
35 // Library ( http://www.moshier.net/ [zeta.c]).
36 //
37 // More specifically, hsl_sf_hzeta_e is a slightly modifed clone of
38 // gsl_sf_hzeta_e as found in GSL 2.4; the remaining is `inspired by'.
39 // [Sooner or later a _Working_Note_ may be deposited at ResearchGate
40 // ( https://www.researchgate.net/profile/Jerome_Benoit )]
41 */
42
43 /* Author: Jerome G. Benoit < jgmbenoit _at_ rezozer _dot_ net > */
44
45 #ifdef _MSC_VER
46 #define _USE_MATH_DEFINES
47 #endif
48
49 #include <math.h>
50 #include <stdio.h>
51 #include "hzeta.h"
52 #include "error.h"
53 #include "platform.h"
54
55 /* imported from gsl_machine.h */
56
57 #define GSL_LOG_DBL_MIN (-7.0839641853226408e+02)
58 #define GSL_LOG_DBL_MAX 7.0978271289338397e+02
59 #define GSL_DBL_EPSILON 2.2204460492503131e-16
60
61 /* imported from gsl_math.h */
62
63 #ifndef M_LOG2E
64 #define M_LOG2E 1.44269504088896340735992468100 /* log_2 (e) */
65 #endif
66
67 /* imported from gsl_sf_result.h */
68
69 struct gsl_sf_result_struct {
70 double val;
71 double err;
72 };
73 typedef struct gsl_sf_result_struct gsl_sf_result;
74
75 /* imported and adapted from hsl/specfunc/specfunc_def.h */
76
77 #define HSL_SF_EVAL_RESULT(FnE) \
78 gsl_sf_result result; \
79 FnE ; \
80 return (result.val);
81
82 #define HSL_SF_EVAL_TUPLE_RESULT(FnET) \
83 gsl_sf_result result0; \
84 gsl_sf_result result1; \
85 FnET ; \
86 *tuple1=result1.val; \
87 *tuple0=result0.val; \
88 return (result0.val);
89
90 /* */
91
92
93 #define HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT 10
94 #define HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER 32
95
96 #define HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX 256
97
98 // B_{2j}/(2j)
99 static
100 double hsl_sf_hzeta_eulermaclaurin_series_coeffs[HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={
101 +1.0,
102 +1.0/12.0,
103 -1.0/720.0,
104 +1.0/30240.0,
105 -1.0/1209600.0,
106 +1.0/47900160.0,
107 -691.0/1307674368000.0,
108 +1.0/74724249600.0,
109 -3.38968029632258286683019539125e-13,
110 +8.58606205627784456413590545043e-15,
111 -2.17486869855806187304151642387e-16,
112 +5.50900282836022951520265260890e-18,
113 -1.39544646858125233407076862641e-19,
114 +3.53470703962946747169322997780e-21,
115 -8.95351742703754685040261131811e-23,
116 +2.26795245233768306031095073887e-24,
117 -5.74479066887220244526388198761e-26,
118 +1.45517247561486490186626486727e-27,
119 -3.68599494066531017818178247991e-29,
120 +9.33673425709504467203255515279e-31,
121 -2.36502241570062993455963519637e-32,
122 +5.99067176248213430465991239682e-34,
123 -1.51745488446829026171081313586e-35,
124 +3.84375812545418823222944529099e-37,
125 -9.73635307264669103526762127925e-39,
126 +2.46624704420068095710640028029e-40,
127 -6.24707674182074369314875679472e-42,
128 +1.58240302446449142975108170683e-43,
129 -4.00827368594893596853001219052e-45,
130 +1.01530758555695563116307139454e-46,
131 -2.57180415824187174992481940976e-48,
132 +6.51445603523381493155843485864e-50,
133 -1.65013099068965245550609878048e-51
134 }; // hsl_sf_hzeta_eulermaclaurin_series_coeffs
135
136 // 4\zeta(2j)/(2\pi)^(2j)
137 static
138 double hsl_sf_hzeta_eulermaclaurin_series_majorantratios[HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={
139 -2.0,
140 +1.0/6.0,
141 +1.0/360.0,
142 +1.0/15120.0,
143 +1.0/604800.0,
144 +1.0/23950080.0,
145 +691.0/653837184000.0,
146 +1.0/37362124800.0,
147 +3617.0/5335311421440000.0,
148 +1.71721241125556891282718109009e-14,
149 +4.34973739711612374608303284773e-16,
150 +1.10180056567204590304053052178e-17,
151 +2.79089293716250466814153725281e-19,
152 +7.06941407925893494338645995561e-21,
153 +1.79070348540750937008052226362e-22,
154 +4.53590490467536612062190147774e-24,
155 +1.14895813377444048905277639752e-25,
156 +2.91034495122972980373252973454e-27,
157 +7.37198988133062035636356495982e-29,
158 +1.86734685141900893440651103056e-30,
159 +4.73004483140125986911927039274e-32,
160 +1.19813435249642686093198247936e-33,
161 +3.03490976893658052342162627173e-35,
162 +7.68751625090837646445889058198e-37,
163 +1.94727061452933820705352425585e-38,
164 +4.93249408840136191421280056051e-40,
165 +1.24941534836414873862975135893e-41,
166 +3.16480604892898285950216341362e-43,
167 +8.01654737189787193706002438098e-45,
168 +2.03061517111391126232614278906e-46,
169 +5.14360831648374349984963881946e-48,
170 +1.30289120704676298631168697172e-49,
171 +3.30026198137930491101219756091e-51
172 }; // hsl_sf_hzeta_eulermaclaurin_series_majorantratios
173
174
175 extern
hsl_sf_hzeta_e(const double s,const double q,gsl_sf_result * result)176 int hsl_sf_hzeta_e(const double s, const double q, gsl_sf_result * result) {
177
178 /* CHECK_POINTER(result) */
179
180 if ((s <= 1.0) || (q <= 0.0)) {
181 PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
182 }
183 else {
184 const double max_bits=54.0; // max_bits=\lceil{s}\rceil with \zeta(s,2)=\zeta(s)-1=GSL_DBL_EPSILON
185 const double ln_term0=-s*log(q);
186 if (ln_term0 < GSL_LOG_DBL_MIN+1.0) {
187 PLFIT_ERROR("underflow", PLFIT_UNDRFLOW);
188 }
189 else if (GSL_LOG_DBL_MAX-1.0 < ln_term0) {
190 PLFIT_ERROR("overflow", PLFIT_OVERFLOW);
191 }
192 #if 1
193 else if (((max_bits < s) && (q < 1.0)) || ((0.5*max_bits < s) && (q < 0.25))) {
194 result->val=pow(q,-s);
195 result->err=2.0*GSL_DBL_EPSILON*fabs(result->val);
196 return (PLFIT_SUCCESS);
197 }
198 else if ((0.5*max_bits < s) && (q < 1.0)) {
199 const double a0=pow(q,-s);
200 const double p1=pow(q/(1.0+q),s);
201 const double p2=pow(q/(2.0+q),s);
202 const double ans=a0*(1.0+p1+p2);
203 result->val=ans;
204 result->err=GSL_DBL_EPSILON*(2.0+0.5*s)*fabs(result->val);
205 return (PLFIT_SUCCESS);
206 }
207 #endif
208 else { // Euler-Maclaurin summation formula
209 const double qshift=HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+q;
210 const double inv_qshift=1.0/qshift;
211 const double sqr_inv_qshift=inv_qshift*inv_qshift;
212 const double inv_sm1=1.0/(s-1.0);
213 const double pmax=pow(qshift,-s);
214 double terms[HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
215 double delta=NAN;
216 double tscp=s;
217 double scp=tscp;
218 double pcp=pmax*inv_qshift;
219 double ratio=scp*pcp;
220 size_t n=0;
221 size_t j=0;
222 double ans=0.0;
223 double mjr=NAN;
224
225 for(j=0;j<HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT;++j) ans+=(terms[n++]=pow(j+q,-s));
226 ans+=(terms[n++]=0.5*pmax);
227 ans+=(terms[n++]=pmax*qshift*inv_sm1);
228 for(j=1;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
229 delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
230 ans+=(terms[n++]=delta);
231 scp*=++tscp;
232 scp*=++tscp;
233 pcp*=sqr_inv_qshift;
234 ratio=scp*pcp;
235 if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
236 }
237
238 ans=0.0; while (n) ans+=terms[--n];
239 mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
240
241 result->val=+ans;
242 result->err=2.0*((HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+1.0)*GSL_DBL_EPSILON*fabs(ans)+mjr);
243 return (PLFIT_SUCCESS);
244 }
245 }
246
247 return (PLFIT_SUCCESS); }
248
249 extern
hsl_sf_hzeta(const double s,const double q)250 double hsl_sf_hzeta(const double s, const double q) {
251 HSL_SF_EVAL_RESULT(hsl_sf_hzeta_e(s,q,&result)); }
252
253 extern
hsl_sf_hzeta_deriv_e(const double s,const double q,gsl_sf_result * result)254 int hsl_sf_hzeta_deriv_e(const double s, const double q, gsl_sf_result * result) {
255
256 /* CHECK_POINTER(result) */
257
258 if ((s <= 1.0) || (q <= 0.0)) {
259 PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
260 }
261 else {
262 const double ln_hz_term0=-s*log(q);
263 if (ln_hz_term0 < GSL_LOG_DBL_MIN+1.0) {
264 PLFIT_ERROR("underflow", PLFIT_UNDRFLOW);
265 }
266 else if (GSL_LOG_DBL_MAX-1.0 < ln_hz_term0) {
267 PLFIT_ERROR("overflow", PLFIT_OVERFLOW);
268 }
269 else { // Euler-Maclaurin summation formula
270 const double qshift=HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+q;
271 const double inv_qshift=1.0/qshift;
272 const double sqr_inv_qshift=inv_qshift*inv_qshift;
273 const double inv_sm1=1.0/(s-1.0);
274 const double pmax=pow(qshift,-s);
275 const double lmax=log(qshift);
276 double terms[HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
277 double delta=NAN;
278 double tscp=s;
279 double scp=tscp;
280 double pcp=pmax*inv_qshift;
281 double lcp=lmax-1.0/s;
282 double ratio=scp*pcp*lcp;
283 double qs=NAN;
284 size_t n=0;
285 size_t j=0;
286 double ans=0.0;
287 double mjr=NAN;
288
289 for(j=0,qs=q;j<HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT;++qs,++j) ans+=(terms[n++]=log(qs)*pow(qs,-s));
290 ans+=(terms[n++]=0.5*lmax*pmax);
291 ans+=(terms[n++]=pmax*qshift*inv_sm1*(lmax+inv_sm1));
292 for(j=1;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
293 delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
294 ans+=(terms[n++]=delta);
295 scp*=++tscp; lcp-=1.0/tscp;
296 scp*=++tscp; lcp-=1.0/tscp;
297 pcp*=sqr_inv_qshift;
298 ratio=scp*pcp*lcp;
299 if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
300 }
301
302 ans=0.0; while (n) ans+=terms[--n];
303 mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
304
305 result->val=-ans;
306 result->err=2.0*((HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+1.0)*GSL_DBL_EPSILON*fabs(ans)+mjr);
307 return (PLFIT_SUCCESS);
308 }
309 }
310
311 return (PLFIT_SUCCESS); }
312
313 extern
hsl_sf_hzeta_deriv(const double s,const double q)314 double hsl_sf_hzeta_deriv(const double s, const double q) {
315 HSL_SF_EVAL_RESULT(hsl_sf_hzeta_deriv_e(s,q,&result)); }
316
317 extern
hsl_sf_hzeta_deriv2_e(const double s,const double q,gsl_sf_result * result)318 int hsl_sf_hzeta_deriv2_e(const double s, const double q, gsl_sf_result * result) {
319
320 /* CHECK_POINTER(result) */
321
322 if ((s <= 1.0) || (q <= 0.0)) {
323 PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
324 }
325 else {
326 const double ln_hz_term0=-s*log(q);
327 if (ln_hz_term0 < GSL_LOG_DBL_MIN+1.0) {
328 PLFIT_ERROR("underflow", PLFIT_UNDRFLOW);
329 }
330 else if (GSL_LOG_DBL_MAX-1.0 < ln_hz_term0) {
331 PLFIT_ERROR("overflow", PLFIT_OVERFLOW);
332 }
333 else { // Euler-Maclaurin summation formula
334 const double qshift=HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+q;
335 const double inv_qshift=1.0/qshift;
336 const double sqr_inv_qshift=inv_qshift*inv_qshift;
337 const double inv_sm1=1.0/(s-1.0);
338 const double pmax=pow(qshift,-s);
339 const double lmax=log(qshift);
340 const double lmax_p_inv_sm1=lmax+inv_sm1;
341 const double sqr_inv_sm1=inv_sm1*inv_sm1;
342 const double sqr_lmax=lmax*lmax;
343 const double sqr_lmax_p_inv_sm1=lmax_p_inv_sm1*lmax_p_inv_sm1;
344 double terms[HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
345 double delta=NAN;
346 double tscp=s;
347 double slcp=NAN;
348 double plcp=NAN;
349 double scp=tscp;
350 double pcp=pmax*inv_qshift;
351 double lcp=1.0/s-lmax;
352 double sqr_lcp=lmax*(lmax-2.0/s);
353 double ratio=scp*pcp*sqr_lcp;
354 double qs=NAN;
355 double lqs=NAN;
356 size_t n=0;
357 size_t j=0;
358 double ans=0.0;
359 double mjr=NAN;
360
361 for(j=0,qs=q;j<HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT;++qs,++j) {
362 lqs=log(qs);
363 ans+=(terms[n++]=lqs*lqs*pow(qs,-s));
364 }
365 ans+=(terms[n++]=0.5*sqr_lmax*pmax);
366 ans+=(terms[n++]=pmax*qshift*inv_sm1*(sqr_lmax_p_inv_sm1+sqr_inv_sm1));
367 for(j=1;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
368 delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
369 ans+=(terms[n++]=delta);
370 scp*=++tscp; slcp=plcp=1.0/tscp;
371 scp*=++tscp; slcp+=1.0/tscp; plcp/=tscp;
372 pcp*=sqr_inv_qshift;
373 sqr_lcp+=2.0*(plcp+slcp*lcp);
374 ratio=scp*pcp*sqr_lcp;
375 if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
376 lcp+=slcp;
377 }
378
379 ans=0.0; while (n) ans+=terms[--n];
380 mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
381
382 result->val=+ans;
383 result->err=2.0*((HSL_SF_HZETA_EULERMACLAURIN_SERIES_SHIFT+1.0)*GSL_DBL_EPSILON*fabs(ans)+mjr);
384 return (PLFIT_SUCCESS);
385 }
386 }
387
388 return (PLFIT_SUCCESS); }
389
390 extern
hsl_sf_hzeta_deriv2(const double s,const double q)391 double hsl_sf_hzeta_deriv2(const double s, const double q) {
392 HSL_SF_EVAL_RESULT(hsl_sf_hzeta_deriv2_e(s,q,&result)); }
393
394 static inline
hsl_sf_hZeta0_zed(const double s,const double q)395 double hsl_sf_hZeta0_zed(const double s, const double q) {
396 #if 1
397 const long double ld_q=(long double)(q);
398 const long double ld_s=(long double)(s);
399 const long double ld_log1prq=log1pl(1.0L/ld_q);
400 const long double ld_epsilon=expm1l(-ld_s*ld_log1prq);
401 const long double ld_z=ld_s+(ld_q+0.5L*ld_s+0.5L)*ld_epsilon;
402 const double z=(double)(ld_z);
403 #else
404 double z=s+(q+0.5*s+0.5)*expm1(-s*log1p(1.0/q));
405 #endif
406 return (z); }
407
408 // Z_{0}(s,a) = a^s \left(\frac{1}{2}+\frac{a}{s-1}\right)^{-1} \zeta(s,a) - 1
409 // Z_{0}(s,a) = O\left(\frac{(s-1)s}{6a^{2}}\right)
410 static
hsl_sf_hZeta0(const double s,const double q,double * value,double * abserror)411 int hsl_sf_hZeta0(const double s, const double q, double * value, double * abserror) {
412 const double criterion=ceil(10.0*s-q);
413 const size_t shift=(criterion<0.0)?0:
414 (criterion<HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX)?(size_t)(llrint(criterion)):
415 HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX;
416 const double qshift=(double)(shift)+q;
417 const double inv_qshift=1.0/qshift;
418 const double sqr_inv_qshift=inv_qshift*inv_qshift;
419 const double sm1=s-1.0;
420 double terms[HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
421 double delta=NAN;
422 double tscp=s;
423 double scp=s*sm1;
424 double pcp=inv_qshift/(2.0*qshift+sm1);
425 double ratio=NAN;
426 size_t n=0;
427 size_t j=0;
428 double ans=0.0;
429 double mjr=NAN;
430
431 if (shift) {
432 const double hsm1=0.5*sm1;
433 const double inv_q=1.0/q;
434 const double qphsm1=q+hsm1;
435 const double inv_qphsm1=1.0/qphsm1;
436 const double qshiftphsm1=qshift+hsm1;
437 double qs=q;
438 double a=1.0;
439 for(j=0;j<shift;) {
440 ans+=(terms[n++]=a*hsl_sf_hZeta0_zed(s,qs++)*inv_qphsm1);
441 a=exp(-s*log1p((++j)*inv_q));
442 }
443 pcp*=a*qshiftphsm1*inv_qphsm1;
444 }
445 ratio=scp*pcp;
446 ans+=terms[n++]=ratio/6.0;
447 scp*=++tscp;
448 scp*=++tscp;
449 pcp*=2.0*sqr_inv_qshift;
450 ratio=scp*pcp;
451 for(j=2;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
452 delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
453 ans+=(terms[n++]=delta);
454 scp*=++tscp;
455 scp*=++tscp;
456 pcp*=sqr_inv_qshift;
457 ratio=scp*pcp;
458 if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
459 }
460
461 ans=0.0; while (n) ans+=terms[--n];
462 mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
463
464 *value=ans;
465 *abserror=2.0*((shift+1)*GSL_DBL_EPSILON*fabs(ans)+mjr);
466
467 return (PLFIT_SUCCESS); }
468
469 static inline
hsl_sf_hZeta1_zed(const double s,const double q)470 double hsl_sf_hZeta1_zed(const double s, const double q) {
471 #if 1
472 const long double ld_q=(long double)(q);
473 const long double ld_s=(long double)(s);
474 const long double ld_sm1=ld_s-1.0L;
475 const long double ld_logq=logl(ld_q);
476 const long double ld_log1prq=log1pl(1.0L/ld_q);
477 const long double ld_inv_logq=1.0L/ld_logq;
478 const long double ld_logratiom1=ld_log1prq*ld_inv_logq;
479 const long double ld_powratiom1=expm1l(-ld_s*ld_log1prq);
480 const long double ld_varepsilon=expm1l(-ld_sm1*ld_log1prq);
481 const long double ld_epsilon=ld_logratiom1+ld_powratiom1+ld_logratiom1*ld_powratiom1;
482 const long double ld_z=ld_s+(ld_q+0.5L*ld_s+0.5L)*ld_epsilon+ld_q/ld_sm1*ld_inv_logq*ld_varepsilon;
483 const double z=(double)(ld_z);
484 #else
485 const double sm1=s-1.0;
486 const double inv_ln_q=1.0/log(q);
487 const double log1prq=log1p(1.0/q);
488 const double logratiom1=log1prq*inv_ln_q;
489 const double powratiom1=expm1(-s*log1prq);
490 const double epsilon=logratiom1+powratiom1+logratiom1*powratiom1;
491 const double z=s+(q+0.5*s+0.5)*epsilon+q/sm1*inv_ln_q*expm1(-sm1*log1prq);
492 #endif
493 return (z); }
494
495 // Z_{1}(s,a) = -\frac{a^s}{\ln(a)} \left(\frac{1}{2}+\frac{a}{s-1}\,\left[1+\frac{1}{(s-1)\,\ln(a)}\right]\right)^{-1} \zeta^{\prime}(s,a) - 1
496 // Z_{1}(s,a) = O\left(\frac{(s-1)s}{6a^{2}}\right)
497 static
hsl_sf_hZeta1(const double s,const double q,const double ln_q,double * value,double * abserror,double * coeff1)498 int hsl_sf_hZeta1(const double s, const double q, const double ln_q, double * value, double * abserror, double * coeff1) {
499 const double criterion=ceil(10.0*s-q);
500 const size_t shift=(criterion<0.0)?0:
501 (criterion<HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX)?(size_t)(llrint(criterion)):
502 HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX;
503 const double qshift=(double)(shift)+q;
504 const double ln_qshift=log(qshift);
505 const double inv_qshift=1.0/qshift;
506 const double inv_ln_q=1.0/ln_q;
507 const double inv_ln_qshift=1.0/ln_qshift;
508 const double sqr_inv_qshift=inv_qshift*inv_qshift;
509 const double sm1=s-1.0;
510 const double hsm1=0.5*sm1;
511 const double q_over_ln_q=q*inv_ln_q;
512 const double qshift_over_ln_qshift=qshift*inv_ln_qshift;
513 const double qphsm1=q+hsm1;
514 const double qshiftphsm1=qshift+hsm1;
515 double terms[HSL_SF_LNHZETA_EULERMACLAURIN_SERIES_SHIFT_MAX+HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER+1]={NAN};
516 double delta=NAN;
517 double tscp=s;
518 double scp=s*sm1;
519 double pcp=inv_qshift*sm1/(qshift_over_ln_qshift+sm1*qshiftphsm1);
520 double lcp=1.0-inv_ln_qshift/s;
521 double ratio=NAN;
522 size_t n=0;
523 size_t j=0;
524 double ans=0.0;
525 double mjr=NAN;
526
527 if (shift) {
528 const double inv_q=1.0/q;
529 const double inv_sm1=1.0/sm1;
530 const double w=1.0+inv_sm1*inv_ln_q;
531 const double wshift=1.0+inv_sm1*inv_ln_qshift;
532 const double qwphsm1=q*w+hsm1;
533 const double inv_qwphsm1=1.0/qwphsm1;
534 const double qshiftwshiftphsm1=qshift*wshift+hsm1;
535 double qs=q;
536 double a=1.0;
537 for(j=0;j<shift;) {
538 ans+=(terms[n++]=a*hsl_sf_hZeta1_zed(s,qs++)*inv_qwphsm1);
539 a=log(qs)*inv_ln_q*exp(-s*log1p((++j)*inv_q));
540 }
541 pcp*=a*qshiftwshiftphsm1*inv_qwphsm1;
542 }
543 ratio=scp*pcp*lcp;
544 ans+=terms[n++]=ratio/12.0;
545 scp*=++tscp; lcp-=inv_ln_qshift/tscp;
546 scp*=++tscp; lcp-=inv_ln_qshift/tscp;
547 pcp*=sqr_inv_qshift;
548 ratio=scp*pcp*lcp;
549 for(j=2;j<=HSL_SF_HZETA_EULERMACLAURIN_SERIES_ORDER;++j) {
550 delta=hsl_sf_hzeta_eulermaclaurin_series_coeffs[j]*ratio;
551 ans+=(terms[n++]=delta);
552 scp*=++tscp; lcp-=inv_ln_qshift/tscp;
553 scp*=++tscp; lcp-=inv_ln_qshift/tscp;
554 pcp*=sqr_inv_qshift;
555 ratio=scp*pcp*lcp;
556 if ((fabs(delta/ans)) < (0.5*GSL_DBL_EPSILON)) break;
557 }
558
559 ans=0.0; while (n) ans+=terms[--n];
560 mjr=hsl_sf_hzeta_eulermaclaurin_series_majorantratios[j]*ratio;
561
562 *value=ans;
563 *abserror=2.0*((shift+1)*GSL_DBL_EPSILON*fabs(ans)+mjr);
564
565 if (coeff1) *coeff1=1.0+q_over_ln_q/qphsm1/sm1;
566
567 return (PLFIT_SUCCESS); }
568
569 extern
hsl_sf_lnhzeta_deriv_tuple_e(const double s,const double q,gsl_sf_result * result,gsl_sf_result * result_deriv)570 int hsl_sf_lnhzeta_deriv_tuple_e(const double s, const double q, gsl_sf_result * result, gsl_sf_result * result_deriv) {
571
572 /* CHECK_POINTER(result) */
573
574 if ((s <= 1.0) || (q <= 0.0)) {
575 PLFIT_ERROR("s must be larger than 1.0 and q must be larger than zero", PLFIT_EINVAL);
576 }
577 else if (q == 1.0) {
578 const double inv_sm1=1.0/(s-1.0);
579 const double inv_qsm1=4.0*inv_sm1;
580 const double hz_coeff0=exp2(s+1.0);
581 const double hz_coeff1=1.0+inv_qsm1;
582 double hZeta0_value=NAN;
583 double hZeta0_abserror=NAN;
584 hsl_sf_hZeta0(s,2.0,&hZeta0_value,&hZeta0_abserror);
585 hZeta0_value+=1.0;
586 if (result) {
587 const double ln_hz_coeff=hz_coeff1/hz_coeff0;
588 const double ln_hZeta0_value=ln_hz_coeff*hZeta0_value;
589 result->val=log1p(ln_hZeta0_value);
590 result->err=(2.0*GSL_DBL_EPSILON*ln_hz_coeff+hZeta0_abserror)/(1.0+ln_hZeta0_value);
591 }
592
593 if (result_deriv) {
594 const double ld_hz_coeff2=1.0+inv_sm1*M_LOG2E;
595 const double ld_hz_coeff1=1.0+inv_qsm1*ld_hz_coeff2;
596 double hZeta1_value=NAN;
597 double hZeta1_abserror=NAN;
598 hsl_sf_hZeta1(s,2.0,M_LN2,&hZeta1_value,&hZeta1_abserror,NULL);
599 hZeta0_value*=hz_coeff1;
600 hZeta0_value+=hz_coeff0;
601 hZeta1_value+=1.0;
602 hZeta1_value*=-M_LN2*ld_hz_coeff1;
603 result_deriv->val=hZeta1_value/hZeta0_value;
604 result_deriv->err=2.0*GSL_DBL_EPSILON*fabs(result_deriv->val)+(hZeta0_abserror+hZeta1_abserror);
605 }
606 }
607 else {
608 const double ln_q=log(q);
609 double hZeta0_value=NAN;
610 double hZeta0_abserror=NAN;
611 hsl_sf_hZeta0(s,q,&hZeta0_value,&hZeta0_abserror);
612 if (result) {
613 const double ln_hz_term0=-s*ln_q;
614 const double ln_hz_term1=log(0.5+q/(s-1.0));
615 result->val=ln_hz_term0+ln_hz_term1+log1p(hZeta0_value);
616 result->err=2.0*GSL_DBL_EPSILON*(fabs(ln_hz_term0)+fabs(ln_hz_term1))+hZeta0_abserror/(1.0+hZeta0_value);
617 }
618 if (result_deriv) {
619 double hZeta1_value=NAN;
620 double hZeta1_abserror=NAN;
621 double ld_hz_coeff1=NAN;
622 hsl_sf_hZeta1(s,q,ln_q,&hZeta1_value,&hZeta1_abserror,&ld_hz_coeff1);
623 result_deriv->val=-ln_q*ld_hz_coeff1*(1.0+hZeta1_value)/(1.0+hZeta0_value);
624 result_deriv->err=2.0*GSL_DBL_EPSILON*fabs(result_deriv->val)+(hZeta0_abserror+hZeta1_abserror);
625 }
626 }
627
628 return (PLFIT_SUCCESS); }
629
630 extern
hsl_sf_lnhzeta_deriv_tuple(const double s,const double q,double * tuple0,double * tuple1)631 double hsl_sf_lnhzeta_deriv_tuple(const double s, const double q, double * tuple0, double * tuple1) {
632 HSL_SF_EVAL_TUPLE_RESULT(hsl_sf_lnhzeta_deriv_tuple_e(s,q,&result0,&result1)); }
633
634 extern
hsl_sf_lnhzeta_e(const double s,const double q,gsl_sf_result * result)635 int hsl_sf_lnhzeta_e(const double s, const double q, gsl_sf_result * result) {
636 return (hsl_sf_lnhzeta_deriv_tuple_e(s,q,result,NULL)); }
637
638 extern
hsl_sf_lnhzeta(const double s,const double q)639 double hsl_sf_lnhzeta(const double s, const double q) {
640 HSL_SF_EVAL_RESULT(hsl_sf_lnhzeta_e(s,q,&result)); }
641
642 extern
hsl_sf_lnhzeta_deriv_e(const double s,const double q,gsl_sf_result * result)643 int hsl_sf_lnhzeta_deriv_e(const double s, const double q, gsl_sf_result * result) {
644 return (hsl_sf_lnhzeta_deriv_tuple_e(s,q,NULL,result)); }
645
646 extern
hsl_sf_lnhzeta_deriv(const double s,const double q)647 double hsl_sf_lnhzeta_deriv(const double s, const double q) {
648 HSL_SF_EVAL_RESULT(hsl_sf_lnhzeta_deriv_e(s,q,&result)); }
649
650 //
651 // End of file `hsl/specfunc/hzeta.c'.
652