Números duffinianos

José A. Alonso

12-01-2021

Los números duffinianos, llamados así por Richard Duffy, son los números compuestos n que son coprimos con la suma de sus divisores; es decir, n y la suma de los divisores de n no tienen ningún factor primo común.

Por ejemplo, 35 es un número duffiniano ya que la suma de sus divisores es 1 + 5 + 7 + 35 = 48 que es coprimo con 35.

Definir las funciones

esDuffiniano :: Integer -> Bool
duffinianos :: [Integer]

tales que

(esDuffiniano n) se verifica si n es duffiniano. Por ejemplo,

esDuffiniano 35    ==  True
esDuffiniano 2021  ==  True
esDuffiniano 11    ==  False
esDuffiniano 12    ==  False
esDuffiniano (product [1..2*10^4])  ==  False

duffinianos es la sucesión de los números duffinianos. Por ejemplo,

take 12 duffinianos  ==  [4,8,9,16,21,25,27,32,35,36,39,49]
length (takeWhile (<10^5) duffinianos)  ==  24434

Comprobar con QuickCheck que los números de la forma p^k, con p un primo mayor que 2 y k un entero mayor que 1, son duffinianos.

Soluciones

module Numeros_duffinianos where

import Data.List (genericLength, group)
import Data.Numbers.Primes (isPrime, primeFactors, primes)
import Test.QuickCheck (Property, (==>), quickCheck)

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

duffinianos :: [Integer]
duffinianos = filter esDuffiniano [1..]

esDuffiniano :: Integer -> Bool
esDuffiniano n =
  n > 1 &&
  not (isPrime n) &&
  gcd n (sumaDivisores n) == 1

-- (sumaDivisores n) es la suma de los divisores de n. Por ejemplo.
--      sumaDivisores 35  ==  48
sumaDivisores :: Integer -> Integer
sumaDivisores = sum . divisores

-- (divisores n) es la lista de los divisores de n. Por ejemplo,
--      divisores 35  ==  [1,5,7,35]
divisores :: Integer -> [Integer]
divisores n = [x | x <- [1..n]
                 , n `mod` x == 0]

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

duffinianos2 :: [Integer]
duffinianos2 = filter esDuffiniano2 [1..]

esDuffiniano2 :: Integer -> Bool
esDuffiniano2 n =
  n > 1 &&
  not (isPrime n) &&
  gcd n (sumaDivisores2 n) == 1

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

-- (sumaDivisores2 n) es la suma de los divisores de n. Por ejemplo.
--      sumaDivisores2 35  ==  48
sumaDivisores2 :: Integer -> Integer
sumaDivisores2 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)
primeroYlongitud _      = error "No tiene elementos"

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

-- La comparación es
--    λ> esDuffiniano (product [1..11])
--    False
--    (14.09 secs, 7,983,535,608 bytes)
--    λ> esDuffiniano2 (product [1..11])
--    False
--    (0.01 secs, 125,760 bytes)
--
--    λ> head (dropWhile (<10^4) duffinianos)
--    10000
--    (13.45 secs, 8,872,967,976 bytes)
--    λ> head (dropWhile (<10^4) duffinianos2)
--    10000
--    (0.15 secs, 280,668,240 bytes)
--
--    λ> length (takeWhile (<10^4) duffinianos)
--    2370
--    (13.43 secs, 8,966,138,016 bytes)
--    λ> length (takeWhile (<10^4) duffinianos2)
--    2370
--    (0.15 secs, 286,120,048 bytes)

-- Propiedad
-- =========

-- La propiedad es
prop_duffinianos :: Int -> Int -> Property
prop_duffinianos n k =
  n >= 0 && k > 1 ==> esDuffiniano2 (p^k)
  where p = primes !! n

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

-- La función de verificación es
verifica :: IO ()
verifica = quickCheck prop_duffinianos

-- La verificación es
--    λ> verifica
--    +++ OK, passed 100 tests.