Gráficas – Página 2

Las sucesiones de Loomis

PorJosé A. Alonso 30-abril-201825-abril-2021

La sucesión de Loomis generada por un número entero positivo x es la sucesión cuyos términos se definen por

f(0) es x
f(n) es la suma de f(n-1) y el producto de los dígitos no nulos de f(n-1)

Los primeros términos de las primeras sucesiones de Loomis son

Generada por 1: 1, 2, 4, 8, 16, 22, 26, 38, 62, 74, 102, 104, 108, 116, 122, …
Generada por 2: 2, 4, 8, 16, 22, 26, 38, 62, 74, 102, 104, 108, 116, 122, 126, …
Generada por 3: 3, 6, 12, 14, 18, 26, 38, 62, 74, 102, 104, 108, 116, 122, 126, …
Generada por 4: 4, 8, 16, 22, 26, 38, 62, 74, 102, 104, 108, 116, 122, 126, 138, …
Generada por 5: 5, 10, 11, 12, 14, 18, 26, 38, 62, 74, 102, 104, 108, 116, 122, …

Se observa que a partir de un término todas coinciden con la generada por 1. Dicho término se llama el punto de convergencia. Por ejemplo,

la generada por 2 converge a 2
la generada por 3 converge a 26
la generada por 4 converge a 4
la generada por 5 converge a 26

Definir las siguientes funciones

   sucLoomis           :: Integer -> [Integer]
   convergencia        :: Integer -> Integer
   graficaConvergencia :: [Integer] -> IO ()

sucLoomis :: Integer -> [Integer]

convergencia :: Integer -> Integer

graficaConvergencia :: [Integer] -> IO ()

tales que

(sucLoomis x) es la sucesión de Loomis generada por x. Por ejemplo,

     λ> take 15 (sucLoomis 1)
     [1,2,4,8,16,22,26,38,62,74,102,104,108,116,122]
     λ> take 15 (sucLoomis 2)
     [2,4,8,16,22,26,38,62,74,102,104,108,116,122,126]
     λ> take 15 (sucLoomis 3)
     [3,6,12,14,18,26,38,62,74,102,104,108,116,122,126]
     λ> take 15 (sucLoomis 4)
     [4,8,16,22,26,38,62,74,102,104,108,116,122,126,138]
     λ> take 15 (sucLoomis 5)
     [5,10,11,12,14,18,26,38,62,74,102,104,108,116,122]
     λ> take 15 (sucLoomis 20)
     [20,22,26,38,62,74,102,104,108,116,122,126,138,162,174]
     λ> take 15 (sucLoomis 100)
     [100,101,102,104,108,116,122,126,138,162,174,202,206,218,234]
     λ> sucLoomis 1 !! (2*10^5)
     235180736652

λ> take 15 (sucLoomis 1)

[1,2,4,8,16,22,26,38,62,74,102,104,108,116,122]

λ> take 15 (sucLoomis 2)

[2,4,8,16,22,26,38,62,74,102,104,108,116,122,126]

λ> take 15 (sucLoomis 3)

[3,6,12,14,18,26,38,62,74,102,104,108,116,122,126]

λ> take 15 (sucLoomis 4)

[4,8,16,22,26,38,62,74,102,104,108,116,122,126,138]

λ> take 15 (sucLoomis 5)

[5,10,11,12,14,18,26,38,62,74,102,104,108,116,122]

λ> take 15 (sucLoomis 20)

[20,22,26,38,62,74,102,104,108,116,122,126,138,162,174]

λ> take 15 (sucLoomis 100)

[100,101,102,104,108,116,122,126,138,162,174,202,206,218,234]

λ> sucLoomis 1 !! (2*10^5)

235180736652

(convergencia x) es el término de convergencia de la sucesioń de Loomis generada por x xon la geerada por 1. Por ejemplo,

     convergencia  2      ==  2
     convergencia  3      ==  26
     convergencia  4      ==  4
     convergencia 17      ==  38
     convergencia 19      ==  102
     convergencia 43      ==  162
     convergencia 27      ==  202
     convergencia 58      ==  474
     convergencia 63      ==  150056
     convergencia 81      ==  150056
     convergencia 89      ==  150056
     convergencia (10^12) ==  1000101125092

convergencia 2 == 2

convergencia 3 == 26

convergencia 4 == 4

convergencia 17 == 38

convergencia 19 == 102

convergencia 43 == 162

convergencia 27 == 202

convergencia 58 == 474

convergencia 63 == 150056

convergencia 81 == 150056

convergencia 89 == 150056

convergencia (10^12) == 1000101125092

