fromIntegral – Exercitium

Operaciones con series de potencias

PorJosé A. Alonso 29-mayo-202025-marzo-2022

Una serie de potencias es una serie de la forma

   a₀ + a₁x + a₂x² + a₃x³ + ...

1	a₀ + a₁x + a₂x² + a₃x³ + ...

Las series de potencias se pueden representar mediante listas infinitas. Por ejemplo, la serie de la función exponencial es

   e^x = 1 + x + x²/2! + x³/3! + ...

1	e^x = 1 + x + x²/2! + x³/3! + ...

y se puede representar por [1, 1, 1/2, 1/6, 1/24, 1/120, …]

Las operaciones con series se pueden ver como una generalización de las de los polinomios.

En lo que sigue, usaremos el tipo (Serie a) para representar las series de potencias con coeficientes en a y su definición es

   type Serie a = [a]

1	type Serie a = [a]

Definir las siguientes funciones

   opuesta      :: Num a => Serie a -> Serie a
   suma         :: Num a => Serie a -> Serie a -> Serie a
   resta        :: Num a => Serie a -> Serie a -> Serie a
   producto     :: Num a => Serie a -> Serie a -> Serie a
   cociente     :: Fractional a => Serie a -> Serie a -> Serie a
   derivada     :: (Num a, Enum a) => Serie a -> Serie a
   integral     :: (Fractional a, Enum a) => Serie a -> Serie a
   expx         :: Serie Rational

opuesta :: Num a => Serie a -> Serie a

suma :: Num a => Serie a -> Serie a -> Serie a

resta :: Num a => Serie a -> Serie a -> Serie a

producto :: Num a => Serie a -> Serie a -> Serie a

cociente :: Fractional a => Serie a -> Serie a -> Serie a

derivada :: (Num a, Enum a) => Serie a -> Serie a

integral :: (Fractional a, Enum a) => Serie a -> Serie a

expx :: Serie Rational

tales que

(opuesta xs) es la opuesta de la serie xs. Por ejemplo,

     λ> take 7 (opuesta [-6,-4..])
     [6,4,2,0,-2,-4,-6]

1 2	λ> take 7 (opuesta [-6,-4..]) [6,4,2,0,-2,-4,-6]

(suma xs ys) es la suma de las series xs e ys. Por ejemplo,

     λ> take 7 (suma [1,3..] [2,4..])
     [3,7,11,15,19,23,27]

1 2	λ> take 7 (suma [1,3..] [2,4..]) [3,7,11,15,19,23,27]

(resta xs ys) es la resta de las series xs es ys. Por ejemplo,

     λ> take 7 (resta [3,5..] [2,4..])
     [1,1,1,1,1,1,1]
     λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])
     [2,4,6,8,10,12,14]

λ> take 7 (resta [3,5..] [2,4..])

[1,1,1,1,1,1,1]

λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])

[2,4,6,8,10,12,14]

(producto xs ys) es el producto de las series xs e ys. Por ejemplo,

     λ> take 7 (producto [3,5..] [2,4..])
     [6,22,52,100,170,266,392]

1 2	λ> take 7 (producto [3,5..] [2,4..]) [6,22,52,100,170,266,392]

(cociente xs ys) es el cociente de las series xs e ys. Por ejemplo,

     λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..])
     [2.0,4.0,6.0,8.0,10.0,12.0,14.0]

1 2	λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..]) [2.0,4.0,6.0,8.0,10.0,12.0,14.0]

(derivada xs) es la derivada de la serie xs. Por ejemplo,

     λ> take 7 (derivada [2,4..])
     [4,12,24,40,60,84,112]

1 2	λ> take 7 (derivada [2,4..]) [4,12,24,40,60,84,112]

(integral xs) es la integral de la serie xs. Por ejemplo,

     λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0))
     [0.0,4.0,6.0,8.0,10.0,12.0,14.0]

1 2	λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0)) [0.0,4.0,6.0,8.0,10.0,12.0,14.0]

expx es la serie de la función exponencial. Por ejemplo,

     λ> take 8 expx
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (derivada expx)
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (integral expx)
     [0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 expx

[1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 (derivada expx)

[1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 (integral expx)

[0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

Soluciones

type Serie a = [a] 

opuesta :: Num a => Serie a -> Serie a
opuesta = map negate

suma :: Num a => Serie a -> Serie a -> Serie a
suma = zipWith (+)

resta :: Num a => Serie a -> Serie a -> Serie a
resta xs ys = suma xs (opuesta ys)

producto :: Num a => Serie a -> Serie a -> Serie a
producto (x:xs) zs@(y:ys) = 
    x*y : suma (producto xs zs) (map (x*) ys)

cociente :: Fractional a => Serie a -> Serie a -> Serie a
cociente (x:xs) (y:ys) = zs 
    where zs = x/y : map (/y) (resta xs (producto zs ys))  

derivada :: (Num a, Enum a) => Serie a -> Serie a
derivada (_:xs) = zipWith (*) xs [1..]

integral :: (Fractional a, Enum a) => Serie a -> Serie a
integral xs = 0 : zipWith (/) xs [1..]

expx :: Serie Rational
expx = map (1/) (map fromIntegral factoriales)

-- factoriales es la lista de los factoriales. Por ejemplo, 
--    take 7 factoriales  ==  [1,1,2,6,24,120,720]
factoriales :: [Integer]
factoriales = 1 : scanl1 (*) [1..]

type Serie a = [a]

opuesta :: Num a => Serie a -> Serie a

opuesta = map negate

suma :: Num a => Serie a -> Serie a -> Serie a

suma = zipWith (+)

resta :: Num a => Serie a -> Serie a -> Serie a

resta xs ys = suma xs (opuesta ys)

producto :: Num a => Serie a -> Serie a -> Serie a

producto (x:xs) zs@(y:ys) =

x*y : suma (producto xs zs) (map (x*) ys)

cociente :: Fractional a => Serie a -> Serie a -> Serie a

cociente (x:xs) (y:ys) = zs

where zs = x/y : map (/y) (resta xs (producto zs ys))

derivada :: (Num a, Enum a) => Serie a -> Serie a

derivada (_:xs) = zipWith (*) xs [1..]

integral :: (Fractional a, Enum a) => Serie a -> Serie a

integral xs = 0 : zipWith (/) xs [1..]

expx :: Serie Rational

expx = map (1/) (map fromIntegral factoriales)

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

-- take 7 factoriales == [1,1,2,6,24,120,720]

factoriales :: [Integer]

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

Cálculo de pi mediante la fracción continua de Lange

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

En 1999, L.J. Lange publicó el artículo An elegant new continued fraction for π.

En el primer teorema del artículo se demuestra la siguiente expresión de π mediante una fracción continua
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange$

La primeras aproximaciones son

   a(1) = 3+1                = 4.0
   a(2) = 3+(1/(6+9))        = 3.066666666666667
   a(3) = 3+(1/(6+9/(6+25))) = 3.158974358974359

a(1) = 3+1 = 4.0

a(2) = 3+(1/(6+9)) = 3.066666666666667

a(3) = 3+(1/(6+9/(6+25))) = 3.158974358974359

Definir las funciones

   aproximacionPi :: Int -> Double
   grafica        :: [Int] -> IO ()

1 2	aproximacionPi :: Int -> Double grafica :: [Int] -> IO ()

tales que

(aproximacionPi n) es la n-ésima aproximación de pi con la fracción continua de Lange. Por ejemplo,

     aproximacionPi 1     ==  4.0
     aproximacionPi 2     ==  3.066666666666667
     aproximacionPi 3     ==  3.158974358974359
     aproximacionPi 10    ==  3.141287132741557
     aproximacionPi 100   ==  3.141592398533554
     aproximacionPi 1000  ==  3.1415926533392926

aproximacionPi 1 == 4.0

aproximacionPi 2 == 3.066666666666667

aproximacionPi 3 == 3.158974358974359

aproximacionPi 10 == 3.141287132741557

aproximacionPi 100 == 3.141592398533554

aproximacionPi 1000 == 3.1415926533392926

(grafica xs) dibuja la gráfica de las k-ésimas aproximaciones de pi donde k toma los valores de la lista xs. Por ejemplo, (grafica [1..10]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_2$
(grafica [10..100]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_3$
y (grafica [100..200]) dibuja
$Calculo_de_pi_mediante_la_fraccion_continua_de_Lange_4$

Soluciones

import Graphics.Gnuplot.Simple
 
-- fraccionPi es la representación de la fracción continua de pi como un
-- par de listas infinitas.
fraccionPi :: [(Integer, Integer)]
fraccionPi = zip (3 : [6,6..]) (map (^2) [1,3..])

-- (aproximacionFC n fc) es la n-ésima aproximación de la fracción
-- continua fc (como un par de listas).  
aproximacionFC :: Int -> [(Integer, Integer)] -> Double
aproximacionFC n =
  foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take n

aproximacionPi :: Int -> Double
aproximacionPi n =
  aproximacionFC n fraccionPi
 
grafica :: [Int] -> IO ()
grafica xs = 
  plotList [Key Nothing]
           [(k,aproximacionPi k) | k <- xs]

import Graphics.Gnuplot.Simple

-- fraccionPi es la representación de la fracción continua de pi como un

-- par de listas infinitas.

fraccionPi :: [(Integer, Integer)]

fraccionPi = zip (3 : [6,6..]) (map (^2) [1,3..])

-- (aproximacionFC n fc) es la n-ésima aproximación de la fracción

-- continua fc (como un par de listas).

aproximacionFC :: Int -> [(Integer, Integer)] -> Double

aproximacionFC n =

foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take n

aproximacionPi :: Int -> Double

aproximacionPi n =

aproximacionFC n fraccionPi

grafica :: [Int] -> IO ()

grafica xs =

plotList [Key Nothing]

[(k,aproximacionPi k) | k <- xs]

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>

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

PorJosé A. Alonso 9-abril-202026-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

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

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>

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>

Cálculo de pi mediante los métodos de Gregory-Leibniz y de Beeler

PorJosé A. Alonso 7-abril-202026-marzo-2022

La fórmula de Gregory-Leibniz para calcular pi es

y la de Beeler es

Definir las funciones

   aproximaPiGL     :: Int -> Double
   aproximaPiBeeler :: Int -> Double
   graficas         :: [Int] -> IO ()

aproximaPiGL :: Int -> Double

aproximaPiBeeler :: Int -> Double

graficas :: [Int] -> IO ()

tales que

(aproximaPiGL n) es la aproximación de pi con los primeros n términos de la fórmula de Gregory-Leibniz. Por ejemplo,

     aproximaPiGL 1       ==  4.0
     aproximaPiGL 2       ==  2.666666666666667
     aproximaPiGL 3       ==  3.466666666666667
     aproximaPiGL 10      ==  3.0418396189294032
     aproximaPiGL 100     ==  3.1315929035585537
     aproximaPiGL 1000    ==  3.140592653839794
     aproximaPiGL 10000   ==  3.1414926535900345
     aproximaPiGL 100000  ==  3.1415826535897198

aproximaPiGL 1 == 4.0

aproximaPiGL 2 == 2.666666666666667

aproximaPiGL 3 == 3.466666666666667

aproximaPiGL 10 == 3.0418396189294032

aproximaPiGL 100 == 3.1315929035585537

aproximaPiGL 1000 == 3.140592653839794

aproximaPiGL 10000 == 3.1414926535900345

aproximaPiGL 100000 == 3.1415826535897198

(aproximaPiBeeler n) es la aproximación de pi con los primeros n términos de la fórmula de Beeler. Por ejemplo,

     aproximaPiBeeler 1   ==  2.0
     aproximaPiBeeler 2   ==  2.6666666666666665
     aproximaPiBeeler 3   ==  2.933333333333333
     aproximaPiBeeler 10  ==  3.140578169680337
     aproximaPiBeeler 60  ==  3.141592653589793
     pi                   ==  3.141592653589793

aproximaPiBeeler 1 == 2.0

aproximaPiBeeler 2 == 2.6666666666666665

aproximaPiBeeler 3 == 2.933333333333333

aproximaPiBeeler 10 == 3.140578169680337

aproximaPiBeeler 60 == 3.141592653589793

pi == 3.141592653589793

(graficas xs) dibuja la gráfica de las k-ésimas aproximaciones de pi, donde k toma los valores de la lista xs, con las fórmulas de Gregory-Leibniz y de Beeler. Por ejemplo, (graficas [1..25]) dibuja

donde la línea morada corresponde a la aproximación de Gregory-Leibniz y la verde a la de Beeler.

Soluciones

import Graphics.Gnuplot.Simple

-- Definiciones de aproximaPiGL
-- ============================

-- 1ª definición de aproximaPiGL
aproximaPiGL :: Int -> Double
aproximaPiGL n = 4 * (sum . take n . sumaA . zipWith (/) [1,1..]) [1,3..]
  where sumaA (x:y:xs) = x:(-y):sumaA xs

-- 2ª definición de aproximaPiGL
aproximaPiGL2 :: Int -> Double
aproximaPiGL2 n =
  4 * (sum (take n (zipWith (/) (cycle [1,-1]) [1,3..])))

-- 3ª definición de aproximaPiGL
aproximaPiGL3 :: Int -> Double
aproximaPiGL3 n =
  4 * (sum . take n . zipWith (/) (cycle [1,-1])) [1,3..]

-- 4ª definición de aproximaPiGL
aproximaPiGL4 :: Int -> Double
aproximaPiGL4 n = serieGL !! (n-1)

serieGL :: [Double]
serieGL = scanl1 (+) (zipWith (/) numeradores denominadores)
  where numeradores   = cycle [4,-4]
        denominadores = [1,3..]

-- Definición de aproximaPiBeeler
aproximaPiBeeler :: Int -> Double
aproximaPiBeeler n = 2 * aux (fromIntegral n) 1
  where
    aux :: Double -> Double -> Double 
    aux n k | n == k    = 1
            | otherwise = 1 + (k/(2*k+1)) * aux n (1+k)

-- Definición de graficas
graficas :: [Int] -> IO ()
graficas xs = 
    plotLists [Key Nothing]
             [[(k,aproximaPiGL k)     | k <- xs],
              [(k,aproximaPiBeeler k) | k <- xs]]

import Graphics.Gnuplot.Simple

-- Definiciones de aproximaPiGL

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

-- 1ª definición de aproximaPiGL

aproximaPiGL :: Int -> Double

aproximaPiGL n = 4 * (sum . take n . sumaA . zipWith (/) [1,1..]) [1,3..]

where sumaA (x:y:xs) = x:(-y):sumaA xs

-- 2ª definición de aproximaPiGL

aproximaPiGL2 :: Int -> Double

aproximaPiGL2 n =

4 * (sum (take n (zipWith (/) (cycle [1,-1]) [1,3..])))

-- 3ª definición de aproximaPiGL

aproximaPiGL3 :: Int -> Double

aproximaPiGL3 n =

4 * (sum . take n . zipWith (/) (cycle [1,-1])) [1,3..]

-- 4ª definición de aproximaPiGL

aproximaPiGL4 :: Int -> Double

aproximaPiGL4 n = serieGL !! (n-1)

serieGL :: [Double]

serieGL = scanl1 (+) (zipWith (/) numeradores denominadores)

where numeradores = cycle [4,-4]

denominadores = [1,3..]

-- Definición de aproximaPiBeeler

aproximaPiBeeler :: Int -> Double

aproximaPiBeeler n = 2 * aux (fromIntegral n) 1

where

aux :: Double -> Double -> Double

aux n k | n == k = 1

| otherwise = 1 + (k/(2*k+1)) * aux n (1+k)

-- Definición de graficas

graficas :: [Int] -> IO ()

graficas xs =

plotLists [Key Nothing]

[[(k,aproximaPiGL k) | k <- xs],

[(k,aproximaPiBeeler k) | k <- xs]]

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>

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 menos conocida de las conjeturas de Goldbach

PorJosé A. Alonso 7-febrero-202017-febrero-2020

Goldbach, el de la famosa conjetura, hizo por lo menos otra conjetura que finalmente resultó ser falsa.

Esta última decía que todo número compuesto impar puede expresarse como la suma de un número primo más dos veces la suma de un cuadrado. Así por ejemplo,

    9 =  7 + 2×1^2
   15 =  7 + 2×2^2
   21 =  3 + 2×3^2
   25 =  7 + 2×3^2
   27 = 19 + 2×2^2
   33 = 31 + 2×1^2

9 = 7 + 2×1^2

15 = 7 + 2×2^2

21 = 3 + 2×3^2

25 = 7 + 2×3^2

27 = 19 + 2×2^2

33 = 31 + 2×1^2

Definir las sucesiones

   imparesCompuestos :: [Integer]
   descomposiciones :: Integer -> [(Integer,Integer)]
   contraejemplosGoldbach :: [Integer]

imparesCompuestos :: [Integer]

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

contraejemplosGoldbach :: [Integer]

tales que

imparesCompuestos es la lista de los números impares compuestos. Por ejemplo,

     take 9 imparesCompuestos  ==  [9,15,21,25,27,33,35,39,45]

1	take 9 imparesCompuestos == [9,15,21,25,27,33,35,39,45]

(descomposiciones n) es la lista de las descomposiciones de n de n como la suma de un número primo más dos veces la suma de un cuadrado. Por ejemplo,

     descomposiciones 9     ==  [(7,1)]
     descomposiciones 21    ==  [(3,9),(13,4),(19,1)]
     descomposiciones 5777  ==  []

descomposiciones 9 == [(7,1)]

descomposiciones 21 == [(3,9),(13,4),(19,1)]

descomposiciones 5777 == []

Las 3 descomposiciones de 21 son

     21 =  3 + 2*9 = 21 + 2*3^2
     21 = 13 + 2*4 = 13 + 2*3^2
     21 = 19 + 2*1 = 19 + 2*1^2

21 = 3 + 2*9 = 21 + 2*3^2

21 = 13 + 2*4 = 13 + 2*3^2

21 = 19 + 2*1 = 19 + 2*1^2

contraejemplosGoldbach es la lista de los contraejemplos de la anterior conjetura de Goldbach; es decir, los números impares compuestos que no pueden expresarse como la suma de un número primo más dos veces la suma de un cuadrado. Por ejemplo,

   take 2 contraejemplosGoldbach  ==  [5777,5993]

1	take 2 contraejemplosGoldbach == [5777,5993]

Comprobar con QuickCheck que la conjetura de Golbach se verifica a partir de 5993; es decir, todo número compuesto impar mayor que 5993 puede expresarse como la suma de un número primo más dos veces la suma de un cuadrado.

Nota: Basado en el artículo La menos conocida de las conjeturas de Goldbach de Claudio Meller en el blog Números y algo más.

Soluciones

import Data.Numbers.Primes
import Test.QuickCheck

imparesCompuestos :: [Integer]
imparesCompuestos = filter esCompuesto [3,5..]

-- (esCompuesto x) se verifica si x es un número compuesto. Por ejemplo,
--    esCompuesto 6  ==  True
--    esCompuesto 7  ==  False
esCompuesto :: Integer -> Bool
esCompuesto = not . isPrime

contraejemplosGoldbach :: [Integer]
contraejemplosGoldbach = filter esContraejemplo imparesCompuestos

-- (esContraejemplo x) es verifica si el número impar compuesto x es un
-- contraejemplo de la conjetura de Goldbach. Por ejemplo,
--    esContraejemplo 5777  ==  True
--    esContraejemplo 15    ==  False
esContraejemplo :: Integer -> Bool
esContraejemplo = null . descomposiciones

descomposiciones :: Integer -> [(Integer,Integer)]
descomposiciones n =
  [(p,x) | p <- takeWhile (<=n) primes
         , (n - p) `mod` 2 == 0
         , let x = (n - p) `div` 2
         , esCuadrado x]

-- (esCuadrado x) es verifica si x es un cuadrado perfecto. Por ejemplo, 
--    esCuadrado 16  ==  True
--    esCuadrado 27  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = y^2 == x
  where y = ceiling (sqrt (fromIntegral x))

-- La propiedad es
prop_conjetura :: Int -> Property
prop_conjetura n =
  n >= 0 ==> not (esContraejemplo (imparesCompuestosMayore5993 !! n))
  where imparesCompuestosMayore5993 = dropWhile (<=5993) imparesCompuestos

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

import Data.Numbers.Primes

import Test.QuickCheck

imparesCompuestos :: [Integer]

imparesCompuestos = filter esCompuesto [3,5..]

-- (esCompuesto x) se verifica si x es un número compuesto. Por ejemplo,

-- esCompuesto 6 == True

-- esCompuesto 7 == False

esCompuesto :: Integer -> Bool

esCompuesto = not . isPrime

contraejemplosGoldbach :: [Integer]

contraejemplosGoldbach = filter esContraejemplo imparesCompuestos

-- (esContraejemplo x) es verifica si el número impar compuesto x es un

-- contraejemplo de la conjetura de Goldbach. Por ejemplo,

-- esContraejemplo 5777 == True

-- esContraejemplo 15 == False

esContraejemplo :: Integer -> Bool

esContraejemplo = null . descomposiciones

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

descomposiciones n =

[(p,x) | p <- takeWhile (<=n) primes

, (n - p) `mod` 2 == 0

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

, esCuadrado x]

-- (esCuadrado x) es verifica si x es un cuadrado perfecto. Por ejemplo,

-- esCuadrado 16 == True

-- esCuadrado 27 == False

esCuadrado :: Integer -> Bool

esCuadrado x = y^2 == x

where y = ceiling (sqrt (fromIntegral x))

-- La propiedad es

prop_conjetura :: Int -> Property

prop_conjetura n =

n >= 0 ==> not (esContraejemplo (imparesCompuestosMayore5993 !! n))

where imparesCompuestosMayore5993 = dropWhile (<=5993) imparesCompuestos

-- La comprobación es

-- λ> quickCheck prop_conjetura

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

«Obvio es la palabra más peligrosa de las matemáticas.»

Eric Temple Bell

Entre dos potencias sucesivas

PorJosé A. Alonso 3-febrero-202010-febrero-2020

Se dice que un número entero está entre potencias sucesivas de n si x-1 es una potencia n-ésima y x+1 es una potencia (n+1)-ésima; es decir, si existen a y b tales que x-1 es a^n y x+1 es b^(n+1). Por ejemplo,

 2 está entre potencias sucesivas de 0, ya que  1 =  1^0 y  3 = 3^1
15 está entre potencias sucesivas de 1, ya que 14 = 14^1 y 16 = 4^2
26 está entre potencias sucesivas de 2, ya que 25 =  5^2 y 27 = 3^3

2 está entre potencias sucesivas de 0, ya que 1 = 1^0 y 3 = 3^1

15 está entre potencias sucesivas de 1, ya que 14 = 14^1 y 16 = 4^2

26 está entre potencias sucesivas de 2, ya que 25 = 5^2 y 27 = 3^3

Definir las funciones

   entrePotencias :: Integer -> Integer -> Bool
   pares :: [(Integer,Integer)]
   paresEntrePotencias :: [(Integer,Integer)]

entrePotencias :: Integer -> Integer -> Bool

pares :: [(Integer,Integer)]

paresEntrePotencias :: [(Integer,Integer)]

tales que

(entrePotencias n x) se verifica si x está entre potencias sucesivas de n. Por ejemplo,

     entrePotencias 0 2   ==  True
     entrePotencias 1 15  ==  True
     entrePotencias 2 26  ==  True

entrePotencias 0 2 == True

entrePotencias 1 15 == True

entrePotencias 2 26 == True

pares es la lista de los números enteros ordenados por su suma y primer elemento. Por ejemplo,

     λ> take 11 pares
     [(0,0),(0,1),(1,0),(0,2),(1,1),(2,0),(0,3),(1,2),(2,1),(3,0),(0,4)]

1 2	λ> take 11 pares [(0,0),(0,1),(1,0),(0,2),(1,1),(2,0),(0,3),(1,2),(2,1),(3,0),(0,4)]

paresEntrePotencias es la lista de los pares (n,x) tales que x está entre potencias sucesivas de n. Por ejemplo,

     λ> take 10 paresEntrePotencias
     [(1,0),(0,2),(1,3),(1,8),(1,15),(1,24),(2,26),(1,35),(1,48),(1,63)]

1 2	λ> take 10 paresEntrePotencias [(1,0),(0,2),(1,3),(1,8),(1,15),(1,24),(2,26),(1,35),(1,48),(1,63)]

Comprobar con QuickCheck que 26 es el único número que está entre potencias sucesivas con exponentes mayor que 1; es decir, que el único par (n,x) tal que x está entre potencias sucesivas de n con n mayor que uno es el (2,26).

Nota: Este ejercicio ha sido propuesto por Rebeca Isabel González Gordillo y está basado en el artículo El número 26 … ¡un número especial! de Amadeo Artacho en MatematicasCercanas.

Soluciones

import Test.QuickCheck

entrePotencias :: Integer -> Integer -> Bool
entrePotencias n x = esPotencia n (x-1) && esPotencia (n+1) (x+1)

-- (esPotencia n x) se verifica si x es una potencia n-ésima. Por
-- ejemplo, 
--    esPotencia 3 27  ==  True
--    esPotencia 3 25  ==  False
esPotencia :: Integer -> Integer -> Bool
esPotencia n x = x == r^n
  where r = ceiling ((fromIntegral x)**(1/fromIntegral n))

pares :: [(Integer,Integer)]
pares = [(x,n-x) | n <- [0..], x <- [0..n]]

paresEntrePotencias :: [(Integer,Integer)]
paresEntrePotencias =
  [(n,x) | (n,x) <- pares
         , entrePotencias n x]

prop_entrePotencias :: Integer -> Integer -> Property
prop_entrePotencias n x =
  n > 1 ==> entrePotencias n x == ((n,x) == (2,26))

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

import Test.QuickCheck

entrePotencias :: Integer -> Integer -> Bool

entrePotencias n x = esPotencia n (x-1) && esPotencia (n+1) (x+1)

-- (esPotencia n x) se verifica si x es una potencia n-ésima. Por

-- ejemplo,

-- esPotencia 3 27 == True

-- esPotencia 3 25 == False

esPotencia :: Integer -> Integer -> Bool

esPotencia n x = x == r^n

where r = ceiling ((fromIntegral x)**(1/fromIntegral n))

pares :: [(Integer,Integer)]

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

paresEntrePotencias :: [(Integer,Integer)]

paresEntrePotencias =

[(n,x) | (n,x) <- pares

, entrePotencias n x]

prop_entrePotencias :: Integer -> Integer -> Property

prop_entrePotencias n x =

n > 1 ==> entrePotencias n x == ((n,x) == (2,26))

-- La comprobación es

-- λ> quickCheck prop_entrePotencias

-- +++ 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 verdadero objetivo de la ciencia es el honor de la mente humana.»

Carl Gustav Jacob Jacobi

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.

Productos de sumas de cuatro cuadrados

PorJosé A. Alonso 24-enero-202031-enero-2020

Definir la función

   productoSuma4Cuadrados :: Integral a => [a] -> [a] -> [a] -> [a] -> a

1	productoSuma4Cuadrados :: Integral a => [a] -> [a] -> [a] -> [a] -> a

tal que (productoSuma4Cuadrados as bs cs ds) es el producto de las sumas de los cuadrados de cada una de las listas que ocupan la misma posición (hasta que alguna se acaba). Por ejemplo,

   productoSuma4Cuadrados [2,3] [1,5] [4,6] [0,3,9]
   = (2² + 1² + 4² + 0²) * (3² + 5² + 6² + 3²)
   = (4 +  1 + 16  + 0)  * (9 + 25 + 36  + 9)
   = 1659

productoSuma4Cuadrados [2,3] [1,5] [4,6] [0,3,9]

= (2² + 1² + 4² + 0²) * (3² + 5² + 6² + 3²)

= (4 + 1 + 16 + 0) * (9 + 25 + 36 + 9)

= 1659

Comprobar con QuickCheckWith que si as, bs cs y ds son listas no vacías de enteros positivos, entonces (productoSuma4Cuadrados as bs cs ds) se puede escribir como la suma de los cuadrados de cuatro enteros positivos.

Soluciones

import Data.List (zip4)
import Test.QuickCheck

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

productoSuma4Cuadrados :: Integral a => [a] -> [a] -> [a] -> [a] -> a
productoSuma4Cuadrados (a:as) (b:bs) (c:cs) (d:ds) =
  (a^2+b^2+c^2+d^2) * productoSuma4Cuadrados as bs cs ds
productoSuma4Cuadrados _ _ _ _ = 1

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

productoSuma4Cuadrados2 :: Integral a => [a] -> [a] -> [a] -> [a] -> a
productoSuma4Cuadrados2 as bs cs ds =
  product [a^2 + b^2 + c^2 + d^2 | (a,b,c,d) <- zip4 as bs cs ds]

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

-- La propiedad es
prop_productoSuma4Cuadrados ::
  [Integer] -> [Integer] -> [Integer] -> [Integer] -> Property
prop_productoSuma4Cuadrados as bs cs ds =
  all (not . null) [as, bs, cs, ds]
  ==> 
  esSuma4Cuadrados (productoSuma4Cuadrados as' bs' cs' ds')
  where as' = [1 + abs a | a <- as]
        bs' = [1 + abs b | b <- bs]
        cs' = [1 + abs c | c <- cs]
        ds' = [1 + abs d | d <- ds]

-- (esSuma4Cuadrados n) se verifica si n es la suma de 4 cuadrados. Por
-- ejemplo, 
--    esSuma4Cuadrados 42  ==  True
--    esSuma4Cuadrados 11  ==  False
--    esSuma4Cuadrados 41  ==  False
esSuma4Cuadrados :: Integer -> Bool
esSuma4Cuadrados = not . null . sumas4Cuadrados

-- (sumas4Cuadrados n) es la lista de las descomposiciones de n como
-- sumas de 4 cuadrados. Por ejemplo,
--    sumas4Cuadrados 42  ==  [(16,16,9,1),(25,9,4,4),(36,4,1,1)]
sumas4Cuadrados :: Integer -> [(Integer,Integer,Integer,Integer)]
sumas4Cuadrados n =
  [(a^2,b^2,c^2,d) | a <- [1 .. floor (sqrt (fromIntegral n / 4))]
                   , b <- [a .. floor (sqrt (fromIntegral (n-a^2) / 3))]
                   , c <- [b .. floor (sqrt (fromIntegral (n-a^2-b^2) / 2))]
                   , let d = n - a^2 - b^2 - c^2
                   , c^2 <= d 
                   , esCuadrado d]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer
raiz x = floor (sqrt (fromIntegral x))

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=5}) prop_productoSuma4Cuadrados
--    +++ OK, passed 100 tests.

import Data.List (zip4)

import Test.QuickCheck

-- 1ª solución

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

productoSuma4Cuadrados :: Integral a => [a] -> [a] -> [a] -> [a] -> a

productoSuma4Cuadrados (a:as) (b:bs) (c:cs) (d:ds) =

(a^2+b^2+c^2+d^2) * productoSuma4Cuadrados as bs cs ds

productoSuma4Cuadrados _ _ _ _ = 1

-- 2ª solución

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

productoSuma4Cuadrados2 :: Integral a => [a] -> [a] -> [a] -> [a] -> a

productoSuma4Cuadrados2 as bs cs ds =

product [a^2 + b^2 + c^2 + d^2 | (a,b,c,d) <- zip4 as bs cs ds]

-- Propiedad

-- =========

-- La propiedad es

prop_productoSuma4Cuadrados ::

[Integer] -> [Integer] -> [Integer] -> [Integer] -> Property

prop_productoSuma4Cuadrados as bs cs ds =

all (not . null) [as, bs, cs, ds]

==>

esSuma4Cuadrados (productoSuma4Cuadrados as' bs' cs' ds')

where as' = [1 + abs a | a <- as]

bs' = [1 + abs b | b <- bs]

cs' = [1 + abs c | c <- cs]

ds' = [1 + abs d | d <- ds]

-- (esSuma4Cuadrados n) se verifica si n es la suma de 4 cuadrados. Por

-- ejemplo,

-- esSuma4Cuadrados 42 == True

-- esSuma4Cuadrados 11 == False

-- esSuma4Cuadrados 41 == False

esSuma4Cuadrados :: Integer -> Bool

esSuma4Cuadrados = not . null . sumas4Cuadrados

-- (sumas4Cuadrados n) es la lista de las descomposiciones de n como

-- sumas de 4 cuadrados. Por ejemplo,

-- sumas4Cuadrados 42 == [(16,16,9,1),(25,9,4,4),(36,4,1,1)]

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

sumas4Cuadrados n =

[(a^2,b^2,c^2,d) | a <- [1 .. floor (sqrt (fromIntegral n / 4))]

, b <- [a .. floor (sqrt (fromIntegral (n-a^2) / 3))]

, c <- [b .. floor (sqrt (fromIntegral (n-a^2-b^2) / 2))]

, let d = n - a^2 - b^2 - c^2

, c^2 <= d

, esCuadrado d]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz 25 == 5

-- raiz 24 == 4

-- raiz 26 == 5

raiz :: Integer -> Integer

raiz x = floor (sqrt (fromIntegral x))

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=5}) prop_productoSuma4Cuadrados

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

¿Vivir? Sencillamente:
la sed y el agua cerca …
o el agua lejos, más, la sed y el agua,
un poco de cansancio ¡y a beberla!.

Antonio Machado

Sumas de cuatro cuadrados

PorJosé A. Alonso 23-enero-202030-enero-2020

El número 42 es una suma de cuatro cuadrados de números enteros positivos ya que

   42 = 16 + 16 + 9 + 1 = 4² + 4² + 3² + 1²

1	42 = 16 + 16 + 9 + 1 = 4² + 4² + 3² + 1²

Definir las funciones

   sumas4Cuadrados :: Integer -> [(Integer,Integer,Integer,Integer)]
   graficaNumeroSumas4Cuadrados :: Integer -> IO ()

1 2	sumas4Cuadrados :: Integer -> [(Integer,Integer,Integer,Integer)] graficaNumeroSumas4Cuadrados :: Integer -> IO ()

tales que

(sumas4Cuadrados n) es la lista de las descompociones de n como suma de cuatro cuadrados. Por ejemplo,

     sumas4Cuadrados 42  ==  [(16,16,9,1),(25,9,4,4),(36,4,1,1)]
     sumas4Cuadrados 14  ==  []
     length (sumas4Cuadrados (5*10^4))  ==  260

sumas4Cuadrados 42 == [(16,16,9,1),(25,9,4,4),(36,4,1,1)]

sumas4Cuadrados 14 == []

length (sumas4Cuadrados (5*10^4)) == 260

(graficaNumeroSumas4Cuadrados n) dibuja la gráfica del número de descomposiciones en sumas de 4 cuadrados de los n primeros. Por ejemplo, (graficaNumeroSumas4Cuadrados 600) dibuja

Soluciones

Pensamiento

import Graphics.Gnuplot.Simple

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

sumas4Cuadrados :: Integer -> [(Integer,Integer,Integer,Integer)]
sumas4Cuadrados n =
  [(a^2,b^2,c^2,d) | a <- [1..n]
                   , b <- [a..n]
                   , c <- [b..n]
                   , let d = n - a^2 - b^2 - c^2
                   , c^2 <= d 
                   , esCuadrado d]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer
raiz x = floor (sqrt (fromIntegral x))

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

-- Los intervalos de búqueda en la definición anterior se pueden reducir
-- teniendo en cuenta las siguientes restricciones
--    1 <= a <= b <= c <= d 
--    n = a² + b² + c² + d² >= 4a² ==> a <= sqrt (n/4)
--    n - a² = b² + c² + d² >= 3b² ==> b <= sqrt ((n-a²)/3)
--    n - a² - b² = c² + d² >= 2c² ==> c <= sqrt ((n-a²-b²)/2)

sumas4Cuadrados2 :: Integer -> [(Integer,Integer,Integer,Integer)]
sumas4Cuadrados2 n =
  [(a^2,b^2,c^2,d) | a <- [1 .. floor (sqrt (fromIntegral n / 4))]
                   , b <- [a .. floor (sqrt (fromIntegral (n-a^2) / 3))]
                   , c <- [b .. floor (sqrt (fromIntegral (n-a^2-b^2) / 2))]
                   , let d = n - a^2 - b^2 - c^2
                   , c^2 <= d 
                   , esCuadrado d]

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

-- La comprobación es
--    λ> length (sumas4Cuadrados 300)
--    11
--    (7.93 secs, 11,280,814,312 bytes)
--    λ> length (sumas4Cuadrados2 300)
--    11
--    (0.01 secs, 901,520 bytes)

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

graficaNumeroSumas4Cuadrados :: Integer -> IO ()
graficaNumeroSumas4Cuadrados n =
  plotList [ Key Nothing
           , Title "Numero de sumas como 4 cuadrados"
           , PNG "Sumas_de_cuatro_cuadrados.png"
           ]
           [length (sumas4Cuadrados2 k) | k <- [0..n]] 

-- Definición de esSuma4Cuadrados
-- ==============================

-- (esSuma4Cuadrados n) se verifica si n es la suma de 4 cuadrados. Por
-- ejemplo, 
--    esSuma4Cuadrados 42  ==  True
--    esSuma4Cuadrados 11  ==  False
--    esSuma4Cuadrados 41  ==  False
esSuma4Cuadrados :: Integer -> Bool
esSuma4Cuadrados = not . null . sumas4Cuadrados2

import Graphics.Gnuplot.Simple

-- 1ª definición de sumas4Cuadrados

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

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

sumas4Cuadrados n =

[(a^2,b^2,c^2,d) | a <- [1..n]

, b <- [a..n]

, c <- [b..n]

, let d = n - a^2 - b^2 - c^2

, c^2 <= d

, esCuadrado d]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz 25 == 5

-- raiz 24 == 4

-- raiz 26 == 5

raiz :: Integer -> Integer

raiz x = floor (sqrt (fromIntegral x))

-- 2ª definición de sumas4Cuadrados

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

-- Los intervalos de búqueda en la definición anterior se pueden reducir

-- teniendo en cuenta las siguientes restricciones

-- 1 <= a <= b <= c <= d

-- n = a² + b² + c² + d² >= 4a² ==> a <= sqrt (n/4)

-- n - a² = b² + c² + d² >= 3b² ==> b <= sqrt ((n-a²)/3)

-- n - a² - b² = c² + d² >= 2c² ==> c <= sqrt ((n-a²-b²)/2)

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

sumas4Cuadrados2 n =

[(a^2,b^2,c^2,d) | a <- [1 .. floor (sqrt (fromIntegral n / 4))]

, b <- [a .. floor (sqrt (fromIntegral (n-a^2) / 3))]

, c <- [b .. floor (sqrt (fromIntegral (n-a^2-b^2) / 2))]

, let d = n - a^2 - b^2 - c^2

, c^2 <= d

, esCuadrado d]

-- Comparación de eficiencia

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

-- La comprobación es

-- λ> length (sumas4Cuadrados 300)

-- 11

-- (7.93 secs, 11,280,814,312 bytes)

-- λ> length (sumas4Cuadrados2 300)

-- 11

-- (0.01 secs, 901,520 bytes)

-- Definición de graficaConvergencia

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

graficaNumeroSumas4Cuadrados :: Integer -> IO ()

graficaNumeroSumas4Cuadrados n =

plotList [ Key Nothing

, Title "Numero de sumas como 4 cuadrados"

, PNG "Sumas_de_cuatro_cuadrados.png"

]

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

-- Definición de esSuma4Cuadrados

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

-- (esSuma4Cuadrados n) se verifica si n es la suma de 4 cuadrados. Por

-- ejemplo,

-- esSuma4Cuadrados 42 == True

-- esSuma4Cuadrados 11 == False

-- esSuma4Cuadrados 41 == False

esSuma4Cuadrados :: Integer -> Bool

esSuma4Cuadrados = not . null . sumas4Cuadrados2

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>

¿Cuál es el peor de todos
los afanes? Preguntar.
¿Y el mejor? – Hacer camino
sin volver la vista atrás.

Antonio Machado

Elementos múltiplos de la longitud de la lista

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

Definir las funciones

   multiplosDeLaLongitud :: [Int] -> [Int]
   multiplosDeLaLongitudDeConsecutivos :: Int -> Int -> [Int]

1 2	multiplosDeLaLongitud :: [Int] -> [Int] multiplosDeLaLongitudDeConsecutivos :: Int -> Int -> [Int]

tales que

(multiplosDeLaLongitud xs) es la lista de los elementos de xs que son múltiplos de la longitud de xs. Por ejemplo,

     multiplosDeLaLongitud [2,4,6,8] == [4,8]

1	multiplosDeLaLongitud [2,4,6,8] == [4,8]

(multiplosDeLaLongitudDeConsecutivos n m) es la lista de elementos de [n..n+m-1] que son múltiplos de n. Por ejemplo,

     multiplosDeLaLongitudDeConsecutivos 9 2  ==  [10]
     multiplosDeLaLongitudDeConsecutivos 9 12 ==  [12]

1 2	multiplosDeLaLongitudDeConsecutivos 9 2 == [10] multiplosDeLaLongitudDeConsecutivos 9 12 == [12]

Comprobar con QuickCheck si se verifican las siguientes propiedades

En cualquier conjunto de m elementos consecutivos, m divide exactamente a uno de dichos elementos. En otras palabras, si n y m son enteros positivos, entonces (multiplosDeLaLongitudDeConsecutivos n m) tiene exactamente un elemento.
Si n es un entero positivo y m >= n, entonces (multiplosDeLaLongitudDeConsecutivos n m) es igual a [m]
Si n y n son enteros positivos y m < n, entonces (multiplosDeLaLongitudDeConsecutivos n m) es igual a [m * ceiling (n’ / m’)] donde n’ y m’ son las formas decimales de n y m respectivamente.

Soluciones

import Test.QuickCheck

multiplosDeLaLongitud :: [Int] -> [Int]
multiplosDeLaLongitud xs =
  [x | x <- xs
     , x `mod`  n == 0]
  where n = length xs

multiplosDeLaLongitudDeConsecutivos :: Int -> Int -> [Int]
multiplosDeLaLongitudDeConsecutivos n m =
  multiplosDeLaLongitud [n..n+m-1]

-- La 1ª propiedad es
prop_multiplosDeLaLongitud :: Int -> Int -> Property
prop_multiplosDeLaLongitud n m =
  n > 0 && m > 0
  ==>
  length (multiplosDeLaLongitudDeConsecutivos n m) == 1

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

-- La 2ª propiedad es
prop_multiplosDeLaLongitud2 :: Int -> Int -> Property
prop_multiplosDeLaLongitud2 n m =
  n > 0 && m >= n
  ==>
  multiplosDeLaLongitudDeConsecutivos n m == [m]

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

-- La 3ª propiedad es
prop_multiplosDeLaLongitud3 :: Int -> Int -> Property
prop_multiplosDeLaLongitud3 n m =
  n > 0 && 0 < m && m < n
  ==>
  multiplosDeLaLongitudDeConsecutivos n m == [m * ceiling (n' / m')]
  where n' = fromIntegral n
        m' = fromIntegral m 
  
-- La comprobación es
--    λ> quickCheck prop_multiplosDeLaLongitud3
--    +++ OK, passed 100 tests.

-- Con las propiedades anteriores se puede redefinir multiplos
multiplosDeLaLongitudDeConsecutivos2 :: Int -> Int -> [Int]
multiplosDeLaLongitudDeConsecutivos2 n m
  | m < n     = [m * ceiling (n' / m')]
  | otherwise = [m]
  where n' = fromIntegral n
        m' = fromIntegral m 

-- Comprobación de la equivalencia
prop_equiv :: Int -> Int -> Property
prop_equiv n m =
  n > 0 && m > 0
  ==>
  multiplosDeLaLongitudDeConsecutivos  n m ==
  multiplosDeLaLongitudDeConsecutivos2 n m

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

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

--    λ> xs1 = [multiplosDeLaLongitudDeConsecutivos n (10^3) | n <- [1..10^4]]
--    (0.01 secs, 0 bytes)
--    λ> xs2 = [multiplosDeLaLongitudDeConsecutivos2 n (10^3) | n <- [1..10^4]]
--    (0.01 secs, 0 bytes)
--    λ> maximum xs1
--    [10000]
--    (2.51 secs, 2,093,241,168 bytes)
--    λ> maximum xs2
--    [10000]
--    (0.02 secs, 16,742,408 bytes)

import Test.QuickCheck

multiplosDeLaLongitud :: [Int] -> [Int]

multiplosDeLaLongitud xs =

[x | x <- xs

, x `mod` n == 0]

where n = length xs

multiplosDeLaLongitudDeConsecutivos :: Int -> Int -> [Int]

multiplosDeLaLongitudDeConsecutivos n m =

multiplosDeLaLongitud [n..n+m-1]

-- La 1ª propiedad es

prop_multiplosDeLaLongitud :: Int -> Int -> Property

prop_multiplosDeLaLongitud n m =

n > 0 && m > 0

==>

length (multiplosDeLaLongitudDeConsecutivos n m) == 1

-- La comprobación es

-- λ> quickCheck prop_multiplosDeLaLongitud

-- +++ OK, passed 100 tests.

-- La 2ª propiedad es

prop_multiplosDeLaLongitud2 :: Int -> Int -> Property

prop_multiplosDeLaLongitud2 n m =

n > 0 && m >= n

==>

multiplosDeLaLongitudDeConsecutivos n m == [m]

-- La comprobación es

-- λ> quickCheck prop_multiplosDeLaLongitud2

-- +++ OK, passed 100 tests.

-- La 3ª propiedad es

prop_multiplosDeLaLongitud3 :: Int -> Int -> Property

prop_multiplosDeLaLongitud3 n m =

n > 0 && 0 < m && m < n

==>

multiplosDeLaLongitudDeConsecutivos n m == [m * ceiling (n' / m')]

where n' = fromIntegral n

m' = fromIntegral m

-- La comprobación es

-- λ> quickCheck prop_multiplosDeLaLongitud3

-- +++ OK, passed 100 tests.

-- Con las propiedades anteriores se puede redefinir multiplos

multiplosDeLaLongitudDeConsecutivos2 :: Int -> Int -> [Int]

multiplosDeLaLongitudDeConsecutivos2 n m

| m < n = [m * ceiling (n' / m')]

| otherwise = [m]

where n' = fromIntegral n

m' = fromIntegral m

-- Comprobación de la equivalencia

prop_equiv :: Int -> Int -> Property

prop_equiv n m =

n > 0 && m > 0

==>

multiplosDeLaLongitudDeConsecutivos n m ==

multiplosDeLaLongitudDeConsecutivos2 n m

-- La comprobación es

-- λ> quickCheck prop_multiplosDeLaLongitud

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> xs1 = [multiplosDeLaLongitudDeConsecutivos n (10^3) | n <- [1..10^4]]

-- (0.01 secs, 0 bytes)

-- λ> xs2 = [multiplosDeLaLongitudDeConsecutivos2 n (10^3) | n <- [1..10^4]]

-- (0.01 secs, 0 bytes)

-- λ> maximum xs1

-- [10000]

-- (2.51 secs, 2,093,241,168 bytes)

-- λ> maximum xs2

-- [10000]

-- (0.02 secs, 16,742,408 bytes)

Referencia

Set of successive numbers contains unique multiple en ProofWiki.

Pensamiento

Pensando que no veía
porque Dios no le miraba,
dijo Abel cuando moría:
Se acabó lo que se daba.

Antonio Machado

Teorema de existencia de divisores

PorJosé A. Alonso 10-enero-20206-junio-2022

El teorema de existencia de divisores afirma que

En cualquier subconjunto de {1, 2, …, 2m} con al menos m+1 elementos existen números distintos a, b tales que a divide a b.

Un conjunto de números naturales xs es mayoritario si existe un m tal que la lista de xs es un subconjunto de {1,2,…,2m} con al menos m+1 elementos. Por ejemplo, {2,3,5,6} porque es un subconjunto de {1,2,…,6} con más de 3 elementos.

Definir las funciones

   divisoresMultiplos :: [Integer] -> [(Integer,Integer)]
   esMayoritario :: [Integer] -> Bool

1 2	divisoresMultiplos :: [Integer] -> [(Integer,Integer)] esMayoritario :: [Integer] -> Bool

tales que

(divisores xs) es la lista de pares de elementos distintos de (a,b) tales que a divide a b. Por ejemplo,

     divisoresMultiplos [2,3,5,6]  ==  [(2,6),(3,6)]
     divisoresMultiplos [2,3,5]    ==  []
     divisoresMultiplos [4..8]     ==  [(4,8)]

divisoresMultiplos [2,3,5,6] == [(2,6),(3,6)]

divisoresMultiplos [2,3,5] == []

divisoresMultiplos [4..8] == [(4,8)]

(esMayoritario xs) se verifica xs es mayoritario. Por ejemplo,

     esMayoritario [2,3,5,6]  ==  True
     esMayoritario [2,3,5]    ==  False

1 2	esMayoritario [2,3,5,6] == True esMayoritario [2,3,5] == False

Comprobar con QuickCheck el teorema de existencia de divisores; es decir, en cualquier conjunto mayoritario existen números distintos a, b tales que a divide a b. Para la comprobación se puede usar el siguiente generador de conjuntos mayoritarios

   mayoritario :: Gen [Integer]
   mayoritario = do
     m' <- arbitrary
     let m = 1 + abs m'
     xs' <- sublistOf [1..2*m] `suchThat` (\ys -> genericLength ys > m)
     return xs'

mayoritario :: Gen [Integer]

mayoritario = do

m' <- arbitrary

let m = 1 + abs m'

xs' <- sublistOf [1..2*m] `suchThat` (\ys -> genericLength ys > m)

return xs'

con lo que la propiedad que hay que comprobar con QuickCheck es

   teorema_de_existencia_de_divisores :: Property
   teorema_de_existencia_de_divisores =
     forAll mayoritario (not . null . divisoresMultiplos)

teorema_de_existencia_de_divisores :: Property

teorema_de_existencia_de_divisores =

forAll mayoritario (not . null . divisoresMultiplos)

Soluciones

import Data.List (genericLength)
import Test.QuickCheck

divisoresMultiplos :: [Integer] -> [(Integer,Integer)]
divisoresMultiplos xs =
  [(x,y) | x <- xs
         , y <- xs
         , y /= x
         , y `mod` x == 0]

esMayoritario :: [Integer] -> Bool
esMayoritario xs =
  not (null xs) && length xs > ceiling (n / 2) 
  where n = fromIntegral (maximum xs)

-- Comprobación del teorema
-- ========================

-- La propiedad es
teorema_de_existencia_de_divisores :: Property
teorema_de_existencia_de_divisores =
  forAll mayoritario (not . null . divisoresMultiplos)

-- mayoritario es un generador de conjuntos mayoritarios. Por ejemplo, 
--    λ> sample mayoritario
--    [1,2]
--    [2,5,7,8]
--    [1,2,8,10,14]
--    [3,8,11,12,13,15,18,19,22,23,25,26]
--    [1,3,4,6]
--    [3,6,9,11,12,14,17,19]
mayoritario :: Gen [Integer]
mayoritario = do
  m' <- arbitrary
  let m = 1 + abs m'
  xs' <- sublistOf [1..2*m] `suchThat` (\ys -> genericLength ys > m)
  return xs'

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

import Data.List (genericLength)

import Test.QuickCheck

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

divisoresMultiplos xs =

[(x,y) | x <- xs

, y <- xs

, y /= x

, y `mod` x == 0]

esMayoritario :: [Integer] -> Bool

esMayoritario xs =

not (null xs) && length xs > ceiling (n / 2)

where n = fromIntegral (maximum xs)

-- Comprobación del teorema

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

-- La propiedad es

teorema_de_existencia_de_divisores :: Property

teorema_de_existencia_de_divisores =

forAll mayoritario (not . null . divisoresMultiplos)

-- mayoritario es un generador de conjuntos mayoritarios. Por ejemplo,

-- λ> sample mayoritario

-- [1,2]

-- [2,5,7,8]

-- [1,2,8,10,14]

-- [3,8,11,12,13,15,18,19,22,23,25,26]

-- [1,3,4,6]

-- [3,6,9,11,12,14,17,19]

mayoritario :: Gen [Integer]

mayoritario = do

m' <- arbitrary

let m = 1 + abs m'

xs' <- sublistOf [1..2*m] `suchThat` (\ys -> genericLength ys > m)

return xs'

-- La comprobación es

-- λ> quickCheck teorema_de_existencia_de_divisores

-- +++ OK, passed 100 tests.

Pensamiento

Guiomar, Guiomar,
mírame en ti castigado:
reo de haberte creado,
ya no te puedo olvidar.

Antonio Machado

El teorema de Navidad de Fermat

PorJosé A. Alonso 20-diciembre-201927-diciembre-2019

El 25 de diciembre de 1640, en una carta a Mersenne, Fermat demostró la conjetura de Girard: todo primo de la forma 4n+1 puede expresarse de manera única como suma de dos cuadrados. Por eso es conocido como el Teorema de Navidad de Fermat.

Definir las funciones

   representaciones :: Integer -> [(Integer,Integer)]
   primosImparesConRepresentacionUnica :: [Integer]
   primos4nM1 :: [Integer]

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

primosImparesConRepresentacionUnica :: [Integer]

primos4nM1 :: [Integer]

tales que

(representaciones n) es la lista de pares de números naturales (x,y) tales que n = x^2 + y^2 con x <= y. Por ejemplo,

     representaciones  20           ==  [(2,4)]
     representaciones  25           ==  [(0,5),(3,4)]
     representaciones 325           ==  [(1,18),(6,17),(10,15)]
     representaciones 100000147984  ==  [(0,316228)]
     length (representaciones (10^10))    ==  6
     length (representaciones (4*10^12))  ==  7

representaciones 20 == [(2,4)]

representaciones 25 == [(0,5),(3,4)]

representaciones 325 == [(1,18),(6,17),(10,15)]

representaciones 100000147984 == [(0,316228)]

length (representaciones (10^10)) == 6

length (representaciones (4*10^12)) == 7

primosImparesConRepresentacionUnica es la lista de los números primos impares que se pueden escribir exactamente de una manera como suma de cuadrados de pares de números naturales (x,y) con x <= y. Por ejemplo,

     λ> take 20 primosImparesConRepresentacionUnica
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

1 2	λ> take 20 primosImparesConRepresentacionUnica [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

primos4nM1 es la lista de los números primos que se pueden escribir como uno más un múltiplo de 4 (es decir, que son congruentes con 1 módulo 4). Por ejemplo,

     λ> take 20 primos4nM1
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

1 2	λ> take 20 primos4nM1 [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

El teorema de Navidad de Fermat afirma que un número primo impar p se puede escribir exactamente de una manera como suma de dos cuadrados de números naturales p = x² + y^2 (con x <= y) si, y sólo si, p se puede escribir como uno más un múltiplo de 4 (es decir, que es congruente con 1 módulo 4).

Comprobar con QuickCheck el teorema de Navidad de Fermat; es decir, que para todo número n, los n-ésimos elementos de primosImparesConRepresentacionUnica y de primos4nM1 son iguales.

Soluciones

import Data.Numbers.Primes (primes)
import Test.QuickCheck

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

representaciones :: Integer -> [(Integer,Integer)]
representaciones n =
  [(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]

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

representaciones2 :: Integer -> [(Integer,Integer)]
representaciones2 n =
  [(x,raiz z) | x <- [0..raiz (n `div` 2)] 
              , let z = n - x*x
              , esCuadrado z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer 
raiz 0 = 0
raiz 1 = 1
raiz x = aux (0,x)
    where aux (a,b) | d == x    = c
                    | c == a    = a
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c) 
              where c = (a+b) `div` 2
                    d = c^2

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

representaciones3 :: Integer -> [(Integer,Integer)]
representaciones3 n =
  [(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)] 
               , let z = n - x*x
               , esCuadrado3 z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado3 25  ==  True
--    esCuadrado3 26  ==  False
esCuadrado3 :: Integer -> Bool
esCuadrado3 x = x == y * y
  where y = raiz3 x

-- (raiz3 x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz3 25  ==  5
--    raiz3 24  ==  4
--    raiz3 26  ==  5
raiz3 :: Integer -> Integer
raiz3 x = floor (sqrt (fromIntegral x))
          
-- 4ª definición de representaciones
-- =================================
 
representaciones4 :: Integer -> [(Integer, Integer)]
representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))
  where aux x y
          | x > y     = [] 
          | otherwise = case compare (x*x + y*y) n of
                          LT -> aux (x + 1) y
                          EQ -> (x, y) : aux (x + 1) (y - 1)
                          GT -> aux x (y - 1)

-- Equivalencia de las definiciones de representaciones
-- ====================================================

-- La propiedad es
prop_representaciones_equiv :: (Positive Integer) -> Bool
prop_representaciones_equiv (Positive n) =
  representaciones  n == representaciones2 n &&
  representaciones2 n == representaciones3 n &&
  representaciones3 n == representaciones4 n

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

-- Comparación de eficiencia de las definiciones de representaciones
-- =================================================================

--    λ> representaciones 3025
--    [(0,55),(33,44)]
--    (2.86 secs, 1,393,133,528 bytes)
--    λ> representaciones2 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 867,944 bytes)
--    λ> representaciones3 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 173,512 bytes)
--    λ> representaciones4 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 423,424 bytes)
--    
--    λ> length (representaciones2 (10^10))
--    6
--    (3.38 secs, 2,188,903,544 bytes)
--    λ> length (representaciones3 (10^10))
--    6
--    (0.10 secs, 62,349,048 bytes)
--    λ> length (representaciones4 (10^10))
--    6
--    (0.11 secs, 48,052,360 bytes)
--
--    λ> length (representaciones3 (4*10^12))
--    7
--    (1.85 secs, 1,222,007,176 bytes)
--    λ> length (representaciones4 (4*10^12))
--    7
--    (1.79 secs, 953,497,480 bytes)

-- Definición de primosImparesConRepresentacionUnica
-- =================================================

primosImparesConRepresentacionUnica :: [Integer]
primosImparesConRepresentacionUnica =
  [x | x <- tail primes
     , length (representaciones4 x) == 1]

-- Definición de primos4nM1
-- ========================

primos4nM1 :: [Integer]
primos4nM1 = [x | x <- primes
                , x `mod` 4 == 1]

-- Teorema de Navidad de Fermat
-- ============================

-- La propiedad es
prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool
prop_teoremaDeNavidadDeFermat (Positive n) =
  primosImparesConRepresentacionUnica !! n == primos4nM1 !! n
             
-- La comprobación es
--    λ> quickCheck prop_teoremaDeNavidadDeFermat
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

import Data.Numbers.Primes (primes)

import Test.QuickCheck

-- 1ª definición de representaciones

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

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

representaciones n =

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

-- 2ª definición de representaciones

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

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

representaciones2 n =

[(x,raiz z) | x <- [0..raiz (n `div` 2)]

, let z = n - x*x

, esCuadrado z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz 25 == 5

-- raiz 24 == 4

-- raiz 26 == 5

raiz :: Integer -> Integer

raiz 0 = 0

raiz 1 = 1

raiz x = aux (0,x)

where aux (a,b) | d == x = c

| c == a = a

| d < x = aux (c,b)

| otherwise = aux (a,c)

where c = (a+b) `div` 2

d = c^2

-- 3ª definición de representaciones

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

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

representaciones3 n =

[(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)]

, let z = n - x*x

, esCuadrado3 z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado3 25 == True

-- esCuadrado3 26 == False

esCuadrado3 :: Integer -> Bool

esCuadrado3 x = x == y * y

where y = raiz3 x

-- (raiz3 x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz3 25 == 5

-- raiz3 24 == 4

-- raiz3 26 == 5

raiz3 :: Integer -> Integer

raiz3 x = floor (sqrt (fromIntegral x))

-- 4ª definición de representaciones

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

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

representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))

where aux x y

| x > y = []

| otherwise = case compare (x*x + y*y) n of

LT -> aux (x + 1) y

EQ -> (x, y) : aux (x + 1) (y - 1)

GT -> aux x (y - 1)

-- Equivalencia de las definiciones de representaciones

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

-- La propiedad es

prop_representaciones_equiv :: (Positive Integer) -> Bool

prop_representaciones_equiv (Positive n) =

representaciones n == representaciones2 n &&

representaciones2 n == representaciones3 n &&

representaciones3 n == representaciones4 n

-- La comprobación es

-- λ> quickCheck prop_representaciones_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de las definiciones de representaciones

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

-- λ> representaciones 3025

-- [(0,55),(33,44)]

-- (2.86 secs, 1,393,133,528 bytes)

-- λ> representaciones2 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 867,944 bytes)

