isPrime

La conjetura de Levy

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

    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.

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

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

Números somirp

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

Un número omirp es un número primo que forma un primo distinto al invertir el orden de sus dígitos.

Definir las funciones

  esOmirp            :: Integer -> Bool
  omirps             :: [Integer]
  nOmirpsIntermedios :: Int -> Int

esOmirp :: Integer -> Bool

omirps :: [Integer]

nOmirpsIntermedios :: Int -> Int

tales que

(esOmirp n) se verifica si n es un número omirp. Por ejemplo,

     esOmirp 13      ==  True
     esOmirp 11      ==  False
     esOmirp 112207  ==  True

esOmirp 13 == True

esOmirp 11 == False

esOmirp 112207 == True

omirps es la lista de los números omirps. Por ejemplo,

     take 15 omirps  ==  [13,17,31,37,71,73,79,97,107,113]
     omirps !! 2000  ==  112207

1 2	take 15 omirps == [13,17,31,37,71,73,79,97,107,113] omirps !! 2000 == 112207

(nOmirpsIntermedios n) es la cantidad de números omirps entre el n-ésimo número omirp y el obtenido al invertir el orden de sus dígitos. Por ejemplo,

     nOmirpsIntermedios 2000  ==  4750

1	nOmirpsIntermedios 2000 == 4750

Nota: Este ejercicio ha sido propuesto por Ángel Ruiz Campos.

Soluciones

import Data.Numbers.Primes (isPrime,primes)

esOmirp :: Integer -> Bool
esOmirp n = n /= rn && isPrime rn
  where rn = read . reverse . show $ n

omirps :: [Integer]
omirps = filter esOmirp primes

nOmirpsIntermedios :: Int -> Int
nOmirpsIntermedios n =
  length
  . filter esOmirp
  . takeWhile (< rx)
  . dropWhile (<= x) $ primes
  where x  = omirps !! n
        rx = read . reverse . show $ x

import Data.Numbers.Primes (isPrime,primes)

esOmirp :: Integer -> Bool

esOmirp n = n /= rn && isPrime rn

where rn = read . reverse . show $ n

omirps :: [Integer]

omirps = filter esOmirp primes

nOmirpsIntermedios :: Int -> Int

nOmirpsIntermedios n =

length

. filter esOmirp

. takeWhile (< rx)

. dropWhile (<= x) $ primes

where x = omirps !! n

rx = read . reverse . show $ x

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

La conjetura de Levy

Soluciones

Sucesión de Lichtenberg

Soluciones

Números somirp

Soluciones

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

Soluciones

Entradas recientes

Buscar

Correo electrónico