map

Teorema de Liouville sobre listas CuCu

PorJosé A. Alonso 3-enero-202010-enero-2020

Una lista CuCu es una lista de números enteros positivos tales que la suma de sus Cubos es igual al Cuadrado de su suma. Por ejemplo, [1, 2, 3, 2, 4, 6] es una lista CuCu ya que

   1³ + 2³ + 3³ + 2³ + 4³ + 6³ = (1 + 2 + 3 + 2 + 4 + 6)²

1	1³ + 2³ + 3³ + 2³ + 4³ + 6³ = (1 + 2 + 3 + 2 + 4 + 6)²

La lista de Liouville correspondiente al número entero positivo n es la lista formada por el número de divisores de cada divisor de n. Por ejemplo, para el número 20 se tiene que sus divisores son

   1, 2, 4, 5, 10, 20

1	1, 2, 4, 5, 10, 20

puesto que el número de sus divisores es

El 1 tiene 1 divisor (el 1 solamente).
El 2 tiene 2 divisores (el 1 y el 2).
El 4 tiene 3 divisores (el 1, el 2 y el 4).
El 5 tiene 2 divisores (el 1 y el 5).
El 10 tiene 4 divisores (el 1, el 2, el 5 y el 10).
El 20 tiene 6 divisores (el 1, el 2, el 4, el 5, el 10 y el 20).

la lista de Liouville de 20 es [1, 2, 3, 2, 4, 6] que, como se comentó anteriormente, es una lista CuCu.

El teorema de Lioville afirma que todas las lista de Lioville son CuCu.

Definir las funciones

   esCuCu :: [Integer] -> Bool
   liouville :: Integer -> [Integer]

1 2	esCuCu :: [Integer] -> Bool liouville :: Integer -> [Integer]

tales que

(esCuCu xs) se verifica si la lista xs es CuCu; es decir, la suma de los cubos de sus elementos es igual al cuadrado de su suma. Por ejemplo,

     esCuCu [1,2,3]        ==  True
     esCuCu [1,2,3,2]      ==  False
     esCuCu [1,2,3,2,4,6]  ==  True

esCuCu [1,2,3] == True

esCuCu [1,2,3,2] == False

esCuCu [1,2,3,2,4,6] == True

(liouville n) es la lista de Lioville correspondiente al número n. Por ejemplo,

     liouville 20  ==  [1,2,3,2,4,6]
     liouville 60  ==  [1,2,2,3,2,4,4,6,4,6,8,12]
     length (liouville (product [1..25]))  ==  340032

liouville 20 == [1,2,3,2,4,6]

liouville 60 == [1,2,2,3,2,4,4,6,4,6,8,12]

length (liouville (product [1..25])) == 340032

Comprobar con QuickCheck

que para todo entero positivo n, (liouville (2^n)) es la lista [1,2,3,…,n+1] y
el teorema de Lioville; es decir, para todo entero positivo n, (liouville n) es una lista CuCu.

Nota: Este ejercicio está basado en Cómo generar conjuntos CuCu de Gaussianos.

Soluciones

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

esCuCu :: [Integer] -> Bool
esCuCu xs = sum (map (^3) xs) == (sum xs)^2

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

liouville :: Integer -> [Integer]
liouville n = map numeroDivisores (divisores n)

-- (divisores x) es el conjunto de divisores de los x. Por ejemplo, 
--   divisores 30  ==  [1,2,3,5,6,10,15,30]
divisores :: Integer -> [Integer]
divisores n = [x | x <- [1..n], n `mod` x == 0]

-- (numeroDivisores x) es el número de divisores de x. Por ejemplo, 
--    numeroDivisores 12  ==  6
--    numeroDivisores 25  ==  3
numeroDivisores :: Integer -> Integer
numeroDivisores n = genericLength (divisores n) 

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

liouville2 :: Integer -> [Integer]
liouville2 n = map numeroDivisores2 (divisores2 n)

