<jupyter_start><jupyter_text>Computational Experiments on Maass Forms<jupyter_code>T = matrix([[1, 1], [0, 1]])
S = matrix([[0, -1], [1, 0]])
Id = matrix([[1, 0], [0, 1]])

def act(gamma, z):
    """
    Computes the action of matrix gamma on point z
    """
    a, b = gamma[0]
    c, d = gamma[1]
    return (a*z + b)/ (c*z + d)

def is_in_fund_domain(z):
    """
    True if z is in fundamental domain, False otherwise.
    """
    x = z.real()
    if x < -1/2 or x > 1/2:
        return False
    if z.norm() < 1:
        return False
    return True


def pullback(z):
    """
    Returns matrix gamma, point z^*, where gamma z = z^* and z^* is in the fundamental domain
    """
    gamma = Id
    while not is_in_fund_domain(z):
        x, y = z.real(), z.imag()
        xshift = -floor(x + .5)
        gamma = T**xshift * gamma
        z = x + xshift + i*y
        if z.norm() < 1:
            z = act(S, z)
            gamma = S * gamma
    return gamma, z
pullback(CC(3.75 + .25*I))
pullback(CDF(3.75 + .25*I))
bessel_K(CC(1.0), CC(3.0))
bessel_K(CDF(1.0), CDF(3.0))
def kappa(r, u):
    return (exp(CC(pi) * CC(r) / 2) * bessel_K(CC(i*r), CC(u)))
#def kappa(r, u):
#    return bessel_K(CDF(i*r), CDF(u))
def make_zj_list(Y, Q=50):
    """
    Returns list of z_j = x_j + iY
    """
    ret = []
    for j in range(1-Q, Q+1):
        ret.append((1/(2*Q)*(j - 0.5) + i*Y))
    return ret

def make_zjstar_list(zjlist):
    """
    Returns list of z_j^* in same order as zjlist
    """
    return list(map(lambda zj: pullback(zj)[1], zjlist))
test_zj_list = make_zj_list(0.5, Q=4)
for elem in test_zj_list:
    print(elem)
test_zjstar_list = make_zjstar_list(test_zj_list)
for elem in test_zjstar_list:
    print(elem)
def perform_sanity_check_zjstar(zjlist, zjstarlist):
    for zj, zjstar in zip(zjlist, zjstarlist):
        mat, _ = pullback(zj)
        assert (act(mat, CC(zj)) - zjstar).norm() < 0.00000001
    print("Passed.")
perform_sanity_check_zjstar(test_zj_list, test_zjstar_list)
def V(n, l, zjlist, r, Q, zjstarlist=None):
    """
    Return the constant V(n, l) --- which also depends on zj, r, and Q --- as
    defined above.
    """
    if zjstarlist is None:
        zjstarlist = make_zjstar_list(zjlist)
    ret = CDF(0)
    assert len(zjstarlist) == 2*Q
    for zj, zjs in zip(zjlist, zjstarlist):
        xj = zj.real()
        xjs = zjs.real()
        yjs = zjs.imag()
        ret += (sqrt(yjs)) * kappa(r, 2 * pi * abs(l) * yjs) * (exp(2*pi*i*(l * xjs - n * xj)))
    ret = ret / (2*Q)
    return CDF(ret)<jupyter_output><empty_output><jupyter_text>Large Speedup TODOEssentially all computation time is spent currently computing $V(n, l)$ in order to make the matrix. There are $Q^2$ matrix entries, each formed from a sum of $Q$ values. And right now, I'm computing all $Q^3$ K-Bessel functions. But this is excessive. In fact, there are at most $Q^2$ K-Bessel functions that need to be computed:$$ K_{ir}(2 \pi \lvert \ell \rvert y_j^*),$$where $1 \leq \lvert \ell \rvert \leq Q$ and $j \in [-Q, Q]$.Rewriting the matrix-generation code to take advantage of this with caching would make it a $Q^2$ operation (assuming that the Bessel functions are where all of the computation time is being spent).Section 3.4 of [this paper](http://www2.math.uu.se/~astrombe/papers/kbessel.pdf) gives a further improvement for repeated honing in, as only the values of $r$ are changing. Sections 3.2 and 3.3 might give improvements for an individual run --- I don't quite understand the implied algorithms.The matrix-generation code should also be vectorized as much as possible.<jupyter_code>V(3, 1, test_zj_list, 9.5**.5, 4)
