diciembre 2019 – Exercitium

Derivada aritmética

PorJosé A. Alonso 31-diciembre-20197-enero-2020

La derivada aritmética es una función definida sobre los números naturales por analogía con la regla del producto para el cálculo de las derivadas usada en análisis.

Para un número natural n su derivada D(n) se define por

   D(0)  = 0
   D(1)  = 0
   D(p)  = 1, si p es primo
   D(ab) = D(a)b + aD(b) (regla de Leibniz para derivar productos)

D(0) = 0

D(1) = 0

D(p) = 1, si p es primo

D(ab) = D(a)b + aD(b) (regla de Leibniz para derivar productos)

Por ejemplo,

   D(6)  = D(2*3) = D(2)*3 + 2*D(3) = 1*3 + 2*1 =  5
   D(12) = D(2*6) = D(2)*6 + 2*D(6) = 1*6 + 2*5 = 16

1 2	D(6) = D(23) = D(2)3 + 2D(3) = 13 + 21 = 5 D(12) = D(26) = D(2)6 + 2D(6) = 16 + 25 = 16

Definir la función

   derivada :: Integer -> Integer

1	derivada :: Integer -> Integer

tal que (derivada n) es la derivada aritmética de n. Por ejemplo,

   derivada  6  ==  5
   derivada 12  ==  16
   maximum [derivada n | n <- [1..60000]]  ==  380928

derivada 6 == 5

derivada 12 == 16

maximum [derivada n | n <- [1..60000]] == 380928

Comprobar con QuickCheck que si x es un número entero positivo y su descomposición en factores primos es

   x = p(1)^e(1) + p(2)^e(2) +...+ p(n)^e(n)

1	x = p(1)^e(1) + p(2)^e(2) +...+ p(n)^e(n)

entonces la derivada de x es

   x * [e(1)/p(1) + e(2)/p(2) +...+ e(n)/p(n)]

1	x * [e(1)/p(1) + e(2)/p(2) +...+ e(n)/p(n)]

Nota: No usar en la definición la propiedad que hay que comprobar.

Soluciones

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

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

derivada :: Integer -> Integer
derivada 0 = 0
derivada 1 = 0
derivada n | esPrimo n = 1
           | otherwise = (derivada a) * b + a * (derivada b)
  where a = menorFactor n
        b = n `div` a

-- (esPrimo n) se verifica si n es primo. Por ejemplo,
--    esPrimo 5  ==  True
--    esPrimo 6  ==  False
esPrimo :: Integer -> Bool
esPrimo 0 = False
esPrimo 1 = False
esPrimo n = n == menorFactor n

-- (menorFactor n) es el menor divisor primo de n (con n >= 2). Por
-- ejemplo, 
--    menorFactor 6   ==  2
--    menorFactor 7   ==  7
--    menorFactor 15  ==  3
menorFactor :: Integer -> Integer
menorFactor n
  | even n = 2
  | otherwise = head [x | x <- [3,5..]
                        , n `mod` x == 0]

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

derivada2 :: Integer -> Integer
derivada2 0 = 0
derivada2 1 = 0
derivada2 n | isPrime n = 1
            | otherwise = (derivada2 a) * b + a * (derivada2 b)
  where (a:_) = primeFactors n
        b     = n `div` a

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

--    λ> maximum [derivada n | n <- [1..10000]]
--    53248
--    (1.59 secs, 1,091,452,552 bytes)
--    λ> maximum [derivada2 n | n <- [1..10000]]
--    53248
--    (0.17 secs, 457,819,120 bytes)

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

-- La propiedad es
prop_derivada :: Integer -> Property
prop_derivada x =
  x > 0 ==>
  derivada x == sum [(x * e) `div` p | (p,e) <- factorizacion x]

-- (factorizacion x) es la lista de las bases y exponentes de
-- la descomposición prima de x. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- Su comprobación es
--    λ> quickCheck prop_derivada
--    +++ OK, passed 100 tests.

import Data.List (genericLength, group)

import Data.Numbers.Primes (isPrime, primeFactors)

import Test.QuickCheck

-- 1ª solución

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

derivada :: Integer -> Integer

derivada 0 = 0

derivada 1 = 0

derivada n | esPrimo n = 1

| otherwise = (derivada a) * b + a * (derivada b)

where a = menorFactor n

b = n `div` a

-- (esPrimo n) se verifica si n es primo. Por ejemplo,

-- esPrimo 5 == True

-- esPrimo 6 == False

esPrimo :: Integer -> Bool

esPrimo 0 = False

esPrimo 1 = False

esPrimo n = n == menorFactor n

-- (menorFactor n) es el menor divisor primo de n (con n >= 2). Por

-- ejemplo,

-- menorFactor 6 == 2

-- menorFactor 7 == 7

-- menorFactor 15 == 3

menorFactor :: Integer -> Integer

menorFactor n

| even n = 2

| otherwise = head [x | x <- [3,5..]

, n `mod` x == 0]

-- 2ª solución

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

derivada2 :: Integer -> Integer

derivada2 0 = 0

derivada2 1 = 0

derivada2 n | isPrime n = 1

| otherwise = (derivada2 a) * b + a * (derivada2 b)

where (a:_) = primeFactors n

b = n `div` a

-- Comparación de eficiencia

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

-- λ> maximum [derivada n | n <- [1..10000]]

-- 53248

-- (1.59 secs, 1,091,452,552 bytes)

-- λ> maximum [derivada2 n | n <- [1..10000]]

-- 53248

-- (0.17 secs, 457,819,120 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_derivada :: Integer -> Property

prop_derivada x =

x > 0 ==>

derivada x == sum [(x * e) `div` p | (p,e) <- factorizacion x]

-- (factorizacion x) es la lista de las bases y exponentes de

-- la descomposición prima de x. Por ejemplo,

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

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

factorizacion n =

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

-- Su comprobación es

-- λ> quickCheck prop_derivada

-- +++ OK, passed 100 tests.

Referencias

Arithmetic derivative en Wikipedia.
Sucesión A003415 de OEIS.

Pensamiento

En ese jardín, Guiomar,
el mutuo jardín que inventan
dos corazones al par,
se funden y complementan
nuestras horas.

Antonio Machado

Postulado de Bertrand

PorJosé A. Alonso 30-diciembre-20196-enero-2020

El postulado de Bertrand afirma que para cualquier número entero n > 1, existe al menos un número primo p con n < p < 2n.

Definir la función

   siguientePrimo :: Integer -> Integer

1	siguientePrimo :: Integer -> Integer

tal que (siguientePrimo n) es el menor primo mayor que n. Por ejemplo,

   siguientePrimo 8   ==  11
   siguientePrimo 11  ==  13

1 2	siguientePrimo 8 == 11 siguientePrimo 11 == 13

Comprobar con QuickCheck el postulado de Bertrand; es decir, para todo entero n > 1, se verifica que n < p < 2n, donde p es (siguientePrimo n).

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck

siguientePrimo :: Integer -> Integer
siguientePrimo n = head (dropWhile (<= n) primes)

-- La propiedad es
postuladoDeBertrand :: Integer -> Property
postuladoDeBertrand n =
  n > 1 ==> n < p && p < 2 * n
  where p = siguientePrimo n

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

import Data.Numbers.Primes

import Test.QuickCheck

siguientePrimo :: Integer -> Integer

siguientePrimo n = head (dropWhile (<= n) primes)

-- La propiedad es

postuladoDeBertrand :: Integer -> Property

postuladoDeBertrand n =

n > 1 ==> n < p && p < 2 * n

where p = siguientePrimo n

-- La comprobación es

-- λ> quickCheck postuladoDeBertrand

-- +++ OK, passed 100 tests.

Referencias

Bertrand’s postulate en Wikipedia.

Pensamiento

Pero caer de cabeza,
en esta noche sin luna,
en medio de esta maleza,
junto a la negra laguna.

Antonio Machado

Posiciones de conjuntos finitos de naturales

PorJosé A. Alonso 27-diciembre-20193-enero-2020

En un ejercicio anterior se mostró que 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.

Además, se definió 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]]

Definir la función

   posicion :: [Integer] -> Integer

1	posicion :: [Integer] -> Integer

tal que (posicion xs) es la posición del conjunto finito de números naturales xs, representado por una lista decreciente, en enumeracionCFN. Por ejemplo,

   posicion [2,0]          ==  5
   posicion [2,1]          ==  6
   posicion [2,1,0]        ==  7
   posicion [0]            ==  1
   posicion [1,0]          ==  3
   posicion [2,1,0]        ==  7
   posicion [3,2,1,0]      ==  15
   posicion [4,3,2,1,0]    ==  31
   posicion [5,4,3,2,1,0]  ==  63
   length (show (posicion [3*10^7])) == 9030900

posicion [2,0] == 5

posicion [2,1] == 6

