# AUTHOR = ShreevatsaR
# https://math.stackexchange.com/users/205/shreevatsar
# Published on:
# https://math.stackexchange.com/q/3198978
# Software: SageMath, version 10.3
# Licence: CC BY-SA 4.0

# Modified by: Peter Luschny, 2024/06/29
# By convention we set
#    is_quadratic_residue(a, 0) = True
# to be able to generate the regular triangle A373748 like this:
#
# def A373748_row(n):
#    return [k if is_quadratic_residue(k, n) else -k for k in range(n + 1)]
# for n in range(11): print([n], A373748_row(n))
#
# If instead you prefer that an "ArithmeticError" is thrown
# to remind you that the "factorization of 0 is not defined"
# then comment out the first line in 'is_quadratic_residue',
# which restors the original form as written by ShreevatsaR.


def peel(a, p):
    # Returns (k, b) such that a = p^k * b and b is minimal.
    if a == 0:
        return (1, 0)
    k = 0
    while a % p == 0:
        k += 1
        a = a // p
    return (k, a)

def is_quadratic_residue(a, n):
    # by convention:
    if n == 0: return True

    for p, e in factor(n):
        (k, b) = peel(a % p**e, p)
        if b == 0:
            continue
        if k % 2:
            return False
        if p == 2:
            if e == 1:
                continue
            if b % 4 != 1:
                return False
            if e >= 3 and b % 8 != 1:
                return False
        else:  # Euler's criterion
            if power_mod(b, (p - 1) // 2, p) != 1:
                return False
    return True