-- λ> representaciones3 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 173,512 bytes)

-- λ> representaciones4 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 423,424 bytes)

-- λ> length (representaciones2 (10^10))

-- 6

-- (3.38 secs, 2,188,903,544 bytes)

-- λ> length (representaciones3 (10^10))

-- 6

-- (0.10 secs, 62,349,048 bytes)

-- λ> length (representaciones4 (10^10))

-- 6

-- (0.11 secs, 48,052,360 bytes)

-- λ> length (representaciones3 (4*10^12))

-- 7

-- (1.85 secs, 1,222,007,176 bytes)

-- λ> length (representaciones4 (4*10^12))

-- 7

-- (1.79 secs, 953,497,480 bytes)

-- Definición de primosImparesConRepresentacionUnica

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

primosImparesConRepresentacionUnica :: [Integer]

primosImparesConRepresentacionUnica =

[x | x <- tail primes

, length (representaciones4 x) == 1]

-- Definición de primos4nM1

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

primos4nM1 :: [Integer]

primos4nM1 = [x | x <- primes

, x `mod` 4 == 1]

-- Teorema de Navidad de Fermat

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

-- La propiedad es

prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool

prop_teoremaDeNavidadDeFermat (Positive n) =

primosImparesConRepresentacionUnica !! n == primos4nM1 !! n

