junio 2015 – Exercitium

Pandigitales tridivisibles

PorJosé A. Alonso 26-junio-201528-octubre-2015

El número 4106357289 tiene la siguientes dos propiedades:

es pandigital, porque tiene todos los dígitos del 0 al 9 exactamente una vez y
es tridivisible, porque los sucesivos subnúmeros de tres dígitos (a partir del segundo) son divisibles por los sucesivos números primos; es decir, representado por d(i) el i-ésimo dígito, se tiene

     d(2)d(3)d(4)  = 106 es divisible por 2
     d(3)d(4)d(5)  = 063 es divisible por 3
     d(4)d(5)d(6)  = 635 es divisible por 5
     d(5)d(6)d(7)  = 357 es divisible por 7
     d(6)d(7)d(8)  = 572 es divisible por 11
     d(7)d(8)d(9)  = 728 es divisible por 13
     d(8)d(9)d(10) = 289 es divisible por 17

d(2)d(3)d(4) = 106 es divisible por 2

d(3)d(4)d(5) = 063 es divisible por 3

d(4)d(5)d(6) = 635 es divisible por 5

d(5)d(6)d(7) = 357 es divisible por 7

d(6)d(7)d(8) = 572 es divisible por 11

d(7)d(8)d(9) = 728 es divisible por 13

d(8)d(9)d(10) = 289 es divisible por 17

Definir la constante

   pandigitalesTridivisibles :: [Integer]

1	pandigitalesTridivisibles :: [Integer]

cuyos elementos son los números pandigitales tridivisibles. Por ejemplo,

   head pandigitalesTridivisibles  ==  4106357289
   sum pandigitalesTridivisibles   ==  16695334890

1 2	head pandigitalesTridivisibles == 4106357289 sum pandigitalesTridivisibles == 16695334890

Soluciones

import Data.List ((\\))

pandigitalesTridivisibles :: [Integer]
pandigitalesTridivisibles = 
    [entero [d1,d2,d3,d4,d5,d6,d7,d8,d9,d10] 
     | d2 <- [0..9]
     , d3 <- [0..9] \\ [d2]
     , d4 <- [0..9] \\ [d2,d3]
     , divisible [d2,d3,d4] 2
     , d5 <- [0..9] \\ [d2,d3,d4]  
     , divisible [d3,d4,d5] 3
     , d6 <- [0..9] \\ [d2,d3,d4,d5]  
     , divisible [d4,d5,d6] 5
     , d7 <- [0..9] \\ [d2,d3,d4,d5,d6]  
     , d8 <- [0..9] \\ [d2,d3,d4,d5,d6,d7]  
     , divisible [d6,d7,d8] 11
     , divisible [d5,d6,d7] 7
     , d9 <- [0..9] \\ [d2,d3,d4,d5,d6,d7,d8]  
     , divisible [d7,d8,d9] 13
     , d10 <- [0..9] \\ [d2,d3,d4,d5,d6,d7,d8,d9]  
     , divisible [d8,d9,d10] 17
     , d1 <- [1..9] \\ [d2,d3,d4,d5,d6,d7,d8,d9,d10]  
    ]

-- (entero ns) es el número cuyos dígitos son ns. Por ejemplo,
--    entero [3,2,5]  ==  325
entero :: [Integer] -> Integer
entero = read . concatMap show

divisible :: [Integer] -> Integer -> Bool
divisible xs n =
    entero xs `mod` n == 0

import Data.List ((\\))

pandigitalesTridivisibles :: [Integer]

pandigitalesTridivisibles =

[entero [d1,d2,d3,d4,d5,d6,d7,d8,d9,d10]

| d2 <- [0..9]

, d3 <- [0..9] \\ [d2]

, d4 <- [0..9] \\ [d2,d3]

, divisible [d2,d3,d4] 2

, d5 <- [0..9] \\ [d2,d3,d4]

, divisible [d3,d4,d5] 3

, d6 <- [0..9] \\ [d2,d3,d4,d5]

, divisible [d4,d5,d6] 5

, d7 <- [0..9] \\ [d2,d3,d4,d5,d6]

, d8 <- [0..9] \\ [d2,d3,d4,d5,d6,d7]

, divisible [d6,d7,d8] 11

, divisible [d5,d6,d7] 7

, d9 <- [0..9] \\ [d2,d3,d4,d5,d6,d7,d8]

, divisible [d7,d8,d9] 13

, d10 <- [0..9] \\ [d2,d3,d4,d5,d6,d7,d8,d9]

, divisible [d8,d9,d10] 17

, d1 <- [1..9] \\ [d2,d3,d4,d5,d6,d7,d8,d9,d10]

]

-- (entero ns) es el número cuyos dígitos son ns. Por ejemplo,

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

entero :: [Integer] -> Integer

entero = read . concatMap show

divisible :: [Integer] -> Integer -> Bool

divisible xs n =

entero xs `mod` n == 0

Partición por longitudes

PorJosé A. Alonso 23-junio-201528-octubre-2015

Definir la función

   particion :: [a] -> [Int] -> [[a]]

1	particion :: [a] -> [Int] -> [[a]]

tal que (particion xs ns) es la partición de xs donde la longitud de cada parte está determinada por los elementos de ns. Por ejemplo,

   particion [1..10] [2,5,0,3]  ==  [[1,2],[3,4,5,6,7],[],[8,9,10]] 
   particion [1..10] [1,4,2,3]  ==  [[1],[2,3,4,5],[6,7],[8,9,10]]

1 2	particion [1..10] [2,5,0,3] == [[1,2],[3,4,5,6,7],[],[8,9,10]] particion [1..10] [1,4,2,3] == [[1],[2,3,4,5],[6,7],[8,9,10]]

Soluciones

particion :: [a] -> [Int] -> [[a]]
particion [] []     = []
particion xs (n:ns) = take n xs : particion (drop n xs) ns

particion :: [a] -> [Int] -> [[a]]

particion [] [] = []

particion xs (n:ns) = take n xs : particion (drop n xs) ns

Máximos de una lista

PorJosé A. Alonso 22-junio-201528-octubre-2015

Definir la función

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

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

tal que (maximos xs) es la lista de los elementos de xs que son mayores que todos sus anteriores. Por ejemplo,

   maximos [1,-3,5,2,3,4,7,6,7]                         ==  [1,5,7]
   maximos "bafcdegag"                                  ==  "bfg"
   maximos (concat (replicate (10^6) "adxbcde")++"yz")  ==  "adxyz"
   length (maximos [1..10^6])                           ==  1000000

maximos [1,-3,5,2,3,4,7,6,7] == [1,5,7]

maximos "bafcdegag" == "bfg"

maximos (concat (replicate (10^6) "adxbcde")++"yz") == "adxyz"

length (maximos [1..10^6]) == 1000000

Soluciones

import Data.List (inits, nub)

-- 1ª definición (por comprensión)
maximos1 :: Ord a => [a] -> [a]
maximos1 xs =
    [x | (ys,x) <- zip (inits xs) xs, all (<x) ys]

-- 2ª definición (por recursión)
maximos2 :: Ord a => [a] -> [a]
maximos2 [] = []
maximos2 (x:xs) = x : maximos2 (filter (>x) xs)

-- 3ª definición (por recursión con acumulador)
maximos3 :: Ord a => [a] -> [a]
maximos3 [] = []
maximos3 (x:xs) = aux xs [x] x
    where aux [] zs _ = reverse zs
          aux (y:ys) zs m | y > m     = aux ys (y:zs) y
                          | otherwise = aux ys zs m 

-- 4ª definición (con scanl1):
maximos4 :: Ord a => [a] -> [a]
maximos4 = nub . scanl1 max 

