#!/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))