-- La comprobación es

-- λ> quickCheck prop_teoremaDeNavidadDeFermat

-- +++ OK, passed 100 tests.

Pensamiento

Dijo Dios: brote la nada
Y alzó su mano derecha,
hasta ocultar su mirada.
Y quedó la nada hecha.

Antonio Machado

Suma de primos menores

PorJosé A. Alonso 9-diciembre-201917-diciembre-2019

La suma de los primos menores que 10 es 2 + 3 + 5 + 7 = 17.

Definir la función

   sumaPrimosMenores :: Integer -> Integer

1	sumaPrimosMenores :: Integer -> Integer

tal que (sumaPrimosMenores n) es la suma de los primos menores que n. Por ejemplo,

   sumaPrimosMenores 10        ==  17
   sumaPrimosMenores (5*10^5)  ==  9914236195

1 2	sumaPrimosMenores 10 == 17 sumaPrimosMenores (5*10^5) == 9914236195

Nota: Este ejercicio está basado en el problema 10 del Proyecto Euler

Soluciones

import Data.Numbers.Primes (primes)

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

sumaPrimosMenores :: Integer -> Integer
sumaPrimosMenores n = sum (takeWhile (<n) primos)

-- primos es la lista de los números primos. Por ejemplo,
--    λ> take 12 primos
--    [2,3,5,7,11,13,17,19,23,29,31,37]
primos :: [Integer]
primos = 2 : filter esPrimo [3,5..]

esPrimo :: Integer -> Bool
esPrimo n = null [x | x <- [2..(ceiling . sqrt . fromIntegral) n]
                    , n `mod` x == 0]

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

sumaPrimosMenores2 :: Integer -> Integer
sumaPrimosMenores2 n = sum (takeWhile (<n) primos2)

primos2 :: [Integer]
primos2 = 2 : filter esPrimo2 [3,5..]
 
esPrimo2 :: Integer -> Bool
esPrimo2 x =
  all ((/= 0) . mod x)
  (takeWhile (<= floor (sqrt (fromIntegral x))) primos2)

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

sumaPrimosMenores3 :: Integer -> Integer
sumaPrimosMenores3 n = sum (takeWhile (<n) primes)

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

--    λ> sumaPrimosMenores (2*10^5)
--    1709600813
--    (2.56 secs, 1,522,015,240 bytes)
--    λ> sumaPrimosMenores2 (2*10^5)
--    1709600813
--    (0.56 secs, 376,951,456 bytes)
--    λ> sumaPrimosMenores3 (2*10^5)
--    1709600813
--    (0.07 secs, 62,321,888 bytes)

import Data.Numbers.Primes (primes)

-- 1ª solución

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

sumaPrimosMenores :: Integer -> Integer

sumaPrimosMenores n = sum (takeWhile (<n) primos)

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

-- λ> take 12 primos

-- [2,3,5,7,11,13,17,19,23,29,31,37]

primos :: [Integer]

primos = 2 : filter esPrimo [3,5..]

esPrimo :: Integer -> Bool

esPrimo n = null [x | x <- [2..(ceiling . sqrt . fromIntegral) n]

, n `mod` x == 0]

-- 2ª solución

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

sumaPrimosMenores2 :: Integer -> Integer

sumaPrimosMenores2 n = sum (takeWhile (<n) primos2)

primos2 :: [Integer]

primos2 = 2 : filter esPrimo2 [3,5..]

esPrimo2 :: Integer -> Bool

esPrimo2 x =

all ((/= 0) . mod x)

(takeWhile (<= floor (sqrt (fromIntegral x))) primos2)

-- 3ª solución

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

sumaPrimosMenores3 :: Integer -> Integer

sumaPrimosMenores3 n = sum (takeWhile (<n) primes)

-- Comparación de eficiencia

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

-- λ> sumaPrimosMenores (2*10^5)

-- 1709600813

-- (2.56 secs, 1,522,015,240 bytes)

-- λ> sumaPrimosMenores2 (2*10^5)

-- 1709600813

-- (0.56 secs, 376,951,456 bytes)

-- λ> sumaPrimosMenores3 (2*10^5)

-- 1709600813

-- (0.07 secs, 62,321,888 bytes)

Pensamiento

El movimiento no es nada esencial. La fuerza puede ser inmóvil (lo es en su estado de pureza); mas no por ello deja de ser activa.

Antonio Machado

Mayor producto de n dígitos consecutivos de un número

PorJosé A. Alonso 29-noviembre-20196-diciembre-2019

Definir la función

   mayorProducto :: Int -> Integer -> Integer

1	mayorProducto :: Int -> Integer -> Integer

tal que (mayorProducto n x) es el mayor producto de n dígitos consecutivos del número x (suponiendo que x tiene al menos n dígitos). Por ejemplo,

   mayorProducto 2 325                  ==  10
   mayorProducto 5 11111                ==  1
   mayorProducto 5 113111               ==  3
   mayorProducto 5 110111               ==  0
   mayorProducto 5 10151112             ==  10
   mayorProducto 5 101511124            ==  10
   mayorProducto 5 (product [1..1000])  ==  41472

mayorProducto 2 325 == 10

mayorProducto 5 11111 == 1

mayorProducto 5 113111 == 3

mayorProducto 5 110111 == 0

mayorProducto 5 10151112 == 10

mayorProducto 5 101511124 == 10

mayorProducto 5 (product [1..1000]) == 41472

Nota: Este ejercicio está basado en el problema 8 del Proyecto Euler

Soluciones

import Data.List (inits, tails)
import Data.Char (digitToInt)

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

mayorProducto :: Int -> Integer -> Integer
mayorProducto n x =
  maximum [product xs | xs <- segmentos n (digitos x)]

-- (digitos x) es la lista de las digitos del número x. Por ejemplo, 
--    digitos 325  ==  [3,2,5]
digitos :: Integer -> [Integer]
digitos x = map (toInteger . digitToInt) (show x)

-- (segmentos n xs) es la lista de los segmentos de longitud n de la
-- lista xs. Por ejemplo,
--    segmentos 2 [3,5,4,6]  ==  [[3,5],[5,4],[4,6]]
segmentos :: Int -> [Integer] -> [[Integer]]
segmentos n xs = take (length xs - n + 1) (map (take n) (tails xs))

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

mayorProducto2 :: Int -> Integer -> Integer
mayorProducto2 n x = maximum (aux ns)
    where ns     = [read [d] | d <- show x]
          aux xs | length xs < n = []
                 | otherwise     = product (take n xs) : aux (tail xs)

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

mayorProducto3 :: Int -> Integer -> Integer
mayorProducto3 n = maximum
                 . map (product . take n)
                 . filter ((>=n) . length) 
                 . tails
                 . digitos

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

mayorProducto4 :: Int -> Integer -> Integer
mayorProducto4 n = maximum  
                 . map (product . map (fromIntegral . digitToInt)) 
                 . filter ((==n) . length) 
                 . concatMap inits
                 . tails 
                 . show

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

--    λ> mayorProducto 5 (product [1..500])
--    28224
--    (0.01 secs, 1,645,256 bytes)
--    λ> mayorProducto2 5 (product [1..500])
--    28224
--    (0.03 secs, 5,848,416 bytes)
--    λ> mayorProducto3 5 (product [1..500])
--    28224
--    (0.03 secs, 1,510,640 bytes)
--    λ> mayorProducto4 5 (product [1..500])
--    28224
--    (1.85 secs, 10,932,551,216 bytes)
--    
--    λ> mayorProducto 5 (product [1..7000])
--    46656
--    (0.10 secs, 68,590,808 bytes)
--    λ> mayorProducto2 5 (product [1..7000])
--    46656
--    (1.63 secs, 157,031,432 bytes)
--    λ> mayorProducto3 5 (product [1..7000])
--    46656
--    (1.55 secs, 65,727,176 bytes)

import Data.List (inits, tails)

import Data.Char (digitToInt)

-- 1ª solución

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

