representaciones :: Integer -> [(Integer,Integer)]

   representaciones  20              ==  [(2,4)]
   representaciones  25              ==  [(0,5),(3,4)]
   representaciones 325              ==  [(1,18),(6,17),(10,15)]
   length (representaciones (10^14)) == 8

module Representaciones_de_un_numero_como_suma_de_dos_cuadrados where

import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
import Test.Hspec (Spec, describe, hspec, it, shouldBe)
import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

representaciones1 :: Integer -> [(Integer,Integer)]
representaciones1 n =
  [(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]

-- 2ª solución
-- ===========

representaciones2 :: Integer -> [(Integer,Integer)]
representaciones2 n =
  [(x,raiz z) | x <- [0..raiz (n `div` 2)],
                let z = n - x*x,
                esCuadrado z]

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer
raiz 0 = 0
raiz 1 = 1
raiz x = aux (0,x)
    where aux (a,b) | d == x    = c
                    | c == a    = a
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c)
              where c = (a+b) `div` 2
                    d = c^2

-- Nota: La siguiente definición de raíz cuadrada falla para números
-- grandes.
--    raiz' :: Integer -> Integer
--    raiz' = floor . sqrt . fromIntegral
-- Por ejemplo,
--    λ> raiz' (10^50)
--    9999999999999998758486016
--    λ> raiz (10^50)
--    10000000000000000000000000

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

-- 3ª solución
-- ===========

representaciones3 :: Integer -> [(Integer, Integer)]
representaciones3 n = aux 0 (raiz n)
  where aux x y
          | x > y     = []
          | otherwise = case compare (x*x + y*y) n of
                          LT -> aux (x + 1) y
                          EQ -> (x, y) : aux (x + 1) (y - 1)
                          GT -> aux x (y - 1)

-- Verificación                                                     --
-- ============

verifica :: IO ()
verifica = hspec spec

specG :: (Integer -> [(Integer, Integer)]) -> Spec
specG representaciones = do
  it "e1" $
    representaciones  20 `shouldBe` [(2,4)]
  it "e2" $
    representaciones  25 `shouldBe` [(0,5),(3,4)]
  it "e3" $
    representaciones 325 `shouldBe` [(1,18),(6,17),(10,15)]


spec :: Spec
spec = do
  describe "def. 1" $ specG representaciones1
  describe "def. 2" $ specG representaciones2
  describe "def. 3" $ specG representaciones3

-- La verificación es
--    λ> verifica
--
--    9 examples, 0 failures

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_representaciones :: Positive Integer -> Bool
prop_representaciones (Positive n) =
  all (== representaciones1 n)
      [representaciones2 n,
       representaciones3 n]

-- La comprobación es
--    λ> quickCheck prop_representaciones
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> representaciones1 4000
--    [(20,60),(36,52)]
--    (4.95 secs, 2,434,929,624 bytes)
--    λ> representaciones2 4000
--    [(20,60),(36,52)]
--    (0.00 secs, 599,800 bytes)
--    λ> representaciones3 4000
--    [(20,60),(36,52)]
--    (0.01 secs, 591,184 bytes)
--
--    λ> length (representaciones2 (10^14))
--    8
--    (6.64 secs, 5,600,837,088 bytes)
--    λ> length (representaciones3 (10^14))
--    8
--    (9.37 secs, 4,720,548,264 bytes)

-- Comprobación de la propiedad
-- ============================

-- La propiedad es
prop_representacion :: Positive Integer -> Bool
prop_representacion (Positive n) =
  not (null (representaciones2 n)) ==
  all (\(p,e) -> p `mod` 4 /= 3 || even e) (factorizacion n)

-- (factorizacion n) es la factorización prima de n. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

-- La comprobación es
--    λ> quickCheck prop_representacion
--    +++ OK, passed 100 tests.

from math import floor, sqrt
from timeit import Timer, default_timer

from hypothesis import given
from hypothesis import strategies as st
from sympy import factorint

# 1ª solución
# ===========