(graficaConvergencia xs) dibuja la gráfica de los términos de convergencia de las sucesiones de Loomis generadas por los elementos de xs. Por ejemplo, (graficaConvergencia ([1..50]) dibuja

y graficaConvergencia ([1..148] \ [63,81,89,137]) dibuja

Soluciones

import Data.List               ((\\))
import Data.Char               (digitToInt)
import Graphics.Gnuplot.Simple (plotList, Attribute (Key, Title, XRange, PNG))

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

sucLoomis :: Integer -> [Integer]
sucLoomis x = map (loomis x) [0..]

loomis :: Integer -> Integer -> Integer
loomis x 0 = x
loomis x n = y + productoDigitosNoNulos y
  where y = loomis x (n-1)

productoDigitosNoNulos :: Integer -> Integer
productoDigitosNoNulos = product . digitosNoNulos

digitosNoNulos :: Integer -> [Integer]
digitosNoNulos x =
  [read [c] | c <- show x, c /= '0']

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

sucLoomis2 :: Integer -> [Integer]
sucLoomis2 = iterate siguienteLoomis 

siguienteLoomis :: Integer -> Integer
siguienteLoomis y = y + productoDigitosNoNulos y

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

sucLoomis3 :: Integer -> [Integer]
sucLoomis3 =
  iterate ((+) <*> product .
           map (toInteger . digitToInt) .
           filter (/= '0') . show)


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

--    λ> sucLoomis 1 !! 30000
--    6571272766
--    (2.45 secs, 987,955,944 bytes)
--    λ> sucLoomis2 1 !! 30000
--    6571272766
--    (2.26 secs, 979,543,328 bytes)
--    λ> sucLoomis3 1 !! 30000
--    6571272766
--    (0.31 secs, 88,323,832 bytes)

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

convergencia1 :: Integer -> Integer
convergencia1 x =
  head (dropWhile noEnSucLoomisDe1 (sucLoomis x))

noEnSucLoomisDe1 :: Integer -> Bool
noEnSucLoomisDe1 x = not (pertenece x sucLoomisDe1)

sucLoomisDe1 :: [Integer]
sucLoomisDe1 = sucLoomis 1

pertenece :: Integer -> [Integer] -> Bool
pertenece x ys =
  x == head (dropWhile (<x) ys)

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

convergencia2 :: Integer -> Integer
convergencia2 = aux (sucLoomis3 1) . sucLoomis3
 where aux as@(x:xs) bs@(y:ys) | x == y    = x
                               | x < y     = aux xs bs
                               | otherwise = aux as ys

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

convergencia3 :: Integer -> Integer
convergencia3 = head . interseccion (sucLoomis3 1) . sucLoomis3
 
-- (interseccion xs ys) es la intersección entre las listas ordenadas xs
-- e ys. Por ejemplo,
--    λ> take 10 (interseccion (sucLoomis3 1) (sucLoomis3 2))
--    [2,4,8,16,22,26,38,62,74,102]
interseccion :: Ord a => [a] -> [a] -> [a]
interseccion = aux
  where aux as@(x:xs) bs@(y:ys) = case compare x y of
                                    LT ->     aux xs bs
                                    EQ -> x : aux xs ys
                                    GT ->     aux as ys
        aux _         _         = []                           

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

convergencia4 :: Integer -> Integer
convergencia4 x = perteneceA (sucLoomis3 x) 1
  where perteneceA (y:ys) n | y == c    = y
                            | otherwise = perteneceA ys c
          where c = head $ dropWhile (< y) $ sucLoomis3 n

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

--    λ> convergencia1 (10^4)
--    150056
--    (2.94 secs, 1,260,809,808 bytes)
--    λ> convergencia2 (10^4)
--    150056
--    (0.03 secs, 700,240 bytes)
--    λ> convergencia3 (10^4)
--    150056
--    (0.03 secs, 1,165,496 bytes)
--    λ> convergencia4 (10^4)
--    150056
--    (0.02 secs, 1,119,648 bytes)
--    
--    λ> convergencia2 (10^12)
--    1000101125092
--    (1.81 secs, 714,901,080 bytes)
--    λ> convergencia3 (10^12)
--    1000101125092
--    (1.92 secs, 744,932,184 bytes)
--    λ> convergencia4 (10^12)
--    1000101125092
--    (1.82 secs, 941,053,328 bytes)

-- Definición de graficaConvergencia
-- ==================================

graficaConvergencia :: [Integer] -> IO ()
graficaConvergencia xs =
  plotList [ Key Nothing
           , Title "Convergencia de sucesiones de Loomis"
           , XRange (fromIntegral (minimum xs),fromIntegral (maximum xs))
           , PNG "Las_sucesiones_de_Loomis_2.png"
           ]
           [(x,convergencia2 x) | x <- xs]

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

import Data.List ((\\))

import Data.Char (digitToInt)

import Graphics.Gnuplot.Simple (plotList, Attribute (Key, Title, XRange, PNG))

-- 1ª definición de sucLoomis

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

sucLoomis :: Integer -> [Integer]

sucLoomis x = map (loomis x) [0..]

loomis :: Integer -> Integer -> Integer

loomis x 0 = x

loomis x n = y + productoDigitosNoNulos y

where y = loomis x (n-1)

productoDigitosNoNulos :: Integer -> Integer

productoDigitosNoNulos = product . digitosNoNulos

digitosNoNulos :: Integer -> [Integer]

digitosNoNulos x =

[read [c] | c <- show x, c /= '0']

-- 2ª definición de sucLoomis

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

sucLoomis2 :: Integer -> [Integer]

sucLoomis2 = iterate siguienteLoomis

siguienteLoomis :: Integer -> Integer

siguienteLoomis y = y + productoDigitosNoNulos y

-- 3ª definición de sucLoomis

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

sucLoomis3 :: Integer -> [Integer]

sucLoomis3 =

iterate ((+) <*> product .

map (toInteger . digitToInt) .

filter (/= '0') . show)

-- Comparación de eficiencia

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

-- λ> sucLoomis 1 !! 30000

-- 6571272766

-- (2.45 secs, 987,955,944 bytes)

-- λ> sucLoomis2 1 !! 30000

-- 6571272766

-- (2.26 secs, 979,543,328 bytes)

-- λ> sucLoomis3 1 !! 30000

-- 6571272766

-- (0.31 secs, 88,323,832 bytes)

-- 1ª definición de convergencia

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

convergencia1 :: Integer -> Integer

convergencia1 x =

head (dropWhile noEnSucLoomisDe1 (sucLoomis x))

noEnSucLoomisDe1 :: Integer -> Bool

noEnSucLoomisDe1 x = not (pertenece x sucLoomisDe1)

sucLoomisDe1 :: [Integer]

sucLoomisDe1 = sucLoomis 1

pertenece :: Integer -> [Integer] -> Bool

pertenece x ys =

x == head (dropWhile (<x) ys)

-- 2ª definición de convergencia

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

convergencia2 :: Integer -> Integer

convergencia2 = aux (sucLoomis3 1) . sucLoomis3

where aux as@(x:xs) bs@(y:ys) | x == y = x

| x < y = aux xs bs

| otherwise = aux as ys

-- 3ª definición de convergencia

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

convergencia3 :: Integer -> Integer

convergencia3 = head . interseccion (sucLoomis3 1) . sucLoomis3

-- (interseccion xs ys) es la intersección entre las listas ordenadas xs

-- e ys. Por ejemplo,

-- λ> take 10 (interseccion (sucLoomis3 1) (sucLoomis3 2))

-- [2,4,8,16,22,26,38,62,74,102]

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

interseccion = aux

where aux as@(x:xs) bs@(y:ys) = case compare x y of

LT -> aux xs bs

EQ -> x : aux xs ys

GT -> aux as ys

aux _ _ = []

-- 4ª definición de convergencia

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

convergencia4 :: Integer -> Integer

convergencia4 x = perteneceA (sucLoomis3 x) 1

where perteneceA (y:ys) n | y == c = y

| otherwise = perteneceA ys c

where c = head $ dropWhile (< y) $ sucLoomis3 n

-- Comparación de eficiencia

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

-- λ> convergencia1 (10^4)

-- 150056

-- (2.94 secs, 1,260,809,808 bytes)

-- λ> convergencia2 (10^4)

-- 150056

-- (0.03 secs, 700,240 bytes)

-- λ> convergencia3 (10^4)

-- 150056

-- (0.03 secs, 1,165,496 bytes)

-- λ> convergencia4 (10^4)

-- 150056

-- (0.02 secs, 1,119,648 bytes)

-- λ> convergencia2 (10^12)

-- 1000101125092

-- (1.81 secs, 714,901,080 bytes)

-- λ> convergencia3 (10^12)

-- 1000101125092

-- (1.92 secs, 744,932,184 bytes)

-- λ> convergencia4 (10^12)

-- 1000101125092

-- (1.82 secs, 941,053,328 bytes)

-- Definición de graficaConvergencia

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

graficaConvergencia :: [Integer] -> IO ()

graficaConvergencia xs =

plotList [ Key Nothing

, Title "Convergencia de sucesiones de Loomis"

, XRange (fromIntegral (minimum xs),fromIntegral (maximum xs))

, PNG "Las_sucesiones_de_Loomis_2.png"

]

[(x,convergencia2 x) | x <- xs]

Alturas primas

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

Se considera una enumeración de los números primos:

    p(1)=2,  p(2)=3, p(3)=5, p(4)=7, p(5)=11, p(6)=13, p(7)=17,...

1	p(1)=2, p(2)=3, p(3)=5, p(4)=7, p(5)=11, p(6)=13, p(7)=17,...

Dado un entero x > 1, su altura prima es el mayor i tal que el primo p(i) aparece en la factorización de x en números primos. Por ejemplo, la altura prima de 3500 tiene longitud 4, pues 3500=2^2×5^3×7^1 y la de 34 tiene es 7, pues 34 = 2×17. Además, se define la altura prima de 1 como 0.

Definir las funciones

   alturaPrima        :: Integer -> Integer
   alturasPrimas      :: Integer -> [Integer]
   graficaAlturaPrima :: Integer -> IO ()

alturaPrima :: Integer -> Integer

alturasPrimas :: Integer -> [Integer]

graficaAlturaPrima :: Integer -> IO ()

tales que

(alturaPrima x) es la altura prima de x. Por ejemplo,

     alturaPrima 3500  ==  4
     alturaPrima 34    ==  7

1 2	alturaPrima 3500 == 4 alturaPrima 34 == 7

(alturasPrimas n) es la lista de las altura prima de los primeros n números enteros positivos. Por ejemplo,

     alturasPrimas 15  ==  [0,1,2,1,3,2,4,1,2,3,5,2,6,4,3]
     maximum (alturasPrimas 10000)  ==  1229
     maximum (alturasPrimas 20000)  ==  2262
     maximum (alturasPrimas 30000)  ==  3245
     maximum (alturasPrimas 40000)  ==  4203

alturasPrimas 15 == [0,1,2,1,3,2,4,1,2,3,5,2,6,4,3]

maximum (alturasPrimas 10000) == 1229

maximum (alturasPrimas 20000) == 2262

maximum (alturasPrimas 30000) == 3245

maximum (alturasPrimas 40000) == 4203

(graficaAlturaPrima n) dibuja las alturas primas de los números entre 2 y n. Por ejemplo, (graficaAlturaPrima 500) dibuja

Soluciones

import Data.List (genericLength)
import Data.Numbers.Primes (isPrime, primes, primeFactors)
import Data.Array
import Graphics.Gnuplot.Simple

-- 1ª definicioń de alturaPrima
-- ============================

alturaPrima :: Integer -> Integer
alturaPrima 1 = 0
alturaPrima n = indice (mayorFactorPrimo n)

-- (mayorFactorPrimo n) es el mayor factor primo de n. Por ejemplo,
--    mayorFactorPrimo 3500  ==  7
--    mayorFactorPrimo 34    ==  17
mayorFactorPrimo :: Integer -> Integer
mayorFactorPrimo = last . primeFactors

-- (indice p) es el índice de p en la sucesión de los números
-- primos. Por ejemplo,
--    indice 7   ==  4
--    indice 17  ==  7
indice :: Integer -> Integer
indice p = genericLength (takeWhile (<=p) primes)

-- 2ª definicioń de alturaPrima
-- ============================

alturaPrima2 :: Integer -> Integer
alturaPrima2 n = v ! n
  where v = array (1,n) [(i,f i) | i <- [1..n]]
        f 1 = 0
        f k | isPrime k = indice2 k
            | otherwise = v ! (k `div` (head (primeFactors k)))

indice2 :: Integer -> Integer
indice2 p = head [n | (x,n) <- indicesPrimos, x == p]

-- indicesPrimos es la suceción formada por los números primos y sus
-- índices. Por ejemplo,
--    λ> take 10 indicesPrimos
--    [(2,1),(3,2),(5,3),(7,4),(11,5),(13,6),(17,7),(19,8),(23,9),(29,10)]
indicesPrimos :: [(Integer,Integer)]
indicesPrimos = zip primes [1..]

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

alturasPrimas :: Integer -> [Integer]
alturasPrimas n = map alturaPrima [1..n]

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

alturasPrimas2 :: Integer -> [Integer]
alturasPrimas2 n = elems v 
  where v = array (1,n) [(i,f i) | i <- [1..n]]
        f 1 = 0
        f k | isPrime k = indice2 k
            | otherwise = v ! (k `div` (head (primeFactors k)))

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

--    λ> maximum (alturasPrimas 20000)
--    2262
--    (29.97 secs, 13,179,984,536 bytes)
--    λ> maximum (alturasPrimas2 20000)
--    2262
--    (2.11 secs, 455,259,448 bytes)

-- Definición de graficaAlturaPrima
-- ================================

graficaAlturaPrima :: Integer -> IO ()
graficaAlturaPrima n =
  plotList [ Key Nothing
           , PNG "Alturas_primas.png"
           ]
           (alturasPrimas2 n)

import Data.List (genericLength)

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

import Data.Array

import Graphics.Gnuplot.Simple

-- 1ª definicioń de alturaPrima

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

alturaPrima :: Integer -> Integer

alturaPrima 1 = 0

alturaPrima n = indice (mayorFactorPrimo n)

-- (mayorFactorPrimo n) es el mayor factor primo de n. Por ejemplo,

-- mayorFactorPrimo 3500 == 7

-- mayorFactorPrimo 34 == 17

mayorFactorPrimo :: Integer -> Integer

mayorFactorPrimo = last . primeFactors

-- (indice p) es el índice de p en la sucesión de los números

-- primos. Por ejemplo,

-- indice 7 == 4

-- indice 17 == 7

indice :: Integer -> Integer

indice p = genericLength (takeWhile (<=p) primes)

-- 2ª definicioń de alturaPrima

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

alturaPrima2 :: Integer -> Integer

alturaPrima2 n = v ! n

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

f 1 = 0

f k | isPrime k = indice2 k

| otherwise = v ! (k `div` (head (primeFactors k)))

indice2 :: Integer -> Integer

indice2 p = head [n | (x,n) <- indicesPrimos, x == p]

-- indicesPrimos es la suceción formada por los números primos y sus

-- índices. Por ejemplo,

-- λ> take 10 indicesPrimos

-- [(2,1),(3,2),(5,3),(7,4),(11,5),(13,6),(17,7),(19,8),(23,9),(29,10)]

indicesPrimos :: [(Integer,Integer)]

indicesPrimos = zip primes [1..]

-- 1ª definición de alturasPrimas

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

alturasPrimas :: Integer -> [Integer]

alturasPrimas n = map alturaPrima [1..n]

-- 2ª definición de alturasPrimas

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

alturasPrimas2 :: Integer -> [Integer]

alturasPrimas2 n = elems v

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

f 1 = 0

f k | isPrime k = indice2 k

| otherwise = v ! (k `div` (head (primeFactors k)))

-- Comparación de eficiencia

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

-- λ> maximum (alturasPrimas 20000)

-- 2262

-- (29.97 secs, 13,179,984,536 bytes)

-- λ> maximum (alturasPrimas2 20000)

-- 2262

-- (2.11 secs, 455,259,448 bytes)

-- Definición de graficaAlturaPrima

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

graficaAlturaPrima :: Integer -> IO ()

graficaAlturaPrima n =

plotList [ Key Nothing

, PNG "Alturas_primas.png"

]

(alturasPrimas2 n)

Sucesiones de números consecutivos con suma dada

PorJosé A. Alonso 29-marzo-20187-abril-2018

El número 15 se puede escribir de 5 formas como suma de números naturales consecutivos:

   15 = 0+1+2+3+4+5 = 1+2+3+4+5 = 4+5+6 = 7+8 = 15.

1	15 = 0+1+2+3+4+5 = 1+2+3+4+5 = 4+5+6 = 7+8 = 15.

Definir las funciones

   sucesionesConSuma        :: Int -> [(Int,Int)]
   graficaSucesionesConSuma :: Int -> IO ()

1 2	sucesionesConSuma :: Int -> [(Int,Int)] graficaSucesionesConSuma :: Int -> IO ()

tales que

(sucesionesConSuma n) es la lista de los pares formados por el primero y por el último elemento de las sucesiones de números naturales consecutivos con suma n. Por ejemplo,

     sucesionesConSuma 15             ==  [(0,5),(1,5),(4,6),(7,8),(15,15)]
     length (sucesionesConSuma 3000)  ==  8

1 2	sucesionesConSuma 15 == [(0,5),(1,5),(4,6),(7,8),(15,15)] length (sucesionesConSuma 3000) == 8

(graficaSucesionesConSuma n) dibuja la gráfica del número de formas de escribir los n primeros números como suma de números naturales consecutivos. Por ejemplo, (graficaSucesionesConSuma 100) dibuja

Soluciones

import Graphics.Gnuplot.Simple

-- 1ª solución
sucesionesConSuma :: Int -> [(Int,Int)]
sucesionesConSuma n =
  [(x,y) | y <- [0..n], x <- [0..y], sum [x..y] == n]

-- 2ª solución
sucesionesConSuma2 :: Int -> [(Int,Int)]
sucesionesConSuma2 n = 
    [(x,y) | y <- [0..n], x <- [0..y], (x+y)*(y-x+1) == 2*n]

-- Comparación de eficiencia
--    λ> sucesionesConSuma 700
--    [(3,37),(16,40),(84,91),(97,103),(138,142),(700,700)]
--    (2.71 secs, 7,894,718,264 bytes)
--    λ> sucesionesConSuma2 700
--    [(3,37),(16,40),(84,91),(97,103),(138,142),(700,700)]
--    (0.22 secs, 118,104,176 bytes)

-- Gráfica
graficaSucesionesConSuma :: Int -> IO ()
graficaSucesionesConSuma n =
  plotList [ Key Nothing
           , PNG "Sucesiones_de_numeros_consecutivos_con_suma_dada.png"]
           [length (sucesionesConSuma2 k) | k <- [0..n]]

import Graphics.Gnuplot.Simple

-- 1ª solución

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

sucesionesConSuma n =

[(x,y) | y <- [0..n], x <- [0..y], sum [x..y] == n]

-- 2ª solución

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

sucesionesConSuma2 n =

[(x,y) | y <- [0..n], x <- [0..y], (x+y)*(y-x+1) == 2*n]

-- Comparación de eficiencia

-- λ> sucesionesConSuma 700

-- [(3,37),(16,40),(84,91),(97,103),(138,142),(700,700)]

-- (2.71 secs, 7,894,718,264 bytes)

-- λ> sucesionesConSuma2 700

-- [(3,37),(16,40),(84,91),(97,103),(138,142),(700,700)]

-- (0.22 secs, 118,104,176 bytes)

-- Gráfica

graficaSucesionesConSuma :: Int -> IO ()

graficaSucesionesConSuma n =

plotList [ Key Nothing

, PNG "Sucesiones_de_numeros_consecutivos_con_suma_dada.png"]

[length (sucesionesConSuma2 k) | k <- [0..n]]

Números tetranacci

PorJosé A. Alonso 27-marzo-201826-marzo-2022

Los números tetranacci son una generalización de los números de Fibonacci definidos por

   T(0) = 0,
   T(1) = 1,
   T(2) = 1,
   T(3) = 2, 
   T(n) = T(n-1) + T(n-2) + T(n-3) + T(n-4), para n > 3.

T(0) = 0,

T(1) = 1,

T(2) = 1,

T(3) = 2,

T(n) = T(n-1) + T(n-2) + T(n-3) + T(n-4), para n > 3.

Los primeros números tetranacci son

   0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208

1	0, 1, 1, 2, 4, 8, 15, 29, 56, 108, 208

Definir las funciones

   tetranacci        :: Int -> Integer
   graficaTetranacci :: Int -> IO ()

1 2	tetranacci :: Int -> Integer graficaTetranacci :: Int -> IO ()

tales que

(tetranacci n) es el n-ésimo número tetranacci. Por ejemplo,

     λ> tetranacci 10
     208
     λ> map tetranacci [0..10]
     [0,1,1,2,4,8,15,29,56,108,208]
     λ> length (show (tetranacci5 (10^5)))
     28501

λ> tetranacci 10

208

λ> map tetranacci [0..10]

[0,1,1,2,4,8,15,29,56,108,208]

λ> length (show (tetranacci5 (10^5)))

28501

(graficaTetranacci n) dibuja la gráfica de los cocientes de n primeros pares de número tetranacci. Por ejemplo, (graficaTetranacci 300) dibuja

Soluciones

import Data.List (zipWith4)
import Data.Array
import Graphics.Gnuplot.Simple

-- 1ª solución (por recursión) 
-- ===========================

tetranacci :: Int -> Integer
tetranacci 0 = 0
tetranacci 1 = 1
tetranacci 2 = 1
tetranacci 3 = 2
tetranacci n =
  tetranacci (n-1) + tetranacci (n-2) + tetranacci (n-3) + tetranacci (n-4) 

-- 2ª solución (programación dinámica con zipWith4)
-- ================================================

tetranacci2 :: Int -> Integer
tetranacci2 n = tetranaccis2 !! n

tetranaccis2 :: [Integer]
tetranaccis2 = 
    0 : 1 : 1 : 2 : zipWith4 f (r 0) (r 1) (r 2) (r 3)
    where f a b c d = a+b+c+d
          r n       = drop n tetranaccis2

-- 3ª solución (con acumuladores)
-- ==============================

tetranacci3 :: Int -> Integer
tetranacci3 n = tetranaccis3 !! n

tetranaccis3 :: [Integer]
tetranaccis3 = p (0, 1, 1, 2)
    where p (a, b, c, d) = a : p (b, c, d, a + b + c + d)

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

tetranacci4 :: Int -> Integer
tetranacci4 n = tetranaccis4 !! n

tetranaccis4 :: [Integer]
tetranaccis4 = 0 : 1 : 1 : 2 : p tetranaccis4
   where p (a:b:c:d:xs) = (a+b+c+d): p (b:c:d:xs)

-- 5ª solución (programación dinámica con vectores)
-- ================================================

tetranacci5 :: Int -> Integer
tetranacci5 n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 0
  f 1 = 1
  f 2 = 1
  f 3 = 2
  f k = v!(k-1) + v!(k-2) + v!(k-3) + v!(k-4) 

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

--    λ> tetranacci 26
--    7555935
--    (3.04 secs, 1,649,520,064 bytes)
--    λ> tetranacci2 26
--    7555935
--    (0.00 secs, 148,064 bytes)
-- 
--    λ> length (show (tetranacci2 (10^5)))
--    28501
--    (1.22 secs, 1,844,457,288 bytes)
--    λ> length (show (tetranacci3 (10^5)))
--    28501
--    (0.88 secs, 1,860,453,968 bytes)
--    λ> length (show (tetranacci4 (10^5)))
--    28501
--    (0.77 secs, 1,882,852,168 bytes)
--    λ> length (show (tetranacci5 (10^5)))
--    28501
--    (0.72 secs, 1,905,707,408 bytes)

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

graficaTetranacci :: Int -> IO ()
graficaTetranacci n =
  plotList [ Key Nothing
           , Title "Tasa de crecimiento de los numeros tetranacci"
           , PNG ("Numeros_tetranacci_" ++ show n ++ ".png")
           ]
           (take n (zipWith (/) (tail xs) xs))
  where xs = (map fromIntegral tetranaccis4) :: [Double]

import Data.List (zipWith4)

import Data.Array

import Graphics.Gnuplot.Simple

-- 1ª solución (por recursión)

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

tetranacci :: Int -> Integer

tetranacci 0 = 0

tetranacci 1 = 1

tetranacci 2 = 1

tetranacci 3 = 2

tetranacci n =

tetranacci (n-1) + tetranacci (n-2) + tetranacci (n-3) + tetranacci (n-4)

-- 2ª solución (programación dinámica con zipWith4)

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

tetranacci2 :: Int -> Integer

tetranacci2 n = tetranaccis2 !! n

tetranaccis2 :: [Integer]

tetranaccis2 =

0 : 1 : 1 : 2 : zipWith4 f (r 0) (r 1) (r 2) (r 3)

where f a b c d = a+b+c+d

r n = drop n tetranaccis2

-- 3ª solución (con acumuladores)

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

tetranacci3 :: Int -> Integer

tetranacci3 n = tetranaccis3 !! n

tetranaccis3 :: [Integer]

tetranaccis3 = p (0, 1, 1, 2)

where p (a, b, c, d) = a : p (b, c, d, a + b + c + d)

-- 4ª solución

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

tetranacci4 :: Int -> Integer

tetranacci4 n = tetranaccis4 !! n

tetranaccis4 :: [Integer]

tetranaccis4 = 0 : 1 : 1 : 2 : p tetranaccis4

where p (a:b:c:d:xs) = (a+b+c+d): p (b:c:d:xs)

-- 5ª solución (programación dinámica con vectores)

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

tetranacci5 :: Int -> Integer

tetranacci5 n = v ! n where

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

f 0 = 0

f 1 = 1

f 2 = 1

f 3 = 2

f k = v!(k-1) + v!(k-2) + v!(k-3) + v!(k-4)

-- Comparación de eficiencia

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

-- λ> tetranacci 26

-- 7555935

-- (3.04 secs, 1,649,520,064 bytes)

-- λ> tetranacci2 26

-- 7555935

-- (0.00 secs, 148,064 bytes)

-- λ> length (show (tetranacci2 (10^5)))

-- 28501

-- (1.22 secs, 1,844,457,288 bytes)

-- λ> length (show (tetranacci3 (10^5)))

-- 28501

-- (0.88 secs, 1,860,453,968 bytes)

-- λ> length (show (tetranacci4 (10^5)))

-- 28501

-- (0.77 secs, 1,882,852,168 bytes)

-- λ> length (show (tetranacci5 (10^5)))

-- 28501

-- (0.72 secs, 1,905,707,408 bytes)

-- Gráfica

-- =======

graficaTetranacci :: Int -> IO ()

graficaTetranacci n =

plotList [ Key Nothing

, Title "Tasa de crecimiento de los numeros tetranacci"

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

]

(take n (zipWith (/) (tail xs) xs))

where xs = (map fromIntegral tetranaccis4) :: [Double]

Conjetura de Goldbach

PorJosé A. Alonso 21-marzo-201826-marzo-2022

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

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

Definir las siguientes funciones

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

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

tales que

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

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

indiceGoldbach 2 == 1

indiceGoldbach 4 == 2

indiceGoldbach 27 == 3

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

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

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

Comprobar con QuickCheck la conjetura de Goldbach anterior.

Soluciones

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

100

101

102

103

104

105

106

107

import Data.Array

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

import Test.QuickCheck

-- 1ª definición

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

indiceGoldbach :: Int -> Int

indiceGoldbach n =

minimum (map length (particiones n))

particiones :: Int -> [[Int]]

particiones n = v ! n where

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

where f 0 = [[]]

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

y <- v ! (m-x),

[x] >= take 1 y]

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

-- 2ª definición

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

indiceGoldbach2 :: Int -> Int

indiceGoldbach2 x =

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

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

-- primos. Por ejemplo,

-- esSumaDe 2 1 == True

-- esSumaDe 4 1 == False

-- esSumaDe 4 2 == True

-- esSumaDe 27 2 == False

-- esSumaDe 27 3 == True

esSumaDe :: Int -> Int -> Bool

esSumaDe x 1 = isPrime x

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

-- 3ª definición

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

indiceGoldbach3 :: Int -> Int

indiceGoldbach3 x =

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

esSumaDe3 :: Int -> Int -> Bool

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

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

f i 1 = isPrime i

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

-- 4ª definición

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

indiceGoldbach4 :: Int -> Int

indiceGoldbach4 n = v ! n where

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

f i | isPrime i = 1

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

-- Comparación de eficiencia

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

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

-- 142

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

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

-- 142

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

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

-- 142

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

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

-- 142

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

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

-- 2030

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

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

-- 2030

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

-- Gráfica

-- =======

graficaGoldbach :: Int -> IO ()

graficaGoldbach n =

plotList [ Key Nothing

, XRange (2,fromIntegral n)

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

]

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

-- Comprobación de la conjetura de Goldbach

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

-- La propiedad es

prop_Goldbach :: Int -> Property

prop_Goldbach x =

x >= 2 ==> indiceGoldbach2 x < 4

-- La comprobación es

-- λ> quickCheck prop_Goldbach

-- +++ OK, passed 100 tests.

La sucesión de Sylvester

PorJosé A. Alonso 19-marzo-201825-abril-2021

La sucesión de Sylvester es la sucesión que comienza en 2 y sus restantes términos se obtienen multiplicando los anteriores y sumándole 1.

Definir las funciones

   sylvester        :: Integer -> Integer
   graficaSylvester :: Integer -> Integer -> IO ()

1 2	sylvester :: Integer -> Integer graficaSylvester :: Integer -> Integer -> IO ()

tales que

(sylvester n) es el n-ésimo término de la sucesión de Sylvester. Por ejemplo,

     λ> [sylvester n | n <- [0..7]]
     [2,3,7,43,1807,3263443,10650056950807,113423713055421844361000443]
     λ> length (show (sylvester 25))
     6830085

λ> [sylvester n | n <- [0..7]]

[2,3,7,43,1807,3263443,10650056950807,113423713055421844361000443]

λ> length (show (sylvester 25))

6830085

(graficaSylvester d n) dibuja la gráfica de los d últimos dígitos de los n primeros términos de la sucesión de Sylvester. Por ejemplo,
- (graficaSylvester 3 30) dibuja
- (graficaSylvester 4 30) dibuja
- (graficaSylvester 5 30) dibuja

