Created
July 28, 2020 18:07
-
-
Save fredrik-johansson/7c2711887811ef9f2d7038b8451a4e63 to your computer and use it in GitHub Desktop.
Extended Rademacher series
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from flint import fmpq, fmpz | |
from cmath import * | |
from mpmath import plot | |
def dedekind_sum(r, k, _cache={}): | |
key = (r << 24) | k | |
if key in _cache: | |
return _cache[key] | |
if fmpz(r).gcd(k) != 1: | |
v = None | |
else: | |
u = fmpq.dedekind_sum(r, k) | |
p = u.p | |
q = u.q | |
v = float(int(p)) / float(int(q)) | |
_cache[key] = v | |
return v | |
def rademacher_A(n, k, exponential=False): | |
from cmath import exp, pi, cos | |
s = 0.0 | |
for r in range(k): | |
v = dedekind_sum(r, k) | |
if v is not None: | |
if exponential: | |
s += exp(pi * 1j * (v - 2.0*n*r/k)) | |
else: | |
s += cos(pi * (v - 2.0*n*r/k)) | |
return s | |
def bessel32(x): | |
from cmath import sinh, cosh, pi, sqrt | |
s = sinh(x) | |
c = cosh(x) | |
return (2*c-2*s/x)/(sqrt(2*pi*x)) | |
def rademacher_p(n, tol=1e-5, consecutive=20, exponential=False, verbose=False): | |
from cmath import sinh, cosh, pi, sqrt | |
n = complex(n) | |
k = 1 | |
s = 0.0 | |
N = 10**9 | |
c = sqrt(2*(n-1.0/24)/3) | |
print(n) | |
while k < N: | |
a = rademacher_A(n, k, exponential=exponential) | |
x = (pi / k) * c | |
t = a * bessel32(x) / k | |
if verbose: | |
print(k, abs(t)) | |
if abs(t) < tol: | |
N = min(N, k + consecutive) | |
else: | |
N = 10**9 | |
s += t | |
k += 1 | |
return 2.0*s*pi/complex(24*n-1)**0.75 | |
plot(lambda x: rademacher_p(x), [-5,5]) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment