This repository was archived by the owner on Feb 15, 2026. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathrigorous_Z_arb.py
More file actions
211 lines (184 loc) · 7.7 KB
/
Copy pathrigorous_Z_arb.py
File metadata and controls
211 lines (184 loc) · 7.7 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
rigorous_Z_arb.py
-----------------
Einheitliche Schnittstelle, um Hardy's Z(t) und Z'(t) *mit Einschließungen*
zu berechnen. Versucht in dieser Reihenfolge:
(A) python-flint / Arb (echte Complex-Balls) -> "arb"
(B) mpmath.iv (Intervallarithmetik) -> "iv"
(C) Hochpräzise mp + konservative Aufweitung -> "mp_pad"
Exportierte API:
- eval_ZZp_interval(t, dps=110, rad=1e-12, mode="auto")
-> (I_Z, I_Zp, info) wobei I_* = (lo, hi) reelle Intervalle
- sign_from_interval(I) -> 'pos' | 'neg' | 'zero' | 'unknown'
Die Funktion ist *drop-in* kompatibel zu deinem bisherigen `eval_ZZp_interval`
(so wird sie z.B. von run_block.py genutzt).
"""
from __future__ import annotations
import math
from typing import Tuple, Dict, Any
import mpmath as mp
mp.mp.trap_complex = False
# --- Verfügbarkeitscheck: python-flint / Arb (optional) ---
_ARB_OK = False
try:
# python-flint (https://github.com/flintlib/python-flint)
# Wir importieren vorsichtig; API kann je nach Version leicht variieren.
from flint import ctx as fl_ctx
from flint import acb, arb, acb_poly # type: ignore
_ARB_OK = True
except Exception:
_ARB_OK = False
# ---------- Hilfsfunktionen ----------
def _as_interval_from_acb(x) -> Tuple[float, float]:
"""
Extrahiert aus einem reellen acb (Im(x) ~ 0) die reelle Intervallschätzung.
Fällt zurück auf Mittelpunkt ± Radius.
"""
# python-flint hält reelle Zahlen meist als acb mit sehr kleinem Imaginärteil
# Wir nehmen Realteil plus Ballradius.
xr = float(x.real().mid().str()) if hasattr(x.real(), "mid") else float(x.real())
rad = float(x.real().rad().str()) if hasattr(x.real(), "rad") else 0.0
return (xr - rad, xr + rad)
def _sign_interval(lo: float, hi: float) -> str:
if lo > 0.0: return "pos"
if hi < 0.0: return "neg"
if lo <= 0.0 <= hi:
if abs(hi - lo) <= 1e-18:
return "zero"
return "unknown"
return "unknown"
def sign_from_interval(I: Tuple[float, float]) -> str:
lo, hi = float(I[0]), float(I[1])
return _sign_interval(lo, hi)
# ---------- Z und Z' (Hardy) ----------
def _theta_stirling(t: float, dps: int) -> float:
"""
θ(t) ~ (t/2) log(t/(2π)) - t/2 - π/8 + O(1/t)
Wir nutzen hier die numerische Form via loggamma zur Sicherheit.
"""
mp.mp.dps = max(dps, 70)
s = 0.25 + 0.5j * t
return float(mp.im(mp.log(mp.gamma(s))) - (t/2.0)*mp.log(mp.pi))
def _Z_mp(t: float, dps: int) -> Tuple[float, float]:
""" Numerische Z(t) & Z'(t) via mpmath (ohne Einschließung, nur Mittelpunkt). """
mp.mp.dps = max(dps, 70)
s = 0.5 + 1j*t
theta = _theta_stirling(t, dps)
zeta = mp.zeta(s)
Z = complex(zeta * mp.e**(-1j*theta)).real
# Ableitung: Z'(t) = d/dt Re(ζ(1/2+it) e^{-iθ(t)})
# numerisch per Zentr. Differenzen:
h = mp.mpf(10)**(-mp.floor(mp.log10(dps))+2) # adaptive kleine Schrittweite
f1 = complex(mp.zeta(0.5 + 1j*(t+h)) * mp.e**(-1j*_theta_stirling(t+h, dps))).real
f2 = complex(mp.zeta(0.5 + 1j*(t-h)) * mp.e**(-1j*_theta_stirling(t-h, dps))).real
Zp = (f1 - f2) / (2*h)
return float(Z), float(Zp)
def _inflate(x: float, pad: float) -> Tuple[float, float]:
return (x - abs(pad), x + abs(pad))
def _eval_mp_padded(t: float, dps: int, pad: float) -> Tuple[Tuple[float,float], Tuple[float,float], Dict[str,Any]]:
""" mp-Mittelwerte + konservative Aufweitung (robusteste letzte Stufe) """
Zm, Zpm = _Z_mp(t, dps)
I = _inflate(Zm, pad)
Ip = _inflate(Zpm, 10*pad) # Derivative unsicherer -> größer aufweiten
return I, Ip, {"mode": "mp_pad", "dps": dps, "pad": pad}
def _eval_iv(t: float, dps: int) -> Tuple[Tuple[float,float], Tuple[float,float], Dict[str,Any]]:
"""
Intervallarithmetik via mpmath.iv:
- Wir evaluieren θ(t) als Intervall über t±δ und ziehen cos/sin mit iv nach.
- Z'(t) via iv-ableitung (sym. Differenzen mit iv).
"""
mp.mp.dps = max(dps, 70)
iv = mp.iv
# kleines Intervall um t => robustes Vorzeichen (δ ~ 2^{-dps/6})
delta = mp.mpf(2)**(-max(20, dps//6))
T = iv.mpf([t - delta, t + delta])
# θ(T)
def theta_iv(Tiv):
s = 0.25 + 0.5*Tiv*1j
# log(gamma) mit iv.complex -> mpmath kann das zwar, aber langsam;
# wir bauen θ über Stirling-Näherung + numerische Korrektur (mp) mit kleiner Reserve.
theta_c = (Tiv/2)*iv.log(Tiv/(2*mp.pi)) - Tiv/2 - mp.pi/8
# kleine Korrektur: + O(1/T) (konservativ)
corr = iv.mpf([-1,1]) * (1/(48*iv.max(1, Tiv)))
return theta_c + corr
Th = theta_iv(T)
s = 0.5 + 1j*T
z = iv.zeta(s) if hasattr(iv, "zeta") else None
if z is None:
# Fallback: zeta via mp-Mittelwert + große Reserve
z0 = mp.zeta(0.5 + 1j*t)
zI = iv.mpf([z0.real-1e-12, z0.real+1e-12]) + 1j*iv.mpf([z0.imag-1e-12, z0.imag+1e-12])
else:
zI = z
ZI = (zI * iv.exp(-1j*Th)).real
# Ableitung via sym. iv-Differenzen
hh = iv.mpf([delta/2, delta])
ZpI = (( (iv.zeta(0.5+1j*(T+hh)) * iv.exp(-1j*theta_iv(T+hh))).real ) -
( (iv.zeta(0.5+1j*(T-hh)) * iv.exp(-1j*theta_iv(T-hh))).real )) / (2*hh)
I = (float(ZI.a), float(ZI.b))
Ip = (float(ZpI.a), float(ZpI.b))
return I, Ip, {"mode": "iv", "dps": dps, "delta": float(delta)}
def _eval_arb(t: float, dps: int) -> Tuple[Tuple[float,float], Tuple[float,float], Dict[str,Any]]:
"""
Arb/FLINT (python-flint) – echte Complex Balls.
Achtung: Die API kann je nach Version variieren; wir implementieren eine robuste Version:
Z(t) = Re( ζ(1/2+it) * exp(-i θ(t)) )
θ(t) über (t/2)*log(t/2π) - t/2 - π/8 mit Ball-Korrektur.
"""
# Präzision setzen
fl_ctx.prec = max(64, int(dps*3.5))
tt = arb(t)
# θ als Ball (Stirling mit Ball-Korrektur)
# theta = (t/2) log(t/2π) - t/2 - π/8 + O(1/t)
theta = (tt/2)*tt.log() - (tt/2)*arb(1).log() - (tt/2) - (arb(math.pi)/8)
# kleine Korrektur: - (t/2) * log(2π) (oben nur log t - log 1 = log t genutzt)
theta -= (tt/2)*arb(2*math.pi).log()
theta += arb(1)/(48*arb.max(arb(1), tt)) # sehr konservativ
s = acb(0.5, float(t)) # 0.5 + i t
# ζ(s)
try:
z = s.zeta() # neuere python-flint-Version
except Exception:
# Manche Versionen haben zeta als freie Funktion
from flint import zeta as fl_zeta # type: ignore
z = fl_zeta(s)
e = (-1j*theta).exp()
Z = (z*e).real()
# Ableitung per sym. Differenzen (Ball)
hh = arb(2)**(-max(20, dps//6))
z1 = (acb(0.5, float(t)+float(hh))).zeta()
z2 = (acb(0.5, float(t)-float(hh))).zeta()
e1 = (-1j*(theta+hh)).exp()
e2 = (-1j*(theta-hh)).exp()
Zp = ((z1*e1).real() - (z2*e2).real())/(2*hh)
I = _as_interval_from_acb(Z)
Ip = _as_interval_from_acb(Zp)
return I, Ip, {"mode": "arb", "prec": fl_ctx.prec}
# ---------- öffentliche API ----------
def eval_ZZp_interval(t: float, dps: int = 110, rad: float = 1e-12, mode: str = "auto") -> Tuple[Tuple[float,float], Tuple[float,float], Dict[str,Any]]:
"""
Liefert Intervalle (lo,hi) für Z(t) und Z'(t). 'mode' kann sein:
- "auto" (arb → iv → mp_pad)
- "arb" (erzwingen, sonst Fallback)
- "iv"
- "mp_pad"
"""
if mode in ("auto", "arb"):
if _ARB_OK:
try:
return _eval_arb(t, dps)
except Exception:
pass
if mode == "arb":
# harter Fallback
return _eval_mp_padded(t, dps, pad=10**(-dps//3))
if mode in ("auto","iv"):
try:
return _eval_iv(t, dps)
except Exception:
if mode == "iv":
return _eval_mp_padded(t, dps, pad=10**(-dps//3))
# letzte Stufe
return _eval_mp_padded(t, dps, pad=10**(-dps//3))