-- Comparación de eficiencia
--    ghci> maximos1 (concat (replicate (10^3) "adxbcde") ++ "yz")
--    "adxyz"
--    (5.82 secs, 2859603952 bytes)
--    ghci> maximos2 (concat (replicate (10^3) "adxbcde") ++ "yz")
--    "adxyz"
--    (0.02 secs, 5879664 bytes)
--    ghci> maximos3 (concat (replicate (10^3) "adxbcde") ++ "yz")
--    "adxyz"
--    (0.03 secs, 4153680 bytes)
--    ghci> maximos4 (concat (replicate (10^3) "adxbcde") ++ "yz")
--    "adxyz"
--    (0.02 secs, 4163296 bytes)
--    
--    ghci> maximos2 (concat (replicate (10^6) "adxbcde") ++ "yz")
--    "adxyz"
--    (3.64 secs, 1314485320 bytes)
--    ghci> maximos3 (concat (replicate (10^6) "adxbcde") ++ "yz")
--    "adxyz"
--    (6.59 secs, 1706434544 bytes)
--    ghci> maximos4 (concat (replicate (10^6) "adxbcde") ++ "yz")
--    "adxyz"
--    (0.89 secs, 1594567000 bytes)
--    
--    ghci> length (maximos2 [1..10^4])
--    10000
--    (17.34 secs, 4302913816 bytes)
--    ghci> length (maximos3 [1..10^4])
--    10000
--    (0.03 secs, 6602488 bytes)
--    ghci> length (maximos4 [1..10^4])
--    10000
--    (1.37 secs, 6820408 bytes)

import Data.List (inits, nub)

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

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

maximos1 xs =

[x | (ys,x) <- zip (inits xs) xs, all (<x) ys]

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

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

maximos2 [] = []

maximos2 (x:xs) = x : maximos2 (filter (>x) xs)

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

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

maximos3 [] = []

maximos3 (x:xs) = aux xs [x] x

where aux [] zs _ = reverse zs

aux (y:ys) zs m | y > m = aux ys (y:zs) y

| otherwise = aux ys zs m

-- 4ª definición (con scanl1):

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

maximos4 = nub . scanl1 max

-- Comparación de eficiencia

-- ghci> maximos1 (concat (replicate (10^3) "adxbcde") ++ "yz")

-- "adxyz"

-- (5.82 secs, 2859603952 bytes)

-- ghci> maximos2 (concat (replicate (10^3) "adxbcde") ++ "yz")

-- "adxyz"

-- (0.02 secs, 5879664 bytes)

-- ghci> maximos3 (concat (replicate (10^3) "adxbcde") ++ "yz")

-- "adxyz"

-- (0.03 secs, 4153680 bytes)

-- ghci> maximos4 (concat (replicate (10^3) "adxbcde") ++ "yz")

-- "adxyz"

-- (0.02 secs, 4163296 bytes)

-- ghci> maximos2 (concat (replicate (10^6) "adxbcde") ++ "yz")

-- "adxyz"

-- (3.64 secs, 1314485320 bytes)

-- ghci> maximos3 (concat (replicate (10^6) "adxbcde") ++ "yz")

-- "adxyz"

-- (6.59 secs, 1706434544 bytes)

-- ghci> maximos4 (concat (replicate (10^6) "adxbcde") ++ "yz")

-- "adxyz"

-- (0.89 secs, 1594567000 bytes)

-- ghci> length (maximos2 [1..10^4])

-- 10000

-- (17.34 secs, 4302913816 bytes)

-- ghci> length (maximos3 [1..10^4])

-- 10000

-- (0.03 secs, 6602488 bytes)

-- ghci> length (maximos4 [1..10^4])

-- 10000

-- (1.37 secs, 6820408 bytes)

Menor n tal que el primo n-ésimo cumple una propiedad

PorJosé A. Alonso 18-junio-201528-octubre-2015

Sea p(n) el n-ésimo primo y sea r el resto de dividir (p(n)-1)^n + (p(n)+1)^n por p(n)^2. Por ejemplo,

   si n = 3, entonces p(3) =  5 y r = ( 4^3 +  6^3) mod  (5^2) =   5
   si n = 7, entonces p(7) = 17 y r = (16^7 + 18^7) mod (17^2) = 238

1 2	si n = 3, entonces p(3) = 5 y r = ( 4^3 + 6^3) mod (5^2) = 5 si n = 7, entonces p(7) = 17 y r = (16^7 + 18^7) mod (17^2) = 238

Definir la función

   menorPR :: Integer -> Integer

1	menorPR :: Integer -> Integer

tal que (menorPR x) es el menor n tal que el resto de dividir (p(n)-1)^n + (p(n)+1)^n por p(n)^2 es mayor que x. Por ejemplo,

   menorPR 100     == 5
   menorPR 345     == 9
   menorPR 1000    == 13
   menorPR (10^9)  == 7037
   menorPR (10^10) == 21035
   menorPR (10^12) == 191041

menorPR 100 == 5

menorPR 345 == 9

menorPR 1000 == 13

menorPR (10^9) == 7037

menorPR (10^10) == 21035

menorPR (10^12) == 191041

Soluciones

import Data.Numbers.Primes (primes)

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

menorPR1 :: Integer -> Integer
menorPR1 x = head [n | (n,p) <- zip [1..] primes
                     , (((p-1)^n + (p+1)^n) `mod` (p^2)) > x]

-- Segunda solución (usando el binomio de Newton)
-- ==============================================

-- Desarrollando por el binomio de Newton
--    (p+1)^n = C(n,0)p^n + C(n,1)p^(n-1) +...+ C(n,n-1)p + 1
--    (p-1)^n = C(n,0)p^n - C(n,1)p^(n-1) +...+ C(n,n-1)p + (-1)^n
-- Sumando se obtiene (según n sea par o impar)
--    2*C(n,0)p^n + 2*C(n,n-2)p^(n-1) +...+ 2*C(n,2)p^2 + 2
--    2*C(n,0)p^n + 2*C(n,n-2)p^(n-1) +...+ 2*C(n,1)p^1
-- Al dividir por p^2, el resto es (según n sea par o impar) 2 ó 2*C(n,1)p

-- (restoM n p) es el resto de de dividir (p-1)^n + (p+1)^n por p^2.
restoM :: Integer -> Integer -> Integer
restoM n p | even n    = 2
           | otherwise = 2*n*p `mod`(p^2)

menorPR2 :: Integer -> Integer
menorPR2 x = head [n | (n,p) <- zip [1..] primes, restoM n p > x]

-- Comparación de eficiencia
--    ghci> menorPR1 (3*10^8)
--    3987
--    (2.44 secs, 120291676 bytes)
--    ghci> menorPR2 (3*10^8)
--    3987
--    (0.04 secs, 8073900 bytes)

import Data.Numbers.Primes (primes)

-- 1ª solución

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

menorPR1 :: Integer -> Integer

menorPR1 x = head [n | (n,p) <- zip [1..] primes

, (((p-1)^n + (p+1)^n) `mod` (p^2)) > x]

-- Segunda solución (usando el binomio de Newton)

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

-- Desarrollando por el binomio de Newton

-- (p+1)^n = C(n,0)p^n + C(n,1)p^(n-1) +...+ C(n,n-1)p + 1

-- (p-1)^n = C(n,0)p^n - C(n,1)p^(n-1) +...+ C(n,n-1)p + (-1)^n

-- Sumando se obtiene (según n sea par o impar)

-- 2*C(n,0)p^n + 2*C(n,n-2)p^(n-1) +...+ 2*C(n,2)p^2 + 2

-- 2*C(n,0)p^n + 2*C(n,n-2)p^(n-1) +...+ 2*C(n,1)p^1

-- Al dividir por p^2, el resto es (según n sea par o impar) 2 ó 2*C(n,1)p

-- (restoM n p) es el resto de de dividir (p-1)^n + (p+1)^n por p^2.

restoM :: Integer -> Integer -> Integer

restoM n p | even n = 2

| otherwise = 2*n*p `mod`(p^2)

menorPR2 :: Integer -> Integer

menorPR2 x = head [n | (n,p) <- zip [1..] primes, restoM n p > x]

-- Comparación de eficiencia

-- ghci> menorPR1 (3*10^8)

-- 3987

