sum

Cálculo de pi mediante el método de Newton

PorJosé A. Alonso 18-marzo-202025-marzo-2020

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

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]

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

«Mi trabajo siempre trató de unir lo verdadero con lo bello; pero cuando tuve que elegir uno u otro, generalmente elegí lo bello.»

Hermann Weyl.

Medias de dígitos de pi

PorJosé A. Alonso 16-marzo-202023-marzo-2020

El fichero Digitos_de_pi.txt contiene el número pi con un millón de decimales; es decir,

   3.1415926535897932384626433832 ... 83996346460422090106105779458151

1	3.1415926535897932384626433832 ... 83996346460422090106105779458151

Definir las funciones

   mediasDigitosDePi        :: IO [Double]
   graficaMediasDigitosDePi :: Int -> IO ()

1 2	mediasDigitosDePi :: IO [Double] graficaMediasDigitosDePi :: Int -> IO ()

tales que

mediasDigitosDePi es la sucesión cuyo n-ésimo elemento es la media de los n primeros dígitos de pi. Por ejemplo,

     λ> xs <- mediasDigitosDePi
     λ> take 5 xs
     [1.0,2.5,2.0,2.75,4.0]
     λ> take 10 xs
     [1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]
     λ> take 10 <$> mediasDigitosDePi
     [1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]

λ> xs <- mediasDigitosDePi

λ> take 5 xs

[1.0,2.5,2.0,2.75,4.0]

λ> take 10 xs

[1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]

λ> take 10 <$> mediasDigitosDePi

[1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]

(graficaMediasDigitosDePi n) dibuja la gráfica de los n primeros términos de mediasDigitosDePi. Por ejemplo,
- (graficaMediasDigitosDePi 20) dibuja
- (graficaMediasDigitosDePi 200) dibuja
- (graficaMediasDigitosDePi 2000) dibuja

Soluciones

import Data.Char (digitToInt)
import Data.List (genericLength, inits, tails)
import Graphics.Gnuplot.Simple ( Attribute (Key, PNG)
                               , plotList )

-- Definición de mediasDigitosDePi
-- ===============================

mediasDigitosDePi :: IO [Double]
mediasDigitosDePi = do
  (_:_:ds) <- readFile "Digitos_de_pi.txt"
  return (medias (digitos ds))

-- (digitos cs) es la lista de los digitos de cs. Por ejempplo,
--    digitos "1415926535"  ==  [1,4,1,5,9,2,6,5,3,5]
digitos :: String -> [Int]
digitos = map digitToInt

-- (medias xs) es la lista de las medias de los segmentos iniciales de
-- xs. Por ejemplo,
--    λ> medias [1,4,1,5,9,2,6,5,3,5]
--    [1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]
medias :: [Int] -> [Double]
medias xs = map media (tail (inits xs))

-- (media xs) es la media aritmética de xs. Por ejemplo,
--    media [1,4,1,5,9]  ==  4.0
media :: [Int] -> Double
media xs = fromIntegral (sum xs) / genericLength xs

-- Definición de graficaMediasDigitosDePi
-- ======================================

graficaMediasDigitosDePi :: Int -> IO ()
graficaMediasDigitosDePi n = do
  xs <- mediasDigitosDePi
  plotList [ Key Nothing
           , PNG ("Medias_de_digitos_de_pi_" ++ show n ++ ".png")
           ]
           (take n xs)

import Data.Char (digitToInt)

import Data.List (genericLength, inits, tails)

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

, plotList )

-- Definición de mediasDigitosDePi

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

mediasDigitosDePi :: IO [Double]

mediasDigitosDePi = do

(_:_:ds) <- readFile "Digitos_de_pi.txt"

return (medias (digitos ds))

-- (digitos cs) es la lista de los digitos de cs. Por ejempplo,

-- digitos "1415926535" == [1,4,1,5,9,2,6,5,3,5]

digitos :: String -> [Int]

digitos = map digitToInt

-- (medias xs) es la lista de las medias de los segmentos iniciales de

-- xs. Por ejemplo,

-- λ> medias [1,4,1,5,9,2,6,5,3,5]

-- [1.0,2.5,2.0,2.75,4.0,3.6666666666666665,4.0,4.125,4.0,4.1]

medias :: [Int] -> [Double]

medias xs = map media (tail (inits xs))

-- (media xs) es la media aritmética de xs. Por ejemplo,

-- media [1,4,1,5,9] == 4.0

media :: [Int] -> Double

media xs = fromIntegral (sum xs) / genericLength xs

-- Definición de graficaMediasDigitosDePi

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

graficaMediasDigitosDePi :: Int -> IO ()