mayorProducto :: Int -> Integer -> Integer

mayorProducto n x =

maximum [product xs | xs <- segmentos n (digitos x)]

-- (digitos x) es la lista de las digitos del número x. Por ejemplo,

-- digitos 325 == [3,2,5]

digitos :: Integer -> [Integer]

digitos x = map (toInteger . digitToInt) (show x)

-- (segmentos n xs) es la lista de los segmentos de longitud n de la

-- lista xs. Por ejemplo,

-- segmentos 2 [3,5,4,6] == [[3,5],[5,4],[4,6]]

segmentos :: Int -> [Integer] -> [[Integer]]

segmentos n xs = take (length xs - n + 1) (map (take n) (tails xs))

-- 2ª solución

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

mayorProducto2 :: Int -> Integer -> Integer

mayorProducto2 n x = maximum (aux ns)

where ns = [read [d] | d <- show x]

aux xs | length xs < n = []

| otherwise = product (take n xs) : aux (tail xs)

-- 3ª solución

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

mayorProducto3 :: Int -> Integer -> Integer

mayorProducto3 n = maximum

. map (product . take n)

. filter ((>=n) . length)

. tails

. digitos

-- 4ª solución

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

mayorProducto4 :: Int -> Integer -> Integer

mayorProducto4 n = maximum

. map (product . map (fromIntegral . digitToInt))

. filter ((==n) . length)

. concatMap inits

. tails

. show

-- Comparación de eficiencia

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

-- λ> mayorProducto 5 (product [1..500])

-- 28224

-- (0.01 secs, 1,645,256 bytes)

-- λ> mayorProducto2 5 (product [1..500])

-- 28224

-- (0.03 secs, 5,848,416 bytes)

-- λ> mayorProducto3 5 (product [1..500])

-- 28224

-- (0.03 secs, 1,510,640 bytes)

-- λ> mayorProducto4 5 (product [1..500])

-- 28224

-- (1.85 secs, 10,932,551,216 bytes)

-- λ> mayorProducto 5 (product [1..7000])

-- 46656

-- (0.10 secs, 68,590,808 bytes)

-- λ> mayorProducto2 5 (product [1..7000])

-- 46656

-- (1.63 secs, 157,031,432 bytes)

-- λ> mayorProducto3 5 (product [1..7000])

-- 46656

-- (1.55 secs, 65,727,176 bytes)

Pensamiento

«El control de la complejidad es la esencia de la programación.» ~ B.W. Kernigan

Mayor divisor primo

PorJosé A. Alonso 26-noviembre-20193-diciembre-2019

Los divisores primos de 13195 son 5, 7, 13 y 29. Por tanto, el mayor divisor primo de 13195 es 29.

Definir la función

   mayorDivisorPrimo :: Integer -> Integer

1	mayorDivisorPrimo :: Integer -> Integer

tal que (mayorDivisorPrimo n) es el mayor divisor primo de n. Por ejemplo,

   mayorDivisorPrimo 13195            ==  29
   mayorDivisorPrimo 152416333181401  ==  12345701

1 2	mayorDivisorPrimo 13195 == 29 mayorDivisorPrimo 152416333181401 == 12345701

Nota: Este ejercicio está basado en el problema 3 del Proyecto Euler

Soluciones

import Data.Numbers.Primes (primeFactors)

-- 1ª solución  (sin librerías auxiliares)
-- =======================================

mayorDivisorPrimo :: Integer -> Integer
mayorDivisorPrimo = last . divisoresPrimos 

-- (divisoresPrimos n) es la lista de los divisores primos de n. Por
-- ejemplo, 
--    divisoresPrimos 13195  ==  [5,7,13,29]
divisoresPrimos :: Integer -> [Integer]
divisoresPrimos 0 = []
divisoresPrimos 1 = []
divisoresPrimos n = m : divisoresPrimos (n `div` m)
  where m = menorDivisorPrimo n 

-- (menorDivisorPrimo n) es el menor divisor primo de n. Por ejemplo, 
--    menorDivisorPrimo 24  ==  2
--    menorDivisorPrimo 25  ==  5
--    menorDivisorPrimo 29  ==  29
menorDivisorPrimo :: Integer -> Integer
menorDivisorPrimo x =
  head [y | y <- 2 : [3,5..(ceiling . sqrt . fromIntegral) x] ++ [x]
          , x `mod` y == 0]

-- 2ª solución (con la librería Data.Numbers.Primes)
-- =================================================

mayorDivisorPrimo2 :: Integer -> Integer
mayorDivisorPrimo2 = last . primeFactors

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

--   λ> mayorDivisorPrimo 152416333181401
--   12345701
--   (1.96 secs, 1,630,201,856 bytes)
--   λ> mayorDivisorPrimo2 152416333181401
--   12345701
--   (2.01 secs, 5,445,284,432 bytes)

import Data.Numbers.Primes (primeFactors)

-- 1ª solución (sin librerías auxiliares)

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

mayorDivisorPrimo :: Integer -> Integer

mayorDivisorPrimo = last . divisoresPrimos

-- (divisoresPrimos n) es la lista de los divisores primos de n. Por

-- ejemplo,

-- divisoresPrimos 13195 == [5,7,13,29]

divisoresPrimos :: Integer -> [Integer]

divisoresPrimos 0 = []

divisoresPrimos 1 = []

divisoresPrimos n = m : divisoresPrimos (n `div` m)

where m = menorDivisorPrimo n

-- (menorDivisorPrimo n) es el menor divisor primo de n. Por ejemplo,

-- menorDivisorPrimo 24 == 2

-- menorDivisorPrimo 25 == 5

-- menorDivisorPrimo 29 == 29

menorDivisorPrimo :: Integer -> Integer

menorDivisorPrimo x =

head [y | y <- 2 : [3,5..(ceiling . sqrt . fromIntegral) x] ++ [x]

, x `mod` y == 0]

-- 2ª solución (con la librería Data.Numbers.Primes)

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

mayorDivisorPrimo2 :: Integer -> Integer

mayorDivisorPrimo2 = last . primeFactors

-- Comparación de eficiencia

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

-- λ> mayorDivisorPrimo 152416333181401

-- 12345701

-- (1.96 secs, 1,630,201,856 bytes)

-- λ> mayorDivisorPrimo2 152416333181401

-- 12345701

-- (2.01 secs, 5,445,284,432 bytes)

Pensamiento

«Un programa de ordenador es una demostración.» ~ Igor Rivin

Operaciones con series de potencias

PorJosé A. Alonso 26-abril-201925-marzo-2022

Una serie de potencias es una serie de la forma

   a₀ + a₁x + a₂x² + a₃x³ + ...

1	a₀ + a₁x + a₂x² + a₃x³ + ...

Las series de potencias se pueden representar mediante listas infinitas. Por ejemplo, la serie de la función exponencial es

   e^x = 1 + x + x²/2! + x³/3! + ...

1	e^x = 1 + x + x²/2! + x³/3! + ...

y se puede representar por [1, 1, 1/2, 1/6, 1/24, 1/120, …]

Las operaciones con series se pueden ver como una generalización de las de los polinomios.

En lo que sigue, usaremos el tipo (Serie a) para representar las series de potencias con coeficientes en a y su definición es

   type Serie a = [a]

1	type Serie a = [a]

Definir las siguientes funciones

   opuesta      :: Num a => Serie a -> Serie a
   suma         :: Num a => Serie a -> Serie a -> Serie a
   resta        :: Num a => Serie a -> Serie a -> Serie a
   producto     :: Num a => Serie a -> Serie a -> Serie a
   cociente     :: Fractional a => Serie a -> Serie a -> Serie a
   derivada     :: (Num a, Enum a) => Serie a -> Serie a
   integral     :: (Fractional a, Enum a) => Serie a -> Serie a
   expx         :: Serie Rational

opuesta :: Num a => Serie a -> Serie a

suma :: Num a => Serie a -> Serie a -> Serie a

resta :: Num a => Serie a -> Serie a -> Serie a

producto :: Num a => Serie a -> Serie a -> Serie a

cociente :: Fractional a => Serie a -> Serie a -> Serie a

derivada :: (Num a, Enum a) => Serie a -> Serie a

integral :: (Fractional a, Enum a) => Serie a -> Serie a

expx :: Serie Rational

tales que

(opuesta xs) es la opuesta de la serie xs. Por ejemplo,

     λ> take 7 (opuesta [-6,-4..])
     [6,4,2,0,-2,-4,-6]

1 2	λ> take 7 (opuesta [-6,-4..]) [6,4,2,0,-2,-4,-6]

(suma xs ys) es la suma de las series xs e ys. Por ejemplo,

     λ> take 7 (suma [1,3..] [2,4..])
     [3,7,11,15,19,23,27]

1 2	λ> take 7 (suma [1,3..] [2,4..]) [3,7,11,15,19,23,27]

(resta xs ys) es la resta de las series xs es ys. Por ejemplo,

     λ> take 7 (resta [3,5..] [2,4..])
     [1,1,1,1,1,1,1]
     λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])
     [2,4,6,8,10,12,14]

λ> take 7 (resta [3,5..] [2,4..])

[1,1,1,1,1,1,1]

λ> take 7 (resta ([3,7,11,15,19,23,27] ++ repeat 0) [1,3..])

[2,4,6,8,10,12,14]

(producto xs ys) es el producto de las series xs e ys. Por ejemplo,

     λ> take 7 (producto [3,5..] [2,4..])
     [6,22,52,100,170,266,392]

1 2	λ> take 7 (producto [3,5..] [2,4..]) [6,22,52,100,170,266,392]

(cociente xs ys) es el cociente de las series xs e ys. Por ejemplo,

     λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..])
     [2.0,4.0,6.0,8.0,10.0,12.0,14.0]

1 2	λ> take 7 (cociente ([6,22,52,100,170,266,392] ++ repeat 0) [3,5..]) [2.0,4.0,6.0,8.0,10.0,12.0,14.0]

(derivada xs) es la derivada de la serie xs. Por ejemplo,

     λ> take 7 (derivada [2,4..])
     [4,12,24,40,60,84,112]

1 2	λ> take 7 (derivada [2,4..]) [4,12,24,40,60,84,112]

(integral xs) es la integral de la serie xs. Por ejemplo,

     λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0))
     [0.0,4.0,6.0,8.0,10.0,12.0,14.0]

1 2	λ> take 7 (integral ([4,12,24,40,60,84,112] ++ repeat 0)) [0.0,4.0,6.0,8.0,10.0,12.0,14.0]