Soluciones

import Data.List               (genericIndex)
import Data.Array              ((!), array)
import Graphics.Gnuplot.Simple (plotList, Attribute (Key, PNG))

-- 1ª solución (por recursión)
-- ===========================

sylvester1 :: Integer -> Integer
sylvester1 0 = 2
sylvester1 n = 1 + product [sylvester1 k | k <- [0..n-1]]

-- 2ª solución (con programación dinámica)
-- =======================================

sylvester2 :: Integer -> Integer
sylvester2 n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 2
  f m = 1 + product [v!k | k <- [0..m-1]]

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

sylvester3 :: Integer -> Integer
sylvester3 0 = 2
sylvester3 n = 1 + x^2 - x
  where x = sylvester3 (n-1)

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

sylvester4 :: Integer -> Integer
sylvester4 n = v ! n where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f 0 = 2
  f m = 1 + x^2 - x
    where x = v ! (m-1)

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

sylvester5 :: Integer -> Integer
sylvester5 0 = 2
sylvester5 n = 1 + (productosSylvester `genericIndex` (n-1))

sucSylvester5 :: [Integer]
sucSylvester5 = map sylvester5 [0..]

productosSylvester :: [Integer]
productosSylvester = scanl1 (*) sucSylvester5

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

--    λ> length (show (sylvester1 20))
--    213441
--    (3.40 secs, 519,249,840 bytes)
--    λ> length (show (sylvester2 20))
--    213441
--    (0.10 secs, 13,716,024 bytes)
--    λ> length (show (sylvester3 20))
--    213441
--    (0.16 secs, 13,646,472 bytes)
--    λ> length (show (sylvester4 20))
--    213441
--    (0.18 secs, 13,654,064 bytes)
--    λ> length (show (sylvester5 20))
--    213441
--    (0.12 secs, 13,497,192 bytes)

graficaSylvester :: Integer -> Integer -> IO ()
graficaSylvester d n =
  plotList [ Key Nothing
           , PNG ("La_sucesion_de_Sylvester_" ++ show (d,n) ++ ".png")
           ]
           [sylvester5 k `mod` (10^d) | k <- [0..n]]

import Data.List (genericIndex)

import Data.Array ((!), array)

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

-- 1ª solución (por recursión)

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

sylvester1 :: Integer -> Integer

sylvester1 0 = 2

sylvester1 n = 1 + product [sylvester1 k | k <- [0..n-1]]

-- 2ª solución (con programación dinámica)

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

sylvester2 :: Integer -> Integer

sylvester2 n = v ! n where

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

f 0 = 2

f m = 1 + product [v!k | k <- [0..m-1]]

-- 3ª solución

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

sylvester3 :: Integer -> Integer

sylvester3 0 = 2

sylvester3 n = 1 + x^2 - x

where x = sylvester3 (n-1)

-- 4ª solución

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

sylvester4 :: Integer -> Integer

sylvester4 n = v ! n where

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

f 0 = 2

f m = 1 + x^2 - x

where x = v ! (m-1)

-- 4ª solución

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

sylvester5 :: Integer -> Integer

sylvester5 0 = 2

sylvester5 n = 1 + (productosSylvester `genericIndex` (n-1))

sucSylvester5 :: [Integer]

sucSylvester5 = map sylvester5 [0..]

productosSylvester :: [Integer]

productosSylvester = scanl1 (*) sucSylvester5

-- Comparación de eficiencia

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

-- λ> length (show (sylvester1 20))

-- 213441

-- (3.40 secs, 519,249,840 bytes)

-- λ> length (show (sylvester2 20))

-- 213441

-- (0.10 secs, 13,716,024 bytes)

-- λ> length (show (sylvester3 20))

-- 213441

-- (0.16 secs, 13,646,472 bytes)

-- λ> length (show (sylvester4 20))

-- 213441

-- (0.18 secs, 13,654,064 bytes)

-- λ> length (show (sylvester5 20))

-- 213441

-- (0.12 secs, 13,497,192 bytes)

graficaSylvester :: Integer -> Integer -> IO ()

graficaSylvester d n =

plotList [ Key Nothing

, PNG ("La_sucesion_de_Sylvester_" ++ show (d,n) ++ ".png")

]

[sylvester5 k `mod` (10^d) | k <- [0..n]]

La conjetura de Levy

PorJosé A. Alonso 12-marzo-201825-abril-2021

Hyman Levy observó que

    7 = 3 + 2 x 2
    9 = 3 + 2 x 3 =  5 + 2 x 2
   11 = 5 + 2 x 3 =  7 + 2 x 2
   13 = 3 + 2 x 5 =  7 + 2 x 3
   15 = 3 + 2 x 5 = 11 + 2 x 2
   17 = 3 + 2 x 7 =  7 + 2 x 5 = 11 + 2 x 3 = 13 + 2 x 2
   19 = 5 + 2 x 7 = 13 + 2 x 3

7 = 3 + 2 x 2

9 = 3 + 2 x 3 = 5 + 2 x 2

11 = 5 + 2 x 3 = 7 + 2 x 2

13 = 3 + 2 x 5 = 7 + 2 x 3

15 = 3 + 2 x 5 = 11 + 2 x 2

17 = 3 + 2 x 7 = 7 + 2 x 5 = 11 + 2 x 3 = 13 + 2 x 2

19 = 5 + 2 x 7 = 13 + 2 x 3

y conjeturó que todos los número impares mayores o iguales que 7 se pueden escribir como la suma de un primo y el doble de un primo. El objetivo de los siguientes ejercicios es comprobar la conjetura de Levy.

Definir las siguientes funciones

   descomposicionesLevy :: Integer -> [(Integer,Integer)]
   graficaLevy          :: Integer -> IO ()

1 2	descomposicionesLevy :: Integer -> [(Integer,Integer)] graficaLevy :: Integer -> IO ()

tales que

(descomposicionesLevy x) es la lista de pares de primos (p,q) tales que x = p + 2q. Por ejemplo,

     descomposicionesLevy  7  ==  [(3,2)]
     descomposicionesLevy  9  ==  [(3,3),(5,2)]
     descomposicionesLevy 17  ==  [(3,7),(7,5),(11,3),(13,2)]

descomposicionesLevy 7 == [(3,2)]

descomposicionesLevy 9 == [(3,3),(5,2)]

descomposicionesLevy 17 == [(3,7),(7,5),(11,3),(13,2)]

(graficaLevy n) dibuja los puntos (x,y) tales que x pertenece a [7,9..7+2x(n-1)] e y es el número de descomposiciones de Levy de x. Por ejemplo, (graficaLevy 200) dibuja

Comprobar con QuickCheck la conjetura de Levy.

Soluciones

[schedule on=’2018-03-19′ at=»06:00″]

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

descomposicionesLevy :: Integer -> [(Integer,Integer)]
descomposicionesLevy x =
  [(p,q) | p <- takeWhile (< x) (tail primes)
         , let q = (x - p) `div` 2
         , isPrime q]

-- graficaLevy 300
graficaLevy :: Integer -> IO ()
graficaLevy n =
  plotList [ Key Nothing
           , XRange (7,fromIntegral (7+2*(n-1)))
           , PNG ("La_conjetura_de_Levy-" ++ show n ++ ".png")
           ]
           [(x, length (descomposicionesLevy x)) | x <- [7,9..7+2*(n-1)]] 

-- La propiedad es
prop_Levy :: Integer -> Bool
prop_Levy x =
  not (null (descomposicionesLevy (7 + 2 * abs x)))

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

import Data.Numbers.Primes

import Test.QuickCheck

import Graphics.Gnuplot.Simple

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

descomposicionesLevy x =

[(p,q) | p <- takeWhile (< x) (tail primes)

, let q = (x - p) `div` 2

, isPrime q]

-- graficaLevy 300

graficaLevy :: Integer -> IO ()

graficaLevy n =

plotList [ Key Nothing

, XRange (7,fromIntegral (7+2*(n-1)))

, PNG ("La_conjetura_de_Levy-" ++ show n ++ ".png")

]

[(x, length (descomposicionesLevy x)) | x <- [7,9..7+2*(n-1)]]

-- La propiedad es

prop_Levy :: Integer -> Bool

prop_Levy x =

not (null (descomposicionesLevy (7 + 2 * abs x)))

-- La comprobación es

-- λ> quickCheck prop_Levy

-- +++ OK, passed 100 tests.

La función de Smarandache

PorJosé A. Alonso 27-febrero-20186-marzo-2018

La función de Smarandache, también conocida como la función de Kempner, es la función que asigna a cada número entero positivo n el menor número cuyo factorial es divisible por n y se representa por S(n). Por ejemplo, el número 8 no divide a 1!, 2!, 3!, pero sí divide 4!; por tanto, S(8) = 4.

Definir las funciones

   smarandache        :: Integer -> Integer
   graficaSmarandache :: Integer -> IO ()

1 2	smarandache :: Integer -> Integer graficaSmarandache :: Integer -> IO ()

tales que

(smarandache n) es el menor número cuyo factorial es divisible por n. Por ejemplo,

     smarandache 8   ==  4
     smarandache 10  ==  5
     smarandache 16  ==  6

smarandache 8 == 4

smarandache 10 == 5

smarandache 16 == 6

(graficaSmarandache n) dibuja la gráfica de los n primeros términos de la sucesión de Smarandache. Por ejemplo, (graficaSmarandache 100) dibuja

(graficaSmarandache 500) dibuja

Soluciones

import Data.List (genericLength)
import Graphics.Gnuplot.Simple

smarandache :: Integer -> Integer
smarandache x =
  head [n | (n,y) <- zip [0..] factoriales
          , y `mod` x == 0]

-- factoriales es la lista de los factoriales. Por ejemplo, 
--    λ> take 12 factoriales
--    [1,1,2,6,24,120,720,5040,40320,362880,3628800,39916800]
factoriales :: [Integer]
factoriales = 1 : scanl1 (*) [1..]

graficaSmarandache :: Integer -> IO ()
graficaSmarandache n =
  plotList [Key Nothing
           , PNG ("La_funcion_de_Smarandache_" ++ show n ++ ".png")
           ]
           (map smarandache [1..n])

import Data.List (genericLength)

import Graphics.Gnuplot.Simple

smarandache :: Integer -> Integer

smarandache x =

head [n | (n,y) <- zip [0..] factoriales

, y `mod` x == 0]

-- factoriales es la lista de los factoriales. Por ejemplo,

-- λ> take 12 factoriales

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

factoriales :: [Integer]

factoriales = 1 : scanl1 (*) [1..]

graficaSmarandache :: Integer -> IO ()

graficaSmarandache n =

plotList [Key Nothing

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

]

(map smarandache [1..n])

Menor potencia de 2 que comienza por n

PorJosé A. Alonso 22-febrero-201826-marzo-2022

Definir las funciones

   menorPotencia            :: Integer -> (Integer,Integer)
   graficaMenoresExponentes :: Integer -> IO ()

1 2	menorPotencia :: Integer -> (Integer,Integer) graficaMenoresExponentes :: Integer -> IO ()

tales que

(menorPotencia n) es el par (k,m) donde m es la menor potencia de 2 que empieza por n y k es su exponentes (es decir, 2^k = m). Por ejemplo,

     menorPotencia 3             ==  (5,32)
     menorPotencia 7             ==  (46,70368744177664)
     fst (menorPotencia 982)     ==  3973
     fst (menorPotencia 32627)   ==  28557
     fst (menorPotencia 158426)  ==  40000

menorPotencia 3 == (5,32)

menorPotencia 7 == (46,70368744177664)

fst (menorPotencia 982) == 3973

fst (menorPotencia 32627) == 28557

fst (menorPotencia 158426) == 40000

(graficaMenoresExponentes n) dibuja la gráfica de los exponentes de 2 en las menores potencias de los n primeros números enteros positivos. Por ejemplo, (graficaMenoresExponentes 200) dibuja

Soluciones

import Data.List               (isPrefixOf)
import Graphics.Gnuplot.Simple (Attribute (Key, PNG), plotList)

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

menorPotencia :: Integer -> (Integer,Integer)
menorPotencia n =
  head [(k,m) | (k,m) <- zip [0..] potenciasDe2
              , cs `isPrefixOf` show m]
  where cs = show n

-- potenciasDe 2 es la lista de las potencias de dos. Por ejemplo,
--    take 12 potenciasDe2  ==  [1,2,4,8,16,32,64,128,256,512,1024,2048]
potenciasDe2 :: [Integer]
potenciasDe2 = iterate (*2) 1

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

menorPotencia2 :: Integer -> (Integer,Integer)
menorPotencia2 n = aux (0,1)
  where aux (k,m) | cs `isPrefixOf` show m = (k,m)
                  | otherwise              = aux (k+1,2*m)
        cs = show n

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

menorPotencia3 :: Integer -> (Integer,Integer)
menorPotencia3 n =
  until (isPrefixOf n1 . show . snd) (\(x,y) -> (x+1,2*y)) (0,1)
  where n1 = show n

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

--    λ> maximum [fst (menorPotencia n) | n <- [1..1000]]
--    3973
--    (3.69 secs, 1,094,923,696 bytes)
--    λ> maximum [fst (menorPotencia2 n) | n <- [1..1000]]
--    3973
--    (5.13 secs, 1,326,382,872 bytes)
--    λ> maximum [fst (menorPotencia3 n) | n <- [1..1000]]
--    3973
--    (4.71 secs, 1,240,498,128 bytes)

graficaMenoresExponentes :: Integer -> IO ()
graficaMenoresExponentes n =
  plotList [ Key Nothing
           , PNG "Menor_potencia_de_2_que_comienza_por_n.png"
           ]
           (map (fst . menorPotencia) [1..n])

import Data.List (isPrefixOf)

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

-- 1ª definición

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

menorPotencia :: Integer -> (Integer,Integer)

menorPotencia n =

head [(k,m) | (k,m) <- zip [0..] potenciasDe2

, cs `isPrefixOf` show m]

where cs = show n

-- potenciasDe 2 es la lista de las potencias de dos. Por ejemplo,

-- take 12 potenciasDe2 == [1,2,4,8,16,32,64,128,256,512,1024,2048]

potenciasDe2 :: [Integer]

potenciasDe2 = iterate (*2) 1

-- 2ª definición

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

menorPotencia2 :: Integer -> (Integer,Integer)

menorPotencia2 n = aux (0,1)

where aux (k,m) | cs `isPrefixOf` show m = (k,m)

| otherwise = aux (k+1,2*m)

cs = show n

-- 3ª definición

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

menorPotencia3 :: Integer -> (Integer,Integer)

menorPotencia3 n =

until (isPrefixOf n1 . show . snd) (\(x,y) -> (x+1,2*y)) (0,1)

where n1 = show n

-- Comparación de eficiencia

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

-- λ> maximum [fst (menorPotencia n) | n <- [1..1000]]

-- 3973

-- (3.69 secs, 1,094,923,696 bytes)

-- λ> maximum [fst (menorPotencia2 n) | n <- [1..1000]]

-- 3973

-- (5.13 secs, 1,326,382,872 bytes)

-- λ> maximum [fst (menorPotencia3 n) | n <- [1..1000]]

-- 3973

-- (4.71 secs, 1,240,498,128 bytes)

graficaMenoresExponentes :: Integer -> IO ()

graficaMenoresExponentes n =

plotList [ Key Nothing

, PNG "Menor_potencia_de_2_que_comienza_por_n.png"

]

(map (fst . menorPotencia) [1..n])

Períodos de Fibonacci

PorJosé A. Alonso 9-febrero-201819-febrero-2018

Los primeros términos de la sucesión de Fibonacci son

   0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610

1	0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610

Al calcular sus restos módulo 3 se obtiene

   0,1,1,2,0,2,2,1, 0,1,1,2,0,2,2,1

1	0,1,1,2,0,2,2,1, 0,1,1,2,0,2,2,1

Se observa que es periódica y su período es

   0,1,1,2,0,2,2,1

1	0,1,1,2,0,2,2,1

Definir las funciones

   fibsMod                   :: Integer -> [Integer]
   periodoFibMod             :: Integer -> [Integer]
   longPeriodosFibMod        :: [Int]
   graficaLongPeriodosFibMod :: Int -> IO ()

fibsMod :: Integer -> [Integer]

periodoFibMod :: Integer -> [Integer]

longPeriodosFibMod :: [Int]

graficaLongPeriodosFibMod :: Int -> IO ()

tales que

(fibsMod n) es la lista de los términos de la sucesión de Fibonacci módulo n. Por ejemplo,

     λ> take 24 (fibsMod 3)
     [0,1,1,2,0,2,2,1, 0,1,1,2,0,2,2,1, 0,1,1,2,0,2,2,1]
     λ> take 24 (fibsMod 4)
     [0,1,1,2,3,1, 0,1,1,2,3,1, 0,1,1,2,3,1, 0,1,1,2,3,1]

λ> take 24 (fibsMod 3)

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

λ> take 24 (fibsMod 4)

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

(periodoFibMod n) es la parte perioica de la sucesión de Fibonacci módulo n. Por ejemplo,

     periodoFibMod 3  ==  [0,1,1,2,0,2,2,1]
     periodoFibMod 4  ==  [0,1,1,2,3,1]
     periodoFibMod 7  ==  [0,1,1,2,3,5,1,6,0,6,6,5,4,2,6,1]

periodoFibMod 3 == [0,1,1,2,0,2,2,1]