graficaMediasDigitosDePi n = do

xs <- mediasDigitosDePi

plotList [ Key Nothing

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

]

(take n xs)

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

Es el mejor de los buenos
quien sabe que en esta vida
todo es cuestión de medida:
un poco más, algo menos.

Antonio Machado

Pandemia

PorJosé A. Alonso 2-marzo-20209-marzo-2020

¡El mundo está en cuarentena! Hay una nueva pandemia que lucha contra la humanidad. Cada continente está aislado de los demás, pero las personas infectadas se han propagado antes de la advertencia.

En este problema se representará el mundo por una cadena como la siguiente

   "01000000X000X011X0X"

1	"01000000X000X011X0X"

donde 0 representa no infectado, 1 representa infectado y X representa un océano

Las reglas de propagación son:

El virus no puede propagarse al otro lado de un océano.
Si una persona se infecta, todas las personas de este continente se infectan también.
El primer y el último continente no están conectados.

El problema consiste en encontrar el porcentaje de la población humana que se infectó al final. Por ejemplo,

   inicio:     "01000000X000X011X0X"
   final:      "11111111X000X111X0X"
   total:      15
   infectados: 11
   porcentaje: 100*11/15 = 73.33333333333333

inicio: "01000000X000X011X0X"

final: "11111111X000X111X0X"

total: 15

infectados: 11

porcentaje: 100*11/15 = 73.33333333333333

Definir la función

   porcentajeInfectados :: String -> Double

1	porcentajeInfectados :: String -> Double

tal que (porcentajeInfectados xs) es el porcentaje final de infectados para el mapa inicial xs. Por ejemplo,

   porcentajeInfectados "01000000X000X011X0X"  == 73.33333333333333
   porcentajeInfectados "01X000X010X011XX"     == 72.72727272727273
   porcentajeInfectados "XXXXX"                == 0.0
   porcentajeInfectados "0000000010"           == 100.0
   porcentajeInfectados "X00X000000X10X0100"   == 42.857142857142854

porcentajeInfectados "01000000X000X011X0X" == 73.33333333333333

porcentajeInfectados "01X000X010X011XX" == 72.72727272727273

porcentajeInfectados "XXXXX" == 0.0

porcentajeInfectados "0000000010" == 100.0

porcentajeInfectados "X00X000000X10X0100" == 42.857142857142854

Soluciones

import Data.List (genericLength)
import Data.List.Split (splitOn)

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

porcentajeInfectados :: String -> Double
porcentajeInfectados xs
  | nh == 0   = 0
  | otherwise = 100 * ni / nh
  where ni = fromIntegral (numeroInfectados xs)
        nh = fromIntegral (numeroHabitantes xs)

-- (continentes xs) es la lista de las poblaciones de los continentes
-- del mapa xs. Por ejemplo,
--    continentes "01000000X000X011X0X" == ["01000000","000","011","0"]
--    continentes "01X000X010X011XX"    == ["01","000","010","011"]
--    continentes "XXXXX"               == [""]
--    continentes "0000000010"          == ["0000000010"]
--    continentes "X00X000000X10X0100"  == ["","00","000000","10","0100"]
continentes :: String -> [String]
continentes [] = []
continentes xs = as : continentes (dropWhile (=='X') bs)
  where (as,bs) = break (=='X') xs

-- (numeroInfectados xs) es el número final de infectados a partir del
-- mapa xs. Por ejemplo,
--    numeroInfectados "01000000X000X011X0X"  ==  11
numeroInfectados :: String -> Int
numeroInfectados xs =
  sum [length ys | ys <- continentes xs
                 , '1' `elem` ys]
  
-- (numeroHabitantes xs) es el número final de habitantes del mapa
-- xs. Por ejemplo, 
--    numeroHabitantes "01000000X000X011X0X"  ==  15
numeroHabitantes :: String -> Int
numeroHabitantes xs = length (filter (/='X') xs)

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

porcentajeInfectados2 :: String -> Double
porcentajeInfectados2 xs
  | nh == 0   = 0
  | otherwise = 100 * ni / nh
  where ni = sum [genericLength ys | ys <- splitOn "X" xs, '1' `elem` ys]
        nh = genericLength (filter (/='X') xs)

import Data.List (genericLength)

import Data.List.Split (splitOn)

-- 1ª solución

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

porcentajeInfectados :: String -> Double

porcentajeInfectados xs

| nh == 0 = 0

| otherwise = 100 * ni / nh

where ni = fromIntegral (numeroInfectados xs)

nh = fromIntegral (numeroHabitantes xs)

-- (continentes xs) es la lista de las poblaciones de los continentes

