from sympy import factorint from operator import mul def P(n): f = factorint(n) return sorted([f[i] for i in f]) def a046523(n): x=1 while True: if P(n) == P(x): return x else: x+=1 def A(n, k): f = factorint(n) return 1 if n == 1 else reduce(mul, [1 if i%4==k else i**f[i] for i in f]) def a001511(n): return bin(n)[2:][::-1].index("1") + 1 def T(n, m): return ((n + m)**2 - n - 3*m + 2)/2 def a286364(n): return T(a046523(n/A(n, 1)), a046523(n/A(n, 3))) def a(n): return T(a001511(n), a286364(2*n - 1))