-
Notifications
You must be signed in to change notification settings - Fork 0
/
A153731.py
59 lines (49 loc) · 2.02 KB
/
A153731.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#! /usr/bin/env python3
from labmath import primegen # Available via pip (https://pypi.org/project/labmath/)
from math import comb
"""
This program computes the coefficients of the Swinnerton-Dyer polynomials for the sets {2}, {2,3}, {2,3,5}, {2,3,5,7}, ....
We start by observing that the SDP for {2} is x^2 - 2:
"""
a = {0:-2, 2:1}
n = 2
"""
n stores the degree of the polynomial stored in a.
"""
#print(2)
print(-2, 1)
for q in primegen():
if q == 2: continue
#print(q)
"""
We have the coefficients of the SDP for {2, 3, 5, 7, ..., prevprime(q)} stored in the dictionary a.
All odd-index coefficients are 0, so we do not store them; instead, we will access the dictionary by a.get(z, 0).
We proceed to compute the coefficients of the SDP for {2,3,5,7,...,q} in order of increasing degree.
The coefficients will be stored in the dictionary newa.
"""
newa = {}
for l in range(0, 2*n+1, 2):
# We skip all odd l because those turn out to be 0.
total = 0
for k in range(l, 2*n+1):
g = 0
for j in range(k+1):
multiplier = a.get(j,0) * a.get(k-j,0)
if multiplier == 0: continue
subtotal = sum(comb(j, l-m) * comb(k-j, m) * int((-1)**(k-j-m)) for m in range(l+1))
# We do not need the int() there for theory reasons;
# we use it because (-1)**(negative) always returns a float, even when the exponent is an integer.
g += subtotal * multiplier
if k % 2 == 1: assert g == 0
else:
assert (k - l) % 2 == 0
# Because of the if-statement, we have even k; because of that and the fact that l loops over even numbers,
# the floor division in this exponent is in fact a true division.
total += q**((k-l)//2) * g
if total != 0:
newa[l] = total
print(newa[l], end=' ', flush=True)
print()
a = newa
n *= 2
assert n == max(a.keys())