José A. Alonso

Transitividad de una relación

PorJosé A. Alonso 13-julio-20226-julio-2022

Una relación binaria R sobre un conjunto A es transitiva cuando se cumple que siempre que un elemento se relaciona con otro y éste último con un tercero, entonces el primero se relaciona con el tercero.

Definir la función

   transitiva :: Ord a => [(a,a)] -> Bool

1	transitiva :: Ord a => [(a,a)] -> Bool

tal que (transitiva r) se verifica si la relación r es transitiva. Por ejemplo,

   transitiva [(1,1),(1,3),(3,1),(3,3),(5,5)]  ==  True
   transitiva [(1,1),(1,3),(3,1),(5,5)]        ==  False
   transitiva [(n,n) | n <- [1..10^4]]         ==  True

transitiva [(1,1),(1,3),(3,1),(3,3),(5,5)] == True

transitiva [(1,1),(1,3),(3,1),(5,5)] == False

transitiva [(n,n) | n <- [1..10^4]] == True

Soluciones

import Data.List (nub)
import Data.Maybe (mapMaybe)
import qualified Data.Map as M (Map, assocs, empty, insertWith, lookup, map)
import Test.QuickCheck (maxSize, quickCheckWith, stdArgs)

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

transitiva1 :: Ord a => [(a,a)] -> Bool
transitiva1 r = 
  and [(x,y) `elem` r | (x,u) <- r, (v,y) <- r, u == v]

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

transitiva2 :: Ord a => [(a,a)] -> Bool
transitiva2 r = 
  all (`elem` r) [(x,y) | (x,u) <- r, (v,y) <- r, u == v]

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

transitiva3 :: Ord a => [(a,a)] -> Bool
transitiva3 r =
  subconjunto (composicion r r) r

-- (subconjunto xs ys) se verifica si xs es un subconjunto de xs. Por
-- ejemplo, 
--    subconjunto [1,3] [3,1,5]  ==  True
--    subconjunto [3,1,5] [1,3]  ==  False
subconjunto :: Ord a => [a] -> [a] -> Bool
subconjunto xs ys =
  all (`elem`ys) xs

-- (composicion r s) es la composición de las relaciones binarias r y
-- s. Por ejemplo, 
--    λ> composicion [(1,2)] [(2,3),(2,4)]
--    [(1,3),(1,4)]
--    λ> composicion [(1,2),(5,2)] [(2,3),(2,4)]
--    [(1,3),(1,4),(5,3),(5,4)]
--    λ> composicion [(1,2),(1,4),(1,5)] [(2,3),(4,3)]
--    [(1,3)]
composicion :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion r s =
  nub [(x,y) | (x,u) <- r, (v,y) <- s, u == v] 

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

transitiva4 :: Ord a => [(a,a)] -> Bool
transitiva4 r =
  subconjunto (composicion4 r r) r

composicion4 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion4 r s =
  relAlista (composicionRel (listaArel r) (listaArel s))

-- Una relación se puede representar por un diccionario donde las claves
-- son los elementos y los valores son las listas de los elementos con
-- los que se relaciona.
type Rel a = M.Map a [a]

-- (listaArel xys) es la relación correspondiente a la lista de pares
-- xys. Por ejemplo.
--    λ> listaArel [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]
--    fromList [(1,[1,2,3,6]),(2,[2,6]),(3,[3,6]),(6,[6])]
listaArel :: Ord a => [(a,a)] -> Rel a
listaArel []          = M.empty
listaArel ((x,y):xys) = M.insertWith (++) x [y] (listaArel xys)

-- (composicionRel r s) es la composición de las relaciones r y s. Por
-- ejemplo,
--    λ> r = listaArel [(1,2),(5,2)]
--    λ> s = listaArel [(2,3),(2,4)]
--    λ> composicionRel r s
--    fromList [(1,[3,4]),(5,[3,4])]
composicionRel :: Ord a => Rel a -> Rel a -> Rel a
composicionRel r s =
  M.map f r 
  where f xs = concat (mapMaybe (`M.lookup` s) xs)

-- (relAlista r) es la lista de pares correspondientes a la relación
-- r. Por ejemplo,
--    λ> relAlista (M.fromList [(1,[3,4]),(5,[3,4])])
--    [(1,3),(1,4),(5,3),(5,4)]
relAlista :: Ord a => Rel a -> [(a,a)]
relAlista r =
  nub [(x,y) | (x,ys) <- M.assocs r, y <- ys] 

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

-- La propiedad es
prop_transitiva :: [(Int,Int)] -> Bool
prop_transitiva r =
  all (== transitiva1 r)
      [transitiva2 r,
       transitiva3 r,
       transitiva4 r] 

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

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