-- (2.44 secs, 120291676 bytes)

-- ghci> menorPR2 (3*10^8)

-- 3987

-- (0.04 secs, 8073900 bytes)

Polinomio cromático de un grafo

PorJosé A. Alonso 17-junio-201526-marzo-2022

El polinomio cromático de un grafo calcula el número de maneras en las cuales puede ser coloreado el grafo usando un número de colores dado, de forma que dos vértices adyacentes no tengan el mismo color.

En el caso del grafo completo de n vértices, su polinomio cromático es

   P(n,x) = x(x-1)(x-2) ... (x-(n-1))

1	P(n,x) = x(x-1)(x-2) ... (x-(n-1))

Por ejemplo,

   P(3,x) = x(x-1)(x-2)      = x^3 - 3*x^2 + 2*x
   P(4,x) = x(x-1)(x-2)(x-3) = x^4 - 6*x^3 + 11*x^2 - 6*x

1 2	P(3,x) = x(x-1)(x-2) = x^3 - 3x^2 + 2x P(4,x) = x(x-1)(x-2)(x-3) = x^4 - 6x^3 + 11x^2 - 6*x

Lo que significa que P(4)(x) es el número de formas de colorear el grafo completo de 4 vértices con x colores. Por tanto,

   P(4,2) =  0 (no se puede colorear con 2 colores)
   P(4,4) = 24 (hay 24 formas de colorearlo con 4 colores)

1 2	P(4,2) = 0 (no se puede colorear con 2 colores) P(4,4) = 24 (hay 24 formas de colorearlo con 4 colores)

Definir la función

   polGC:: Int -> Polinomio Int

1	polGC:: Int -> Polinomio Int

tal que (polGC n) es el polinomio cromático del grafo completo de n vértices. Por ejemplo,

   polGC 4  ==  x^4 + -6*x^3 + 11*x^2 + -6*x
   polGC 5  ==  x^5 + -10*x^4 + 35*x^3 + -50*x^2 + 24*x

1 2	polGC 4 == x^4 + -6x^3 + 11x^2 + -6x polGC 5 == x^5 + -10x^4 + 35x^3 + -50x^2 + 24*x

Comprobar con QuickCheck que si el número de colores (x) coincide con el número de vértices del grafo (n), el número de maneras de colorear el grafo es n!.

Nota 1. Al hacer la comprobación limitar el tamaño de las pruebas como se indica a continuación

   ghci> quickCheckWith (stdArgs {maxSize=7}) prop_polGC
   +++ OK, passed 100 tests.

1 2	ghci> quickCheckWith (stdArgs {maxSize=7}) prop_polGC +++ OK, passed 100 tests.

Nota 2: Este ejercicio debe realizarse usando únicamente las funciones de la librería de polinomios (I1M.PolOperaciones) que se describe aquí y se encuentra aquí.

Soluciones

import I1M.PolOperaciones
import Test.QuickCheck

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

polGC :: Int -> Polinomio Int
polGC 0 = consPol 0 1 polCero
polGC n = polGC (n-1) `multPol` consPol 1 1 (consPol 0 (-n+1) polCero)

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

polGC2 :: Int -> Polinomio Int
polGC2 n = multLista (map polMon [0..n-1])

-- (polMon n) es el monomio x-n. Por ejemplo,
--    polMon 3  ==  1*x + -3
polMon:: Int -> Polinomio Int
polMon n = consPol 1 1 (consPol 0 (-n) polCero)

-- (multLista ps) es el producto de la lista de polinomios ps.
multLista :: [Polinomio Int] -> Polinomio Int
multLista []     = polUnidad
multLista (p:ps) = multPol p (multLista ps)

-- 3ª solución (por plegado)
-- =========================

polGC3 :: Int -> Polinomio Int
polGC3 n = foldl multPol polUnidad
           [consPol 1 1 (consPol 0 (-i) polCero) | i <- [0..n-1]]

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

-- La propiedad es
prop_polGC :: Int -> Property
prop_polGC n = 
    n > 0 ==> valor (polGC n) n == product [1..n]

-- La comprobación es
--    ghci> quickCheckWith (stdArgs {maxSize=7}) prop_polGC
--    +++ OK, passed 100 tests.
--    (0.04 secs, 7785800 bytes)

import I1M.PolOperaciones

import Test.QuickCheck

-- 1ª solución

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

polGC :: Int -> Polinomio Int

polGC 0 = consPol 0 1 polCero

polGC n = polGC (n-1) `multPol` consPol 1 1 (consPol 0 (-n+1) polCero)

-- 2ª solución

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

polGC2 :: Int -> Polinomio Int

polGC2 n = multLista (map polMon [0..n-1])

-- (polMon n) es el monomio x-n. Por ejemplo,

-- polMon 3 == 1*x + -3

polMon:: Int -> Polinomio Int

polMon n = consPol 1 1 (consPol 0 (-n) polCero)

-- (multLista ps) es el producto de la lista de polinomios ps.

multLista :: [Polinomio Int] -> Polinomio Int

multLista [] = polUnidad

multLista (p:ps) = multPol p (multLista ps)

-- 3ª solución (por plegado)

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

polGC3 :: Int -> Polinomio Int

polGC3 n = foldl multPol polUnidad

[consPol 1 1 (consPol 0 (-i) polCero) | i <- [0..n-1]]

-- Comprobación

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

-- La propiedad es

prop_polGC :: Int -> Property

prop_polGC n =

n > 0 ==> valor (polGC n) n == product [1..n]

-- La comprobación es

-- ghci> quickCheckWith (stdArgs {maxSize=7}) prop_polGC

-- +++ OK, passed 100 tests.

-- (0.04 secs, 7785800 bytes)

Suma de posteriores

PorJosé A. Alonso 16-junio-201528-octubre-2015

Definir la función

   sumaPosteriores :: [Int] -> [Int]

1	sumaPosteriores :: [Int] -> [Int]

tal que (sumaPosteriores xs) es la lista obtenida sustituyendo cada elemento de xs por la suma de los elementos posteriores. Por ejemplo,

   sumaPosteriores [1..8]        == [35,33,30,26,21,15,8,0]
   sumaPosteriores [1,-3,2,5,-8] == [-4,-1,-3,-8,0]

1 2	sumaPosteriores [1..8] == [35,33,30,26,21,15,8,0] sumaPosteriores [1,-3,2,5,-8] == [-4,-1,-3,-8,0]

Comprobar con QuickCheck que el último elemento de la lista (sumaPosteriores xs) siempre es 0.

Soluciones

import Data.List (tails)
import Test.QuickCheck

-- 1ª definición (por recursión):
sumaPosteriores1 :: [Int] -> [Int]
sumaPosteriores1 []     = []
sumaPosteriores1 (x:xs) = sum xs : sumaPosteriores1 xs

-- 2ª definición (sin argumentos)
sumaPosteriores2 :: [Int] -> [Int]
sumaPosteriores2 = map (sum . tail) . init . tails

-- La propiedad es
propSumaP:: [Int] -> Property
propSumaP xs = not (null xs) ==> last (sumaPosteriores1 xs) == 0

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

import Data.List (tails)

import Test.QuickCheck

-- 1ª definición (por recursión):

sumaPosteriores1 :: [Int] -> [Int]

sumaPosteriores1 [] = []

sumaPosteriores1 (x:xs) = sum xs : sumaPosteriores1 xs

-- 2ª definición (sin argumentos)

sumaPosteriores2 :: [Int] -> [Int]

sumaPosteriores2 = map (sum . tail) . init . tails

-- La propiedad es

propSumaP:: [Int] -> Property

propSumaP xs = not (null xs) ==> last (sumaPosteriores1 xs) == 0

-- La comprobación es

-- ghci> quickCheck propSumaP

-- +++ OK, passed 100 tests.

Algoritmo de bajada para resolver un sistema triangular inferior

PorJosé A. Alonso 15-junio-201528-octubre-2015

