Menu Close

Etiqueta: product

Matrices de Pascal

El triángulo de Pascal es un triángulo de números

         1
        1 1
       1 2 1
     1  3 3  1
    1 4  6  4 1
   1 5 10 10 5 1
  ...............

construido de la siguiente forma

  • la primera fila está formada por el número 1;
  • las filas siguientes se construyen sumando los números adyacentes de la fila superior y añadiendo un 1 al principio y al final de la fila.

La matriz de Pascal es la matriz cuyas filas son los elementos de la
correspondiente fila del triángulo de Pascal completadas con ceros. Por ejemplo, la matriz de Pascal de orden 6 es

   |1 0  0  0 0 0|
   |1 1  0  0 0 0|
   |1 2  1  0 0 0|
   |1 3  3  1 0 0|
   |1 4  6  4 1 0|
   |1 5 10 10 5 1|

Definir la función

   matrizPascal :: Int -> Matrix Integer

tal que (matrizPascal n) es la matriz de Pascal de orden n. Por ejemplo,

   λ> matrizPascal 6
   (  1  0  0  0  0  0 )
   (  1  1  0  0  0  0 )
   (  1  2  1  0  0  0 )
   (  1  3  3  1  0  0 )
   (  1  4  6  4  1  0 )
   (  1  5 10 10  5  1 )

Soluciones

import Data.Matrix
 
-- 1ª solución
-- ===========
 
matrizPascal :: Int -> Matrix Integer 
matrizPascal 1 = fromList 1 1 [1]
matrizPascal n = matrix n n f 
  where f (i,j) | i < n && j <  n  = p!(i,j)
                | i < n && j == n  = 0
                | j == 1 || j == n = 1
                | otherwise        = p!(i-1,j-1) + p!(i-1,j)
        p = matrizPascal (n-1)
 
-- 2ª solución
-- ===========
 
matrizPascal2 :: Int -> Matrix Integer
matrizPascal2 n = fromLists xss
  where yss = take n pascal
        xss = map (take n) (map (++ repeat 0) yss)
 
pascal :: [[Integer]]
pascal = [1] : map f pascal
    where f xs = zipWith (+) (0:xs) (xs ++ [0])
 
-- 3ª solución
-- ===========
 
matrizPascal3 :: Int -> Matrix Integer
matrizPascal3 n =  matrix n n f
  where f (i,j) | i >=  j   = comb (i-1) (j-1)
                | otherwise = 0
 