periodoFibMod 4 == [0,1,1,2,3,1]

periodoFibMod 7 == [0,1,1,2,3,5,1,6,0,6,6,5,4,2,6,1]

longPeriodosFibMod es la sucesión de las longitudes de los períodos de las sucesiones de Fibonacci módulo n, para n > 0. Por ejemplo,

     λ> take 20 longPeriodosFibMod
     [1,3,8,6,20,24,16,12,24,60,10,24,28,48,40,24,36,24,18,60]

1 2	λ> take 20 longPeriodosFibMod [1,3,8,6,20,24,16,12,24,60,10,24,28,48,40,24,36,24,18,60]

(graficaLongPeriodosFibMod n) dibuja la gráfica de los n primeros términos de la sucesión longPeriodosFibMod. Por ejemplo, (graficaLongPeriodosFibMod n) dibuja

Soluciones

import Graphics.Gnuplot.Simple

fibsMod :: Integer -> [Integer]
fibsMod n = map (`mod` n) fibs

-- fibs es la sucesión de Fibonacci. Por ejemplo,
--    λ> take 20 fibs
--    [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]
fibs :: [Integer]
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

periodoFibMod :: Integer -> [Integer]
periodoFibMod 1 = [0]
periodoFibMod n = 0 : 1 : aux (drop 2 (fibsMod n))
  where aux (0:1:xs) = []
        aux (a:b:xs) = a : aux (b:xs)

longPeriodosFibMod :: [Int]
longPeriodosFibMod =
  [length (periodoFibMod n) | n <- [1..]]

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

longPeriodosFibMod2 :: [Int]
longPeriodosFibMod2 = map longPeriodoFibMod [1..]

longPeriodoFibMod :: Integer -> Int
longPeriodoFibMod 1 = 1
longPeriodoFibMod n = aux 1 (tail (fibsMod n)) 0
  where aux 0 (1 : xs) k = k
        aux _ (x : xs) k = aux x xs (k + 1)

graficaLongPeriodosFibMod :: Int -> IO ()
graficaLongPeriodosFibMod n =
  plotList [ Key Nothing
           , Title ("graficaLongPeriodosFibMod " ++ show n)
           , PNG ("Periodos_de_Fibonacci " ++ show n ++ ".png")]
           (take n longPeriodosFibMod)

import Graphics.Gnuplot.Simple

fibsMod :: Integer -> [Integer]

fibsMod n = map (`mod` n) fibs

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

-- λ> take 20 fibs

-- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181]

fibs :: [Integer]

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

periodoFibMod :: Integer -> [Integer]

periodoFibMod 1 = [0]

periodoFibMod n = 0 : 1 : aux (drop 2 (fibsMod n))

where aux (0:1:xs) = []

aux (a:b:xs) = a : aux (b:xs)

longPeriodosFibMod :: [Int]

longPeriodosFibMod =

[length (periodoFibMod n) | n <- [1..]]

-- 2ª definición de longPeriodosFibMod

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

longPeriodosFibMod2 :: [Int]

longPeriodosFibMod2 = map longPeriodoFibMod [1..]

longPeriodoFibMod :: Integer -> Int

longPeriodoFibMod 1 = 1

longPeriodoFibMod n = aux 1 (tail (fibsMod n)) 0

where aux 0 (1 : xs) k = k

aux _ (x : xs) k = aux x xs (k + 1)

graficaLongPeriodosFibMod :: Int -> IO ()

graficaLongPeriodosFibMod n =

plotList [ Key Nothing

, Title ("graficaLongPeriodosFibMod " ++ show n)

, PNG ("Periodos_de_Fibonacci " ++ show n ++ ".png")]

(take n longPeriodosFibMod)

Huecos binarios

PorJosé A. Alonso 7-febrero-201814-febrero-2018

Los huecos binarios de un número natural n son las listas de cer0 entre dos unos en la representación binaria de n. Por ejemplo, puesto que la representación binaria de 20 es 10100 tiene dos huecos binarios de longitudes 1 y 2. La longitud del mayor hueco binario de 529 es 4 ya que la representación binaria de 529 es 1000010001.

Definir las funciones

   longMayorHuecoBinario        :: Int -> Int
   graficaLongMayorHuecoBinario :: Int -> IO ()

1 2	longMayorHuecoBinario :: Int -> Int graficaLongMayorHuecoBinario :: Int -> IO ()

tales que

(longMayorHuecoBinario n) es la longitud del mayor hueco binario de n. Por ejemplo,

     longMayorHuecoBinario 20    ==  2
     longMayorHuecoBinario 529   ==  4
     longMayorHuecoBinario 2018  ==  3

longMayorHuecoBinario 20 == 2

longMayorHuecoBinario 529 == 4

longMayorHuecoBinario 2018 == 3

(graficaLongMayorHuecoBinario n) dibuja la gráfica de las longitudes de los mayores huecos binarios de los n primeros números naturales. Por ejemplo, (graficaLongMayorHuecoBinario 200) dibuja

Soluciones

import Data.List (group)
import Graphics.Gnuplot.Simple

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

longMayorHuecoBinario :: Int -> Int
longMayorHuecoBinario n =
  maximum (0 : map length (huecosBinarios (decimalAbinario n)))

-- (decimalAbinario x) es la representación binaria del número c. Por
-- ejemplo, 
--    decimalAbinario 20   ==  [0,0,1,0,1]
--    decimalAbinario 529  ==  [1,0,0,0,1,0,0,0,0,1] 
decimalAbinario :: Int -> [Int]
decimalAbinario x
  | x < 2     = [x]
  | otherwise = x `mod` 2 : decimalAbinario (x `div` 2)

-- (huecosBinarios xs) es la lista de los ceros consecutivos de xs entre
-- dos elementos distintos de cero. Por ejemplo,
--    huecosBinarios [0,0,1,0,1]            ==  [[0,0],[0]]
--    huecosBinarios [1,0,0,0,1,0,0,0,0,1]  ==  [[0,0,0],[0,0,0,0]]
huecosBinarios :: [Int] -> [[Int]]
huecosBinarios xs =
  [ys | ys <- group xs, 0 `elem` ys]

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

longMayorHuecoBinario2 :: Int -> Int
longMayorHuecoBinario2 =
  maximum . (0 :) . map length . huecosBinarios . decimalAbinario

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

longMayorHuecoBinario3 :: Int -> Int
longMayorHuecoBinario3 = maximum . longHuecosBinarios

-- (longHuecosBinarios n) es la lista de las longitudes de los huecos
-- binarios de n. Por ejemplo,
--    longHuecosBinarios 20  ==  [2,1,0]
--    longHuecosBinarios 529  ==  [3,4,0]
longHuecosBinarios :: Int -> [Int]
longHuecosBinarios 1 = [0]
longHuecosBinarios n | mod n 2 == 1 = longHuecosBinarios (div n 2)
                     | mod n 4 == 2 = 1 : h : t
                     | otherwise = 1 + h : t
  where (h : t) = longHuecosBinarios (div n 2)

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

--    λ> maximum [longMayorHuecoBinario x | x <- [1..100000]]
--    16
--    (5.51 secs, 941,090,272 bytes)
--    λ> maximum [longMayorHuecoBinario2 x | x <- [1..100000]]
--    16
--    (5.54 secs, 933,093,168 bytes)
--    λ> maximum [longMayorHuecoBinario3 x | x <- [1..100000]]
--    16
--    (6.00 secs, 966,584,848 bytes)

graficaLongMayorHuecoBinario :: Int -> IO ()
graficaLongMayorHuecoBinario n =
  plotList [ Key Nothing
           , Title ("graficaLongMayorHuecoBinario " ++ show n)
           , PNG ("Huecos_binarios_" ++ show n ++ ".png")
           ]
           [longMayorHuecoBinario k | k <- [0..n-1]]

import Data.List (group)

import Graphics.Gnuplot.Simple

-- 1ª solución

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

longMayorHuecoBinario :: Int -> Int

longMayorHuecoBinario n =

maximum (0 : map length (huecosBinarios (decimalAbinario n)))

-- (decimalAbinario x) es la representación binaria del número c. Por

-- ejemplo,

-- decimalAbinario 20 == [0,0,1,0,1]

-- decimalAbinario 529 == [1,0,0,0,1,0,0,0,0,1]

decimalAbinario :: Int -> [Int]

decimalAbinario x

| x < 2 = [x]

| otherwise = x `mod` 2 : decimalAbinario (x `div` 2)

-- (huecosBinarios xs) es la lista de los ceros consecutivos de xs entre

-- dos elementos distintos de cero. Por ejemplo,

-- huecosBinarios [0,0,1,0,1] == [[0,0],[0]]

-- huecosBinarios [1,0,0,0,1,0,0,0,0,1] == [[0,0,0],[0,0,0,0]]

huecosBinarios :: [Int] -> [[Int]]

huecosBinarios xs =

[ys | ys <- group xs, 0 `elem` ys]

-- 2ª solución

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

longMayorHuecoBinario2 :: Int -> Int

longMayorHuecoBinario2 =

maximum . (0 :) . map length . huecosBinarios . decimalAbinario

-- 3ª solución

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

longMayorHuecoBinario3 :: Int -> Int

longMayorHuecoBinario3 = maximum . longHuecosBinarios

-- (longHuecosBinarios n) es la lista de las longitudes de los huecos

-- binarios de n. Por ejemplo,

-- longHuecosBinarios 20 == [2,1,0]

-- longHuecosBinarios 529 == [3,4,0]

longHuecosBinarios :: Int -> [Int]

longHuecosBinarios 1 = [0]

longHuecosBinarios n | mod n 2 == 1 = longHuecosBinarios (div n 2)

| mod n 4 == 2 = 1 : h : t

| otherwise = 1 + h : t

where (h : t) = longHuecosBinarios (div n 2)

-- Comparación de eficiencia

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

-- λ> maximum [longMayorHuecoBinario x | x <- [1..100000]]

-- 16

-- (5.51 secs, 941,090,272 bytes)

-- λ> maximum [longMayorHuecoBinario2 x | x <- [1..100000]]

-- 16

-- (5.54 secs, 933,093,168 bytes)

-- λ> maximum [longMayorHuecoBinario3 x | x <- [1..100000]]

-- 16

-- (6.00 secs, 966,584,848 bytes)

graficaLongMayorHuecoBinario :: Int -> IO ()

graficaLongMayorHuecoBinario n =

plotList [ Key Nothing

, Title ("graficaLongMayorHuecoBinario " ++ show n)

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

]

[longMayorHuecoBinario k | k <- [0..n-1]]

Dígitos iniciales

PorJosé A. Alonso 2-febrero-201819-febrero-2018

Definir las funciones

   digitosIniciales        :: [Int]
   graficaDigitosIniciales :: Int -> IO ()

1 2	digitosIniciales :: [Int] graficaDigitosIniciales :: Int -> IO ()

tales que

digitosIniciales es la lista de los dígitos iniciales de los números naturales. Por ejemplo,

     λ> take 100 digitosIniciales
     [0,1,2,3,4,5,6,7,8,9,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,
      3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,
      6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,
      9,9,9,9,9,9,9,9,9,9]

λ> take 100 digitosIniciales

[0,1,2,3,4,5,6,7,8,9,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,

3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,

6,6,6,6,6,6,6,6,6,6,7,7,7,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,

9,9,9,9,9,9,9,9,9,9]

(graficaDigitosIniciales n) dibuja la gráfica de los primeros n términos de la sucesión digitosIniciales. Por ejemplo, (graficaDigitosIniciales 100) dibuja

y (graficaDigitosIniciales 1000) dibuja

Soluciones

import Graphics.Gnuplot.Simple

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

digitosIniciales :: [Int]
digitosIniciales = map digitoInicial [0..]

digitoInicial :: Integer -> Int
digitoInicial n = read [head (show n)]

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

digitosIniciales2 :: [Int]
digitosIniciales2 = map (read . return . head . show) [0..]

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

digitosIniciales3 :: [Int]
digitosIniciales3 = map digitoInicial3 [0..]

digitoInicial3 :: Integer -> Int
digitoInicial3 = fromInteger . until (< 10) (`div` 10)

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

digitosIniciales4 :: [Int]
digitosIniciales4 = map (fromInteger . until (< 10) (`div` 10)) [0..]

-- 5ª definición
-- =============

digitosIniciales5 :: [Int]
digitosIniciales5 =
  0 : concat [replicate k x | k <- [10^n | n <- [0..]]
                            , x <- [1..9]]

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

--    λ> digitosIniciales !! (2*10^6)
--    2
--    (0.46 secs, 320,145,984 bytes)
--    λ> digitosIniciales2 !! (2*10^6)
--    2
--    (0.46 secs, 320,143,288 bytes)
--    λ> digitosIniciales3 !! (2*10^6)
--    2
--    (0.17 secs, 320,139,216 bytes)
--    λ> digitosIniciales4 !! (2*10^6)
--    2
--    (0.55 secs, 320,139,248 bytes)
--    λ> digitosIniciales5 !! (2*10^6)
--    2
--    (0.12 secs, 224,158,992 bytes)

graficaDigitosIniciales :: Int -> IO ()
graficaDigitosIniciales n =
  plotList [ Key Nothing
           , Title ("graficaDigitosIniciales " ++ show n)
           , PNG ("Digitos_iniciales_" ++ show n ++ ".png" )
           ]
           (take n digitosIniciales)

import Graphics.Gnuplot.Simple

-- 1ª definición

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

digitosIniciales :: [Int]

digitosIniciales = map digitoInicial [0..]

digitoInicial :: Integer -> Int

digitoInicial n = read [head (show n)]

-- 2ª definición

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

digitosIniciales2 :: [Int]

digitosIniciales2 = map (read . return . head . show) [0..]

-- 3ª definición

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

digitosIniciales3 :: [Int]

digitosIniciales3 = map digitoInicial3 [0..]

digitoInicial3 :: Integer -> Int

digitoInicial3 = fromInteger . until (< 10) (`div` 10)

-- 4ª definición

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

digitosIniciales4 :: [Int]

digitosIniciales4 = map (fromInteger . until (< 10) (`div` 10)) [0..]

-- 5ª definición

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

digitosIniciales5 :: [Int]

digitosIniciales5 =

0 : concat [replicate k x | k <- [10^n | n <- [0..]]

, x <- [1..9]]

-- Comparación de eficiencia

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

-- λ> digitosIniciales !! (2*10^6)

-- 2

-- (0.46 secs, 320,145,984 bytes)

-- λ> digitosIniciales2 !! (2*10^6)

-- 2

-- (0.46 secs, 320,143,288 bytes)

-- λ> digitosIniciales3 !! (2*10^6)

-- 2

-- (0.17 secs, 320,139,216 bytes)

-- λ> digitosIniciales4 !! (2*10^6)

-- 2

-- (0.55 secs, 320,139,248 bytes)

-- λ> digitosIniciales5 !! (2*10^6)

-- 2

-- (0.12 secs, 224,158,992 bytes)

graficaDigitosIniciales :: Int -> IO ()

graficaDigitosIniciales n =

plotList [ Key Nothing

, Title ("graficaDigitosIniciales " ++ show n)

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

]

(take n digitosIniciales)

Sucesión de Lichtenberg

PorJosé A. Alonso 25-enero-20181-febrero-2018

La sucesión de Lichtenberg esta formada por la representación decimal de los números binarios de la sucesión de dígitos 0 y 1 alternados Los primeros términos de ambas sucesiones son

   Alternada ..... Lichtenberg
   0 ....................... 0
   1 ....................... 1
   10 ...................... 2
   101 ..................... 5
   1010 ................... 10
   10101 .................. 21
   101010 ................. 42
   1010101 ................ 85
   10101010 .............. 170
   101010101 ............. 341
   1010101010 ............ 682
   10101010101 .......... 1365
   101010101010 ......... 2730

Alternada ..... Lichtenberg

0 ....................... 0

1 ....................... 1

10 ...................... 2

101 ..................... 5

1010 ................... 10

10101 .................. 21

101010 ................. 42

1010101 ................ 85

10101010 .............. 170

101010101 ............. 341

1010101010 ............ 682

10101010101 .......... 1365

101010101010 ......... 2730

Definir las funciones

   lichtenberg        :: [Integer]
   graficaLichtenberg :: Int -> IO ()

1 2	lichtenberg :: [Integer] graficaLichtenberg :: Int -> IO ()

tales que

lichtenberg es la lista cuyos elementos son los términos de la sucesión de Lichtenberg. Por ejemplo,

     λ> take 17 lichtenberg
     [0,1,2,5,10,21,42,85,170,341,682,1365,2730,5461,10922,21845,43690]

1 2	λ> take 17 lichtenberg [0,1,2,5,10,21,42,85,170,341,682,1365,2730,5461,10922,21845,43690]

(graficaLichtenberg n) dibuja la gráfica del número de dígitos de los n primeros términos de la sucesión de Lichtenberg. Por ejemlo, (graficaLichtenberg 100) dibuja

Comprobar con QuickCheck que todos los términos de la sucesión de Lichtenberg, a partir del 4º, son números compuestos.

Soluciones

import Data.Char (digitToInt)
import Graphics.Gnuplot.Simple
import Test.QuickCheck
import Data.Numbers.Primes (isPrime)

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

lichtenberg1 :: [Integer]
lichtenberg1 = map binarioAdecimal sucAlternada

-- sucAlternada es la lista cuyos elementos son los términos de la
-- sucesión de los dígitos 0 y 1 alternados. Por ejemplo,
--    λ> take 7 sucAlternada
--    ["0","1","10","101","1010","10101","101010"]
sucAlternada :: [String]
sucAlternada =
  ['0'] : [take n cadenaAlternada | n <- [1..]]

