Kernel: SageMath (latest)

Example of loading fast modular symbols code in a Jupyter notebook

WRONG

Exacty this same code should work in a .py file...

In [2]:

load("modular_symbol_map.pyx")

Compiling ./modular_symbol_map.pyx...

In [3]:

A = ModularSymbols(389,sign=1).cuspidal_subspace().new_subspace().decomposition()[0]
f = ModularSymbolMap(A)

In [4]:

%timeit f._eval1(-3,7)

The slowest run took 21.72 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 1.11 µs per loop

In [0]:

f._denom

Now try $d=29$ as in 11a.ipynb

In [5]:

d = 29
# much more ms, since this code is massively faster...
ms = [m for m in prime_range(3,100000) if gcd(m, 11) == 1 and euler_phi(m) % d == 0]
print(ms)

[59, 233, 349, 523, 929, 1103, 1277, 1451, 1567, 1741, 1973, 2089, 2437, 2843, 3191, 3307, 3539, 4003, 4177, 4409, 4583, 4931, 5279, 5569, 5743, 5801, 6091, 6323, 6961, 7019, 7193, 7309, 7541, 8179, 8237, 8353, 8527, 8933, 9049, 9281, 9397, 9629, 9803, 10093, 10151, 10267, 10499, 10789, 10847, 11311, 11369, 11717, 11833, 12007, 12239, 12413, 12703, 13109, 13399, 13457, 13921, 14153, 14327, 15139, 15313, 15661, 16067, 16183, 16763, 16879, 16937, 17053, 17401, 17749, 17807, 17923, 17981, 18097, 18329, 18503, 18793, 19141, 19373, 19489, 20011, 20359, 20533, 20707, 20939, 21577, 21751, 22157, 22273, 22447, 22621, 22679, 22853, 23027, 23143, 23201, 23549, 24071, 24419, 24593, 24709, 24767, 25057, 25463, 25579, 26449, 26681, 27551, 28537, 28711, 29059, 29581, 30103, 30161, 30509, 31147, 31321, 31379, 31727, 32191, 32423, 32713, 32771, 32887, 33119, 33409, 33641, 33757, 33931, 34337, 34511, 35149, 35323, 35381, 35671, 35729, 36251, 36541, 36599, 36947, 37643, 37991, 38281, 38629, 38803, 38861, 38977, 39209, 39383, 39499, 39847, 40253, 40427, 40543, 40949, 41413, 41761, 42283, 42457, 42689, 42863, 42979, 43037, 44371, 44777, 44893, 45589, 45763, 45821, 46633, 46691, 46807, 47387, 48779, 48953, 49069, 49417, 49823, 49939, 50287, 50461, 50867, 51157, 51563, 51679, 51853, 52027, 52201, 52259, 52433, 53129, 53419, 53593, 54347, 54521, 54869, 55217, 55333, 55681, 56087, 56377, 56783, 56957, 57073, 57131, 57653, 57943, 59219, 59393, 59509, 59567, 60089, 60611, 60727, 60901, 62003, 62119, 62351, 62467, 62873, 62989, 63337, 63743, 64033, 64091, 64381, 64439, 64613, 65309, 65599, 65657, 65831, 66179, 66643, 66701, 67049, 67339, 68209, 68963, 69427, 69833, 70123, 70181, 70297, 70529, 70877, 71167, 71341, 71399, 72211, 72269, 72559, 72617, 72733, 72907, 73951, 74357, 74531, 74821, 75169, 75227, 75401, 76039, 76213, 76387, 76561, 77141, 77431, 77489, 78301, 78649, 78707, 78823, 79229, 79693, 79867, 80273, 80447, 80621, 80737, 80911, 81839, 82013, 82129, 82361, 82651, 82883, 83231, 83579, 83869, 84391, 84449, 84913, 85087, 85667, 86131, 86711, 87407, 87523, 87697, 88741, 88799, 89611, 89669, 89959, 90017, 90191, 90481, 90887, 91583, 91757, 91873, 92221, 92569, 92627, 92801, 93323, 93497, 93787, 94309, 94483, 94541, 94889, 95063, 95527, 95701, 96223, 96281, 97151, 97441, 97499, 97673, 97789, 97847, 98369, 98543, 98717, 99181, 99529, 99761, 99877]

In [6]:

M = ModularSymbols(11,sign=1).cuspidal_submodule()
ms_map = ModularSymbolMap(M)
ms_denom = ZZ(ms_map.denom)
def f(a,b):
    return ms_map._eval1(a,b)[0] / ms_denom

In [7]:

f(1,11)

-1

In [8]:

def alphas(m, d, normalize=True):
    gen = Integers(m)(primitive_root(m))
    n = euler_phi(m)//d
    b = gen^n
    h = gen^d
    if normalize:
        denom = float(sqrt(m*euler_phi(m)*log(m)))
    else:
        denom = 1
    alphas = []
    for i in range(d):
        s = 0
        for j in range(n):
            period = f((b^i * h^j).lift(), m)
            s += period
        alphas.append(s / denom)
    return alphas

In [9]:

print alphas(ms[0], d)

