#!/usr/bin/env python3
"""
Comprobaciones reproducibles del escrito
«Extracción de raíces cúbicas enteras mediante congruencias.
Reconstrucción y desarrollo de una nota de 2014» (M. Cerdá Bennassar, 2026),
apartado 4.6.

Implementa con aritmética entera exacta las dos variantes del método:
  A) con el grupo inicial G (sección 3);
  B) sin G, mediante el filtro módulo 48 y una reconstrucción modular
     más amplia (apartado 4.4).

Comprobaciones:
  1. Ambas variantes recuperan las raíces de los cubos de 0^3 a 100000^3.
  2. Ambas rechazan los valores n^3 - 1 y n^3 + 1 para 2 <= n <= 30000.
  3. Ambas recuperan 600 raíces aleatorias: diez para cada longitud de
     1 a 60 cifras (semilla fija 20261018).
  4. Fórmulas de los filtros y de las columnas para todas las raíces de 0
     a 100000 (lemas 3 y 5, proposiciones 1, 13 y 14, apartados 3.3 y 4.1).
  5. Criterio de existencia y unicidad de la proposición 10 para todos los
     módulos 1 <= M <= 20, longitudes 1 <= W <= M, restos módulo M y
     extremos inferiores L = -7, 0, 11.

Uso:  python3 verificar_extension.py
Requiere Python 3.8 o posterior; solo usa la biblioteca estándar.
"""

import random
import sys
import time

SEMILLA = 20261018

# ---------------------------------------------------------------------------
# Utilidades
# ---------------------------------------------------------------------------

_primos = []


def primos_2mod3(k):
    """Los k primeros primos mayores que 5 y congruentes con 2 módulo 3."""
    p = _primos[-1] + 1 if _primos else 11
    while len(_primos) < k:
        if p % 3 == 2 and all(p % d for d in range(2, int(p ** 0.5) + 1)):
            _primos.append(p)
        p += 1
    return _primos[:k]


