Coverage for /usr/lib/python3/dist-packages/mpmath/functions/rszeta.py: 2%
893 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2---------------------------------------------------------------------
3.. sectionauthor:: Juan Arias de Reyna <arias@us.es>
5This module implements zeta-related functions using the Riemann-Siegel
6expansion: zeta_offline(s,k=0)
8* coef(J, eps): Need in the computation of Rzeta(s,k)
10* Rzeta_simul(s, der=0) computes Rzeta^(k)(s) and Rzeta^(k)(1-s) simultaneously
11 for 0 <= k <= der. Used by zeta_offline and z_offline
13* Rzeta_set(s, derivatives) computes Rzeta^(k)(s) for given derivatives, used by
14 z_half(t,k) and zeta_half
16* z_offline(w,k): Z(w) and its derivatives of order k <= 4
17* z_half(t,k): Z(t) (Riemann Siegel function) and its derivatives of order k <= 4
18* zeta_offline(s): zeta(s) and its derivatives of order k<= 4
19* zeta_half(1/2+it,k): zeta(s) and its derivatives of order k<= 4
21* rs_zeta(s,k=0) Computes zeta^(k)(s) Unifies zeta_half and zeta_offline
22* rs_z(w,k=0) Computes Z^(k)(w) Unifies z_offline and z_half
23----------------------------------------------------------------------
25This program uses Riemann-Siegel expansion even to compute
26zeta(s) on points s = sigma + i t with sigma arbitrary not
27necessarily equal to 1/2.
29It is founded on a new deduction of the formula, with rigorous
30and sharp bounds for the terms and rest of this expansion.
32More information on the papers:
34 J. Arias de Reyna, High Precision Computation of Riemann's
35 Zeta Function by the Riemann-Siegel Formula I, II
37 We refer to them as I, II.
39 In them we shall find detailed explanation of all the
40 procedure.
42The program uses Riemann-Siegel expansion.
43This is useful when t is big, ( say t > 10000 ).
44The precision is limited, roughly it can compute zeta(sigma+it)
45with an error less than exp(-c t) for some constant c depending
46on sigma. The program gives an error when the Riemann-Siegel
47formula can not compute to the wanted precision.
49"""
51import math
53class RSCache(object):
54 def __init__(ctx):
55 ctx._rs_cache = [0, 10, {}, {}]
57from .functions import defun
59#-------------------------------------------------------------------------------#
60# #
61# coef(ctx, J, eps, _cache=[0, 10, {} ] ) #
62# #
63#-------------------------------------------------------------------------------#
65# This function computes the coefficients c[n] defined on (I, equation (47))
66# but see also (II, section 3.14).
67#
68# Since these coefficients are very difficult to compute we save the values
69# in a cache. So if we compute several values of the functions Rzeta(s) for
70# near values of s, we do not recompute these coefficients.
71#
72# c[n] are the Taylor coefficients of the function:
73#
74# F(z):= (exp(pi*j*(z*z/2+3/8))-j* sqrt(2) cos(pi*z/2))/(2*cos(pi *z))
75#
76#
78def _coef(ctx, J, eps):
79 r"""
80 Computes the coefficients `c_n` for `0\le n\le 2J` with error less than eps
82 **Definition**
84 The coefficients c_n are defined by
86 .. math ::
88 \begin{equation}
89 F(z)=\frac{e^{\pi i
90 \bigl(\frac{z^2}{2}+\frac38\bigr)}-i\sqrt{2}\cos\frac{\pi}{2}z}{2\cos\pi
91 z}=\sum_{n=0}^\infty c_{2n} z^{2n}
92 \end{equation}
94 they are computed applying the relation
96 .. math ::
98 \begin{multline}
99 c_{2n}=-\frac{i}{\sqrt{2}}\Bigl(\frac{\pi}{2}\Bigr)^{2n}
100 \sum_{k=0}^n\frac{(-1)^k}{(2k)!}
101 2^{2n-2k}\frac{(-1)^{n-k}E_{2n-2k}}{(2n-2k)!}+\\
102 +e^{3\pi i/8}\sum_{j=0}^n(-1)^j\frac{
103 E_{2j}}{(2j)!}\frac{i^{n-j}\pi^{n+j}}{(n-j)!2^{n-j+1}}.
104 \end{multline}
105 """
107 newJ = J+2 # compute more coefficients that are needed
108 neweps6 = eps/2. # compute with a slight more precision that are needed
110 # PREPARATION FOR THE COMPUTATION OF V(N) AND W(N)
111 # See II Section 3.16
112 #
113 # Computing the exponent wpvw of the error II equation (81)
114 wpvw = max(ctx.mag(10*(newJ+3)), 4*newJ+5-ctx.mag(neweps6))
116 # Preparation of Euler numbers (we need until the 2*RS_NEWJ)
117 E = ctx._eulernum(2*newJ)
119 # Now we have in the cache all the needed Euler numbers.
120 #
121 # Computing the powers of pi
122 #
123 # We need to compute the powers pi**n for 1<= n <= 2*J
124 # with relative error less than 2**(-wpvw)
125 # it is easy to show that this is obtained
126 # taking wppi as the least d with
127 # 2**d>40*J and 2**d> 4.24 *newJ + 2**wpvw
128 # In II Section 3.9 we need also that
129 # wppi > wptcoef[0], and that the powers
130 # here computed 0<= k <= 2*newJ are more
131 # than those needed there that are 2*L-2.
132 # so we need J >= L this will be checked
133 # before computing tcoef[]
134 wppi = max(ctx.mag(40*newJ), ctx.mag(newJ)+3 +wpvw)
135 ctx.prec = wppi
136 pipower = {}
137 pipower[0] = ctx.one
138 pipower[1] = ctx.pi
139 for n in range(2,2*newJ+1):
140 pipower[n] = pipower[n-1]*ctx.pi
142 # COMPUTING THE COEFFICIENTS v(n) AND w(n)
143 # see II equation (61) and equations (81) and (82)
144 ctx.prec = wpvw+2
145 v={}
146 w={}
147 for n in range(0,newJ+1):
148 va = (-1)**n * ctx._eulernum(2*n)
149 va = ctx.mpf(va)/ctx.fac(2*n)
150 v[n]=va*pipower[2*n]
151 for n in range(0,2*newJ+1):
152 wa = ctx.one/ctx.fac(n)
153 wa=wa/(2**n)
154 w[n]=wa*pipower[n]
156 # COMPUTATION OF THE CONVOLUTIONS RS_P1 AND RS_P2
157 # See II Section 3.16
158 ctx.prec = 15
159 wpp1a = 9 - ctx.mag(neweps6)
160 P1 = {}
161 for n in range(0,newJ+1):
162 ctx.prec = 15
163 wpp1 = max(ctx.mag(10*(n+4)),4*n+wpp1a)
164 ctx.prec = wpp1
165 sump = 0
166 for k in range(0,n+1):
167 sump += ((-1)**k) * v[k]*w[2*n-2*k]
168 P1[n]=((-1)**(n+1))*ctx.j*sump
169 P2={}
170 for n in range(0,newJ+1):
171 ctx.prec = 15
172 wpp2 = max(ctx.mag(10*(n+4)),4*n+wpp1a)
173 ctx.prec = wpp2
174 sump = 0
175 for k in range(0,n+1):
176 sump += (ctx.j**(n-k)) * v[k]*w[n-k]
177 P2[n]=sump
178 # COMPUTING THE COEFFICIENTS c[2n]
179 # See II Section 3.14
180 ctx.prec = 15
181 wpc0 = 5 - ctx.mag(neweps6)
182 wpc = max(6,4*newJ+wpc0)
183 ctx.prec = wpc
184 mu = ctx.sqrt(ctx.mpf('2'))/2
185 nu = ctx.expjpi(3./8)/2
186 c={}
187 for n in range(0,newJ):
188 ctx.prec = 15
189 wpc = max(6,4*n+wpc0)
190 ctx.prec = wpc
191 c[2*n] = mu*P1[n]+nu*P2[n]
192 for n in range(1,2*newJ,2):
193 c[n] = 0
194 return [newJ, neweps6, c, pipower]
196def coef(ctx, J, eps):
197 _cache = ctx._rs_cache
198 if J <= _cache[0] and eps >= _cache[1]:
199 return _cache[2], _cache[3]
200 orig = ctx._mp.prec
201 try:
202 data = _coef(ctx._mp, J, eps)
203 finally:
204 ctx._mp.prec = orig
205 if ctx is not ctx._mp:
206 data[2] = dict((k,ctx.convert(v)) for (k,v) in data[2].items())
207 data[3] = dict((k,ctx.convert(v)) for (k,v) in data[3].items())
208 ctx._rs_cache[:] = data
209 return ctx._rs_cache[2], ctx._rs_cache[3]
211#-------------------------------------------------------------------------------#
212# #
213# Rzeta_simul(s,k=0) #
214# #
215#-------------------------------------------------------------------------------#
216# This function return a list with the values:
217# Rzeta(sigma+it), conj(Rzeta(1-sigma+it)),Rzeta'(sigma+it), conj(Rzeta'(1-sigma+it)),
218# .... , Rzeta^{(k)}(sigma+it), conj(Rzeta^{(k)}(1-sigma+it))
219#
220# Useful to compute the function zeta(s) and Z(w) or its derivatives.
221#
223def aux_M_Fp(ctx, xA, xeps4, a, xB1, xL):
224 # COMPUTING M NUMBER OF DERIVATIVES Fp[m] TO COMPUTE
225 # See II Section 3.11 equations (47) and (48)
226 aux1 = 126.0657606*xA/xeps4 # 126.06.. = 316/sqrt(2*pi)
227 aux1 = ctx.ln(aux1)
228 aux2 = (2*ctx.ln(ctx.pi)+ctx.ln(xB1)+ctx.ln(a))/3 -ctx.ln(2*ctx.pi)/2
229 m = 3*xL-3
230 aux3= (ctx.loggamma(m+1)-ctx.loggamma(m/3.0+2))/2 -ctx.loggamma((m+1)/2.)
231 while((aux1 < m*aux2+ aux3)and (m>1)):
232 m = m - 1
233 aux3 = (ctx.loggamma(m+1)-ctx.loggamma(m/3.0+2))/2 -ctx.loggamma((m+1)/2.)
234 xM = m
235 return xM
237def aux_J_needed(ctx, xA, xeps4, a, xB1, xM):
238 # DETERMINATION OF J THE NUMBER OF TERMS NEEDED
239 # IN THE TAYLOR SERIES OF F.
240 # See II Section 3.11 equation (49))
241 # Only determine one
242 h1 = xeps4/(632*xA)
243 h2 = xB1*a * 126.31337419529260248 # = pi^2*e^2*sqrt(3)
244 h2 = h1 * ctx.power((h2/xM**2),(xM-1)/3) / xM
245 h3 = min(h1,h2)
246 return h3
248def Rzeta_simul(ctx, s, der=0):
249 # First we take the value of ctx.prec
250 wpinitial = ctx.prec
252 # INITIALIZATION
253 # Take the real and imaginary part of s
254 t = ctx._im(s)
255 xsigma = ctx._re(s)
256 ysigma = 1 - xsigma
258 # Now compute several parameter that appear on the program
259 ctx.prec = 15
260 a = ctx.sqrt(t/(2*ctx.pi))
261 xasigma = a ** xsigma
262 yasigma = a ** ysigma
264 # We need a simple bound A1 < asigma (see II Section 3.1 and 3.3)
265 xA1=ctx.power(2, ctx.mag(xasigma)-1)
266 yA1=ctx.power(2, ctx.mag(yasigma)-1)
268 # We compute various epsilon's (see II end of Section 3.1)
269 eps = ctx.power(2, -wpinitial)
270 eps1 = eps/6.
271 xeps2 = eps * xA1/3.
272 yeps2 = eps * yA1/3.
274 # COMPUTING SOME COEFFICIENTS THAT DEPENDS
275 # ON sigma
276 # constant b and c (see I Theorem 2 formula (26) )
277 # coefficients A and B1 (see I Section 6.1 equation (50))
278 #
279 # here we not need high precision
280 ctx.prec = 15
281 if xsigma > 0:
282 xb = 2.
283 xc = math.pow(9,xsigma)/4.44288
284 # 4.44288 =(math.sqrt(2)*math.pi)
285 xA = math.pow(9,xsigma)
286 xB1 = 1
287 else:
288 xb = 2.25158 # math.sqrt( (3-2* math.log(2))*math.pi )
289 xc = math.pow(2,-xsigma)/4.44288
290 xA = math.pow(2,-xsigma)
291 xB1 = 1.10789 # = 2*sqrt(1-log(2))
293 if(ysigma > 0):
294 yb = 2.
295 yc = math.pow(9,ysigma)/4.44288
296 # 4.44288 =(math.sqrt(2)*math.pi)
297 yA = math.pow(9,ysigma)
298 yB1 = 1
299 else:
300 yb = 2.25158 # math.sqrt( (3-2* math.log(2))*math.pi )
301 yc = math.pow(2,-ysigma)/4.44288
302 yA = math.pow(2,-ysigma)
303 yB1 = 1.10789 # = 2*sqrt(1-log(2))
305 # COMPUTING L THE NUMBER OF TERMS NEEDED IN THE RIEMANN-SIEGEL
306 # CORRECTION
307 # See II Section 3.2
308 ctx.prec = 15
309 xL = 1
310 while 3*xc*ctx.gamma(xL*0.5) * ctx.power(xb*a,-xL) >= xeps2:
311 xL = xL+1
312 xL = max(2,xL)
313 yL = 1
314 while 3*yc*ctx.gamma(yL*0.5) * ctx.power(yb*a,-yL) >= yeps2:
315 yL = yL+1
316 yL = max(2,yL)
318 # The number L has to satify some conditions.
319 # If not RS can not compute Rzeta(s) with the prescribed precision
320 # (see II, Section 3.2 condition (20) ) and
321 # (II, Section 3.3 condition (22) ). Also we have added
322 # an additional technical condition in Section 3.17 Proposition 17
323 if ((3*xL >= 2*a*a/25.) or (3*xL+2+xsigma<0) or (abs(xsigma) > a/2.) or \
324 (3*yL >= 2*a*a/25.) or (3*yL+2+ysigma<0) or (abs(ysigma) > a/2.)):
325 ctx.prec = wpinitial
326 raise NotImplementedError("Riemann-Siegel can not compute with such precision")
328 # We take the maximum of the two values
329 L = max(xL, yL)
331 # INITIALIZATION (CONTINUATION)
332 #
333 # eps3 is the constant defined on (II, Section 3.5 equation (27) )
334 # each term of the RS correction must be computed with error <= eps3
335 xeps3 = xeps2/(4*xL)
336 yeps3 = yeps2/(4*yL)
338 # eps4 is defined on (II Section 3.6 equation (30) )
339 # each component of the formula (II Section 3.6 equation (29) )
340 # must be computed with error <= eps4
341 xeps4 = xeps3/(3*xL)
342 yeps4 = yeps3/(3*yL)
344 # COMPUTING M NUMBER OF DERIVATIVES Fp[m] TO COMPUTE
345 xM = aux_M_Fp(ctx, xA, xeps4, a, xB1, xL)
346 yM = aux_M_Fp(ctx, yA, yeps4, a, yB1, yL)
347 M = max(xM, yM)
349 # COMPUTING NUMBER OF TERMS J NEEDED
350 h3 = aux_J_needed(ctx, xA, xeps4, a, xB1, xM)
351 h4 = aux_J_needed(ctx, yA, yeps4, a, yB1, yM)
352 h3 = min(h3,h4)
353 J = 12
354 jvalue = (2*ctx.pi)**J / ctx.gamma(J+1)
355 while jvalue > h3:
356 J = J+1
357 jvalue = (2*ctx.pi)*jvalue/J
359 # COMPUTING eps5[m] for 1 <= m <= 21
360 # See II Section 10 equation (43)
361 # We choose the minimum of the two possibilities
362 eps5={}
363 xforeps5 = math.pi*math.pi*xB1*a
364 yforeps5 = math.pi*math.pi*yB1*a
365 for m in range(0,22):
366 xaux1 = math.pow(xforeps5, m/3)/(316.*xA)
367 yaux1 = math.pow(yforeps5, m/3)/(316.*yA)
368 aux1 = min(xaux1, yaux1)
369 aux2 = ctx.gamma(m+1)/ctx.gamma(m/3.0+0.5)
370 aux2 = math.sqrt(aux2)
371 eps5[m] = (aux1*aux2*min(xeps4,yeps4))
373 # COMPUTING wpfp
374 # See II Section 3.13 equation (59)
375 twenty = min(3*L-3, 21)+1
376 aux = 6812*J
377 wpfp = ctx.mag(44*J)
378 for m in range(0,twenty):
379 wpfp = max(wpfp, ctx.mag(aux*ctx.gamma(m+1)/eps5[m]))
381 # COMPUTING N AND p
382 # See II Section
383 ctx.prec = wpfp + ctx.mag(t)+20
384 a = ctx.sqrt(t/(2*ctx.pi))
385 N = ctx.floor(a)
386 p = 1-2*(a-N)
388 # now we get a rounded version of p
389 # to the precision wpfp
390 # this possibly is not necessary
391 num=ctx.floor(p*(ctx.mpf('2')**wpfp))
392 difference = p * (ctx.mpf('2')**wpfp)-num
393 if (difference < 0.5):
394 num = num
395 else:
396 num = num+1
397 p = ctx.convert(num * (ctx.mpf('2')**(-wpfp)))
399 # COMPUTING THE COEFFICIENTS c[n] = cc[n]
400 # We shall use the notation cc[n], since there is
401 # a constant that is called c
402 # See II Section 3.14
403 # We compute the coefficients and also save then in a
404 # cache. The bulk of the computation is passed to
405 # the function coef()
406 #
407 # eps6 is defined in II Section 3.13 equation (58)
408 eps6 = ctx.power(ctx.convert(2*ctx.pi), J)/(ctx.gamma(J+1)*3*J)
410 # Now we compute the coefficients
411 cc = {}
412 cont = {}
413 cont, pipowers = coef(ctx, J, eps6)
414 cc=cont.copy() # we need a copy since we have to change his values.
415 Fp={} # this is the adequate locus of this
416 for n in range(M, 3*L-2):
417 Fp[n] = 0
418 Fp={}
419 ctx.prec = wpfp
420 for m in range(0,M+1):
421 sumP = 0
422 for k in range(2*J-m-1,-1,-1):
423 sumP = (sumP * p)+ cc[k]
424 Fp[m] = sumP
425 # preparation of the new coefficients
426 for k in range(0,2*J-m-1):
427 cc[k] = (k+1)* cc[k+1]
429 # COMPUTING THE NUMBERS xd[u,n,k], yd[u,n,k]
430 # See II Section 3.17
431 #
432 # First we compute the working precisions xwpd[k]
433 # Se II equation (92)
434 xwpd={}
435 d1 = max(6,ctx.mag(40*L*L))
436 xd2 = 13+ctx.mag((1+abs(xsigma))*xA)-ctx.mag(xeps4)-1
437 xconst = ctx.ln(8/(ctx.pi*ctx.pi*a*a*xB1*xB1)) /2
438 for n in range(0,L):
439 xd3 = ctx.mag(ctx.sqrt(ctx.gamma(n-0.5)))-ctx.floor(n*xconst)+xd2
440 xwpd[n]=max(xd3,d1)
442 # procedure of II Section 3.17
443 ctx.prec = xwpd[1]+10
444 xpsigma = 1-(2*xsigma)
445 xd = {}
446 xd[0,0,-2]=0; xd[0,0,-1]=0; xd[0,0,0]=1; xd[0,0,1]=0
447 xd[0,-1,-2]=0; xd[0,-1,-1]=0; xd[0,-1,0]=1; xd[0,-1,1]=0
448 for n in range(1,L):
449 ctx.prec = xwpd[n]+10
450 for k in range(0,3*n//2+1):
451 m = 3*n-2*k
452 if(m!=0):
453 m1 = ctx.one/m
454 c1= m1/4
455 c2=(xpsigma*m1)/2
456 c3=-(m+1)
457 xd[0,n,k]=c3*xd[0,n-1,k-2]+c1*xd[0,n-1,k]+c2*xd[0,n-1,k-1]
458 else:
459 xd[0,n,k]=0
460 for r in range(0,k):
461 add=xd[0,n,r]*(ctx.mpf('1.0')*ctx.fac(2*k-2*r)/ctx.fac(k-r))
462 xd[0,n,k] -= ((-1)**(k-r))*add
463 xd[0,n,-2]=0; xd[0,n,-1]=0; xd[0,n,3*n//2+1]=0
464 for mu in range(-2,der+1):
465 for n in range(-2,L):
466 for k in range(-3,max(1,3*n//2+2)):
467 if( (mu<0)or (n<0) or(k<0)or (k>3*n//2)):
468 xd[mu,n,k] = 0
469 for mu in range(1,der+1):
470 for n in range(0,L):
471 ctx.prec = xwpd[n]+10
472 for k in range(0,3*n//2+1):
473 aux=(2*mu-2)*xd[mu-2,n-2,k-3]+2*(xsigma+n-2)*xd[mu-1,n-2,k-3]
474 xd[mu,n,k] = aux - xd[mu-1,n-1,k-1]
476 # Now we compute the working precisions ywpd[k]
477 # Se II equation (92)
478 ywpd={}
479 d1 = max(6,ctx.mag(40*L*L))
480 yd2 = 13+ctx.mag((1+abs(ysigma))*yA)-ctx.mag(yeps4)-1
481 yconst = ctx.ln(8/(ctx.pi*ctx.pi*a*a*yB1*yB1)) /2
482 for n in range(0,L):
483 yd3 = ctx.mag(ctx.sqrt(ctx.gamma(n-0.5)))-ctx.floor(n*yconst)+yd2
484 ywpd[n]=max(yd3,d1)
486 # procedure of II Section 3.17
487 ctx.prec = ywpd[1]+10
488 ypsigma = 1-(2*ysigma)
489 yd = {}
490 yd[0,0,-2]=0; yd[0,0,-1]=0; yd[0,0,0]=1; yd[0,0,1]=0
491 yd[0,-1,-2]=0; yd[0,-1,-1]=0; yd[0,-1,0]=1; yd[0,-1,1]=0
492 for n in range(1,L):
493 ctx.prec = ywpd[n]+10
494 for k in range(0,3*n//2+1):
495 m = 3*n-2*k
496 if(m!=0):
497 m1 = ctx.one/m
498 c1= m1/4
499 c2=(ypsigma*m1)/2
500 c3=-(m+1)
501 yd[0,n,k]=c3*yd[0,n-1,k-2]+c1*yd[0,n-1,k]+c2*yd[0,n-1,k-1]
502 else:
503 yd[0,n,k]=0
504 for r in range(0,k):
505 add=yd[0,n,r]*(ctx.mpf('1.0')*ctx.fac(2*k-2*r)/ctx.fac(k-r))
506 yd[0,n,k] -= ((-1)**(k-r))*add
507 yd[0,n,-2]=0; yd[0,n,-1]=0; yd[0,n,3*n//2+1]=0
509 for mu in range(-2,der+1):
510 for n in range(-2,L):
511 for k in range(-3,max(1,3*n//2+2)):
512 if( (mu<0)or (n<0) or(k<0)or (k>3*n//2)):
513 yd[mu,n,k] = 0
514 for mu in range(1,der+1):
515 for n in range(0,L):
516 ctx.prec = ywpd[n]+10
517 for k in range(0,3*n//2+1):
518 aux=(2*mu-2)*yd[mu-2,n-2,k-3]+2*(ysigma+n-2)*yd[mu-1,n-2,k-3]
519 yd[mu,n,k] = aux - yd[mu-1,n-1,k-1]
521 # COMPUTING THE COEFFICIENTS xtcoef[k,l]
522 # See II Section 3.9
523 #
524 # computing the needed wp
525 xwptcoef={}
526 xwpterm={}
527 ctx.prec = 15
528 c1 = ctx.mag(40*(L+2))
529 xc2 = ctx.mag(68*(L+2)*xA)
530 xc4 = ctx.mag(xB1*a*math.sqrt(ctx.pi))-1
531 for k in range(0,L):
532 xc3 = xc2 - k*xc4+ctx.mag(ctx.fac(k+0.5))/2.
533 xwptcoef[k] = (max(c1,xc3-ctx.mag(xeps4)+1)+1 +20)*1.5
534 xwpterm[k] = (max(c1,ctx.mag(L+2)+xc3-ctx.mag(xeps3)+1)+1 +20)
535 ywptcoef={}
536 ywpterm={}
537 ctx.prec = 15
538 c1 = ctx.mag(40*(L+2))
539 yc2 = ctx.mag(68*(L+2)*yA)
540 yc4 = ctx.mag(yB1*a*math.sqrt(ctx.pi))-1
541 for k in range(0,L):
542 yc3 = yc2 - k*yc4+ctx.mag(ctx.fac(k+0.5))/2.
543 ywptcoef[k] = ((max(c1,yc3-ctx.mag(yeps4)+1))+10)*1.5
544 ywpterm[k] = (max(c1,ctx.mag(L+2)+yc3-ctx.mag(yeps3)+1)+1)+10
546 # check of power of pi
547 # computing the fortcoef[mu,k,ell]
548 xfortcoef={}
549 for mu in range(0,der+1):
550 for k in range(0,L):
551 for ell in range(-2,3*k//2+1):
552 xfortcoef[mu,k,ell]=0
553 for mu in range(0,der+1):
554 for k in range(0,L):
555 ctx.prec = xwptcoef[k]
556 for ell in range(0,3*k//2+1):
557 xfortcoef[mu,k,ell]=xd[mu,k,ell]*Fp[3*k-2*ell]/pipowers[2*k-ell]
558 xfortcoef[mu,k,ell]=xfortcoef[mu,k,ell]/((2*ctx.j)**ell)
560 def trunc_a(t):
561 wp = ctx.prec
562 ctx.prec = wp + 2
563 aa = ctx.sqrt(t/(2*ctx.pi))
564 ctx.prec = wp
565 return aa
567 # computing the tcoef[k,ell]
568 xtcoef={}
569 for mu in range(0,der+1):
570 for k in range(0,L):
571 for ell in range(-2,3*k//2+1):
572 xtcoef[mu,k,ell]=0
573 ctx.prec = max(xwptcoef[0],ywptcoef[0])+3
574 aa= trunc_a(t)
575 la = -ctx.ln(aa)
577 for chi in range(0,der+1):
578 for k in range(0,L):
579 ctx.prec = xwptcoef[k]
580 for ell in range(0,3*k//2+1):
581 xtcoef[chi,k,ell] =0
582 for mu in range(0, chi+1):
583 tcoefter=ctx.binomial(chi,mu)*ctx.power(la,mu)*xfortcoef[chi-mu,k,ell]
584 xtcoef[chi,k,ell] += tcoefter
586 # COMPUTING THE COEFFICIENTS ytcoef[k,l]
587 # See II Section 3.9
588 #
589 # computing the needed wp
590 # check of power of pi
591 # computing the fortcoef[mu,k,ell]
592 yfortcoef={}
593 for mu in range(0,der+1):
594 for k in range(0,L):
595 for ell in range(-2,3*k//2+1):
596 yfortcoef[mu,k,ell]=0
597 for mu in range(0,der+1):
598 for k in range(0,L):
599 ctx.prec = ywptcoef[k]
600 for ell in range(0,3*k//2+1):
601 yfortcoef[mu,k,ell]=yd[mu,k,ell]*Fp[3*k-2*ell]/pipowers[2*k-ell]
602 yfortcoef[mu,k,ell]=yfortcoef[mu,k,ell]/((2*ctx.j)**ell)
603 # computing the tcoef[k,ell]
604 ytcoef={}
605 for chi in range(0,der+1):
606 for k in range(0,L):
607 for ell in range(-2,3*k//2+1):
608 ytcoef[chi,k,ell]=0
609 for chi in range(0,der+1):
610 for k in range(0,L):
611 ctx.prec = ywptcoef[k]
612 for ell in range(0,3*k//2+1):
613 ytcoef[chi,k,ell] =0
614 for mu in range(0, chi+1):
615 tcoefter=ctx.binomial(chi,mu)*ctx.power(la,mu)*yfortcoef[chi-mu,k,ell]
616 ytcoef[chi,k,ell] += tcoefter
618 # COMPUTING tv[k,ell]
619 # See II Section 3.8
620 #
621 # a has a good value
622 ctx.prec = max(xwptcoef[0], ywptcoef[0])+2
623 av = {}
624 av[0] = 1
625 av[1] = av[0]/a
627 ctx.prec = max(xwptcoef[0],ywptcoef[0])
628 for k in range(2,L):
629 av[k] = av[k-1] * av[1]
631 # Computing the quotients
632 xtv = {}
633 for chi in range(0,der+1):
634 for k in range(0,L):
635 ctx.prec = xwptcoef[k]
636 for ell in range(0,3*k//2+1):
637 xtv[chi,k,ell] = xtcoef[chi,k,ell]* av[k]
638 # Computing the quotients
639 ytv = {}
640 for chi in range(0,der+1):
641 for k in range(0,L):
642 ctx.prec = ywptcoef[k]
643 for ell in range(0,3*k//2+1):
644 ytv[chi,k,ell] = ytcoef[chi,k,ell]* av[k]
646 # COMPUTING THE TERMS xterm[k]
647 # See II Section 3.6
648 xterm = {}
649 for chi in range(0,der+1):
650 for n in range(0,L):
651 ctx.prec = xwpterm[n]
652 te = 0
653 for k in range(0, 3*n//2+1):
654 te += xtv[chi,n,k]
655 xterm[chi,n] = te
657 # COMPUTING THE TERMS yterm[k]
658 # See II Section 3.6
659 yterm = {}
660 for chi in range(0,der+1):
661 for n in range(0,L):
662 ctx.prec = ywpterm[n]
663 te = 0
664 for k in range(0, 3*n//2+1):
665 te += ytv[chi,n,k]
666 yterm[chi,n] = te
668 # COMPUTING rssum
669 # See II Section 3.5
670 xrssum={}
671 ctx.prec=15
672 xrsbound = math.sqrt(ctx.pi) * xc /(xb*a)
673 ctx.prec=15
674 xwprssum = ctx.mag(4.4*((L+3)**2)*xrsbound / xeps2)
675 xwprssum = max(xwprssum, ctx.mag(10*(L+1)))
676 ctx.prec = xwprssum
677 for chi in range(0,der+1):
678 xrssum[chi] = 0
679 for k in range(1,L+1):
680 xrssum[chi] += xterm[chi,L-k]
681 yrssum={}
682 ctx.prec=15
683 yrsbound = math.sqrt(ctx.pi) * yc /(yb*a)
684 ctx.prec=15
685 ywprssum = ctx.mag(4.4*((L+3)**2)*yrsbound / yeps2)
686 ywprssum = max(ywprssum, ctx.mag(10*(L+1)))
687 ctx.prec = ywprssum
688 for chi in range(0,der+1):
689 yrssum[chi] = 0
690 for k in range(1,L+1):
691 yrssum[chi] += yterm[chi,L-k]
693 # COMPUTING S3
694 # See II Section 3.19
695 ctx.prec = 15
696 A2 = 2**(max(ctx.mag(abs(xrssum[0])), ctx.mag(abs(yrssum[0]))))
697 eps8 = eps/(3*A2)
698 T = t *ctx.ln(t/(2*ctx.pi))
699 xwps3 = 5 + ctx.mag((1+(2/eps8)*ctx.power(a,-xsigma))*T)
700 ywps3 = 5 + ctx.mag((1+(2/eps8)*ctx.power(a,-ysigma))*T)
702 ctx.prec = max(xwps3, ywps3)
704 tpi = t/(2*ctx.pi)
705 arg = (t/2)*ctx.ln(tpi)-(t/2)-ctx.pi/8
706 U = ctx.expj(-arg)
707 a = trunc_a(t)
708 xasigma = ctx.power(a, -xsigma)
709 yasigma = ctx.power(a, -ysigma)
710 xS3 = ((-1)**(N-1)) * xasigma * U
711 yS3 = ((-1)**(N-1)) * yasigma * U
713 # COMPUTING S1 the zetasum
714 # See II Section 3.18
715 ctx.prec = 15
716 xwpsum = 4+ ctx.mag((N+ctx.power(N,1-xsigma))*ctx.ln(N) /eps1)
717 ywpsum = 4+ ctx.mag((N+ctx.power(N,1-ysigma))*ctx.ln(N) /eps1)
718 wpsum = max(xwpsum, ywpsum)
720 ctx.prec = wpsum +10
721 '''
722 # This can be improved
723 xS1={}
724 yS1={}
725 for chi in range(0,der+1):
726 xS1[chi] = 0
727 yS1[chi] = 0
728 for n in range(1,int(N)+1):
729 ln = ctx.ln(n)
730 xexpn = ctx.exp(-ln*(xsigma+ctx.j*t))
731 yexpn = ctx.conj(1/(n*xexpn))
732 for chi in range(0,der+1):
733 pown = ctx.power(-ln, chi)
734 xterm = pown*xexpn
735 yterm = pown*yexpn
736 xS1[chi] += xterm
737 yS1[chi] += yterm
738 '''
739 xS1, yS1 = ctx._zetasum(s, 1, int(N)-1, range(0,der+1), True)
741 # END OF COMPUTATION of xrz, yrz
742 # See II Section 3.1
743 ctx.prec = 15
744 xabsS1 = abs(xS1[der])
745 xabsS2 = abs(xrssum[der] * xS3)
746 xwpend = max(6, wpinitial+ctx.mag(6*(3*xabsS1+7*xabsS2) ) )
748 ctx.prec = xwpend
749 xrz={}
750 for chi in range(0,der+1):
751 xrz[chi] = xS1[chi]+xrssum[chi]*xS3
753 ctx.prec = 15
754 yabsS1 = abs(yS1[der])
755 yabsS2 = abs(yrssum[der] * yS3)
756 ywpend = max(6, wpinitial+ctx.mag(6*(3*yabsS1+7*yabsS2) ) )
758 ctx.prec = ywpend
759 yrz={}
760 for chi in range(0,der+1):
761 yrz[chi] = yS1[chi]+yrssum[chi]*yS3
762 yrz[chi] = ctx.conj(yrz[chi])
763 ctx.prec = wpinitial
764 return xrz, yrz
766def Rzeta_set(ctx, s, derivatives=[0]):
767 r"""
768 Computes several derivatives of the auxiliary function of Riemann `R(s)`.
770 **Definition**
772 The function is defined by
774 .. math ::
776 \begin{equation}
777 {\mathop{\mathcal R }\nolimits}(s)=
778 \int_{0\swarrow1}\frac{x^{-s} e^{\pi i x^2}}{e^{\pi i x}-
779 e^{-\pi i x}}\,dx
780 \end{equation}
782 To this function we apply the Riemann-Siegel expansion.
783 """
784 der = max(derivatives)
785 # First we take the value of ctx.prec
786 # During the computation we will change ctx.prec, and finally we will
787 # restaurate the initial value
788 wpinitial = ctx.prec
789 # Take the real and imaginary part of s
790 t = ctx._im(s)
791 sigma = ctx._re(s)
792 # Now compute several parameter that appear on the program
793 ctx.prec = 15
794 a = ctx.sqrt(t/(2*ctx.pi)) # Careful
795 asigma = ctx.power(a, sigma) # Careful
796 # We need a simple bound A1 < asigma (see II Section 3.1 and 3.3)
797 A1 = ctx.power(2, ctx.mag(asigma)-1)
798 # We compute various epsilon's (see II end of Section 3.1)
799 eps = ctx.power(2, -wpinitial)
800 eps1 = eps/6.
801 eps2 = eps * A1/3.
802 # COMPUTING SOME COEFFICIENTS THAT DEPENDS
803 # ON sigma
804 # constant b and c (see I Theorem 2 formula (26) )
805 # coefficients A and B1 (see I Section 6.1 equation (50))
806 # here we not need high precision
807 ctx.prec = 15
808 if sigma > 0:
809 b = 2.
810 c = math.pow(9,sigma)/4.44288
811 # 4.44288 =(math.sqrt(2)*math.pi)
812 A = math.pow(9,sigma)
813 B1 = 1
814 else:
815 b = 2.25158 # math.sqrt( (3-2* math.log(2))*math.pi )
816 c = math.pow(2,-sigma)/4.44288
817 A = math.pow(2,-sigma)
818 B1 = 1.10789 # = 2*sqrt(1-log(2))
819 # COMPUTING L THE NUMBER OF TERMS NEEDED IN THE RIEMANN-SIEGEL
820 # CORRECTION
821 # See II Section 3.2
822 ctx.prec = 15
823 L = 1
824 while 3*c*ctx.gamma(L*0.5) * ctx.power(b*a,-L) >= eps2:
825 L = L+1
826 L = max(2,L)
827 # The number L has to satify some conditions.
828 # If not RS can not compute Rzeta(s) with the prescribed precision
829 # (see II, Section 3.2 condition (20) ) and
830 # (II, Section 3.3 condition (22) ). Also we have added
831 # an additional technical condition in Section 3.17 Proposition 17
832 if ((3*L >= 2*a*a/25.) or (3*L+2+sigma<0) or (abs(sigma)> a/2.)):
833 #print 'Error Riemann-Siegel can not compute with such precision'
834 ctx.prec = wpinitial
835 raise NotImplementedError("Riemann-Siegel can not compute with such precision")
837 # INITIALIZATION (CONTINUATION)
838 #
839 # eps3 is the constant defined on (II, Section 3.5 equation (27) )
840 # each term of the RS correction must be computed with error <= eps3
841 eps3 = eps2/(4*L)
843 # eps4 is defined on (II Section 3.6 equation (30) )
844 # each component of the formula (II Section 3.6 equation (29) )
845 # must be computed with error <= eps4
846 eps4 = eps3/(3*L)
848 # COMPUTING M. NUMBER OF DERIVATIVES Fp[m] TO COMPUTE
849 M = aux_M_Fp(ctx, A, eps4, a, B1, L)
850 Fp = {}
851 for n in range(M, 3*L-2):
852 Fp[n] = 0
854 # But I have not seen an instance of M != 3*L-3
855 #
856 # DETERMINATION OF J THE NUMBER OF TERMS NEEDED
857 # IN THE TAYLOR SERIES OF F.
858 # See II Section 3.11 equation (49))
859 h1 = eps4/(632*A)
860 h2 = ctx.pi*ctx.pi*B1*a *ctx.sqrt(3)*math.e*math.e
861 h2 = h1 * ctx.power((h2/M**2),(M-1)/3) / M
862 h3 = min(h1,h2)
863 J=12
864 jvalue = (2*ctx.pi)**J / ctx.gamma(J+1)
865 while jvalue > h3:
866 J = J+1
867 jvalue = (2*ctx.pi)*jvalue/J
869 # COMPUTING eps5[m] for 1 <= m <= 21
870 # See II Section 10 equation (43)
871 eps5={}
872 foreps5 = math.pi*math.pi*B1*a
873 for m in range(0,22):
874 aux1 = math.pow(foreps5, m/3)/(316.*A)
875 aux2 = ctx.gamma(m+1)/ctx.gamma(m/3.0+0.5)
876 aux2 = math.sqrt(aux2)
877 eps5[m] = aux1*aux2*eps4
879 # COMPUTING wpfp
880 # See II Section 3.13 equation (59)
881 twenty = min(3*L-3, 21)+1
882 aux = 6812*J
883 wpfp = ctx.mag(44*J)
884 for m in range(0, twenty):
885 wpfp = max(wpfp, ctx.mag(aux*ctx.gamma(m+1)/eps5[m]))
886 # COMPUTING N AND p
887 # See II Section
888 ctx.prec = wpfp + ctx.mag(t) + 20
889 a = ctx.sqrt(t/(2*ctx.pi))
890 N = ctx.floor(a)
891 p = 1-2*(a-N)
893 # now we get a rounded version of p to the precision wpfp
894 # this possibly is not necessary
895 num = ctx.floor(p*(ctx.mpf(2)**wpfp))
896 difference = p * (ctx.mpf(2)**wpfp)-num
897 if difference < 0.5:
898 num = num
899 else:
900 num = num+1
901 p = ctx.convert(num * (ctx.mpf(2)**(-wpfp)))
903 # COMPUTING THE COEFFICIENTS c[n] = cc[n]
904 # We shall use the notation cc[n], since there is
905 # a constant that is called c
906 # See II Section 3.14
907 # We compute the coefficients and also save then in a
908 # cache. The bulk of the computation is passed to
909 # the function coef()
910 #
911 # eps6 is defined in II Section 3.13 equation (58)
912 eps6 = ctx.power(2*ctx.pi, J)/(ctx.gamma(J+1)*3*J)
914 # Now we compute the coefficients
915 cc={}
916 cont={}
917 cont, pipowers = coef(ctx, J, eps6)
918 cc = cont.copy() # we need a copy since we have
919 Fp={}
920 for n in range(M, 3*L-2):
921 Fp[n] = 0
922 ctx.prec = wpfp
923 for m in range(0,M+1):
924 sumP = 0
925 for k in range(2*J-m-1,-1,-1):
926 sumP = (sumP * p) + cc[k]
927 Fp[m] = sumP
928 # preparation of the new coefficients
929 for k in range(0, 2*J-m-1):
930 cc[k] = (k+1) * cc[k+1]
932 # COMPUTING THE NUMBERS d[n,k]
933 # See II Section 3.17
935 # First we compute the working precisions wpd[k]
936 # Se II equation (92)
937 wpd = {}
938 d1 = max(6, ctx.mag(40*L*L))
939 d2 = 13+ctx.mag((1+abs(sigma))*A)-ctx.mag(eps4)-1
940 const = ctx.ln(8/(ctx.pi*ctx.pi*a*a*B1*B1)) /2
941 for n in range(0,L):
942 d3 = ctx.mag(ctx.sqrt(ctx.gamma(n-0.5)))-ctx.floor(n*const)+d2
943 wpd[n] = max(d3,d1)
945 # procedure of II Section 3.17
946 ctx.prec = wpd[1]+10
947 psigma = 1-(2*sigma)
948 d = {}
949 d[0,0,-2]=0; d[0,0,-1]=0; d[0,0,0]=1; d[0,0,1]=0
950 d[0,-1,-2]=0; d[0,-1,-1]=0; d[0,-1,0]=1; d[0,-1,1]=0
951 for n in range(1,L):
952 ctx.prec = wpd[n]+10
953 for k in range(0,3*n//2+1):
954 m = 3*n-2*k
955 if (m!=0):
956 m1 = ctx.one/m
957 c1 = m1/4
958 c2 = (psigma*m1)/2
959 c3 = -(m+1)
960 d[0,n,k] = c3*d[0,n-1,k-2]+c1*d[0,n-1,k]+c2*d[0,n-1,k-1]
961 else:
962 d[0,n,k]=0
963 for r in range(0,k):
964 add = d[0,n,r]*(ctx.one*ctx.fac(2*k-2*r)/ctx.fac(k-r))
965 d[0,n,k] -= ((-1)**(k-r))*add
966 d[0,n,-2]=0; d[0,n,-1]=0; d[0,n,3*n//2+1]=0
968 for mu in range(-2,der+1):
969 for n in range(-2,L):
970 for k in range(-3,max(1,3*n//2+2)):
971 if ((mu<0)or (n<0) or(k<0)or (k>3*n//2)):
972 d[mu,n,k] = 0
974 for mu in range(1,der+1):
975 for n in range(0,L):
976 ctx.prec = wpd[n]+10
977 for k in range(0,3*n//2+1):
978 aux=(2*mu-2)*d[mu-2,n-2,k-3]+2*(sigma+n-2)*d[mu-1,n-2,k-3]
979 d[mu,n,k] = aux - d[mu-1,n-1,k-1]
981 # COMPUTING THE COEFFICIENTS t[k,l]
982 # See II Section 3.9
983 #
984 # computing the needed wp
985 wptcoef = {}
986 wpterm = {}
987 ctx.prec = 15
988 c1 = ctx.mag(40*(L+2))
989 c2 = ctx.mag(68*(L+2)*A)
990 c4 = ctx.mag(B1*a*math.sqrt(ctx.pi))-1
991 for k in range(0,L):
992 c3 = c2 - k*c4+ctx.mag(ctx.fac(k+0.5))/2.
993 wptcoef[k] = max(c1,c3-ctx.mag(eps4)+1)+1 +10
994 wpterm[k] = max(c1,ctx.mag(L+2)+c3-ctx.mag(eps3)+1)+1 +10
996 # check of power of pi
998 # computing the fortcoef[mu,k,ell]
999 fortcoef={}
1000 for mu in derivatives:
1001 for k in range(0,L):
1002 for ell in range(-2,3*k//2+1):
1003 fortcoef[mu,k,ell]=0
1005 for mu in derivatives:
1006 for k in range(0,L):
1007 ctx.prec = wptcoef[k]
1008 for ell in range(0,3*k//2+1):
1009 fortcoef[mu,k,ell]=d[mu,k,ell]*Fp[3*k-2*ell]/pipowers[2*k-ell]
1010 fortcoef[mu,k,ell]=fortcoef[mu,k,ell]/((2*ctx.j)**ell)
1012 def trunc_a(t):
1013 wp = ctx.prec
1014 ctx.prec = wp + 2
1015 aa = ctx.sqrt(t/(2*ctx.pi))
1016 ctx.prec = wp
1017 return aa
1019 # computing the tcoef[chi,k,ell]
1020 tcoef={}
1021 for chi in derivatives:
1022 for k in range(0,L):
1023 for ell in range(-2,3*k//2+1):
1024 tcoef[chi,k,ell]=0
1025 ctx.prec = wptcoef[0]+3
1026 aa = trunc_a(t)
1027 la = -ctx.ln(aa)
1029 for chi in derivatives:
1030 for k in range(0,L):
1031 ctx.prec = wptcoef[k]
1032 for ell in range(0,3*k//2+1):
1033 tcoef[chi,k,ell] = 0
1034 for mu in range(0, chi+1):
1035 tcoefter = ctx.binomial(chi,mu) * la**mu * \
1036 fortcoef[chi-mu,k,ell]
1037 tcoef[chi,k,ell] += tcoefter
1039 # COMPUTING tv[k,ell]
1040 # See II Section 3.8
1042 # Computing the powers av[k] = a**(-k)
1043 ctx.prec = wptcoef[0] + 2
1045 # a has a good value of a.
1046 # See II Section 3.6
1047 av = {}
1048 av[0] = 1
1049 av[1] = av[0]/a
1051 ctx.prec = wptcoef[0]
1052 for k in range(2,L):
1053 av[k] = av[k-1] * av[1]
1055 # Computing the quotients
1056 tv = {}
1057 for chi in derivatives:
1058 for k in range(0,L):
1059 ctx.prec = wptcoef[k]
1060 for ell in range(0,3*k//2+1):
1061 tv[chi,k,ell] = tcoef[chi,k,ell]* av[k]
1063 # COMPUTING THE TERMS term[k]
1064 # See II Section 3.6
1065 term = {}
1066 for chi in derivatives:
1067 for n in range(0,L):
1068 ctx.prec = wpterm[n]
1069 te = 0
1070 for k in range(0, 3*n//2+1):
1071 te += tv[chi,n,k]
1072 term[chi,n] = te
1074 # COMPUTING rssum
1075 # See II Section 3.5
1076 rssum={}
1077 ctx.prec=15
1078 rsbound = math.sqrt(ctx.pi) * c /(b*a)
1079 ctx.prec=15
1080 wprssum = ctx.mag(4.4*((L+3)**2)*rsbound / eps2)
1081 wprssum = max(wprssum, ctx.mag(10*(L+1)))
1082 ctx.prec = wprssum
1083 for chi in derivatives:
1084 rssum[chi] = 0
1085 for k in range(1,L+1):
1086 rssum[chi] += term[chi,L-k]
1088 # COMPUTING S3
1089 # See II Section 3.19
1090 ctx.prec = 15
1091 A2 = 2**(ctx.mag(rssum[0]))
1092 eps8 = eps/(3* A2)
1093 T = t * ctx.ln(t/(2*ctx.pi))
1094 wps3 = 5 + ctx.mag((1+(2/eps8)*ctx.power(a,-sigma))*T)
1096 ctx.prec = wps3
1097 tpi = t/(2*ctx.pi)
1098 arg = (t/2)*ctx.ln(tpi)-(t/2)-ctx.pi/8
1099 U = ctx.expj(-arg)
1100 a = trunc_a(t)
1101 asigma = ctx.power(a, -sigma)
1102 S3 = ((-1)**(N-1)) * asigma * U
1104 # COMPUTING S1 the zetasum
1105 # See II Section 3.18
1106 ctx.prec = 15
1107 wpsum = 4 + ctx.mag((N+ctx.power(N,1-sigma))*ctx.ln(N)/eps1)
1109 ctx.prec = wpsum + 10
1110 '''
1111 # This can be improved
1112 S1 = {}
1113 for chi in derivatives:
1114 S1[chi] = 0
1115 for n in range(1,int(N)+1):
1116 ln = ctx.ln(n)
1117 expn = ctx.exp(-ln*(sigma+ctx.j*t))
1118 for chi in derivatives:
1119 term = ctx.power(-ln, chi)*expn
1120 S1[chi] += term
1121 '''
1122 S1 = ctx._zetasum(s, 1, int(N)-1, derivatives)[0]
1124 # END OF COMPUTATION
1125 # See II Section 3.1
1126 ctx.prec = 15
1127 absS1 = abs(S1[der])
1128 absS2 = abs(rssum[der] * S3)
1129 wpend = max(6, wpinitial + ctx.mag(6*(3*absS1+7*absS2)))
1130 ctx.prec = wpend
1131 rz = {}
1132 for chi in derivatives:
1133 rz[chi] = S1[chi]+rssum[chi]*S3
1134 ctx.prec = wpinitial
1135 return rz
1138def z_half(ctx,t,der=0):
1139 r"""
1140 z_half(t,der=0) Computes Z^(der)(t)
1141 """
1142 s=ctx.mpf('0.5')+ctx.j*t
1143 wpinitial = ctx.prec
1144 ctx.prec = 15
1145 tt = t/(2*ctx.pi)
1146 wptheta = wpinitial +1 + ctx.mag(3*(tt**1.5)*ctx.ln(tt))
1147 wpz = wpinitial + 1 + ctx.mag(12*tt*ctx.ln(tt))
1148 ctx.prec = wptheta
1149 theta = ctx.siegeltheta(t)
1150 ctx.prec = wpz
1151 rz = Rzeta_set(ctx,s, range(der+1))
1152 if der > 0: ps1 = ctx._re(ctx.psi(0,s/2)/2 - ctx.ln(ctx.pi)/2)
1153 if der > 1: ps2 = ctx._re(ctx.j*ctx.psi(1,s/2)/4)
1154 if der > 2: ps3 = ctx._re(-ctx.psi(2,s/2)/8)
1155 if der > 3: ps4 = ctx._re(-ctx.j*ctx.psi(3,s/2)/16)
1156 exptheta = ctx.expj(theta)
1157 if der == 0:
1158 z = 2*exptheta*rz[0]
1159 if der == 1:
1160 zf = 2j*exptheta
1161 z = zf*(ps1*rz[0]+rz[1])
1162 if der == 2:
1163 zf = 2 * exptheta
1164 z = -zf*(2*rz[1]*ps1+rz[0]*ps1**2+rz[2]-ctx.j*rz[0]*ps2)
1165 if der == 3:
1166 zf = -2j*exptheta
1167 z = 3*rz[1]*ps1**2+rz[0]*ps1**3+3*ps1*rz[2]
1168 z = zf*(z-3j*rz[1]*ps2-3j*rz[0]*ps1*ps2+rz[3]-rz[0]*ps3)
1169 if der == 4:
1170 zf = 2*exptheta
1171 z = 4*rz[1]*ps1**3+rz[0]*ps1**4+6*ps1**2*rz[2]
1172 z = z-12j*rz[1]*ps1*ps2-6j*rz[0]*ps1**2*ps2-6j*rz[2]*ps2-3*rz[0]*ps2*ps2
1173 z = z + 4*ps1*rz[3]-4*rz[1]*ps3-4*rz[0]*ps1*ps3+rz[4]+ctx.j*rz[0]*ps4
1174 z = zf*z
1175 ctx.prec = wpinitial
1176 return ctx._re(z)
1178def zeta_half(ctx, s, k=0):
1179 """
1180 zeta_half(s,k=0) Computes zeta^(k)(s) when Re s = 0.5
1181 """
1182 wpinitial = ctx.prec
1183 sigma = ctx._re(s)
1184 t = ctx._im(s)
1185 #--- compute wptheta, wpR, wpbasic ---
1186 ctx.prec = 53
1187 # X see II Section 3.21 (109) and (110)
1188 if sigma > 0:
1189 X = ctx.sqrt(abs(s))
1190 else:
1191 X = (2*ctx.pi)**(sigma-1) * abs(1-s)**(0.5-sigma)
1192 # M1 see II Section 3.21 (111) and (112)
1193 if sigma > 0:
1194 M1 = 2*ctx.sqrt(t/(2*ctx.pi))
1195 else:
1196 M1 = 4 * t * X
1197 # T see II Section 3.21 (113)
1198 abst = abs(0.5-s)
1199 T = 2* abst*math.log(abst)
1200 # computing wpbasic, wptheta, wpR see II Section 3.21
1201 wpbasic = max(6,3+ctx.mag(t))
1202 wpbasic2 = 2+ctx.mag(2.12*M1+21.2*M1*X+1.3*M1*X*T)+wpinitial+1
1203 wpbasic = max(wpbasic, wpbasic2)
1204 wptheta = max(4, 3+ctx.mag(2.7*M1*X)+wpinitial+1)
1205 wpR = 3+ctx.mag(1.1+2*X)+wpinitial+1
1206 ctx.prec = wptheta
1207 theta = ctx.siegeltheta(t-ctx.j*(sigma-ctx.mpf('0.5')))
1208 if k > 0: ps1 = (ctx._re(ctx.psi(0,s/2)))/2 - ctx.ln(ctx.pi)/2
1209 if k > 1: ps2 = -(ctx._im(ctx.psi(1,s/2)))/4
1210 if k > 2: ps3 = -(ctx._re(ctx.psi(2,s/2)))/8
1211 if k > 3: ps4 = (ctx._im(ctx.psi(3,s/2)))/16
1212 ctx.prec = wpR
1213 xrz = Rzeta_set(ctx,s,range(k+1))
1214 yrz={}
1215 for chi in range(0,k+1):
1216 yrz[chi] = ctx.conj(xrz[chi])
1217 ctx.prec = wpbasic
1218 exptheta = ctx.expj(-2*theta)
1219 if k==0:
1220 zv = xrz[0]+exptheta*yrz[0]
1221 if k==1:
1222 zv1 = -yrz[1] - 2*yrz[0]*ps1
1223 zv = xrz[1] + exptheta*zv1
1224 if k==2:
1225 zv1 = 4*yrz[1]*ps1+4*yrz[0]*(ps1**2)+yrz[2]+2j*yrz[0]*ps2
1226 zv = xrz[2]+exptheta*zv1
1227 if k==3:
1228 zv1 = -12*yrz[1]*ps1**2-8*yrz[0]*ps1**3-6*yrz[2]*ps1-6j*yrz[1]*ps2
1229 zv1 = zv1 - 12j*yrz[0]*ps1*ps2-yrz[3]+2*yrz[0]*ps3
1230 zv = xrz[3]+exptheta*zv1
1231 if k == 4:
1232 zv1 = 32*yrz[1]*ps1**3 +16*yrz[0]*ps1**4+24*yrz[2]*ps1**2
1233 zv1 = zv1 +48j*yrz[1]*ps1*ps2+48j*yrz[0]*(ps1**2)*ps2
1234 zv1 = zv1+12j*yrz[2]*ps2-12*yrz[0]*ps2**2+8*yrz[3]*ps1-8*yrz[1]*ps3
1235 zv1 = zv1-16*yrz[0]*ps1*ps3+yrz[4]-2j*yrz[0]*ps4
1236 zv = xrz[4]+exptheta*zv1
1237 ctx.prec = wpinitial
1238 return zv
1240def zeta_offline(ctx, s, k=0):
1241 """
1242 Computes zeta^(k)(s) off the line
1243 """
1244 wpinitial = ctx.prec
1245 sigma = ctx._re(s)
1246 t = ctx._im(s)
1247 #--- compute wptheta, wpR, wpbasic ---
1248 ctx.prec = 53
1249 # X see II Section 3.21 (109) and (110)
1250 if sigma > 0:
1251 X = ctx.power(abs(s), 0.5)
1252 else:
1253 X = ctx.power(2*ctx.pi, sigma-1)*ctx.power(abs(1-s),0.5-sigma)
1254 # M1 see II Section 3.21 (111) and (112)
1255 if (sigma > 0):
1256 M1 = 2*ctx.sqrt(t/(2*ctx.pi))
1257 else:
1258 M1 = 4 * t * X
1259 # M2 see II Section 3.21 (111) and (112)
1260 if (1-sigma > 0):
1261 M2 = 2*ctx.sqrt(t/(2*ctx.pi))
1262 else:
1263 M2 = 4*t*ctx.power(2*ctx.pi, -sigma)*ctx.power(abs(s),sigma-0.5)
1264 # T see II Section 3.21 (113)
1265 abst = abs(0.5-s)
1266 T = 2* abst*math.log(abst)
1267 # computing wpbasic, wptheta, wpR see II Section 3.21
1268 wpbasic = max(6,3+ctx.mag(t))
1269 wpbasic2 = 2+ctx.mag(2.12*M1+21.2*M2*X+1.3*M2*X*T)+wpinitial+1
1270 wpbasic = max(wpbasic, wpbasic2)
1271 wptheta = max(4, 3+ctx.mag(2.7*M2*X)+wpinitial+1)
1272 wpR = 3+ctx.mag(1.1+2*X)+wpinitial+1
1273 ctx.prec = wptheta
1274 theta = ctx.siegeltheta(t-ctx.j*(sigma-ctx.mpf('0.5')))
1275 s1 = s
1276 s2 = ctx.conj(1-s1)
1277 ctx.prec = wpR
1278 xrz, yrz = Rzeta_simul(ctx, s, k)
1279 if k > 0: ps1 = (ctx.psi(0,s1/2)+ctx.psi(0,(1-s1)/2))/4 - ctx.ln(ctx.pi)/2
1280 if k > 1: ps2 = ctx.j*(ctx.psi(1,s1/2)-ctx.psi(1,(1-s1)/2))/8
1281 if k > 2: ps3 = -(ctx.psi(2,s1/2)+ctx.psi(2,(1-s1)/2))/16
1282 if k > 3: ps4 = -ctx.j*(ctx.psi(3,s1/2)-ctx.psi(3,(1-s1)/2))/32
1283 ctx.prec = wpbasic
1284 exptheta = ctx.expj(-2*theta)
1285 if k == 0:
1286 zv = xrz[0]+exptheta*yrz[0]
1287 if k == 1:
1288 zv1 = -yrz[1]-2*yrz[0]*ps1
1289 zv = xrz[1]+exptheta*zv1
1290 if k == 2:
1291 zv1 = 4*yrz[1]*ps1+4*yrz[0]*(ps1**2) +yrz[2]+2j*yrz[0]*ps2
1292 zv = xrz[2]+exptheta*zv1
1293 if k == 3:
1294 zv1 = -12*yrz[1]*ps1**2 -8*yrz[0]*ps1**3-6*yrz[2]*ps1-6j*yrz[1]*ps2
1295 zv1 = zv1 - 12j*yrz[0]*ps1*ps2-yrz[3]+2*yrz[0]*ps3
1296 zv = xrz[3]+exptheta*zv1
1297 if k == 4:
1298 zv1 = 32*yrz[1]*ps1**3 +16*yrz[0]*ps1**4+24*yrz[2]*ps1**2
1299 zv1 = zv1 +48j*yrz[1]*ps1*ps2+48j*yrz[0]*(ps1**2)*ps2
1300 zv1 = zv1+12j*yrz[2]*ps2-12*yrz[0]*ps2**2+8*yrz[3]*ps1-8*yrz[1]*ps3
1301 zv1 = zv1-16*yrz[0]*ps1*ps3+yrz[4]-2j*yrz[0]*ps4
1302 zv = xrz[4]+exptheta*zv1
1303 ctx.prec = wpinitial
1304 return zv
1306def z_offline(ctx, w, k=0):
1307 r"""
1308 Computes Z(w) and its derivatives off the line
1309 """
1310 s = ctx.mpf('0.5')+ctx.j*w
1311 s1 = s
1312 s2 = ctx.conj(1-s1)
1313 wpinitial = ctx.prec
1314 ctx.prec = 35
1315 # X see II Section 3.21 (109) and (110)
1316 # M1 see II Section 3.21 (111) and (112)
1317 if (ctx._re(s1) >= 0):
1318 M1 = 2*ctx.sqrt(ctx._im(s1)/(2 * ctx.pi))
1319 X = ctx.sqrt(abs(s1))
1320 else:
1321 X = (2*ctx.pi)**(ctx._re(s1)-1) * abs(1-s1)**(0.5-ctx._re(s1))
1322 M1 = 4 * ctx._im(s1)*X
1323 # M2 see II Section 3.21 (111) and (112)
1324 if (ctx._re(s2) >= 0):
1325 M2 = 2*ctx.sqrt(ctx._im(s2)/(2 * ctx.pi))
1326 else:
1327 M2 = 4 * ctx._im(s2)*(2*ctx.pi)**(ctx._re(s2)-1)*abs(1-s2)**(0.5-ctx._re(s2))
1328 # T see II Section 3.21 Prop. 27
1329 T = 2*abs(ctx.siegeltheta(w))
1330 # defining some precisions
1331 # see II Section 3.22 (115), (116), (117)
1332 aux1 = ctx.sqrt(X)
1333 aux2 = aux1*(M1+M2)
1334 aux3 = 3 +wpinitial
1335 wpbasic = max(6, 3+ctx.mag(T), ctx.mag(aux2*(26+2*T))+aux3)
1336 wptheta = max(4,ctx.mag(2.04*aux2)+aux3)
1337 wpR = ctx.mag(4*aux1)+aux3
1338 # now the computations
1339 ctx.prec = wptheta
1340 theta = ctx.siegeltheta(w)
1341 ctx.prec = wpR
1342 xrz, yrz = Rzeta_simul(ctx,s,k)
1343 pta = 0.25 + 0.5j*w
1344 ptb = 0.25 - 0.5j*w
1345 if k > 0: ps1 = 0.25*(ctx.psi(0,pta)+ctx.psi(0,ptb)) - ctx.ln(ctx.pi)/2
1346 if k > 1: ps2 = (1j/8)*(ctx.psi(1,pta)-ctx.psi(1,ptb))
1347 if k > 2: ps3 = (-1./16)*(ctx.psi(2,pta)+ctx.psi(2,ptb))
1348 if k > 3: ps4 = (-1j/32)*(ctx.psi(3,pta)-ctx.psi(3,ptb))
1349 ctx.prec = wpbasic
1350 exptheta = ctx.expj(theta)
1351 if k == 0:
1352 zv = exptheta*xrz[0]+yrz[0]/exptheta
1353 j = ctx.j
1354 if k == 1:
1355 zv = j*exptheta*(xrz[1]+xrz[0]*ps1)-j*(yrz[1]+yrz[0]*ps1)/exptheta
1356 if k == 2:
1357 zv = exptheta*(-2*xrz[1]*ps1-xrz[0]*ps1**2-xrz[2]+j*xrz[0]*ps2)
1358 zv =zv + (-2*yrz[1]*ps1-yrz[0]*ps1**2-yrz[2]-j*yrz[0]*ps2)/exptheta
1359 if k == 3:
1360 zv1 = -3*xrz[1]*ps1**2-xrz[0]*ps1**3-3*xrz[2]*ps1+j*3*xrz[1]*ps2
1361 zv1 = (zv1+ 3j*xrz[0]*ps1*ps2-xrz[3]+xrz[0]*ps3)*j*exptheta
1362 zv2 = 3*yrz[1]*ps1**2+yrz[0]*ps1**3+3*yrz[2]*ps1+j*3*yrz[1]*ps2
1363 zv2 = j*(zv2 + 3j*yrz[0]*ps1*ps2+ yrz[3]-yrz[0]*ps3)/exptheta
1364 zv = zv1+zv2
1365 if k == 4:
1366 zv1 = 4*xrz[1]*ps1**3+xrz[0]*ps1**4 + 6*xrz[2]*ps1**2
1367 zv1 = zv1-12j*xrz[1]*ps1*ps2-6j*xrz[0]*ps1**2*ps2-6j*xrz[2]*ps2
1368 zv1 = zv1-3*xrz[0]*ps2*ps2+4*xrz[3]*ps1-4*xrz[1]*ps3-4*xrz[0]*ps1*ps3
1369 zv1 = zv1+xrz[4]+j*xrz[0]*ps4
1370 zv2 = 4*yrz[1]*ps1**3+yrz[0]*ps1**4 + 6*yrz[2]*ps1**2
1371 zv2 = zv2+12j*yrz[1]*ps1*ps2+6j*yrz[0]*ps1**2*ps2+6j*yrz[2]*ps2
1372 zv2 = zv2-3*yrz[0]*ps2*ps2+4*yrz[3]*ps1-4*yrz[1]*ps3-4*yrz[0]*ps1*ps3
1373 zv2 = zv2+yrz[4]-j*yrz[0]*ps4
1374 zv = exptheta*zv1+zv2/exptheta
1375 ctx.prec = wpinitial
1376 return zv
1378@defun
1379def rs_zeta(ctx, s, derivative=0, **kwargs):
1380 if derivative > 4:
1381 raise NotImplementedError
1382 s = ctx.convert(s)
1383 re = ctx._re(s); im = ctx._im(s)
1384 if im < 0:
1385 z = ctx.conj(ctx.rs_zeta(ctx.conj(s), derivative))
1386 return z
1387 critical_line = (re == 0.5)
1388 if critical_line:
1389 return zeta_half(ctx, s, derivative)
1390 else:
1391 return zeta_offline(ctx, s, derivative)
1393@defun
1394def rs_z(ctx, w, derivative=0):
1395 w = ctx.convert(w)
1396 re = ctx._re(w); im = ctx._im(w)
1397 if re < 0:
1398 return rs_z(ctx, -w, derivative)
1399 critical_line = (im == 0)
1400 if critical_line :
1401 return z_half(ctx, w, derivative)
1402 else:
1403 return z_offline(ctx, w, derivative)