Ejercicio

Potencias perfectas

PorJosé A. Alonso 20-julio-202213-julio-2022

Un número natural n es una potencia perfecta si existen dos números naturales m > 1 y k > 1 tales que n = m^k. Las primeras potencias perfectas son

   4 = 2², 8 = 2³, 9 = 3², 16 = 2⁴, 25 = 5², 27 = 3³, 32 = 2⁵, 
   36 = 6², 49 = 7², 64 = 2⁶, ...

1 2	4 = 2², 8 = 2³, 9 = 3², 16 = 2⁴, 25 = 5², 27 = 3³, 32 = 2⁵, 36 = 6², 49 = 7², 64 = 2⁶, ...

Definir la sucesión

   potenciasPerfectas :: [Integer]

1	potenciasPerfectas :: [Integer]

cuyos términos son las potencias perfectas. Por ejemplo,

   take 10 potenciasPerfectas  ==  [4,8,9,16,25,27,32,36,49,64]
   potenciasPerfectas !! 3000  ==  7778521

1 2	take 10 potenciasPerfectas == [4,8,9,16,25,27,32,36,49,64] potenciasPerfectas !! 3000 == 7778521

Definir el procedimiento

   grafica :: Int -> IO ()

1	grafica :: Int -> IO ()

tal que (grafica n) es la representación gráfica de las n primeras potencias perfectas. Por ejemplo, para (grafica 30) dibuja

Soluciones

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

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

potenciasPerfectas1 :: [Integer]
potenciasPerfectas1 = filter esPotenciaPerfecta [4..]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta 36  ==  True
--    esPotenciaPerfecta 72  ==  False
esPotenciaPerfecta :: Integer -> Bool
esPotenciaPerfecta = not . null. potenciasPerfectasDe 

-- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que 
-- x = a^b. Por ejemplo,
--    potenciasPerfectasDe 64  ==  [(2,6),(4,3),(8,2)]
--    potenciasPerfectasDe 72  ==  []
potenciasPerfectasDe :: Integer -> [(Integer,Integer)]
potenciasPerfectasDe n = 
  [(m,k) | m <- takeWhile (\x -> x*x <= n) [2..]
         , k <- takeWhile (\x -> m^x <= n) [2..]
         , m^k == n]

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

potenciasPerfectas2 :: [Integer]
potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x]