expx es la serie de la función exponencial. Por ejemplo,

     λ> take 8 expx
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (derivada expx)
     [1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]
     λ> take 8 (integral expx)
     [0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 expx

[1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 (derivada expx)

[1 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

λ> take 8 (integral expx)

[0 % 1,1 % 1,1 % 2,1 % 6,1 % 24,1 % 120,1 % 720,1 % 5040]

Soluciones

type Serie a = [a] 

opuesta :: Num a => Serie a -> Serie a
opuesta = map negate

suma :: Num a => Serie a -> Serie a -> Serie a
suma = zipWith (+)

resta :: Num a => Serie a -> Serie a -> Serie a
resta xs ys = suma xs (opuesta ys)

producto :: Num a => Serie a -> Serie a -> Serie a
producto (x:xs) zs@(y:ys) = 
    x*y : suma (producto xs zs) (map (x*) ys)

cociente :: Fractional a => Serie a -> Serie a -> Serie a
cociente (x:xs) (y:ys) = zs 
    where zs = x/y : map (/y) (resta xs (producto zs ys))  

derivada :: (Num a, Enum a) => Serie a -> Serie a
derivada (_:xs) = zipWith (*) xs [1..]

integral :: (Fractional a, Enum a) => Serie a -> Serie a
integral xs = 0 : zipWith (/) xs [1..]

expx :: Serie Rational
expx = map (1/) (map fromIntegral factoriales)

-- factoriales es la lista de los factoriales. Por ejemplo, 
--    take 7 factoriales  ==  [1,1,2,6,24,120,720]
factoriales :: [Integer]
factoriales = 1 : scanl1 (*) [1..]

type Serie a = [a]

opuesta :: Num a => Serie a -> Serie a

opuesta = map negate

suma :: Num a => Serie a -> Serie a -> Serie a

suma = zipWith (+)

resta :: Num a => Serie a -> Serie a -> Serie a

resta xs ys = suma xs (opuesta ys)

producto :: Num a => Serie a -> Serie a -> Serie a

producto (x:xs) zs@(y:ys) =

x*y : suma (producto xs zs) (map (x*) ys)

cociente :: Fractional a => Serie a -> Serie a -> Serie a

cociente (x:xs) (y:ys) = zs

where zs = x/y : map (/y) (resta xs (producto zs ys))

derivada :: (Num a, Enum a) => Serie a -> Serie a

derivada (_:xs) = zipWith (*) xs [1..]

integral :: (Fractional a, Enum a) => Serie a -> Serie a

integral xs = 0 : zipWith (/) xs [1..]

expx :: Serie Rational

expx = map (1/) (map fromIntegral factoriales)

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

-- take 7 factoriales == [1,1,2,6,24,120,720]

factoriales :: [Integer]

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

Pensamiento

Ni mármol duro y eterno,
ni música ni pintura,
sino palabra en el tiempo.

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

Caminos minimales en un árbol numérico

PorJosé A. Alonso 19-marzo-201930-marzo-2019

En la librería Data.Tree se definen los tipos de árboles y bosques como sigue

   data Tree a   = Node a (Forest a)
   type Forest a = [Tree a]

1 2	data Tree a = Node a (Forest a) type Forest a = [Tree a]

Se pueden definir árboles. Por ejemplo,

   ej = Node 3 [Node 5 [Node 9 []], Node 7 []]

1	ej = Node 3 [Node 5 [Node 9 []], Node 7 []]

Y se pueden dibujar con la función drawTree. Por ejemplo,

   λ> putStrLn (drawTree (fmap show ej))
   3
   |
   +- 5
   |  |
   |  `- 9
   |
   `- 7

λ> putStrLn (drawTree (fmap show ej))

+- 5

| |

| `- 9

`- 7

Los mayores divisores de un número x son los divisores u tales que u > 1 y existe un v tal que 1 < v < u y u.v = x. Por ejemplo, los mayores divisores de 24 son 12, 8 y 6.

El árbol de los predecesores y mayores divisores de un número x es el árbol cuya raíz es x y los sucesores de cada nodo y > 1 es el conjunto formado por y-1 junto con los mayores divisores de y. Los nodos con valor 1 no tienen sucesores. Por ejemplo, el árbol de los predecesores y mayores divisores del número 6 es

/ \

5 3

| |

4 2

/ \ |

3 2 1

| |

2 1

Definir las siguientes funciones

   mayoresDivisores :: Int -> [Int]
   arbol            :: Int -> Tree Int
   caminos          :: Int -> [[Int]]
   caminosMinimales :: Int -> [[Int]]

mayoresDivisores :: Int -> [Int]

arbol :: Int -> Tree Int

caminos :: Int -> [[Int]]

caminosMinimales :: Int -> [[Int]]

tales que
+ (mayoresDivisores x) es la lista de los mayores divisores de x. Por ejemplo,

     mayoresDivisores 24  ==  [12,8,6]
     mayoresDivisores 16  ==  [8,4]
     mayoresDivisores 10  ==  [5]
     mayoresDivisores 17  ==  []

mayoresDivisores 24 == [12,8,6]

mayoresDivisores 16 == [8,4]

mayoresDivisores 10 == [5]

mayoresDivisores 17 == []

(arbol x) es el árbol de los predecesores y mayores divisores del número x. Por ejemplo,

     λ> putStrLn (drawTree (fmap show (arbol 6)))
     6
     |
     +- 5
     |  |
     |  `- 4
     |     |
     |     +- 3
     |     |  |
     |     |  `- 2
     |     |     |
     |     |     `- 1
     |     |
     |     `- 2
     |        |
     |        `- 1
     |
     `- 3
        |
        `- 2
           |
           `- 1

λ> putStrLn (drawTree (fmap show (arbol 6)))

+- 5

| |

| `- 4

| |

| +- 3

| | |

| | `- 2

| | |

| | `- 1

| |

| `- 2

| |

| `- 1

`- 3

`- 2

`- 1

(caminos x) es la lista de los caminos en el árbol de los predecesores y mayores divisores del número x. Por ejemplo,

     λ> caminos 6
     [[6,5,4,3,2,1],[6,5,4,2,1],[6,3,2,1]]

1 2	λ> caminos 6 [[6,5,4,3,2,1],[6,5,4,2,1],[6,3,2,1]]

(caminosMinimales x) es la lista de los caminos en de menor longitud en el árbol de los predecesores y mayores divisores del número x. Por ejemplo,

     λ> caminosMinimales 6
     [[6,3,2,1]]
     λ> caminosMinimales 17
     [[17,16,4,2,1]]
     λ> caminosMinimales 50
     [[50,25,5,4,2,1],[50,10,9,3,2,1],[50,10,5,4,2,1]]

λ> caminosMinimales 6

[[6,3,2,1]]

λ> caminosMinimales 17

[[17,16,4,2,1]]

λ> caminosMinimales 50

[[50,25,5,4,2,1],[50,10,9,3,2,1],[50,10,5,4,2,1]]

Soluciones

import Data.Tree
import Test.QuickCheck

mayoresDivisores :: Int -> [Int]
mayoresDivisores x =
  [max u v | u <- [2..floor (sqrt (fromIntegral x))]
           , x `mod` u == 0
           , let v = x `div` u]  

arbol :: Int -> Tree Int
arbol 1 = Node 1 []
arbol x = Node x (arbol (x-1) : [arbol y | y <- mayoresDivisores x])

caminos :: Int -> [[Int]]
caminos = caminosArbol . arbol

--    λ> caminosArbol (arbol 6)
--    [[6,5,4,3,2,1],[6,5,4,2,1],[6,3,2,1]]
caminosArbol :: Tree a -> [[a]]
caminosArbol (Node x []) = [[x]]
caminosArbol (Node x as) = [x:ys | ys <- caminosBosque as]

caminosBosque :: Forest a -> [[a]]
caminosBosque = concatMap caminosArbol

caminosMinimales :: Int -> [[Int]]
caminosMinimales x = [ys | ys <- yss, length ys == m]
  where yss = caminos x
        m   = minimum (map length yss)

import Data.Tree

import Test.QuickCheck

mayoresDivisores :: Int -> [Int]

mayoresDivisores x =

[max u v | u <- [2..floor (sqrt (fromIntegral x))]

, x `mod` u == 0

, let v = x `div` u]

arbol :: Int -> Tree Int

arbol 1 = Node 1 []

arbol x = Node x (arbol (x-1) : [arbol y | y <- mayoresDivisores x])

caminos :: Int -> [[Int]]

caminos = caminosArbol . arbol

-- λ> caminosArbol (arbol 6)

-- [[6,5,4,3,2,1],[6,5,4,2,1],[6,3,2,1]]

caminosArbol :: Tree a -> [[a]]

caminosArbol (Node x []) = [[x]]

caminosArbol (Node x as) = [x:ys | ys <- caminosBosque as]

caminosBosque :: Forest a -> [[a]]

caminosBosque = concatMap caminosArbol

caminosMinimales :: Int -> [[Int]]

caminosMinimales x = [ys | ys <- yss, length ys == m]

where yss = caminos x

m = minimum (map length yss)

Pensamiento

Tras el vivir y el soñar,
está lo que más importa:
despertar.

Antonio Machado

Número de descomposiciones en sumas de cuatro cuadrados

PorJosé A. Alonso 1-marzo-201926-marzo-2022

Definir la función

   nDescomposiciones       :: Int -> Int
   graficaDescomposiciones :: Int -> IO ()

1 2	nDescomposiciones :: Int -> Int graficaDescomposiciones :: Int -> IO ()

tales que

(nDescomposiciones x) es el número de listas de los cuadrados de cuatro números enteros positivos cuya suma es x. Por ejemplo.

     nDescomposiciones 4      ==  1
     nDescomposiciones 5      ==  0
     nDescomposiciones 7      ==  4
     nDescomposiciones 10     ==  6
     nDescomposiciones 15     ==  12
     nDescomposiciones 50000  ==  5682

nDescomposiciones 4 == 1

nDescomposiciones 5 == 0

nDescomposiciones 7 == 4

nDescomposiciones 10 == 6

nDescomposiciones 15 == 12

nDescomposiciones 50000 == 5682

(graficaDescomposiciones n) dibuja la gráfica del número de descomposiciones de los n primeros números naturales. Por ejemplo, (graficaDescomposiciones 500) dibuja

Soluciones

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

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

nDescomposiciones :: Int -> Int
nDescomposiciones = length . descomposiciones

-- (descomposiciones x) es la lista de las listas de los cuadrados de
-- cuatro números enteros positivos cuya suma es x. Por  ejemplo. 
--    λ> descomposiciones 4
--    [[1,1,1,1]]
--    λ> descomposiciones 5
--    []
--    λ> descomposiciones 7
--    [[1,1,1,4],[1,1,4,1],[1,4,1,1],[4,1,1,1]]
--    λ> descomposiciones 10
--    [[1,1,4,4],[1,4,1,4],[1,4,4,1],[4,1,1,4],[4,1,4,1],[4,4,1,1]]
--    λ> descomposiciones 15
--    [[1,1,4,9],[1,1,9,4],[1,4,1,9],[1,4,9,1],[1,9,1,4],[1,9,4,1],
--     [4,1,1,9],[4,1,9,1],[4,9,1,1],[9,1,1,4],[9,1,4,1],[9,4,1,1]]
descomposiciones :: Int -> [[Int]]
descomposiciones x = aux x 4
  where 
    aux 0 1 = []
    aux 1 1 = [[1]]
    aux 2 1 = []
    aux 3 1 = []
    aux y 1 | esCuadrado y = [[y]]
            | otherwise    = []
    aux y n = [x^2 : zs | x <- [1..raizEntera y]
                        , zs <- aux (y - x^2) (n-1)]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Int -> Bool
esCuadrado x = (raizEntera x)^2 == x

-- (raizEntera n) es el mayor entero cuya raíz cuadrada es menor o igual
-- que n. Por ejemplo,
--    raizEntera 15  ==  3
--    raizEntera 16  ==  4
--    raizEntera 17  ==  4
raizEntera :: Int -> Int
raizEntera = floor . sqrt . fromIntegral 

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

nDescomposiciones2 :: Int -> Int
nDescomposiciones2 = length . descomposiciones2

descomposiciones2 :: Int -> [[Int]]
descomposiciones2 x = a ! (x,4)
  where
    a = array ((0,1),(x,4)) [((i,j), f i j) | i <- [0..x], j <- [1..4]]
    f 0 1 = []
    f 1 1 = [[1]]
    f 2 1 = []
    f 3 1 = []
    f i 1 | esCuadrado i = [[i]]
          | otherwise    = []
    f i j = [x^2 : zs | x <- [1..raizEntera i]
                      , zs <- a ! (i - x^2,j-1)]

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

nDescomposiciones3 :: Int -> Int
nDescomposiciones3 x = aux x 4
  where
    aux 0 1 = 0
    aux 1 1 = 1
    aux 2 1 = 0
    aux 3 1 = 0
    aux y 1 | esCuadrado y = 1
            | otherwise    = 0
    aux y n = sum [aux (y - x^2) (n-1) | x <- [1..raizEntera y]]

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

nDescomposiciones4 :: Int -> Int
nDescomposiciones4 x = a ! (x,4)
  where
    a = array ((0,1),(x,4)) [((i,j), f i j) | i <- [0..x], j <- [1..4]]
    f 0 1 = 0
    f 1 1 = 1
    f 2 1 = 0
    f 3 1 = 0
    f i 1 | esCuadrado i = 1
          | otherwise    = 0
    f i j = sum [a ! (i- x^2,j-1) | x <- [1..raizEntera i]]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_nDescomposiciones :: Positive Int -> Bool
prop_nDescomposiciones (Positive x) =
  all (== nDescomposiciones x) [f x | f <- [ nDescomposiciones2
                                           , nDescomposiciones3
                                           , nDescomposiciones4]]

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

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

--    λ> nDescomposiciones 20000
--    1068
--    (3.69 secs, 3,307,250,128 bytes)
--    λ> nDescomposiciones2 20000
--    1068
--    (0.72 secs, 678,419,328 bytes)
--    λ> nDescomposiciones3 20000
--    1068
--    (3.94 secs, 3,485,725,552 bytes)
--    λ> nDescomposiciones4 20000
--    1068
--    (0.74 secs, 716,022,456 bytes)
--    
--    λ> nDescomposiciones2 50000
--    5682
--    (2.64 secs, 2,444,206,000 bytes)
--    λ> nDescomposiciones4 50000
--    5682
--    (2.77 secs, 2,582,443,448 bytes)

-- Definición de graficaDescomposiciones
-- =====================================

graficaDescomposiciones :: Int -> IO ()
graficaDescomposiciones n =
  plotList [ Key Nothing
           , PNG ("Numero_de_descomposiciones_en_sumas_de_cuadrados.png")
           ]
           (map nDescomposiciones3 [0..n])

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

import Data.Array

import Graphics.Gnuplot.Simple

import Test.QuickCheck

-- 1ª solución

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

nDescomposiciones :: Int -> Int

nDescomposiciones = length . descomposiciones

-- (descomposiciones x) es la lista de las listas de los cuadrados de

-- cuatro números enteros positivos cuya suma es x. Por ejemplo.

-- λ> descomposiciones 4

-- [[1,1,1,1]]

-- λ> descomposiciones 5

-- []

-- λ> descomposiciones 7

-- [[1,1,1,4],[1,1,4,1],[1,4,1,1],[4,1,1,1]]

-- λ> descomposiciones 10

-- [[1,1,4,4],[1,4,1,4],[1,4,4,1],[4,1,1,4],[4,1,4,1],[4,4,1,1]]

-- λ> descomposiciones 15

-- [[1,1,4,9],[1,1,9,4],[1,4,1,9],[1,4,9,1],[1,9,1,4],[1,9,4,1],

-- [4,1,1,9],[4,1,9,1],[4,9,1,1],[9,1,1,4],[9,1,4,1],[9,4,1,1]]

descomposiciones :: Int -> [[Int]]

descomposiciones x = aux x 4

where

aux 0 1 = []

aux 1 1 = [[1]]

aux 2 1 = []

aux 3 1 = []

aux y 1 | esCuadrado y = [[y]]

| otherwise = []

aux y n = [x^2 : zs | x <- [1..raizEntera y]

, zs <- aux (y - x^2) (n-1)]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Int -> Bool

esCuadrado x = (raizEntera x)^2 == x

-- (raizEntera n) es el mayor entero cuya raíz cuadrada es menor o igual

-- que n. Por ejemplo,

-- raizEntera 15 == 3

-- raizEntera 16 == 4

-- raizEntera 17 == 4

raizEntera :: Int -> Int

raizEntera = floor . sqrt . fromIntegral

-- 2ª solución

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

nDescomposiciones2 :: Int -> Int

nDescomposiciones2 = length . descomposiciones2

descomposiciones2 :: Int -> [[Int]]

descomposiciones2 x = a ! (x,4)

where

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

f 0 1 = []

f 1 1 = [[1]]

f 2 1 = []

f 3 1 = []

f i 1 | esCuadrado i = [[i]]

| otherwise = []

f i j = [x^2 : zs | x <- [1..raizEntera i]

, zs <- a ! (i - x^2,j-1)]

-- 3ª solución

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

nDescomposiciones3 :: Int -> Int

nDescomposiciones3 x = aux x 4

where

aux 0 1 = 0

aux 1 1 = 1

aux 2 1 = 0

aux 3 1 = 0

aux y 1 | esCuadrado y = 1

| otherwise = 0

aux y n = sum [aux (y - x^2) (n-1) | x <- [1..raizEntera y]]

-- 4ª solución

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

nDescomposiciones4 :: Int -> Int

nDescomposiciones4 x = a ! (x,4)

where

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

f 0 1 = 0

f 1 1 = 1

f 2 1 = 0

f 3 1 = 0

f i 1 | esCuadrado i = 1

| otherwise = 0

f i j = sum [a ! (i- x^2,j-1) | x <- [1..raizEntera i]]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_nDescomposiciones :: Positive Int -> Bool

prop_nDescomposiciones (Positive x) =

all (== nDescomposiciones x) [f x | f <- [ nDescomposiciones2

, nDescomposiciones3

, nDescomposiciones4]]

-- La comprobación es

-- λ> quickCheck prop_nDescomposiciones

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- λ> nDescomposiciones 20000

-- 1068

-- (3.69 secs, 3,307,250,128 bytes)

-- λ> nDescomposiciones2 20000

-- 1068

-- (0.72 secs, 678,419,328 bytes)

-- λ> nDescomposiciones3 20000

-- 1068

-- (3.94 secs, 3,485,725,552 bytes)

-- λ> nDescomposiciones4 20000

-- 1068

-- (0.74 secs, 716,022,456 bytes)

-- λ> nDescomposiciones2 50000

-- 5682

-- (2.64 secs, 2,444,206,000 bytes)

-- λ> nDescomposiciones4 50000

-- 5682

-- (2.77 secs, 2,582,443,448 bytes)

-- Definición de graficaDescomposiciones

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

graficaDescomposiciones :: Int -> IO ()

graficaDescomposiciones n =

plotList [ Key Nothing

, PNG ("Numero_de_descomposiciones_en_sumas_de_cuadrados.png")

]

(map nDescomposiciones3 [0..n])

Pensamiento

Ya habrá cigüeñas al sol,
mirando la tarde roja,
entre Moncayo y Urbión.

Antonio Machado

Descomposiciones en sumas de cuatro cuadrados

PorJosé A. Alonso 28-febrero-20197-marzo-2019

Definir la función

   descomposiciones :: Int -> [[Int]]

1	descomposiciones :: Int -> [[Int]]

tal que (descomposiciones x) es la lista de las listas de los cuadrados de cuatro números enteros positivos cuya suma es x. Por ejemplo.

   λ> descomposiciones 4
   [[1,1,1,1]]
   λ> descomposiciones 5
   []
   λ> descomposiciones 7
   [[1,1,1,4],[1,1,4,1],[1,4,1,1],[4,1,1,1]]
   λ> descomposiciones 10
   [[1,1,4,4],[1,4,1,4],[1,4,4,1],[4,1,1,4],[4,1,4,1],[4,4,1,1]]
   λ> descomposiciones 15
   [[1,1,4,9],[1,1,9,4],[1,4,1,9],[1,4,9,1],[1,9,1,4],[1,9,4,1],
    [4,1,1,9],[4,1,9,1],[4,9,1,1],[9,1,1,4],[9,1,4,1],[9,4,1,1]]
   λ> length (descomposiciones 50000)
   5682

λ> descomposiciones 4

[[1,1,1,1]]

λ> descomposiciones 5

[]

λ> descomposiciones 7

[[1,1,1,4],[1,1,4,1],[1,4,1,1],[4,1,1,1]]

λ> descomposiciones 10

[[1,1,4,4],[1,4,1,4],[1,4,4,1],[4,1,1,4],[4,1,4,1],[4,4,1,1]]

λ> descomposiciones 15

[[1,1,4,9],[1,1,9,4],[1,4,1,9],[1,4,9,1],[1,9,1,4],[1,9,4,1],

[4,1,1,9],[4,1,9,1],[4,9,1,1],[9,1,1,4],[9,1,4,1],[9,4,1,1]]

λ> length (descomposiciones 50000)

5682

Soluciones

import Data.Array
import Test.QuickCheck

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

descomposiciones :: Int -> [[Int]]
descomposiciones x = aux x 4
  where 
    aux 0 1 = []
    aux 1 1 = [[1]]
    aux 2 1 = []
    aux 3 1 = []
    aux y 1 | esCuadrado y = [[y]]
            | otherwise    = []
    aux y n = [x^2 : zs | x <- [1..raizEntera y]
                        , zs <- aux (y - x^2) (n-1)]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Int -> Bool
esCuadrado x = (raizEntera x)^2 == x

-- (raizEntera n) es el mayor entero cuya raíz cuadrada es menor o igual
-- que n. Por ejemplo,
--    raizEntera 15  ==  3
--    raizEntera 16  ==  4
--    raizEntera 17  ==  4
raizEntera :: Int -> Int
raizEntera = floor . sqrt . fromIntegral 

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

descomposiciones2 :: Int -> [[Int]]
descomposiciones2 x = a ! (x,4)
  where
    a = array ((0,1),(x,4)) [((i,j), f i j) | i <- [0..x], j <- [1..4]]
    f 0 1 = []
    f 1 1 = [[1]]
    f 2 1 = []
    f 3 1 = []
    f i 1 | esCuadrado i = [[i]]
          | otherwise    = []
    f i j = [x^2 : zs | x <- [1..raizEntera i]
                      , zs <- a ! (i - x^2,j-1)]

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_descomposiciones :: Positive Int -> Bool
prop_descomposiciones (Positive x) =
  descomposiciones x == descomposiciones2 x

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

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

-- La comparación es
--    λ> length (descomposiciones (2*10^4))
--    1068
--    (3.70 secs, 3,307,251,704 bytes)
--    λ> length (descomposiciones2 (2*10^4))
--    1068
--    (0.72 secs, 678,416,144 bytes)

import Data.Array

import Test.QuickCheck

-- 1ª definición

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

descomposiciones :: Int -> [[Int]]

descomposiciones x = aux x 4

where

aux 0 1 = []

aux 1 1 = [[1]]

aux 2 1 = []

aux 3 1 = []

aux y 1 | esCuadrado y = [[y]]

| otherwise = []

aux y n = [x^2 : zs | x <- [1..raizEntera y]

, zs <- aux (y - x^2) (n-1)]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Int -> Bool

esCuadrado x = (raizEntera x)^2 == x

-- (raizEntera n) es el mayor entero cuya raíz cuadrada es menor o igual

-- que n. Por ejemplo,

-- raizEntera 15 == 3

-- raizEntera 16 == 4

-- raizEntera 17 == 4

raizEntera :: Int -> Int

raizEntera = floor . sqrt . fromIntegral

-- 2ª definición

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

descomposiciones2 :: Int -> [[Int]]

descomposiciones2 x = a ! (x,4)

where

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

f 0 1 = []

f 1 1 = [[1]]

f 2 1 = []

f 3 1 = []

f i 1 | esCuadrado i = [[i]]

| otherwise = []

f i j = [x^2 : zs | x <- [1..raizEntera i]

, zs <- a ! (i - x^2,j-1)]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_descomposiciones :: Positive Int -> Bool

prop_descomposiciones (Positive x) =

descomposiciones x == descomposiciones2 x

-- La comprobación es

-- λ> quickCheck prop_descomposiciones

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (descomposiciones (2*10^4))

-- 1068

-- (3.70 secs, 3,307,251,704 bytes)

-- λ> length (descomposiciones2 (2*10^4))

-- 1068

-- (0.72 secs, 678,416,144 bytes)

Pensamiento

No extrañéis, dulces amigos,
que esté mi frente arrugada;
yo vivo en paz con los hombres
y en guerra con mis entrañas.

Antonio Machado

Número de sumandos en suma de cuadrados

PorJosé A. Alonso 21-febrero-201926-marzo-2022

El teorema de Lagrange de los cuatro cuadrados asegura que cualquier número entero positivo es la suma de, como máximo,cuatro cuadrados de números enteros. Por ejemplo,

   16 = 4²
   29 = 2² + 5²  
   14 = 1² + 2² + 3²
   15 = 1² + 1² + 2² + 3²

16 = 4²

29 = 2² + 5²

14 = 1² + 2² + 3²

15 = 1² + 1² + 2² + 3²

Definir las funciones

   ordenLagrange        :: Integer -> Int
   graficaOrdenLagrange :: Integer -> IO ()

1 2	ordenLagrange :: Integer -> Int graficaOrdenLagrange :: Integer -> IO ()

tales que

(ordenLagrange n) es el menor número de cuadrados necesarios para escribir n como suma de cuadrados. Por ejemplo.

     ordenLagrange 16     ==  1
     ordenLagrange 29     ==  2
     ordenLagrange 14     ==  3
     ordenLagrange 15     ==  4
     ordenLagrange 10000  ==  1
     ordenLagrange 10001  ==  2
     ordenLagrange 10002  ==  3
     ordenLagrange 10007  ==  4

ordenLagrange 16 == 1

ordenLagrange 29 == 2

ordenLagrange 14 == 3

ordenLagrange 15 == 4

ordenLagrange 10000 == 1

ordenLagrange 10001 == 2

ordenLagrange 10002 == 3

ordenLagrange 10007 == 4

(graficaOrdenLagrange n) dibuja la gráfica de los órdenes de Lagrange de los n primeros números naturales. Por ejemplo, (graficaOrdenLagrange 100) dibuja

Comprobar con QuickCheck que. para todo entero positivo k, el orden de Lagrange de k es menos o igual que 4, el de 4k+3 es distinto de 2 y el de 8k+7 es distinto de 3.

Soluciones

import Data.Array (Array, (!), array)
import Graphics.Gnuplot.Simple

import Test.QuickCheck

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

ordenLagrange :: Integer -> Int
ordenLagrange n
  | esCuadrado n = 1
  | otherwise    = 1 + minimum [ ordenLagrange (n - x^2)
                               | x <- [1..raizEntera n]]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = (raizEntera x)^2 == x

-- (raizEntera n) es el mayor entero cuya raíz cuadrada es menor o igual
-- que n. Por ejemplo,
--    raizEntera 15  ==  3
--    raizEntera 16  ==  4
--    raizEntera 17  ==  4
raizEntera :: Integer -> Integer
raizEntera = floor . sqrt . fromIntegral 

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

ordenLagrange2 :: Integer -> Int
ordenLagrange2 n = (vectorOrdenLagrange n) ! n

vectorOrdenLagrange :: Integer -> Array Integer Int
vectorOrdenLagrange n = v where
  v = array (0,n) [(i,f i) | i <- [0..n]]
  f i | esCuadrado i = 1
      | otherwise    = 1 + minimum [ v ! (i - j^2)
                                   | j <- [1..raizEntera i]]

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

-- La comparación es
--    λ> ordenLagrange 50
--    2
--    (10.39 secs, 1,704,144,464 bytes)
--    λ> ordenLagrange2 50
--    2
--    (0.01 secs, 341,920 bytes)

-- Definición de graficaOrdenLagrange
-- ==================================

graficaOrdenLagrange :: Integer -> IO ()
graficaOrdenLagrange n = 
  plotList [ Key Nothing
           , PNG ("Numero_de_sumandos_en_suma_de_cuadrados.png")
           ]
           (map ordenLagrange2 [0..n-1])

-- Comprobación de la propiedad
-- ============================

-- La propiedad es
prop_OrdenLagrange :: Positive Integer -> Bool
prop_OrdenLagrange (Positive k) =
  ordenLagrange2 k <= 4 &&
  ordenLagrange2 (4*k+3) /= 2 &&
  ordenLagrange2 (8*k+7) /= 3

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

import Data.Array (Array, (!), array)

import Graphics.Gnuplot.Simple

import Test.QuickCheck

-- 1ª definición

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

ordenLagrange :: Integer -> Int

ordenLagrange n

| esCuadrado n = 1

| otherwise = 1 + minimum [ ordenLagrange (n - x^2)

| x <- [1..raizEntera n]]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = (raizEntera x)^2 == x

-- (raizEntera n) es el mayor entero cuya raíz cuadrada es menor o igual

-- que n. Por ejemplo,

-- raizEntera 15 == 3

-- raizEntera 16 == 4

-- raizEntera 17 == 4

raizEntera :: Integer -> Integer

raizEntera = floor . sqrt . fromIntegral

-- 2ª definición

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

ordenLagrange2 :: Integer -> Int

ordenLagrange2 n = (vectorOrdenLagrange n) ! n

vectorOrdenLagrange :: Integer -> Array Integer Int

vectorOrdenLagrange n = v where

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

f i | esCuadrado i = 1

| otherwise = 1 + minimum [ v ! (i - j^2)

| j <- [1..raizEntera i]]

-- Comparación de eficiencia

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

-- La comparación es

-- λ> ordenLagrange 50

-- 2

-- (10.39 secs, 1,704,144,464 bytes)

-- λ> ordenLagrange2 50

-- 2

-- (0.01 secs, 341,920 bytes)

-- Definición de graficaOrdenLagrange

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

graficaOrdenLagrange :: Integer -> IO ()

graficaOrdenLagrange n =

plotList [ Key Nothing

, PNG ("Numero_de_sumandos_en_suma_de_cuadrados.png")

]

(map ordenLagrange2 [0..n-1])

-- Comprobación de la propiedad

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

-- La propiedad es

prop_OrdenLagrange :: Positive Integer -> Bool

prop_OrdenLagrange (Positive k) =

ordenLagrange2 k <= 4 &&

ordenLagrange2 (4*k+3) /= 2 &&

ordenLagrange2 (8*k+7) /= 3

-- La comprobación es

-- λ> quickCheck prop_OrdenLagrange

-- +++ OK, passed 100 tests.

Pensamiento

— Nuestro español bosteza.
¿Es hambre? ¿Sueño? ¿Hastío?
Doctor, ¿tendrá el estómago vacío?
— El vacío es más bien en la cabeza.

Antonio Machado

Ternas euclídeas

PorJosé A. Alonso 18-febrero-201925-febrero-2019

Uno de los problemas planteados por Euclides en los Elementos consiste en encontrar tres números tales que cada uno de sus productos, dos a dos, aumentados en la unidad sea un cuadrado perfecto.

Diremos que (x,y,z) es una terna euclídea si es una solución del problema; es decir, si x <= y <= z y xy+1, yz+1 y zx+1 son cuadrados. Por ejemplo, (4,6,20) es una terna euclídea ya que

   4x6+1 = 5^2, 6x20+1 = 11^2 y 20*4+1 = 9^2

1	4x6+1 = 5^2, 6x20+1 = 11^2 y 20*4+1 = 9^2

Definir la funciones

   ternasEuclideas        :: [(Integer,Integer,Integer)]
   esMayorDeTernaEuclidea :: Integer -> Bool

1 2	ternasEuclideas :: [(Integer,Integer,Integer)] esMayorDeTernaEuclidea :: Integer -> Bool

tales que

ternasEuclideas es la lista de las ternas euclídeas. Por ejemplo,

     λ> take 7 ternasEuclideas
     [(1,3,8),(2,4,12),(1,8,15),(3,5,16),(4,6,20),(3,8,21),(5,7,24)]

1 2	λ> take 7 ternasEuclideas [(1,3,8),(2,4,12),(1,8,15),(3,5,16),(4,6,20),(3,8,21),(5,7,24)]

(esMayorDeTernaEuclidea z) se verifica si existen x, y tales que (x,y,z) es una terna euclídea. Por ejemplo,

     esMayorDeTernaEuclidea 20  ==  True
     esMayorDeTernaEuclidea 22  ==  False

1 2	esMayorDeTernaEuclidea 20 == True esMayorDeTernaEuclidea 22 == False

Comprobar con QuickCheck que z es el mayor de una terna euclídea si, y sólo si, existe un número natural x tal que 1 < x < z – 1 y x^2 es congruente con 1 módulo z.

Soluciones

import Test.QuickCheck

ternasEuclideas :: [(Integer,Integer,Integer)]
ternasEuclideas =
  [(x,y,z) | z <- [1..]
           , y <- [1..z]
           , esCuadrado (y * z + 1)
           , x <- [1..y]
           , esCuadrado (x * y + 1)
           , esCuadrado (z * x + 1)]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = (raizEntera x)^2 == x
  where raizEntera :: Integer -> Integer
        raizEntera = floor . sqrt . fromIntegral 

esMayorDeTernaEuclidea :: Integer -> Bool
esMayorDeTernaEuclidea z =
  not (null [(x,y) | y <- [1..z]
                   , esCuadrado (y * z + 1)
                   , x <- [1..y]
                   , esCuadrado (x * y + 1)
                   , esCuadrado (z * x + 1)])


-- La propiedad es
prop_esMayorDeTernaEuclidea :: Positive Integer -> Bool
prop_esMayorDeTernaEuclidea (Positive z) =
  esMayorDeTernaEuclidea z == any (\x -> (x^2) `mod` z == 1) [2..z-2]

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

import Test.QuickCheck

ternasEuclideas :: [(Integer,Integer,Integer)]

ternasEuclideas =

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

, y <- [1..z]

, esCuadrado (y * z + 1)

, x <- [1..y]

, esCuadrado (x * y + 1)

, esCuadrado (z * x + 1)]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = (raizEntera x)^2 == x

where raizEntera :: Integer -> Integer

raizEntera = floor . sqrt . fromIntegral

esMayorDeTernaEuclidea :: Integer -> Bool

esMayorDeTernaEuclidea z =

not (null [(x,y) | y <- [1..z]

, esCuadrado (y * z + 1)

, x <- [1..y]

, esCuadrado (x * y + 1)

, esCuadrado (z * x + 1)])

-- La propiedad es

prop_esMayorDeTernaEuclidea :: Positive Integer -> Bool

prop_esMayorDeTernaEuclidea (Positive z) =

esMayorDeTernaEuclidea z == any (\x -> (x^2) `mod` z == 1) [2..z-2]

-- La comprobación es

-- λ> quickCheck prop_esMayorDeTernaEuclidea

-- +++ OK, passed 100 tests.

Pensamiento

Todo pasa y todo queda,
pero lo nuestro es pasar,
pasar haciendo caminos,
caminos sobre la mar.

Antonio Machado

Medias de dígitos de pi

PorJosé A. Alonso 5-febrero-201925-abril-2021

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)

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

Dígitos en las posiciones pares de cuadrados

PorJosé A. Alonso 1-febrero-20198-febrero-2019

Definir las funciones

   digitosPosParesCuadrado    :: Integer -> ([Integer],Int)
   invDigitosPosParesCuadrado :: ([Integer],Int) -> [Integer]

1 2	digitosPosParesCuadrado :: Integer -> ([Integer],Int) invDigitosPosParesCuadrado :: ([Integer],Int) -> [Integer]

tales que

(digitosPosParesCuadrado n) es el par formados por los dígitos de n² en la posiciones pares y por el número de dígitos de n². Por ejemplo,

     digitosPosParesCuadrado 8     ==  ([6],2)
     digitosPosParesCuadrado 14    ==  ([1,6],3)
     digitosPosParesCuadrado 36    ==  ([1,9],4)
     digitosPosParesCuadrado 116   ==  ([1,4,6],5)
     digitosPosParesCuadrado 2019  ==  ([4,7,3,1],7)

digitosPosParesCuadrado 8 == ([6],2)

digitosPosParesCuadrado 14 == ([1,6],3)

digitosPosParesCuadrado 36 == ([1,9],4)

digitosPosParesCuadrado 116 == ([1,4,6],5)

digitosPosParesCuadrado 2019 == ([4,7,3,1],7)

(invDigitosPosParesCuadrado (xs,k)) es la lista de los números n tales que xs es la lista de los dígitos de n² en la posiciones pares y k es el número de dígitos de n². Por ejemplo,

     invDigitosPosParesCuadrado ([6],2)             ==  [8]
     invDigitosPosParesCuadrado ([1,6],3)           ==  [14]
     invDigitosPosParesCuadrado ([1,9],4)           ==  [36]
     invDigitosPosParesCuadrado ([1,4,6],5)         ==  [116,136]
     invDigitosPosParesCuadrado ([4,7,3,1],7)       ==  [2019,2139,2231]
     invDigitosPosParesCuadrado ([1,2],3)           ==  []
     invDigitosPosParesCuadrado ([1,2],4)           ==  [32,35,39]
     invDigitosPosParesCuadrado ([1,2,3,4,5,6],11)  ==  [115256,127334,135254]

invDigitosPosParesCuadrado ([6],2) == [8]

invDigitosPosParesCuadrado ([1,6],3) == [14]

invDigitosPosParesCuadrado ([1,9],4) == [36]

invDigitosPosParesCuadrado ([1,4,6],5) == [116,136]

invDigitosPosParesCuadrado ([4,7,3,1],7) == [2019,2139,2231]

invDigitosPosParesCuadrado ([1,2],3) == []

invDigitosPosParesCuadrado ([1,2],4) == [32,35,39]

invDigitosPosParesCuadrado ([1,2,3,4,5,6],11) == [115256,127334,135254]

Comprobar con QuickCheck que para todo entero positivo n se verifica que para todo entero positivo m, m pertenece a (invDigitosPosParesCuadrado (digitosPosParesCuadrado n)) si, y sólo si, (digitosPosParesCuadrado m) es igual a (digitosPosParesCuadrado n)

Soluciones

import Test.QuickCheck

-- Definición de digitosPosParesCuadrado
-- =====================================

digitosPosParesCuadrado :: Integer -> ([Integer],Int)
digitosPosParesCuadrado n =
  (digitosPosPares (n^2),length (show (n^2)))

-- (digitosPosPares n) es la lista de los dígitos de n en posiciones
-- pares. Por ejemplo,
--    digitosPosPares 24012019  ==  [2,0,2,1]
digitosPosPares :: Integer -> [Integer]
digitosPosPares n = elementosPosPares (digitos n)

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

-- (elementosPosPares xs) es la lista de los elementos de xs en
-- posiciones pares. Por ejemplo,
--    elementosPosPares [3,2,5,7,6,4]  ==  [3,5,6]
elementosPosPares :: [a] -> [a]
elementosPosPares []       = []
elementosPosPares [x]      = [x]
elementosPosPares (x:_:zs) = x : elementosPosPares zs

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

invDigitosPosParesCuadrado :: ([Integer],Int) -> [Integer]
invDigitosPosParesCuadrado (xs, a) =
  [x | x <- [ceiling (sqrt 10^(a-1))..ceiling (sqrt 10^a)]
     , digitosPosParesCuadrado x == (xs,a)]

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

invDigitosPosParesCuadrado2 :: ([Integer],Int) -> [Integer]
invDigitosPosParesCuadrado2 x =
  [n | n <- [a..b], digitosPosParesCuadrado n == x]
  where a = floor (sqrt (fromIntegral (completaNum x 0)))
        b = ceiling (sqrt (fromIntegral (completaNum x 9)))

-- (completaNum (xs,k) n) es el número cuyos dígitos en las posiciones
-- pares son los de xs y los de las posiciones impares son iguales a n
-- (se supone que k es igual al doble de la longitud de xs o un
-- menos). Por ejemplo, 
--    completaNum ([1,3,8],5) 4  ==  14348
--    completaNum ([1,3,8],6) 4  ==  143484
completaNum :: ([Integer],Int) -> Integer -> Integer
completaNum x n = digitosAnumero (completa x n)

-- (completa (xs,k) n) es la lista cuyos elementos en las posiciones
-- pares son los de xs y los de las posiciones impares son iguales a n
-- (se supone que k es igual al doble de la longitud de xs o un
-- menos). Por ejemplo, 
--    completa ([1,3,8],5) 4  ==  [1,4,3,4,8]
--    completa ([1,3,8],6) 4  ==  [1,4,3,4,8,4]
completa :: ([Integer],Int) -> Integer -> [Integer]
completa (xs,k) n
  | even k    = ys
  | otherwise = init ys
  where ys = concat [[x,n] | x <- xs]

-- (digitosAnumero ds) es el número cuyos dígitos son ds. Por ejemplo,
--    digitosAnumero [2,0,1,9]  ==  2019
digitosAnumero :: [Integer] -> Integer
digitosAnumero = read . concatMap show

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

--    λ> invDigitosPosParesCuadrado ([1,2,1,5,7,4,9],13)
--    [1106393,1234567,1314597]
--    (7.55 secs, 13,764,850,536 bytes)
--    λ> invDigitosPosParesCuadrado2 ([1,2,1,5,7,4,9],13)
--    [1106393,1234567,1314597]
--    (1.96 secs, 3,780,368,816 bytes)

-- Comprobación de la propiedad
-- ============================

-- La propiedad es  
prop_digitosPosParesCuadrado :: Positive Integer -> Positive Integer -> Bool
prop_digitosPosParesCuadrado (Positive n) (Positive m) =
  (digitosPosParesCuadrado m == x)
  == (m `elem` invDigitosPosParesCuadrado x)
  where x = digitosPosParesCuadrado n

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

import Test.QuickCheck

-- Definición de digitosPosParesCuadrado

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

digitosPosParesCuadrado :: Integer -> ([Integer],Int)

digitosPosParesCuadrado n =

(digitosPosPares (n^2),length (show (n^2)))

-- (digitosPosPares n) es la lista de los dígitos de n en posiciones

-- pares. Por ejemplo,

-- digitosPosPares 24012019 == [2,0,2,1]

digitosPosPares :: Integer -> [Integer]

digitosPosPares n = elementosPosPares (digitos n)

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

-- digitos 325 == [3,2,5]

digitos :: Integer -> [Integer]

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

-- (elementosPosPares xs) es la lista de los elementos de xs en

-- posiciones pares. Por ejemplo,

-- elementosPosPares [3,2,5,7,6,4] == [3,5,6]

elementosPosPares :: [a] -> [a]

elementosPosPares [] = []

elementosPosPares [x] = [x]

elementosPosPares (x:_:zs) = x : elementosPosPares zs

-- 1ª definición de invDigitosPosParesCuadrado

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

invDigitosPosParesCuadrado :: ([Integer],Int) -> [Integer]

invDigitosPosParesCuadrado (xs, a) =

[x | x <- [ceiling (sqrt 10^(a-1))..ceiling (sqrt 10^a)]

, digitosPosParesCuadrado x == (xs,a)]

-- 2ª definición de invDigitosPosParesCuadrado

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

invDigitosPosParesCuadrado2 :: ([Integer],Int) -> [Integer]

invDigitosPosParesCuadrado2 x =

[n | n <- [a..b], digitosPosParesCuadrado n == x]

where a = floor (sqrt (fromIntegral (completaNum x 0)))

b = ceiling (sqrt (fromIntegral (completaNum x 9)))

-- (completaNum (xs,k) n) es el número cuyos dígitos en las posiciones

-- pares son los de xs y los de las posiciones impares son iguales a n

-- (se supone que k es igual al doble de la longitud de xs o un

-- menos). Por ejemplo,

-- completaNum ([1,3,8],5) 4 == 14348

-- completaNum ([1,3,8],6) 4 == 143484

completaNum :: ([Integer],Int) -> Integer -> Integer

completaNum x n = digitosAnumero (completa x n)

-- (completa (xs,k) n) es la lista cuyos elementos en las posiciones

-- pares son los de xs y los de las posiciones impares son iguales a n

-- (se supone que k es igual al doble de la longitud de xs o un

-- menos). Por ejemplo,

-- completa ([1,3,8],5) 4 == [1,4,3,4,8]

-- completa ([1,3,8],6) 4 == [1,4,3,4,8,4]

completa :: ([Integer],Int) -> Integer -> [Integer]

completa (xs,k) n

| even k = ys

| otherwise = init ys

where ys = concat [[x,n] | x <- xs]

-- (digitosAnumero ds) es el número cuyos dígitos son ds. Por ejemplo,

-- digitosAnumero [2,0,1,9] == 2019

digitosAnumero :: [Integer] -> Integer

digitosAnumero = read . concatMap show

-- Comparación de eficiencia

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

-- λ> invDigitosPosParesCuadrado ([1,2,1,5,7,4,9],13)

-- [1106393,1234567,1314597]

-- (7.55 secs, 13,764,850,536 bytes)

-- λ> invDigitosPosParesCuadrado2 ([1,2,1,5,7,4,9],13)

-- [1106393,1234567,1314597]

-- (1.96 secs, 3,780,368,816 bytes)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_digitosPosParesCuadrado :: Positive Integer -> Positive Integer -> Bool

prop_digitosPosParesCuadrado (Positive n) (Positive m) =

(digitosPosParesCuadrado m == x)

== (m `elem` invDigitosPosParesCuadrado x)

where x = digitosPosParesCuadrado n

-- La comprobación es

-- λ> quickCheck prop_digitosPosParesCuadrado

-- +++ OK, passed 100 tests.

Pensamiento

¡Ojos que a la luz se abrieron
un día para, después,
ciegos tornar a la tierra,
hartos de mirar sin ver.

Antonio Machado

El teorema de Navidad de Fermat

PorJosé A. Alonso 31-diciembre-201825-abril-2021

Definir las funciones

   representaciones :: Integer -> [(Integer,Integer)]
   primosImparesConRepresentacionUnica :: [Integer]
   primos4nM1 :: [Integer]

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

primosImparesConRepresentacionUnica :: [Integer]

primos4nM1 :: [Integer]

tales que

(representaciones n) es la lista de pares de números naturales (x,y) tales que n = x^2 + y^2 con x <= y. Por ejemplo.

     representaciones  20           ==  [(2,4)]
     representaciones  25           ==  [(0,5),(3,4)]
     representaciones 325           ==  [(1,18),(6,17),(10,15)]
     representaciones 100000147984  ==  [(0,316228)]
     length (representaciones (10^10))    ==  6
     length (representaciones (4*10^12))  ==  7

representaciones 20 == [(2,4)]

representaciones 25 == [(0,5),(3,4)]

representaciones 325 == [(1,18),(6,17),(10,15)]

representaciones 100000147984 == [(0,316228)]

length (representaciones (10^10)) == 6

length (representaciones (4*10^12)) == 7

primosImparesConRepresentacionUnica es la lista de los números primos impares que se pueden escribir exactamente de una manera como suma de cuadrados de pares de números naturales (x,y) con x <= y. Por ejemplo,

     λ> take 20 primosImparesConRepresentacionUnica
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

1 2	λ> take 20 primosImparesConRepresentacionUnica [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

primos4nM1 es la lista de los números primos que se pueden escribir como uno más un múltiplo de 4 (es decir, que son congruentes con 1 módulo 4). Por ejemplo,

     λ> take 20 primos4nM1
     [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

1 2	λ> take 20 primos4nM1 [5,13,17,29,37,41,53,61,73,89,97,101,109,113,137,149,157,173,181,193]

Comprobar con QuickCheck el torema de Navidad de Fermat; es decir, que para todo número n, los n-ésimos elementos de primosImparesConRepresentacionUnica y de primos4nM1 son iguales.

Soluciones

import Data.Numbers.Primes (primes)
import Test.QuickCheck

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

representaciones :: Integer -> [(Integer,Integer)]
representaciones n =
  [(x,y) | x <- [0..n], y <- [x..n], n == x*x + y*y]

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

representaciones2 :: Integer -> [(Integer,Integer)]
representaciones2 n =
  [(x,raiz z) | x <- [0..raiz (n `div` 2)] 
              , let z = n - x*x
              , esCuadrado z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado 25  ==  True
--    esCuadrado 26  ==  False
esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz 25  ==  5
--    raiz 24  ==  4
--    raiz 26  ==  5
raiz :: Integer -> Integer 
raiz 0 = 0
raiz 1 = 1
raiz x = aux (0,x)
    where aux (a,b) | d == x    = c
                    | c == a    = a
                    | d < x     = aux (c,b)
                    | otherwise = aux (a,c) 
              where c = (a+b) `div` 2
                    d = c^2

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

representaciones3 :: Integer -> [(Integer,Integer)]
representaciones3 n =
  [(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)] 
               , let z = n - x*x
               , esCuadrado3 z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por
-- ejemplo,
--    esCuadrado3 25  ==  True
--    esCuadrado3 26  ==  False
esCuadrado3 :: Integer -> Bool
esCuadrado3 x = x == y * y
  where y = raiz3 x

-- (raiz3 x) es la raíz cuadrada entera de x. Por ejemplo,
--    raiz3 25  ==  5
--    raiz3 24  ==  4
--    raiz3 26  ==  5
raiz3 :: Integer -> Integer
raiz3 x = floor (sqrt (fromIntegral x))
          
-- 4ª definición de representaciones
-- =================================
 
representaciones4 :: Integer -> [(Integer, Integer)]
representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))
  where aux x y
          | x > y     = [] 
          | otherwise = case compare (x*x + y*y) n of
                          LT -> aux (x + 1) y
                          EQ -> (x, y) : aux (x + 1) (y - 1)
                          GT -> aux x (y - 1)

-- Equivalencia de las definiciones de representaciones
-- ====================================================

-- La propiedad es
prop_representaciones_equiv :: (Positive Integer) -> Bool
prop_representaciones_equiv (Positive n) =
  representaciones  n == representaciones2 n &&
  representaciones2 n == representaciones3 n &&
  representaciones3 n == representaciones4 n

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

-- Comparación de eficiencia de las definiciones de representaciones
-- =================================================================

--    λ> representaciones 3025
--    [(0,55),(33,44)]
--    (2.86 secs, 1,393,133,528 bytes)
--    λ> representaciones2 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 867,944 bytes)
--    λ> representaciones3 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 173,512 bytes)
--    λ> representaciones4 3025
--    [(0,55),(33,44)]
--    (0.00 secs, 423,424 bytes)
--    
--    λ> length (representaciones2 (10^10))
--    6
--    (3.38 secs, 2,188,903,544 bytes)
--    λ> length (representaciones3 (10^10))
--    6
--    (0.10 secs, 62,349,048 bytes)
--    λ> length (representaciones4 (10^10))
--    6
--    (0.11 secs, 48,052,360 bytes)
--
--    λ> length (representaciones3 (4*10^12))
--    7
--    (1.85 secs, 1,222,007,176 bytes)
--    λ> length (representaciones4 (4*10^12))
--    7
--    (1.79 secs, 953,497,480 bytes)