-- (comb n k) es el número de combinaciones (o coeficiente binomial) de
-- n sobre k. Por ejemplo,
comb :: Int -> Int -> Integer
comb n k = product [n',n'-1..n'-k'+1] `div` product [1..k']
  where n' = fromIntegral n
        k' = fromIntegral k
 
-- 4ª solución
-- ===========
 
matrizPascal4 :: Int -> Matrix Integer
matrizPascal4 n = p
  where p = matrix n n (\(i,j) -> f i j)
        f i 1 = 1
        f i j
          | j >  i    = 0
          | i == j    = 1
          | otherwise = p!(i-1,j) + p!(i-1,j-1)
 
-- Comparación de eficiencia
-- =========================
 
--    λ> maximum (matrizPascal 150)
--    46413034868354394849492907436302560970058760
--    (2.58 secs, 394,030,504 bytes)
--    λ> maximum (matrizPascal2 150)
--    46413034868354394849492907436302560970058760
--    (0.03 secs, 8,326,784 bytes)
--    λ> maximum (matrizPascal3 150)
--    46413034868354394849492907436302560970058760
--    (0.38 secs, 250,072,360 bytes)
--    λ> maximum (matrizPascal4 150)
--    46413034868354394849492907436302560970058760
--    (0.10 secs, 13,356,360 bytes)
--    
--    λ> length (show (maximum (matrizPascal2 300)))
--    89
--    (0.06 secs, 27,286,296 bytes)
--    λ> length (show (maximum (matrizPascal3 300)))
--    89
--    (2.74 secs, 2,367,037,536 bytes)
--    λ> length (show (maximum (matrizPascal4 300)))
--    89
--    (0.36 secs, 53,934,792 bytes)
--    
--    λ> length (show (maximum (matrizPascal2 700)))
--    209
--    (0.83 secs, 207,241,080 bytes)
--    λ> length (show (maximum (matrizPascal4 700)))
--    209
--    (2.22 secs, 311,413,008 bytes)

Números cuyos factoriales son divisibles por x pero no por y

Hay 3 números (el 2, 3 y 4) cuyos factoriales son divisibles por 2 pero no por 5. Análogamente, hay números 5 (el 5, 6, 7, 8, 9) cuyos factoriales son divisibles por 15 pero no por 25.

Definir la función

   nNumerosConFactorialesDivisibles :: Integer -> Integer -> Integer

tal que (nNumerosConFactorialesDivisibles x y) es la cantidad de números cuyo factorial es divisible por x pero no por y. Por ejemplo,

  nNumerosConFactorialesDivisibles 2   5     ==  3
  nNumerosConFactorialesDivisibles 15  25    ==  5
  nNumerosConFactorialesDivisibles 100 2000  ==  5

Soluciones

import Data.List (genericLength)
 
-- 1ª solución
-- ===========
 
nNumerosConFactorialesDivisibles :: Integer -> Integer -> Integer
nNumerosConFactorialesDivisibles x y =
  genericLength (numerosConFactorialesDivisibles x y)
 
-- (numerosConFactorialesDivisibles x y) es la lista de números
-- divisibles por el factorial de x pero no divisibles por el 
-- factorial de y. Por ejemplo,
--   numerosConFactorialesDivisibles 2  5   ==  [2,3,4]
--   numerosConFactorialesDivisibles 15 25  ==  [5,6,7,8,9]
numerosConFactorialesDivisibles :: Integer -> Integer -> [Integer]
numerosConFactorialesDivisibles x y =
  [z | z <- [0..y-1]
     , factorial z `mod` x == 0
     , factorial z `mod` y /= 0]
 
-- (factorial n) es el factorial de n. Por ejemplo, 
--   factorial 4  ==  24
factorial :: Integer -> Integer
factorial n = product [1..n]
 
-- 2ª solución (usando la función de Smarandache)
-- ==============================================
 
nNumerosConFactorialesDivisibles2 :: Integer -> Integer -> Integer
nNumerosConFactorialesDivisibles2 x y =
  max 0 (smarandache y - smarandache x)
 
--(smarandache n) es el menor número cuyo factorial es divisible por
-- n. Por ejemplo,   
--    smarandache 8   ==  4
--    smarandache 10  ==  5
--    smarandache 16  ==  6
smarandache :: Integer -> Integer
smarandache x =
  head [n | (n,y) <- zip [0..] factoriales
          , y `mod` x == 0]
 
-- factoriales es la lista de los factoriales. Por ejemplo, 
--    λ> take 12 factoriales
--    [1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800]
factoriales :: [Integer]
factoriales = 1 : scanl1 (*) [1..]
 
-- Comparación de eficiencia
--    λ> nNumerosConFactorialesDivisibles 100 2000
--    5
--    (2.70 secs, 3,933,938,648 bytes)
--    λ> nNumerosConFactorialesDivisibles2 100 2000
--    5
--    (0.01 secs, 148,200 bytes)

Factorial módulo

Definir la función

   factorialMod :: Integer -> Integer -> Integer

tal que (factorialMod n x) es el factorial de x módulo n. Por ejemplo,

   factorialMod (7+10^9) 100       ==  437918130
   factorialMod (7+10^9) (5*10^6)  ==  974067448

Soluciones

import Data.List (foldl')
 
-- 1ª definición
-- =============
 
factorialMod :: Integer -> Integer -> Integer
factorialMod n x =
  factorial x `mod` n
 
-- (factorial x) es el factorial de x. Por ejemplo,
--    factorial 3  ==  6
factorial :: Integer -> Integer
factorial x = product [1..x]
 
-- 2ª definición
-- =============
 
factorialMod2 :: Integer -> Integer -> Integer
factorialMod2 n x =
  foldl' (prodMod n) 1 [1..x]
 
-- (prodMod n x y) es el producto de x e y módulo n. Por ejemplo,
--    prodMod 3 4 7  ==  1
--    prodMod 5 4 7  ==  -- 3
prodMod :: Integer -> Integer -> Integer -> Integer
prodMod n x y =
  ((x `mod` n) * (y `mod` n)) `mod` n
 
-- Comparación de eficiencia
-- =========================
 
--    λ> factorialMod (7+10^9) (5*10^4)
--    737935835
--    (2.62 secs, 2,601,358,640 bytes)
--    λ> factorialMod2 (7+10^9) (5*10^4)
--    737935835
--    (0.07 secs, 23,471,880 bytes)

Aplicaciones biyectivas

Definir las funciones

   biyectivas  :: (Ord a, Ord b) => [a] -> [b] -> [[(a,b)]]
   nBiyectivas :: (Ord a, Ord b) => [a] -> [b] -> Integer

tales que

  • (biyectivas xs ys) es el conjunto de las aplicaciones biyectivas del conjunto xs en el conjunto ys. Por ejemplo,
     λ> biyectivas [1,3] [2,4]
     [[(1,2),(3,4)],[(1,4),(3,2)]]
     λ> biyectivas [1,3,5] [2,4,6]
     [[(1,2),(3,4),(5,6)],[(1,4),(3,2),(5,6)],[(1,6),(3,4),(5,2)],
      [(1,4),(3,6),(5,2)],[(1,6),(3,2),(5,4)],[(1,2),(3,6),(5,4)]]
     λ> biyectivas [1,3] [2,4,6]
     []
  • (nBiyectivas xs ys) es el número de aplicaciones biyectivas del conjunto xs en el conjunto ys. Por ejemplo,
     nBiyectivas [1,3] [2,4]      ==  2
     nBiyectivas [1,3,5] [2,4,6]  ==  6
     λ> nBiyectivas [1,3] [2,4,6] ==  0
     length (show (nBiyectivas2 [1..2*10^4] [1..2*10^4]))  ==  77338

Nota: En este ejercicio los conjuntos se representan mediante listas ordenadas de elementos distintos.

Soluciones

import Data.List (genericLength, permutations)
 
-- 1ª definición de biyectivas
biyectivas :: (Ord a, Ord b) => [a] -> [b] -> [[(a,b)]]
biyectivas xs ys
  | length xs == length ys = [zip xs zs | zs <- permutations ys]
  | otherwise              = []
 
-- 1ª definición de nBiyectivas
nBiyectivas :: (Ord a, Ord b) => [a] -> [b] -> Integer
nBiyectivas xs ys
  | length xs == length ys = genericLength (biyectivas xs ys)
  | otherwise              = 0
 
-- 2ª definición de nBiyectivas
nBiyectivas2 :: (Ord a, Ord b) => [a] -> [b] -> Integer
nBiyectivas2 xs ys 
  | n == m    = product [1..n]
  | otherwise = 0
  where n = genericLength xs
        m = genericLength ys

Cadenas de sumas de factoriales de los dígitos

Dado un número n se considera la sucesión cuyo primer término es n y los restantes se obtienen sumando los factoriales de los dígitos del anterior. Por ejemplo, la sucesión que empieza en 69 es

         69
     363600  (porque 6! + 9! = 363600)  
       1454  (porque 3! + 6! + 3! + 6! + 0! + 0! = 1454)
        169  (porque 1! + 4! + 5! + 4! = 169)
     363601  (porque 1! + 6! + 9! = 363601)
       1454  (porque 3! + 6! + 3! + 6! + 0! + 1! = 1454)
     ......

La cadena correspondiente a un número n son los términos de la sucesión que empieza en n hasta la primera repetición de un elemento en la sucesión. Por ejemplo, la cadena de 69 es

   [69,363600,1454,169,363601,1454]

Consta de una parte no periódica ([69,363600]) y de una periódica ([1454,169,363601]).

Definir las funciones

   cadena  :: Integer -> [Integer]
   periodo :: Integer -> [Integer]

tales que

  • (cadena n es la cadena correspondiente al número n. Por ejemplo,
      cadena 69    ==  [69,363600,1454,169,363601,1454]
      cadena 145   ==  [145,145]
      cadena 78    ==  [78,45360,871,45361,871]
      cadena 569   ==  [569,363720,5775,10320,11,2,2]
      cadena 3888  ==  [3888,120966,364324,782,45362,872,45362]
      maximum [length (cadena n) | n <- [1..5000]]  ==  61
      length (cadena 1479)                          ==  61
  • (periodo n) es la parte periódica de la cadena de n. Por ejemplo,
      periodo 69    ==  [169,363601,1454]
      periodo 145   ==  [145]
      periodo 78    ==  [45361,871]
      periodo 569   ==  [2]
      periodo 3888  ==  [872,45362]
      maximum [length (periodo n) | n <- [1..5000]]  ==  3
      length (periodo 1479)                          ==  3

Soluciones

-- Definición de cadena
-- ====================
 
cadena :: Integer -> [Integer]
cadena n = reverse (extension [n])
 
extension :: [Integer] -> [Integer]
extension (n:ns)
  | m `elem` (n:ns) = m : n : ns
  | otherwise       = extension (m : n : ns)
  where m = siguiente n
 
siguiente :: Integer -> Integer
siguiente n =
  sum [factorial d | d <- digitos n]
 
factorial :: Integer -> Integer
factorial n =
  product [1..n]
 
digitos :: Integer -> [Integer]
digitos n =
  [read [c] | c <- show n]
 
-- Definición de periodo
-- =====================
 
periodo :: Integer -> [Integer]
periodo n =
  reverse (x : takeWhile (/= x) xs)
  where (x:xs) = reverse (cadena n)

Pares definidos por su MCD y su MCM

Definir las siguientes funciones

   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
  • (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

Soluciones

import Data.Numbers.Primes (primeFactors)
import Data.List (genericLength, group, nub, sort, subsequences)
import Test.QuickCheck
 
-- 1ª definición de pares
-- ======================
 
pares1 :: Integer -> Integer -> [(Integer,Integer)]
pares1 a b = [(x,y) | x <- [1..b]
                    , y <- [1..b]
                    , gcd x y == a
                    , lcm x y == b]
 
-- 2ª definición de pares
-- ======================
 
pares2 :: Integer -> Integer -> [(Integer,Integer)]
pares2 a b = [(x,y) | x <- [a,a+a..b]
                    , y <- [a,a+a..b]
                    , gcd x y == a
                    , lcm x y == b]
 
-- Comparacioń de eficiencia
--    λ> length (pares1 3 (product [3,5..11]))
--    16
--    (95.12 secs, 86,534,165,528 bytes)
--    λ> length (pares2 3 (product [3,5..11]))
--    16
--    (15.80 secs, 14,808,762,128 bytes)
 
-- 3ª definición de pares
-- ======================
 
pares3 :: Integer -> Integer -> [(Integer,Integer)]
pares3 a b = [(x,y) | x <- [a,a+a..b]
                    , c `rem` x == 0
                    , let y = c `div` x
                    , gcd x y == a
                    ]
  where c = a * b
 
-- Comparacioń de eficiencia
--    λ> length (pares2 3 (product [3,5..11]))
--    16
--    (15.80 secs, 14,808,762,128 bytes)
--    λ> length (pares3 3 (product [3,5..11]))
--    16
--    (0.01 secs, 878,104 bytes)
 
-- 4ª definición de pares
-- ======================
 
-- Para la cuarta definición de pares se observa la relación con los
-- factores primos
--    λ> [(primeFactors x, primeFactors y) | (x,y) <- pares1 2 12]
--    [([2],[2,2,3]),([2,2],[2,3]),([2,3],[2,2]),([2,2,3],[2])]
--    λ> [primeFactors x | (x,y) <- pares1 2 12]
--    [[2],[2,2],[2,3],[2,2,3]]
--    λ> [primeFactors x | (x,y) <- pares1 2 60]
--    [[2],[2,2],[2,3],[2,5],[2,2,3],[2,2,5],[2,3,5],[2,2,3,5]]
--    λ> [primeFactors x | (x,y) <- pares1 6 60]
--    [[2,3],[2,2,3],[2,3,5],[2,2,3,5]]
--    λ> [primeFactors x | (x,y) <- pares1 2 24]
--    [[2],[2,3],[2,2,2],[2,2,2,3]]
-- Se observa que cada pare se obtiene de uno de los subconjuntos de los
-- divisores primos de b/a. Por ejemplo,
--    λ> (a,b) = (2,24)
--    λ> b `div` a
--    12
--    λ> primeFactors it
--    [2,2,3]
--    λ> group it
--    [[2,2],[3]]
--    λ> subsequences it
--    [[],[[2,2]],[[3]],[[2,2],[3]]]
--    λ> map concat it
--    [[],[2,2],[3],[2,2,3]]
--    λ> map product it
--    [1,4,3,12]
--    λ> [(a * x, b `div` x) | x <- it]
--    [(2,24),(8,6),(6,8),(24,2)]
-- A partir de la observación se construye la siguiente definición
 
pares4 :: Integer -> Integer -> [(Integer,Integer)]
pares4 a b
  | b `mod` a /= 0 = []
  | otherwise =
    [(a * x, b `div` x)
    | x <- map (product . concat)
               ((subsequences . group . primeFactors) (b `div` a))]
 
-- Nota. La función pares4 calcula el mismo conjunto que las anteriores,
-- pero no necesariamente en el mismo orde. Por ejemplo,
--    λ> pares3 2 60 
--    [(2,60),(4,30),(6,20),(10,12),(12,10),(20,6),(30,4),(60,2)]
--    λ> pares4 2 60 
--    [(2,60),(4,30),(6,20),(12,10),(10,12),(20,6),(30,4),(60,2)]
--    λ> pares3 2 60 == sort (pares4 2 60)
--    True
 
-- Comparacioń de eficiencia
--    λ> length (pares3 3 (product [3,5..17]))
--    64
--    (4.44 secs, 2,389,486,440 bytes)
--    λ> length (pares4 3 (product [3,5..17]))
--    64
--    (0.00 secs, 177,704 bytes)
 
-- Propiedades de equivalencia de pares
-- ====================================
 
prop_pares :: Integer -> Integer -> Property
prop_pares a b =
  a > 0 && b > 0 ==>
  all (== pares1 a b)
      [sort (f a b) | f <- [ pares2
                           , pares3
                           , pares4
                           ]]
 
prop_pares2 :: Integer -> Integer -> Property
prop_pares2 a b =
  a > 0 && b > 0 ==>
  all (== pares1 a (a * b))
      [sort (f a (a * b)) | f <- [ pares2
                                 , pares3
                                 , pares4
                                 ]]
 
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares
--    +++ OK, passed 100 tests.
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares
--    +++ OK, passed 100 tests.
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_pares2
--    +++ OK, passed 100 tests.
 
-- 1ª definición de nPares
-- =======================
 
nPares1 :: Integer -> Integer -> Integer
nPares1 a b = genericLength (pares4 a b)
 
-- 2ª definición de nPares
-- =======================
 
nPares2 :: Integer -> Integer -> Integer
nPares2 a b = 2^(length (nub (primeFactors (b `div` a))))
 
-- Comparación de eficiencia
--    λ> nPares1 3 (product [3,5..91])
--    8388608
--    (4.68 secs, 4,178,295,920 bytes)
--    λ> nPares2 3 (product [3,5..91])
--    8388608
--    (0.00 secs, 234,688 bytes)
 
-- Propiedad de equivalencia de nPares
-- ===================================
 
prop_nPares :: Integer -> Integer -> Property
prop_nPares a b =
  a > 0 && b > 0 ==>
  nPares1 a (a * b) == nPares2 a (a * b)
 
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_nPares
--    +++ OK, passed 100 tests.

Números libres de cuadrados

Un número entero positivo es libre de cuadrados si no es divisible el cuadrado de ningún entero mayor que 1. Por ejemplo, 70 es libre de cuadrado porque sólo es divisible por 1, 2, 5, 7 y 70; en cambio, 40 no es libre de cuadrados porque es divisible por 2^2.

Definir la función

   libreDeCuadrados :: Integer -> Bool

tal que (libreDeCuadrados x) se verifica si x es libre de cuadrados. Por ejemplo,

   libreDeCuadrados 70                    ==  True
   libreDeCuadrados 40                    ==  False
   libreDeCuadrados 510510                ==  True
   libreDeCuadrados (((10^10)^10)^10)     ==  False

Otro ejemplo,

   λ> filter (not . libreDeCuadrados) [1..50]
   [4,8,9,12,16,18,20,24,25,27,28,32,36,40,44,45,48,49,50]

Soluciones

import Data.Numbers.Primes (primeFactors, primes)
import Data.List (nub)
import Test.QuickCheck
 
import Data.List (genericLength) -- Para OS
 
-- 1ª definición:
libreDeCuadrados :: Integer -> Bool
libreDeCuadrados x = x == product (divisoresPrimos x)
 
-- (divisoresPrimos x) es la lista de los divisores primos de x. Por
-- ejemplo,  
--    divisoresPrimos 40  ==  [2,5]
--    divisoresPrimos 70  ==  [2,5,7]
divisoresPrimos :: Integer -> [Integer]
divisoresPrimos x = [n | n <- divisores x, primo n]
 
-- (divisores n) es la lista de los divisores del número n. 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]
 
-- (primo n) se verifica si n es primo. Por ejemplo,
--    primo 30  == False
--    primo 31  == True  
primo :: Integer -> Bool
primo n = divisores n == [1, n]
 
-- 2ª definición
libreDeCuadrados2 :: Integer -> Bool
libreDeCuadrados2 n = 
    null [x | x <- [2..n], rem n (x^2) == 0]
 
-- 3ª definición
libreDeCuadrados3 :: Integer -> Bool
libreDeCuadrados3 n = 
    null [x | x <- [2..floor (sqrt (fromIntegral n))], 
              rem n (x^2) == 0]
 
-- 4ª definición
libreDeCuadrados4 :: Integer -> Bool
libreDeCuadrados4 n = 
  xs == nub xs
  where xs = primeFactors n
 
-- 5ª definición
libreDeCuadrados5 :: Integer -> Bool
libreDeCuadrados5 n = 
  all (\(x,y) -> x /= y) (zip xs (tail xs))
  where xs = primeFactors n
 
-- Equivalencia de las definiciones
-- ================================
 
prop_equivalencia :: Integer -> Property
prop_equivalencia n =
  n > 0 ==>
  all (== libreDeCuadrados n)
      [f n | f <- [ libreDeCuadrados2
                  , libreDeCuadrados3
                  , libreDeCuadrados4
                  , libreDeCuadrados5
                  ]]
 
-- La comprobación es
--    λ> quickCheck prop_equivalencia
--    +++ OK, passed 100 tests.
 
-- Comparación de eficiencia
-- =========================
 
--    λ> libreDeCuadrados 510510
--    True
--    (0.76 secs, 89,522,360 bytes)
--    λ> libreDeCuadrados2 510510
--    True
--    (1.78 secs, 371,826,320 bytes)
--    λ> libreDeCuadrados3 510510
--    True
--    (0.01 secs, 0 bytes)
--    λ> libreDeCuadrados4 510510
--    True
--    (0.00 secs, 152,712 bytes)
--
--    λ> filter libreDeCuadrados3 [1..] !! (2*10^4)
--    32906
--    (2.24 secs, 1,812,139,456 bytes)
--    λ> filter libreDeCuadrados4 [1..] !! (2*10^4)
--    32906
--    (0.51 secs, 936,216,664 bytes)
--    λ> filter libreDeCuadrados5 [1..] !! (2*10^4)
--    32906
--    (0.38 secs, 806,833,952 bytes)
--
--    λ> filter libreDeCuadrados4 [1..] !! (10^5)
--    164499
--    (3.28 secs, 7,470,629,888 bytes)
--    λ> filter libreDeCuadrados5 [1..] !! (10^5)
--    164499
--    (2.88 secs, 6,390,072,384 bytes)

El problema de las N torres

El problema de las N torres consiste en colocar N torres en un tablero con N filas y N columnas de forma que no haya dos torres en la misma fila ni en la misma columna.

Cada solución del problema de puede representar mediante una matriz con ceros y unos donde los unos representan las posiciones ocupadas por las torres y los ceros las posiciones libres. Por ejemplo,

   ( 0 1 0 )
   ( 1 0 0 )
   ( 0 0 1 )

representa una solución del problema de las 3 torres.

Definir las funciones

   torres  :: Int -> [Matrix Int]
   nTorres :: Int -> Integer

tales que
+ (torres n) es la lista de las soluciones del problema de las n torres. Por ejemplo,

      λ> torres 3
      [( 1 0 0 )
       ( 0 1 0 )
       ( 0 0 1 )
      ,( 1 0 0 )
       ( 0 0 1 )
       ( 0 1 0 )
      ,( 0 1 0 )
       ( 1 0 0 )
       ( 0 0 1 )
      ,( 0 1 0 )
       ( 0 0 1 )
       ( 1 0 0 )
      ,( 0 0 1 )
       ( 1 0 0 )
       ( 0 1 0 )
      ,( 0 0 1 )
       ( 0 1 0 )
       ( 1 0 0 )
      ]
  • (nTorres n) es el número de soluciones del problema de las n torres. Por ejemplo,
      λ> nTorres 3
      6
      λ> length (show (nTorres (10^4)))
      35660

Soluciones

import Data.List (genericLength, sort, permutations)
import Data.Matrix 
 
-- 1ª definición de torres
-- =======================
 
torres1 :: Int -> [Matrix Int]
torres1 n = 
    [permutacionAmatriz n p | p <- sort (permutations [1..n])]
 
permutacionAmatriz :: Int -> [Int] -> Matrix Int
permutacionAmatriz n p =
    matrix n n f
    where f (i,j) | (i,j) `elem` posiciones = 1
                  | otherwise               = 0
          posiciones = zip [1..n] p    
 
-- 2ª definición de torres
-- =======================
 
torres2 :: Int -> [Matrix Int]
torres2 = map fromLists . permutations . toLists . identity
 
-- El cálculo con la definición anterior es:
--    λ> identity 3
--    ( 1 0 0 )
--    ( 0 1 0 )
--    ( 0 0 1 )
--    
--    λ> toLists it
--    [[1,0,0],[0,1,0],[0,0,1]]
--    λ> permutations it
--    [[[1,0,0],[0,1,0],[0,0,1]],
--     [[0,1,0],[1,0,0],[0,0,1]],
--     [[0,0,1],[0,1,0],[1,0,0]],
--     [[0,1,0],[0,0,1],[1,0,0]],
--     [[0,0,1],[1,0,0],[0,1,0]],
--     [[1,0,0],[0,0,1],[0,1,0]]]
--    λ> map fromLists it
--    [( 1 0 0 )
--     ( 0 1 0 )
--     ( 0 0 1 )
--    ,( 0 1 0 )
--     ( 1 0 0 )
--     ( 0 0 1 )
--    ,( 0 0 1 )
--     ( 0 1 0 )
--     ( 1 0 0 )
--    ,( 0 1 0 )
--     ( 0 0 1 )
--     ( 1 0 0 )
--    ,( 0 0 1 )
--     ( 1 0 0 )
--     ( 0 1 0 )
--    ,( 1 0 0 )
--     ( 0 0 1 )
--     ( 0 1 0 )
--    ]
 
-- 1ª definición de nTorres
-- ========================
 
nTorres1 :: Int -> Integer
nTorres1 = genericLength . torres1
 
-- 2ª definición de nTorres
-- ========================
 
nTorres2 :: Int -> Integer
nTorres2 n = product [1..fromIntegral n]
 
-- Comparación de eficiencia
-- =========================
 
--    λ> nTorres1 9
--    362880
--    (4.22 secs, 693,596,128 bytes)
--    λ> nTorres2 9
--    362880
--    (0.00 secs, 0 bytes)

Número de dígitos del factorial

Definir las funciones

   nDigitosFact :: Integer -> Integer
   graficas     :: [Integer] -> IO ()

tales que

  • (nDigitosFact n) es el número de dígitos de n!. Por ejemplo,
     nDigitosFact 0        ==  1
     nDigitosFact 4        ==  2
     nDigitosFact 5        ==  3
     nDigitosFact 10       ==  7
     nDigitosFact 100      ==  158
     nDigitosFact 1000     ==  2568
     nDigitosFact 10000    ==  35660
     nDigitosFact 100000   ==  456574
     nDigitosFact 1000000  ==  5565709
  • (graficas xs) dibuja las gráficas de los números de dígitos del factorial de k (para k en xs) y de la recta y = 5.5 x. Por ejemplo, (graficas [0,500..10^6]) dibuja
    Numero_de_digitos_del_factorial

Nota: Este ejercicio está basado en el problema How many digits? de Kattis en donde se impone la restricción de calcular, en menos de 1 segundo, el número de dígitos de los factoriales de 10.000 números del rango [0,1.000.000].

Se puede simular como sigue

   λ> import System.Random 
   λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^3]]
   λ> maximum (map nDigitosFact xs)
   5561492
   λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^3]]
   λ> maximum (map nDigitosFact xs)
   5563303