-- La comparación es
--    λ> transitiva1 [(n,n) | n <- [1..5*10^3]]
--    True
--    (4.27 secs, 1,403,139,032 bytes)
--    λ> transitiva2 [(n,n) | n <- [1..5*10^3]]
--    True
--    (4.24 secs, 1,402,739,056 bytes)
--    λ> transitiva3 [(n,n) | n <- [1..5*10^3]]
--    True
--    (4.41 secs, 1,403,219,880 bytes)
--    λ> transitiva4 [(n,n) | n <- [1..5*10^3]]
--    True
--    (0.46 secs, 16,206,624 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

import Data.List (nub)

import Data.Maybe (mapMaybe)

import qualified Data.Map as M (Map, assocs, empty, insertWith, lookup, map)

import Test.QuickCheck (maxSize, quickCheckWith, stdArgs)

-- 1ª solución

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

transitiva1 :: Ord a => [(a,a)] -> Bool

transitiva1 r =

and [(x,y) `elem` r | (x,u) <- r, (v,y) <- r, u == v]

-- 2ª solución

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

transitiva2 :: Ord a => [(a,a)] -> Bool

transitiva2 r =

all (`elem` r) [(x,y) | (x,u) <- r, (v,y) <- r, u == v]

-- 3ª solución

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

transitiva3 :: Ord a => [(a,a)] -> Bool

transitiva3 r =

subconjunto (composicion r r) r

-- (subconjunto xs ys) se verifica si xs es un subconjunto de xs. Por

-- ejemplo,

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

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

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

subconjunto xs ys =

all (`elem`ys) xs

-- (composicion r s) es la composición de las relaciones binarias r y

-- s. Por ejemplo,

-- λ> composicion [(1,2)] [(2,3),(2,4)]

-- [(1,3),(1,4)]

-- λ> composicion [(1,2),(5,2)] [(2,3),(2,4)]

-- [(1,3),(1,4),(5,3),(5,4)]

-- λ> composicion [(1,2),(1,4),(1,5)] [(2,3),(4,3)]

-- [(1,3)]

composicion :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion r s =

nub [(x,y) | (x,u) <- r, (v,y) <- s, u == v]

-- 3ª solución

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

transitiva4 :: Ord a => [(a,a)] -> Bool

transitiva4 r =

subconjunto (composicion4 r r) r

composicion4 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion4 r s =

relAlista (composicionRel (listaArel r) (listaArel s))

-- Una relación se puede representar por un diccionario donde las claves

-- son los elementos y los valores son las listas de los elementos con

-- los que se relaciona.

type Rel a = M.Map a [a]

-- (listaArel xys) es la relación correspondiente a la lista de pares

-- xys. Por ejemplo.

-- λ> listaArel [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]

-- fromList [(1,[1,2,3,6]),(2,[2,6]),(3,[3,6]),(6,[6])]

listaArel :: Ord a => [(a,a)] -> Rel a

listaArel [] = M.empty

listaArel ((x,y):xys) = M.insertWith (++) x [y] (listaArel xys)

-- (composicionRel r s) es la composición de las relaciones r y s. Por

-- ejemplo,

-- λ> r = listaArel [(1,2),(5,2)]

-- λ> s = listaArel [(2,3),(2,4)]

-- λ> composicionRel r s

-- fromList [(1,[3,4]),(5,[3,4])]

composicionRel :: Ord a => Rel a -> Rel a -> Rel a

composicionRel r s =

M.map f r

where f xs = concat (mapMaybe (`M.lookup` s) xs)

-- (relAlista r) es la lista de pares correspondientes a la relación

-- r. Por ejemplo,

-- λ> relAlista (M.fromList [(1,[3,4]),(5,[3,4])])

-- [(1,3),(1,4),(5,3),(5,4)]

relAlista :: Ord a => Rel a -> [(a,a)]

relAlista r =

nub [(x,y) | (x,ys) <- M.assocs r, y <- ys]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_transitiva :: [(Int,Int)] -> Bool

prop_transitiva r =

all (== transitiva1 r)

[transitiva2 r,

transitiva3 r,

transitiva4 r]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> transitiva1 [(n,n) | n <- [1..5*10^3]]

-- True

-- (4.27 secs, 1,403,139,032 bytes)

-- λ> transitiva2 [(n,n) | n <- [1..5*10^3]]

-- True

-- (4.24 secs, 1,402,739,056 bytes)

-- λ> transitiva3 [(n,n) | n <- [1..5*10^3]]

-- True

-- (4.41 secs, 1,403,219,880 bytes)

-- λ> transitiva4 [(n,n) | n <- [1..5*10^3]]

-- True

-- (0.46 secs, 16,206,624 bytes)

El código se encuentra en GitHub.

Composición de relaciones binarias

PorJosé A. Alonso 12-julio-20225-julio-2022

Las relaciones binarias en un conjunto A se pueden representar mediante conjuntos de pares de elementos de A. Por ejemplo, la relación de divisibilidad en el conjunto {1,2,3,6} se representa por

   [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]

1	[(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]

La composición de dos relaciones binarias R y S en el conjunto A es la relación binaria formada por los pares (x,y) para los que existe un z tal que (x,z) ∈ R y (z,y) ∈ S.

Definir la función

   composicion :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

1	composicion :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

tal que (composicion r s) es la composición de las relaciones binarias r y s. Por ejemplo,

   λ> composicion [(1,2)] [(2,3),(2,4)]
   [(1,3),(1,4)]
   λ> composicion [(1,2),(5,2)] [(2,3),(2,4)]
   [(1,3),(1,4),(5,3),(5,4)]
   λ> composicion [(1,2),(1,4),(1,5)] [(2,3),(4,3)]
   [(1,3)]

λ> composicion [(1,2)] [(2,3),(2,4)]

[(1,3),(1,4)]

λ> composicion [(1,2),(5,2)] [(2,3),(2,4)]

[(1,3),(1,4),(5,3),(5,4)]

λ> composicion [(1,2),(1,4),(1,5)] [(2,3),(4,3)]

[(1,3)]

Nota: Se supone que las relaciones binarias son listas sin elementos repetidos.

Soluciones

import Data.List (nub, sort)
import Data.Maybe (mapMaybe)
import qualified Data.Set.Monad as S (Set, fromList, toList)
import qualified Data.Map as M (Map, assocs, empty, insertWith, lookup, map)
import Test.QuickCheck (quickCheck)

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

composicion1 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion1 r s =
  nub [(x,y) | (x,u) <- r, (v,y) <- s, u == v] 

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

composicion2 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion2 r s =
  S.toList (composicionS (S.fromList r) (S.fromList s))

composicionS :: Ord a => S.Set (a,a) -> S.Set (a,a) -> S.Set (a,a)
composicionS r s =  
  [(x,y) | (x,u) <- r, (v,y) <- s, u == v] 

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

composicion3 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion3 r s =
  relAlista (composicionRel (listaArel r) (listaArel s))

-- Una relación se puede representar por un diccionario donde las claves
-- son los elementos y los valores son las listas de los elementos con
-- los que se relaciona.
type Rel a = M.Map a [a]

-- (listaArel xys) es la relación correspondiente a la lista de pares
-- xys. Por ejemplo.
--    λ> listaArel [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]
--    fromList [(1,[1,2,3,6]),(2,[2,6]),(3,[3,6]),(6,[6])]
listaArel :: Ord a => [(a,a)] -> Rel a
listaArel []          = M.empty
listaArel ((x,y):xys) = M.insertWith (++) x [y] (listaArel xys)

-- (composicionRel r s) es la composición de las relaciones r y s. Por
-- ejemplo,
--    λ> r = listaArel [(1,2),(5,2)]
--    λ> s = listaArel [(2,3),(2,4)]
--    λ> composicionRel r s
--    fromList [(1,[3,4]),(5,[3,4])]
composicionRel :: Ord a => Rel a -> Rel a -> Rel a
composicionRel r s =
  M.map f r 
  where f xs = concat (mapMaybe (`M.lookup` s) xs)

-- (relAlista r) es la lista de pares correspondientes a la relación
-- r. Por ejemplo,
--    λ> relAlista (M.fromList [(1,[3,4]),(5,[3,4])])
--    [(1,3),(1,4),(5,3),(5,4)]
relAlista :: Ord a => Rel a -> [(a,a)]
relAlista r =
  nub [(x,y) | (x,ys) <- M.assocs r, y <- ys] 

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

-- La propiedad es
prop_composicion :: [(Int,Int)] -> [(Int,Int)] -> Bool
prop_composicion r s =
  all (== sort (composicion1 r' s'))
      [sort (composicion2 r' s'),
       sort (composicion3 r' s')]
  where r' = nub r
        s' = nub s

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

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

-- La comparación es
--    λ> length (composicion1 [(n,n+1) | n <- [1..2000]] [(n,n+1) | n <- [1..2000]])
--    1999
--    (1.54 secs, 770,410,352 bytes)
--    λ> length (composicion2 [(n,n+1) | n <- [1..2000]] [(n,n+1) | n <- [1..2000]])
--    1999
--    (1.66 secs, 1,348,948,096 bytes)
--    λ> length (composicion3 [(n,n+1) | n <- [1..2000]] [(n,n+1) | n <- [1..2000]])
--    1999
--    (0.07 secs, 6,960,552 bytes)
--
--    λ> r100 = [(n,k) | n <- [1..100], k <- [1..n]]
--    λ> length (composicion1 r100 r100)
--    5050
--    (9.91 secs, 4,946,556,944 bytes)
--    λ> length (composicion2 r100 r100)
--    5050
--    (13.59 secs, 15,241,357,536 bytes)
--    λ> length (composicion3 r100 r100)
--    5050
--    (0.35 secs, 52,015,544 bytes)

100

101

102

103

import Data.List (nub, sort)

import Data.Maybe (mapMaybe)

import qualified Data.Set.Monad as S (Set, fromList, toList)

import qualified Data.Map as M (Map, assocs, empty, insertWith, lookup, map)

import Test.QuickCheck (quickCheck)

-- 1ª solución

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

composicion1 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion1 r s =

nub [(x,y) | (x,u) <- r, (v,y) <- s, u == v]

-- 2ª solución

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

composicion2 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion2 r s =

S.toList (composicionS (S.fromList r) (S.fromList s))

composicionS :: Ord a => S.Set (a,a) -> S.Set (a,a) -> S.Set (a,a)

composicionS r s =

[(x,y) | (x,u) <- r, (v,y) <- s, u == v]

-- 3ª solución

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

composicion3 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion3 r s =

relAlista (composicionRel (listaArel r) (listaArel s))

-- Una relación se puede representar por un diccionario donde las claves

-- son los elementos y los valores son las listas de los elementos con

-- los que se relaciona.

type Rel a = M.Map a [a]

-- (listaArel xys) es la relación correspondiente a la lista de pares

-- xys. Por ejemplo.

-- λ> listaArel [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]

-- fromList [(1,[1,2,3,6]),(2,[2,6]),(3,[3,6]),(6,[6])]

listaArel :: Ord a => [(a,a)] -> Rel a

listaArel [] = M.empty

listaArel ((x,y):xys) = M.insertWith (++) x [y] (listaArel xys)

-- (composicionRel r s) es la composición de las relaciones r y s. Por

-- ejemplo,

-- λ> r = listaArel [(1,2),(5,2)]

-- λ> s = listaArel [(2,3),(2,4)]

-- λ> composicionRel r s

-- fromList [(1,[3,4]),(5,[3,4])]

composicionRel :: Ord a => Rel a -> Rel a -> Rel a

composicionRel r s =

M.map f r

where f xs = concat (mapMaybe (`M.lookup` s) xs)

-- (relAlista r) es la lista de pares correspondientes a la relación

-- r. Por ejemplo,

-- λ> relAlista (M.fromList [(1,[3,4]),(5,[3,4])])

-- [(1,3),(1,4),(5,3),(5,4)]

relAlista :: Ord a => Rel a -> [(a,a)]

relAlista r =

nub [(x,y) | (x,ys) <- M.assocs r, y <- ys]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_composicion :: [(Int,Int)] -> [(Int,Int)] -> Bool

prop_composicion r s =

all (== sort (composicion1 r' s'))

[sort (composicion2 r' s'),

sort (composicion3 r' s')]

where r' = nub r

s' = nub s

-- La comprobación es

-- λ> quickCheck prop_composicion

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (composicion1 [(n,n+1) | n <- [1..2000]] [(n,n+1) | n <- [1..2000]])

-- 1999

-- (1.54 secs, 770,410,352 bytes)

-- λ> length (composicion2 [(n,n+1) | n <- [1..2000]] [(n,n+1) | n <- [1..2000]])

-- 1999

-- (1.66 secs, 1,348,948,096 bytes)

-- λ> length (composicion3 [(n,n+1) | n <- [1..2000]] [(n,n+1) | n <- [1..2000]])

-- 1999

-- (0.07 secs, 6,960,552 bytes)

-- λ> r100 = [(n,k) | n <- [1..100], k <- [1..n]]

-- λ> length (composicion1 r100 r100)

-- 5050

-- (9.91 secs, 4,946,556,944 bytes)

-- λ> length (composicion2 r100 r100)

-- 5050

-- (13.59 secs, 15,241,357,536 bytes)

-- λ> length (composicion3 r100 r100)

-- 5050

-- (0.35 secs, 52,015,544 bytes)

El código se encuentra en GitHub.

Número de particiones en k subconjuntos

PorJosé A. Alonso 11-julio-202229-junio-2022

Definir la función

   numeroParticiones :: Int -> Int -> Int

1	numeroParticiones :: Int -> Int -> Int

tal que (numeroParticiones n k) es el número de particiones de conjunto de n elementos en k subconjuntos disjuntos. Por ejemplo,

   numeroParticiones 3 2    ==  3
   numeroParticiones 3 3    ==  1
   numeroParticiones 4 3    ==  6
   numeroParticiones 4 1    ==  1
   numeroParticiones 4 4    ==  1
   numeroParticiones 91 89  ==  8139495

numeroParticiones 3 2 == 3

numeroParticiones 3 3 == 1

numeroParticiones 4 3 == 6

numeroParticiones 4 1 == 1

numeroParticiones 4 4 == 1

numeroParticiones 91 89 == 8139495

Soluciones

import Data.Array (Array, (!), array)
import Test.QuickCheck (Positive (Positive), quickCheckWith)

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

numeroParticiones1 :: Int -> Int -> Int
numeroParticiones1 n k =
  length (particiones1 [1..n] k)

particiones1 :: [a] -> Int -> [[[a]]]
particiones1 [] _     = []
particiones1 _  0     = []
particiones1 xs 1     = [[xs]]
particiones1 (x:xs) k = [[x]:ys | ys <- particiones1 xs (k-1)] ++ 
                        concat [inserta x ys | ys <- particiones1 xs k]

-- (inserta x yss) es la lista obtenida insertando x en cada uno de los
-- conjuntos de yss. Por ejemplo,
--    inserta 4 [[3],[2,5]]  ==  [[[4,3],[2,5]],[[3],[4,2,5]]]
inserta :: a -> [[a]] -> [[[a]]]
inserta _ []       = []
inserta x (ys:yss) = ((x:ys):yss) : [ys:zss | zss <- inserta x yss]
  
-- 2ª solución
-- ===========

numeroParticiones2 :: Int -> Int -> Int
numeroParticiones2 0 _ = 0
numeroParticiones2 _ 0 = 0
numeroParticiones2 _ 1 = 1
numeroParticiones2 n k = k * numeroParticiones2 (n-1) k +
                         numeroParticiones2 (n-1) (k-1) 

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

numeroParticiones3 :: Int -> Int -> Int
numeroParticiones3 n k = matrizNumeroParticiones n k ! (n,k)

matrizNumeroParticiones :: Int -> Int -> Array (Int,Int) Int
matrizNumeroParticiones n k = q where
  q = array ((0,0),(n,k)) [((i,j), f i j) | i <- [0..n], j <- [0..k]]
  f _ 0 = 0
  f 0 _ = 0
  f _ 1 = 1
  f i j | i == j    = 1
        | otherwise = j * f (i-1) j + f (i-1) (j-1)

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

-- La propiedad es
prop_numeroParticiones :: Positive Int -> Positive Int -> Bool
prop_numeroParticiones (Positive n) (Positive k) =
  all (== numeroParticiones1 n k)
      [numeroParticiones2 n k,
       numeroParticiones3 n k]

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

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

-- La comparación es
--    λ> numeroParticiones1 12 6
--    1323652
--    (1.22 secs, 1,152,945,608 bytes)
--    λ> numeroParticiones2 12 6
--    1323652
--    (0.00 secs, 1,283,152 bytes)
--    λ> numeroParticiones3 12 6
--    1323652
--    (0.01 secs, 1,155,864 bytes)
--
--    λ> numeroParticiones2 21 19
--    19285
--    (2.04 secs, 990,274,976 bytes)
--    λ> numeroParticiones3 21 19
--    19285
--    (0.00 secs, 940,736 bytes)

import Data.Array (Array, (!), array)

import Test.QuickCheck (Positive (Positive), quickCheckWith)

-- 1ª solución

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

numeroParticiones1 :: Int -> Int -> Int

numeroParticiones1 n k =

length (particiones1 [1..n] k)

particiones1 :: [a] -> Int -> [[[a]]]

particiones1 [] _ = []

particiones1 _ 0 = []

particiones1 xs 1 = [[xs]]

particiones1 (x:xs) k = [[x]:ys | ys <- particiones1 xs (k-1)] ++

concat [inserta x ys | ys <- particiones1 xs k]

-- (inserta x yss) es la lista obtenida insertando x en cada uno de los

-- conjuntos de yss. Por ejemplo,

-- inserta 4 [[3],[2,5]] == [[[4,3],[2,5]],[[3],[4,2,5]]]

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

inserta _ [] = []

inserta x (ys:yss) = ((x:ys):yss) : [ys:zss | zss <- inserta x yss]

-- 2ª solución

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

numeroParticiones2 :: Int -> Int -> Int

numeroParticiones2 0 _ = 0

numeroParticiones2 _ 0 = 0

numeroParticiones2 _ 1 = 1

numeroParticiones2 n k = k * numeroParticiones2 (n-1) k +

numeroParticiones2 (n-1) (k-1)

-- 3ª solución

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

numeroParticiones3 :: Int -> Int -> Int

numeroParticiones3 n k = matrizNumeroParticiones n k ! (n,k)

matrizNumeroParticiones :: Int -> Int -> Array (Int,Int) Int

matrizNumeroParticiones n k = q where

q = array ((0,0),(n,k)) [((i,j), f i j) | i <- [0..n], j <- [0..k]]

f _ 0 = 0

f 0 _ = 0

f _ 1 = 1

f i j | i == j = 1

| otherwise = j * f (i-1) j + f (i-1) (j-1)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_numeroParticiones :: Positive Int -> Positive Int -> Bool

prop_numeroParticiones (Positive n) (Positive k) =

all (== numeroParticiones1 n k)

[numeroParticiones2 n k,

numeroParticiones3 n k]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_numeroParticiones

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> numeroParticiones1 12 6

-- 1323652

-- (1.22 secs, 1,152,945,608 bytes)

-- λ> numeroParticiones2 12 6

-- 1323652

-- (0.00 secs, 1,283,152 bytes)

-- λ> numeroParticiones3 12 6

-- 1323652

-- (0.01 secs, 1,155,864 bytes)

-- λ> numeroParticiones2 21 19

-- 19285

-- (2.04 secs, 990,274,976 bytes)

-- λ> numeroParticiones3 21 19

-- 19285

-- (0.00 secs, 940,736 bytes)

El código se encuentra en GitHub.

Particiones en k subconjuntos

PorJosé A. Alonso 8-julio-20226-julio-2022

Definir la función

   particiones :: [a] -> Int -> [[[a]]]

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

tal que (particiones xs k) es la lista de las particiones de xs en k subconjuntos disjuntos. Por ejemplo,

   λ> particiones [2,3,6] 2
   [[[2],[3,6]],[[2,3],[6]],[[3],[2,6]]]
   λ> particiones [2,3,6] 3
   [[[2],[3],[6]]]
   λ> particiones [4,2,3,6] 3
   [[[4],[2],[3,6]],[[4],[2,3],[6]],[[4],[3],[2,6]],
    [[4,2],[3],[6]],[[2],[4,3],[6]],[[2],[3],[4,6]]]
   λ> particiones [4,2,3,6] 1
   [[[4,2,3,6]]]
   λ> particiones [4,2,3,6] 4
   [[[4],[2],[3],[6]]]

λ> particiones [2,3,6] 2

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

λ> particiones [2,3,6] 3

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

λ> particiones [4,2,3,6] 3

[[[4],[2],[3,6]],[[4],[2,3],[6]],[[4],[3],[2,6]],

[[4,2],[3],[6]],[[2],[4,3],[6]],[[2],[3],[4,6]]]

λ> particiones [4,2,3,6] 1

[[[4,2,3,6]]]

λ> particiones [4,2,3,6] 4

[[[4],[2],[3],[6]]]

Soluciones

import Data.List (nub, sort)
import Data.Array (Array, (!), array, listArray)
import Test.QuickCheck (Positive (Positive), quickCheckWith)

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

particiones1 :: [a] -> Int -> [[[a]]]
particiones1 [] _     = []
particiones1 _  0     = []
particiones1 xs 1     = [[xs]]
particiones1 (x:xs) k = [[x]:ys | ys <- particiones1 xs (k-1)] ++ 
                        concat [inserta x ys | ys <- particiones1 xs k]

-- (inserta x yss) es la lista obtenida insertando x en cada uno de los
-- conjuntos de yss. Por ejemplo,
--    inserta 4 [[3],[2,5]]  ==  [[[4,3],[2,5]],[[3],[4,2,5]]]
inserta :: a -> [[a]] -> [[[a]]]
inserta _ []       = []
inserta x (ys:yss) = ((x:ys):yss) : [ys:zss | zss <- inserta x yss]

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

particiones2 :: [a] -> Int -> [[[a]]]
particiones2 [] _     = []
particiones2 _  0     = []
particiones2 xs 1     = [[xs]]
particiones2 (x:xs) k = map ([x]:) (particiones2 xs (k-1)) ++ 
                        concatMap (inserta x) (particiones2 xs k)

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

particiones3 :: [a] -> Int -> [[[a]]]
particiones3 xs k = matrizParticiones xs k ! (length xs, k)

matrizParticiones :: [a] -> Int -> Array (Int,Int) [[[a]]]
matrizParticiones xs k = q where
  q = array ((0,0),(n,k)) [((i,j), f i j) | i <- [0..n], j <- [0..k]]
  n = length xs
  v = listArray (1,n) xs
  f _ 0 = []
  f 0 _ = []
  f m 1 = [[take m xs]]
  f i j | i == j = [[[x] | x <- take i xs]]
        | otherwise = map ([v!i] :) (q!(i-1,j-1)) ++
                      concatMap (inserta (v!i)) (q!(i-1,j))

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

-- La propiedad es
prop_particiones :: [Int] -> Positive Int -> Bool
prop_particiones xs (Positive k) =
  all (iguales (particiones1 xs' k))
      [particiones2 xs' k,
       particiones3 xs' k]
  where
    xs' = nub xs
    iguales xss yss = sort (map sort [map sort x | x <- xss]) ==
                      sort (map sort [map sort y | y <- yss])

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

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

-- La comparación es
--    λ> length (particiones1 [1..12] 6)
--    1323652
--    (1.33 secs, 1,152,945,584 bytes)
--    λ> length (particiones2 [1..12] 6)
--    1323652
--    (1.07 secs, 1,104,960,360 bytes)
--    λ> length (particiones3 [1..12] 6)
--    1323652
--    (1.68 secs, 1,047,004,368 bytes)

import Data.List (nub, sort)

import Data.Array (Array, (!), array, listArray)

import Test.QuickCheck (Positive (Positive), quickCheckWith)

-- 1ª solución

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

particiones1 :: [a] -> Int -> [[[a]]]

particiones1 [] _ = []

particiones1 _ 0 = []

particiones1 xs 1 = [[xs]]

particiones1 (x:xs) k = [[x]:ys | ys <- particiones1 xs (k-1)] ++

concat [inserta x ys | ys <- particiones1 xs k]

-- (inserta x yss) es la lista obtenida insertando x en cada uno de los

-- conjuntos de yss. Por ejemplo,

-- inserta 4 [[3],[2,5]] == [[[4,3],[2,5]],[[3],[4,2,5]]]

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

inserta _ [] = []

inserta x (ys:yss) = ((x:ys):yss) : [ys:zss | zss <- inserta x yss]

-- 2ª solución

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

particiones2 :: [a] -> Int -> [[[a]]]

particiones2 [] _ = []

particiones2 _ 0 = []

particiones2 xs 1 = [[xs]]

particiones2 (x:xs) k = map ([x]:) (particiones2 xs (k-1)) ++

concatMap (inserta x) (particiones2 xs k)

-- 3ª solución

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

particiones3 :: [a] -> Int -> [[[a]]]

particiones3 xs k = matrizParticiones xs k ! (length xs, k)

matrizParticiones :: [a] -> Int -> Array (Int,Int) [[[a]]]

matrizParticiones xs k = q where

q = array ((0,0),(n,k)) [((i,j), f i j) | i <- [0..n], j <- [0..k]]

n = length xs

v = listArray (1,n) xs

f _ 0 = []

f 0 _ = []

f m 1 = [[take m xs]]

f i j | i == j = [[[x] | x <- take i xs]]

| otherwise = map ([v!i] :) (q!(i-1,j-1)) ++

concatMap (inserta (v!i)) (q!(i-1,j))

-- Comprobación de equivalencia

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

-- La propiedad es

prop_particiones :: [Int] -> Positive Int -> Bool

prop_particiones xs (Positive k) =

all (iguales (particiones1 xs' k))

[particiones2 xs' k,

particiones3 xs' k]

where

xs' = nub xs

iguales xss yss = sort (map sort [map sort x | x <- xss]) ==

sort (map sort [map sort y | y <- yss])

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_particiones

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (particiones1 [1..12] 6)

-- 1323652

-- (1.33 secs, 1,152,945,584 bytes)

-- λ> length (particiones2 [1..12] 6)

-- 1323652

-- (1.07 secs, 1,104,960,360 bytes)

-- λ> length (particiones3 [1..12] 6)

-- 1323652

-- (1.68 secs, 1,047,004,368 bytes)

El código se encuentra en GitHub.

Mayor semiprimo menor que n

PorJosé A. Alonso 7-julio-20226-julio-2022

Un número semiprimo es un número natural es producto de dos números primos no necesariamente distintos. Por ejemplo, 26 es semiprimo (porque 26 = 2·13) y 49 también lo es (porque 49 = 7·7).

Definir la función

   mayorSemiprimoMenor :: Integer -> Integer

1	mayorSemiprimoMenor :: Integer -> Integer

tal que (mayorSemiprimoMenor n) es el mayor semiprimo menor que n (suponiendo que n > 4). Por ejemplo,

   mayorSemiprimoMenor 27      ==  26
   mayorSemiprimoMenor 50      ==  49
   mayorSemiprimoMenor 49      ==  46
   mayorSemiprimoMenor (10^15) == 999999999999998

mayorSemiprimoMenor 27 == 26

mayorSemiprimoMenor 50 == 49

mayorSemiprimoMenor 49 == 46

mayorSemiprimoMenor (10^15) == 999999999999998

Soluciones

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

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

mayorSemiprimoMenor1 :: Integer -> Integer
mayorSemiprimoMenor1 n =
  head [x | x <- [n-1,n-2..2], semiPrimo x]

semiPrimo :: Integer -> Bool
semiPrimo n =
  not (null [x | x <- [n,n-1..2], 
                 primo x,
                 n `mod` x == 0,
                 primo (n `div` x)])

primo :: Integer -> Bool
primo n = [x | x <- [1..n], n `mod` x == 0] == [1,n] 

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

mayorSemiprimoMenor2 :: Integer -> Integer
mayorSemiprimoMenor2 n =
  head [x | x <- [n-1,n-2..2], semiPrimo2 x]

semiPrimo2 :: Integer -> Bool
semiPrimo2 n =
  not (null [x | x <- [n-1,n-2..2], 
                 isPrime x,
                 n `mod` x == 0,
                 isPrime (n `div` x)])

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

mayorSemiprimoMenor3 :: Integer -> Integer
mayorSemiprimoMenor3 n =
  head [x | x <- [n-1,n-2..2], semiPrimo3 x]

semiPrimo3 :: Integer -> Bool
semiPrimo3 n =
  not (null [x | x <- reverse (takeWhile (<n) primes),
                 n `mod` x == 0,
                 isPrime (n `div` x)])

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

mayorSemiprimoMenor4 :: Integer -> Integer
mayorSemiprimoMenor4 n =
  head [ p | p <- [n-1,n-2..2]
           , (length . primeFactors) p == 2]

-- 5ª solución
-- ===========

mayorSemiprimoMenor5 :: Integer -> Integer
mayorSemiprimoMenor5 n
  | semiPrimo5 (n-1) = n-1
  | otherwise        = mayorSemiprimoMenor5 (n-1)

semiPrimo5 :: Integer -> Bool
semiPrimo5 x = length (primeFactors x) == 2

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

-- La propiedad es
prop_mayorSemiprimoMenor :: Integer -> Property
prop_mayorSemiprimoMenor n =
  n > 4 ==>
  all (== mayorSemiprimoMenor1 n)
      [mayorSemiprimoMenor2 n,
       mayorSemiprimoMenor3 n,
       mayorSemiprimoMenor4 n,
       mayorSemiprimoMenor5 n]

-- La comprobación es
--    λ> quickCheck prop_mayorSemiprimoMenor
--    +++ OK, passed 100 tests; 353 discarded.

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

-- La comparación es
--    λ> mayorSemiprimoMenor1 5000
--    4997
--    (1.92 secs, 945,507,880 bytes)
--    λ> mayorSemiprimoMenor2 5000
--    4997
--    (0.05 secs, 123,031,264 bytes)
--    λ> mayorSemiprimoMenor3 5000
--    4997
--    (0.01 secs, 5,865,120 bytes)
--    λ> mayorSemiprimoMenor4 5000
--    4997
--    (0.00 secs, 593,528 bytes)
--    λ> mayorSemiprimoMenor5 5000
--    4997
--    (0.00 secs, 593,200 bytes)
--
--    λ> mayorSemiprimoMenor3 (3*10^6)
--    2999995
--    (2.34 secs, 6,713,620,000 bytes)
--    λ> mayorSemiprimoMenor4 (2*10^6)
--    1999997
--    (0.01 secs, 728,936 bytes)
--    λ> mayorSemiprimoMenor5 (2*10^6)
--    1999997
--    (0.01 secs, 728,608 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

import Data.Numbers.Primes (primeFactors, isPrime, primes)

import Test.QuickCheck

-- 1ª solución

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

mayorSemiprimoMenor1 :: Integer -> Integer

mayorSemiprimoMenor1 n =

head [x | x <- [n-1,n-2..2], semiPrimo x]

semiPrimo :: Integer -> Bool

semiPrimo n =

not (null [x | x <- [n,n-1..2],

primo x,

n `mod` x == 0,

primo (n `div` x)])

primo :: Integer -> Bool

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

-- 2ª solución

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

mayorSemiprimoMenor2 :: Integer -> Integer

mayorSemiprimoMenor2 n =

head [x | x <- [n-1,n-2..2], semiPrimo2 x]

semiPrimo2 :: Integer -> Bool

semiPrimo2 n =

not (null [x | x <- [n-1,n-2..2],

isPrime x,

n `mod` x == 0,

isPrime (n `div` x)])

-- 3ª solución

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

mayorSemiprimoMenor3 :: Integer -> Integer

mayorSemiprimoMenor3 n =

head [x | x <- [n-1,n-2..2], semiPrimo3 x]

semiPrimo3 :: Integer -> Bool

semiPrimo3 n =

not (null [x | x <- reverse (takeWhile (<n) primes),

n `mod` x == 0,

isPrime (n `div` x)])

-- 4ª solución

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

mayorSemiprimoMenor4 :: Integer -> Integer

mayorSemiprimoMenor4 n =

head [ p | p <- [n-1,n-2..2]

, (length . primeFactors) p == 2]

-- 5ª solución

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

mayorSemiprimoMenor5 :: Integer -> Integer

mayorSemiprimoMenor5 n

| semiPrimo5 (n-1) = n-1

| otherwise = mayorSemiprimoMenor5 (n-1)

semiPrimo5 :: Integer -> Bool

semiPrimo5 x = length (primeFactors x) == 2

-- Comprobación de equivalencia

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

-- La propiedad es

prop_mayorSemiprimoMenor :: Integer -> Property

prop_mayorSemiprimoMenor n =

n > 4 ==>

all (== mayorSemiprimoMenor1 n)

[mayorSemiprimoMenor2 n,

mayorSemiprimoMenor3 n,

mayorSemiprimoMenor4 n,

mayorSemiprimoMenor5 n]

-- La comprobación es

-- λ> quickCheck prop_mayorSemiprimoMenor

-- +++ OK, passed 100 tests; 353 discarded.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> mayorSemiprimoMenor1 5000

-- 4997

-- (1.92 secs, 945,507,880 bytes)

-- λ> mayorSemiprimoMenor2 5000

-- 4997

-- (0.05 secs, 123,031,264 bytes)

-- λ> mayorSemiprimoMenor3 5000

-- 4997

-- (0.01 secs, 5,865,120 bytes)

-- λ> mayorSemiprimoMenor4 5000

-- 4997

-- (0.00 secs, 593,528 bytes)

-- λ> mayorSemiprimoMenor5 5000

-- 4997

-- (0.00 secs, 593,200 bytes)

-- λ> mayorSemiprimoMenor3 (3*10^6)

-- 2999995

-- (2.34 secs, 6,713,620,000 bytes)

-- λ> mayorSemiprimoMenor4 (2*10^6)

-- 1999997

-- (0.01 secs, 728,936 bytes)

-- λ> mayorSemiprimoMenor5 (2*10^6)

-- 1999997

-- (0.01 secs, 728,608 bytes)

El código se encuentra en GitHub.

Intersecciones parciales

PorJosé A. Alonso 6-julio-20226-julio-2022

Definir la función

   interseccionParcial :: Ord a => Int -> [[a]] -> [a]

1	interseccionParcial :: Ord a => Int -> [[a]] -> [a]

tal que (interseccionParcial n xss) es la lista de los elementos que pertenecen al menos a n conjuntos de xss. Por ejemplo,

   interseccionParcial 1 [[3,4],[4,5,9],[5,4,7]]  == [3,4,5,9,7]
   interseccionParcial 2 [[3,4],[4,5,9],[5,4,7]]  == [4,5]
   interseccionParcial 3 [[3,4],[4,5,9],[5,4,7]]  == [4]
   interseccionParcial 4 [[3,4],[4,5,9],[5,4,7]]  == []

interseccionParcial 1 [[3,4],[4,5,9],[5,4,7]] == [3,4,5,9,7]

interseccionParcial 2 [[3,4],[4,5,9],[5,4,7]] == [4,5]

interseccionParcial 3 [[3,4],[4,5,9],[5,4,7]] == [4]

interseccionParcial 4 [[3,4],[4,5,9],[5,4,7]] == []

Soluciones

import Data.List (foldl', nub, union, sort)
import Test.QuickCheck

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

interseccionParcial1 :: Ord a => Int -> [[a]] -> [a]
interseccionParcial1 n xss = 
  [x | x <- sort (elementos xss)
     , pertenecenAlMenos n xss x]

elementos :: Ord a => [[a]] -> [a]
elementos []       = []
elementos (xs:xss) = xs `union` elementos xss

pertenecenAlMenos :: Ord a => Int -> [[a]] -> a -> Bool
pertenecenAlMenos n xss x =
  length [xs | xs <- xss, x `elem` xs] >= n
  
-- 2ª solución
-- ===========

interseccionParcial2 :: Ord a => Int -> [[a]] -> [a]
interseccionParcial2 n xss = 
  [x | x <- sort (elementos2 xss)
     , pertenecenAlMenos2 n xss x]

elementos2 :: Ord a => [[a]] -> [a]
elementos2 = foldl' union []
  
pertenecenAlMenos2 :: Ord a => Int -> [[a]] -> a -> Bool
pertenecenAlMenos2 n xss x =
  length (filter (x `elem`) xss) >= n
  
-- 3ª solución
-- ===========

interseccionParcial3 :: Ord a => Int -> [[a]] -> [a]
interseccionParcial3 n xss = 
  [x | x <- sort (nub (concat xss))
     , length [xs | xs <- xss, x `elem` xs] >= n]

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

-- La propiedad es
prop_interseccionParcial :: Positive Int -> [[Int]] -> Bool
prop_interseccionParcial (Positive n) xss =
  all (== interseccionParcial1 n yss)
      [interseccionParcial2 n yss,
       interseccionParcial3 n yss]
  where yss = map nub xss

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

import Data.List (foldl', nub, union, sort)

import Test.QuickCheck

-- 1ª solución

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

interseccionParcial1 :: Ord a => Int -> [[a]] -> [a]

interseccionParcial1 n xss =

[x | x <- sort (elementos xss)

, pertenecenAlMenos n xss x]

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

elementos [] = []

elementos (xs:xss) = xs `union` elementos xss

pertenecenAlMenos :: Ord a => Int -> [[a]] -> a -> Bool

pertenecenAlMenos n xss x =

length [xs | xs <- xss, x `elem` xs] >= n

-- 2ª solución

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

interseccionParcial2 :: Ord a => Int -> [[a]] -> [a]

interseccionParcial2 n xss =

[x | x <- sort (elementos2 xss)

, pertenecenAlMenos2 n xss x]

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

elementos2 = foldl' union []

pertenecenAlMenos2 :: Ord a => Int -> [[a]] -> a -> Bool

pertenecenAlMenos2 n xss x =

length (filter (x `elem`) xss) >= n

-- 3ª solución

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

interseccionParcial3 :: Ord a => Int -> [[a]] -> [a]

interseccionParcial3 n xss =

[x | x <- sort (nub (concat xss))

, length [xs | xs <- xss, x `elem` xs] >= n]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_interseccionParcial :: Positive Int -> [[Int]] -> Bool

prop_interseccionParcial (Positive n) xss =

all (== interseccionParcial1 n yss)

[interseccionParcial2 n yss,

interseccionParcial3 n yss]

where yss = map nub xss

-- La comprobación es

-- λ> quickCheck prop_interseccionParcial

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Unión e intersección general de conjuntos

PorJosé A. Alonso 5-julio-20226-julio-2022

Definir las funciones

   unionGeneral        :: Eq a => [[a]] -> [a]
   interseccionGeneral :: Eq a => [[a]] -> [a]

1 2	unionGeneral :: Eq a => [[a]] -> [a] interseccionGeneral :: Eq a => [[a]] -> [a]

tales que

(unionGeneral xs) es la unión de los conjuntos de la lista de conjuntos xs (es decir, el conjunto de los elementos que pertenecen a alguno de los elementos de xs). Por ejemplo,

     unionGeneral []                    ==  []
     unionGeneral [[1]]                 ==  [1]
     unionGeneral [[1],[1,2],[2,3]]     ==  [1,2,3]
     unionGeneral ([[x] | x <- [1..9]]) ==  [1,2,3,4,5,6,7,8,9]

unionGeneral [] == []

unionGeneral [[1]] == [1]

unionGeneral [[1],[1,2],[2,3]] == [1,2,3]

unionGeneral ([[x] | x <- [1..9]]) == [1,2,3,4,5,6,7,8,9]

(interseccionGeneral xs) es la intersección de los conjuntos de la lista de conjuntos xs (es decir, el conjunto de los elementos que pertenecen a todos los elementos de xs). Por ejemplo,

     interseccionGeneral [[1]]                      ==  [1]
     interseccionGeneral [[2],[1,2],[2,3]]          ==  [2]
     interseccionGeneral [[2,7,5],[1,5,2],[5,2,3]]  ==  [2,5]
     interseccionGeneral ([[x] | x <- [1..9]])      ==  []
     interseccionGeneral (replicate (10^6) [1..5])  ==  [1,2,3,4,5]

interseccionGeneral [[1]] == [1]

interseccionGeneral [[2],[1,2],[2,3]] == [2]

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

interseccionGeneral ([[x] | x <- [1..9]]) == []

interseccionGeneral (replicate (10^6) [1..5]) == [1,2,3,4,5]

Soluciones

import Data.List (foldl', foldl1', intersect, nub, union)
import Test.QuickCheck (NonEmptyList (NonEmpty), quickCheck)

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

unionGeneral1 :: Eq a => [[a]] -> [a]
unionGeneral1 []     = []
unionGeneral1 (x:xs) = x `union` unionGeneral1 xs 

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

unionGeneral2 :: Eq a => [[a]] -> [a]
unionGeneral2 = foldr union []

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

unionGeneral3 :: Eq a => [[a]] -> [a]
unionGeneral3 = foldl' union []

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

-- La propiedad es
prop_unionGeneral :: [[Int]] -> Bool
prop_unionGeneral xss =
  all (== unionGeneral1 xss')
      [unionGeneral2 xss',
       unionGeneral3 xss']
  where xss' = nub (map nub xss)
  
-- La comprobación es
--    λ> quickCheck prop_unionGeneral
--    +++ OK, passed 100 tests.
--    (0.85 secs, 1,017,807,600 bytes)

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

-- La comparación es
--    λ> length (unionGeneral1 ([[x] | x <- [1..10^3]]))
--    1000
--    (1.56 secs, 107,478,456 bytes)
--    λ> length (unionGeneral2 ([[x] | x <- [1..10^3]]))
--    1000
--    (1.50 secs, 107,406,560 bytes)
--    λ> length (unionGeneral3 ([[x] | x <- [1..10^3]]))
--    1000
--    (0.07 secs, 92,874,024 bytes)

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

interseccionGeneral1 :: Eq a => [[a]] -> [a]
interseccionGeneral1 [x]    = x
interseccionGeneral1 (x:xs) = x `intersect` interseccionGeneral1 xs 

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

interseccionGeneral2 :: Eq a => [[a]] -> [a]
interseccionGeneral2 = foldr1 intersect

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

interseccionGeneral3 :: Eq a => [[a]] -> [a]
interseccionGeneral3 = foldl1' intersect

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

-- La propiedad es
prop_interseccionGeneral :: NonEmptyList [Int] -> Bool
prop_interseccionGeneral (NonEmpty xss) =
  all (== interseccionGeneral1 xss')
      [interseccionGeneral2 xss',
       interseccionGeneral3 xss']
  where xss' = nub (map nub xss)

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

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

-- La comparación es
--    λ> interseccionGeneral1 (replicate (10^6) [1..5])
--    [1,2,3,4,5]
--    (2.02 secs, 1,173,618,400 bytes)
--    λ> interseccionGeneral2 (replicate (10^6) [1..5])
--    [1,2,3,4,5]
--    (1.83 secs, 1,092,120,224 bytes)
--    λ> interseccionGeneral3 (replicate (10^6) [1..5])
--    [1,2,3,4,5]
--    (1.33 secs, 985,896,136 bytes)

import Data.List (foldl', foldl1', intersect, nub, union)

import Test.QuickCheck (NonEmptyList (NonEmpty), quickCheck)

-- 1ª definición de unionGeneral

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

unionGeneral1 :: Eq a => [[a]] -> [a]

unionGeneral1 [] = []

unionGeneral1 (x:xs) = x `union` unionGeneral1 xs

-- 2ª definición de unionGeneral

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

unionGeneral2 :: Eq a => [[a]] -> [a]

unionGeneral2 = foldr union []

-- 3ª definición de unionGeneral

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

unionGeneral3 :: Eq a => [[a]] -> [a]

unionGeneral3 = foldl' union []

-- Comprobación de equivalencia de unionGeneral

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

-- La propiedad es

prop_unionGeneral :: [[Int]] -> Bool

prop_unionGeneral xss =

all (== unionGeneral1 xss')

[unionGeneral2 xss',

unionGeneral3 xss']

where xss' = nub (map nub xss)

-- La comprobación es

-- λ> quickCheck prop_unionGeneral

-- +++ OK, passed 100 tests.

-- (0.85 secs, 1,017,807,600 bytes)

-- Comparación de eficiencia de unionGeneral

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

-- La comparación es

-- λ> length (unionGeneral1 ([[x] | x <- [1..10^3]]))

-- 1000

-- (1.56 secs, 107,478,456 bytes)

-- λ> length (unionGeneral2 ([[x] | x <- [1..10^3]]))

-- 1000

-- (1.50 secs, 107,406,560 bytes)

-- λ> length (unionGeneral3 ([[x] | x <- [1..10^3]]))

-- 1000

-- (0.07 secs, 92,874,024 bytes)

-- 1ª definición de interseccionGeneral

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

interseccionGeneral1 :: Eq a => [[a]] -> [a]

interseccionGeneral1 [x] = x

interseccionGeneral1 (x:xs) = x `intersect` interseccionGeneral1 xs

-- 2ª definición de interseccionGeneral

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

interseccionGeneral2 :: Eq a => [[a]] -> [a]

interseccionGeneral2 = foldr1 intersect

-- 3ª definición de interseccionGeneral

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

interseccionGeneral3 :: Eq a => [[a]] -> [a]

interseccionGeneral3 = foldl1' intersect

-- Comprobación de equivalencia

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

-- La propiedad es

prop_interseccionGeneral :: NonEmptyList [Int] -> Bool

prop_interseccionGeneral (NonEmpty xss) =

all (== interseccionGeneral1 xss')

[interseccionGeneral2 xss',

interseccionGeneral3 xss']

where xss' = nub (map nub xss)

-- La comprobación es

-- λ> quickCheck prop_interseccionGeneral

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> interseccionGeneral1 (replicate (10^6) [1..5])

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

-- (2.02 secs, 1,173,618,400 bytes)

-- λ> interseccionGeneral2 (replicate (10^6) [1..5])

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

-- (1.83 secs, 1,092,120,224 bytes)

-- λ> interseccionGeneral3 (replicate (10^6) [1..5])

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

-- (1.33 secs, 985,896,136 bytes)

El código se encuentra en GitHub.

Ceros finales del factorial

PorJosé A. Alonso 4-julio-20226-julio-2022

Definir la función

   cerosDelFactorial :: Integer -> Integer

1	cerosDelFactorial :: Integer -> Integer

tal que (cerosDelFactorial n) es el número de ceros en que termina el factorial de n. Por ejemplo,

   cerosDelFactorial 24                         == 4
   cerosDelFactorial 25                         == 6
   length (show (cerosDelFactorial (10^70000))) == 70000

cerosDelFactorial 24 == 4

cerosDelFactorial 25 == 6

length (show (cerosDelFactorial (10^70000))) == 70000

Soluciones

import Data.List (genericLength)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

cerosDelFactorial1 :: Integer -> Integer
cerosDelFactorial1 n = ceros (factorial n)

-- (factorial n) es el factorial n. Por ejemplo,
--    factorial 3  ==  6
factorial :: Integer -> Integer
factorial n = product [1..n]

-- (ceros n) es el número de ceros en los que termina el número n. Por
-- ejemplo, 
--    ceros 320000  ==  4
ceros :: Integer -> Integer
ceros n | rem n 10 /= 0 = 0
        | otherwise     = 1 + ceros (div n 10)

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

cerosDelFactorial2 :: Integer -> Integer
cerosDelFactorial2 = ceros2 . factorial 

ceros2 :: Integer -> Integer
ceros2 n = genericLength (takeWhile (=='0') (reverse (show n)))

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

cerosDelFactorial3 :: Integer -> Integer
cerosDelFactorial3 n
  | n < 5     = 0
  | otherwise = m + cerosDelFactorial3 m
  where m = n `div` 5

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

-- La propiedad es
prop_cerosDelFactorial :: Positive Integer -> Bool
prop_cerosDelFactorial (Positive n) =
  all (== cerosDelFactorial1 n)
      [cerosDelFactorial2 n,
       cerosDelFactorial3 n]

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

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

-- La comparación es
--    λ> cerosDelFactorial1 (4*10^4)
--    9998
--    (1.93 secs, 2,296,317,904 bytes)
--    λ> cerosDelFactorial2 (4*10^4)
--    9998
--    (1.57 secs, 1,636,242,040 bytes)
--    λ> cerosDelFactorial3 (4*10^4)
--    9998
--    (0.02 secs, 527,584 bytes)

import Data.List (genericLength)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

cerosDelFactorial1 :: Integer -> Integer

cerosDelFactorial1 n = ceros (factorial n)

-- (factorial n) es el factorial n. Por ejemplo,

-- factorial 3 == 6

factorial :: Integer -> Integer

factorial n = product [1..n]

-- (ceros n) es el número de ceros en los que termina el número n. Por

-- ejemplo,

-- ceros 320000 == 4

ceros :: Integer -> Integer

ceros n | rem n 10 /= 0 = 0

| otherwise = 1 + ceros (div n 10)

-- 2ª solución

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

cerosDelFactorial2 :: Integer -> Integer

cerosDelFactorial2 = ceros2 . factorial

ceros2 :: Integer -> Integer

ceros2 n = genericLength (takeWhile (=='0') (reverse (show n)))

-- 3ª solución

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

cerosDelFactorial3 :: Integer -> Integer

cerosDelFactorial3 n

| n < 5 = 0

| otherwise = m + cerosDelFactorial3 m

where m = n `div` 5

-- Comprobación de equivalencia

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

-- La propiedad es

prop_cerosDelFactorial :: Positive Integer -> Bool

prop_cerosDelFactorial (Positive n) =

all (== cerosDelFactorial1 n)

[cerosDelFactorial2 n,

cerosDelFactorial3 n]

-- La comprobación es

-- λ> quickCheck prop_cerosDelFactorial

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> cerosDelFactorial1 (4*10^4)

-- 9998

-- (1.93 secs, 2,296,317,904 bytes)

-- λ> cerosDelFactorial2 (4*10^4)

-- 9998

-- (1.57 secs, 1,636,242,040 bytes)

-- λ> cerosDelFactorial3 (4*10^4)

-- 9998

-- (0.02 secs, 527,584 bytes)

El código se encuentra en GitHub.

Números autodescriptivos

PorJosé A. Alonso 1-julio-202221-junio-2022

Un número n es autodescriptivo cuando para cada posición k de n (empezando a contar las posiciones a partir de 0), el dígito en la posición k es igual al número de veces que ocurre k en n. Por ejemplo, 1210 es autodescriptivo porque tiene 1 dígito igual a «0», 2 dígitos iguales a «1», 1 dígito igual a «2» y ningún dígito igual a «3».

Definir la función

   autodescriptivo :: Integer -> Bool

1	autodescriptivo :: Integer -> Bool

tal que (autodescriptivo n) se verifica si n es autodescriptivo. Por ejemplo,

   λ> autodescriptivo 1210
   True
   λ> [x | x <- [1..100000], autodescriptivo x]
   [1210,2020,21200]
   λ> autodescriptivo 9210000001000
   True

λ> autodescriptivo 1210

True

λ> [x | x <- [1..100000], autodescriptivo x]

[1210,2020,21200]

λ> autodescriptivo 9210000001000

True

Soluciones

import Data.Char (digitToInt)
import Data.List (genericLength)
import Test.QuickCheck

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

autodescriptivo1 :: Integer -> Bool
autodescriptivo1 n = autodescriptiva (digitos n)

-- (digitos n) es la lista de los dígitos de n. Por ejemplo.
--    digitos 325 == [3,2,5]
digitos :: Integer -> [Integer]
digitos n = [read [d] | d <- show n]

-- (autodescriptiva ns) se verifica si la lista de dígitos ns es
-- autodescriptiva; es decir, si para cada posición k de ns
-- (empezando a contar las posiciones a partir de 0), el dígito en la
-- posición k es igual al número de veces que ocurre k en ns. Por
-- ejemplo, 
--    autodescriptiva [1,2,1,0] == True
--    autodescriptiva [1,2,1,1] == False
autodescriptiva :: [Integer] -> Bool
autodescriptiva ns = 
  and [x == ocurrencias k ns | (k,x) <- zip [0..] ns]

-- (ocurrencias x ys) es el número de veces que ocurre x en ys. Por
-- ejemplo, 
--    ocurrencias 1 [1,2,1,0,1] == 3
ocurrencias :: Integer -> [Integer] -> Integer
ocurrencias x ys = genericLength (filter (==x) ys)

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

autodescriptivo2 :: Integer -> Bool
autodescriptivo2 n =
  and (zipWith (==) (map digitToInt xs)
                    [length (filter (==c) xs) |  c <- ['0'..'9']])
  where xs = show n

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

-- La propiedad es
prop_autodescriptivo :: Positive Integer -> Bool
prop_autodescriptivo (Positive n) =
  autodescriptivo1 n == autodescriptivo2 n

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

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

-- La comparación es
--    λ> [x | x <- [1..3*10^5], autodescriptivo1 x]
--    [1210,2020,21200]
--    (2.59 secs, 6,560,244,696 bytes)
--    λ> [x | x <- [1..3*10^5], autodescriptivo2 x]
--    [1210,2020,21200]
--    (0.67 secs, 425,262,848 bytes)

import Data.Char (digitToInt)

import Data.List (genericLength)

import Test.QuickCheck

-- 1ª solución

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

autodescriptivo1 :: Integer -> Bool

autodescriptivo1 n = autodescriptiva (digitos n)

-- (digitos n) es la lista de los dígitos de n. Por ejemplo.

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

digitos :: Integer -> [Integer]

digitos n = [read [d] | d <- show n]

-- (autodescriptiva ns) se verifica si la lista de dígitos ns es

-- autodescriptiva; es decir, si para cada posición k de ns

-- (empezando a contar las posiciones a partir de 0), el dígito en la

-- posición k es igual al número de veces que ocurre k en ns. Por

-- ejemplo,

-- autodescriptiva [1,2,1,0] == True

-- autodescriptiva [1,2,1,1] == False

autodescriptiva :: [Integer] -> Bool

autodescriptiva ns =

and [x == ocurrencias k ns | (k,x) <- zip [0..] ns]

-- (ocurrencias x ys) es el número de veces que ocurre x en ys. Por

-- ejemplo,

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

ocurrencias :: Integer -> [Integer] -> Integer

ocurrencias x ys = genericLength (filter (==x) ys)

-- 2ª solución

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

autodescriptivo2 :: Integer -> Bool

autodescriptivo2 n =

and (zipWith (==) (map digitToInt xs)

[length (filter (==c) xs) | c <- ['0'..'9']])

where xs = show n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_autodescriptivo :: Positive Integer -> Bool

prop_autodescriptivo (Positive n) =

autodescriptivo1 n == autodescriptivo2 n

-- La comprobación es

-- λ> quickCheck prop_autodescriptivo

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> [x | x <- [1..3*10^5], autodescriptivo1 x]

-- [1210,2020,21200]

-- (2.59 secs, 6,560,244,696 bytes)

-- λ> [x | x <- [1..3*10^5], autodescriptivo2 x]

-- [1210,2020,21200]

-- (0.67 secs, 425,262,848 bytes)

El código se encuentra en GitHub.

Aproximación del número pi

PorJosé A. Alonso 30-junio-20222-julio-2022

Una forma de aproximar el número π es usando la siguiente igualdad:

    π         1     1·2     1·2·3     1·2·3·4     
   --- = 1 + --- + ----- + ------- + --------- + ....
    2         3     3·5     3·5·7     3·5·7·9

π 1 1·2 1·2·3 1·2·3·4

--- = 1 + --- + ----- + ------- + --------- + ....

2 3 3·5 3·5·7 3·5·7·9

Es decir, la serie cuyo término general n-ésimo es el cociente entre el producto de los primeros n números y los primeros n números impares:

               Π i   
   s(n) =  -----------
            Π (2*i+1)

Π i

s(n) = -----------

Π (2*i+1)

Definir la función

   aproximaPi :: Integer -> Double

1	aproximaPi :: Integer -> Double

tal que (aproximaPi n) es la aproximación del número π calculada con la serie anterior hasta el término n-ésimo. Por ejemplo,

   aproximaPi 10     == 3.1411060206
   aproximaPi 20     == 3.1415922987403397
   aproximaPi 30     == 3.1415926533011596
   aproximaPi 40     == 3.1415926535895466
   aproximaPi 50     == 3.141592653589793
   aproximaPi (10^4) == 3.141592653589793
   pi                == 3.141592653589793

aproximaPi 10 == 3.1411060206

aproximaPi 20 == 3.1415922987403397

aproximaPi 30 == 3.1415926533011596

aproximaPi 40 == 3.1415926535895466

aproximaPi 50 == 3.141592653589793

aproximaPi (10^4) == 3.141592653589793

pi == 3.141592653589793

Soluciones

import Data.Ratio ((%))
import Data.List (genericTake)
import Test.QuickCheck (Property, arbitrary, forAll, suchThat, quickCheck)

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

aproximaPi1 :: Integer -> Double
aproximaPi1 n = 
  fromRational (2 * sum [product [1..i] % product [1,3..2*i+1] | i <- [0..n]])

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

aproximaPi2 :: Integer -> Double
aproximaPi2 0 = 2
aproximaPi2 n = 
  aproximaPi2 (n-1) + fromRational (2 * product [1..n] % product [3,5..2*n+1])

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

aproximaPi3 :: Integer -> Double
aproximaPi3 n = 
  fromRational (2 * (1 + sum (zipWith (%) numeradores (genericTake n denominadores))))

-- numeradores es la sucesión de los numeradores. Por ejemplo,
--    λ> take 10 numeradores
--    [1,2,6,24,120,720,5040,40320,362880,3628800]
numeradores :: [Integer]
numeradores = scanl (*) 1 [2..]

-- denominadores es la sucesión de los denominadores. Por ejemplo,
--    λ> take 10 denominadores
--    [3,15,105,945,10395,135135,2027025,34459425,654729075,13749310575]
denominadores :: [Integer]
denominadores = scanl (*) 3 [5, 7..]
 
-- 4ª solución
-- ===========

aproximaPi4 :: Integer -> Double
aproximaPi4 n = 
  read (x : "." ++ xs)
  where (x:xs) = show (aproximaPi4' n)

aproximaPi4' :: Integer -> Integer
aproximaPi4' n = 
  2 * (p + sum (zipWith div (map (*p) numeradores) (genericTake n denominadores))) 
  where p = 10^n

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

-- La propiedad es
prop_aproximaPi :: Property
prop_aproximaPi =
  forAll (arbitrary `suchThat` (> 3)) $ \n ->
  all (=~ aproximaPi1 n)
      [aproximaPi2 n,
       aproximaPi3 n,
       aproximaPi4 n]

(=~) :: Double -> Double -> Bool
x =~ y = abs (x - y) < 0.001

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

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

-- La comparación es
--    λ> aproximaPi1 3000 
--    3.141592653589793
--    (4.96 secs, 27,681,824,408 bytes)
--    λ> aproximaPi2 3000 
--    3.1415926535897922
--    (3.00 secs, 20,496,194,496 bytes)
--    λ> aproximaPi3 3000 
--    3.141592653589793
--    (3.13 secs, 13,439,528,432 bytes)
--    λ> aproximaPi4 3000 
--    3.141592653589793
--    (0.09 secs, 23,142,144 bytes)

import Data.Ratio ((%))

import Data.List (genericTake)

import Test.QuickCheck (Property, arbitrary, forAll, suchThat, quickCheck)

-- 1ª solución

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

aproximaPi1 :: Integer -> Double

aproximaPi1 n =

fromRational (2 * sum [product [1..i] % product [1,3..2*i+1] | i <- [0..n]])

-- 2ª solución

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

aproximaPi2 :: Integer -> Double

aproximaPi2 0 = 2

aproximaPi2 n =

aproximaPi2 (n-1) + fromRational (2 * product [1..n] % product [3,5..2*n+1])

-- 3ª solución

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

aproximaPi3 :: Integer -> Double

aproximaPi3 n =

fromRational (2 * (1 + sum (zipWith (%) numeradores (genericTake n denominadores))))

-- numeradores es la sucesión de los numeradores. Por ejemplo,

-- λ> take 10 numeradores

-- [1,2,6,24,120,720,5040,40320,362880,3628800]

numeradores :: [Integer]

numeradores = scanl (*) 1 [2..]

-- denominadores es la sucesión de los denominadores. Por ejemplo,

-- λ> take 10 denominadores

-- [3,15,105,945,10395,135135,2027025,34459425,654729075,13749310575]

denominadores :: [Integer]

denominadores = scanl (*) 3 [5, 7..]

-- 4ª solución

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

aproximaPi4 :: Integer -> Double

aproximaPi4 n =

read (x : "." ++ xs)

where (x:xs) = show (aproximaPi4' n)

aproximaPi4' :: Integer -> Integer

aproximaPi4' n =

2 * (p + sum (zipWith div (map (*p) numeradores) (genericTake n denominadores)))

where p = 10^n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_aproximaPi :: Property

prop_aproximaPi =

forAll (arbitrary `suchThat` (> 3)) $ \n ->

all (=~ aproximaPi1 n)

[aproximaPi2 n,

aproximaPi3 n,

aproximaPi4 n]

(=~) :: Double -> Double -> Bool

x =~ y = abs (x - y) < 0.001

-- La comprobación es

-- λ> quickCheck prop_aproximaPi

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> aproximaPi1 3000

-- 3.141592653589793

-- (4.96 secs, 27,681,824,408 bytes)

-- λ> aproximaPi2 3000

-- 3.1415926535897922

-- (3.00 secs, 20,496,194,496 bytes)

-- λ> aproximaPi3 3000

-- 3.141592653589793

-- (3.13 secs, 13,439,528,432 bytes)

-- λ> aproximaPi4 3000

-- 3.141592653589793

-- (0.09 secs, 23,142,144 bytes)

El código se encuentra en GitHub.

Pandigitales primos

PorJosé A. Alonso 29-junio-202221-junio-2022

Un número con n dígitos es pandigital si contiene todos los dígitos del 1 a n exactamente una vez. Por ejemplo, 2143 es un pandigital con 4 dígitos y, además, es primo.

Definir la lista

   pandigitalesPrimos :: [Int]

1	pandigitalesPrimos :: [Int]

tal que sus elementos son los números pandigitales primos, ordenados de mayor a menor. Por ejemplo,

   take 3 pandigitalesPrimos       ==  [7652413,7642513,7641253]
   2143 `elem` pandigitalesPrimos  ==  True
   length pandigitalesPrimos       ==  538

take 3 pandigitalesPrimos == [7652413,7642513,7641253]

2143 `elem` pandigitalesPrimos == True

length pandigitalesPrimos == 538

Soluciones

import Data.List (permutations, sort)
import Data.Char (intToDigit)
import Data.Numbers.Primes (isPrime)

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

pandigitalesPrimos1 :: [Int]
pandigitalesPrimos1 =
  concatMap nPandigitalesPrimos1 [9,8..1]

-- (nPandigitalesPrimos n) es la lista de los números pandigitales con n
-- dígitos, ordenada de mayor a menor. Por ejemplo,
--    nPandigitalesPrimos 4  ==  [4231,2341,2143,1423]
--    nPandigitalesPrimos 5  ==  []
nPandigitalesPrimos1 :: Int -> [Int]
nPandigitalesPrimos1 n = filter isPrime (pandigitales n)

-- (pandigitales n) es la lista de los números pandigitales de n dígitos
-- ordenada de mayor a menor. Por ejemplo,
--    pandigitales 3  ==  [321,312,231,213,132,123]
pandigitales :: Int -> [Int]
pandigitales n = 
    reverse $ sort $ map digitosAentero (permutations [1..n])

-- (digitosAentero ns) es el número cuyos dígitos son ns. Por ejemplo,
--    digitosAentero [3,2,5]  ==  325
digitosAentero :: [Int] -> Int
digitosAentero = read . map intToDigit

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

pandigitalesPrimos2 :: [Int]
pandigitalesPrimos2 =
  concatMap nPandigitalesPrimos2 [9,8..1]

-- Nota. La definición de nPandigitalesPrimos1 se puede simplificar, ya
-- que la suma de los números de 1 a n es divisible por 3, entonces los
-- números  pandigitales con n dígitos también lo son y, por tanto, no
-- son primos. 
nPandigitalesPrimos2 :: Int -> [Int]
nPandigitalesPrimos2 n 
  | sum [1..n] `mod` 3 == 0 = []
  | otherwise               = filter isPrime (pandigitales n)

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

pandigitalesPrimos3 :: [Int]
pandigitalesPrimos3 =
  concatMap nPandigitalesPrimos3 [9,8..1]

-- La definición de nPandigitales se puede simplificar, ya que
--    λ> [n | n <- [1..9], sum [1..n] `mod` 3 /= 0]
--    [1,4,7]
nPandigitalesPrimos3 :: Int -> [Int]
nPandigitalesPrimos3 n 
  | n `elem` [4,7] = filter isPrime (pandigitales n)
  | otherwise      = []

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

-- La propiedad es
prop_pandigitalesPrimos :: Bool
prop_pandigitalesPrimos =
  all (== pandigitalesPrimos1)
      [pandigitalesPrimos2,
       pandigitalesPrimos3]

-- La comprobación es
--    λ> prop_pandigitalesPrimos
--    True

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

-- La comparación es
--    λ> length (pandigitalesPrimos1)
--    538
--    (1.44 secs, 5,249,850,032 bytes)
--    λ> length (pandigitalesPrimos2)
--    538
--    (0.14 secs, 619,249,632 bytes)
--    λ> length (pandigitalesPrimos3)
--    538
--    (0.14 secs, 619,237,464 bytes)

import Data.List (permutations, sort)

import Data.Char (intToDigit)

import Data.Numbers.Primes (isPrime)

-- 1ª solución

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

pandigitalesPrimos1 :: [Int]

pandigitalesPrimos1 =

concatMap nPandigitalesPrimos1 [9,8..1]

-- (nPandigitalesPrimos n) es la lista de los números pandigitales con n

-- dígitos, ordenada de mayor a menor. Por ejemplo,

-- nPandigitalesPrimos 4 == [4231,2341,2143,1423]

-- nPandigitalesPrimos 5 == []

nPandigitalesPrimos1 :: Int -> [Int]

nPandigitalesPrimos1 n = filter isPrime (pandigitales n)

-- (pandigitales n) es la lista de los números pandigitales de n dígitos

-- ordenada de mayor a menor. Por ejemplo,

-- pandigitales 3 == [321,312,231,213,132,123]

pandigitales :: Int -> [Int]

pandigitales n =

reverse $ sort $ map digitosAentero (permutations [1..n])

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

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

digitosAentero :: [Int] -> Int

digitosAentero = read . map intToDigit

-- 2ª solución

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

pandigitalesPrimos2 :: [Int]

pandigitalesPrimos2 =

concatMap nPandigitalesPrimos2 [9,8..1]

-- Nota. La definición de nPandigitalesPrimos1 se puede simplificar, ya

-- que la suma de los números de 1 a n es divisible por 3, entonces los

-- números pandigitales con n dígitos también lo son y, por tanto, no

-- son primos.

nPandigitalesPrimos2 :: Int -> [Int]

nPandigitalesPrimos2 n

| sum [1..n] `mod` 3 == 0 = []

| otherwise = filter isPrime (pandigitales n)

-- 2ª solución

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

pandigitalesPrimos3 :: [Int]

pandigitalesPrimos3 =

concatMap nPandigitalesPrimos3 [9,8..1]

-- La definición de nPandigitales se puede simplificar, ya que

-- λ> [n | n <- [1..9], sum [1..n] `mod` 3 /= 0]

-- [1,4,7]

nPandigitalesPrimos3 :: Int -> [Int]

nPandigitalesPrimos3 n

| n `elem` [4,7] = filter isPrime (pandigitales n)

| otherwise = []

-- Comprobación de equivalencia

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

-- La propiedad es

prop_pandigitalesPrimos :: Bool

prop_pandigitalesPrimos =

all (== pandigitalesPrimos1)

[pandigitalesPrimos2,

pandigitalesPrimos3]

-- La comprobación es

-- λ> prop_pandigitalesPrimos

-- True

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (pandigitalesPrimos1)

-- 538

-- (1.44 secs, 5,249,850,032 bytes)

-- λ> length (pandigitalesPrimos2)

-- 538

-- (0.14 secs, 619,249,632 bytes)

-- λ> length (pandigitalesPrimos3)

-- 538

-- (0.14 secs, 619,237,464 bytes)

El código se encuentra en GitHub.

Codificación de Fibonacci

PorJosé A. Alonso 28-junio-202221-junio-2022

La codificación de Fibonacci http://bit.ly/1Lllqjv de un número n es una cadena d = d(0)d(1)…d(k-1)d(k) de ceros y unos tal que

   n = d(0)·F(2) + d(1)·F(3) +...+ d(k-1)·F(k+1) 
   d(k-1) = d(k) = 1

1 2	n = d(0)·F(2) + d(1)·F(3) +...+ d(k-1)·F(k+1) d(k-1) = d(k) = 1

donde F(i) es el i-ésimo término de la sucesión de Fibonacci

   0, 1, 1, 2, 3, 5, 8, 13, 21, 34, ...

1	0, 1, 1, 2, 3, 5, 8, 13, 21, 34, ...

Por ejemplo, la codificación de Fibonacci de 4 es «1011» ya que los dos últimos elementos son iguales a 1 y

   1·F(2) + 0·F(3) + 1·F(4) = 1·1 + 0·2 + 1·3 = 4

1	1·F(2) + 0·F(3) + 1·F(4) = 1·1 + 0·2 + 1·3 = 4

La codificación de Fibonacci de los primeros números se muestra en la siguiente tabla

    1  = 1     = F(2)           ≡       11
    2  = 2     = F(3)           ≡      011
    3  = 3     = F(4)           ≡     0011
    4  = 1+3   = F(2)+F(4)      ≡     1011
    5  = 5     = F(5)           ≡    00011
    6  = 1+5   = F(2)+F(5)      ≡    10011
    7  = 2+5   = F(3)+F(5)      ≡    01011
    8  = 8     = F(6)           ≡   000011
    9  = 1+8   = F(2)+F(6)      ≡   100011
   10  = 2+8   = F(3)+F(6)      ≡   010011
   11  = 3+8   = F(4)+F(6)      ≡   001011
   12  = 1+3+8 = F(2)+F(4)+F(6) ≡   101011
   13  = 13    = F(7)           ≡  0000011
   14  = 1+13  = F(2)+F(7)      ≡  1000011

1 = 1 = F(2) ≡ 11

2 = 2 = F(3) ≡ 011

3 = 3 = F(4) ≡ 0011

4 = 1+3 = F(2)+F(4) ≡ 1011

5 = 5 = F(5) ≡ 00011

6 = 1+5 = F(2)+F(5) ≡ 10011

7 = 2+5 = F(3)+F(5) ≡ 01011

8 = 8 = F(6) ≡ 000011

9 = 1+8 = F(2)+F(6) ≡ 100011

10 = 2+8 = F(3)+F(6) ≡ 010011

11 = 3+8 = F(4)+F(6) ≡ 001011

12 = 1+3+8 = F(2)+F(4)+F(6) ≡ 101011

13 = 13 = F(7) ≡ 0000011

14 = 1+13 = F(2)+F(7) ≡ 1000011

Definir la función

   codigoFib :: Integer -> String

1	codigoFib :: Integer -> String

tal que (codigoFib n) es la codificación de Fibonacci del número n. Por ejemplo,

   λ> codigoFib 65
   "0100100011"
   λ> [codigoFib n | n <- [1..7]]
   ["11","011","0011","1011","00011","10011","01011"]

λ> codigoFib 65

"0100100011"

λ> [codigoFib n | n <- [1..7]]

["11","011","0011","1011","00011","10011","01011"]

Comprobar con QuickCheck las siguientes propiedades:

Todo entero positivo se puede descomponer en suma de números de la sucesión de Fibonacci.
Las codificaciones de Fibonacci tienen como mínimo 2 elementos.
En las codificaciones de Fibonacci, la cadena «11» sólo aparece una vez y la única vez que aparece es al final.

Soluciones

import Data.List (isInfixOf)
import Data.Array (Array, accumArray, elems)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

codigoFib1 :: Integer -> String
codigoFib1 = concatMap show . codificaFibLista

-- (codificaFibLista n) es la lista correspondiente a la codificación de
-- Fibonacci del número n. Por ejemplo,
--    λ> codificaFibLista 65
--    [0,1,0,0,1,0,0,0,1,1]
--    λ> [codificaFibLista n | n <- [1..7]]
--    [[1,1],[0,1,1],[0,0,1,1],[1,0,1,1],[0,0,0,1,1],[1,0,0,1,1],[0,1,0,1,1]]
codificaFibLista :: Integer -> [Integer]
codificaFibLista n = map f [2..head xs] ++ [1]
  where xs = map fst (descomposicion n)
        f i | i `elem` xs = 1
            | otherwise = 0

-- (descomposicion n) es la lista de pares (i,f) tales que f es el
-- i-ésimo número de Fibonacci y las segundas componentes es una
-- sucesión decreciente de números de Fibonacci cuya suma es n. Por
-- ejemplo, 
--    descomposicion 65  ==  [(10,55),(6,8),(3,2)]
--    descomposicion 66  ==  [(10,55),(6,8),(4,3)]
descomposicion :: Integer -> [(Integer, Integer)]
descomposicion 0 = []
descomposicion 1 = [(2,1)]
descomposicion n = (i,x) : descomposicion (n-x)
  where (i,x) = fibAnterior n

-- (fibAnterior n) es el mayor número de Fibonacci menor o igual que
-- n. Por ejemplo,
--    fibAnterior 33  ==  (8,21)
--    fibAnterior 34  ==  (9,34)
fibAnterior :: Integer -> (Integer, Integer)
fibAnterior n = last (takeWhile p fibsConIndice)
  where p (_,x) = x <= n

-- fibsConIndice es la sucesión de los números de Fibonacci junto con
-- sus índices. Por ejemplo,
--    λ> take 10 fibsConIndice
--    [(0,0),(1,1),(2,1),(3,2),(4,3),(5,5),(6,8),(7,13),(8,21),(9,34)]
fibsConIndice :: [(Integer, Integer)]
fibsConIndice = zip [0..] fibs

-- fibs es la sucesión de Fibonacci. Por ejemplo, 
--    take 10 fibs  ==  [0,1,1,2,3,5,8,13,21,34]
fibs :: [Integer]
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

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

codigoFib2 :: Integer -> String
codigoFib2 = concatMap show . elems . codificaFibVec

-- (codificaFibVec n) es el vector correspondiente a la codificación de
-- Fibonacci del número n. Por ejemplo,
--    λ> codificaFibVec 65
--    array (0,9) [(0,0),(1,1),(2,0),(3,0),(4,1),(5,0),(6,0),(7,0),(8,1),(9,1)]
--    λ> [elems (codificaFibVec n) | n <- [1..7]]
--    [[1,1],[0,1,1],[0,0,1,1],[1,0,1,1],[0,0,0,1,1],[1,0,0,1,1],[0,1,0,1,1]]
codificaFibVec :: Integer -> Array Integer Integer
codificaFibVec n = accumArray (+) 0 (0,a+1) ((a+1,1):is) 
  where is = [(i-2,1) | (i,_) <- descomposicion n]
        a  = fst (head is)

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

-- La propiedad es
prop_codigoFib :: Positive Integer -> Bool 
prop_codigoFib (Positive n) =
  codigoFib1 n == codigoFib2 n
  
-- La comprobación es
--    λ> quickCheck prop_codigoFib
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> head [n | n <- [1..], length (codigoFib1 n) > 25]
--    121393
--    (4.30 secs, 3,031,108,104 bytes)
--    λ> head [n | n <- [1..], length (codigoFib2 n) > 25]
--    121393
--    (3.46 secs, 2,505,869,616 bytes)

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

-- Usaremos la 2ª definición
codigoFib :: Integer -> String
codigoFib = codigoFib2

-- Prop.: La función descomposicion es correcta:
prop_descomposicion_correcta :: Positive Integer -> Bool
prop_descomposicion_correcta (Positive n) =
  n == sum (map snd (descomposicion n))

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

-- Prop.: Todo entero positivo se puede descomponer en suma de números de
-- la sucesión de Fibonacci.
prop_descomposicion :: Positive Integer -> Bool
prop_descomposicion (Positive n) =
  not (null (descomposicion n))

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

-- Prop.: Las codificaciones de Fibonacci tienen como mínimo 2 elementos.
prop_length_codigoFib :: Positive Integer -> Bool
prop_length_codigoFib (Positive n) =
  length (codigoFib n) >= 2

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

-- Prop.: En las codificaciones de Fibonacci, la cadena "11" sólo
-- aparece una vez y la única vez que aparece es al final.
prop3_cadena_11_en_codigoFib :: Positive Integer -> Bool
prop3_cadena_11_en_codigoFib (Positive n) = 
  take 2 xs == "11" && not ("11" `isInfixOf` drop 2 xs)
  where xs = reverse (codigoFib n)

-- La comprobación es
--    λ> quickCheck prop3_cadena_11_en_codigoFib
--    +++ 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

import Data.List (isInfixOf)

import Data.Array (Array, accumArray, elems)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

codigoFib1 :: Integer -> String

codigoFib1 = concatMap show . codificaFibLista

-- (codificaFibLista n) es la lista correspondiente a la codificación de

-- Fibonacci del número n. Por ejemplo,

-- λ> codificaFibLista 65

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

-- λ> [codificaFibLista n | n <- [1..7]]

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

codificaFibLista :: Integer -> [Integer]

codificaFibLista n = map f [2..head xs] ++ [1]

where xs = map fst (descomposicion n)

f i | i `elem` xs = 1

| otherwise = 0

-- (descomposicion n) es la lista de pares (i,f) tales que f es el

-- i-ésimo número de Fibonacci y las segundas componentes es una

-- sucesión decreciente de números de Fibonacci cuya suma es n. Por

-- ejemplo,

-- descomposicion 65 == [(10,55),(6,8),(3,2)]

-- descomposicion 66 == [(10,55),(6,8),(4,3)]

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

descomposicion 0 = []

descomposicion 1 = [(2,1)]

descomposicion n = (i,x) : descomposicion (n-x)

where (i,x) = fibAnterior n

-- (fibAnterior n) es el mayor número de Fibonacci menor o igual que

-- n. Por ejemplo,

-- fibAnterior 33 == (8,21)

-- fibAnterior 34 == (9,34)

fibAnterior :: Integer -> (Integer, Integer)

fibAnterior n = last (takeWhile p fibsConIndice)

where p (_,x) = x <= n

-- fibsConIndice es la sucesión de los números de Fibonacci junto con

-- sus índices. Por ejemplo,

-- λ> take 10 fibsConIndice

-- [(0,0),(1,1),(2,1),(3,2),(4,3),(5,5),(6,8),(7,13),(8,21),(9,34)]

fibsConIndice :: [(Integer, Integer)]

fibsConIndice = zip [0..] fibs

-- fibs es la sucesión de Fibonacci. Por ejemplo,

-- take 10 fibs == [0,1,1,2,3,5,8,13,21,34]

fibs :: [Integer]

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

--- 2ª solución

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

codigoFib2 :: Integer -> String

codigoFib2 = concatMap show . elems . codificaFibVec

-- (codificaFibVec n) es el vector correspondiente a la codificación de

-- Fibonacci del número n. Por ejemplo,

-- λ> codificaFibVec 65

-- array (0,9) [(0,0),(1,1),(2,0),(3,0),(4,1),(5,0),(6,0),(7,0),(8,1),(9,1)]

-- λ> [elems (codificaFibVec n) | n <- [1..7]]

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

codificaFibVec :: Integer -> Array Integer Integer

codificaFibVec n = accumArray (+) 0 (0,a+1) ((a+1,1):is)

where is = [(i-2,1) | (i,_) <- descomposicion n]

a = fst (head is)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_codigoFib :: Positive Integer -> Bool

prop_codigoFib (Positive n) =

codigoFib1 n == codigoFib2 n

-- La comprobación es

-- λ> quickCheck prop_codigoFib

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> head [n | n <- [1..], length (codigoFib1 n) > 25]

-- 121393

-- (4.30 secs, 3,031,108,104 bytes)

-- λ> head [n | n <- [1..], length (codigoFib2 n) > 25]

-- 121393

-- (3.46 secs, 2,505,869,616 bytes)

-- Propiedades

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

-- Usaremos la 2ª definición

codigoFib :: Integer -> String

codigoFib = codigoFib2

-- Prop.: La función descomposicion es correcta:

prop_descomposicion_correcta :: Positive Integer -> Bool

prop_descomposicion_correcta (Positive n) =

n == sum (map snd (descomposicion n))

-- La comprobación es

-- λ> quickCheck prop_descomposicion_correcta

-- +++ OK, passed 100 tests.

-- Prop.: Todo entero positivo se puede descomponer en suma de números de

-- la sucesión de Fibonacci.

prop_descomposicion :: Positive Integer -> Bool

prop_descomposicion (Positive n) =

not (null (descomposicion n))

-- La comprobación es

-- λ> quickCheck prop_descomposicion

-- +++ OK, passed 100 tests.

-- Prop.: Las codificaciones de Fibonacci tienen como mínimo 2 elementos.

prop_length_codigoFib :: Positive Integer -> Bool

prop_length_codigoFib (Positive n) =

length (codigoFib n) >= 2

-- La comprobación es

-- λ> quickCheck prop_length_codigoFib

-- +++ OK, passed 100 tests.

-- Prop.: En las codificaciones de Fibonacci, la cadena "11" sólo

-- aparece una vez y la única vez que aparece es al final.

prop3_cadena_11_en_codigoFib :: Positive Integer -> Bool

prop3_cadena_11_en_codigoFib (Positive n) =

take 2 xs == "11" && not ("11" `isInfixOf` drop 2 xs)

where xs = reverse (codigoFib n)

-- La comprobación es

-- λ> quickCheck prop3_cadena_11_en_codigoFib

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

La sucesión del reloj astronómico de Praga

PorJosé A. Alonso 27-junio-202221-junio-2022

La cadena infinita «1234321234321234321…», formada por la repetición de los dígitos 123432, tiene una propiedad (en la que se basa el funcionamiento del reloj astronómico de Praga: la cadena se puede partir en una sucesión de números, de forma que la suma de los dígitos de dichos números es la sucesión de los números naturales, como se observa a continuación:

    1, 2, 3, 4, 32, 123, 43, 2123, 432, 1234, 32123, ...
    1, 2, 3, 4,  5,   6,  7,    8,   9,   10,    11, ...

1 2	1, 2, 3, 4, 32, 123, 43, 2123, 432, 1234, 32123, ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ...

Definir la lista

   reloj :: [Integer]

1	reloj :: [Integer]

cuyos elementos son los términos de la sucesión anterior. Por ejemplo,

   λ> take 11 reloj
   [1,2,3,4,32,123,43,2123,432,1234,32123]
   λ> (reloj !! 1000) `mod` (10^50)
   23432123432123432123432123432123432123432123432123

λ> take 11 reloj

[1,2,3,4,32,123,43,2123,432,1234,32123]

λ> (reloj !! 1000) `mod` (10^50)

23432123432123432123432123432123432123432123432123

Soluciones

import Data.List (inits, tails)
import Data.Char (digitToInt)
import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

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

reloj1 :: [Integer]
reloj1 = aux [1..] (cycle "123432")
  where aux (n:ns) xs = read ys : aux ns zs
          where (ys,zs) = prefijoSuma n xs

-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo
-- suma es n y el resto de xs. Por ejemplo,
--    prefijoSuma 6 "12343"  ==  ("123","43")
prefijoSuma :: Int -> String -> (String,String)
prefijoSuma n xs = 
  head [(us,vs) | (us,vs) <- zip (inits xs) (tails xs)
                , sumaD us == n]

-- (sumaD xs) es la suma de los dígitos de xs. Por ejemplo,
--    sumaD "123"  ==  6
sumaD :: String -> Int
sumaD = sum . map digitToInt

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

reloj2 :: [Integer]
reloj2 = aux [1..] (cycle "123432")
  where aux (n:ns) xs = read ys : aux ns zs
          where (ys,zs) = prefijoSuma2 n xs

-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo
-- suma es n y el resto de xs. Por ejemplo,
--    prefijoSuma2 6 "12343"  ==  ("123","43")
prefijoSuma2 :: Int -> String -> (String,String)
prefijoSuma2 n (x:xs) 
  | y == n = ([x],xs)
  | otherwise = (x:ys,zs) 
  where y       = read [x]
        (ys,zs) = prefijoSuma2 (n-y) xs

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

-- La propiedad es
prop_reloj :: NonNegative Int -> Bool
prop_reloj (NonNegative n) =
  reloj1 !! n == reloj2 !! n
  
-- La comprobación es
--    λ> quickCheck prop_reloj
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> (reloj1 !! 1000) `mod` (10^9)
--    123432123
--    (2.47 secs, 5,797,620,784 bytes)
--    λ> (reloj2 !! 1000) `mod` (10^9)
--    123432123
--    (0.44 secs, 798,841,528 bytes)

import Data.List (inits, tails)

import Data.Char (digitToInt)

import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

-- 1ª solución

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

reloj1 :: [Integer]

reloj1 = aux [1..] (cycle "123432")

where aux (n:ns) xs = read ys : aux ns zs

where (ys,zs) = prefijoSuma n xs

-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo

-- suma es n y el resto de xs. Por ejemplo,

-- prefijoSuma 6 "12343" == ("123","43")

prefijoSuma :: Int -> String -> (String,String)

prefijoSuma n xs =

head [(us,vs) | (us,vs) <- zip (inits xs) (tails xs)

, sumaD us == n]

-- (sumaD xs) es la suma de los dígitos de xs. Por ejemplo,

-- sumaD "123" == 6

sumaD :: String -> Int

sumaD = sum . map digitToInt

-- 2ª solución

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

reloj2 :: [Integer]

reloj2 = aux [1..] (cycle "123432")

where aux (n:ns) xs = read ys : aux ns zs

where (ys,zs) = prefijoSuma2 n xs

-- (prefijoSuma n xs) es el par formado por el primer prefijo de xs cuyo

-- suma es n y el resto de xs. Por ejemplo,

-- prefijoSuma2 6 "12343" == ("123","43")

prefijoSuma2 :: Int -> String -> (String,String)

prefijoSuma2 n (x:xs)

| y == n = ([x],xs)

| otherwise = (x:ys,zs)

where y = read [x]

(ys,zs) = prefijoSuma2 (n-y) xs

-- Comprobación de equivalencia

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

-- La propiedad es

prop_reloj :: NonNegative Int -> Bool

prop_reloj (NonNegative n) =

reloj1 !! n == reloj2 !! n

-- La comprobación es

-- λ> quickCheck prop_reloj

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> (reloj1 !! 1000) `mod` (10^9)

-- 123432123

-- (2.47 secs, 5,797,620,784 bytes)

-- λ> (reloj2 !! 1000) `mod` (10^9)

-- 123432123

-- (0.44 secs, 798,841,528 bytes)

El código se encuentra en GitHub.

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

PorJosé A. Alonso 24-junio-202221-junio-2022

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

import Test.QuickCheck (Arbitrary, Gen, Property, (==>), arbitrary, choose, sized, quickCheck)

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

biseccion1 :: (Double -> Double) -> Double -> Double -> Double -> Double
biseccion1 f a b e  
    | abs (f c) < e = c
    | f a * f c < 0 = biseccion1 f a c e
    | otherwise     = biseccion1 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

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

newtype Polinomio = P [Int]
  deriving Show

valorPolinomio :: Polinomio -> Double -> Double
valorPolinomio (P cs) x =
  sum [fromIntegral c * x^n | (c,n) <- zip cs [0..]]

polinomioArbitrario :: Int -> Gen Polinomio
polinomioArbitrario 0 = return (P [])
polinomioArbitrario n = do
  c <- choose (-10,10)
  (P xs) <- polinomioArbitrario (n `div` 2)
  return (P (c:xs))

instance Arbitrary Polinomio where
  arbitrary = sized polinomioArbitrario

-- La propiedad es
prop_biseccion :: Polinomio -> Double -> Double -> Double -> Property
prop_biseccion p a b e =
  f a * f b < 0 && e > 0 ==>
  biseccion1 f a b e =~ biseccion2 f a b e
  where
    f = valorPolinomio p 
    x =~ y = abs (x - y) < 0.001

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxDiscardRatio=30}) prop_biseccion
--    +++ OK, passed 100 tests; 2156 discarded.

import Test.QuickCheck (Arbitrary, Gen, Property, (==>), arbitrary, choose, sized, quickCheck)

-- 1ª solución

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

biseccion1 :: (Double -> Double) -> Double -> Double -> Double -> Double

biseccion1 f a b e

| abs (f c) < e = c

| f a * f c < 0 = biseccion1 f a c e

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

-- Comprobación de equivalencia

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

newtype Polinomio = P [Int]

deriving Show

valorPolinomio :: Polinomio -> Double -> Double

valorPolinomio (P cs) x =

sum [fromIntegral c * x^n | (c,n) <- zip cs [0..]]

polinomioArbitrario :: Int -> Gen Polinomio

polinomioArbitrario 0 = return (P [])

polinomioArbitrario n = do

c <- choose (-10,10)

(P xs) <- polinomioArbitrario (n `div` 2)

return (P (c:xs))

instance Arbitrary Polinomio where

arbitrary = sized polinomioArbitrario

-- La propiedad es

prop_biseccion :: Polinomio -> Double -> Double -> Double -> Property

prop_biseccion p a b e =

f a * f b < 0 && e > 0 ==>

biseccion1 f a b e =~ biseccion2 f a b e

where

f = valorPolinomio p

x =~ y = abs (x - y) < 0.001

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxDiscardRatio=30}) prop_biseccion

-- +++ OK, passed 100 tests; 2156 discarded.

El código se encuentra en GitHub.

Números para los que mcm(1,2,…n-1) = mcm(1,2,…,n)

PorJosé A. Alonso 23-junio-202221-junio-2022

Un número n es especial si mcm(1,2,…,n-1) = mcm(1,2,…,n). Por ejemplo, el 6 es especial ya que

   mcm(1,2,3,4,5) = 60 = mcm(1,2,3,4,5,6)

1	mcm(1,2,3,4,5) = 60 = mcm(1,2,3,4,5,6)

Definir la sucesión

   especiales :: [Integer]

1	especiales :: [Integer]

cuyos términos son los números especiales. Por ejemplo,

   take 10 especiales     ==  [1,6,10,12,14,15,18,20,21,22]
   especiales !! 50       ==  84
   especiales !! 500      ==  638
   especiales !! 5000     ==  5806
   especiales !! 50000    ==  55746

take 10 especiales == [1,6,10,12,14,15,18,20,21,22]

especiales !! 50 == 84

especiales !! 500 == 638

especiales !! 5000 == 5806

especiales !! 50000 == 55746

Soluciones

import Test.QuickCheck (NonNegative(NonNegative), quickCheck)

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

especiales1 :: [Integer]
especiales1 = filter especial1 [1..]

especial1 :: Integer -> Bool
especial1 n = mcm1 [1..n-1] == mcm1 [1..n]

mcm1 :: [Integer] -> Integer
mcm1 []     = 1
mcm1 (x:xs) = lcm x (mcm1 xs)

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

especiales2 :: [Integer]
especiales2 = filter especial2 [1..]

especial2 :: Integer -> Bool
especial2 n = mcm2 [1..n-1] == mcm2 [1..n]

mcm2 :: [Integer] -> Integer
mcm2 = foldr lcm 1

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

especiales3 :: [Integer]
especiales3 = [n | ((n,x),(_,y)) <- zip mcms (tail mcms)
                 , x == y]

mcms :: [(Integer,Integer)]
mcms = zip [1..] (scanl lcm 1 [1..])

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

-- La propiedad es
prop_especiales :: NonNegative Int -> Bool
prop_especiales (NonNegative n) =
  all (== especiales1 !! n)
      [especiales2 !! n,
       especiales3 !! n]

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

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

-- La comparación es

-- Comparación 
--    λ> especiales1 !! 2000
--    2390
--    (3.38 secs, 4,724,497,192 bytes)
--    λ> especiales2 !! 2000
--    2390
--    (1.91 secs, 4,303,415,512 bytes)
--    λ> especiales3 !! 2000
--    2390
--    (0.01 secs, 4,209,664 bytes)

import Test.QuickCheck (NonNegative(NonNegative), quickCheck)

-- 1ª solución

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

especiales1 :: [Integer]

especiales1 = filter especial1 [1..]

especial1 :: Integer -> Bool

especial1 n = mcm1 [1..n-1] == mcm1 [1..n]

mcm1 :: [Integer] -> Integer

mcm1 [] = 1

mcm1 (x:xs) = lcm x (mcm1 xs)

-- 2ª solución

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

especiales2 :: [Integer]

especiales2 = filter especial2 [1..]

especial2 :: Integer -> Bool

especial2 n = mcm2 [1..n-1] == mcm2 [1..n]

mcm2 :: [Integer] -> Integer

mcm2 = foldr lcm 1

-- 3ª solución

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

especiales3 :: [Integer]

especiales3 = [n | ((n,x),(_,y)) <- zip mcms (tail mcms)

, x == y]

mcms :: [(Integer,Integer)]

mcms = zip [1..] (scanl lcm 1 [1..])

-- Comprobación de equivalencia

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

-- La propiedad es

prop_especiales :: NonNegative Int -> Bool

prop_especiales (NonNegative n) =

all (== especiales1 !! n)

[especiales2 !! n,

especiales3 !! n]

-- La comprobación es

-- λ> quickCheck prop_especiales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- Comparación

-- λ> especiales1 !! 2000

-- 2390

-- (3.38 secs, 4,724,497,192 bytes)

-- λ> especiales2 !! 2000

-- 2390

-- (1.91 secs, 4,303,415,512 bytes)

-- λ> especiales3 !! 2000

-- 2390

-- (0.01 secs, 4,209,664 bytes)

El código se encuentra en GitHub.

Cálculo de la suma 11! + 22! + 33! + … + nn!

PorJosé A. Alonso 22-junio-202221-junio-2022

Definir la función

   suma :: Integer -> Integer

1	suma :: Integer -> Integer

tal que (suma n) es la suma 1·1! + 2·2! + 3·3! + ... + n·n!. Por ejemplo,

   suma 1  ==  1
   suma 2  ==  5
   suma 3  ==  23
   suma 4  ==  119
   suma 5  ==  719
   take 9 (show (suma 70000))  ==  "823780458"

suma 1 == 1

suma 2 == 5

suma 3 == 23

suma 4 == 119

suma 5 == 719

take 9 (show (suma 70000)) == "823780458"

Soluciones

import Test.QuickCheck (Positive (Positive), quickCheck)

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

suma1 :: Integer -> Integer
suma1 n = sum [k * factorial k | k <- [1..n]]

factorial :: Integer -> Integer
factorial n = product [1..n]

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

suma2 :: Integer -> Integer
suma2 n = sum (zipWith (*) [1..n] factoriales)

factoriales :: [Integer]
factoriales = scanl (*) 1 [2..]

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

-- Basada en los siguientes cálculos
--    λ> [suma1 n | n <- [0..10]]
--    [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]
--    λ> [factorial n | n <- [0..10]]
--    [1,1,2,6,24,120,720,5040,40320,362880,3628800]
--    λ> [factorial n | n <- [1..11]]
--    [1,2,6,24,120,720,5040,40320,362880,3628800,39916800]
--    λ> [factorial n - 1 | n <- [1..11]]
--    [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

suma3 :: Integer -> Integer
suma3 n = factorial (n+1) - 1

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

-- La propiedad es
prop_suma :: Positive Integer -> Bool
prop_suma (Positive n) =
  all (== suma1 n)
      [suma2 n,
       suma3 n]
  
-- La comprobación es
--    λ> quickCheck prop_suma
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> take 5 (show (suma1 4000))
--    "73170"
--    (5.04 secs, 16,225,195,448 bytes)
--    λ> take 5 (show (suma2 4000))
--    "73170"
--    (0.08 secs, 35,862,152 bytes)
--    λ> take 5 (show (suma3 4000))
--    "73170"
--    (0.01 secs, 12,896,968 bytes)
--    
--    
--    λ> take 5 (show (suma2 40000))
--    "83669"
--    (1.82 secs, 4,549,612,264 bytes)
--    λ> take 5 (show (suma3 40000))
--    "83669"
--    (0.24 secs, 1,620,976,984 bytes)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

suma1 :: Integer -> Integer

suma1 n = sum [k * factorial k | k <- [1..n]]

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª solución

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

suma2 :: Integer -> Integer

suma2 n = sum (zipWith (*) [1..n] factoriales)

factoriales :: [Integer]

factoriales = scanl (*) 1 [2..]

-- 3ª solución

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

-- Basada en los siguientes cálculos

-- λ> [suma1 n | n <- [0..10]]

-- [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

-- λ> [factorial n | n <- [0..10]]

-- [1,1,2,6,24,120,720,5040,40320,362880,3628800]

-- λ> [factorial n | n <- [1..11]]

-- [1,2,6,24,120,720,5040,40320,362880,3628800,39916800]

-- λ> [factorial n - 1 | n <- [1..11]]

-- [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

suma3 :: Integer -> Integer

suma3 n = factorial (n+1) - 1

-- Comprobación de equivalencia

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

-- La propiedad es

prop_suma :: Positive Integer -> Bool

prop_suma (Positive n) =

all (== suma1 n)

[suma2 n,

suma3 n]

-- La comprobación es

-- λ> quickCheck prop_suma

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> take 5 (show (suma1 4000))

-- "73170"

-- (5.04 secs, 16,225,195,448 bytes)

-- λ> take 5 (show (suma2 4000))

-- "73170"

-- (0.08 secs, 35,862,152 bytes)

-- λ> take 5 (show (suma3 4000))

-- "73170"

-- (0.01 secs, 12,896,968 bytes)

-- λ> take 5 (show (suma2 40000))

-- "83669"

-- (1.82 secs, 4,549,612,264 bytes)

-- λ> take 5 (show (suma3 40000))

-- "83669"

-- (0.24 secs, 1,620,976,984 bytes)

El código se encuentra en GitHub.

Cálculo aproximado de integrales definidas

PorJosé A. Alonso 21-junio-202221-junio-2022

La integral definida de una función f entre los límites a y b puede calcularse mediante la regla del rectángulo usando la fórmula

   h * (f(a+h/2) + f(a+h+h/2) + f(a+2h+h/2) + ... + f(a+n*h+h/2))

1	h * (f(a+h/2) + f(a+h+h/2) + f(a+2h+h/2) + ... + f(a+n*h+h/2))

con a+n*h+h/2 <= b < a+(n+1)*h+h/2 y usando valores pequeños para h.

Definir la función

   integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

1	integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

tal que (integral a b f h) es el valor de dicha expresión. Por ejemplo, el cálculo de la integral de f(x) = x^3 entre 0 y 1, con paso 0.01, es

   integral 0 1 (^3) 0.01  ==  0.24998750000000042

1	integral 0 1 (^3) 0.01 == 0.24998750000000042

Otros ejemplos son

   integral 0 1 (^4) 0.01                   ==  0.19998333362500048
   integral 0 1 (\x -> 3*x^2 + 4*x^3) 0.01  ==  1.9999250000000026
   log 2 - integral 1 2 (\x -> 1/x) 0.01         ==  3.124931644782336e-6
   pi - 4 * integral 0 1 (\x -> 1/(x^2+1)) 0.01  ==  -8.333333331389525e-6

integral 0 1 (^4) 0.01 == 0.19998333362500048

integral 0 1 (\x -> 3*x^2 + 4*x^3) 0.01 == 1.9999250000000026

log 2 - integral 1 2 (\x -> 1/x) 0.01 == 3.124931644782336e-6

pi - 4 * integral 0 1 (\x -> 1/(x^2+1)) 0.01 == -8.333333331389525e-6

Soluciones

import Test.QuickCheck.HigherOrder (quickCheck')
import Test.QuickCheck (Property, (==>), quickCheck)

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

integral1 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral1 a b f h
  | a+h/2 > b = 0
  | otherwise = h * f (a+h/2) + integral1 (a+h) b f h

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

integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral2 a b f h = aux a where
  aux x | x+h/2 > b = 0
        | otherwise = h * f (x+h/2) + aux (x+h)

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

integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral3 a b f h = h * suma (a+h/2) b (+h) f

-- (suma a b s f) es l valor de
--    f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))
-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,
--    suma 2 5 (1+) (^3)  ==  224
suma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a
suma a b s f = sum [f x | x <- sucesion a b s]

-- (sucesion x y s) es la lista
--    [a, s(a), s(s(a), ..., s(...(s(a))...)]
-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,
--    sucesion 3 20 (+2)  ==  [3,5,7,9,11,13,15,17,19]
sucesion :: Ord a => a -> a -> (a -> a) -> [a]
sucesion a b s = takeWhile (<=b) (iterate s a)

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

-- La propiedad es
prop_integral :: Int -> Int -> (Int -> Int) -> Int -> Property
prop_integral a b f h =
  a < b && h > 0 ==> 
  all (=~ integral1 a' b' f' h')
      [integral2 a' b' f' h',
       integral3 a' b' f' h']
  where
    a' = fromIntegral a
    b' = fromIntegral b
    h' = fromIntegral h
    f' = fromIntegral . f. round
    x =~ y = abs (x - y) < 0.001
        
-- La comprobación es
--    λ> quickCheck' prop_integral
--    +++ OK, passed 100 tests; 385 discarded.

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

-- La comparación es
--    λ> integral1 0 10 (^3) 0.00001
--    2499.999999881125
--    (2.63 secs, 1,491,006,744 bytes)
--    λ> integral2 0 10 (^3) 0.00001
--    2499.999999881125
--    (1.93 secs, 1,419,006,696 bytes)
--    λ> integral3 0 10 (^3) 0.00001
--    2499.9999998811422
--    (1.28 secs, 817,772,216 bytes)

import Test.QuickCheck.HigherOrder (quickCheck')

import Test.QuickCheck (Property, (==>), quickCheck)

-- 1ª solución

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

integral1 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral1 a b f h

| a+h/2 > b = 0

| otherwise = h * f (a+h/2) + integral1 (a+h) b f h

-- 2ª solución

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

integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral2 a b f h = aux a where

aux x | x+h/2 > b = 0

| otherwise = h * f (x+h/2) + aux (x+h)

-- 3ª solución

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

integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral3 a b f h = h * suma (a+h/2) b (+h) f

-- (suma a b s f) es l valor de

-- f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))

-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,

-- suma 2 5 (1+) (^3) == 224

suma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a

suma a b s f = sum [f x | x <- sucesion a b s]

-- (sucesion x y s) es la lista

-- [a, s(a), s(s(a), ..., s(...(s(a))...)]

-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,

-- sucesion 3 20 (+2) == [3,5,7,9,11,13,15,17,19]

sucesion :: Ord a => a -> a -> (a -> a) -> [a]

sucesion a b s = takeWhile (<=b) (iterate s a)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_integral :: Int -> Int -> (Int -> Int) -> Int -> Property

prop_integral a b f h =

a < b && h > 0 ==>

all (=~ integral1 a' b' f' h')

[integral2 a' b' f' h',

integral3 a' b' f' h']

where

a' = fromIntegral a

b' = fromIntegral b

h' = fromIntegral h

f' = fromIntegral . f. round

x =~ y = abs (x - y) < 0.001

-- La comprobación es

-- λ> quickCheck' prop_integral

-- +++ OK, passed 100 tests; 385 discarded.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> integral1 0 10 (^3) 0.00001

-- 2499.999999881125

-- (2.63 secs, 1,491,006,744 bytes)

-- λ> integral2 0 10 (^3) 0.00001

-- 2499.999999881125

-- (1.93 secs, 1,419,006,696 bytes)

-- λ> integral3 0 10 (^3) 0.00001

-- 2499.9999998811422

-- (1.28 secs, 817,772,216 bytes)

El código se encuentra en GitHub.

Menor número con una cantidad dada de divisores

PorJosé A. Alonso 20-junio-202221-junio-2022

El menor número con 2 divisores es el 2, ya que tiene 2 divisores (el 1 y el 2) y el anterior al 2 (el 1) sólo tiene 1 divisor (el 1).

El menor número con 4 divisores es el 6, ya que tiene 4 divisores (el 1, 2, 3 y 6) y sus anteriores (el 1, 2, 3, 4 y 5) tienen menos de 4 divisores (tienen 1, 1, 1, 3 y 1, respectivamente).

El menor número con 8 divisores es el 24, ya que tiene 8 divisores (el 1, 2, 3, 4, 6, 8, 12 y 24) y sus anteriores (del 1 al 23) tienen menos de 8 divisores.

El menor número con 16 divisores es el 120, ya que tiene 16 divisores (el 1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 24, 30, 40, 60 y 120) y sus anteriores (del 1 al 119) tienen menos de 16 divisores.

Definir la función

   menor :: Integer -> Integer

1	menor :: Integer -> Integer

tal que (menor n) es el menor número con 2^n divisores. Por ejemplo,

   menor 1  ==  2
   menor 2  ==  6
   menor 3  ==  24
   menor 4  ==  120
   length (show (menor (4*10^4)))  ==  207945

menor 1 == 2

menor 2 == 6

menor 3 == 24

menor 4 == 120

length (show (menor (4*10^4))) == 207945

Comprobar con QuickCheck que, para todo k >=0, (menor (2^k)) es un divisor de (menor (2^(k+1))).

Nota: Este ejercicio está basado en el problema N1 de la Olimpíada Internacional de Matemáticas (IMO) del 2011.

Soluciones

import Data.List           (genericLength, genericTake, group)
import Data.Numbers.Primes (primeFactors, primes)
import Test.QuickCheck     (Positive (Positive), maxSize, stdArgs,
                            quickCheck, quickCheckWith)

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

menor1 :: Integer -> Integer
menor1 n =
  head [x | x <- [1..], numeroDivisores x == 2^n]

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

-- (divisores x) es la lista de los divisores de x. Por ejemplo,
--    divisores 12  ==  [1,3,2,6,4,12]
--    divisores 25  ==  [1,5,25]
divisores :: Integer -> [Integer]
divisores n =
  [x | x <- [1..n], n `rem` x == 0]

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

menor2 :: Integer -> Integer
menor2 n =
  head [x | x <- [1..], numeroDivisores2 x == 2^n]

numeroDivisores2 :: Integer -> Integer
numeroDivisores2 =
  product . map ((+1) . genericLength) . group . primeFactors

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

menor3 :: Integer -> Integer
menor3 n = product (genericTake n potencias)

-- potencias es la sucesión de las potencias de la forma p^(2^k),
-- donde p es un número primo y k es un número natural, ordenadas de
-- menor a mayor. Por ejemplo,
--    take 14 potencias    ==  [2,3,4,5,7,9,11,13,16,17,19,23,25,29]
potencias :: [Integer]
potencias = 2 : mezcla (tail primes) (map (^2) potencias)

-- (mezcla xs ys) es la lista obtenida mezclando las dos listas xs e ys,
-- que se suponen ordenadas y disjuntas. Por ejemplo,
--    λ> take 15 (mezcla [2^n | n <- [1..]] [3^n | n <- [1..]])
--    [2,3,4,8,9,16,27,32,64,81,128,243,256,512,729]
mezcla :: Ord a => [a] -> [a] -> [a]
mezcla (x:xs) (y:ys) | x < y = x : mezcla xs (y:ys)
                     | x > y = y : mezcla (x:xs) ys

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

-- La propiedad es
prop_menor :: Positive Integer -> Bool
prop_menor (Positive n) =
  all (== menor1 n)
      [menor2 n,
       menor3 n]

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

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

-- La comparación es
--   λ> menor 8
--   1081080
--   (47.69 secs, 94,764,856,352 bytes)
--   λ> menor2 8
--   1081080
--   (36.17 secs, 94,764,856,368 bytes)
--   λ> menor3 8
--   1081080
--   (0.00 secs, 116,960 bytes)

-- Definición de menor
-- ===================

-- En lo que sigue, usaremos menor3 como menor.
menor :: Integer -> Integer
menor = menor3

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

-- La propiedad es
prop_menor_divide :: Positive Integer -> Bool
prop_menor_divide (Positive n) =
  menor (n+1) `mod` menor n == 0

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

100

101

102

103

import Data.List (genericLength, genericTake, group)

import Data.Numbers.Primes (primeFactors, primes)

import Test.QuickCheck (Positive (Positive), maxSize, stdArgs,

quickCheck, quickCheckWith)

-- 1ª solución

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

menor1 :: Integer -> Integer

menor1 n =

head [x | x <- [1..], numeroDivisores x == 2^n]

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

-- numeroDivisores 12 == 6

numeroDivisores :: Integer -> Integer

numeroDivisores =

genericLength . divisores

-- (divisores x) es la lista de los divisores de x. Por ejemplo,

-- divisores 12 == [1,3,2,6,4,12]

-- divisores 25 == [1,5,25]

divisores :: Integer -> [Integer]

divisores n =

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

-- 2ª solución

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

menor2 :: Integer -> Integer

menor2 n =

head [x | x <- [1..], numeroDivisores2 x == 2^n]

numeroDivisores2 :: Integer -> Integer

numeroDivisores2 =

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

-- 3ª solución

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

menor3 :: Integer -> Integer

menor3 n = product (genericTake n potencias)

-- potencias es la sucesión de las potencias de la forma p^(2^k),

-- donde p es un número primo y k es un número natural, ordenadas de

-- menor a mayor. Por ejemplo,

-- take 14 potencias == [2,3,4,5,7,9,11,13,16,17,19,23,25,29]

potencias :: [Integer]

potencias = 2 : mezcla (tail primes) (map (^2) potencias)

-- (mezcla xs ys) es la lista obtenida mezclando las dos listas xs e ys,

-- que se suponen ordenadas y disjuntas. Por ejemplo,

-- λ> take 15 (mezcla [2^n | n <- [1..]] [3^n | n <- [1..]])

-- [2,3,4,8,9,16,27,32,64,81,128,243,256,512,729]

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

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

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

-- Comprobación de equivalencia

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

-- La propiedad es

prop_menor :: Positive Integer -> Bool

prop_menor (Positive n) =

all (== menor1 n)

[menor2 n,

menor3 n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=5}) prop_menor

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> menor 8

-- 1081080

-- (47.69 secs, 94,764,856,352 bytes)

-- λ> menor2 8

-- 1081080

-- (36.17 secs, 94,764,856,368 bytes)

-- λ> menor3 8

-- 1081080

-- (0.00 secs, 116,960 bytes)

-- Definición de menor

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

-- En lo que sigue, usaremos menor3 como menor.

menor :: Integer -> Integer

menor = menor3

-- Propiedad

-- =========

-- La propiedad es

prop_menor_divide :: Positive Integer -> Bool

prop_menor_divide (Positive n) =

menor (n+1) `mod` menor n == 0

-- La comprobación es

-- λ> quickCheck prop_menor_divide

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Distancia esperada entre dos puntos de un cuadrado unitario

PorJosé A. Alonso 10-junio-20227-junio-2022

Definir, por simulación, la función

   distanciaEsperada :: Int -> IO Double

1	distanciaEsperada :: Int -> IO Double

tal que (distanciaEsperada n) es la distancia esperada entre n puntos del cuadrado unitario de vértices opuestos (0,0) y (1,1), elegidos aleatoriamente. Por ejemplo,

   distanciaEsperada 10     ==  0.43903617921423593
   distanciaEsperada 10     ==  0.6342350621260004
   distanciaEsperada 100    ==  0.5180418995364429
   distanciaEsperada 100    ==  0.5288261085653962
   distanciaEsperada 1000   ==  0.5143804432569616
   distanciaEsperada 10000  ==  0.5208360147922616

distanciaEsperada 10 == 0.43903617921423593

distanciaEsperada 10 == 0.6342350621260004

distanciaEsperada 100 == 0.5180418995364429

distanciaEsperada 100 == 0.5288261085653962

distanciaEsperada 1000 == 0.5143804432569616

distanciaEsperada 10000 == 0.5208360147922616

El valor exacto de la distancia esperada es

   ve = (sqrt(2) + 2 + 5*log(1+sqrt(2)))/15 = 0.5214054331647207

1	ve = (sqrt(2) + 2 + 5*log(1+sqrt(2)))/15 = 0.5214054331647207

Definir la función

   graficaDistanciaEsperada :: [Int] -> IO ()

1	graficaDistanciaEsperada :: [Int] -> IO ()

tal que (graficaDistanciaEsperada ns) dibuja las gráficas de los pares (n, distanciaEsperada n) para n en la lista creciente ns junto con la recta y = ve, donde ve es el valor exacto. Por ejemplo, (graficaDistanciaEsperada [10,30..4000]) dibuja

Soluciones

import Data.List     (genericLength)
import System.Random (newStdGen, randomRIO, randomRs)
import Control.Monad (replicateM)
import Graphics.Gnuplot.Simple

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

-- Un punto es un par de números reales.
type Punto = (Double, Double)

-- (puntosDelCuadrado n) es una lista de n puntos del cuadrado
-- unitario de vértices opuestos (0,0) y (1,1). Por ejemplo, 
--    λ> puntosDelCuadrado 3
--    [(0.6067427807212623,0.24785843546479303),
--     (0.9579158098726746,8.047408846191773e-2),
--     (0.856758357789639,0.9814972717003113)]
--    λ> puntosDelCuadrado 3
--    [(1.9785720974027532e-2,0.6343219201012211),
--     (0.21903717179861604,0.20947986189590784),
--     (0.4739903340716357,1.2262474491489095e-2)]
puntosDelCuadrado :: Int -> IO [Punto]
puntosDelCuadrado n = do
  gen <- newStdGen
  let xs = randomRs (0,1) gen
      (as, ys) = splitAt n xs
      (bs, _)  = splitAt n ys
  return (zip as bs)
  
-- (distancia p1 p2) es la distancia entre los puntos p1 y p2. Por
-- ejemplo,
--    distancia (0,0) (3,4)  ==  5.0
distancia :: Punto -> Punto -> Double
distancia (x1,y1) (x2,y2) = sqrt ((x1-x2)^2+(y1-y2)^2)

-- (distancias ps) es la lista de las distancias entre los elementos 1º
-- y 2º, 3º y 4º, ... de ps. Por ejemplo,
--    distancias [(0,0),(3,4),(1,1),(7,9)]  ==  [5.0,10.0]
distancias :: [Punto] -> [Double]
distancias []         = []
distancias (p1:p2:ps) = distancia p1 p2 : distancias ps

-- (media xs) es la media aritmética de los elementos de xs. Por ejemplo,
--    media [1,7,1]  ==  3.0
media :: [Double] -> Double
media xs =
  sum xs / genericLength xs

-- (distanciaEsperada n) es la distancia esperada entre n puntos
-- aleatorios en el cuadrado unitario. Por ejemplo,
--    λ> distanciaEsperada 100
--    0.4712858421448363
--    λ> distanciaEsperada 100
--    0.4947745206856711
distanciaEsperada :: Int -> IO Double
distanciaEsperada n = do
  ps <- puntosDelCuadrado (2*n)
  return (media (distancias ps))

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

distanciaEsperada2 :: Int -> IO Double
distanciaEsperada2 n = do
  ps <- puntosDelCuadrado2 (2*n)
  return (media (distancias ps))

-- (puntosDelCuadrado2 n) es una lista de n puntos del cuadrado
-- unitario de vértices opuestos (0,0) y (1,1). Por ejemplo, 
--    λ> puntosDelCuadrado2 3
--    [(0.9836699352638695,0.5143414844876929),
--     (0.8715237339877027,0.9905157772823782),
--     (0.29502946161912935,0.16889248111565192)]
--    λ> puntosDelCuadrado2 3
--    [(0.20405570457106392,0.47574116941605116),
--     (0.7128182811364226,3.201419787777959e-2),
--     (0.5576891231675457,0.9994474730919443)]
puntosDelCuadrado2 :: Int -> IO [Punto]
puntosDelCuadrado2 n =
  replicateM n puntoDelCuadrado2

-- (puntoDelCuadrado2 n) es un punto del cuadrado unitario de vértices
-- opuestos (0,0) y (1,1). Por ejemplo,  
--    λ> puntoDelCuadrado2
--    (0.7512991739803923,0.966436016138578)
--    λ> puntoDelCuadrado2
--    (0.7306826194847795,0.8984574498515252)
puntoDelCuadrado2 :: IO Punto
puntoDelCuadrado2 = do
  x <- randomRIO (0, 1.0)
  y <- randomRIO (0, 1.0)
  return (x, y)

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

distanciaEsperada3 :: Int -> IO Double
distanciaEsperada3 n = do
  ds <- distanciasAleatorias n
  return (media ds)

-- (distanciasAleatorias n) es la lista de las distancias aleatorias
-- entre n pares de puntos del cuadrado unitario. Por ejemplo, 
--    λ> distanciasAleatorias 3
--    [0.8325589110989705,0.6803336613847881,0.1690051224111662]
--    λ> distanciasAleatorias 3
--    [0.3470124940889039,0.459002678562019,0.7665623634969365]
distanciasAleatorias :: Int -> IO [Double]
distanciasAleatorias n = 
  replicateM n distanciaAleatoria

-- distanciaAleatoria es la distancia de un par de punto del cuadrado
-- unitario elegidos aleatoriamente. Por ejemplo,
--    λ> distanciaAleatoria
--    0.8982361685460913
--    λ> distanciaAleatoria
--    0.9777207485571939
--    λ> distanciaAleatoria
--    0.6042223512347842
distanciaAleatoria :: IO Double
distanciaAleatoria = do 
  p1 <- puntoDelCuadrado2
  distancia p1 <$> puntoDelCuadrado2

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

distanciaEsperada4 :: Int -> IO Double
distanciaEsperada4 n =
  media <$> distanciasAleatorias n

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

graficaDistanciaEsperada :: [Int] -> IO ()
graficaDistanciaEsperada ns = do
  ys <- mapM distanciaEsperada ns
  let e = (sqrt 2 + 2 + 5 * log (1 + sqrt 2)) / 15
  plotLists [ Key Nothing
            -- , PNG "Distancia_esperada_entre_dos_puntos.png"
            ]
            [ zip ns ys
            , zip ns (repeat e)]

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 (genericLength)

import System.Random (newStdGen, randomRIO, randomRs)

import Control.Monad (replicateM)

import Graphics.Gnuplot.Simple

-- 1ª solución

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

-- Un punto es un par de números reales.

type Punto = (Double, Double)

-- (puntosDelCuadrado n) es una lista de n puntos del cuadrado

-- unitario de vértices opuestos (0,0) y (1,1). Por ejemplo,

-- λ> puntosDelCuadrado 3

-- [(0.6067427807212623,0.24785843546479303),

-- (0.9579158098726746,8.047408846191773e-2),

-- (0.856758357789639,0.9814972717003113)]

-- λ> puntosDelCuadrado 3

-- [(1.9785720974027532e-2,0.6343219201012211),

-- (0.21903717179861604,0.20947986189590784),

-- (0.4739903340716357,1.2262474491489095e-2)]

puntosDelCuadrado :: Int -> IO [Punto]

puntosDelCuadrado n = do

gen <- newStdGen

let xs = randomRs (0,1) gen

(as, ys) = splitAt n xs

(bs, _) = splitAt n ys

return (zip as bs)

-- (distancia p1 p2) es la distancia entre los puntos p1 y p2. Por

-- ejemplo,

-- distancia (0,0) (3,4) == 5.0

distancia :: Punto -> Punto -> Double

distancia (x1,y1) (x2,y2) = sqrt ((x1-x2)^2+(y1-y2)^2)

-- (distancias ps) es la lista de las distancias entre los elementos 1º

-- y 2º, 3º y 4º, ... de ps. Por ejemplo,

-- distancias [(0,0),(3,4),(1,1),(7,9)] == [5.0,10.0]

distancias :: [Punto] -> [Double]

distancias [] = []

distancias (p1:p2:ps) = distancia p1 p2 : distancias ps

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

-- media [1,7,1] == 3.0

media :: [Double] -> Double

media xs =

sum xs / genericLength xs

-- (distanciaEsperada n) es la distancia esperada entre n puntos

-- aleatorios en el cuadrado unitario. Por ejemplo,

-- λ> distanciaEsperada 100

-- 0.4712858421448363

-- λ> distanciaEsperada 100

-- 0.4947745206856711

distanciaEsperada :: Int -> IO Double

distanciaEsperada n = do

ps <- puntosDelCuadrado (2*n)

return (media (distancias ps))

-- 2ª solución

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

distanciaEsperada2 :: Int -> IO Double

distanciaEsperada2 n = do

ps <- puntosDelCuadrado2 (2*n)

return (media (distancias ps))

-- (puntosDelCuadrado2 n) es una lista de n puntos del cuadrado

-- unitario de vértices opuestos (0,0) y (1,1). Por ejemplo,

-- λ> puntosDelCuadrado2 3

-- [(0.9836699352638695,0.5143414844876929),

-- (0.8715237339877027,0.9905157772823782),

-- (0.29502946161912935,0.16889248111565192)]

-- λ> puntosDelCuadrado2 3

-- [(0.20405570457106392,0.47574116941605116),

-- (0.7128182811364226,3.201419787777959e-2),

-- (0.5576891231675457,0.9994474730919443)]

puntosDelCuadrado2 :: Int -> IO [Punto]

puntosDelCuadrado2 n =

replicateM n puntoDelCuadrado2

-- (puntoDelCuadrado2 n) es un punto del cuadrado unitario de vértices

-- opuestos (0,0) y (1,1). Por ejemplo,

-- λ> puntoDelCuadrado2

-- (0.7512991739803923,0.966436016138578)

-- λ> puntoDelCuadrado2

-- (0.7306826194847795,0.8984574498515252)

puntoDelCuadrado2 :: IO Punto

puntoDelCuadrado2 = do

x <- randomRIO (0, 1.0)

y <- randomRIO (0, 1.0)

return (x, y)

-- 3ª solución

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

distanciaEsperada3 :: Int -> IO Double

distanciaEsperada3 n = do

ds <- distanciasAleatorias n

return (media ds)

-- (distanciasAleatorias n) es la lista de las distancias aleatorias

-- entre n pares de puntos del cuadrado unitario. Por ejemplo,

-- λ> distanciasAleatorias 3

-- [0.8325589110989705,0.6803336613847881,0.1690051224111662]

-- λ> distanciasAleatorias 3

-- [0.3470124940889039,0.459002678562019,0.7665623634969365]

distanciasAleatorias :: Int -> IO [Double]

distanciasAleatorias n =

replicateM n distanciaAleatoria

-- distanciaAleatoria es la distancia de un par de punto del cuadrado

-- unitario elegidos aleatoriamente. Por ejemplo,

-- λ> distanciaAleatoria

-- 0.8982361685460913

-- λ> distanciaAleatoria

-- 0.9777207485571939

-- λ> distanciaAleatoria

-- 0.6042223512347842

distanciaAleatoria :: IO Double

distanciaAleatoria = do

p1 <- puntoDelCuadrado2

distancia p1 <$> puntoDelCuadrado2

-- 4ª solución

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

distanciaEsperada4 :: Int -> IO Double

distanciaEsperada4 n =

media <$> distanciasAleatorias n

-- Gráfica

-- =======

graficaDistanciaEsperada :: [Int] -> IO ()

graficaDistanciaEsperada ns = do

ys <- mapM distanciaEsperada ns

let e = (sqrt 2 + 2 + 5 * log (1 + sqrt 2)) / 15

plotLists [ Key Nothing

-- , PNG "Distancia_esperada_entre_dos_puntos.png"

]

[ zip ns ys

, zip ns (repeat e)]

El código se encuentra en GitHub.

Representación matricial de relaciones binarias

PorJosé A. Alonso 9-junio-20226-junio-2022

Dada una relación r sobre un conjunto de números enteros, la matriz asociada a r es una matriz booleana p (cuyos elementos son True o False), tal que p(i,j) = True si y sólo si i está relacionado con j mediante la relación r.

Las relaciones binarias homogéneas y las matrices booleanas se pueden representar por

   type Relacion = ([Int],[(Int,Int)])
   type Matriz = Array (Int,Int) Bool

1 2	type Relacion = ([Int],[(Int,Int)]) type Matriz = Array (Int,Int) Bool

Definir la función

   matrizRB:: Relacion -> Matriz

1	matrizRB:: Relacion -> Matriz

tal que (matrizRB r) es la matriz booleana asociada a r. Por ejemplo,

   λ> matrizRB ([1..3],[(1,1), (1,3), (3,1), (3,3)])
   array ((1,1),(3,3)) [((1,1),True) ,((1,2),False),((1,3),True),
                        ((2,1),False),((2,2),False),((2,3),False),
                        ((3,1),True) ,((3,2),False),((3,3),True)]
   λ> matrizRB ([1..3],[(1,3), (3,1)])
   array ((1,1),(3,3)) [((1,1),False),((1,2),False),((1,3),True),
                        ((2,1),False),((2,2),False),((2,3),False),
                        ((3,1),True) ,((3,2),False),((3,3),False)]
   λ> let n = 10^4 in matrizRB3 ([1..n],[(1,n),(n,1)]) ! (n,n)
   False

λ> matrizRB ([1..3],[(1,1), (1,3), (3,1), (3,3)])

array ((1,1),(3,3)) [((1,1),True) ,((1,2),False),((1,3),True),

((2,1),False),((2,2),False),((2,3),False),

((3,1),True) ,((3,2),False),((3,3),True)]

λ> matrizRB ([1..3],[(1,3), (3,1)])

array ((1,1),(3,3)) [((1,1),False),((1,2),False),((1,3),True),

((2,1),False),((2,2),False),((2,3),False),

((3,1),True) ,((3,2),False),((3,3),False)]

λ> let n = 10^4 in matrizRB3 ([1..n],[(1,n),(n,1)]) ! (n,n)

False

Soluciones

import Data.Array      (Array, accumArray, array, listArray)
import Test.QuickCheck (Arbitrary, Gen, arbitrary, sublistOf, suchThat,
                        quickCheck)

type Relacion = ([Int],[(Int,Int)])
type Matriz   = Array (Int,Int) Bool

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

matrizRB1 :: Relacion -> Matriz 
matrizRB1 r = 
    array ((1,1),(n,n)) 
          [((a,b), (a,b) `elem` grafo r) | a <- [1..n], b <- [1..n]]
    where n = maximum (universo r)
          universo (us,_) = us
          grafo (_,ps)    = ps

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

matrizRB2 :: Relacion -> Matriz
matrizRB2 r = 
    listArray ((1,1),(n,n)) 
              [(a,b) `elem` snd r | a <- [1..n], b <- [1..n]]
    where n = maximum (fst r)

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

matrizRB3 :: Relacion -> Matriz
matrizRB3 r = 
    accumArray (||) False ((1,1),(n,n)) (zip (snd r) (repeat True))
    where n = maximum (fst r)

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

-- Tipo de relaciones binarias
newtype RB = RB Relacion
  deriving Show

-- relacionArbitraria genera una relación arbitraria. Por ejemplo, 
--    λ> generate relacionArbitraria
--    RB ([1,2,3,4,5],[(1,4),(1,5),(2,3),(2,4),(4,2),(4,3),(4,4),(5,1),(5,2),(5,3),(5,4)])
relacionArbitraria :: Gen RB
relacionArbitraria = do
  n <- arbitrary `suchThat` (> 1)
  xs <- sublistOf [(x,y) | x <- [1..n], y <- [1..n]]
  return (RB ([1..n], xs))

-- RB es una subclase de Arbitrary
instance Arbitrary RB where
  arbitrary = relacionArbitraria

-- La propiedad es
prop_matrizRB :: RB -> Bool
prop_matrizRB (RB r) =
  all (== matrizRB1 r)
      [matrizRB2 r,
       matrizRB3 r]

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

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

-- La comparación es
--    λ> let n = 2000 in matrizRB1 ([1..n],[(1,n),(n,1)]) ! (n,n)
--    False
--    (2.02 secs, 1,505,248,912 bytes)
--    λ> let n = 2000 in matrizRB2 ([1..n],[(1,n),(n,1)]) ! (n,n)
--    False
--    (1.92 secs, 833,232,360 bytes)
--    λ> let n = 2000 in matrizRB3 ([1..n],[(1,n),(n,1)]) ! (n,n)
--    False
--    (0.05 secs, 32,848,696 bytes)

import Data.Array (Array, accumArray, array, listArray)

import Test.QuickCheck (Arbitrary, Gen, arbitrary, sublistOf, suchThat,

quickCheck)

type Relacion = ([Int],[(Int,Int)])

type Matriz = Array (Int,Int) Bool

-- 1ª solución

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

matrizRB1 :: Relacion -> Matriz

matrizRB1 r =

array ((1,1),(n,n))

[((a,b), (a,b) `elem` grafo r) | a <- [1..n], b <- [1..n]]

where n = maximum (universo r)

universo (us,_) = us

grafo (_,ps) = ps

-- 2ª solución

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

matrizRB2 :: Relacion -> Matriz

matrizRB2 r =

listArray ((1,1),(n,n))

[(a,b) `elem` snd r | a <- [1..n], b <- [1..n]]

where n = maximum (fst r)

-- 3ª solución

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

matrizRB3 :: Relacion -> Matriz

matrizRB3 r =

accumArray (||) False ((1,1),(n,n)) (zip (snd r) (repeat True))

where n = maximum (fst r)

-- Comprobación de equivalencia

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

-- Tipo de relaciones binarias

newtype RB = RB Relacion

deriving Show

-- relacionArbitraria genera una relación arbitraria. Por ejemplo,

-- λ> generate relacionArbitraria

-- RB ([1,2,3,4,5],[(1,4),(1,5),(2,3),(2,4),(4,2),(4,3),(4,4),(5,1),(5,2),(5,3),(5,4)])

relacionArbitraria :: Gen RB

relacionArbitraria = do

n <- arbitrary `suchThat` (> 1)

xs <- sublistOf [(x,y) | x <- [1..n], y <- [1..n]]

return (RB ([1..n], xs))

-- RB es una subclase de Arbitrary

instance Arbitrary RB where

arbitrary = relacionArbitraria

-- La propiedad es

prop_matrizRB :: RB -> Bool

prop_matrizRB (RB r) =

all (== matrizRB1 r)

[matrizRB2 r,

matrizRB3 r]

-- La comprobación es

-- λ> quickCheck prop_matrixzB

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> let n = 2000 in matrizRB1 ([1..n],[(1,n),(n,1)]) ! (n,n)

-- False

-- (2.02 secs, 1,505,248,912 bytes)

-- λ> let n = 2000 in matrizRB2 ([1..n],[(1,n),(n,1)]) ! (n,n)

-- False

-- (1.92 secs, 833,232,360 bytes)

-- λ> let n = 2000 in matrizRB3 ([1..n],[(1,n),(n,1)]) ! (n,n)

-- False

-- (0.05 secs, 32,848,696 bytes)

El código se encuentra en GitHub.

Codificación de Gödel

PorJosé A. Alonso 8-junio-20225-junio-2022

Dada una lista de números naturales xs, codificación de Gödel de xs se obtiene multiplicando las potencias de los primos sucesivos, siendo los exponentes los sucesores de los elementos de xs. Por ejemplo, si xs = [6,0,4], la codificación de xs es

   2^7 * 3^1 * 5^5 = 1200000

1	2^7 * 3^1 * 5^5 = 1200000

Definir las funciones

   codificaG   :: [Integer] -> Integer
   decodificaG :: Integer -> [Integer]

1 2	codificaG :: [Integer] -> Integer decodificaG :: Integer -> [Integer]

tales que

(codificaG xs) es la codificación de Gödel de xs. Por ejemplo,

     codificaG [6,0,4]            ==  1200000
     codificaG [3,1,1]            ==  3600
     codificaG [3,1,0,0,0,0,0,1]  ==  4423058640
     codificaG [1..6]             ==  126111168580452537982500

codificaG [6,0,4] == 1200000

codificaG [3,1,1] == 3600

codificaG [3,1,0,0,0,0,0,1] == 4423058640

codificaG [1..6] == 126111168580452537982500

(decodificaG n) es la lista xs cuya codificación es n. Por ejemplo,

     decodificaG 1200000                   ==  [6,0,4]
     decodificaG 3600                      ==  [3,1,1]
     decodificaG 4423058640                ==  [3,1,0,0,0,0,0,1]
     decodificaG 126111168580452537982500  ==  [1,2,3,4,5,6]

decodificaG 1200000 == [6,0,4]

decodificaG 3600 == [3,1,1]

decodificaG 4423058640 == [3,1,0,0,0,0,0,1]

decodificaG 126111168580452537982500 == [1,2,3,4,5,6]

Comprobar con QuickCheck que ambas funciones son inversas; es decir,

   decodificaG (codificaG xs) = xs
   codificaG (decodificaG n) = n

1 2	decodificaG (codificaG xs) = xs codificaG (decodificaG n) = n

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (primes, primeFactors)
import Test.QuickCheck (NonNegative (getNonNegative),
                        NonEmptyList (NonEmpty),
                        Positive (Positive, getPositive), quickCheck)

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

codificaG1 :: [Integer] -> Integer
codificaG1 = codificaG1' . codificaAux

codificaAux :: [Integer] -> [Integer]
codificaAux = map succ

codificaG1' :: [Integer] -> Integer
codificaG1' = aux primes
  where aux _ []          = 1
        aux (p:ps) (x:xs) = p^x * aux ps xs

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

codificaG2 :: [Integer] -> Integer
codificaG2 = codificaG2' . codificaAux

codificaG2' :: [Integer] -> Integer
codificaG2' xs = product [p^x | (p, x) <- zip primes xs]

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

codificaG3 :: [Integer] -> Integer
codificaG3 = codificaG3' . codificaAux

codificaG3' :: [Integer] -> Integer
codificaG3' xs = product (zipWith (^) primes xs)

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

codificaG4 :: [Integer] -> Integer
codificaG4 = codificaG4' . codificaAux

codificaG4' :: [Integer] -> Integer
codificaG4' = product . zipWith (^) primes

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

-- La propiedad es
prop_codificaG :: [NonNegative Integer] -> Bool
prop_codificaG xs =
  all (== codificaG1 ys)
      [codificaG2 ys,
       codificaG3 ys,
       codificaG4 ys]
  where ys = map getNonNegative xs

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

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

-- La comparación es
--    λ> length (show (codificaG1 (replicate (4*10^4) 1)))
--    208100
--    (1.73 secs, 2,312,100,136 bytes)
--    λ> length (show (codificaG2 (replicate (4*10^4) 1)))
--    208100
--    (1.63 secs, 2,327,676,832 bytes)
--    λ> length (show (codificaG3 (replicate (4*10^4) 1)))
--    208100
--    (1.62 secs, 2,323,836,832 bytes)
--    λ> length (show (codificaG4 (replicate (4*10^4) 1)))
--    208100
--    (1.54 secs, 2,147,635,680 bytes)

-- Definición de codificaG
-- =======================

-- Usaremos la 4ª
codificaG :: [Integer] -> Integer
codificaG = codificaG4

-- Definición de decodificaG
-- =========================

decodificaG :: Integer -> [Integer]
decodificaG = decodificaAux . decodificaG'

decodificaAux :: [Integer] -> [Integer]
decodificaAux = map pred

decodificaG' :: Integer -> [Integer]
decodificaG' 1 = [0]
decodificaG' n = aux primes (group (primeFactors n))
  where aux _ [] = []
        aux (x:xs) ((y:ys):yss) | x == y    = 1 + genericLength ys : aux xs yss
                                | otherwise = 0 : aux xs ((y:ys):yss)

-- Comprobación de propiedades
-- ===========================

-- La primera propiedad es
prop_decodificaG_codificaG :: NonEmptyList (NonNegative Integer) -> Bool
prop_decodificaG_codificaG (NonEmpty xs) = 
  decodificaG (codificaG ys) == ys
  where ys = map getNonNegative xs

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

-- La 2ª propiedad es
prop_codificaG_decodificaG :: Positive Integer -> Bool
prop_codificaG_decodificaG (Positive n) = 
  codificaG (decodificaG n) == n

-- la comprobación es
--    λ> quickCheck prop_codificaG_decodificaG
--    +++ 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

import Data.List (genericLength, group)

import Data.Numbers.Primes (primes, primeFactors)

import Test.QuickCheck (NonNegative (getNonNegative),

NonEmptyList (NonEmpty),

Positive (Positive, getPositive), quickCheck)

-- 1ª definición de codificaG

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

codificaG1 :: [Integer] -> Integer

codificaG1 = codificaG1' . codificaAux

codificaAux :: [Integer] -> [Integer]

codificaAux = map succ

codificaG1' :: [Integer] -> Integer

codificaG1' = aux primes

where aux _ [] = 1

aux (p:ps) (x:xs) = p^x * aux ps xs

-- 2ª definición de codificaG

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

codificaG2 :: [Integer] -> Integer

codificaG2 = codificaG2' . codificaAux

codificaG2' :: [Integer] -> Integer

codificaG2' xs = product [p^x | (p, x) <- zip primes xs]

-- 3ª definición de codificaG

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

codificaG3 :: [Integer] -> Integer

codificaG3 = codificaG3' . codificaAux

codificaG3' :: [Integer] -> Integer

codificaG3' xs = product (zipWith (^) primes xs)

-- 4ª definición de codificaG

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

codificaG4 :: [Integer] -> Integer

codificaG4 = codificaG4' . codificaAux

codificaG4' :: [Integer] -> Integer

codificaG4' = product . zipWith (^) primes

-- Comprobación de equivalencia de codificaG

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

-- La propiedad es

prop_codificaG :: [NonNegative Integer] -> Bool

prop_codificaG xs =

all (== codificaG1 ys)

[codificaG2 ys,

codificaG3 ys,

codificaG4 ys]

where ys = map getNonNegative xs

-- La comprobación es

-- λ> quickCheck prop_codificaG

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de codificaG

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

-- La comparación es

-- λ> length (show (codificaG1 (replicate (4*10^4) 1)))

-- 208100

-- (1.73 secs, 2,312,100,136 bytes)

-- λ> length (show (codificaG2 (replicate (4*10^4) 1)))

-- 208100

-- (1.63 secs, 2,327,676,832 bytes)

-- λ> length (show (codificaG3 (replicate (4*10^4) 1)))

-- 208100

-- (1.62 secs, 2,323,836,832 bytes)

-- λ> length (show (codificaG4 (replicate (4*10^4) 1)))

-- 208100

-- (1.54 secs, 2,147,635,680 bytes)

-- Definición de codificaG

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

-- Usaremos la 4ª

codificaG :: [Integer] -> Integer

codificaG = codificaG4

-- Definición de decodificaG

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

decodificaG :: Integer -> [Integer]

decodificaG = decodificaAux . decodificaG'

decodificaAux :: [Integer] -> [Integer]

decodificaAux = map pred

decodificaG' :: Integer -> [Integer]

decodificaG' 1 = [0]

decodificaG' n = aux primes (group (primeFactors n))

where aux _ [] = []

aux (x:xs) ((y:ys):yss) | x == y = 1 + genericLength ys : aux xs yss

| otherwise = 0 : aux xs ((y:ys):yss)

-- Comprobación de propiedades

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

-- La primera propiedad es

prop_decodificaG_codificaG :: NonEmptyList (NonNegative Integer) -> Bool

prop_decodificaG_codificaG (NonEmpty xs) =

decodificaG (codificaG ys) == ys

where ys = map getNonNegative xs

-- La comprobación es

-- λ> quickCheck prop_decodificaG_codificaG

-- +++ OK, passed 100 tests.

-- La 2ª propiedad es

prop_codificaG_decodificaG :: Positive Integer -> Bool

prop_codificaG_decodificaG (Positive n) =

codificaG (decodificaG n) == n

-- la comprobación es

-- λ> quickCheck prop_codificaG_decodificaG

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Primos circulares

PorJosé A. Alonso 7-junio-20223-junio-2022

Un primo circular es un número tal que todas las rotaciones de dígitos producen números primos. Por ejemplo, 195 es un primo circular ya que las rotaciones de sus dígitos son 197, 971 y 719 y los tres números son primos.

Definir la lista

   circulares :: [Integer]

1	circulares :: [Integer]

cuyo valor es la lista de los números primos circulares. Por ejemplo,

   take 16 circulares == [2,3,5,7,11,13,17,31,37,71,73,79,97,113,131,197]
   circulares !! 50   == 933199

1 2	take 16 circulares == [2,3,5,7,11,13,17,31,37,71,73,79,97,113,131,197] circulares !! 50 == 933199

Soluciones

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

-- 1ª solución
-- ===========
  
circulares1 :: [Integer]
circulares1 = filter esCircular1 primes

-- (esCircular1 n) se verifica si n es un número circular. Por ejemplo, 
--    esCircular1 197  ==  True
--    esCircular1 157  ==  False
esCircular1 :: Integer -> Bool
esCircular1 = all (isPrime . read) . rotaciones1 . show 

-- (rotaciones1 xs) es la lista de las rotaciones obtenidas desplazando
-- el primer elemento xs al final. Por ejemplo,
--    rotaciones1 [2,3,5]  ==  [[2,3,5],[3,5,2],[5,2,3]]
rotaciones1 :: [a] -> [[a]]
rotaciones1 xs = reverse (aux (length xs) [xs])
    where aux 1 yss      = yss
          aux n (ys:yss) = aux (n-1) (rota ys : ys :yss)

-- 2ª solución
-- ===========
  
circulares2 :: [Integer]
circulares2 = filter esCircular2 primes

esCircular2 :: Integer -> Bool
esCircular2 = all (isPrime . read) . rotaciones2 . show 

rotaciones2 :: [a] -> [[a]]
rotaciones2 xs = take (length xs) (iterate rota xs)

-- (rota xs) es la lista añadiendo el primer elemento de xs al
-- final. Por ejemplo, 
--    rota [3,2,5,7]  ==  [2,5,7,3]
rota :: [a] -> [a]
rota (x:xs) = xs ++ [x]

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

circulares3 :: [Integer]
circulares3 = filter (all isPrime . rotaciones3) primes 

rotaciones3 :: Integer -> [Integer]
rotaciones3 n = [read (take m (drop i (cycle s))) | i <- [1..m]]
    where s = show n
          m = length s

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

-- Nota. La 4ª definición es una mejora observando que para que n sea un
-- número primo circular es necesario que todos los dígitos de n sean
-- impares, salvo para n = 2. 

circulares4 :: [Integer]
circulares4 = 2 : filter esCircular4 primes

-- (esCircular4 n) se verifica si n es un número circular. Por ejemplo, 
--    esCircular4 197  ==  True
--    esCircular4 157  ==  False
esCircular4 :: Integer -> Bool
esCircular4 n = digitosImpares n && 
                all (isPrime . read) (rotaciones2 (show n))

-- (digitosImpares n) se verifica si todos los dígitos de n son
-- impares. Por ejemplo,
--    digitosImpares 7351  ==  True
--    digitosImpares 7341  ==  False
digitosImpares :: Integer -> Bool
digitosImpares = all (`elem` "135679") . show

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

-- La propiedad es
prop_circulares :: Positive Int -> Bool
prop_circulares (Positive n) =
  all (== circulares1 !! n)
      [circulares2 !! n,
       circulares3 !! n,
       circulares4 !! n]

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

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

-- La comparación es
--    λ> circulares1 !! 46
--    331999
--    (2.08 secs, 7,229,208,200 bytes)
--    λ> circulares2 !! 46
--    331999
--    (1.93 secs, 7,165,043,992 bytes)
--    λ> circulares3 !! 46
--    331999
--    (0.74 secs, 2,469,098,648 bytes)
--    λ> circulares4 !! 46
--    331999
--    (0.28 secs, 917,501,600 bytes)

100

101

102

103

104

105

106

107

import Data.Numbers.Primes (isPrime, primes)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

circulares1 :: [Integer]

circulares1 = filter esCircular1 primes

-- (esCircular1 n) se verifica si n es un número circular. Por ejemplo,

-- esCircular1 197 == True

-- esCircular1 157 == False

esCircular1 :: Integer -> Bool

esCircular1 = all (isPrime . read) . rotaciones1 . show

-- (rotaciones1 xs) es la lista de las rotaciones obtenidas desplazando

-- el primer elemento xs al final. Por ejemplo,

-- rotaciones1 [2,3,5] == [[2,3,5],[3,5,2],[5,2,3]]

rotaciones1 :: [a] -> [[a]]

rotaciones1 xs = reverse (aux (length xs) [xs])

where aux 1 yss = yss

aux n (ys:yss) = aux (n-1) (rota ys : ys :yss)

-- 2ª solución

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

circulares2 :: [Integer]

circulares2 = filter esCircular2 primes

esCircular2 :: Integer -> Bool

esCircular2 = all (isPrime . read) . rotaciones2 . show

rotaciones2 :: [a] -> [[a]]

rotaciones2 xs = take (length xs) (iterate rota xs)

-- (rota xs) es la lista añadiendo el primer elemento de xs al

-- final. Por ejemplo,

-- rota [3,2,5,7] == [2,5,7,3]

rota :: [a] -> [a]

rota (x:xs) = xs ++ [x]

-- 3ª solución

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

circulares3 :: [Integer]

circulares3 = filter (all isPrime . rotaciones3) primes

rotaciones3 :: Integer -> [Integer]

rotaciones3 n = [read (take m (drop i (cycle s))) | i <- [1..m]]

where s = show n

m = length s

-- 4ª definición

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

-- Nota. La 4ª definición es una mejora observando que para que n sea un

-- número primo circular es necesario que todos los dígitos de n sean

-- impares, salvo para n = 2.

circulares4 :: [Integer]

circulares4 = 2 : filter esCircular4 primes

-- (esCircular4 n) se verifica si n es un número circular. Por ejemplo,

-- esCircular4 197 == True

-- esCircular4 157 == False

esCircular4 :: Integer -> Bool

esCircular4 n = digitosImpares n &&

all (isPrime . read) (rotaciones2 (show n))

-- (digitosImpares n) se verifica si todos los dígitos de n son

-- impares. Por ejemplo,

-- digitosImpares 7351 == True

-- digitosImpares 7341 == False

digitosImpares :: Integer -> Bool

digitosImpares = all (`elem` "135679") . show

-- Comprobación de equivalencia

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

-- La propiedad es

prop_circulares :: Positive Int -> Bool

prop_circulares (Positive n) =

all (== circulares1 !! n)

[circulares2 !! n,

circulares3 !! n,

circulares4 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=50}) prop_circulares

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> circulares1 !! 46

-- 331999

-- (2.08 secs, 7,229,208,200 bytes)

-- λ> circulares2 !! 46

-- 331999

-- (1.93 secs, 7,165,043,992 bytes)

-- λ> circulares3 !! 46

-- 331999

-- (0.74 secs, 2,469,098,648 bytes)

-- λ> circulares4 !! 46

-- 331999

-- (0.28 secs, 917,501,600 bytes)

El código se encuentra en GitHub.

Diccionario de frecuencias

PorJosé A. Alonso 6-junio-20222-junio-2022

Definir la función

   frecuencias :: Ord a => [a] -> Map a Int

1	frecuencias :: Ord a => [a] -> Map a Int

tal que (frecuencias xs) es el diccionario formado por los elementos de xs junto con el número de veces que aparecen en xs. Por ejemplo,

   λ> frecuencias "sosos"
   fromList [('o',2),('s',3)]
   λ> frecuencias (show (10^100))
   fromList [('0',100),('1',1)]
   λ> frecuencias (take (10^6) (cycle "abc"))
   fromList [('a',333334),('b',333333),('c',333333)]
   λ> size (frecuencias (take (10^6) (cycle [1..10^6])))
   1000000

λ> frecuencias "sosos"

fromList [('o',2),('s',3)]

λ> frecuencias (show (10^100))

fromList [('0',100),('1',1)]

λ> frecuencias (take (10^6) (cycle "abc"))

fromList [('a',333334),('b',333333),('c',333333)]

λ> size (frecuencias (take (10^6) (cycle [1..10^6])))

1000000

Soluciones

import Data.List (foldl')
import Data.Map (Map, empty, insertWith, fromList, fromListWith, size)
import Test.QuickCheck

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

frecuencias1 :: Ord a => [a] -> Map a Int
frecuencias1 []     = empty
frecuencias1 (x:xs) = insertWith (+) x 1 (frecuencias1 xs)

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

frecuencias2 :: Ord a => [a] -> Map a Int
frecuencias2 = foldl' (\d x-> insertWith (+) x 1 d) empty

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

frecuencias3 :: Ord a => [a] -> Map a Int
frecuencias3 xs = fromListWith (+) (zip xs (repeat 1))

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

-- La propiedad es
prop_frecuencias :: [Int] -> Bool
prop_frecuencias xs =
  all (== frecuencias1 xs)
      [ frecuencias2 xs
      , frecuencias3 xs]

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

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

-- La comparación es
--    λ> frecuencias1 (take (10^6) (cycle "abc"))
--    fromList [('a',333334),('b',333333),('c',333333)]
--    (0.89 secs, 453,842,448 bytes)
--    λ> frecuencias2 (take (10^6) (cycle "abc"))
--    fromList [('a',333334),('b',333333),('c',333333)]
--    (0.54 secs, 274,181,128 bytes)
--    λ> frecuencias3 (take (10^6) (cycle "abc"))
--    fromList [('a',333334),('b',333333),('c',333333)]
--    (0.29 secs, 313,787,976 bytes)
     
--    λ> size (frecuencias1 (take (10^6) (cycle [1..10^6])))
--    1000000
--    (3.76 secs, 2,651,926,024 bytes)
--    λ> size (frecuencias2 (take (10^6) (cycle [1..10^6])))
--    1000000
--    (1.03 secs, 1,640,678,448 bytes)
--    λ> size (frecuencias3 (take (10^6) (cycle [1..10^6])))
--    1000000
--    (0.88 secs, 1,672,678,536 bytes)

import Data.List (foldl')

import Data.Map (Map, empty, insertWith, fromList, fromListWith, size)

import Test.QuickCheck

-- 1ª solución

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

frecuencias1 :: Ord a => [a] -> Map a Int

frecuencias1 [] = empty

frecuencias1 (x:xs) = insertWith (+) x 1 (frecuencias1 xs)

-- 2ª solución

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

frecuencias2 :: Ord a => [a] -> Map a Int

frecuencias2 = foldl' (\d x-> insertWith (+) x 1 d) empty

-- 3ª solución

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

frecuencias3 :: Ord a => [a] -> Map a Int

frecuencias3 xs = fromListWith (+) (zip xs (repeat 1))

-- Equivalencia de las definiciones

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

-- La propiedad es

prop_frecuencias :: [Int] -> Bool

prop_frecuencias xs =

all (== frecuencias1 xs)

[ frecuencias2 xs

, frecuencias3 xs]

-- La comprobación es

-- λ> quickCheck prop_frecuencias

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> frecuencias1 (take (10^6) (cycle "abc"))

-- fromList [('a',333334),('b',333333),('c',333333)]

-- (0.89 secs, 453,842,448 bytes)

-- λ> frecuencias2 (take (10^6) (cycle "abc"))

-- fromList [('a',333334),('b',333333),('c',333333)]

-- (0.54 secs, 274,181,128 bytes)

-- λ> frecuencias3 (take (10^6) (cycle "abc"))

-- fromList [('a',333334),('b',333333),('c',333333)]

-- (0.29 secs, 313,787,976 bytes)

-- λ> size (frecuencias1 (take (10^6) (cycle [1..10^6])))

-- 1000000

-- (3.76 secs, 2,651,926,024 bytes)

-- λ> size (frecuencias2 (take (10^6) (cycle [1..10^6])))

-- 1000000

-- (1.03 secs, 1,640,678,448 bytes)

-- λ> size (frecuencias3 (take (10^6) (cycle [1..10^6])))

-- 1000000

-- (0.88 secs, 1,672,678,536 bytes)

El código se encuentra en GitHub.

Descomposiciones con sumandos 1 ó 2

PorJosé A. Alonso 3-junio-20225-junio-2022

Definir la funciones

   sumas  :: Int -> [[Int]]
   nSumas :: Int -> Integer

1 2	sumas :: Int -> [[Int]] nSumas :: Int -> Integer

tales que

(sumas n) es la lista de las descomposiciones de n como sumas cuyos sumandos son 1 ó 2. Por ejemplo,

      sumas 1            ==  [[1]]
      sumas 2            ==  [[1,1],[2]]
      sumas 3            ==  [[1,1,1],[1,2],[2,1]]
      sumas 4            ==  [[1,1,1,1],[1,1,2],[1,2,1],[2,1,1],[2,2]]
      length (sumas 26)  ==  196418
      length (sumas 33)  ==  5702887

sumas 1 == [[1]]

sumas 2 == [[1,1],[2]]

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

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

length (sumas 26) == 196418

length (sumas 33) == 5702887

(nSumas n) es el número de descomposiciones de n como sumas cuyos sumandos son 1 ó 2. Por ejemplo,

      nSumas 4                      ==  5
      nSumas 123                    ==  36726740705505779255899443
      length (show (nSumas 123456)) ==  25801

nSumas 4 == 5

nSumas 123 == 36726740705505779255899443

length (show (nSumas 123456)) == 25801

Soluciones

import Data.List  (genericIndex, genericLength)
import Data.Array ((!), array)
import Test.QuickCheck (Positive(Positive), quickCheckWith)

-- 1ª solución de sumas
-- ====================

sumas1 :: Int -> [[Int]]
sumas1 0 = [[]]
sumas1 1 = [[1]]
sumas1 n = [1:xs | xs <- sumas1 (n-1)] ++ [2:xs | xs <- sumas1 (n-2)]

-- 2ª solución de sumas
-- ====================

sumas2 :: Int -> [[Int]]
sumas2 0 = [[]]
sumas2 1 = [[1]]
sumas2 n = map (1:) (sumas2 (n-1)) ++ map (2:) (sumas2 (n-2))

-- 3ª solución de sumas
-- ====================

sumas3 :: Int -> [[Int]]
sumas3 n = v ! n
  where v = array (0,n) [(i, f i) | i <- [0..n]]
        f 0 = [[]]
        f 1 = [[1]]
        f k = map (1:) (v!(k-1)) ++ map (2:) (v!(k-2))
 
-- 4ª solución de sumas
-- ====================

sumas4 :: Int -> [[Int]]
sumas4 n = sucSumas !! n

-- sucSumas es la sucesión cuyo n-ésimo elemento es la lista de las
-- descomposiciones de n como sumas cuyos sumandos son 1 ó 2. Por
-- ejemplo,
--    λ> take 4 sucSumas
--    [[[]],[[1]],[[1,1],[2]],[[1,1,1],[1,2],[2,1]]]
--    λ> mapM_ print (take 5 sucSumas)
--    [[]]
--    [[1]]
--    [[1,1],[2]]
--    [[1,1,1],[1,2],[2,1]]
--    [[1,1,1,1],[1,1,2],[1,2,1],[2,1,1],[2,2]]
sucSumas :: [[[Int]]]
sucSumas = [[]] : [[1]] : zipWith f (tail sucSumas) sucSumas
  where f xs ys = map (1:) xs ++ map (2:) ys

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

-- La propiedad es
prop_sumas :: Positive Int -> Bool
prop_sumas (Positive n) =
  all (== sumas1 n)
      [sumas2 n,
       sumas3 n,
       sumas4 n]

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

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

-- La comparación es
--    λ> length (sumas1 28)
--    514229
--    (2.79 secs, 1,739,784,512 bytes)
--    λ> length (sumas2 28)
--    514229
--    (1.33 secs, 1,512,291,248 bytes)
--    λ> length (sumas3 28)
--    514229
--    (0.20 secs, 165,215,800 bytes)
--    λ> length (sumas4 28)
--    514229
--    (0.17 secs, 165,201,592 bytes)
--
--    λ> length (sumas3 33)
--    5702887
--    (2.16 secs, 1,830,761,864 bytes)
--    λ> length (sumas4 33)
--    5702887
--    (1.44 secs, 1,830,749,832 bytes)

-- Definición de sumas
-- ===================

-- La cuarta solución es más eficiente y es la que usaremos en lo
-- sucesivo:
sumas :: Int -> [[Int]]
sumas = sumas4

-- 1ª solución de nSumas
-- =====================

nSumas1 :: Int -> Integer
nSumas1 = genericLength . sumas2

-- 2ª solución de nSumas
-- =====================

nSumas2 :: Int -> Integer
nSumas2 0 = 1
nSumas2 1 = 1
nSumas2 n = nSumas2 (n-1) + nSumas2 (n-2)

-- 3ª solución de nSumas
-- =====================

nSumas3 :: Int -> Integer
nSumas3 n = v ! n
  where v = array (0,n) [(i,f i) | i <- [0..n]]
        f 0 = 1
        f 1 = 1
        f k = v ! (k-1) + v ! (k-2)

-- 4ª solución de nSumas
-- =====================

nSumas4 :: Int -> Integer
nSumas4 n = aux `genericIndex` n
  where aux = 1 : 1 : zipWith (+) aux (tail aux) 

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

-- La propiedad es
prop_nSumas :: Positive Int -> Bool
prop_nSumas (Positive n) =
  all (== nSumas1 n)
      [nSumas2 n,
       nSumas3 n,
       nSumas4 n]

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

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

-- La comparación es
--    λ> nSumas1 33
--    5702887
--    (17.32 secs, 23,140,562,600 bytes)
--    λ> nSumas2 33
--    5702887
--    (3.48 secs, 1,870,676,904 bytes)
--    λ> nSumas3 33
--    5702887
--    (0.00 secs, 152,960 bytes)
--    λ> nSumas4 33
--    5702887
--    (0.00 secs, 139,456 bytes)
--    
--    λ> length (show (nSumas3 (2*10^5)))
--    41798
--    (1.41 secs, 1,895,295,528 bytes)
--    λ> length (show (nSumas4 (2*10^5)))
--    41798
--    (2.39 secs, 1,834,998,800 bytes)

-- Nota. El valor de (nSumas n) es el n-ésimo término de la sucesión de
-- Fibonacci 1, 1, 2, 3, 5, 8, ...

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

169

170

import Data.List (genericIndex, genericLength)

import Data.Array ((!), array)

import Test.QuickCheck (Positive(Positive), quickCheckWith)

-- 1ª solución de sumas

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

sumas1 :: Int -> [[Int]]

sumas1 0 = [[]]

sumas1 1 = [[1]]

sumas1 n = [1:xs | xs <- sumas1 (n-1)] ++ [2:xs | xs <- sumas1 (n-2)]

-- 2ª solución de sumas

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

sumas2 :: Int -> [[Int]]

sumas2 0 = [[]]

sumas2 1 = [[1]]

sumas2 n = map (1:) (sumas2 (n-1)) ++ map (2:) (sumas2 (n-2))

-- 3ª solución de sumas

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

sumas3 :: Int -> [[Int]]

sumas3 n = v ! n

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

f 0 = [[]]

f 1 = [[1]]

f k = map (1:) (v!(k-1)) ++ map (2:) (v!(k-2))

-- 4ª solución de sumas

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

sumas4 :: Int -> [[Int]]

sumas4 n = sucSumas !! n

-- sucSumas es la sucesión cuyo n-ésimo elemento es la lista de las

-- descomposiciones de n como sumas cuyos sumandos son 1 ó 2. Por

-- ejemplo,

-- λ> take 4 sucSumas

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

-- λ> mapM_ print (take 5 sucSumas)

-- [[]]

-- [[1]]

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

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

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

sucSumas :: [[[Int]]]

sucSumas = [[]] : [[1]] : zipWith f (tail sucSumas) sucSumas

where f xs ys = map (1:) xs ++ map (2:) ys

-- Comprobación de equivalencia de sumas

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

-- La propiedad es

prop_sumas :: Positive Int -> Bool

prop_sumas (Positive n) =

all (== sumas1 n)

[sumas2 n,

sumas3 n,

sumas4 n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=20}) prop_sumas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de sumas

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

-- La comparación es

-- λ> length (sumas1 28)

-- 514229

-- (2.79 secs, 1,739,784,512 bytes)

-- λ> length (sumas2 28)

-- 514229

-- (1.33 secs, 1,512,291,248 bytes)

-- λ> length (sumas3 28)

-- 514229

-- (0.20 secs, 165,215,800 bytes)

-- λ> length (sumas4 28)

-- 514229

-- (0.17 secs, 165,201,592 bytes)

-- λ> length (sumas3 33)

-- 5702887

-- (2.16 secs, 1,830,761,864 bytes)

-- λ> length (sumas4 33)

-- 5702887

-- (1.44 secs, 1,830,749,832 bytes)

-- Definición de sumas

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

-- La cuarta solución es más eficiente y es la que usaremos en lo

-- sucesivo:

sumas :: Int -> [[Int]]

sumas = sumas4

-- 1ª solución de nSumas

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

nSumas1 :: Int -> Integer

nSumas1 = genericLength . sumas2

-- 2ª solución de nSumas

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

nSumas2 :: Int -> Integer

nSumas2 0 = 1

nSumas2 1 = 1

nSumas2 n = nSumas2 (n-1) + nSumas2 (n-2)

-- 3ª solución de nSumas

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

nSumas3 :: Int -> Integer

nSumas3 n = v ! n

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

f 0 = 1

f 1 = 1

f k = v ! (k-1) + v ! (k-2)

-- 4ª solución de nSumas

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

nSumas4 :: Int -> Integer

nSumas4 n = aux `genericIndex` n

where aux = 1 : 1 : zipWith (+) aux (tail aux)

-- Comprobación de equivalencia de nSumas

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

-- La propiedad es

prop_nSumas :: Positive Int -> Bool

prop_nSumas (Positive n) =

all (== nSumas1 n)

[nSumas2 n,

nSumas3 n,

nSumas4 n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=20}) prop_nSumas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de nSumas

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

-- La comparación es

-- λ> nSumas1 33

-- 5702887

-- (17.32 secs, 23,140,562,600 bytes)

-- λ> nSumas2 33

-- 5702887

-- (3.48 secs, 1,870,676,904 bytes)

-- λ> nSumas3 33

-- 5702887

-- (0.00 secs, 152,960 bytes)

-- λ> nSumas4 33

-- 5702887

-- (0.00 secs, 139,456 bytes)

-- λ> length (show (nSumas3 (2*10^5)))

-- 41798

-- (1.41 secs, 1,895,295,528 bytes)

-- λ> length (show (nSumas4 (2*10^5)))

-- 41798

-- (2.39 secs, 1,834,998,800 bytes)

-- Nota. El valor de (nSumas n) es el n-ésimo término de la sucesión de

-- Fibonacci 1, 1, 2, 3, 5, 8, ...

El código se encuentra en GitHub.

Suma de los elementos de las diagonales matrices espirales

PorJosé A. Alonso 2-junio-202231-mayo-2022

Empezando con el número 1 y moviéndose en el sentido de las agujas del reloj se obtienen las matrices espirales

   |1 2|   |7 8 9|   | 7  8  9 10|   |21 22 23 24 25|
   |4 3|   |6 1 2|   | 6  1  2 11|   |20  7  8  9 10|
           |5 4 3|   | 5  4  3 12|   |19  6  1  2 11|
                     |16 15 14 13|   |18  5  4  3 12|
                                     |17 16 15 14 13|

|1 2| |7 8 9| | 7 8 9 10| |21 22 23 24 25|

|4 3| |6 1 2| | 6 1 2 11| |20 7 8 9 10|

|5 4 3| | 5 4 3 12| |19 6 1 2 11|

|16 15 14 13| |18 5 4 3 12|

|17 16 15 14 13|

La suma los elementos de sus diagonales es

   + en la 2x2: 1+3+2+4               =  10
   + en la 3x3: 1+3+5+7+9             =  25
   + en la 4x4: 1+2+3+4+7+10+13+16    =  56
   + en la 5x5: 1+3+5+7+9+13+17+21+25 = 101

+ en la 2x2: 1+3+2+4 = 10

+ en la 3x3: 1+3+5+7+9 = 25

+ en la 4x4: 1+2+3+4+7+10+13+16 = 56

+ en la 5x5: 1+3+5+7+9+13+17+21+25 = 101

Definir la función

   sumaDiagonales :: Integer -> Integer

1	sumaDiagonales :: Integer -> Integer

tal que (sumaDiagonales n) es la suma de los elementos en las diagonales de la matriz espiral de orden nxn. Por ejemplo.

   sumaDiagonales 1         ==  1
   sumaDiagonales 2         ==  10
   sumaDiagonales 3         ==  25
   sumaDiagonales 4         ==  56
   sumaDiagonales 5         ==  101
   sumaDiagonales (10^6)    ==  666667166668000000
   sumaDiagonales (1+10^6)  ==  666669166671000001

   sumaDiagonales (10^2)  ==         671800
   sumaDiagonales (10^3)  ==        667168000
   sumaDiagonales (10^4)  ==       666716680000
   sumaDiagonales (10^5)  ==      666671666800000
   sumaDiagonales (10^6)  ==     666667166668000000
   sumaDiagonales (10^7)  ==    666666716666680000000
   sumaDiagonales (10^8)  ==   666666671666666800000000
   sumaDiagonales (10^9)  ==  666666667166666668000000000

sumaDiagonales 1 == 1

sumaDiagonales 2 == 10

sumaDiagonales 3 == 25

sumaDiagonales 4 == 56

sumaDiagonales 5 == 101

sumaDiagonales (10^6) == 666667166668000000

sumaDiagonales (1+10^6) == 666669166671000001

sumaDiagonales (10^2) == 671800

sumaDiagonales (10^3) == 667168000

sumaDiagonales (10^4) == 666716680000

sumaDiagonales (10^5) == 666671666800000

sumaDiagonales (10^6) == 666667166668000000

sumaDiagonales (10^7) == 666666716666680000000

sumaDiagonales (10^8) == 666666671666666800000000

sumaDiagonales (10^9) == 666666667166666668000000000

Comprobar con QuickCheck que el último dígito de (sumaDiagonales n) es 0, 4 ó 6 si n es par y es 1, 5 ó 7 en caso contrario.

Soluciones

import Test.QuickCheck (Positive (Positive), quickCheck)

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

sumaDiagonales1 :: Integer -> Integer
sumaDiagonales1 = sum . elementosEnDiagonales

-- (elementosEnDiagonales n) es la lista de los elementos en las
-- diagonales de la matriz espiral de orden nxn. Por ejemplo,
--    elementosEnDiagonales 1  ==  [1]
--    elementosEnDiagonales 2  ==  [1,2,3,4]
--    elementosEnDiagonales 3  ==  [1,3,5,7,9]
--    elementosEnDiagonales 4  ==  [1,2,3,4,7,10,13,16]
--    elementosEnDiagonales 5  ==  [1,3,5,7,9,13,17,21,25]
elementosEnDiagonales :: Integer -> [Integer]
elementosEnDiagonales n 
  | even n    = tail (scanl (+) 0 (concatMap (replicate 4) [1,3..n-1]))
  | otherwise = scanl (+) 1 (concatMap (replicate 4) [2,4..n-1])

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

sumaDiagonales2 :: Integer -> Integer
sumaDiagonales2 n
  | even n    = (-1) + n `div` 2 + sum [2*k^2-k+1 | k <- [0..n]]
  | otherwise = 1 + sum [4*k^2-6*k+6 | k <- [3,5..n]]

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

sumaDiagonales3 :: Integer -> Integer
sumaDiagonales3 n
  | even n    = n * (4*n^2 + 3*n + 8) `div` 6
  | otherwise = (4*n^3 + 3*n^2 + 8*n - 9) `div` 6

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

-- La propiedad es
prop_sumaDiagonales :: Positive Integer -> Bool
prop_sumaDiagonales (Positive n) =
  all (== sumaDiagonales1 n)
      [sumaDiagonales2 n,
       sumaDiagonales3 n]

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

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

-- La comparación es
--    λ> sumaDiagonales (2*10^6)
--    5333335333336000000
--    (2.30 secs, 1,521,955,848 bytes)
--    λ> sumaDiagonales2 (2*10^6)
--    5333335333336000000
--    (2.77 secs, 1,971,411,440 bytes)
--    λ> sumaDiagonales3 (2*10^6)
--    5333335333336000000
--    (0.01 secs, 139,520 bytes)

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

-- La propiedad es
prop_sumaDiagonales2 :: Positive Integer -> Bool
prop_sumaDiagonales2 (Positive n) 
  | even n    = x `elem` [0,4,6] 
  | otherwise = x `elem` [1,5,7] 
  where x = sumaDiagonales1 n `mod` 10

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

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

sumaDiagonales1 :: Integer -> Integer

sumaDiagonales1 = sum . elementosEnDiagonales

-- (elementosEnDiagonales n) es la lista de los elementos en las

-- diagonales de la matriz espiral de orden nxn. Por ejemplo,

-- elementosEnDiagonales 1 == [1]

-- elementosEnDiagonales 2 == [1,2,3,4]

-- elementosEnDiagonales 3 == [1,3,5,7,9]

-- elementosEnDiagonales 4 == [1,2,3,4,7,10,13,16]

-- elementosEnDiagonales 5 == [1,3,5,7,9,13,17,21,25]

elementosEnDiagonales :: Integer -> [Integer]

elementosEnDiagonales n

| even n = tail (scanl (+) 0 (concatMap (replicate 4) [1,3..n-1]))

| otherwise = scanl (+) 1 (concatMap (replicate 4) [2,4..n-1])

-- 2ª solución

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

sumaDiagonales2 :: Integer -> Integer

sumaDiagonales2 n

| even n = (-1) + n `div` 2 + sum [2*k^2-k+1 | k <- [0..n]]

| otherwise = 1 + sum [4*k^2-6*k+6 | k <- [3,5..n]]

-- 3ª solución

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

sumaDiagonales3 :: Integer -> Integer

sumaDiagonales3 n

| even n = n * (4*n^2 + 3*n + 8) `div` 6

| otherwise = (4*n^3 + 3*n^2 + 8*n - 9) `div` 6

-- Equivalencia de las definiciones

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

-- La propiedad es

prop_sumaDiagonales :: Positive Integer -> Bool

prop_sumaDiagonales (Positive n) =

all (== sumaDiagonales1 n)

[sumaDiagonales2 n,

sumaDiagonales3 n]

-- La comprobación es

-- λ> quickCheck prop_sumaDiagonales_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> sumaDiagonales (2*10^6)

-- 5333335333336000000

-- (2.30 secs, 1,521,955,848 bytes)

-- λ> sumaDiagonales2 (2*10^6)

-- 5333335333336000000

-- (2.77 secs, 1,971,411,440 bytes)

-- λ> sumaDiagonales3 (2*10^6)

-- 5333335333336000000

-- (0.01 secs, 139,520 bytes)

-- Propiedad

-- =========

-- La propiedad es

prop_sumaDiagonales2 :: Positive Integer -> Bool

prop_sumaDiagonales2 (Positive n)

| even n = x `elem` [0,4,6]

| otherwise = x `elem` [1,5,7]

where x = sumaDiagonales1 n `mod` 10

-- La comprobación es

-- λ> quickCheck prop_sumaDiagonales2

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Término ausente en una progresión aritmética

PorJosé A. Alonso 1-junio-202230-mayo-2022

Una progresión aritmética es una sucesión de números tales que la diferencia de dos términos sucesivos cualesquiera de la sucesión es constante.

Definir la función

   ausente :: Integral a => [a] -> a

1	ausente :: Integral a => [a] -> a

tal que (ausente xs) es el único término ausente de la progresión aritmética xs. Por ejemplo,

   ausente [3,7,9,11]               ==  5
   ausente [3,5,9,11]               ==  7
   ausente [3,5,7,11]               ==  9
   ausente ([1..9]++[11..])         ==  10
   ausente ([1..10^6] ++ [2+10^6])  ==  1000001

ausente [3,7,9,11] == 5

ausente [3,5,9,11] == 7

ausente [3,5,7,11] == 9

ausente ([1..9]++[11..]) == 10

ausente ([1..10^6] ++ [2+10^6]) == 1000001

Nota. Se supone que la lista tiene al menos 3 elementos, que puede ser infinita y que sólo hay un término de la progresión aritmética que no está en la lista.

Soluciones

import Data.List (group, genericLength)
import Test.QuickCheck (Arbitrary, Gen,
                        arbitrary, frequency, suchThat, quickCheck) 

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

ausente1 :: Integral a => [a] -> a
ausente1 (x1:xs@(x2:x3:_))
  | d1 == d2     = ausente1 xs
  | d1 == 2 * d2 = x1 + d2
  | d2 == 2 * d1 = x2 + d1
  where d1 = x2 - x1
        d2 = x3 - x2          

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

ausente2 :: Integral a => [a] -> a
ausente2 s@(x1:x2:x3:_) 
  | x1 + x3 /= 2 * x2 = x1 + (x3 - x2)
  | otherwise         = head [a | (a,b) <- zip [x1,x2..] s
                                , a /= b]

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

ausente3 :: Integral a => [a] -> a
ausente3  xs@(x1:x2:_) 
  | null us   = x1 + v
  | otherwise = x2 + u * genericLength (u:us) 
  where ((u:us):(v:_):_) = group (zipWith (-) (tail xs) xs)

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

-- Tipo de progresiones aritméticas con un término ausente.
newtype PAconAusente = PA [Integer]
  deriving Show

-- Generación de progresiones aritméticas con un término ausente. 
progresionConAusenteArbitraria :: Gen PAconAusente
progresionConAusenteArbitraria = do
  x <- arbitrary
  d <- arbitrary `suchThat` (/= 0)
  n <- arbitrary `suchThat` (> 2)
  k <- arbitrary `suchThat` (> 0)
  let (as,_:bs) = splitAt k [x,x+d..]
  frequency [(1,return (PA (as ++ bs))),
             (1,return (PA (take (length as + n) (as ++ bs))))]

-- Inclusión del tipo PAconAusente en Arbitrary.
instance Arbitrary PAconAusente where
  arbitrary = progresionConAusenteArbitraria

-- La propiedad es
prop_ausente :: PAconAusente -> Bool 
prop_ausente (PA xs) =
  all (== ausente1 xs)
      [ausente2 xs,
       ausente3 xs]

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

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

-- La comparación es
--    λ> let n = 10^6 in ausente1 ([1..n] ++ [n+2])
--    1000001
--    (1.15 secs, 560,529,520 bytes)
--    λ> let n = 10^6 in ausente2 ([1..n] ++ [n+2])
--    1000001
--    (0.33 secs, 336,530,680 bytes)
--    λ> let n = 10^6 in ausente3 ([1..n] ++ [n+2])
--    1000001
--    (0.50 secs, 498,047,584 bytes)

import Data.List (group, genericLength)

import Test.QuickCheck (Arbitrary, Gen,

arbitrary, frequency, suchThat, quickCheck)

-- 1ª solución

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

ausente1 :: Integral a => [a] -> a

ausente1 (x1:xs@(x2:x3:_))

| d1 == d2 = ausente1 xs

| d1 == 2 * d2 = x1 + d2

| d2 == 2 * d1 = x2 + d1

where d1 = x2 - x1

d2 = x3 - x2

-- 2ª solución

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

ausente2 :: Integral a => [a] -> a

ausente2 s@(x1:x2:x3:_)

| x1 + x3 /= 2 * x2 = x1 + (x3 - x2)

| otherwise = head [a | (a,b) <- zip [x1,x2..] s

, a /= b]

-- 3ª solución

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

ausente3 :: Integral a => [a] -> a

ausente3 xs@(x1:x2:_)

| null us = x1 + v

| otherwise = x2 + u * genericLength (u:us)

where ((u:us):(v:_):_) = group (zipWith (-) (tail xs) xs)

-- Comprobación de equivalencia

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

-- Tipo de progresiones aritméticas con un término ausente.

newtype PAconAusente = PA [Integer]

deriving Show

-- Generación de progresiones aritméticas con un término ausente.

progresionConAusenteArbitraria :: Gen PAconAusente

progresionConAusenteArbitraria = do

x <- arbitrary

d <- arbitrary `suchThat` (/= 0)

n <- arbitrary `suchThat` (> 2)

k <- arbitrary `suchThat` (> 0)

let (as,_:bs) = splitAt k [x,x+d..]

frequency [(1,return (PA (as ++ bs))),

(1,return (PA (take (length as + n) (as ++ bs))))]

-- Inclusión del tipo PAconAusente en Arbitrary.

instance Arbitrary PAconAusente where

arbitrary = progresionConAusenteArbitraria

-- La propiedad es

prop_ausente :: PAconAusente -> Bool

prop_ausente (PA xs) =

all (== ausente1 xs)

[ausente2 xs,

ausente3 xs]

-- La comprobación es

-- λ> quickCheck prop_ausente

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> let n = 10^6 in ausente1 ([1..n] ++ [n+2])

-- 1000001

-- (1.15 secs, 560,529,520 bytes)

-- λ> let n = 10^6 in ausente2 ([1..n] ++ [n+2])

-- 1000001

-- (0.33 secs, 336,530,680 bytes)

-- λ> let n = 10^6 in ausente3 ([1..n] ++ [n+2])

-- 1000001

-- (0.50 secs, 498,047,584 bytes)

El código se encuentra en GitHub.

Polinomios de Bell

PorJosé A. Alonso 31-mayo-202229-mayo-2022

Los polinomios de Bell forman una sucesión de polinomios, definida como sigue:

B₀(x) = 1 (polinomio unidad)
Bₙ(x) = x·[Bₙ(x) + Bₙ'(x)] (donde Bₙ'(x) es la derivada de Bₙ(x))

Por ejemplo,

   B₀(x) = 1                     = 1
   B₁(x) = x·(1+0)               = x     
   B₂(x) = x·(x+1)               = x²+x         
   B₃(x) = x·(x²+x+2x+1)         = x³+3x²+x    
   B₄(x) = x·(x³+3x²+x+3x²+6x+1) = x⁴+6x³+7x²+x

B₀(x) = 1 = 1

B₁(x) = x·(1+0) = x

B₂(x) = x·(x+1) = x²+x

B₃(x) = x·(x²+x+2x+1) = x³+3x²+x

B₄(x) = x·(x³+3x²+x+3x²+6x+1) = x⁴+6x³+7x²+x

Definir la función

   polBell :: Integer -> Polinomio Integer

1	polBell :: Integer -> Polinomio Integer

tal que (polBell n) es el polinomio de Bell de grado n. Por ejemplo,

   polBell 4                    ==  x^4 + 6*x^3 + 7*x^2 + 1*x
   coeficiente 2 (polBell 4)    ==  7
   coeficiente 2 (polBell 30)   ==  536870911
   coeficiente 1 (polBell 1000) == 1
   length (show (coeficiente 9 (polBell 2000)))  ==  1903

polBell 4 == x^4 + 6*x^3 + 7*x^2 + 1*x

coeficiente 2 (polBell 4) == 7

coeficiente 2 (polBell 30) == 536870911

coeficiente 1 (polBell 1000) == 1

length (show (coeficiente 9 (polBell 2000))) == 1903

Notas: Se usa la librería I1M.PolOperaciones que se encuentra aquí y se describe aquí. Además, en el último ejemplo se usa la función coeficiente tal que (coeficiente k p) es el coeficiente del término de grado k en el polinomio p definida por

   coeficiente :: Num a => Int -> Polinomio a -> a
   coeficiente k p | k == n                 = coefLider p
                   | k > grado (restoPol p) = 0
                   | otherwise              = coeficiente k (restoPol p)
                   where n = grado p

coeficiente :: Num a => Int -> Polinomio a -> a

coeficiente k p | k == n = coefLider p

| k > grado (restoPol p) = 0

| otherwise = coeficiente k (restoPol p)

where n = grado p

Soluciones

import Data.List          (genericIndex)
import I1M.PolOperaciones (Polinomio, coefLider, consPol, derivada,
                           grado, multPol, polCero, polUnidad, restoPol,
                           sumaPol) 
import Test.QuickCheck    (Positive (Positive), quickCheck)

-- Función auxiliar
-- ================

-- (coeficiente k p) es el coeficiente del término de grado k en el
-- polinomio p.
coeficiente :: Num a => Int -> Polinomio a -> a
coeficiente k p | k == n                 = coefLider p
                | k > grado (restoPol p) = 0
                | otherwise              = coeficiente k (restoPol p)
                where n = grado p

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

polBell1 :: Integer -> Polinomio Integer
polBell1 0 = polUnidad
polBell1 n = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))
  where p = polBell1 (n-1)

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

polBell2 :: Integer -> Polinomio Integer
polBell2 n = sucPolinomiosBell `genericIndex` n

sucPolinomiosBell :: [Polinomio Integer]
sucPolinomiosBell = iterate f polUnidad
  where f p = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))

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

-- La propiedad es
prop_polBell :: Positive Integer -> Bool 
prop_polBell (Positive n) =
  polBell1 n == polBell2 n

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

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

-- La comparación es
--    λ> length (show (coeficiente 9 (polBell1 2000)))
--    1903
--    (5.37 secs, 4,829,322,368 bytes)
--    λ> length (show (coeficiente 9 (polBell2 2000)))
--    1903
--    (4.03 secs, 4,825,094,064 bytes)

import Data.List (genericIndex)

import I1M.PolOperaciones (Polinomio, coefLider, consPol, derivada,

grado, multPol, polCero, polUnidad, restoPol,

sumaPol)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- Función auxiliar

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

-- (coeficiente k p) es el coeficiente del término de grado k en el

-- polinomio p.

coeficiente :: Num a => Int -> Polinomio a -> a

coeficiente k p | k == n = coefLider p

| k > grado (restoPol p) = 0

| otherwise = coeficiente k (restoPol p)

where n = grado p

-- 1ª solución

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

polBell1 :: Integer -> Polinomio Integer

polBell1 0 = polUnidad

polBell1 n = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))

where p = polBell1 (n-1)

-- 2ª solución

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

polBell2 :: Integer -> Polinomio Integer

polBell2 n = sucPolinomiosBell `genericIndex` n

sucPolinomiosBell :: [Polinomio Integer]

sucPolinomiosBell = iterate f polUnidad

where f p = multPol (consPol 1 1 polCero) (sumaPol p (derivada p))

-- Comprobación de equivalencia

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

-- La propiedad es

prop_polBell :: Positive Integer -> Bool

prop_polBell (Positive n) =

polBell1 n == polBell2 n

-- La comprobación es

-- λ> quickCheck prop_polBell

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (show (coeficiente 9 (polBell1 2000)))

-- 1903

-- (5.37 secs, 4,829,322,368 bytes)

-- λ> length (show (coeficiente 9 (polBell2 2000)))

-- 1903

-- (4.03 secs, 4,825,094,064 bytes)

El código se encuentra en GitHub.

Ordenación de los racionales

PorJosé A. Alonso 30-mayo-202229-mayo-2022

En este ejercicio, representamos las fracciones mediante pares de números de enteros.

Definir la función

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

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

tal que (fraccionesOrd n) es la lista con las fracciones propias positivas ordenadas, con denominador menor o igual que n. Por ejemplo,

   λ> fraccionesOrd 4
   [(1,4),(1,3),(1,2),(2,3),(3,4)]
   λ> fraccionesOrd 5
   [(1,5),(1,4),(1,3),(2,5),(1,2),(3,5),(2,3),(3,4),(4,5)]

λ> fraccionesOrd 4

[(1,4),(1,3),(1,2),(2,3),(3,4)]

λ> fraccionesOrd 5

[(1,5),(1,4),(1,3),(2,5),(1,2),(3,5),(2,3),(3,4),(4,5)]

Soluciones

import Data.List (sort, sortBy)
import Data.Ratio ((%), numerator, denominator)
import Test.QuickCheck

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

fraccionesOrd1 :: Integer -> [(Integer,Integer)]
fraccionesOrd1 n = 
  [(x,y) | (_,(x,y)) <- sort [(fromIntegral x/fromIntegral y,(x,y))
                              | y <- [2..n], 
                                x <- [1..y-1], 
                                gcd x y == 1]]

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

fraccionesOrd2 :: Integer -> [(Integer,Integer)]
fraccionesOrd2 n = 
  map snd (sort [(fromIntegral x/fromIntegral y,(x,y))
                 | y <- [2..n], 
                   x <- [1..y-1], 
                   gcd x y == 1])

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

fraccionesOrd3 :: Integer -> [(Integer,Integer)]
fraccionesOrd3 n = 
  sortBy comp [(x,y) | y <- [2..n], x <- [1..y-1], gcd x y == 1]
  where comp (a,b) (c,d) = compare (a*d) (b*c)

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

fraccionesOrd4 :: Integer -> [(Integer,Integer)]
fraccionesOrd4 n = 
  [(numerator x, denominator x) | x <- racionalesOrd4 n]
  
-- (racionalesOrd4 n) es la lista con los racionales ordenados, con
-- denominador menor o igual que n. Por ejemplo,
--    λ> racionalesOrd4 4
--    [1 % 4,1 % 3,1 % 2,2 % 3,3 % 4]
racionalesOrd4 :: Integer -> [Rational]
racionalesOrd4 n =
  sort [x % y | y <- [2..n], x <- [1..y-1], gcd x y == 1]

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

-- La propiedad es
prop_fraccionesOrd :: Positive Integer -> Bool 
prop_fraccionesOrd (Positive n) =
  all (== fraccionesOrd1 n)
      [fraccionesOrd2 n,
       fraccionesOrd3 n,
       fraccionesOrd4 n]

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

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

-- La comparación es
--    λ> length (fraccionesOrd1 2000)
--    1216587
--    (3.65 secs, 2,879,842,368 bytes)
--    λ> length (fraccionesOrd2 2000)
--    1216587
--    (3.36 secs, 2,870,109,640 bytes)
--    λ> length (fraccionesOrd3 2000)
--    1216587
--    (8.83 secs, 5,700,519,584 bytes)
--    λ> length (fraccionesOrd4 2000)
--    1216587
--    (4.12 secs, 5,181,904,336 bytes)

import Data.List (sort, sortBy)

import Data.Ratio ((%), numerator, denominator)

import Test.QuickCheck

-- 1ª solución

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

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

fraccionesOrd1 n =

[(x,y) | (_,(x,y)) <- sort [(fromIntegral x/fromIntegral y,(x,y))

| y <- [2..n],

x <- [1..y-1],

gcd x y == 1]]

-- 2ª solución

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

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

fraccionesOrd2 n =

map snd (sort [(fromIntegral x/fromIntegral y,(x,y))

| y <- [2..n],

x <- [1..y-1],

gcd x y == 1])

-- 3ª solución

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

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

fraccionesOrd3 n =

sortBy comp [(x,y) | y <- [2..n], x <- [1..y-1], gcd x y == 1]

where comp (a,b) (c,d) = compare (a*d) (b*c)

-- 4ª solución

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

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

fraccionesOrd4 n =

[(numerator x, denominator x) | x <- racionalesOrd4 n]

-- (racionalesOrd4 n) es la lista con los racionales ordenados, con

-- denominador menor o igual que n. Por ejemplo,

-- λ> racionalesOrd4 4

-- [1 % 4,1 % 3,1 % 2,2 % 3,3 % 4]

racionalesOrd4 :: Integer -> [Rational]

racionalesOrd4 n =

sort [x % y | y <- [2..n], x <- [1..y-1], gcd x y == 1]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_fraccionesOrd :: Positive Integer -> Bool

prop_fraccionesOrd (Positive n) =

all (== fraccionesOrd1 n)

[fraccionesOrd2 n,

fraccionesOrd3 n,

fraccionesOrd4 n]

-- La comprobación es

-- λ> quickCheck prop_fraccionesOrd

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (fraccionesOrd1 2000)

-- 1216587

-- (3.65 secs, 2,879,842,368 bytes)

-- λ> length (fraccionesOrd2 2000)

-- 1216587

-- (3.36 secs, 2,870,109,640 bytes)

-- λ> length (fraccionesOrd3 2000)

-- 1216587

-- (8.83 secs, 5,700,519,584 bytes)

-- λ> length (fraccionesOrd4 2000)

-- 1216587

-- (4.12 secs, 5,181,904,336 bytes)

El código se encuentra en GitHub.

Polinomios cuadráticos generadores de primos

PorJosé A. Alonso 27-mayo-2022

En 1772, Euler publicó que el polinomio n² + n + 41 genera 40 números primos para todos los valores de n entre 0 y 39. Sin embargo, cuando n = 40, 40²+40+41 = 40(40+1)+41 es divisible por 41.

Usando ordenadores, se descubrió que el polinomio n² – 79n + 1601 genera 80 números primos para todos los valores de n entre 0 y 79.

Definir la función

   generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

1	generadoresMaximales :: Integer -> (Int,[(Integer,Integer)])

tal que (generadoresMaximales n) es el par (m,xs) donde

xs es la lista de pares (x,y) tales que n²+xn+y es uno de polinomios que genera un número máximo de números primos consecutivos a partir de cero entre todos los polinomios de la forma n²+an+b, con |a| ≤ n y |b| ≤ n y
m es dicho número máximo.

Por ejemplo,

   generadoresMaximales    4  ==  ( 3,[(-2,3),(-1,3),(3,3)])
   generadoresMaximales    6  ==  ( 5,[(-1,5),(5,5)])
   generadoresMaximales   41  ==  (41,[(-1,41)])
   generadoresMaximales   50  ==  (43,[(-5,47)])
   generadoresMaximales  100  ==  (48,[(-15,97)])
   generadoresMaximales  200  ==  (53,[(-25,197)])
   generadoresMaximales 1650  ==  (80,[(-79,1601)])

generadoresMaximales 4 == ( 3,[(-2,3),(-1,3),(3,3)])

generadoresMaximales 6 == ( 5,[(-1,5),(5,5)])

generadoresMaximales 41 == (41,[(-1,41)])

generadoresMaximales 50 == (43,[(-5,47)])

generadoresMaximales 100 == (48,[(-15,97)])

generadoresMaximales 200 == (53,[(-25,197)])

generadoresMaximales 1650 == (80,[(-79,1601)])

Soluciones

import Data.List (sort)
import Data.Numbers.Primes (primes, isPrime)
import I1M.PolOperaciones (valor, consPol, polCero)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales1 n =
  (m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])
  where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]

-- (nPrimos a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos (-1) 41     ==  41
--    nPrimos (-79) 1601  ==  80
nPrimos :: Integer -> Integer -> Int
nPrimos a b =
  length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])

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

-- Notas:
-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que
--    0²+a*0+b = b es primo.
-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para
--    n = 1.

generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales2 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    a <- [-n..n],
                                    isPrime(1+a+b)]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

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

generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales3 n = (m,map snd zs)
  where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,
                                    p <- takeWhile (<=1+2*n) primes,
                                    let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- 4ª solución (con la librería de polinomios)
-- ===========================================

generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])
generadoresMaximales4 n = (m,map snd zs)
  where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,
                                     p <- takeWhile (<=1+2*n) primes,
                                     let a = p-b-1]
        ys = reverse (sort xs)
        m  = fst (head ys)
        zs = takeWhile (\(k,_) -> k == m) ys

-- (nPrimos2 a b) es el número de primos consecutivos generados por el
-- polinomio n² + an + b a partir de n=0. Por ejemplo,
--    nPrimos2 (-1) 41     ==  41
--    nPrimos2 (-79) 1601  ==  80
nPrimos2 :: Integer -> Integer -> Int
nPrimos2 a b =
  length (takeWhile isPrime [valor p n | n <- [0..]])
  where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))

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

-- La propiedad es
prop_generadoresMaximales :: Positive Integer -> Bool
prop_generadoresMaximales (Positive n) =
  all (equivalentes (generadoresMaximales1 n'))
      [generadoresMaximales2 n',
       generadoresMaximales3 n',
       generadoresMaximales4 n']
  where n' = n+1

equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool
equivalentes (n,xs) (m,ys) =
  n == m && sort xs == sort ys

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

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

-- La comparación es
--    λ> generadoresMaximales1 300
--    (56,[(-31,281)])
--    (2.10 secs, 2,744,382,760 bytes)
--    λ> generadoresMaximales2 300
--    (56,[(-31,281)])
--    (0.17 secs, 382,103,656 bytes)
--    λ> generadoresMaximales3 300
--    (56,[(-31,281)])
--    (0.19 secs, 346,725,872 bytes)
--    λ> generadoresMaximales4 300
--    (56,[(-31,281)])
--    (0.20 secs, 388,509,808 bytes)

100

101

102

103

104

105

106

107

108

import Data.List (sort)

import Data.Numbers.Primes (primes, isPrime)

import I1M.PolOperaciones (valor, consPol, polCero)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

generadoresMaximales1 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales1 n =

(m,[(a,b) | a <- [-n..n], b <- [-n..n], nPrimos a b == m])

where m = maximum [nPrimos a b | a <- [-n..n], b <- [-n..n]]

-- (nPrimos a b) es el número de primos consecutivos generados por el

-- polinomio n² + an + b a partir de n=0. Por ejemplo,

-- nPrimos (-1) 41 == 41

-- nPrimos (-79) 1601 == 80

nPrimos :: Integer -> Integer -> Int

nPrimos a b =

length (takeWhile isPrime [n*n+a*n+b | n <- [0..]])

-- 2ª solución

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

-- Notas:

-- 1. Se tiene que b es primo, ya que para n = 0, se tiene que

-- 0²+a*0+b = b es primo.

-- 2. Se tiene que 1+a+b es primo, ya que es el valor del polinomio para

-- n = 1.

generadoresMaximales2 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales2 n = (m,map snd zs)

where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,

a <- [-n..n],

isPrime(1+a+b)]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- 3ª solución

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

generadoresMaximales3 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales3 n = (m,map snd zs)

where xs = [(nPrimos a b,(a,b)) | b <- takeWhile (<=n) primes,

p <- takeWhile (<=1+2*n) primes,

let a = p-b-1]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- 4ª solución (con la librería de polinomios)

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

generadoresMaximales4 :: Integer -> (Int,[(Integer,Integer)])

generadoresMaximales4 n = (m,map snd zs)

where xs = [(nPrimos2 a b,(a,b)) | b <- takeWhile (<=n) primes,

p <- takeWhile (<=1+2*n) primes,

let a = p-b-1]

ys = reverse (sort xs)

m = fst (head ys)

zs = takeWhile (\(k,_) -> k == m) ys

-- (nPrimos2 a b) es el número de primos consecutivos generados por el

-- polinomio n² + an + b a partir de n=0. Por ejemplo,

-- nPrimos2 (-1) 41 == 41

-- nPrimos2 (-79) 1601 == 80

nPrimos2 :: Integer -> Integer -> Int

nPrimos2 a b =

length (takeWhile isPrime [valor p n | n <- [0..]])

where p = consPol 2 1 (consPol 1 a (consPol 0 b polCero))

-- Comprobación de equivalencia

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

-- La propiedad es

prop_generadoresMaximales :: Positive Integer -> Bool

prop_generadoresMaximales (Positive n) =

all (equivalentes (generadoresMaximales1 n'))

[generadoresMaximales2 n',

generadoresMaximales3 n',

generadoresMaximales4 n']

where n' = n+1

equivalentes :: (Int,[(Integer,Integer)]) -> (Int,[(Integer,Integer)]) -> Bool

equivalentes (n,xs) (m,ys) =

n == m && sort xs == sort ys

-- La comprobación es

-- λ> quickCheck prop_generadoresMaximales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> generadoresMaximales1 300

-- (56,[(-31,281)])

-- (2.10 secs, 2,744,382,760 bytes)

-- λ> generadoresMaximales2 300

-- (56,[(-31,281)])

-- (0.17 secs, 382,103,656 bytes)

-- λ> generadoresMaximales3 300

-- (56,[(-31,281)])

-- (0.19 secs, 346,725,872 bytes)

-- λ> generadoresMaximales4 300

-- (56,[(-31,281)])

-- (0.20 secs, 388,509,808 bytes)

El código se encuentra en GitHub.

El triángulo de Floyd

PorJosé A. Alonso 26-mayo-2022

El triángulo de Floyd, llamado así en honor a Robert Floyd, es un triángulo rectángulo formado con números naturales. Para crear un triángulo de Floyd, se comienza con un 1 en la esquina superior izquierda, y se continúa escribiendo la secuencia de los números naturales de manera que cada línea contenga un número más que la anterior. Las 5 primeras líneas del triángulo de Floyd son

    1
    2   3
    4   5   6
    7   8   9  10
   11  12  13  14  15

2 3

4 5 6

7 8 9 10

11 12 13 14 15

Definir la función

   trianguloFloyd :: [[Integer]]

1	trianguloFloyd :: [[Integer]]

tal que trianguloFloyd es el triángulo de Floyd. Por ejemplo,

   λ> take 4 trianguloFloyd
   [[1],
    [2,3],
    [4,5,6],
    [7,8,9,10]]
  (trianguloFloyd !! (10^5)) !! 0  ==  5000050001
  (trianguloFloyd !! (10^6)) !! 0  ==  500000500001
  (trianguloFloyd !! (10^7)) !! 0  ==  50000005000001

λ> take 4 trianguloFloyd

[[1],

[2,3],

[4,5,6],

[7,8,9,10]]

(trianguloFloyd !! (10^5)) !! 0 == 5000050001

(trianguloFloyd !! (10^6)) !! 0 == 500000500001

(trianguloFloyd !! (10^7)) !! 0 == 50000005000001

Soluciones

import Data.List (genericLength)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

trianguloFloyd1 :: [[Integer]]
trianguloFloyd1 = floyd 1 [1..]
  where floyd n xs = i : floyd (n+1) r
          where (i,r) = splitAt n xs

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

trianguloFloyd2 :: [[Integer]]
trianguloFloyd2 = iterate siguienteF [1]

-- (siguienteF xs) es la lista de los elementos de la línea xs en el
-- triángulo de Floyd. Por ejemplo,
--    siguienteF [2,3]    ==  [4,5,6]
--    siguienteF [4,5,6]  ==  [7,8,9,10]
siguienteF :: [Integer] -> [Integer]
siguienteF xs = [a..a+n]
    where a = 1 + last xs
          n = genericLength xs

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

trianguloFloyd3 :: [[Integer]]
trianguloFloyd3 =
  [[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]

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

trianguloFloyd4 :: [[Integer]]
trianguloFloyd4 =
  scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]

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

-- La propiedad es
prop_trianguloFloyd :: Positive Int -> Bool
prop_trianguloFloyd (Positive n) =
  all (== (trianguloFloyd1 !! n))
      [trianguloFloyd2 !! n,
       trianguloFloyd3 !! n,
       trianguloFloyd4 !! n]

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

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

-- La comparación es
--    λ> (trianguloFloyd1 !! 5000) !! 5000
--    12507501
--    (1.47 secs, 2,505,005,752 bytes)
--    λ> (trianguloFloyd2 !! 5000) !! 5000
--    12507501
--    (0.79 secs, 2,416,259,176 bytes)
--    λ> (trianguloFloyd3 !! 5000) !! 5000
--    12507501
--    (0.00 secs, 1,809,152 bytes)
--    λ> (trianguloFloyd4 !! 5000) !! 5000
--    12507501
--    (0.01 secs, 3,517,896 bytes)
--
--    λ> (trianguloFloyd3 !! (10^7)) !! 0
--    50000005000001
--    (2.45 secs, 1,656,534,080 bytes)
--    λ> (trianguloFloyd4 !! (10^7)) !! 0
--    50000005000001
--    (10.86 secs, 5,302,760,752 bytes)

import Data.List (genericLength)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

trianguloFloyd1 :: [[Integer]]

trianguloFloyd1 = floyd 1 [1..]

where floyd n xs = i : floyd (n+1) r

where (i,r) = splitAt n xs

-- 2ª solución

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

trianguloFloyd2 :: [[Integer]]

trianguloFloyd2 = iterate siguienteF [1]

-- (siguienteF xs) es la lista de los elementos de la línea xs en el

-- triángulo de Floyd. Por ejemplo,

-- siguienteF [2,3] == [4,5,6]

-- siguienteF [4,5,6] == [7,8,9,10]

siguienteF :: [Integer] -> [Integer]

siguienteF xs = [a..a+n]

where a = 1 + last xs

n = genericLength xs

-- 3ª solución

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

trianguloFloyd3 :: [[Integer]]

trianguloFloyd3 =

[[(n*(n-1) `div` 2) + 1 .. (n*(n+1) `div` 2)] | n <- [1..]]

-- 4ª solución

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

trianguloFloyd4 :: [[Integer]]

trianguloFloyd4 =

scanl (\(x:_) y -> [x+y..x+2*y]) [1] [1..]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_trianguloFloyd :: Positive Int -> Bool

prop_trianguloFloyd (Positive n) =

all (== (trianguloFloyd1 !! n))

[trianguloFloyd2 !! n,

trianguloFloyd3 !! n,

trianguloFloyd4 !! n]

-- La comprobación es

-- λ> quickCheck prop_trianguloFloyd

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> (trianguloFloyd1 !! 5000) !! 5000

-- 12507501

-- (1.47 secs, 2,505,005,752 bytes)

-- λ> (trianguloFloyd2 !! 5000) !! 5000

-- 12507501

-- (0.79 secs, 2,416,259,176 bytes)

-- λ> (trianguloFloyd3 !! 5000) !! 5000

-- 12507501

-- (0.00 secs, 1,809,152 bytes)

-- λ> (trianguloFloyd4 !! 5000) !! 5000

-- 12507501

-- (0.01 secs, 3,517,896 bytes)

-- λ> (trianguloFloyd3 !! (10^7)) !! 0

-- 50000005000001

-- (2.45 secs, 1,656,534,080 bytes)

-- λ> (trianguloFloyd4 !! (10^7)) !! 0

-- 50000005000001

-- (10.86 secs, 5,302,760,752 bytes)

El código se encuentra en GitHub.

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Entradas recientes

Buscar

Correo electrónico