import time
now = time.time()
print(V(3, 1, test_zj_list, 9.5**.5, 4))
print("Took {}".format(time.time() - now))
def make_matrix(size, Y, r, withtime=True):
    if withtime:
        now = time.time()
    Q = size
    zjlist = make_zj_list(Y, Q)
    zjstarlist = make_zjstar_list(zjlist)
    if withtime:
        print("lists made in {}".format(time.time() - now))
    matrix = []
    B = []
    for n in range(-Q, Q+1):
        if withtime:
            now = time.time()
        if n == 0:
            continue
        row = []
        for l in range(-Q, Q+1):
            if l == 0:
                continue
            entry = CDF(V(n, l, zjlist, CDF(r), CDF(Q), zjstarlist=zjstarlist))
            if n == l:
                #entry = entry - sqrt(Y) * bessel_K(i*r, 2 * pi * abs(n) * Y)
                entry = entry - CDF(sqrt(Y)) * kappa(r, 2 * pi * abs(n) * Y)
            if l != 1:
                row.append(entry.n())
            else:
                B.append([-entry.n()])
                #if n == 1:
                #    print(V(n, l, zjlist, r, Q, zjstarlist=zjstarlist))
                #    print(sqrt(Y) * kappa(r, 2 * pi * abs(n) * Y))
        matrix.append(row)
        if withtime:
            print("row {} computed in {} seconds".format(n, time.time() - now))
    return matrix, B
Q = 8
Y = 0.28
r = 9.53369526135
now = time.time()
#testmat, testB = make_matrix(20, 0.5, 3)
testmat, testB = make_matrix(Q, Y, r)
print("Done in {}".format(time.time() - now))<jupyter_output>lists made in 0.03257441520690918
row -8 computed in 0.6628572940826416 seconds
row -7 computed in 0.6433212757110596 seconds
row -6 computed in 0.6426317691802979 seconds
row -5 computed in 0.6682097911834717 seconds
row -4 computed in 0.6700131893157959 seconds
row -3 computed in 0.6402966976165771 seconds
row -2 computed in 0.663346529006958 seconds
row -1 computed in 0.6863448619842529 seconds
row 1 computed in 0.7095780372619629 seconds
row 2 computed in 0.9151451587677002 seconds
row 3 computed in 0.798314094543457 seconds
row 4 computed in 0.9088175296783447 seconds
row 5 computed in 0.8387832641601562 seconds
row 6 computed in 0.8383357524871826 seconds
row 7 computed in 0.8320372104644775 seconds
row 8 computed in 0.8334674835205078 seconds
Done in 11.985921859741211<jupyter_text>Contents of tiny cython file```pythonimport sage.libs.mpmath.all as mpmathmpmath_ctx = mpmath.fpdef cykappa(r, u):    return mpmath_ctx.besselk(1j*r, u)def cykappa_with_mult(r, u):    pi = mpmath_ctx.pi    cdef double exp_Pih_R=mpmath_ctx.exp(pi * r / 2.)    return mpmath_ctx.besselk(1j*r, u) * exp_Pih_Rdef cyV(n, l, r, Q, zjlist, zjstarlist):    pi = mpmath_ctx.pi    cdef double exp_Pih_R=mpmath_ctx.exp(pi * r / 2.)    cdef double complex ret = 0    for zj, zjs in zip(zjlist, zjstarlist):        xj = zj.real()        xjs = zjs.real()        yjs = zjs.imag()        ret += mpmath_ctx.sqrt(yjs) * exp_Pih_R * cykappa(r, 2*pi*abs(l)*yjs) \               * mpmath_ctx.exp(2*pi*1j*(l * xjs - n * xj))    ret = ret / (2*Q)    return ret```<jupyter_code>load("maass-cython.spyx")
cykappa_with_mult(1, 3+i)
kappa(1, 3+i)
test_zj_starlist = make_zjstar_list(test_zj_list)
cyV(3, 1, 9.5**.5, 4, test_zj_list, test_zj_starlist)
V(3, 1, test_zj_list, 9.5**.5, 4)
%time V(3, 1, test_zj_list, 9.5**.5, 4)
%time cyV(3, 1, 9.5**.5, 4, test_zj_list, test_zj_starlist)
def make_matrix_with_some_cy(size, Y, r, withtime=True):
    if withtime:
        now = time.time()
    Q = size
    zjlist = make_zj_list(Y, Q)
    zjstarlist = make_zjstar_list(zjlist)
    if withtime:
        print("lists made in {}".format(time.time() - now))
    matrix = []
    B = []
    for n in range(-Q, Q+1):
        if withtime:
            now = time.time()
        if n == 0:
            continue
        row = []
        for l in range(-Q, Q+1):
            if l == 0:
                continue
            entry = cyV(n, l, CDF(r), CDF(Q), zjlist, zjstarlist=zjstarlist)
            if n == l:
                #entry = entry - sqrt(Y) * bessel_K(i*r, 2 * pi * abs(n) * Y)
                entry = entry - CDF(sqrt(Y)) * cykappa_with_mult(r, 2 * pi * abs(n) * Y)
            if l != 1:
                row.append(entry)
            else:
                B.append([-entry])
                #if n == 1:
                #    print(V(n, l, zjlist, r, Q, zjstarlist=zjstarlist))
                #    print(sqrt(Y) * kappa(r, 2 * pi * abs(n) * Y))
        matrix.append(row)
        if withtime:
            print("row {} computed in {} seconds".format(n, time.time() - now))
    return matrix, B