Un sistema de ecuaciones lineales Ax = b es triangular inferior si todos los elementos de la matriz A que están por encima de la diagonal principal son nulos; es decir, es de la forma

   a(1,1)*x(1)                                               = b(1)
   a(2,1)*x(1) + a(2,2)*x(2)                                 = b(2)
   a(3,1)*x(1) + a(3,2)*x(2) + a(3,3)*x(3)                   = b(3)
   ...
   a(n,1)*x(1) + a(n,2)*x(2) + a(n,3)*x(3) +...+ a(x,x)*x(n) = b(n)

a(1,1)*x(1) = b(1)

a(2,1)*x(1) + a(2,2)*x(2) = b(2)

a(3,1)*x(1) + a(3,2)*x(2) + a(3,3)*x(3) = b(3)

...

a(n,1)*x(1) + a(n,2)*x(2) + a(n,3)*x(3) +...+ a(x,x)*x(n) = b(n)

El sistema es compatible si, y sólo si, el producto de los elementos de la diagonal principal es distinto de cero. En este caso, la solución se puede calcular mediante el algoritmo de bajada:

   x(1) = b(1) / a(1,1)
   x(2) = (b(2) - a(2,1)*x(1)) / a(2,2)
   x(3) = (b(3) - a(3,1)*x(1) - a(3,2)*x(2)) / a(3,3)
   ...
   x(n) = (b(n) - a(n,1)*x(1) - a(n,2)*x(2) -...- a(n,n-1)*x(n-1)) / a(n,n)

x(1) = b(1) / a(1,1)

x(2) = (b(2) - a(2,1)*x(1)) / a(2,2)

x(3) = (b(3) - a(3,1)*x(1) - a(3,2)*x(2)) / a(3,3)

...

x(n) = (b(n) - a(n,1)*x(1) - a(n,2)*x(2) -...- a(n,n-1)*x(n-1)) / a(n,n)

Definir la función

   bajada :: Matrix Double -> Matrix Double -> Matrix Double

1	bajada :: Matrix Double -> Matrix Double -> Matrix Double

tal que (bajada a b) es la solución, mediante el algoritmo de bajada, del sistema compatible triangular superior ax = b. Por ejemplo,

   ghci> let a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]
   ghci> let b = fromLists [[3],[6.5],[10]]
   ghci> bajada a b
   ( 1.5 )
   ( 2.0 )
   ( 0.0 )

ghci> let a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]

ghci> let b = fromLists [[3],[6.5],[10]]

ghci> bajada a b

( 1.5 )

( 2.0 )

( 0.0 )

Es decir, la solución del sistema

   2x            = 3
   3x + y        = 6.5
   4x + 2y + 5 z = 10

2x = 3

3x + y = 6.5

4x + 2y + 5 z = 10

es x=1.5, y=2 y z=0.

Soluciones

import Data.Matrix

bajada :: Matrix Double -> Matrix Double -> Matrix Double
bajada a b = fromLists [[x i] | i <- [1..m]]
    where m = nrows a
          x k = (b!(k,1) - sum [a!(k,j) * x j | j <- [1..k-1]]) / a!(k,k)

import Data.Matrix

bajada :: Matrix Double -> Matrix Double -> Matrix Double

bajada a b = fromLists [[x i] | i <- [1..m]]

where m = nrows a

x k = (b!(k,1) - sum [a!(k,j) * x j | j <- [1..k-1]]) / a!(k,k)

Mayor producto de n números adyacentes en una matriz

PorJosé A. Alonso 12-junio-201528-octubre-2015

Definir la función

   mayorProductoAdyacentes :: (Num a, Ord a) => Int -> Matriz a -> [[a]]

1	mayorProductoAdyacentes :: (Num a, Ord a) => Int -> Matriz a -> [[a]]

tal que (mayorProductoAdyacentes n p) es la lista de los segmentos formados por n elementos adyacentes en la misma fila, columna o diagonal de la matriz p cuyo productos son máximo. Por ejemplo,

   ghci> mayorProductoAdyacentes 3 (listArray ((1,1),(3,4)) [1..12])
   [[10,11,12]]
   ghci> mayorProductoAdyacentes 3 (listArray ((1,1),(3,4)) [1,3,4,5, 0,7,2,1, 3,9,2,1])
   [[3,7,9]]
   ghci> mayorProductoAdyacentes 2 (listArray ((1,1),(2,3)) [1,3,4, 0,3,2])
   [[3,4],[4,3]]
   ghci> mayorProductoAdyacentes 2 (listArray ((1,1),(2,3)) [1,2,1, 3,0,3])
   [[2,3],[2,3]]
   ghci> mayorProductoAdyacentes 2 (listArray ((1,1),(2,3)) [1,2,1, 3,4,3])
   [[3,4],[4,3]]
   ghci> mayorProductoAdyacentes 2 (listArray ((1,1),(2,3)) [1,5,1, 3,4,3])
   [[5,4]]
   ghci> mayorProductoAdyacentes 3 (listArray ((1,1),(3,4)) [1,3,4,5, 0,7,2,1, 3,9,2,1])
   [[3,7,9]]

ghci> mayorProductoAdyacentes 3 (listArray ((1,1),(3,4)) [1..12])

[[10,11,12]]

ghci> mayorProductoAdyacentes 3 (listArray ((1,1),(3,4)) [1,3,4,5, 0,7,2,1, 3,9,2,1])

[[3,7,9]]

ghci> mayorProductoAdyacentes 2 (listArray ((1,1),(2,3)) [1,3,4, 0,3,2])

[[3,4],[4,3]]

ghci> mayorProductoAdyacentes 2 (listArray ((1,1),(2,3)) [1,2,1, 3,0,3])

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

ghci> mayorProductoAdyacentes 2 (listArray ((1,1),(2,3)) [1,2,1, 3,4,3])

[[3,4],[4,3]]

ghci> mayorProductoAdyacentes 2 (listArray ((1,1),(2,3)) [1,5,1, 3,4,3])

[[5,4]]

ghci> mayorProductoAdyacentes 3 (listArray ((1,1),(3,4)) [1,3,4,5, 0,7,2,1, 3,9,2,1])

[[3,7,9]]

Soluciones

import Data.Array

type Matriz a = Array (Int,Int) a

mayorProductoAdyacentes :: (Num a, Ord a) => Int -> Matriz a -> [[a]]
mayorProductoAdyacentes n p = 
    [xs | xs <- segmentos, product xs == m]
    where segmentos = adyacentes n p 
          m         = maximum (map product segmentos)

-- (adyacentes n p) es la lista de los segmentos de longitud n de las
-- líneas de la matriz p. Por ejemplo,
--    ghci> adyacentes 3 (listArray ((1,1),(3,4)) [1..12])
--    [[1,2,3],[2,3,4],[5,6,7],[6,7,8],[9,10,11],[10,11,12],[1,5,9],
--     [2,6,10],[3,7,11],[4,8,12],[1,6,11],[2,7,12],[3,6,9],[4,7,10]]
adyacentes :: Num a => Int -> Matriz a -> [[a]]
adyacentes n p = concatMap (segmentos n) (lineas p)

-- (lineas p) es la lista de las líneas de la matriz p. Por ejemplo, 
--    ghci> lineas (listArray ((1,1),(3,4)) [1..12])
--    [[1,2,3,4],[5,6,7,8],[9,10,11,12],
--     [1,5,9],[2,6,10],[3,7,11],[4,8,12],
--     [9],[5,10],[1,6,11],[2,7,12],[3,8],[4],
--     [1],[2,5],[3,6,9],[4,7,10],[8,11],[12]]
lineas :: Num a => Matriz a -> [[a]]
lineas p = filas p ++ 
           columnas p ++ 
           diagonalesPrincipales p ++
           diagonalesSecundarias p

-- (filas p) es la lista de las filas de la matriz p. Por ejemplo,
--    ghci> filas (listArray ((1,1),(3,4)) [1..12])
--    [[1,2,3,4],[5,6,7,8],[9,10,11,12]]
filas :: Num a => Matriz a -> [[a]]
filas p = [fila i p | i <- [1..m]]
    where (_,(m,_)) = bounds p