-- cadenaAltenada es la cadena formada alternando los caracteres 1 y
-- 0. Por ejemplo,
--    take 20 cadenaAlternada  ==  "10101010101010101010"
cadenaAlternada :: String
cadenaAlternada = cycle ['1','0']

-- (binarioAdecimal cs) es el número decimal correspondiente al número
-- binario cuya cadena de dígitos es cs. Por ejemplo,
--    binarioAdecimal "11101"  ==  29
binarioAdecimal :: String -> Integer
binarioAdecimal =
  foldl (\acc x -> acc * 2 + (toInteger . digitToInt) x) 0
  
-- 2ª solución
lichtenberg2 :: [Integer]
lichtenberg2 = map a [0..]
  where a 0 = 0
        a 1 = 1
        a n = a (n-1) + 2 * a (n-2) + 1

-- 3ª solución
lichtenberg3 :: [Integer]
lichtenberg3 =
  0 : 1 : map (+1) (zipWith (+) (tail lichtenberg3) (map (*2) lichtenberg3)) 

-- Comprobación de eficiencia
--    λ> length (show (lichtenberg1 !! 27))
--    8
--    (0.02 secs, 155,384 bytes)
--    λ> length (show (lichtenberg2 !! 27))
--    8
--    (2.22 secs, 311,157,760 bytes)
--    
--    λ> length (show (lichtenberg1 !! (8*10^4)))
--    24083
--    (1.28 secs, 664,207,040 bytes)
--    λ> length (show (lichtenberg3 !! (8*10^4)))
--    24083
--    (2.59 secs, 1,253,328,200 bytes)

-- La propiedad es
propLichtenberg :: Int -> Property
propLichtenberg n =
  n > 4 ==> not (isPrime (lichtenberg1 !! n))

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

graficaLichtenberg :: Int -> IO ()
graficaLichtenberg n =
  plotList [ Key Nothing
           , Title "Numero de digitos de la sucesion de Lichtenberg"
           , PNG "Sucesion_de_Lichtenberg.png"
           ]
           (take n (map (length . show) lichtenberg1))

import Data.Char (digitToInt)

import Graphics.Gnuplot.Simple

import Test.QuickCheck

import Data.Numbers.Primes (isPrime)

-- 1ª solución

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

lichtenberg1 :: [Integer]

lichtenberg1 = map binarioAdecimal sucAlternada

-- sucAlternada es la lista cuyos elementos son los términos de la

-- sucesión de los dígitos 0 y 1 alternados. Por ejemplo,

-- λ> take 7 sucAlternada

-- ["0","1","10","101","1010","10101","101010"]

sucAlternada :: [String]

sucAlternada =

['0'] : [take n cadenaAlternada | n <- [1..]]

-- cadenaAltenada es la cadena formada alternando los caracteres 1 y

-- 0. Por ejemplo,

-- take 20 cadenaAlternada == "10101010101010101010"

cadenaAlternada :: String

cadenaAlternada = cycle ['1','0']

-- (binarioAdecimal cs) es el número decimal correspondiente al número

-- binario cuya cadena de dígitos es cs. Por ejemplo,

-- binarioAdecimal "11101" == 29

binarioAdecimal :: String -> Integer

binarioAdecimal =

foldl (\acc x -> acc * 2 + (toInteger . digitToInt) x) 0

-- 2ª solución

lichtenberg2 :: [Integer]

lichtenberg2 = map a [0..]

where a 0 = 0

a 1 = 1

a n = a (n-1) + 2 * a (n-2) + 1

-- 3ª solución

lichtenberg3 :: [Integer]

lichtenberg3 =

0 : 1 : map (+1) (zipWith (+) (tail lichtenberg3) (map (*2) lichtenberg3))

-- Comprobación de eficiencia

-- λ> length (show (lichtenberg1 !! 27))

-- 8

-- (0.02 secs, 155,384 bytes)

-- λ> length (show (lichtenberg2 !! 27))

-- 8

-- (2.22 secs, 311,157,760 bytes)

-- λ> length (show (lichtenberg1 !! (8*10^4)))

-- 24083

-- (1.28 secs, 664,207,040 bytes)

-- λ> length (show (lichtenberg3 !! (8*10^4)))

-- 24083

-- (2.59 secs, 1,253,328,200 bytes)

-- La propiedad es

propLichtenberg :: Int -> Property

propLichtenberg n =

n > 4 ==> not (isPrime (lichtenberg1 !! n))

-- La comprobación es

-- λ> quickCheck propLichtenberg

-- +++ OK, passed 100 tests.

graficaLichtenberg :: Int -> IO ()

graficaLichtenberg n =

plotList [ Key Nothing

, Title "Numero de digitos de la sucesion de Lichtenberg"

, PNG "Sucesion_de_Lichtenberg.png"

]

(take n (map (length . show) lichtenberg1))

Complemento potencial

PorJosé A. Alonso 22-enero-201829-enero-2018

El complemento potencial de un número entero positivo x es el menor número y tal que el producto de x por y es un una potencia perfecta. Por ejemplo,

el complemento potencial de 12 es 3 ya que 12 y 24 no son potencias perfectas pero 36 sí lo es;
el complemento potencial de 54 es 4 ya que 54, 108 y 162 no son potencias perfectas pero 216 = 6^3 sí lo es.

Definir las funciones

   complemento                 :: Integer -> Integer
   graficaComplementoPotencial :: Integer -> IO ()

1 2	complemento :: Integer -> Integer graficaComplementoPotencial :: Integer -> IO ()

tales que

(complemento x) es el complemento potencial de x; por ejemplo,

     complemento 12     ==  3
     complemento 54     ==  4
     complemento 720    ==  5
     complemento 24000  ==  9
     complemento 2018   ==  2018

complemento 12 == 3

complemento 54 == 4

complemento 720 == 5

complemento 24000 == 9

complemento 2018 == 2018

(graficaComplementoPotencial n) dibuja la gráfica de los complementos potenciales de los n primeros números enteros positivos. Por ejemplo, (graficaComplementoPotencial 100) dibuja

y (graficaComplementoPotencial 500) dibuja

Comprobar con QuickCheck que (complemento x) es menor o igual que x.

Soluciones

import Data.Numbers.Primes     (primeFactors)
import Data.List               (genericLength, group)
import Graphics.Gnuplot.Simple (plotList, Attribute (Key, PNG, Title))
import Test.QuickCheck

complemento :: Integer -> Integer
complemento 1 = 1
complemento x =
  head [y | y <- [1..]
          , esPotenciaPerfecta (x*y)]
  
-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta 36  ==  True
--    esPotenciaPerfecta 72  ==  False
esPotenciaPerfecta :: Integer -> Bool
esPotenciaPerfecta x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de la factorización prima
-- de x. Por ejemplos,
--    exponentes 36  ==  [2,2]
--    exponentes 72  ==  [3,2]
exponentes :: Integer -> [Integer]
exponentes = map snd . factorizacion
  
-- (factorizacion n) es la factorizacion prima de n. Por ejemplo,
--    factorizacion 1400  ==  [(2,3),(5,2),(7,1)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(x,genericLength xs) | xs@(x:_) <- group (primeFactors n)]

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,
--    mcd [4,6,10]  ==  2
--    mcd [4,5,10]  ==  1
mcd :: [Integer] -> Integer
mcd = foldl1 gcd

-- La propiedad es
prop_complemento :: (Positive Integer) -> Bool
prop_complemento (Positive x) =
  complemento x <= x

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

graficaComplementoPotencial :: Integer -> IO ()
graficaComplementoPotencial n =
  plotList [Key Nothing
           , PNG ("Complemento_potencial_"  ++ show n ++ ".png")
           , Title ("(graficaComplementoPotencial " ++ show n ++ ")") 
           ]
           (map complemento [1..n])

import Data.Numbers.Primes (primeFactors)

import Data.List (genericLength, group)

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

import Test.QuickCheck

complemento :: Integer -> Integer

complemento 1 = 1

complemento x =

head [y | y <- [1..]

, esPotenciaPerfecta (x*y)]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por

-- ejemplo,

-- esPotenciaPerfecta 36 == True

-- esPotenciaPerfecta 72 == False

esPotenciaPerfecta :: Integer -> Bool

esPotenciaPerfecta x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de la factorización prima

-- de x. Por ejemplos,

-- exponentes 36 == [2,2]

-- exponentes 72 == [3,2]

exponentes :: Integer -> [Integer]

exponentes = map snd . factorizacion

-- (factorizacion n) es la factorizacion prima de n. Por ejemplo,

-- factorizacion 1400 == [(2,3),(5,2),(7,1)]

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

factorizacion n =

[(x,genericLength xs) | xs@(x:_) <- group (primeFactors n)]

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,

-- mcd [4,6,10] == 2

-- mcd [4,5,10] == 1

mcd :: [Integer] -> Integer

mcd = foldl1 gcd

-- La propiedad es

prop_complemento :: (Positive Integer) -> Bool

prop_complemento (Positive x) =

complemento x <= x

-- La comprobación es

-- λ> quickCheck prop_complemento

-- +++ OK, passed 100 tests.

graficaComplementoPotencial :: Integer -> IO ()

graficaComplementoPotencial n =

plotList [Key Nothing

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

, Title ("(graficaComplementoPotencial " ++ show n ++ ")")

]

(map complemento [1..n])

Terna pitagórica a partir de un lado

PorJosé A. Alonso 19-enero-201819-febrero-2018

Una terna pitagórica con primer lado x es una terna (x,y,z) tal que x^2 + y^2 = z^2. Por ejemplo, las ternas pitagóricas con primer lado 16 son (16,12,20), (16,30,34) y (16,63,65).

Definir las funciones

   ternasPitagoricas      :: Integer -> [(Integer,Integer,Integer)]
   mayorTernaPitagorica   :: Integer -> (Integer,Integer,Integer)
   graficaMayorHipotenusa :: Integer -> IO ()

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

mayorTernaPitagorica :: Integer -> (Integer,Integer,Integer)

graficaMayorHipotenusa :: Integer -> IO ()

tales que

(ternasPitgoricas x) es la lista de las ternas pitagóricas con primer lado x. Por ejemplo,

     ternasPitagoricas 16 == [(16,12,20),(16,30,34),(16,63,65)]
     ternasPitagoricas 20 == [(20,15,25),(20,21,29),(20,48,52),(20,99,101)]
     ternasPitagoricas 25 == [(25,60,65),(25,312,313)]
     ternasPitagoricas 26 == [(26,168,170)]

ternasPitagoricas 16 == [(16,12,20),(16,30,34),(16,63,65)]

ternasPitagoricas 20 == [(20,15,25),(20,21,29),(20,48,52),(20,99,101)]

ternasPitagoricas 25 == [(25,60,65),(25,312,313)]

ternasPitagoricas 26 == [(26,168,170)]

(mayorTernaPitagorica x) es la mayor de las ternas pitagóricas con primer lado x. Por ejemplo,

     mayorTernaPitagorica 16     ==  (16,63,65)
     mayorTernaPitagorica 20     ==  (20,99,101)
     mayorTernaPitagorica 25     ==  (25,312,313)
     mayorTernaPitagorica 26     ==  (26,168,170)
     mayorTernaPitagorica 2018   ==  (2018,1018080,1018082)
     mayorTernaPitagorica 2019   ==  (2019,2038180,2038181)

mayorTernaPitagorica 16 == (16,63,65)

mayorTernaPitagorica 20 == (20,99,101)

mayorTernaPitagorica 25 == (25,312,313)

mayorTernaPitagorica 26 == (26,168,170)

mayorTernaPitagorica 2018 == (2018,1018080,1018082)

mayorTernaPitagorica 2019 == (2019,2038180,2038181)

(graficaMayorHipotenusa n) dibuja la gráfica de las sucesión de las mayores hipotenusas de las ternas pitagóricas con primer lado x, para x entre 3 y n. Por ejemplo, (graficaMayorHipotenusa 100) dibuja

Soluciones

import Graphics.Gnuplot.Simple

-- Definición de ternasPitagoricas
-- ===============================

ternasPitagoricas :: Integer -> [(Integer,Integer,Integer)]
ternasPitagoricas x =
  [(x,y,z) | y <- [1..(x^ 2 - 1) `div` 2 ]
           , z <- raizCuadrada (x^2 + y^2)]

-- La justificación de la cota es
--    x > 2
--    x^2 + y^2 >= (y+1)^2
--    x^2 + y^2 >= y^2 + 2*y + 1
--    y =< (x^ 2 - 1) `div` 2 

-- (raizCuadrada x) es la lista formada por la raíz cuadrada entera de
-- x, si existe y la lista vacía, en caso contrario. Por ejemplo, 
--    raizCuadrada 25  ==  [5]
--    raizCuadrada 26  ==  []
raizCuadrada :: Integer -> [Integer]
raizCuadrada x =
  [y | y <- [(round . sqrt . fromIntegral) x]
     , y^2 == x]


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

mayorTernaPitagorica :: Integer -> (Integer,Integer,Integer)
mayorTernaPitagorica =
  last . ternasPitagoricas

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

mayorTernaPitagorica2 :: Integer -> (Integer,Integer,Integer)
mayorTernaPitagorica2 x =
  head [(x,y,z) | y <- [k, k-1 .. 1]
                , z <- raizCuadrada (x^2 + y^2)]
  where k = (x^2 - 1) `div` 2

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

-- Se supone que x > 2. Se consideran dos casos:
-- 
-- Primer caso: Supongamos que x es par. Entonces x^2 > 4 y es divisible
-- por 4. Por tanto, existe un y tal que x^2 = 4*y + 4; luego,
--    x^2 + y^2 = 4*y + 4 + y^2
--              = (y + 2)^2
-- La terna es (x,y,y+2) donde y = (x^2 - 4) / 4.
--
-- Segundo caso: Supongamos que x es impar. Entonces x^2 es impar. Por
-- tanto, existe un y tal que x^2 = 2*y + 1; luego,
--    x^2 + y^2 = 2*y + 1 + y^2
--              = (y+1)^2
-- La terna es (x,y,y+1) donde y = (x^2 - 1) / 2.

mayorTernaPitagorica3 :: Integer -> (Integer,Integer,Integer)
mayorTernaPitagorica3 x
  | even x    = (x, y1, y1 + 2)
  | otherwise = (x, y2, y2 + 1)
    where y1 = (x^2 - 4) `div` 4
          y2 = (x^2 - 1) `div` 2 

-- Comparación de eficiencia
--    λ> mayorTernaPitagorica 1006
--    (1006,253008,253010)
--    (7.36 secs, 1,407,793,992 bytes)
--    λ> mayorTernaPitagorica2 1006
--    (1006,253008,253010)
--    (3.76 secs, 704,007,456 bytes)
--    λ> mayorTernaPitagorica3 1006
--    (1006,253008,253010)
--    (0.01 secs, 157,328 bytes)

graficaMayorHipotenusa :: Integer -> IO ()
graficaMayorHipotenusa n =
  plotList [ Key Nothing
           , PNG "Terna_pitagorica_a_partir_de_un_lado.png"
           ]
           [(x,z) | x <- [3..n]
                  , let (_,_,z) = mayorTernaPitagorica3 x]

import Graphics.Gnuplot.Simple

-- Definición de ternasPitagoricas

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

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

ternasPitagoricas x =

[(x,y,z) | y <- [1..(x^ 2 - 1) `div` 2 ]

, z <- raizCuadrada (x^2 + y^2)]

-- La justificación de la cota es

-- x > 2

-- x^2 + y^2 >= (y+1)^2

-- x^2 + y^2 >= y^2 + 2*y + 1

-- y =< (x^ 2 - 1) `div` 2

-- (raizCuadrada x) es la lista formada por la raíz cuadrada entera de

-- x, si existe y la lista vacía, en caso contrario. Por ejemplo,

-- raizCuadrada 25 == [5]

-- raizCuadrada 26 == []

raizCuadrada :: Integer -> [Integer]

raizCuadrada x =

[y | y <- [(round . sqrt . fromIntegral) x]

, y^2 == x]

-- 1ª definición de mayorTernaPitagorica

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

mayorTernaPitagorica :: Integer -> (Integer,Integer,Integer)

mayorTernaPitagorica =

last . ternasPitagoricas

-- 2ª definición de mayorTernaPitagorica

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

mayorTernaPitagorica2 :: Integer -> (Integer,Integer,Integer)

mayorTernaPitagorica2 x =

head [(x,y,z) | y <- [k, k-1 .. 1]

, z <- raizCuadrada (x^2 + y^2)]

where k = (x^2 - 1) `div` 2

-- 3ª definición de mayorTernaPitagorica

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

-- Se supone que x > 2. Se consideran dos casos:

-- Primer caso: Supongamos que x es par. Entonces x^2 > 4 y es divisible

-- por 4. Por tanto, existe un y tal que x^2 = 4*y + 4; luego,

-- x^2 + y^2 = 4*y + 4 + y^2

-- = (y + 2)^2

-- La terna es (x,y,y+2) donde y = (x^2 - 4) / 4.

-- Segundo caso: Supongamos que x es impar. Entonces x^2 es impar. Por

-- tanto, existe un y tal que x^2 = 2*y + 1; luego,

-- x^2 + y^2 = 2*y + 1 + y^2

-- = (y+1)^2

-- La terna es (x,y,y+1) donde y = (x^2 - 1) / 2.

mayorTernaPitagorica3 :: Integer -> (Integer,Integer,Integer)