-- Se usan las funciones
-- + divisores de "Conjunto de divisores" http://bit.ly/2OtbFIj
-- + numeroDivisores de "Número de divisores" http://bit.ly/2DgVh74

-- (divisores2 x) es el conjunto de divisores de los x. Por ejemplo, 
--   divisores2 30  ==  [1,2,3,5,6,10,15,30]
divisores2 :: Integer -> [Integer]
divisores2 = sort
           . map (product . concat)
           . sequence
           . map inits
           . group
           . primeFactors

-- (numeroDivisores2 x) es el número de divisores de x. Por ejemplo, 
--    numeroDivisores2 12  ==  6
--    numeroDivisores2 25  ==  3
numeroDivisores2 :: Integer -> Integer
numeroDivisores2 =
  product . map ((+1) . genericLength) . group . primeFactors

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

-- La comparación es
--    λ> length (liouville (product [1..11]))
--    540
--    (13.66 secs, 7,983,550,640 bytes)
--    λ> length (liouville2 (product [1..11]))
--    540
--    (0.01 secs, 1,255,328 bytes)

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

-- La propiedad es
prop_Liouville :: Integer -> Property
prop_Liouville n =
  n > 0 ==> liouville2 (2^n) == [1..n+1]

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

-- Teorema de Liouville
-- ====================

-- La propiedad es
teorema_Liouville :: Integer -> Property
teorema_Liouville n =
  n > 0 ==> esCuCu (liouville n)

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

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

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

esCuCu :: [Integer] -> Bool

esCuCu xs = sum (map (^3) xs) == (sum xs)^2

-- 1ª definición de liouville

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

liouville :: Integer -> [Integer]

liouville n = map numeroDivisores (divisores n)

-- (divisores x) es el conjunto de divisores de los x. Por ejemplo,

-- divisores 30 == [1,2,3,5,6,10,15,30]

divisores :: Integer -> [Integer]

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

-- (numeroDivisores x) es el número de divisores de x. Por ejemplo,

-- numeroDivisores 12 == 6

-- numeroDivisores 25 == 3

numeroDivisores :: Integer -> Integer

numeroDivisores n = genericLength (divisores n)

-- 2ª definición de liouville

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

liouville2 :: Integer -> [Integer]

liouville2 n = map numeroDivisores2 (divisores2 n)

-- Se usan las funciones

-- + divisores de "Conjunto de divisores" http://bit.ly/2OtbFIj

-- + numeroDivisores de "Número de divisores" http://bit.ly/2DgVh74

-- (divisores2 x) es el conjunto de divisores de los x. Por ejemplo,

-- divisores2 30 == [1,2,3,5,6,10,15,30]

divisores2 :: Integer -> [Integer]

divisores2 = sort

. map (product . concat)

. sequence

. map inits

. group

. primeFactors

-- (numeroDivisores2 x) es el número de divisores de x. Por ejemplo,

-- numeroDivisores2 12 == 6

-- numeroDivisores2 25 == 3

numeroDivisores2 :: Integer -> Integer

numeroDivisores2 =

product . map ((+1) . genericLength) . group . primeFactors

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (liouville (product [1..11]))

-- 540

-- (13.66 secs, 7,983,550,640 bytes)

-- λ> length (liouville2 (product [1..11]))

-- 540