def restos_raiz(N, primos):
    """Resto de la raíz módulo 10 y módulo cada primo (lemas 7 y 8),
    combinados por el teorema chino de los restos. Devuelve (a, M)."""
    a, M = pow(N % 10, 3, 10), 10
    for p in primos:
        ap = pow(N % p, (2 * p - 1) // 3, p)
        t = ((ap - a) * pow(M, -1, p)) % p
        a, M = a + M * t, M * p
    return a, M


ADM12 = {0, 1, 3, 4, 5, 7, 8, 9, 11}
ADM9 = {0, 1, 8}
ADM7 = {0, 1, 6}

# ---------------------------------------------------------------------------
# Variante A: con el grupo inicial G (sección 3)
# ---------------------------------------------------------------------------


def raiz_con_G(N):
    """Devuelve la raíz cúbica entera de N, o None si N no es un cubo."""
    if N == 0:
        return 0
    if N % 12 not in ADM12 or N % 9 not in ADM9 or N % 7 not in ADM7:
        return None
    D = len(str(N))
    g = (D + 2) // 3
    W = 10 ** (g - 1)
    G = N // 10 ** (3 * (g - 1))
    d = max(c for c in range(1, 10) if c ** 3 <= G)
    L = d * W
    if g == 1:
        n0 = d
    else:
        a, M = restos_raiz(N, primos_2mod3(g - 2))
        n0 = L + (a - L) % M
        if n0 >= L + W:
            return None
    return n0 if n0 ** 3 == N else None


# ---------------------------------------------------------------------------
# Variante B: sin G (apartado 4.4)
# ---------------------------------------------------------------------------


def raiz_sin_G(N):
    """Devuelve la raíz cúbica entera de N, o None si N no es un cubo."""
    if N == 0:
        return 0
    D = len(str(N))
    g = (D + 2) // 3
    W = 10 ** (g - 1)
    # Clase de la raíz módulo 12 (proposición 14).
    if N % 2:
        s = N % 12
    else:
        u = N % 48
        if u not in {0, 8, 16, 24, 32, 40}:
            return None
        s = u // 4
    # Menor producto M = 10 p1 ... pt con H = 6M >= 9W.
    t = 0
    while True:
        prod = 10
        for p in primos_2mod3(t):
            prod *= p
        if 6 * prod >= 9 * W:
            break
        t += 1
    a, M = restos_raiz(N, primos_2mod3(t))
    if (a - s) % 2:
        return None
    t0 = ((s - a) // 2 * pow(M // 2, -1, 6)) % 6
    A = a + M * t0
    H = 6 * M
    n0 = W + (A - W) % H
    if n0 >= 10 * W:
        return None
    return n0 if n0 ** 3 == N else None


# ---------------------------------------------------------------------------
# Comprobaciones
# ---------------------------------------------------------------------------

fallos = []


def informe(texto, ok):
    print(("  correcto   " if ok else "  FALLO      ") + texto)
    if not ok:
        fallos.append(texto)


def comprobacion_1():
    print("1. Raíces de los cubos de 0^3 a 100000^3")
    malos_A = malos_B = 0
    for n in range(100001):
        N = n ** 3
        malos_A += raiz_con_G(N) != n
        malos_B += raiz_sin_G(N) != n
    informe("variante con G: 100001 raíces recuperadas", malos_A == 0)
    informe("variante sin G: 100001 raíces recuperadas", malos_B == 0)


def comprobacion_2():
    print("2. Rechazo de n^3 - 1 y n^3 + 1 para 2 <= n <= 30000")
    total = malos_A = malos_B = 0
    for n in range(2, 30001):
        for N in (n ** 3 - 1, n ** 3 + 1):
            total += 1
            malos_A += raiz_con_G(N) is not None
            malos_B += raiz_sin_G(N) is not None
    informe("variante con G: %d valores rechazados" % total, malos_A == 0 and total == 59998)
    informe("variante sin G: %d valores rechazados" % total, malos_B == 0 and total == 59998)


def comprobacion_3():
    print("3. Raíces aleatorias de 1 a 60 cifras (semilla %d)" % SEMILLA)
    rng = random.Random(SEMILLA)
    total = malos_A = malos_B = 0
    for cifras in range(1, 61):
        for _ in range(10):
            n = rng.randrange(10 ** (cifras - 1) if cifras > 1 else 1, 10 ** cifras)
            N = n ** 3
            total += 1
            malos_A += raiz_con_G(N) != n
            malos_B += raiz_sin_G(N) != n
    informe("variante con G: %d raíces recuperadas" % total, malos_A == 0 and total == 600)
    informe("variante sin G: %d raíces recuperadas" % total, malos_B == 0 and total == 600)


def comprobacion_4():
    print("4. Filtros y columnas para todas las raíces de 0 a 100000")
    S9 = {0, 1, 3, 4, 5, 7, 8, 9, 11}
    ok = {k: True for k in ("p1", "l3", "l5", "p13", "p14", "c33", "c41")}
    for n in range(100001):
        N = n ** 3
        r = N % 12
        m = n - r
        q, s = divmod(n, 12)
        # Proposición 1.
        if not (N % 12 in ADM12 and N % 9 in ADM9 and N % 7 in ADM7):
            ok["p1"] = False
        if n >= 1:
            # Lema 3.
            D, gn = len(str(N)), len(str(n))
            if gn != (D + 2) // 3:
                ok["l3"] = False
            # Lema 5.
            G = N // 10 ** (3 * (gn - 1))
            d = max(c for c in range(1, 10) if c ** 3 <= G)
            if d != int(str(n)[0]):
                ok["l5"] = False
        # Proposición 13.
        if not ((m % 12 == 6) == (n % 4 == 2) == (N % 16 == 8)):
            ok["p13"] = False
        if m % 12 not in (0, 6):
            ok["p13"] = False
        # Proposición 14 y regla de la clase módulo 12.
        if n % 2 == 0:
            if (N - 4 * n) % 48 or N % 48 not in {0, 8, 16, 24, 32, 40} \
                    or (N % 48) // 4 != s:
                ok["p14"] = False
        elif N % 12 != s:
            ok["p14"] = False
        # Apartado 3.3.
        esperado = 12 * q - 6 if s == 2 else 12 * q + 6 if s in (6, 10) else 12 * q
        if m != esperado:
            ok["c33"] = False
        # Apartado 4.1.
        if m % 12 == 0 and (n - m) not in S9:
            ok["c41"] = False
        if m % 12 == 6 and (n - m) not in (0, 4, 8):
            ok["c41"] = False
    informe("proposición 1: restos de los cubos módulo 12, 9 y 7", ok["p1"])
    informe("lema 3: número de cifras de la raíz", ok["l3"])
    informe("lema 5: primera cifra de la raíz", ok["l5"])
    informe("proposición 13: tipo de columna y resto módulo 16", ok["p13"])
    informe("proposición 14: clase de la raíz módulo 12", ok["p14"])
    informe("apartado 3.3: corrección m según n mod 12", ok["c33"])
    informe("apartado 4.1: columnas de nueve y tres valores", ok["c41"])


def comprobacion_5():
    print("5. Proposición 10: existencia y unicidad en un intervalo")
    casos = 0
    bien = True
    for M in range(1, 21):
        for W in range(1, M + 1):
            for a in range(M):
                for L in (-7, 0, 11):
                    casos += 1
                    enteros = [x for x in range(L, L + W) if (x - a) % M == 0]
                    delta = (a - L) % M
                    if len(enteros) > 1:
                        bien = False
                    if (delta < W) != (len(enteros) == 1):
                        bien = False
                    if enteros and enteros[0] != L + delta:
                        bien = False
    informe("%d casos comprobados" % casos, bien)


def main():
    inicio = time.time()
    print("Verificación del apartado 4.6\n")
    comprobacion_1()
    comprobacion_2()
    comprobacion_3()
    comprobacion_4()
    comprobacion_5()
    print("\nTiempo: %.1f s" % (time.time() - inicio))
    if fallos:
        print("Han fallado %d comprobaciones." % len(fallos))
        sys.exit(1)
    print("Todas las comprobaciones se han superado.")


if __name__ == "__main__":
    main()