[0.04232831912577409, 0.04232831912577409, 0.04232831912577409, 0.04232831912577409, 0.04232831912577409, 0.04232831912577409, -0.16931327650309635, 0.04232831912577409, 0.04232831912577409, -0.16931327650309635, 0.04232831912577409, 0.04232831912577409, -0.16931327650309635, 0.04232831912577409, -0.16931327650309635, 0.04232831912577409, 0.04232831912577409, -0.16931327650309635, 0.04232831912577409, -0.16931327650309635, 0.04232831912577409, 0.04232831912577409, -0.16931327650309635, 0.04232831912577409, 0.04232831912577409, -0.16931327650309635, 0.04232831912577409, 0.04232831912577409, 0.04232831912577409]

In [10]:

%%time
data = []
for m in ms:
    data += alphas(m, d)

CPU times: user 1min 19s, sys: 96 ms, total: 1min 19s
Wall time: 1min 19s

In [11]:

print len(data)
t = stats.TimeSeries(data)
print t.mean()
t.plot_histogram(bins=200)

10005
-0.0107567619125

In [12]:

t.plot()

In [13]:

stats.TimeSeries(t[:i].mean() for i in range(5,len(t))).plot(
      gridlines='minor', ymax=0)

see https://docs.scipy.org/doc/scipy-0.7.x/reference/generated/scipy.stats.kurtosis.html

In [14]:

import scipy.stats
scipy.stats.kurtosis(t.numpy(), fisher=False)

90.78104263656937

In [15]:

import matplotlib.pyplot as plt
plt.figure(figsize=(14,6))
plt.hist(t.numpy(), bins=150)
plt.show()

Now try $d=97$

In [16]:

d = 97
# much more ms, since this code is massively faster...
ms = [m for m in prime_range(3,100000) if gcd(m, 11) == 1 and euler_phi(m) % d == 0]
print(ms)

[389, 971, 1553, 1747, 3299, 3881, 4463, 4657, 5821, 6791, 8537, 8731, 10477, 11059, 11447, 12611, 14551, 14939, 16103, 16879, 18043, 19013, 19207, 20177, 20759, 21341, 22699, 23087, 23669, 24251, 25609, 25997, 27743, 29101, 29683, 30071, 31817, 33563, 33757, 36473, 37831, 38219, 39383, 42293, 42487, 43457, 43651, 44621, 45979, 47143, 48889, 49277, 50053, 50441, 51217, 52769, 52963, 54709, 55291, 56843, 57037, 59753, 60917, 62081, 63439, 66931, 67901, 68483, 69259, 70229, 70423, 72169, 73721, 75079, 76243, 76631, 77213, 78571, 79153, 80317, 81869, 83227, 84391, 85361, 86137, 86719, 87107, 88853, 90017, 90599, 90793, 91957, 92927, 93703, 96419, 97001, 97583, 97777, 99523]

In [17]:

%%time
data = []
for m in ms:
    data += alphas(m, d)

CPU times: user 25.4 s, sys: 32 ms, total: 25.4 s
Wall time: 25.4 s

In [18]:

print len(data)
t = stats.TimeSeries(data)
print t.mean()
t.plot_histogram(bins=200)

9603
-0.00321230474805

In [19]:

t.variance()

1.6244063751141795e-05

In [20]:

save(t, 't-11a-97')

Now try $d=997$

In [21]:

d = 997
# much more ms, since this code is massively faster...
ms = [m for m in prime_range(3,100000) if gcd(m, 11) == 1 and euler_phi(m) % d == 0]
print(ms)

[3989, 23929, 27917, 45863, 47857, 63809, 75773, 93719, 95713]

In [22]:

%%time
data = []
for m in ms:
    data += alphas(m, d)

CPU times: user 2.26 s, sys: 8 ms, total: 2.27 s
Wall time: 2.26 s

In [25]:

print len(data)
t = stats.TimeSeries(data)
print t.mean()
t.plot_histogram(bins=30)

8973
-0.000309958617831

In [24]:

t.variance()

9.57983832643914e-07

In [27]:

show(t.plot_histogram(bins=100), xmin=-0.01, xmax=0.01)

Now try $d=2017$

In [34]:

d = 2017
# much more ms, since this code is massively faster...
ms = [m for m in prime_range(3,200000) if gcd(m, 11) == 1 and euler_phi(m) % d == 0]
print(ms)

[8069, 36307, 48409, 56477, 72613, 80681, 121021, 129089, 157327, 189599]

In [44]:

%%time
data = []
for m in ms:
    data += alphas(m, d)

CPU times: user 5.18 s, sys: 8 ms, total: 5.19 s
Wall time: 5.17 s

In [45]:

print len(data)
t = stats.TimeSeries(data)
print t.mean()
t.plot_histogram(bins=10)

20170
-0.000148878246339

In [46]:

import scipy.stats
scipy.stats.kurtosis(t.numpy(), fisher=False)

8.863885397374593

In [0]:

Example of loading fast modular symbols code in a Jupyter notebook

WRONG

Now try d=29d=29d=29 as in 11a.ipynb

Now try d=97d=97d=97

Now try d=997d=997d=997

Now try d=2017d=2017d=2017

Now try $d=29$ as in 11a.ipynb

Now try $d=97$

Now try $d=997$

Now try $d=2017$