iterate

Números sin 2 en base 3

PorJosé A. Alonso 22-enero-202029-enero-2020

Definir la sucesión

   numerosSin2EnBase3 :: [Integer]

1	numerosSin2EnBase3 :: [Integer]

cuyos términos son los números cuya representación en base 3 no contiene el dígito 2. Por ejemplo,

   λ> take 20 numerosSin2EnBase3
   [0,1,3,4,9,10,12,13,27,28,30,31,36,37,39,40,81,82,84,85]

1 2	λ> take 20 numerosSin2EnBase3 [0,1,3,4,9,10,12,13,27,28,30,31,36,37,39,40,81,82,84,85]

Se observa que

12 está en la sucesión ya que su representación en base 3 es 110 (porque 1·3² + 1·3¹ + 0.3⁰ = 12) y no contiene a 2.
14 no está en la sucesión ya que su representación en base 3 es 112 (porque 1·3² + 1·3¹ + 2.3⁰ = 14) y contiene a 2.

Comprobar con QuickCheck que las sucesiones numerosSin2EnBase3 y sucesionSin3enPA (del ejercicio anterior) son iguales; es decir, para todo número natural n, el n-ésimo término de numerosSin2EnBase3 es igual al n-ésimo término de sucesionSin3enPA.

Soluciones

import Data.List ((\\))
import Test.QuickCheck

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

numerosSin2EnBase3a :: [Integer]
numerosSin2EnBase3a =
  [n | n <- [0..]
     , 2 `notElem` (enBase3 n)]

-- (enBase3 n) es la representación de n en base 3. Por ejemplo,
--    enBase3 7   ==  [1,2]
--    enBase3 9   ==  [0,0,1]
--    enBase3 10  ==  [1,0,1]
--    enBase3 11  ==  [2,0,1]
enBase3 :: Integer -> [Integer]
enBase3 n | n < 3     = [n]
          | otherwise = r : enBase3 q
  where (q,r) = quotRem n 3

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

-- Se puede construir como un triángulo:
--    0
--    1
--    3 4
--    9 10 12 13
--    27 28 30 31 36 37 39 40
--    ....

numerosSin2EnBase3b :: [Integer]
numerosSin2EnBase3b = 0 : concat (iterate siguientes [1])
  where siguientes xs = concatMap (\x -> [3*x,3*x+1]) xs

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

-- La compración es
--    λ> numerosSin2EnBase3a !! (10^4)
--    1679697
--    (4.06 secs, 2,245,415,808 bytes)
--    λ> numerosSin2EnBase3b !! (10^4)
--    1679697
--    (0.03 secs, 2,109,912 bytes)

-- Definición
-- ==========

-- En lo que sigue usaremos la 2ª definición.
numerosSin2EnBase3 :: [Integer]
numerosSin2EnBase3 = numerosSin2EnBase3b

-- Propiedad
-- =========

-- La propiedad es
prop_equivalencia :: Int -> Property
prop_equivalencia n =
  n > 0 ==> sucesionSin3enPA !! n == numerosSin2EnBase3 !! n

-- sucesionSin3enPA donde cada uno de sus términos es el menor número
-- natural tal que no está en PA con cualesquiera dos términos
-- anteriores de la sucesión. 
sucesionSin3enPA :: [Integer]
sucesionSin3enPA = aux [] [0..] 
  where
    aux xs (y:ys) = y : aux (y:xs) (ys \\ [2 * y - x | x <- xs])

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

import Data.List ((\\))

import Test.QuickCheck

-- 1ª solución

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

numerosSin2EnBase3a :: [Integer]

numerosSin2EnBase3a =

[n | n <- [0..]

, 2 `notElem` (enBase3 n)]

-- (enBase3 n) es la representación de n en base 3. Por ejemplo,

-- enBase3 7 == [1,2]

-- enBase3 9 == [0,0,1]

-- enBase3 10 == [1,0,1]

-- enBase3 11 == [2,0,1]

enBase3 :: Integer -> [Integer]

enBase3 n | n < 3 = [n]

| otherwise = r : enBase3 q

where (q,r) = quotRem n 3

-- 2ª definición

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

-- Se puede construir como un triángulo:

-- 0

-- 1

-- 3 4

-- 9 10 12 13

-- 27 28 30 31 36 37 39 40

-- ....

numerosSin2EnBase3b :: [Integer]

numerosSin2EnBase3b = 0 : concat (iterate siguientes [1])

where siguientes xs = concatMap (\x -> [3*x,3*x+1]) xs

-- Comparación de eficiencia

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

-- La compración es

-- λ> numerosSin2EnBase3a !! (10^4)

-- 1679697

-- (4.06 secs, 2,245,415,808 bytes)