-- del mapa xs. Por ejemplo,

-- continentes "01000000X000X011X0X" == ["01000000","000","011","0"]

-- continentes "01X000X010X011XX" == ["01","000","010","011"]

-- continentes "XXXXX" == [""]

-- continentes "0000000010" == ["0000000010"]

-- continentes "X00X000000X10X0100" == ["","00","000000","10","0100"]

continentes :: String -> [String]

continentes [] = []

continentes xs = as : continentes (dropWhile (=='X') bs)

where (as,bs) = break (=='X') xs

-- (numeroInfectados xs) es el número final de infectados a partir del

-- mapa xs. Por ejemplo,

-- numeroInfectados "01000000X000X011X0X" == 11

numeroInfectados :: String -> Int

numeroInfectados xs =

sum [length ys | ys <- continentes xs

, '1' `elem` ys]

-- (numeroHabitantes xs) es el número final de habitantes del mapa

-- xs. Por ejemplo,

-- numeroHabitantes "01000000X000X011X0X" == 15

numeroHabitantes :: String -> Int

numeroHabitantes xs = length (filter (/='X') xs)

-- 2ª solución

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

porcentajeInfectados2 :: String -> Double

porcentajeInfectados2 xs

| nh == 0 = 0

| otherwise = 100 * ni / nh

where ni = sum [genericLength ys | ys <- splitOn "X" xs, '1' `elem` ys]

nh = genericLength (filter (/='X') xs)

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

«El avance de las matemáticas puede ser visto como un progreso de lo infinito a lo finito.»

Gian-Carlo Rota.

La conjetura de Mertens

PorJosé A. Alonso 27-enero-20203-febrero-2020

Un número entero n es libre de cuadrados si no existe un número primo p tal que p² divide a n; es decir, los factores primos de n son todos distintos.

La función de Möbius μ(n) está definida para todos los enteros positivos como sigue:

μ(n) = 1 si n es libre de cuadrados y tiene un número par de factores primos.
μ(n) = -1 si n es libre de cuadrados y tiene un número impar de factores primos.
μ(n) = 0 si n no es libre de cuadrados.

Sus primeros valores son 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, …

La función de Mertens M(n) está definida para todos los enteros positivos como la suma de μ(k) para 1 ≤ k ≤ n. Sus primeros valores son 1, 0, -1, -1, -2, -1, -2, -2, …

La conjetura de Mertens afirma que

Para todo entero x mayor que 1, el valor absoluto de la función de Mertens en x es menor que la raíz cuadrada de x.

La conjetura fue planteada por Franz Mertens en 1897. Riele Odlyzko, demostraronen 1985 que la conjetura de Mertens deja de ser cierta más o menos a partir de $10^{10^{64}}$ , cifra que luego de algunos refinamientos se redujo a $10^{10^{40}}$ .

Definir las funciones

   mobius :: Integer -> Integer
   mertens :: Integer -> Integer
   graficaMertens :: Integer -> IO ()

mobius :: Integer -> Integer

mertens :: Integer -> Integer

graficaMertens :: Integer -> IO ()

tales que

(mobius n) es el valor de la función de Möbius en n. Por ejemplo,

     mobius 6   ==  1
     mobius 30  ==  -1
     mobius 12  ==  0

mobius 6 == 1

mobius 30 == -1

mobius 12 == 0

(mertens n) es el valor de la función de Mertens en n. Por ejemplo,

     mertens 1     ==  1
     mertens 2     ==  0
     mertens 3     ==  -1
     mertens 5     ==  -2
     mertens 661   ==  -11
     mertens 1403  ==  11

mertens 1 == 1

mertens 2 == 0

mertens 3 == -1

mertens 5 == -2

mertens 661 == -11

mertens 1403 == 11

(graficaMertens n) dibuja la gráfica de la función de Mertens, la raíz cuadrada y el opuestos de la raíz cuadrada para los n primeros n enteros positivos. Por ejemplo, (graficaMertens 1000) dibuja

Comprobar con QuickCheck la conjetura de Mertens.

Nota: El ejercicio está basado en La conjetura de Merterns y su relación con un número tan raro como extremada y colosalmente grande publicado por @Alvy la semana pasada en Microsiervos.

Soluciones

import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck
import Graphics.Gnuplot.Simple

mobius :: Integer -> Integer
mobius n | tieneRepetidos xs = 0
         | otherwise         = (-1)^(length xs)
  where xs = primeFactors n

tieneRepetidos :: [Integer] -> Bool
tieneRepetidos xs =
  or [x == y | (x,y) <- zip xs (tail xs)]
        