posicion [2,1,0] == 7

posicion [0] == 1

posicion [1,0] == 3

posicion [2,1,0] == 7

posicion [3,2,1,0] == 15

posicion [4,3,2,1,0] == 31

posicion [5,4,3,2,1,0] == 63

length (show (posicion [3*10^7])) == 9030900

Comprobar con QuickCheck que para todo número natural n,

   posicion [n,n-1..0] == 2^(n+1) - 1.

1	posicion [n,n-1..0] == 2^(n+1) - 1.

Soluciones

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

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

posicion :: [Integer] -> Integer
posicion xs =
  genericLength (takeWhile (< xs) enumeracionCFN)

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

posicion2 :: [Integer] -> Integer
posicion2 []     = 0
posicion2 (x:xs) = 2^x + posicion2 xs

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

posicion3 :: [Integer] -> Integer
posicion3 = foldr (\x -> (+) (2^x)) 0 

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

posicion4 :: [Integer] -> Integer
posicion4 xs = sum [2^x | x <- xs ]

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

-- La propiedad es
prop_equiv :: [Integer] -> Bool
prop_equiv xs =
  all (== posicion xs') [f xs' | f <- [ posicion2
                                      , posicion3
                                      , posicion4]]
  where xs' = reverse (sort (nub (map abs xs)))

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

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

-- La comparación es
--    λ> posicion [19,18,17,16,14,9,6]
--    1000000
--    (2.61 secs, 1,754,265,776 bytes)
--    λ> posicion2 [19,18,17,16,14,9,6]
--    1000000
--    (0.01 secs, 111,808 bytes)
--    λ> posicion3 [19,18,17,16,14,9,6]
--    1000000
--    λ> posicion4 [19,18,17,16,14,9,6]
--    1000000
--    (0.01 secs, 111,704 bytes)

--    λ> length (show (posicion2 [3*10^7]))
--    9030900
--    (2.06 secs, 571,911,304 bytes)
--    λ> length (show (posicion3 [3*10^7]))
--    9030900
--    (2.06 secs, 571,911,544 bytes)
--    λ> length (show (posicion4 [3*10^7]))
--    9030900
--    (2.05 secs, 571,911,464 bytes)

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

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

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

import Data.List (genericLength, nub, sort)

import Test.QuickCheck

-- 1ª solución

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

posicion :: [Integer] -> Integer

posicion xs =

genericLength (takeWhile (< xs) enumeracionCFN)

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

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

posicion2 :: [Integer] -> Integer

posicion2 [] = 0

posicion2 (x:xs) = 2^x + posicion2 xs

-- 3ª solución

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

posicion3 :: [Integer] -> Integer

posicion3 = foldr (\x -> (+) (2^x)) 0

-- 4ª solución

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

posicion4 :: [Integer] -> Integer

posicion4 xs = sum [2^x | x <- xs ]

-- Equivalencia de las definiciones

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

-- La propiedad es

prop_equiv :: [Integer] -> Bool

prop_equiv xs =

all (== posicion xs') [f xs' | f <- [ posicion2

, posicion3

, posicion4]]

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

-- La comprobación es

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

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> posicion [19,18,17,16,14,9,6]

-- 1000000

-- (2.61 secs, 1,754,265,776 bytes)

-- λ> posicion2 [19,18,17,16,14,9,6]

-- 1000000

-- (0.01 secs, 111,808 bytes)

-- λ> posicion3 [19,18,17,16,14,9,6]

-- 1000000

-- λ> posicion4 [19,18,17,16,14,9,6]

-- 1000000

-- (0.01 secs, 111,704 bytes)

-- λ> length (show (posicion2 [3*10^7]))

-- 9030900

-- (2.06 secs, 571,911,304 bytes)

-- λ> length (show (posicion3 [3*10^7]))

-- 9030900

-- (2.06 secs, 571,911,544 bytes)

-- λ> length (show (posicion4 [3*10^7]))

-- 9030900

-- (2.05 secs, 571,911,464 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_posicion :: Integer -> Property

prop_posicion n =

n >= 0 ==> posicion3 [n,n-1..0] == 2^(n+1) - 1

-- La comprobación es

-- λ> quickCheck prop_posicion

-- +++ OK, passed 100 tests.

Pensamiento

¡Volar sin alas donde todo es cielo!

Antonio Machado

Conjuntos con más sumas que restas

PorJosé A. Alonso 26-diciembre-20192-enero-2020

Dado un conjunto de números naturales, por ejemplo A = {0, 2, 3, 4}, calculamos las sumas de todos los pares de elementos de A. Como A tiene 4 elementos hay 16 pares, pero no todas sus sumas son distintas. En este caso solo hay 8 sumas distintas: {0, 2, 3, 4, 5, 6, 7, 8}. Procediendo análogamente hay 9 diferencias distinatas entre los pares de A: {-4, -3, -2, -1, 0, 1, 2, 3, 4}.

Experimentando con más conjuntos, se puede conjeturar que el número de restas es mayor que el de sumas y argumentar que que mientras que con dos números distintos sólo se produce una suma distints sin embargo se producen dos restas distintas. Por ejemplo, con 5 y 7 sólo se produce una suma (ya que 5+7 y 7+5 ambos dan 12) pero dos restas (ya que 5-7 y 7-5 dan -2 y 2, respectivamente).

Sin embargo, la conjetura es falsa. Un contraejemplo en el conjunto {0, 2, 3, 4, 7, 11, 12, 14}, que tiene 26 sumas distintas con sus pares de elementos pero sólo 25 restas.

Los conjuntos con más sumas distintas con sus pares de elementos que restas se llaman conjuntos MSQR (por «más sumas que restas»).

El objetivo de este ejercicio es calcular los conjuntos MSQR.

Definir las funciones

   tieneMSQR :: [Integer] -> Bool
   conjuntosMSQR :: [[Integer]]

1 2	tieneMSQR :: [Integer] -> Bool conjuntosMSQR :: [[Integer]]

tales que

(tieneMSQR xs) se verifica si el conjunto xs tiene más sumas que restas. Por ejemplo,

     tieneMSQR [0, 2, 3, 4]                 ==  False
     tieneMSQR [0, 2, 3, 4, 7, 11, 12, 14]  ==  True

1 2	tieneMSQR [0, 2, 3, 4] == False tieneMSQR [0, 2, 3, 4, 7, 11, 12, 14] == True

conjuntosMSQR es la lista de los conjuntos MSQR. Por ejemplo,

     λ> take 5 conjuntosMSQR
     [[14,12,11,7,4,3,2,0],
      [14,12,11,10,7,3,2,0],
      [14,13,12,9,5,4,2,1,0],
      [14,13,12,10,9,5,2,1,0],
      [15,13,12,8,5,4,3,1]]
   
      length (takeWhile (< [14]) conjuntosMSQR)   ==  0
      length (takeWhile (< [15]) conjuntosMSQR)   ==  4
      length (takeWhile (< [16]) conjuntosMSQR)   ==  10
      length (takeWhile (< [17]) conjuntosMSQR)   ==  30

λ> take 5 conjuntosMSQR

[[14,12,11,7,4,3,2,0],

[14,12,11,10,7,3,2,0],

[14,13,12,9,5,4,2,1,0],

[14,13,12,10,9,5,2,1,0],

[15,13,12,8,5,4,3,1]]

length (takeWhile (< [14]) conjuntosMSQR) == 0

length (takeWhile (< [15]) conjuntosMSQR) == 4

length (takeWhile (< [16]) conjuntosMSQR) == 10

length (takeWhile (< [17]) conjuntosMSQR) == 30

Soluciones

import Data.List (tails, nub, sort)

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

-- (sumas xs) es el conjunto de las sumas de pares de elementos de
-- xs. Por ejemplo,
--    sumas2 [0,2,3,4]  ==  [0,2,3,4,5,6,7,8]
sumas :: [Integer] -> [Integer]
sumas xs = nub [x + y | x <- xs, y <- xs]

-- (restas xs) es el conjunto de las restas de pares de elementos de
-- xs. Por ejemplo,
--    restas [0,2,3,4]  ==  [0,-2,-3,-4,2,-1,3,1,4]
restas :: [Integer] -> [Integer]
restas xs = nub [x - y | x <- xs, y <- xs]

tieneMSQR :: [Integer] -> Bool
tieneMSQR xs = length (sumas xs) > length (restas xs)

conjuntosMSQR :: [[Integer]]
conjuntosMSQR = [xs | xs <- enumeracionCFN, tieneMSQR xs]

-- enumeracionCFN es la enumeración de los conjuntos finitos de números
-- naturales del ejercicio anterior.
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
-- ===========

-- (sumas2 xs) es el conjunto de las sumas de pares de elementos de
-- xs. Por ejemplo,
--    sumas2 [0,2,3,4]  ==  [0,2,3,4,5,6,7,8]
--    sumas2 [0,2,3,4]  ==  [0,2,3,4,5,6,7,8]
sumas2 :: [Integer] -> [Integer]
sumas2 xs = nub [x + y | (x:ys) <- tails xs, y <- (x:ys)]

-- (restas2 xs) es el conjunto de las restas de pares de elementos de
-- xs. Por ejemplo,
--    sumas2 [0,2,3,4]  ==  [0,2,3,4,5,6,7,8]
--    restas2 [0,2,3,4]  ==  [0,-2,-3,-4,2,-1,3,1,4]
restas2 :: [Integer] -> [Integer]
restas2 xs = 0 : ys ++ map negate ys
  where ys = nub [x - y | (x:ys) <- tails (sort xs), y <- ys]

tieneMSQR2 :: [Integer] -> Bool
tieneMSQR2 xs = length (sumas2 xs) > length (restas2 xs)

conjuntosMSQR2 :: [[Integer]]
conjuntosMSQR2 = [xs | xs <- enumeracionCFN, tieneMSQR2 xs]

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

--    λ> length (takeWhile (< [17,16..0]) conjuntosMSQR)
--    66
--    (21.36 secs, 10,301,222,168 bytes)
--    λ> length (takeWhile (< [17,16..0]) conjuntosMSQR2)
--    66
--    (10.13 secs, 7,088,969,752 bytes)

import Data.List (tails, nub, sort)

-- 1ª solución

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

-- (sumas xs) es el conjunto de las sumas de pares de elementos de

-- xs. Por ejemplo,

-- sumas2 [0,2,3,4] == [0,2,3,4,5,6,7,8]

sumas :: [Integer] -> [Integer]

sumas xs = nub [x + y | x <- xs, y <- xs]

-- (restas xs) es el conjunto de las restas de pares de elementos de

-- xs. Por ejemplo,

-- restas [0,2,3,4] == [0,-2,-3,-4,2,-1,3,1,4]

restas :: [Integer] -> [Integer]

restas xs = nub [x - y | x <- xs, y <- xs]

tieneMSQR :: [Integer] -> Bool

tieneMSQR xs = length (sumas xs) > length (restas xs)

conjuntosMSQR :: [[Integer]]

conjuntosMSQR = [xs | xs <- enumeracionCFN, tieneMSQR xs]

-- enumeracionCFN es la enumeración de los conjuntos finitos de números

-- naturales del ejercicio anterior.

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

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

-- (sumas2 xs) es el conjunto de las sumas de pares de elementos de

-- xs. Por ejemplo,

-- sumas2 [0,2,3,4] == [0,2,3,4,5,6,7,8]

sumas2 :: [Integer] -> [Integer]

sumas2 xs = nub [x + y | (x:ys) <- tails xs, y <- (x:ys)]

-- (restas2 xs) es el conjunto de las restas de pares de elementos de

-- xs. Por ejemplo,

-- sumas2 [0,2,3,4] == [0,2,3,4,5,6,7,8]

-- restas2 [0,2,3,4] == [0,-2,-3,-4,2,-1,3,1,4]

restas2 :: [Integer] -> [Integer]

restas2 xs = 0 : ys ++ map negate ys

where ys = nub [x - y | (x:ys) <- tails (sort xs), y <- ys]

tieneMSQR2 :: [Integer] -> Bool

tieneMSQR2 xs = length (sumas2 xs) > length (restas2 xs)

conjuntosMSQR2 :: [[Integer]]

conjuntosMSQR2 = [xs | xs <- enumeracionCFN, tieneMSQR2 xs]

-- Comparación de eficiencia

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

-- λ> length (takeWhile (< [17,16..0]) conjuntosMSQR)

-- 66

-- (21.36 secs, 10,301,222,168 bytes)

-- λ> length (takeWhile (< [17,16..0]) conjuntosMSQR2)

-- 66

-- (10.13 secs, 7,088,969,752 bytes)

Pensamiento

¡Qué fácil es volar, qué fácil es!
Todo consiste en no dejar que el suelo
se acerque a nuestros pies.

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

Infinitud de primos gemelos

PorJosé A. Alonso 24-diciembre-201931-diciembre-2019

Un par de números primos (p,q) es un par de números primos gemelos si su distancia de 2; es decir, si q = p+2. Por ejemplo, (17,19) es una par de números primos gemelos.

La conjetura de los primos gemelos postula la existencia de infinitos pares de primos gemelos.

Definir la constante

   primosGemelos :: [(Integer,Integer)]

1	primosGemelos :: [(Integer,Integer)]

tal que sus elementos son los pares de primos gemelos. Por ejemplo,

   λ> take 7 primosGemelos
   [(3,5),(5,7),(11,13),(17,19),(29,31),(41,43),(59,61)]
   λ> primosGemelos !! (4*10^4)
   (6381911,6381913)

λ> take 7 primosGemelos

[(3,5),(5,7),(11,13),(17,19),(29,31),(41,43),(59,61)]

λ> primosGemelos !! (4*10^4)

(6381911,6381913)

Comprobar con QuickCheck la conjetura de los primos gemelos.

Soluciones

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

-- 1ª solución
primosGemelos :: [(Integer,Integer)]
primosGemelos = [(x,x+2) | x <- primes, isPrime (x+2)]

-- 2ª solución
primosGemelos2 :: [(Integer,Integer)]
primosGemelos2 = [(x,y) | (x,y) <- zip primes (tail primes)
                        , y - x == 2]

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

--    λ> primosGemelos !! (2*10^4)
--    (2840447,2840449)
--    (3.93 secs, 12,230,474,952 bytes)
--    λ> primosGemelos2 !! (2*10^4)
--    (2840447,2840449)
--    (0.77 secs, 2,202,822,456 bytes)

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

-- La propiedad es
prop_primosGemelos :: Integer -> Property
prop_primosGemelos n =
  n >= 0 ==> not (null [(x,y) | (x,y) <- primosGemelos2, x > n])

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

import Data.Numbers.Primes (primes, isPrime)

import Test.QuickCheck

-- 1ª solución

primosGemelos :: [(Integer,Integer)]

primosGemelos = [(x,x+2) | x <- primes, isPrime (x+2)]

-- 2ª solución

primosGemelos2 :: [(Integer,Integer)]

primosGemelos2 = [(x,y) | (x,y) <- zip primes (tail primes)

, y - x == 2]

-- Comparación de eficiencia

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

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

-- (2840447,2840449)

-- (3.93 secs, 12,230,474,952 bytes)

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

-- (2840447,2840449)

-- (0.77 secs, 2,202,822,456 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_primosGemelos :: Integer -> Property

prop_primosGemelos n =

n >= 0 ==> not (null [(x,y) | (x,y) <- primosGemelos2, x > n])

-- La comprobación es

-- λ> quickCheck prop_primosGemelos

-- +++ OK, passed 100 tests.

Pensamiento

El sentimiento ha de tener tanto de individual como de genérico; debe orientarse hacia valores universales, o que pretenden serlo.

Antonio Machado

Suma de números de Fibonacci con índice impar

PorJosé A. Alonso 23-diciembre-201930-diciembre-2019

La sucesión de Fibonacci, F(n), es la siguiente sucesión infinita de números naturales:

   0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, ...

1	0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, ...

La sucesión comienza con los números 0 y 1. A partir de estos, cada término es la suma de los dos anteriores.

Definir la función

   sumaFibsIndiceImpar :: Int -> Integer

1	sumaFibsIndiceImpar :: Int -> Integer

tal que (sumaFibsIndiceImpar n) es la suma de los n primeros términos de la sucesión de Fibonacci no índice impar; es decir,

   sumaFibsIndiceImpar n = F(1) + F(3) + ... + F(2*n-1)

1	sumaFibsIndiceImpar n = F(1) + F(3) + ... + F(2*n-1)

Por ejemplo,

   sumaFibsIndiceImpar 1  ==  1
   sumaFibsIndiceImpar 2  ==  3
   sumaFibsIndiceImpar 3  ==  8
   sumaFibsIndiceImpar 4  ==  21
   sumaFibsIndiceImpar 5  ==  55
   sumaFibsIndiceImpar (10^4) `rem` (10^9)  ==  213093125

sumaFibsIndiceImpar 1 == 1

sumaFibsIndiceImpar 2 == 3

sumaFibsIndiceImpar 3 == 8

sumaFibsIndiceImpar 4 == 21

sumaFibsIndiceImpar 5 == 55

sumaFibsIndiceImpar (10^4) `rem` (10^9) == 213093125

En los ejemplos anteriores se observa que

   sumaFibsIndiceImpar 1  ==  F(2)
   sumaFibsIndiceImpar 2  ==  F(4)
   sumaFibsIndiceImpar 3  ==  F(6)
   sumaFibsIndiceImpar 4  ==  F(8)
   sumaFibsIndiceImpar 5  ==  F(10)

sumaFibsIndiceImpar 1 == F(2)

sumaFibsIndiceImpar 2 == F(4)

sumaFibsIndiceImpar 3 == F(6)

sumaFibsIndiceImpar 4 == F(8)

sumaFibsIndiceImpar 5 == F(10)

Comprobar con QuickCheck que (sumaFibsIndiceImpar n) es F(2n); es decir, el 2n-ésimo número de Fibonacci

Soluciones

import Test.QuickCheck

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

sumaFibsIndiceImpar :: Int -> Integer
sumaFibsIndiceImpar n = sum [fib (2*k-1) | k <- [1..n]]

-- (fib n) es el n-ésimo término de la sucesión de Fibonacci. Por
-- ejemplo,
--    fib 6  ==  8
fib :: Int -> Integer
fib n = fibs !! n

-- fibs es la lista de términos de la sucesión de Fibonacci. Por ejemplo,
--    λ> take 20 fibs
--    [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]
fibs :: [Integer]
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

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

sumaFibsIndiceImpar2 :: Int -> Integer
sumaFibsIndiceImpar2 n =
  sum [a | (a,b) <- zip fibs [0..2*n], odd b]

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

--    λ> sumaFibsIndiceImpar (10^4) `rem` (10^9)
--    213093125
--    (0.98 secs, 13,889,312 bytes)
--    λ> sumaFibsIndiceImpar2 (10^4) `rem` (10^9)
--    213093125
--    (0.05 secs, 18,047,720 bytes)

-- Comprobación
-- ============

-- La propiedad es
prop_sumaFibsIndiceImpar :: Int -> Property
prop_sumaFibsIndiceImpar n =
  n >= 0 ==> sumaFibsIndiceImpar n == fib (2*n)

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

import Test.QuickCheck

-- 1ª solución

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

sumaFibsIndiceImpar :: Int -> Integer

sumaFibsIndiceImpar n = sum [fib (2*k-1) | k <- [1..n]]

-- (fib n) es el n-ésimo término de la sucesión de Fibonacci. Por

-- ejemplo,

-- fib 6 == 8

fib :: Int -> Integer

fib n = fibs !! n

-- fibs es la lista de términos de la sucesión de Fibonacci. Por ejemplo,

-- λ> take 20 fibs

-- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]

fibs :: [Integer]

fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

-- 2ª solución

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

sumaFibsIndiceImpar2 :: Int -> Integer

sumaFibsIndiceImpar2 n =

sum [a | (a,b) <- zip fibs [0..2*n], odd b]

-- Comparación de eficiencia

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

-- λ> sumaFibsIndiceImpar (10^4) `rem` (10^9)

-- 213093125

-- (0.98 secs, 13,889,312 bytes)

-- λ> sumaFibsIndiceImpar2 (10^4) `rem` (10^9)

-- 213093125

-- (0.05 secs, 18,047,720 bytes)

-- Comprobación

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

-- La propiedad es

prop_sumaFibsIndiceImpar :: Int -> Property

prop_sumaFibsIndiceImpar n =

n >= 0 ==> sumaFibsIndiceImpar n == fib (2*n)

-- La comprobación es

-- λ> quickCheck prop_sumaFibsIndiceImpar

-- +++ OK, passed 100 tests.

Referencia

Sum of sequence of odd index Fibonacci numbers en ProofWiki.

Pensamiento

El corazón del poeta, tan rico en sonoridades, es casi un insulto a la afonía cordial de la masa.

Antonio Machado

El teorema de Navidad de Fermat

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

El 25 de diciembre de 1640, en una carta a Mersenne, Fermat demostró la conjetura de Girard: todo primo de la forma 4n+1 puede expresarse de manera única como suma de dos cuadrados. Por eso es conocido como el Teorema de Navidad de Fermat.

Definir las funciones

   representaciones :: Integer -> [(Integer,Integer)]
   primosImparesConRepresentacionUnica :: [Integer]
   primos4nM1 :: [Integer]

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

primosImparesConRepresentacionUnica :: [Integer]

primos4nM1 :: [Integer]

tales que

(representaciones n) es la lista de pares de números naturales (x,y) tales que n = x^2 + y^2 con x <= y. Por ejemplo,

     representaciones  20           ==  [(2,4)]
     representaciones  25           ==  [(0,5),(3,4)]
     representaciones 325           ==  [(1,18),(6,17),(10,15)]
     representaciones 100000147984  ==  [(0,316228)]
     length (representaciones (10^10))    ==  6
     length (representaciones (4*10^12))  ==  7

representaciones 20 == [(2,4)]

representaciones 25 == [(0,5),(3,4)]

representaciones 325 == [(1,18),(6,17),(10,15)]

representaciones 100000147984 == [(0,316228)]

length (representaciones (10^10)) == 6

length (representaciones (4*10^12)) == 7

primosImparesConRepresentacionUnica es la lista de los números primos impares que se pueden escribir exactamente de una manera como suma de cuadrados de pares de números naturales (x,y) con x <= y. Por ejemplo,

     λ> take 20 primosImparesConRepresentacionUnica
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

1 2	λ> take 20 primosImparesConRepresentacionUnica [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

primos4nM1 es la lista de los números primos que se pueden escribir como uno más un múltiplo de 4 (es decir, que son congruentes con 1 módulo 4). Por ejemplo,

     λ> take 20 primos4nM1
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

1 2	λ> take 20 primos4nM1 [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

El teorema de Navidad de Fermat afirma que un número primo impar p se puede escribir exactamente de una manera como suma de dos cuadrados de números naturales p = x² + y^2 (con x <= y) si, y sólo si, p se puede escribir como uno más un múltiplo de 4 (es decir, que es congruente con 1 módulo 4).

Comprobar con QuickCheck el teorema de Navidad de Fermat; es decir, que para todo número n, los n-ésimos elementos de primosImparesConRepresentacionUnica y de primos4nM1 son iguales.

Soluciones

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

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

representaciones :: Integer -> [(Integer,Integer)]
representaciones n =
  [(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]

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

representaciones2 :: Integer -> [(Integer,Integer)]
representaciones2 n =
  [(x,raiz z) | x <- [0..raiz (n `div` 2)] 
              , let z = n - x*x
              , esCuadrado z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer 
raiz 0 = 0
raiz 1 = 1
raiz x = aux (0,x)
    where aux (a,b) | d == x    = c
                    | c == a    = a
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c) 
              where c = (a+b) `div` 2
                    d = c^2

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

representaciones3 :: Integer -> [(Integer,Integer)]
representaciones3 n =
  [(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)] 
               , let z = n - x*x
               , esCuadrado3 z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado3 25  ==  True
--    esCuadrado3 26  ==  False
esCuadrado3 :: Integer -> Bool
esCuadrado3 x = x == y * y
  where y = raiz3 x

-- (raiz3 x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz3 25  ==  5
--    raiz3 24  ==  4
--    raiz3 26  ==  5
raiz3 :: Integer -> Integer
raiz3 x = floor (sqrt (fromIntegral x))
          
-- 4ª definición de representaciones
-- =================================
 
representaciones4 :: Integer -> [(Integer, Integer)]
representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))
  where aux x y
          | x > y     = [] 
          | otherwise = case compare (x*x + y*y) n of
                          LT -> aux (x + 1) y
                          EQ -> (x, y) : aux (x + 1) (y - 1)
                          GT -> aux x (y - 1)

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

-- La propiedad es
prop_representaciones_equiv :: (Positive Integer) -> Bool
prop_representaciones_equiv (Positive n) =
  representaciones  n == representaciones2 n &&
  representaciones2 n == representaciones3 n &&
  representaciones3 n == representaciones4 n

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

-- Comparación de eficiencia de las definiciones de representaciones
-- =================================================================

--    λ> representaciones 3025
--    [(0,55),(33,44)]
--    (2.86 secs, 1,393,133,528 bytes)
--    λ> representaciones2 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 867,944 bytes)
--    λ> representaciones3 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 173,512 bytes)
--    λ> representaciones4 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 423,424 bytes)
--    
--    λ> length (representaciones2 (10^10))
--    6
--    (3.38 secs, 2,188,903,544 bytes)
--    λ> length (representaciones3 (10^10))
--    6
--    (0.10 secs, 62,349,048 bytes)
--    λ> length (representaciones4 (10^10))
--    6
--    (0.11 secs, 48,052,360 bytes)
--
--    λ> length (representaciones3 (4*10^12))
--    7
--    (1.85 secs, 1,222,007,176 bytes)
--    λ> length (representaciones4 (4*10^12))
--    7
--    (1.79 secs, 953,497,480 bytes)

-- Definición de primosImparesConRepresentacionUnica
-- =================================================

primosImparesConRepresentacionUnica :: [Integer]
primosImparesConRepresentacionUnica =
  [x | x <- tail primes
     , length (representaciones4 x) == 1]

-- Definición de primos4nM1
-- ========================

primos4nM1 :: [Integer]
primos4nM1 = [x | x <- primes
                , x `mod` 4 == 1]

-- Teorema de Navidad de Fermat
-- ============================

-- La propiedad es
prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool
prop_teoremaDeNavidadDeFermat (Positive n) =
  primosImparesConRepresentacionUnica !! n == primos4nM1 !! n
             
-- La comprobación es
--    λ> quickCheck prop_teoremaDeNavidadDeFermat
--    +++ 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

import Data.Numbers.Primes (primes)

import Test.QuickCheck

-- 1ª definición de representaciones

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

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

representaciones n =

[(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]

-- 2ª definición de representaciones

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

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

representaciones2 n =

[(x,raiz z) | x <- [0..raiz (n `div` 2)]

, let z = n - x*x

, esCuadrado z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz 25 == 5

-- raiz 24 == 4

-- raiz 26 == 5

raiz :: Integer -> Integer

raiz 0 = 0

raiz 1 = 1

raiz x = aux (0,x)

where aux (a,b) | d == x = c

| c == a = a

| d < x = aux (c,b)

| otherwise = aux (a,c)

where c = (a+b) `div` 2

d = c^2

-- 3ª definición de representaciones

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

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

representaciones3 n =

[(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)]

, let z = n - x*x

, esCuadrado3 z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado3 25 == True

-- esCuadrado3 26 == False

esCuadrado3 :: Integer -> Bool

esCuadrado3 x = x == y * y

where y = raiz3 x

-- (raiz3 x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz3 25 == 5

-- raiz3 24 == 4

-- raiz3 26 == 5

raiz3 :: Integer -> Integer

raiz3 x = floor (sqrt (fromIntegral x))

-- 4ª definición de representaciones

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

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

representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))

where aux x y

| x > y = []

| otherwise = case compare (x*x + y*y) n of

LT -> aux (x + 1) y

EQ -> (x, y) : aux (x + 1) (y - 1)

GT -> aux x (y - 1)

-- Equivalencia de las definiciones de representaciones

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

-- La propiedad es

prop_representaciones_equiv :: (Positive Integer) -> Bool

prop_representaciones_equiv (Positive n) =

representaciones n == representaciones2 n &&

representaciones2 n == representaciones3 n &&

representaciones3 n == representaciones4 n

-- La comprobación es

-- λ> quickCheck prop_representaciones_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de las definiciones de representaciones

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

-- λ> representaciones 3025

-- [(0,55),(33,44)]

-- (2.86 secs, 1,393,133,528 bytes)

-- λ> representaciones2 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 867,944 bytes)

-- λ> representaciones3 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 173,512 bytes)

-- λ> representaciones4 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 423,424 bytes)

-- λ> length (representaciones2 (10^10))

-- 6

-- (3.38 secs, 2,188,903,544 bytes)

-- λ> length (representaciones3 (10^10))

-- 6

-- (0.10 secs, 62,349,048 bytes)

-- λ> length (representaciones4 (10^10))

-- 6

-- (0.11 secs, 48,052,360 bytes)

-- λ> length (representaciones3 (4*10^12))

-- 7

-- (1.85 secs, 1,222,007,176 bytes)

-- λ> length (representaciones4 (4*10^12))

-- 7

-- (1.79 secs, 953,497,480 bytes)

-- Definición de primosImparesConRepresentacionUnica

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

primosImparesConRepresentacionUnica :: [Integer]

primosImparesConRepresentacionUnica =

[x | x <- tail primes

, length (representaciones4 x) == 1]

-- Definición de primos4nM1

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

primos4nM1 :: [Integer]

primos4nM1 = [x | x <- primes

, x `mod` 4 == 1]

-- Teorema de Navidad de Fermat

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

-- La propiedad es

prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool

prop_teoremaDeNavidadDeFermat (Positive n) =

primosImparesConRepresentacionUnica !! n == primos4nM1 !! n

-- La comprobación es

-- λ> quickCheck prop_teoremaDeNavidadDeFermat

-- +++ OK, passed 100 tests.

Pensamiento

Dijo Dios: brote la nada
Y alzó su mano derecha,
hasta ocultar su mirada.
Y quedó la nada hecha.

Antonio Machado

Árbol binario de divisores

PorJosé A. Alonso 18-diciembre-201925-diciembre-2019

El árbol binario de los divisores de 24 es

2 45

3 15

3 5

Se puede representar por

   N 90 (H 2) (N 45 (H 3) (N 15 (H 3) (H 5)))

1	N 90 (H 2) (N 45 (H 3) (N 15 (H 3) (H 5)))

usando el tipo de dato definido por

   data Arbol = H Int
              | N Int Arbol Arbol
     deriving (Eq, Show)

data Arbol = H Int

| N Int Arbol Arbol

deriving (Eq, Show)

Análogamente se obtiene el árbol binario de cualquier número x: se comienza en x y en cada paso se tiene dos hijos (su menor divisor y su cociente) hasta obtener números primos en las hojas.

Definir las funciones

   arbolDivisores      :: Int -> Arbol
   hojasArbolDivisores :: Int -> [Int]

1 2	arbolDivisores :: Int -> Arbol hojasArbolDivisores :: Int -> [Int]

tales que

(arbolDivisores x) es el árbol binario de los divisores de x. Por ejemplo,

     λ> arbolDivisores 90
     N 90 (H 2) (N 45 (H 3) (N 15 (H 3) (H 5)))
     λ> arbolDivisores 24
     N 24 (H 2) (N 12 (H 2) (N 6 (H 2) (H 3)))
     λ> arbolDivisores 300
     N 300 (H 2) (N 150 (H 2) (N 75 (H 3) (N 25 (H 5) (H 5))))

λ> arbolDivisores 90

N 90 (H 2) (N 45 (H 3) (N 15 (H 3) (H 5)))

λ> arbolDivisores 24

N 24 (H 2) (N 12 (H 2) (N 6 (H 2) (H 3)))

λ> arbolDivisores 300

N 300 (H 2) (N 150 (H 2) (N 75 (H 3) (N 25 (H 5) (H 5))))

(hojasArbolDivisores x) es la lista de las hohas del árbol binario de los divisores de x. Por ejemplo

     hojasArbolDivisores 90   ==  [2,3,3,5]
     hojasArbolDivisores 24   ==  [2,2,2,3]
     hojasArbolDivisores 300  ==  [2,2,3,5,5]

hojasArbolDivisores 90 == [2,3,3,5]

hojasArbolDivisores 24 == [2,2,2,3]

hojasArbolDivisores 300 == [2,2,3,5,5]

Soluciones

import Data.Numbers.Primes (isPrime, primeFactors)

data Arbol = H Int
           | N Int Arbol Arbol
  deriving (Eq, Show)

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

arbolDivisores :: Int -> Arbol
arbolDivisores x
  | y == x    = H x
  | otherwise = N x (H y) (arbolDivisores (x `div` y))
  where y = menorDivisor x

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

hojasArbolDivisores :: Int -> [Int]
hojasArbolDivisores = hojas . arbolDivisores

-- (hojas a) es la lista de las hojas del árbol a. Por ejemplo,
--    hojas (N 3 (H 4) (N 5 (H 7) (H 2)))  ==  [4,7,2]
hojas :: Arbol -> [Int]
hojas (H x)     = [x]
hojas (N _ i d) = hojas i ++ hojas d

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

arbolDivisores2 :: Int -> Arbol
arbolDivisores2 x
  | y == x    = H x
  | otherwise = N x (H y) (arbolDivisores (x `div` y))
  where (y:_) = primeFactors x

hojasArbolDivisores2 :: Int -> [Int]
hojasArbolDivisores2 = primeFactors

import Data.Numbers.Primes (isPrime, primeFactors)

data Arbol = H Int

| N Int Arbol Arbol

deriving (Eq, Show)

-- 1ª solución

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

arbolDivisores :: Int -> Arbol

arbolDivisores x

| y == x = H x

| otherwise = N x (H y) (arbolDivisores (x `div` y))

where y = menorDivisor x

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

-- menorDivisor 45 == 3

-- menorDivisor 5 == 5

menorDivisor :: Int -> Int

menorDivisor x =

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

hojasArbolDivisores :: Int -> [Int]

hojasArbolDivisores = hojas . arbolDivisores

-- (hojas a) es la lista de las hojas del árbol a. Por ejemplo,

-- hojas (N 3 (H 4) (N 5 (H 7) (H 2))) == [4,7,2]

hojas :: Arbol -> [Int]

hojas (H x) = [x]

hojas (N _ i d) = hojas i ++ hojas d

-- 2ª solución

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

arbolDivisores2 :: Int -> Arbol

arbolDivisores2 x

| y == x = H x

| otherwise = N x (H y) (arbolDivisores (x `div` y))

where (y:_) = primeFactors x

hojasArbolDivisores2 :: Int -> [Int]

hojasArbolDivisores2 = primeFactors

Pensamiento

Cuando el Ser que se es hizo la nada
y reposó que bien lo merecía,
ya tuvo el día noche, y compañía
tuvo el amante en la ausencia de la amada.

Antonio Machado

Intersección de listas infinitas crecientes

PorJosé A. Alonso 17-diciembre-201924-diciembre-2019

Definir la función

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

1	interseccion :: Ord a => [[a]] -> [a]

tal que (interseccion xss) es la intersección de la lista no vacía de listas infinitas crecientes xss; es decir, la lista de los elementos que pertenecen a todas las listas de xss. Por ejemplo,

   λ> take 10 (interseccion [[2,4..],[3,6..],[5,10..]])
   [30,60,90,120,150,180,210,240,270,300]
   λ> take 10 (interseccion [[2,5..],[3,5..],[5,7..]])
   [5,11,17,23,29,35,41,47,53,59]

λ> take 10 (interseccion [[2,4..],[3,6..],[5,10..]])

[30,60,90,120,150,180,210,240,270,300]

λ> take 10 (interseccion [[2,5..],[3,5..],[5,7..]])

[5,11,17,23,29,35,41,47,53,59]

Soluciones

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

interseccion :: Ord a => [[a]] -> [a]
interseccion [xs]        = xs
interseccion (xs:ys:zss) = interseccionDos xs (interseccion (ys:zss))

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

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

interseccion2 :: Ord a => [[a]] -> [a]
interseccion2 = foldl1 interseccionDos

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

interseccion3 :: Ord a => [[a]] -> [a]
interseccion3 (xs:xss) =
  [x | x <- xs, all (x `pertenece`) xss]
 
pertenece :: Ord a => a -> [a] -> Bool
pertenece x xs = x == head (dropWhile (<x) xs)

-- 1ª solución

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

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

interseccion [xs] = xs

interseccion (xs:ys:zss) = interseccionDos xs (interseccion (ys:zss))

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

interseccionDos (x:xs) (y:ys)

| x == y = x : interseccionDos xs ys

| x < y = interseccionDos (dropWhile (<y) xs) (y:ys)

| otherwise = interseccionDos (x:xs) (dropWhile (<x) ys)

-- 2ª solución

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

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

interseccion2 = foldl1 interseccionDos

-- 3ª solución

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

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

interseccion3 (xs:xss) =

[x | x <- xs, all (x `pertenece`) xss]

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

pertenece x xs = x == head (dropWhile (<x) xs)

Pensamiento

Dios no es el creador del mundo (según Martín), sino el creador de la nada.

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

Evaluación de árboles de expresiones aritméticas

PorJosé A. Alonso 12-diciembre-201919-diciembre-2019

Las expresiones aritméticas se pueden representar como árboles con números en las hojas y operaciones en los nodos. Por ejemplo, la expresión «9-2*4» se puede representar por el árbol

/ \

9 *

/ \

2 4

Definiendo el tipo de dato Arbol por

   data Arbol = H Int | N (Int -> Int -> Int) Arbol Arbol

1	data Arbol = H Int \| N (Int -> Int -> Int) Arbol Arbol

la representación del árbol anterior es

   N (-) (H 9) (N (*) (H 2) (H 4))

1	N (-) (H 9) (N (*) (H 2) (H 4))

Definir la función

   valor :: Arbol -> Int

1	valor :: Arbol -> Int

tal que (valor a) es el valor de la expresión aritmética correspondiente al árbol a. Por ejemplo,

   valor (N (-) (H 9) (N (*) (H 2) (H 4)))    ==  1
   valor (N (+) (H 9) (N (*) (H 2) (H 4)))    ==  17
   valor (N (+) (H 9) (N (div) (H 4) (H 2)))  ==  11
   valor (N (+) (H 9) (N (max) (H 4) (H 2)))  ==  13

valor (N (-) (H 9) (N (*) (H 2) (H 4))) == 1

valor (N (+) (H 9) (N (*) (H 2) (H 4))) == 17

valor (N (+) (H 9) (N (div) (H 4) (H 2))) == 11

valor (N (+) (H 9) (N (max) (H 4) (H 2))) == 13

Soluciones

data Arbol = H Int | N (Int -> Int -> Int) Arbol Arbol

valor :: Arbol -> Int
valor (H x)     = x
valor (N f i d) = f (valor i) (valor d)

data Arbol = H Int | N (Int -> Int -> Int) Arbol Arbol

valor :: Arbol -> Int

valor (H x) = x

valor (N f i d) = f (valor i) (valor d)

Pensamiento

En cierto modo, las matemáticas no son el arte de responder preguntas matemáticas, es el arte de hacer las preguntas correctas, las preguntas que te dan una idea, las que te guían en direcciones interesantes, las que se conectan con muchas otras preguntas interesantes, las que tienen hermosas respuestas. ~ Gregory Chaitin

Suma de primos menores

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

La suma de los primos menores que 10 es 2 + 3 + 5 + 7 = 17.

Definir la función

   sumaPrimosMenores :: Integer -> Integer

1	sumaPrimosMenores :: Integer -> Integer

tal que (sumaPrimosMenores n) es la suma de los primos menores que n. Por ejemplo,

   sumaPrimosMenores 10        ==  17
   sumaPrimosMenores (5*10^5)  ==  9914236195

1 2	sumaPrimosMenores 10 == 17 sumaPrimosMenores (5*10^5) == 9914236195

Nota: Este ejercicio está basado en el problema 10 del Proyecto Euler

Soluciones

import Data.Numbers.Primes (primes)

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

sumaPrimosMenores :: Integer -> Integer
sumaPrimosMenores n = sum (takeWhile (<n) primos)

-- primos es la lista de los números primos. Por ejemplo,
--    λ> take 12 primos
--    [2,3,5,7,11,13,17,19,23,29,31,37]
primos :: [Integer]
primos = 2 : filter esPrimo [3,5..]

esPrimo :: Integer -> Bool
esPrimo n = null [x | x <- [2..(ceiling . sqrt . fromIntegral) n]
                    , n `mod` x == 0]

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

sumaPrimosMenores2 :: Integer -> Integer
sumaPrimosMenores2 n = sum (takeWhile (<n) primos2)

primos2 :: [Integer]
primos2 = 2 : filter esPrimo2 [3,5..]
 
esPrimo2 :: Integer -> Bool
esPrimo2 x =
  all ((/= 0) . mod x)
  (takeWhile (<= floor (sqrt (fromIntegral x))) primos2)

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

sumaPrimosMenores3 :: Integer -> Integer
sumaPrimosMenores3 n = sum (takeWhile (<n) primes)

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

--    λ> sumaPrimosMenores (2*10^5)
--    1709600813
--    (2.56 secs, 1,522,015,240 bytes)
--    λ> sumaPrimosMenores2 (2*10^5)
--    1709600813
--    (0.56 secs, 376,951,456 bytes)
--    λ> sumaPrimosMenores3 (2*10^5)
--    1709600813
--    (0.07 secs, 62,321,888 bytes)

import Data.Numbers.Primes (primes)

-- 1ª solución

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

sumaPrimosMenores :: Integer -> Integer

sumaPrimosMenores n = sum (takeWhile (<n) primos)

-- primos es la lista de los números primos. Por ejemplo,

-- λ> take 12 primos

-- [2,3,5,7,11,13,17,19,23,29,31,37]

primos :: [Integer]

primos = 2 : filter esPrimo [3,5..]

esPrimo :: Integer -> Bool

esPrimo n = null [x | x <- [2..(ceiling . sqrt . fromIntegral) n]

, n `mod` x == 0]

-- 2ª solución

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

sumaPrimosMenores2 :: Integer -> Integer

sumaPrimosMenores2 n = sum (takeWhile (<n) primos2)

primos2 :: [Integer]

primos2 = 2 : filter esPrimo2 [3,5..]

esPrimo2 :: Integer -> Bool

esPrimo2 x =

all ((/= 0) . mod x)

(takeWhile (<= floor (sqrt (fromIntegral x))) primos2)

-- 3ª solución

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

sumaPrimosMenores3 :: Integer -> Integer

sumaPrimosMenores3 n = sum (takeWhile (<n) primes)

-- Comparación de eficiencia

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

-- λ> sumaPrimosMenores (2*10^5)

-- 1709600813

-- (2.56 secs, 1,522,015,240 bytes)

-- λ> sumaPrimosMenores2 (2*10^5)

-- 1709600813

-- (0.56 secs, 376,951,456 bytes)

-- λ> sumaPrimosMenores3 (2*10^5)

-- 1709600813

-- (0.07 secs, 62,321,888 bytes)

Pensamiento

El movimiento no es nada esencial. La fuerza puede ser inmóvil (lo es en su estado de pureza); mas no por ello deja de ser activa.

Antonio Machado

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

Sublistas con producto dado

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

Definir las funciones

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

1 2	sublistasConProducto :: Integer -> [Integer] -> [[Integer]] unifactorizables :: [Integer]

tales que

(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]
     []
     λ> length (sublistasConProducto 1234567 [1..1234567])
     4

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

[]

λ> length (sublistasConProducto 1234567 [1..1234567])

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

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

Soluciones

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

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

sublistasConProducto :: Integer -> [Integer] -> [[Integer]]
sublistasConProducto n xs =
  [ys | ys <- subsequences xs
      , product ys == n]

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

sublistasConProducto2 :: Integer -> [Integer] -> [[Integer]]
sublistasConProducto2 _ [] = []
sublistasConProducto2 n (x:xs)
  | x > n     = []
  | x == n    = [[x]]
  | r == 0    = map (x:) (sublistasConProducto2 q xs)
                ++ sublistasConProducto2 n xs
  | otherwise = sublistasConProducto2 n xs
  where (q,r) = quotRem n x

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

-- La propiedad es
prop_sublistasConProducto :: Integer -> [Integer] -> Bool
prop_sublistasConProducto n xs =
  sort (sublistasConProducto n' xs') == sublistasConProducto2 n' xs'
  where n'  = 2 + abs n
        xs' = (nub . sort . map ((+2) . abs)) xs

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

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

--    λ> sublistasConProducto 15 [1..23]
--    [[3,5],[1,3,5],[15],[1,15]]
--    (3.44 secs, 7,885,411,472 bytes)
--    λ> sublistasConProducto2 15 [1..23]
--    [[1,3,5],[1,15],[3,5],[15]]
--    (0.01 secs, 135,056 bytes)
--
--    λ> length (sublistasConProducto2 1234567 [1..1234567])
--    4
--    (1.49 secs, 1,054,380,480 bytes)

-- Definición de unifactorizables
-- ==============================

unifactorizables :: [Integer]
unifactorizables =
  [n | n <- [2..]
     , length (sublistasConProducto2 n [2..n]) == 1]

import Test.QuickCheck

import Data.List (nub, sort, subsequences)

-- 1ª solución

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

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

sublistasConProducto n xs =

[ys | ys <- subsequences xs

, product ys == n]

-- 2ª solución

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

sublistasConProducto2 :: Integer -> [Integer] -> [[Integer]]

sublistasConProducto2 _ [] = []

sublistasConProducto2 n (x:xs)

| x > n = []

| x == n = [[x]]

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

++ sublistasConProducto2 n xs

| otherwise = sublistasConProducto2 n xs

where (q,r) = quotRem n x

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sublistasConProducto :: Integer -> [Integer] -> Bool

prop_sublistasConProducto n xs =

sort (sublistasConProducto n' xs') == sublistasConProducto2 n' xs'

where n' = 2 + abs n

xs' = (nub . sort . map ((+2) . abs)) xs

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=30}) prop_sublistasConProducto

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> sublistasConProducto 15 [1..23]

-- [[3,5],[1,3,5],[15],[1,15]]

-- (3.44 secs, 7,885,411,472 bytes)

-- λ> sublistasConProducto2 15 [1..23]

-- [[1,3,5],[1,15],[3,5],[15]]

-- (0.01 secs, 135,056 bytes)

-- λ> length (sublistasConProducto2 1234567 [1..1234567])

-- 4

-- (1.49 secs, 1,054,380,480 bytes)

-- Definición de unifactorizables

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

unifactorizables :: [Integer]

unifactorizables =