-- (esPotenciaPerfecta2 x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta2 36  ==  True
--    esPotenciaPerfecta2 72  ==  False
esPotenciaPerfecta2 :: Integer -> Bool
esPotenciaPerfecta2 x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de l factorización prima
-- de x. Por ejemplos,
--    exponentes 36  ==  [2,2]
--    exponentes 72  ==  [3,2]
exponentes :: Integer -> [Int]
exponentes x = [length ys | ys <- group (primeFactors x)] 

-- (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 :: [Int] -> Int
mcd = foldl1 gcd

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

potenciasPerfectas3 :: [Integer]
potenciasPerfectas3 = mezclaTodas potencias

-- potencias es la lista las listas de potencias de todos los números
-- mayores que 1 con exponentes mayores que 1. Por ejemplo,
--    λ> map (take 3) (take 4 potencias)
--    [[4,8,16],[9,27,81],[16,64,256],[25,125,625]]
potencias:: [[Integer]]
potencias = [[n^k | k <- [2..]] | n <- [2..]]

-- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las
-- listas ordenadas xss. Por ejemplo,
--    take 7 (mezclaTodas potencias)  ==  [4,8,9,16,25,27,32]
mezclaTodas :: Ord a => [[a]] -> [a]
mezclaTodas = foldr1 xmezcla
  where xmezcla (x:xs) ys = x : mezcla xs ys

-- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las
-- listas ordenadas xs e ys. Por ejemplo,
--    take 7 (mezcla [2,5..] [4,6..])  ==  [2,4,5,6,8,10,11]
mezcla :: Ord a => [a] -> [a] -> [a]
mezcla (x:xs) (y:ys) | x < y  = x : mezcla xs (y:ys)
                     | x == y = x : mezcla xs ys
                     | x > y  = y : mezcla (x:xs) ys

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

-- La propiedad es
prop_potenciasPerfectas :: NonNegative Int -> Bool
prop_potenciasPerfectas (NonNegative n) =
  all (== potenciasPerfectas1 !! n)
      [potenciasPerfectas2 !! n,
       potenciasPerfectas3 !! n]

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

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

-- La comparación es
--    λ> potenciasPerfectas1 !! 200
--    28224
--    (10.56 secs, 8,434,647,368 bytes)
--    λ> potenciasPerfectas2 !! 200
--    28224
--    (0.36 secs, 825,040,416 bytes)
--    λ> potenciasPerfectas3 !! 200
--    28224
--    (0.05 secs, 7,474,280 bytes)
--    
--    λ> potenciasPerfectas2 !! 500
--    191844
--    (4.16 secs, 9,899,367,112 bytes)
--    λ> potenciasPerfectas3 !! 500
--    191844
--    (0.09 secs, 51,275,464 bytes)

-- En lo que sigue se usa la 3ª solución
potenciasPerfectas :: [Integer]
potenciasPerfectas = potenciasPerfectas3

-- Representación gráfica
-- ======================

grafica :: Int -> IO ()
grafica n = 
  plotList [ Key Nothing
           , PNG "Potencias_perfectas.png"
           ]
           (take n potenciasPerfectas)

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

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

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

import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

-- 1ª solución

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

potenciasPerfectas1 :: [Integer]

potenciasPerfectas1 = filter esPotenciaPerfecta [4..]

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

-- ejemplo,

-- esPotenciaPerfecta 36 == True

-- esPotenciaPerfecta 72 == False

esPotenciaPerfecta :: Integer -> Bool

esPotenciaPerfecta = not . null. potenciasPerfectasDe

-- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que

-- x = a^b. Por ejemplo,

-- potenciasPerfectasDe 64 == [(2,6),(4,3),(8,2)]

-- potenciasPerfectasDe 72 == []

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

potenciasPerfectasDe n =

[(m,k) | m <- takeWhile (\x -> x*x <= n) [2..]

, k <- takeWhile (\x -> m^x <= n) [2..]

, m^k == n]

-- 2ª solución

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

potenciasPerfectas2 :: [Integer]

potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x]

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

-- ejemplo,

-- esPotenciaPerfecta2 36 == True

-- esPotenciaPerfecta2 72 == False

esPotenciaPerfecta2 :: Integer -> Bool

esPotenciaPerfecta2 x = mcd (exponentes x) > 1

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

-- de x. Por ejemplos,

-- exponentes 36 == [2,2]

-- exponentes 72 == [3,2]

exponentes :: Integer -> [Int]

exponentes x = [length ys | ys <- group (primeFactors x)]

-- (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 :: [Int] -> Int

mcd = foldl1 gcd

-- 3ª solución

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

potenciasPerfectas3 :: [Integer]

potenciasPerfectas3 = mezclaTodas potencias

-- potencias es la lista las listas de potencias de todos los números

-- mayores que 1 con exponentes mayores que 1. Por ejemplo,

-- λ> map (take 3) (take 4 potencias)

-- [[4,8,16],[9,27,81],[16,64,256],[25,125,625]]

potencias:: [[Integer]]

potencias = [[n^k | k <- [2..]] | n <- [2..]]

-- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las

-- listas ordenadas xss. Por ejemplo,

-- take 7 (mezclaTodas potencias) == [4,8,9,16,25,27,32]

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

mezclaTodas = foldr1 xmezcla

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

-- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las

-- listas ordenadas xs e ys. Por ejemplo,

-- take 7 (mezcla [2,5..] [4,6..]) == [2,4,5,6,8,10,11]

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

mezcla (x:xs) (y:ys) | x < y = x : mezcla xs (y:ys)

| x == y = x : mezcla xs ys

| x > y = y : mezcla (x:xs) ys

-- Comprobación de equivalencia

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

-- La propiedad es

prop_potenciasPerfectas :: NonNegative Int -> Bool

prop_potenciasPerfectas (NonNegative n) =

all (== potenciasPerfectas1 !! n)

[potenciasPerfectas2 !! n,

potenciasPerfectas3 !! n]

-- La comprobación es

-- λ> quickCheck prop_potenciasPerfectas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> potenciasPerfectas1 !! 200

-- 28224

-- (10.56 secs, 8,434,647,368 bytes)

-- λ> potenciasPerfectas2 !! 200

-- 28224

-- (0.36 secs, 825,040,416 bytes)

-- λ> potenciasPerfectas3 !! 200

-- 28224

-- (0.05 secs, 7,474,280 bytes)

-- λ> potenciasPerfectas2 !! 500

-- 191844

-- (4.16 secs, 9,899,367,112 bytes)

-- λ> potenciasPerfectas3 !! 500

-- 191844

-- (0.09 secs, 51,275,464 bytes)

-- En lo que sigue se usa la 3ª solución

potenciasPerfectas :: [Integer]

potenciasPerfectas = potenciasPerfectas3

-- Representación gráfica

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

grafica :: Int -> IO ()

grafica n =

plotList [ Key Nothing

, PNG "Potencias_perfectas.png"

]

(take n potenciasPerfectas)

El código se encuentra en GitHub.

Sumas alternas de factoriales

PorJosé A. Alonso 19-julio-202212-julio-2022

Las primeras sumas alternas de los factoriales son números primos; en efecto,

   3! - 2! + 1! = 5
   4! - 3! + 2! - 1! = 19
   5! - 4! + 3! - 2! + 1! = 101
   6! - 5! + 4! - 3! + 2! - 1! = 619
   7! - 6! + 5! - 4! + 3! - 2! + 1! = 4421
   8! - 7! + 6! - 5! + 4! - 3! + 2! - 1! = 35899

3! - 2! + 1! = 5

4! - 3! + 2! - 1! = 19

5! - 4! + 3! - 2! + 1! = 101

6! - 5! + 4! - 3! + 2! - 1! = 619

7! - 6! + 5! - 4! + 3! - 2! + 1! = 4421

8! - 7! + 6! - 5! + 4! - 3! + 2! - 1! = 35899

son primos, pero

   9! - 8! + 7! - 6! + 5! - 4! + 3! - 2! + 1! = 326981

1	9! - 8! + 7! - 6! + 5! - 4! + 3! - 2! + 1! = 326981

no es primo.

Definir las funciones

   sumaAlterna         :: Integer -> Integer
   sumasAlternas       :: [Integer]
   conSumaAlternaPrima :: [Integer]

sumaAlterna :: Integer -> Integer

sumasAlternas :: [Integer]

conSumaAlternaPrima :: [Integer]

tales que

(sumaAlterna n) es la suma alterna de los factoriales desde n hasta 1. Por ejemplo,

     sumaAlterna 3  ==  5
     sumaAlterna 4  ==  19
     sumaAlterna 5  ==  101
     sumaAlterna 6  ==  619
     sumaAlterna 7  ==  4421
     sumaAlterna 8  ==  35899
     sumaAlterna 9  ==  326981
     sumaAlterna (5*10^4) `mod` (10^6) == 577019

sumaAlterna 3 == 5

sumaAlterna 4 == 19

sumaAlterna 5 == 101

sumaAlterna 6 == 619

sumaAlterna 7 == 4421

sumaAlterna 8 == 35899

sumaAlterna 9 == 326981

sumaAlterna (5*10^4) `mod` (10^6) == 577019

sumasAlternas es la sucesión de las sumas alternas de factoriales. Por ejemplo,

     λ> take 10 sumasAlternas1
     [0,1,1,5,19,101,619,4421,35899,326981]

1 2	λ> take 10 sumasAlternas1 [0,1,1,5,19,101,619,4421,35899,326981]

conSumaAlternaPrima es la sucesión de los números cuya suma alterna de factoriales es prima. Por ejemplo,

     λ> take 8 conSumaAlternaPrima
     [3,4,5,6,7,8,10,15]

1 2	λ> take 8 conSumaAlternaPrima [3,4,5,6,7,8,10,15]

Soluciones

import Data.List (genericTake)
import Data.Numbers.Primes (isPrime)
import Test.QuickCheck

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

sumaAlterna1 :: Integer -> Integer
sumaAlterna1 1 = 1
sumaAlterna1 n = factorial n - sumaAlterna1 (n-1)

factorial :: Integer -> Integer
factorial n = product [1..n]

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

sumaAlterna2 :: Integer -> Integer
sumaAlterna2 n =
  sum (genericTake n (zipWith (*) signos (tail factoriales)))
  where
    signos | odd n     = concat (repeat [1,-1])
           | otherwise = concat (repeat [-1,1])

-- factoriales es la lista de los factoriales. Por ejemplo,
--    take 7 factoriales  ==  [1,1,2,6,24,120,720]
factoriales :: [Integer]
factoriales = 1 : scanl1 (*) [1..]

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

sumaAlterna3 :: Integer -> Integer
sumaAlterna3 n = 
  sum (genericTake n (zipWith (*) signos (tail factoriales)))
  where signos | odd n     = cycle [1,-1]
               | otherwise = cycle [-1,1]

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

sumaAlterna4 :: Integer -> Integer
sumaAlterna4 n =
  foldl (flip (-)) 0 (scanl1 (*) [1..n])

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

-- La propiedad es
prop_sumaAlterna :: Positive Integer -> Bool 
prop_sumaAlterna (Positive n) =
  all (== sumaAlterna1 n)
      [sumaAlterna2 n,
       sumaAlterna3 n,
       sumaAlterna4 n]

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

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

-- La comparación es
--    λ> sumaAlterna1 4000 `mod` (10^6)
--    577019
--    (6.21 secs, 16,154,113,192 bytes)
--    λ> sumaAlterna2 4000 `mod` (10^6)
--    577019
--    (0.01 secs, 24,844,664 bytes)
--    
--    λ> sumaAlterna2 (5*10^4) `mod` (10^6)
--    577019
--    (1.81 secs, 4,729,583,864 bytes)
--    λ> sumaAlterna3 (5*10^4) `mod` (10^6)
--    577019
--    (0.89 secs, 4,725,983,928 bytes)
--    λ> sumaAlterna4 (5*10^4) `mod` (10^6)
--    577019
--    (0.70 secs, 4,710,770,592 bytes)

-- En lo que sigue se usa la 3ª definición
sumaAlterna :: Integer -> Integer
sumaAlterna = sumaAlterna3

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

sumasAlternas1 :: [Integer]
sumasAlternas1 =
  map sumaAlterna [0..]

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

sumasAlternas2 :: [Integer]
sumasAlternas2 =
  0 : zipWith (-) (tail factoriales) sumasAlternas2

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

sumasAlternas3 :: [Integer]
sumasAlternas3 =
  scanl (flip (-)) 0 $ scanl1 (*) [1..]

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

-- La propiedad es
prop_sumasAlternas :: NonNegative Int -> Bool
prop_sumasAlternas (NonNegative n) =
  all (== sumasAlternas1 !! n)
      [sumasAlternas2 !! n,
       sumasAlternas3 !! n]
  
-- La comprobación es
--    λ> quickCheck prop_sumasAlternas
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> length (show (sumasAlternas1 !! (5*10^4)))
--    213237
--    (4.90 secs, 4,731,620,600 bytes)
--    λ> length (show (sumasAlternas2 !! (5*10^4)))
--    213237
--    (2.39 secs, 4,726,820,456 bytes)
--    λ> length (show (sumasAlternas3 !! (5*10^4)))
--    213237
--    (1.78 secs, 4,726,820,280 bytes)

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

conSumaAlternaPrima1 :: [Integer]
conSumaAlternaPrima1 =
  [n | n <- [0..], isPrime (sumaAlterna n)]

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

conSumaAlternaPrima2 :: [Integer]
conSumaAlternaPrima2 =
  [x | (x,y) <- zip [0..] sumasAlternas2, isPrime y]

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

conSumaAlternaPrima3 :: [Integer]
conSumaAlternaPrima3 =
  filter (isPrime . sumaAlterna) [0..]

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

-- La propiedad es
prop_conSumaAlternaPrima :: NonNegative Int -> Bool
prop_conSumaAlternaPrima (NonNegative n) =
  all (== conSumaAlternaPrima1 !! n)
      [conSumaAlternaPrima2 !! n,
       conSumaAlternaPrima3 !! n]

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=5}) prop_conSumaAlternaPrima
--    +++ 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.List (genericTake)

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck

-- 1ª definición de sumaAlterna

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

sumaAlterna1 :: Integer -> Integer

sumaAlterna1 1 = 1

sumaAlterna1 n = factorial n - sumaAlterna1 (n-1)

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª definición de sumaAlterna

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

sumaAlterna2 :: Integer -> Integer

sumaAlterna2 n =

sum (genericTake n (zipWith (*) signos (tail factoriales)))

where

signos | odd n = concat (repeat [1,-1])

| otherwise = concat (repeat [-1,1])

-- factoriales es la lista de los factoriales. Por ejemplo,

-- take 7 factoriales == [1,1,2,6,24,120,720]

factoriales :: [Integer]

factoriales = 1 : scanl1 (*) [1..]

-- 3ª definición de sumaAlterna

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

sumaAlterna3 :: Integer -> Integer

sumaAlterna3 n =

sum (genericTake n (zipWith (*) signos (tail factoriales)))

where signos | odd n = cycle [1,-1]

| otherwise = cycle [-1,1]

-- 3ª definición de sumaAlterna

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

sumaAlterna4 :: Integer -> Integer

sumaAlterna4 n =

foldl (flip (-)) 0 (scanl1 (*) [1..n])

-- Comprobación de equivalencia de sumaAlterna

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

-- La propiedad es

prop_sumaAlterna :: Positive Integer -> Bool

prop_sumaAlterna (Positive n) =

all (== sumaAlterna1 n)

[sumaAlterna2 n,

sumaAlterna3 n,

sumaAlterna4 n]