Q = 8
Y = 0.28
r = 9.53369526135
#testmat, testB = make_matrix(20, 0.5, 3)
now = time.time()
testmat, testB = make_matrix_with_some_cy(Q, Y, r)
print("Done in {}".format(time.time() - now))
def solve_for_coeffs(Q, Y, r):
    mat, B = make_matrix_with_some_cy(Q, Y, r, withtime=False)
    preinversiontime = time.time()
    mm = matrix(CC, mat[1:])
    Bm = matrix(CC, B[1:])
    res = mm.inverse() * Bm
    print("Inversion done in {}".format(time.time() - preinversiontime))
    indices = [i for i in range(-Q, Q) if i != 0 and i != 1]
    resdict = dict()
    for index, re in zip(indices, res):
        resdict[index] = re[0].real()
    return resdict
%time coeffdict = solve_for_coeffs(Q, Y, r)
coeffdict
def error_func(r, Y1=0.28, Y2=0.3, Q=8, power=2):
    coeffdict1 = solve_for_coeffs(Q, Y1, r)
    coeffdict2 = solve_for_coeffs(Q, Y2, r)
    return ((coeffdict1[2] - coeffdict2[2])**power
            + (coeffdict1[3] - coeffdict2[3]) **power
            + (coeffdict1[5] - coeffdict2[5])**power)
%time error_func(r)
def secant_method_single_iteration(r1, r2, Q=8, power=2):
    e1 = error_func(r1, Q=Q, power=power)
    e2 = error_func(r2, Q=Q, power=power)
    return r1 - e1*(r1 - r2)/(e1 - e2)
r1 = 9.5
r2 = 9.55
rlist = [r1, r2]
NOW = time.time()
now = time.time()
for idx in range(10):
    rnew = secant_method_single_iteration(rlist[-2], rlist[-1], power=2)
    rlist.append(rnew)
    print("computed r: {}      time: {}".format(rnew, time.time() - now))
    now = time.time()
print("Done in {}".format(time.time() - NOW))
r1 = 9.5
r2 = 9.55
rlist = [r1, r2]
NOW = time.time()
now = time.time()
for idx in range(10):
    rnew = secant_method_single_iteration(rlist[-2], rlist[-1], power=1)
    rlist.append(rnew)
    print("computed r: {}      time: {}".format(rnew, time.time() - now))
    now = time.time()