-- (0.01 secs, 1,255,328 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_Liouville :: Integer -> Property

prop_Liouville n =

n > 0 ==> liouville2 (2^n) == [1..n+1]

-- La comprobación es

-- λ> quickCheck prop_Liouville

-- +++ OK, passed 100 tests.

-- Teorema de Liouville

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

-- La propiedad es

teorema_Liouville :: Integer -> Property

teorema_Liouville n =

n > 0 ==> esCuCu (liouville n)

-- La comprobación es

-- λ> quickCheck teorema_Liouville

-- +++ OK, passed 100 tests.

Pensamiento

¡Oh, tarde viva y quieta
que opuso al panta rhei su nada corre.

Antonio Machado

Enumeración de conjuntos finitos de naturales

PorJosé A. Alonso 25-diciembre-20191-enero-2020

Los conjuntos finitos de números naturales se pueden enumerar como sigue

    0: []
    1: [0]
    2: [1]
    3: [1,0]
    4: [2]
    5: [2,0]
    6: [2,1]
    7: [2,1,0]
    8: [3]
    9: [3,0]
   10: [3,1]
   11: [3,1,0]
   12: [3,2]
   13: [3,2,0]
   14: [3,2,1]
   15: [3,2,1,0]
   16: [4]
   17: [4,0]
   18: [4,1]
   19: [4,1,0]

0: []

1: [0]

2: [1]

3: [1,0]

4: [2]

5: [2,0]

6: [2,1]

7: [2,1,0]

8: [3]

9: [3,0]

10: [3,1]

11: [3,1,0]

12: [3,2]

13: [3,2,0]

14: [3,2,1]

15: [3,2,1,0]

16: [4]

17: [4,0]

18: [4,1]

19: [4,1,0]

en la que los elementos están ordenados de manera decreciente.

Definir la constante

   enumeracionCFN :: [[Integer]]

1	enumeracionCFN :: [[Integer]]

tal que sus elementos son los conjuntos de los números naturales con la ordenación descrita anteriormente. Por ejemplo,

   λ> take 20 enumeracionCFN
   [[],
    [0],
    [1],[1,0],
    [2],[2,0],[2,1],[2,1,0],
    [3],[3,0],[3,1],[3,1,0],[3,2],[3,2,0],[3,2,1],[3,2,1,0],
    [4],[4,0],[4,1],[4,1,0]]

λ> take 20 enumeracionCFN

[[],

[0],

[1],[1,0],

[2],[2,0],[2,1],[2,1,0],

[3],[3,0],[3,1],[3,1,0],[3,2],[3,2,0],[3,2,1],[3,2,1,0],

[4],[4,0],[4,1],[4,1,0]]

Comprobar con QuickCheck que

si (xs,ys) es un par de elementos consecutivos de enumeracionCFN, entonces xs < ys;
todo conjunto finito de números naturales, representado por una lista decreciente, está en enumeracionCFN.

Soluciones

import Data.List (genericLength, nub, sort)
import Test.QuickCheck

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

enumeracionCFN :: [[Integer]]
enumeracionCFN = concatMap enumeracionCFNHasta [0..]

-- (enumeracionCFNHasta n) es la lista de conjuntos con la enumeración
-- anterior cuyo primer elemento es n. Por ejemplo,
--    λ> enumeracionCFNHasta 1
--    [[1],[1,0]]
--    λ> enumeracionCFNHasta 2
--    [[2],[2,0],[2,1],[2,1,0]]
--    λ> enumeracionCFNHasta 3
--    [[3],[3,0],[3,1],[3,1,0],[3,2],[3,2,0],[3,2,1],[3,2,1,0]]
enumeracionCFNHasta :: Integer -> [[Integer]]
enumeracionCFNHasta 0 = [[],[0]]
enumeracionCFNHasta n =
  [n:xs | k <- [0..n-1], xs <- enumeracionCFNHasta k]

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

enumeracionCFN2 :: [[Integer]]
enumeracionCFN2 = [] : aux 0 [[]]
  where aux n xs = yss ++ aux (n+1) (xs ++ yss)
          where yss = map (n:) xs

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

enumeracionCFN3 :: [[Integer]]
enumeracionCFN3 = map conjunto [0..]

-- (conjunto n) es el conjunto en la posición n. Por ejemplo,
--   conjunto 6  ==  [2,1]
conjunto :: Integer -> [Integer]
conjunto n = reverse [x | (x,y) <- zip [0..] (binario n), y == 1]

-- (binario n) es la representación binarioa del número n (en orden
-- inverso). Por ejemplo,
--   binario 6  ==  [0,1,1]
binario :: Integer -> [Integer]
binario 0 = [0]
binario 1 = [1]
binario n = n `mod` 2 : binario (n `div` 2)

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

--    λ> enumeracionCFN !! (4*10^5)
--    [18,17,12,11,9,7]
--    (1.18 secs, 576,924,344 bytes)
--    λ> enumeracionCFN2 !! (4*10^5)
--    [18,17,12,11,9,7]
--    (0.10 secs, 72,399,784 bytes)
--    λ> enumeracionCFN3 !! (4*10^5)
--    [18,17,12,11,9,7]
--    (0.07 secs, 64,123,952 bytes)
--
--    λ> enumeracionCFN2 !! (6*10^6)
--    [22,20,19,17,16,15,11,10,8,7]
--    (1.25 secs, 1,082,690,216 bytes)
--    λ> enumeracionCFN3 !! (6*10^6)
--    [22,20,19,17,16,15,11,10,8,7]
--    (0.38 secs, 960,134,256 bytes)

-- Propiedades
-- ===========

-- La primera propiedad es
prop_enumeracionCFN :: Int -> Property
prop_enumeracionCFN n =
  n >= 0 ==> xs < ys
  where (xs:ys:_) = drop n enumeracionCFN

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

-- La segunda propiedad es
prop_completa :: [Integer] -> Bool
prop_completa xs =
  xs' `elem` enumeracionCFN
  where xs' = reverse (sort (nub (map abs xs)))

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=15}) prop_completa
--    +++ OK, passed 100 tests.

