Representaciones de un número como suma de dos cuadrados

José A. Alonso

08-02-2016

Definir la función

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

tal que (representaciones n) es la lista de pares de números naturales (x,y) tales que n = x^2 + y^2. Por ejemplo.

representaciones  20           ==  [(2,4)]
representaciones  25           ==  [(0,5),(3,4)]
representaciones 325           ==  [(1,18),(6,17),(10,15)]
representaciones 100000147984  ==  [(0,316228)]

Comprobar con QuickCheck que un número natural n se puede representar como suma de dos cuadrados si, y sólo si, en la factorización prima de n todos los exponentes de sus factores primos congruentes con 3 módulo 4 son pares.

Soluciones

import Test.QuickCheck
import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)

-- 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]

-- (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

-- (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 x = floor (sqrt (fromIntegral x))

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

representaciones3 :: Integer -> [(Integer, Integer)]
representaciones3 n = aux 0 (floor (sqrt (fromIntegral 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)

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

--    λ> representaciones1 2000
--    [(8,44),(20,40)]
--    (4.08 secs, 613,941,832 bytes)
--    λ> representaciones2 2000
--    [(8,44),(20,40)]
--    (0.01 secs, 0 bytes)
--    λ> representaciones3 2000
--    [(8,44),(20,40)]
--    (0.01 secs, 0 bytes)
--
--    λ> length (representaciones2 1000000000000)
--    7
--    (3.95 secs, 760,612,080 bytes)
--    λ> length (representaciones3 1000000000000)
--    7
--    (3.41 secs, 470,403,832 bytes)

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

-- La propiedad es
prop_representacion :: Integer -> Property
prop_representacion n =
    n > 0 ==>
      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
--    λ> quickCheckWith (stdArgs {maxSize=7}) prop_representacion
--    +++ OK, passed 100 tests.