gcd

Conjunto de primos relativos

PorJosé A. Alonso 30-marzo-20226-abril-2022

Dos números enteros positivos son primos relativos si no tienen ningún factor primo en común; es decit, si 1 es su único divisor común. Por ejemplo, 6 y 35 son primos entre sí, pero 6 y 27 no lo son porque ambos son divisibles por 3.

Definir la función

   primosRelativos :: [Int] -> Bool

1	primosRelativos :: [Int] -> Bool

tal que (primosRelativos xs) se verifica si los elementos de xs son primos relativos dos a dos. Por ejemplo,

   primosRelativos [6,35]         ==  True
   primosRelativos [6,27]         ==  False
   primosRelativos [2,3,4]        ==  False
   primosRelativos [6,35,11]      ==  True
   primosRelativos [6,35,11,221]  ==  True
   primosRelativos [6,35,11,231]  ==  False

primosRelativos [6,35] == True

primosRelativos [6,27] == False

primosRelativos [2,3,4] == False

primosRelativos [6,35,11] == True

primosRelativos [6,35,11,221] == True

primosRelativos [6,35,11,231] == False

Soluciones

import Test.QuickCheck
import Data.List (delete, intersect)
import Data.Numbers.Primes (primeFactors, primes)
import qualified Data.Set as S (disjoint, fromList)

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

primosRelativos1 :: [Int] -> Bool
primosRelativos1 []     = True
primosRelativos1 (x:xs) =
  and [sonPrimosRelativos1 x y | y <- xs] && primosRelativos1 xs

-- (sonPrimosRelativos x y) se verifica si x e y son primos
-- relativos. Por ejemplo,
--    sonPrimosRelativos1 6 35  ==  True
--    sonPrimosRelativos1 6 27  ==  False
sonPrimosRelativos1 :: Int -> Int -> Bool
sonPrimosRelativos1 x y =
  null (divisoresPrimos x `intersect` divisoresPrimos y)

-- (divisoresPrimos x) es la lista de los divisores primos de x. Por
-- ejemplo,
--    divisoresPrimos 600  ==  [2,2,2,3,5,5]
divisoresPrimos :: Int -> [Int]
divisoresPrimos 1 = []
divisoresPrimos x =
  y : divisoresPrimos (x `div` y)
  where y = menorDivisorPrimo x

-- (menorDivisorPrimo x) es el menor divisor primo de x. Por ejemplo,
--    menorDivisorPrimo 15  ==  3
--    menorDivisorPrimo 11  ==  11
menorDivisorPrimo :: Int -> Int
menorDivisorPrimo x =
  head [y | y <- [2..], x `mod` y == 0]

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

primosRelativos2 :: [Int] -> Bool
primosRelativos2 []     = True
primosRelativos2 (x:xs) =
  all (sonPrimosRelativos1 x) xs && primosRelativos2 xs

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

primosRelativos3 :: [Int] -> Bool
primosRelativos3 []     = True
primosRelativos3 (x:xs) =
  all (sonPrimosRelativos2 x) xs && primosRelativos3 xs

sonPrimosRelativos2 :: Int -> Int -> Bool
sonPrimosRelativos2 x y =
  null (primeFactors x `intersect` primeFactors y)

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

primosRelativos4 :: [Int] -> Bool
primosRelativos4 []     = True
primosRelativos4 (x:xs) =
  all (sonPrimosRelativos3 x) xs && primosRelativos4 xs

sonPrimosRelativos3 :: Int -> Int -> Bool
sonPrimosRelativos3 x y =
  S.fromList (primeFactors x) `S.disjoint` S.fromList (primeFactors y)

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

primosRelativos5 :: [Int] -> Bool
primosRelativos5 []     = True
primosRelativos5 (x:xs) =
  all (sonPrimosRelativos5 x) xs && primosRelativos5 xs

sonPrimosRelativos5 :: Int -> Int -> Bool
sonPrimosRelativos5 x y =
  gcd x y == 1

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

-- La propiedad es
prop_primosRelativos :: [Positive Int] -> Bool
prop_primosRelativos xs =
  all (== primosRelativos1 ys)
      [primosRelativos2 ys,
       primosRelativos3 ys,
       primosRelativos4 ys,
       primosRelativos5 ys]
  where ys = getPositive <$> xs

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

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

