Mayor órbita de la sucesión de Collatz
Publicado el 24 de marzo de 2024 por José A. Alonso

13, 40, 20, 10, 5, 16, 8, 4, 2, 1, 4, 2, 1,...

mayoresGeneradores :: Integer -> [Integer]

mayoresGeneradores 20      ==  [18,19]
mayoresGeneradores (10^6)  ==  [837799]

module Mayor_orbita_de_la_sucesion_de_Collatz where

import qualified Data.MemoCombinators as Memo (integral)
import Data.List (genericLength, genericTake)
import Test.Hspec (Spec, describe, hspec, it, shouldBe)
import Test.QuickCheck

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

mayoresGeneradores1 :: Integer -> [Integer]
mayoresGeneradores1 n =
  [x | (x,y) <- ps, y == m]
  where ps = genericTake n longitudesOrbitas
        m  = maximum (map snd ps)

-- longitudesOrbita es la lista de los números junto a las longitudes de
-- las órbitas de Collatz que generan. Por ejemplo,
--    λ> take 10 longitudesOrbitas
--    [(1,1),(2,2),(3,8),(4,3),(5,6),(6,9),(7,17),(8,4),(9,20),(10,7)]
longitudesOrbitas :: [(Integer, Integer)]
longitudesOrbitas =
  [(n, genericLength (collatz n)) | n <- [1..]]

-- (siguiente n) es el siguiente de n en la sucesión de Collatz. Por
-- ejemplo,
--    siguiente 13  ==  40
--    siguiente 40  ==  20
siguiente :: Integer -> Integer
siguiente n | even n    = n `div` 2
            | otherwise = 3*n+1

-- (collatz1 n) es la órbita de Collatz de n hasta alcanzar el
-- 1. Por ejemplo,
--    collatz 13  ==  [13,40,20,10,5,16,8,4,2,1]

-- 1ª definición de collatz
collatz1 :: Integer -> [Integer]
collatz1 1 = [1]
collatz1 n = n : collatz1 (siguiente n)

-- 2ª definición de collatz
collatz2 :: Integer -> [Integer]
collatz2 n = takeWhile (/=1) (iterate siguiente n) ++ [1]

-- Usaremos la 2ª definición de collatz
collatz :: Integer -> [Integer]
collatz = collatz2

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

mayoresGeneradores2 :: Integer -> [Integer]
mayoresGeneradores2 n =
  [x | (x,y) <- ps, y == m]
  where ps = [(x, longitudOrbita x) | x <- [1..n]]
        m  = maximum (map snd ps)

-- (longitudOrbita x) es la longitud de la órbita de x. Por ejemplo,
--    longitudOrbita 13  ==  10
longitudOrbita :: Integer -> Integer
longitudOrbita 1 = 1
longitudOrbita x = 1 + longitudOrbita (siguiente x)

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

mayoresGeneradores3 :: Integer -> [Integer]
mayoresGeneradores3 n =
  [x | (x,y) <- ps, y == m]
  where ps = [(x, longitudOrbita2 x) | x <- [1..n]]
        m  = maximum (map snd ps)

longitudOrbita2 :: Integer -> Integer
longitudOrbita2 = Memo.integral longitudOrbita2'
  where
    longitudOrbita2' 1 = 1
    longitudOrbita2' x = 1 + longitudOrbita2 (siguiente x)

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

verifica :: IO ()
verifica = hspec spec

specG :: (Integer -> [Integer]) -> Spec
specG mayoresGeneradores = do
  it "e1" $
    mayoresGeneradores 20 `shouldBe` [18,19]

spec :: Spec
spec = do
  describe "def. 1" $ specG mayoresGeneradores1
  describe "def. 2" $ specG mayoresGeneradores2
  describe "def. 3" $ specG mayoresGeneradores3

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

-- Equivalencia de definiciones
-- ============================

-- La propiedad es
prop_mayoresGeneradores :: Positive Integer -> Bool
prop_mayoresGeneradores (Positive n) =
  all (== mayoresGeneradores1 n)
      [mayoresGeneradores2 n,
       mayoresGeneradores3 n]

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

-- Comprobación de eficiencia
-- ==========================

-- La comprobación es
--    λ> mayoresGeneradores (10^5)
--    [77031]
--    (5.43 secs, 6,232,320,064 bytes)
--    λ> mayoresGeneradores2 (10^5)
--    [77031]
--    (7.68 secs, 5,238,991,616 bytes)
--    λ> mayoresGeneradores3 (10^5)
--    [77031]
--    (0.88 secs, 571,788,736 bytes)