import Data.List (genericLength, nub, sort)

import Test.QuickCheck

-- 1ª solución

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

enumeracionCFN :: [[Integer]]

enumeracionCFN = concatMap enumeracionCFNHasta [0..]

-- (enumeracionCFNHasta n) es la lista de conjuntos con la enumeración

-- anterior cuyo primer elemento es n. Por ejemplo,

-- λ> enumeracionCFNHasta 1

-- [[1],[1,0]]

-- λ> enumeracionCFNHasta 2

-- [[2],[2,0],[2,1],[2,1,0]]

-- λ> enumeracionCFNHasta 3

-- [[3],[3,0],[3,1],[3,1,0],[3,2],[3,2,0],[3,2,1],[3,2,1,0]]

enumeracionCFNHasta :: Integer -> [[Integer]]

enumeracionCFNHasta 0 = [[],[0]]

enumeracionCFNHasta n =

[n:xs | k <- [0..n-1], xs <- enumeracionCFNHasta k]

-- 2ª solución

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

enumeracionCFN2 :: [[Integer]]

enumeracionCFN2 = [] : aux 0 [[]]

where aux n xs = yss ++ aux (n+1) (xs ++ yss)

where yss = map (n:) xs

-- 3ª solución

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

enumeracionCFN3 :: [[Integer]]

enumeracionCFN3 = map conjunto [0..]

-- (conjunto n) es el conjunto en la posición n. Por ejemplo,

-- conjunto 6 == [2,1]

conjunto :: Integer -> [Integer]

conjunto n = reverse [x | (x,y) <- zip [0..] (binario n), y == 1]

-- (binario n) es la representación binarioa del número n (en orden

-- inverso). Por ejemplo,

-- binario 6 == [0,1,1]

binario :: Integer -> [Integer]

binario 0 = [0]

binario 1 = [1]

binario n = n `mod` 2 : binario (n `div` 2)

-- Comparación de eficiencia

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

-- λ> enumeracionCFN !! (4*10^5)

-- [18,17,12,11,9,7]

-- (1.18 secs, 576,924,344 bytes)

-- λ> enumeracionCFN2 !! (4*10^5)

-- [18,17,12,11,9,7]

-- (0.10 secs, 72,399,784 bytes)

-- λ> enumeracionCFN3 !! (4*10^5)

-- [18,17,12,11,9,7]

-- (0.07 secs, 64,123,952 bytes)

-- λ> enumeracionCFN2 !! (6*10^6)

-- [22,20,19,17,16,15,11,10,8,7]

-- (1.25 secs, 1,082,690,216 bytes)

-- λ> enumeracionCFN3 !! (6*10^6)

-- [22,20,19,17,16,15,11,10,8,7]

