accumArray

Cálculo de dígitos de pi y su distribución

PorJosé A. Alonso 8-abril-202021-abril-2020

Se pueden generar los dígitos de Pi, como se explica en el artículo Unbounded spigot algorithms for the digits of pi c0on la función digitosPi definida por

   digitosPi :: [Integer]
   digitosPi = g(1,0,1,1,3,3) where
     g (q,r,t,k,n,l) = 
       if 4*q+r-t < n*t
       then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
       else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

digitosPi :: [Integer]

digitosPi = g(1,0,1,1,3,3) where

g (q,r,t,k,n,l) =

if 4*q+r-t < n*t

then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)

else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

Por ejemplo,

   λ> take 25 digitosPi
   [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

1 2	λ> take 25 digitosPi [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

La distribución de los primeros 25 dígitos de pi es [0,2,3,5,3,3,3,1,2,3] ya que el 0 no aparece, el 1 ocurre 2 veces, el 3 ocurre 3 veces, el 4 ocurre 5 veces, …

Usando digitosPi, definir las siguientes funciones

   distribucionDigitosPi :: Int -> [Int]
   frecuenciaDigitosPi   :: Int -> [Double]

1 2	distribucionDigitosPi :: Int -> [Int] frecuenciaDigitosPi :: Int -> [Double]

tales que

(distribucionDigitosPi n) es la distribución de los n primeros dígitos de pi. Por ejemplo,

     λ> distribucionDigitosPi 10
     [0,2,1,2,1,2,1,0,0,1]
     λ> distribucionDigitosPi 100
     [8,8,12,12,10,8,9,8,12,13]
     λ> distribucionDigitosPi 1000
     [93,116,103,103,93,97,94,95,101,105]
     λ> distribucionDigitosPi 5000
     [466,531,496,460,508,525,513,488,492,521]

λ> distribucionDigitosPi 10

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

λ> distribucionDigitosPi 100

[8,8,12,12,10,8,9,8,12,13]

λ> distribucionDigitosPi 1000

[93,116,103,103,93,97,94,95,101,105]

λ> distribucionDigitosPi 5000

[466,531,496,460,508,525,513,488,492,521]

(frecuenciaDigitosPi n) es la frecuencia de los n primeros dígitos de pi. Por ejemplo,

   λ> frecuenciaDigitosPi 10
   [0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]
   λ> frecuenciaDigitosPi 100
   [8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]
   λ> frecuenciaDigitosPi 1000
   [9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]
   λ> frecuenciaDigitosPi 5000
   [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

λ> frecuenciaDigitosPi 10

[0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]

λ> frecuenciaDigitosPi 100

[8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]

λ> frecuenciaDigitosPi 1000

[9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]

λ> frecuenciaDigitosPi 5000

[9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

Soluciones

import Data.Array
import Data.List (group, sort)

digitosPi :: [Integer]
digitosPi = g(1,0,1,1,3,3) where
  g (q,r,t,k,n,l) = 
    if 4*q+r-t < n*t
    then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
    else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

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

distribucionDigitosPi :: Int -> [Int]
distribucionDigitosPi n =
  elems (accumArray (+) 0 (0,9) [ (i,1)
                                | i <- take n digitosPi]) 

frecuenciaDigitosPi :: Int -> [Double]
frecuenciaDigitosPi n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]
  where m = fromIntegral n

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

distribucionDigitosPi2 :: Int -> [Int]
distribucionDigitosPi2 n =
  [length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]

frecuenciaDigitosPi2 :: Int -> [Double]
frecuenciaDigitosPi2 n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]
  where m = fromIntegral n

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

--    λ> last (take 5000 digitosPi)
--    2
--    (4.47 secs, 3,927,848,448 bytes)
--    λ> frecuenciaDigitosPi 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.01 secs, 0 bytes)
--    λ> frecuenciaDigitosPi2 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.02 secs, 0 bytes)

import Data.Array

import Data.List (group, sort)

digitosPi :: [Integer]

digitosPi = g(1,0,1,1,3,3) where

g (q,r,t,k,n,l) =

if 4*q+r-t < n*t

then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)

else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

-- 1ª definición

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

distribucionDigitosPi :: Int -> [Int]

distribucionDigitosPi n =

elems (accumArray (+) 0 (0,9) [ (i,1)

| i <- take n digitosPi])

frecuenciaDigitosPi :: Int -> [Double]

frecuenciaDigitosPi n =

[100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]

where m = fromIntegral n

-- 2ª definición

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

distribucionDigitosPi2 :: Int -> [Int]

distribucionDigitosPi2 n =

[length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]

frecuenciaDigitosPi2 :: Int -> [Double]

frecuenciaDigitosPi2 n =

[100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]

where m = fromIntegral n

-- Comparación de eficiencia

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

-- λ> last (take 5000 digitosPi)

-- 2

-- (4.47 secs, 3,927,848,448 bytes)

-- λ> frecuenciaDigitosPi 5000