def representaciones1(n: int) -> list[tuple[int, int]]:
    return [(x,y) for x in range(n+1) for y in range(x,n+1) if n == x*x + y*y]

# 2ª solución
# ===========

# raiz(x) es la raíz cuadrada entera de x. Por ejemplo,
#    raiz(25)     ==  5
#    raiz(24)     ==  4
#    raiz(26)     ==  5
#    raiz(10**46) == 100000000000000000000000
def raiz(x: int) -> int:
    def aux(a: int, b: int) -> int:
        c = (a+b) // 2
        d = c**2
        if d == x:
            return c
        if c == a:
            return a
        if d < x:
            return aux (c,b)
        return aux (a,c)
    if x == 0:
        return 0
    if x == 1:
        return 1
    return aux(0, x)

# Nota. La siguiente definición de raíz cuadrada entera falla para
# números grandes. Por ejemplo,
#    >>> raiz2(10**46)
#    99999999999999991611392
def raiz2(x: int) -> int:
    return floor(sqrt(x))

# esCuadrado(x) se verifica si x es un número al cuadrado. Por
# ejemplo,
#    esCuadrado(25)     == True
#    esCuadrado(26)     == False
#    esCuadrado(10**46) == True
#    esCuadrado(10**47) == False
def esCuadrado(x: int) -> bool:
    y = raiz(x)
    return x == y * y

def representaciones2(n: int) -> list[tuple[int, int]]:
    r: list[tuple[int, int]] = []
    for x in range(1 + raiz(n // 2)):
        z = n - x*x
        if esCuadrado(z):
            r.append((x, raiz(z)))
    return r

# 3ª solución
# ===========

def representaciones3(n: int) -> list[tuple[int, int]]:
    r: list[tuple[int, int]] = []
    for x in range(1 + raiz(n // 2)):
        y = n - x*x
        z = raiz(y)
        if y == z * z:
            r.append((x, z))
    return r

# Verificación
# ============

def test_representaciones() -> None:
    assert representaciones1(20) == [(2,4)]
    assert representaciones1(25) == [(0,5),(3,4)]
    assert representaciones1(325) == [(1,18),(6,17),(10,15)]
    assert representaciones2(20) == [(2,4)]
    assert representaciones2(25) == [(0,5),(3,4)]
    assert representaciones2(325) == [(1,18),(6,17),(10,15)]
    assert representaciones3(20) == [(2,4)]
    assert representaciones3(25) == [(0,5),(3,4)]
    assert representaciones3(325) == [(1,18),(6,17),(10,15)]
    print("Verificado")

# La comprobación es
#    >>> test_representaciones()
#    Verificado

# Equivalencia de las definiciones de representaciones
# ====================================================

# La propiedad es
@given(st.integers(min_value=1, max_value=1000))
def test_representaciones_equiv(x: int) -> None:
    xs = representaciones1(x)
    assert representaciones2(x) == xs
    assert representaciones3(x) == xs

# La comprobación es
#    >>> test_representaciones_equiv()
#    >>>

# Comparación de eficiencia de las definiciones de representaciones
# =================================================================

def tiempo(e: str) -> None:
    """Tiempo (en segundos) de evaluar la expresión e."""
    t = Timer(e, "", default_timer, globals()).timeit(1)
    print(f"{t:0.2f} segundos")

# La comparación es
#    >>> tiempo('representaciones1(5000)')
#    1.27 segundos
#    >>> tiempo('representaciones2(5000)')
#    0.00 segundos
#    >>> tiempo('representaciones3(5000)')
#    0.00 segundos
#
#    >>> tiempo('len(representaciones2(10**12))')
#    11.54 segundos
#    >>> tiempo('len(representaciones3(10**12))')
#    12.08 segundos

# Comprobación de la propiedad
# ============================

# La propiedad es
@given(st.integers(min_value=1, max_value=1000))
def test_representaciones_prop(n: int) -> None:
    factores = factorint(n)
    assert (len(representaciones2(n)) > 0) == \
        all(p % 4 != 3 or e % 2 == 0 for p, e in factores.items())

# La comprobación es
#    >>> test_representaciones_prop()
#    >>>