-- (fila i p) es la fila i-ésima de la matriz p. Por ejemplo, 
--    fila 2 (listArray ((1,1),(3,4)) [1..12])  ==  [5,6,7,8]
fila :: Num a => Int -> Matriz a -> [a]
fila i p = [p!(i,j) | j <- [1..n]]
    where (_,(_,n)) = bounds p

-- (columnas p) es la lista de las columnas de la matriz p. Por ejemplo,
--    ghci> columnas (listArray ((1,1),(3,4)) [1..12])
--    [[1,5,9],[2,6,10],[3,7,11],[4,8,12]]
columnas :: Num a => Matriz a -> [[a]]
columnas p = [columna j p | j <- [1..n]]
    where (_,(_,n)) = bounds p

-- (columna j p) es la columna j-ésima de la matriz p. Por ejemplo, 
--    columna 2 (listArray ((1,1),(3,4)) [1..12])  ==  [2,6,10]
columna :: Num a => Int -> Matriz a -> [a]
columna j p = [p!(i,j) | i <- [1..m]]
    where (_,(m,_)) = bounds p

-- (diagonalesPrincipales p) es la lista de las diagonales principales
-- de p. Por ejemplo, 
--    ghci> diagonalesPrincipales (listArray ((1,1),(3,4)) [1..12])
--    [[9],[5,10],[1,6,11],[2,7,12],[3,8],[4]]
diagonalesPrincipales :: Matriz a -> [[a]]
diagonalesPrincipales p = 
    [[p!ij1 | ij1 <- extension ij] | ij <- iniciales] 
    where (_,(m,n)) = bounds p
          iniciales = [(i,1) | i <- [m,m-1..2]] ++ [(1,j) | j <- [1..n]] 
          extension (i,j) = [(i+k,j+k) | k <- [0..min (m-i) (n-j)]]

-- (diagonalesSecundarias p) es la lista de las diagonales secundarias
-- de p. Por ejemplo, 
--    ghci> diagonalesSecundarias (listArray ((1,1),(3,4)) [1..12])
--    [[1],[2,5],[3,6,9],[4,7,10],[8,11],[12]]
diagonalesSecundarias p = 
    [[p!ij1 | ij1 <- extension ij] | ij <- iniciales]
    where (_,(m,n)) = bounds p
          iniciales = [(1,j) | j <- [1..n]] ++ [(i,n) | i <- [2..m]] 
          extension (i,j) = [(i+k,j-k) | k <- [0..min (j-1) (m-i)]]

-- (segmentos n xs) es la lista de los segmentos de longitud n de la
-- lista xs. Por ejemplo, 
--    segmentos 3 [1..5]  ==  [[1,2,3],[2,3,4],[3,4,5]]
segmentos :: Int -> [a] -> [[a]]
segmentos n xs 
    | length xs < n = []
    | otherwise     = take n xs : segmentos n (tail xs)

import Data.Array

type Matriz a = Array (Int,Int) a

mayorProductoAdyacentes :: (Num a, Ord a) => Int -> Matriz a -> [[a]]

mayorProductoAdyacentes n p =

[xs | xs <- segmentos, product xs == m]

where segmentos = adyacentes n p

m = maximum (map product segmentos)

-- (adyacentes n p) es la lista de los segmentos de longitud n de las

-- líneas de la matriz p. Por ejemplo,

-- ghci> adyacentes 3 (listArray ((1,1),(3,4)) [1..12])

-- [[1,2,3],[2,3,4],[5,6,7],[6,7,8],[9,10,11],[10,11,12],[1,5,9],

-- [2,6,10],[3,7,11],[4,8,12],[1,6,11],[2,7,12],[3,6,9],[4,7,10]]

adyacentes :: Num a => Int -> Matriz a -> [[a]]

adyacentes n p = concatMap (segmentos n) (lineas p)

-- (lineas p) es la lista de las líneas de la matriz p. Por ejemplo,

-- ghci> lineas (listArray ((1,1),(3,4)) [1..12])

-- [[1,2,3,4],[5,6,7,8],[9,10,11,12],

-- [1,5,9],[2,6,10],[3,7,11],[4,8,12],

-- [9],[5,10],[1,6,11],[2,7,12],[3,8],[4],

-- [1],[2,5],[3,6,9],[4,7,10],[8,11],[12]]

lineas :: Num a => Matriz a -> [[a]]

lineas p = filas p ++

columnas p ++

diagonalesPrincipales p ++

diagonalesSecundarias p

-- (filas p) es la lista de las filas de la matriz p. Por ejemplo,

-- ghci> filas (listArray ((1,1),(3,4)) [1..12])

-- [[1,2,3,4],[5,6,7,8],[9,10,11,12]]

filas :: Num a => Matriz a -> [[a]]

filas p = [fila i p | i <- [1..m]]

where (_,(m,_)) = bounds p

-- (fila i p) es la fila i-ésima de la matriz p. Por ejemplo,

-- fila 2 (listArray ((1,1),(3,4)) [1..12]) == [5,6,7,8]

fila :: Num a => Int -> Matriz a -> [a]

fila i p = [p!(i,j) | j <- [1..n]]

where (_,(_,n)) = bounds p

-- (columnas p) es la lista de las columnas de la matriz p. Por ejemplo,

-- ghci> columnas (listArray ((1,1),(3,4)) [1..12])

-- [[1,5,9],[2,6,10],[3,7,11],[4,8,12]]

columnas :: Num a => Matriz a -> [[a]]

columnas p = [columna j p | j <- [1..n]]

where (_,(_,n)) = bounds p

-- (columna j p) es la columna j-ésima de la matriz p. Por ejemplo,

-- columna 2 (listArray ((1,1),(3,4)) [1..12]) == [2,6,10]

columna :: Num a => Int -> Matriz a -> [a]

columna j p = [p!(i,j) | i <- [1..m]]

where (_,(m,_)) = bounds p

-- (diagonalesPrincipales p) es la lista de las diagonales principales

-- de p. Por ejemplo,

-- ghci> diagonalesPrincipales (listArray ((1,1),(3,4)) [1..12])

-- [[9],[5,10],[1,6,11],[2,7,12],[3,8],[4]]

diagonalesPrincipales :: Matriz a -> [[a]]

diagonalesPrincipales p =

[[p!ij1 | ij1 <- extension ij] | ij <- iniciales]

where (_,(m,n)) = bounds p

iniciales = [(i,1) | i <- [m,m-1..2]] ++ [(1,j) | j <- [1..n]]

extension (i,j) = [(i+k,j+k) | k <- [0..min (m-i) (n-j)]]

-- (diagonalesSecundarias p) es la lista de las diagonales secundarias

-- de p. Por ejemplo,

-- ghci> diagonalesSecundarias (listArray ((1,1),(3,4)) [1..12])

-- [[1],[2,5],[3,6,9],[4,7,10],[8,11],[12]]

diagonalesSecundarias p =

[[p!ij1 | ij1 <- extension ij] | ij <- iniciales]

where (_,(m,n)) = bounds p

iniciales = [(1,j) | j <- [1..n]] ++ [(i,n) | i <- [2..m]]

extension (i,j) = [(i+k,j-k) | k <- [0..min (j-1) (m-i)]]

-- (segmentos n xs) es la lista de los segmentos de longitud n de la

-- lista xs. Por ejemplo,

-- segmentos 3 [1..5] == [[1,2,3],[2,3,4],[3,4,5]]

segmentos :: Int -> [a] -> [[a]]

segmentos n xs

| length xs < n = []

| otherwise = take n xs : segmentos n (tail xs)

Método de bisección para encontrar raíces de funciones

PorJosé A. Alonso 11-junio-201528-octubre-2015

El método de bisección para calcular un cero de una función en el intervalo
[a,b] se basa en el teorema de Bolzano:

