primeFactors

Números primos de Pierpont

PorJosé A. Alonso 18-diciembre-201825-diciembre-2018

Un número primo de Pierpont es un número primo de la forma $2^{u}3^{v}+1$ , para u y v enteros no negativos.

Definir la sucesión

   primosPierpont :: [Integer]

1	primosPierpont :: [Integer]

tal que sus elementos son los números primos de Pierpont. Por ejemplo,

   λ> take 20 primosPierpont
   [2,3,5,7,13,17,19,37,73,97,109,163,193,257,433,487,577,769,1153,1297]
   λ> primosPierpont !! 49
   8503057

λ> take 20 primosPierpont

[2,3,5,7,13,17,19,37,73,97,109,163,193,257,433,487,577,769,1153,1297]

λ> primosPierpont !! 49

8503057

Soluciones

import Data.Numbers.Primes (primes, primeFactors)

primosPierpont :: [Integer]
primosPierpont =
  [n | n <- primes
     , primoPierpont n]

primoPierpont :: Integer -> Bool
primoPierpont n =
  primeFactors (n-1) `contenidoEn` [2,3]

-- (contenidoEn xs ys) se verifica si xs está contenido en ys. Por
-- ejemplo,
--    contenidoEn [2,3,2,2,3] [2,3]  ==  True
--    contenidoEn [2,3,2,2,1] [2,3]  ==  False
contenidoEn :: [Integer] -> [Integer] -> Bool
contenidoEn xs ys =
  all (`elem` ys) xs

import Data.Numbers.Primes (primes, primeFactors)

primosPierpont :: [Integer]

primosPierpont =

[n | n <- primes

, primoPierpont n]

primoPierpont :: Integer -> Bool

primoPierpont n =

primeFactors (n-1) `contenidoEn` [2,3]

-- (contenidoEn xs ys) se verifica si xs está contenido en ys. Por

-- ejemplo,

-- contenidoEn [2,3,2,2,3] [2,3] == True

-- contenidoEn [2,3,2,2,1] [2,3] == False

contenidoEn :: [Integer] -> [Integer] -> Bool

contenidoEn xs ys =

all (`elem` ys) xs

Pensamiento

«La memoria es infiel: no sólo borra y confunde, sino que, a veces, inventa, para desorientarnos.»

Antonio Machado

Grafo de divisibilidad

PorJosé A. Alonso 14-mayo-201825-abril-2021

El grafo de divisibilidad de orden n es el grafo cuyos nodos son los números naturales entre 1 y n, cuyas aristas son los pares (x,y) tales que x divide a y o y divide a x. El coste de cada arista es el cociente entre su mayor y menor elemento.

Definir las siguientes funciones:

   grafoDivisibilidad :: Int -> Grafo Int Int
   coste              :: Int -> Int

1 2	grafoDivisibilidad :: Int -> Grafo Int Int coste :: Int -> Int

tales que

(grafoDivisibilidad n) es el grafo de divisibilidad de orden n. Por ejemplo,

      λ> grafoDivisibilidad 12
      G ND (array (1,12)
                  [(1,[(2,2),(3,3),(4,4),(5,5),(6,6),(7,7),
                       (8,8),(9,9),(10,10),(11,11),(12,12)]),
                   (2,[(1,2),(4,2),(6,3),(8,4),(10,5),(12,6)]),
                   (3,[(1,3),(6,2),(9,3),(12,4)]),
                   (4,[(1,4),(2,2),(8,2),(12,3)]),
                   (5,[(1,5),(10,2)]),
                   (6,[(1,6),(2,3),(3,2),(12,2)]),
                   (7,[(1,7)]),
                   (8,[(1,8),(2,4),(4,2)]),
                   (9,[(1,9),(3,3)]),
                   (10,[(1,10),(2,5),(5,2)]),
                   (11,[(1,11)]),
                   (12,[(1,12),(2,6),(3,4),(4,3),(6,2)])])

λ> grafoDivisibilidad 12

G ND (array (1,12)

[(1,[(2,2),(3,3),(4,4),(5,5),(6,6),(7,7),

(8,8),(9,9),(10,10),(11,11),(12,12)]),

(2,[(1,2),(4,2),(6,3),(8,4),(10,5),(12,6)]),

(3,[(1,3),(6,2),(9,3),(12,4)]),

(4,[(1,4),(2,2),(8,2),(12,3)]),

(5,[(1,5),(10,2)]),

(6,[(1,6),(2,3),(3,2),(12,2)]),

(7,[(1,7)]),

(8,[(1,8),(2,4),(4,2)]),

(9,[(1,9),(3,3)]),

(10,[(1,10),(2,5),(5,2)]),

(11,[(1,11)]),

(12,[(1,12),(2,6),(3,4),(4,3),(6,2)])])

(coste n) es el coste del árbol de expansión mínimo del grafo de divisibilidad de orden n. Por ejemplo,

      coste 12        ==  41
      coste 3000      ==  605305
      coste (2*10^5)  ==  1711798835

coste 12 == 41

coste 3000 == 605305

coste (2*10^5) == 1711798835

Soluciones

import Data.Ix
import Data.List (delete, sort)
import qualified Data.Set as S
import I1M.Grafo
import I1M.Tabla
import Data.Numbers.Primes (primeFactors)

grafoDivisibilidad :: Int -> Grafo Int Int
grafoDivisibilidad n =
  creaGrafo ND (1,n) [(x,y,y `div` x) | y <- [1..n]
                                      , x <- [1..y-1]
                                      , y `mod` x == 0]

-- 1ª solución (con el algoritmo de Kruskal)
-- =========================================

coste1 :: Int -> Int
coste1 n = sum [p | (p,x,y) <- kruskal (grafoDivisibilidad n)]

-- (kruskal g) es el árbol de expansión mínimo del grafo g calculado
-- mediante el algoritmo de Kruskal. Por ejemplo,
--    λ> kruskal (grafoDivisibilidad 12)
--    [(11,1,11),(7,1,7),(5,1,5),(3,3,9),(3,1,3),(2,6,12),(2,5,10),
--     (2,4,8),(2,3,6),(2,2,4),(2,1,2)]
kruskal :: (Ix v, Num p, Ord p) => Grafo v p -> [(p,v,v)]
kruskal g = kruskal' cola                           -- Cola ordenada
                     (tabla [(x,x) | x <- nodos g]) -- Tabla de raices
                     []                             -- Árbol de expansión
                     (length (nodos g) - 1)         -- Aristas por
                                                    -- colocar
    where cola = sort [(p,x,y) | (x,y,p) <- aristas g]