-- La comprobación es

-- λ> quickCheck prop_sumaAlterna

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de sumaAlterna

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

-- La comparación es

-- λ> sumaAlterna1 4000 `mod` (10^6)

-- 577019

-- (6.21 secs, 16,154,113,192 bytes)

-- λ> sumaAlterna2 4000 `mod` (10^6)

-- 577019

-- (0.01 secs, 24,844,664 bytes)

-- λ> sumaAlterna2 (5*10^4) `mod` (10^6)

-- 577019

-- (1.81 secs, 4,729,583,864 bytes)

-- λ> sumaAlterna3 (5*10^4) `mod` (10^6)

-- 577019

-- (0.89 secs, 4,725,983,928 bytes)

-- λ> sumaAlterna4 (5*10^4) `mod` (10^6)

-- 577019

-- (0.70 secs, 4,710,770,592 bytes)

-- En lo que sigue se usa la 3ª definición

sumaAlterna :: Integer -> Integer

sumaAlterna = sumaAlterna3

-- 1ª definición de sumasAlternas

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

sumasAlternas1 :: [Integer]

sumasAlternas1 =

map sumaAlterna [0..]

-- 2ª definición de sumasAlternas

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

sumasAlternas2 :: [Integer]

sumasAlternas2 =

0 : zipWith (-) (tail factoriales) sumasAlternas2

-- 3ª definición de sumasAlternas

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

sumasAlternas3 :: [Integer]

sumasAlternas3 =

scanl (flip (-)) 0 $ scanl1 (*) [1..]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sumasAlternas :: NonNegative Int -> Bool

prop_sumasAlternas (NonNegative n) =

all (== sumasAlternas1 !! n)

[sumasAlternas2 !! n,

sumasAlternas3 !! n]

-- La comprobación es

-- λ> quickCheck prop_sumasAlternas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (show (sumasAlternas1 !! (5*10^4)))

-- 213237

-- (4.90 secs, 4,731,620,600 bytes)

-- λ> length (show (sumasAlternas2 !! (5*10^4)))

-- 213237

-- (2.39 secs, 4,726,820,456 bytes)

-- λ> length (show (sumasAlternas3 !! (5*10^4)))

-- 213237

-- (1.78 secs, 4,726,820,280 bytes)

-- 1ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima1 :: [Integer]

conSumaAlternaPrima1 =

[n | n <- [0..], isPrime (sumaAlterna n)]

-- 2ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima2 :: [Integer]

conSumaAlternaPrima2 =

[x | (x,y) <- zip [0..] sumasAlternas2, isPrime y]

-- 3ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima3 :: [Integer]

conSumaAlternaPrima3 =

filter (isPrime . sumaAlterna) [0..]

-- Comprobación de equivalencia de conSumaAlternaPrima

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

-- La propiedad es

prop_conSumaAlternaPrima :: NonNegative Int -> Bool

prop_conSumaAlternaPrima (NonNegative n) =

all (== conSumaAlternaPrima1 !! n)

[conSumaAlternaPrima2 !! n,

conSumaAlternaPrima3 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=5}) prop_conSumaAlternaPrima

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Primos con cubos

PorJosé A. Alonso 18-julio-202211-julio-2022

Un primo con cubo es un número primo p para el que existe algún entero positivo n tal que la expresión n^3 + n^2p es un cubo perfecto. Por ejemplo, 19 es un primo con cubo ya que 8^3 + 8^2×19 = 12^3.

Definir la sucesión

   primosConCubos :: [Integer]

1	primosConCubos :: [Integer]

tal que sus elementos son los primos con cubo. Por ejemplo,

   λ> take 6 primosConCubos
   [7,19,37,61,127,271]
   λ> length (takeWhile (< 1000000) primosConCubos)
   173

λ> take 6 primosConCubos

[7,19,37,61,127,271]

λ> length (takeWhile (< 1000000) primosConCubos)

173

Soluciones

import Data.Numbers.Primes (isPrime)
import Test.QuickCheck (NonNegative (NonNegative), maxSize, quickCheckWith, stdArgs) 

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

primosConCubos1 :: [Integer]
primosConCubos1 =
  [p | x <- [1..],
       n <- [1..x],
       (x^3 - n^3) `mod` (n^2) == 0,
       let p = (x^3 - n^3) `div` (n^2),
       isPrime p]

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

-- Para analizar la respuesta, en esta solución se calculan los pares
-- (p,n) tales que p es un primo con cubo y n es un número positivo tal
-- que n^3 + n^2p es un cubo

primosConCubos2' :: [(Integer,Integer)]
primosConCubos2' =
  [(p,n) | x <- [1..],
           n <- [1..x],
           (x^3 - n^3) `mod` (n^2) == 0,
           let p = (x^3 - n^3) `div` (n^2),
           isPrime p]

-- El cálculo es
--    λ> take 7 primosConCubos2
--    [(7,1),(19,8),(37,27),(61,64),(127,216),(271,729),(331,1000)]

-- Se observa que la sucesión de los segundos elementos [1,8,27,64,...]
-- es la de los cubos y que los primeros elementos se obtienen restando
-- los segundos elementos consecutivos; es decir,
--     7 =  8 -  1 = 2^3 - 1^3  
--    19 = 27 -  8 = 3^3 - 2^3
--    37 = 64 - 27 = 4^3 - 3^3
-- Continuando el patrón,
--     61 =  5^3 - 4^3 =  125 -   64
--     91 =  6^3 - 5^3 =  216 -  125
--    127 =  7^3 - 6^3 =  343 -  216
--    271 = 10^3 - 9^3 = 1000 -  729
--    331 = 11^3 -10^3 = 1331 - 1000
-- Por tanto, los primos con cubos son diferencias de dos cubos
-- consecutivos; es decir, coinciden con los números cubanos del
-- ejercicio anterior. A partir de la conjetura se obtienen las
-- siguientes definiciones

-- Basado en las anteriores observaciones se obtiene la siguiente
-- definición 
primosConCubos2 :: [Integer]
primosConCubos2 = 
  filter isPrime [(x+1)^3 - x^3 | x <- [1..]] 

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

