#!/usr/bin/python3 from functools import lru_cache from math import sqrt from sympy.ntheory import jacobi_symbol from sympy.ntheory.factor_ import core, primefactors from sympy.ntheory.multinomial import binomial_coefficients_list # Cache binomial coefficients @lru_cache(maxsize=None) def binomial_row(n): return binomial_coefficients_list(n) def binomial(n, k): return binomial_row(n)[k] # stirling2(n, k) * k! @lru_cache(maxsize=None) def stirling2_fact(n, k): return sum((-1) ** i * binomial(k, i) * (k - i) ** n for i in range(0, k + 1)) def A000182(n): # NB This uses the second formula given by Vladimir Kruchinin, slightly reworked and with offset 1 return (-1) ** (n + 1) * sum((-2) ** (2 * n - 1 - j) * stirling2_fact(2 * n - 1, j) for j in range(1, 2 * n)) @lru_cache(maxsize=None) def d_squarefree(b, n): if b == 1: return A000182(n) # Eqn 22 of the Shanks paper if b % 4 == 1: D = (-1) ** (n - 1) * sum(jacobi_symbol(k, b) * (b - 4 * k) ** (2 * n - 1) for k in range(1, (b + 1) // 2)) else: D = (-1) ** (n - 1) * sum(jacobi_symbol(b, m) * (b - m) ** (2 * n - 1) for m in range(1, b, 2)) # Eqn 21 of the Shanks paper return D - sum((-b * b) ** i * binomial(2 * n - 1, 2 * i) * d(b, n-i) for i in range(1, n)) @lru_cache(maxsize=None) def d(a, n): # Decompose into squarefree part and square part b = core(a) m = int(sqrt(a // b)) result = d_squarefree(b, n) * m ** (4 * n - 1) for p in primefactors(m): if p > 2: p2n = p ** (2 * n) result = result * (p2n - jacobi_symbol(b, p)) // p2n if b == 1 and m > 1: result //= 2 return result def A000061(n): return d(n, 1) def A000176(n): return d(n, 2) def A000488(n): return d(n, 3) def A000518(n): return d(n, 4) # Print A000061 in B-file format if __name__ == "__main__": for n in range(1, 64): print(n, A000061(n))