mayorTernaPitagorica3 x

| even x = (x, y1, y1 + 2)

| otherwise = (x, y2, y2 + 1)

where y1 = (x^2 - 4) `div` 4

y2 = (x^2 - 1) `div` 2

-- Comparación de eficiencia

-- λ> mayorTernaPitagorica 1006

-- (1006,253008,253010)

-- (7.36 secs, 1,407,793,992 bytes)

-- λ> mayorTernaPitagorica2 1006

-- (1006,253008,253010)

-- (3.76 secs, 704,007,456 bytes)

-- λ> mayorTernaPitagorica3 1006

-- (1006,253008,253010)

-- (0.01 secs, 157,328 bytes)

graficaMayorHipotenusa :: Integer -> IO ()

graficaMayorHipotenusa n =

plotList [ Key Nothing

, PNG "Terna_pitagorica_a_partir_de_un_lado.png"

]

[(x,z) | x <- [3..n]

, let (_,_,z) = mayorTernaPitagorica3 x]

Escalada hasta un primo

PorJosé A. Alonso 18-enero-201819-febrero-2018

Este ejercicio está basado en el artículo La conjetura de la «escalada hasta un primo» publicado esta semana por Miguel Ángel Morales en su blog Gaussianos.

La conjetura de escalada hasta un primo trata, propuesta por John Horton Conway, es sencilla de plantear, pero primero vamos a ver qué es eso de escalar hasta un primo. Tomamos un número cualquiera y lo descomponemos en factores primos (colocados en orden ascendente). Si el número era primo, ya hemos acabado; si no era primo, construimos el número formado por los factores primos y los exponentes de los mismos colocados tal cual salen en la factorización. Con el número obtenido hacemos lo mismo que antes. La escalada finaliza cuando obtengamos un número primo. Por ejemplo, para obtener la escalada prima de 1400, como no es primo, se factoriza (obteniéndose 2^3 * 5^2 * 7) y se unen bases y exponentes (obteniéndose 23527). Con el 23527 se repite el proceso obteniéndose la factorización (7 * 3361) y su unión (73361). Como el 73361 es primo, termina la escalada. Por tanto, la escalada de 1400 es [1400,23527,73361].

La conjetura de Conway sobre «escalada hasta un primo» dice que todo número natural mayor o igual que 2 termina su escalada en un número primo.

Definir las funciones

   escaladaPrima                :: Integer -> [Integer]
   longitudEscaladaPrima        :: Integer -> Integer
   longitudEscaladaPrimaAcotada :: Integer -> Integer -> Integer
   graficaEscalada              :: Integer -> Integer -> IO ()

escaladaPrima :: Integer -> [Integer]

longitudEscaladaPrima :: Integer -> Integer

longitudEscaladaPrimaAcotada :: Integer -> Integer -> Integer

graficaEscalada :: Integer -> Integer -> IO ()

tales que

(escaladaPrima n) es la escalada prima de n. Por ejemplo,

     λ> escaladaPrima 1400
     [1400,23527,73361]
     λ> escaladaPrima 333
     [333,3237,31383,3211317,3217139151,39722974813]
     λ> take 10 (escaladaPrima 20)
     [20,225,3252,223271,297699,399233,715623,3263907,32347303,160720129]
     λ> take 3 (escaladaPrima 13532385396179)
     [13532385396179,13532385396179,13532385396179]
     λ> take 2 (escaladaPrima 45214884853168941713016664887087462487)
     [45214884853168941713016664887087462487,13532385396179]

λ> escaladaPrima 1400

[1400,23527,73361]

λ> escaladaPrima 333

[333,3237,31383,3211317,3217139151,39722974813]

λ> take 10 (escaladaPrima 20)

[20,225,3252,223271,297699,399233,715623,3263907,32347303,160720129]

λ> take 3 (escaladaPrima 13532385396179)

[13532385396179,13532385396179,13532385396179]

λ> take 2 (escaladaPrima 45214884853168941713016664887087462487)

[45214884853168941713016664887087462487,13532385396179]

(longitudEscaladaPrima n) es la longitud de la escalada prima de n. Por ejemplo,

     longitudEscaladaPrima 1400  ==  3
     longitudEscaladaPrima 333   ==  6

1 2	longitudEscaladaPrima 1400 == 3 longitudEscaladaPrima 333 == 6

(longitudEscaladaPrimaAcotada n k) es el mínimo entre la longitud de la escalada prima de n y k. Por ejemplo,

     longitudEscaladaPrimaAcotada 333 10  ==  6
     longitudEscaladaPrimaAcotada 333 4   ==  4
     longitudEscaladaPrimaAcotada 20 4    ==  4

longitudEscaladaPrimaAcotada 333 10 == 6

longitudEscaladaPrimaAcotada 333 4 == 4

longitudEscaladaPrimaAcotada 20 4 == 4

(graficaEscalada n k) dibuja la gráfica de (longitudEscaladaPrimaAcotada x k) para x entre 2 y n. Por ejemplo, (graficaEscalada 120 15) dibuja

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes
import Graphics.Gnuplot.Simple

-- Definición de escaladaPrima 
escaladaPrima :: Integer -> [Integer]
escaladaPrima n
  | isPrime n = [n]
  | otherwise = n : escaladaPrima (siguiente n)

-- (siguiente n) es el siguiente de n en la escalada hacia un primo. Por
-- ejemplo, 
--    siguiente 1400  ==  23527
siguiente :: Integer -> Integer
siguiente = paresAnumero . factorizacion

-- (factorizacion n) es la factorizacion prima de n. Por ejemplo,
--    factorizacion 1400  ==  [(2,3),(5,2),(7,1)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(x,genericLength xs) | xs@(x:_) <- group (primeFactors n)]

-- (paAnumero p) es el número obtenido pegando las dos componentes de p
-- (si la segunda es mayor que 1) o sólo la primera componente, en caso
-- contrario. Por ejemplo, 
--    parAnumero (7,1)  ==  7
--    parAnumero (23,5)  ==  235
parAnumero :: (Integer,Integer) -> Integer
parAnumero (n,1) = n
parAnumero (n,k) = read (show n ++ show k)

-- (paresA numeros ps) es el número obtenido aplicando parAnumero a cada
-- uno de los pares de ps. Por ejemplo,
--    paresAnumero [(2,3),(5,2),(7,1)]  ==  23527
paresAnumero :: [(Integer,Integer)] -> Integer
paresAnumero = read . concatMap (show . parAnumero)

-- Definición de longitudEscaladaPrima
longitudEscaladaPrima :: Integer -> Integer
longitudEscaladaPrima n
  | isPrime n = 1
  | otherwise = 1 + longitudEscaladaPrima (siguiente n)

-- Definición de longitudEscaladaPrimaAcotada
longitudEscaladaPrimaAcotada :: Integer -> Integer -> Integer
longitudEscaladaPrimaAcotada _ 0 = 0
longitudEscaladaPrimaAcotada n k
  | isPrime n = 1
  | otherwise = 1 + longitudEscaladaPrimaAcotada (siguiente n) (k-1)

-- Definición de graficaEscalada 
graficaEscalada :: Integer -> Integer -> IO ()
graficaEscalada n k =
  plotList [ Key Nothing
           , PNG "Escalada_hasta_un_primo.png"
           , Title ("graficaEscalada " ++ show n ++ " " ++ show k)
           , XLabel "Numeros"
           , YLabel "Longitud de escalada"
           ]
           [(x,longitudEscaladaPrimaAcotada x k) | x <- [2..n]]

import Data.List (genericLength, group)

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

-- Definición de escaladaPrima

escaladaPrima :: Integer -> [Integer]

escaladaPrima n

| isPrime n = [n]

| otherwise = n : escaladaPrima (siguiente n)

-- (siguiente n) es el siguiente de n en la escalada hacia un primo. Por

-- ejemplo,

-- siguiente 1400 == 23527

siguiente :: Integer -> Integer

siguiente = paresAnumero . factorizacion

-- (factorizacion n) es la factorizacion prima de n. Por ejemplo,

-- factorizacion 1400 == [(2,3),(5,2),(7,1)]

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

factorizacion n =

[(x,genericLength xs) | xs@(x:_) <- group (primeFactors n)]

-- (paAnumero p) es el número obtenido pegando las dos componentes de p

-- (si la segunda es mayor que 1) o sólo la primera componente, en caso

-- contrario. Por ejemplo,

-- parAnumero (7,1) == 7

-- parAnumero (23,5) == 235

parAnumero :: (Integer,Integer) -> Integer

parAnumero (n,1) = n

parAnumero (n,k) = read (show n ++ show k)

-- (paresA numeros ps) es el número obtenido aplicando parAnumero a cada

-- uno de los pares de ps. Por ejemplo,

-- paresAnumero [(2,3),(5,2),(7,1)] == 23527

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

paresAnumero = read . concatMap (show . parAnumero)

-- Definición de longitudEscaladaPrima

longitudEscaladaPrima :: Integer -> Integer

longitudEscaladaPrima n

| isPrime n = 1

| otherwise = 1 + longitudEscaladaPrima (siguiente n)

-- Definición de longitudEscaladaPrimaAcotada

longitudEscaladaPrimaAcotada :: Integer -> Integer -> Integer

longitudEscaladaPrimaAcotada _ 0 = 0

longitudEscaladaPrimaAcotada n k

| isPrime n = 1

| otherwise = 1 + longitudEscaladaPrimaAcotada (siguiente n) (k-1)

-- Definición de graficaEscalada

graficaEscalada :: Integer -> Integer -> IO ()

graficaEscalada n k =

plotList [ Key Nothing

, PNG "Escalada_hasta_un_primo.png"

, Title ("graficaEscalada " ++ show n ++ " " ++ show k)

, XLabel "Numeros"

, YLabel "Longitud de escalada"

]

[(x,longitudEscaladaPrimaAcotada x k) | x <- [2..n]]

Sucesión de Recamán

PorJosé A. Alonso 24-abril-201726-marzo-2022

La sucesión de Recamán está definida como sigue:

   a(0) = 0
   a(n) = a(n-1) - n, si a(n-1) > n y no figura ya en la sucesión
   a(n) = a(n-1) + n, en caso contrario.

a(0) = 0

a(n) = a(n-1) - n, si a(n-1) > n y no figura ya en la sucesión

a(n) = a(n-1) + n, en caso contrario.

Definir las funciones

   sucRecaman :: [Int]
   invRecaman :: Int -> Int
   graficaSucRecaman :: Int -> IO ()
   graficaInvRecaman :: Int -> IO ()

sucRecaman :: [Int]

invRecaman :: Int -> Int

graficaSucRecaman :: Int -> IO ()

graficaInvRecaman :: Int -> IO ()

tales que

sucRecaman es la lista de los términos de la sucesión de Recamám. Por ejemplo,

      λ> take 25 sucRecaman3
      [0,1,3,6,2,7,13,20,12,21,11,22,10,23,9,24,8,25,43,62,42,63,41,18,42]
      λ> sucRecaman !! 1000
      3686
      λ> sucRecaman !! 1000001
      1057163

λ> take 25 sucRecaman3

[0,1,3,6,2,7,13,20,12,21,11,22,10,23,9,24,8,25,43,62,42,63,41,18,42]

λ> sucRecaman !! 1000

3686

λ> sucRecaman !! 1000001

1057163

(invRecaman n) es la primera posición de n en la sucesión de Recamán. Por ejemplo,

      invRecaman 10       ==  12
      invRecaman 3686     ==  1000
      invRecaman 1057163  ==  1000001

invRecaman 10 == 12

invRecaman 3686 == 1000

invRecaman 1057163 == 1000001

(graficaSucRecaman n) dibuja los n primeros términos de la sucesión de Recamán. Por ejemplo, (graficaSucRecaman 300) dibuja
(graficaInvRecaman n) dibuja los valores de (invRecaman k) para k entre 0 y n. Por ejemplo, (graficaInvRecaman 17) dibuja

y (graficaInvRecaman 100) dibuja

Soluciones

import qualified Data.Set as S

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

sucRecaman1 :: [Int]
sucRecaman1 = map suc1 [0..]

suc1 :: Int -> Int
suc1 0 = 0
suc1 n | y > n && y - n `notElem` ys = y - n
       | otherwise                   = y + n
  where y  = suc1 (n - 1)
        ys = [suc1 k | k <- [0..n - 1]]

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

sucRecaman2 :: [Int]
sucRecaman2 = 0:zipWith3 f sucRecaman2 [1..] (repeat sucRecaman2)
  where f y n ys | y > n && y - n `notElem` take n ys = y - n
                 | otherwise                          = y + n

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

sucRecaman3 :: [Int]
sucRecaman3 = 0 : recaman (S.singleton 0) 1 0

recaman :: S.Set Int -> Int -> Int -> [Int]
recaman s n x
  | x > n && (x-n) `S.notMember` s =
    (x-n) : recaman (S.insert (x-n) s) (n+1) (x-n)
  | otherwise =
    (x+n):recaman (S.insert (x+n) s) (n+1) (x+n) 

-- Comparación de eficiencia:
--    λ> sucRecaman1 !! 25
--    17
--    (3.76 secs, 2,394,593,952 bytes)
--    λ> sucRecaman2 !! 25
--    17
--    (0.00 secs, 0 bytes)
--    λ> sucRecaman3 !! 25
--    17
--    (0.00 secs, 0 bytes)
--
--    λ> sucRecaman2 !! (2*10^4)
--    14358
--    (2.69 secs, 6,927,559,784 bytes)
--    λ> sucRecaman3 !! (2*10^4)
--    14358
--    (0.04 secs, 0 bytes)

-- Definición de invRecaman
invRecaman :: Int -> Int
invRecaman n =
  length (takeWhile (/=n) sucRecaman3)

graficaSucRecaman :: Int -> IO ()
graficaSucRecaman n =
  plotList [Key Nothing]
           (take n sucRecaman3)

graficaInvRecaman :: Int -> IO ()
graficaInvRecaman n =
  plotList [Key Nothing]
           [invRecaman k | k <- [0..n]]

import qualified Data.Set as S

-- 1ª solución

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

sucRecaman1 :: [Int]

sucRecaman1 = map suc1 [0..]

suc1 :: Int -> Int

suc1 0 = 0

suc1 n | y > n && y - n `notElem` ys = y - n

| otherwise = y + n

where y = suc1 (n - 1)

ys = [suc1 k | k <- [0..n - 1]]

-- 2ª solución

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

sucRecaman2 :: [Int]

sucRecaman2 = 0:zipWith3 f sucRecaman2 [1..] (repeat sucRecaman2)

where f y n ys | y > n && y - n `notElem` take n ys = y - n

| otherwise = y + n

-- 3ª solución

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

sucRecaman3 :: [Int]

sucRecaman3 = 0 : recaman (S.singleton 0) 1 0

recaman :: S.Set Int -> Int -> Int -> [Int]

recaman s n x

| x > n && (x-n) `S.notMember` s =

(x-n) : recaman (S.insert (x-n) s) (n+1) (x-n)

| otherwise =

(x+n):recaman (S.insert (x+n) s) (n+1) (x+n)

-- Comparación de eficiencia:

-- λ> sucRecaman1 !! 25

-- 17

-- (3.76 secs, 2,394,593,952 bytes)

-- λ> sucRecaman2 !! 25

-- 17

-- (0.00 secs, 0 bytes)

-- λ> sucRecaman3 !! 25

-- 17

-- (0.00 secs, 0 bytes)

-- λ> sucRecaman2 !! (2*10^4)

-- 14358

-- (2.69 secs, 6,927,559,784 bytes)

-- λ> sucRecaman3 !! (2*10^4)

-- 14358

-- (0.04 secs, 0 bytes)

-- Definición de invRecaman

invRecaman :: Int -> Int

invRecaman n =

length (takeWhile (/=n) sucRecaman3)

graficaSucRecaman :: Int -> IO ()

graficaSucRecaman n =

plotList [Key Nothing]

(take n sucRecaman3)

graficaInvRecaman :: Int -> IO ()

graficaInvRecaman n =

plotList [Key Nothing]

[invRecaman k | k <- [0..n]]

Distancias entre primos consecutivos

PorJosé A. Alonso 4-abril-201711-abril-2017

Los 15 primeros números primos son

   2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47

1	2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47

Las distancias entre los elementos consecutivos son

   1, 2, 2, 4, 2,  4,  2,  4,  6,  2,  6,  4,  2,  4

1	1, 2, 2, 4, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4

La distribución de las distancias es

   (1,1), (2,6), (4,5), (6,2)

1	(1,1), (2,6), (4,5), (6,2)

(es decir, el 1 aparece una vez, el 2 aparece 6 veces, etc.) La frecuencia de las distancias es

   (1,7.142857), (2,42.857143), (4,35.714287), (6,14.285714)

1	(1,7.142857), (2,42.857143), (4,35.714287), (6,14.285714)

(es decir, el 1 aparece el 7.142857%, el 2 el 42.857143% etc.)

Definir las funciones

  cuentaDistancias        :: Int -> [(Int,Int)]
  frecuenciasDistancias   :: Int -> [(Int,Float)]
  graficas                :: [Int] -> IO ()
  distanciasMasFrecuentes :: Int -> [Int]

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

frecuenciasDistancias :: Int -> [(Int,Float)]

graficas :: [Int] -> IO ()

distanciasMasFrecuentes :: Int -> [Int]

tales que