-- [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

-- (0.01 secs, 0 bytes)

-- λ> frecuenciaDigitosPi2 5000

-- [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

-- (0.02 secs, 0 bytes)

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>

Distribución de diferencias de dígitos consecutivos de pi

PorJosé A. Alonso 4-abril-201911-abril-2019

Usando la librería Data.Number.CReal, que se instala con

   cabal install number

1	cabal install number

se pueden calcular el número pi con la precisión que se desee. Por ejemplo,

   λ> import Data.Number.CReal
   λ> showCReal 60 pi
   "3.141592653589793238462643383279502884197169399375105820974945"

λ> import Data.Number.CReal

λ> showCReal 60 pi

"3.141592653589793238462643383279502884197169399375105820974945"

importa la librería y calcula el número pi con 60 decimales.

La distribución de las diferencias de los dígitos consecutivos para los 18 primeros n 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

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)

Pensamiento

Doy consejo, a fuer de viejo:
nunca sigas mi consejo.

Antonio Machado

Cálculo de dígitos de pi y su distribución

PorJosé A. Alonso 3-abril-201925-abril-2021

Se pueden generar los dígitos de Pi, como se explica en el artículo Unbounded spigot algorithms for the digits of pi c0on la función digitosPi definida por

   digitosPi :: [Integer]
   digitosPi = g(1,0,1,1,3,3) where
     g (q,r,t,k,n,l) = 
       if 4*q+r-t < n*t
       then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
       else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

digitosPi :: [Integer]

digitosPi = g(1,0,1,1,3,3) where

g (q,r,t,k,n,l) =

if 4*q+r-t < n*t

then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)

else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

Por ejemplo,

   λ> take 25 digitosPi
   [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

1 2	λ> take 25 digitosPi [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

La distribución de los primeros 25 dígitos de pi es [0,2,3,5,3,3,3,1,2,3] ya que el 0 no aparece, el 1 ocurre 2 veces, el 3 ocurre 3 veces, el 4 ocurre 5 veces, …

Usando digitosPi, definir las siguientes funciones

   distribucionDigitosPi :: Int -> [Int]
   frecuenciaDigitosPi   :: Int -> [Double]

1 2	distribucionDigitosPi :: Int -> [Int] frecuenciaDigitosPi :: Int -> [Double]

tales que

(distribucionDigitosPi n) es la distribución de los n primeros dígitos de pi. Por ejemplo,

     λ> distribucionDigitosPi 10
     [0,2,1,2,1,2,1,0,0,1]
     λ> distribucionDigitosPi 100
     [8,8,12,12,10,8,9,8,12,13]
     λ> distribucionDigitosPi 1000
     [93,116,103,103,93,97,94,95,101,105]
     λ> distribucionDigitosPi 5000
     [466,531,496,460,508,525,513,488,492,521]

λ> distribucionDigitosPi 10

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

λ> distribucionDigitosPi 100

[8,8,12,12,10,8,9,8,12,13]

λ> distribucionDigitosPi 1000

[93,116,103,103,93,97,94,95,101,105]

λ> distribucionDigitosPi 5000

[466,531,496,460,508,525,513,488,492,521]

(frecuenciaDigitosPi n) es la frecuencia de los n primeros dígitos de pi. Por ejemplo,

   λ> frecuenciaDigitosPi 10
   [0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]
   λ> frecuenciaDigitosPi 100
   [8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]
   λ> frecuenciaDigitosPi 1000
   [9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]
   λ> frecuenciaDigitosPi 5000
   [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

λ> frecuenciaDigitosPi 10

[0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]

λ> frecuenciaDigitosPi 100

[8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]

λ> frecuenciaDigitosPi 1000

[9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]

λ> frecuenciaDigitosPi 5000

[9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

Soluciones

import Data.Array
import Data.List (group, sort)

digitosPi :: [Integer]
digitosPi = g(1,0,1,1,3,3) where
  g (q,r,t,k,n,l) = 
    if 4*q+r-t < n*t
    then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
    else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

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

distribucionDigitosPi :: Int -> [Int]
distribucionDigitosPi n =
    elems (accumArray (+) 0 (0,9) [(i,1)
                                  | i <- take n digitosPi]) 

frecuenciaDigitosPi :: Int -> [Double]
frecuenciaDigitosPi n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]
  where m = fromIntegral n

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

distribucionDigitosPi2 :: Int -> [Int]
distribucionDigitosPi2 n =
  [length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]

frecuenciaDigitosPi2 :: Int -> [Double]
frecuenciaDigitosPi2 n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]
  where m = fromIntegral n

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

--    λ> last (take 5000 digitosPi)
--    2
--    (4.47 secs, 3,927,848,448 bytes)
--    λ> frecuenciaDigitosPi 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.01 secs, 0 bytes)
--    λ> frecuenciaDigitosPi2 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.02 secs, 0 bytes)

import Data.Array

import Data.List (group, sort)

digitosPi :: [Integer]

digitosPi = g(1,0,1,1,3,3) where

g (q,r,t,k,n,l) =

if 4*q+r-t < n*t

then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)

else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

-- 1ª definición

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

distribucionDigitosPi :: Int -> [Int]

distribucionDigitosPi n =

elems (accumArray (+) 0 (0,9) [(i,1)

| i <- take n digitosPi])

frecuenciaDigitosPi :: Int -> [Double]

frecuenciaDigitosPi n =

[100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]

where m = fromIntegral n

-- 2ª definición

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

distribucionDigitosPi2 :: Int -> [Int]

distribucionDigitosPi2 n =

[length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]

frecuenciaDigitosPi2 :: Int -> [Double]

frecuenciaDigitosPi2 n =

[100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]

where m = fromIntegral n

-- Comparación de eficiencia

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

-- λ> last (take 5000 digitosPi)

-- 2

-- (4.47 secs, 3,927,848,448 bytes)

-- λ> frecuenciaDigitosPi 5000

-- [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

-- (0.01 secs, 0 bytes)

-- λ> frecuenciaDigitosPi2 5000

-- [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

-- (0.02 secs, 0 bytes)

Pensamiento

¿Cuál es la verdad? ¿El río
que fluye y pasa
donde el barco y el barquero
son también ondas de agua?
¿O este soñar del marino
siempre con ribera y ancla?

Antonio Machado

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)

Distribución de dígitos de pi

PorJosé A. Alonso 16-febrero-201725-abril-2021

Se pueden generar los dígitos de Pi, como se explica en el artículo Unbounded spigot algorithms for the digits of pi, con la función digitosPi definida por

   digitosPi :: [Integer]
   digitosPi = g(1,0,1,1,3,3) where
     g (q,r,t,k,n,l) = 
       if 4*q+r-t < n*t
       then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
       else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

digitosPi :: [Integer]

digitosPi = g(1,0,1,1,3,3) where

g (q,r,t,k,n,l) =

if 4*q+r-t < n*t

then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)