mertens :: Integer -> Integer
mertens n = sum (map mobius [1..n])

-- Definición de graficaMertens
-- ============================

graficaMertens :: Integer -> IO ()
graficaMertens n = do
  plotLists [ Key Nothing
            , Title "Conjetura de Mertens"
            , PNG "La_conjetura_de_Mertens.png"
            ]
            [ [mertens k | k <- [1..n]]
            , raices
            , map negate raices
            ]
    
  where
    raices = [ceiling (sqrt k) | k <- [1..fromIntegral n]]

-- Conjetura de Mertens
-- ====================

-- La conjetura es
conjeturaDeMertens :: Integer -> Property
conjeturaDeMertens n =
  n > 1
  ==>
  abs (mertens n) < ceiling (sqrt n')
  where n' = fromIntegral n

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

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck

import Graphics.Gnuplot.Simple

mobius :: Integer -> Integer

mobius n | tieneRepetidos xs = 0

| otherwise = (-1)^(length xs)

where xs = primeFactors n

tieneRepetidos :: [Integer] -> Bool

tieneRepetidos xs =

or [x == y | (x,y) <- zip xs (tail xs)]

mertens :: Integer -> Integer

mertens n = sum (map mobius [1..n])

-- Definición de graficaMertens

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

graficaMertens :: Integer -> IO ()

graficaMertens n = do

plotLists [ Key Nothing

, Title "Conjetura de Mertens"

, PNG "La_conjetura_de_Mertens.png"

]

[ [mertens k | k <- [1..n]]

, raices

, map negate raices

]

where

raices = [ceiling (sqrt k) | k <- [1..fromIntegral n]]

-- Conjetura de Mertens

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

-- La conjetura es

conjeturaDeMertens :: Integer -> Property

conjeturaDeMertens n =

n > 1

==>

abs (mertens n) < ceiling (sqrt n')

where n' = fromIntegral n

-- La comprobación es

-- λ> quickCheck conjeturaDeMertens

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

«El control de la complejidad es la esencia de la programación informática.»

Brian Kernighan.

Números de Munchausen

PorJosé A. Alonso 13-enero-202020-enero-2020

Un número de Munchausen es un número entero positivo tal que es igual a la suma de sus dígitos elevados a sí mismo. Por ejemplo, 3435 es un número de Munchausen ya que

   3³ + 4⁴ + 3³ + 5⁵ = 27 + 256 + 27 + 3125 = 3435

1	3³ + 4⁴ + 3³ + 5⁵ = 27 + 256 + 27 + 3125 = 3435

Definir la función

   esMunchausen :: Integer -> Bool

1	esMunchausen :: Integer -> Bool

tal que (esMunchausen n) se verifica si n es un número de Munchausen. Por ejemplo,

   esMunchausen 3435  ==  True
   esMunchausen 2020  ==  False

1 2	esMunchausen 3435 == True esMunchausen 2020 == False

Comprobar con QuickCheck que que los únicos números de Munchausen son 1 y 3435.

Nota 1: No usar la propiedad en la definición.

Nota 2: El ejercicio está basado en el artículo ¿Por qué 3435 es uno de mis números favoritos? de Miguel Ángel Morales en El Aleph.

Soluciones

import Test.QuickCheck

esMunchausen :: Integer -> Bool
esMunchausen n =
  n == sum [x^x | x <- digitos n]

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

-- La propiedad es
prop_Munchausen :: Integer -> Property
prop_Munchausen n =
  n > 0
  ==>
  esMunchausen n == elem n [1, 3435]

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

import Test.QuickCheck

esMunchausen :: Integer -> Bool

esMunchausen n =

n == sum [x^x | x <- digitos n]

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

-- digitos 3435 == [3,4,3,5]

digitos :: Integer -> [Integer]

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

-- La propiedad es

prop_Munchausen :: Integer -> Property

prop_Munchausen n =

n > 0

==>

esMunchausen n == elem n [1, 3435]

-- La comprobación es

-- λ> quickCheck prop_Munchausen

-- +++ OK, passed 100 tests.

Pensamiento

Escribiré en tu abanico:
te quiero para olvidarte,
para quererte te olvido.

Antonio Machado

Cálculo de pi mediante el método de Newton

Soluciones

Otras soluciones

Pensamiento

Medias de dígitos de pi

Soluciones

Otras soluciones

Pensamiento

Pandemia

Soluciones

Otras soluciones

Pensamiento

La conjetura de Mertens

Soluciones

Otras soluciones

Pensamiento

Números de Munchausen

Soluciones

Pensamiento

Entradas recientes

Buscar

Correo electrónico