primosConCubos3 :: [Integer]
primosConCubos3 =
  filter isPrime diferenciasCubosConsecutivos

diferenciasCubosConsecutivos :: [Integer]
diferenciasCubosConsecutivos =
  zipWith (-) (tail cubos) cubos 

cubos :: [Integer]
cubos = map (^3) [0..]

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

-- Simplificando la expresión (x+1)^3 - x^3 se obtiene 3*x^2+ 3*x + 1,
-- con lo que la 3ª definición se reduce a

primosConCubos4 :: [Integer]
primosConCubos4 =
  [p | x <- [1..],
       let p = 3*x^2+ 3*x + 1,
       isPrime p]

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

-- La propiedad es
prop_primosConCubos :: NonNegative Int -> Bool
prop_primosConCubos (NonNegative n) =
  all (== primosConCubos1 !! n)
      [primosConCubos2 !! n,
       primosConCubos3 !! n,
       primosConCubos4 !! n]
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_primosConCubos
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> primosConCubos1 !! 6
--    331
--    (1.15 secs, 1,105,612,336 bytes)
--    λ> primosConCubos1 !! 7
--    397
--    (1.96 secs, 1,909,369,592 bytes)
--    λ> primosConCubos2 !! 7
--    397
--    
--    (0.01 secs, 648,840 bytes)
--    λ> primosConCubos2 !! (10^3)
--    65580901
--    (0.53 secs, 1,726,837,688 bytes)
--    λ> primosConCubos3 !! (10^3)
--    65580901
--    (0.49 secs, 1,724,258,632 bytes)
--    λ> primosConCubos4 !! (10^3)
--    65580901
--    (0.47 secs, 1,724,833,992 bytes)

-- Demostración de la conjetura
-- ============================

-- Vamos a demostrar que los primos con cubos son diferencias de dos
-- cubos consecutivos.
--
-- Sea p un primo con cubo. Por la definición, existe un entero
-- positivo n tal que n³ + n²p es un cubo. 
--
-- Lema 1: Los números n y p son coprimos (es decir, mcd(n,p) = 1).
-- Dem.: En caso contrario, puesto que p es primo, existe un a tal que 
-- n = ap. Luego n³ + n²p = (a³+a²)p³ es un cubo y, por tanto,
-- a³+a² es un cubo lo que es imposible ya que el siguiente cubo de
-- a³ es a³+3a²+3a+1.
--
-- Lema 2: Los números n² y n+p son coprimos.
-- Dem.: Sea k = mcd(n^2,n+p). Por k divide n², luego k divide a n;
-- además, k divide a n+p y (usando el lema 1 y el ser p primo), se
-- tiene que k = 1.
--
-- Puesto que n³+n²p = n²(n+p) es un cubo, usando el lema 2, se tiene
-- que n² y n+p son cubos y, por serlo n², n también es un cubo. Es
-- decir, existen enteros positivos x e y tales que n = x³ y 
-- n+p = y³. Por tanto, p = y³-x³. Sea k = y-x. Se tiene que k = 1 ya
-- que 
--    p = y³-x³
--      = (n+k)³-n³
--      = 3k+3k²+k³
-- no es primo para k > 1.
--
-- Por consiguiente, p = (x+1)³-x³.

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

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck (NonNegative (NonNegative), maxSize, quickCheckWith, stdArgs)

-- 1ª solución

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

primosConCubos1 :: [Integer]

primosConCubos1 =

[p | x <- [1..],

n <- [1..x],

(x^3 - n^3) `mod` (n^2) == 0,

let p = (x^3 - n^3) `div` (n^2),

isPrime p]

-- 2ª solución

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

-- Para analizar la respuesta, en esta solución se calculan los pares

-- (p,n) tales que p es un primo con cubo y n es un número positivo tal

-- que n^3 + n^2p es un cubo

primosConCubos2' :: [(Integer,Integer)]

primosConCubos2' =

[(p,n) | x <- [1..],

n <- [1..x],

(x^3 - n^3) `mod` (n^2) == 0,

let p = (x^3 - n^3) `div` (n^2),

isPrime p]

-- El cálculo es

-- λ> take 7 primosConCubos2

-- [(7,1),(19,8),(37,27),(61,64),(127,216),(271,729),(331,1000)]

-- Se observa que la sucesión de los segundos elementos [1,8,27,64,...]

-- es la de los cubos y que los primeros elementos se obtienen restando

-- los segundos elementos consecutivos; es decir,

-- 7 = 8 - 1 = 2^3 - 1^3

-- 19 = 27 - 8 = 3^3 - 2^3

-- 37 = 64 - 27 = 4^3 - 3^3

-- Continuando el patrón,

-- 61 = 5^3 - 4^3 = 125 - 64

-- 91 = 6^3 - 5^3 = 216 - 125

-- 127 = 7^3 - 6^3 = 343 - 216

-- 271 = 10^3 - 9^3 = 1000 - 729

-- 331 = 11^3 -10^3 = 1331 - 1000

-- Por tanto, los primos con cubos son diferencias de dos cubos

-- consecutivos; es decir, coinciden con los números cubanos del

-- ejercicio anterior. A partir de la conjetura se obtienen las

-- siguientes definiciones

-- Basado en las anteriores observaciones se obtiene la siguiente

-- definición

primosConCubos2 :: [Integer]

primosConCubos2 =

filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

-- 3ª definición

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

primosConCubos3 :: [Integer]

primosConCubos3 =

filter isPrime diferenciasCubosConsecutivos

diferenciasCubosConsecutivos :: [Integer]

diferenciasCubosConsecutivos =

zipWith (-) (tail cubos) cubos

cubos :: [Integer]

cubos = map (^3) [0..]

-- 4ª solución

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

-- Simplificando la expresión (x+1)^3 - x^3 se obtiene 3*x^2+ 3*x + 1,

-- con lo que la 3ª definición se reduce a

primosConCubos4 :: [Integer]

primosConCubos4 =