else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

Por ejemplo,

   λ> take 25 digitosPi
   [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

1 2	λ> take 25 digitosPi [3,1,4,1,5,9,2,6,5,3,5,8,9,7,9,3,2,3,8,4,6,2,6,4,3]

La distribución de los primeros 25 dígitos de pi es [0,2,3,5,3,3,3,1,2,3] ya que el 0 no aparece, el 1 ocurre 2 veces, el 3 ocurre 3 veces, el 4 ocurre 5 veces, …

Usando digitosPi, definir las siguientes funciones

   distribucionDigitosPi :: Int -> [Int]
   frecuenciaDigitosPi   :: Int -> [Double]

1 2	distribucionDigitosPi :: Int -> [Int] frecuenciaDigitosPi :: Int -> [Double]

tales que

(distribucionDigitosPi n) es la distribución de los n primeros dígitos de pi. Por ejemplo,

     λ> distribucionDigitosPi 10
     [0,2,1,2,1,2,1,0,0,1]
     λ> distribucionDigitosPi 100
     [8,8,12,12,10,8,9,8,12,13]
     λ> distribucionDigitosPi 1000
     [93,116,103,103,93,97,94,95,101,105]
     λ> distribucionDigitosPi 5000
     [466,531,496,460,508,525,513,488,492,521]

λ> distribucionDigitosPi 10

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

λ> distribucionDigitosPi 100

[8,8,12,12,10,8,9,8,12,13]

λ> distribucionDigitosPi 1000

[93,116,103,103,93,97,94,95,101,105]

λ> distribucionDigitosPi 5000

[466,531,496,460,508,525,513,488,492,521]

(frecuenciaDigitosPi n) es la frecuencia de los n primeros dígitos de pi. Por ejemplo,

   λ> frecuenciaDigitosPi 10
   [0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]
   λ> frecuenciaDigitosPi 100
   [8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]
   λ> frecuenciaDigitosPi 1000
   [9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]
   λ> frecuenciaDigitosPi 5000
   [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

λ> frecuenciaDigitosPi 10

[0.0,20.0,10.0,20.0,10.0,20.0,10.0,0.0,0.0,10.0]

λ> frecuenciaDigitosPi 100

[8.0,8.0,12.0,12.0,10.0,8.0,9.0,8.0,12.0,13.0]

λ> frecuenciaDigitosPi 1000

[9.3,11.6,10.3,10.3,9.3,9.7,9.4,9.5,10.1,10.5]

λ> frecuenciaDigitosPi 5000

[9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

Soluciones

import Data.Array
import Data.List (group, sort)

digitosPi :: [Integer]
digitosPi = g(1,0,1,1,3,3) where
  g (q,r,t,k,n,l) = 
    if 4*q+r-t < n*t
    then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)
    else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

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

distribucionDigitosPi :: Int -> [Int]
distribucionDigitosPi n =
    elems (accumArray (+) 0 (0,9) [(i,1)
                                  | i <- take n digitosPi]) 

frecuenciaDigitosPi :: Int -> [Double]
frecuenciaDigitosPi n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]
  where m = fromIntegral n

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

distribucionDigitosPi2 :: Int -> [Int]
distribucionDigitosPi2 n =
  [length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]

frecuenciaDigitosPi2 :: Int -> [Double]
frecuenciaDigitosPi2 n =
  [100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]
  where m = fromIntegral n

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

--    λ> last (take 5000 digitosPi)
--    2
--    (4.47 secs, 3,927,848,448 bytes)
--    λ> frecuenciaDigitosPi 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.01 secs, 0 bytes)
--    λ> frecuenciaDigitosPi2 5000
--    [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]
--    (0.02 secs, 0 bytes)

import Data.Array

import Data.List (group, sort)

digitosPi :: [Integer]

digitosPi = g(1,0,1,1,3,3) where

g (q,r,t,k,n,l) =

if 4*q+r-t < n*t

then n : g (10*q, 10*(r-n*t), t, k, div (10*(3*q+r)) t - 10*n, l)

else g (q*k, (2*q+r)*l, t*l, k+1, div (q*(7*k+2)+r*l) (t*l), l+2)

-- 1ª definición

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

distribucionDigitosPi :: Int -> [Int]

distribucionDigitosPi n =

elems (accumArray (+) 0 (0,9) [(i,1)

| i <- take n digitosPi])

frecuenciaDigitosPi :: Int -> [Double]

frecuenciaDigitosPi n =

[100 * (fromIntegral x / m) | x <- distribucionDigitosPi n]

where m = fromIntegral n

-- 2ª definición

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

distribucionDigitosPi2 :: Int -> [Int]

distribucionDigitosPi2 n =

[length xs - 1 | xs <- group (sort (take n digitosPi ++ [0..9]))]

frecuenciaDigitosPi2 :: Int -> [Double]

frecuenciaDigitosPi2 n =

[100 * (fromIntegral x / m) | x <- distribucionDigitosPi2 n]

where m = fromIntegral n

-- Comparación de eficiencia

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

-- λ> last (take 5000 digitosPi)

-- 2

-- (4.47 secs, 3,927,848,448 bytes)

-- λ> frecuenciaDigitosPi 5000