-- Definición de primosImparesConRepresentacionUnica
-- =================================================

primosImparesConRepresentacionUnica :: [Integer]
primosImparesConRepresentacionUnica =
  [x | x <- tail primes
     , length (representaciones4 x) == 1]

-- Definición de primos4nM1
-- ========================

primos4nM1 :: [Integer]
primos4nM1 = [x | x <- primes
                , x `mod` 4 == 1]

-- Teorema de Navidad de Fermat
-- ============================

-- La propiedad es
prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool
prop_teoremaDeNavidadDeFermat (Positive n) =
  primosImparesConRepresentacionUnica !! n == primos4nM1 !! n
             
-- La comprobación es
--    λ> quickCheck prop_teoremaDeNavidadDeFermat
--    +++ OK, passed 100 tests.

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

import Data.Numbers.Primes (primes)

import Test.QuickCheck

-- 1ª definición de representaciones

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

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

representaciones n =

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

-- 2ª definición de representaciones

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

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

representaciones2 n =

[(x,raiz z) | x <- [0..raiz (n `div` 2)]

, let z = n - x*x

, esCuadrado z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado 25 == True

-- esCuadrado 26 == False

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

-- (raiz x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz 25 == 5

-- raiz 24 == 4

-- raiz 26 == 5

raiz :: Integer -> Integer

raiz 0 = 0

raiz 1 = 1

raiz x = aux (0,x)

where aux (a,b) | d == x = c

| c == a = a

| d < x = aux (c,b)

| otherwise = aux (a,c)

where c = (a+b) `div` 2

d = c^2

-- 3ª definición de representaciones

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

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

representaciones3 n =

[(x,raiz3 z) | x <- [0..raiz3 (n `div` 2)]

, let z = n - x*x

, esCuadrado3 z]

-- (esCuadrado x) se verifica si x es un número al cuadrado. Por

-- ejemplo,

-- esCuadrado3 25 == True

-- esCuadrado3 26 == False

esCuadrado3 :: Integer -> Bool

esCuadrado3 x = x == y * y

where y = raiz3 x

-- (raiz3 x) es la raíz cuadrada entera de x. Por ejemplo,

-- raiz3 25 == 5

-- raiz3 24 == 4

-- raiz3 26 == 5

raiz3 :: Integer -> Integer

raiz3 x = floor (sqrt (fromIntegral x))

-- 4ª definición de representaciones

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

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

representaciones4 n = aux 0 (floor (sqrt (fromIntegral n)))

where aux x y

| x > y = []

| otherwise = case compare (x*x + y*y) n of

LT -> aux (x + 1) y

EQ -> (x, y) : aux (x + 1) (y - 1)

GT -> aux x (y - 1)

-- Equivalencia de las definiciones de representaciones

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

-- La propiedad es

prop_representaciones_equiv :: (Positive Integer) -> Bool

prop_representaciones_equiv (Positive n) =

representaciones n == representaciones2 n &&

representaciones2 n == representaciones3 n &&

representaciones3 n == representaciones4 n

-- La comprobación es

-- λ> quickCheck prop_representaciones_equiv

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de las definiciones de representaciones

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

-- λ> representaciones 3025

-- [(0,55),(33,44)]

-- (2.86 secs, 1,393,133,528 bytes)

-- λ> representaciones2 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 867,944 bytes)

-- λ> representaciones3 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 173,512 bytes)