[p | x <- [1..],

let p = 3*x^2+ 3*x + 1,

isPrime p]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosConCubos :: NonNegative Int -> Bool

prop_primosConCubos (NonNegative n) =

all (== primosConCubos1 !! n)

[primosConCubos2 !! n,

primosConCubos3 !! n,

primosConCubos4 !! n]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosConCubos1 !! 6

-- 331

-- (1.15 secs, 1,105,612,336 bytes)

-- λ> primosConCubos1 !! 7

-- 397

-- (1.96 secs, 1,909,369,592 bytes)

-- λ> primosConCubos2 !! 7

-- 397

-- (0.01 secs, 648,840 bytes)

-- λ> primosConCubos2 !! (10^3)

-- 65580901

-- (0.53 secs, 1,726,837,688 bytes)

-- λ> primosConCubos3 !! (10^3)

-- 65580901

-- (0.49 secs, 1,724,258,632 bytes)

-- λ> primosConCubos4 !! (10^3)

-- 65580901

-- (0.47 secs, 1,724,833,992 bytes)

-- Demostración de la conjetura

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

-- Vamos a demostrar que los primos con cubos son diferencias de dos

-- cubos consecutivos.

-- Sea p un primo con cubo. Por la definición, existe un entero

-- positivo n tal que n³ + n²p es un cubo.

-- Lema 1: Los números n y p son coprimos (es decir, mcd(n,p) = 1).

-- Dem.: En caso contrario, puesto que p es primo, existe un a tal que

-- n = ap. Luego n³ + n²p = (a³+a²)p³ es un cubo y, por tanto,

-- a³+a² es un cubo lo que es imposible ya que el siguiente cubo de

-- a³ es a³+3a²+3a+1.

-- Lema 2: Los números n² y n+p son coprimos.

-- Dem.: Sea k = mcd(n^2,n+p). Por k divide n², luego k divide a n;

-- además, k divide a n+p y (usando el lema 1 y el ser p primo), se

-- tiene que k = 1.

-- Puesto que n³+n²p = n²(n+p) es un cubo, usando el lema 2, se tiene

-- que n² y n+p son cubos y, por serlo n², n también es un cubo. Es

-- decir, existen enteros positivos x e y tales que n = x³ y

-- n+p = y³. Por tanto, p = y³-x³. Sea k = y-x. Se tiene que k = 1 ya

-- que

-- p = y³-x³

-- = (n+k)³-n³

-- = 3k+3k²+k³

-- no es primo para k > 1.

-- Por consiguiente, p = (x+1)³-x³.

El código se encuentra en GitHub.

Primos cubanos

PorJosé A. Alonso 15-julio-20229-julio-2022

Un primo cubano es un número primo que se puede escribir como diferencia de dos cubos consecutivos. Por ejemplo, el 61 es un primo cubano porque es primo y 61 = 5³-4³.

Definir la sucesión

   cubanos :: [Integer]

1	cubanos :: [Integer]

tal que sus elementos son los números cubanos. Por ejemplo,

   λ> take 15 cubanos
   [7,19,37,61,127,271,331,397,547,631,919,1657,1801,1951,2269]

1 2	λ> take 15 cubanos [7,19,37,61,127,271,331,397,547,631,919,1657,1801,1951,2269]

Soluciones

import Data.Numbers.Primes (isPrime)
import Test.QuickCheck 

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

cubanos1 :: [Integer]
cubanos1 = filter isPrime (zipWith (-) (tail cubos) cubos) 

-- cubos es la lista de los cubos. Por ejemplo,
--    λ> take 10 cubos
--    [1,8,27,64,125,216,343,512,729,1000]
cubos :: [Integer]
cubos = map (^3) [1..]

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

cubanos2 :: [Integer]
cubanos2 = filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

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

cubanos3 :: [Integer]
cubanos3 = filter isPrime [3*x^2 + 3*x + 1 | x <- [1..]]

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

-- La propiedad es
prop_cubanos :: NonNegative Int -> Bool
prop_cubanos (NonNegative n) =
  all (== cubanos1 !! n)
      [cubanos2 !! n,
       cubanos3 !! n]

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

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

-- La comparación es
--    λ> cubanos1 !! 3000
--    795066361
--    (4.21 secs, 16,953,612,192 bytes)
--    λ> cubanos2 !! 3000
--    795066361
--    (4.27 secs, 16,962,597,288 bytes)
--    λ> cubanos3 !! 3000
--    795066361
--    (4.29 secs, 16,956,085,672 bytes)

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck

-- 1ª solución

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

cubanos1 :: [Integer]

cubanos1 = filter isPrime (zipWith (-) (tail cubos) cubos)

-- cubos es la lista de los cubos. Por ejemplo,

-- λ> take 10 cubos

-- [1,8,27,64,125,216,343,512,729,1000]

cubos :: [Integer]

cubos = map (^3) [1..]

-- 2ª solución

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

cubanos2 :: [Integer]

cubanos2 = filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

-- 3ª solución

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

cubanos3 :: [Integer]

cubanos3 = filter isPrime [3*x^2 + 3*x + 1 | x <- [1..]]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_cubanos :: NonNegative Int -> Bool

prop_cubanos (NonNegative n) =

all (== cubanos1 !! n)

[cubanos2 !! n,

cubanos3 !! n]

-- La comprobación es

-- λ> quickCheck prop_cubanos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> cubanos1 !! 3000

-- 795066361

-- (4.21 secs, 16,953,612,192 bytes)

-- λ> cubanos2 !! 3000

-- 795066361

-- (4.27 secs, 16,962,597,288 bytes)

-- λ> cubanos3 !! 3000

-- 795066361

-- (4.29 secs, 16,956,085,672 bytes)

El código se encuentra en GitHub.

Clausura transitiva de una relación binaria

PorJosé A. Alonso 14-julio-20228-julio-2022

La clausura transitiva de una relación binaria R es la relación transitiva que contiene a R. Se puede calcular
usando la composición de relaciones. Veamos un ejemplo, en el que (R ∘ S) representa la composición de R y S: sea

   R = [(1,2),(2,5),(5,6)]