-- [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

-- (0.01 secs, 0 bytes)

-- λ> frecuenciaDigitosPi2 5000

-- [9.32,10.62,9.92,9.2,10.16,10.5,10.26,9.76,9.84,10.42]

-- (0.02 secs, 0 bytes)

Comportamiento del último dígito en primos consecutivos

PorJosé A. Alonso 18-marzo-20162-mayo-2016

El pasado 11 de marzo se ha publicado el artículo Unexpected biases in the distribution of consecutive primes en el que muestra que los números primos repelen a otros primos que terminan en el mismo dígito.

La lista de los últimos dígitos de los 30 primeros números es

   [2,3,5,7,1,3,7,9,3,9,1,7,1,3,7,3,9,1,7,1,3,9,3,9,7,1,3,7,9,3]

1	[2,3,5,7,1,3,7,9,3,9,1,7,1,3,7,3,9,1,7,1,3,9,3,9,7,1,3,7,9,3]

Se observa que hay 6 números que su último dígito es un 1 y de sus consecutivos 4 terminan en 3 y 2 terminan en 7.

Definir la función

   distribucionUltimos :: Int -> M.Matrix Integer

1	distribucionUltimos :: Int -> M.Matrix Integer

tal que (distribucionUltimos n) es la matriz cuyo elemento (i,j) indica cuántos de los n primeros números primos terminan en i y su siguiente número primo termina en j. Por ejemplo,

   λ> distribucionUltimos 30
   ( 0 0 4 0 0 0 2 0 0 )
   ( 0 0 1 0 0 0 0 0 0 )
   ( 0 0 0 0 1 0 4 0 4 )
   ( 0 0 0 0 0 0 0 0 0 )
   ( 0 0 0 0 0 0 1 0 0 )
   ( 0 0 0 0 0 0 0 0 0 )
   ( 4 0 1 0 0 0 0 0 2 )
   ( 0 0 0 0 0 0 0 0 0 )
   ( 2 0 3 0 0 0 1 0 0 )
   
   λ> distribucionUltimos (10^5)
   ( 4104    0 7961    0    0    0 8297    0 4605 )
   (    0    0    1    0    0    0    0    0    0 )
   ( 5596    0 3604    0    1    0 7419    0 8387 )
   (    0    0    0    0    0    0    0    0    0 )
   (    0    0    0    0    0    0    1    0    0 )
   (    0    0    0    0    0    0    0    0    0 )
   ( 6438    0 6928    0    0    0 3627    0 8022 )
   (    0    0    0    0    0    0    0    0    0 )
   ( 8830    0 6513    0    0    0 5671    0 3995 )

λ> distribucionUltimos 30

( 0 0 4 0 0 0 2 0 0 )

( 0 0 1 0 0 0 0 0 0 )

( 0 0 0 0 1 0 4 0 4 )

( 0 0 0 0 0 0 0 0 0 )

( 0 0 0 0 0 0 1 0 0 )

( 0 0 0 0 0 0 0 0 0 )

( 4 0 1 0 0 0 0 0 2 )

( 0 0 0 0 0 0 0 0 0 )

( 2 0 3 0 0 0 1 0 0 )

λ> distribucionUltimos (10^5)

( 4104 0 7961 0 0 0 8297 0 4605 )

( 0 0 1 0 0 0 0 0 0 )

( 5596 0 3604 0 1 0 7419 0 8387 )

( 0 0 0 0 0 0 0 0 0 )

( 0 0 0 0 0 0 1 0 0 )

( 0 0 0 0 0 0 0 0 0 )

( 6438 0 6928 0 0 0 3627 0 8022 )

( 0 0 0 0 0 0 0 0 0 )

( 8830 0 6513 0 0 0 5671 0 3995 )

Nota: Se observa cómo se «repelen» ya que en las filas del 1, 3, 7 y 9 el menor elemento es el de la diagonal.

Soluciones

import Data.Numbers.Primes
import Data.Array
import qualified Data.Matrix as M

-- (ultimo n) es el último dígito de n.    
ultimo :: Integer -> Integer
ultimo n = n `mod` 10

-- ultimos es la lista de el último dígito de los primos.
--    λ> take 20 ultimos
--    [2,3,5,7,1,3,7,9,3,9,1,7,1,3,7,3,9,1,7,1]
ultimos :: [Integer]
ultimos = map ultimo primes

-- ultimosConsecutivos es la lista de los últimos dígitos de los primos
-- consecutivos. 
--    λ> take 10 ultimosConsecutivos
--    [(2,3),(3,5),(5,7),(7,1),(1,3),(3,7),(7,9),(9,3),(3,9),(9,1)]
ultimosConsecutivos :: [(Integer,Integer)]
ultimosConsecutivos = zip ultimos (tail ultimos)

-- (histograma r is) es el vector formado contando cuantas veces
-- aparecen los elementos del rango r en la lista de índices is. Por
-- ejemplo, 
--    ghci> histograma (0,5) [3,1,4,1,5,4,2,7]
--    array (0,5) [(0,0),(1,2),(2,1),(3,1),(4,2),(5,1)]
histograma :: (Ix a, Num b) => (a,a) -> [a] -> Array a b
histograma r is = 
    accumArray (+) 0 r [(i,1) | i <- is, inRange r i]

distribucionUltimos :: Int -> M.Matrix Integer
distribucionUltimos n =
    M.fromList 9 9 (elems (histograma ((1,1),(9,9)) (take n ultimosConsecutivos)))

import Data.Numbers.Primes

import Data.Array

import qualified Data.Matrix as M

-- (ultimo n) es el último dígito de n.

ultimo :: Integer -> Integer

ultimo n = n `mod` 10

-- ultimos es la lista de el último dígito de los primos.

-- λ> take 20 ultimos

-- [2,3,5,7,1,3,7,9,3,9,1,7,1,3,7,3,9,1,7,1]

ultimos :: [Integer]

ultimos = map ultimo primes

-- ultimosConsecutivos es la lista de los últimos dígitos de los primos

-- consecutivos.

-- λ> take 10 ultimosConsecutivos

-- [(2,3),(3,5),(5,7),(7,1),(1,3),(3,7),(7,9),(9,3),(3,9),(9,1)]

ultimosConsecutivos :: [(Integer,Integer)]

ultimosConsecutivos = zip ultimos (tail ultimos)

-- (histograma r is) es el vector formado contando cuantas veces

-- aparecen los elementos del rango r en la lista de índices is. Por

-- ejemplo,

-- ghci> histograma (0,5) [3,1,4,1,5,4,2,7]

-- array (0,5) [(0,0),(1,2),(2,1),(3,1),(4,2),(5,1)]

histograma :: (Ix a, Num b) => (a,a) -> [a] -> Array a b

histograma r is =

accumArray (+) 0 r [(i,1) | i <- is, inRange r i]

distribucionUltimos :: Int -> M.Matrix Integer

distribucionUltimos n =

M.fromList 9 9 (elems (histograma ((1,1),(9,9)) (take n ultimosConsecutivos)))

Solución en Maxima

distribucionUltimos (n) := block (
  [r : zeromatrix (9,9),
   xs : ultimos (n),
   i, j],
  unless length (xs) < 2 do
    ( i : first (xs),
      j : second (xs),
      r[i,j] : 1 + r[i,j],
      xs : rest (xs)),
  r)$

/* ultimos(n) es la lista del último dígito de los n primeros primos.
      (%i5) ultimos (30);
      (%o5) [2,3,5,7,1,3,7,9,3,9,1,7,1,3,7,3,9,1,7,1,3,9,3,9,7,1,3,7,9,3,7]
*/
ultimos (n) := block ([r:[], p:2],
  for k from 0 thru n do
    ( r : cons (mod (p,10), r),
      p : next_prime (p) ),
  reverse (r))$

distribucionUltimos (n) := block (

[r : zeromatrix (9,9),

xs : ultimos (n),

i, j],

unless length (xs) < 2 do

( i : first (xs),

j : second (xs),

r[i,j] : 1 + r[i,j],

xs : rest (xs)),

r)$

/* ultimos(n) es la lista del último dígito de los n primeros primos.

(%i5) ultimos (30);

(%o5) [2,3,5,7,1,3,7,9,3,9,1,7,1,3,7,3,9,1,7,1,3,9,3,9,7,1,3,7,9,3,7]

ultimos (n) := block ([r:[], p:2],

for k from 0 thru n do

( r : cons (mod (p,10), r),

p : next_prime (p) ),

reverse (r))$

Perímetro más frecuente de triángulos rectángulos

PorJosé A. Alonso 7-abril-20151-mayo-2015

El grado perimetral de un número p es la cantidad de tres triángulos rectángulos de lados enteros cuyo perímetro es p. Por ejemplo, el grado perimetral de 120 es 3 ya que sólo hay 3 triángulos rectángulos de lados enteros cuyo perímetro es 120: {20,48,52}, {24,45,51} y {30,40,50}.

Definir la función

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

1	maxGradoPerimetral :: Int -> (Int,[Int])

tal que (maxGradoPerimetral n) es el par (m,ps) tal que m es el máximo grado perimetral de los números menores o iguales que n y ps son los perímetros, menores o iguales que n, cuyo grado perimetral es m. Por ejemplo,

   maxGradoPerimetral   50  ==  (1,[12,24,30,36,40,48])
   maxGradoPerimetral  100  ==  (2,[60,84,90])
   maxGradoPerimetral  200  ==  (3,[120,168,180])
   maxGradoPerimetral  400  ==  (4,[240,360])
   maxGradoPerimetral  500  ==  (5,[420])
   maxGradoPerimetral  750  ==  (6,[720])
   maxGradoPerimetral  839  ==  (6,[720])
   maxGradoPerimetral  840  ==  (8,[840])
   maxGradoPerimetral 1500  ==  (8,[840,1260])
   maxGradoPerimetral 2000  ==  (10,[1680])
   maxGradoPerimetral 3000  ==  (12,[2520])

maxGradoPerimetral 50 == (1,[12,24,30,36,40,48])

maxGradoPerimetral 100 == (2,[60,84,90])

maxGradoPerimetral 200 == (3,[120,168,180])

maxGradoPerimetral 400 == (4,[240,360])

maxGradoPerimetral 500 == (5,[420])

maxGradoPerimetral 750 == (6,[720])

maxGradoPerimetral 839 == (6,[720])

maxGradoPerimetral 840 == (8,[840])

maxGradoPerimetral 1500 == (8,[840,1260])

maxGradoPerimetral 2000 == (10,[1680])

maxGradoPerimetral 3000 == (12,[2520])

Soluciones

import Data.List
import Data.Array (accumArray, assocs)

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

maxGradoPerimetral1 :: Int -> (Int,[Int])
maxGradoPerimetral1 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(length (triangulos x),x) | x <- [1..p]] 
          (m,_) = maximum ts 

-- (triangulos p) es el conjunto de triángulos rectángulos de perímetro
-- p. Por ejemplo,
--    triangulos 120  ==  [(20,48,52),(24,45,51),(30,40,50)]
triangulos :: Int -> [(Int,Int,Int)]
triangulos p = 
    [(a,b,c) | a <- [1..q],
               b <- [a..q],
               let c = p-a-b,
               a*a+b*b == c*c]
    where q = p `div` 2

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

maxGradoPerimetral2 :: Int -> (Int,[Int])
maxGradoPerimetral2 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(n,x) | (x,n) <- numeroPerimetrosTriangulos p, n > 0]
          (m,_) = maximum ts 

-- (numeroPerimetrosTriangulos p) es la lista formado por los números de
-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro
-- es dicho número. Por ejemplo,
--    ghci>  [(p,n) | (p,n) <- numeroPerimetrosTriangulos 70, n > 0]
--    [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]
numeroPerimetrosTriangulos :: Int -> [(Int,Int)] 
numeroPerimetrosTriangulos p = 
    assocs (accumArray (\x _ -> 1+x) 0 (1,p) (perimetrosTriangulos p))

-- (perimetrosTriangulos p) es la lista formada por los perímetros y los
-- lados de los triángulos rectángulos enteros cuyo perímetro es menor o
-- igual que p. Por ejemplo,
--    ghci> perimetrosTriangulos 70
--    [(12,(3,4,5)),   (30,(5,12,13)),(24,(6,8,10)),  (56,(7,24,25)),
--     (40,(8,15,17)), (36,(9,12,15)),(60,(10,24,26)),(48,(12,16,20)),
--     (60,(15,20,25)),(70,(20,21,29))]
perimetrosTriangulos :: Int -> [(Int,(Int,Int,Int))]
perimetrosTriangulos p =
    [(q,(a,b,c)) | a <- [1..p1],
                   b <- [a..p1],
                   esCuadrado (a*a+b*b), 
                   let c = raizCuadradaE (a*a+b*b), 
                   let q = a+b+c,
                   q <= p]
    where p1 = p `div` 2

-- (esCuadrado n) se verifica si n es un cuadrado. Por ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 27  ==  False
esCuadrado :: Int -> Bool
esCuadrado n = a*a == n
    where a = raizCuadradaE n

-- (raizCuadradaE n) es la raíz cuadrada entera de n. Por ejemplo,
--    raizCuadradaE 25  ==  5
--    raizCuadradaE 27  ==  5
--    raizCuadradaE 35  ==  5
--    raizCuadradaE 36  ==  6
raizCuadradaE :: Int -> Int
raizCuadradaE = floor . sqrt . fromIntegral

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

maxGradoPerimetral3 :: Int -> (Int,[Int])
maxGradoPerimetral3 p = (m,[x | (n,x) <- ts, n == m])
    where ts    = [(n,x) | (x,n) <- numeroPerimetrosTriangulos2 p, n > 0]
          (m,_) = maximum ts 

-- (numeroPerimetrosTriangulos2 p) es la lista formado por los números de
-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro
-- es dicho número. Por ejemplo,
--    ghci>  [(p,n) | (p,n) <- numeroPerimetrosTriangulos2 70, n > 0]
--    [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]
numeroPerimetrosTriangulos2 :: Int -> [(Int,Int)] 
numeroPerimetrosTriangulos2 p = 
    [(head xs, length xs) | xs <- group (sort (perimetrosTriangulos2 p))]

-- (perimetrosTriangulos2 p) es la lista formada por los perímetros de
-- los triángulos rectángulos enteros cuyo perímetro es menor o igual
-- que p. Por ejemplo, 
--    perimetrosTriangulos2 70  ==  [12,30,24,56,40,36,60,48,60,70]
perimetrosTriangulos2 :: Int -> [Int]
perimetrosTriangulos2 p =
    [q | a <- [1..p1],
         b <- [a..p1],
         esCuadrado (a*a+b*b), 
         let c = raizCuadradaE (a*a+b*b), 
         let q = a+b+c,
         q <= p]
    where p1 = p `div` 2

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

-- La comparación es
--    ghci> maxGradoPerimetral1 1000
--    (8,[840])
--    (120.08 secs, 21116625136 bytes)
--    ghci> maxGradoPerimetral2 1000
--    (8,[840])
--    (0.66 secs, 132959056 bytes)
--    ghci> maxGradoPerimetral3 1000
--    (1000,[1])
--    (0.66 secs, 133443816 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

import Data.List

import Data.Array (accumArray, assocs)

-- 1ª solución --

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

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

maxGradoPerimetral1 p = (m,[x | (n,x) <- ts, n == m])

where ts = [(length (triangulos x),x) | x <- [1..p]]

(m,_) = maximum ts

-- (triangulos p) es el conjunto de triángulos rectángulos de perímetro

-- p. Por ejemplo,

-- triangulos 120 == [(20,48,52),(24,45,51),(30,40,50)]

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

triangulos p =

[(a,b,c) | a <- [1..q],

b <- [a..q],

let c = p-a-b,

a*a+b*b == c*c]

where q = p `div` 2

-- 2ª solución --

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

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

maxGradoPerimetral2 p = (m,[x | (n,x) <- ts, n == m])

where ts = [(n,x) | (x,n) <- numeroPerimetrosTriangulos p, n > 0]

(m,_) = maximum ts

-- (numeroPerimetrosTriangulos p) es la lista formado por los números de

-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro

-- es dicho número. Por ejemplo,

-- ghci> [(p,n) | (p,n) <- numeroPerimetrosTriangulos 70, n > 0]

-- [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]

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

numeroPerimetrosTriangulos p =

assocs (accumArray (\x _ -> 1+x) 0 (1,p) (perimetrosTriangulos p))

-- (perimetrosTriangulos p) es la lista formada por los perímetros y los

-- lados de los triángulos rectángulos enteros cuyo perímetro es menor o

-- igual que p. Por ejemplo,

-- ghci> perimetrosTriangulos 70

-- [(12,(3,4,5)), (30,(5,12,13)),(24,(6,8,10)), (56,(7,24,25)),

-- (40,(8,15,17)), (36,(9,12,15)),(60,(10,24,26)),(48,(12,16,20)),

-- (60,(15,20,25)),(70,(20,21,29))]

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

perimetrosTriangulos p =

[(q,(a,b,c)) | a <- [1..p1],

b <- [a..p1],

esCuadrado (a*a+b*b),

let c = raizCuadradaE (a*a+b*b),

let q = a+b+c,

q <= p]

where p1 = p `div` 2

-- (esCuadrado n) se verifica si n es un cuadrado. Por ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 27 == False

esCuadrado :: Int -> Bool

esCuadrado n = a*a == n

where a = raizCuadradaE n

-- (raizCuadradaE n) es la raíz cuadrada entera de n. Por ejemplo,

-- raizCuadradaE 25 == 5

-- raizCuadradaE 27 == 5

-- raizCuadradaE 35 == 5

-- raizCuadradaE 36 == 6

raizCuadradaE :: Int -> Int

raizCuadradaE = floor . sqrt . fromIntegral

-- 3ª solución --

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

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

maxGradoPerimetral3 p = (m,[x | (n,x) <- ts, n == m])

where ts = [(n,x) | (x,n) <- numeroPerimetrosTriangulos2 p, n > 0]

(m,_) = maximum ts

-- (numeroPerimetrosTriangulos2 p) es la lista formado por los números de

-- 1 a p y la cantidad de triángulos rectángulos enteros cuyo perímetro

-- es dicho número. Por ejemplo,

-- ghci> [(p,n) | (p,n) <- numeroPerimetrosTriangulos2 70, n > 0]

-- [(12,1),(24,1),(30,1),(36,1),(40,1),(48,1),(56,1),(60,2),(70,1)]

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

numeroPerimetrosTriangulos2 p =

[(head xs, length xs) | xs <- group (sort (perimetrosTriangulos2 p))]

-- (perimetrosTriangulos2 p) es la lista formada por los perímetros de

-- los triángulos rectángulos enteros cuyo perímetro es menor o igual

-- que p. Por ejemplo,

-- perimetrosTriangulos2 70 == [12,30,24,56,40,36,60,48,60,70]

perimetrosTriangulos2 :: Int -> [Int]

perimetrosTriangulos2 p =

[q | a <- [1..p1],

b <- [a..p1],

esCuadrado (a*a+b*b),

let c = raizCuadradaE (a*a+b*b),

let q = a+b+c,

q <= p]

where p1 = p `div` 2

-- Comparación de eficiencia

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

-- La comparación es

-- ghci> maxGradoPerimetral1 1000

-- (8,[840])

-- (120.08 secs, 21116625136 bytes)

-- ghci> maxGradoPerimetral2 1000

-- (8,[840])

-- (0.66 secs, 132959056 bytes)

-- ghci> maxGradoPerimetral3 1000

-- (1000,[1])

-- (0.66 secs, 133443816 bytes)

Suma de conjuntos de polinomios

PorJosé A. Alonso 25-febrero-20151-mayo-2016

Los conjuntos de polinomios se pueden representar por listas de listas de la misma longitud. Por ejemplo, los polinomios 3x²+5x+9, 10x³+9 y 8x³+5x²+x-1 se pueden representar por las listas [0,3,5,9], [10,0,0,9] y [8,5,1,-1].

Definir la función

   sumaPolinomios :: Num a => [[a]] -> [a]

1	sumaPolinomios :: Num a => [[a]] -> [a]

tal que (sumaPolinomios ps) es la suma de los polinomios ps. Por ejemplo,

   ghci> sumaPolinomios1 [[0,3,5,9],[10,0,0,9],[8,5,1,-1]]
   [18,8,6,17]
   ghci> sumaPolinomios6 (replicate 1000000 (replicate 3 1))
   [1000000,1000000,1000000]

ghci> sumaPolinomios1 [[0,3,5,9],[10,0,0,9],[8,5,1,-1]]

[18,8,6,17]

ghci> sumaPolinomios6 (replicate 1000000 (replicate 3 1))

[1000000,1000000,1000000]

Soluciones

import Data.List (transpose)
import Data.Array ((!),accumArray,elems,listArray)

-- 1ª definición (por recursión):
sumaPolinomios1 :: Num a => [[a]] -> [a]
sumaPolinomios1 []          = []
sumaPolinomios1 [xs]        = xs
sumaPolinomios1 (xs:ys:zss) = suma xs (sumaPolinomios1 (ys:zss))

-- (suma xs ys) es la suma de los vectores xs e ys. Por ejemplo,
--    suma [4,7,3] [1,2,5]  == [5,9,8]
suma :: Num a => [a] -> [a] -> [a]
suma [] []         = []
suma (x:xs) (y:ys) = x+y : suma xs ys

-- 2ª definición (por recursión con zipWith): 
sumaPolinomios2 :: Num a => [[a]] -> [a]
sumaPolinomios2 []       = []
sumaPolinomios2 [xs]     = xs
sumaPolinomios2 (xs:xss) = zipWith (+) xs (sumaPolinomios2 xss)

-- 3ª definición (por plegado)
sumaPolinomios3 :: Num a => [[a]] -> [a]
sumaPolinomios3 = foldr1 (zipWith (+))

-- 4ª definición (por comprensión con transpose):
sumaPolinomios4 :: Num a => [[a]] -> [a]
sumaPolinomios4 xss = [sum xs | xs <- transpose xss]

-- 5ª definición (con map y transpose):
sumaPolinomios5 :: Num a => [[a]] -> [a]
sumaPolinomios5 = map sum . transpose 

-- 6ª definición (con array)
sumaPolinomios6 :: Num a => [[a]] -> [a]
sumaPolinomios6 xss = [sum [p!(i,j) | i <- [1..m]] | j <- [1..n]] 
    where m = length xss
          n = length (head xss)
          p = listArray ((1,1),(m,n)) (concat xss) 

-- 7ª definición (con accumArray)
sumaPolinomios7 :: Num a => [[a]] -> [a]
sumaPolinomios7 xss = 
    elems $ accumArray (+) 0 (1,n) (concat [zip [1..] xs | xs <- xss])
    where n = length (head xss)

-- Comparación de eficiencia
--    ghci> sumaPolinomios1 (replicate 300000 (replicate 5 1))
--    [300000,300000,300000,300000,300000]
--    (3.94 secs, 354713532 bytes)
--    
--    ghci> sumaPolinomios2 (replicate 300000 (replicate 5 1))
--    [300000,300000,300000,300000,300000]
--    (2.08 secs, 185506908 bytes)
--    
--    ghci> sumaPolinomios3 (replicate 300000 (replicate 5 1))
--    [300000,300000,300000,300000,300000]
--    (1.48 secs, 167026728 bytes)
--    
--    ghci> sumaPolinomios4 (replicate 300000 (replicate 5 1))
--    [300000,300000,300000,300000,300000]
--    (1.08 secs, 148564752 bytes)
--    
--    ghci> sumaPolinomios5 (replicate 300000 (replicate 5 1))
--    [300000,300000,300000,300000,300000]
--    (1.02 secs, 148062764 bytes)
--    
--    ghci> sumaPolinomios6 (replicate 300000 (replicate 5 1))
--    [300000,300000,300000,300000,300000]
--    (3.17 secs, 463756028 bytes)
--    
--    ghci> sumaPolinomios7 (replicate 300000 (replicate 5 1))
--    [300000,300000,300000,300000,300000]
--    (1.50 secs, 291699548 bytes)

import Data.List (transpose)

import Data.Array ((!),accumArray,elems,listArray)

-- 1ª definición (por recursión):

sumaPolinomios1 :: Num a => [[a]] -> [a]

sumaPolinomios1 [] = []

sumaPolinomios1 [xs] = xs

sumaPolinomios1 (xs:ys:zss) = suma xs (sumaPolinomios1 (ys:zss))

-- (suma xs ys) es la suma de los vectores xs e ys. Por ejemplo,

-- suma [4,7,3] [1,2,5] == [5,9,8]

suma :: Num a => [a] -> [a] -> [a]

suma [] [] = []

suma (x:xs) (y:ys) = x+y : suma xs ys

-- 2ª definición (por recursión con zipWith):

sumaPolinomios2 :: Num a => [[a]] -> [a]

sumaPolinomios2 [] = []

sumaPolinomios2 [xs] = xs

sumaPolinomios2 (xs:xss) = zipWith (+) xs (sumaPolinomios2 xss)

-- 3ª definición (por plegado)

sumaPolinomios3 :: Num a => [[a]] -> [a]

sumaPolinomios3 = foldr1 (zipWith (+))

-- 4ª definición (por comprensión con transpose):

sumaPolinomios4 :: Num a => [[a]] -> [a]

sumaPolinomios4 xss = [sum xs | xs <- transpose xss]

-- 5ª definición (con map y transpose):

sumaPolinomios5 :: Num a => [[a]] -> [a]

sumaPolinomios5 = map sum . transpose

-- 6ª definición (con array)

sumaPolinomios6 :: Num a => [[a]] -> [a]

sumaPolinomios6 xss = [sum [p!(i,j) | i <- [1..m]] | j <- [1..n]]

where m = length xss

n = length (head xss)

p = listArray ((1,1),(m,n)) (concat xss)

-- 7ª definición (con accumArray)

sumaPolinomios7 :: Num a => [[a]] -> [a]

sumaPolinomios7 xss =

elems $ accumArray (+) 0 (1,n) (concat [zip [1..] xs | xs <- xss])

where n = length (head xss)

-- Comparación de eficiencia

-- ghci> sumaPolinomios1 (replicate 300000 (replicate 5 1))

-- [300000,300000,300000,300000,300000]

-- (3.94 secs, 354713532 bytes)

-- ghci> sumaPolinomios2 (replicate 300000 (replicate 5 1))

-- [300000,300000,300000,300000,300000]

-- (2.08 secs, 185506908 bytes)

-- ghci> sumaPolinomios3 (replicate 300000 (replicate 5 1))

-- [300000,300000,300000,300000,300000]

-- (1.48 secs, 167026728 bytes)

-- ghci> sumaPolinomios4 (replicate 300000 (replicate 5 1))

-- [300000,300000,300000,300000,300000]

-- (1.08 secs, 148564752 bytes)

-- ghci> sumaPolinomios5 (replicate 300000 (replicate 5 1))

-- [300000,300000,300000,300000,300000]

-- (1.02 secs, 148062764 bytes)

-- ghci> sumaPolinomios6 (replicate 300000 (replicate 5 1))

-- [300000,300000,300000,300000,300000]

-- (3.17 secs, 463756028 bytes)

-- ghci> sumaPolinomios7 (replicate 300000 (replicate 5 1))

-- [300000,300000,300000,300000,300000]

-- (1.50 secs, 291699548 bytes)

Distribución de diferencias de dígitos consecutivos de pi

Soluciones

Pensamiento

Distribución de diferencias de dígitos consecutivos de pi

Soluciones

Comportamiento del último dígito en primos consecutivos

Soluciones

Solución en Maxima

Suma de conjuntos de polinomios

Soluciones

Entradas recientes

Buscar

Correo electrónico