Si f(x) es una función continua en el intervalo [a, b], y si, además, en los extremos del intervalo la función f(x) toma valores de signo opuesto (f(a) * f(b) < 0), entonces existe al menos un valor c en (a, b) para el que f(c) = 0″.

El método para calcular un cero de la función f en el intervalo [a,b] con un error menor que e consiste en tomar el punto medio del intervalo c = (a+b)/2 y considerar los siguientes casos:

Si |f(c)| < e, hemos encontrado una aproximación del punto que anula f en el intervalo con un error aceptable.
Si f(c) tiene signo distinto de f(a), repetir el proceso en el intervalo [a,c].
Si no, repetir el proceso en el intervalo [c,b].

Definir la función

   biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double

1	biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double

tal que (biseccion f a b e) es una aproximación del punto del intervalo [a,b] en el que se anula la función f, con un error menor que e, calculada mediante el método de la bisección. Por ejemplo,

   biseccion (\x -> x^2 - 3) 0 5 0.01             ==  1.7333984375
   biseccion (\x -> x^3 - x - 2) 0 4 0.01         ==  1.521484375
   biseccion cos 0 2 0.01                         ==  1.5625
   biseccion (\x -> log (50-x) - 4) (-10) 3 0.01  ==  -5.125

biseccion (\x -> x^2 - 3) 0 5 0.01 == 1.7333984375

biseccion (\x -> x^3 - x - 2) 0 4 0.01 == 1.521484375

biseccion cos 0 2 0.01 == 1.5625

biseccion (\x -> log (50-x) - 4) (-10) 3 0.01 == -5.125

Soluciones

-- 1ª solución
biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double
biseccion f a b e  
    | abs (f c) < e   = c
    | (f a)*(f c) < 0 = biseccion f a c e
    | otherwise       = biseccion f c b e
    where c = (a+b)/2

-- 2ª solución
biseccion2 :: (Double -> Double) -> Double -> Double -> Double -> Double
biseccion2 f a b e = aux a b
    where aux a b | abs (f c) < e   = c
                  | (f a)*(f c) < 0 = aux a c 
                  | otherwise       = aux c b
                  where c = (a+b)/2

-- 1ª solución

biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double

biseccion f a b e

| abs (f c) < e = c

| (f a)*(f c) < 0 = biseccion f a c e

| otherwise = biseccion f c b e

where c = (a+b)/2

-- 2ª solución

biseccion2 :: (Double -> Double) -> Double -> Double -> Double -> Double

biseccion2 f a b e = aux a b

where aux a b | abs (f c) < e = c

| (f a)*(f c) < 0 = aux a c

| otherwise = aux c b

where c = (a+b)/2

Descomposiciones en sumas de primos

PorJosé A. Alonso 4-junio-201513-junio-2015

Definir la función

   sumaDePrimos :: Int -> [[Int]]

1	sumaDePrimos :: Int -> [[Int]]

tal que (sumaDePrimos x) es la lista de las listas no crecientes de números primos que suman x. Por ejemplo,

   sumaDePrimos 10  ==  [[7,3],[5,5],[5,3,2],[3,3,2,2],[2,2,2,2,2]]
   sumaDePrimos 0   ==  []
   sumaDePrimos 1   ==  []

sumaDePrimos 10 == [[7,3],[5,5],[5,3,2],[3,3,2,2],[2,2,2,2,2]]

sumaDePrimos 0 == []

sumaDePrimos 1 == []

Soluciones

import Data.Numbers.Primes (primes)

sumaDePrimos :: Integer -> [[Integer]]
sumaDePrimos n = aux n (reverse (takeWhile (<=n) primes))
    where aux _ [] = []
          aux n (x:xs) | x > n     = aux n xs
                       | x == n    = [n] : aux n xs
                       | otherwise = map (x:) (aux (n-x) (x:xs)) ++ aux n xs

import Data.Numbers.Primes (primes)

sumaDePrimos :: Integer -> [[Integer]]

sumaDePrimos n = aux n (reverse (takeWhile (<=n) primes))

where aux _ [] = []

aux n (x:xs) | x > n = aux n xs

| x == n = [n] : aux n xs

| otherwise = map (x:) (aux (n-x) (x:xs)) ++ aux n xs

Caminos en un grafo

PorJosé A. Alonso 2-junio-201526-marzo-2022

Definir las funciones

   grafo   :: [(Int,Int)] -> Grafo Int Int
   caminos :: Grafo Int Int -> Int -> Int -> [[Int]]

1 2	grafo :: [(Int,Int)] -> Grafo Int Int caminos :: Grafo Int Int -> Int -> Int -> [[Int]]

tales que

(grafo as) es el grafo no dirigido definido cuyas aristas son as. Por ejemplo,

     ghci> grafo [(2,4),(4,5)]
     G ND (array (2,5) [(2,[(4,0)]),(3,[]),(4,[(2,0),(5,0)]),(5,[(4,0)])])

1 2	ghci> grafo [(2,4),(4,5)] G ND (array (2,5) [(2,[(4,0)]),(3,[]),(4,[(2,0),(5,0)]),(5,[(4,0)])])

(caminos g a b) es la lista los caminos en el grafo g desde a hasta b sin pasar dos veces por el mismo nodo. Por ejemplo,

     ghci> sort (caminos (grafo [(1,3),(2,5),(3,5),(3,7),(5,7)]) 1 7)
     [[1,3,5,7],[1,3,7]]
     ghci> sort (caminos (grafo [(1,3),(2,5),(3,5),(3,7),(5,7)]) 2 7)
     [[2,5,3,7],[2,5,7]]
     ghci> sort (caminos (grafo [(1,3),(2,5),(3,5),(3,7),(5,7)]) 1 2)
     [[1,3,5,2],[1,3,7,5,2]]
     ghci> caminos (grafo [(1,3),(2,5),(3,5),(3,7),(5,7)]) 1 4
     []
     ghci> length (caminos (grafo [(i,j) | i <- [1..10], j <- [i..10]]) 1 10)
     109601

ghci> sort (caminos (grafo [(1,3),(2,5),(3,5),(3,7),(5,7)]) 1 7)

[[1,3,5,7],[1,3,7]]

ghci> sort (caminos (grafo [(1,3),(2,5),(3,5),(3,7),(5,7)]) 2 7)

[[2,5,3,7],[2,5,7]]

ghci> sort (caminos (grafo [(1,3),(2,5),(3,5),(3,7),(5,7)]) 1 2)

[[1,3,5,2],[1,3,7,5,2]]

ghci> caminos (grafo [(1,3),(2,5),(3,5),(3,7),(5,7)]) 1 4

[]

ghci> length (caminos (grafo [(i,j) | i <- [1..10], j <- [i..10]]) 1 10)

109601

Nota: Este ejercicio debe realizarse usando únicamente las funciones de la librería de grafos (I1M.Grafo) que se describe aquí y se encuentra aquí.

Soluciones

import Data.List (sort)
import I1M.Grafo
import I1M.BusquedaEnEspaciosDeEstados

grafo :: [(Int,Int)] -> Grafo Int Int
grafo as = creaGrafo ND (m,n) [(x,y,0) | (x,y) <- as]
    where ns = map fst as ++ map snd as
          m  = minimum ns
          n  = maximum ns

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

caminos :: Grafo Int Int -> Int -> Int -> [[Int]]
caminos g a b = aux [[b]] where 
    aux [] = []
    aux ((x:xs):yss)
        | x == a    = (x:xs) : aux yss
        | otherwise = aux ([z:x:xs | z <- adyacentes g x
                                   , z `notElem` (x:xs)] 
                           ++ yss) 

-- 2ª solución (mediante espacio de estados)
-- =========================================

caminos2 :: Grafo Int Int -> Int -> Int -> [[Int]]
caminos2 g a b = buscaEE sucesores esFinal inicial
    where inicial          = [b]
          sucesores (x:xs) = [z:x:xs | z <- adyacentes g x
                                     , z `notElem` (x:xs)] 
          esFinal (x:xs)   = x == a

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