print("Done in {}".format(time.time() - NOW))
r1 = 9.5
r2 = 9.55
rlist = [r1, r2]
NOW = time.time()
now = time.time()
for idx in range(10):
    rnew = secant_method_single_iteration(rlist[-2], rlist[-1], Q=6)
    rlist.append(rnew)
    print("computed r: {}      time: {}".format(rnew, time.time() - now))
    now = time.time()
print("Done in {}".format(time.time() - NOW))
error_func(r, Q=6)<jupyter_output>Inversion done in 0.00270843505859375
Inversion done in 0.0026862621307373047<jupyter_text>I don't quite understand the required bounds on $Q$ at the moment. But surprisingly small values of $Q$ work surprisingly well.<jupyter_code>r1 = 12.0
r2 = 12.5
rlist = [r1, r2]
NOW = time.time()
now = time.time()
for idx in range(20):
    rnew = secant_method_single_iteration(rlist[-2], rlist[-1], Q=8)
    rlist.append(rnew)
    print("computed r: {}      time: {}".format(rnew, time.time() - now))
    now = time.time()
print("Done in {}".format(time.time() - NOW))
r1 = 12.15
r2 = 12.2
rlist = [r1, r2]
NOW = time.time()
now = time.time()
for idx in range(10):
    rnew = secant_method_single_iteration(rlist[-2], rlist[-1], Q=6)
    rlist.append(rnew)
    print("computed r: {}      time: {}".format(rnew, time.time() - now))
    now = time.time()
print("Done in {}".format(time.time() - NOW))<jupyter_output>Inversion done in 0.002987384796142578
Inversion done in 0.003139019012451172
Inversion done in 0.002670764923095703
Inversion done in 0.003651857376098633
computed r: 12.2686437411483      time: 2.285745143890381
Inversion done in 0.002882719039916992
Inversion done in 0.003132343292236328
Inversion done in 0.0028274059295654297
Inversion done in 0.0029087066650390625
computed r: 12.1999812501054      time: 2.3465421199798584
Inversion done in 0.002807140350341797
Inversion done in 0.0027446746826171875
Inversion done in 0.0028433799743652344
Inversion done in 0.002790212631225586
computed r: 12.1999625422140      time: 2.2490718364715576
Inversion done in 0.0027322769165039062
Inversion done in 0.0028107166290283203
Inversion done in 0.002768278121948242
Inversion done in 0.003839254379272461
computed r: 12.1925157918053      time: 2.279813051223755
Inversion done in 0.0028047561645507812
Inversion done in 0.0028917789459228516
Inversion done in 0.0028030872344970703
Inversion done i[...]<jupyter_text>I am rapidly running into precision errors. Many of these are caused, I think, by the larger-than-necessary K-bessel functions. I'm using $K_{ir}(u)$ instead of Andy's $\kappa_{ir}(u)$, simply because the libraries I have can compute these quickly. I should try to find a better K-bessel computation.<jupyter_code>r1 = 12.15
r2 = 12.2
rlist = [r1, r2]
NOW = time.time()
now = time.time()
for idx in range(10):
    rnew = secant_method_single_iteration(rlist[-2], rlist[-1], Q=10)
    rlist.append(rnew)
    print("computed r: {}      time: {}".format(rnew, time.time() - now))
    now = time.time()
print("Done in {}".format(time.time() - NOW))<jupyter_output>Inversion done in 0.010133981704711914
Inversion done in 0.009264945983886719
Inversion done in 0.009323835372924805
Inversion done in 0.011097908020019531
computed r: 12.1495446393491      time: 8.355088472366333
Inversion done in 0.00950479507446289
Inversion done in 0.010741949081420898
Inversion done in 0.009484529495239258
Inversion done in 0.00935220718383789
computed r: 12.1490742084784      time: 8.347170114517212
Inversion done in 0.009395122528076172
Inversion done in 0.009287595748901367
Inversion done in 0.009331941604614258
Inversion done in 0.011080503463745117
computed r: 12.1694270709510      time: 8.331717491149902
Inversion done in 0.009650468826293945
Inversion done in 0.009494304656982422
Inversion done in 0.009322166442871094
Inversion done in 0.009551525115966797
computed r: 12.1705659328399      time: 8.360275268554688
Inversion done in 0.010320663452148438
Inversion done in 0.011347293853759766
Inversion done in 0.010271310806274414
Inversion done in 0.009859323[...]<jupyter_text>There are apparently precision problems; see the following computation. $Q = 20$ is much higher, but the actual convergence is not any faster.<jupyter_code>r1 = 12.15
r2 = 12.2
rlist = [r1, r2]
NOW = time.time()
now = time.time()
for idx in range(10):
    rnew = secant_method_single_iteration(rlist[-2], rlist[-1], Q=20)
    rlist.append(rnew)
    print("computed r: {}      time: {}".format(rnew, time.time() - now))
    now = time.time()