kruskal' ((p,x,y):as) t ae n 
  | n == 0      = ae
  | actualizado = kruskal' as t' ((p,x,y):ae) (n-1)
  | otherwise   = kruskal' as t  ae           n
  where (actualizado,t') = buscaActualiza (x,y) t

-- (raiz t n) es la raíz de n en la tabla t. Por ejemplo,
--    raiz (crea [(1,1),(3,1),(4,3),(5,4),(2,6),(6,6)]) 5  == 1
--    raiz (crea [(1,1),(3,1),(4,3),(5,4),(2,6),(6,6)]) 2  == 6
raiz:: Eq n => Tabla n n -> n -> n
raiz t x | v == x    = v
         | otherwise = raiz t v
         where v = valor t x

-- (buscaActualiza a t) es el par formado por False y la tabla t, si los
-- dos vértices de la arista a tienen la misma raíz en t y el par
-- formado por True y la tabla obtenida añadiéndole a t la arista
-- formada por el vértice de a de mayor raíz y la raíz del vértice de
-- a de menor raíz. Por ejemplo,
--    ghci> let t = crea [(1,1),(2,2),(3,1),(4,1)]
--    ghci> buscaActualiza (2,3) t
--    (True,Tbl [(1,1),(2,1),(3,1),(4,1)])
--    ghci> buscaActualiza (3,4) t
--    (False,Tbl [(1,1),(2,2),(3,1),(4,1)])
buscaActualiza :: (Eq n, Ord n) => (n,n) -> Tabla n n -> (Bool,Tabla n n)
buscaActualiza (x,y) t 
  | x' == y'  = (False, t) 
  | y' <  x'  = (True, modifica (x,y') t)
  | otherwise = (True, modifica (y,x') t)
  where x' = raiz t x 
        y' = raiz t y

-- 2ª solución (con el algoritmo de Prim)
-- ======================================

coste2 :: Int -> Int
coste2 n = sum [p | (p,x,y) <- prim (grafoDivisibilidad n)]

-- (prim g) es el árbol de expansión mínimo del grafo g calculado
-- mediante el algoritmo de Prim. Por ejemplo,
--    λ> prim (grafoDivisibilidad 12)
--    [(11,1,11),(7,1,7),(2,5,10),(5,1,5),(3,3,9),(2,6,12),(2,3,6),
--     (3,1,3),(2,4,8),(2,2,4),(2,1,2)]
prim :: (Ix v, Num p, Ord p) => Grafo v p -> [(p,v,v)]
prim g = prim' [n]              -- Nodos colocados
               ns               -- Nodos por colocar 
               []               -- Árbol de expansión
               (aristas g)      -- Aristas del grafo
         where (n:ns) = nodos g

prim' t [] ae as = ae
prim' t r  ae as = prim' (v':t) (delete v' r) (e:ae) as
  where e@(c,u', v') = minimum [(c,u,v)| (u,v,c) <- as,
                                         u `elem` t, 
                                         v `elem` r]

-- 3ª solución (con el algoritmo de Prim con conjuntos)
-- ====================================================

coste3 :: Int -> Int
coste3 n = sum [p | (p,x,y) <- prim2 (grafoDivisibilidad n)]

-- (prim2 g) es el árbol de expansión mínimo del grafo g calculado
-- mediante el algoritmo de Prim. Por ejemplo,
--    λ> prim2 (grafoDivisibilidad 12)
--    [(11,1,11),(7,1,7),(2,5,10),(5,1,5),(3,3,9),(2,6,12),(2,3,6),
--     (3,1,3),(2,4,8),(2,2,4),(2,1,2)]
prim2 :: (Ix v, Num p, Ord p) => Grafo v p -> [(p,v,v)]
prim2 g = prim2' (S.singleton n)  -- Nodos colocados
                 (S.fromList ns)  -- Nodos por colocar 
                 []               -- Árbol de expansión
                 (aristas g)      -- Aristas del grafo
  where (n:ns) = nodos g

prim2' t r ae as
  | S.null r  = ae
  | otherwise = prim2' (S.insert v' t)
                       (S.delete v' r)
                       (e:ae)
                       as
  where e@(c,u', v') = minimum [(c,u,v)| (u,v,c) <- as,
                                         S.member u t, 
                                         S.member v r]

-- 4ª solución
-- ===========

coste4 :: Int -> Int
coste4 n = sum [head (primeFactors x) | x <- [2..n]]

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

--    λ> coste1 400
--    14923
--    (0.08 secs, 31,336,440 bytes)
--    λ> coste2 400
--    14923
--    (4.54 secs, 220,745,608 bytes)
--    λ> coste3 400
--    14923
--    (0.69 secs, 217,031,144 bytes)
--    λ> coste4 400
--    14923
--    (0.01 secs, 2,192,336 bytes)
--    
--    λ> coste1 2000
--    284105
--    (2.09 secs, 842,601,904 bytes)
--    λ> coste4 2000
--    284105
--    (0.02 secs, 14,586,888 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

import Data.Ix

import Data.List (delete, sort)

import qualified Data.Set as S

import I1M.Grafo

import I1M.Tabla

import Data.Numbers.Primes (primeFactors)

grafoDivisibilidad :: Int -> Grafo Int Int

grafoDivisibilidad n =

creaGrafo ND (1,n) [(x,y,y `div` x) | y <- [1..n]

, x <- [1..y-1]

, y `mod` x == 0]

-- 1ª solución (con el algoritmo de Kruskal)

-- =========================================

coste1 :: Int -> Int

coste1 n = sum [p | (p,x,y) <- kruskal (grafoDivisibilidad n)]

-- (kruskal g) es el árbol de expansión mínimo del grafo g calculado

-- mediante el algoritmo de Kruskal. Por ejemplo,

-- λ> kruskal (grafoDivisibilidad 12)

-- [(11,1,11),(7,1,7),(5,1,5),(3,3,9),(3,1,3),(2,6,12),(2,5,10),

-- (2,4,8),(2,3,6),(2,2,4),(2,1,2)]

kruskal :: (Ix v, Num p, Ord p) => Grafo v p -> [(p,v,v)]

kruskal g = kruskal' cola -- Cola ordenada

(tabla [(x,x) | x <- nodos g]) -- Tabla de raices

[] -- Árbol de expansión

(length (nodos g) - 1) -- Aristas por

-- colocar

where cola = sort [(p,x,y) | (x,y,p) <- aristas g]

kruskal' ((p,x,y):as) t ae n

| n == 0 = ae

| actualizado = kruskal' as t' ((p,x,y):ae) (n-1)

| otherwise = kruskal' as t ae n

where (actualizado,t') = buscaActualiza (x,y) t

-- (raiz t n) es la raíz de n en la tabla t. Por ejemplo,

-- raiz (crea [(1,1),(3,1),(4,3),(5,4),(2,6),(6,6)]) 5 == 1

-- raiz (crea [(1,1),(3,1),(4,3),(5,4),(2,6),(6,6)]) 2 == 6

raiz:: Eq n => Tabla n n -> n -> n

raiz t x | v == x = v

| otherwise = raiz t v

where v = valor t x

-- (buscaActualiza a t) es el par formado por False y la tabla t, si los

-- dos vértices de la arista a tienen la misma raíz en t y el par

-- formado por True y la tabla obtenida añadiéndole a t la arista

-- formada por el vértice de a de mayor raíz y la raíz del vértice de

-- a de menor raíz. Por ejemplo,

-- ghci> let t = crea [(1,1),(2,2),(3,1),(4,1)]

-- ghci> buscaActualiza (2,3) t

-- (True,Tbl [(1,1),(2,1),(3,1),(4,1)])

-- ghci> buscaActualiza (3,4) t

-- (False,Tbl [(1,1),(2,2),(3,1),(4,1)])

buscaActualiza :: (Eq n, Ord n) => (n,n) -> Tabla n n -> (Bool,Tabla n n)

buscaActualiza (x,y) t

| x' == y' = (False, t)

| y' < x' = (True, modifica (x,y') t)

| otherwise = (True, modifica (y,x') t)

where x' = raiz t x

y' = raiz t y

-- 2ª solución (con el algoritmo de Prim)

-- ======================================

coste2 :: Int -> Int

coste2 n = sum [p | (p,x,y) <- prim (grafoDivisibilidad n)]

-- (prim g) es el árbol de expansión mínimo del grafo g calculado

-- mediante el algoritmo de Prim. Por ejemplo,

-- λ> prim (grafoDivisibilidad 12)

-- [(11,1,11),(7,1,7),(2,5,10),(5,1,5),(3,3,9),(2,6,12),(2,3,6),

-- (3,1,3),(2,4,8),(2,2,4),(2,1,2)]

prim :: (Ix v, Num p, Ord p) => Grafo v p -> [(p,v,v)]

prim g = prim' [n] -- Nodos colocados

ns -- Nodos por colocar

[] -- Árbol de expansión

(aristas g) -- Aristas del grafo

where (n:ns) = nodos g

prim' t [] ae as = ae

prim' t r ae as = prim' (v':t) (delete v' r) (e:ae) as

where e@(c,u', v') = minimum [(c,u,v)| (u,v,c) <- as,

u `elem` t,

v `elem` r]

-- 3ª solución (con el algoritmo de Prim con conjuntos)

-- ====================================================

coste3 :: Int -> Int

coste3 n = sum [p | (p,x,y) <- prim2 (grafoDivisibilidad n)]

-- (prim2 g) es el árbol de expansión mínimo del grafo g calculado

-- mediante el algoritmo de Prim. Por ejemplo,

-- λ> prim2 (grafoDivisibilidad 12)

-- [(11,1,11),(7,1,7),(2,5,10),(5,1,5),(3,3,9),(2,6,12),(2,3,6),

-- (3,1,3),(2,4,8),(2,2,4),(2,1,2)]

prim2 :: (Ix v, Num p, Ord p) => Grafo v p -> [(p,v,v)]

prim2 g = prim2' (S.singleton n) -- Nodos colocados

(S.fromList ns) -- Nodos por colocar

[] -- Árbol de expansión

(aristas g) -- Aristas del grafo

where (n:ns) = nodos g

prim2' t r ae as

| S.null r = ae

| otherwise = prim2' (S.insert v' t)

(S.delete v' r)

(e:ae)

where e@(c,u', v') = minimum [(c,u,v)| (u,v,c) <- as,

S.member u t,

S.member v r]

-- 4ª solución

-- ===========

coste4 :: Int -> Int

coste4 n = sum [head (primeFactors x) | x <- [2..n]]

-- Comparación de eficiencia

-- =========================

-- λ> coste1 400

-- 14923

-- (0.08 secs, 31,336,440 bytes)

-- λ> coste2 400

-- 14923

-- (4.54 secs, 220,745,608 bytes)

-- λ> coste3 400

-- 14923

-- (0.69 secs, 217,031,144 bytes)

-- λ> coste4 400

-- 14923

-- (0.01 secs, 2,192,336 bytes)

-- λ> coste1 2000

-- 284105

-- (2.09 secs, 842,601,904 bytes)

-- λ> coste4 2000

-- 284105

-- (0.02 secs, 14,586,888 bytes)

Números compuestos por un conjunto de primos

PorJosé A. Alonso 20-abril-201830-abril-2018

Los números compuestos por un conjunto de primos son los números cuyos factores primos pertenecen al conjunto. Por ejemplo, los primeros números compuestos por [2,5,7] son

   1,2,4,5,7,8,10,14,16,20,25,28,32,35,40,49,50,56,64,70,...

1	1,2,4,5,7,8,10,14,16,20,25,28,32,35,40,49,50,56,64,70,...

El 28 es compuesto ya que sus divisores primos son 2 y 7 que están en [2,5,7].

Definir la función

   compuestos :: [Integer] -> [Integer]

1	compuestos :: [Integer] -> [Integer]

tal que (compuesto ps) es la lista de los números compuestos por el conjunto de primos ps. Por ejemplo,

   λ> take 20 (compuestos [2,5,7])
   [1,2,4,5,7,8,10,14,16,20,25,28,32,35,40,49,50,56,64,70]
   λ> take 20 (compuestos [2,5])
   [1,2,4,5,8,10,16,20,25,32,40,50,64,80,100,125,128,160,200,250]
   λ> take 20 (compuestos [2,3,5])
   [1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36]
   λ> take 20 (compuestos [3,5,7,11,13])
   [1,3,5,7,9,11,13,15,21,25,27,33,35,39,45,49,55,63,65,75]
   λ> take 15 (compuestos [2])
   [1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384]
   λ> compuestos [2,7] !! (10^4)
   57399514149595471961908157955229677377312712667508119466382354072731648
   λ> compuestos [2,3,5] !! (10^5)
   290237644800000000000000000000000000000

λ> take 20 (compuestos [2,5,7])

[1,2,4,5,7,8,10,14,16,20,25,28,32,35,40,49,50,56,64,70]

λ> take 20 (compuestos [2,5])

[1,2,4,5,8,10,16,20,25,32,40,50,64,80,100,125,128,160,200,250]

λ> take 20 (compuestos [2,3,5])

[1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36]

λ> take 20 (compuestos [3,5,7,11,13])

[1,3,5,7,9,11,13,15,21,25,27,33,35,39,45,49,55,63,65,75]

λ> take 15 (compuestos [2])

[1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384]

λ> compuestos [2,7] !! (10^4)

57399514149595471961908157955229677377312712667508119466382354072731648

λ> compuestos [2,3,5] !! (10^5)

290237644800000000000000000000000000000

Soluciones

import Data.Numbers.Primes (primeFactors)

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

compuestos1 :: [Integer] -> [Integer]
compuestos1 ps =
  [n | n <- [1..], esCompuesto ps n]

-- (esCompuesto ps n) se verifica si los factores primos de n pertenecen
-- a ps. Por ejemplo, 
--    esCompuesto [2,3,7]    28  ==  True
--    esCompuesto [2,3,7]   140  ==  False
--    esCompuesto [2,3,5,7] 140  ==  True
esCompuesto :: [Integer] -> Integer -> Bool
esCompuesto ps n =
  subconjunto (primeFactors n) ps

-- (subconjunto xs ys) se verifica si todos los elementos de xs
-- pertenecen a ys. Por ejemplo, 
--    subconjunto [2,7,2] [7,5,2]  ==  True
--    subconjunto [2,7,3] [7,5,2]  ==  False
subconjunto :: Eq a => [a] -> [a] -> Bool
subconjunto xs ys =
  all (`elem` ys) xs

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

compuestos2 :: [Integer] -> [Integer]
compuestos2 ps =
   1 : mezclaTodas (combinaciones ps)

-- (combinaciones ps) es la lista de los productos de cada elemento de
-- ps por los números compuestos con ps. Por ejemplo,
--    λ> take 8 (compuestos4 [2,5,7])
--    [1,2,4,5,7,8,10,14]
--    λ> map (take 6) (combinaciones [2,5,7])
--    [[2,4,8,10,14,16],[5,10,20,25,35,40],[7,14,28,35,49,56]]
combinaciones :: [Integer] -> [[Integer]]
combinaciones ps =
  [[p * q | q <- compuestos2 ps] | p <- ps]

-- (mezclaTodas xss) es la mezcla ordenada de xss, donde tanto xss como
-- sus elementos son listas infinitas ordenadas. Por ejemplo, 
--    λ> take 10 (mezclaTodas [[n,2*n..] | n <- [2..]])
--    [2,3,4,5,6,7,8,9,10,11]
--    λ> take 10 (mezclaTodas [[n,2*n..] | n <- [2,9..]])
--    [2,4,6,8,9,10,12,14,16,18]
mezclaTodas :: [[Integer]] -> [Integer]
mezclaTodas = foldr1 xmezcla
  where xmezcla (x:xs) ys = x : mezcla xs ys

-- (mezcla xs ys) es la mezcla, eliminando repetidos, de las lista
-- ordenadas xs e ys. Por ejemplo,  
mezcla :: [Integer] -> [Integer] -> [Integer]
mezcla []     ys              = ys
mezcla xs     []              = xs
mezcla us@(x:xs) vs@(y:ys) | x == y     = x : mezcla xs ys
                           | x < y      = x : mezcla xs vs
                           | otherwise  = y : mezcla us ys

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

compuestos3 :: [Integer] -> [Integer]
compuestos3 [] = [1]
compuestos3 (p:ps) =
  mezclaTodas [map (*y) (compuestos3 ps) | y <- [p^k | k <- [0..]]]

-- 4ª solución
-- ===========

compuestos4 :: [Integer] -> [Integer]
compuestos4 ps = foldl aux xs (tail ps)
  where p        = head ps
        xs       = [p^k | k <- [0..]]
        aux xs p = mezclaTodas [map (*y) xs | y <- [p^k | k <- [0..]]]

-- 5ª solución
-- ===========

compuestos5 :: [Integer] -> [Integer]
compuestos5 = foldl aux [1] 
  where aux xs p = mezclaTodas [map (*y) xs | y <- [p^k | k <- [0..]]]

-- 6ª solución
-- ===========

compuestos6 :: [Integer] -> [Integer]
compuestos6 xs = aux
  where aux = 1 : mezclas xs aux
        mezclas []     _  = []
        mezclas (x:xs) zs = mezcla (map (x*) zs) (mezclas xs zs)

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

--    λ> compuestos1 [2,3,5] !! 300
--    84375
--    (5.85 secs, 2,961,101,088 bytes)
--    λ> compuestos2 [2,3,5] !! 300
--    84375
--    (3.54 secs, 311,137,952 bytes)
--    λ> compuestos2 [2,3,5] !! 400
--    312500
--    (13.01 secs, 1,229,801,184 bytes)
--    λ> compuestos3 [2,3,5] !! 400
--    312500
--    (0.02 secs, 2,066,152 bytes)
--    λ> compuestos3 [2,3,5] !! 20000
--    15441834907098675000000
--    (1.57 secs, 203,061,864 bytes)
--    λ> compuestos4 [2,3,5] !! 20000
--    15441834907098675000000
--    (0.40 secs, 53,335,080 bytes)
--    λ> compuestos4 [2,3,5] !! 50000
--    2379528690747474604574166220800
--    (1.25 secs, 170,058,496 bytes)
--    λ> compuestos5 [2,3,5] !! 50000
--    2379528690747474604574166220800
--    (1.26 secs, 170,104,648 bytes)
--    λ> compuestos6 [2,3,5] !! 50000
--    2379528690747474604574166220800
--    (0.26 secs, 40,490,280 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

import Data.Numbers.Primes (primeFactors)

-- 1ª solución

-- ===========

compuestos1 :: [Integer] -> [Integer]

compuestos1 ps =

[n | n <- [1..], esCompuesto ps n]

-- (esCompuesto ps n) se verifica si los factores primos de n pertenecen

-- a ps. Por ejemplo,

-- esCompuesto [2,3,7] 28 == True

-- esCompuesto [2,3,7] 140 == False

-- esCompuesto [2,3,5,7] 140 == True

esCompuesto :: [Integer] -> Integer -> Bool

esCompuesto ps n =

subconjunto (primeFactors n) ps

-- (subconjunto xs ys) se verifica si todos los elementos de xs

-- pertenecen a ys. Por ejemplo,

-- subconjunto [2,7,2] [7,5,2] == True

-- subconjunto [2,7,3] [7,5,2] == False

subconjunto :: Eq a => [a] -> [a] -> Bool

subconjunto xs ys =

all (`elem` ys) xs

-- 2ª solución

-- ===========

compuestos2 :: [Integer] -> [Integer]

compuestos2 ps =

1 : mezclaTodas (combinaciones ps)

-- (combinaciones ps) es la lista de los productos de cada elemento de

-- ps por los números compuestos con ps. Por ejemplo,

-- λ> take 8 (compuestos4 [2,5,7])

-- [1,2,4,5,7,8,10,14]

-- λ> map (take 6) (combinaciones [2,5,7])

-- [[2,4,8,10,14,16],[5,10,20,25,35,40],[7,14,28,35,49,56]]

combinaciones :: [Integer] -> [[Integer]]

combinaciones ps =

[[p * q | q <- compuestos2 ps] | p <- ps]

-- (mezclaTodas xss) es la mezcla ordenada de xss, donde tanto xss como

-- sus elementos son listas infinitas ordenadas. Por ejemplo,

-- λ> take 10 (mezclaTodas [[n,2*n..] | n <- [2..]])

-- [2,3,4,5,6,7,8,9,10,11]

-- λ> take 10 (mezclaTodas [[n,2*n..] | n <- [2,9..]])

-- [2,4,6,8,9,10,12,14,16,18]

mezclaTodas :: [[Integer]] -> [Integer]

mezclaTodas = foldr1 xmezcla

where xmezcla (x:xs) ys = x : mezcla xs ys

-- (mezcla xs ys) es la mezcla, eliminando repetidos, de las lista

-- ordenadas xs e ys. Por ejemplo,

mezcla :: [Integer] -> [Integer] -> [Integer]

mezcla [] ys = ys

mezcla xs [] = xs

mezcla us@(x:xs) vs@(y:ys) | x == y = x : mezcla xs ys

| x < y = x : mezcla xs vs

| otherwise = y : mezcla us ys

-- 3ª solución

-- ===========

compuestos3 :: [Integer] -> [Integer]

compuestos3 [] = [1]

compuestos3 (p:ps) =

mezclaTodas [map (*y) (compuestos3 ps) | y <- [p^k | k <- [0..]]]

-- 4ª solución

-- ===========

compuestos4 :: [Integer] -> [Integer]

compuestos4 ps = foldl aux xs (tail ps)

where p = head ps

xs = [p^k | k <- [0..]]

aux xs p = mezclaTodas [map (*y) xs | y <- [p^k | k <- [0..]]]

-- 5ª solución

-- ===========

compuestos5 :: [Integer] -> [Integer]

compuestos5 = foldl aux [1]

where aux xs p = mezclaTodas [map (*y) xs | y <- [p^k | k <- [0..]]]

-- 6ª solución

-- ===========

compuestos6 :: [Integer] -> [Integer]

compuestos6 xs = aux

where aux = 1 : mezclas xs aux

mezclas [] _ = []

mezclas (x:xs) zs = mezcla (map (x*) zs) (mezclas xs zs)

-- Comparación de eficiencia

-- =========================

-- λ> compuestos1 [2,3,5] !! 300

-- 84375

-- (5.85 secs, 2,961,101,088 bytes)

-- λ> compuestos2 [2,3,5] !! 300

-- 84375

-- (3.54 secs, 311,137,952 bytes)

-- λ> compuestos2 [2,3,5] !! 400

-- 312500

-- (13.01 secs, 1,229,801,184 bytes)

-- λ> compuestos3 [2,3,5] !! 400

-- 312500

-- (0.02 secs, 2,066,152 bytes)

-- λ> compuestos3 [2,3,5] !! 20000

-- 15441834907098675000000

-- (1.57 secs, 203,061,864 bytes)

-- λ> compuestos4 [2,3,5] !! 20000

-- 15441834907098675000000

-- (0.40 secs, 53,335,080 bytes)

-- λ> compuestos4 [2,3,5] !! 50000

-- 2379528690747474604574166220800

-- (1.25 secs, 170,058,496 bytes)

-- λ> compuestos5 [2,3,5] !! 50000

-- 2379528690747474604574166220800

-- (1.26 secs, 170,104,648 bytes)

-- λ> compuestos6 [2,3,5] !! 50000

-- 2379528690747474604574166220800

-- (0.26 secs, 40,490,280 bytes)

Alturas primas

PorJosé A. Alonso 10-abril-201826-marzo-2022

Se considera una enumeración de los números primos:

    p(1)=2,  p(2)=3, p(3)=5, p(4)=7, p(5)=11, p(6)=13, p(7)=17,...

1	p(1)=2, p(2)=3, p(3)=5, p(4)=7, p(5)=11, p(6)=13, p(7)=17,...

Dado un entero x > 1, su altura prima es el mayor i tal que el primo p(i) aparece en la factorización de x en números primos. Por ejemplo, la altura prima de 3500 tiene longitud 4, pues 3500=2^2×5^3×7^1 y la de 34 tiene es 7, pues 34 = 2×17. Además, se define la altura prima de 1 como 0.

Definir las funciones

   alturaPrima        :: Integer -> Integer
   alturasPrimas      :: Integer -> [Integer]
   graficaAlturaPrima :: Integer -> IO ()

alturaPrima :: Integer -> Integer

alturasPrimas :: Integer -> [Integer]

graficaAlturaPrima :: Integer -> IO ()

tales que

(alturaPrima x) es la altura prima de x. Por ejemplo,

     alturaPrima 3500  ==  4
     alturaPrima 34    ==  7

1 2	alturaPrima 3500 == 4 alturaPrima 34 == 7

(alturasPrimas n) es la lista de las altura prima de los primeros n números enteros positivos. Por ejemplo,

     alturasPrimas 15  ==  [0,1,2,1,3,2,4,1,2,3,5,2,6,4,3]
     maximum (alturasPrimas 10000)  ==  1229
     maximum (alturasPrimas 20000)  ==  2262
     maximum (alturasPrimas 30000)  ==  3245
     maximum (alturasPrimas 40000)  ==  4203

alturasPrimas 15 == [0,1,2,1,3,2,4,1,2,3,5,2,6,4,3]

maximum (alturasPrimas 10000) == 1229

maximum (alturasPrimas 20000) == 2262

maximum (alturasPrimas 30000) == 3245

maximum (alturasPrimas 40000) == 4203

(graficaAlturaPrima n) dibuja las alturas primas de los números entre 2 y n. Por ejemplo, (graficaAlturaPrima 500) dibuja

Soluciones

import Data.List (genericLength)
import Data.Numbers.Primes (isPrime, primes, primeFactors)
import Data.Array
import Graphics.Gnuplot.Simple

-- 1ª definicioń de alturaPrima
-- ============================

alturaPrima :: Integer -> Integer
alturaPrima 1 = 0
alturaPrima n = indice (mayorFactorPrimo n)

-- (mayorFactorPrimo n) es el mayor factor primo de n. Por ejemplo,
--    mayorFactorPrimo 3500  ==  7
--    mayorFactorPrimo 34    ==  17
mayorFactorPrimo :: Integer -> Integer
mayorFactorPrimo = last . primeFactors

-- (indice p) es el índice de p en la sucesión de los números
-- primos. Por ejemplo,
--    indice 7   ==  4
--    indice 17  ==  7
indice :: Integer -> Integer
indice p = genericLength (takeWhile (<=p) primes)

-- 2ª definicioń de alturaPrima
-- ============================

alturaPrima2 :: Integer -> Integer
alturaPrima2 n = v ! n
  where v = array (1,n) [(i,f i) | i <- [1..n]]
        f 1 = 0
        f k | isPrime k = indice2 k
            | otherwise = v ! (k `div` (head (primeFactors k)))

indice2 :: Integer -> Integer
indice2 p = head [n | (x,n) <- indicesPrimos, x == p]

-- indicesPrimos es la suceción formada por los números primos y sus
-- índices. Por ejemplo,
--    λ> take 10 indicesPrimos
--    [(2,1),(3,2),(5,3),(7,4),(11,5),(13,6),(17,7),(19,8),(23,9),(29,10)]
indicesPrimos :: [(Integer,Integer)]
indicesPrimos = zip primes [1..]

-- 1ª definición de alturasPrimas
-- ==============================

alturasPrimas :: Integer -> [Integer]
alturasPrimas n = map alturaPrima [1..n]

-- 2ª definición de alturasPrimas
-- ==============================

alturasPrimas2 :: Integer -> [Integer]
alturasPrimas2 n = elems v 
  where v = array (1,n) [(i,f i) | i <- [1..n]]
        f 1 = 0
        f k | isPrime k = indice2 k
            | otherwise = v ! (k `div` (head (primeFactors k)))

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

--    λ> maximum (alturasPrimas 20000)
--    2262
--    (29.97 secs, 13,179,984,536 bytes)
--    λ> maximum (alturasPrimas2 20000)
--    2262
--    (2.11 secs, 455,259,448 bytes)

-- Definición de graficaAlturaPrima
-- ================================

graficaAlturaPrima :: Integer -> IO ()
graficaAlturaPrima n =
  plotList [ Key Nothing
           , PNG "Alturas_primas.png"
           ]
           (alturasPrimas2 n)

import Data.List (genericLength)

import Data.Numbers.Primes (isPrime, primes, primeFactors)

import Data.Array

import Graphics.Gnuplot.Simple

-- 1ª definicioń de alturaPrima

-- ============================

alturaPrima :: Integer -> Integer

alturaPrima 1 = 0

alturaPrima n = indice (mayorFactorPrimo n)

-- (mayorFactorPrimo n) es el mayor factor primo de n. Por ejemplo,

-- mayorFactorPrimo 3500 == 7

-- mayorFactorPrimo 34 == 17

mayorFactorPrimo :: Integer -> Integer

mayorFactorPrimo = last . primeFactors

-- (indice p) es el índice de p en la sucesión de los números

-- primos. Por ejemplo,

-- indice 7 == 4

-- indice 17 == 7

indice :: Integer -> Integer

indice p = genericLength (takeWhile (<=p) primes)

-- 2ª definicioń de alturaPrima

-- ============================

alturaPrima2 :: Integer -> Integer

alturaPrima2 n = v ! n

where v = array (1,n) [(i,f i) | i <- [1..n]]

f 1 = 0

f k | isPrime k = indice2 k

| otherwise = v ! (k `div` (head (primeFactors k)))

indice2 :: Integer -> Integer

indice2 p = head [n | (x,n) <- indicesPrimos, x == p]

-- indicesPrimos es la suceción formada por los números primos y sus

-- índices. Por ejemplo,

-- λ> take 10 indicesPrimos

-- [(2,1),(3,2),(5,3),(7,4),(11,5),(13,6),(17,7),(19,8),(23,9),(29,10)]

indicesPrimos :: [(Integer,Integer)]

indicesPrimos = zip primes [1..]

-- 1ª definición de alturasPrimas

-- ==============================

alturasPrimas :: Integer -> [Integer]

alturasPrimas n = map alturaPrima [1..n]

-- 2ª definición de alturasPrimas

-- ==============================

alturasPrimas2 :: Integer -> [Integer]

alturasPrimas2 n = elems v

where v = array (1,n) [(i,f i) | i <- [1..n]]

f 1 = 0

f k | isPrime k = indice2 k

| otherwise = v ! (k `div` (head (primeFactors k)))

-- Comparación de eficiencia

-- =========================

-- λ> maximum (alturasPrimas 20000)

-- 2262

-- (29.97 secs, 13,179,984,536 bytes)

-- λ> maximum (alturasPrimas2 20000)

-- 2262

-- (2.11 secs, 455,259,448 bytes)

-- Definición de graficaAlturaPrima

-- ================================

graficaAlturaPrima :: Integer -> IO ()

graficaAlturaPrima n =

plotList [ Key Nothing

, PNG "Alturas_primas.png"

]

(alturasPrimas2 n)

Exponentes de Hamming

PorJosé A. Alonso 1-febrero-20188-febrero-2018

Los números de Hamming forman una sucesión estrictamente creciente de números que cumplen las siguientes condiciones:

El número 1 está en la sucesión.
Si x está en la sucesión, entonces 2x, 3x y 5x también están.
Ningún otro número está en la sucesión.

Los primeros números de Hamming son 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, …

Los exponentes de un número de Hamming n es una terna (x,y,z) tal que n = 2^x*3^y*5^z. Por ejemplo, los exponentes de 600 son (3,1,2) ya que 600 = 2^x*3^1*5^z.

Definir la sucesión

   sucExponentesHamming :: [(Int,Int,Int)]

1	sucExponentesHamming :: [(Int,Int,Int)]

cuyos elementos son los exponentes de los números de Hamming. Por ejemplo,

   λ> take 21 sucExponentesHamming
   [(0,0,0),(1,0,0),(0,1,0),(2,0,0),(0,0,1),(1,1,0),(3,0,0),
    (0,2,0),(1,0,1),(2,1,0),(0,1,1),(4,0,0),(1,2,0),(2,0,1),
    (3,1,0),(0,0,2),(0,3,0),(1,1,1),(5,0,0),(2,2,0),(3,0,1)]
   λ> sucExponentesHamming !! (5*10^5)
   (74,82,7)

λ> take 21 sucExponentesHamming

[(0,0,0),(1,0,0),(0,1,0),(2,0,0),(0,0,1),(1,1,0),(3,0,0),

(0,2,0),(1,0,1),(2,1,0),(0,1,1),(4,0,0),(1,2,0),(2,0,1),

(3,1,0),(0,0,2),(0,3,0),(1,1,1),(5,0,0),(2,2,0),(3,0,1)]

λ> sucExponentesHamming !! (5*10^5)

(74,82,7)

Soluciones

import Data.Numbers.Primes (primeFactors)

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

sucExponentesHamming :: [(Int,Int,Int)]
sucExponentesHamming = map exponentes hamming

-- (exponentes n) es la terna de exponentes del número de Hamming n. Por
-- ejemplo, 
--    exponentes 600  ==  (3,1,2)
exponentes :: Integer -> (Int,Int,Int)
exponentes x = (length as, length cs, length ds)
  where xs = primeFactors x
        (as,bs) = span (==2) xs
        (cs,ds) = span (==3) bs

-- hamming es la sucesión de los números de Hamming. Por ejemplo,
--    λ> take 21 hamming
--    [1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40]
hamming :: [Integer]
hamming = 1 : mezcla3 [2*i | i <- hamming]  
                      [3*i | i <- hamming]  
                      [5*i | i <- hamming]  

-- mezcla3 xs ys zs es la lista obtenida mezclando las listas ordenadas
-- xs, ys y zs y eliminando los elementos duplicados. Por ejemplo, 
--    mezcla3 [2,4,6,8,10] [3,6,9,12] [5,10]  ==  [2,3,4,5,6,8,9,10,12]
mezcla3 :: Ord a => [a] -> [a] -> [a] -> [a]
mezcla3 xs ys zs = mezcla2 xs (mezcla2 ys zs)  

-- mezcla2 xs ys zs es la lista obtenida mezclando las listas ordenadas
-- xs e ys y eliminando los elementos duplicados. Por ejemplo, 
--    mezcla2 [2,4,6,8,10,12] [3,6,9,12]  ==  [2,3,4,6,8,9,10,12]
mezcla2 :: Ord a => [a] -> [a] -> [a] 
mezcla2 p@(x:xs) q@(y:ys) | x < y     = x:mezcla2 xs q
                          | x > y     = y:mezcla2 p  ys  
                          | otherwise = x:mezcla2 xs ys
mezcla2 []       ys                   = ys
mezcla2 xs       []                   = xs

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

sucExponentesHamming2 :: [(Int,Int,Int)]
sucExponentesHamming2 = map exponentes2 hamming

exponentes2 :: Integer -> (Int,Int,Int)
exponentes2 = aux (0,0,0)
  where aux (a,b,c) 1 = (a,b,c)
        aux (a,b,c) x | mod x 2 == 0 = aux (a+1,b,c) (div x 2)
                      | mod x 3 == 0 = aux (a,b+1,c) (div x 3)
                      | otherwise    = aux (a,b,c+1) (div x 5)

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

sucExponentesHamming3 :: [(Int,Int,Int)]
sucExponentesHamming3 = map exponentes3 hamming

exponentes3 :: Integer -> (Int,Int,Int)
exponentes3 1 = (0,0,0)
exponentes3 x
  | x `mod` 2 == 0 = suma (1,0,0) (descomposicion (x `div` 2))
  | x `mod` 3 == 0 = suma (0,1,0) (descomposicion (x `div` 3))
  | otherwise      = suma (0,0,1) (descomposicion (x `div` 5))
  where suma (x,y,z) (a,b,c) = (x+a,y+b,z+c)

-- 4ª solución
-- ===========

type Terna = (Int,Int,Int)

sucExponentesHamming4 :: [Terna]
sucExponentesHamming4 =
  (0,0,0) : mezclaT3 [(x+1,y,z) | (x,y,z) <- sucExponentesHamming4]
                     [(x,y+1,z) | (x,y,z) <- sucExponentesHamming4]
                     [(x,y,z+1) | (x,y,z) <- sucExponentesHamming4]
 
mezclaT3 :: [Terna] -> [Terna] -> [Terna] -> [Terna]
mezclaT3 t1 t2 t3 = mezclaT2 t1 (mezclaT2 t2 t3)

mezclaT2 :: [Terna] -> [Terna] -> [Terna]
mezclaT2 ts1@((i,j,k):xs) ts2@((a,b,c):ys)
  | x < y     = (i,j,k) : mezclaT2 xs ts2
  | x > y     = (a,b,c) : mezclaT2 ts1 ys
  | otherwise = (i,j,k) : mezclaT2 xs ys
  where x = 2^i*3^j*5^k
        y = 2^a*3^b*5^c

import Data.Numbers.Primes (primeFactors)

-- 1ª solución

-- ===========

sucExponentesHamming :: [(Int,Int,Int)]

sucExponentesHamming = map exponentes hamming

-- (exponentes n) es la terna de exponentes del número de Hamming n. Por

-- ejemplo,

-- exponentes 600 == (3,1,2)

exponentes :: Integer -> (Int,Int,Int)

exponentes x = (length as, length cs, length ds)

where xs = primeFactors x

(as,bs) = span (==2) xs

(cs,ds) = span (==3) bs

-- hamming es la sucesión de los números de Hamming. Por ejemplo,

-- λ> take 21 hamming

-- [1,2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40]

hamming :: [Integer]

hamming = 1 : mezcla3 [2*i | i <- hamming]

[3*i | i <- hamming]

[5*i | i <- hamming]

-- mezcla3 xs ys zs es la lista obtenida mezclando las listas ordenadas

-- xs, ys y zs y eliminando los elementos duplicados. Por ejemplo,

-- mezcla3 [2,4,6,8,10] [3,6,9,12] [5,10] == [2,3,4,5,6,8,9,10,12]

mezcla3 :: Ord a => [a] -> [a] -> [a] -> [a]

mezcla3 xs ys zs = mezcla2 xs (mezcla2 ys zs)

-- mezcla2 xs ys zs es la lista obtenida mezclando las listas ordenadas

-- xs e ys y eliminando los elementos duplicados. Por ejemplo,

-- mezcla2 [2,4,6,8,10,12] [3,6,9,12] == [2,3,4,6,8,9,10,12]

mezcla2 :: Ord a => [a] -> [a] -> [a]

mezcla2 p@(x:xs) q@(y:ys) | x < y = x:mezcla2 xs q

| x > y = y:mezcla2 p ys

| otherwise = x:mezcla2 xs ys

mezcla2 [] ys = ys

mezcla2 xs [] = xs

-- 2ª solución

-- ===========

sucExponentesHamming2 :: [(Int,Int,Int)]

sucExponentesHamming2 = map exponentes2 hamming

exponentes2 :: Integer -> (Int,Int,Int)

exponentes2 = aux (0,0,0)

where aux (a,b,c) 1 = (a,b,c)

aux (a,b,c) x | mod x 2 == 0 = aux (a+1,b,c) (div x 2)

| mod x 3 == 0 = aux (a,b+1,c) (div x 3)

| otherwise = aux (a,b,c+1) (div x 5)

-- 3ª solución

-- ===========

sucExponentesHamming3 :: [(Int,Int,Int)]

sucExponentesHamming3 = map exponentes3 hamming

exponentes3 :: Integer -> (Int,Int,Int)

exponentes3 1 = (0,0,0)

exponentes3 x

| x `mod` 2 == 0 = suma (1,0,0) (descomposicion (x `div` 2))

| x `mod` 3 == 0 = suma (0,1,0) (descomposicion (x `div` 3))

| otherwise = suma (0,0,1) (descomposicion (x `div` 5))

where suma (x,y,z) (a,b,c) = (x+a,y+b,z+c)

-- 4ª solución

-- ===========

type Terna = (Int,Int,Int)

sucExponentesHamming4 :: [Terna]

sucExponentesHamming4 =

(0,0,0) : mezclaT3 [(x+1,y,z) | (x,y,z) <- sucExponentesHamming4]

[(x,y+1,z) | (x,y,z) <- sucExponentesHamming4]

[(x,y,z+1) | (x,y,z) <- sucExponentesHamming4]

mezclaT3 :: [Terna] -> [Terna] -> [Terna] -> [Terna]

mezclaT3 t1 t2 t3 = mezclaT2 t1 (mezclaT2 t2 t3)

mezclaT2 :: [Terna] -> [Terna] -> [Terna]

mezclaT2 ts1@((i,j,k):xs) ts2@((a,b,c):ys)

| x < y = (i,j,k) : mezclaT2 xs ts2

| x > y = (a,b,c) : mezclaT2 ts1 ys

| otherwise = (i,j,k) : mezclaT2 xs ys

where x = 2^i*3^j*5^k

y = 2^a*3^b*5^c

Complemento potencial

PorJosé A. Alonso 22-enero-201829-enero-2018

El complemento potencial de un número entero positivo x es el menor número y tal que el producto de x por y es un una potencia perfecta. Por ejemplo,

el complemento potencial de 12 es 3 ya que 12 y 24 no son potencias perfectas pero 36 sí lo es;
el complemento potencial de 54 es 4 ya que 54, 108 y 162 no son potencias perfectas pero 216 = 6^3 sí lo es.

Definir las funciones

   complemento                 :: Integer -> Integer
   graficaComplementoPotencial :: Integer -> IO ()

1 2	complemento :: Integer -> Integer graficaComplementoPotencial :: Integer -> IO ()

tales que

(complemento x) es el complemento potencial de x; por ejemplo,

     complemento 12     ==  3
     complemento 54     ==  4
     complemento 720    ==  5
     complemento 24000  ==  9
     complemento 2018   ==  2018

complemento 12 == 3

complemento 54 == 4

complemento 720 == 5

complemento 24000 == 9

complemento 2018 == 2018

(graficaComplementoPotencial n) dibuja la gráfica de los complementos potenciales de los n primeros números enteros positivos. Por ejemplo, (graficaComplementoPotencial 100) dibuja

y (graficaComplementoPotencial 500) dibuja

Comprobar con QuickCheck que (complemento x) es menor o igual que x.

Soluciones

import Data.Numbers.Primes     (primeFactors)
import Data.List               (genericLength, group)
import Graphics.Gnuplot.Simple (plotList, Attribute (Key, PNG, Title))
import Test.QuickCheck

complemento :: Integer -> Integer
complemento 1 = 1
complemento x =
  head [y | y <- [1..]
          , esPotenciaPerfecta (x*y)]
  
-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta 36  ==  True
--    esPotenciaPerfecta 72  ==  False
esPotenciaPerfecta :: Integer -> Bool
esPotenciaPerfecta x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de la factorización prima
-- de x. Por ejemplos,
--    exponentes 36  ==  [2,2]
--    exponentes 72  ==  [3,2]
exponentes :: Integer -> [Integer]
exponentes = map snd . factorizacion
  
-- (factorizacion n) es la factorizacion prima de n. Por ejemplo,
--    factorizacion 1400  ==  [(2,3),(5,2),(7,1)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(x,genericLength xs) | xs@(x:_) <- group (primeFactors n)]

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,
--    mcd [4,6,10]  ==  2
--    mcd [4,5,10]  ==  1
mcd :: [Integer] -> Integer
mcd = foldl1 gcd

-- La propiedad es
prop_complemento :: (Positive Integer) -> Bool
prop_complemento (Positive x) =
  complemento x <= x

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

graficaComplementoPotencial :: Integer -> IO ()
graficaComplementoPotencial n =
  plotList [Key Nothing
           , PNG ("Complemento_potencial_"  ++ show n ++ ".png")
           , Title ("(graficaComplementoPotencial " ++ show n ++ ")") 
           ]
           (map complemento [1..n])

import Data.Numbers.Primes (primeFactors)

import Data.List (genericLength, group)

import Graphics.Gnuplot.Simple (plotList, Attribute (Key, PNG, Title))

import Test.QuickCheck

complemento :: Integer -> Integer

complemento 1 = 1

complemento x =

head [y | y <- [1..]

, esPotenciaPerfecta (x*y)]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por

-- ejemplo,

-- esPotenciaPerfecta 36 == True

-- esPotenciaPerfecta 72 == False

esPotenciaPerfecta :: Integer -> Bool

esPotenciaPerfecta x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de la factorización prima

-- de x. Por ejemplos,

-- exponentes 36 == [2,2]

-- exponentes 72 == [3,2]

exponentes :: Integer -> [Integer]

exponentes = map snd . factorizacion

-- (factorizacion n) es la factorizacion prima de n. Por ejemplo,

-- factorizacion 1400 == [(2,3),(5,2),(7,1)]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion n =

[(x,genericLength xs) | xs@(x:_) <- group (primeFactors n)]

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,

-- mcd [4,6,10] == 2

-- mcd [4,5,10] == 1

mcd :: [Integer] -> Integer

mcd = foldl1 gcd

-- La propiedad es

prop_complemento :: (Positive Integer) -> Bool

prop_complemento (Positive x) =

complemento x <= x

-- La comprobación es

-- λ> quickCheck prop_complemento

-- +++ OK, passed 100 tests.

graficaComplementoPotencial :: Integer -> IO ()

graficaComplementoPotencial n =

plotList [Key Nothing

, PNG ("Complemento_potencial_" ++ show n ++ ".png")

, Title ("(graficaComplementoPotencial " ++ show n ++ ")")

]

(map complemento [1..n])

Escalada hasta un primo

PorJosé A. Alonso 18-enero-201819-febrero-2018

Este ejercicio está basado en el artículo La conjetura de la «escalada hasta un primo» publicado esta semana por Miguel Ángel Morales en su blog Gaussianos.

La conjetura de escalada hasta un primo trata, propuesta por John Horton Conway, es sencilla de plantear, pero primero vamos a ver qué es eso de escalar hasta un primo. Tomamos un número cualquiera y lo descomponemos en factores primos (colocados en orden ascendente). Si el número era primo, ya hemos acabado; si no era primo, construimos el número formado por los factores primos y los exponentes de los mismos colocados tal cual salen en la factorización. Con el número obtenido hacemos lo mismo que antes. La escalada finaliza cuando obtengamos un número primo. Por ejemplo, para obtener la escalada prima de 1400, como no es primo, se factoriza (obteniéndose 2^3 * 5^2 * 7) y se unen bases y exponentes (obteniéndose 23527). Con el 23527 se repite el proceso obteniéndose la factorización (7 * 3361) y su unión (73361). Como el 73361 es primo, termina la escalada. Por tanto, la escalada de 1400 es [1400,23527,73361].

La conjetura de Conway sobre «escalada hasta un primo» dice que todo número natural mayor o igual que 2 termina su escalada en un número primo.

Definir las funciones

   escaladaPrima                :: Integer -> [Integer]
   longitudEscaladaPrima        :: Integer -> Integer
   longitudEscaladaPrimaAcotada :: Integer -> Integer -> Integer
   graficaEscalada              :: Integer -> Integer -> IO ()

escaladaPrima :: Integer -> [Integer]

longitudEscaladaPrima :: Integer -> Integer

longitudEscaladaPrimaAcotada :: Integer -> Integer -> Integer

graficaEscalada :: Integer -> Integer -> IO ()

tales que

(escaladaPrima n) es la escalada prima de n. Por ejemplo,

     λ> escaladaPrima 1400
     [1400,23527,73361]
     λ> escaladaPrima 333
     [333,3237,31383,3211317,3217139151,39722974813]
     λ> take 10 (escaladaPrima 20)
     [20,225,3252,223271,297699,399233,715623,3263907,32347303,160720129]
     λ> take 3 (escaladaPrima 13532385396179)
     [13532385396179,13532385396179,13532385396179]
     λ> take 2 (escaladaPrima 45214884853168941713016664887087462487)
     [45214884853168941713016664887087462487,13532385396179]

λ> escaladaPrima 1400

[1400,23527,73361]

λ> escaladaPrima 333

[333,3237,31383,3211317,3217139151,39722974813]

λ> take 10 (escaladaPrima 20)

[20,225,3252,223271,297699,399233,715623,3263907,32347303,160720129]

λ> take 3 (escaladaPrima 13532385396179)

[13532385396179,13532385396179,13532385396179]

λ> take 2 (escaladaPrima 45214884853168941713016664887087462487)

[45214884853168941713016664887087462487,13532385396179]

(longitudEscaladaPrima n) es la longitud de la escalada prima de n. Por ejemplo,

     longitudEscaladaPrima 1400  ==  3
     longitudEscaladaPrima 333   ==  6

1 2	longitudEscaladaPrima 1400 == 3 longitudEscaladaPrima 333 == 6

(longitudEscaladaPrimaAcotada n k) es el mínimo entre la longitud de la escalada prima de n y k. Por ejemplo,

     longitudEscaladaPrimaAcotada 333 10  ==  6
     longitudEscaladaPrimaAcotada 333 4   ==  4
     longitudEscaladaPrimaAcotada 20 4    ==  4

longitudEscaladaPrimaAcotada 333 10 == 6

longitudEscaladaPrimaAcotada 333 4 == 4

longitudEscaladaPrimaAcotada 20 4 == 4

(graficaEscalada n k) dibuja la gráfica de (longitudEscaladaPrimaAcotada x k) para x entre 2 y n. Por ejemplo, (graficaEscalada 120 15) dibuja

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes
import Graphics.Gnuplot.Simple

-- Definición de escaladaPrima 
escaladaPrima :: Integer -> [Integer]
escaladaPrima n
  | isPrime n = [n]
  | otherwise = n : escaladaPrima (siguiente n)

-- (siguiente n) es el siguiente de n en la escalada hacia un primo. Por
-- ejemplo, 
--    siguiente 1400  ==  23527
siguiente :: Integer -> Integer
siguiente = paresAnumero . factorizacion

-- (factorizacion n) es la factorizacion prima de n. Por ejemplo,
--    factorizacion 1400  ==  [(2,3),(5,2),(7,1)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(x,genericLength xs) | xs@(x:_) <- group (primeFactors n)]

-- (paAnumero p) es el número obtenido pegando las dos componentes de p
-- (si la segunda es mayor que 1) o sólo la primera componente, en caso
-- contrario. Por ejemplo, 
--    parAnumero (7,1)  ==  7
--    parAnumero (23,5)  ==  235
parAnumero :: (Integer,Integer) -> Integer
parAnumero (n,1) = n
parAnumero (n,k) = read (show n ++ show k)

-- (paresA numeros ps) es el número obtenido aplicando parAnumero a cada
-- uno de los pares de ps. Por ejemplo,
--    paresAnumero [(2,3),(5,2),(7,1)]  ==  23527
paresAnumero :: [(Integer,Integer)] -> Integer
paresAnumero = read . concatMap (show . parAnumero)

-- Definición de longitudEscaladaPrima
longitudEscaladaPrima :: Integer -> Integer
longitudEscaladaPrima n
  | isPrime n = 1
  | otherwise = 1 + longitudEscaladaPrima (siguiente n)

-- Definición de longitudEscaladaPrimaAcotada
longitudEscaladaPrimaAcotada :: Integer -> Integer -> Integer
longitudEscaladaPrimaAcotada _ 0 = 0
longitudEscaladaPrimaAcotada n k
  | isPrime n = 1
  | otherwise = 1 + longitudEscaladaPrimaAcotada (siguiente n) (k-1)

-- Definición de graficaEscalada 
graficaEscalada :: Integer -> Integer -> IO ()
graficaEscalada n k =
  plotList [ Key Nothing
           , PNG "Escalada_hasta_un_primo.png"
           , Title ("graficaEscalada " ++ show n ++ " " ++ show k)
           , XLabel "Numeros"
           , YLabel "Longitud de escalada"
           ]
           [(x,longitudEscaladaPrimaAcotada x k) | x <- [2..n]]

import Data.List (genericLength, group)

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

-- Definición de escaladaPrima

escaladaPrima :: Integer -> [Integer]

escaladaPrima n

| isPrime n = [n]

| otherwise = n : escaladaPrima (siguiente n)

-- (siguiente n) es el siguiente de n en la escalada hacia un primo. Por

-- ejemplo,

-- siguiente 1400 == 23527

siguiente :: Integer -> Integer

siguiente = paresAnumero . factorizacion

-- (factorizacion n) es la factorizacion prima de n. Por ejemplo,

-- factorizacion 1400 == [(2,3),(5,2),(7,1)]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion n =

[(x,genericLength xs) | xs@(x:_) <- group (primeFactors n)]

-- (paAnumero p) es el número obtenido pegando las dos componentes de p

-- (si la segunda es mayor que 1) o sólo la primera componente, en caso

-- contrario. Por ejemplo,

-- parAnumero (7,1) == 7

-- parAnumero (23,5) == 235

parAnumero :: (Integer,Integer) -> Integer

parAnumero (n,1) = n

parAnumero (n,k) = read (show n ++ show k)

-- (paresA numeros ps) es el número obtenido aplicando parAnumero a cada

-- uno de los pares de ps. Por ejemplo,

-- paresAnumero [(2,3),(5,2),(7,1)] == 23527

paresAnumero :: [(Integer,Integer)] -> Integer

paresAnumero = read . concatMap (show . parAnumero)

-- Definición de longitudEscaladaPrima

longitudEscaladaPrima :: Integer -> Integer

longitudEscaladaPrima n

| isPrime n = 1

| otherwise = 1 + longitudEscaladaPrima (siguiente n)

-- Definición de longitudEscaladaPrimaAcotada

longitudEscaladaPrimaAcotada :: Integer -> Integer -> Integer

longitudEscaladaPrimaAcotada _ 0 = 0

longitudEscaladaPrimaAcotada n k

| isPrime n = 1

| otherwise = 1 + longitudEscaladaPrimaAcotada (siguiente n) (k-1)

-- Definición de graficaEscalada

graficaEscalada :: Integer -> Integer -> IO ()

graficaEscalada n k =

plotList [ Key Nothing

, PNG "Escalada_hasta_un_primo.png"

, Title ("graficaEscalada " ++ show n ++ " " ++ show k)

, XLabel "Numeros"

, YLabel "Longitud de escalada"

]

[(x,longitudEscaladaPrimaAcotada x k) | x <- [2..n]]

Pares definidos por su MCD y su MCM

PorJosé A. Alonso 10-noviembre-201725-abril-2021

Definir las siguientes funciones

   pares  :: Integer -> Integer -> [(Integer,Integer)]
   nPares :: Integer -> Integer -> Integer

1 2	pares :: Integer -> Integer -> [(Integer,Integer)] nPares :: Integer -> Integer -> Integer

tales que

(pares a b) es la lista de los pares de números enteros positivos tales que su máximo común divisor es a y su mínimo común múltiplo es b. Por ejemplo,

     pares 3 3  == [(3,3)]
     pares 4 12 == [(4,12),(12,4)]
     pares 2 12 == [(2,12),(4,6),(6,4),(12,2)]
     pares 2 60 == [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]
     pares 2 7  == []
     pares 12 3 == []
     length (pares 3 (product [3,5..91]))  ==  8388608

pares 3 3 == [(3,3)]

pares 4 12 == [(4,12),(12,4)]

pares 2 12 == [(2,12),(4,6),(6,4),(12,2)]

pares 2 60 == [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]

pares 2 7 == []

pares 12 3 == []

length (pares 3 (product [3,5..91])) == 8388608

(nPares a b) es el número de pares de enteros positivos tales que su máximo común divisor es a y su mínimo común múltiplo es b. Por ejemplo,

     nPares 3 3   ==  1
     nPares 4 12  ==  2
     nPares 2 12  ==  4
     nPares 2 60  ==  8
     nPares 2 7   ==  0
     nPares 12 3  ==  0
     nPares 3 (product [3..3*10^4]) `mod` (10^12)  ==  477999992832
     length (show (nPares 3 (product [3..3*10^4])))  ==  977

nPares 3 3 == 1

nPares 4 12 == 2

nPares 2 12 == 4

nPares 2 60 == 8

nPares 2 7 == 0

nPares 12 3 == 0

nPares 3 (product [3..3*10^4]) `mod` (10^12) == 477999992832

length (show (nPares 3 (product [3..3*10^4]))) == 977

Soluciones

import Data.Numbers.Primes (primeFactors)
import Data.List (genericLength, group, nub, sort, subsequences)
import Test.QuickCheck

-- 1ª definición de pares
-- ======================

pares1 :: Integer -> Integer -> [(Integer,Integer)]
pares1 a b = [(x,y) | x <- [1..b]
                    , y <- [1..b]
                    , gcd x y == a
                    , lcm x y == b]

-- 2ª definición de pares
-- ======================

pares2 :: Integer -> Integer -> [(Integer,Integer)]
pares2 a b = [(x,y) | x <- [a,a+a..b]
                    , y <- [a,a+a..b]
                    , gcd x y == a
                    , lcm x y == b]

-- Comparacioń de eficiencia
--    λ> length (pares1 3 (product [3,5..11]))
--    16
--    (95.12 secs, 86,534,165,528 bytes)
--    λ> length (pares2 3 (product [3,5..11]))
--    16
--    (15.80 secs, 14,808,762,128 bytes)

-- 3ª definición de pares
-- ======================

pares3 :: Integer -> Integer -> [(Integer,Integer)]
pares3 a b = [(x,y) | x <- [a,a+a..b]
                    , c `rem` x == 0
                    , let y = c `div` x
                    , gcd x y == a
                    ]
  where c = a * b

-- Comparacioń de eficiencia
--    λ> length (pares2 3 (product [3,5..11]))
--    16
--    (15.80 secs, 14,808,762,128 bytes)
--    λ> length (pares3 3 (product [3,5..11]))
--    16
--    (0.01 secs, 878,104 bytes)

-- 4ª definición de pares
-- ======================

-- Para la cuarta definición de pares se observa la relación con los
-- factores primos
--    λ> [(primeFactors x, primeFactors y) | (x,y) <- pares1 2 12]
--    [([2],[2,2,3]),([2,2],[2,3]),([2,3],[2,2]),([2,2,3],[2])]
--    λ> [primeFactors x | (x,y) <- pares1 2 12]
--    [[2],[2,2],[2,3],[2,2,3]]
--    λ> [primeFactors x | (x,y) <- pares1 2 60]
--    [[2],[2,2],[2,3],[2,5],[2,2,3],[2,2,5],[2,3,5],[2,2,3,5]]
--    λ> [primeFactors x | (x,y) <- pares1 6 60]
--    [[2,3],[2,2,3],[2,3,5],[2,2,3,5]]
--    λ> [primeFactors x | (x,y) <- pares1 2 24]
--    [[2],[2,3],[2,2,2],[2,2,2,3]]
-- Se observa que cada pare se obtiene de uno de los subconjuntos de los
-- divisores primos de b/a. Por ejemplo,
--    λ> (a,b) = (2,24)
--    λ> b `div` a
--    12
--    λ> primeFactors it
--    [2,2,3]
--    λ> group it
--    [[2,2],[3]]
--    λ> subsequences it
--    [[],[[2,2]],[[3]],[[2,2],[3]]]
--    λ> map concat it
--    [[],[2,2],[3],[2,2,3]]
--    λ> map product it
--    [1,4,3,12]
--    λ> [(a * x, b `div` x) | x <- it]
--    [(2,24),(8,6),(6,8),(24,2)]
-- A partir de la observación se construye la siguiente definición

pares4 :: Integer -> Integer -> [(Integer,Integer)]
pares4 a b
  | b `mod` a /= 0 = []
  | otherwise =
    [(a * x, b `div` x)
    | x <- map (product . concat)
               ((subsequences . group . primeFactors) (b `div` a))]

-- Nota. La función pares4 calcula el mismo conjunto que las anteriores,
-- pero no necesariamente en el mismo orde. Por ejemplo,
--    λ> pares3 2 60 
--    [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]
--    λ> pares4 2 60 
--    [(2,60),(4,30),(6,20),(12,10),(10,12),(20,6),(30,4),(60,2)]
--    λ> pares3 2 60 == sort (pares4 2 60)
--    True

-- Comparacioń de eficiencia
--    λ> length (pares3 3 (product [3,5..17]))
--    64
--    (4.44 secs, 2,389,486,440 bytes)
--    λ> length (pares4 3 (product [3,5..17]))
--    64
--    (0.00 secs, 177,704 bytes)
  
-- Propiedades de equivalencia de pares
-- ====================================

prop_pares :: Integer -> Integer -> Property
prop_pares a b =
  a > 0 && b > 0 ==>
  all (== pares1 a b)
      [sort (f a b) | f <- [ pares2
                           , pares3
                           , pares4
                           ]]
  
prop_pares2 :: Integer -> Integer -> Property
prop_pares2 a b =
  a > 0 && b > 0 ==>
  all (== pares1 a (a * b))
      [sort (f a (a * b)) | f <- [ pares2
                                 , pares3
                                 , pares4
                                 ]]
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares
--    +++ OK, passed 100 tests.
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares
--    +++ OK, passed 100 tests.
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares2
--    +++ OK, passed 100 tests.
             
-- 1ª definición de nPares
-- =======================

nPares1 :: Integer -> Integer -> Integer
nPares1 a b = genericLength (pares4 a b)

-- 2ª definición de nPares
-- =======================

nPares2 :: Integer -> Integer -> Integer
nPares2 a b = 2^(length (nub (primeFactors (b `div` a))))

-- Comparación de eficiencia
--    λ> nPares1 3 (product [3,5..91])
--    8388608
--    (4.68 secs, 4,178,295,920 bytes)
--    λ> nPares2 3 (product [3,5..91])
--    8388608
--    (0.00 secs, 234,688 bytes)

-- Propiedad de equivalencia de nPares
-- ===================================

prop_nPares :: Integer -> Integer -> Property
prop_nPares a b =
  a > 0 && b > 0 ==>
  nPares1 a (a * b) == nPares2 a (a * b)
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_nPares
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

import Data.Numbers.Primes (primeFactors)

import Data.List (genericLength, group, nub, sort, subsequences)

import Test.QuickCheck

-- 1ª definición de pares

-- ======================

pares1 :: Integer -> Integer -> [(Integer,Integer)]

pares1 a b = [(x,y) | x <- [1..b]

, y <- [1..b]

, gcd x y == a

, lcm x y == b]

-- 2ª definición de pares

-- ======================

pares2 :: Integer -> Integer -> [(Integer,Integer)]

pares2 a b = [(x,y) | x <- [a,a+a..b]

, y <- [a,a+a..b]

, gcd x y == a

, lcm x y == b]

-- Comparacioń de eficiencia

-- λ> length (pares1 3 (product [3,5..11]))

-- 16

-- (95.12 secs, 86,534,165,528 bytes)

-- λ> length (pares2 3 (product [3,5..11]))

-- 16

-- (15.80 secs, 14,808,762,128 bytes)

-- 3ª definición de pares

-- ======================

pares3 :: Integer -> Integer -> [(Integer,Integer)]

pares3 a b = [(x,y) | x <- [a,a+a..b]

, c `rem` x == 0

, let y = c `div` x

, gcd x y == a

]

where c = a * b

-- Comparacioń de eficiencia

-- λ> length (pares2 3 (product [3,5..11]))

-- 16

-- (15.80 secs, 14,808,762,128 bytes)

-- λ> length (pares3 3 (product [3,5..11]))

-- 16

-- (0.01 secs, 878,104 bytes)

-- 4ª definición de pares

-- ======================

-- Para la cuarta definición de pares se observa la relación con los

-- factores primos

-- λ> [(primeFactors x, primeFactors y) | (x,y) <- pares1 2 12]

-- [([2],[2,2,3]),([2,2],[2,3]),([2,3],[2,2]),([2,2,3],[2])]

-- λ> [primeFactors x | (x,y) <- pares1 2 12]

-- [[2],[2,2],[2,3],[2,2,3]]

-- λ> [primeFactors x | (x,y) <- pares1 2 60]

-- [[2],[2,2],[2,3],[2,5],[2,2,3],[2,2,5],[2,3,5],[2,2,3,5]]

-- λ> [primeFactors x | (x,y) <- pares1 6 60]

-- [[2,3],[2,2,3],[2,3,5],[2,2,3,5]]

-- λ> [primeFactors x | (x,y) <- pares1 2 24]

-- [[2],[2,3],[2,2,2],[2,2,2,3]]

-- Se observa que cada pare se obtiene de uno de los subconjuntos de los

-- divisores primos de b/a. Por ejemplo,

-- λ> (a,b) = (2,24)

-- λ> b `div` a

-- 12

-- λ> primeFactors it

-- [2,2,3]

-- λ> group it

-- [[2,2],[3]]

-- λ> subsequences it

-- [[],[[2,2]],[[3]],[[2,2],[3]]]

-- λ> map concat it

-- [[],[2,2],[3],[2,2,3]]

-- λ> map product it

-- [1,4,3,12]

-- λ> [(a * x, b `div` x) | x <- it]

-- [(2,24),(8,6),(6,8),(24,2)]

-- A partir de la observación se construye la siguiente definición

pares4 :: Integer -> Integer -> [(Integer,Integer)]

pares4 a b

| b `mod` a /= 0 = []

| otherwise =

[(a * x, b `div` x)

| x <- map (product . concat)

((subsequences . group . primeFactors) (b `div` a))]

-- Nota. La función pares4 calcula el mismo conjunto que las anteriores,

-- pero no necesariamente en el mismo orde. Por ejemplo,

-- λ> pares3 2 60

-- [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]

-- λ> pares4 2 60

-- [(2,60),(4,30),(6,20),(12,10),(10,12),(20,6),(30,4),(60,2)]

-- λ> pares3 2 60 == sort (pares4 2 60)

-- True

-- Comparacioń de eficiencia

-- λ> length (pares3 3 (product [3,5..17]))

-- 64

-- (4.44 secs, 2,389,486,440 bytes)

-- λ> length (pares4 3 (product [3,5..17]))

-- 64

-- (0.00 secs, 177,704 bytes)

-- Propiedades de equivalencia de pares

-- ====================================

prop_pares :: Integer -> Integer -> Property

prop_pares a b =

a > 0 && b > 0 ==>

all (== pares1 a b)

[sort (f a b) | f <- [ pares2

, pares3

, pares4

]]

prop_pares2 :: Integer -> Integer -> Property

prop_pares2 a b =

a > 0 && b > 0 ==>

all (== pares1 a (a * b))

[sort (f a (a * b)) | f <- [ pares2

, pares3

, pares4

]]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares

-- +++ OK, passed 100 tests.

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares

-- +++ OK, passed 100 tests.

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares2

-- +++ OK, passed 100 tests.

-- 1ª definición de nPares

-- =======================

nPares1 :: Integer -> Integer -> Integer

nPares1 a b = genericLength (pares4 a b)

-- 2ª definición de nPares

-- =======================

nPares2 :: Integer -> Integer -> Integer

nPares2 a b = 2^(length (nub (primeFactors (b `div` a))))

-- Comparación de eficiencia

-- λ> nPares1 3 (product [3,5..91])

-- 8388608

-- (4.68 secs, 4,178,295,920 bytes)

-- λ> nPares2 3 (product [3,5..91])

-- 8388608

-- (0.00 secs, 234,688 bytes)

-- Propiedad de equivalencia de nPares

-- ===================================

prop_nPares :: Integer -> Integer -> Property

prop_nPares a b =

a > 0 && b > 0 ==>

nPares1 a (a * b) == nPares2 a (a * b)

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_nPares

-- +++ OK, passed 100 tests.

Cálculo de pi mediante la variante de Euler de la serie armónica

PorJosé A. Alonso 21-febrero-201726-marzo-2022

En el artículo El desarrollo más bello de Pi como suma infinita, Miguel Ángel Morales comenta el desarrollo de pi publicado por Leonhard Euler en su libro «Introductio in Analysis Infinitorum» (1748).

El desarrollo es el siguiente

y se obtiene a partir de la serie armónica

modificando sólo el signo de algunos términos según el siguiente criterio:

Dejamos un + cuando el denominador de la fracción sea un 2 o un primo de la forma 4m-1.
Cambiamos a – si el denominador de la fracción es un primo de la forma 4m+1.
Si el número es compuesto ponemos el signo que quede al multiplicar los signos correspondientes a cada factor.

Por ejemplo,

la de denominador 3 = 4×1-1 lleva un +,
la de denominador 5 = 4×1+1 lleva un -,
la de denominador 13 = 4×3+1 lleva un -,
la de denominador 6 = 2×3 lleva un + (porque los dos llevan un +),
la de denominador 10 = 2×5 lleva un – (porque el 2 lleva un + y el 5 lleva un -) y
la de denominador 50 = 5x5x2 lleva un + (un – por el primer 5, otro – por el segundo 5 y un + por el 2).

Definir las funciones

  aproximacionPi :: Int -> Double
  grafica        :: Int -> IO ()

1 2	aproximacionPi :: Int -> Double grafica :: Int -> IO ()

tales que

(aproximacionPi n) es la aproximación de pi obtenida sumando los n primeros términos de la serie de Euler. Por ejemplo.

     aproximacionPi 1        ==  1.0
     aproximacionPi 10       ==  2.3289682539682537
     aproximacionPi 100      ==  2.934318000847734
     aproximacionPi 1000     ==  3.0603246224585128
     aproximacionPi 10000    ==  3.1105295744825403
     aproximacionPi 100000   ==  3.134308801935256
     aproximacionPi 1000000  ==  3.1395057903490806

aproximacionPi 1 == 1.0

aproximacionPi 10 == 2.3289682539682537

aproximacionPi 100 == 2.934318000847734

aproximacionPi 1000 == 3.0603246224585128

aproximacionPi 10000 == 3.1105295744825403

aproximacionPi 100000 == 3.134308801935256

aproximacionPi 1000000 == 3.1395057903490806

(grafica n) dibuja la gráfica de las aproximaciones de pi usando k sumando donde k toma los valores de la lista [100,110..n]. Por ejemplo, al evaluar (grafica 4000) se obtiene

Nota: Este ejercicio ha sido propuesto por Paula Macías.

Soluciones

import Data.Numbers.Primes
import Graphics.Gnuplot.Simple

-- 1ª definición
-- =============

aproximacionPi :: Int -> Double
aproximacionPi n =
  sum [1 / fromIntegral (k * signo k) | k <- [1..n]] 

signoPrimo :: Int -> Int
signoPrimo 2 = 1
signoPrimo p | p `mod` 4 == 3 = 1
             | otherwise      = -1

signo :: Int -> Int
signo n | isPrime n = signoPrimo n
        | otherwise = product (map signoPrimo (primeFactors n))

-- 2ª definición
-- =============

aproximacionPi2 :: Int -> Double
aproximacionPi2 n = serieEuler !! (n-1)

serieEuler :: [Double]
serieEuler =
  scanl1 (+) [1 / fromIntegral (n * signo n) | n <- [1..]]

-- Definición de grafica
-- =====================

grafica :: Int -> IO ()
grafica n = 
    plotList [Key Nothing]
             [(k,aproximacionPi2 k) | k <- [100,110..n]]

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

-- 1ª definición

-- =============

aproximacionPi :: Int -> Double

aproximacionPi n =

sum [1 / fromIntegral (k * signo k) | k <- [1..n]]

signoPrimo :: Int -> Int

signoPrimo 2 = 1

signoPrimo p | p `mod` 4 == 3 = 1

| otherwise = -1

signo :: Int -> Int

signo n | isPrime n = signoPrimo n

| otherwise = product (map signoPrimo (primeFactors n))

-- 2ª definición

-- =============

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = serieEuler !! (n-1)

serieEuler :: [Double]

serieEuler =

scanl1 (+) [1 / fromIntegral (n * signo n) | n <- [1..]]

-- Definición de grafica

-- =====================

grafica :: Int -> IO ()

grafica n =

plotList [Key Nothing]

[(k,aproximacionPi2 k) | k <- [100,110..n]]

Números cuyos dígitos coinciden con los de sus factores primos

PorJosé A. Alonso 30-marzo-201625-marzo-2022

Un número n es especial si al unir los dígitos de sus factores primos, se obtienen exactamente los dígitos de n, aunque puede ser en otro orden. Por ejemplo, 1255 es especial, pues los factores primos de 1255 son 5 y 251.

Definir la función

   esEspecial :: Integer -> Bool

1	esEspecial :: Integer -> Bool

tal que (esEspecial n) se verifica si un número n es especial. Por ejemplo,

   esEspecial 1255 == True
   esEspecial 125  == False
   esEspecial 132  == False

esEspecial 1255 == True

esEspecial 125 == False

esEspecial 132 == False

Comprobar con QuickCheck que todo número primo es especial.

Calcular los 5 primeros números especiales que no son primos.

Soluciones

import Data.List (nub, sort)
import Data.Numbers.Primes (isPrime, primeFactors)
import Test.QuickCheck

esEspecial :: Integer -> Bool
esEspecial n = 
    sort (show n) == sort (concatMap show (nub (primeFactors n)))

-- La propiedad es
prop_primos:: Integer -> Property
prop_primos n = 
    isPrime (abs n) ==> esEspecial (abs n)

-- La comprobación es
--    ghci> quickCheck prop_primos
--    +++ OK, passed 100 tests.

-- El cálculo es
--    ghci> take 5 [n | n <- [2..], esEspecial n, not (isPrime n)]
--    [735,1255,3792,7236,11913]

import Data.List (nub, sort)

import Data.Numbers.Primes (isPrime, primeFactors)

import Test.QuickCheck

esEspecial :: Integer -> Bool

esEspecial n =

sort (show n) == sort (concatMap show (nub (primeFactors n)))

-- La propiedad es

prop_primos:: Integer -> Property

prop_primos n =

isPrime (abs n) ==> esEspecial (abs n)

-- La comprobación es

-- ghci> quickCheck prop_primos

-- +++ OK, passed 100 tests.

-- El cálculo es

-- ghci> take 5 [n | n <- [2..], esEspecial n, not (isPrime n)]

-- [735,1255,3792,7236,11913]

Paridad del número de divisores

PorJosé A. Alonso 10-marzo-201617-marzo-2016

Definir la función

   nDivisoresPar :: Integer -> Bool

1	nDivisoresPar :: Integer -> Bool

tal que (nDivisoresPar n) se verifica si n tiene un número par de divisores. Por ejemplo,

   nDivisoresPar 12                     ==  True
   nDivisoresPar 16                     ==  False
   nDivisoresPar (product [1..2*10^4])  ==  True

nDivisoresPar 12 == True

nDivisoresPar 16 == False

nDivisoresPar (product [1..2*10^4]) == True

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primeFactors)

-- 1ª definición
-- =============

nDivisoresPar1 :: Integer -> Bool
nDivisoresPar1 = even . length . divisores

-- (divisores n) es la lista de los divisores de n. Por ejemplo,
--    divisores 12  ==  [1,2,3,4,6,12]
--    divisores 16  ==  [1,2,4,8,16]
divisores :: Integer -> [Integer]
divisores n =
    [x | x <- [1..n], n `mod` x == 0]

-- 2ª definición
-- =============

nDivisoresPar2 :: Integer -> Bool
nDivisoresPar2 n =
    any odd (map length (group (primeFactors n)))

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

--    λ> nDivisoresPar1 (10^7)
--    True
--    (14.18 secs, 2,078,736,992 bytes)
--    λ> nDivisoresPar2 (10^7)
--    True
--    (0.01 secs, 0 bytes)
--
--    λ> nDivisoresPar2 (product [1..2*10^4])
--    True
--    (2.30 secs, 1,003,013,112 bytes)

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

-- 1ª definición

-- =============

nDivisoresPar1 :: Integer -> Bool

nDivisoresPar1 = even . length . divisores

-- (divisores n) es la lista de los divisores de n. Por ejemplo,

-- divisores 12 == [1,2,3,4,6,12]

-- divisores 16 == [1,2,4,8,16]

divisores :: Integer -> [Integer]

divisores n =

[x | x <- [1..n], n `mod` x == 0]

-- 2ª definición

-- =============

nDivisoresPar2 :: Integer -> Bool

nDivisoresPar2 n =

any odd (map length (group (primeFactors n)))

-- Comparación de eficiencia

-- =========================

-- λ> nDivisoresPar1 (10^7)

-- True

-- (14.18 secs, 2,078,736,992 bytes)

-- λ> nDivisoresPar2 (10^7)

-- True

-- (0.01 secs, 0 bytes)

-- λ> nDivisoresPar2 (product [1..2*10^4])

-- True

-- (2.30 secs, 1,003,013,112 bytes)

Solución en Maxima

/* concat(xs) es la lista obtenida concatenado los elementos de xss. Por
   ejemplo,
      concat([[1,3],[2,4],[5,7]])  ==  [1, 3, 2, 4, 5, 7]             */
concat (xs) := block ([r:[]],
  for x in reverse (xs) do r : append (x,r),
  r)$

nDivisoresPares (n) := block(
  [xs : concat (map (rest, ifactors (n))), p : false],
  for x in xs do p: oddp (x) or p,
  p)$

/* concat(xs) es la lista obtenida concatenado los elementos de xss. Por

ejemplo,

concat([[1,3],[2,4],[5,7]]) == [1, 3, 2, 4, 5, 7] */

concat (xs) := block ([r:[]],

for x in reverse (xs) do r : append (x,r),

r)$

nDivisoresPares (n) := block(

[xs : concat (map (rest, ifactors (n))), p : false],

for x in xs do p: oddp (x) or p,

p)$

Dígitos en la factorización

PorJosé A. Alonso 19-febrero-20162-marzo-2016

El enunciado del problema 652 de Números y algo más es el siguiente

Si factorizamos los factoriales de un número en función de sus divisores primos y sus potencias, ¿Cuál es el menor número n tal que entre los factores primos y los exponentes de estos, n! contiene los dígitos del cero al nueve? Por ejemplo

6! = 2⁴x3²x5¹, le faltan los dígitos 0,6,7,8 y 9

12! = 2¹⁰x3⁵x5²x7¹x11¹, le faltan los dígitos 4,6,8 y 9

Definir la función

   digitosDeFactorizacion :: Integer -> [Integer]

1	digitosDeFactorizacion :: Integer -> [Integer]

tal que (digitosDeFactorizacion n) es el conjunto de los dígitos que aparecen en la factorización de n. Por ejemplo,

   digitosDeFactorizacion (factorial 6)   ==  [1,2,3,4,5]
   digitosDeFactorizacion (factorial 12)  ==  [0,1,2,3,5,7]

1 2	digitosDeFactorizacion (factorial 6) == [1,2,3,4,5] digitosDeFactorizacion (factorial 12) == [0,1,2,3,5,7]

Usando la función anterior, calcular la solución del problema.

Comprobar con QuickCheck que si n es mayor que 100, entonces

   digitosDeFactorizacion (factorial n) == [0..9]

1	digitosDeFactorizacion (factorial n) == [0..9]

Soluciones

import Data.List (genericLength, group, nub, sort)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck

-- 1ª definición
-- =============

digitosDeFactorizacion1 :: Integer -> [Integer]
digitosDeFactorizacion1 n =
   sort (nub (concat [digitos x | x <- numerosDeFactorizacion n]))

-- (digitos n) es la lista de los digitos del número n. Por ejemplo, 
--    digitos 320274  ==  [3,2,0,2,7,4]
digitos :: Integer -> [Integer]
digitos n = [read [x] | x <- show n]

-- (numerosDeFactorizacion n) es el conjunto de los números en la
-- factorización de n. Por ejemplo,
--    numerosDeFactorizacion 60  ==  [1,2,3,5]
numerosDeFactorizacion :: Integer -> [Integer]
numerosDeFactorizacion n = 
   sort (nub (aux (factorizacion n)))
   where aux [] = []
         aux ((x,y):zs) = x : y : aux zs

-- (factorización n) es la factorización de n. Por ejemplo,
--    factorizacion 300  ==  [(2,2),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n = 
    [(head xs, genericLength xs) | xs <- group (factorizacion' n)]

-- (factorizacion' n) es la lista de todos los factores primos de n; es
-- decir, es una lista de números primos cuyo producto es n. Por ejemplo,
--    factorizacion 300  ==  [2,2,3,5,5]
factorizacion' :: Integer -> [Integer]
factorizacion' n | n == 1    = []
                 | otherwise = x : factorizacion' (div n x)
                 where x = menorFactor n

-- (menorFactor n) es el menor factor primo de n. Por ejemplo,
--    menorFactor 15  ==  3
menorFactor :: Integer -> Integer
menorFactor n = head [x | x <- [2..], rem n x == 0]

-- (factorial n) es el factorial de n. Por ejemplo,
--    factorial 5  ==  120
factorial :: Integer -> Integer
factorial n = product [1..n]

-- 2ª definición
-- =============

digitosDeFactorizacion2 :: Integer -> [Integer]
digitosDeFactorizacion2 n =
   sort (nub (concat [digitos x | x <- numerosDeFactorizacion2 n]))

-- (numerosDeFactorizacion2 n) es el conjunto de los números en la
-- factorización de n. Por ejemplo,
--    numerosDeFactorizacion2 60  ==  [1,2,3,5]
numerosDeFactorizacion2 :: Integer -> [Integer]
numerosDeFactorizacion2 n = 
    sort (nub (aux (factorizacion2 n)))
    where aux xs = concat [[a,b] | (a,b) <- xs]

-- (factorización2 n) es la factorización de n. Por ejemplo,
--    factorizacion2 300  ==  [(2,2),(3,1),(5,2)]
factorizacion2 :: Integer -> [(Integer,Integer)]
factorizacion2 n =
    [(head xs, genericLength xs) | xs <- group (primeFactors n)]

-- 3ª definición
-- =============

digitosDeFactorizacion3 :: Integer -> [Integer]
digitosDeFactorizacion3 n =
    sort (nub (concat [digitos x | x <- aux (group (primeFactors n))]))
    where aux  []            = []
          aux (ys@(y:_):xss) = y : genericLength ys : aux xss

-- Definición
-- ==========

digitosDeFactorizacion :: Integer -> [Integer]
digitosDeFactorizacion = digitosDeFactorizacion3

-- Solución
-- ========

-- Para calcular la solución, se define la constante
solucion :: Integer
solucion = 
    head [n | n <- [1..], digitosDeFactorizacion (factorial n) == [0..9]]

-- El cálculo de la solución es
--    ghci> solucion2
--    49

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

-- La propiedad es
prop_completa :: Integer -> Bool
prop_completa n =
    digitosDeFactorizacion (factorial n1) == [0..9]
    where n1 = 101 + abs n

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

100

101

102

103

104

105

106

107

108

109

import Data.List (genericLength, group, nub, sort)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

-- 1ª definición

-- =============

digitosDeFactorizacion1 :: Integer -> [Integer]

digitosDeFactorizacion1 n =

sort (nub (concat [digitos x | x <- numerosDeFactorizacion n]))

-- (digitos n) es la lista de los digitos del número n. Por ejemplo,

-- digitos 320274 == [3,2,0,2,7,4]

digitos :: Integer -> [Integer]

digitos n = [read [x] | x <- show n]

-- (numerosDeFactorizacion n) es el conjunto de los números en la

-- factorización de n. Por ejemplo,

-- numerosDeFactorizacion 60 == [1,2,3,5]

numerosDeFactorizacion :: Integer -> [Integer]

numerosDeFactorizacion n =

sort (nub (aux (factorizacion n)))

where aux [] = []

aux ((x,y):zs) = x : y : aux zs

-- (factorización n) es la factorización de n. Por ejemplo,

-- factorizacion 300 == [(2,2),(3,1),(5,2)]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion n =

[(head xs, genericLength xs) | xs <- group (factorizacion' n)]

-- (factorizacion' n) es la lista de todos los factores primos de n; es

-- decir, es una lista de números primos cuyo producto es n. Por ejemplo,

-- factorizacion 300 == [2,2,3,5,5]

factorizacion' :: Integer -> [Integer]

factorizacion' n | n == 1 = []

| otherwise = x : factorizacion' (div n x)

where x = menorFactor n

-- (menorFactor n) es el menor factor primo de n. Por ejemplo,

-- menorFactor 15 == 3

menorFactor :: Integer -> Integer

menorFactor n = head [x | x <- [2..], rem n x == 0]

-- (factorial n) es el factorial de n. Por ejemplo,

-- factorial 5 == 120

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª definición

-- =============

digitosDeFactorizacion2 :: Integer -> [Integer]

digitosDeFactorizacion2 n =

sort (nub (concat [digitos x | x <- numerosDeFactorizacion2 n]))

-- (numerosDeFactorizacion2 n) es el conjunto de los números en la

-- factorización de n. Por ejemplo,

-- numerosDeFactorizacion2 60 == [1,2,3,5]

numerosDeFactorizacion2 :: Integer -> [Integer]

numerosDeFactorizacion2 n =

sort (nub (aux (factorizacion2 n)))

where aux xs = concat [[a,b] | (a,b) <- xs]

-- (factorización2 n) es la factorización de n. Por ejemplo,

-- factorizacion2 300 == [(2,2),(3,1),(5,2)]

factorizacion2 :: Integer -> [(Integer,Integer)]

factorizacion2 n =

[(head xs, genericLength xs) | xs <- group (primeFactors n)]

-- 3ª definición

-- =============

digitosDeFactorizacion3 :: Integer -> [Integer]

digitosDeFactorizacion3 n =

sort (nub (concat [digitos x | x <- aux (group (primeFactors n))]))

where aux [] = []

aux (ys@(y:_):xss) = y : genericLength ys : aux xss

-- Definición

-- ==========

digitosDeFactorizacion :: Integer -> [Integer]

digitosDeFactorizacion = digitosDeFactorizacion3

-- Solución

-- ========

-- Para calcular la solución, se define la constante

solucion :: Integer

solucion =

head [n | n <- [1..], digitosDeFactorizacion (factorial n) == [0..9]]

-- El cálculo de la solución es

-- ghci> solucion2

-- 49

-- Propiedad

-- =========

-- La propiedad es

prop_completa :: Integer -> Bool

prop_completa n =

digitosDeFactorizacion (factorial n1) == [0..9]

where n1 = 101 + abs n

-- La comprobación es

-- λ> quickCheck prop_completa

-- +++ OK, passed 100 tests.

La solución en Maxima

/* digitos(n) es la lista de los digitos del número n. Por ejemplo, 
      digitos (320274) == [3,2,0,2,7,4]
*/      
digitos (n) := charlist (string (n))$

digitosDeFactorizacion (n) :=
   unique (apply (append, map (digitos, apply (append, ifactors (n))))) $

solucion () := block ([n:1],
  while length (digitosDeFactorizacion (n!)) < 10 do
     n : n+1,
  n)$   

/*
   (%i5) digitosDeFactorizacion (6!);
   (%o5) [1, 2, 3, 4, 5]

   (%i6) solucion ();
   (%o6) 49
*/

/* digitos(n) es la lista de los digitos del número n. Por ejemplo,

digitos (320274) == [3,2,0,2,7,4]

digitos (n) := charlist (string (n))$

digitosDeFactorizacion (n) :=

unique (apply (append, map (digitos, apply (append, ifactors (n))))) $

solucion () := block ([n:1],

while length (digitosDeFactorizacion (n!)) < 10 do

n : n+1,

n)$

(%i5) digitosDeFactorizacion (6!);

(%o5) [1, 2, 3, 4, 5]

(%i6) solucion ();

(%o6) 49

Puntos visibles en la cuadrícula de un plano

PorJosé A. Alonso 13-enero-20161-mayo-2016

La cuadrícula entera de lado n, Cₙ, es el conjunto de los puntos (x,y) donde x e y son números enteros tales que 1 ≤ x, y ≤ n.

Un punto (x,y) de Cₙ es visible desde el origen si el máximo común divisor de x e y es 1. Por ejemplo, el punto (4,6) no es visible porque está ocultado por el (2,3); en cambio, el (2,3) sí es visible.

El conjunto de los puntos visibles en la cuadrícula entera de lado 6 son (1,1), (1,2), (1,3), (1,4), (1,5), (1,6), (2,1), (2,3), (2,5), (3,1), (3,2), (3,4), (3,5), (4,1), (4,3), (4,5), (5,1), (5,2), (5,3), (5,4), (5,6), (6,1) y (6,5).

Definir la función

   nVisibles :: Integer -> Integer

1	nVisibles :: Integer -> Integer

tal que (nVisibles n) es el número de los puntos visibles en la cuadrícula de lado n.Por ejemplo,

   nVisibles 6       ==  23
   nVisibles 10      ==  63
   nVisibles 100     ==  6087
   nVisibles 1000    ==  608383
   nVisibles 10000   ==  60794971
   nVisibles 100000  ==  6079301507

nVisibles 6 == 23

nVisibles 10 == 63

nVisibles 100 == 6087

nVisibles 1000 == 608383

nVisibles 10000 == 60794971

nVisibles 100000 == 6079301507

Soluciones

import Data.List (genericLength, group, sort)
import Data.Numbers.Primes (primeFactors)

type Punto = (Integer,Integer)

-- (visible p) se verifica si el punto p es visible. Por ejemplo,
--    visible (4,5)  ==  True
--    visible (4,6)  ==  False
visible :: Punto -> Bool
visible (x,y) = gcd x y == 1

-- (visibles n) es el conjunto de los puntos visibles en la cuadrícula
-- de lado 1. Por ejemplo,
-- cuadrado 1 <= x, y <= n. Por ejemplo,
--    λ> sort (visibles 6)
--    [(1,1),(1,2),(1,3),(1,4),(1,5),(1,6),(2,1),(2,3),(2,5),(3,1),(3,2),(3,4),
--     (3,5),(4,1),(4,3),(4,5),(5,1),(5,2),(5,3),(5,4),(5,6),(6,1),(6,5)]

-- 1ª definición de visibles
visibles1 :: Integer -> [Punto]
visibles1 n = [(x,y) | x <- [1..n], y <- [1..n], visible (x,y)]

-- 2ª definición de visibles
visibles2 :: Integer -> [Punto]
visibles2 n = (1,1) : ps ++ [(y,x) | (x,y) <- ps] 
    where ps = [(x,y) | x <- [1..n], y <- [x+1..n], visible (x,y)]

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

--    λ> length (visibles1 1000)
--    608383
--    (3.40 secs, 758,346,208 bytes)
--    λ> length (visibles2 1000)
--    608383
--    (2.22 secs, 464,276,856 bytes)

-- Usaremos la 2ª definición de visibles
visibles :: Integer -> [Punto]
visibles = visibles2

-- 1ª definición de nVisibles
nVisibles1 :: Integer -> Integer
nVisibles1 = genericLength . visibles


-- 2ª definición de nVisibles
nVisibles2 :: Integer -> Integer
nVisibles2 n = 
    1 + 2 * genericLength [(x,y) | x <- [1..n], 
                                   y <- [x+1..n], 
                                   visible (x,y)]

-- 3ª definición de nVisibles
nVisibles3 :: Integer -> Integer
nVisibles3 1 = 1
nVisibles3 n = nVisibles3 (n-1) + (2 * phi n)

-- La función φ de Euler (del ejercicio anterior).
phi :: Integer -> Integer
phi n = product [(p-1)*p^(e-1) | (p,e) <- factorizacion n] 
    
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
    [(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- 4ª definición de nVisibles
nVisibles4 :: Integer -> Integer
nVisibles4 n = 2 * sum [phi x | x <- [1..n]] - 1

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

--    λ> nVisibles1 1000
--    608383
--    (2.69 secs, 521,599,512 bytes)
--    λ> nVisibles2 1000
--    608383
--    λ> nVisibles3 1000
--    608383
--    (0.05 secs, 0 bytes)
--    λ> nVisibles4 1000
--    608383
--    (0.04 secs, 0 bytes)

import Data.List (genericLength, group, sort)

import Data.Numbers.Primes (primeFactors)

type Punto = (Integer,Integer)

-- (visible p) se verifica si el punto p es visible. Por ejemplo,

-- visible (4,5) == True

-- visible (4,6) == False

visible :: Punto -> Bool

visible (x,y) = gcd x y == 1

-- (visibles n) es el conjunto de los puntos visibles en la cuadrícula

-- de lado 1. Por ejemplo,

-- cuadrado 1 <= x, y <= n. Por ejemplo,

-- λ> sort (visibles 6)

-- [(1,1),(1,2),(1,3),(1,4),(1,5),(1,6),(2,1),(2,3),(2,5),(3,1),(3,2),(3,4),

-- (3,5),(4,1),(4,3),(4,5),(5,1),(5,2),(5,3),(5,4),(5,6),(6,1),(6,5)]

-- 1ª definición de visibles

visibles1 :: Integer -> [Punto]

visibles1 n = [(x,y) | x <- [1..n], y <- [1..n], visible (x,y)]

-- 2ª definición de visibles

visibles2 :: Integer -> [Punto]

visibles2 n = (1,1) : ps ++ [(y,x) | (x,y) <- ps]

where ps = [(x,y) | x <- [1..n], y <- [x+1..n], visible (x,y)]

-- Comparación de eficiencia

-- =========================

-- λ> length (visibles1 1000)

-- 608383

-- (3.40 secs, 758,346,208 bytes)

-- λ> length (visibles2 1000)

-- 608383

-- (2.22 secs, 464,276,856 bytes)

-- Usaremos la 2ª definición de visibles

visibles :: Integer -> [Punto]

visibles = visibles2

-- 1ª definición de nVisibles

nVisibles1 :: Integer -> Integer

nVisibles1 = genericLength . visibles

-- 2ª definición de nVisibles

nVisibles2 :: Integer -> Integer

nVisibles2 n =

1 + 2 * genericLength [(x,y) | x <- [1..n],

y <- [x+1..n],

visible (x,y)]

-- 3ª definición de nVisibles

nVisibles3 :: Integer -> Integer

nVisibles3 1 = 1

nVisibles3 n = nVisibles3 (n-1) + (2 * phi n)

-- La función φ de Euler (del ejercicio anterior).

phi :: Integer -> Integer

phi n = product [(p-1)*p^(e-1) | (p,e) <- factorizacion n]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion n =

[(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- 4ª definición de nVisibles

nVisibles4 :: Integer -> Integer

nVisibles4 n = 2 * sum [phi x | x <- [1..n]] - 1

-- Comparación de eficiencia

-- =========================

-- λ> nVisibles1 1000

-- 608383

-- (2.69 secs, 521,599,512 bytes)

-- λ> nVisibles2 1000

-- 608383

-- λ> nVisibles3 1000

-- 608383

-- (0.05 secs, 0 bytes)

-- λ> nVisibles4 1000

-- 608383

-- (0.04 secs, 0 bytes)

Referencias

N.J.A. Sloane, Sucesión A018805 en OEIS.

Parte libre de cuadrados y parte cuadrada de un número

PorJosé A. Alonso 11-enero-201618-enero-2016

La parte libre de cuadrados de un número n es el producto de todos sus divisores primos con exponente impar en la factorización prima de n. Por ejemplo, la parte libre de cuadrados de 360 es 10 ya que 360 = 2³3²5 y 2.5 = 10; además, 360 = 10.6²

La parte cuadrada de un número n es el mayor número cuadrado que divide a n. Por ejemplo, la parte cuadrada de 360 es 6.

Definir las funciones

   parteLibre    :: Integer -> Integer
   parteCuadrada :: Integer -> Integer

1 2	parteLibre :: Integer -> Integer parteCuadrada :: Integer -> Integer

tales que

(parteLibre x) es la parte libre de x. Por ejemplo,

     parteLibre 360                 ==  10
     parteLibre 1800                ==  2
     [parteLibre n | n <- [1..14]]  ==  [1,2,3,1,5,6,7,2,1,10,11,3,13,14]

parteLibre 360 == 10

parteLibre 1800 == 2

[parteLibre n | n <- [1..14]] == [1,2,3,1,5,6,7,2,1,10,11,3,13,14]

(parteCuadrada x) es la parte cuadrada de x. Por ejemplo,

     parteCuadrada 360                 ==  36
     parteCuadrada 1800                ==  900
     [parteCuadrada n | n <- [1..14]]  ==  [1,1,1,4,1,1,1,4,9,1,1,4,1,1]

parteCuadrada 360 == 36

parteCuadrada 1800 == 900

[parteCuadrada n | n <- [1..14]] == [1,1,1,4,1,1,1,4,9,1,1,4,1,1]

Soluciones

import Data.List (group, genericLength)
import Data.Numbers.Primes (primeFactors)

parteLibre :: Integer -> Integer
parteLibre n = product [p | (p,e) <- factorizacion n, odd e]

-- (factorizacion n) es la factorización de n. Por ejemplo, 
--    factorizacion  600  ==  [(2,3),(3,1),(5,2)]
--    factorizacion 1400  ==  [(2,3),(5,2),(7,1)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
    [(head xs,genericLength xs) | xs <- group (primeFactors n)]

parteCuadrada :: Integer -> Integer
parteCuadrada n = n `div` parteLibre n

import Data.List (group, genericLength)

import Data.Numbers.Primes (primeFactors)

parteLibre :: Integer -> Integer

parteLibre n = product [p | (p,e) <- factorizacion n, odd e]

-- (factorizacion n) es la factorización de n. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

-- factorizacion 1400 == [(2,3),(5,2),(7,1)]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion n =

[(head xs,genericLength xs) | xs <- group (primeFactors n)]

parteCuadrada :: Integer -> Integer

parteCuadrada n = n `div` parteLibre n

Referencias

N.J.A. Sloane, Sequence A008833 en OEIS.
N.J.A. Sloane, Sequence A007913 en OEIS.
E.W. Weisstein, Square part en MathWorld.
E.W. Weisstein, Squarefree part en MathWorld.

2016 es un número práctico

PorJosé A. Alonso 4-enero-20161-mayo-2016

Un entero positivo n es un número práctico si todos los enteros positivos menores que él se pueden expresar como suma de distintos divisores de n. Por ejemplo, el 12 es un número práctico, ya que todos los enteros positivos menores que 12 se pueden expresar como suma de divisores de 12 (1, 2, 3, 4 y 6) sin usar ningún divisor más de una vez en cada suma:

    1 = 1
    2 = 2
    3 = 3
    4 = 4
    5 = 2 + 3
    6 = 6
    7 = 1 + 6
    8 = 2 + 6
    9 = 3 + 6
   10 = 4 + 6
   11 = 1 + 4 + 6

1 = 1

2 = 2

3 = 3

4 = 4

5 = 2 + 3

6 = 6

7 = 1 + 6

8 = 2 + 6

9 = 3 + 6

10 = 4 + 6

11 = 1 + 4 + 6

En cambio, 14 no es un número práctico ya que 6 no se puede escribir como suma, con sumandos distintos, de divisores de 14.

Definir la función

   esPractico :: Integer -> Bool

1	esPractico :: Integer -> Bool

tal que (esPractico n) se verifica si n es un número práctico. Por ejemplo,

   esPractico 12                                      ==  True
   esPractico 14                                      ==  False
   esPractico 2016                                    ==  True
   esPractico 42535295865117307932921825928971026432  ==  True

esPractico 12 == True

esPractico 14 == False

esPractico 2016 == True

esPractico 42535295865117307932921825928971026432 == True

Soluciones

import Data.List (genericLength, group, nub, sort, subsequences)
import Data.Numbers.Primes (primeFactors)
import Graphics.Gnuplot.Simple

-- 1ª definición
-- =============

esPractico1 :: Integer -> Bool
esPractico1 n =
    takeWhile (<n) (sumas (divisores n)) == [0..n-1]

-- (divisores n) es la lista de los divisores de n. Por ejemplo,
--    divisores 12  ==  [1,2,3,4,6]
--    divisores 14  ==  [1,2,7]
divisores :: Integer -> [Integer]
divisores n = [k | k <- [1..n-1], n `mod` k == 0]

-- (sumas xs) es la lista ordenada de números que se pueden obtener como
-- sumas de elementos de xs sin usar ningún elemento más de una vez en
-- cada suma. Por ejemplo,  
--    sumas [1,2,3]  ==  [0,1,2,3,4,5,6]
--    sumas [1,2,7]  ==  [0,1,2,3,7,8,9,10]
sumas :: [Integer] -> [Integer]
sumas xs = sort (nub (map sum (subsequences xs)))

-- 2ª definición
-- =============

esPractico2 :: Integer -> Bool
esPractico2 n = all (esSumable (divisores n)) [1..n-1]

-- (esSumable xs n) se verifica si n se puede escribir como una suma de
-- elementos distintos de la lista creciente xs. Por ejemplo,
--    esSumable [1,2,7] 8  ==  True
--    esSumable [1,2,7] 6  ==  False
--    esSumable [1,2,7] 4  ==  False
--    esSumable [1,2,7] 2  ==  True
--    esSumable [1,2,7] 0  ==  True
esSumable :: [Integer] -> Integer -> Bool
esSumable _ 0  = True
esSumable [] _ = False
esSumable (x:xs) n = x <= n && (esSumable xs (n-x) || esSumable xs n)

-- 3ª definición
-- =============

-- Usando la caracterización de Stewart y Sierpiński: un entero n >= 2
-- es práctico syss para su factorización prima
--    n = p(1)^e(1) * p(2)*e(2) *...* p(k)^e(k)
-- se cumple que p(1) = 2 y, para cada i de 2 a k se cumple que
--                         1+e(j) 
--                i-1  p(j)       - 1
--    p(i) <= 1 +  ∏  ----------------
--                j=1     p(j) - 1

esPractico3 :: Integer -> Bool
esPractico3 1 = True
esPractico3 n = 
    x == 2 &&
    and [p <= 1 + c | (p,c) <- zip bases cotas]
    where xss       = factorizacion n
          (x:bases) = map fst xss
          cotas     = scanl1 (*) [(p^(1+e)-1) `div` (p-1) | (p,e) <- xss]

-- (factorizacion n) es la factorización de n. Por ejemplo, 
--    factorizacion  600  ==  [(2,3),(3,1),(5,2)]
--    factorizacion 1400  ==  [(2,3),(5,2),(7,1)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
    [(head xs,genericLength xs) | xs <- group (primeFactors n)]

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

--    λ> length [n | n <- [1..400], esPractico1 n]
--    92
--    (40.21 secs, 8,378,539,464 bytes)
--    λ> length [n | n <- [1..400], esPractico2 n]
--    92
--    (8.29 secs, 1,109,669,760 bytes)
--    λ> length [n | n <- [1..400], esPractico3 n]
--    92
--    (0.02 secs, 0 bytes)

import Data.List (genericLength, group, nub, sort, subsequences)

import Data.Numbers.Primes (primeFactors)

import Graphics.Gnuplot.Simple

-- 1ª definición

-- =============

esPractico1 :: Integer -> Bool

esPractico1 n =

takeWhile (<n) (sumas (divisores n)) == [0..n-1]

-- (divisores n) es la lista de los divisores de n. Por ejemplo,

-- divisores 12 == [1,2,3,4,6]

-- divisores 14 == [1,2,7]

divisores :: Integer -> [Integer]

divisores n = [k | k <- [1..n-1], n `mod` k == 0]

-- (sumas xs) es la lista ordenada de números que se pueden obtener como

-- sumas de elementos de xs sin usar ningún elemento más de una vez en

-- cada suma. Por ejemplo,

-- sumas [1,2,3] == [0,1,2,3,4,5,6]

-- sumas [1,2,7] == [0,1,2,3,7,8,9,10]

sumas :: [Integer] -> [Integer]

sumas xs = sort (nub (map sum (subsequences xs)))

-- 2ª definición

-- =============

esPractico2 :: Integer -> Bool

esPractico2 n = all (esSumable (divisores n)) [1..n-1]

-- (esSumable xs n) se verifica si n se puede escribir como una suma de

-- elementos distintos de la lista creciente xs. Por ejemplo,

-- esSumable [1,2,7] 8 == True

-- esSumable [1,2,7] 6 == False

-- esSumable [1,2,7] 4 == False

-- esSumable [1,2,7] 2 == True

-- esSumable [1,2,7] 0 == True

esSumable :: [Integer] -> Integer -> Bool

esSumable _ 0 = True

esSumable [] _ = False

esSumable (x:xs) n = x <= n && (esSumable xs (n-x) || esSumable xs n)

-- 3ª definición

-- =============

-- Usando la caracterización de Stewart y Sierpiński: un entero n >= 2

-- es práctico syss para su factorización prima

-- n = p(1)^e(1) * p(2)*e(2) *...* p(k)^e(k)

-- se cumple que p(1) = 2 y, para cada i de 2 a k se cumple que

-- 1+e(j)

-- i-1 p(j) - 1

-- p(i) <= 1 + ∏ ----------------

-- j=1 p(j) - 1

esPractico3 :: Integer -> Bool

esPractico3 1 = True

esPractico3 n =

x == 2 &&

and [p <= 1 + c | (p,c) <- zip bases cotas]

where xss = factorizacion n

(x:bases) = map fst xss

cotas = scanl1 (*) [(p^(1+e)-1) `div` (p-1) | (p,e) <- xss]

-- (factorizacion n) es la factorización de n. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

-- factorizacion 1400 == [(2,3),(5,2),(7,1)]

factorizacion :: Integer -> [(Integer,Integer)]

factorizacion n =

[(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- Comparación de eficiencia

-- =========================

-- λ> length [n | n <- [1..400], esPractico1 n]

-- 92

-- (40.21 secs, 8,378,539,464 bytes)

-- λ> length [n | n <- [1..400], esPractico2 n]

-- 92

-- (8.29 secs, 1,109,669,760 bytes)

-- λ> length [n | n <- [1..400], esPractico3 n]

-- 92

-- (0.02 secs, 0 bytes)

Referencias

Basado en el artículo de Gaussianos Feliz Navidad y Feliz Año (número práctico) 2016.

Otras referencias

M. Margenstern (1991), Les nombres pratiques: théorie, observations et conjectures.
G. Melfi, Practical numbers.
G. Melfi (2008), A survey on practical numbers.
OEIS, Sucesión A005153.
PlanetMath.org, Practical number.
E.W. Weisstein, Practical number.
A. Weingartner (2015), Practical numbers and the distribution of divisors.
Wikipedia, Practical number.

Los números de Smith

PorJosé A. Alonso 4-diciembre-201511-diciembre-2015

Un número de Smith es un número natural compuesto que cumple que la suma de sus dígitos es igual a la suma de los dígitos de todos sus factores primos (si tenemos algún factor primo repetido lo sumamos tantas veces como aparezca). Por ejemplo, el 22 es un número de Smith ya que

    22 = 2*11 y
   2+2 = 2+(1+1)

1 2	22 = 2*11 y 2+2 = 2+(1+1)

y el 4937775 también lo es ya que

   4937775       = 3*5*5*65837 y 
   4+9+3+7+7+7+5 = 3+5+5+(6+5+8+3+7)

1 2	4937775 = 355*65837 y 4+9+3+7+7+7+5 = 3+5+5+(6+5+8+3+7)

Definir las funciones

   esSmith :: Integer -> Bool
   smith :: [Integer]

1 2	esSmith :: Integer -> Bool smith :: [Integer]

tales que

(esSmith x) se verifica si x es un número de Smith. Por ejemplo,

     esSmith 22          ==  True
     esSmith 29          ==  False
     esSmith 2015        ==  False
     esSmith 4937775     ==  True
     esSmith 4567597056  ==  True

esSmith 22 == True

esSmith 29 == False

esSmith 2015 == False

esSmith 4937775 == True

esSmith 4567597056 == True

smith es la lista cuyos elementos son los números de Smith. Por ejemplo,

     λ> take 17 smith
     [4,22,27,58,85,94,121,166,202,265,274,319,346,355,378,382,391]
     λ> smith !! 2000
     62158

λ> take 17 smith

[4,22,27,58,85,94,121,166,202,265,274,319,346,355,378,382,391]

λ> smith !! 2000

62158

Soluciones

import Data.Numbers.Primes

esSmith :: Integer -> Bool
esSmith x = 
    not (isPrime x) && 
    sumaDigitos x == sum (map sumaDigitos (primeFactors x))

sumaDigitos :: Integer -> Integer
sumaDigitos x | x < 10 = x
              | otherwise = x `mod` 10 + sumaDigitos (x `div` 10)

smith :: [Integer]
smith = [x | x <- [1..], esSmith x]

import Data.Numbers.Primes

esSmith :: Integer -> Bool

esSmith x =

not (isPrime x) &&

sumaDigitos x == sum (map sumaDigitos (primeFactors x))

sumaDigitos :: Integer -> Integer

sumaDigitos x | x < 10 = x

| otherwise = x `mod` 10 + sumaDigitos (x `div` 10)

smith :: [Integer]

smith = [x | x <- [1..], esSmith x]

Múltiplos especiales

PorJosé A. Alonso 7-mayo-201515-mayo-2015

Dado dos números n y m, decimos que m es un múltiplo especial de n si m es un múltiplo de n y m no tiene ningún factor primo que sea congruente con 1 módulo 3.

Definir la función

   multiplosEspecialesCota :: Int -> Int -> [Int]

1	multiplosEspecialesCota :: Int -> Int -> [Int]

tal que (multiplosEspecialesCota n k) es la lista ordenada de todos los múltiplos especiales de n que son menores o iguales que k. Por ejemplo,

   multiplosEspecialesCota 5 50  ==  [5,10,15,20,25,30,40,45,50]
   multiplosEspecialesCota 7 50  ==  []

1 2	multiplosEspecialesCota 5 50 == [5,10,15,20,25,30,40,45,50] multiplosEspecialesCota 7 50 == []

Soluciones

import Data.Numbers.Primes (primeFactors)

multiplosEspecialesCota :: Int -> Int -> [Int]
multiplosEspecialesCota n k =
    [m | m <- [n,2*n..k], 
         all (\p -> p `mod` 3 /= 1) (primeFactors m)]

import Data.Numbers.Primes (primeFactors)

multiplosEspecialesCota :: Int -> Int -> [Int]

multiplosEspecialesCota n k =

[m | m <- [n,2*n..k],

all (\p -> p `mod` 3 /= 1) (primeFactors m)]

Números cuyas cifras coinciden con las de sus factores primos

PorJosé A. Alonso 4-mayo-201525-marzo-2022

Un número n es especial si al unir las cifras de sus factores primos, se obtienen exactamente las cifras de n, aunque puede ser en otro orden. Por ejemplo, 1255 es especial, pues los factores primos de 1255 son 5 y 251.

Definir la función

   esEspecial :: Integer -> Bool

1	esEspecial :: Integer -> Bool

tal que (esEspecial n) se verifica si un número n es especial. Por ejemplo,

   esEspecial 1255 == True
   esEspecial 125  == False

1 2	esEspecial 1255 == True esEspecial 125 == False

Comprobar con QuickCheck que todo número primo es especial.

Calcular los 5 primeros números especiales que no son primos.

Soluciones

import Data.List (nub, sort)
import Data.Numbers.Primes (isPrime, primeFactors)
import Test.QuickCheck

esEspecial :: Integer -> Bool
esEspecial n = 
    sort (show n) == sort (concatMap show (nub (primeFactors n)))

-- La propiedad es
prop_primos:: Integer -> Property
prop_primos n = 
    isPrime (abs n) ==> esEspecial (abs n)

-- La comprobación es
--    ghci> quickCheck prop_primos
--    +++ OK, passed 100 tests.

-- El cálculo es
--    ghci> take 5 [n | n <- [2..], esEspecial n, not (isPrime n)]
--    [735,1255,3792,7236,11913]

import Data.List (nub, sort)

import Data.Numbers.Primes (isPrime, primeFactors)

import Test.QuickCheck

esEspecial :: Integer -> Bool

esEspecial n =

sort (show n) == sort (concatMap show (nub (primeFactors n)))

-- La propiedad es

prop_primos:: Integer -> Property

prop_primos n =

isPrime (abs n) ==> esEspecial (abs n)

-- La comprobación es

-- ghci> quickCheck prop_primos

-- +++ OK, passed 100 tests.

-- El cálculo es

-- ghci> take 5 [n | n <- [2..], esEspecial n, not (isPrime n)]

-- [735,1255,3792,7236,11913]

Pares de enteros con sólo un factor primo común

PorJosé A. Alonso 31-marzo-201525-marzo-2022

Dos enteros positivos a y b se dirán relacionados si poseen, exactamente, un factor primo en común. Por ejemplo, 12 y 20 están relacionados, pero 6 y 30 no lo están.

Definir la lista infinita

   paresRel :: [(Int,Int)]

1	paresRel :: [(Int,Int)]

tal que paresRel enumera todos los pares (a,b), con 1 ≤ a < b, tal que a y b están relacionados. Por ejemplo,

   ghci> take 10 paresRel
   [(2,4),(2,6),(3,6),(4,6),(2,8),(4,8),(6,8),(3,9),(6,9),(2,10)]

1 2	ghci> take 10 paresRel [(2,4),(2,6),(3,6),(4,6),(2,8),(4,8),(6,8),(3,9),(6,9),(2,10)]

¿Qué lugar ocupa el par (51,111) en la lista infinita paresRel?

Soluciones

import Data.List (group, intersect, nub)
import Data.Numbers.Primes (primeFactors)

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

paresRel1 :: [(Int,Int)]
paresRel1 = [(a,b) | b <- [1..], a <- [1..b-1], relacionados a b]

relacionados :: Int -> Int -> Bool
relacionados a b = 
    length (nub (primeFactors a `intersect` primeFactors b)) == 1

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

paresRel2 :: [(Int,Int)]
paresRel2 = [(x,y) | y <- [4..], x <- [2..y-2], rel x y]
    where rel x y = m /= 1 && all (== head ps) ps
              where m  = gcd x y
                    ps = primeFactors m

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

paresRel3 :: [(Int,Int)]
paresRel3 =
    [(x,y) | y <- [2..], x <- [2..y-1], relacionados3 x y]

relacionados3 :: Int -> Int -> Bool
relacionados3 x y =
    length (group (primeFactors (gcd x y))) == 1

-- Comparación de eficiencia
--    λ> paresRel1 !! 40000
--    (216,489)
--    (3.19 secs, 1,825,423,056 bytes)
--    λ> paresRel2 !! 40000
--    (216,489)
--    (0.96 secs, 287,174,864 bytes)
--    λ> paresRel3 !! 40000
--    (216,489)
--    (0.70 secs, 264,137,928 bytes)

-- Cálculo
-- =======

-- El cálculo es
--    ghci> 1 + length (takeWhile (/=(51,111)) paresRel1)
--    2016

import Data.List (group, intersect, nub)

import Data.Numbers.Primes (primeFactors)

-- 1ª solución

-- ===========

paresRel1 :: [(Int,Int)]

paresRel1 = [(a,b) | b <- [1..], a <- [1..b-1], relacionados a b]

relacionados :: Int -> Int -> Bool

relacionados a b =

length (nub (primeFactors a `intersect` primeFactors b)) == 1

-- 2ª solución

-- ===========

paresRel2 :: [(Int,Int)]

paresRel2 = [(x,y) | y <- [4..], x <- [2..y-2], rel x y]

where rel x y = m /= 1 && all (== head ps) ps

where m = gcd x y

ps = primeFactors m

-- 3ª solución

-- ===========

paresRel3 :: [(Int,Int)]

paresRel3 =

[(x,y) | y <- [2..], x <- [2..y-1], relacionados3 x y]

relacionados3 :: Int -> Int -> Bool

relacionados3 x y =

length (group (primeFactors (gcd x y))) == 1

-- Comparación de eficiencia

-- λ> paresRel1 !! 40000

-- (216,489)

-- (3.19 secs, 1,825,423,056 bytes)

-- λ> paresRel2 !! 40000

-- (216,489)

-- (0.96 secs, 287,174,864 bytes)

-- λ> paresRel3 !! 40000

-- (216,489)

-- (0.70 secs, 264,137,928 bytes)

-- Cálculo

-- =======

-- El cálculo es

-- ghci> 1 + length (takeWhile (/=(51,111)) paresRel1)

-- 2016

Sucesión de números parientes

PorJosé A. Alonso 17-marzo-20151-mayo-2016

Se dice que dos números naturales son parientes sitienen exactamente un factor primo en común, independientemente de su multiplicidad. Por ejemplo,

Los números 12 (2²·3) y 40 (2³·5) son parientes, pues tienen al 2 como único factor primo en común.
Los números 49 (7²) y 63 (3²·7) son parientes, pues tienen al 7 como único factor primo en común.
Los números 12 (2²·3) y 30 (2·3·5) no son parientes, pues tienen dos factores primos en común.
Los números 49 (7²) y 25 (5²) no son parientes, pues no tienen factores primos en común.

Se dice que una lista de números naturales es una secuencia de parientes si cada par de números consecutivos son parientes. Por ejemplo,

La lista [12,40,35,28] es una secuencia de parientes.
La lista [12,30,21,49] no es una secuencia de parientes.

Definir la función

   secuenciaParientes :: [Integer] -> Bool

1	secuenciaParientes :: [Integer] -> Bool

tal que (secuenciaParientes xs) se verifica si xs es una secuencia de parientes. Por ejemplo,

   secuenciaParientes [12,40,35,28]           ==  True
   secuenciaParientes [12,30,21,49]           ==  False
   secuenciaParientes [2^n | n <- [1..2000]]  ==  True

secuenciaParientes [12,40,35,28] == True

secuenciaParientes [12,30,21,49] == False

secuenciaParientes [2^n | n <- [1..2000]] == True

Soluciones

import Data.List (intersect, nub)
import Data.Numbers.Primes (primes, primeFactors)

-- (parientes x y) se verifica si x e y son parientes. Por ejemplo,
--    parientes 12 40  ==  True
--    parientes 49 63  ==  True
--    parientes 12 30  ==  False
--    parientes 49 25  ==  False

-- 1ª definición (con gcd)
parientes1 :: Integer -> Integer -> Bool
parientes1 x y =
    length [p | p <- takeWhile (<= d) primes, d `mod` p == 0] == 1 
    where d = gcd x y

-- 2ª definición (con primeFactors)
parientes2 :: Integer -> Integer -> Bool
parientes2 0 0 = False
parientes2 x y = 
    length (nub (primeFactors x `intersect` primeFactors y)) == 1

-- Comparación de eficiencia
--    ghci> parientes1 (2^25) (2^25)
--    True
--    (34.34 secs, 15974866184 bytes)
--    ghci> parientes2 (2^25) (2^25)
--    True
--    (0.01 secs, 3093024 bytes)

-- Usaremos la 2ª definición
parientes :: Integer -> Integer -> Bool
parientes = parientes2

-- Definiciones de secuenciaParientes 
-- ==================================

-- 1ª definición (por recursión)
secuenciaParientes :: [Integer] -> Bool
secuenciaParientes []         = True
secuenciaParientes [x]        = True
secuenciaParientes (x1:x2:xs) =
    parientes x1 x2 && secuenciaParientes (x2:xs)

-- 2ª definición (por recursión con 2 ecuaciones)
secuenciaParientes2 :: [Integer] -> Bool
secuenciaParientes2 (x1:x2:xs) =
    parientes x1 x2 && secuenciaParientes2 (x2:xs)
secuenciaParientes2 _         = True

-- 3ª definición (sin recursión):
secuenciaParientes3 :: [Integer] -> Bool
secuenciaParientes3 xs = all (\(x,y) -> parientes x y) (zip xs (tail xs)) 

-- 4ª definición
secuenciaParientes4 :: [Integer] -> Bool
secuenciaParientes4 xs = all (uncurry parientes) (zip xs (tail xs))

import Data.List (intersect, nub)

import Data.Numbers.Primes (primes, primeFactors)

-- (parientes x y) se verifica si x e y son parientes. Por ejemplo,

-- parientes 12 40 == True

-- parientes 49 63 == True

-- parientes 12 30 == False

-- parientes 49 25 == False

-- 1ª definición (con gcd)

parientes1 :: Integer -> Integer -> Bool

parientes1 x y =

length [p | p <- takeWhile (<= d) primes, d `mod` p == 0] == 1

where d = gcd x y

-- 2ª definición (con primeFactors)

parientes2 :: Integer -> Integer -> Bool

parientes2 0 0 = False

parientes2 x y =

length (nub (primeFactors x `intersect` primeFactors y)) == 1

-- Comparación de eficiencia

-- ghci> parientes1 (2^25) (2^25)

-- True

-- (34.34 secs, 15974866184 bytes)

-- ghci> parientes2 (2^25) (2^25)

-- True

-- (0.01 secs, 3093024 bytes)

-- Usaremos la 2ª definición

parientes :: Integer -> Integer -> Bool

parientes = parientes2

-- Definiciones de secuenciaParientes

-- ==================================

-- 1ª definición (por recursión)

secuenciaParientes :: [Integer] -> Bool

secuenciaParientes [] = True

secuenciaParientes [x] = True

secuenciaParientes (x1:x2:xs) =

parientes x1 x2 && secuenciaParientes (x2:xs)

-- 2ª definición (por recursión con 2 ecuaciones)

secuenciaParientes2 :: [Integer] -> Bool

secuenciaParientes2 (x1:x2:xs) =

parientes x1 x2 && secuenciaParientes2 (x2:xs)

secuenciaParientes2 _ = True

-- 3ª definición (sin recursión):

secuenciaParientes3 :: [Integer] -> Bool

secuenciaParientes3 xs = all (\(x,y) -> parientes x y) (zip xs (tail xs))

-- 4ª definición

secuenciaParientes4 :: [Integer] -> Bool

secuenciaParientes4 xs = all (uncurry parientes) (zip xs (tail xs))

Menor número triangular con más de n divisores

PorJosé A. Alonso 6-febrero-201525-abril-2021

La sucesión de los números triangulares se obtiene sumando los números naturales.

Así, el 7º número triangular es

   1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.

1	1 + 2 + 3 + 4 + 5 + 6 + 7 = 28.

Los primeros 10 números triangulares son

   1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

1	1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Los divisores de los primeros 7 números triangulares son:

    1: 1
    3: 1,3
    6: 1,2,3,6
   10: 1,2,5,10
   15: 1,3,5,15
   21: 1,3,7,21
   28: 1,2,4,7,14,28

1: 1

3: 1,3

6: 1,2,3,6

10: 1,2,5,10

15: 1,3,5,15

21: 1,3,7,21

28: 1,2,4,7,14,28

Como se puede observar, 28 es el menor número triangular con más de 5 divisores.

Definir la función

   menorTriangularConAlMenosNDivisores :: Int -> Integer

1	menorTriangularConAlMenosNDivisores :: Int -> Integer

tal que (menorTriangularConAlMenosNDivisores n) es el menor número triangular que tiene al menos n divisores. Por ejemplo,

   menorTriangularConAlMenosNDivisores 5    ==  28
   menorTriangularConAlMenosNDivisores 50   ==  25200
   menorTriangularConAlMenosNDivisores 500  ==  76576500

menorTriangularConAlMenosNDivisores 5 == 28

menorTriangularConAlMenosNDivisores 50 == 25200

menorTriangularConAlMenosNDivisores 500 == 76576500

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primes,primeFactors)

-- 1ª definición
-- =============

menorTriangularConAlMenosNDivisores1 :: Int -> Integer
menorTriangularConAlMenosNDivisores1 n = 
    head [x | x <- triangulares, nDivisores x >= n]

-- triangulares es la sucesión de los números triangulares. Por ejemplo,
--    take 10 triangulares  ==  [1,3,6,10,15,21,28,36,45,55]
triangulares :: [Integer]
triangulares = scanl (+) 1 [2..]

-- (nDivisores x) es el número de divisores de x. Por ejemplo,
--    nDivisores 28  ==  6
nDivisores :: Integer -> Int
nDivisores x = 
    1 + length [y | y <- [1..x `div` 2], mod x y == 0]

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

menorTriangularConAlMenosNDivisores2 :: Int -> Integer
menorTriangularConAlMenosNDivisores2 n = 
    head [x | x <- triangulares, nDivisores2 x >= n]

-- (nDivisores2 x) es el número de divisores de x. Por ejemplo,
--    nDivisores2 28  ==  6
nDivisores2 :: Integer -> Int
nDivisores2 n = product [1 + length xs | xs <- group (factoresPrimos n)]    

-- (factoresPrimos n) es la lista de los factores primos de n. Por
-- ejemplo, 
--    factoresPrimos 28  ==  [2,2,7]
factoresPrimos :: Integer -> [Integer]
factoresPrimos n = aux n primos
    where aux n (p:ps) 
              | p*p > n        = [n]
              | n `mod` p == 0 = p : aux (n `div` p) (p:ps)
              | otherwise      = aux n ps

-- primos es la lista de los números primos. Por ejemplo,
--    take 10 primos  ==  [2,3,5,7,11,13,17,19,23,29]
primos :: [Integer]
primos = 2 : filter ((==1) . length . factoresPrimos) [3,5..]

-- 3ª solución (usando primes)
-- ===========================

menorTriangularConAlMenosNDivisores3 :: Int -> Integer
menorTriangularConAlMenosNDivisores3 n = 
    head [x | x <- triangulares, nDivisores3 x >= n]

-- (nDivisores3 x) es el número de divisores de x. Por ejemplo,
--    nDivisores3 28  ==  6
nDivisores3 :: Integer -> Int
nDivisores3 n = product [1 + length xs | xs <- group (factoresPrimos3 n)]    

-- (factoresPrimos3 n) es la lista de los factores primos de n. Por
-- ejemplo, 
--    factoresPrimos3 28  ==  [2,2,7]
factoresPrimos3 n = aux n primes
  where
    aux n (p:ps) 
        | p*p > n        = [n]
        | n `mod` p == 0 = p : aux (n `div` p) (p:ps)
        | otherwise      = aux n ps

-- 4ª solución (usando primeFactors)
-- =================================

menorTriangularConAlMenosNDivisores4 :: Int -> Integer
menorTriangularConAlMenosNDivisores4 n = 
    head [x | x <- triangulares, nDivisores4 x >= n]

-- (nDivisores4 x) es el número de divisores de x. Por ejemplo,
--    nDivisores4 28  ==  6
nDivisores4 :: Integer -> Int
nDivisores4 n = product [1 + length xs | xs <- group (primeFactors n)]    

-- ---------------------------------------------------------------------
-- § Comparación de eficiencia                                        --
-- ---------------------------------------------------------------------

-- La comparación es
--    ghci> menorTriangularConAlMenosNDivisores1 50
--    25200
--    (1.25 secs, 200236512 bytes)
--    
--    ghci> menorTriangularConAlMenosNDivisores2 50
--    25200
--    (0.02 secs, 4199904 bytes)
--    
--    ghci> menorTriangularConAlMenosNDivisores3 50
--    25200
--    (0.03 secs, 6265128 bytes)
--    
--    ghci> menorTriangularConAlMenosNDivisores4 50
--    25200
--    (0.01 secs, 5753048 bytes)

100

101

102

import Data.List (group)

import Data.Numbers.Primes (primes,primeFactors)

-- 1ª definición

-- =============

menorTriangularConAlMenosNDivisores1 :: Int -> Integer

menorTriangularConAlMenosNDivisores1 n =

head [x | x <- triangulares, nDivisores x >= n]

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

-- take 10 triangulares == [1,3,6,10,15,21,28,36,45,55]

triangulares :: [Integer]

triangulares = scanl (+) 1 [2..]

-- (nDivisores x) es el número de divisores de x. Por ejemplo,

-- nDivisores 28 == 6

nDivisores :: Integer -> Int

nDivisores x =

1 + length [y | y <- [1..x `div` 2], mod x y == 0]

-- 2ª solución

-- ===========

menorTriangularConAlMenosNDivisores2 :: Int -> Integer

menorTriangularConAlMenosNDivisores2 n =

head [x | x <- triangulares, nDivisores2 x >= n]

-- (nDivisores2 x) es el número de divisores de x. Por ejemplo,

-- nDivisores2 28 == 6

nDivisores2 :: Integer -> Int

nDivisores2 n = product [1 + length xs | xs <- group (factoresPrimos n)]

-- (factoresPrimos n) es la lista de los factores primos de n. Por

-- ejemplo,

-- factoresPrimos 28 == [2,2,7]

factoresPrimos :: Integer -> [Integer]

factoresPrimos n = aux n primos

where aux n (p:ps)

| p*p > n = [n]

| n `mod` p == 0 = p : aux (n `div` p) (p:ps)

| otherwise = aux n ps

-- primos es la lista de los números primos. Por ejemplo,

-- take 10 primos == [2,3,5,7,11,13,17,19,23,29]

primos :: [Integer]

primos = 2 : filter ((==1) . length . factoresPrimos) [3,5..]

-- 3ª solución (usando primes)

-- ===========================

menorTriangularConAlMenosNDivisores3 :: Int -> Integer

menorTriangularConAlMenosNDivisores3 n =

head [x | x <- triangulares, nDivisores3 x >= n]

-- (nDivisores3 x) es el número de divisores de x. Por ejemplo,

-- nDivisores3 28 == 6

nDivisores3 :: Integer -> Int

nDivisores3 n = product [1 + length xs | xs <- group (factoresPrimos3 n)]

-- (factoresPrimos3 n) es la lista de los factores primos de n. Por

-- ejemplo,

-- factoresPrimos3 28 == [2,2,7]

factoresPrimos3 n = aux n primes

where

aux n (p:ps)

| p*p > n = [n]

| n `mod` p == 0 = p : aux (n `div` p) (p:ps)

| otherwise = aux n ps

-- 4ª solución (usando primeFactors)

-- =================================

menorTriangularConAlMenosNDivisores4 :: Int -> Integer

menorTriangularConAlMenosNDivisores4 n =

head [x | x <- triangulares, nDivisores4 x >= n]

-- (nDivisores4 x) es el número de divisores de x. Por ejemplo,

-- nDivisores4 28 == 6

nDivisores4 :: Integer -> Int

nDivisores4 n = product [1 + length xs | xs <- group (primeFactors n)]

-- ---------------------------------------------------------------------

-- § Comparación de eficiencia --

-- ---------------------------------------------------------------------

-- La comparación es

-- ghci> menorTriangularConAlMenosNDivisores1 50

-- 25200

-- (1.25 secs, 200236512 bytes)

-- ghci> menorTriangularConAlMenosNDivisores2 50

-- 25200

-- (0.02 secs, 4199904 bytes)

-- ghci> menorTriangularConAlMenosNDivisores3 50

-- 25200

-- (0.03 secs, 6265128 bytes)

-- ghci> menorTriangularConAlMenosNDivisores4 50

-- 25200

-- (0.01 secs, 5753048 bytes)

Números primos de Pierpont

Soluciones

Pensamiento

Grafo de divisibilidad

Soluciones

Números compuestos por un conjunto de primos

Soluciones

Alturas primas

Soluciones

Exponentes de Hamming

Soluciones

Cálculo de pi mediante la variante de Euler de la serie armónica

Soluciones

Números cuyos dígitos coinciden con los de sus factores primos

Soluciones

Los números de Smith

Soluciones

Múltiplos especiales

Soluciones

Números cuyas cifras coinciden con las de sus factores primos

Soluciones

Sucesión de números parientes

Soluciones

Entradas recientes

Buscar

Correo electrónico