-- λ> numerosSin2EnBase3b !! (10^4)

-- 1679697

-- (0.03 secs, 2,109,912 bytes)

-- Definición

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

-- En lo que sigue usaremos la 2ª definición.

numerosSin2EnBase3 :: [Integer]

numerosSin2EnBase3 = numerosSin2EnBase3b

-- Propiedad

-- =========

-- La propiedad es

prop_equivalencia :: Int -> Property

prop_equivalencia n =

n > 0 ==> sucesionSin3enPA !! n == numerosSin2EnBase3 !! n

-- sucesionSin3enPA donde cada uno de sus términos es el menor número

-- natural tal que no está en PA con cualesquiera dos términos

-- anteriores de la sucesión.

sucesionSin3enPA :: [Integer]

sucesionSin3enPA = aux [] [0..]

where

aux xs (y:ys) = y : aux (y:xs) (ys \\ [2 * y - x | x <- xs])

-- La comprobación es

-- λ> quickCheck prop_equivalencia

-- +++ OK, passed 100 tests.

Otras soluciones

Se pueden escribir otras soluciones en los comentarios.
El código se debe escribir entre una línea con <pre lang=»haskell»> y otra con </pre>

Pensamiento

O que yo pueda asesinar un día
en mi alma, al despertar, esa persona
que me hizo el mundo mientras yo dormía.

Antonio Machado

La mayor potencia de dos no es divisor

PorJosé A. Alonso 7-enero-202014-enero-2020

Para cada número entero positivo n, se define el conjunto

   S(n) = {1, 2, 3, ..., n}

1	S(n) = {1, 2, 3, ..., n}

de los números desde el 1 hasta n.

Definir la función

   mayorPotenciaDeDosEnS :: Integer -> Integer

1	mayorPotenciaDeDosEnS :: Integer -> Integer

tal que (mayorPotenciaDeDosEnS n) es la mayor potencia de 2 en S(n). Por ejemplo,

   mayorPotenciaDeDosEnS 3  ==  2
   mayorPotenciaDeDosEnS 4  ==  4

1 2	mayorPotenciaDeDosEnS 3 == 2 mayorPotenciaDeDosEnS 4 == 4

Comprobar con QuickCheck que la mayor potencia de 2 en S(n) no divide a ningún otro elemento de S(n).

Soluciones

import Data.List (delete)
import Test.QuickCheck

mayorPotenciaDeDosEnS :: Integer -> Integer
mayorPotenciaDeDosEnS n =
  last (takeWhile (<=n) potenciasDeDos)

-- potenciasDeDos es la lista de las potencias de 2. Por ejemplo,
--    take 5 potenciasDeDos  ==  [1,2,4,8,16]
potenciasDeDos :: [Integer]
potenciasDeDos = iterate (*2) 1

-- La propiedad es
prop_mayorPotenciaDeDosEnS :: Integer -> Property
prop_mayorPotenciaDeDosEnS n =
  n > 0 ==> all (x `noDivide`) (delete x [1..n])
  where x = mayorPotenciaDeDosEnS n
        x `noDivide` y = y `mod` x /= 0

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

import Data.List (delete)

import Test.QuickCheck

mayorPotenciaDeDosEnS :: Integer -> Integer

mayorPotenciaDeDosEnS n =

last (takeWhile (<=n) potenciasDeDos)

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

-- take 5 potenciasDeDos == [1,2,4,8,16]

potenciasDeDos :: [Integer]

potenciasDeDos = iterate (*2) 1

-- La propiedad es

prop_mayorPotenciaDeDosEnS :: Integer -> Property

prop_mayorPotenciaDeDosEnS n =

n > 0 ==> all (x `noDivide`) (delete x [1..n])

where x = mayorPotenciaDeDosEnS n

x `noDivide` y = y `mod` x /= 0

-- La comprobación es

-- λ> quickCheck prop_mayorPotenciaDeDosEnS

-- +++ OK, passed 100 tests.

Referencia

Basado en Greatest power of two not divisor de ProofWiki.

Pensamiento

¡Sólo tu figura,
como una centella blanca,
en mi noche oscura.

Antonio Machado

Cálculo de pi usando la fórmula de Vieta

PorJosé A. Alonso 16-diciembre-201923-diciembre-2019

La fórmula de Vieta para el cálculo de pi es la siguiente

Definir las funciones

   aproximacionPi :: Int -> Double
   errorPi :: Double -> Int

1 2	aproximacionPi :: Int -> Double errorPi :: Double -> Int

tales que

(aproximacionPi n) es la aproximación de pi usando n factores de la fórmula de Vieta. Por ejemplo,

     aproximacionPi  5  ==  3.140331156954753
     aproximacionPi 10  ==  3.1415914215112
     aproximacionPi 15  ==  3.141592652386592
     aproximacionPi 20  ==  3.1415926535886207
     aproximacionPi 25  ==  3.141592653589795