print("Done in {}".format(time.time() - NOW))<jupyter_output>Inversion done in 0.06637048721313477
Inversion done in 0.12478208541870117
Inversion done in 0.06429123878479004
Inversion done in 0.06825780868530273
computed r: 12.1013321704435      time: 55.27284812927246
Inversion done in 0.08541512489318848
Inversion done in 0.06437945365905762
Inversion done in 0.07300662994384766
Inversion done in 0.06478047370910645
computed r: 12.2206422540283      time: 54.621328830718994
Inversion done in 0.06634902954101562
Inversion done in 0.06354093551635742
Inversion done in 0.07012414932250977
Inversion done in 0.06394267082214355
computed r: 10.8987101117853      time: 54.14669919013977
Inversion done in 0.06333732604980469
Inversion done in 0.06509089469909668
Inversion done in 0.06349635124206543
Inversion done in 0.06322932243347168
computed r: 12.2228958502830      time: 53.6555712223053
Inversion done in 0.06446576118469238
Inversion done in 0.061765432357788086
Inversion done in 0.06392073631286621
Inversion done in 0.06449246406555176
compute[...]<jupyter_text>Ah, or maybe everything is improved if we use the power 1 in the error function, instead of the power 2. I had thought finding a minimum would be superior to finding a zero, but not with the secant method setup I suppose.<jupyter_code>r1 = 12.15
r2 = 12.2
rlist = [r1, r2]
NOW = time.time()
now = time.time()
for idx in range(10):
    rnew = secant_method_single_iteration(rlist[-2], rlist[-1], Q=20, power=1)
    rlist.append(rnew)
    print("computed r: {}      time: {}".format(rnew, time.time() - now))
    now = time.time()
print("Done in {}".format(time.time() - NOW))
def make_single_matrix(size, Y, r, withtime=True):
    if withtime:
        now = time.time()
    Q = size
    zjlist = make_zj_list(Y, Q)
    zjstarlist = make_zjstar_list(zjlist)
    if withtime:
        print("lists made in {}".format(time.time() - now))
    matrix = []
    B = []
    for n in range(-Q, Q+1):
        if withtime:
            now = time.time()
        if n == 0:
            continue
        row = []
        for l in range(-Q, Q+1):
            if l == 0:
                continue
            entry = CDF(V(n, l, zjlist, CDF(r), CDF(Q), zjstarlist=zjstarlist))
            if n == l:
                #entry = entry - sqrt(Y) * bessel_K(i*r, 2 * pi * abs(n) * Y)
                entry = entry - CDF(sqrt(Y)) * kappa(r, 2 * pi * abs(n) * Y)
            if l != 1:
                row.append(entry.n())
            else:
                B.append([-entry.n()])
                #if n == 1:
                #    print(V(n, l, zjlist, r, Q, zjstarlist=zjstarlist))
                #    print(sqrt(Y) * kappa(r, 2 * pi * abs(n) * Y))
        matrix.append(row)
        if withtime:
            print("row {} computed in {} seconds".format(n, time.time() - now))
    return matrix, B
make_single_matrixe_matrixl<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13233264320220842, 'fraction_numerical': 0.1267540832758224, 'mean_word_length': 3.043953488372093, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '518014', 'license': 'MIT', 'n_tokens_mistral': 7986, 'n_tokens_neox': 6635, 'n_words': 1797, 'path': 'Experimentation/maass_computational_experiments.ipynb', 'repo_name': 'davidlowryduda/notebooks', 'size': '17388'}