1	R = [(1,2),(2,5),(5,6)]

la relación R no es transitiva ya que (1,2) y (1,5) pertenecen a R pero (1,5) no pertenece; sea

   R1 = R ∪ (R ∘ R)
      = [(1,2),(2,5),(5,6),(1,5),(2,6)]

1 2	R1 = R ∪ (R ∘ R) = [(1,2),(2,5),(5,6),(1,5),(2,6)]

la relación R1 tampoco es transitiva ya que (1,2) y (2,6) pertenecen a R pero (1,6) no pertenece; sea

   R2 = R1 ∪ (R1 ∘ R1)
      = [(1,2),(2,5),(5,6),(1,5),(2,6),(1,6)]

1 2	R2 = R1 ∪ (R1 ∘ R1) = [(1,2),(2,5),(5,6),(1,5),(2,6),(1,6)]

La relación R2 es transitiva y contiene a R. Además, R2 es la clausura transitiva de R.

Definir la función

   clausuraTransitiva :: Ord a => [(a,a)] -> [(a,a)]

1	clausuraTransitiva :: Ord a => [(a,a)] -> [(a,a)]

tal que (clausuraTransitiva r) es la clausura transitiva de r; es decir, la menor relación transitiva que contiene a r. Por ejemplo,

   λ> clausuraTransitiva [(1,2),(2,5),(5,6)]
   [(1,2),(2,5),(5,6),(1,5),(2,6),(1,6)]
   λ> clausuraTransitiva [(1,2),(2,5),(5,6),(6,3)]
   [(1,2),(2,5),(5,6),(6,3),(1,5),(2,6),(5,3),(1,6),(2,3),(1,3)]
   λ> length (clausuraTransitiva [(n,n+1) | n <- [1..100]])
   5050

λ> clausuraTransitiva [(1,2),(2,5),(5,6)]

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

λ> clausuraTransitiva [(1,2),(2,5),(5,6),(6,3)]

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

λ> length (clausuraTransitiva [(n,n+1) | n <- [1..100]])

5050

Soluciones

import Data.List (union, nub, sort)
import Data.Maybe (mapMaybe)
import qualified Data.Map as M (Map, assocs, empty, insertWith, lookup, map)
import Test.QuickCheck (quickCheck)

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

clausuraTransitiva1 :: Ord a => [(a,a)] -> [(a,a)]  
clausuraTransitiva1 r
  | transitiva r = r
  | otherwise    = clausuraTransitiva1 r1
  where r1 = r `union` composicion r r

-- (transitiva r) se verifica si la relación r es transitiva. Por
-- ejemplo, 
--    transitiva [(1,1),(1,3),(3,1),(3,3),(5,5)]  ==  True
--    transitiva [(1,1),(1,3),(3,1),(5,5)]        ==  False
transitiva :: Ord a => [(a,a)] -> Bool
transitiva r = subconjunto (composicion r r) r

-- (composicion r s) es la composición de las relaciones binarias r y
-- s. Por ejemplo, 
--    λ> composicion [(1,2)] [(2,3),(2,4)]
--    [(1,3),(1,4)]
--    λ> composicion [(1,2),(5,2)] [(2,3),(2,4)]
--    [(1,3),(1,4),(5,3),(5,4)]
--    λ> composicion [(1,2),(1,4),(1,5)] [(2,3),(4,3)]
--    [(1,3)]
composicion :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion r s = nub [(x,y) | (x,u) <- r, (v,y) <- s, u == v] 

-- (subconjunto xs ys) se verifica si xs es un subconjunto de xs. Por
-- ejemplo, 
--    subconjunto [1,3] [3,1,5]  ==  True
--    subconjunto [3,1,5] [1,3]  ==  False
subconjunto :: Ord a => [a] -> [a] -> Bool
subconjunto xs ys = all (`elem` ys) xs

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

clausuraTransitiva2 :: Ord a => [(a,a)] -> [(a,a)]  
clausuraTransitiva2 r
  | r1 == r   = r
  | otherwise = clausuraTransitiva2 r1
  where r1 = r `union` composicion r r

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

-- La propiedad es
prop_clausuraTransitiva :: [(Int,Int)] -> Bool
prop_clausuraTransitiva r =
  all (== sort (clausuraTransitiva1 r))
      [sort (clausuraTransitiva2 r),
       sort (clausuraTransitiva3 r)]
  
-- La comprobación es
--    λ> quickCheck prop_clausuraTransitiva
--    +++ OK, passed 100 tests.

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

clausuraTransitiva3 :: Ord a => [(a,a)] -> [(a,a)]  
clausuraTransitiva3 r
  | transitiva3 r = r
  | otherwise     = clausuraTransitiva3 r1
  where r1 = r `union` composicion3 r r

transitiva3 :: Ord a => [(a,a)] -> Bool
transitiva3 r = subconjunto (composicion3 r r) r

composicion3 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion3 r s =
  relAlista (composicionRel (listaArel r) (listaArel s))

-- Una relación se puede representar por un diccionario donde las claves
-- son los elementos y los valores son las listas de los elementos con
-- los que se relaciona.
type Rel a = M.Map a [a]

-- (listaArel xys) es la relación correspondiente a la lista de pares
-- xys. Por ejemplo.
--    λ> listaArel [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]
--    fromList [(1,[1,2,3,6]),(2,[2,6]),(3,[3,6]),(6,[6])]
listaArel :: Ord a => [(a,a)] -> Rel a
listaArel []          = M.empty
listaArel ((x,y):xys) = M.insertWith (++) x [y] (listaArel xys)

-- (composicionRel r s) es la composición de las relaciones r y s. Por
-- ejemplo,
--    λ> r = listaArel [(1,2),(5,2)]
--    λ> s = listaArel [(2,3),(2,4)]
--    λ> composicionRel r s
--    fromList [(1,[3,4]),(5,[3,4])]
composicionRel :: Ord a => Rel a -> Rel a -> Rel a
composicionRel r s =
  M.map f r 
  where f xs = concat (mapMaybe (`M.lookup` s) xs)