aproximacionPi 5 == 3.140331156954753

aproximacionPi 10 == 3.1415914215112

aproximacionPi 15 == 3.141592652386592

aproximacionPi 20 == 3.1415926535886207

aproximacionPi 25 == 3.141592653589795

(errorPi x) es el menor número de factores de la fórmula de Vieta necesarios para obtener pi con un error menor que x. Por ejemplo,

     errorPi 0.1        ==  2
     errorPi 0.01       ==  4
     errorPi 0.001      ==  6
     errorPi 0.0001     ==  7
     errorPi 1e-4       ==  7
     errorPi 1e-14      ==  24
     pi                 ==  3.141592653589793
     aproximacionPi 24  ==  3.1415926535897913

errorPi 0.1 == 2

errorPi 0.01 == 4

errorPi 0.001 == 6

errorPi 0.0001 == 7

errorPi 1e-4 == 7

errorPi 1e-14 == 24

pi == 3.141592653589793

aproximacionPi 24 == 3.1415926535897913

Soluciones

-- 1ª definición de aproximacionPi
aproximacionPi :: Int -> Double
aproximacionPi n = product [2 / aux x | x <- [0..n]]
  where
    aux 0 = 1
    aux 1 = sqrt 2
    aux n = sqrt (2 + aux (n-1))

-- 2ª definición de aproximacionPi
aproximacionPi2 :: Int -> Double
aproximacionPi2 n = product [2/x | x <- 1 : xs] 
  where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)

-- 3ª definición de aproximaxionPi
aproximacionPi3 :: Int -> Double
aproximacionPi3 n =  product (2 : take n (map (2/) xs))
  where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]

-- 1ª definición de errorPi
errorPi :: Double -> Int
errorPi x = head [n | n <- [1..]
                    , abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi
errorPi2 :: Double -> Int
errorPi2 x = until aceptable (+1) 1
  where aceptable n = abs (pi - aproximacionPi n) < x

-- 1ª definición de aproximacionPi

aproximacionPi :: Int -> Double

aproximacionPi n = product [2 / aux x | x <- [0..n]]

where

aux 0 = 1

aux 1 = sqrt 2

aux n = sqrt (2 + aux (n-1))

-- 2ª definición de aproximacionPi

aproximacionPi2 :: Int -> Double

aproximacionPi2 n = product [2/x | x <- 1 : xs]

where xs = take n $ iterate (\x -> sqrt (2+x)) (sqrt 2)

-- 3ª definición de aproximaxionPi

aproximacionPi3 :: Int -> Double

aproximacionPi3 n = product (2 : take n (map (2/) xs))

where xs = sqrt 2 : [sqrt (2 + x) | x <- xs]

-- 1ª definición de errorPi

errorPi :: Double -> Int

errorPi x = head [n | n <- [1..]

, abs (pi - aproximacionPi n) < x]

-- 2ª definición de errorPi

errorPi2 :: Double -> Int

errorPi2 x = until aceptable (+1) 1

where aceptable n = abs (pi - aproximacionPi n) < x

Pensamiento

El tiempo que la barba me platea,
cavó mis ojos y agrandó mi frente,
va siendo en mi recuerdo transparente,
y mientras más al fondo, más clarea.

Antonio Machado

Las sucesiones de Loomis

PorJosé A. Alonso 15-mayo-201925-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

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]

Pensamiento

Era una noche del mes
de mayo, azul y serena.
Sobre el agudo ciprés
brillaba la luna llena.

Antonio Machado

La conjetura de Collatz

PorJosé A. Alonso 9-abril-201926-marzo-2022

La conjetura de Collatz, conocida también como conjetura 3n+1, fue enunciada por Lothar Collatz en 1937 y, hasta la fecha, no se ha resuelto.

La conjetura hace referencia a una propiedad de las sucesiones de Siracusa. La sucesión de Siracusa de un número entero positivo x es la sucesión cuyo primer término es x y el siguiente de un término se obtiene dividiéndolo entre 2, si es par o multiplicándolo por 3 y sumándole 1, si es impar. Por ejemplo, la sucesión de Siracusa de 12 es

   12, 6, 3, 10, 5, 16, 8, 4, 2, 1, 4, 2, 1, 4, 2, 1, 4, 2, 1, 4, ....

1	12, 6, 3, 10, 5, 16, 8, 4, 2, 1, 4, 2, 1, 4, 2, 1, 4, 2, 1, 4, ....

La conjetura de Collatz afirma que para todo número entero positivo x, el 1 pertenece a la sucesión de Siracusa de x.