--    ghci> length (caminos (grafo [(i,j) | i <- [1..10], j <- [i..10]]) 1 10)
--    109601
--    (3.57 secs, 500533816 bytes)
--    ghci> length (caminos2 (grafo [(i,j) | i <- [1..10], j <- [i..10]]) 1 10)
--    109601
--    (3.53 secs, 470814096 bytes)

import Data.List (sort)

import I1M.Grafo

import I1M.BusquedaEnEspaciosDeEstados

grafo :: [(Int,Int)] -> Grafo Int Int

grafo as = creaGrafo ND (m,n) [(x,y,0) | (x,y) <- as]

where ns = map fst as ++ map snd as

m = minimum ns

n = maximum ns

-- 1ª solución

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

caminos :: Grafo Int Int -> Int -> Int -> [[Int]]

caminos g a b = aux [[b]] where

aux [] = []

aux ((x:xs):yss)

| x == a = (x:xs) : aux yss

| otherwise = aux ([z:x:xs | z <- adyacentes g x

, z `notElem` (x:xs)]

++ yss)

-- 2ª solución (mediante espacio de estados)

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

caminos2 :: Grafo Int Int -> Int -> Int -> [[Int]]

caminos2 g a b = buscaEE sucesores esFinal inicial

where inicial = [b]

sucesores (x:xs) = [z:x:xs | z <- adyacentes g x

, z `notElem` (x:xs)]

esFinal (x:xs) = x == a

-- Comparación de eficiencia

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

-- ghci> length (caminos (grafo [(i,j) | i <- [1..10], j <- [i..10]]) 1 10)

-- 109601

-- (3.57 secs, 500533816 bytes)

-- ghci> length (caminos2 (grafo [(i,j) | i <- [1..10], j <- [i..10]]) 1 10)

-- 109601

-- (3.53 secs, 470814096 bytes)

Problema de las bolas de Dijkstra

PorJosé A. Alonso 1-junio-201513-junio-2015

En el juego de las bolas de Dijkstra se dispone de una bolsa con bolas blancas y negras. El juego consiste en elegir al azar dos bolas de la bolsa y añadir una bola negra si las dos bolas elegidas son del mismo color o una bola blanca en caso contrario. El juego termina cuando queda sólo una bola en la bolsa.

Vamos a representar las bolas blancas por 0, las negras por 1 y la bolsa la representaremos por una lista cuyos elementos son 0 ó 1.

Definir las funciones

   juego  :: [Int] -> [[Int]]
   ultima :: [Int] -> Int

1 2	juego :: [Int] -> [[Int]] ultima :: [Int] -> Int

tales que

(juego xs) es la lista de los pasos aleatorios de un juego de Dijkstra a partir de la lista xs. Por ejemplo,

     juego [1,1,0,0,1]  ==  [[1,1,0,0,1],[1,1,0,0],[1,1,1],[1,1],[1]]
     juego [1,1,0,0,1]  ==  [[1,1,0,0,1],[0,1,1,0],[0,0,1],[1,1],[1]]
     juego [1,0,0,0,1]  ==  [[1,0,0,0,1],[0,0,0,1],[0,1,1],[1,0],[0]]
     juego [1,0,1,1,1]  ==  [[1,0,1,1,1],[1,1,0,1],[1,0,1],[0,1],[0]]

juego [1,1,0,0,1] == [[1,1,0,0,1],[1,1,0,0],[1,1,1],[1,1],[1]]

juego [1,1,0,0,1] == [[1,1,0,0,1],[0,1,1,0],[0,0,1],[1,1],[1]]

juego [1,0,0,0,1] == [[1,0,0,0,1],[0,0,0,1],[0,1,1],[1,0],[0]]

juego [1,0,1,1,1] == [[1,0,1,1,1],[1,1,0,1],[1,0,1],[0,1],[0]]

(ultima xs) es la bola que queda en la bolsa al final del juego de Dijkstra a partir de xs. Por ejemplo,

     ultima [1,1,0,0,1]  ==  1
     ultima [1,0,0,0,1]  ==  0
     ultima [1,0,1,1,1]  ==  0

ultima [1,1,0,0,1] == 1

ultima [1,0,0,0,1] == 0

ultima [1,0,1,1,1] == 0

Comprobar con QuickCheck que la bola que queda en la bolsa al final del juego de Dijkstra es blanca si, y sólo si, el número de bolas blancas en la bolsa inicial es impar.

Soluciones

import System.IO.Unsafe
import System.Random
import Test.QuickCheck

-- (aleatorio a b) es un número aleatorio entre a y b. Por ejemplo, 
--    ghci> aleatorio 0 1000
--    681
--    ghci> aleatorio 0 1000
--    66
aleatorio :: Random t => t -> t -> t
aleatorio a b = unsafePerformIO $ 
                getStdRandom (randomR (a,b))

-- (aleatorios m n) es una lista infinita de números aleatorios entre m y
-- n. Por ejemplo, 
--    ghci> take 20 (aleatorios 2 9)
--    [6,5,3,9,6,3,6,6,2,7,9,6,8,6,2,4,2,6,9,4]
--    ghci> take 20 (aleatorios 2 9)
--    [3,7,7,5,7,7,5,8,6,4,7,2,8,8,2,8,7,6,5,5]
aleatorios :: Random t => t -> t -> [t]
aleatorios m n = aleatorio m n : aleatorios m n

-- (selecciona n x) es el par formado por el n-ésimo elemento de xs y
-- los restantes elementos. Por ejemplo,
--    selecciona 0 "abc"  ==  ('a',"bc")
--    selecciona 1 "abc"  ==  ('b',"ac")
--    selecciona 2 "abc"  ==  ('c',"ab")
selecciona :: Int -> [a] -> (a,[a])
selecciona n xs = (z,ys++zs)
    where (ys,z:zs) = splitAt n xs

-- (extraeAleatorio xs) es el par formado aleatoriamente por un elemento
-- de xs y las restantes. Por ejemplo, 
--    extraeAleatorio [1,0,0,1,1]  ==  (0,[1,0,1,1])
--    extraeAleatorio [1,0,0,1,1]  ==  (1,[0,0,1,1])
--    extraeAleatorio [1,0,0,1,1]  ==  (0,[1,0,1,1])
--    extraeAleatorio [1,0,0,1,1]  ==  (1,[1,0,0,1])
extraeAleatorio :: [a] -> (a,[a])
extraeAleatorio xs = selecciona n xs
    where n = aleatorio 0 (length xs - 1)

-- (inserta n x xs) inserta, en la posición n, el elemento x en la lista
-- xs. Por ejemplo,  
--    inserta 0 'a' "bcd"  ==  "abcd"
--    inserta 1 'a' "bcd"  ==  "bacd"
--    inserta 2 'a' "bcd"  ==  "bcad"
--    inserta 3 'a' "bcd"  ==  "bcda"
inserta :: Int -> a -> [a] -> [a]
inserta n x xs = us ++ (x:vs)
    where (us,vs) = splitAt n xs

-- (insertaAleatorio x xs) inserta aleatoriamente el elemento x en la
-- lista xs. Por ejemplo,
--    insertaAleatorio 9 [0..8]  ==  [0,1,2,3,4,9,5,6,7,8]
--    insertaAleatorio 9 [0..8]  ==  [0,1,9,2,3,4,5,6,7,8]
--    insertaAleatorio 9 [0..8]  ==  [0,1,2,3,4,5,6,9,7,8]
insertaAleatorio :: a -> [a] -> [a]
insertaAleatorio x xs = inserta n x xs
    where n = aleatorio 0 (length xs - 1)

-- (paso xs) elige aleatoriamente dos bolas de xs e inserta
-- aleatoriamente una bola negra si las dos extraídas son del mismo
-- color o una blanca, en caso contrario. Por ejemplo, 
--    paso [1,0,1,1,0,0,1,1,1,0,0,0]  ==  [0,1,1,0,0,1,1,0,1,0,0]
--    paso [1,0,1,1,0,0,1,1,1,0,0,0]  ==  [1,0,0,0,1,1,1,0,0,1,0]
paso :: [Int] -> [Int]
paso [x] = []
paso xs | x == y    = insertaAleatorio 1 zs
        | otherwise = insertaAleatorio 0 zs 
    where (x,ys) = extraeAleatorio xs
          (y,zs) = extraeAleatorio ys