-- La comparación es
--    λ> primosRelativos1 (take 120 primes)
--    True
--    (1.92 secs, 869,909,416 bytes)
--    λ> primosRelativos2 (take 120 primes)
--    True
--    (1.99 secs, 869,045,656 bytes)
--    λ> primosRelativos3 (take 120 primes)
--    True
--    (0.09 secs, 221,183,200 bytes)
--
--    λ> primosRelativos3 (take 600 primes)
--    True
--    (2.62 secs, 11,196,690,856 bytes)
--    λ> primosRelativos4 (take 600 primes)
--    True
--    (2.66 secs, 11,190,940,456 bytes)
--    λ> primosRelativos5 (take 600 primes)
--    True
--    (0.14 secs, 123,673,648 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

import Test.QuickCheck

import Data.List (delete, intersect)

import Data.Numbers.Primes (primeFactors, primes)

import qualified Data.Set as S (disjoint, fromList)

-- 1ª solución

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

primosRelativos1 :: [Int] -> Bool

primosRelativos1 [] = True

primosRelativos1 (x:xs) =

and [sonPrimosRelativos1 x y | y <- xs] && primosRelativos1 xs

-- (sonPrimosRelativos x y) se verifica si x e y son primos

-- relativos. Por ejemplo,

-- sonPrimosRelativos1 6 35 == True

-- sonPrimosRelativos1 6 27 == False

sonPrimosRelativos1 :: Int -> Int -> Bool

sonPrimosRelativos1 x y =

null (divisoresPrimos x `intersect` divisoresPrimos y)

-- (divisoresPrimos x) es la lista de los divisores primos de x. Por

-- ejemplo,

-- divisoresPrimos 600 == [2,2,2,3,5,5]

divisoresPrimos :: Int -> [Int]

divisoresPrimos 1 = []

divisoresPrimos x =

y : divisoresPrimos (x `div` y)

where y = menorDivisorPrimo x

-- (menorDivisorPrimo x) es el menor divisor primo de x. Por ejemplo,

-- menorDivisorPrimo 15 == 3

-- menorDivisorPrimo 11 == 11

menorDivisorPrimo :: Int -> Int

menorDivisorPrimo x =

head [y | y <- [2..], x `mod` y == 0]

-- 2ª solución

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

primosRelativos2 :: [Int] -> Bool

primosRelativos2 [] = True

primosRelativos2 (x:xs) =

all (sonPrimosRelativos1 x) xs && primosRelativos2 xs

-- 3ª solución

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

primosRelativos3 :: [Int] -> Bool

primosRelativos3 [] = True

primosRelativos3 (x:xs) =

all (sonPrimosRelativos2 x) xs && primosRelativos3 xs

sonPrimosRelativos2 :: Int -> Int -> Bool

sonPrimosRelativos2 x y =

null (primeFactors x `intersect` primeFactors y)

-- 4ª solución

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

primosRelativos4 :: [Int] -> Bool

primosRelativos4 [] = True

primosRelativos4 (x:xs) =

all (sonPrimosRelativos3 x) xs && primosRelativos4 xs

sonPrimosRelativos3 :: Int -> Int -> Bool

sonPrimosRelativos3 x y =

S.fromList (primeFactors x) `S.disjoint` S.fromList (primeFactors y)

-- 5ª solución

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

primosRelativos5 :: [Int] -> Bool

primosRelativos5 [] = True

primosRelativos5 (x:xs) =

all (sonPrimosRelativos5 x) xs && primosRelativos5 xs

sonPrimosRelativos5 :: Int -> Int -> Bool

sonPrimosRelativos5 x y =

gcd x y == 1

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosRelativos :: [Positive Int] -> Bool

prop_primosRelativos xs =

all (== primosRelativos1 ys)

[primosRelativos2 ys,

primosRelativos3 ys,

primosRelativos4 ys,

primosRelativos5 ys]

where ys = getPositive <$> xs

-- La comprobación es

-- λ> quickCheck prop_primosRelativos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosRelativos1 (take 120 primes)

-- True

-- (1.92 secs, 869,909,416 bytes)

-- λ> primosRelativos2 (take 120 primes)

-- True

-- (1.99 secs, 869,045,656 bytes)

-- λ> primosRelativos3 (take 120 primes)

-- True

-- (0.09 secs, 221,183,200 bytes)

-- λ> primosRelativos3 (take 600 primes)

-- True

-- (2.62 secs, 11,196,690,856 bytes)

-- λ> primosRelativos4 (take 600 primes)

-- True

-- (2.66 secs, 11,190,940,456 bytes)

-- λ> primosRelativos5 (take 600 primes)

-- True

-- (0.14 secs, 123,673,648 bytes)

El código se encuentra en GitHub.

La elaboración de las soluciones se describe en el siguiente vídeo

Nuevas soluciones

En los comentarios se pueden escribir nuevas soluciones.
El código se debe escribir entre una línea con <pre lang="haskell"> y otra con </pre>

La sucesión ECG

PorJosé A. Alonso 1-junio-20202-enero-2021

La sucesión ECG estás definida por a(1) = 1, a(2) = 2 y, para n >= 3, a(n) es el menor natural que aún no está en la sucesión tal que a(n) tiene algún divisor común con a(n-1).

Los primeros términos de la sucesión son 1, 2, 4, 6, 3, 9, 12, 8, 10, 5, 15, …

Al dibujar su gráfica, se parece a la de los electrocardiogramas (abreviadamente, ECG). Por ello, la sucesión se conoce como la sucesión ECG.

Definir las funciones

   sucECG :: [Integer]
   graficaSucECG :: Int -> IO ()

1 2	sucECG :: [Integer] graficaSucECG :: Int -> IO ()

tales que

sucECG es la lista de los términos de la sucesión ECG. Por ejemplo,

     λ> take 20 sucECG
     [1,2,4,6,3,9,12,8,10,5,15,18,14,7,21,24,16,20,22,11]
     λ> sucECG !! 6000
     6237

λ> take 20 sucECG

[1,2,4,6,3,9,12,8,10,5,15,18,14,7,21,24,16,20,22,11]

λ> sucECG !! 6000

6237

(graficaSucECG n) dibuja la gráfica de los n primeros términos de la sucesión ECG. Por ejemplo, (graficaSucECG 160) dibuja

Soluciones

import Data.List (delete)
import Graphics.Gnuplot.Simple

sucECG :: [Integer]
sucECG = 1 : ecg 2 [2..]
  where ecg x zs = f zs
          where f (y:ys) | gcd x y > 1 = y : ecg y (delete y zs)
                         | otherwise   = f ys

graficaSucECG :: Int -> IO ()
graficaSucECG n =
  plotList [ Key Nothing
           , PNG "La_sucesion_ECG.png" 
           ]
           (take n sucECG)

import Data.List (delete)

import Graphics.Gnuplot.Simple

sucECG :: [Integer]

sucECG = 1 : ecg 2 [2..]

where ecg x zs = f zs

where f (y:ys) | gcd x y > 1 = y : ecg y (delete y zs)

| otherwise = f ys

graficaSucECG :: Int -> IO ()

graficaSucECG n =

plotList [ Key Nothing

, PNG "La_sucesion_ECG.png"

]

(take n sucECG)

Huecos de Aquiles

PorJosé A. Alonso 23-abril-202026-marzo-2022

Un número de Aquiles es un número natural n que es potente (es decir, si p es un divisor primo de n, entonces p² también lo es) y no es una potencia perfecta (es decir, no existen números naturales m y k tales que n es igual a m^k). Por ejemplo,

108 es un número de Aquiles proque es un número potente (ya que su factorización es 2^2 · 3^3, sus divisores primos son 2 and 3 y sus cuadrados (2^2 = 4 y 3^2 = 9) son divisores de 108. Además, 108 no es una potencia perfecta.
360 no es un número de Aquiles ya que 5 es un divisor primo de 360, pero 5^2 = 15 no lo es.
784 no es un número de Aquiles porque, aunque es potente, es una potencia perfecta ya que 784 = 28^2.

Los primeros números de Aquiles son

   72, 108, 200, 288, 392, 432, 500, 648, 675, 800, 864, 968, 972, ...

1	72, 108, 200, 288, 392, 432, 500, 648, 675, 800, 864, 968, 972, ...

Definir las funciones

   esAquiles              :: Integer -> Bool
   huecosDeAquiles        :: [Integer]
   graficaHuecosDeAquiles :: Int -> IO ()

esAquiles :: Integer -> Bool

huecosDeAquiles :: [Integer]

graficaHuecosDeAquiles :: Int -> IO ()

tales que

(esAquiles x) se verifica si x es un número de Aquiles. Por ejemplo,

     esAquiles 108         ==  True
     esAquiles 360         ==  False
     esAquiles 784         ==  False
     esAquiles 5425069447  ==  True
     esAquiles 5425069448  ==  True

esAquiles 108 == True

esAquiles 360 == False

esAquiles 784 == False

esAquiles 5425069447 == True

esAquiles 5425069448 == True

huecosDeAquiles es la sucesión de la diferencias entre los números de Aquiles consecutivos. Por ejemplo,

     λ> take 15 huecosDeAquiles
     [36,92,88,104,40,68,148,27,125,64,104,4,153,27,171]

1 2	λ> take 15 huecosDeAquiles [36,92,88,104,40,68,148,27,125,64,104,4,153,27,171]

(graficaHuecosDeAquiles n) dibuja la gráfica de los n primeros huecos de Aquiles. Por ejemplo, (graficaHuecosDeAquiles 160) dibuja

Soluciones

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

-- Definición de esAquiles
-- =======================

esAquiles :: Integer -> Bool
esAquiles x = esPotente x && noEsPotenciaPerfecta x

-- (esPotente x) se verifica si x es potente. Por ejemplo,
--    esPotente 108  ==  True
--    esPotente 360  ==  False
--    esPotente 784  ==  True
esPotente :: Integer -> Bool
esPotente x = all (>1) (exponentes x)

-- (exponentes x) es la lista de los exponentes en la factorización de
-- x. Por ejemplo,
--    exponentes 108  ==  [2,3]
--    exponentes 360  ==  [3,2,1]
--    exponentes 784  ==  [4,2]
exponentes :: Integer -> [Int]
exponentes x = map length (group (primeFactors x))

-- (noEsPotenciaPerfecta x) se verifica si x no es una potencia
-- perfecta. Por ejemplo,
--    noEsPotenciaPerfecta 108  ==  True
--    noEsPotenciaPerfecta 360  ==  True
--    noEsPotenciaPerfecta 784  ==  False
noEsPotenciaPerfecta :: Integer -> Bool
noEsPotenciaPerfecta x = foldl1 gcd (exponentes x) == 1 

-- Definición de huecosDeAquiles
-- =============================

huecosDeAquiles :: [Integer]
huecosDeAquiles = zipWith (-) (tail aquiles) aquiles

-- aquiles es la sucesión de los números de Aquiles. Por ejemplo, 
--    λ> take 15 aquiles
--    [72,108,200,288,392,432,500,648,675,800,864,968,972,1125,1152]
aquiles :: [Integer]
aquiles = filter esAquiles [2..]

-- Definición de graficaHuecosDeAquiles
-- ====================================

graficaHuecosDeAquiles :: Int -> IO ()
graficaHuecosDeAquiles n =
  plotList [ Key Nothing
           , PNG "Huecos_de_Aquiles.png"
           ]
           (take n huecosDeAquiles)

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

import Graphics.Gnuplot.Simple

-- Definición de esAquiles

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

esAquiles :: Integer -> Bool

esAquiles x = esPotente x && noEsPotenciaPerfecta x

-- (esPotente x) se verifica si x es potente. Por ejemplo,

-- esPotente 108 == True

-- esPotente 360 == False

-- esPotente 784 == True

esPotente :: Integer -> Bool

esPotente x = all (>1) (exponentes x)

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

-- x. Por ejemplo,

-- exponentes 108 == [2,3]

-- exponentes 360 == [3,2,1]

-- exponentes 784 == [4,2]

exponentes :: Integer -> [Int]

exponentes x = map length (group (primeFactors x))

-- (noEsPotenciaPerfecta x) se verifica si x no es una potencia

-- perfecta. Por ejemplo,

-- noEsPotenciaPerfecta 108 == True

-- noEsPotenciaPerfecta 360 == True

-- noEsPotenciaPerfecta 784 == False

noEsPotenciaPerfecta :: Integer -> Bool

noEsPotenciaPerfecta x = foldl1 gcd (exponentes x) == 1

-- Definición de huecosDeAquiles

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

huecosDeAquiles :: [Integer]

huecosDeAquiles = zipWith (-) (tail aquiles) aquiles

-- aquiles es la sucesión de los números de Aquiles. Por ejemplo,

-- λ> take 15 aquiles

-- [72,108,200,288,392,432,500,648,675,800,864,968,972,1125,1152]

aquiles :: [Integer]

aquiles = filter esAquiles [2..]

-- Definición de graficaHuecosDeAquiles

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

graficaHuecosDeAquiles :: Int -> IO ()

graficaHuecosDeAquiles n =

plotList [ Key Nothing

, PNG "Huecos_de_Aquiles.png"

]

(take n huecosDeAquiles)

Otras soluciones

Se pueden escribir otras soluciones en los comentarios.
El código se debe escribir entre una línea con <pre lang="haskell"> y otra con </pre>

Acotación del primorial

PorJosé A. Alonso 31-enero-20207-febrero-2020

El primorial de un número natural n es el producto de todos los números primos menores o iguales a n. Por ejemplo, el primorial de 5 es 30 porque el producto de los primos menores o iguales que 5 es

   2 * 3 * 5 = 30

1	2 * 3 * 5 = 30

La propiedad de Erdös de acotación de los primoriales afirma que

Para todo número natural n, su primorial es menor o igual que 4ⁿ.

Definir las funciones

   primorial :: Integer -> Integer
   primoriales :: [Integer]

1 2	primorial :: Integer -> Integer primoriales :: [Integer]

tales que

(primorial n) es el primorial de n. Por ejemplo,

     primorial 3  ==  6
     primorial 5  ==  30
     primorial 8  ==  210

primorial 3 == 6

primorial 5 == 30

primorial 8 == 210

primoriales es la sucesión de los primoriales. Por ejemplo,

   λ> take 15 primoriales
   [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

1 2	λ> take 15 primoriales [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

Comprobar con QuickCheck la propiedad de Erdös de acotación de los primoriales.

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck

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

primorial :: Integer -> Integer
primorial n = product (takeWhile (<= n) primes)

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

primorial2 :: Integer -> Integer
primorial2 0 = 1
primorial2 n | gcd n x == 1 = n*x
             | otherwise    = x
  where x = primorial2 (n-1)

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

--    λ> length (show (primorial (5*10^5)))
--    216852
--    (1.65 secs, 2,472,977,584 bytes)
--    λ> length (show (primorial2 (5*10^5)))
--    216852
--    (3.56 secs, 2,719,162,272 bytes)

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

--    λ> take 15 primoriales
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales :: [Integer]
primoriales = map primorial [0..]

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

--    λ> take 15 primoriales2
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales2 :: [Integer]
primoriales2 = map primorial2 [0..]

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

--    λ> take 15 primoriales3
--    [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]
primoriales3 :: [Integer]
primoriales3 = scanl1 f [1..]
  where f x n | gcd n x == 1 = n*x
              | otherwise    = x

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

--    λ> minimum (take 5000 primoriales)
--    1
--    (1.56 secs, 4,857,760,464 bytes)
--    λ> minimum (take 5000 primoriales2)
--    1
--    (9.39 secs, 10,942,848,240 bytes)
--    λ> minimum (take 5000 primoriales3)
--    1
--    (0.01 secs, 5,575,024 bytes)
--    
--    λ> minimum (take 6000 primoriales)
--    1
--    (2.22 secs, 7,013,937,248 bytes)
--    λ> minimum (take 6000 primoriales3)
--    1
--    (0.01 secs, 6,737,328 bytes)

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

prop_primorial :: Integer -> Property
prop_primorial n =
  n >= 0 ==> primorial n <= 4^n

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

import Data.Numbers.Primes

import Test.QuickCheck

-- 1ª definición de primorial

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

primorial :: Integer -> Integer

primorial n = product (takeWhile (<= n) primes)

-- 2ª definición de primorial

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

primorial2 :: Integer -> Integer

primorial2 0 = 1

primorial2 n | gcd n x == 1 = n*x

| otherwise = x

where x = primorial2 (n-1)

-- Comparación de eficiencia

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

-- λ> length (show (primorial (5*10^5)))

-- 216852

-- (1.65 secs, 2,472,977,584 bytes)

-- λ> length (show (primorial2 (5*10^5)))

-- 216852

-- (3.56 secs, 2,719,162,272 bytes)

-- 1ª definición de primoriales

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

-- λ> take 15 primoriales

-- [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

primoriales :: [Integer]

primoriales = map primorial [0..]

-- 2ª definición de primoriales

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

-- λ> take 15 primoriales2

-- [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

primoriales2 :: [Integer]

primoriales2 = map primorial2 [0..]

-- 3ª definición de primoriales

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

-- λ> take 15 primoriales3

-- [1,1,2,6,6,30,30,210,210,210,210,2310,2310,30030,30030]

primoriales3 :: [Integer]

primoriales3 = scanl1 f [1..]

where f x n | gcd n x == 1 = n*x

| otherwise = x

-- Comparación de eficiencia

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

-- λ> minimum (take 5000 primoriales)

-- 1

-- (1.56 secs, 4,857,760,464 bytes)

-- λ> minimum (take 5000 primoriales2)

-- 1

-- (9.39 secs, 10,942,848,240 bytes)

-- λ> minimum (take 5000 primoriales3)

-- 1

-- (0.01 secs, 5,575,024 bytes)

-- λ> minimum (take 6000 primoriales)

-- 1

-- (2.22 secs, 7,013,937,248 bytes)

-- λ> minimum (take 6000 primoriales3)

-- 1

-- (0.01 secs, 6,737,328 bytes)

-- Propiedad

-- =========

prop_primorial :: Integer -> Property

prop_primorial n =

n >= 0 ==> primorial n <= 4^n

-- La comprobación es

-- λ> quickCheck prop_primorial

-- +++ OK, passed 100 tests.

Otras soluciones

Se pueden escribir otras soluciones en los comentarios.
El código se debe escribir entre una línea con <pre lang=»haskell»> y otra con </pre>

Pensamiento

«Las matemáticas son la reina de las ciencias y la teoría de los números es la reina de las matemáticas.»

Carl Friedrich Gauss.

Las conjeturas de Catalan y de Pillai

PorJosé A. Alonso 15-enero-202022-enero-2020

La conjetura de Catalan, enunciada en 1844 por Eugène Charles Catalan y demostrada 2002 por Preda Mihăilescu1, afirma que

Las únicas dos potencias de números enteros consecutivos son 8 y 9 (que son respectivamente 2³ y 3²).

En otras palabras, la única solución entera de la ecuación

   x^a - y^b = 1

1	x^a - y^b = 1

para x, a, y, b > 1 es x = 3, a = 2, y = 2, b = 3.

La conjetura de Pillai, propuesta por S.S. Pillai en 1942, generaliza este resultado y es un problema abierto. Afirma que cada entero se puede escribir sólo un número finito de veces como una diferencia de dos potencias perfectas. En otras palabras, para todo entero positivo n, el conjunto de soluciones de

   x^a - y^b = n

1	x^a - y^b = n

para x, a, y, b > 1 es finito.

Por ejemplo, para n = 4, hay 3 soluciones

   (2,3, 2,2) ya que 2³ -  2² =   8 -   4 = 4
   (6,2, 2,5) ya que 6² -  2⁵ =  36 -  32 = 4
   (5,3,11,2) ya que 5³ - 11² = 125 - 121 = 4

(2,3, 2,2) ya que 2³ - 2² = 8 - 4 = 4

(6,2, 2,5) ya que 6² - 2⁵ = 36 - 32 = 4

(5,3,11,2) ya que 5³ - 11² = 125 - 121 = 4

Las soluciones se pueden representar por la menor potencia (en el caso anterior, por 4, 32 y 121) ya que dado n (en el caso anterior es 4), la potencia mayor es la menor más n.

Definir las funciones

   potenciasPerfectas :: [Integer]
   solucionesPillati :: Integer -> [Integer]
   solucionesPillatiAcotadas :: Integer -> Integer -> [Integer]

potenciasPerfectas :: [Integer]

solucionesPillati :: Integer -> [Integer]

solucionesPillatiAcotadas :: Integer -> Integer -> [Integer]

tales que

potenciasPerfectas es la lista de las potencias perfectas (es decir, de los números de la forma x^a con x y a mayores que 1). Por ejemplo,

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

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

(solucionesPillati n) es la lista de las menores potencias de las soluciones de la ecuación de Pillati x^a – y^b = n; es decir, es la lista de los u tales que u y u+n son potencias perfectas. Por ejemplo,

     take 3 (solucionesPillati 4)  ==  [4,32,121]
     take 2 (solucionesPillati 5)  ==  [4,27]
     take 4 (solucionesPillati 7)  ==  [9,25,121,32761]

take 3 (solucionesPillati 4) == [4,32,121]

take 2 (solucionesPillati 5) == [4,27]

take 4 (solucionesPillati 7) == [9,25,121,32761]

(solucionesPillatiAcotadas c n) es la lista de elementos de (solucionesPillati n) menores que n. Por ejemplo,

     solucionesPillatiAcotadas (10^3) 1  ==  [8]
     solucionesPillatiAcotadas (10^3) 2  ==  [25]
     solucionesPillatiAcotadas (10^3) 3  ==  [125]
     solucionesPillatiAcotadas (10^3) 4  ==  [4,32,121]
     solucionesPillatiAcotadas (10^3) 5  ==  [4,27]
     solucionesPillatiAcotadas (10^3) 6  ==  []
     solucionesPillatiAcotadas (10^3) 7  ==  [9,25,121]
     solucionesPillatiAcotadas (10^5) 7  ==  [9,25,121,32761]

solucionesPillatiAcotadas (10^3) 1 == [8]

solucionesPillatiAcotadas (10^3) 2 == [25]

solucionesPillatiAcotadas (10^3) 3 == [125]

solucionesPillatiAcotadas (10^3) 4 == [4,32,121]

solucionesPillatiAcotadas (10^3) 5 == [4,27]

solucionesPillatiAcotadas (10^3) 6 == []

solucionesPillatiAcotadas (10^3) 7 == [9,25,121]

solucionesPillatiAcotadas (10^5) 7 == [9,25,121,32761]

Soluciones

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

-- Definiciones de potenciasPerfectas
-- ==================================

-- 1ª definició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ª definició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ª definició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

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

--    λ> potenciasPerfectas1 !! 200
--    28224
--    (7.24 secs, 9,245,991,160 bytes)
--    λ> potenciasPerfectas2 !! 200
--    28224
--    (0.30 secs, 814,597,152 bytes)
--    λ> potenciasPerfectas3 !! 200
--    28224
--    (0.01 secs, 7,061,120 bytes)

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

-- Definición de solucionesPillati
-- ===============================

solucionesPillati :: Integer -> [Integer]
solucionesPillati n =
  [x | x <- potenciasPerfectas
     , esPotenciaPerfecta2 (x+n)]

-- Definición de solucionesPillatiAcotadas
-- =======================================

solucionesPillatiAcotadas :: Integer -> Integer -> [Integer]
solucionesPillatiAcotadas c n =
  [x | x <- takeWhile (< (c-n)) potenciasPerfectas
     , esPotenciaPerfecta2 (x+n)]

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

-- Definiciones de potenciasPerfectas

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

-- 1ª definició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ª definició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ª definició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

-- Comparación de eficiencia

-- -------------------------

-- λ> potenciasPerfectas1 !! 200

-- 28224

-- (7.24 secs, 9,245,991,160 bytes)

-- λ> potenciasPerfectas2 !! 200

-- 28224

-- (0.30 secs, 814,597,152 bytes)

-- λ> potenciasPerfectas3 !! 200

-- 28224

-- (0.01 secs, 7,061,120 bytes)

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

potenciasPerfectas :: [Integer]

potenciasPerfectas = potenciasPerfectas3

-- Definición de solucionesPillati

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

solucionesPillati :: Integer -> [Integer]

solucionesPillati n =

[x | x <- potenciasPerfectas

, esPotenciaPerfecta2 (x+n)]

-- Definición de solucionesPillatiAcotadas

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

solucionesPillatiAcotadas :: Integer -> Integer -> [Integer]

solucionesPillatiAcotadas c n =

[x | x <- takeWhile (< (c-n)) potenciasPerfectas

, esPotenciaPerfecta2 (x+n)]

Referencia

Catalan’s conjecture en Wikipedia.

Pensamiento

Y te enviaré mi canción:
«Se canta lo que se pierde»,
con un papagayo verde
que la diga en tu balcón.

Antonio Machado

Enteros como sumas de tres coprimos.

PorJosé A. Alonso 8-enero-202015-enero-2020

Dos números enteros son coprimos (o primos entre sí) si no tienen ningún factor primo en común. Por ejemplo, 4 y 15 son coprimos.

Una terna coprima es una terna (a,b,c) tal que

a y b son coprimos,
a y c son coprimos y
b y c son coprimos.

Por ejemplo, (3,4,5) es una terna coprima.

Definir la función

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

1	sumas3coprimos :: Integer -> [(Integer,Integer,Integer)]

tal que (sumas3coprimos n) es la lista de las ternas coprimas cuya suma es n. Por ejemplo,

   sumas3coprimos 10  ==  [(2,3,5)]
   sumas3coprimos 11  ==  []
   sumas3coprimos 12  ==  [(2,3,7),(3,4,5)]
   length (sumas3coprimos 4000)  ==  546146

sumas3coprimos 10 == [(2,3,5)]

sumas3coprimos 11 == []

sumas3coprimos 12 == [(2,3,7),(3,4,5)]

length (sumas3coprimos 4000) == 546146

Comprobar con QuickCheck que todo número entero mayor que 17 se puede escribir como suma de alguna terna coprima; es decir, para todo entero n, (sumas3coprimos2 (18 + abs n)) tiene algún elemento.

Soluciones

import Test.QuickCheck

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

sumas3coprimos :: Integer -> [(Integer,Integer,Integer)]
sumas3coprimos n =
  [(a,b,c) | a <- [2..n]
           , b <- [a+1..n]
           , c <- [b+1..n]
           , a + b + c == n
           , gcd a b == 1
           , gcd a c == 1
           , gcd b c == 1]

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

sumas3coprimos2 :: Integer -> [(Integer,Integer,Integer)]
sumas3coprimos2 n =
  [(a,b,c) | a <- [2..n `div` 3]
           , b <- [a+1..(n - a) `div` 2]
           , gcd a b == 1
           , let c = n - a - b  
           , gcd a c == 1
           , gcd b c == 1]

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

-- La propiedad de equivalencia es
prop_sumas3coprimos_equiv :: Integer -> Property
prop_sumas3coprimos_equiv n =
  n > 0 ==> sumas3coprimos n == sumas3coprimos2 n

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

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

--    λ> length (sumas3coprimos 400)
--    5345
--    (4.16 secs, 2,894,799,744 bytes)
--    λ> length (sumas3coprimos2 400)
--    5345
--    (0.06 secs, 16,565,136 bytes)

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

-- La propiedad
prop_sumas3coprimos :: Integer -> Bool
prop_sumas3coprimos n =
  not (null (sumas3coprimos2 (18 + abs n)))

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

import Test.QuickCheck

-- 1ª solución

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

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

sumas3coprimos n =

[(a,b,c) | a <- [2..n]

, b <- [a+1..n]

, c <- [b+1..n]

, a + b + c == n

, gcd a b == 1

, gcd a c == 1

, gcd b c == 1]

-- 2ª solución

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

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

sumas3coprimos2 n =

[(a,b,c) | a <- [2..n `div` 3]

, b <- [a+1..(n - a) `div` 2]

, gcd a b == 1

, let c = n - a - b

, gcd a c == 1

, gcd b c == 1]

-- Equivalencia de las definiciones

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

-- La propiedad de equivalencia es

prop_sumas3coprimos_equiv :: Integer -> Property

prop_sumas3coprimos_equiv n =

n > 0 ==> sumas3coprimos n == sumas3coprimos2 n

-- La comprobación es

-- λ> quickCheck prop_sumas3coprimos_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> length (sumas3coprimos 400)

-- 5345

-- (4.16 secs, 2,894,799,744 bytes)

-- λ> length (sumas3coprimos2 400)

-- 5345

-- (0.06 secs, 16,565,136 bytes)

-- Propiedad

-- =========

-- La propiedad

prop_sumas3coprimos :: Integer -> Bool

prop_sumas3coprimos n =

not (null (sumas3coprimos2 (18 + abs n)))

-- La comprobación es

-- λ> quickCheck prop_sumas3coprimos

-- +++ OK, passed 100 tests.

Referencias

Basado en Integers as sum of three pairwise coprime integers de ProofWiki.

Pensamiento

Todo amor es fantasía;
él inventa el año, el día,
la hora y su melodía;
inventa el amante y, más
la amada. No prueba nada,
contra el amor, que la amada
no haya existido jamás.

Antonio Machado

Pares definidos por su MCD y su MCM

PorJosé A. Alonso 13-diciembre-201920-diciembre-2019

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]

-- Comparación 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 pares 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 orden. 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]

-- Comparación 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 pares 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 orden. 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.

Pensamiento

Largo es el camino de la enseñanza por medio de teorías; breve y eficaz por medio de ejemplos. ~ Séneca

Huecos de Aquiles

PorJosé A. Alonso 2-abril-201926-marzo-2022

108 es un número de Aquiles proque es un número potente (ya que su factorización es 2^2 · 3^3, sus divisores primos son 2 and 3 y sus cuadrados (2^2 = 4 y 3^2 = 9) son divisores de 108. Además, 108 no es una potencia perfecta.
360 no es un número de Aquiles ya que 5 es un divisor primo de 360, pero 5^2 = 15 no lo es.
784 no es un número de Aquiles porque, aunque es potente, es una potencia perfecta ya que 784 = 28^2.

Los primeros números de Aquiles son

   72, 108, 200, 288, 392, 432, 500, 648, 675, 800, 864, 968, 972, ...

1	72, 108, 200, 288, 392, 432, 500, 648, 675, 800, 864, 968, 972, ...

Definir las funciones

   esAquiles              :: Integer -> Bool
   huecosDeAquiles        :: [Integer]
   graficaHuecosDeAquiles :: Int -> IO ()

esAquiles :: Integer -> Bool

huecosDeAquiles :: [Integer]

graficaHuecosDeAquiles :: Int -> IO ()

tales que

(esAquiles x) se verifica si x es un número de Aquiles. Por ejemplo,

     esAquiles 108         ==  True
     esAquiles 360         ==  False
     esAquiles 784         ==  False
     esAquiles 5425069447  ==  True
     esAquiles 5425069448  ==  True

esAquiles 108 == True

esAquiles 360 == False

esAquiles 784 == False

esAquiles 5425069447 == True

esAquiles 5425069448 == True

huecosDeAquiles es la sucesión de la diferencias entre los números de Aquiles consecutivos. Por ejemplo,

     λ> take 15 huecosDeAquiles
     [36,92,88,104,40,68,148,27,125,64,104,4,153,27,171]

1 2	λ> take 15 huecosDeAquiles [36,92,88,104,40,68,148,27,125,64,104,4,153,27,171]

(graficaHuecosDeAquiles n) dibuja la gráfica de los n primeros huecos de Aquiles. Por ejemplo, (graficaHuecosDeAquiles 160) dibuja

Soluciones

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

-- Definición de esAquiles
-- =======================

esAquiles :: Integer -> Bool
esAquiles x = esPotente x && noEsPotenciaPerfecta x

-- (esPotente x) se verifica si x es potente. Por ejemplo,
--    esPotente 108  ==  True
--    esPotente 360  ==  False
--    esPotente 784  ==  True
esPotente :: Integer -> Bool
esPotente x = all (>1) (exponentes x)

-- (exponentes x) es la lista de los exponentes en la factorización de
-- x. Por ejemplo,
--    exponentes 108  ==  [2,3]
--    exponentes 360  ==  [3,2,1]
--    exponentes 784  ==  [4,2]
exponentes :: Integer -> [Int]
exponentes x = map length (group (primeFactors x))

-- (noEsPotenciaPerfecta x) se verifica si x no es una potencia
-- perfecta. Por ejemplo,
--    noEsPotenciaPerfecta 108  ==  True
--    noEsPotenciaPerfecta 360  ==  True
--    noEsPotenciaPerfecta 784  ==  False
noEsPotenciaPerfecta :: Integer -> Bool
noEsPotenciaPerfecta x = foldl1 gcd (exponentes x) == 1 

-- Definición de huecosDeAquiles
-- =============================

huecosDeAquiles :: [Integer]
huecosDeAquiles = zipWith (-) (tail aquiles) aquiles

-- aquiles es la sucesión de los números de Aquiles. Por ejemplo, 
--    λ> take 15 aquiles
--    [72,108,200,288,392,432,500,648,675,800,864,968,972,1125,1152]
aquiles :: [Integer]
aquiles = filter esAquiles [2..]

-- Definición de graficaHuecosDeAquiles
-- ====================================

graficaHuecosDeAquiles :: Int -> IO ()
graficaHuecosDeAquiles n =
  plotList [ Key Nothing
           , PNG "Huecos_de_Aquiles.png"
           ]
           (take n huecosDeAquiles)

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

import Graphics.Gnuplot.Simple

-- Definición de esAquiles

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

esAquiles :: Integer -> Bool

esAquiles x = esPotente x && noEsPotenciaPerfecta x

-- (esPotente x) se verifica si x es potente. Por ejemplo,

-- esPotente 108 == True

-- esPotente 360 == False

-- esPotente 784 == True

esPotente :: Integer -> Bool

esPotente x = all (>1) (exponentes x)

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

-- x. Por ejemplo,

-- exponentes 108 == [2,3]

-- exponentes 360 == [3,2,1]

-- exponentes 784 == [4,2]

exponentes :: Integer -> [Int]

exponentes x = map length (group (primeFactors x))

-- (noEsPotenciaPerfecta x) se verifica si x no es una potencia

-- perfecta. Por ejemplo,

-- noEsPotenciaPerfecta 108 == True

-- noEsPotenciaPerfecta 360 == True

-- noEsPotenciaPerfecta 784 == False

noEsPotenciaPerfecta :: Integer -> Bool

noEsPotenciaPerfecta x = foldl1 gcd (exponentes x) == 1

-- Definición de huecosDeAquiles

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

huecosDeAquiles :: [Integer]

huecosDeAquiles = zipWith (-) (tail aquiles) aquiles

-- aquiles es la sucesión de los números de Aquiles. Por ejemplo,

-- λ> take 15 aquiles

-- [72,108,200,288,392,432,500,648,675,800,864,968,972,1125,1152]

aquiles :: [Integer]

aquiles = filter esAquiles [2..]

-- Definición de graficaHuecosDeAquiles

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

graficaHuecosDeAquiles :: Int -> IO ()

graficaHuecosDeAquiles n =

plotList [ Key Nothing

, PNG "Huecos_de_Aquiles.png"

]

(take n huecosDeAquiles)

Pensamiento

Tengo a mis amigos
en mi soledad;
cuando estoy con ellos
¡qué lejos están!

Antonio Machado

La sucesión ECG

PorJosé A. Alonso 21-marzo-201925-abril-2021

La sucesión ECG estás definida por a(1) = 1, a(2) = 2 y, para n >= 3, a(n) es el menor natural que aún no está en la sucesión tal que a(n) tiene algún divisor común con a(n-1).

Los primeros términos de la sucesión son 1, 2, 4, 6, 3, 9, 12, 8, 10, 5, 15, …

Al dibujar su gráfica, se parece a la de los electrocardiogramas (abreviadamente, ECG). Por ello, la sucesión se conoce como la sucesión ECG.

Definir las funciones

   sucECG :: [Integer]
   graficaSucECG :: Int -> IO ()

1 2	sucECG :: [Integer] graficaSucECG :: Int -> IO ()

tales que

sucECG es la lista de los términos de la sucesión ECG. Por ejemplo,

     λ> take 20 sucECG
     [1,2,4,6,3,9,12,8,10,5,15,18,14,7,21,24,16,20,22,11]
     λ> sucECG !! 6000
     6237

λ> take 20 sucECG

[1,2,4,6,3,9,12,8,10,5,15,18,14,7,21,24,16,20,22,11]

λ> sucECG !! 6000

6237

(graficaSucECG n) dibuja la gráfica de los n primeros términos de la sucesión ECG. Por ejemplo, (graficaSucECG 160) dibuja

Soluciones

import Data.List (delete)
import Graphics.Gnuplot.Simple

sucECG :: [Integer]
sucECG = 1 : ecg 2 [2..]
  where ecg x zs = f zs
          where f (y:ys) | gcd x y > 1 = y : ecg y (delete y zs)
                         | otherwise   = f ys

graficaSucECG :: Int -> IO ()
graficaSucECG n =
  plotList [ Key Nothing
           , PNG "La_sucesion_ECG.png" 
           ]
           (take n sucECG)

import Data.List (delete)

import Graphics.Gnuplot.Simple

sucECG :: [Integer]

sucECG = 1 : ecg 2 [2..]

where ecg x zs = f zs

where f (y:ys) | gcd x y > 1 = y : ecg y (delete y zs)

| otherwise = f ys

graficaSucECG :: Int -> IO ()

graficaSucECG n =

plotList [ Key Nothing

, PNG "La_sucesion_ECG.png"

]

(take n sucECG)

Pensamiento

Algunos desesperados
sólo se curan con soga;
otros, con siete palabras:
la fe se ha puesto de moda.

Antonio Machado

Mayor exponente

PorJosé A. Alonso 20-febrero-201927-febrero-2019

Definir las funciones

   mayorExponente        :: Integer -> Integer
   graficaMayorExponente :: Integer -> IO ()

1 2	mayorExponente :: Integer -> Integer graficaMayorExponente :: Integer -> IO ()

tales que

(mayorExponente n) es el mayor número b para el que existe un a tal que n = a^b. Se supone que n > 1. Por ejemplo,

     mayorExponente 9   ==  2
     mayorExponente 8   ==  3
     mayorExponente 7   ==  1
     mayorExponente 18  ==  1
     mayorExponente 36  ==  2
     mayorExponente (10^(10^5))  ==  100000

mayorExponente 9 == 2

mayorExponente 8 == 3

mayorExponente 7 == 1

mayorExponente 18 == 1

mayorExponente 36 == 2

mayorExponente (10^(10^5)) == 100000

(graficaMayorExponente n) dibuja la gráfica de los mayores exponentes de los números entre 2 y n. Por ejemplo, (graficaMayorExponente 50) dibuja

Soluciones

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


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

mayorExponente :: Integer -> Integer
mayorExponente x =
  last [b | b <- [1..x]
          , a <- [1..x]
          , a^b == x]

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

mayorExponente2 :: Integer -> Integer
mayorExponente2 x =
  head [b | b <- [x,x-1..1]
          , a <- [1..x]
          , a^b == x]

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

mayorExponente3 :: Integer -> Integer
mayorExponente3 x = aux x
  where aux 1 = 1
        aux b | any (\a -> a^b == x) [1..x] = b
              | otherwise                   = aux (b-1)

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

mayorExponente4 :: Integer -> Integer
mayorExponente4 x =
  mcd (exponentes x)

-- (exponentes x) es la lista de los exponentes en la factorizacioń de
-- x. por ejemplo.
--    exponentes 720  ==  [4,2,1]
exponentes :: Integer -> [Integer]
exponentes x =
  map genericLength (group (primeFactors x))

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

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

-- La propiedad es
prop_mayorExponente :: Integer -> Property
prop_mayorExponente n =
  n >= 0 ==>
  mayorExponente  n == mayorExponente2 n &&
  mayorExponente2 n == mayorExponente3 n

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

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

--    λ> mayorExponente (10^3)
--    3
--    (3.96 secs, 4,671,928,464 bytes)
--    λ> mayorExponente2 (10^3)
--    3
--    (3.99 secs, 4,670,107,024 bytes)
--    λ> mayorExponente3 (10^3)
--    3
--    (3.90 secs, 4,686,383,952 bytes)
--    λ> mayorExponente4 (10^3)
--    3
--    (0.02 secs, 131,272 bytes)

-- Definición de graficaMayorExponente
-- ======================================

graficaMayorExponente :: Integer -> IO ()
graficaMayorExponente n = 
  plotList [ Key Nothing
           , PNG ("MayorExponente.png")
           ]
           (map mayorExponente3 [2..n])

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

import Graphics.Gnuplot.Simple

-- 1ª solución

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

mayorExponente :: Integer -> Integer

mayorExponente x =

last [b | b <- [1..x]

, a <- [1..x]

, a^b == x]

-- 2ª solución

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

mayorExponente2 :: Integer -> Integer

mayorExponente2 x =

head [b | b <- [x,x-1..1]

, a <- [1..x]

, a^b == x]

-- 3ª solución

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

mayorExponente3 :: Integer -> Integer

mayorExponente3 x = aux x

where aux 1 = 1

aux b | any (\a -> a^b == x) [1..x] = b

| otherwise = aux (b-1)

-- 4ª solución

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

mayorExponente4 :: Integer -> Integer

mayorExponente4 x =

mcd (exponentes x)

-- (exponentes x) es la lista de los exponentes en la factorizacioń de

-- x. por ejemplo.

-- exponentes 720 == [4,2,1]

exponentes :: Integer -> [Integer]

exponentes x =

map genericLength (group (primeFactors x))

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

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

mcd :: [Integer] -> Integer

mcd = foldr1 gcd

-- Comprobación de equivalencia

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

-- La propiedad es

prop_mayorExponente :: Integer -> Property

prop_mayorExponente n =

n >= 0 ==>

mayorExponente n == mayorExponente2 n &&

mayorExponente2 n == mayorExponente3 n

-- La comprobación es

-- λ> quickCheck prop_mayorExponente

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> mayorExponente (10^3)

-- 3

-- (3.96 secs, 4,671,928,464 bytes)

-- λ> mayorExponente2 (10^3)

-- 3

-- (3.99 secs, 4,670,107,024 bytes)

-- λ> mayorExponente3 (10^3)

-- 3

-- (3.90 secs, 4,686,383,952 bytes)

-- λ> mayorExponente4 (10^3)

-- 3

-- (0.02 secs, 131,272 bytes)

-- Definición de graficaMayorExponente

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

graficaMayorExponente :: Integer -> IO ()

graficaMayorExponente n =

plotList [ Key Nothing

, PNG ("MayorExponente.png")

]

(map mayorExponente3 [2..n])

Pensamiento

Mirando mi calavera
un nuevo Hamlet dirá:
He aquí un lindo fósil de una
careta de carnaval.

Antonio Machado

Entre dos conjuntos

PorJosé A. Alonso 12-diciembre-201817-enero-2019

Se dice que un x número se encuentra entre dos conjuntos xs e ys si x es divisible por todos los elementos de xs y todos los elementos de zs son divisibles por x. Por ejemplo, 12 se encuentra entre los conjuntos {2, 6} y {24, 36}.

Definir la función

   entreDosConjuntos :: [Int] -> [Int] -> [Int]

1	entreDosConjuntos :: [Int] -> [Int] -> [Int]

tal que (entreDosConjuntos xs ys) es la lista de elementos entre xs e ys (se supone que xs e ys son listas no vacías de números enteros positivos). Por ejemplo,

   entreDosConjuntos [2,6] [24,36]     ==  [6,12]
   entreDosConjuntos [2,4] [32,16,96]  ==  [4,8,16]

1 2	entreDosConjuntos [2,6] [24,36] == [6,12] entreDosConjuntos [2,4] [32,16,96] == [4,8,16]

Otros ejemplos

   λ> (xs,a) = ([1..15],product xs) 
   λ> length (entreDosConjuntos5 xs [a,2*a..10*a])
   270
   λ> (xs,a) = ([1..16],product xs) 
   λ> length (entreDosConjuntos5 xs [a,2*a..10*a])
   360

λ> (xs,a) = ([1..15],product xs)

λ> length (entreDosConjuntos5 xs [a,2*a..10*a])

270

λ> (xs,a) = ([1..16],product xs)

λ> length (entreDosConjuntos5 xs [a,2*a..10*a])

360

Soluciones

import Test.QuickCheck 

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

entreDosConjuntos :: [Int] -> [Int] -> [Int]
entreDosConjuntos xs ys =
  [z | z <- [a..b]
     , and [z `mod` x == 0 | x <- xs]
     , and [y `mod` z == 0 | y <- ys]]
  where a = maximum xs
        b = minimum ys

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

entreDosConjuntos2 :: [Int] -> [Int] -> [Int]
entreDosConjuntos2 xs ys =
  [z | z <- [a..b]
     , all (`divideA` z) xs
     , all (z `divideA`) ys]
  where a = mcmL xs
        b = mcdL ys

--    mcmL [2,3,18]  ==  18
--    mcmL [2,3,15]  ==  30
mcdL :: [Int] -> Int
mcdL [x]    = x
mcdL (x:xs) = gcd x (mcdL xs)

--    mcmL [12,30,18]  ==  6
--    mcmL [12,30,14]  ==  2
mcmL :: [Int] -> Int
mcmL [x]    = x
mcmL (x:xs) = lcm x (mcmL xs)

divideA :: Int -> Int -> Bool
divideA x y = y `mod` x == 0

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

entreDosConjuntos3 :: [Int] -> [Int] -> [Int]
entreDosConjuntos3 xs ys =
  [z | z <- [a..b]
     , all (`divideA` z) xs
     , all (z `divideA`) ys]
  where a = mcmL2 xs
        b = mcdL2 ys

-- Definición equivalente
mcdL2 :: [Int] -> Int
mcdL2 = foldl1 gcd

-- Definición equivalente
mcmL2 :: [Int] -> Int
mcmL2 = foldl1 lcm

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

entreDosConjuntos4 :: [Int] -> [Int] -> [Int]
entreDosConjuntos4 xs ys =
  [z | z <- [a,a+a..b]
     , z `divideA` b] 
  where a = mcmL2 xs
        b = mcdL2 ys

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

entreDosConjuntos5 :: [Int] -> [Int] -> [Int]
entreDosConjuntos5 xs ys =
  filter (`divideA` b) [a,a+a..b]
  where a = mcmL2 xs
        b = mcdL2 ys

-- Equivalencia
-- ============

-- Para comprobar la equivalencia se define el tipo de listas no vacías
-- de números enteros positivos:
newtype ListaNoVaciaDePositivos = L [Int]
  deriving Show

-- genListaNoVaciaDePositivos es un generador de listas no vacióas de
-- enteros positivos. Por ejemplo,
--    λ> sample genListaNoVaciaDePositivos
--    L [1]
--    L [1,2,2]
--    L [4,3,4]
--    L [1,6,5,2,4]
--    L [2,8]
--    L [11]
--    L [13,2,3]
--    L [7,3,9,15,11,12,13,3,9,6,13,3]
--    L [16,2,11,10,6,5,16,4,1,15,9,11,8,15,2,15,7]
--    L [5,4,9,13,5,6,7]
--    L [7,4,6,12,2,11,6,14,14,13,14,11,6,2,18,8,16,2,13,9]
genListaNoVaciaDePositivos :: Gen ListaNoVaciaDePositivos
genListaNoVaciaDePositivos = do
  x  <- arbitrary
  xs <- arbitrary
  return (L (map ((+1) . abs) (x:xs)))

-- Generación arbitraria de listas no vacías de enteros positivos.
instance Arbitrary ListaNoVaciaDePositivos where
  arbitrary = genListaNoVaciaDePositivos

-- La propiedad es
prop_entreDosConjuntos_equiv ::
     ListaNoVaciaDePositivos
  -> ListaNoVaciaDePositivos
  -> Bool
prop_entreDosConjuntos_equiv (L xs) (L ys) =
  entreDosConjuntos xs ys == entreDosConjuntos2 xs ys &&
  entreDosConjuntos xs ys == entreDosConjuntos3 xs ys &&
  entreDosConjuntos xs ys == entreDosConjuntos4 xs ys &&
  entreDosConjuntos xs ys == entreDosConjuntos5 xs ys 

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

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

--    λ> (xs,a) = ([1..10],product xs) 
--    λ> length (entreDosConjuntos xs [a,2*a..10*a])
--    36
--    (5.08 secs, 4,035,689,200 bytes)
--    λ> length (entreDosConjuntos2 xs [a,2*a..10*a])
--    36
--    (3.75 secs, 2,471,534,072 bytes)
--    λ> length (entreDosConjuntos3 xs [a,2*a..10*a])
--    36
--    (3.73 secs, 2,471,528,664 bytes)
--    λ> length (entreDosConjuntos4 xs [a,2*a..10*a])
--    36
--    (0.01 secs, 442,152 bytes)
--    λ> length (entreDosConjuntos5 xs [a,2*a..10*a])
--    36
--    (0.00 secs, 374,824 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

import Test.QuickCheck

-- 1ª solución

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

entreDosConjuntos :: [Int] -> [Int] -> [Int]

entreDosConjuntos xs ys =

[z | z <- [a..b]

, and [z `mod` x == 0 | x <- xs]

, and [y `mod` z == 0 | y <- ys]]

where a = maximum xs

b = minimum ys

-- 2ª solución

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

entreDosConjuntos2 :: [Int] -> [Int] -> [Int]

entreDosConjuntos2 xs ys =

[z | z <- [a..b]

, all (`divideA` z) xs

, all (z `divideA`) ys]

where a = mcmL xs

b = mcdL ys

-- mcmL [2,3,18] == 18

-- mcmL [2,3,15] == 30

mcdL :: [Int] -> Int

mcdL [x] = x

mcdL (x:xs) = gcd x (mcdL xs)

-- mcmL [12,30,18] == 6

-- mcmL [12,30,14] == 2

mcmL :: [Int] -> Int

mcmL [x] = x

mcmL (x:xs) = lcm x (mcmL xs)

divideA :: Int -> Int -> Bool

divideA x y = y `mod` x == 0

-- 3ª solución

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

entreDosConjuntos3 :: [Int] -> [Int] -> [Int]

entreDosConjuntos3 xs ys =

[z | z <- [a..b]

, all (`divideA` z) xs

, all (z `divideA`) ys]

where a = mcmL2 xs

b = mcdL2 ys

-- Definición equivalente

mcdL2 :: [Int] -> Int

mcdL2 = foldl1 gcd

-- Definición equivalente

mcmL2 :: [Int] -> Int

mcmL2 = foldl1 lcm

-- 4ª solución

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

entreDosConjuntos4 :: [Int] -> [Int] -> [Int]

entreDosConjuntos4 xs ys =

[z | z <- [a,a+a..b]

, z `divideA` b]

where a = mcmL2 xs

b = mcdL2 ys

-- 5ª solución

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

entreDosConjuntos5 :: [Int] -> [Int] -> [Int]

entreDosConjuntos5 xs ys =

filter (`divideA` b) [a,a+a..b]

where a = mcmL2 xs

b = mcdL2 ys

-- Equivalencia

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

-- Para comprobar la equivalencia se define el tipo de listas no vacías

-- de números enteros positivos:

newtype ListaNoVaciaDePositivos = L [Int]

deriving Show

-- genListaNoVaciaDePositivos es un generador de listas no vacióas de

-- enteros positivos. Por ejemplo,

-- λ> sample genListaNoVaciaDePositivos

-- L [1]

-- L [1,2,2]

-- L [4,3,4]

-- L [1,6,5,2,4]

-- L [2,8]

-- L [11]

-- L [13,2,3]

-- L [7,3,9,15,11,12,13,3,9,6,13,3]

-- L [16,2,11,10,6,5,16,4,1,15,9,11,8,15,2,15,7]

-- L [5,4,9,13,5,6,7]

-- L [7,4,6,12,2,11,6,14,14,13,14,11,6,2,18,8,16,2,13,9]

genListaNoVaciaDePositivos :: Gen ListaNoVaciaDePositivos

genListaNoVaciaDePositivos = do

x <- arbitrary

xs <- arbitrary

return (L (map ((+1) . abs) (x:xs)))

-- Generación arbitraria de listas no vacías de enteros positivos.

instance Arbitrary ListaNoVaciaDePositivos where

arbitrary = genListaNoVaciaDePositivos

-- La propiedad es

prop_entreDosConjuntos_equiv ::

ListaNoVaciaDePositivos

-> ListaNoVaciaDePositivos

-> Bool

prop_entreDosConjuntos_equiv (L xs) (L ys) =

entreDosConjuntos xs ys == entreDosConjuntos2 xs ys &&

entreDosConjuntos xs ys == entreDosConjuntos3 xs ys &&

entreDosConjuntos xs ys == entreDosConjuntos4 xs ys &&

entreDosConjuntos xs ys == entreDosConjuntos5 xs ys

-- La comprobación es

-- λ> quickCheck prop_entreDosConjuntos_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> (xs,a) = ([1..10],product xs)

-- λ> length (entreDosConjuntos xs [a,2*a..10*a])

-- 36

-- (5.08 secs, 4,035,689,200 bytes)

-- λ> length (entreDosConjuntos2 xs [a,2*a..10*a])

-- 36

-- (3.75 secs, 2,471,534,072 bytes)

-- λ> length (entreDosConjuntos3 xs [a,2*a..10*a])

-- 36

-- (3.73 secs, 2,471,528,664 bytes)

-- λ> length (entreDosConjuntos4 xs [a,2*a..10*a])

-- 36

-- (0.01 secs, 442,152 bytes)

-- λ> length (entreDosConjuntos5 xs [a,2*a..10*a])

-- 36

-- (0.00 secs, 374,824 bytes)

Referencia

Este ejercicio está basado en el problema Between two sets de HackerRank.

Pensamiento

Las razones no se transmiten, se engendran, por cooperación, en el diálogo.

Antonio Machado

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

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.

Máximo común divisor de x e y veces n

PorJosé A. Alonso 9-febrero-201716-febrero-2017

Definir las funciones

   repite :: Int -> Integer -> Integer
   mcdR   :: Integer -> Int -> Int -> Integer

1 2	repite :: Int -> Integer -> Integer mcdR :: Integer -> Int -> Int -> Integer

tales que

(repite x n) es el número obtenido repitiendo x veces el número n. Por ejemplo.

     repite 3 123  ==  123123123

1	repite 3 123 == 123123123

(mcdR n x y) es el máximo común divisor de los números obtenidos repitiendo x veces e y veces el número n. Por ejemplo.

     mcdR 123 2 3                     ==  123
     mcdR 4 4 6                       ==  44
     mcdR 2017 (10^1000) (2+10^1000)  ==  20172017

mcdR 123 2 3 == 123

mcdR 4 4 6 == 44

mcdR 2017 (10^1000) (2+10^1000) == 20172017

Soluciones

import Test.QuickCheck

repite :: Int -> Integer -> Integer
repite x n =
  read (concat (replicate x (show n)))

-- 1ª definición
mcdR :: Integer -> Int -> Int -> Integer
mcdR n x y =
  gcd (repite x n) (repite y n)

-- 2ª definición
mcdR2 :: Integer -> Int -> Int -> Integer
mcdR2 n x y =
  repite (gcd x y) n

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

-- La propiedad es
prop_equivalencia :: (Positive Integer) -> (Positive Int) -> (Positive Int)
                     -> Bool
prop_equivalencia (Positive n) (Positive x) (Positive y) =
  mcdR n x y == mcdR2 n x y

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

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

--    λ> mcdR 2017 (10^6) (2+10^6)
--    20172017
--    (18.92 secs, 6,012,285,224 bytes)
--    λ> mcdR2 2017 (10^6) (2+10^6)
--    20172017
--    (0.01 secs, 0 bytes)

import Test.QuickCheck

repite :: Int -> Integer -> Integer

repite x n =

read (concat (replicate x (show n)))

-- 1ª definición

mcdR :: Integer -> Int -> Int -> Integer

mcdR n x y =

gcd (repite x n) (repite y n)

-- 2ª definición

mcdR2 :: Integer -> Int -> Int -> Integer

mcdR2 n x y =

repite (gcd x y) n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_equivalencia :: (Positive Integer) -> (Positive Int) -> (Positive Int)

-> Bool

prop_equivalencia (Positive n) (Positive x) (Positive y) =

mcdR n x y == mcdR2 n x y

-- La comprobación es

-- λ> quickCheck prop_equivalencia

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> mcdR 2017 (10^6) (2+10^6)

-- 20172017

-- (18.92 secs, 6,012,285,224 bytes)

-- λ> mcdR2 2017 (10^6) (2+10^6)

-- 20172017

-- (0.01 secs, 0 bytes)

Inverso multiplicativo modular

PorJosé A. Alonso 6-abril-201626-marzo-2022

El inverso multiplicativo modular de un entero n módulo p es el número m, entre 1 y p-1, tal que

   mn = 1 (mod p)

1	mn = 1 (mod p)

Por ejemplo, el inverso multiplicativo de 2 módulo 5 es 3, ya que 1 <= 3 <= 4 y 2×3 = 1 (mod 5).

El inverso multipicativo de n módulo p existe si y sólo si n y p son coprimos; es decir, si mcd(n,p) = 1.

Definir la función

   invMod :: Integer -> Integer -> Maybe Integer

1	invMod :: Integer -> Integer -> Maybe Integer

tal que (invMod n p) es justo el inverso multiplicativo de n módulo p, si existe y Nothing en caso contrario. Por ejemplo,

   λ> invMod 2 5
   Just 3
   λ> invMod 2 6
   Nothing
   λ> [(x,invMod x 5) | x <- [0..4]]
   [(0,Nothing),(1,Just 1),(2,Just 3),(3,Just 2),(4,Just 4)]
   λ> [(x,invMod x 6) | x <- [0..5]]
   [(0,Nothing),(1,Just 1),(2,Nothing),(3,Nothing),(4,Nothing),(5,Just 5)]
   λ> let n = 10^7 in invMod (10^n) (1+10^n) == Just (10^n)
   True

λ> invMod 2 5

Just 3

λ> invMod 2 6

Nothing

λ> [(x,invMod x 5) | x <- [0..4]]

[(0,Nothing),(1,Just 1),(2,Just 3),(3,Just 2),(4,Just 4)]

λ> [(x,invMod x 6) | x <- [0..5]]

[(0,Nothing),(1,Just 1),(2,Nothing),(3,Nothing),(4,Nothing),(5,Just 5)]

λ> let n = 10^7 in invMod (10^n) (1+10^n) == Just (10^n)

True

Soluciones

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

invMod1 :: Integer -> Integer -> Maybe Integer
invMod1 n p | gcd n p == 1 = Just (head [m | m <- [1..p-1], m*n `mod` p == 1])
            | otherwise    = Nothing

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

invMod2 :: Integer -> Integer -> Maybe Integer
invMod2 n p | m /= 1    = Nothing
            | x < 0     = Just (x + p)
            | otherwise = Just x
    where (x,_,m) = mcdExt n p

-- [Algoritmo extendido de Euclides]
-- (mcd a b) es la terna (x,y,g) tal que g es el máximo común divisor
-- de a y b y se cumple que ax + by = g. Por ejemplo,
--    mcdExt  2  5  ==  (-2,1,1)
--    mcdExt  2  6  ==  ( 1,0,2)
--    mcdExt 12 15  ==  (-1,1,3)
mcdExt :: Integer -> Integer -> (Integer,Integer,Integer)
mcdExt a 0 = (1, 0, a)
mcdExt a b = (t, s - q * t, g)
  where (q, r)    = a `quotRem` b
        (s, t, g) = mcdExt b r

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

invMod3 :: Integer -> Integer -> Maybe Integer
invMod3 n p | gcd n p == 1 = Just (aux n p)
            | otherwise    = Nothing
  where aux 1 p = 1
        aux n p = (x * p + 1) `div` n
          where x = n - aux (p `mod` n) n                    

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

--    λ> invMod1 (10^7) (1+10^7)
--    Just 10000000
--    (5.05 secs, 2,872,350,896 bytes)
--    λ> invMod2 (10^7) (1+10^7)
--    Just 10000000
--    (0.00 secs, 0 bytes)
--    λ> invMod3 (10^7) (1+10^7)
--    Just 10000000
--    (0.00 secs, 0 bytes)
--    
--    λ> let n = 10^7 in invMod2 (10^n) (1+10^n) == Just (10^n)
--    True
--    (0.88 secs, 55,728,120 bytes)
--    λ> let n = 10^7 in invMod3 (10^n) (1+10^n) == Just (10^n)
--    True
--    (1.99 secs, 80,644,352 bytes)

-- 1ª definición

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

invMod1 :: Integer -> Integer -> Maybe Integer

invMod1 n p | gcd n p == 1 = Just (head [m | m <- [1..p-1], m*n `mod` p == 1])

| otherwise = Nothing

-- 2ª definición

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

invMod2 :: Integer -> Integer -> Maybe Integer

invMod2 n p | m /= 1 = Nothing

| x < 0 = Just (x + p)

| otherwise = Just x

where (x,_,m) = mcdExt n p

-- [Algoritmo extendido de Euclides]

-- (mcd a b) es la terna (x,y,g) tal que g es el máximo común divisor

-- de a y b y se cumple que ax + by = g. Por ejemplo,

-- mcdExt 2 5 == (-2,1,1)

-- mcdExt 2 6 == ( 1,0,2)

-- mcdExt 12 15 == (-1,1,3)

mcdExt :: Integer -> Integer -> (Integer,Integer,Integer)

mcdExt a 0 = (1, 0, a)

mcdExt a b = (t, s - q * t, g)

where (q, r) = a `quotRem` b

(s, t, g) = mcdExt b r

-- 3ª definición

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

invMod3 :: Integer -> Integer -> Maybe Integer

invMod3 n p | gcd n p == 1 = Just (aux n p)

| otherwise = Nothing

where aux 1 p = 1

aux n p = (x * p + 1) `div` n

where x = n - aux (p `mod` n) n

-- Comparación de eficiencia

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

-- λ> invMod1 (10^7) (1+10^7)

-- Just 10000000

-- (5.05 secs, 2,872,350,896 bytes)

-- λ> invMod2 (10^7) (1+10^7)

-- Just 10000000

-- (0.00 secs, 0 bytes)

-- λ> invMod3 (10^7) (1+10^7)

-- Just 10000000

-- (0.00 secs, 0 bytes)

-- λ> let n = 10^7 in invMod2 (10^n) (1+10^n) == Just (10^n)

-- True

-- (0.88 secs, 55,728,120 bytes)

-- λ> let n = 10^7 in invMod3 (10^n) (1+10^n) == Just (10^n)

-- True

-- (1.99 secs, 80,644,352 bytes)

Solución en Maxima

invMod (n,p) := block ([r],
  r : inv_mod (n,p),
  if integerp (r)
  then Just (r)
  else Nothing)$

invMod (n,p) := block ([r],

r : inv_mod (n,p),

if integerp (r)

then Just (r)

else Nothing)$

La evaluación de los ejemplos es

(%i2) invMod (2,5);
(%o2) Just(3)
(%i3) invMod (2,6);
(%o3) Nothing
(%i4) makelist([x,invMod(x,5)],x,0,4);
(%o4) [[0,Nothing],[1,Just(1)],[2,Just(3)],[3,Just(2)],[4,Just(4)]]
(%i5) makelist([x,invMod(x,6)],x,0,5);
(%o5) [[0,Nothing],[1,Just(1)],[2,Nothing],[3,Nothing],[4,Nothing],[5,Just(5)]]
(%i6) (n : 10^7, is (invMod (10^n,1+10^n) = Just (10^n)));
(%o6) true

(%i2) invMod (2,5);

(%o2) Just(3)

(%i3) invMod (2,6);

(%o3) Nothing

(%i4) makelist([x,invMod(x,5)],x,0,4);

(%o4) [[0,Nothing],[1,Just(1)],[2,Just(3)],[3,Just(2)],[4,Just(4)]]

(%i5) makelist([x,invMod(x,6)],x,0,5);

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

(%i6) (n : 10^7, is (invMod (10^n,1+10^n) = Just (10^n)));

(%o6) true

Referencia

Calculating multiplicative inverses in modular arithmetic

Representación decimal de números racionales

PorJosé A. Alonso 1-abril-20168-abril-2016

Los números decimales se representan por ternas, donde el primer elemento es la parte entera, el segundo es el anteperíodo y el tercero es el período. Por ejemplo,

    6/2  = 3                  se representa por (3,[],[])
    1/2  = 0.5                se representa por (0,[5],[])
    1/3  = 0.333333...        se representa por (0,[],[3])  
   23/14 = 1.6428571428571... se representa por (1,[6],[4,2,8,5,7,1])

6/2 = 3 se representa por (3,[],[])

1/2 = 0.5 se representa por (0,[5],[])

1/3 = 0.333333... se representa por (0,[],[3])

23/14 = 1.6428571428571... se representa por (1,[6],[4,2,8,5,7,1])

Su tipo es

   type Decimal = (Integer,[Integer],[Integer])

1	type Decimal = (Integer,[Integer],[Integer])

Los números racionales se representan por un par de enteros, donde el primer elemento es el numerador y el segundo el denominador. Por ejemplo, el número 2/3 se representa por (2,3). Su tipo es

   type Racional = (Integer,Integer)

1	type Racional = (Integer,Integer)

Definir las funciones

   decimal  :: Racional -> Decimal
   racional :: Decimal -> Racional

1 2	decimal :: Racional -> Decimal racional :: Decimal -> Racional

tales que

(decimal r) es la representación decimal del número racional r. Por ejemplo,

        decimal (1,4)    ==  (0,[2,5],[])
        decimal (1,3)    ==  (0,[],[3])
        decimal (23,14)  ==  (1,[6],[4,2,8,5,7,1])

decimal (1,4) == (0,[2,5],[])

decimal (1,3) == (0,[],[3])

decimal (23,14) == (1,[6],[4,2,8,5,7,1])

(racional d) es el número racional cuya representación decimal es d. Por ejemplo,

        racional (0,[2,5],[])           ==  (1,4)
        racional (0,[],[3])             ==  (1,3)
        racional (1,[6],[4,2,8,5,7,1])  ==  (23,14)

racional (0,[2,5],[]) == (1,4)

racional (0,[],[3]) == (1,3)

racional (1,[6],[4,2,8,5,7,1]) == (23,14)

Con la función decimal se puede calcular los períodos de los números racionales. Por ejemplo,

   ghci> let (_,_,p) = decimal (23,14) in concatMap show p
   "428571"
   ghci> let (_,_,p) = decimal (1,47) in concatMap show p
   "0212765957446808510638297872340425531914893617"
   ghci> let (_,_,p) = decimal (1,541) in length (concatMap show p)
   540

ghci> let (_,_,p) = decimal (23,14) in concatMap show p

"428571"

ghci> let (_,_,p) = decimal (1,47) in concatMap show p

"0212765957446808510638297872340425531914893617"

ghci> let (_,_,p) = decimal (1,541) in length (concatMap show p)

540

Comprobar con QuickCheck si las funciones decimal y racional son inversas.

Soluciones

-- ---------------------------------------------------------------------
-- § Librerías auxiliares                                             --
-- ---------------------------------------------------------------------

import Data.List
import Test.QuickCheck

-- ---------------------------------------------------------------------
-- § Tipos                                                            --
-- ---------------------------------------------------------------------

-- Un número racional se representa por un par de enteros, donde el
-- primer elemento es el numerador y el segundo el denominador.
type Racional = (Integer,Integer)

-- Un número decimal es una terna donde el primer elemento es la parte
-- entera, el segundo es el anteperíodo y el tercero es el período.
type Decimal = (Integer,[Integer],[Integer])

-- ---------------------------------------------------------------------
-- § De racional a decimal                                            --
-- ---------------------------------------------------------------------

-- (mayorExponente a x) es el mayor n tal que a^n divide a x. Por ejemplo,
--    mayorExponente 2 40  ==  3
--    mayorExponente 5 40  ==  1
--    mayorExponente 5 41  ==  0
mayorExponente :: Integer -> Integer -> Integer
mayorExponente a x = head [n | n <- [1..], mod x (a^n) /= 0] - 1

-- (longitudAnteperiodo n) es la longitud del anteperíodo de 1/n (y de
-- cualquier fracción irreducible de denominador n). Por ejemplo, 
--    longitudAnteperiodo 2   ==  1
--    longitudAnteperiodo 3   ==  0
--    longitudAnteperiodo 14  ==  1
longitudAnteperiodo :: Integer -> Integer
longitudAnteperiodo n = max (mayorExponente 2 n) (mayorExponente 5 n)

-- (longitudPeriodo n) es la longitud del período de 1/n (y de
-- cualquier fracción irreducible de denominador n). Por ejemplo, 
--    longitudPeriodo 2  ==  0
--    longitudPeriodo 3  ==  1
--    longitudPeriodo 7  ==  6
longitudPeriodo :: Integer -> Integer
longitudPeriodo n
    | m == 1    = 0 
    | otherwise = head [k | k <- [1..], (10^k) `mod` m == 1]
    where m = n `div` (2^(mayorExponente 2 n) * 5^(mayorExponente 5 n))

-- (expansionDec (x,y)) es la expansión decimal de x/y. Por ejemplo, 
--    take 10 (expansionDec (1,4))    ==  [0,2,5]
--    take 10 (expansionDec (1,7))    ==  [0,1,4,2,8,5,7,1,4,2]
--    take 12 (expansionDec (90,7))   ==  [12,8,5,7,1,4,2,8,5,7,1,4]
--    take 12 (expansionDec (23,14))  ==  [1,6,4,2,8,5,7,1,4,2,8,5]
expansionDec :: Racional -> [Integer]
expansionDec (x,y) 
    | r == 0    = [q] 
    | otherwise = q : expansionDec (r*10,y)
    where (q,r) = quotRem x y

-- (parteEntera (a,b)) es la parte entera de a/b. Por ejemplo,
--    parteEntera (125,3)  ==  41
parteEntera :: Racional -> Integer
parteEntera (a,b) = a `div` b

-- (antePeriodo (a,b)) es el anteperíodo de a/b; es decir, la lista de
-- cifras que se encuentran entre la parte entera y el primer período de
-- a/b. Por ejemplo,
--    antePeriodo (23,14)  ==  [6]
--    antePeriodo (1,5)    ==  [2]
antePeriodo :: Racional -> [Integer]
antePeriodo (a,b) = genericTake s (tail xs)
    where xs = expansionDec (a,b)
          s  = longitudAnteperiodo b

-- (periodo (a,b)) es el período de a/b; es decir, la lista de cifras que
-- se encuentran entre la parte entera y el primer período de a/b. Por
-- ejemplo, 
--    periodo (1,3)                   ==  [3]
--    periodo (1,5)                   ==  []
--    periodo (1,7)                   ==  [1,4,2,8,5,7]
--    periodo (23,14)                 ==  [4,2,8,5,7,1]
--    concatMap show $ periodo (1,29) == "0344827586206896551724137931"
periodo :: Racional -> [Integer]
periodo (a,b) = genericTake t (genericDrop (1+s) xs)
    where xs = expansionDec (a,b)
          s  = longitudAnteperiodo b
          t  = longitudPeriodo b

-- (reducido (a,b)) es la forma reducida del número racional a/b. Por
-- ejemplo, 
--    reducido (40,80)  ==  (1,2)
reducido :: Racional -> Racional
reducido (a,b) = (a `div` m, b `div` m) 
    where m = gcd a b

-- (decimal (x,y)) es la forma decimal de x/y; es decir, la terna
-- formada por la parte entera, la parte decimal pura y la parte decimal
-- periódica. Por ejemplo, 
--    decimal (1,4)    ==  (0,[2,5],[])
--    decimal (1,3)    ==  (0,[],[3])
--    decimal (23,14)  ==  (1,[6],[4,2,8,5,7,1])
decimal :: Racional -> Decimal
decimal (a,b) = (parteEntera r, antePeriodo r, periodo r)
    where r = reducido (a,b)

-- ---------------------------------------------------------------------
-- § De decimal a racional                                            --
-- ---------------------------------------------------------------------

-- (digitosAnumero xs) es el número correspondiente a la lista de
-- dígitos xs. Por ejemplo,
--    digitosAnumero [3,2,5]  ==  325
digitosAnumero :: [Integer] -> Integer
digitosAnumero = foldl (\a y -> 10*a+y) 0

-- 2ª definición de digitosAnumero (por comprensión)
digitosAnumero2 :: [Integer] -> Integer
digitosAnumero2 xs = sum [x*(10^i) | (x,i) <- zip xs [n-1,n-2..0]]
    where n = length xs

-- (racional (x,ys,zs)) es el número racional cuya representación
-- decimal es (x,ys,zs). Por ejemplo,
--    racional (0,[2,5],[])           ==  (1,4)
--    racional (0,[],[3])             ==  (1,3)
--    racional (1,[6],[4,2,8,5,7,1])  ==  (23,14)
racional :: Decimal -> Racional
racional (x,ys,[]) = reducido (a,b) where
    a = digitosAnumero (x:ys)
    b = 10^(length ys)
racional (x,ys,zs) = reducido (a,b) where 
    a = digitosAnumero (x:ys++zs) - digitosAnumero (x:ys)
    b = digitosAnumero (replicate (length zs) 9 ++ replicate (length ys) 0)

-- ---------------------------------------------------------------------
-- § Propiedades                                                      --
-- ---------------------------------------------------------------------

-- La 1ª propiedad es
prop_RacionalDecimal :: Racional -> Property
prop_RacionalDecimal (a,b) =
    a >= 0 && b > 0  ==>
    racional (decimal (a,b)) == reducido (a,b)

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

-- En lugar de reducido se puede usar un generador de números racionales
numeroRacional :: Gen Racional
numeroRacional = do a <- arbitrary
                    b <- arbitrary
                    return (reducido (abs a, 1+ abs b))

-- La propiedad es
prop_RacionalDecimal2 :: Property
prop_RacionalDecimal2 =
    forAll numeroRacional
           (\(a,b) -> racional (decimal (a,b)) == (a,b))

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

-- Para la 2ª propiedad se define un generador de números decimales
numeroDecimal :: Gen Decimal
numeroDecimal = do a <- arbitrary
                   b <- arbitrary
                   return (decimal (abs a, 1+ abs b))

-- La 2ª propiedad es
prop_DecimalRacional :: Property
prop_DecimalRacional =
    forAll numeroDecimal 
           (\(x,ys,zs) -> decimal (racional (x,ys,zs)) == (x,ys,zs))

-- La comprobación es
--   ghci> quickCheck prop_DecimalRacional
--   +++ 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

169

170

171

172

173

174

175

176

177

178

179

-- ---------------------------------------------------------------------

-- § Librerías auxiliares --

-- ---------------------------------------------------------------------

import Data.List

import Test.QuickCheck

-- ---------------------------------------------------------------------

-- § Tipos --

-- ---------------------------------------------------------------------

-- Un número racional se representa por un par de enteros, donde el

-- primer elemento es el numerador y el segundo el denominador.

type Racional = (Integer,Integer)

-- Un número decimal es una terna donde el primer elemento es la parte

-- entera, el segundo es el anteperíodo y el tercero es el período.

type Decimal = (Integer,[Integer],[Integer])

-- ---------------------------------------------------------------------

-- § De racional a decimal --

-- ---------------------------------------------------------------------

-- (mayorExponente a x) es el mayor n tal que a^n divide a x. Por ejemplo,

-- mayorExponente 2 40 == 3

-- mayorExponente 5 40 == 1

-- mayorExponente 5 41 == 0

mayorExponente :: Integer -> Integer -> Integer

mayorExponente a x = head [n | n <- [1..], mod x (a^n) /= 0] - 1

-- (longitudAnteperiodo n) es la longitud del anteperíodo de 1/n (y de

-- cualquier fracción irreducible de denominador n). Por ejemplo,

-- longitudAnteperiodo 2 == 1

-- longitudAnteperiodo 3 == 0

-- longitudAnteperiodo 14 == 1

longitudAnteperiodo :: Integer -> Integer

longitudAnteperiodo n = max (mayorExponente 2 n) (mayorExponente 5 n)

-- (longitudPeriodo n) es la longitud del período de 1/n (y de

-- cualquier fracción irreducible de denominador n). Por ejemplo,

-- longitudPeriodo 2 == 0

-- longitudPeriodo 3 == 1

-- longitudPeriodo 7 == 6

longitudPeriodo :: Integer -> Integer

longitudPeriodo n

| m == 1 = 0

| otherwise = head [k | k <- [1..], (10^k) `mod` m == 1]

where m = n `div` (2^(mayorExponente 2 n) * 5^(mayorExponente 5 n))

-- (expansionDec (x,y)) es la expansión decimal de x/y. Por ejemplo,

-- take 10 (expansionDec (1,4)) == [0,2,5]

-- take 10 (expansionDec (1,7)) == [0,1,4,2,8,5,7,1,4,2]

-- take 12 (expansionDec (90,7)) == [12,8,5,7,1,4,2,8,5,7,1,4]

-- take 12 (expansionDec (23,14)) == [1,6,4,2,8,5,7,1,4,2,8,5]

expansionDec :: Racional -> [Integer]

expansionDec (x,y)

| r == 0 = [q]

| otherwise = q : expansionDec (r*10,y)

where (q,r) = quotRem x y

-- (parteEntera (a,b)) es la parte entera de a/b. Por ejemplo,

-- parteEntera (125,3) == 41

parteEntera :: Racional -> Integer

parteEntera (a,b) = a `div` b

-- (antePeriodo (a,b)) es el anteperíodo de a/b; es decir, la lista de

-- cifras que se encuentran entre la parte entera y el primer período de

-- a/b. Por ejemplo,

-- antePeriodo (23,14) == [6]

-- antePeriodo (1,5) == [2]

antePeriodo :: Racional -> [Integer]

antePeriodo (a,b) = genericTake s (tail xs)

where xs = expansionDec (a,b)

s = longitudAnteperiodo b

-- (periodo (a,b)) es el período de a/b; es decir, la lista de cifras que

-- se encuentran entre la parte entera y el primer período de a/b. Por

-- ejemplo,

-- periodo (1,3) == [3]

-- periodo (1,5) == []

-- periodo (1,7) == [1,4,2,8,5,7]

-- periodo (23,14) == [4,2,8,5,7,1]

-- concatMap show $ periodo (1,29) == "0344827586206896551724137931"

periodo :: Racional -> [Integer]

periodo (a,b) = genericTake t (genericDrop (1+s) xs)

where xs = expansionDec (a,b)

s = longitudAnteperiodo b

t = longitudPeriodo b

-- (reducido (a,b)) es la forma reducida del número racional a/b. Por

-- ejemplo,

-- reducido (40,80) == (1,2)

reducido :: Racional -> Racional

reducido (a,b) = (a `div` m, b `div` m)

where m = gcd a b

-- (decimal (x,y)) es la forma decimal de x/y; es decir, la terna

-- formada por la parte entera, la parte decimal pura y la parte decimal

-- periódica. Por ejemplo,

-- decimal (1,4) == (0,[2,5],[])

-- decimal (1,3) == (0,[],[3])

-- decimal (23,14) == (1,[6],[4,2,8,5,7,1])

decimal :: Racional -> Decimal

decimal (a,b) = (parteEntera r, antePeriodo r, periodo r)

where r = reducido (a,b)

-- ---------------------------------------------------------------------

-- § De decimal a racional --

-- ---------------------------------------------------------------------

-- (digitosAnumero xs) es el número correspondiente a la lista de

-- dígitos xs. Por ejemplo,

-- digitosAnumero [3,2,5] == 325

digitosAnumero :: [Integer] -> Integer

digitosAnumero = foldl (\a y -> 10*a+y) 0

-- 2ª definición de digitosAnumero (por comprensión)

digitosAnumero2 :: [Integer] -> Integer

digitosAnumero2 xs = sum [x*(10^i) | (x,i) <- zip xs [n-1,n-2..0]]

where n = length xs

-- (racional (x,ys,zs)) es el número racional cuya representación

-- decimal es (x,ys,zs). Por ejemplo,

-- racional (0,[2,5],[]) == (1,4)

-- racional (0,[],[3]) == (1,3)

-- racional (1,[6],[4,2,8,5,7,1]) == (23,14)

racional :: Decimal -> Racional

racional (x,ys,[]) = reducido (a,b) where

a = digitosAnumero (x:ys)

b = 10^(length ys)

racional (x,ys,zs) = reducido (a,b) where

a = digitosAnumero (x:ys++zs) - digitosAnumero (x:ys)

b = digitosAnumero (replicate (length zs) 9 ++ replicate (length ys) 0)

-- ---------------------------------------------------------------------

-- § Propiedades --

-- ---------------------------------------------------------------------

-- La 1ª propiedad es

prop_RacionalDecimal :: Racional -> Property

prop_RacionalDecimal (a,b) =

a >= 0 && b > 0 ==>

racional (decimal (a,b)) == reducido (a,b)

-- La comprobación es

-- ghci> quickCheck prop_RacionalDecimal

-- +++ OK, passed 100 tests.

-- En lugar de reducido se puede usar un generador de números racionales

numeroRacional :: Gen Racional

numeroRacional = do a <- arbitrary

b <- arbitrary

return (reducido (abs a, 1+ abs b))

-- La propiedad es

prop_RacionalDecimal2 :: Property

prop_RacionalDecimal2 =

forAll numeroRacional

(\(a,b) -> racional (decimal (a,b)) == (a,b))

-- La comprobación es

-- ghci> quickCheck prop_RacionalDecimal2

-- +++ OK, passed 100 tests.

-- Para la 2ª propiedad se define un generador de números decimales

numeroDecimal :: Gen Decimal

numeroDecimal = do a <- arbitrary

b <- arbitrary

return (decimal (abs a, 1+ abs b))

-- La 2ª propiedad es

prop_DecimalRacional :: Property

prop_DecimalRacional =

forAll numeroDecimal

(\(x,ys,zs) -> decimal (racional (x,ys,zs)) == (x,ys,zs))

-- La comprobación es

-- ghci> quickCheck prop_DecimalRacional

-- +++ OK, passed 100 tests.

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.

Listas hermanadas

PorJosé A. Alonso 21-diciembre-201528-diciembre-2015

Una lista hermanada es una lista de números estrictamente positivos en la que cada elemento tiene algún factor primo en común con el siguiente, en caso de que exista, o alguno de los dos es un 1. Por ejemplo,

[2,6,3,9,1,5] es una lista hermanada pues 2 y 6 tienen un factor en común (2); 6 y 3 tienen un factor en común (3); 3 y 9 tienen un factor en común (3); de 9 y 1 uno es el número 1; y de 1 y 5 uno es el número 1.
[2,3,5] no es una lista hermanada pues 2 y 3 no tienen ningún factor primo en común.

Definir la función

   hermanada :: [Int] -> Bool

1	hermanada :: [Int] -> Bool

tal que (hermanada xs) se verifica si la lista xs es hermanada según la definición anterior. Por ejemplo,

   hermanada [2,6,3,9,1,5]   ==  True
   hermanada [2,3,5]         ==  False
   hermanada [2,4..1000000]  ==  True

hermanada [2,6,3,9,1,5] == True

hermanada [2,3,5] == False

hermanada [2,4..1000000] == True

Nota: Este ejercicio es parte del examen del grupo 3 del 2 de diciembre.

Soluciones

-- 1ª definición (por comprensión)
hermanada1 :: [Int] -> Bool
hermanada1 xs = and [hermanos p | p <- zip xs (tail xs)]

-- (hermanos (x,y)) se verifica si x e y son hermanos; es decir, alguno es
-- igual a 1 o tienen algún factor primo en común
hermanos :: (Int, Int) -> Bool
hermanos (x,y) = x == 1 || y == 1 || gcd x y /= 1

-- 2ª definición (con all)
hermanada2 :: [Int] -> Bool
hermanada2 xs = all hermanos (zip xs (tail xs))

-- 3ª definición (por recursión)
hermanada3 :: [Int] -> Bool
hermanada3 (x1:x:xs) = hermanos (x1,x) && hermanada3 (x:xs)
hermanada3 _          = True

-- 4ª definición (por plegado)
hermanada4 :: [Int] -> Bool
hermanada4 xs =
    foldl (\ws p -> hermanos p && ws) True (zip xs (tail xs))

-- Comparación de eficiencia
--    λ> hermanada1 [2,4..1000000]
--    True
--    (2.33 secs, 476,586,552 bytes)
--    λ> hermanada2 [2,4..1000000]
--    True
--    (1.80 secs, 422,879,072 bytes)
--    λ> hermanada3 [2,4..1000000]
--    True
--    (2.58 secs, 477,251,896 bytes)
--    λ> hermanada4 [2,4..1000000]
--    True
--    (2.36 secs, 440,047,520 bytes)

-- 1ª definición (por comprensión)

hermanada1 :: [Int] -> Bool

hermanada1 xs = and [hermanos p | p <- zip xs (tail xs)]

-- (hermanos (x,y)) se verifica si x e y son hermanos; es decir, alguno es

-- igual a 1 o tienen algún factor primo en común

hermanos :: (Int, Int) -> Bool

hermanos (x,y) = x == 1 || y == 1 || gcd x y /= 1

-- 2ª definición (con all)

hermanada2 :: [Int] -> Bool

hermanada2 xs = all hermanos (zip xs (tail xs))

-- 3ª definición (por recursión)

hermanada3 :: [Int] -> Bool

hermanada3 (x1:x:xs) = hermanos (x1,x) && hermanada3 (x:xs)

hermanada3 _ = True

-- 4ª definición (por plegado)

hermanada4 :: [Int] -> Bool

hermanada4 xs =

foldl (\ws p -> hermanos p && ws) True (zip xs (tail xs))

-- Comparación de eficiencia

-- λ> hermanada1 [2,4..1000000]

-- True

-- (2.33 secs, 476,586,552 bytes)

-- λ> hermanada2 [2,4..1000000]

-- True

-- (1.80 secs, 422,879,072 bytes)

-- λ> hermanada3 [2,4..1000000]

-- True

-- (2.58 secs, 477,251,896 bytes)

-- λ> hermanada4 [2,4..1000000]

-- True

-- (2.36 secs, 440,047,520 bytes)

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

Conjunto de primos relativos

Soluciones

Nuevas soluciones

La sucesión ECG

Soluciones

Acotación del primorial

Soluciones

Otras soluciones

Pensamiento

Enteros como sumas de tres coprimos.

Soluciones

Referencias

Pensamiento

La sucesión ECG

Soluciones

Pensamiento

Entre dos conjuntos

Soluciones

Referencia

Pensamiento

Máximo común divisor de x e y veces n

Soluciones

Inverso multiplicativo modular

Soluciones

Solución en Maxima

Referencia

Representación decimal de números racionales

Soluciones

Listas hermanadas

Soluciones

Sucesión de números parientes

Soluciones

Entradas recientes

Buscar

Correo electrónico