# Python code for OEIS A260550 # Michael S. Branicky, Jan 15 2021 from itertools import product from numba import njit @njit def decomposable(A): M = max(A) if M == 1: return False a, b, c, d = A for e in range(1, min(a, b, M)): for f in range(1, min(a-e+1, b-e+1, M)): for g in range(1, min(c, d, M)): for h in range(1, min(c-g+1, d-g+1, M)): for i in range(1, min(a//e, c//g)+1): for j in range(1, min(b//e, d//g)+1): for k in range(1, min((a-e*i)//f, (c-g*i)//h)+1): for l in range(1, min((b-e*j)//f, (d-g*j)//h)+1): if a != e*i + f*k: continue if b != e*j + f*l: continue if c != g*i + h*k: continue if d != g*j + h*l: continue return True return False def a(n): c = n**4 - (n - 1)**4 # count matrices with at least one 1 return c + sum(not decomposable(A) for A in product(range(2, n+1), repeat=4)) def a2(n): # don't sum over matrices with 1's return sum(not decomposable(A) for A in product(range(2, n+1), repeat=4)) def a2_incremental(n): # incrementally, only sum matrices with an n entry return sum(not decomposable(A) for A in product(range(2, n+1), repeat=4) if any(A[i]==n for i in range(4))) # RUN CODE, PRINT/COMPARE/SAVE RESULTS from time import time time0 = time() data = [1, 15, 75, 237, 559, 1157, 2055, 3471, 5449, 8131, 11633, 16361, 22041, 29349, 38329, 48839, 61325, 76479, 93957, 114717, 138041, 164153, 194505, 229625, 268259, 311031, 359719, 413245, 472145, 537835, 608837, 688121, 774877, 867549, 971403, 1080637, 1198233, 1326059, 1467029, 1617451, 1777881, 1948219, 2132381, 2329081, 2539351] prevrhomb = 0 an = 0 for n in range(1, 1001): # an = a(n) rhomb = n**4 - (n-1)**4 # Rhombic dodecahedral numbers, A005917 # an = a2(n) + rhomb an += (rhomb - prevrhomb) + a2_incremental(n) prevrhomb = rhomb if n < len(data): print(n, an, data[n-1]==an, time()-time0) else: print(n, an, time()-time0) with open('b260550.txt', 'a') as f: f.write(f"{n} {an}\n")