-- (0.38 secs, 960,134,256 bytes)

-- Propiedades

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

-- La primera propiedad es

prop_enumeracionCFN :: Int -> Property

prop_enumeracionCFN n =

n >= 0 ==> xs < ys

where (xs:ys:_) = drop n enumeracionCFN

-- La comprobación es

-- λ> quickCheck prop_enumeracionCFN

-- +++ OK, passed 100 tests.

-- La segunda propiedad es

prop_completa :: [Integer] -> Bool

prop_completa xs =

xs' `elem` enumeracionCFN

where xs' = reverse (sort (nub (map abs xs)))

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=15}) prop_completa

-- +++ OK, passed 100 tests.

Pensamiento

Junto al agua fría,
en la senda clara,
sombra dará algún día,
ese arbolillo en que nadie repara.

Antonio Machado

Cálculo de pi usando la fórmula de Vieta

PorJosé A. Alonso 16-diciembre-201923-diciembre-2019

La fórmula de Vieta para el cálculo de pi es la siguiente

Definir las funciones

   aproximacionPi :: Int -> Double
   errorPi :: Double -> Int

1 2	aproximacionPi :: Int -> Double errorPi :: Double -> Int

tales que

(aproximacionPi n) es la aproximación de pi usando n factores de la fórmula de Vieta. Por ejemplo,

     aproximacionPi  5  ==  3.140331156954753
     aproximacionPi 10  ==  3.1415914215112
     aproximacionPi 15  ==  3.141592652386592
     aproximacionPi 20  ==  3.1415926535886207
     aproximacionPi 25  ==  3.141592653589795

aproximacionPi 5 == 3.140331156954753

aproximacionPi 10 == 3.1415914215112

aproximacionPi 15 == 3.141592652386592

aproximacionPi 20 == 3.1415926535886207

aproximacionPi 25 == 3.141592653589795

(errorPi x) es el menor número de factores de la fórmula de Vieta necesarios para obtener pi con un error menor que x. Por ejemplo,

     errorPi 0.1        ==  2
     errorPi 0.01       ==  4
     errorPi 0.001      ==  6
     errorPi 0.0001     ==  7
     errorPi 1e-4       ==  7
     errorPi 1e-14      ==  24
     pi                 ==  3.141592653589793
     aproximacionPi 24  ==  3.1415926535897913

errorPi 0.1 == 2

errorPi 0.01 == 4

errorPi 0.001 == 6

errorPi 0.0001 == 7

errorPi 1e-4 == 7

errorPi 1e-14 == 24

pi == 3.141592653589793

aproximacionPi 24 == 3.1415926535897913

Soluciones

-- 1ª definición de aproximacionPi
aproximacionPi :: Int -> Double
aproximacionPi n = product [2 / aux x | x <- [0..n]]
  where
    aux 0 = 1
    aux 1 = sqrt 2
    aux n = sqrt (2 + aux (n-1))

-- 2ª definición de aproximacionPi
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = product [2/x | x <- 1 : xs] 
  where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)

-- 3ª definición de aproximaxionPi
aproximacionPi3 :: Int -> Double
aproximacionPi3 n =  product (2 : take n (map (2/) xs))
  where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]

-- 1ª definición de errorPi
errorPi :: Double -> Int
errorPi x = head [n | n <- [1..]
                    , abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi
errorPi2 :: Double -> Int
errorPi2 x = until aceptable (+1) 1
  where aceptable n = abs (pi - aproximacionPi n) < x

-- 1ª definición de aproximacionPi

aproximacionPi :: Int -> Double

aproximacionPi n = product [2 / aux x | x <- [0..n]]

where

aux 0 = 1

aux 1 = sqrt 2

aux n = sqrt (2 + aux (n-1))

-- 2ª definición de aproximacionPi

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = product [2/x | x <- 1 : xs]

where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)

-- 3ª definición de aproximaxionPi

aproximacionPi3 :: Int -> Double

aproximacionPi3 n = product (2 : take n (map (2/) xs))