(cuentaDistancias n) es la distribución de distancias entre los n primeros primos consecutivos. Por ejemplo,

     λ> cuentaDistancias 15
     [(1,1),(2,6),(4,5),(6,2)]
     λ> cuentaDistancias 100
     [(1,1),(2,25),(4,26),(6,25),(8,7),(10,7),(12,4),(14,3),(18,1)]

λ> cuentaDistancias 15

[(1,1),(2,6),(4,5),(6,2)]

λ> cuentaDistancias 100

[(1,1),(2,25),(4,26),(6,25),(8,7),(10,7),(12,4),(14,3),(18,1)]

(frecuenciasDistancias n) es la frecuencia de distancias entre los n primeros primos consecutivos. Por ejemplo,

     λ> frecuenciasDistancias 15
     [(1,7.142857),(2,42.857143),(4,35.714287),(6,14.285714)]
     λ> frecuenciasDistancias 30
     [(1,3.4482758),(2,34.482758),(4,34.482758),(6,24.137932),(8,3.4482758)]

λ> frecuenciasDistancias 15

[(1,7.142857),(2,42.857143),(4,35.714287),(6,14.285714)]

λ> frecuenciasDistancias 30

[(1,3.4482758),(2,34.482758),(4,34.482758),(6,24.137932),(8,3.4482758)]

(graficas ns) dibuja las gráficas de (frecuenciasDistancias k) para k en ns. Por ejemplo, (graficas [10,20,30]) dibuja

(graficas [1000,2000,3000]) dibuja

y (graficas [100000,200000,300000]) dibuja
(distanciasMasFrecuentes n) es la lista de las distancias más frecuentes entre los elementos consecutivos de la lista de los n primeros primos. Por ejemplo,

     distanciasMasFrecuentes 15     ==  [2]
     distanciasMasFrecuentes 26     ==  [2,4]
     distanciasMasFrecuentes 32     ==  [4]
     distanciasMasFrecuentes 41     ==  [2,4,6]
     distanciasMasFrecuentes 77     ==  [6]
     distanciasMasFrecuentes 160    ==  [4,6]
     distanciasMasFrecuentes (10^6) == [6]

distanciasMasFrecuentes 15 == [2]

distanciasMasFrecuentes 26 == [2,4]

distanciasMasFrecuentes 32 == [4]

distanciasMasFrecuentes 41 == [2,4,6]

distanciasMasFrecuentes 77 == [6]

distanciasMasFrecuentes 160 == [4,6]

distanciasMasFrecuentes (10^6) == [6]

Comprobar con QuickCheck si para todo n > 160 se verifica que (distanciasMasFrecuentes n) es [6].

Soluciones

import Data.Numbers.Primes
import qualified Data.Map as M 
import Graphics.Gnuplot.Simple
import Test.QuickCheck

cuentaDistancias :: Int -> [(Int,Int)]
cuentaDistancias = M.toList . dicDistancias

dicDistancias :: Int -> M.Map Int Int
dicDistancias n =
  M.fromListWith (+) (zip (distancias n) (repeat 1))

distancias :: Int -> [Int]
distancias n =
  zipWith (-) (tail xs) xs
  where xs = take n primes

frecuenciasDistancias :: Int -> [(Int,Float)]
frecuenciasDistancias n =
  [(k,(100 * fromIntegral x) / n1) | (k,x) <- cuentaDistancias n]
  where n1 = fromIntegral (n-1)

graficas :: [Int] -> IO ()
graficas ns =
  plotLists [Key Nothing]
            (map frecuenciasDistancias ns)

distanciasMasFrecuentes :: Int -> [Int]
distanciasMasFrecuentes n =
  M.keys (M.filter (==m) d)
  where d = dicDistancias n
        m = maximum (M.elems d)

-- La propiedad es
prop_distanciasMasFrecuentes :: Int -> Bool
prop_distanciasMasFrecuentes n =
  distanciasMasFrecuentes (161 + abs n) == [6]

import Data.Numbers.Primes

import qualified Data.Map as M

import Graphics.Gnuplot.Simple

import Test.QuickCheck

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

cuentaDistancias = M.toList . dicDistancias

dicDistancias :: Int -> M.Map Int Int

dicDistancias n =

M.fromListWith (+) (zip (distancias n) (repeat 1))

distancias :: Int -> [Int]

distancias n =

zipWith (-) (tail xs) xs

where xs = take n primes

frecuenciasDistancias :: Int -> [(Int,Float)]

frecuenciasDistancias n =

[(k,(100 * fromIntegral x) / n1) | (k,x) <- cuentaDistancias n]

where n1 = fromIntegral (n-1)

graficas :: [Int] -> IO ()

graficas ns =

plotLists [Key Nothing]

(map frecuenciasDistancias ns)

distanciasMasFrecuentes :: Int -> [Int]

distanciasMasFrecuentes n =

M.keys (M.filter (==m) d)

where d = dicDistancias n

m = maximum (M.elems d)

-- La propiedad es

prop_distanciasMasFrecuentes :: Int -> Bool

prop_distanciasMasFrecuentes n =

distanciasMasFrecuentes (161 + abs n) == [6]

Distribución de diferencias de dígitos consecutivos de pi

PorJosé A. Alonso 14-marzo-201725-abril-2021

La distribución de las diferencias de los dígitos consecutivos para los 18 primeros dígitos de pi se calcula como sigue: los primeros 18 dígitos de pi son

   3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3

1	3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 3, 2, 3

Las diferencias de sus elementos consecutivos es

   2, -3, 3, -4, -4, 7, -4, 1, 2, -2, -3, -1, 2, -2, 6, 1, -1

1	2, -3, 3, -4, -4, 7, -4, 1, 2, -2, -3, -1, 2, -2, 6, 1, -1

y la distribución de sus frecuencias en el intervalo [-9,9] es

   0, 0, 0, 0, 0, 3, 2, 2, 2, 0, 2, 3, 1, 0, 0, 1, 1, 0, 0

1	0, 0, 0, 0, 0, 3, 2, 2, 2, 0, 2, 3, 1, 0, 0, 1, 1, 0, 0

es decir, el desde el -9 a -5 no aparecen, el -4 aparece 3 veces, el -2 aparece 2 veces y así sucesivamente.

Definir las funciones

   distribucionDDCpi :: Int -> [Int]
   graficas :: [Int] -> FilePath -> IO ()

1 2	distribucionDDCpi :: Int -> [Int] graficas :: [Int] -> FilePath -> IO ()

tales que

(distribucionDDCpi n) es la distribución de las diferencias de los dígitos consecutivos para los primeros n dígitos de pi. Por ejemplo,

   λ> distribucionDDCpi 18
   [0,0,0,0,0,3,2,2,2,0,2,3,1,0,0,1,1,0,0]
   λ> distribucionDDCpi 100
   [1,2,1,7,7,7,6,5,8,6,7,14,4,9,3,6,4,1,0]
   λ> distribucionDDCpi 200
   [3,6,2,13,14,12,11,12,15,17,15,19,11,17,8,13,9,2,0]
   λ> distribucionDDCpi 1000
   [16,25,23,44,57,61,55,75,92,98,80,88,64,65,42,54,39,14,8]
   λ> distribucionDDCpi 5000
   [67,99,130,196,245,314,361,391,453,468,447,407,377,304,242,221,134,97,47]

λ> distribucionDDCpi 18

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

λ> distribucionDDCpi 100

[1,2,1,7,7,7,6,5,8,6,7,14,4,9,3,6,4,1,0]

λ> distribucionDDCpi 200

[3,6,2,13,14,12,11,12,15,17,15,19,11,17,8,13,9,2,0]

λ> distribucionDDCpi 1000

[16,25,23,44,57,61,55,75,92,98,80,88,64,65,42,54,39,14,8]

λ> distribucionDDCpi 5000

[67,99,130,196,245,314,361,391,453,468,447,407,377,304,242,221,134,97,47]

