CoCalc Shared Filescode / fast-modsym-notebook.ipynb

# 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 [ ]:
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)

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 [ ]: