map – Página 2 – Exercitium

Hojas con caminos no decrecientes

PorJosé A. Alonso 17-abril-202029-abril-2020

Los árboles se pueden representar mediante el siguiente tipo de datos

   data Arbol = N Int [Arbol]
     deriving Show

1 2	data Arbol = N Int [Arbol] deriving Show

Por ejemplo, los árboles

         1             1             1  
        /  \          / \           / \ 
       /    \        8   3         8   3
      2      6          /|\       /|\  |
     / \    / \        4 2 6     4 5 6 2
    4   5  5   7

1 1 1

/ \ / \ / \

/ \ 8 3 8 3

2 6 /|\ /|\ |

/ \ / \ 4 2 6 4 5 6 2

4 5 5 7

se representan por

   ej1, ej2, ej3 :: Arbol
   ej1 = N 1 [N 2 [N 4 [], N 5 []], N 6 [N 5 [], N 7 []]]
   ej2 = N 1 [N 8 [], N 3 [N 4 [], N 2 [], N 6 []]]
   ej3 = N 1 [N 8 [N 4 [], N 5 [], N 6 []], N 3 [N 2 []]]

ej1, ej2, ej3 :: Arbol

ej1 = N 1 [N 2 [N 4 [], N 5 []], N 6 [N 5 [], N 7 []]]

ej2 = N 1 [N 8 [], N 3 [N 4 [], N 2 [], N 6 []]]

ej3 = N 1 [N 8 [N 4 [], N 5 [], N 6 []], N 3 [N 2 []]]

Definir la función

   hojasEnNoDecreciente :: Arbol -> [Int]

1	hojasEnNoDecreciente :: Arbol -> [Int]

tal que (hojasEnNoDecreciente a) es el conjunto de las hojas de a que se encuentran en alguna rama no decreciente. Por ejemplo,

   hojasEnNoDecreciente ej1  ==  [4,5,7]
   hojasEnNoDecreciente ej2  ==  [4,6,8]
   hojasEnNoDecreciente ej3  ==  []

hojasEnNoDecreciente ej1 == [4,5,7]

hojasEnNoDecreciente ej2 == [4,6,8]

hojasEnNoDecreciente ej3 == []

Soluciones

import Data.List (sort, nub)

data Arbol = N Int [Arbol]
  deriving Show
           
ej1, ej2, ej3 :: Arbol
ej1 = N 1 [N 2 [N 4 [], N 5 []], N 6 [N 5 [], N 7 []]]
ej2 = N 1 [N 8 [], N 3 [N 4 [], N 2 [], N 6 []]]
ej3 = N 1 [N 8 [N 4 [], N 5 [], N 6 []], N 3 [N 2 []]]

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

hojasEnNoDecreciente :: Arbol -> [Int]
hojasEnNoDecreciente a =
  sort (nub (map last (ramasNoDecrecientes a)))

--    ramasNoDecrecientes ej1  ==  [[1,2,4],[1,2,5],[1,6,7]]
--    ramasNoDecrecientes ej2  ==  [[1,8],[1,3,4],[1,3,6]]
--    ramasNoDecrecientes ej3  ==  []
ramasNoDecrecientes :: Arbol -> [[Int]]
ramasNoDecrecientes a =
  filter esNoDecreciente (ramas a)

-- (ramas a) es la lista de las ramas del árbol a. Por ejemplo,
--    λ> ramas ej1
--    [[1,2,4],[1,2,5],[1,6,5],[1,6,7]]
--    λ> ramas ej2
--    [[1,8],[1,3,4],[1,3,2],[1,3,6]]
--    λ> ramas ej3
--    [[1,8,4],[1,8,5],[1,8,6],[1,3,2]]
ramas :: Arbol -> [[Int]]
ramas (N x []) = [[x]]
ramas (N x as) = map (x:) (concatMap ramas as)

-- (esNoDecreciente xs) se verifica si la lista xs es no
-- decreciente. Por ejemplo, 
--    esNoDecreciente [1,3,3,5]  ==  True
--    esNoDecreciente [1,3,5,3]  ==  False
esNoDecreciente :: [Int] -> Bool
esNoDecreciente xs =
  and (zipWith (<=) xs (tail xs))

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

--    hojasEnNoDecreciente ej1  ==  [4,5,7]
--    hojasEnNoDecreciente ej2  ==  [4,6,8]
--    hojasEnNoDecreciente ej3  ==  []
hojasEnNoDecreciente2 :: Arbol -> [Int]
hojasEnNoDecreciente2 = sort . nub . aux
  where
    aux (N x []) = [x]
    aux (N x as) = concat [aux (N y bs) | (N y bs) <- as, x <= y]

import Data.List (sort, nub)

data Arbol = N Int [Arbol]

deriving Show

ej1, ej2, ej3 :: Arbol

ej1 = N 1 [N 2 [N 4 [], N 5 []], N 6 [N 5 [], N 7 []]]

ej2 = N 1 [N 8 [], N 3 [N 4 [], N 2 [], N 6 []]]

ej3 = N 1 [N 8 [N 4 [], N 5 [], N 6 []], N 3 [N 2 []]]

-- 1ª solución

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

hojasEnNoDecreciente :: Arbol -> [Int]

hojasEnNoDecreciente a =

sort (nub (map last (ramasNoDecrecientes a)))

-- ramasNoDecrecientes ej1 == [[1,2,4],[1,2,5],[1,6,7]]

-- ramasNoDecrecientes ej2 == [[1,8],[1,3,4],[1,3,6]]

-- ramasNoDecrecientes ej3 == []

ramasNoDecrecientes :: Arbol -> [[Int]]

ramasNoDecrecientes a =

filter esNoDecreciente (ramas a)

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

-- λ> ramas ej1

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

-- λ> ramas ej2

-- [[1,8],[1,3,4],[1,3,2],[1,3,6]]

-- λ> ramas ej3

-- [[1,8,4],[1,8,5],[1,8,6],[1,3,2]]

ramas :: Arbol -> [[Int]]

ramas (N x []) = [[x]]

ramas (N x as) = map (x:) (concatMap ramas as)

-- (esNoDecreciente xs) se verifica si la lista xs es no

-- decreciente. Por ejemplo,

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

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

esNoDecreciente :: [Int] -> Bool

esNoDecreciente xs =

and (zipWith (<=) xs (tail xs))

-- 2ª solución

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

-- hojasEnNoDecreciente ej1 == [4,5,7]

-- hojasEnNoDecreciente ej2 == [4,6,8]

-- hojasEnNoDecreciente ej3 == []

hojasEnNoDecreciente2 :: Arbol -> [Int]

hojasEnNoDecreciente2 = sort . nub . aux

where

aux (N x []) = [x]

aux (N x as) = concat [aux (N y bs) | (N y bs) <- as, x <= y]

Otras soluciones

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

Orden de divisibilidad

PorJosé A. Alonso 14-abril-202021-abril-2020

El orden de divisibilidad de un número x es el mayor n tal que para todo i menor o igual que n, los i primeros dígitos de n es divisible por i. Por ejemplo, el orden de divisibilidad de 74156 es 3 porque

   7       es divisible por 1
   74      es divisible por 2
   741     es divisible por 3
   7415 no es divisible por 4

7 es divisible por 1

74 es divisible por 2

741 es divisible por 3

7415 no es divisible por 4

Definir la función

   ordenDeDivisibilidad :: Integer -> Int

1	ordenDeDivisibilidad :: Integer -> Int

tal que (ordenDeDivisibilidad x) es el orden de divisibilidad de x. Por ejemplo,

   ordenDeDivisibilidad 74156                      ==  3
   ordenDeDivisibilidad 12                         ==  2
   ordenDeDivisibilidad 7                          ==  1
   ordenDeDivisibilidad 3608528850368400786036725  ==  25

ordenDeDivisibilidad 74156 == 3

ordenDeDivisibilidad 12 == 2

ordenDeDivisibilidad 7 == 1

ordenDeDivisibilidad 3608528850368400786036725 == 25

Soluciones

import Data.List (inits)

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

ordenDeDivisibilidad :: Integer -> Int
ordenDeDivisibilidad n = 
  length (takeWhile (\(x,k) -> x `mod` k == 0) (zip (sucDigitos n) [1..]))

-- (sucDigitos x) es la sucesión de los dígitos de x. Por ejemplo,
--    sucDigitos 325    ==  [3,32,325]
--    sucDigitos 32050  ==  [3,32,320,3205,32050]
sucDigitos :: Integer -> [Integer]
sucDigitos n = 
    [n `div` (10^i) | i <- [k-1,k-2..0]]
    where k = length (show n)

-- 2ª definición de sucDigitos
sucDigitos2 :: Integer -> [Integer]
sucDigitos2 n = [read xs | xs <- aux (show n)]
  where aux []     = []
        aux (d:ds) = [d] : map (d:) (aux ds)

-- 3ª definición de sucDigitos
sucDigitos3 :: Integer -> [Integer]
sucDigitos3 n = 
  [read (take k ds) | k <- [1..length ds]]
  where ds = show n

-- 4ª definición de sucDigitos
sucDigitos4 :: Integer -> [Integer]
sucDigitos4 n = [read xs | xs <- tail (inits (show n))]

-- 5ª definición de sucDigitos
sucDigitos5 :: Integer -> [Integer]
sucDigitos5 n = map read (tail (inits (show n)))

-- 6ª definición de sucDigitos
sucDigitos6 :: Integer -> [Integer]
sucDigitos6 = map read . (tail . inits . show)

-- Eficiencia de las definiciones de sucDigitos
--    ghci> length (sucDigitos (10^5000))
--    5001
--    (0.01 secs, 1550688 bytes)
--    ghci> length (sucDigitos2 (10^5000))
--    5001
--    (1.25 secs, 729411872 bytes)
--    ghci> length (sucDigitos3 (10^5000))
--    5001
--    (0.02 secs, 2265120 bytes)
--    ghci> length (sucDigitos4 (10^5000))
--    5001
--    (1.10 secs, 728366872 bytes)
--    ghci> length (sucDigitos5 (10^5000))
--    5001
--    (1.12 secs, 728393864 bytes)
--    ghci> length (sucDigitos6 (10^5000))
--    5001
--    (1.20 secs, 728403052 bytes)
-- 
--    ghci> length (sucDigitos (10^3000000))
--    3000001
--    (2.73 secs, 820042696 bytes)
--    ghci> length (sucDigitos3 (10^3000000))
--    3000001
--    (3.69 secs, 820043688 bytes)

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

ordenDeDivisibilidad2 :: Integer -> Int
ordenDeDivisibilidad2 x =
  length
  $ takeWhile (==0)
  $ zipWith (mod . read) (tail $ inits $ show x) [1..]

import Data.List (inits)

-- 1ª definición de ordenDeDivisibilidad

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

ordenDeDivisibilidad :: Integer -> Int

ordenDeDivisibilidad n =

length (takeWhile (\(x,k) -> x `mod` k == 0) (zip (sucDigitos n) [1..]))

-- (sucDigitos x) es la sucesión de los dígitos de x. Por ejemplo,

-- sucDigitos 325 == [3,32,325]

-- sucDigitos 32050 == [3,32,320,3205,32050]

sucDigitos :: Integer -> [Integer]

sucDigitos n =

[n `div` (10^i) | i <- [k-1,k-2..0]]

where k = length (show n)

-- 2ª definición de sucDigitos

sucDigitos2 :: Integer -> [Integer]

sucDigitos2 n = [read xs | xs <- aux (show n)]

where aux [] = []

aux (d:ds) = [d] : map (d:) (aux ds)

-- 3ª definición de sucDigitos

sucDigitos3 :: Integer -> [Integer]

sucDigitos3 n =

[read (take k ds) | k <- [1..length ds]]

where ds = show n

-- 4ª definición de sucDigitos

sucDigitos4 :: Integer -> [Integer]

sucDigitos4 n = [read xs | xs <- tail (inits (show n))]

-- 5ª definición de sucDigitos

sucDigitos5 :: Integer -> [Integer]

sucDigitos5 n = map read (tail (inits (show n)))

-- 6ª definición de sucDigitos

sucDigitos6 :: Integer -> [Integer]

sucDigitos6 = map read . (tail . inits . show)

-- Eficiencia de las definiciones de sucDigitos

-- ghci> length (sucDigitos (10^5000))

-- 5001

-- (0.01 secs, 1550688 bytes)

-- ghci> length (sucDigitos2 (10^5000))

-- 5001

-- (1.25 secs, 729411872 bytes)

-- ghci> length (sucDigitos3 (10^5000))

-- 5001

-- (0.02 secs, 2265120 bytes)

-- ghci> length (sucDigitos4 (10^5000))

-- 5001

-- (1.10 secs, 728366872 bytes)

-- ghci> length (sucDigitos5 (10^5000))

-- 5001

-- (1.12 secs, 728393864 bytes)

-- ghci> length (sucDigitos6 (10^5000))

-- 5001

-- (1.20 secs, 728403052 bytes)

-- ghci> length (sucDigitos (10^3000000))

-- 3000001

-- (2.73 secs, 820042696 bytes)

-- ghci> length (sucDigitos3 (10^3000000))

-- 3000001

-- (3.69 secs, 820043688 bytes)

-- 2ª definición de ordenDeDivisibilidad

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

ordenDeDivisibilidad2 :: Integer -> Int

ordenDeDivisibilidad2 x =

length

$ takeWhile (==0)

$ zipWith (mod . read) (tail $ inits $ show x) [1..]

Otras soluciones

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

Cálculo de pi mediante la fracción continua de Lange

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

En 1999, L.J. Lange publicó el artículo An elegant new continued fraction for π.

En el primer teorema del artículo se demuestra la siguiente expresión de π mediante una fracción continua
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange$

La primeras aproximaciones son

   a(1) = 3+1                = 4.0
   a(2) = 3+(1/(6+9))        = 3.066666666666667
   a(3) = 3+(1/(6+9/(6+25))) = 3.158974358974359

a(1) = 3+1 = 4.0

a(2) = 3+(1/(6+9)) = 3.066666666666667

a(3) = 3+(1/(6+9/(6+25))) = 3.158974358974359

Definir las funciones

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

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

tales que

(aproximacionPi n) es la n-ésima aproximación de pi con la fracción continua de Lange. Por ejemplo,

     aproximacionPi 1     ==  4.0
     aproximacionPi 2     ==  3.066666666666667
     aproximacionPi 3     ==  3.158974358974359
     aproximacionPi 10    ==  3.141287132741557
     aproximacionPi 100   ==  3.141592398533554
     aproximacionPi 1000  ==  3.1415926533392926

aproximacionPi 1 == 4.0

aproximacionPi 2 == 3.066666666666667

aproximacionPi 3 == 3.158974358974359

aproximacionPi 10 == 3.141287132741557

aproximacionPi 100 == 3.141592398533554

aproximacionPi 1000 == 3.1415926533392926

(grafica xs) dibuja la gráfica de las k-ésimas aproximaciones de pi donde k toma los valores de la lista xs. Por ejemplo, (grafica [1..10]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_2$
(grafica [10..100]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_3$
y (grafica [100..200]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_4$

Soluciones

import Graphics.Gnuplot.Simple
 
-- fraccionPi es la representación de la fracción continua de pi como un
-- par de listas infinitas.
fraccionPi :: [(Integer, Integer)]
fraccionPi = zip (3 : [6,6..]) (map (^2) [1,3..])

-- (aproximacionFC n fc) es la n-ésima aproximación de la fracción
-- continua fc (como un par de listas).  
aproximacionFC :: Int -> [(Integer, Integer)] -> Double
aproximacionFC n =
  foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take n

aproximacionPi :: Int -> Double
aproximacionPi n =
  aproximacionFC n fraccionPi
 
grafica :: [Int] -> IO ()
grafica xs = 
  plotList [Key Nothing]
           [(k,aproximacionPi k) | k <- xs]

import Graphics.Gnuplot.Simple

-- fraccionPi es la representación de la fracción continua de pi como un

-- par de listas infinitas.

fraccionPi :: [(Integer, Integer)]

fraccionPi = zip (3 : [6,6..]) (map (^2) [1,3..])

-- (aproximacionFC n fc) es la n-ésima aproximación de la fracción

-- continua fc (como un par de listas).

aproximacionFC :: Int -> [(Integer, Integer)] -> Double

aproximacionFC n =

foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take n

aproximacionPi :: Int -> Double

aproximacionPi n =

aproximacionFC n fraccionPi

grafica :: [Int] -> IO ()

grafica xs =

plotList [Key Nothing]

[(k,aproximacionPi k) | k <- xs]

Otras soluciones

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

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

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

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

El desarrollo es el siguiente

y se obtiene a partir de la serie armónica

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

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

Por ejemplo,

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

Definir las funciones

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

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

tales que

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

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

aproximacionPi 1 == 1.0

aproximacionPi 10 == 2.3289682539682537

aproximacionPi 100 == 2.934318000847734

aproximacionPi 1000 == 3.0603246224585128

aproximacionPi 10000 == 3.1105295744825403

aproximacionPi 100000 == 3.134308801935256

aproximacionPi 1000000 == 3.1395057903490806

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

Soluciones

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

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

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

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

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

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

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

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

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

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

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

-- 1ª definición

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

aproximacionPi :: Int -> Double

aproximacionPi n =

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

signoPrimo :: Int -> Int

signoPrimo 2 = 1

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

| otherwise = -1

signo :: Int -> Int

signo n | isPrime n = signoPrimo n

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

-- 2ª definición

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

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = serieEuler !! (n-1)

serieEuler :: [Double]

serieEuler =

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

-- Definición de grafica

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

grafica :: Int -> IO ()

grafica n =

plotList [Key Nothing]

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

Otras soluciones

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

Conjetura de Goldbach

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

Una forma de la conjetura de Golbach afirma que todo entero mayor que 1 se puede escribir como la suma de uno, dos o tres números primos.

Si se define el índice de Goldbach de n > 1 como la mínima cantidad de primos necesarios para que su suma sea n, entonces la conjetura de Goldbach afirma que todos los índices de Goldbach de los enteros mayores que 1 son menores que 4.

Definir las siguientes funciones

   indiceGoldbach  :: Int -> Int
   graficaGoldbach :: Int -> IO ()

1 2	indiceGoldbach :: Int -> Int graficaGoldbach :: Int -> IO ()

tales que

(indiceGoldbach n) es el índice de Goldbach de n. Por ejemplo,

     indiceGoldbach 2                        ==  1
     indiceGoldbach 4                        ==  2
     indiceGoldbach 27                       ==  3
     sum (map indiceGoldbach [2..5000])      ==  10619
     maximum (map indiceGoldbach [2..5000])  ==  3

indiceGoldbach 2 == 1

indiceGoldbach 4 == 2

indiceGoldbach 27 == 3

sum (map indiceGoldbach [2..5000]) == 10619

maximum (map indiceGoldbach [2..5000]) == 3

(graficaGoldbach n) dibuja la gráfica de los índices de Goldbach de los números entre 2 y n. Por ejemplo, (graficaGoldbach 150) dibuja

Comprobar con QuickCheck la conjetura de Goldbach anterior.

Soluciones

import Data.Array
import Data.Numbers.Primes
import Graphics.Gnuplot.Simple
import Test.QuickCheck


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

indiceGoldbach :: Int -> Int
indiceGoldbach n =
  minimum (map length (particiones n))

particiones :: Int -> [[Int]]
particiones n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
    where f 0 = [[]]
          f m = [x:y | x <- xs, 
                       y <- v ! (m-x), 
                       [x] >= take 1 y]
            where xs = reverse (takeWhile (<= m) primes)

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

indiceGoldbach2 :: Int -> Int
indiceGoldbach2 x =
  head [n | n <- [1..], esSumaDe x n]

-- (esSumaDe x n) se verifica si x se puede escribir como la suma de n
-- primos. Por ejemplo,
--    esSumaDe 2  1  ==  True
--    esSumaDe 4  1  ==  False
--    esSumaDe 4  2  ==  True
--    esSumaDe 27 2  ==  False
--    esSumaDe 27 3  ==  True
esSumaDe :: Int -> Int -> Bool
esSumaDe x 1 = isPrime x
esSumaDe x n = or [esSumaDe (x-p) (n-1) | p <- takeWhile (<= x) primes]

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

indiceGoldbach3 :: Int -> Int
indiceGoldbach3 x =
  head [n | n <- [1..], esSumaDe3 x n]

esSumaDe3 :: Int -> Int -> Bool
esSumaDe3 x n = a ! (x,n) where
  a = array ((2,1),(x,9)) [((i,j),f i j) | i <- [2..x], j <- [1..9]]
  f i 1 = isPrime i
  f i j = or [a!(i-k,j-1) | k <- takeWhile (<= i) primes]

-- 4ª definición
-- =============

indiceGoldbach4 :: Int -> Int
indiceGoldbach4 n = v ! n where
  v = array (2,n) [(i,f i) | i <- [2..n]]
  f i | isPrime i = 1
      | otherwise = 1 + minimum [v!(i-p) | p <- takeWhile (< (i-1)) primes]

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

--    λ> sum (map indiceGoldbach [2..80])
--    142
--    (2.66 secs, 1,194,330,496 bytes)
--    λ> sum (map indiceGoldbach2 [2..80])
--    142
--    (0.01 secs, 1,689,944 bytes)
--    λ> sum (map indiceGoldbach3 [2..80])
--    142
--    (0.03 secs, 27,319,296 bytes)
--    λ> sum (map indiceGoldbach4 [2..80])
--    142
--    (0.03 secs, 47,823,656 bytes)
--    
--    λ> sum (map indiceGoldbach2 [2..1000])
--    2030
--    (0.10 secs, 200,140,264 bytes)
--    λ> sum (map indiceGoldbach3 [2..1000])
--    2030
--    (3.10 secs, 4,687,467,664 bytes)

-- Gráfica
-- =======

graficaGoldbach :: Int -> IO ()
graficaGoldbach n =
  plotList [ Key Nothing
           , XRange (2,fromIntegral n)
           , PNG ("Conjetura_de_Goldbach_" ++ show n ++ ".png")
           ]
           [indiceGoldbach2 k | k <- [2..n]]

-- Comprobación de la conjetura de Goldbach
-- ========================================

-- La propiedad es
prop_Goldbach :: Int -> Property
prop_Goldbach x =
  x >= 2 ==> indiceGoldbach2 x < 4

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

100

101

102

103

104

105

106

107

import Data.Array

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

import Test.QuickCheck

-- 1ª definición

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

indiceGoldbach :: Int -> Int

indiceGoldbach n =

minimum (map length (particiones n))

particiones :: Int -> [[Int]]

particiones n = v ! n where

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

where f 0 = [[]]

f m = [x:y | x <- xs,

y <- v ! (m-x),

[x] >= take 1 y]

where xs = reverse (takeWhile (<= m) primes)

-- 2ª definición

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

indiceGoldbach2 :: Int -> Int

indiceGoldbach2 x =

head [n | n <- [1..], esSumaDe x n]

-- (esSumaDe x n) se verifica si x se puede escribir como la suma de n

-- primos. Por ejemplo,

-- esSumaDe 2 1 == True

-- esSumaDe 4 1 == False

-- esSumaDe 4 2 == True

-- esSumaDe 27 2 == False

-- esSumaDe 27 3 == True

esSumaDe :: Int -> Int -> Bool

esSumaDe x 1 = isPrime x

esSumaDe x n = or [esSumaDe (x-p) (n-1) | p <- takeWhile (<= x) primes]

-- 3ª definición

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

indiceGoldbach3 :: Int -> Int

indiceGoldbach3 x =

head [n | n <- [1..], esSumaDe3 x n]

esSumaDe3 :: Int -> Int -> Bool

esSumaDe3 x n = a ! (x,n) where

a = array ((2,1),(x,9)) [((i,j),f i j) | i <- [2..x], j <- [1..9]]

f i 1 = isPrime i

f i j = or [a!(i-k,j-1) | k <- takeWhile (<= i) primes]

-- 4ª definición

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

indiceGoldbach4 :: Int -> Int

indiceGoldbach4 n = v ! n where

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

f i | isPrime i = 1

| otherwise = 1 + minimum [v!(i-p) | p <- takeWhile (< (i-1)) primes]

-- Comparación de eficiencia

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

-- λ> sum (map indiceGoldbach [2..80])

-- 142

-- (2.66 secs, 1,194,330,496 bytes)

-- λ> sum (map indiceGoldbach2 [2..80])

-- 142

-- (0.01 secs, 1,689,944 bytes)

-- λ> sum (map indiceGoldbach3 [2..80])

-- 142

-- (0.03 secs, 27,319,296 bytes)

-- λ> sum (map indiceGoldbach4 [2..80])

-- 142

-- (0.03 secs, 47,823,656 bytes)

-- λ> sum (map indiceGoldbach2 [2..1000])

-- 2030

-- (0.10 secs, 200,140,264 bytes)

-- λ> sum (map indiceGoldbach3 [2..1000])

-- 2030

-- (3.10 secs, 4,687,467,664 bytes)

-- Gráfica

-- =======

graficaGoldbach :: Int -> IO ()

graficaGoldbach n =

plotList [ Key Nothing

, XRange (2,fromIntegral n)

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

]

[indiceGoldbach2 k | k <- [2..n]]

-- Comprobación de la conjetura de Goldbach

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

-- La propiedad es

prop_Goldbach :: Int -> Property

prop_Goldbach x =

x >= 2 ==> indiceGoldbach2 x < 4

-- La comprobación es

-- λ> quickCheck prop_Goldbach

-- +++ OK, passed 100 tests.

Otras soluciones

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

Pensamiento

«La diferencia entre los matemáticos y los físicos es que después de que los físicos prueban un gran resultado piensan que es fantástico, pero después de que los matemáticos prueban un gran resultado piensan que es trivial.»

Lucien Szpiro.

Repeticiones consecutivas

PorJosé A. Alonso 17-marzo-202024-marzo-2020

Se dice que una palabra tiene una repetición en una frase si es igual a una, o más, de las palabras consecutivas sin distinguir mayúsculas de minúsculas.

Definir la función

   nRepeticionesConsecutivas :: String ->Int

1	nRepeticionesConsecutivas :: String ->Int

tal que (nRepeticionesConsecutivas cs) es el número de repeticiones de palabras consecutivas de la cadena cs. Por ejemplo,

   nRepeticionesConsecutivas "oso rana"                    == 0      
   nRepeticionesConsecutivas "oso rana oso"                == 0
   nRepeticionesConsecutivas "oso oSo rana"                == 1
   nRepeticionesConsecutivas "oso oso oso rana"            == 1
   nRepeticionesConsecutivas "coronavirus virus oso rana"  == 0
   nRepeticionesConsecutivas "virus     virus oso rana"    == 1
   nRepeticionesConsecutivas "virus oso virus oso rana"    == 0
   nRepeticionesConsecutivas "oso oso oso oso oso oso"     == 1
   nRepeticionesConsecutivas "oso oso oso oso rana rana"   == 2
   nRepeticionesConsecutivas "rana rana oso oso rana rana" == 3

nRepeticionesConsecutivas "oso rana" == 0

nRepeticionesConsecutivas "oso rana oso" == 0

nRepeticionesConsecutivas "oso oSo rana" == 1

nRepeticionesConsecutivas "oso oso oso rana" == 1

nRepeticionesConsecutivas "coronavirus virus oso rana" == 0

nRepeticionesConsecutivas "virus virus oso rana" == 1

nRepeticionesConsecutivas "virus oso virus oso rana" == 0

nRepeticionesConsecutivas "oso oso oso oso oso oso" == 1

nRepeticionesConsecutivas "oso oso oso oso rana rana" == 2

nRepeticionesConsecutivas "rana rana oso oso rana rana" == 3

Soluciones

import Data.List (group)
import Data.Char (toUpper)

-- 1ª solución
nRepeticionesConsecutivas :: String ->Int
nRepeticionesConsecutivas = aux . words . map toUpper 
  where aux (x:y:zs) | x == y    = 1 + aux (dropWhile (== x) zs)
                     | otherwise = aux (y:zs)
        aux _ = 0

-- 2ª solución
nRepeticionesConsecutivas2 :: String ->Int
nRepeticionesConsecutivas2 cs =
  length [xs | xs <- group (words (map toUpper cs)), length xs > 1]

-- 3ª solución
nRepeticionesConsecutivas3 :: String ->Int
nRepeticionesConsecutivas3 =
  length . filter ((>1) . length) . group . words . map toUpper

import Data.List (group)

import Data.Char (toUpper)

-- 1ª solución

nRepeticionesConsecutivas :: String ->Int

nRepeticionesConsecutivas = aux . words . map toUpper

where aux (x:y:zs) | x == y = 1 + aux (dropWhile (== x) zs)

| otherwise = aux (y:zs)

aux _ = 0

-- 2ª solución

nRepeticionesConsecutivas2 :: String ->Int

nRepeticionesConsecutivas2 cs =

length [xs | xs <- group (words (map toUpper cs)), length xs > 1]

-- 3ª solución

nRepeticionesConsecutivas3 :: String ->Int

nRepeticionesConsecutivas3 =

length . filter ((>1) . length) . group . words . map toUpper

Otras soluciones

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

Pensamiento

«En el campo de la computación, el momento de la verdad es la ejecución de un programa; todo lo demás es profecía.»

Herbert A. Simon.

Medias de dígitos de pi

PorJosé A. Alonso 16-marzo-202023-marzo-2020

El fichero Digitos_de_pi.txt contiene el número pi con un millón de decimales; es decir,

   3.1415926535897932384626433832 ... 83996346460422090106105779458151

1	3.1415926535897932384626433832 ... 83996346460422090106105779458151

Definir las funciones

   mediasDigitosDePi        :: IO [Double]
   graficaMediasDigitosDePi :: Int -> IO ()

1 2	mediasDigitosDePi :: IO [Double] graficaMediasDigitosDePi :: Int -> IO ()

tales que

mediasDigitosDePi es la sucesión cuyo n-ésimo elemento es la media de los n primeros dígitos de pi. Por ejemplo,

     λ> xs <- mediasDigitosDePi
     λ> take 5 xs
     [1.0,2.5,2.0,2.75,4.0]
     λ> take 10 xs
     [1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]
     λ> take 10 <$> mediasDigitosDePi
     [1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]

λ> xs <- mediasDigitosDePi

λ> take 5 xs

[1.0,2.5,2.0,2.75,4.0]

λ> take 10 xs

[1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]

λ> take 10 <$> mediasDigitosDePi

[1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]

(graficaMediasDigitosDePi n) dibuja la gráfica de los n primeros términos de mediasDigitosDePi. Por ejemplo,
- (graficaMediasDigitosDePi 20) dibuja
- (graficaMediasDigitosDePi 200) dibuja
- (graficaMediasDigitosDePi 2000) dibuja

Soluciones

import Data.Char (digitToInt)
import Data.List (genericLength, inits, tails)
import Graphics.Gnuplot.Simple ( Attribute (Key, PNG)
                               , plotList )

-- Definición de mediasDigitosDePi
-- ===============================

mediasDigitosDePi :: IO [Double]
mediasDigitosDePi = do
  (_:_:ds) <- readFile "Digitos_de_pi.txt"
  return (medias (digitos ds))

-- (digitos cs) es la lista de los digitos de cs. Por ejempplo,
--    digitos "1415926535"  ==  [1,4,1,5,9,2,6,5,3,5]
digitos :: String -> [Int]
digitos = map digitToInt

-- (medias xs) es la lista de las medias de los segmentos iniciales de
-- xs. Por ejemplo,
--    λ> medias [1,4,1,5,9,2,6,5,3,5]
--    [1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]
medias :: [Int] -> [Double]
medias xs = map media (tail (inits xs))

-- (media xs) es la media aritmética de xs. Por ejemplo,
--    media [1,4,1,5,9]  ==  4.0
media :: [Int] -> Double
media xs = fromIntegral (sum xs) / genericLength xs

-- Definición de graficaMediasDigitosDePi
-- ======================================

graficaMediasDigitosDePi :: Int -> IO ()
graficaMediasDigitosDePi n = do
  xs <- mediasDigitosDePi
  plotList [ Key Nothing
           , PNG ("Medias_de_digitos_de_pi_" ++ show n ++ ".png")
           ]
           (take n xs)

import Data.Char (digitToInt)

import Data.List (genericLength, inits, tails)

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

, plotList )

-- Definición de mediasDigitosDePi

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

mediasDigitosDePi :: IO [Double]

mediasDigitosDePi = do

(_:_:ds) <- readFile "Digitos_de_pi.txt"

return (medias (digitos ds))

-- (digitos cs) es la lista de los digitos de cs. Por ejempplo,

-- digitos "1415926535" == [1,4,1,5,9,2,6,5,3,5]

digitos :: String -> [Int]

digitos = map digitToInt

-- (medias xs) es la lista de las medias de los segmentos iniciales de

-- xs. Por ejemplo,

-- λ> medias [1,4,1,5,9,2,6,5,3,5]

-- [1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]

medias :: [Int] -> [Double]

medias xs = map media (tail (inits xs))

-- (media xs) es la media aritmética de xs. Por ejemplo,

-- media [1,4,1,5,9] == 4.0

media :: [Int] -> Double

media xs = fromIntegral (sum xs) / genericLength xs

-- Definición de graficaMediasDigitosDePi

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

graficaMediasDigitosDePi :: Int -> IO ()

graficaMediasDigitosDePi n = do

xs <- mediasDigitosDePi

plotList [ Key Nothing

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

]

(take n xs)

Otras soluciones

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

Pensamiento

Es el mejor de los buenos
quien sabe que en esta vida
todo es cuestión de medida:
un poco más, algo menos.

Antonio Machado

Reducción de SAT a Clique

PorJosé A. Alonso 27-febrero-20205-marzo-2020

Nota: En este ejercicio se usa la misma notación que en los anteriores importando los módulos

+ Evaluacion_de_FNC
+ Modelos_de_FNC
+ Problema_SAT_para_FNC
+ Cliques
+ KCliques
+ Grafo_FNC

+ Evaluacion_de_FNC

+ Modelos_de_FNC

+ Problema_SAT_para_FNC

+ Cliques

+ KCliques

+ Grafo_FNC

Definir las funciones

   cliquesFNC :: FNC -> [[(Int,Literal)]]
   cliquesCompletos :: FNC -> [[(Int,Literal)]]
   esSatisfaciblePorClique :: FNC -> Bool

cliquesFNC :: FNC -> [[(Int,Literal)]]

cliquesCompletos :: FNC -> [[(Int,Literal)]]

esSatisfaciblePorClique :: FNC -> Bool

tales que

(cliquesFNCf) es la lista de los cliques del grafo de f. Por ejemplo,

     λ> cliquesFNC [[1,-2,3],[-1,2],[-2,3]]
     [[], [(0,1)], [(1,2)], [(0,1),(1,2)], [(2,-2)],
      [(0,1),(2,-2)], [(2,3)], [(0,1),(2,3)], [(1,2),(2,3)],
      [(0,1),(1,2),(2,3)], [(0,-2)], [(2,-2),(0,-2)], [(2,3),(0,-2)],
      [(1,-1)], [(2,-2),(1,-1)], [(2,3),(1,-1)], [(0,-2),(1,-1)],
      [(2,-2),(0,-2),(1,-1)], [(2,3),(0,-2),(1,-1)], [(0,3)],
      [(1,2),(0,3)], [(2,-2),(0,3)], [(2,3),(0,3)],
      [(1,2),(2,3),(0,3)], [(1,-1),(0,3)],
      [(2,-2),(1,-1),(0,3)], [(2,3),(1,-1),(0,3)]]

λ> cliquesFNC [[1,-2,3],[-1,2],[-2,3]]

[[], [(0,1)], [(1,2)], [(0,1),(1,2)], [(2,-2)],

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

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

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

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

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

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

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

(cliquesCompletos f) es la lista de los cliques del grafo de f que tiene tantos elementos como cláusulas tiene f. Por ejemplo,

     λ> cliquesCompletos [[1,-2,3],[-1,2],[-2,3]]
     [[(0,1),(1,2),(2,3)],   [(2,-2),(0,-2),(1,-1)],
      [(2,3),(0,-2),(1,-1)], [(1,2),(2,3),(0,3)],
      [(2,-2),(1,-1),(0,3)], [(2,3),(1,-1),(0,3)]]
     λ> cliquesCompletos [[1,2],[1,-2],[-1,2],[-1,-2]]
     []

λ> cliquesCompletos [[1,-2,3],[-1,2],[-2,3]]

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

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

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

λ> cliquesCompletos [[1,2],[1,-2],[-1,2],[-1,-2]]

[]

(esSatisfaciblePorClique f) se verifica si f no contiene la cláusula vacía, tiene más de una cláusula y posee algún clique completo. Por ejemplo,

     λ> esSatisfaciblePorClique [[1,-2,3],[-1,2],[-2,3]]
     True
     λ> esSatisfaciblePorClique [[1,2],[1,-2],[-1,2],[-1,-2]]
     False

λ> esSatisfaciblePorClique [[1,-2,3],[-1,2],[-2,3]]

True

λ> esSatisfaciblePorClique [[1,2],[1,-2],[-1,2],[-1,-2]]

False

Comprobar con QuickCheck que toda fórmula en FNC es satisfacible si, y solo si, es satisfacible por Clique.

Soluciones

module Reduccion_de_SAT_a_Clique where

import Evaluacion_de_FNC
import Modelos_de_FNC
import Problema_SAT_para_FNC
import Cliques
import KCliques
import Grafo_FNC
import Data.List (nub, sort)
import Test.QuickCheck

cliquesFNC :: FNC -> [[(Int,Literal)]]
cliquesFNC f = cliques (grafoFNC f)

cliquesCompletos :: FNC -> [[(Int,Literal)]]
cliquesCompletos cs = kCliques (grafoFNC cs) (length cs)

esSatisfaciblePorClique :: FNC -> Bool
esSatisfaciblePorClique f =
     [] `notElem` f'
  && (length f' <= 1 || not (null (cliquesCompletos f')))
  where f' = nub (map (nub . sort) f) 

-- La propiedad es
prop_esSatisfaciblePorClique :: FNC -> Bool
prop_esSatisfaciblePorClique f =
  esSatisfacible f == esSatisfaciblePorClique f

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

module Reduccion_de_SAT_a_Clique where

import Evaluacion_de_FNC

import Modelos_de_FNC

import Problema_SAT_para_FNC

import Cliques

import KCliques

import Grafo_FNC

import Data.List (nub, sort)

import Test.QuickCheck

cliquesFNC :: FNC -> [[(Int,Literal)]]

cliquesFNC f = cliques (grafoFNC f)

cliquesCompletos :: FNC -> [[(Int,Literal)]]

cliquesCompletos cs = kCliques (grafoFNC cs) (length cs)

esSatisfaciblePorClique :: FNC -> Bool

esSatisfaciblePorClique f =

[] `notElem` f'

&& (length f' <= 1 || not (null (cliquesCompletos f')))

where f' = nub (map (nub . sort) f)

-- La propiedad es

prop_esSatisfaciblePorClique :: FNC -> Bool

prop_esSatisfaciblePorClique f =

esSatisfacible f == esSatisfaciblePorClique f

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=7}) prop_esSatisfaciblePorClique

-- +++ OK, passed 100 tests.

Otras soluciones

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

Pensamiento

«La resolución de problemas es una habilidad práctica como, digamos, la natación. Adquirimos cualquier habilidad práctica por imitación y práctica. Tratando de nadar, imitas lo que otras personas hacen con sus manos y pies para mantener sus cabezas sobre el agua, y, finalmente, aprendes a nadar practicando la natación. Al intentar resolver problemas, hay que observar e imitar lo que hacen otras personas al resolver problemas y, finalmente, se aprende a resolver problemas haciéndolos.»

George Pólya.

Átomos de FNC (fórmulas en forma normal conjuntiva)

PorJosé A. Alonso 13-febrero-202020-febrero-2020

Nota: En este ejercicio usaremos las mismas notaciones que en el anterior importando el módulo Evaluacion_de_FNC.

Definir las siguientes funciones

    atomosClausula :: Clausula -> [Atomo]
    atomosFNC      :: FNC -> [Atomo]

1 2	atomosClausula :: Clausula -> [Atomo] atomosFNC :: FNC -> [Atomo]

tales que

(atomosClausula c) es el conjunto de los átomos de la cláusula c. Por ejemplo,

     atomosClausula [3,1,-3] == [1,3]

1	atomosClausula [3,1,-3] == [1,3]

(atomosFNC f) es el conjunto de los átomos de la FNC f. Por ejemplo,

   atomosFNC [[4,5],[1,-2],[-4,-1,-5]]  ==  [1,2,4,5]

1	atomosFNC [[4,5],[1,-2],[-4,-1,-5]] == [1,2,4,5]

Nota: Escribir la solución en el módulo Atomos_de_FNC para poderlo usar en los siguientes ejercicios.

Soluciones

module Atomos_de_FNC where
  
import Evaluacion_de_FNC
import Data.List (sort, nub)

-- 1ª definición de atomosClausula
atomosClausula :: Clausula -> [Atomo]
atomosClausula c = sort (nub (map abs c))

-- 2ª definición de atomosClausula
atomosClausula2 :: Clausula -> [Atomo]
atomosClausula2 = sort . nub . map abs

-- 1ª definición de atomosFNC
atomosFNC :: FNC -> [Atomo]
atomosFNC f = sort (nub (concat [atomosClausula c | c <- f]))

-- 2ª definición
atomosFNC2 :: FNC -> [Atomo]
atomosFNC2 f = sort (nub (concatMap atomosClausula f))

-- 3ª definición
atomosFNC3 :: FNC -> [Atomo]
atomosFNC3 = sort . nub . concatMap atomosClausula

module Atomos_de_FNC where

import Evaluacion_de_FNC

import Data.List (sort, nub)

-- 1ª definición de atomosClausula

atomosClausula :: Clausula -> [Atomo]

atomosClausula c = sort (nub (map abs c))

-- 2ª definición de atomosClausula

atomosClausula2 :: Clausula -> [Atomo]

atomosClausula2 = sort . nub . map abs

-- 1ª definición de atomosFNC

atomosFNC :: FNC -> [Atomo]

atomosFNC f = sort (nub (concat [atomosClausula c | c <- f]))

-- 2ª definición

atomosFNC2 :: FNC -> [Atomo]

atomosFNC2 f = sort (nub (concatMap atomosClausula f))

-- 3ª definición

atomosFNC3 :: FNC -> [Atomo]

atomosFNC3 = sort . nub . concatMap atomosClausula

Otras soluciones

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

Pensamiento

«La esencia de las matemáticas es su libertad.»

Georg Cantor.

Máximo número de consecutivos iguales al dado

PorJosé A. Alonso 4-febrero-202011-febrero-2020

Definir la función

   maximoConsecutivosIguales :: Eq a => a -> [a] -> Int

1	maximoConsecutivosIguales :: Eq a => a -> [a] -> Int

tal que (maximoConsecutivosIguales x xs) es el mayor número de elementos consecutivos en xs iguales a x. Por ejemplo,

   maximoConsecutivosIguales 'b' "abbcccbbbd"    ==  3
   maximoConsecutivosIguales 'b' "abbbbcccbbbd"  ==  4
   maximoConsecutivosIguales 'e' "abbcccbbbd"    ==  0

maximoConsecutivosIguales 'b' "abbcccbbbd" == 3

maximoConsecutivosIguales 'b' "abbbbcccbbbd" == 4

maximoConsecutivosIguales 'e' "abbcccbbbd" == 0

Soluciones

import Data.List (group)

maximoConsecutivosIguales :: Eq a => a -> [a] -> Int
maximoConsecutivosIguales x = maximum
                            . (0:)
                            . map length
                            . filter ((== x) . head)
                            . group

import Data.List (group)

maximoConsecutivosIguales :: Eq a => a -> [a] -> Int

maximoConsecutivosIguales x = maximum

. (0:)

. map length

. filter ((== x) . head)

. group

Otras soluciones

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

Pensamiento

«La programación de computadoras es un arte, porque aplica el conocimiento
acumulado al mundo, porque requiere habilidad e ingenio, y especialmente
porque produce belleza. Un programador que subconscientemente se ve
a sí mismo como un artista disfrutará con lo que hace y lo hará mejor.»

Donald Knuth.

Acotación del primorial

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

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

   2 * 3 * 5 = 30

1	2 * 3 * 5 = 30

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

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

Definir las funciones

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

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

tales que

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

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

primorial 3 == 6

primorial 5 == 30

primorial 8 == 210

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

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

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

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

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

import Data.Numbers.Primes

import Test.QuickCheck

-- 1ª definición de primorial

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

primorial :: Integer -> Integer

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

-- 2ª definición de primorial

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

primorial2 :: Integer -> Integer

primorial2 0 = 1

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

| otherwise = x

where x = primorial2 (n-1)

-- Comparación de eficiencia

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

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

-- 216852

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

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

-- 216852

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

-- 1ª definición de primoriales

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

-- λ> take 15 primoriales

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

primoriales :: [Integer]

primoriales = map primorial [0..]

-- 2ª definición de primoriales

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

-- λ> take 15 primoriales2

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

primoriales2 :: [Integer]

primoriales2 = map primorial2 [0..]

-- 3ª definición de primoriales

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

-- λ> take 15 primoriales3

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

primoriales3 :: [Integer]

primoriales3 = scanl1 f [1..]

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

| otherwise = x

-- Comparación de eficiencia

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

-- λ> minimum (take 5000 primoriales)

-- 1

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

-- λ> minimum (take 5000 primoriales2)

-- 1

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

-- λ> minimum (take 5000 primoriales3)

-- 1

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

-- λ> minimum (take 6000 primoriales)

-- 1

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

-- λ> minimum (take 6000 primoriales3)

-- 1

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

-- Propiedad

-- =========

prop_primorial :: Integer -> Property

prop_primorial n =

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

-- La comprobación es

-- λ> quickCheck prop_primorial

-- +++ OK, passed 100 tests.

Otras soluciones

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

Pensamiento

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

Carl Friedrich Gauss.

Longitud de la parte periódica

PorJosé A. Alonso 30-enero-202026-marzo-2022

La propiedad de la longitud de la parte periódica afirma que

Si p es un número primo distinto de 2 y de 5, entonces la longitud del período de 1/p es el menor entero positivo n tal que p divide a $10^n - 1$ .

El objetivo de este ejercicio es la verificación de dicha propiedad.

Las fracciones se representan por un par de enteros. Por ejemplo, el número 2/3 se representa por (2,3). Su tipo es

   type Fraccion = (Integer,Integer)

1	type Fraccion = (Integer,Integer)

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

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

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

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

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

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

Su tipo es

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

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

Definir, usando las funciones cocientesRestos y primerRepetido de los ejercicios anteriores, las funciones

   decimal :: Fraccion -> Decimal
   longitudPeriodo :: Fraccion -> Integer

1 2	decimal :: Fraccion -> Decimal longitudPeriodo :: Fraccion -> Integer

tales que

(decimal f) es la representación decimal de la fracción f. Por ejemplo,

     decimal (6,2)          ==  (3,[],[])
     decimal (3,4)          ==  (0,[7,5],[])
     decimal (1,3)          ==  (0,[],[3])
     decimal (23,14)        ==  (1,[6],[4,2,8,5,7,1])
     decimal (247813,19980) ==  (12,[4,0],[3,0,5])
     decimal (1,101)        ==  (0,[],[0,0,9,9])

decimal (6,2) == (3,[],[])

decimal (3,4) == (0,[7,5],[])

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

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

decimal (247813,19980) == (12,[4,0],[3,0,5])

decimal (1,101) == (0,[],[0,0,9,9])

(longitudPeriodo f) es la longitud de la parte periódica de la representación decimal de la fracción f. Por ejemplo,

     longitudPeriodo (6,2)           ==  0
     longitudPeriodo (3,4)           ==  0
     longitudPeriodo (1,3)           ==  1
     longitudPeriodo (23,14)         ==  6
     longitudPeriodo (247813,19980)  ==  3
     longitudPeriodo (1,101)         ==  4
     longitudPeriodo (1,1229)        ==  1228

longitudPeriodo (6,2) == 0

longitudPeriodo (3,4) == 0

longitudPeriodo (1,3) == 1

longitudPeriodo (23,14) == 6

longitudPeriodo (247813,19980) == 3

longitudPeriodo (1,101) == 4

longitudPeriodo (1,1229) == 1228

Comprobar con QuickCheck la propiedad de la longitud de la parte periódica; es decir, k es un número natural distinto de 0 y 2 y p es el primo k-ésimo, entonces la longitud del período de 1/p es el menor entero positivo n tal que p divide a $10^n - 1$ ..

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck

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

decimal :: Fraccion -> Decimal
decimal (n,d) 
  | snd y == 0 = (fst x, map fst xs, [])
  | otherwise  = (fst x, map fst xs, map fst (y:zs))
  where
    qrs         = cocientesRestos (n,d)
    Just (q,r)  = primerRepetido qrs
    (x:xs,y:ys) = break (==(q,r)) qrs
    zs          = takeWhile (/=(q,r)) ys

cocientesRestos :: Fraccion -> [(Integer,Integer)]
cocientesRestos (n,d) =
  (q,r) : cocientesRestos (10*r, d)
  where (q,r) = quotRem n d

