Día: 1 octubre, 2022

   sumaDivisores :: Integer -> Integer

   sumaDivisores 12                 ==  28
   sumaDivisores 25                 ==  31
   sumaDivisores (product [1..25])  ==  93383273455325195473152000
   length (show (sumaDivisores (product [1..30000])))  ==  121289
   maximum (map sumaDivisores [1..2*10^6])             ==  8851392

import Data.List (foldl', genericLength, group, inits)
import Data.Set (toList)
import Data.Numbers.Primes (primeFactors)
import Math.NumberTheory.ArithmeticFunctions (divisors, sigma)
import Test.QuickCheck
 
-- 1ª solución
-- ===========
 
sumaDivisores1 :: Integer -> Integer
sumaDivisores1 n = sum (divisores1 n)
 
-- (divisores x) es la lista de los divisores de x. Por ejemplo,
--    divisores 60  ==  [1,5,3,15,2,10,6,30,4,20,12,60]
divisores1 :: Integer -> [Integer]
divisores1 n = [x | x <- [1..n], n `rem` x == 0]
 
-- 2ª solución
-- ===========
 
-- Sustituyendo la definición de divisores de la solución anterior por
-- cada una de las del ejercicio [Divisores de un número](https://bit.ly/3S1HYwi)
-- Se obtiene una nueva definición de sumaDivisores. La usada en la
-- definición anterior es la menos eficiente y la que se usa en la
-- siguiente definición es la más eficiente.
 
sumaDivisores2 :: Integer -> Integer
sumaDivisores2 = sum . divisores2
 
divisores2 :: Integer -> [Integer]
divisores2 = toList . divisors
 
-- 3ª solución
-- ===========
 
-- La solución anterior se puede simplificar
 
sumaDivisores3 :: Integer -> Integer
sumaDivisores3 = sum . divisors
 
-- 4ª solución
-- ===========
 
sumaDivisores4 :: Integer -> Integer
sumaDivisores4 = foldl' (+) 0 . divisores2
 
-- 5ª solución
-- ===========
 
sumaDivisores5 :: Integer -> Integer
sumaDivisores5 n = aux [1..n]
  where aux [] = 0
        aux (x:xs) | n `rem` x == 0 = x + aux xs
                   | otherwise      = aux xs
 
-- 6ª solución
-- ===========
 
sumaDivisores6 :: Integer -> Integer
sumaDivisores6 = sum
               . map (product . concat)
               . mapM inits
               . group
               . primeFactors
 
-- 7ª solución
-- ===========
 
-- Si la descomposición de x en factores primos es
--    x = p(1)^e(1) . p(2)^e(2) . .... . p(n)^e(n)
-- entonces la suma de los divisores de x es
--    p(1)^(e(1)+1) - 1     p(2)^(e(2)+1) - 1       p(n)^(e(2)+1) - 1
--   ------------------- . ------------------- ... -------------------
--        p(1)-1                p(2)-1                  p(n)-1
-- Ver la demostración en http://bit.ly/2zUXZPc
 
sumaDivisores7 :: Integer -> Integer
sumaDivisores7 x =
  product [(p^(e+1)-1) `div` (p-1) | (p,e) <- factorizacion x]
 
-- (factorizacion x) es la lista de las bases y exponentes de la
-- descomposición prima de x. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion = map primeroYlongitud . group . primeFactors
 
-- (primeroYlongitud xs) es el par formado por el primer elemento de xs
-- y la longitud de xs. Por ejemplo,
--    primeroYlongitud [3,2,5,7] == (3,4)
primeroYlongitud :: [a] -> (a,Integer)
primeroYlongitud (x:xs) =
  (x, 1 + genericLength xs)
 
-- 8ª solución
-- ===========
 
sumaDivisores8 :: Integer -> Integer
sumaDivisores8 = sigma 1
 
-- Comprobación de equivalencia
-- ============================
 
-- La propiedad es
prop_sumaDivisores :: Positive Integer -> Bool
prop_sumaDivisores (Positive x) =
  all (== sumaDivisores1 x)
      [ sumaDivisores2 x
      , sumaDivisores3 x
      , sumaDivisores4 x
      , sumaDivisores5 x
      , sumaDivisores6 x
      , sumaDivisores7 x
      , sumaDivisores8 x
      ]
 
-- La comprobación es
--    λ> quickCheck prop_sumaDivisores
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
-- La comparación es
--    λ> sumaDivisores1 5336100
--    21386001
--    (2.25 secs, 1,067,805,248 bytes)
--    λ> sumaDivisores2 5336100
--    21386001
--    (0.01 secs, 659,112 bytes)
--    λ> sumaDivisores3 5336100
--    21386001
--    (0.01 secs, 635,688 bytes)
--    λ> sumaDivisores4 5336100
--    21386001
--    (0.01 secs, 648,992 bytes)
--    λ> sumaDivisores5 5336100
--    21386001
--    (2.44 secs, 1,323,924,176 bytes)
--    λ> sumaDivisores6 5336100
--    21386001
--    (0.01 secs, 832,104 bytes)
--    λ> sumaDivisores7 5336100
--    21386001
--    (0.01 secs, 571,040 bytes)
--    λ> sumaDivisores8 5336100
--    21386001
--    (0.00 secs, 558,296 bytes)
--
--    λ> sumaDivisores2 251888923423315469521109880000000
--    1471072204661054993275791673480320
--    (2.30 secs, 1,130,862,080 bytes)
--    λ> sumaDivisores3 251888923423315469521109880000000
--    1471072204661054993275791673480320
--    (1.83 secs, 896,386,232 bytes)
--    λ> sumaDivisores4 251888923423315469521109880000000
--    1471072204661054993275791673480320
--    (1.52 secs, 997,992,328 bytes)
--    λ> sumaDivisores6 251888923423315469521109880000000
--    1471072204661054993275791673480320
--    (2.35 secs, 5,719,848,600 bytes)
--    λ> sumaDivisores7 251888923423315469521109880000000
--    1471072204661054993275791673480320
--    (0.00 secs, 628,136 bytes)
--    λ> sumaDivisores8 251888923423315469521109880000000
--    1471072204661054993275791673480320
--    (0.00 secs, 591,352 bytes)
--
--    λ> length (show (sumaDivisores7 (product [1..30000])))
--    121289
--    (2.76 secs, 4,864,576,304 bytes)
--    λ> length (show (sumaDivisores8 (product [1..30000])))
--    121289
--    (1.65 secs, 3,173,319,312 bytes)

from operator import mul
from functools import reduce
from timeit import Timer, default_timer
from sys import setrecursionlimit
from sympy import divisors, divisor_sigma, factorint
from hypothesis import given, strategies as st
 
setrecursionlimit(10**6)
 
# 1ª solución
# ===========
 
# divisores(x) es la lista de los divisores de x. Por ejemplo,
#    divisores(60)  ==  [1,5,3,15,2,10,6,30,4,20,12,60]
def divisores(n: int) -> list[int]:
    return [x for x in range(1, n + 1) if n % x == 0]
 
def sumaDivisores1(n: int) -> int:
    return sum(divisores(n))
 
# 2ª solución
# ===========
 
# Sustituyendo la definición de divisores de la solución anterior por
# cada una de las del ejercicio Divisores de un número https://bit.ly/3S1HYwi)
# Se obtiene una nueva definición de sumaDivisores. La usada en la
# definición anterior es la menos eficiente y la que se usa en la
# siguiente definición es la más eficiente.
def sumaDivisores2(n: int) -> int:
    return sum(divisors(n))
 
# 3ª solución
# ===========
 
def sumaDivisores3(n: int) -> int:
    def aux(xs: list[int]) -> int:
        if xs:
            if n % xs[0] == 0:
                return xs[0] + aux(xs[1:])
            return aux(xs[1:])
        return 0
 
    return aux(list(range(1, n + 1)))
 
# 4ª solución
# ===========
 
# Si la descomposición de x en factores primos es
#    x = p(1)^e(1) . p(2)^e(2) . .... . p(n)^e(n)
# entonces la suma de los divisores de x es
#    p(1)^(e(1)+1) - 1     p(2)^(e(2)+1) - 1       p(n)^(e(2)+1) - 1
#   ------------------- . ------------------- ... -------------------
#        p(1)-1                p(2)-1                  p(n)-1
# Ver la demostración en http://bit.ly/2zUXZPc
 
def sumaDivisores4(n: int) -> int:
    return reduce(mul, [(p ** (e + 1) - 1) // (p - 1)
                        for (p, e) in factorint(n).items()])
 
# 5ª solución
# ===========
 
def sumaDivisores5(n: int) -> int:
    x = 1
    r1 = 0
    r2 = 0
    while x * x < n:
        if n % x == 0:
            r1 += x
            r2 += n // x
        x += 1
    if x * x == n:
        r1 += x
    return r1 + r2
 
# 6ª solución
# ===========
 
def sumaDivisores6(n: int) -> int:
    return divisor_sigma(n, 1)
 
# Comprobación de equivalencia
# ============================
 
# La propiedad es
@given(st.integers(min_value=2, max_value=1000))
def test_sumaDivisores(n):
    r = sumaDivisores1(n)
    assert sumaDivisores2(n) == r
    assert sumaDivisores3(n) == r
    assert sumaDivisores4(n) == r
    assert sumaDivisores5(n) == r
    assert sumaDivisores6(n) == r
 
# La comprobación es
#    src> poetry run pytest -q suma_de_divisores.py
#    1 passed in 0.90s
 
# Comparación de eficiencia
# =========================
 
def tiempo(e):
    """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('sumaDivisores1(5336100)')
#    0.29 segundos
#    >>> tiempo('sumaDivisores2(5336100)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores3(5336100)')
#    Process Python terminado (killed)
#    >>> tiempo('sumaDivisores4(5336100)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores5(5336100)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores6(5336100)')
#    0.00 segundos
#
#    >>> tiempo('sumaDivisores1(2**9 * 3**8 * 5**2)')
#    4.52 segundos
#    >>> tiempo('sumaDivisores2(2**9 * 3**8 * 5**2)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores4(2**9 * 3**8 * 5**2)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores5(2**9 * 3**8 * 5**2)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores6(2**9 * 3**8 * 5**2)')
#    0.00 segundos
#
#    >>> tiempo('sumaDivisores2(2**9 * 3**8 * 5**7 * 7**4)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores4(2**9 * 3**8 * 5**7 * 7**4)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores5(2**9 * 3**8 * 5**7 * 7**4)')
#    3.24 segundos
#    >>> tiempo('sumaDivisores6(2**9 * 3**8 * 5**7 * 7**4)')
#    0.00 segundos
#
#    >>> tiempo('sumaDivisores2(251888923423315469521109880000000)')
#    1.13 segundos
#    >>> tiempo('sumaDivisores4(251888923423315469521109880000000)')
#    0.00 segundos
#    >>> tiempo('sumaDivisores6(251888923423315469521109880000000)')
#    0.00 segundos
#
#    >>> tiempo('sumaDivisores4(reduce(mul, list(range(1, 30000))))')
#    1.89 segundos
#    >>> tiempo('sumaDivisores6(reduce(mul, list(range(1, 30000))))')
#    1.88 segundos