from primesieve import primes #primesieve is a third-party library available on pypi from math import gcd from functools import reduce N = 100000 lprimes = primes(N) lpdivs = [list() for i in range(N+1)] for p in lprimes: for i in range(p, N+1, p): lpdivs[i].append(p) def _factor(x): for p in lpdivs[x]: po = 0 while x%p==0: po += 1 x//=p yield (p, po) def factor(x): if x==1: return [1] l = [1] for p, po in _factor(x): l.extend([i*p**poi for poi in range(1, po+1) for i in l]) l.sort() return l def isA331828(x): s = [1] for f in factor(x)[1:]: if f not in (s[i]+s[j] for i in range(len(s)) for j in range(i, len(s))): return False s.append(f) return True def generate_A331828(t): x = 1 for i in range(t): while not isA331828(x): x += 1 yield x x += 1