(graficas ns f) dibuja en el fichero f las gráficas de las distribuciones de las diferencias de los dígitos consecutivos para los primeros n dígitos de pi, para n en ns. Por ejemplo, al evaluar (graficas [100,250..4000] «distribucionDDCpi.png» se escribe en el fichero «distribucionDDCpi.png» la siguiente gráfica

Nota: Se puede usar la librería Data.Number.CReal.

Soluciones

import Data.Number.CReal
import Graphics.Gnuplot.Simple
import Data.Array

--    λ> digitosPi 18
--    [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3]
digitosPi :: Int -> [Int]
digitosPi n = init [read [c] | c <- (x:xs)]
  where (x:_:xs) = showCReal n pi

--    λ> diferenciasConsecutivos (digitosPi 18)
--    [2,-3,3,-4,-4,7,-4,1,2,-2,-3,-1,2,-2,6,1,-1]
diferenciasConsecutivos :: Num a => [a] -> [a]
diferenciasConsecutivos xs =
  zipWith (-) xs (tail xs)

distribucionDDCpi :: Int -> [Int]
distribucionDDCpi =
  distribucion . diferenciasConsecutivos . digitosPi
  where distribucion xs =
          elems (accumArray (+) 0 (-9,9) (zip xs (repeat 1)))

graficas :: [Int] -> FilePath -> IO ()
graficas ns f = 
  plotLists [Key Nothing, PNG f]
            [puntos n | n <- ns]
  where puntos :: Int -> [(Int,Int)]
        puntos n = zip [-9..9] (distribucionDDCpi n)

import Data.Number.CReal

import Graphics.Gnuplot.Simple

import Data.Array

-- λ> digitosPi 18

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

digitosPi :: Int -> [Int]

digitosPi n = init [read [c] | c <- (x:xs)]

where (x:_:xs) = showCReal n pi

-- λ> diferenciasConsecutivos (digitosPi 18)

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

diferenciasConsecutivos :: Num a => [a] -> [a]

diferenciasConsecutivos xs =

zipWith (-) xs (tail xs)

distribucionDDCpi :: Int -> [Int]

distribucionDDCpi =

distribucion . diferenciasConsecutivos . digitosPi

where distribucion xs =

elems (accumArray (+) 0 (-9,9) (zip xs (repeat 1)))

graficas :: [Int] -> FilePath -> IO ()

graficas ns f =

plotLists [Key Nothing, PNG f]

[puntos n | n <- ns]

where puntos :: Int -> [(Int,Int)]

puntos n = zip [-9..9] (distribucionDDCpi n)

Números de Catalan

PorJosé A. Alonso 13-marzo-201726-marzo-2022

Los números de Catalan forman la sucesión cuyo término general es

Los primeros números de Catalan son

   1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786

1	1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786

Los números de Catalan satisfacen la siguiente relación de recurrencia:

Asintóticamente, los números de Catalan crecen como:

considerando que el cociente entre el n-ésimo número de Catalan y la expresión de la derecha tiende hacia 1 cuando n tiende a infinito.

Definir las funciones

   catalan :: [Integer]
   grafica :: Int -> Int -> IO ()

1 2	catalan :: [Integer] grafica :: Int -> Int -> IO ()

tales que

catalan es la lista de términos de la sucesión de Catalan. Por ejemplo,

     λ> take 12 catalan
     [1,1,2,5,14,42,132,429,1430,4862,16796,58786]
     λ> length (show (catalan !! 50000))
     30096

λ> take 12 catalan

[1,1,2,5,14,42,132,429,1430,4862,16796,58786]

λ> length (show (catalan !! 50000))

30096

(grafica a b) dibuja los n-ésimos términos de la sucesión de Catalan, para n entre a y b, junto con los de la expresión de la derecha de

Por ejemplo, (grafica 5 10) dibuja

y (grafica 55 60) dibuja

Soluciones

import Graphics.Gnuplot.Simple

catalan :: [Integer]
catalan = scanl (\c n -> c*2*(2*n-1) `div` (n+1)) 1 [1..]

grafica :: Int -> Int -> IO ()
grafica a b = 
  plotLists [Key Nothing]
            [[(fromIntegral n, fromIntegral (catalan !! n)) | n <- [a..b]]
            ,[(n,4**n/(n**(3/2)*(sqrt pi))) | n <- [c..d]]]
  where c, d :: Double
        c = fromIntegral a
        d = fromIntegral b

import Graphics.Gnuplot.Simple

catalan :: [Integer]

catalan = scanl (\c n -> c*2*(2*n-1) `div` (n+1)) 1 [1..]

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

grafica a b =

plotLists [Key Nothing]

[[(fromIntegral n, fromIntegral (catalan !! n)) | n <- [a..b]]

,[(n,4**n/(n**(3/2)*(sqrt pi))) | n <- [c..d]]]

where c, d :: Double

c = fromIntegral a

d = fromIntegral b

Cálculo de pi mediante el método de Newton

PorJosé A. Alonso 24-febrero-201725-abril-2021

El método de Newton para el cálculo de pi se basa en la relación

y en el desarrollo del arco seno

de donde se obtiene la fórmula

La primeras aproximaciones son

   a(0) = 6*(1/2)                               = 3.0
   a(1) = 6*(1/2+1/(2*3*2^3))                   = 3.125
   a(2) = 6*(1/2+1/(2*3*2^3)+(1*3)/(2*4*5*2^5)) = 3.1390625

a(0) = 6*(1/2) = 3.0

a(1) = 6*(1/2+1/(2*3*2^3)) = 3.125

a(2) = 6*(1/2+1/(2*3*2^3)+(1*3)/(2*4*5*2^5)) = 3.1390625

Definir las funciones

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

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

tales que

(aproximacionPi n) es la n-ésima aproximación de pi con la fórmula de Newton. Por ejemplo,

     aproximacionPi 0   ==  3.0
     aproximacionPi 1   ==  3.125
     aproximacionPi 2   ==  3.1390625
     aproximacionPi 10  ==  3.1415926468755613
     aproximacionPi 21  ==  3.141592653589793
     pi                 ==  3.141592653589793

aproximacionPi 0 == 3.0

aproximacionPi 1 == 3.125

aproximacionPi 2 == 3.1390625

aproximacionPi 10 == 3.1415926468755613

aproximacionPi 21 == 3.141592653589793

pi == 3.141592653589793

(grafica xs) dibuja la gráfica de las k-ésimas aproximaciones de pi donde k toma los valores de la lista xs. Por ejemplo, (grafica [1..30]) dibuja

Nota: Este ejercicio ha sido propuesto por Manuel Herrera.

Soluciones

import Graphics.Gnuplot.Simple

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

aproximacionPi :: Int -> Double
aproximacionPi n = 6 * arcsinX
  where arcsinX = 0.5 + sum (take n factoresN)

factoresN :: [Double]
factoresN = zipWith (*) (potenciasK 3) fraccionesPI

potenciasK :: Double -> [Double]
potenciasK k = (0.5**k)/k : potenciasK (k+2)

fraccionesPI :: [Double]
fraccionesPI =
  scanl (*) (1/2) (tail (zipWith (/) [1,3..] [2,4..]))

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

aproximacionPi2 :: Int -> Double
aproximacionPi2 n = 6 * (serie !! n)

serie :: [Double]
serie = scanl1 (+) (zipWith (/)
                            (map fromIntegral numeradores)
                            (map fromIntegral denominadores))
  where numeradores    = 1 : scanl1 (*) [1,3..]
        denominadores  = zipWith (*) denominadores1 denominadores2
        denominadores1 = 2 : scanl1 (*) [2,4..]
        denominadores2 = 1 : [n * 2^n | n <- [3,5..]]

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

import Graphics.Gnuplot.Simple

-- 1ª definición

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

aproximacionPi :: Int -> Double

aproximacionPi n = 6 * arcsinX

where arcsinX = 0.5 + sum (take n factoresN)

factoresN :: [Double]

factoresN = zipWith (*) (potenciasK 3) fraccionesPI

potenciasK :: Double -> [Double]

potenciasK k = (0.5**k)/k : potenciasK (k+2)

fraccionesPI :: [Double]

fraccionesPI =

scanl (*) (1/2) (tail (zipWith (/) [1,3..] [2,4..]))

-- 2ª definición

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

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = 6 * (serie !! n)

serie :: [Double]

serie = scanl1 (+) (zipWith (/)

(map fromIntegral numeradores)

(map fromIntegral denominadores))

where numeradores = 1 : scanl1 (*) [1,3..]

denominadores = zipWith (*) denominadores1 denominadores2

denominadores1 = 2 : scanl1 (*) [2,4..]

denominadores2 = 1 : [n * 2^n | n <- [3,5..]]

grafica :: [Int] -> IO ()

grafica xs =

plotList [Key Nothing]

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

Cálculo de pi mediante los métodos de Gregory-Leibniz y de Beeler

PorJosé A. Alonso 23-febrero-201726-marzo-2022

La fórmula de Gregory-Leibniz para calcular pi es

y la de Beeler es

Definir las funciones

   aproximaPiGL     :: Int -> Double
   aproximaPiBeeler :: Int -> Double
   graficas         :: [Int] -> IO ()

aproximaPiGL :: Int -> Double

aproximaPiBeeler :: Int -> Double

graficas :: [Int] -> IO ()

tales que

(aproximaPiGL n) es la aproximación de pi con los primeros n términos de la fórmula de Gregory-Leibniz. Por ejemplo,

     aproximaPiGL 1       ==  4.0
     aproximaPiGL 2       ==  2.666666666666667
     aproximaPiGL 3       ==  3.466666666666667
     aproximaPiGL 10      ==  3.0418396189294032
     aproximaPiGL 100     ==  3.1315929035585537
     aproximaPiGL 1000    ==  3.140592653839794
     aproximaPiGL 10000   ==  3.1414926535900345
     aproximaPiGL 100000  ==  3.1415826535897198

aproximaPiGL 1 == 4.0

aproximaPiGL 2 == 2.666666666666667

aproximaPiGL 3 == 3.466666666666667

aproximaPiGL 10 == 3.0418396189294032

aproximaPiGL 100 == 3.1315929035585537

aproximaPiGL 1000 == 3.140592653839794

aproximaPiGL 10000 == 3.1414926535900345

aproximaPiGL 100000 == 3.1415826535897198

(aproximaPiBeeler n) es la aproximación de pi con los primeros n términos de la fórmula de Beeler. Por ejemplo,

     aproximaPiBeeler 1   ==  2.0
     aproximaPiBeeler 2   ==  2.6666666666666665
     aproximaPiBeeler 3   ==  2.933333333333333
     aproximaPiBeeler 10  ==  3.140578169680337
     aproximaPiBeeler 60  ==  3.141592653589793
     pi                   ==  3.141592653589793

aproximaPiBeeler 1 == 2.0

aproximaPiBeeler 2 == 2.6666666666666665

aproximaPiBeeler 3 == 2.933333333333333

aproximaPiBeeler 10 == 3.140578169680337

aproximaPiBeeler 60 == 3.141592653589793

pi == 3.141592653589793

(graficas xs) dibuja la gráfica de las k-ésimas aproximaciones de pi, donde k toma los valores de la lista xs, con las fórmulas de Gregory-Leibniz y de Beeler. Por ejemplo, (graficas [1..25]) dibuja

donde la línea morada corresponde a la aproximación de Gregory-Leibniz y la verde a la de Beeler.

Nota: Este ejercicio ha sido propuesto por Enrique Naranjo.

Soluciones

import Graphics.Gnuplot.Simple

-- Definiciones de aproximaPiGL
-- ============================

-- 1ª definición de aproximaPiGL
aproximaPiGL :: Int -> Double
aproximaPiGL n = 4 * (sum . take n . sumaA . zipWith (/) [1,1..]) [1,3..]
  where sumaA (x:y:xs) = x:(-y):sumaA xs

-- 2ª definición de aproximaPiGL
aproximaPiGL2 :: Int -> Double
aproximaPiGL2 n =
  4 * (sum (take n (zipWith (/) (cycle [1,-1]) [1,3..])))

-- 3ª definición de aproximaPiGL
aproximaPiGL3 :: Int -> Double
aproximaPiGL3 n =
  4 * (sum . take n . zipWith (/) (cycle [1,-1])) [1,3..]

-- 4ª definición de aproximaPiGL
aproximaPiGL4 :: Int -> Double
aproximaPiGL4 n = serieGL !! (n-1)

serieGL :: [Double]
serieGL = scanl1 (+) (zipWith (/) numeradores denominadores)
  where numeradores   = cycle [4,-4]
        denominadores = [1,3..]

-- Definición de aproximaPiBeeler
aproximaPiBeeler :: Int -> Double
aproximaPiBeeler n = 2 * aux (fromIntegral n) 1
  where
    aux :: Double -> Double -> Double 
    aux n k | n == k    = 1
            | otherwise = 1 + (k/(2*k+1)) * aux n (1+k)

-- Definición de graficas
graficas :: [Int] -> IO ()
graficas xs = 
    plotLists [Key Nothing]
             [[(k,aproximaPiGL k)     | k <- xs],
              [(k,aproximaPiBeeler k) | k <- xs]]

import Graphics.Gnuplot.Simple

-- Definiciones de aproximaPiGL

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

-- 1ª definición de aproximaPiGL

aproximaPiGL :: Int -> Double

aproximaPiGL n = 4 * (sum . take n . sumaA . zipWith (/) [1,1..]) [1,3..]

where sumaA (x:y:xs) = x:(-y):sumaA xs

-- 2ª definición de aproximaPiGL

aproximaPiGL2 :: Int -> Double

aproximaPiGL2 n =

4 * (sum (take n (zipWith (/) (cycle [1,-1]) [1,3..])))

-- 3ª definición de aproximaPiGL

aproximaPiGL3 :: Int -> Double

aproximaPiGL3 n =

4 * (sum . take n . zipWith (/) (cycle [1,-1])) [1,3..]

-- 4ª definición de aproximaPiGL

aproximaPiGL4 :: Int -> Double

aproximaPiGL4 n = serieGL !! (n-1)

serieGL :: [Double]

serieGL = scanl1 (+) (zipWith (/) numeradores denominadores)

where numeradores = cycle [4,-4]

denominadores = [1,3..]

-- Definición de aproximaPiBeeler

aproximaPiBeeler :: Int -> Double

aproximaPiBeeler n = 2 * aux (fromIntegral n) 1

where

aux :: Double -> Double -> Double

aux n k | n == k = 1

| otherwise = 1 + (k/(2*k+1)) * aux n (1+k)

-- Definición de graficas

graficas :: [Int] -> IO ()

graficas xs =

plotLists [Key Nothing]

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

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

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

PorJosé A. Alonso 22-febrero-201726-marzo-2022

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

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

La primeras aproximaciones son

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

a(1) = 3+1 = 4.0

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

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

Definir las funciones

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

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

tales que

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

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

aproximacionPi 1 == 4.0

aproximacionPi 2 == 3.066666666666667

aproximacionPi 3 == 3.158974358974359

aproximacionPi 10 == 3.141287132741557

aproximacionPi 100 == 3.141592398533554

aproximacionPi 1000 == 3.1415926533392926

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

Nota: Este ejercicio ha sido propuesto por Antonio Morales.

Soluciones

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

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

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

import Graphics.Gnuplot.Simple

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

-- par de listas infinitas.

fraccionPi :: [(Integer, Integer)]

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

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

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

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

aproximacionFC n =

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

aproximacionPi :: Int -> Double

aproximacionPi n =

aproximacionFC n fraccionPi

grafica :: [Int] -> IO ()

grafica xs =

plotList [Key Nothing]

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

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

PorJosé A. Alonso 21-febrero-201726-marzo-2022

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

El desarrollo es el siguiente

y se obtiene a partir de la serie armónica

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

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

Por ejemplo,

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

Definir las funciones

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

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

tales que

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

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

aproximacionPi 1 == 1.0

aproximacionPi 10 == 2.3289682539682537

aproximacionPi 100 == 2.934318000847734

aproximacionPi 1000 == 3.0603246224585128

aproximacionPi 10000 == 3.1105295744825403

aproximacionPi 100000 == 3.134308801935256

aproximacionPi 1000000 == 3.1395057903490806

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

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

Soluciones

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

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

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

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

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

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

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

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

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

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

import Data.Numbers.Primes

import Graphics.Gnuplot.Simple

-- 1ª definición

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

aproximacionPi :: Int -> Double

aproximacionPi n =

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

signoPrimo :: Int -> Int

signoPrimo 2 = 1

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

| otherwise = -1

signo :: Int -> Int

signo n | isPrime n = signoPrimo n

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

-- 2ª definición

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

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = serieEuler !! (n-1)

serieEuler :: [Double]

serieEuler =

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

-- Definición de grafica

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

grafica :: Int -> IO ()

grafica n =

plotList [Key Nothing]

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

Regresión lineal

PorJosé A. Alonso 20-febrero-201726-marzo-2022

Dadas dos listas de valores

   xs = [x(1), x(2), ..., x(n)]
   ys = [y(1), y(2), ..., y(n)]

1 2	xs = [x(1), x(2), ..., x(n)] ys = [y(1), y(2), ..., y(n)]

la ecuación de la recta de regresión de ys sobre xs es y = a+bx, donde

   b = (nΣx(i)y(i) - Σx(i)Σy(i)) / (nΣx(i)² - (Σx(i))²)
   a = (Σy(i) - bΣx(i)) / n

1 2	b = (nΣx(i)y(i) - Σx(i)Σy(i)) / (nΣx(i)² - (Σx(i))²) a = (Σy(i) - bΣx(i)) / n

Definir la función

   regresionLineal :: [Double] -> [Double] -> (Double,Double)

1	regresionLineal :: [Double] -> [Double] -> (Double,Double)

tal que (regresionLineal xs ys) es el par (a,b) de los coeficientes de la recta de regresión de ys sobre xs. Por ejemplo, para los valores

   ejX, ejY :: [Double]
   ejX = [5,  7, 10, 12, 16, 20, 23, 27, 19, 14]
   ejY = [9, 11, 15, 16, 20, 24, 27, 29, 22, 20]

ejX, ejY :: [Double]

ejX = [5, 7, 10, 12, 16, 20, 23, 27, 19, 14]

ejY = [9, 11, 15, 16, 20, 24, 27, 29, 22, 20]

se tiene

   λ> regresionLineal ejX ejY
   (5.195045748716805,0.9218924347243919)

1 2	λ> regresionLineal ejX ejY (5.195045748716805,0.9218924347243919)

Para comprobar la definición, se importa la librería Graphics.Gnuplot.Simple y se define el procedimiento

   grafica :: [Double] -> [Double] -> IO ()
   grafica xs ys = 
       plotPathsStyle 
         [YRange (0,10+mY)] 
         [(defaultStyle {plotType = Points,
                         lineSpec = CustomStyle [LineTitle "Datos",
                                                 PointType 2,
                                                 PointSize 2.5]},
                        zip xs ys),
          (defaultStyle {plotType = Lines,
                         lineSpec = CustomStyle [LineTitle "Ajuste",
                                                 LineWidth 2]},
                        [(x,a+b*x) | x <- [0..mX]])]
       where (a,b) = regresionLineal xs ys                      
             mX    = maximum xs
             mY    = maximum ys

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

grafica xs ys =

plotPathsStyle

[YRange (0,10+mY)]

[(defaultStyle {plotType = Points,

lineSpec = CustomStyle [LineTitle "Datos",

PointType 2,

PointSize 2.5]},

zip xs ys),

(defaultStyle {plotType = Lines,

lineSpec = CustomStyle [LineTitle "Ajuste",

LineWidth 2]},

[(x,a+b*x) | x <- [0..mX]])]

where (a,b) = regresionLineal xs ys

mX = maximum xs

mY = maximum ys

tal que (grafica xs ys) pinta los puntos correspondientes a las listas de valores xs e ys y su recta de regresión. Por ejemplo, con (grafica ejX ejY) se obtiene el siguiente dibujo

Soluciones

import Data.List (genericLength)
import Graphics.Gnuplot.Simple

ejX, ejY :: [Double]
ejX = [5,  7, 10, 12, 16, 20, 23, 27, 19, 14]
ejY = [9, 11, 15, 16, 20, 24, 27, 29, 22, 20]

regresionLineal :: [Double] -> [Double] -> (Double,Double)
regresionLineal xs ys = (a,b)
    where n     = genericLength xs
          sumX  = sum xs
          sumY  = sum ys
          sumX2 = sum (zipWith (*) xs xs)
          sumY2 = sum (zipWith (*) ys ys)
          sumXY = sum (zipWith (*) xs ys)
          b     = (n*sumXY - sumX*sumY) / (n*sumX2 - sumX^2)
          a     = (sumY - b*sumX) / n

grafica :: [Double] -> [Double] -> IO ()
grafica xs ys = 
    plotPathsStyle 
      [YRange (0,10+mY)] 
      [(defaultStyle {plotType = Points,
                      lineSpec = CustomStyle [LineTitle "Datos",
                                              PointType 2,
                                              PointSize 2.5]},
                     zip xs ys),
       (defaultStyle {plotType = Lines,
                      lineSpec = CustomStyle [LineTitle "Ajuste",
                                              LineWidth 2]},
                     [(x,a+b*x) | x <- [0..mX]])]
    where (a,b) = regresionLineal xs ys                      
          mX    = maximum xs
          mY    = maximum ys

import Data.List (genericLength)

import Graphics.Gnuplot.Simple

ejX, ejY :: [Double]

ejX = [5, 7, 10, 12, 16, 20, 23, 27, 19, 14]

ejY = [9, 11, 15, 16, 20, 24, 27, 29, 22, 20]

regresionLineal :: [Double] -> [Double] -> (Double,Double)

regresionLineal xs ys = (a,b)

where n = genericLength xs

sumX = sum xs

sumY = sum ys

sumX2 = sum (zipWith (*) xs xs)

sumY2 = sum (zipWith (*) ys ys)

sumXY = sum (zipWith (*) xs ys)

b = (n*sumXY - sumX*sumY) / (n*sumX2 - sumX^2)

a = (sumY - b*sumX) / n

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

grafica xs ys =

plotPathsStyle

[YRange (0,10+mY)]

[(defaultStyle {plotType = Points,

lineSpec = CustomStyle [LineTitle "Datos",

PointType 2,

PointSize 2.5]},

zip xs ys),

(defaultStyle {plotType = Lines,

lineSpec = CustomStyle [LineTitle "Ajuste",

LineWidth 2]},

[(x,a+b*x) | x <- [0..mX]])]

where (a,b) = regresionLineal xs ys

mX = maximum xs

mY = maximum ys

Regresión lineal

PorJosé A. Alonso 22-febrero-201626-marzo-2022

Dadas dos listas de valores

   xs = [x(1), x(2), ..., x(n)]
   ys = [y(1), y(2), ..., y(n)]

1 2	xs = [x(1), x(2), ..., x(n)] ys = [y(1), y(2), ..., y(n)]

la ecuación de la recta de regresión de ys sobre xs es y = a+bx, donde

   b = (nΣx(i)y(i) - Σx(i)Σy(i)) / (nΣx(i)² - (Σx(i))²)
   a = (Σy(i) - bΣx(i)) / n

1 2	b = (nΣx(i)y(i) - Σx(i)Σy(i)) / (nΣx(i)² - (Σx(i))²) a = (Σy(i) - bΣx(i)) / n

Definir la función

   regresionLineal :: [Double] -> [Double] -> (Double,Double)

1	regresionLineal :: [Double] -> [Double] -> (Double,Double)

tal que (regresionLineal xs ys) es el par (a,b) de los coeficientes de la recta de regresión de ys sobre xs. Por ejemplo, para los valores

   ejX, ejY :: [Double]
   ejX = [5,  7, 10, 12, 16, 20, 23, 27, 19, 14]
   ejY = [9, 11, 15, 16, 20, 24, 27, 29, 22, 20]

ejX, ejY :: [Double]

ejX = [5, 7, 10, 12, 16, 20, 23, 27, 19, 14]

ejY = [9, 11, 15, 16, 20, 24, 27, 29, 22, 20]

se tiene

   λ> regresionLineal ejX ejY
   (5.195045748716805,0.9218924347243919)

1 2	λ> regresionLineal ejX ejY (5.195045748716805,0.9218924347243919)

Definir el procedimiento

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

1	grafica :: [Double] -> [Double] -> IO ()

tal que (grafica xs ys) pinte los puntos correspondientes a las listas de valores xs e ys y su recta de regresión. Por ejemplo, con (grafica ejX ejY) se obtiene el siguiente dibujo

Soluciones

import Graphics.Gnuplot.Simple

ejX, ejY :: [Double]
ejX = [5,  7, 10, 12, 16, 20, 23, 27, 19, 14]
ejY = [9, 11, 15, 16, 20, 24, 27, 29, 22, 20]

regresionLineal :: [Double] -> [Double] -> (Double,Double)
regresionLineal xs ys = (a,b)
    where n     = fromIntegral (length xs)
          sumX  = sum xs
          sumY  = sum ys
          sumX2 = sum (zipWith (*) xs xs)
          sumY2 = sum (zipWith (*) ys ys)
          sumXY = sum (zipWith (*) xs ys)
          b     = (n*sumXY - sumX*sumY) / (n*sumX2 - sumX^2)
          a     = (sumY - b*sumX) / n

grafica :: [Double] -> [Double] -> IO ()
grafica xs ys = 
    plotPathsStyle 
      [YRange (0,10+mY)] 
      [(defaultStyle {plotType = Points,
                      lineSpec = CustomStyle [LineTitle "Datos",
                                              PointType 2,
                                              PointSize 2.5]},
                     zip xs ys),
       (defaultStyle {plotType = Lines,
                      lineSpec = CustomStyle [LineTitle "Ajuste",
                                              LineWidth 2]},
                     [(x,a+b*x) | x <- [0..mX]])]
    where (a,b) = regresionLineal xs ys                      
          mX    = maximum xs 
          mY    = maximum ys

import Graphics.Gnuplot.Simple

ejX, ejY :: [Double]

ejX = [5, 7, 10, 12, 16, 20, 23, 27, 19, 14]

ejY = [9, 11, 15, 16, 20, 24, 27, 29, 22, 20]

regresionLineal :: [Double] -> [Double] -> (Double,Double)

regresionLineal xs ys = (a,b)

where n = fromIntegral (length xs)

sumX = sum xs

sumY = sum ys

sumX2 = sum (zipWith (*) xs xs)

sumY2 = sum (zipWith (*) ys ys)

sumXY = sum (zipWith (*) xs ys)

b = (n*sumXY - sumX*sumY) / (n*sumX2 - sumX^2)

a = (sumY - b*sumX) / n

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

grafica xs ys =

plotPathsStyle

[YRange (0,10+mY)]

[(defaultStyle {plotType = Points,

lineSpec = CustomStyle [LineTitle "Datos",

PointType 2,

PointSize 2.5]},

zip xs ys),

(defaultStyle {plotType = Lines,

lineSpec = CustomStyle [LineTitle "Ajuste",

LineWidth 2]},

[(x,a+b*x) | x <- [0..mX]])]

where (a,b) = regresionLineal xs ys

mX = maximum xs

mY = maximum ys

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Entradas recientes

Buscar

Correo electrónico