Soluciones

import Data.List (genericLength, genericIndex)
import Graphics.Gnuplot.Simple
 
-- 1ª definición
nDigitosFact1 :: Integer -> Integer
nDigitosFact1 n =
  genericLength (show (product [1..n]))
 
-- 2ª definición (usando logaritmos)
-- =================================
 
-- Basado en
--    nDigitos (n!) = 1 + floor (log (n!))
--                  = 1 + floor (log (1*2*3*...*n))
--                  = 1 + floor (log(1) + log(2) + log(3) +...+ log(n))
 
nDigitosFact2 :: Integer -> Integer
nDigitosFact2 n =
  1 + floor (sum [logBase 10 (fromIntegral k) | k <- [1..n]])
 
-- 3ª definición (usando logaritmos y programación dinámica)
-- =========================================================
 
nDigitosFact3 :: Integer -> Integer
nDigitosFact3 0 = 1
nDigitosFact3 n = 1 + floor ((sumLogs `genericIndex` (n-1)) / log 10)
 
logs :: [Double]
logs = map log [1..]
 
sumLogs :: [Double]
sumLogs = scanl1 (+) logs
 
-- 4ª definición (Usando la fórmula de Kamenetsky)
-- ===============================================
 
-- La fórmula se encuentra en https://oeis.org/A034886
 