-- (relAlista r) es la lista de pares correspondientes a la relación
-- r. Por ejemplo,
--    λ> relAlista (M.fromList [(1,[3,4]),(5,[3,4])])
--    [(1,3),(1,4),(5,3),(5,4)]
relAlista :: Ord a => Rel a -> [(a,a)]
relAlista r =
  nub [(x,y) | (x,ys) <- M.assocs r, y <- ys] 

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

-- La comparación es
--    λ> length (clausuraTransitiva1 [(n,n+1) | n <- [1..60]])
--    1830
--    (2.15 secs, 453,533,992 bytes)
--    λ> length (clausuraTransitiva2 [(n,n+1) | n <- [1..60]])
--    1830
--    (2.23 secs, 558,571,904 bytes)
--    λ> length (clausuraTransitiva3 [(n,n+1) | n <- [1..60]])
--    1830
--    (0.25 secs, 207,168,552 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

import Data.List (union, nub, sort)

import Data.Maybe (mapMaybe)

import qualified Data.Map as M (Map, assocs, empty, insertWith, lookup, map)

import Test.QuickCheck (quickCheck)

-- 1ª solución

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

clausuraTransitiva1 :: Ord a => [(a,a)] -> [(a,a)]

clausuraTransitiva1 r

| transitiva r = r

| otherwise = clausuraTransitiva1 r1

where r1 = r `union` composicion r r

-- (transitiva r) se verifica si la relación r es transitiva. Por

-- ejemplo,

-- transitiva [(1,1),(1,3),(3,1),(3,3),(5,5)] == True

-- transitiva [(1,1),(1,3),(3,1),(5,5)] == False

transitiva :: Ord a => [(a,a)] -> Bool

transitiva r = subconjunto (composicion r r) r

-- (composicion r s) es la composición de las relaciones binarias r y

-- s. Por ejemplo,

-- λ> composicion [(1,2)] [(2,3),(2,4)]

-- [(1,3),(1,4)]

-- λ> composicion [(1,2),(5,2)] [(2,3),(2,4)]

-- [(1,3),(1,4),(5,3),(5,4)]

-- λ> composicion [(1,2),(1,4),(1,5)] [(2,3),(4,3)]

-- [(1,3)]

composicion :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion r s = nub [(x,y) | (x,u) <- r, (v,y) <- s, u == v]

-- (subconjunto xs ys) se verifica si xs es un subconjunto de xs. Por

-- ejemplo,

-- subconjunto [1,3] [3,1,5] == True

-- subconjunto [3,1,5] [1,3] == False

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

subconjunto xs ys = all (`elem` ys) xs

-- 2ª solución

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

clausuraTransitiva2 :: Ord a => [(a,a)] -> [(a,a)]

clausuraTransitiva2 r

| r1 == r = r

| otherwise = clausuraTransitiva2 r1

where r1 = r `union` composicion r r

-- Comprobación de equivalencia

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

-- La propiedad es

prop_clausuraTransitiva :: [(Int,Int)] -> Bool

prop_clausuraTransitiva r =

all (== sort (clausuraTransitiva1 r))

[sort (clausuraTransitiva2 r),

sort (clausuraTransitiva3 r)]

-- La comprobación es

-- λ> quickCheck prop_clausuraTransitiva

-- +++ OK, passed 100 tests.

-- 3ª solución

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

clausuraTransitiva3 :: Ord a => [(a,a)] -> [(a,a)]

clausuraTransitiva3 r

| transitiva3 r = r

| otherwise = clausuraTransitiva3 r1

where r1 = r `union` composicion3 r r

transitiva3 :: Ord a => [(a,a)] -> Bool

transitiva3 r = subconjunto (composicion3 r r) r

composicion3 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion3 r s =

relAlista (composicionRel (listaArel r) (listaArel s))

-- Una relación se puede representar por un diccionario donde las claves

-- son los elementos y los valores son las listas de los elementos con

-- los que se relaciona.

type Rel a = M.Map a [a]

-- (listaArel xys) es la relación correspondiente a la lista de pares

-- xys. Por ejemplo.

-- λ> listaArel [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]

-- fromList [(1,[1,2,3,6]),(2,[2,6]),(3,[3,6]),(6,[6])]

listaArel :: Ord a => [(a,a)] -> Rel a

listaArel [] = M.empty

listaArel ((x,y):xys) = M.insertWith (++) x [y] (listaArel xys)

-- (composicionRel r s) es la composición de las relaciones r y s. Por

-- ejemplo,

-- λ> r = listaArel [(1,2),(5,2)]

-- λ> s = listaArel [(2,3),(2,4)]

-- λ> composicionRel r s

-- fromList [(1,[3,4]),(5,[3,4])]

composicionRel :: Ord a => Rel a -> Rel a -> Rel a

composicionRel r s =

M.map f r

where f xs = concat (mapMaybe (`M.lookup` s) xs)

-- (relAlista r) es la lista de pares correspondientes a la relación

-- r. Por ejemplo,

-- λ> relAlista (M.fromList [(1,[3,4]),(5,[3,4])])

-- [(1,3),(1,4),(5,3),(5,4)]

relAlista :: Ord a => Rel a -> [(a,a)]

relAlista r =

nub [(x,y) | (x,ys) <- M.assocs r, y <- ys]

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (clausuraTransitiva1 [(n,n+1) | n <- [1..60]])

-- 1830

-- (2.15 secs, 453,533,992 bytes)

-- λ> length (clausuraTransitiva2 [(n,n+1) | n <- [1..60]])

-- 1830

-- (2.23 secs, 558,571,904 bytes)

-- λ> length (clausuraTransitiva3 [(n,n+1) | n <- [1..60]])

-- 1830

-- (0.25 secs, 207,168,552 bytes)

El código se encuentra en GitHub.

Potencias perfectas

Soluciones

Sumas alternas de factoriales

Soluciones

Primos con cubos

Soluciones

Primos cubanos

Soluciones

Clausura transitiva de una relación binaria

Soluciones

Entradas recientes

Buscar

Correo electrónico