primosH :: [Integer]

   take 15 primosH     == [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]
   primosH !! (3*10^4) == 313661

module Numeros_primos_de_Hilbert where

import Data.Numbers.Primes (isPrime, primeFactors)
import Test.QuickCheck (NonNegative (NonNegative), quickCheck)
import Test.Hspec (Spec, hspec, it, shouldBe)

-- 1ª solució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
-- ===========

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 ).

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

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

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "e1" $
    take 15 primosH1 `shouldBe` [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]
  it "e2" $
    take 15 primosH2 `shouldBe` [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]
  it "e3" $
    take 15 primosH3 `shouldBe` [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]

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

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

-- La propiedad es
prop_primosH :: NonNegative Int -> Bool
prop_primosH (NonNegative n) =
  all (== primosH1 !! n)
      [primosH2 !! n,
       primosH3 !! n]

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

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

-- La comparación es
--    λ> primosH1 !! 2000
--    16957
--    (2.16 secs, 752,085,752 bytes)
--    λ> primosH2 !! 2000
--    16957
--    (0.03 secs, 19,771,008 bytes)
--    λ> primosH3 !! 2000
--    16957
--    (0.07 secs, 152,029,168 bytes)
--
--    λ> primosH2 !! (3*10^4)
--    313661
--    (1.44 secs, 989,761,888 bytes)
--    λ> primosH3 !! (3*10^4)
--    313661
--    (2.06 secs, 6,554,068,992 bytes)

from itertools import count, islice, takewhile
from math import ceil, sqrt
from timeit import Timer, default_timer
from typing import Iterator

from hypothesis import given
from hypothesis import strategies as st

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

# numerosH es la sucesión de los números de Hilbert. Por ejemplo,
#    >>> list(islice(numerosH(), 15))
#    [1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, 49, 53, 57]
def numerosH() -> Iterator[int]:
    return count(1, 4)

# (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]
def divisoresH(n: int) -> list[int]:
    return [x for x in takewhile(lambda x: x <= n, numerosH())
            if n % x == 0]

def primosH1() -> Iterator[int]:
    return (n for n in islice(numerosH(), 1, None)
            if divisoresH(n) == [1, n])

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

def primosH2() -> Iterator[int]:
    def esPrimoH(n: int) -> bool:
        m = ceil(sqrt(n))
        def noDivideAn(x: int) -> bool:
            return n % x != 0
        return all(noDivideAn(x) for x in range(5, m+1, 4))
    return filter(esPrimoH, islice(numerosH(), 1, None))

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

def test_primosH() -> None:
    assert list(islice(primosH1(), 15)) == \
        [5, 9, 13, 17, 21, 29, 33, 37, 41, 49, 53, 57, 61, 69, 73]
    assert list(islice(primosH2(), 15)) == \
        [5, 9, 13, 17, 21, 29, 33, 37, 41, 49, 53, 57, 61, 69, 73]
    print ("Verificado")

# La verificación es
#    >>> test_primosH()
#    Verificado

# Comprobación de equivalencia
# ============================

# nth(i, n) es el n-ésimo elemento del iterador i. Por ejemplo,
#    nth(primos(), 4) == 11
def nth(i: Iterator[int], n: int) -> int:
    return list(islice(i, n, n+1))[0]

# La propiedad es
@given(st.integers(min_value=1, max_value=1000))
def test_primosH_equiv(n: int) -> None:
    assert nth(primosH1(), n) == nth(primosH2(), n)

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

# Comparación de eficiencia
# =========================

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('nth(primosH1(), 3000)')
#    2.51 segundos
#    >>> tiempo('nth(primosH2(), 3000)')
#    0.05 segundos