junio 2022 – Exercitium

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.

Meses: junio 2022

Aproximación del número pi

Soluciones

Pandigitales primos

Soluciones

Codificación de Fibonacci

Soluciones

La sucesión del reloj astronómico de Praga

Soluciones

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

Soluciones

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

Soluciones

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

Soluciones

Cálculo aproximado de integrales definidas

Soluciones

Menor número con una cantidad dada de divisores

Soluciones

Distancia esperada entre dos puntos de un cuadrado unitario

Soluciones

Representación matricial de relaciones binarias

Soluciones

Codificación de Gödel

Soluciones

Primos circulares

Soluciones

Diccionario de frecuencias

Soluciones

Descomposiciones con sumandos 1 ó 2

Soluciones

Suma de los elementos de las diagonales matrices espirales

Soluciones

Término ausente en una progresión aritmética

Soluciones

Entradas recientes

Buscar

Correo electrónico