from functools import lru_cache
from itertools import count, islice
from timeit import Timer, default_timer
from typing import Iterator

from hypothesis import given
from hypothesis import strategies as st

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

# siguiente(n) es el siguiente de n en la sucesión de Collatz. Por
# ejemplo,
#    siguiente(13)  ==  40
#    siguiente(40)  ==  20
def siguiente (n: int) -> int:
    if n % 2 == 0:
        return n // 2
    return 3*n+1

# collatz1(n) es la órbita de Collatz de n hasta alcanzar el
# 1. Por ejemplo,
#    collatz1(13)  ==  [13,40,20,10,5,16,8,4,2,1]
def collatz1(n: int) -> list[int]:
    if n == 1:
        return [1]
    return [n]+ collatz1(siguiente(n))

# longitudesOrbita() es la lista de los números junto a las longitudes de
# las órbitas de Collatz que generan. Por ejemplo,
#    >>> list(islice(longitudesOrbitas(), 10))
#    [(1,1),(2,2),(3,8),(4,3),(5,6),(6,9),(7,17),(8,4),(9,20),(10,7)]
def longitudesOrbitas() -> Iterator[tuple[int, int]]:
    return ((n, len(collatz1(n))) for n in count(1))

def mayoresGeneradores1(n: int) -> list[int]:
    ps = list(islice(longitudesOrbitas(), n))
    m = max((y for (_, y) in ps))
    return [x for (x,y) in ps if y == m]

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

def collatz2(n: int) -> list[int]:
    r = [n]
    while n != 1:
        n = siguiente(n)
        r.append(n)
    return r

def longitudesOrbitas2() -> Iterator[tuple[int, int]]:
    return ((n, len(collatz2(n))) for n in count(1))

def mayoresGeneradores2(n: int) -> list[int]:
    ps = list(islice(longitudesOrbitas2(), n))
    m = max((y for (_, y) in ps))
    return [x for (x,y) in ps if y == m]

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

# longitudOrbita(x) es la longitud de la órbita de x. Por ejemplo,
#    longitudOrbita(13)  ==  10
def longitudOrbita(x: int) -> int:
    if x == 1:
        return 1
    return 1 + longitudOrbita(siguiente(x))

def mayoresGeneradores3(n: int) -> list[int]:
    ps = [(x, longitudOrbita(x)) for x in range(1, n+1)]
    m = max((y for (_, y) in ps))
    return [x for (x,y) in ps if y == m]

# 4ª solución
# ===========

# longitudOrbita2(x) es la longitud de la órbita de x. Por ejemplo,
#    longitudOrbita2(13)  ==  10
def longitudOrbita2(x: int) -> int:
    r = 0
    while x != 1:
        x = siguiente(x)
        r += 1
    return r + 1

def mayoresGeneradores4(n: int) -> list[int]:
    ps = [(x, longitudOrbita2(x)) for x in range(1, n+1)]
    m = max((y for (_, y) in ps))
    return [x for (x,y) in ps if y == m]

# 5ª solución
# ===========

@lru_cache(maxsize=None)
def longitudOrbita3(x: int) -> int:
    if x == 1:
        return 1
    return 1 + longitudOrbita3(siguiente(x))

def mayoresGeneradores5(n: int) -> list[int]:
    ps = [(x, longitudOrbita3(x)) for x in range(1, n+1)]
    m = max((y for (_, y) in ps))
    return [x for (x,y) in ps if y == m]

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

def test_mayoresGeneradores() -> None:
    for mayoresGeneradores in [mayoresGeneradores1,
                               mayoresGeneradores2,
                               mayoresGeneradores3,
                               mayoresGeneradores4,
                               mayoresGeneradores5]:
        assert mayoresGeneradores(20) == [18,19]
    print("Verificado")

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

# Equivalencia de definiciones
# ============================

# La propiedad es
@given(st.integers(min_value=1, max_value=1000))
def test_mayoresGeneradores_equiv(n: int) -> None:
    r = mayoresGeneradores1(n)
    assert mayoresGeneradores2(n) == r
    assert mayoresGeneradores3(n) == r

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

# Comprobació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 comprobación es
#    >>> tiempo('mayoresGeneradores1(10**5)')
#    4.08 segundos
#    >>> tiempo('mayoresGeneradores2(10**5)')
#    1.95 segundos
#    >>> tiempo('mayoresGeneradores3(10**5)')
#    2.16 segundos
#    >>> tiempo('mayoresGeneradores4(10**5)')
#    1.71 segundos
#    >>> tiempo('mayoresGeneradores5(10**5)')
#    0.14 segundos