primerRepetido :: Eq a => [a] -> Maybe a
primerRepetido xs = aux xs []
  where
    aux [] _                     = Nothing
    aux (x:xs') ys | x `elem` ys = Just x
                   | otherwise   = aux xs' (x:ys) 

longitudPeriodo :: Fraccion -> Int
longitudPeriodo (n,d) = length xs
  where (_,_,xs) = decimal (n,d)

-- La propiedad es
prop_LongitudPeriodo :: Int -> Property
prop_LongitudPeriodo k =
  k > 0 && k /= 2 
  ==>
  longitudPeriodo (1,p) ==
  head [n | n <- [1..], (10^n-1) `mod` p == 0]
  where p = primes !! k

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

import Data.Numbers.Primes

import Test.QuickCheck

type Fraccion = (Integer,Integer)

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

decimal :: Fraccion -> Decimal

decimal (n,d)

| snd y == 0 = (fst x, map fst xs, [])

| otherwise = (fst x, map fst xs, map fst (y:zs))

where

qrs = cocientesRestos (n,d)

Just (q,r) = primerRepetido qrs

(x:xs,y:ys) = break (==(q,r)) qrs

zs = takeWhile (/=(q,r)) ys

cocientesRestos :: Fraccion -> [(Integer,Integer)]

cocientesRestos (n,d) =

(q,r) : cocientesRestos (10*r, d)

where (q,r) = quotRem n d

primerRepetido :: Eq a => [a] -> Maybe a

primerRepetido xs = aux xs []

where

aux [] _ = Nothing

aux (x:xs') ys | x `elem` ys = Just x

| otherwise = aux xs' (x:ys)

longitudPeriodo :: Fraccion -> Int

longitudPeriodo (n,d) = length xs

where (_,_,xs) = decimal (n,d)

-- La propiedad es

prop_LongitudPeriodo :: Int -> Property

prop_LongitudPeriodo k =

k > 0 && k /= 2

==>

longitudPeriodo (1,p) ==

head [n | n <- [1..], (10^n-1) `mod` p == 0]

where p = primes !! k

-- La comprobación es

-- λ> quickCheck prop_LongitudPeriodo

-- +++ OK, passed 100 tests.

Otras soluciones

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

Pensamiento

«En el desarrollo de la comprensión de los fenómenos complejos, la herramienta más poderosa de que dispone el intelecto humano es la abstracción. La abstracción surge del reconocimiento de las similitudes entre ciertos objetos, situaciones o procesos en el mundo real y de la decisión de concentrarse en estas similitudes e ignorar, por el momento, sus diferencias.»

Tony Hoare

La conjetura de Mertens

PorJosé A. Alonso 27-enero-20203-febrero-2020

Un número entero n es libre de cuadrados si no existe un número primo p tal que p² divide a n; es decir, los factores primos de n son todos distintos.

La función de Möbius μ(n) está definida para todos los enteros positivos como sigue:

μ(n) = 1 si n es libre de cuadrados y tiene un número par de factores primos.
μ(n) = -1 si n es libre de cuadrados y tiene un número impar de factores primos.
μ(n) = 0 si n no es libre de cuadrados.

Sus primeros valores son 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, …

La función de Mertens M(n) está definida para todos los enteros positivos como la suma de μ(k) para 1 ≤ k ≤ n. Sus primeros valores son 1, 0, -1, -1, -2, -1, -2, -2, …

La conjetura de Mertens afirma que

Para todo entero x mayor que 1, el valor absoluto de la función de Mertens en x es menor que la raíz cuadrada de x.

La conjetura fue planteada por Franz Mertens en 1897. Riele Odlyzko, demostraronen 1985 que la conjetura de Mertens deja de ser cierta más o menos a partir de $10^{10^{64}}$ , cifra que luego de algunos refinamientos se redujo a $10^{10^{40}}$ .

Definir las funciones

   mobius :: Integer -> Integer
   mertens :: Integer -> Integer
   graficaMertens :: Integer -> IO ()

mobius :: Integer -> Integer

mertens :: Integer -> Integer

graficaMertens :: Integer -> IO ()

tales que

(mobius n) es el valor de la función de Möbius en n. Por ejemplo,

     mobius 6   ==  1
     mobius 30  ==  -1
     mobius 12  ==  0

mobius 6 == 1

mobius 30 == -1

mobius 12 == 0

(mertens n) es el valor de la función de Mertens en n. Por ejemplo,

     mertens 1     ==  1
     mertens 2     ==  0
     mertens 3     ==  -1
     mertens 5     ==  -2
     mertens 661   ==  -11
     mertens 1403  ==  11

mertens 1 == 1

mertens 2 == 0

mertens 3 == -1

mertens 5 == -2

mertens 661 == -11

mertens 1403 == 11

(graficaMertens n) dibuja la gráfica de la función de Mertens, la raíz cuadrada y el opuestos de la raíz cuadrada para los n primeros n enteros positivos. Por ejemplo, (graficaMertens 1000) dibuja

Comprobar con QuickCheck la conjetura de Mertens.

Nota: El ejercicio está basado en La conjetura de Merterns y su relación con un número tan raro como extremada y colosalmente grande publicado por @Alvy la semana pasada en Microsiervos.

Soluciones

import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck
import Graphics.Gnuplot.Simple

mobius :: Integer -> Integer
mobius n | tieneRepetidos xs = 0
         | otherwise         = (-1)^(length xs)
  where xs = primeFactors n

tieneRepetidos :: [Integer] -> Bool
tieneRepetidos xs =
  or [x == y | (x,y) <- zip xs (tail xs)]
        
mertens :: Integer -> Integer
mertens n = sum (map mobius [1..n])

-- Definición de graficaMertens
-- ============================

graficaMertens :: Integer -> IO ()
graficaMertens n = do
  plotLists [ Key Nothing
            , Title "Conjetura de Mertens"
            , PNG "La_conjetura_de_Mertens.png"
            ]
            [ [mertens k | k <- [1..n]]
            , raices
            , map negate raices
            ]
    
  where
    raices = [ceiling (sqrt k) | k <- [1..fromIntegral n]]

-- Conjetura de Mertens
-- ====================

-- La conjetura es
conjeturaDeMertens :: Integer -> Property
conjeturaDeMertens n =
  n > 1
  ==>
  abs (mertens n) < ceiling (sqrt n')
  where n' = fromIntegral n

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

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

import Graphics.Gnuplot.Simple

mobius :: Integer -> Integer

mobius n | tieneRepetidos xs = 0

| otherwise = (-1)^(length xs)

where xs = primeFactors n

tieneRepetidos :: [Integer] -> Bool

tieneRepetidos xs =

or [x == y | (x,y) <- zip xs (tail xs)]

mertens :: Integer -> Integer

mertens n = sum (map mobius [1..n])

-- Definición de graficaMertens

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

graficaMertens :: Integer -> IO ()

graficaMertens n = do

plotLists [ Key Nothing

, Title "Conjetura de Mertens"

, PNG "La_conjetura_de_Mertens.png"

]

[ [mertens k | k <- [1..n]]

, raices

, map negate raices

]

where

raices = [ceiling (sqrt k) | k <- [1..fromIntegral n]]

-- Conjetura de Mertens

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

-- La conjetura es

conjeturaDeMertens :: Integer -> Property

conjeturaDeMertens n =

n > 1

==>

abs (mertens n) < ceiling (sqrt n')

where n' = fromIntegral n

-- La comprobación es

-- λ> quickCheck conjeturaDeMertens

-- +++ OK, passed 100 tests.

Otras soluciones

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

Pensamiento

«El control de la complejidad es la esencia de la programación informática.»

Brian Kernighan.

Teorema de la amistad

PorJosé A. Alonso 17-enero-202024-enero-2020

El teorema de la amistad afirma que

En cualquier reunión de n personas hay al menos dos personas que tienen el mismo número de amigos (suponiendo que la relación de amistad es simétrica).

Se pueden usar las siguientes representaciones:

números enteros para representar a las personas,
pares de enteros (x,y), con x < y, para representar que la persona x e y son amigas y
lista de pares de enteros para representar la reunión junto con las relaciones de amistad.

Por ejemplo, [(2,3),(3,5)] representa una reunión de tres personas
(2, 3 y 5) donde

2 es amiga de 3,
3 es amiga de 2 y 5 y
5 es amiga de 3.
Si clasificamos las personas poniendo en la misma clase las que tienen el mismo número de amigos, se obtiene [[2,5],[3]] ya que 2 y 5 tienen 1 amigo y 3 tiene 2 amigos.

Definir la función

   clasesAmigos :: [(Int,Int)] -> [[Int]]

1	clasesAmigos :: [(Int,Int)] -> [[Int]]

tal que (clasesAmigos r) es la clasificación según el número de amigos de las personas de la reunión r; es decir, la lista cuyos elementos son las listas de personas con 1 amigo, con 2 amigos y así hasta que se completa todas las personas de la reunión r. Por ejemplo,

   clasesAmigos [(2,3),(3,5)]            ==  [[2,5],[3]]
   clasesAmigos [(2,3),(4,5)]            ==  [[2,3,4,5]]
   clasesAmigos [(2,3),(2,5),(3,5)]      ==  [[2,3,5]]
   clasesAmigos [(2,3),(3,4),(2,5)]      ==  [[4,5],[2,3]]
   clasesAmigos [(x,x+1) | x <- [1..5]]  ==  [[1,6],[2,3,4,5]]
   length (clasesAmigos [(x,x+1) | x <- [1..2020]]) == 2

clasesAmigos [(2,3),(3,5)] == [[2,5],[3]]

clasesAmigos [(2,3),(4,5)] == [[2,3,4,5]]

clasesAmigos [(2,3),(2,5),(3,5)] == [[2,3,5]]

clasesAmigos [(2,3),(3,4),(2,5)] == [[4,5],[2,3]]

clasesAmigos [(x,x+1) | x <- [1..5]] == [[1,6],[2,3,4,5]]

length (clasesAmigos [(x,x+1) | x <- [1..2020]]) == 2

Comprobar con QuickCheck el teorema de la amistad; es decir, si r es una lista de pares de enteros, entonces (clasesAmigos r’) donde r’ es la lista de los pares (x,y) de r con x < y y se supone que r’ es no vacía.

Soluciones

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

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

clasesAmigos :: [(Int,Int)] -> [[Int]]
clasesAmigos ps =
  filter (not . null)
         [[x | x <- xs, numeroDeAmigos ps x == n] | n <- [1..length xs]] 
  where xs = personas ps
    
-- (personas ps) es la lista de personas en la reunión ps. Por ejemplo,
--    personas [(2,3),(3,5)]  ==  [2,3,5]
personas :: [(Int,Int)] -> [Int]
personas ps = sort (nub (map fst ps ++ map snd ps))

-- (numeroDeAmigos ps x) es el número de amigos de x en la reunión
-- ps. Por ejemplo, 
--    numeroDeAmigos [(2,3),(3,5)] 2  ==  1
--    numeroDeAmigos [(2,3),(3,5)] 3  ==  2
--    numeroDeAmigos [(2,3),(3,5)] 5  ==  1
numeroDeAmigos :: [(Int,Int)] -> Int -> Int
numeroDeAmigos ps x = length (amigos ps x)

-- (amigos ps x) es la lista de los amigos de x en la reunión ps. Por
-- ejemplo, 
--    amigos [(2,3),(3,5)] 2  ==  [3]
--    amigos [(2,3),(3,5)] 3  ==  [5,2]
--    amigos [(2,3),(3,5)] 5  ==  [3]
amigos :: [(Int,Int)] -> Int -> [Int]
amigos ps x =
  nub ([b | (a,b) <- ps, a == x] ++ [a | (a,b) <- ps, b == x])

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

clasesAmigos2 :: [(Int,Int)] -> [[Int]]
clasesAmigos2 = clases . sort . tablaAmigos
  where
    clases [] = []
    clases ps@((x,y):ps') = (map snd (takeWhile (\(a,b) -> a == x) ps)) :
                            clases (dropWhile (\(a,b) -> a == x) ps')

-- (tablaAmigos ps) es la lista de pares (a,b) tales que b es una
-- persona de la reunión ps y a es su número de amigos. Por ejemplo,
--    tablaAmigos [(2,3),(3,5)]   ==  [(1,2),(2,3),(1,5)]
tablaAmigos :: [(Int,Int)] -> [(Int,Int)]
tablaAmigos ps = [(numeroDeAmigos ps x,x) | x <- personas ps]

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

-- La propiedad es
prop_equivalencia :: [(Int,Int)] -> Property
prop_equivalencia ps =
  not (null ps')
  ==> 
  clasesAmigos ps' == clasesAmigos2 ps'
  where ps' = [(x,y) | (x,y) <- ps, x < y]

-- La comprobación es
--    λ> quickCheck prop_equivalencia
--    +++ OK, passed 100 tests.
--    (1.06 secs, 337,106,752 bytes)

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

--    λ> length (clasesAmigos [(x,x+1) | x <- [1..200]]) 
--    2
--    (2.37 secs, 804,402,848 bytes)
--    λ> length (clasesAmigos2 [(x,x+1) | x <- [1..200]]) 
--    2
--    (0.02 secs, 4,287,256 bytes)

-- El teorema de la amistad
-- ========================

-- La propiedad es
teoremaDeLaAmistad :: [(Int,Int)] -> Property
teoremaDeLaAmistad ps =
  not (null ps')
  ==> 
  not (null [xs | xs <- clasesAmigos2 ps', length xs > 1])
  where ps' = [(x,y) | (x,y) <- ps, x < y]

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

import Data.List (nub, sort)

import Test.QuickCheck

-- 1ª solución

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

clasesAmigos :: [(Int,Int)] -> [[Int]]

clasesAmigos ps =

filter (not . null)

[[x | x <- xs, numeroDeAmigos ps x == n] | n <- [1..length xs]]

where xs = personas ps

-- (personas ps) es la lista de personas en la reunión ps. Por ejemplo,

-- personas [(2,3),(3,5)] == [2,3,5]

personas :: [(Int,Int)] -> [Int]

personas ps = sort (nub (map fst ps ++ map snd ps))

-- (numeroDeAmigos ps x) es el número de amigos de x en la reunión

-- ps. Por ejemplo,

-- numeroDeAmigos [(2,3),(3,5)] 2 == 1

-- numeroDeAmigos [(2,3),(3,5)] 3 == 2

-- numeroDeAmigos [(2,3),(3,5)] 5 == 1

numeroDeAmigos :: [(Int,Int)] -> Int -> Int

numeroDeAmigos ps x = length (amigos ps x)

-- (amigos ps x) es la lista de los amigos de x en la reunión ps. Por

-- ejemplo,

-- amigos [(2,3),(3,5)] 2 == [3]

-- amigos [(2,3),(3,5)] 3 == [5,2]

-- amigos [(2,3),(3,5)] 5 == [3]

amigos :: [(Int,Int)] -> Int -> [Int]

amigos ps x =

nub ([b | (a,b) <- ps, a == x] ++ [a | (a,b) <- ps, b == x])

-- 2ª solución

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

clasesAmigos2 :: [(Int,Int)] -> [[Int]]

clasesAmigos2 = clases . sort . tablaAmigos

where

clases [] = []

clases ps@((x,y):ps') = (map snd (takeWhile (\(a,b) -> a == x) ps)) :

clases (dropWhile (\(a,b) -> a == x) ps')

-- (tablaAmigos ps) es la lista de pares (a,b) tales que b es una

-- persona de la reunión ps y a es su número de amigos. Por ejemplo,

-- tablaAmigos [(2,3),(3,5)] == [(1,2),(2,3),(1,5)]

tablaAmigos :: [(Int,Int)] -> [(Int,Int)]

tablaAmigos ps = [(numeroDeAmigos ps x,x) | x <- personas ps]

-- Equivalencia de las definiciones

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

-- La propiedad es

prop_equivalencia :: [(Int,Int)] -> Property

prop_equivalencia ps =

not (null ps')

==>

clasesAmigos ps' == clasesAmigos2 ps'

where ps' = [(x,y) | (x,y) <- ps, x < y]

-- La comprobación es

-- λ> quickCheck prop_equivalencia

-- +++ OK, passed 100 tests.

-- (1.06 secs, 337,106,752 bytes)

-- Comparación de eficiencia

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

-- λ> length (clasesAmigos [(x,x+1) | x <- [1..200]])

-- 2

-- (2.37 secs, 804,402,848 bytes)

-- λ> length (clasesAmigos2 [(x,x+1) | x <- [1..200]])

-- 2

-- (0.02 secs, 4,287,256 bytes)

-- El teorema de la amistad

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

-- La propiedad es

teoremaDeLaAmistad :: [(Int,Int)] -> Property

teoremaDeLaAmistad ps =

not (null ps')

==>

not (null [xs | xs <- clasesAmigos2 ps', length xs > 1])

where ps' = [(x,y) | (x,y) <- ps, x < y]

-- La comprobación es

-- λ> quickCheck teoremaDeLaAmistad

-- +++ OK, passed 100 tests.

Referencia

Pigeonhole among friends

Pensamiento

Me dijo el agua clara que reía,
bajo el sol, sobre el mármol de la fuente:
si te inquieta el enigma del presente
aprende el son de la salmodia mía.

Antonio Machado

Teorema de Hilbert-Waring

PorJosé A. Alonso 9-enero-202016-enero-2020

El problema de Waring, propuesto por Edward Waring consiste en déterminar si, para cada número entero k mayor que 1, existe un número n tal que todo entero positivo se puede escribir como una suma de k-potencias de números positivos con n sumandos como máximo.

La respuesta afirmativa al problema, aportada por David Hilbert, se conoce como el teorema de Hilbert-Waring. Su enunciado es

Para cada número entero k, con k ≥ 2, existe un entero positivo g(k) tal que todo entero positivo se puede expresar como una suma de a lo más g(k) k-ésimas potencias.

Definir las funciones

   descomposiciones :: Integer -> Integer -> Integer -> [[Integer]]
   orden :: Integer -> Integer -> Integer

1 2	descomposiciones :: Integer -> Integer -> Integer -> [[Integer]] orden :: Integer -> Integer -> Integer

tales que

(descomposiciones x k n) es la lista de descomposiciones de x como suma de n potencias con exponente k de números enteros positivos.

     descomposiciones 9   2 1  ==  [[9]]
     descomposiciones 9   3 1  ==  []
     descomposiciones 9   3 2  ==  [[1,8]]
     descomposiciones 9   4 9  ==  [[1,1,1,1,1,1,1,1,1]]
     descomposiciones 25  2 2  ==  [[9,16]]
     descomposiciones 133 2 3  ==  [[8,125]]
     descomposiciones 38  3 2  ==  [[1,1,36],[4,9,25]]

descomposiciones 9 2 1 == [[9]]

descomposiciones 9 3 1 == []

descomposiciones 9 3 2 == [[1,8]]

descomposiciones 9 4 9 == [[1,1,1,1,1,1,1,1,1]]

descomposiciones 25 2 2 == [[9,16]]

descomposiciones 133 2 3 == [[8,125]]

descomposiciones 38 3 2 == [[1,1,36],[4,9,25]]

(orden x k) es el menor número de sumandos necesario para expresar x como suma de k-ésimas potencias. Por ejemplo,

     orden 9  2  ==  1
     orden 9  3  ==  2
     orden 9  4  ==  9
     orden 10 2  ==  2
     orden 10 3  ==  3
     orden 10 4  ==  10
     [maximum [orden x k | x <- [1..1000]] | k <- [1..6]] == [1,4,9,19,37,73]

orden 9 2 == 1

orden 9 3 == 2

orden 9 4 == 9

orden 10 2 == 2

orden 10 3 == 3

orden 10 4 == 10

[maximum [orden x k | x <- [1..1000]] | k <- [1..6]] == [1,4,9,19,37,73]

Comprobar el teorema de Hilbert-Waring para k hasta 7; es decir, para todo número x positivo se verifica que

   orden x 2 <= 4   
   orden x 3 <= 9   
   orden x 4 <= 19  
   orden x 5 <= 37  
   orden x 6 <= 73  
   orden x 7 <= 143

orden x 2 <= 4

orden x 3 <= 9

orden x 4 <= 19

orden x 5 <= 37

orden x 6 <= 73

orden x 7 <= 143

y, en general,

   orden x k <= 2^k + floor ((3/2)^k) - 2

1	orden x k <= 2^k + floor ((3/2)^k) - 2

Soluciones

import Test.QuickCheck

descomposiciones :: Integer -> Integer -> Integer -> [[Integer]]
descomposiciones x k n =
  sumas x (takeWhile (<= x) (potencias k)) n

-- (potencias n) es la lista de las potencias de n
--    take 7 (potencias 2)  ==  [1,4,9,16,25,36,49]
--    take 7 (potencias 3)  ==  [1,8,27,64,125,216,343]
potencias :: Integer -> [Integer]
potencias n = map (^n) [1..]

-- (sumas n ys x) es la lista de las descomposiciones de x como
-- sumas de n sumandos de la lista creciente ys. Por ejemplo,
--    sumas 3 [1,2] 2  ==  [[1,2]]
--    sumas 4 [1,2] 2  ==  [[2,2]]
--    sumas 5 [1,2] 2  ==  []
--    sumas 5 [1,2] 3  ==  [[1,2,2]]
--    sumas 6 [1,2] 3  ==  [[2,2,2]]
--    sumas 6 [1,2,5] 2  ==  [[1,5]]
sumas :: Integer -> [Integer] -> Integer -> [[Integer]]
sumas _ [] _                   = []
sumas x ys 1     | x `elem` ys = [[x]]
                 | otherwise   = []
sumas x (y:ys) n | y > x       = []
                 | otherwise   = map (y:) (sumas (x-y) (y:ys) (n-1)) ++
                                 sumas x ys n

orden :: Integer -> Integer -> Integer
orden x k = head [n | n <- [1..]
                    , not (null (descomposiciones x k n))]

-- El teorema es                                 
teorema_Hilbert_Waring :: Integer -> Integer -> Property
teorema_Hilbert_Waring x k =
  x > 0 && k >= 2
  ==>
  orden x 2 <= 4   &&
  orden x 3 <= 9   &&
  orden x 4 <= 19  &&
  orden x 5 <= 37  &&
  orden x 6 <= 73  &&
  orden x 7 <= 143 &&
  orden x k <= 2^k + floor ((3/2)^k) - 2

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

import Test.QuickCheck

descomposiciones :: Integer -> Integer -> Integer -> [[Integer]]

descomposiciones x k n =

sumas x (takeWhile (<= x) (potencias k)) n

-- (potencias n) es la lista de las potencias de n

-- take 7 (potencias 2) == [1,4,9,16,25,36,49]

-- take 7 (potencias 3) == [1,8,27,64,125,216,343]

potencias :: Integer -> [Integer]

potencias n = map (^n) [1..]

-- (sumas n ys x) es la lista de las descomposiciones de x como

-- sumas de n sumandos de la lista creciente ys. Por ejemplo,

-- sumas 3 [1,2] 2 == [[1,2]]

-- sumas 4 [1,2] 2 == [[2,2]]

-- sumas 5 [1,2] 2 == []

-- sumas 5 [1,2] 3 == [[1,2,2]]

-- sumas 6 [1,2] 3 == [[2,2,2]]

-- sumas 6 [1,2,5] 2 == [[1,5]]

sumas :: Integer -> [Integer] -> Integer -> [[Integer]]

sumas _ [] _ = []

sumas x ys 1 | x `elem` ys = [[x]]

| otherwise = []

sumas x (y:ys) n | y > x = []

| otherwise = map (y:) (sumas (x-y) (y:ys) (n-1)) ++

sumas x ys n

orden :: Integer -> Integer -> Integer

orden x k = head [n | n <- [1..]

, not (null (descomposiciones x k n))]

-- El teorema es

teorema_Hilbert_Waring :: Integer -> Integer -> Property

teorema_Hilbert_Waring x k =

x > 0 && k >= 2

==>

orden x 2 <= 4 &&

orden x 3 <= 9 &&

orden x 4 <= 19 &&

orden x 5 <= 37 &&

orden x 6 <= 73 &&

orden x 7 <= 143 &&

orden x k <= 2^k + floor ((3/2)^k) - 2

-- La comprobación es

-- λ> quickCheck teorema_Hilbert_Waring

-- +++ OK, passed 100 tests.

Referencia

Basado en Hilbert-Waring theorem de ProofWiki.

Pensamiento

¡Y en la tersa arena,
cerca de la mar,
tu carne rosa y morena,
súbitamente Guiomar!

Antonio Machado

Teorema de Liouville sobre listas CuCu

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

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

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

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

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

   1, 2, 4, 5, 10, 20

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

puesto que el número de sus divisores es

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

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

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

Definir las funciones

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

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

tales que

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

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

esCuCu [1,2,3] == True

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

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

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

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

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

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

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

Comprobar con QuickCheck

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

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

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

esCuCu :: [Integer] -> Bool

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

-- 1ª definición de liouville

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

liouville :: Integer -> [Integer]

liouville n = map numeroDivisores (divisores n)

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

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

divisores :: Integer -> [Integer]

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

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

-- numeroDivisores 12 == 6

-- numeroDivisores 25 == 3

numeroDivisores :: Integer -> Integer

numeroDivisores n = genericLength (divisores n)

-- 2ª definición de liouville

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

liouville2 :: Integer -> [Integer]

liouville2 n = map numeroDivisores2 (divisores2 n)

-- Se usan las funciones

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

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

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

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

divisores2 :: Integer -> [Integer]

divisores2 = sort

. map (product . concat)

. sequence

. map inits

. group

. primeFactors

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

-- numeroDivisores2 12 == 6

-- numeroDivisores2 25 == 3

numeroDivisores2 :: Integer -> Integer

numeroDivisores2 =

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

-- Comparación de eficiencia

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

-- La comparación es

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

-- 540

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

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

-- 540

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

-- Propiedad

-- =========

-- La propiedad es

prop_Liouville :: Integer -> Property

prop_Liouville n =

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

-- La comprobación es

-- λ> quickCheck prop_Liouville

-- +++ OK, passed 100 tests.

-- Teorema de Liouville

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

-- La propiedad es

teorema_Liouville :: Integer -> Property

teorema_Liouville n =

n > 0 ==> esCuCu (liouville n)

-- La comprobación es

-- λ> quickCheck teorema_Liouville

-- +++ OK, passed 100 tests.

Pensamiento

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

Antonio Machado

Enumeración de conjuntos finitos de naturales

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

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

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

0: []

1: [0]

2: [1]

3: [1,0]

4: [2]

5: [2,0]

6: [2,1]

7: [2,1,0]

8: [3]

9: [3,0]

10: [3,1]

11: [3,1,0]

12: [3,2]

13: [3,2,0]

14: [3,2,1]

15: [3,2,1,0]

16: [4]

17: [4,0]

18: [4,1]

19: [4,1,0]

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

Definir la constante

   enumeracionCFN :: [[Integer]]

1	enumeracionCFN :: [[Integer]]

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

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

λ> take 20 enumeracionCFN

[[],

[0],

[1],[1,0],

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

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

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

Comprobar con QuickCheck que

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

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

import Data.List (genericLength, nub, sort)

import Test.QuickCheck

-- 1ª solución

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

enumeracionCFN :: [[Integer]]

enumeracionCFN = concatMap enumeracionCFNHasta [0..]

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

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

-- λ> enumeracionCFNHasta 1

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

-- λ> enumeracionCFNHasta 2

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

-- λ> enumeracionCFNHasta 3

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

enumeracionCFNHasta :: Integer -> [[Integer]]

enumeracionCFNHasta 0 = [[],[0]]

enumeracionCFNHasta n =

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

-- 2ª solución

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

enumeracionCFN2 :: [[Integer]]

enumeracionCFN2 = [] : aux 0 [[]]

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

where yss = map (n:) xs

-- 3ª solución

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

enumeracionCFN3 :: [[Integer]]

enumeracionCFN3 = map conjunto [0..]

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

-- conjunto 6 == [2,1]

conjunto :: Integer -> [Integer]

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

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

-- inverso). Por ejemplo,

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

binario :: Integer -> [Integer]

binario 0 = [0]

binario 1 = [1]

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

-- Comparación de eficiencia

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

-- Propiedades

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

-- La primera propiedad es

prop_enumeracionCFN :: Int -> Property

prop_enumeracionCFN n =

n >= 0 ==> xs < ys

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

-- La comprobación es

-- λ> quickCheck prop_enumeracionCFN

-- +++ OK, passed 100 tests.

-- La segunda propiedad es

prop_completa :: [Integer] -> Bool

prop_completa xs =

xs' `elem` enumeracionCFN

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

-- La comprobación es

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

-- +++ OK, passed 100 tests.

Pensamiento

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

Antonio Machado

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

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

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

Definir las funciones

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

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

tales que

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

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

aproximacionPi 5 == 3.140331156954753

aproximacionPi 10 == 3.1415914215112

aproximacionPi 15 == 3.141592652386592

aproximacionPi 20 == 3.1415926535886207

aproximacionPi 25 == 3.141592653589795

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

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

errorPi 0.1 == 2

errorPi 0.01 == 4

errorPi 0.001 == 6

errorPi 0.0001 == 7

errorPi 1e-4 == 7

errorPi 1e-14 == 24

pi == 3.141592653589793

aproximacionPi 24 == 3.1415926535897913

Soluciones

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

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

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

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

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

-- 1ª definición de aproximacionPi

aproximacionPi :: Int -> Double

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

where

aux 0 = 1

aux 1 = sqrt 2

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

-- 2ª definición de aproximacionPi

aproximacionPi2 :: Int -> Double

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

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

-- 3ª definición de aproximaxionPi

aproximacionPi3 :: Int -> Double

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

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

-- 1ª definición de errorPi

errorPi :: Double -> Int

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

, abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi

errorPi2 :: Double -> Int

errorPi2 x = until aceptable (+1) 1

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

Pensamiento

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

Antonio Machado

Pares definidos por su MCD y su MCM

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

Definir las siguientes funciones

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

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

tales que

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

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

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

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

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

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

pares 2 7 == []

pares 12 3 == []

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

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

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

nPares 3 3 == 1

nPares 4 12 == 2

nPares 2 12 == 4

nPares 2 60 == 8

nPares 2 7 == 0

nPares 12 3 == 0

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

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

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

import Data.Numbers.Primes (primeFactors)

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

import Test.QuickCheck

-- 1ª definición de pares

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

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

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

, y <- [1..b]

, gcd x y == a

, lcm x y == b]

-- 2ª definición de pares

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

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

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

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

, gcd x y == a

, lcm x y == b]

-- Comparación de eficiencia

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

-- 16

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

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

-- 16

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

-- 3ª definición de pares

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

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

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

, c `rem` x == 0

, let y = c `div` x

, gcd x y == a

]

where c = a * b

-- Comparacioń de eficiencia

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

-- 16

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

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

-- 16

-- (0.01 secs, 878,104 bytes)

-- 4ª definición de pares

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

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

-- factores primos

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

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

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

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

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

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

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

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

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

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

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

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

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

-- λ> b `div` a

-- 12

-- λ> primeFactors it

-- [2,2,3]

-- λ> group it

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

-- λ> subsequences it

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

-- λ> map concat it

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

-- λ> map product it

-- [1,4,3,12]

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

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

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

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

pares4 a b

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

| otherwise =

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

| x <- map (product . concat)

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

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

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

-- λ> pares3 2 60

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

-- λ> pares4 2 60

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

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

-- True

-- Comparacioń de eficiencia

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

-- 64

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

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

-- 64

-- (0.00 secs, 177,704 bytes)

-- Propiedades de equivalencia de pares

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

prop_pares :: Integer -> Integer -> Property

prop_pares a b =

a > 0 && b > 0 ==>

all (== pares1 a b)

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

, pares3

, pares4

]]

prop_pares2 :: Integer -> Integer -> Property

prop_pares2 a b =

a > 0 && b > 0 ==>

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

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

, pares3

, pares4

]]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

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

-- +++ OK, passed 100 tests.

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

-- +++ OK, passed 100 tests.

-- 1ª definición de nPares

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

nPares1 :: Integer -> Integer -> Integer

nPares1 a b = genericLength (pares4 a b)

-- 2ª definición de nPares

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

nPares2 :: Integer -> Integer -> Integer

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

-- Comparación de eficiencia

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

-- 8388608

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

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

-- 8388608

-- (0.00 secs, 234,688 bytes)

-- Propiedad de equivalencia de nPares

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

prop_nPares :: Integer -> Integer -> Property

prop_nPares a b =

a > 0 && b > 0 ==>

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

-- La comprobación es

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

-- +++ OK, passed 100 tests.

Pensamiento

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

Primos o cuadrados de primos

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

Definir la constante

   primosOcuadradosDePrimos :: [Integer]

1	primosOcuadradosDePrimos :: [Integer]

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

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

λ> take 20 primosOcuadradosDePrimos

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

λ> primosOcuadradosDePrimos !! (10^6)

15476729

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

Soluciones

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

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

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

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

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

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

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

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

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

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

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

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

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

import Test.QuickCheck

import Data.Numbers.Primes (primeFactors, primes)

-- 1ª solución

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

primosOcuadradosDePrimos :: [Integer]

primosOcuadradosDePrimos =

filter esPrimoOcuadradoDePrimo [2..]

esPrimoOcuadradoDePrimo :: Integer -> Bool

esPrimoOcuadradoDePrimo n = aux xs

where xs = primeFactors n

aux [_] = True

aux [x,y] = x == y

aux _ = False

-- 2ª solución

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

primosOcuadradosDePrimos2 :: [Integer]

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

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

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

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

-- Comparación de eficiencia

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

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

-- 223589

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

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

-- 223589

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

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

-- 7362497

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

-- Propiedad de equivalencia

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

-- La propiedad es

prop_equivalencia :: Int -> Property

prop_equivalencia n =

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

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

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

-- enteros distintos mayores que uno. Por ejemplo,

-- λ> take 20 unifactorizables

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

-- λ> unifactorizables !! 300

-- 1873

unifactorizables :: [Integer]

unifactorizables =

[n | n <- [2..]

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

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

-- lista ordenada estrictamente creciente xs (cuyos elementos son

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

-- mayor que 1). Por ejemplo,

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

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

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

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

-- λ> sublistasConProducto 2 [4,7]

-- []

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

sublistasConProducto _ [] = []

sublistasConProducto n (x:xs)

| x > n = []

| x == n = [[x]]

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

++ sublistasConProducto n xs

| otherwise = sublistasConProducto n xs

where (q,r) = quotRem n x

-- La comprobación es

-- λ> quickCheck prop_equivalencia

-- +++ OK, passed 100 tests.

Pensamiento

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

Antonio Machado

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

Mayor producto de n dígitos consecutivos de un número

PorJosé A. Alonso 29-noviembre-20196-diciembre-2019

Definir la función

   mayorProducto :: Int -> Integer -> Integer

1	mayorProducto :: Int -> Integer -> Integer

tal que (mayorProducto n x) es el mayor producto de n dígitos consecutivos del número x (suponiendo que x tiene al menos n dígitos). Por ejemplo,

   mayorProducto 2 325                  ==  10
   mayorProducto 5 11111                ==  1
   mayorProducto 5 113111               ==  3
   mayorProducto 5 110111               ==  0
   mayorProducto 5 10151112             ==  10
   mayorProducto 5 101511124            ==  10
   mayorProducto 5 (product [1..1000])  ==  41472

mayorProducto 2 325 == 10

mayorProducto 5 11111 == 1

mayorProducto 5 113111 == 3

mayorProducto 5 110111 == 0

mayorProducto 5 10151112 == 10

mayorProducto 5 101511124 == 10

mayorProducto 5 (product [1..1000]) == 41472

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

Soluciones

import Data.List (inits, tails)
import Data.Char (digitToInt)

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

mayorProducto :: Int -> Integer -> Integer
mayorProducto n x =
  maximum [product xs | xs <- segmentos n (digitos x)]

-- (digitos x) es la lista de las digitos del número x. Por ejemplo, 
--    digitos 325  ==  [3,2,5]
digitos :: Integer -> [Integer]
digitos x = map (toInteger . digitToInt) (show x)

-- (segmentos n xs) es la lista de los segmentos de longitud n de la
-- lista xs. Por ejemplo,
--    segmentos 2 [3,5,4,6]  ==  [[3,5],[5,4],[4,6]]
segmentos :: Int -> [Integer] -> [[Integer]]
segmentos n xs = take (length xs - n + 1) (map (take n) (tails xs))

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

mayorProducto2 :: Int -> Integer -> Integer
mayorProducto2 n x = maximum (aux ns)
    where ns     = [read [d] | d <- show x]
          aux xs | length xs < n = []
                 | otherwise     = product (take n xs) : aux (tail xs)

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

mayorProducto3 :: Int -> Integer -> Integer
mayorProducto3 n = maximum
                 . map (product . take n)
                 . filter ((>=n) . length) 
                 . tails
                 . digitos

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

mayorProducto4 :: Int -> Integer -> Integer
mayorProducto4 n = maximum  
                 . map (product . map (fromIntegral . digitToInt)) 
                 . filter ((==n) . length) 
                 . concatMap inits
                 . tails 
                 . show

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

--    λ> mayorProducto 5 (product [1..500])
--    28224
--    (0.01 secs, 1,645,256 bytes)
--    λ> mayorProducto2 5 (product [1..500])
--    28224
--    (0.03 secs, 5,848,416 bytes)
--    λ> mayorProducto3 5 (product [1..500])
--    28224
--    (0.03 secs, 1,510,640 bytes)
--    λ> mayorProducto4 5 (product [1..500])
--    28224
--    (1.85 secs, 10,932,551,216 bytes)
--    
--    λ> mayorProducto 5 (product [1..7000])
--    46656
--    (0.10 secs, 68,590,808 bytes)
--    λ> mayorProducto2 5 (product [1..7000])
--    46656
--    (1.63 secs, 157,031,432 bytes)
--    λ> mayorProducto3 5 (product [1..7000])
--    46656
--    (1.55 secs, 65,727,176 bytes)

import Data.List (inits, tails)

import Data.Char (digitToInt)

-- 1ª solución

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

mayorProducto :: Int -> Integer -> Integer

mayorProducto n x =

maximum [product xs | xs <- segmentos n (digitos x)]

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

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

digitos :: Integer -> [Integer]

digitos x = map (toInteger . digitToInt) (show x)

-- (segmentos n xs) es la lista de los segmentos de longitud n de la

-- lista xs. Por ejemplo,

-- segmentos 2 [3,5,4,6] == [[3,5],[5,4],[4,6]]

segmentos :: Int -> [Integer] -> [[Integer]]

segmentos n xs = take (length xs - n + 1) (map (take n) (tails xs))

-- 2ª solución

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

mayorProducto2 :: Int -> Integer -> Integer

mayorProducto2 n x = maximum (aux ns)

where ns = [read [d] | d <- show x]

aux xs | length xs < n = []

| otherwise = product (take n xs) : aux (tail xs)

-- 3ª solución

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

mayorProducto3 :: Int -> Integer -> Integer

mayorProducto3 n = maximum

. map (product . take n)

. filter ((>=n) . length)

. tails

. digitos

-- 4ª solución

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

mayorProducto4 :: Int -> Integer -> Integer

mayorProducto4 n = maximum

. map (product . map (fromIntegral . digitToInt))

. filter ((==n) . length)

. concatMap inits

. tails

. show

-- Comparación de eficiencia

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

-- λ> mayorProducto 5 (product [1..500])

-- 28224

-- (0.01 secs, 1,645,256 bytes)

-- λ> mayorProducto2 5 (product [1..500])

-- 28224

-- (0.03 secs, 5,848,416 bytes)

-- λ> mayorProducto3 5 (product [1..500])

-- 28224

-- (0.03 secs, 1,510,640 bytes)

-- λ> mayorProducto4 5 (product [1..500])

-- 28224

-- (1.85 secs, 10,932,551,216 bytes)

-- λ> mayorProducto 5 (product [1..7000])

-- 46656

-- (0.10 secs, 68,590,808 bytes)

-- λ> mayorProducto2 5 (product [1..7000])

-- 46656

-- (1.63 secs, 157,031,432 bytes)

-- λ> mayorProducto3 5 (product [1..7000])

-- 46656

-- (1.55 secs, 65,727,176 bytes)

Pensamiento

«El control de la complejidad es la esencia de la programación.» ~ B.W. Kernigan

Menor número triangular con más de n divisores

PorJosé A. Alonso 27-noviembre-20194-diciembre-2019

La sucesión de los números triangulares se obtiene sumando los números naturales.

   *     *      *        *         *   
        * *    * *      * *       * *  
              * * *    * * *     * * * 
                      * * * *   * * * *
                               * * * * * 
   1     3      6        10        15

* * * * *

* * * * * * * *

* * * * * * * * *

* * * * * * * *

* * * * *

1 3 6 10 15

Así, el 7º número triangular es

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

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

Los primeros 10 números triangulares son

   1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

1	1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...

Los divisores de los primeros 7 números triangulares son:

    1: 1
    3: 1,3
    6: 1,2,3,6
   10: 1,2,5,10
   15: 1,3,5,15
   21: 1,3,7,21
   28: 1,2,4,7,14,28

1: 1

3: 1,3

6: 1,2,3,6

10: 1,2,5,10

15: 1,3,5,15

21: 1,3,7,21

28: 1,2,4,7,14,28

Como se puede observar, 28 es el menor número triangular con más de 5 divisores.

Definir la función

   menorTriangularConAlMenosNDivisores :: Int -> Integer

1	menorTriangularConAlMenosNDivisores :: Int -> Integer

tal que (menorTriangularConAlMenosNDivisores n) es el menor número triangular que tiene al menos n divisores. Por ejemplo,

   menorTriangularConAlMenosNDivisores 5    ==  28
   menorTriangularConAlMenosNDivisores 50   ==  25200
   menorTriangularConAlMenosNDivisores 500  ==  76576500

menorTriangularConAlMenosNDivisores 5 == 28

menorTriangularConAlMenosNDivisores 50 == 25200

menorTriangularConAlMenosNDivisores 500 == 76576500

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

Soluciones

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

menorTriangularConAlMenosNDivisores :: Int -> Integer
menorTriangularConAlMenosNDivisores n = 
  head [x | x <- triangulares, nDivisores x >= n]

-- Nota: Se usarán las funciones
-- + triangulares definida en [Números triangulares](http://bit.ly/2rtr6a3) y
-- + nDivisores definida en [Número de divisores](http://bit.ly/2DgVh74)

-- triangulares es la sucesión de los números triangulares. Por ejemplo,
--    take 10 triangulares  ==  [1,3,6,10,15,21,28,36,45,55]
triangulares :: [Integer]
triangulares = scanl1 (+) [1..]

-- (nDivisores x) es el número de divisores de x. Por ejemplo,
--    nDivisores 28  ==  6
nDivisores :: Integer -> Int
nDivisores = product . map ((+1) . length) . group . primeFactors

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

menorTriangularConAlMenosNDivisores :: Int -> Integer

menorTriangularConAlMenosNDivisores n =

head [x | x <- triangulares, nDivisores x >= n]

-- Nota: Se usarán las funciones

-- + triangulares definida en [Números triangulares](http://bit.ly/2rtr6a3) y

-- + nDivisores definida en [Número de divisores](http://bit.ly/2DgVh74)

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

-- take 10 triangulares == [1,3,6,10,15,21,28,36,45,55]

triangulares :: [Integer]

triangulares = scanl1 (+) [1..]

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

-- nDivisores 28 == 6

nDivisores :: Integer -> Int

nDivisores = product . map ((+1) . length) . group . primeFactors

Pensamiento

«La Matemática es una ciencia experimental y la computación es el experimento.» ~ Rivin

Números triangulares

PorJosé A. Alonso 25-noviembre-201925-marzo-2022

La sucesión de los números triangulares se obtiene sumando los números naturales.

   *     *      *        *         *   
        * *    * *      * *       * *  
              * * *    * * *     * * * 
                      * * * *   * * * *
                               * * * * * 
   1     3      6        10        15

* * * * *

* * * * * * * *

* * * * * * * * *

* * * * * * * *

* * * * *

1 3 6 10 15

Así, los 5 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

1 = 1

3 = 1+2

6 = 1+2+3

10 = 1+2+3+4

15 = 1+2+3+4+5

Definir la función

   triangulares :: [Integer]

1	triangulares :: [Integer]

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

   take 10 triangulares  ==  [1,3,6,10,15,21,28,36,45,55]
   maximum (take (5*10^6) triangulares4)  ==  12500002500000

1 2	take 10 triangulares == [1,3,6,10,15,21,28,36,45,55] maximum (take (5*10^6) triangulares4) == 12500002500000

Comprobar con QuickCheck que entre dos números triangulares consecutivos siempre hay un número primo.

Soluciones

import Test.QuickCheck (Property, (==>), quickCheck)
import Data.Numbers.Primes (primes)

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

triangulares :: [Integer]
triangulares = [sum [1..n] | n <- [1..]]

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

triangulares2 :: [Integer]
triangulares2 = map triangular [1..]

-- (triangular n) es el n-ésimo número triangular. Por ejemplo, 
--    triangular 5  ==  15
triangular :: Integer -> Integer
triangular 1 = 1
triangular n = n + triangular (n-1)

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

triangulares3 :: [Integer]
triangulares3 = 1 : [x+y | (x,y) <- zip [2..] triangulares]

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

triangulares4 :: [Integer]
triangulares4 = scanl1 (+) [1..]

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

triangulares5 :: [Integer]
triangulares5 = [(n*(n+1)) `div` 2 | n <- [1..]]

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

--    λ> maximum (take (10^4) triangulares)
--    50005000
--    (2.10 secs, 8,057,774,104 bytes)
--    λ> maximum (take (10^4) triangulares2)
--    50005000
--    (18.89 secs, 12,142,690,784 bytes)
--    λ> maximum (take (10^4) triangulares3)
--    50005000
--    (0.01 secs, 4,600,976 bytes)
--    λ> maximum (take (10^4) triangulares4)
--    50005000
--    (0.01 secs, 3,643,192 bytes)
--    λ> maximum (take (10^4) triangulares5)
--    50005000
--    (0.02 secs, 5,161,464 bytes)
--    
--    λ> maximum (take (3*10^4) triangulares3)
--    450015000
--    (26.06 secs, 72,546,027,136 bytes)
--    λ> maximum (take (3*10^4) triangulares4)
--    450015000
--    (0.02 secs, 10,711,600 bytes)
--    λ> maximum (take (3*10^4) triangulares5)
--    450015000
--    (0.03 secs, 15,272,320 bytes)
--    
--    λ> maximum (take (5*10^6) triangulares4)
--    12500002500000
--    (1.67 secs, 1,772,410,336 bytes)
--    λ> maximum (take (5*10^6) triangulares5)
--    12500002500000
--    (4.09 secs, 2,532,407,720 bytes)

-- La propiedad es
prop_triangulares :: Int -> Property
prop_triangulares n =
  n >= 0 ==> siguientePrimo x < y
  where (x:y:_) = drop n triangulares4

-- (siguientePrimo n) es el menor primo mayor o igual que n. Por
-- ejemplo, 
--    siguientePrimo 14  ==  17
--    siguientePrimo 17  ==  17
siguientePrimo :: Integer -> Integer
siguientePrimo n = head (dropWhile (< n) primes)

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

import Test.QuickCheck (Property, (==>), quickCheck)

import Data.Numbers.Primes (primes)

-- 1ª solución

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

triangulares :: [Integer]

triangulares = [sum [1..n] | n <- [1..]]

-- 2ª solución

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

triangulares2 :: [Integer]

triangulares2 = map triangular [1..]

-- (triangular n) es el n-ésimo número triangular. Por ejemplo,

-- triangular 5 == 15

triangular :: Integer -> Integer

triangular 1 = 1

triangular n = n + triangular (n-1)

-- 3ª solución

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

triangulares3 :: [Integer]

triangulares3 = 1 : [x+y | (x,y) <- zip [2..] triangulares]

-- 4ª solución

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

triangulares4 :: [Integer]

triangulares4 = scanl1 (+) [1..]

-- 5ª solución

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

triangulares5 :: [Integer]

triangulares5 = [(n*(n+1)) `div` 2 | n <- [1..]]

-- Comparación de eficiencia

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

-- λ> maximum (take (10^4) triangulares)

-- 50005000

-- (2.10 secs, 8,057,774,104 bytes)

-- λ> maximum (take (10^4) triangulares2)

-- 50005000

-- (18.89 secs, 12,142,690,784 bytes)

-- λ> maximum (take (10^4) triangulares3)

-- 50005000

-- (0.01 secs, 4,600,976 bytes)

-- λ> maximum (take (10^4) triangulares4)

-- 50005000

-- (0.01 secs, 3,643,192 bytes)

-- λ> maximum (take (10^4) triangulares5)

-- 50005000

-- (0.02 secs, 5,161,464 bytes)

-- λ> maximum (take (3*10^4) triangulares3)

-- 450015000

-- (26.06 secs, 72,546,027,136 bytes)

-- λ> maximum (take (3*10^4) triangulares4)

-- 450015000

-- (0.02 secs, 10,711,600 bytes)

-- λ> maximum (take (3*10^4) triangulares5)

-- 450015000

-- (0.03 secs, 15,272,320 bytes)

-- λ> maximum (take (5*10^6) triangulares4)

-- 12500002500000

-- (1.67 secs, 1,772,410,336 bytes)

-- λ> maximum (take (5*10^6) triangulares5)

-- 12500002500000

-- (4.09 secs, 2,532,407,720 bytes)

-- La propiedad es

prop_triangulares :: Int -> Property

prop_triangulares n =

n >= 0 ==> siguientePrimo x < y

where (x:y:_) = drop n triangulares4

-- (siguientePrimo n) es el menor primo mayor o igual que n. Por

-- ejemplo,

-- siguientePrimo 14 == 17

-- siguientePrimo 17 == 17

siguientePrimo :: Integer -> Integer

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

-- La comprobación es

-- λ> quickCheck prop_triangulares

-- +++ OK, passed 100 tests.

Pensamiento

Autores, la escena acaba
con un dogma de teatro:
En el principio era la máscara.

Antonio Machado

Las sucesiones de Loomis

PorJosé A. Alonso 15-mayo-201925-abril-2021

La sucesión de Loomis generada por un número entero positivo x es la sucesión cuyos términos se definen por

f(0) es x
f(n) es la suma de f(n-1) y el producto de los dígitos no nulos de f(n-1)

Los primeros términos de las primeras sucesiones de Loomis son

Generada por 1: 1, 2, 4, 8, 16, 22, 26, 38, 62, 74, 102, 104, 108, 116, 122, …
Generada por 2: 2, 4, 8, 16, 22, 26, 38, 62, 74, 102, 104, 108, 116, 122, 126, …
Generada por 3: 3, 6, 12, 14, 18, 26, 38, 62, 74, 102, 104, 108, 116, 122, 126, …
Generada por 4: 4, 8, 16, 22, 26, 38, 62, 74, 102, 104, 108, 116, 122, 126, 138, …
Generada por 5: 5, 10, 11, 12, 14, 18, 26, 38, 62, 74, 102, 104, 108, 116, 122, …

Se observa que a partir de un término todas coinciden con la generada por 1. Dicho término se llama el punto de convergencia. Por ejemplo,

la generada por 2 converge a 2
la generada por 3 converge a 26
la generada por 4 converge a 4
la generada por 5 converge a 26

Definir las siguientes funciones

   sucLoomis           :: Integer -> [Integer]
   convergencia        :: Integer -> Integer
   graficaConvergencia :: [Integer] -> IO ()

sucLoomis :: Integer -> [Integer]

convergencia :: Integer -> Integer

graficaConvergencia :: [Integer] -> IO ()

tales que

(sucLoomis x) es la sucesión de Loomis generada por x. Por ejemplo,

     λ> take 15 (sucLoomis 1)
     [1,2,4,8,16,22,26,38,62,74,102,104,108,116,122]
     λ> take 15 (sucLoomis 2)
     [2,4,8,16,22,26,38,62,74,102,104,108,116,122,126]
     λ> take 15 (sucLoomis 3)
     [3,6,12,14,18,26,38,62,74,102,104,108,116,122,126]
     λ> take 15 (sucLoomis 4)
     [4,8,16,22,26,38,62,74,102,104,108,116,122,126,138]
     λ> take 15 (sucLoomis 5)
     [5,10,11,12,14,18,26,38,62,74,102,104,108,116,122]
     λ> take 15 (sucLoomis 20)
     [20,22,26,38,62,74,102,104,108,116,122,126,138,162,174]
     λ> take 15 (sucLoomis 100)
     [100,101,102,104,108,116,122,126,138,162,174,202,206,218,234]
     λ> sucLoomis 1 !! (2*10^5)
     235180736652

λ> take 15 (sucLoomis 1)

[1,2,4,8,16,22,26,38,62,74,102,104,108,116,122]

λ> take 15 (sucLoomis 2)

[2,4,8,16,22,26,38,62,74,102,104,108,116,122,126]

λ> take 15 (sucLoomis 3)

[3,6,12,14,18,26,38,62,74,102,104,108,116,122,126]

λ> take 15 (sucLoomis 4)

[4,8,16,22,26,38,62,74,102,104,108,116,122,126,138]

λ> take 15 (sucLoomis 5)

[5,10,11,12,14,18,26,38,62,74,102,104,108,116,122]

λ> take 15 (sucLoomis 20)

[20,22,26,38,62,74,102,104,108,116,122,126,138,162,174]

λ> take 15 (sucLoomis 100)

[100,101,102,104,108,116,122,126,138,162,174,202,206,218,234]

λ> sucLoomis 1 !! (2*10^5)

235180736652

(convergencia x) es el término de convergencia de la sucesioń de Loomis generada por x xon la geerada por 1. Por ejemplo,

     convergencia  2      ==  2
     convergencia  3      ==  26
     convergencia  4      ==  4
     convergencia 17      ==  38
     convergencia 19      ==  102
     convergencia 43      ==  162
     convergencia 27      ==  202
     convergencia 58      ==  474
     convergencia 63      ==  150056
     convergencia 81      ==  150056
     convergencia 89      ==  150056
     convergencia (10^12) ==  1000101125092

convergencia 2 == 2

convergencia 3 == 26

convergencia 4 == 4

convergencia 17 == 38

convergencia 19 == 102

convergencia 43 == 162

convergencia 27 == 202

convergencia 58 == 474

convergencia 63 == 150056

convergencia 81 == 150056

convergencia 89 == 150056

convergencia (10^12) == 1000101125092

(graficaConvergencia xs) dibuja la gráfica de los términos de convergencia de las sucesiones de Loomis generadas por los elementos de xs. Por ejemplo, (graficaConvergencia ([1..50]) dibuja

y graficaConvergencia ([1..148] \ [63,81,89,137]) dibuja

Soluciones

import Data.List               ((\\))
import Data.Char               (digitToInt)
import Graphics.Gnuplot.Simple (plotList, Attribute (Key, Title, XRange, PNG))

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

sucLoomis :: Integer -> [Integer]
sucLoomis x = map (loomis x) [0..]

loomis :: Integer -> Integer -> Integer
loomis x 0 = x
loomis x n = y + productoDigitosNoNulos y
  where y = loomis x (n-1)

productoDigitosNoNulos :: Integer -> Integer
productoDigitosNoNulos = product . digitosNoNulos

digitosNoNulos :: Integer -> [Integer]
digitosNoNulos x =
  [read [c] | c <- show x, c /= '0']

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

sucLoomis2 :: Integer -> [Integer]
sucLoomis2 = iterate siguienteLoomis 

siguienteLoomis :: Integer -> Integer
siguienteLoomis y = y + productoDigitosNoNulos y

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

sucLoomis3 :: Integer -> [Integer]
sucLoomis3 =
  iterate ((+) <*> product .
           map (toInteger . digitToInt) .
           filter (/= '0') . show)

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

--    λ> sucLoomis 1 !! 30000
--    6571272766
--    (2.45 secs, 987,955,944 bytes)
--    λ> sucLoomis2 1 !! 30000
--    6571272766
--    (2.26 secs, 979,543,328 bytes)
--    λ> sucLoomis3 1 !! 30000
--    6571272766
--    (0.31 secs, 88,323,832 bytes)

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

convergencia1 :: Integer -> Integer
convergencia1 x =
  head (dropWhile noEnSucLoomisDe1 (sucLoomis x))

noEnSucLoomisDe1 :: Integer -> Bool
noEnSucLoomisDe1 x = not (pertenece x sucLoomisDe1)

sucLoomisDe1 :: [Integer]
sucLoomisDe1 = sucLoomis 1

pertenece :: Integer -> [Integer] -> Bool
pertenece x ys =
  x == head (dropWhile (<x) ys)

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

convergencia2 :: Integer -> Integer
convergencia2 = aux (sucLoomis3 1) . sucLoomis3
 where aux as@(x:xs) bs@(y:ys) | x == y    = x
                               | x < y     = aux xs bs
                               | otherwise = aux as ys

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

convergencia3 :: Integer -> Integer
convergencia3 = head . interseccion (sucLoomis3 1) . sucLoomis3
 
-- (interseccion xs ys) es la intersección entre las listas ordenadas xs
-- e ys. Por ejemplo,
--    λ> take 10 (interseccion (sucLoomis3 1) (sucLoomis3 2))
--    [2,4,8,16,22,26,38,62,74,102]
interseccion :: Ord a => [a] -> [a] -> [a]
interseccion = aux
  where aux as@(x:xs) bs@(y:ys) = case compare x y of
                                    LT ->     aux xs bs
                                    EQ -> x : aux xs ys
                                    GT ->     aux as ys
        aux _         _         = []                           

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

convergencia4 :: Integer -> Integer
convergencia4 x = perteneceA (sucLoomis3 x) 1
  where perteneceA (y:ys) n | y == c    = y
                            | otherwise = perteneceA ys c
          where c = head $ dropWhile (< y) $ sucLoomis3 n

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

--    λ> convergencia1 (10^4)
--    150056
--    (2.94 secs, 1,260,809,808 bytes)
--    λ> convergencia2 (10^4)
--    150056
--    (0.03 secs, 700,240 bytes)
--    λ> convergencia3 (10^4)
--    150056
--    (0.03 secs, 1,165,496 bytes)
--    λ> convergencia4 (10^4)
--    150056
--    (0.02 secs, 1,119,648 bytes)
--    
--    λ> convergencia2 (10^12)
--    1000101125092
--    (1.81 secs, 714,901,080 bytes)
--    λ> convergencia3 (10^12)
--    1000101125092
--    (1.92 secs, 744,932,184 bytes)
--    λ> convergencia4 (10^12)
--    1000101125092
--    (1.82 secs, 941,053,328 bytes)

-- Definición de graficaConvergencia
-- ==================================

graficaConvergencia :: [Integer] -> IO ()
graficaConvergencia xs =
  plotList [ Key Nothing
           , Title "Convergencia de sucesiones de Loomis"
           , XRange (fromIntegral (minimum xs),fromIntegral (maximum xs))
           , PNG "Las_sucesiones_de_Loomis_2.png"
           ]
           [(x,convergencia2 x) | x <- xs]

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

import Data.List ((\\))

import Data.Char (digitToInt)

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

-- 1ª definición de sucLoomis

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

sucLoomis :: Integer -> [Integer]

sucLoomis x = map (loomis x) [0..]

loomis :: Integer -> Integer -> Integer

loomis x 0 = x

loomis x n = y + productoDigitosNoNulos y

where y = loomis x (n-1)

productoDigitosNoNulos :: Integer -> Integer

productoDigitosNoNulos = product . digitosNoNulos

digitosNoNulos :: Integer -> [Integer]

digitosNoNulos x =

[read [c] | c <- show x, c /= '0']

-- 2ª definición de sucLoomis

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

sucLoomis2 :: Integer -> [Integer]

sucLoomis2 = iterate siguienteLoomis

siguienteLoomis :: Integer -> Integer

siguienteLoomis y = y + productoDigitosNoNulos y

-- 3ª definición de sucLoomis

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

sucLoomis3 :: Integer -> [Integer]

sucLoomis3 =

iterate ((+) <*> product .

map (toInteger . digitToInt) .

filter (/= '0') . show)

-- Comparación de eficiencia

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

-- λ> sucLoomis 1 !! 30000

-- 6571272766

-- (2.45 secs, 987,955,944 bytes)

-- λ> sucLoomis2 1 !! 30000

-- 6571272766

-- (2.26 secs, 979,543,328 bytes)

-- λ> sucLoomis3 1 !! 30000

-- 6571272766

-- (0.31 secs, 88,323,832 bytes)

-- 1ª definición de convergencia

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

convergencia1 :: Integer -> Integer

convergencia1 x =

head (dropWhile noEnSucLoomisDe1 (sucLoomis x))

noEnSucLoomisDe1 :: Integer -> Bool

noEnSucLoomisDe1 x = not (pertenece x sucLoomisDe1)

sucLoomisDe1 :: [Integer]

sucLoomisDe1 = sucLoomis 1

pertenece :: Integer -> [Integer] -> Bool

pertenece x ys =

x == head (dropWhile (<x) ys)

-- 2ª definición de convergencia

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

convergencia2 :: Integer -> Integer

convergencia2 = aux (sucLoomis3 1) . sucLoomis3

where aux as@(x:xs) bs@(y:ys) | x == y = x

| x < y = aux xs bs

| otherwise = aux as ys

-- 3ª definición de convergencia

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

convergencia3 :: Integer -> Integer

convergencia3 = head . interseccion (sucLoomis3 1) . sucLoomis3

-- (interseccion xs ys) es la intersección entre las listas ordenadas xs

-- e ys. Por ejemplo,

-- λ> take 10 (interseccion (sucLoomis3 1) (sucLoomis3 2))

-- [2,4,8,16,22,26,38,62,74,102]

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

interseccion = aux

where aux as@(x:xs) bs@(y:ys) = case compare x y of

LT -> aux xs bs

EQ -> x : aux xs ys

GT -> aux as ys

aux _ _ = []

-- 4ª definición de convergencia

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

convergencia4 :: Integer -> Integer

convergencia4 x = perteneceA (sucLoomis3 x) 1

where perteneceA (y:ys) n | y == c = y

| otherwise = perteneceA ys c

where c = head $ dropWhile (< y) $ sucLoomis3 n

-- Comparación de eficiencia

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

-- λ> convergencia1 (10^4)

-- 150056

-- (2.94 secs, 1,260,809,808 bytes)

-- λ> convergencia2 (10^4)

-- 150056

-- (0.03 secs, 700,240 bytes)

-- λ> convergencia3 (10^4)

-- 150056

-- (0.03 secs, 1,165,496 bytes)

-- λ> convergencia4 (10^4)

-- 150056

-- (0.02 secs, 1,119,648 bytes)

-- λ> convergencia2 (10^12)

-- 1000101125092

-- (1.81 secs, 714,901,080 bytes)

-- λ> convergencia3 (10^12)

-- 1000101125092

-- (1.92 secs, 744,932,184 bytes)

-- λ> convergencia4 (10^12)

-- 1000101125092

-- (1.82 secs, 941,053,328 bytes)

-- Definición de graficaConvergencia

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

graficaConvergencia :: [Integer] -> IO ()

graficaConvergencia xs =

plotList [ Key Nothing

, Title "Convergencia de sucesiones de Loomis"

, XRange (fromIntegral (minimum xs),fromIntegral (maximum xs))

, PNG "Las_sucesiones_de_Loomis_2.png"

]

[(x,convergencia2 x) | x <- xs]

Pensamiento

Era una noche del mes
de mayo, azul y serena.
Sobre el agudo ciprés
brillaba la luna llena.

Antonio Machado

Mayor producto de n dígitos consecutivos de un número

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

Definir la función

   mayorProducto :: Int -> Integer -> Integer

1	mayorProducto :: Int -> Integer -> Integer

tal que (mayorProducto n x) es el mayor producto de n dígitos consecutivos del número x (suponiendo que x tiene al menos n dígitos). Por ejemplo,

   mayorProducto 2 325                  ==  10
   mayorProducto 5 11111                ==  1
   mayorProducto 5 113111               ==  3
   mayorProducto 5 110111               ==  0
   mayorProducto 5 10151112             ==  10
   mayorProducto 5 101511124            ==  10
   mayorProducto 5 (product [1..1000])  ==  41472

mayorProducto 2 325 == 10

mayorProducto 5 11111 == 1

mayorProducto 5 113111 == 3

mayorProducto 5 110111 == 0

mayorProducto 5 10151112 == 10

mayorProducto 5 101511124 == 10

mayorProducto 5 (product [1..1000]) == 41472

Soluciones

import Test.QuickCheck
import Data.List (inits, tails)
import Data.Char (digitToInt)

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

mayorProducto1 :: Int -> Integer -> Integer
mayorProducto1 n x =
  maximum [product xs | xs <- segmentos n (cifras x)]

-- (cifras x) es la lista de las cifras del número x. Por ejemplo, 
--    cifras 325  ==  [3,2,5]
cifras :: Integer -> [Integer]
cifras x = map (toInteger . digitToInt) (show x)

-- (segmentos n xs) es la lista de los segmentos de longitud n de la
-- lista xs. Por ejemplo,
--    segmentos 2 [3,5,4,6]  ==  [[3,5],[5,4],[4,6]]
segmentos :: Int -> [Integer] -> [[Integer]]
segmentos n xs = take (length xs - n + 1) (map (take n) (tails xs))

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

mayorProducto2 :: Int -> Integer -> Integer
mayorProducto2 n x = maximum (aux ns)
    where ns     = [read [d] | d <- show x]
          aux xs | length xs < n = []
                 | otherwise     = product (take n xs) : aux (tail xs)

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

mayorProducto3 :: Int -> Integer -> Integer
mayorProducto3 n = maximum
                 . map (product . take n)
                 . filter ((>=n) . length) 
                 . tails
                 . cifras

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

mayorProducto4 :: Int -> Integer -> Integer
mayorProducto4 n = maximum  
                 . map (product . map (fromIntegral . digitToInt)) 
                 . filter ((==n) . length) 
                 . concatMap inits
                 . tails 
                 . show

-- ---------------------------------------------------------------------
-- Comparación de soluciones                                          --
-- ---------------------------------------------------------------------

-- Tiempo (en segundos) del cálculo de (mayorProducto 5 (product [1..n]))
-- 
--    | Def | 500  | 1000 | 5000 | 10000 
--    | 1   | 0.01 | 0.02 | 0.06 | 0.11
--    | 2   | 0.01 | 0.03 | 0.80 | 3.98
--    | 3   | 0.01 | 0.03 | 0.82 | 4.13
--    | 4   | 2.77 |      |      |

import Test.QuickCheck

import Data.List (inits, tails)

import Data.Char (digitToInt)

-- 1ª solución

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

mayorProducto1 :: Int -> Integer -> Integer

mayorProducto1 n x =

maximum [product xs | xs <- segmentos n (cifras x)]

-- (cifras x) es la lista de las cifras del número x. Por ejemplo,

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

cifras :: Integer -> [Integer]

cifras x = map (toInteger . digitToInt) (show x)

-- (segmentos n xs) es la lista de los segmentos de longitud n de la

-- lista xs. Por ejemplo,

-- segmentos 2 [3,5,4,6] == [[3,5],[5,4],[4,6]]

segmentos :: Int -> [Integer] -> [[Integer]]

segmentos n xs = take (length xs - n + 1) (map (take n) (tails xs))

-- 2ª solución

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

mayorProducto2 :: Int -> Integer -> Integer

mayorProducto2 n x = maximum (aux ns)

where ns = [read [d] | d <- show x]

aux xs | length xs < n = []

| otherwise = product (take n xs) : aux (tail xs)

-- 3ª solución

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

mayorProducto3 :: Int -> Integer -> Integer

mayorProducto3 n = maximum

. map (product . take n)

. filter ((>=n) . length)

. tails

. cifras

-- 4ª solución

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

mayorProducto4 :: Int -> Integer -> Integer

mayorProducto4 n = maximum

. map (product . map (fromIntegral . digitToInt))

. filter ((==n) . length)

. concatMap inits

. tails

. show

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

-- Comparación de soluciones --

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

-- Tiempo (en segundos) del cálculo de (mayorProducto 5 (product [1..n]))

-- | Def | 500 | 1000 | 5000 | 10000

-- | 1 | 0.01 | 0.02 | 0.06 | 0.11

-- | 2 | 0.01 | 0.03 | 0.80 | 3.98

-- | 3 | 0.01 | 0.03 | 0.82 | 4.13

-- | 4 | 2.77 | | |

Pensamiento

A las palabras de amor
les sienta bien su poquito
de exageración.

Antonio Machado

Operaciones con series de potencias

PorJosé A. Alonso 26-abril-201925-marzo-2022

Una serie de potencias es una serie de la forma

   a₀ + a₁x + a₂x² + a₃x³ + ...

1	a₀ + a₁x + a₂x² + a₃x³ + ...

Las series de potencias se pueden representar mediante listas infinitas. Por ejemplo, la serie de la función exponencial es

   e^x = 1 + x + x²/2! + x³/3! + ...

1	e^x = 1 + x + x²/2! + x³/3! + ...

y se puede representar por [1, 1, 1/2, 1/6, 1/24, 1/120, …]

Las operaciones con series se pueden ver como una generalización de las de los polinomios.

En lo que sigue, usaremos el tipo (Serie a) para representar las series de potencias con coeficientes en a y su definición es

   type Serie a = [a]

1	type Serie a = [a]

Definir las siguientes funciones

   opuesta      :: Num a => Serie a -> Serie a
   suma         :: Num a => Serie a -> Serie a -> Serie a
   resta        :: Num a => Serie a -> Serie a -> Serie a
   producto     :: Num a => Serie a -> Serie a -> Serie a
   cociente     :: Fractional a => Serie a -> Serie a -> Serie a
   derivada     :: (Num a, Enum a) => Serie a -> Serie a
   integral     :: (Fractional a, Enum a) => Serie a -> Serie a
   expx         :: Serie Rational

opuesta :: Num a => Serie a -> Serie a

suma :: Num a => Serie a -> Serie a -> Serie a

resta :: Num a => Serie a -> Serie a -> Serie a

producto :: Num a => Serie a -> Serie a -> Serie a

cociente :: Fractional a => Serie a -> Serie a -> Serie a

derivada :: (Num a, Enum a) => Serie a -> Serie a

integral :: (Fractional a, Enum a) => Serie a -> Serie a

expx :: Serie Rational

tales que

(opuesta xs) es la opuesta de la serie xs. Por ejemplo,

     λ> take 7 (opuesta [-6,-4..])
     [6,4,2,0,-2,-4,-6]

1 2	λ> take 7 (opuesta [-6,-4..]) [6,4,2,0,-2,-4,-6]

(suma xs ys) es la suma de las series xs e ys. Por ejemplo,

     λ> take 7 (suma [1,3..] [2,4..])
     [3,7,11,15,19,23,27]

1 2	λ> take 7 (suma [1,3..] [2,4..]) [3,7,11,15,19,23,27]

(resta xs ys) es la resta de las series xs es ys. Por ejemplo,

     λ> take 7 (resta [3,5..] [2,4..])
     [1,1,1,1,1,1,1]
     λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])
     [2,4,6,8,10,12,14]

λ> take 7 (resta [3,5..] [2,4..])

[1,1,1,1,1,1,1]

λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])

[2,4,6,8,10,12,14]

(producto xs ys) es el producto de las series xs e ys. Por ejemplo,

     λ> take 7 (producto [3,5..] [2,4..])
     [6,22,52,100,170,266,392]

1 2	λ> take 7 (producto [3,5..] [2,4..]) [6,22,52,100,170,266,392]

(cociente xs ys) es el cociente de las series xs e ys. Por ejemplo,

     λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..])
     [2.0,4.0,6.0,8.0,10.0,12.0,14.0]

1 2	λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..]) [2.0,4.0,6.0,8.0,10.0,12.0,14.0]

(derivada xs) es la derivada de la serie xs. Por ejemplo,

     λ> take 7 (derivada [2,4..])
     [4,12,24,40,60,84,112]

1 2	λ> take 7 (derivada [2,4..]) [4,12,24,40,60,84,112]

(integral xs) es la integral de la serie xs. Por ejemplo,

     λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0))
     [0.0,4.0,6.0,8.0,10.0,12.0,14.0]

1 2	λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0)) [0.0,4.0,6.0,8.0,10.0,12.0,14.0]

expx es la serie de la función exponencial. Por ejemplo,

     λ> take 8 expx
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (derivada expx)
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (integral expx)
     [0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 expx

[1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 (derivada expx)

[1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 (integral expx)

[0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

Soluciones

type Serie a = [a] 

opuesta :: Num a => Serie a -> Serie a
opuesta = map negate

suma :: Num a => Serie a -> Serie a -> Serie a
suma = zipWith (+)

resta :: Num a => Serie a -> Serie a -> Serie a
resta xs ys = suma xs (opuesta ys)

producto :: Num a => Serie a -> Serie a -> Serie a
producto (x:xs) zs@(y:ys) = 
    x*y : suma (producto xs zs) (map (x*) ys)

cociente :: Fractional a => Serie a -> Serie a -> Serie a
cociente (x:xs) (y:ys) = zs 
    where zs = x/y : map (/y) (resta xs (producto zs ys))  

derivada :: (Num a, Enum a) => Serie a -> Serie a
derivada (_:xs) = zipWith (*) xs [1..]

integral :: (Fractional a, Enum a) => Serie a -> Serie a
integral xs = 0 : zipWith (/) xs [1..]

expx :: Serie Rational
expx = map (1/) (map fromIntegral factoriales)

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

type Serie a = [a]

opuesta :: Num a => Serie a -> Serie a

opuesta = map negate

suma :: Num a => Serie a -> Serie a -> Serie a

suma = zipWith (+)

resta :: Num a => Serie a -> Serie a -> Serie a

resta xs ys = suma xs (opuesta ys)

producto :: Num a => Serie a -> Serie a -> Serie a

producto (x:xs) zs@(y:ys) =

x*y : suma (producto xs zs) (map (x*) ys)

cociente :: Fractional a => Serie a -> Serie a -> Serie a

cociente (x:xs) (y:ys) = zs

where zs = x/y : map (/y) (resta xs (producto zs ys))

derivada :: (Num a, Enum a) => Serie a -> Serie a

derivada (_:xs) = zipWith (*) xs [1..]

integral :: (Fractional a, Enum a) => Serie a -> Serie a

integral xs = 0 : zipWith (/) xs [1..]

expx :: Serie Rational

expx = map (1/) (map fromIntegral factoriales)

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

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

factoriales :: [Integer]

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

Pensamiento

Ni mármol duro y eterno,
ni música ni pintura,
sino palabra en el tiempo.

Antonio Machado

Reparto de escaños por la ley d’Hont

PorJosé A. Alonso 23-abril-201925-abril-2021

El sistema D’Hondt es una fórmula creada por Victor d’Hondt, que permite obtener el número de cargos electos asignados a las candidaturas, en proporción a los votos conseguidos.

Tras el recuento de los votos, se calcula una serie de divisores para cada partido. La fórmula de los divisores es V/N, donde V representa el número total de votos recibidos por el partido, y N representa cada uno de los números enteros desde 1 hasta el número de cargos electos de la circunscripción objeto de escrutinio. Una vez realizadas las divisiones de los votos de cada partido por cada uno de los divisores desde 1 hasta N, la asignación de cargos electos se hace ordenando los cocientes de las divisiones de mayor a menor y asignando a cada uno un escaño hasta que éstos se agoten

Definir la función

   reparto :: Int -> [Int] -> [(Int,Int)]

1	reparto :: Int -> [Int] -> [(Int,Int)]

tal que (reparto n vs) es la lista de los pares formados por los números de los partidos y el número de escaño que les corresponden al repartir n escaños en función de la lista de sus votos. Por ejemplo,

   ghci> reparto 7 [340000,280000,160000,60000,15000]
   [(1,3),(2,3),(3,1)]
   ghci> reparto 21 [391000,311000,184000,73000,27000,12000,2000]
   [(1,9),(2,7),(3,4),(4,1)]

ghci> reparto 7 [340000,280000,160000,60000,15000]

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

ghci> reparto 21 [391000,311000,184000,73000,27000,12000,2000]

[(1,9),(2,7),(3,4),(4,1)]

es decir, en el primer ejemplo,

al 1º partido (que obtuvo 340000 votos) le corresponden 3 escaños,
al 2º partido (que obtuvo 280000 votos) le corresponden 3 escaños,
al 3º partido (que obtuvo 160000 votos) le corresponden 1 escaño.

Soluciones

import Data.List (sort, group)

-- Para los ejemplos que siguen, se usará la siguiente ditribución de
-- votos entre 5 partidos.
ejVotos :: [Int]
ejVotos = [340000,280000,160000,60000,15000]

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

reparto :: Int -> [Int] -> [(Int,Int)]
reparto n vs = 
  [(x,1 + length xs) | (x:xs) <- group (sort (repartoAux n vs))] 

-- (repartoAux n vs) es el número de los partidos, cuyos votos son vs, que
-- obtienen los n escaños. Por ejemplo,
--    ghci> repartoAux 7 ejVotos
--    [1,2,1,3,2,1,2]
repartoAux :: Int -> [Int] -> [Int]
repartoAux n vs = map snd (repartoAux' n vs)

-- (repartoAux' n vs) es la lista formada por los n restos mayores
-- correspondientes a la lista de votos vs. Por ejemplo,
--    ghci> repartoAux' 7 ejVotos
--    [(340000,1),(280000,2),(170000,1),(160000,3),(140000,2),(113333,1),
--     (93333,2)]
repartoAux' :: Int -> [Int] -> [(Int,Int)]
repartoAux' n vs = 
  take n (reverse (sort (concatMap (restos n) (votosPartidos vs))))

-- (votosPartidos vs) es la lista con los pares formados por los votos y
-- el número de cada partido. Por ejemplo, 
--    ghci> votosPartidos ejVotos
--    [(340000,1),(280000,2),(160000,3),(60000,4),(15000,5)]
votosPartidos :: [Int] -> [(Int,Int)]
votosPartidos vs = zip vs [1..]

-- (restos n (x,i)) es la lista obtenidas dividiendo n entre 1, 2,..., n.
-- Por ejemplo, 
--    ghci> restos 5 (340000,1)
--    [(340000,1),(170000,1),(113333,1),(85000,1),(68000,1)]
restos :: Int -> (Int,Int) -> [(Int,Int)]
restos n (x,i) = [(x `div` k,i) | k <- [1..n]]

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

reparto2 :: Int -> [Int] -> [(Int,Int)]
reparto2 n xs = 
  ( map (\x -> (head x, length x))  
  . group  
  . sort  
  . map snd  
  . take n  
  . reverse  
  . sort
  ) [(x `div` i, p) | (x,p) <- zip xs [1..], i <- [1..n]]

import Data.List (sort, group)

-- Para los ejemplos que siguen, se usará la siguiente ditribución de

-- votos entre 5 partidos.

ejVotos :: [Int]

ejVotos = [340000,280000,160000,60000,15000]

-- 1ª solución

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

reparto :: Int -> [Int] -> [(Int,Int)]

reparto n vs =

[(x,1 + length xs) | (x:xs) <- group (sort (repartoAux n vs))]

-- (repartoAux n vs) es el número de los partidos, cuyos votos son vs, que

-- obtienen los n escaños. Por ejemplo,

-- ghci> repartoAux 7 ejVotos

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

repartoAux :: Int -> [Int] -> [Int]

repartoAux n vs = map snd (repartoAux' n vs)

-- (repartoAux' n vs) es la lista formada por los n restos mayores

-- correspondientes a la lista de votos vs. Por ejemplo,

-- ghci> repartoAux' 7 ejVotos

-- [(340000,1),(280000,2),(170000,1),(160000,3),(140000,2),(113333,1),

-- (93333,2)]

repartoAux' :: Int -> [Int] -> [(Int,Int)]

repartoAux' n vs =

take n (reverse (sort (concatMap (restos n) (votosPartidos vs))))

-- (votosPartidos vs) es la lista con los pares formados por los votos y

-- el número de cada partido. Por ejemplo,

-- ghci> votosPartidos ejVotos

-- [(340000,1),(280000,2),(160000,3),(60000,4),(15000,5)]

votosPartidos :: [Int] -> [(Int,Int)]

votosPartidos vs = zip vs [1..]

-- (restos n (x,i)) es la lista obtenidas dividiendo n entre 1, 2,..., n.

-- Por ejemplo,

-- ghci> restos 5 (340000,1)

-- [(340000,1),(170000,1),(113333,1),(85000,1),(68000,1)]

restos :: Int -> (Int,Int) -> [(Int,Int)]

restos n (x,i) = [(x `div` k,i) | k <- [1..n]]

-- 2ª solución

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

reparto2 :: Int -> [Int] -> [(Int,Int)]

reparto2 n xs =

( map (\x -> (head x, length x))

. group

. sort

. map snd

. take n

. reverse

. sort

) [(x `div` i, p) | (x,p) <- zip xs [1..], i <- [1..n]]

Pensamiento

Sus cantares llevan
agua de remanso,
que parece quieta.
Y que no lo está;
mas no tiene prisa
por ir a la mar.

Antonio Machado

Matriz dodecafónica

PorJosé A. Alonso 10-abril-201925-abril-2019

Como se explica en Create a Twelve-Tone Melody With a Twelve-Tone Matrix una matriz dodecafónica es una matriz de 12 filas y 12 columnas construidas siguiendo los siguientes pasos:

Se escribe en la primera fila una permutación de los números del 1 al 12. Por ejemplo,

     (  3  1  9  5  4  6  8  7 12 10 11  2 )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )
     (                                     )

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( )

Escribir la primera columna de forma que, para todo i (entre 2 y 12), a(i,1) es el número entre 1 y 12 que verifica la siguiente condición

     (a(1,1) - a(i,1)) = (a(1,i) - a(1,1)) (módulo 12)

1	(a(1,1) - a(i,1)) = (a(1,i) - a(1,1)) (módulo 12)

Siguiendo con el ejemplo anterior, la matriz con la 1ª fila y la 1ª columna es

     (  3  1  9  5  4  6  8  7 12 10 11  2 )
     (  5                                  )
     (  9                                  )
     (  1                                  )
     (  2                                  )
     ( 12                                  )
     ( 10                                  )
     ( 11                                  )
     (  6                                  )
     (  8                                  )
     (  7                                  )
     (  4                                  )

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( 5 )

( 9 )

( 1 )

( 2 )

( 12 )

( 10 )

( 11 )

( 6 )

( 8 )

( 7 )

( 4 )

Escribir la segunda fila de forma que, para todo j (entre 2 y 12), a(j,2) es el número entre 1 y 12 que verifica la siguiente condición

     (a(2,j) - a(1,j)) = (a(2,1) - a(1,1)) (módulo 12)

1	(a(2,j) - a(1,j)) = (a(2,1) - a(1,1)) (módulo 12)

Siguiendo con el ejemplo anterior, la matriz con la 1ª fila, 1ª columna y 2ª fila es

     (  3  1  9  5  4  6  8  7 12 10 11  2 )
     (  5  3 11  7  6  8 10  9  2 12  1  4 )
     (  9                                  )
     (  1                                  )
     (  2                                  )
     ( 12                                  )
     ( 10                                  )
     ( 11                                  )
     (  6                                  )
     (  8                                  )
     (  7                                  )
     (  4                                  )

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( 5 3 11 7 6 8 10 9 2 12 1 4 )

( 9 )

( 1 )

( 2 )

( 12 )

( 10 )

( 11 )

( 6 )

( 8 )

( 7 )

( 4 )

Las restantes filas se completan como la 2ª; es decir, para todo i (entre 3 y 12) y todo j (entre 2 y 12), a(i,j) es el número entre 1 y 12 que verifica la siguiente relación.

     (a(i,j) - a(1,j)) = (a(i,1) - a(1,1)) (módulo 12)

1	(a(i,j) - a(1,j)) = (a(i,1) - a(1,1)) (módulo 12)

Siguiendo con el ejemplo anterior, la matriz dodecafónica es

     (  3  1  9  5  4  6  8  7 12 10 11  2 )
     (  5  3 11  7  6  8 10  9  2 12  1  4 )
     (  9  7  3 11 10 12  2  1  6  4  5  8 )
     (  1 11  7  3  2  4  6  5 10  8  9 12 )
     (  2 12  8  4  3  5  7  6 11  9 10  1 )
     ( 12 10  6  2  1  3  5  4  9  7  8 11 )
     ( 10  8  4 12 11  1  3  2  7  5  6  9 )
     ( 11  9  5  1 12  2  4  3  8  6  7 10 )
     (  6  4 12  8  7  9 11 10  3  1  2  5 )
     (  8  6  2 10  9 11  1 12  5  3  4  7 )
     (  7  5  1  9  8 10 12 11  4  2  3  6 )
     (  4  2 10  6  5  7  9  8  1 11 12  3 )

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( 5 3 11 7 6 8 10 9 2 12 1 4 )

( 9 7 3 11 10 12 2 1 6 4 5 8 )

( 1 11 7 3 2 4 6 5 10 8 9 12 )

( 2 12 8 4 3 5 7 6 11 9 10 1 )

( 12 10 6 2 1 3 5 4 9 7 8 11 )

( 10 8 4 12 11 1 3 2 7 5 6 9 )

( 11 9 5 1 12 2 4 3 8 6 7 10 )

( 6 4 12 8 7 9 11 10 3 1 2 5 )

( 8 6 2 10 9 11 1 12 5 3 4 7 )

( 7 5 1 9 8 10 12 11 4 2 3 6 )

( 4 2 10 6 5 7 9 8 1 11 12 3 )

Definir la función

   matrizDodecafonica :: [Int] -> Matrix Int

1	matrizDodecafonica :: [Int] -> Matrix Int

tal que (matrizDodecafonica xs) es la matriz dodecafónica cuya primera fila es xs (que se supone que es una permutación de los números del 1 al 12). Por ejemplo,

   λ> matrizDodecafonica [3,1,9,5,4,6,8,7,12,10,11,2]
   (  3  1  9  5  4  6  8  7 12 10 11  2 )
   (  5  3 11  7  6  8 10  9  2 12  1  4 )
   (  9  7  3 11 10 12  2  1  6  4  5  8 )
   (  1 11  7  3  2  4  6  5 10  8  9 12 )
   (  2 12  8  4  3  5  7  6 11  9 10  1 )
   ( 12 10  6  2  1  3  5  4  9  7  8 11 )
   ( 10  8  4 12 11  1  3  2  7  5  6  9 )
   ( 11  9  5  1 12  2  4  3  8  6  7 10 )
   (  6  4 12  8  7  9 11 10  3  1  2  5 )
   (  8  6  2 10  9 11  1 12  5  3  4  7 )
   (  7  5  1  9  8 10 12 11  4  2  3  6 )
   (  4  2 10  6  5  7  9  8  1 11 12  3 )

λ> matrizDodecafonica [3,1,9,5,4,6,8,7,12,10,11,2]

( 3 1 9 5 4 6 8 7 12 10 11 2 )

( 5 3 11 7 6 8 10 9 2 12 1 4 )

( 9 7 3 11 10 12 2 1 6 4 5 8 )

( 1 11 7 3 2 4 6 5 10 8 9 12 )

( 2 12 8 4 3 5 7 6 11 9 10 1 )

( 12 10 6 2 1 3 5 4 9 7 8 11 )

( 10 8 4 12 11 1 3 2 7 5 6 9 )

( 11 9 5 1 12 2 4 3 8 6 7 10 )

( 6 4 12 8 7 9 11 10 3 1 2 5 )

( 8 6 2 10 9 11 1 12 5 3 4 7 )

( 7 5 1 9 8 10 12 11 4 2 3 6 )

( 4 2 10 6 5 7 9 8 1 11 12 3 )

Comprobar con QuickCheck para toda matriz dodecafónica D se verifican las siguientes propiedades:

todas las filas de D son permutaciones de los números 1 a 12,
todos los elementos de la diagonal de D son iguales y
la suma de todos los elementos de D es 936.

Nota: Este ejercicio ha sido propuesto por Francisco J. Hidalgo.

Soluciones

import Data.List
import Test.QuickCheck
import Data.Matrix

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

matrizDodecafonica :: [Int] -> Matrix Int
matrizDodecafonica xs = matrix 12 12 f
  where f (1,j) = xs !! (j-1)
        f (i,1) = modulo12 (2 * f (1,1) - f (1,i)) 
        f (i,j) = modulo12 (f (1,j) + f (i,1) - f (1,1)) 
        modulo12 0  = 12
        modulo12 12 = 12
        modulo12 x  = x `mod` 12
  
-- 2ª solución
-- ===========

matrizDodecafonica2 :: [Int] -> Matrix Int
matrizDodecafonica2 xs = fromLists (secuencias xs)

secuencias :: [Int] -> [[Int]]
secuencias xs = [secuencia a xs | a <- inversa xs]

inversa :: [Int] -> [Int]
inversa xs = map conv (map (\x -> (-x) + 2* (abs a)) xs)
  where a = head xs
        
secuencia :: Int -> [Int] -> [Int]
secuencia n xs = [conv (a+(n-b)) | a <- xs] 
  where b = head xs

conv :: Int -> Int
conv n | n == 0 = 12
       | n < 0 = conv (n+12)
       | n > 11 = conv (mod n 12)
       | otherwise = n          

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

-- Las propiedades son
prop_dodecafonica :: Int -> Property
prop_dodecafonica n = 
  n >= 0 ==>
  all esPermutacion (toLists d)
  && all (== d!(1,1)) [d!(i,i) | i <- [2..12]]
  && sum d == 936
  where xss = permutations [1..12]
        k   = n `mod` product [1..12]
        d   = matrizDodecafonica (xss !! k)
        esPermutacion ys = sort ys == [1..12]

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

import Data.List

import Test.QuickCheck

import Data.Matrix

-- 1ª solución

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

matrizDodecafonica :: [Int] -> Matrix Int

matrizDodecafonica xs = matrix 12 12 f

where f (1,j) = xs !! (j-1)

f (i,1) = modulo12 (2 * f (1,1) - f (1,i))

f (i,j) = modulo12 (f (1,j) + f (i,1) - f (1,1))

modulo12 0 = 12

modulo12 12 = 12

modulo12 x = x `mod` 12

-- 2ª solución

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

matrizDodecafonica2 :: [Int] -> Matrix Int

matrizDodecafonica2 xs = fromLists (secuencias xs)

secuencias :: [Int] -> [[Int]]

secuencias xs = [secuencia a xs | a <- inversa xs]

inversa :: [Int] -> [Int]

inversa xs = map conv (map (\x -> (-x) + 2* (abs a)) xs)

where a = head xs

secuencia :: Int -> [Int] -> [Int]

secuencia n xs = [conv (a+(n-b)) | a <- xs]

where b = head xs

conv :: Int -> Int

conv n | n == 0 = 12

| n < 0 = conv (n+12)

| n > 11 = conv (mod n 12)

| otherwise = n

-- Propiedades

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

-- Las propiedades son

prop_dodecafonica :: Int -> Property

prop_dodecafonica n =

n >= 0 ==>

all esPermutacion (toLists d)

&& all (== d!(1,1)) [d!(i,i) | i <- [2..12]]

&& sum d == 936

where xss = permutations [1..12]

k = n `mod` product [1..12]

d = matrizDodecafonica (xss !! k)

esPermutacion ys = sort ys == [1..12]

-- La comprobación es

-- λ> quickCheck prop_dodecafonica

-- +++ OK, passed 100 tests.

Pensamiento

Como el olivar,
mucho fruto lleva,
poca sombra da.

Antonio Machado

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Referencia

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

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Entradas recientes

Buscar

Correo electrónico