Skip to content

Instantly share code, notes, and snippets.

@fredrik-johansson
Created July 28, 2020 18:07
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save fredrik-johansson/7c2711887811ef9f2d7038b8451a4e63 to your computer and use it in GitHub Desktop.
Save fredrik-johansson/7c2711887811ef9f2d7038b8451a4e63 to your computer and use it in GitHub Desktop.
Extended Rademacher series
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