Search for completely multiplicative sequences
From Polymath Wiki
This is a simple Python script to search for all completely multiplicative sequences of successively greater length. At the top of the script, 'N' is the limit of the search, and 'limit' is the constraint on the sums. It prints out lists of primes where the sequence is [math]\displaystyle{ -1 }[/math].
N = 300
limit = 2
# fisrt use Eratosthenes to find the primes up to N
# is_prime[n] = true iff n is prime
is_prime = [None, False] + [True]*(N-1)
p = 2
ok = True
while ok:
for i in range(2*p, N+1, p):
is_prime[i] = False
ok = (True in is_prime[p+1:])
if ok:
p = p+1 + is_prime[p+1:].index(True)
# make list of primes
primes = []
for i in range(1, N+1):
if is_prime[i]:
primes.append(i)
n_primes = len(primes)
# make list a[1..N] where a[n] is a dict with primes as keys and a[n][p] = True
# iff p|n an odd number of times
a = [None]
for n in range(1,N+1):
fact = {}
k = n
for p in primes:
v = False
while k % p == 0:
v = not v
k /= p
fact[p] = v
a.append(fact)
n_seqs = [1]
xx = [[[False]*n_primes, 1]] # 0th item = specification of primes for which x[p]=-1, 1st item = partial sum
c = 1
n = 2
while (n<=N) & (c>0):
print "first sequence:", [primes[i] for i in range(n_primes) if xx[0][0][i]]
new_xs = []
if is_prime[n]:
# ... fork
i = primes.index(n)
for x in xx:
seq = x[0]
tot = x[1]
new_seq = list(seq)
new_seq[i] = True
new_xs.append([new_seq, tot-1])
x[1] = tot+1
else:
# ... factorise
fact = a[n]
for x in xx:
seq = x[0]
v = 1
for i in range(n_primes):
if fact[primes[i]] & seq[i]:
v = -v
x[1] += v
# remove sequences with bad sums
xx = [x for x in xx+new_xs if abs(x[1])<=limit]
c = len(xx)
n_seqs.append(c)
print "\nn =", n
print "number of sequences =", c
n += 1