import Data.Numbers.Primes (isPrime, primeFactors)
import Test.QuickCheck (Positive (Positive), quickCheck)
-- 1ª solución
-- ===========
factorizacionesH1 :: Integer -> [[Integer]]
factorizacionesH1 = aux primosH1
where
aux (x:xs) n
| x == n = [[n]]
| x > n = []
| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n
| otherwise = aux xs n
primosH1 :: [Integer]
primosH1 = [n | n <- tail numerosH,
divisoresH n == [1,n]]
-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,
-- take 15 numerosH == [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]
numerosH :: [Integer]
numerosH = [1,5..]
-- (divisoresH n) es la lista de los números de Hilbert que dividen a
-- n. Por ejemplo,
-- divisoresH 117 == [1,9,13,117]
-- divisoresH 21 == [1,21]
divisoresH :: Integer -> [Integer]
divisoresH n = [x | x <- takeWhile (<=n) numerosH,
n `mod` x == 0]
-- 2ª solución
-- ===========
factorizacionesH2 :: Integer -> [[Integer]]
factorizacionesH2 = aux primosH2
where
aux (x:xs) n
| x == n = [[n]]
| x > n = []
| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n
| otherwise = aux xs n
primosH2 :: [Integer]
primosH2 = filter esPrimoH (tail numerosH)
where esPrimoH n = all noDivideAn [5,9..m]
where noDivideAn x = n `mod` x /= 0
m = ceiling (sqrt (fromIntegral n))
-- 3ª solución
-- ===========
-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo
-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)
-- (ver en https://bit.ly/3zq7h4e ).
factorizacionesH3 :: Integer -> [[Integer]]
factorizacionesH3 = aux primosH3
where
aux (x:xs) n
| x == n = [[n]]
| x > n = []
| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n
| otherwise = aux xs n
primosH3 :: [Integer]
primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]
where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs
where xs = primeFactors n
-- Comprobación de equivalencia
-- ============================
-- La propiedad es
prop_factorizacionesH :: Positive Integer -> Bool
prop_factorizacionesH (Positive n) =
all (== factorizacionesH1 m)
[factorizacionesH2 m,
factorizacionesH3 m]
where m = 1 + 4 * n
-- La comprobación es
-- λ> quickCheck prop_factorizacionesH
-- +++ OK, passed 100 tests.
-- Comparación de eficiencia
-- =========================
-- La comparación es
-- λ> factorizacionesH1 80109
-- [[9,9,989],[9,69,129]]
-- (42.77 secs, 14,899,787,640 bytes)
-- λ> factorizacionesH2 80109
-- [[9,9,989],[9,69,129]]
-- (0.26 secs, 156,051,104 bytes)
-- λ> factorizacionesH3 80109
-- [[9,9,989],[9,69,129]]
-- (0.35 secs, 1,118,236,536 bytes)
-- Propiedad de factorización
-- ==========================
-- La propiedad es
prop_factorizable :: Positive Integer -> Bool
prop_factorizable (Positive n) =
not (null (factorizacionesH1 (1 + 4 * n)))
-- La comprobación es
-- λ> quickCheck prop_factorizable
-- +++ OK, passed 100 tests.