[n | n <- [2..]

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

Pensamiento

Y en el encinar,
¡luna redonda y beata,
siempre conmigo a la par!
Cerca de Úbeda la grande,
cuyos cerros nadie verá,
me iba siguiendo la luna
sobre el olivar.
Una luna jadeante,
siempre conmigo a la par.

Antonio Machado

Transformaciones lineales de números triangulares

PorJosé A. Alonso 4-diciembre-201911-diciembre-2019

La sucesión de los números triangulares se obtiene sumando los números naturales. Así, los 8 primeros números triangulares son

    1 = 1
    3 = 1+2
    6 = 1+2+3
   10 = 1+2+3+4
   15 = 1+2+3+4+5
   21 = 1+2+3+4+5+6
   28 = 1+2+3+4+5+6+7
   36 = 1+2+3+4+5+6+7+8

1 = 1

3 = 1+2

6 = 1+2+3

10 = 1+2+3+4

15 = 1+2+3+4+5

21 = 1+2+3+4+5+6

28 = 1+2+3+4+5+6+7

36 = 1+2+3+4+5+6+7+8

Para cada número triangular n existen números naturales a y b, tales que a . n + b también es triangular. Para n = 6, se tiene que

    6 = 1 * 6 + 0
   15 = 2 * 6 + 3
   21 = 3 * 6 + 3
   28 = 4 * 6 + 4
   36 = 5 * 6 + 6

6 = 1 * 6 + 0

15 = 2 * 6 + 3

21 = 3 * 6 + 3

28 = 4 * 6 + 4

36 = 5 * 6 + 6

son números triangulares

Definir la función

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

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

tal que si n es triangular, (transformaciones n) es la lista de los pares (a,b) tales que a es un entero positivo y b el menor número tal que a . n + b es triangular. Por ejemplo,

   take 5 (transformaciones 6)  == [(1,0),(2,3),(3,3),(4,4),(5,6)]
   take 5 (transformaciones 15) == [(1,0),(2,6),(3,10),(4,6),(5,3)]
   transformaciones 21 !! (7*10^7) == (70000001,39732)

take 5 (transformaciones 6) == [(1,0),(2,3),(3,3),(4,4),(5,6)]

take 5 (transformaciones 15) == [(1,0),(2,6),(3,10),(4,6),(5,3)]

transformaciones 21 !! (7*10^7) == (70000001,39732)

Soluciones

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

transformaciones :: Integer -> [(Integer,Integer)]
transformaciones n = (1,0) : [(a, f a) | a <- [2..]]
  where f a = head (dropWhile (<= a*n) triangulares) - a*n

-- triangulares es la lista de los números triangulares. Por ejemplo,  
--    take 5 triangulares == [1,3,6,10,15]
triangulares :: [Integer]
triangulares = scanl1 (+) [1..]

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

transformaciones2 :: Integer -> [(Integer,Integer)]
transformaciones2 n = (1,0): map g [2..]
  where g a = (a, head (dropWhile (<= a*n) triangulares) - a*n)

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

--    λ> transformaciones 21 !! (2*10^7)
--    (20000001,21615)
--    (3.02 secs, 4,320,111,544 bytes)
--    λ> transformaciones2 21 !! (2*10^7)
--    (20000001,21615)
--    (0.44 secs, 3,200,112,320 bytes)
--
--    λ> transformaciones2 21 !! (7*10^7)
--    (70000001,39732)
--    (1.41 secs, 11,200,885,336 bytes)

-- 1ª solución

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

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

transformaciones n = (1,0) : [(a, f a) | a <- [2..]]

where f a = head (dropWhile (<= a*n) triangulares) - a*n

-- triangulares es la lista de los números triangulares. Por ejemplo,

-- take 5 triangulares == [1,3,6,10,15]

triangulares :: [Integer]

triangulares = scanl1 (+) [1..]

-- 2ª solución

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

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

transformaciones2 n = (1,0): map g [2..]

where g a = (a, head (dropWhile (<= a*n) triangulares) - a*n)

-- Comparación de eficiencia

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

-- λ> transformaciones 21 !! (2*10^7)

-- (20000001,21615)

-- (3.02 secs, 4,320,111,544 bytes)

-- λ> transformaciones2 21 !! (2*10^7)

-- (20000001,21615)

-- (0.44 secs, 3,200,112,320 bytes)

-- λ> transformaciones2 21 !! (7*10^7)

-- (70000001,39732)

-- (1.41 secs, 11,200,885,336 bytes)

Pensamiento

A la hora del rocío,
de la niebla salen
sierra blanca y prado verde.
¡El sol en los encinares!

Antonio Machado

Múltiplos con ceros y unos

PorJosé A. Alonso 3-diciembre-201925-marzo-2022

Se observa que todos los primeros números naturales tienen al menos un múltiplo no nulo que está formado solamente por ceros y unos. Por ejemplo, 1×10=10, 2×5=10, 3×37=111, 4×25=100, 5×2=10, 6×185=1110; 7×143=1001; 8X125=1000; 9×12345679=111111111.

Definir la función

   multiplosCon1y0 :: Integer -> [Integer]

1	multiplosCon1y0 :: Integer -> [Integer]

tal que (multiplosCon1y0 n) es la lista de los múltiplos de n cuyos dígitos son 1 ó 0. Por ejemplo,

   take 4 (multiplosCon1y0 3)      ==  [111,1011,1101,1110]
   take 3 (multiplosCon1y0 23)     ==  [110101,1011011,1101010]
   head (multiplosCon1y0 1234658)  ==  110101101101000000110

take 4 (multiplosCon1y0 3) == [111,1011,1101,1110]

take 3 (multiplosCon1y0 23) == [110101,1011011,1101010]

head (multiplosCon1y0 1234658) == 110101101101000000110

Comprobar con QuickCheck que todo entero positivo tiene algún múltiplo cuyos dígitos son 1 ó 0.

Soluciones

import Test.QuickCheck

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

multiplosCon1y0 :: Integer -> [Integer]
multiplosCon1y0 n = [x | x <- multiplos n
                       , todos1y0 x]

-- (multiplos n) es la lista de los múltiplos de n. Por ejemplo, 
--    take 12 (multiplos 5)  ==  [5,10,15,20,25,30,35,40,45,50,55,60]
multiplos :: Integer -> [Integer]
multiplos n = [n,2*n..]

-- (todos1y0 n) se verifica si todos los dígitos de n son el 1 o el
-- 0. Por ejmplo,
--    todos1y0 1101110  ==  True
--    todos1y0 1102110  ==  False
todos1y0 :: Integer -> Bool
todos1y0 n = all (`elem` "01") (show n)

-- 2ª definición
-- =============
    
multiplosCon1y0b :: Integer -> [Integer] 
multiplosCon1y0b n = 
    [x | x <- numerosCon1y0
       , x `rem` n == 0] 

-- numerosCon1y0 es la lista de los números cuyos dígitos son 1 ó 0. Por
-- ejemplo,  
--    ghci> take 15 numerosCon1y0
--    [1,10,11,100,101,110,111,1000,1001,1010,1011,1100,1101,1110,1111]
numerosCon1y0 :: [Integer]
numerosCon1y0 = 1 : concat [[10*x,10*x+1] | x <- numerosCon1y0]

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

--    λ> head (multiplosCon1y0 9)
--    111111111
--    (7.70 secs, 10,853,320,456 bytes)
--    λ> head (multiplosCon1y0b 9)
--    111111111
--    (0.01 secs, 167,992 bytes)
                        
-- Comprobación de la propiedad
-- ============================
                        
-- La propiedad es
prop_existe_multiplosCon1y0 :: Integer -> Property
prop_existe_multiplosCon1y0 n = 
    n > 0 ==> (not . null) (multiplosCon1y0b n)

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

import Test.QuickCheck

-- 1ª definición

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

multiplosCon1y0 :: Integer -> [Integer]

multiplosCon1y0 n = [x | x <- multiplos n

, todos1y0 x]

-- (multiplos n) es la lista de los múltiplos de n. Por ejemplo,

-- take 12 (multiplos 5) == [5,10,15,20,25,30,35,40,45,50,55,60]

multiplos :: Integer -> [Integer]

multiplos n = [n,2*n..]

-- (todos1y0 n) se verifica si todos los dígitos de n son el 1 o el

-- 0. Por ejmplo,

-- todos1y0 1101110 == True

-- todos1y0 1102110 == False

todos1y0 :: Integer -> Bool

todos1y0 n = all (`elem` "01") (show n)

-- 2ª definición

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

multiplosCon1y0b :: Integer -> [Integer]

multiplosCon1y0b n =

[x | x <- numerosCon1y0

, x `rem` n == 0]

-- numerosCon1y0 es la lista de los números cuyos dígitos son 1 ó 0. Por

-- ejemplo,

-- ghci> take 15 numerosCon1y0

-- [1,10,11,100,101,110,111,1000,1001,1010,1011,1100,1101,1110,1111]

numerosCon1y0 :: [Integer]

numerosCon1y0 = 1 : concat [[10*x,10*x+1] | x <- numerosCon1y0]

-- Comparación de eficiencia

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

-- λ> head (multiplosCon1y0 9)

-- 111111111

-- (7.70 secs, 10,853,320,456 bytes)

-- λ> head (multiplosCon1y0b 9)

-- 111111111

-- (0.01 secs, 167,992 bytes)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_existe_multiplosCon1y0 :: Integer -> Property

prop_existe_multiplosCon1y0 n =

n > 0 ==> (not . null) (multiplosCon1y0b n)

-- La comprobación es

-- λ> quickCheck prop_existe_multiplosCon1y0

-- +++ OK, passed 100 tests.

Pensamiento

Huye del triste amor, amor pacato,
sin peligro, sin venda ni aventura,
que espera del amor prenda segura,
porque en amor locura es lo sensato.

Antonio Machado

Múltiplos palíndromos

PorJosé A. Alonso 2-diciembre-201911-diciembre-2019

Los números 545, 5995 y 15151 son los tres menores palíndromos (capicúas) que son divisibles por 109.

Definir las funciones

   multiplosPalindromos :: Integer -> [Integer]
   multiplosPalindromosMenores :: Integer -> Integer -> [Integer]

1 2	multiplosPalindromos :: Integer -> [Integer] multiplosPalindromosMenores :: Integer -> Integer -> [Integer]

tales que

(multiplosPalindromos n) es la lista de los palíndromos divisibles por n. Por ejemplo,

     take 5 (multiplosPalindromos 109) == [545,5995,15151,64746,74447]

1	take 5 (multiplosPalindromos 109) == [545,5995,15151,64746,74447]

(multiplosPalindromosMenores x n) es la lista de los palíndromos divisibles por n, menores que x. Por ejemplo,

     λ> multiplosPalindromosMenores (10^5) 109
     [545,5995,15151,64746,74447,79897,84148,89598,99299]

1 2	λ> multiplosPalindromosMenores (10^5) 109 [545,5995,15151,64746,74447,79897,84148,89598,99299]

Nota: Este ejercicio está basado en el problema 655 del Proyecto Euler.

Soluciones

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

multiplosPalindromos :: Integer -> [Integer]
multiplosPalindromos n =
  [x | x <- multiplos n
     , esPalindromo x]

-- (esPalindromo n) se verifica si n es palíndromo. Por ejemplo, 
--    esPalindromo 32523  ==  True
--    esPalindromo 32533  ==  False
esPalindromo :: Integer -> Bool
esPalindromo n = reverse xs == xs
  where xs = show n

-- (multiplos n) es la lista de los múltiplos de n. Por ejemplo, 
--    take 12 (multiplos 5)  ==  [5,10,15,20,25,30,35,40,45,50,55,60]
multiplos :: Integer -> [Integer]
multiplos n = [n,2*n..]

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

multiplosPalindromos2 :: Integer -> [Integer]
multiplosPalindromos2 = filter esPalindromo . multiplos 

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

multiplosPalindromosMenores :: Integer -> Integer -> [Integer]
multiplosPalindromosMenores x n =
  takeWhile (<x) (multiplosPalindromos n)

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

multiplosPalindromosMenores2 :: Integer -> Integer -> [Integer]
multiplosPalindromosMenores2 x =
  takeWhile (<x) . multiplosPalindromos

-- 1ª definición de multiplosPalindromos

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

multiplosPalindromos :: Integer -> [Integer]

multiplosPalindromos n =

[x | x <- multiplos n

, esPalindromo x]

-- (esPalindromo n) se verifica si n es palíndromo. Por ejemplo,

-- esPalindromo 32523 == True

-- esPalindromo 32533 == False

esPalindromo :: Integer -> Bool

esPalindromo n = reverse xs == xs

where xs = show n

-- (multiplos n) es la lista de los múltiplos de n. Por ejemplo,

-- take 12 (multiplos 5) == [5,10,15,20,25,30,35,40,45,50,55,60]

multiplos :: Integer -> [Integer]

multiplos n = [n,2*n..]

-- 2ª definición de multiplosPalindromos

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

multiplosPalindromos2 :: Integer -> [Integer]

multiplosPalindromos2 = filter esPalindromo . multiplos

-- 1ª definición de multiplosPalindromosMenores

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

multiplosPalindromosMenores :: Integer -> Integer -> [Integer]

multiplosPalindromosMenores x n =

takeWhile (<x) (multiplosPalindromos n)

-- 2ª definición de multiplosPalindromosMenores

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

multiplosPalindromosMenores2 :: Integer -> Integer -> [Integer]

multiplosPalindromosMenores2 x =

takeWhile (<x) . multiplosPalindromos

Pensamiento

Esta luz de Sevilla… Es el palacio
donde nací, con su rumor de fuente.

Antonio Machado

Soluciones

Referencias

Pensamiento

Soluciones

Referencias

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Referencia

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Entradas recientes

Buscar

Correo electrónico