-- λ> representaciones4 3025

-- [(0,55),(33,44)]

-- (0.00 secs, 423,424 bytes)

-- λ> length (representaciones2 (10^10))

-- 6

-- (3.38 secs, 2,188,903,544 bytes)

-- λ> length (representaciones3 (10^10))

-- 6

-- (0.10 secs, 62,349,048 bytes)

-- λ> length (representaciones4 (10^10))

-- 6

-- (0.11 secs, 48,052,360 bytes)

-- λ> length (representaciones3 (4*10^12))

-- 7

-- (1.85 secs, 1,222,007,176 bytes)

-- λ> length (representaciones4 (4*10^12))

-- 7

-- (1.79 secs, 953,497,480 bytes)

-- Definición de primosImparesConRepresentacionUnica

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

primosImparesConRepresentacionUnica :: [Integer]

primosImparesConRepresentacionUnica =

[x | x <- tail primes

, length (representaciones4 x) == 1]

-- Definición de primos4nM1

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

primos4nM1 :: [Integer]

primos4nM1 = [x | x <- primes

, x `mod` 4 == 1]

-- Teorema de Navidad de Fermat

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

-- La propiedad es

prop_teoremaDeNavidadDeFermat :: Positive Int -> Bool

prop_teoremaDeNavidadDeFermat (Positive n) =

primosImparesConRepresentacionUnica !! n == primos4nM1 !! n

-- La comprobación es

-- λ> quickCheck prop_teoremaDeNavidadDeFermat

-- +++ OK, passed 100 tests.

Pensamiento

– ¡Cuándo llegará otro día!
– Hoy es siempre todavía.

Antonio Machado

Aproximación entre pi y e

PorJosé A. Alonso 7-diciembre-201817-enero-2019

El día 11 de noviembre, se publicó en la cuenta de Twitter de Fermat’s Library la siguiente curiosa identidad que relaciona los números e y pi:

Definir las siguientes funciones:

   sumaTerminos :: Int -> Double
   aproximacion :: Double -> Int

1 2	sumaTerminos :: Int -> Double aproximacion :: Double -> Int

tales que

(sumaTerminos n) es la suma de los primeros n términos de la serie 1/(π²+ 1) + 1/(4π²+1) + 1/(9π²+1) + 1/(16π²+ ) + … Por ejemplo,

     sumaTerminos 10     ==  0.14687821811081034
     sumaTerminos 100    ==  0.15550948345688423
     sumaTerminos 1000   ==  0.15641637221314514
     sumaTerminos 10000  ==  0.15650751113789382

sumaTerminos 10 == 0.14687821811081034

sumaTerminos 100 == 0.15550948345688423

sumaTerminos 1000 == 0.15641637221314514

sumaTerminos 10000 == 0.15650751113789382

(aproximación x) es el menor número de términos que hay que sumar de la serie anterior para que se diferencie (en valor absoluto) de 1/(e²-1) menos que x. Por ejemplo,

     aproximacion 0.1     ==  1
     aproximacion 0.01    ==  10
     aproximacion 0.001   ==  101
     aproximacion 0.0001  ==  1013

aproximacion 0.1 == 1

aproximacion 0.01 == 10

aproximacion 0.001 == 101

aproximacion 0.0001 == 1013

Soluciones

-- 1ª definición de sumaTerminos
sumaTerminos :: Int -> Double
sumaTerminos n =
  sum [1 / (((x ^ 2) * (pi ^ 2)) + 1) | x <- [1 .. fromIntegral n]]

-- 2ª definición de sumaTerminos
sumaTerminos2 :: Int -> Double
sumaTerminos2 0 = 0
sumaTerminos2 n = 1 / (m^2 * pi^2 + 1) + sumaTerminos2 (n-1)
  where m = fromIntegral n

-- Definición de aproximacion
aproximacion :: Double -> Int
aproximacion x =
  head [n | n <- [0..]
          , abs (sumaTerminos n - 1 / (e^2 - 1)) < x]
  where e = exp 1

-- 1ª definición de sumaTerminos

sumaTerminos :: Int -> Double

sumaTerminos n =

sum [1 / (((x ^ 2) * (pi ^ 2)) + 1) | x <- [1 .. fromIntegral n]]

-- 2ª definición de sumaTerminos

sumaTerminos2 :: Int -> Double

sumaTerminos2 0 = 0

sumaTerminos2 n = 1 / (m^2 * pi^2 + 1) + sumaTerminos2 (n-1)

where m = fromIntegral n

-- Definición de aproximacion

aproximacion :: Double -> Int

aproximacion x =

head [n | n <- [0..]

, abs (sumaTerminos n - 1 / (e^2 - 1)) < x]

where e = exp 1

Pensamiento

«Sólo sé que no se nada» contenía la jactancia de un excesivo saber, puesto que olvidó añadir: y aun de esto mismo no estoy completamente seguro.

Antonio Machado

Soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Otras soluciones

Pensamiento

Soluciones

Pensamiento

Otras soluciones

Soluciones

Referencia

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Soluciones

Pensamiento

Entradas recientes

Buscar

Correo electrónico