where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]

-- 1ª definición de errorPi

errorPi :: Double -> Int

errorPi x = head [n | n <- [1..]

, abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi

errorPi2 :: Double -> Int

errorPi2 x = until aceptable (+1) 1

where aceptable n = abs (pi - aproximacionPi n) < x

Pensamiento

El tiempo que la barba me platea,
cavó mis ojos y agrandó mi frente,
va siendo en mi recuerdo transparente,
y mientras más al fondo, más clarea.

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

Primos o cuadrados de primos

PorJosé A. Alonso 6-diciembre-201917-diciembre-2019

Definir la constante

   primosOcuadradosDePrimos :: [Integer]

1	primosOcuadradosDePrimos :: [Integer]

cuyos elementos son los número primos o cuadrados de primos. Por ejemplo,

   λ> take 20 primosOcuadradosDePrimos
   [2,3,4,5,7,9,11,13,17,19,23,25,29,31,37,41,43,47,49,53]
   λ> primosOcuadradosDePrimos !! (10^6)
   15476729

λ> take 20 primosOcuadradosDePrimos

[2,3,4,5,7,9,11,13,17,19,23,25,29,31,37,41,43,47,49,53]

λ> primosOcuadradosDePrimos !! (10^6)

15476729

Comprobar con QuickCheck que las lista primosOcuadradosDePrimos y unifactorizables (definida en el ejercicio anterior) son iguales.

Soluciones

import Test.QuickCheck
import Data.Numbers.Primes (primeFactors, primes)

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

primosOcuadradosDePrimos :: [Integer]
primosOcuadradosDePrimos =
  filter esPrimoOcuadradoDePrimo [2..]

esPrimoOcuadradoDePrimo :: Integer -> Bool
esPrimoOcuadradoDePrimo n = aux xs
  where xs = primeFactors n
        aux [_]   = True
        aux [x,y] = x == y
        aux _     = False

-- 2ª solución
-- ===========
  
primosOcuadradosDePrimos2 :: [Integer]
primosOcuadradosDePrimos2 = mezcla primes (map (^2) primes)

mezcla :: Ord a => [a] -> [a] -> [a]
mezcla (x:xs) (y:ys) | x < y     = x : mezcla xs (y:ys)
                     | otherwise = y : mezcla (x:xs) ys

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

--    λ> primosOcuadradosDePrimos !! (2*10^4)
--    223589
--    (3.08 secs, 9,829,725,096 bytes)
--    λ> primosOcuadradosDePrimos2 !! (2*10^4)
--    223589
--    (0.04 secs, 73,751,888 bytes)
--
--    λ> primosOcuadradosDePrimos2 !! (5*10^5)
--    7362497
--    (1.29 secs, 3,192,803,040 bytes)

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

-- La propiedad es
prop_equivalencia :: Int -> Property
prop_equivalencia n =
  n >= 0 ==> primosOcuadradosDePrimos2 !! n == unifactorizables !! n

--  unifactorizables es la lísta de los números enteros mayores que 1
--  que se pueden escribir sólo de una forma única como producto de
--  enteros distintos mayores que uno. Por ejemplo,
--     λ> take 20 unifactorizables
--     [2,3,4,5,7,9,11,13,17,19,23,25,29,31,37,41,43,47,49,53]
--     λ> unifactorizables !! 300
--     1873
unifactorizables :: [Integer]
unifactorizables =
  [n | n <- [2..]
     , length (sublistasConProducto n [2..n]) == 1]

-- (sublistasConProducto n xs) es la lista de las sublistas de la
-- lista ordenada estrictamente creciente xs (cuyos elementos son
-- enteros mayores que 1) cuyo producto es el número entero n (con n
-- mayor que 1). Por ejemplo,  
--    λ> sublistasConProducto 72 [2,3,4,5,6,7,9,10,16]
--    [[2,4,9],[3,4,6]]
--    λ> sublistasConProducto 720 [2,3,4,5,6,7,9,10,16]
--    [[2,3,4,5,6],[2,4,9,10],[3,4,6,10],[5,9,16]]
--    λ> sublistasConProducto 2 [4,7]
--    []
sublistasConProducto :: Integer -> [Integer] -> [[Integer]]
sublistasConProducto _ [] = []
sublistasConProducto n (x:xs)
  | x > n     = []
  | x == n    = [[x]]
  | r == 0    = map (x:) (sublistasConProducto q xs)
                ++ sublistasConProducto n xs
  | otherwise = sublistasConProducto n xs
  where (q,r) = quotRem n x

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

import Test.QuickCheck

import Data.Numbers.Primes (primeFactors, primes)

-- 1ª solución

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

primosOcuadradosDePrimos :: [Integer]

primosOcuadradosDePrimos =

filter esPrimoOcuadradoDePrimo [2..]

esPrimoOcuadradoDePrimo :: Integer -> Bool

esPrimoOcuadradoDePrimo n = aux xs

where xs = primeFactors n

aux [_] = True

aux [x,y] = x == y

aux _ = False

-- 2ª solución

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

primosOcuadradosDePrimos2 :: [Integer]

primosOcuadradosDePrimos2 = mezcla primes (map (^2) primes)

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

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

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

-- Comparación de eficiencia

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

-- λ> primosOcuadradosDePrimos !! (2*10^4)

-- 223589

-- (3.08 secs, 9,829,725,096 bytes)

-- λ> primosOcuadradosDePrimos2 !! (2*10^4)

-- 223589

-- (0.04 secs, 73,751,888 bytes)

-- λ> primosOcuadradosDePrimos2 !! (5*10^5)

-- 7362497

-- (1.29 secs, 3,192,803,040 bytes)

-- Propiedad de equivalencia

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

-- La propiedad es

prop_equivalencia :: Int -> Property

prop_equivalencia n =

n >= 0 ==> primosOcuadradosDePrimos2 !! n == unifactorizables !! n

-- unifactorizables es la lísta de los números enteros mayores que 1

-- que se pueden escribir sólo de una forma única como producto de

-- enteros distintos mayores que uno. Por ejemplo,

-- λ> take 20 unifactorizables

-- [2,3,4,5,7,9,11,13,17,19,23,25,29,31,37,41,43,47,49,53]

-- λ> unifactorizables !! 300

-- 1873

unifactorizables :: [Integer]

unifactorizables =

[n | n <- [2..]

, length (sublistasConProducto n [2..n]) == 1]

-- (sublistasConProducto n xs) es la lista de las sublistas de la

-- lista ordenada estrictamente creciente xs (cuyos elementos son

-- enteros mayores que 1) cuyo producto es el número entero n (con n

-- mayor que 1). Por ejemplo,

-- λ> sublistasConProducto 72 [2,3,4,5,6,7,9,10,16]

-- [[2,4,9],[3,4,6]]

-- λ> sublistasConProducto 720 [2,3,4,5,6,7,9,10,16]

-- [[2,3,4,5,6],[2,4,9,10],[3,4,6,10],[5,9,16]]

-- λ> sublistasConProducto 2 [4,7]

-- []

sublistasConProducto :: Integer -> [Integer] -> [[Integer]]

sublistasConProducto _ [] = []

sublistasConProducto n (x:xs)

| x > n = []

| x == n = [[x]]

| r == 0 = map (x:) (sublistasConProducto q xs)

++ sublistasConProducto n xs

| otherwise = sublistasConProducto n xs

where (q,r) = quotRem n x

-- La comprobación es

-- λ> quickCheck prop_equivalencia

-- +++ OK, passed 100 tests.

Pensamiento

Despacito y buena letra: el hacer las cosas bien importa más que el hacerlas.

Antonio Machado

Enumeración de conjuntos finitos de naturales

Soluciones

Pensamiento

Cálculo de pi usando la fórmula de Vieta

Soluciones

Pensamiento

Primos o cuadrados de primos

Soluciones

Pensamiento

Entradas recientes

Buscar

Correo electrónico