nDigitosFact4 :: Integer -> Integer
nDigitosFact4 0 = 1
nDigitosFact4 1 = 1
nDigitosFact4 n =
  1 + floor (m * logBase 10 (m/e) + logBase 10 (2*pi*m)/2)
  where m = fromIntegral n
        e = exp 1
 
--    λ> nDigitosFact1 (4*10^4)
--    166714
--    (5.61 secs, 1,649,912,864 bytes)
--    λ> nDigitosFact2 (4*10^4)
--    166714
--    (0.10 secs, 13,741,360 bytes)
--
--    λ> nDigitosFact2 (7*10^5)
--    3787566
--    (1.86 secs, 187,666,328 bytes)
--    λ> nDigitosFact3 (7*10^5)
--    3787566
--    (0.02 secs, 0 bytes)
--    λ> nDigitosFact3 (7*10^5)
--    3787566
--    (1.01 secs, 238,682,064 bytes)
--    λ> nDigitosFact4 (7*10^5)
--    3787566
--    (0.01 secs, 0 bytes)
-- 
--    λ> import System.Random 
--    λ> xs <- sequence [randomRIO (0,10^6) | _ <- [1..10^2]]
--    λ> maximum (map nDigitosFact3 xs)
--    5565097
--    (11.19 secs, 7,336,595,440 bytes)
--    λ> maximum (map nDigitosFact4 xs)
--    5565097
--    (0.01 secs, 0 bytes)
 
graficas :: [Integer] -> IO ()
graficas xs = 
    plotLists [Key Nothing]
              [p1, p2]
    where p1, p2 :: [(Float, Float)]
          p1 = [(fi k, fi (nDigitosFact4 k)) | k <- xs]
          p2 = [(k',5.5*k') | k <- xs, let k' = fi k]
          fi = fromIntegral

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

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
Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_1
y se obtiene a partir de la serie armónica
Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_2
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 ()

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
  • (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
    Calculo_de_pi_mediante_la_variante_de_Euler_de_la_serie_armonica_3.png

Nota: Este ejercicio ha sido propuesto por Paula Macías.

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