Día: 2 agosto, 2022

   factorizacionesH :: Integer -> [[Integer]]

  factorizacionesH  25    ==  [[5,5]]
  factorizacionesH  45    ==  [[5,9]]
  factorizacionesH 441    ==  [[9,49],[21,21]]
  factorizacionesH 80109  ==  [[9,9,989],[9,69,129]]

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.