# Python program for Andersen's algorithm extended to arbitrary base/ordering # Can be used for generating terms in A248915, A259047, A371821, ... # Michael S. Branicky, Apr 12 2024 # NOTES: # 1. uses an adaptation to arbitrary base (using parameter base) and # either ascending (descending = False) or descending (= True) ordering # of the method by J. K. Andersen described in # https://www.primepuzzles.net/puzzles/puzz_472.htm # 2. does not generate all terms, but only a subsequence # 3. can go to higher terms using pre-computed factorizations of # the numbers base^b+1, b = 2..oo, e.g., using factordb.com, # since factoring is the bottleneck # A248915 base = 10 descending = True # use False if ascending # OUTPUT for A248915 # FOUND 1: 12467 (3 base-10 digits) # FOUND 2: 3953318131772867 (9 base-10 digits) # FOUND 3: 32165423814468180271728153293405458790089242105094971116825620121198112217 (39 base-10 digits) # ... # A259047 # base = 10 # descending = False # use True if descending # OUTPUT for A259047 # FOUND 1: 7810053011863508278028459 (13 base-10 digits) # FOUND 2: 243548589812785639501171192374599427569 (21 base-10 digits) # FOUND 3: 156246680967649520709779707010827783165216977094621063970970293 (33 base-10 digits) # ... # A371821 # base = 2 # descending = False # use True if descending # OUTPUT for A371821 # FOUND 1: 85329 (9 base-2 digits) # FOUND 2: 333577497 (15 base-2 digits) # FOUND 3: 55133857902732922904331439521901 (54 base-2 digits) # FOUND 4: 5515275719895841070923935785814145263907960646346230906875207571 (107 base-2 digits) # FOUND 5: 22356328491253661835707832778239168766337726737251132033922408900117132812731486368439225 (150 base-2 digits) # ... verbose = False # use True to print values of b examined from itertools import combinations from math import prod from sympy.ntheory import digits from sympy import factorint, isprime def fromdigits(d, base): n = 0 for di in d: n *= base; n += di return n def ok(n, base, descending): if isprime(n): return False fb = [] fi = list(factorint(n).items()) if descending: fi = fi[::-1] for p, e in fi: fb.extend(digits(p, base)[1:]*e) ic = fromdigits(fb, base) return ic%n == 0 ans = [] for b in range(2, 1001): t = base**b + 1 fb = [] fi = list(factorint(t).items()) if descending: fi = fi[::-1] for p, e in fi: fb.extend([tuple(digits(p, base)[1:])]*e) for m in range(1, len(fb)+1): for c in combinations(fb, m): bP = tuple() for pi in c: bP += pi if len(bP) == b: P = fromdigits(bP, base) if isprime(P): pp = prod([fromdigits(pi, base) for pi in c]) bcsdn = P*pp if bcsdn not in ans: ans.append(bcsdn) ans = sorted(set(ans)) print(f"FOUND {len(ans)}: {bcsdn} ({b} base-{base} digits)") assert ok(bcsdn, base, descending) if verbose: print("...", b) a