(Python) import sympy as sp
import numpy as np
N = 1600
primes = list(sp.primerange(2, N+1))
table = np.zeros((N+1, len(primes)), dtype=np.uint8)
for n in range(2, N+1):
factors = sp.factorint(n)
for j, p in enumerate(primes):
if p in factors:
table[n, j] = factors[p] % 2
def test(e): # assumes 6 < e <= N
fac1 = np.zeros(len(primes), dtype=np.uint8)
for d in reversed(range(1, e)):
if d+1 in primes:
return False
fac1 += table[d+1]
for c in reversed(range(1, d)):
fac2 = fac1 % 2
for b in reversed(range(1, c)):
if b+1 in primes:
if fac1[primes.index(b+1)] == 0:
break
fac2 += table[b+1]
fac3 = fac2 % 2
for a in range(1, b):
fac3 += table[a]
if np.all(fac3 % 2 == 0):
return True
return False
(Python)
from math import prod
from itertools import count, islice
from sympy import factorint, isprime
def
A389148_gen(): # generator of terms
g, p, r = [set()], 0, {}
for c in count(2):
f = factorint(c)
d = g[-1]^set(i for i, j in f.items() if j&1)
for h in g:
if (t:=prod(d^h)) not in r:
r[t] = c-1
g.append(d)
if isprime(c):
p = c-2
elif any(not c%q for q in (2, 3, 5, 7, 11)):
continue
elif max(f.values())>1:
continue
else:
for b in range(c-2, p, -1):
u = d^g[b]
for a in range(b-1, 1, -1):
e = prod(u^g[a])
if e in r and r[e]<a:
break
else:
continue
break
else:
yield c