-- (juego xs) es la lista de los pasos aleatorios de un juego a partir
-- de la lista xs. Por ejemplo,
--    juego [1,1,0,0,1]  ==  [[1,1,0,0,1],[1,1,0,0],[1,1,1],[1,1],[1]]
--    juego [1,1,0,0,1]  ==  [[1,1,0,0,1],[0,1,1,0],[0,0,1],[1,1],[1]]
--    juego [1,1,0,0,0]  ==  [[1,1,0,0,0],[0,1,0,0],[1,1,0],[0,1],[0]]
juego :: [Int] -> [[Int]]
juego xs = takeWhile (not . null) (iterate paso xs)

-- (ultima xs) es la bola que queda en la bolsa al final del juego de
-- Dijkstra. Por ejemplo, 
--    ultima [1,0,0,1]  ==  1
ultima :: [Int] -> Int
ultima = head . last . juego

-- (numeroDeCeros xs) es el número de ceros en xs. Por ejemplo,
--    numeroDeCeros [1,0,0,1,0]  ==  3
numeroDeCeros :: [Int] -> Int
numeroDeCeros xs = length (filter (==0) xs)

-- Propiedad. La bola que queda en la bolsa al final del juego de
-- Dijkstra es blanca si, y sólo si, el número de bolas blancas en la
-- bolsa inicial es impar. 
prop_Dijkstra :: [Int] -> Property
prop_Dijkstra xs = 
    not (null xs) ==> (ultima ys == 0) == odd (numeroDeCeros ys)
    where ys = [x `mod` 2 | x <- xs]

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

100

101

102

import System.IO.Unsafe

import System.Random

import Test.QuickCheck

-- (aleatorio a b) es un número aleatorio entre a y b. Por ejemplo,

-- ghci> aleatorio 0 1000

-- 681

-- ghci> aleatorio 0 1000

-- 66

aleatorio :: Random t => t -> t -> t

aleatorio a b = unsafePerformIO $

getStdRandom (randomR (a,b))

-- (aleatorios m n) es una lista infinita de números aleatorios entre m y

-- n. Por ejemplo,

-- ghci> take 20 (aleatorios 2 9)

-- [6,5,3,9,6,3,6,6,2,7,9,6,8,6,2,4,2,6,9,4]

-- ghci> take 20 (aleatorios 2 9)

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

aleatorios :: Random t => t -> t -> [t]

aleatorios m n = aleatorio m n : aleatorios m n

-- (selecciona n x) es el par formado por el n-ésimo elemento de xs y

-- los restantes elementos. Por ejemplo,

-- selecciona 0 "abc" == ('a',"bc")

-- selecciona 1 "abc" == ('b',"ac")

-- selecciona 2 "abc" == ('c',"ab")

selecciona :: Int -> [a] -> (a,[a])

selecciona n xs = (z,ys++zs)

where (ys,z:zs) = splitAt n xs

-- (extraeAleatorio xs) es el par formado aleatoriamente por un elemento

-- de xs y las restantes. Por ejemplo,

-- extraeAleatorio [1,0,0,1,1] == (0,[1,0,1,1])

-- extraeAleatorio [1,0,0,1,1] == (1,[0,0,1,1])

-- extraeAleatorio [1,0,0,1,1] == (0,[1,0,1,1])

-- extraeAleatorio [1,0,0,1,1] == (1,[1,0,0,1])

extraeAleatorio :: [a] -> (a,[a])

extraeAleatorio xs = selecciona n xs

where n = aleatorio 0 (length xs - 1)

-- (inserta n x xs) inserta, en la posición n, el elemento x en la lista

-- xs. Por ejemplo,

-- inserta 0 'a' "bcd" == "abcd"

-- inserta 1 'a' "bcd" == "bacd"

-- inserta 2 'a' "bcd" == "bcad"

-- inserta 3 'a' "bcd" == "bcda"

inserta :: Int -> a -> [a] -> [a]

inserta n x xs = us ++ (x:vs)

where (us,vs) = splitAt n xs

-- (insertaAleatorio x xs) inserta aleatoriamente el elemento x en la

-- lista xs. Por ejemplo,

-- insertaAleatorio 9 [0..8] == [0,1,2,3,4,9,5,6,7,8]

-- insertaAleatorio 9 [0..8] == [0,1,9,2,3,4,5,6,7,8]

-- insertaAleatorio 9 [0..8] == [0,1,2,3,4,5,6,9,7,8]

insertaAleatorio :: a -> [a] -> [a]

insertaAleatorio x xs = inserta n x xs

where n = aleatorio 0 (length xs - 1)

-- (paso xs) elige aleatoriamente dos bolas de xs e inserta

-- aleatoriamente una bola negra si las dos extraídas son del mismo

-- color o una blanca, en caso contrario. Por ejemplo,

-- paso [1,0,1,1,0,0,1,1,1,0,0,0] == [0,1,1,0,0,1,1,0,1,0,0]

-- paso [1,0,1,1,0,0,1,1,1,0,0,0] == [1,0,0,0,1,1,1,0,0,1,0]

paso :: [Int] -> [Int]

paso [x] = []

paso xs | x == y = insertaAleatorio 1 zs

| otherwise = insertaAleatorio 0 zs

where (x,ys) = extraeAleatorio xs

(y,zs) = extraeAleatorio ys

-- (juego xs) es la lista de los pasos aleatorios de un juego a partir

-- de la lista xs. Por ejemplo,

-- juego [1,1,0,0,1] == [[1,1,0,0,1],[1,1,0,0],[1,1,1],[1,1],[1]]

-- juego [1,1,0,0,1] == [[1,1,0,0,1],[0,1,1,0],[0,0,1],[1,1],[1]]

-- juego [1,1,0,0,0] == [[1,1,0,0,0],[0,1,0,0],[1,1,0],[0,1],[0]]

juego :: [Int] -> [[Int]]

juego xs = takeWhile (not . null) (iterate paso xs)

-- (ultima xs) es la bola que queda en la bolsa al final del juego de

-- Dijkstra. Por ejemplo,

-- ultima [1,0,0,1] == 1

ultima :: [Int] -> Int

ultima = head . last . juego

-- (numeroDeCeros xs) es el número de ceros en xs. Por ejemplo,

-- numeroDeCeros [1,0,0,1,0] == 3

numeroDeCeros :: [Int] -> Int

numeroDeCeros xs = length (filter (==0) xs)

-- Propiedad. La bola que queda en la bolsa al final del juego de

-- Dijkstra es blanca si, y sólo si, el número de bolas blancas en la

-- bolsa inicial es impar.

prop_Dijkstra :: [Int] -> Property

prop_Dijkstra xs =

not (null xs) ==> (ultima ys == 0) == odd (numeroDeCeros ys)

where ys = [x `mod` 2 | x <- xs]

-- La comprobación es

-- ghci> quickCheck prop_Dijkstra

-- +++ OK, passed 100 tests.

Meses: junio 2015

Pandigitales tridivisibles

Soluciones

Partición por longitudes

Soluciones

Máximos de una lista

Soluciones

Menor n tal que el primo n-ésimo cumple una propiedad

Soluciones

Polinomio cromático de un grafo

Soluciones

Suma de posteriores

Soluciones

Algoritmo de bajada para resolver un sistema triangular inferior

Soluciones

Mayor producto de n números adyacentes en una matriz

Soluciones

Método de bisección para encontrar raíces de funciones

Soluciones

Descomposiciones en sumas de primos

Soluciones

Caminos en un grafo

Soluciones

Problema de las bolas de Dijkstra

Soluciones

Entradas recientes

Buscar

Correo electrónico