Definir las funciones

   siracusa        :: Integer -> [Integer]
   graficaSiracusa :: Int -> [Integer] -> IO ()

1 2	siracusa :: Integer -> [Integer] graficaSiracusa :: Int -> [Integer] -> IO ()

tales que

(siracusa x) es la sucesión de Siracusa de x. Por ejemplo,

     λ> take 20 (siracusa 12)
     [12,6,3,10,5,16,8,4,2,1,4,2,1,4,2,1,4,2,1,4]
     λ> take 20 (siracusa 22)
     [22,11,34,17,52,26,13,40,20,10,5,16,8,4,2,1,4,2,1,4]

λ> take 20 (siracusa 12)

[12,6,3,10,5,16,8,4,2,1,4,2,1,4,2,1,4,2,1,4]

λ> take 20 (siracusa 22)

[22,11,34,17,52,26,13,40,20,10,5,16,8,4,2,1,4,2,1,4]

(graficaSiracusa n xs) dibuja los n primeros términos de las sucesiones de Siracusas de los elementos de xs. Por ejemplo, (graficaSiracusa 100 [27]) dibuja

y (graficaSiracusa 150 [1..1000]) dibuja

Comprobar con QuickCheck la conjetura de Collatz.

Soluciones

import Test.QuickCheck
import Graphics.Gnuplot.Simple

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

siracusa :: Integer -> [Integer]
siracusa n | even n    = n : siracusa (n `div` 2)
           | otherwise = n : siracusa (3*n+1)

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

siracusa2 :: Integer -> [Integer]
siracusa2 = iterate siguiente 
  where siguiente x | even x    = x `div` 2
                    | otherwise = 3*x+1

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

-- La comparación es
--    λ> siracusa 12 !! (10^6)
--    4
--    (0.55 secs, 362,791,008 bytes)
--    λ> siracusa 12 !! (2*10^6)
--    2
--    (1.05 secs, 725,456,376 bytes)
--    λ> siracusa2 12 !! (10^6)
--    4
--    (1.66 secs, 647,189,664 bytes)
--    λ> siracusa2 12 !! (2*10^6)
--    2
--    (3.11 secs, 1,294,286,792 bytes)

-- Definición de graficaSiracusa
-- =============================

graficaSiracusa :: Int -> [Integer] -> IO ()
graficaSiracusa n xs =
  plotLists [ Key Nothing
            , PNG "La_conjetura_de_Collatz.png"
            ]
            [take n (siracusa x) | x <- xs]

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

-- La conjetura es
conjeturaCollatz :: Positive Integer -> Bool
conjeturaCollatz (Positive n) =
  1 `elem` siracusa n

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

import Test.QuickCheck

import Graphics.Gnuplot.Simple

-- 1ª definición de siracusa

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

siracusa :: Integer -> [Integer]

siracusa n | even n = n : siracusa (n `div` 2)

| otherwise = n : siracusa (3*n+1)

-- 2ª definición de siracusa

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

siracusa2 :: Integer -> [Integer]

siracusa2 = iterate siguiente

where siguiente x | even x = x `div` 2

| otherwise = 3*x+1

-- Comparación de eficiencia

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

-- La comparación es

-- λ> siracusa 12 !! (10^6)

-- 4

-- (0.55 secs, 362,791,008 bytes)

-- λ> siracusa 12 !! (2*10^6)

-- 2

-- (1.05 secs, 725,456,376 bytes)

-- λ> siracusa2 12 !! (10^6)

-- 4

-- (1.66 secs, 647,189,664 bytes)

-- λ> siracusa2 12 !! (2*10^6)

-- 2

-- (3.11 secs, 1,294,286,792 bytes)

-- Definición de graficaSiracusa

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

graficaSiracusa :: Int -> [Integer] -> IO ()

graficaSiracusa n xs =

plotLists [ Key Nothing

, PNG "La_conjetura_de_Collatz.png"

]

[take n (siracusa x) | x <- xs]

-- Comprobación de la conjetura

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

-- La conjetura es

conjeturaCollatz :: Positive Integer -> Bool

conjeturaCollatz (Positive n) =

1 `elem` siracusa n

-- La comprobación es

-- λ> quickCheck conjeturaCollatz

-- +++ OK, passed 100 tests.

Pensamiento

Que el caminante es suma del camino …

Antonio Machado

Números sin 2 en base 3

Soluciones

Otras soluciones

Pensamiento

La mayor potencia de dos no es divisor

Soluciones

Referencia

Pensamiento

Cálculo de pi usando la fórmula de Vieta

Soluciones

Pensamiento

Las sucesiones de Loomis

Soluciones

Pensamiento

La conjetura de Collatz

Soluciones

Pensamiento

Entradas recientes

Buscar

Correo electrónico