José A. Alonso

Elementos finales

PorJosé A. Alonso 23-agosto-202214-diciembre-2022

Definir la función

   finales :: Int -> [a] -> [a]

1	finales :: Int -> [a] -> [a]

tal que (finales n xs) es la lista formada por los n finales elementos de xs. Por ejemplo,

   finales 3 [2,5,4,7,9,6]  ==  [7,9,6]

1	finales 3 [2,5,4,7,9,6] == [7,9,6]

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

import Test.QuickCheck

-- 1ª definición
finales1 :: Int -> [a] -> [a]
finales1 n xs = drop (length xs - n) xs

-- 2ª definición
finales2 :: Int -> [a] -> [a]
finales2 n xs = reverse (take n (reverse xs))

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

-- La propiedad es
prop_finales :: Int -> [Int] -> Bool
prop_finales n xs =
  finales1 n xs == finales2 n xs

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

import Test.QuickCheck

-- 1ª definición

finales1 :: Int -> [a] -> [a]

finales1 n xs = drop (length xs - n) xs

-- 2ª definición

finales2 :: Int -> [a] -> [a]

finales2 n xs = reverse (take n (reverse xs))

-- Comprobación de equivalencia

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

-- La propiedad es

prop_finales :: Int -> [Int] -> Bool

prop_finales n xs =

finales1 n xs == finales2 n xs

-- La comprobación es

-- λ> quickCheck prop_finales

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub

Soluciones en Python

from typing import TypeVar
from hypothesis import given, strategies as st

A = TypeVar('A')

# 1ª definición
def finales1(n: int, xs: list[A]) -> list[A]:
    if len(xs) <= n:
        return xs
    return xs[len(xs) - n:]

# 2ª definición
def finales2(n: int, xs: list[A]) -> list[A]:
    if n == 0:
        return []
    return xs[-n:]

# 3ª definición
def finales3(n: int, xs: list[A]) -> list[A]:
    ys = list(reversed(xs))
    return list(reversed(ys[:n]))

# La propiedad de equivalencia es
@given(st.integers(min_value=0), st.lists(st.integers()))
def test_equiv_finales(n, xs):
    assert finales1(n, xs) == finales2(n, xs) == finales3(n, xs)

# La comprobación es
#    src> poetry run pytest -q elementos_finales.py
#    1 passed in 0.18s

from typing import TypeVar

from hypothesis import given, strategies as st

A = TypeVar('A')

# 1ª definición

def finales1(n: int, xs: list[A]) -> list[A]:

if len(xs) <= n:

return xs

return xs[len(xs) - n:]

# 2ª definición

def finales2(n: int, xs: list[A]) -> list[A]:

if n == 0:

return []

return xs[-n:]

# 3ª definición

def finales3(n: int, xs: list[A]) -> list[A]:

ys = list(reversed(xs))

return list(reversed(ys[:n]))

# La propiedad de equivalencia es

@given(st.integers(min_value=0), st.lists(st.integers()))

def test_equiv_finales(n, xs):

assert finales1(n, xs) == finales2(n, xs) == finales3(n, xs)

# La comprobación es

# src> poetry run pytest -q elementos_finales.py

# 1 passed in 0.18s

El código se encuentra en GitHub

Comentarios

La longitud de la lista xs se calcula
- en Haskell, con length xs y
- en Python, con len(xs).

Interior de una lista

PorJosé A. Alonso 22-agosto-202214-diciembre-2022

Definir la función

   interior :: [a] -> [a]

1	interior :: [a] -> [a]

tal que (interior xs) es la lista obtenida eliminando los extremos de la lista xs. Por ejemplo,

   interior [2,5,3,7,3]  ==  [5,3,7]
   interior [2..7]       ==  [3,4,5,6]

1 2	interior [2,5,3,7,3] == [5,3,7] interior [2..7] == [3,4,5,6]

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

interior :: [a] -> [a]
interior xs = tail (init xs)

1 2	interior :: [a] -> [a] interior xs = tail (init xs)

El código se encuentra en GitHub.

Soluciones en Python

from typing import TypeVar
from hypothesis import given, strategies as st

A = TypeVar('A')

# 1ª solución
def interior1(xs: list[A]) -> list[A]:
    return xs[1:][:-1]

# 2ª solución
def interior2(xs: list[A]) -> list[A]:
    return xs[1:-1]

# La propiedad de equivalencia es
@given(st.lists(st.integers()))
def test_triangular(xs):
    assert interior1(xs) == interior2(xs)

# La comprobación es
#    src> poetry run pytest -q interior_de_una_lista.py
#    1 passed in 0.21s

from typing import TypeVar

from hypothesis import given, strategies as st

A = TypeVar('A')

# 1ª solución

def interior1(xs: list[A]) -> list[A]:

return xs[1:][:-1]

# 2ª solución

def interior2(xs: list[A]) -> list[A]:

return xs[1:-1]

# La propiedad de equivalencia es

@given(st.lists(st.integers()))

def test_triangular(xs):

assert interior1(xs) == interior2(xs)

# La comprobación es

# src> poetry run pytest -q interior_de_una_lista.py

# 1 passed in 0.21s

El código se encuentra en GitHub.

Comentarios

Los elementos iniciales de una lista xs se calcula
- en Haskell, con init xs y
- en Python, con xs[:-1].

Reconocimiento de palíndromos

PorJosé A. Alonso 19-agosto-202214-diciembre-2022

Definir la función

   palindromo :: Eq a => [a] -> Bool

1	palindromo :: Eq a => [a] -> Bool

tal que (palindromo xs) se verifica si xs es un palíndromo; es decir, es lo mismo leer xs de izquierda a derecha que de derecha a izquierda. Por ejemplo,

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

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

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

palindromo :: Eq a => [a] -> Bool
palindromo xs =
  xs == reverse xs

palindromo :: Eq a => [a] -> Bool

palindromo xs =

xs == reverse xs

El código se encuentra en GitHub.

Soluciones en Python

from typing import TypeVar

A = TypeVar('A')

def palindromo(xs: list[A]) -> bool:
    return xs == list(reversed(xs))

from typing import TypeVar

A = TypeVar('A')

def palindromo(xs: list[A]) -> bool:

return xs == list(reversed(xs))

El código se encuentra en GitHub.

Comentarios

La inversa de la lista xs se calcula
- en Haskell, con reverse xs y
- en Python, con list(reversed(xs)).
Para comparar la igualdad de dos listas xs e ys se escribe igual qh Haskell y en Python: xs == ys.

Rango de una lista

PorJosé A. Alonso 18-agosto-202214-diciembre-2022

Definir la función

   rango :: [Int] -> [Int]

1	rango :: [Int] -> [Int]

tal que (rango xs) es la lista formada por el menor y mayor elemento de xs. Por ejemplo,

   rango [3,2,7,5]  ==  [2,7]

1	rango [3,2,7,5] == [2,7]

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

rango :: [Int] -> [Int]
rango xs = [minimum xs, maximum xs]

1 2	rango :: [Int] -> [Int] rango xs = [minimum xs, maximum xs]

El código se encuentra en GitHub.

Soluciones en Python

def rango(xs: list[int]) -> list[int]:
    return [min(xs), max(xs)]

1 2	def rango(xs: list[int]) -> list[int]: return [min(xs), max(xs)]

El código se encuentra en GitHub.

Comentarios

El menor elemento de la lista xs se calcula
- en Haskell, con minimum xs y
- en Python, con min(xs).
El mayor elemento de la lista xs se calcula
- en Haskell, con maximum xs y
- en Python, con max(xs).

Los primeros al final

PorJosé A. Alonso 17-agosto-202214-diciembre-2022

Definir la función

   rota :: Int -> [a] -> [a]

1	rota :: Int -> [a] -> [a]

tal que (rota n xs) es la lista obtenida poniendo los n primeros elementos de xs al final de la lista. Por ejemplo,

   rota 1 [3,2,5,7]  ==  [2,5,7,3]
   rota 2 [3,2,5,7]  ==  [5,7,3,2]
   rota 3 [3,2,5,7]  ==  [7,3,2,5]

rota 1 [3,2,5,7] == [2,5,7,3]

rota 2 [3,2,5,7] == [5,7,3,2]

rota 3 [3,2,5,7] == [7,3,2,5]

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

rota :: Int -> [a] -> [a]
rota n xs = drop n xs ++ take n xs

1 2	rota :: Int -> [a] -> [a] rota n xs = drop n xs ++ take n xs

El código se encuentra en GitHub.

Soluciones en Python

from typing import TypeVar

A = TypeVar('A')

def rota(n: int, xs: list[A]) -> list[A]:
    return xs[n:] + xs[:n]

from typing import TypeVar

A = TypeVar('A')

def rota(n: int, xs: list[A]) -> list[A]:

return xs[n:] + xs[:n]

El código se encuentra en GitHub.

Comentarios
+ Los n primeros elementos de la lista xs se calcula
+ en Haskell, con take n xs y
+ en Python, con xs[n:].
+ La lista xs sin sus n primeros elementos se calcula
+ en Haskell, con drop n xs y
+ en Python, con xs[:n].

El primero al final

PorJosé A. Alonso 16-agosto-202214-diciembre-2022

Definir la función

   rota1 :: [a] -> [a]

1	rota1 :: [a] -> [a]

tal que (rota1 xs) es la lista obtenida poniendo el primer elemento de xs al final de la lista. Por ejemplo,

   rota1 [3,2,5,7]  ==  [2,5,7,3]

1	rota1 [3,2,5,7] == [2,5,7,3]

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

import Test.QuickCheck

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

rota1a :: [a] -> [a]
rota1a [] = []
rota1a xs = tail xs ++ [head xs]

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

rota1b :: [a] -> [a]
rota1b []     = []
rota1b (x:xs) = xs ++ [x]

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

-- La propiedad es
prop_rota1 :: [Int] -> Bool
prop_rota1 xs =
  rota1a xs == rota1b xs

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

import Test.QuickCheck

-- 1ª solución

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

rota1a :: [a] -> [a]

rota1a [] = []

rota1a xs = tail xs ++ [head xs]

-- 2ª solución

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

rota1b :: [a] -> [a]

rota1b [] = []

rota1b (x:xs) = xs ++ [x]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_rota1 :: [Int] -> Bool

prop_rota1 xs =

rota1a xs == rota1b xs

-- La comprobación es

-- λ> quickCheck prop_rota1

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Soluciones en Python

from typing import TypeVar
from hypothesis import given, strategies as st

A = TypeVar('A')

# 1ª solución
def rota1a(xs: list[A]) -> list[A]:
    if xs == []:
        return []
    return xs[1:] + [xs[0]]

# 2ª solución
def rota1b(xs: list[A]) -> list[A]:
    if xs == []:
        return []
    ys = xs[1:]
    ys.append(xs[0])
    return ys

# 3ª solución
def rota1c(xs: list[A]) -> list[A]:
    if xs == []:
        return []
    y, *ys = xs
    return ys + [y]

# La equivalencia de las definiciones es
@given(st.lists(st.integers()))
def test_rota1(xs: list[int]) -> None:
    assert rota1a(xs) == rota1b(xs) == rota1c(xs)

# La comprobación es
#    src> poetry run pytest -q el_primero_al_final.py
#    1 passed in 0.20s

from typing import TypeVar

from hypothesis import given, strategies as st

A = TypeVar('A')

# 1ª solución

def rota1a(xs: list[A]) -> list[A]:

if xs == []:

return []

return xs[1:] + [xs[0]]

# 2ª solución

def rota1b(xs: list[A]) -> list[A]:

if xs == []:

return []

ys = xs[1:]

ys.append(xs[0])

return ys

# 3ª solución

def rota1c(xs: list[A]) -> list[A]:

if xs == []:

return []

y, *ys = xs

return ys + [y]

# La equivalencia de las definiciones es

@given(st.lists(st.integers()))

def test_rota1(xs: list[int]) -> None:

assert rota1a(xs) == rota1b(xs) == rota1c(xs)

# La comprobación es

# src> poetry run pytest -q el_primero_al_final.py

# 1 passed in 0.20s

El código se encuentra en GitHub.

Comentarios

El primer elemento de la lista xs se calcula
- en Haskell, con head xs
- en Python, con xs[0]
El resto de la lista xs se calcula
- en Haskell, con tail xs
- en Python, con xs[1:]
La concatenación de las listas xse ysse calcula
- en Haskell, con xs ++ ys
- en Python, con xs + ys
En Python. xs.append(x) modifica la lista xs añadiéndole x al final. Por ejemplo,

>>> xs = [3, 2, 5]
>>> xs.append(1)
>>> xs
[3, 2, 5, 1]

>>> xs = [3, 2, 5]

>>> xs.append(1)

>>> xs

[3, 2, 5, 1]

Máximo de tres números

PorJosé A. Alonso 15-agosto-202214-diciembre-2022

Definir la función

   maxTres :: Int -> Int -> Int -> Int

1	maxTres :: Int -> Int -> Int -> Int

tal que (maxTres x y z) es el máximo de x, y y z. Por ejemplo,

   maxTres 6 2 4  ==  6
   maxTres 6 7 4  ==  7
   maxTres 6 7 9  ==  9

maxTres 6 2 4 == 6

maxTres 6 7 4 == 7

maxTres 6 7 9 == 9

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Solución en Haskell

maxTres :: Int -> Int -> Int -> Int
maxTres x y z = max x (max y z)

1 2	maxTres :: Int -> Int -> Int -> Int maxTres x y z = max x (max y z)

El código se encuentra en GitHub.

Soluciones en Python

def maxTres(x: int, y: int, z: int) -> int:
    return max(x, max(y, z))

1 2	def maxTres(x: int, y: int, z: int) -> int: return max(x, max(y, z))

El código se encuentra en GitHub.

Comentarios

El máximo de x e y se escribe
- en Haskell, max x y y
- en Python, max(x, y).

Último dígito

PorJosé A. Alonso 12-agosto-202214-diciembre-2022

Definir la función

   ultimoDigito :: Int -> Int

1	ultimoDigito :: Int -> Int

tal que (ultimoDigito x) es el último dígito del número x. Por ejemplo,

   ultimoDigito 325  ==  5

1	ultimoDigito 325 == 5

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Solución en Haskell

ultimoDigito :: Int -> Int
ultimoDigito x = rem x 10

1 2	ultimoDigito :: Int -> Int ultimoDigito x = rem x 10

El código se encuentra en GitHub.

Solución en Python

def ultimoDigito(x: int) -> int:
    return x % 10

1 2	def ultimoDigito(x: int) -> int: return x % 10

El código se encuentra en GitHub.

Comentarios

El resto de la división entera se x entre y sn ecribe
- en Haskell, rem x y,
- en Python, x % y.

Área de la corona circular

PorJosé A. Alonso 11-agosto-202214-diciembre-2022

Definir la función

   areaDeCoronaCircular :: Double -> Double -> Double

1	areaDeCoronaCircular :: Double -> Double -> Double

tal que (areaDeCoronaCircular r1 r2) es el área de una corona circular de radio interior r1 y radio exterior r2. Por ejemplo,

   areaDeCoronaCircular 1 2 == 9.42477796076938
   areaDeCoronaCircular 2 5 == 65.97344572538566
   areaDeCoronaCircular 3 5 == 50.26548245743669

areaDeCoronaCircular 1 2 == 9.42477796076938

areaDeCoronaCircular 2 5 == 65.97344572538566

areaDeCoronaCircular 3 5 == 50.26548245743669

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Solución en Haskell

areaDeCoronaCircular :: Double -> Double -> Double
areaDeCoronaCircular r1 r2 = pi*(r2^2 -r1^2)

1 2	areaDeCoronaCircular :: Double -> Double -> Double areaDeCoronaCircular r1 r2 = pi*(r2^2 -r1^2)

El código se encuentra en GitHub.

Solución en Python

from math import pi

def areaDeCoronaCircular(r1: float, r2: float) -> float:
    return pi * (r2 ** 2 - r1 ** 2)

from math import pi

def areaDeCoronaCircular(r1: float, r2: float) -> float:

return pi * (r2 ** 2 - r1 ** 2)

El código se encuentra en GitHub.

Comentarios

La diferencia de dos números x e y se escribe en Python igual que en Haskell: x - y.

Volumen de la esfera

PorJosé A. Alonso 10-agosto-202214-diciembre-2022

Definir la función

   volumenEsfera :: Double -> Double

1	volumenEsfera :: Double -> Double

tal que (volumenEsfera r) es el volumen de la esfera de radio r. Por ejemplo,

   volumenEsfera 10  ==  4188.790204786391

1	volumenEsfera 10 == 4188.790204786391

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Solución en Haskell

volumenEsfera :: Double -> Double
volumenEsfera r = (4/3)*pi*r^3

1 2	volumenEsfera :: Double -> Double volumenEsfera r = (4/3)pir^3

El código se encuentra en GitHub.

Solución en Python

from math import pi

def volumenEsfera(r: float) -> float:
    return (4 / 3) * pi * r ** 3

from math import pi

def volumenEsfera(r: float) -> float:

return (4 / 3) * pi * r ** 3

El código se encuentra en GitHub.

Comentarios

El número $\pi$ se representa igual en Python que en Haskell; pero, en Python. para usarlo hay que importarlo de la librería math.
La potencia de número x elevado al entero n se escribe
- en Haskell, x^n y
- en Python, x ** n.

Suma de monedas

PorJosé A. Alonso 9-agosto-202214-diciembre-2022

Definir la función

   sumaMonedas :: Int -> Int -> Int -> Int -> Int -> Int

1	sumaMonedas :: Int -> Int -> Int -> Int -> Int -> Int

tal que (sumaMonedas a b c d e) es la suma de los euros correspondientes a a monedas de 1 euro, b de 2 euros, c de 5 euros, d de 10 euros y e de 20 euros. Por ejemplo,

   sumaMonedas 0 0 0 0 1  ==  20
   sumaMonedas 0 0 8 0 3  == 100
   sumaMonedas 1 1 1 1 1  ==  38

sumaMonedas 0 0 0 0 1 == 20

sumaMonedas 0 0 8 0 3 == 100

sumaMonedas 1 1 1 1 1 == 38

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

sumaMonedas :: Int -> Int -> Int -> Int -> Int -> Int
sumaMonedas a b c d e = 1*a+2*b+5*c+10*d+20*e

1 2	sumaMonedas :: Int -> Int -> Int -> Int -> Int -> Int sumaMonedas a b c d e = 1a+2b+5c+10d+20*e

El código se encuentra en GitHub.

Soluciones en Python

def sumaMonedas(a: int, b: int, c: int, d: int, e: int) -> int:
    return 1 * a + 2 * b + 5 * c + 10 * d + 20 * e

1 2	def sumaMonedas(a: int, b: int, c: int, d: int, e: int) -> int: return 1 * a + 2 * b + 5 * c + 10 * d + 20 * e

El código se encuentra en GitHub.

Comentarios

El producto de dos números x e y se escribe en Python igual que en Haskell: x * y.

Media aritmética de tres números

PorJosé A. Alonso 8-agosto-202214-diciembre-2022

Definir la función

   media3 :: Float -> Float -> Float -> Float

1	media3 :: Float -> Float -> Float -> Float

tal que (media3 x y z) es la media aritmética de los números x, y y z. Por ejemplo,

   media3 1 3 8     ==  4.0
   media3 (-1) 0 7  ==  2.0
   media3 (-3) 0 3  ==  0.0

media3 1 3 8 == 4.0

media3 (-1) 0 7 == 2.0

media3 (-3) 0 3 == 0.0

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Solución en Haskell

media3 :: Float -> Float -> Float -> Float
media3 x y z = (x+y+z)/3

1 2	media3 :: Float -> Float -> Float -> Float media3 x y z = (x+y+z)/3

El código se encuentra en GitHub.

Soluciones en Python

def media3(x: float, y: float, z: float) -> float:
    return (x + y + z)/3

1 2	def media3(x: float, y: float, z: float) -> float: return (x + y + z)/3

El código se encuentra en GitHub.

Comentarios

Los comentarios irán resaltando las diferencias de la solución en Python respecto de la de Haskell (que no se hayan comentado en ningún ejercicio anterior).
La estructura de la definición en Python es

def <nombre de la función>(<argumento 1>,...,<argumento n>):
    # type: <signatura de la función>
    ...
    return <resultado>

def <nombre de la función>(<argumento 1>,...,<argumento n>):

# type: <signatura de la función>

...

return <resultado>

La suma de dos números x e y se escribe en Python igual que en Haskell: x + y.
El cociente de dos números decimales x e y se escribe en Python igual que en Haskell: x / y.

Curso de introducción a la programación con Haskell y Python

PorJosé A. Alonso 7-agosto-202210-agosto-2022

Mañana comenzará en Exercitium un curso práctico de introducción a la programación con Haskell y Python.

Diariamente, se publicará un ejercicio con sus soluciones en Haskell y en Python. El orden de los ejercicios se corresponde con el de los temas del curso de Programación funcional con Haskell. Además, en cada ejercicio se comentarán las diferencias entre ambos lenguajes y se irá extendiendo la tabla de equivalencia entre Haskell y Python.

Las soluciones de los ejercicios con Haskell se irán añadiendo al repositorio Exercitium em GitHub. Están como un proyecto con Stack (en concreto, la LTS-18.24) de forma que se puede verificar los ejemplos con

stack test

1	stack test

Las soluciones de los ejercicios con Python se irán añadiendo al repositorio Exercitium-Python em GitHub. Están como un proyecto con Poetry que usa la versión 3.10 de Python. También se pueden verificar los ejemplos con

poetry run pytest

1	poetry run pytest

He intentado escribir soluciones en Python parecidas a la de Haskell, conservando el estilo funcional. De todas formas, en los comentarios se pueden escribir definiciones alternativas en cualquiera de los lenguajes.

Los libros de Haskell en los que me he basado son

Algorithms: A functional programming approach de F. Rabhi y G. Lapalme.
Algorithm design with Haskell de E. Bird y J. Gibbons.
An introduction to functional programming systems using Haskell de A.J.T. Davie.
An introduction to functional programming through lambda calculus de G. Michaelson.
¡Aprende Haskell por el bien de todos! de M. Lipovača.
Beautiful code, compelling evidence (Functional programming for information visualization) de J.R. Heard.
Beginning Haskell (A project based approach). ~ A. Serrano.
Categories and Haskell de J.W. Buurlage.
Category theory applied to functional programming J.P. Villa.
Computational semantics and type theory de J. van Eijck.
Computational semantics with functional programming de J. van Eijck y C. Unger.
Developing web applications with Haskell and Yesod de M. Snoyman.
Discrete mathematics using a computer de J. O’Donnell, C. Hall y R. Page.
Domain-specific languages of Mathematics de P. Jansson, C. Ionescu y J.P. Bernardy.
Ejercicios de programación funcional con Haskell de J.A. Alonso y M.J. Hidalgo.
Exámenes de programación funcional con Haskell de J.A. Alonso y als.
Functional programming and graph algorithms de D.J. King.
Get programming with Haskell de W. Kurt.
Haskell cookbook de Y. Sajanikar.
Haskell data analysis cookbook de N. Shukla.
Haskell design patterns de R. Lemmer.
Haskell financial data modeling and predictive analytics de P. Ryzhov.
Haskell high performance programming de S. Thomasson.
Haskell (Notes for professionals).
Haskell programming from first principles de C. Allen y J. Moronuki.
Haskell quick syntax reference de S.L. Nita y M. Mihailescu.
Haskell tutorial and cookbook de M. Watson.
Haskell (The craft of functional programming) de S. Thompson.
Introducción a la programación funcional con Haskell de R. Bird.
Introduction to computation Haskell, logic and automata de D. Sannella, M. Fourman, H. Peng y P. Wadler.
Learn you a Haskell for great good! de M. Lipovaca y J. Brock.
Learning Haskell data analysis de J. Church.
Parallel and concurrent programming in Haskell. S. Marlow.
Pearls of functional algorithm design de R. Bird.
Piensa en Haskell de J.A. Alonso y M.J. Hidalgo.
Practical Haskell (A real World guide to programming) de A. Serrano.
Practical Web development with Haskell de E. Putrady.
Practical concurrent Haskell (with big data applications) de S.L. Nita y M. Mihailescu.
Programming in Haskell (2nd edition) de G. Hutton.
Purely functional data structures de D. Okasaki.
Real world Haskell de B. O’Sullivan, J. Goerzen y D. Stewart.
Temas de programación funcional de J.A. Alonso.
The Dao of functional programming de B. Milewski.
The Haskell school of music (from signals to symphonies) de P. Hudak.
The Haskell road to logic, maths and programming de K. Doets y J van Eijck.
Thinking functionally with Haskell de R. Bird.
Thinking with types de S. Maguire.
Yet another Haskell tutorial de Hal Daume III et al.

Los libros de Python en los que me he basado son

Advanced guide to Python 3 programming de J. Hunt.
Algorithmic problem solving with Python de J.B. Schneider, S.L. Broschat y J. Dahmen.
Aprenda a pensar como un programador con Python de A. Downey, J. Elkner y C. Meyers.
Aprende a programar en Python como si estuvieras en el siglo XXI de J.J. Merelo.
Artificial intelligence with Python de P. Joshi.
Automate the boring stuff with Python de A. Sweigart.
A beginners guide to Python 3 programming de J. Hunt.
A concise introduction to programming in Python de M.J. Johnson.
A functional start to computing with Python de T. Herman.
A primer on scientific programming with Python de H.P. Langtangen.
A programmer’s introduction to mathematics de J. Kun.
A very basic introduction to scientific Python programming de H.P. Langtangen.
Basic cheat sheet for Python de Swigartpdf
Beginning Python (From novice to professional) de M.L. Hetland.
Clean Python (Elegant coding in Python) de S. Kapil.
Competitive programming in Python de C. Durr y J.L. Vie.
Computational physics with Python de E. Ayars.
Computational quantum mechanics de J. Izaac y J. Wang.
Data science from scratch (First principles with Python) de J. Grus.
Data structures and algorithms in Python de M.T. Goodrich, R. Tamassia y M.H.Goldwasser.
Doing math with Python de A. Saha.
Effective computation in physics de A. Scopatz y K.D. Huff.
El lenguaje Python de D. Masip.
Fluent Python de L. Ramalho.
Functional Python programming de D.F. Lott.
Functional programming in Python de D. Mertz.
Fundamentals of Python (From first programs through data structures) de K.A. Lambert.
Hacking secret ciphers with Python de A. Sweigart.
Introducción a la programación con Python 3 de A. Marzal. I. Gracia y P. García.
Introduction a l’algorithmique avec Python de P. Olivier.
Introduction to Python for computational science and engineering de H. ~ Fangohr.
Introduction to Python for econometrics, statistics and data analysis. ~ Sheppard-14.pdf
Introduction to Python for science and engineering de D.J. Pine.
Introduction to Python for scientists and engineers de S. Nagar.
Introduction to computation and programming using Python de J.V. Guttag.
Introduction to recursive programming de M. Rubio.
Learning to program using Python de C. Jackson.
Lectures on mathematical computing with Python de J. Gopalakrishnan.
Mastering Python de R. van Hattem.
Problem solving with algorithms and data structures using Python de B. Miller y D. Ranum.
Programmation en Python pour les mathématiques de A. Casamayou P. Chauvin y G. Connan.
Programming and mathematical thinking de A.M. Stavely.
Programming for computations (A gentle introduction to numerical simulations with Python 3.6) de S. Linge y H.P. Langtangen.
Programming with Python de T.R. Padmanabhan.
Python 3 for absolute beginners de T. Hall y J-P Stacey.
Python algorithms (Mastering basic algorithms in the Python language) de M.L. Lie.
Python descriptors de J. Zimmerman.
Python programming (An introduction to computer science) de J. Zelle.
Python scientific lecture notes de V. Haenel, E. Gouillart y G. Varoquaux..
Scipy Lecture Notes (One document to learn numerics, science, and data with Python) de G. Varoquaux et als.
Serious Python de J. Danjou.
The Python workbook (A brief introduction with exercises and solutions) de B. Stephenson.

Números de Pentanacci

PorJosé A. Alonso 5-agosto-202214-diciembre-2022

Los números de Fibonacci se definen mediante las ecuaciones

   F(0) = 0
   F(1) = 1
   F(n) = F(n-1) + F(n-2), si n > 1

F(0) = 0

F(1) = 1

F(n) = F(n-1) + F(n-2), si n > 1

Los primeros números de Fibonacci son

   0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, ...

1	0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, ...

Una generalización de los anteriores son los números de Pentanacci definidos por las siguientes ecuaciones

   P(0) = 0
   P(1) = 1
   P(2) = 1
   P(3) = 2
   P(4) = 4
   P(n) = P(n-1) + P(n-2) + P(n-3) + P(n-4) + P(n-5), si n > 4

P(0) = 0

P(1) = 1

P(2) = 1

P(3) = 2

P(4) = 4

P(n) = P(n-1) + P(n-2) + P(n-3) + P(n-4) + P(n-5), si n > 4

Los primeros números de Pentanacci son

  0, 1, 1, 2, 4, 8, 16, 31, 61, 120, 236, 464, 912, 1793, 3525, ...

1	0, 1, 1, 2, 4, 8, 16, 31, 61, 120, 236, 464, 912, 1793, 3525, ...

Definir la sucesión

   pentanacci :: [Integer]

1	pentanacci :: [Integer]

cuyos elementos son los números de Pentanacci. Por ejemplo,

   λ> take 15 pentanacci
   [0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]
   λ> (pentanacci !! (10^5)) `mod` (10^30) 
   482929150584077921552549215816
   231437922897686901289110700696
   λ> length (show (pentanacci !! (10^5)))
   29357

λ> take 15 pentanacci

[0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]

λ> (pentanacci !! (10^5)) `mod` (10^30)

482929150584077921552549215816

231437922897686901289110700696

λ> length (show (pentanacci !! (10^5)))

29357

Soluciones

import Data.List (zipWith5)
import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

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

pentanacci1 :: [Integer]
pentanacci1 = [pent n | n <- [0..]]

pent :: Integer -> Integer
pent 0 = 0
pent 1 = 1
pent 2 = 1
pent 3 = 2
pent 4 = 4
pent n = pent (n-1) + pent (n-2) + pent (n-3) + pent (n-4) + pent (n-5)

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

pentanacci2 :: [Integer]
pentanacci2 = 
  0 : 1 : 1 : 2 : 4 : zipWith5 f (r 0) (r 1) (r 2) (r 3) (r 4)
  where f a b c d e = a+b+c+d+e
        r n         = drop n pentanacci2

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

pentanacci3 :: [Integer]
pentanacci3 = p (0, 1, 1, 2, 4)
  where p (a, b, c, d, e) = a : p (b, c, d, e, a + b + c + d + e)

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

pentanacci4 :: [Integer]
pentanacci4 = 0: 1: 1: 2: 4: p pentanacci4
  where p (a:b:c:d:e:xs) = (a+b+c+d+e): p (b:c:d:e:xs)

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

-- La propiedad es
prop_pentanacci :: NonNegative Int -> Bool
prop_pentanacci (NonNegative n) =
  all (== pentanacci1 !! n)
      [pentanacci1 !! n,
       pentanacci2 !! n,
       pentanacci3 !! n]

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

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

-- La comparación es
--    λ> pentanacci1 !! 25
--    5976577
--    (3.18 secs, 1,025,263,896 bytes)
--    λ> pentanacci2 !! 25
--    5976577
--    (0.00 secs, 562,360 bytes)
--    
--    λ> length (show (pentanacci2 !! (10^5)))
--    29357
--    (1.04 secs, 2,531,259,408 bytes)
--    λ> length (show (pentanacci3 !! (10^5)))
--    29357
--    (1.00 secs, 2,548,868,384 bytes)
--    λ> length (show (pentanacci4 !! (10^5)))
--    29357
--    (0.96 secs, 2,580,065,520 bytes)

import Data.List (zipWith5)

import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

-- 1ª solución

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

pentanacci1 :: [Integer]

pentanacci1 = [pent n | n <- [0..]]

pent :: Integer -> Integer

pent 0 = 0

pent 1 = 1

pent 2 = 1

pent 3 = 2

pent 4 = 4

pent n = pent (n-1) + pent (n-2) + pent (n-3) + pent (n-4) + pent (n-5)

-- 2ª solución

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

pentanacci2 :: [Integer]

pentanacci2 =

0 : 1 : 1 : 2 : 4 : zipWith5 f (r 0) (r 1) (r 2) (r 3) (r 4)

where f a b c d e = a+b+c+d+e

r n = drop n pentanacci2

-- 3ª solución

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

pentanacci3 :: [Integer]

pentanacci3 = p (0, 1, 1, 2, 4)

where p (a, b, c, d, e) = a : p (b, c, d, e, a + b + c + d + e)

-- 4ª solución

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

pentanacci4 :: [Integer]

pentanacci4 = 0: 1: 1: 2: 4: p pentanacci4

where p (a:b:c:d:e:xs) = (a+b+c+d+e): p (b:c:d:e:xs)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_pentanacci :: NonNegative Int -> Bool

prop_pentanacci (NonNegative n) =

all (== pentanacci1 !! n)

[pentanacci1 !! n,

pentanacci2 !! n,

pentanacci3 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=25}) prop_pentanacci

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> pentanacci1 !! 25

-- 5976577

-- (3.18 secs, 1,025,263,896 bytes)

-- λ> pentanacci2 !! 25

-- 5976577

-- (0.00 secs, 562,360 bytes)

-- λ> length (show (pentanacci2 !! (10^5)))

-- 29357

-- (1.04 secs, 2,531,259,408 bytes)

-- λ> length (show (pentanacci3 !! (10^5)))

-- 29357

-- (1.00 secs, 2,548,868,384 bytes)

-- λ> length (show (pentanacci4 !! (10^5)))

-- 29357

-- (0.96 secs, 2,580,065,520 bytes)

El código se encuentra en GitHub.

Referencias

Tito III Piezas y Eric Weisstein, Pentanacci number.
N. J. A. Sloane, Sucesión A001591 de la OEIS.

Número de representaciones de n como suma de dos cuadrados

PorJosé A. Alonso 4-agosto-202214-diciembre-2022

Sea $n$ un número natural cuya factorización prima es
$$n = 2^{a} \times p(1)^{b(1)} \times \dots \times p(n)^{b(n)} \times q(1)^{c(1)} \times \dots \times q(m)^{c(m)}$$
donde los $p(i)$ son los divisores primos de $n$ congruentes con 3 módulo 4 y los $q(j)$ son los divisores primos de $n$ congruentes con 1 módulo 4. Entonces, el número de forma de descomponer $n$ como suma de dos
cuadrados es 0, si algún $b(i)$ es impar y es el techo (es decir, el número entero más próximo por exceso) de
$$\frac{(1+c(1)) \times \dots \times (1+c(m))}{2}$$
en caso contrario. Por ejemplo, el número
$$2^{3} \times (3^{9} \times 7^{8}) \times (5^{3} \times 13^{6})$$
no se puede descomponer como sumas de dos cuadrados (porque el exponente de 3 es impar) y el número
$$2^{3} \times (3^{2} \times 7^{8}) \times (5^{3} \times 13^{6})$$
tiene 14 descomposiciones como suma de dos cuadrados (porque los exponentes de 3 y 7 son pares y el techo de
$$\frac{(1+3) \times (1+6)}{2}$$
es 14).

Definir la función

   nRepresentaciones :: Integer -> Integer

1	nRepresentaciones :: Integer -> Integer

tal que (nRepresentaciones n) es el número de formas de representar n como suma de dos cuadrados. Por ejemplo,

   nRepresentaciones (2^3*3^9*5^3*7^8*13^6)        ==  0
   nRepresentaciones (2^3*3^2*5^3*7^8*13^6)        ==  14
   head [n | n <- [1..], nRepresentaciones n > 8]  ==  71825

nRepresentaciones (2^3*3^9*5^3*7^8*13^6) == 0

nRepresentaciones (2^3*3^2*5^3*7^8*13^6) == 14

head [n | n <- [1..], nRepresentaciones n > 8] == 71825

Usando la función representaciones del ejercicio anterior, comprobar con QuickCheck la siguiente propiedad

   prop_representacion :: Positive Integer -> Bool
   prop_representacion (Positive n) =
     nRepresentaciones2 n == genericLength (representaciones n)

prop_representacion :: Positive Integer -> Bool

prop_representacion (Positive n) =

nRepresentaciones2 n == genericLength (representaciones n)

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

nRepresentaciones1 :: Integer -> Integer
nRepresentaciones1 n =
  ceiling (fromIntegral (product (map aux (factorizacion n))) / 2)
  where aux (p,e) | p == 2         = 1
                  | p `mod` 4 == 3 = if even e then 1 else 0
                  | otherwise      = e+1

-- (factorizacion n) es la factorización prima de n. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

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

nRepresentaciones2 :: Integer -> Integer
nRepresentaciones2 n =
  (1 + product (map aux (factorizacion n))) `div`  2
  where aux (p,e) | p == 2         = 1
                  | p `mod` 4 == 3 = if even e then 1 else 0
                  | otherwise      = e+1

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

-- La propiedad es
prop_nRepresentaciones :: Positive Integer -> Bool
prop_nRepresentaciones (Positive n) =
  nRepresentaciones1 n == nRepresentaciones2 n
  
-- La comprobación es
--    λ> quickCheck prop_nRepresentaciones
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> head [n | n <- [1..], nRepresentaciones1 n > 8]
--    71825
--    (1.39 secs, 2,970,063,760 bytes)
--    λ> head [n | n <- [1..], nRepresentaciones2 n > 8]
--    71825
--    (1.71 secs, 2,943,788,424 bytes)

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

-- La propiedad es
prop_representacion :: Positive Integer -> Bool
prop_representacion (Positive n) =
  nRepresentaciones2 n == genericLength (representaciones n)

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

esCuadrado :: Integer -> Bool
esCuadrado x = x == y * y
  where y = raiz x

raiz :: Integer -> Integer
raiz x = floor (sqrt (fromIntegral x))
    
-- La comprobación es
--    λ> quickCheck prop_representacion
--    +++ OK, passed 100 tests.

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

nRepresentaciones1 :: Integer -> Integer

nRepresentaciones1 n =

ceiling (fromIntegral (product (map aux (factorizacion n))) / 2)

where aux (p,e) | p == 2 = 1

| p `mod` 4 == 3 = if even e then 1 else 0

| otherwise = e+1

-- (factorizacion n) es la factorización prima de n. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

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

factorizacion n =

map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

-- 2ª solución

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

nRepresentaciones2 :: Integer -> Integer

nRepresentaciones2 n =

(1 + product (map aux (factorizacion n))) `div` 2

where aux (p,e) | p == 2 = 1

| p `mod` 4 == 3 = if even e then 1 else 0

| otherwise = e+1

-- Comprobación de equivalencia

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

-- La propiedad es

prop_nRepresentaciones :: Positive Integer -> Bool

prop_nRepresentaciones (Positive n) =

nRepresentaciones1 n == nRepresentaciones2 n

-- La comprobación es

-- λ> quickCheck prop_nRepresentaciones

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> head [n | n <- [1..], nRepresentaciones1 n > 8]

-- 71825

-- (1.39 secs, 2,970,063,760 bytes)

-- λ> head [n | n <- [1..], nRepresentaciones2 n > 8]

-- 71825

-- (1.71 secs, 2,943,788,424 bytes)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_representacion :: Positive Integer -> Bool

prop_representacion (Positive n) =

nRepresentaciones2 n == genericLength (representaciones n)

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

representaciones n =

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

let z = n - x*x,

esCuadrado z]

esCuadrado :: Integer -> Bool

esCuadrado x = x == y * y

where y = raiz x

raiz :: Integer -> Integer

raiz x = floor (sqrt (fromIntegral x))

-- La comprobación es

-- λ> quickCheck prop_representacion

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Referencias

W. Stein, Which numbers are the sum of two squares?.
E.W. Weisstein, Sum of squares function en MathWorld.
N.J.A. Sloane,Sucesión A004018 de OEIS.
Expressing a number as a sum of two squares.
Sum of squares.

Representaciones de un número como suma de dos cuadrados

PorJosé A. Alonso 3-agosto-202214-diciembre-2022

Definir la función

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

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

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

   representaciones  20              ==  [(2,4)]
   representaciones  25              ==  [(0,5),(3,4)]
   representaciones 325              ==  [(1,18),(6,17),(10,15)]
   length (representaciones (10^14)) == 8

representaciones 20 == [(2,4)]

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

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

length (representaciones (10^14)) == 8

Comprobar con QuickCheck que un número natural n se puede como suma de dos cuadrados si, y sólo si, en la factorización prima de n todos los exponentes de sus factores primos congruentes con 3 módulo 4 son pares.

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

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

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

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

-- (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 = floor . sqrt . fromIntegral
    
-- (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

-- 3ª solución
-- ===========
 
representaciones3 :: Integer -> [(Integer, Integer)]
representaciones3 n = aux 0 (raiz 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)

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

-- La propiedad es
prop_representaciones :: Positive Integer -> Bool
prop_representaciones (Positive n) =
  all (== representaciones1 n)
      [representaciones2 n,
       representaciones3 n]

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

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

-- La comparación es
--    λ> representaciones1 4000
--    [(20,60),(36,52)]
--    (4.95 secs, 2,434,929,624 bytes)
--    λ> representaciones2 4000
--    [(20,60),(36,52)]
--    (0.00 secs, 599,800 bytes)
--    λ> representaciones3 4000
--    [(20,60),(36,52)]
--    (0.01 secs, 591,184 bytes)
--    
--    λ> length (representaciones2 (10^14))
--    8
--    (6.64 secs, 5,600,837,088 bytes)
--    λ> length (representaciones3 (10^14))
--    8
--    (9.37 secs, 4,720,548,264 bytes)

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

-- La propiedad es
prop_representacion :: Positive Integer -> Bool
prop_representacion (Positive n) =
  not (null (representaciones2 n)) == 
  all (\(p,e) -> p `mod` 4 /= 3 || even e) (factorizacion n)

-- (factorizacion n) es la factorización prima de n. Por ejemplo,
--    factorizacion 600  ==  [(2,3),(3,1),(5,2)]
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

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

100

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

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

representaciones1 n =

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

-- 2ª solución

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

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

representaciones2 n =

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

let z = n - x*x,

esCuadrado z]

-- (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 = floor . sqrt . fromIntegral

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

-- 3ª solución

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

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

representaciones3 n = aux 0 (raiz 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)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_representaciones :: Positive Integer -> Bool

prop_representaciones (Positive n) =

all (== representaciones1 n)

[representaciones2 n,

representaciones3 n]

-- La comprobación es

-- λ> quickCheck prop_representaciones

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> representaciones1 4000

-- [(20,60),(36,52)]

-- (4.95 secs, 2,434,929,624 bytes)

-- λ> representaciones2 4000

-- [(20,60),(36,52)]

-- (0.00 secs, 599,800 bytes)

-- λ> representaciones3 4000

-- [(20,60),(36,52)]

-- (0.01 secs, 591,184 bytes)

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

-- 8

-- (6.64 secs, 5,600,837,088 bytes)

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

-- 8

-- (9.37 secs, 4,720,548,264 bytes)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_representacion :: Positive Integer -> Bool

prop_representacion (Positive n) =

not (null (representaciones2 n)) ==

all (\(p,e) -> p `mod` 4 /= 3 || even e) (factorizacion n)

-- (factorizacion n) es la factorización prima de n. Por ejemplo,

-- factorizacion 600 == [(2,3),(3,1),(5,2)]

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

factorizacion n =

map (\xs -> (head xs, genericLength xs)) (group (primeFactors n))

-- La comprobación es

-- λ> quickCheck prop_representacion

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Factorizaciones de números de Hilbert

PorJosé A. Alonso 2-agosto-202224-julio-2022

Un número de Hilbert es un entero positivo de la forma 4n+1. Los primeros números de Hilbert son 1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, 49, 53, 57, 61, 65, 69, …

Un primo de Hilbert es un número de Hilbert n que no es por ningún número de Hilbert menor que n (salvo el 1). Los primeros primos de Hilbert son 5, 9, 13, 17, 21, 29, 33, 37, 41, 49, 53, 57, 61, 69, 73, 77, 89, 93, 97, 101, 109, 113, 121, 129, 133, 137, …

Definir la función

   factorizacionesH :: Integer -> [[Integer]]

1	factorizacionesH :: Integer -> [[Integer]]

tal que (factorizacionesH n) es la listas de primos de Hilbert cuyo producto es el número de Hilbert n. Por ejemplo,

  factorizacionesH  25    ==  [[5,5]]
  factorizacionesH  45    ==  [[5,9]]
  factorizacionesH 441    ==  [[9,49],[21,21]]
  factorizacionesH 80109  ==  [[9,9,989],[9,69,129]]

factorizacionesH 25 == [[5,5]]

factorizacionesH 45 == [[5,9]]

factorizacionesH 441 == [[9,49],[21,21]]

factorizacionesH 80109 == [[9,9,989],[9,69,129]]

Comprobar con QuickCheck que todos los números de Hilbert son factorizables como producto de primos de Hilbert (aunque la factorización, como para el 441, puede no ser única).

Soluciones

import Data.Numbers.Primes (isPrime, primeFactors)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

factorizacionesH1 :: Integer -> [[Integer]]
factorizacionesH1 = aux primosH1
  where
    aux (x:xs) n 
      | x == n         = [[n]]
      | x > n          = []
      | n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n 
      | otherwise      = aux xs n 

primosH1 :: [Integer]
primosH1 = [n | n <- tail numerosH,
                divisoresH n == [1,n]]

-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,
--    take 15 numerosH  ==  [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]
numerosH :: [Integer]
numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a
-- n. Por ejemplo,
--   divisoresH 117  ==  [1,9,13,117]
--   divisoresH  21  ==  [1,21]
divisoresH :: Integer -> [Integer]
divisoresH n = [x | x <- takeWhile (<=n) numerosH,
                    n `mod` x == 0]

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

factorizacionesH2 :: Integer -> [[Integer]]
factorizacionesH2 = aux primosH2
  where
    aux (x:xs) n 
      | x == n         = [[n]]
      | x > n          = []
      | n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n 
      | otherwise      = aux xs n 

primosH2 :: [Integer]
primosH2 = filter esPrimoH (tail numerosH) 
  where esPrimoH n = all noDivideAn [5,9..m]
          where noDivideAn x = n `mod` x /= 0
                m            = ceiling (sqrt (fromIntegral n))

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

-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo 
-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)
-- (ver en https://bit.ly/3zq7h4e ).

factorizacionesH3 :: Integer -> [[Integer]]
factorizacionesH3 = aux primosH3
  where
    aux (x:xs) n 
      | x == n         = [[n]]
      | x > n          = []
      | n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n 
      | otherwise      = aux xs n 

primosH3 :: [Integer]
primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]
  where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs
          where xs = primeFactors n

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

-- La propiedad es
prop_factorizacionesH :: Positive Integer -> Bool
prop_factorizacionesH (Positive n) =
  all (== factorizacionesH1 m)
      [factorizacionesH2 m,
       factorizacionesH3 m]
  where m = 1 + 4 * n
  
-- La comprobación es
--    λ> quickCheck prop_factorizacionesH
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> factorizacionesH1 80109
--    [[9,9,989],[9,69,129]]
--    (42.77 secs, 14,899,787,640 bytes)
--    λ> factorizacionesH2 80109
--    [[9,9,989],[9,69,129]]
--    (0.26 secs, 156,051,104 bytes)
--    λ> factorizacionesH3 80109
--    [[9,9,989],[9,69,129]]
--    (0.35 secs, 1,118,236,536 bytes)

-- Propiedad de factorización
-- ==========================

-- La propiedad es
prop_factorizable :: Positive Integer -> Bool
prop_factorizable (Positive n) =
  not (null (factorizacionesH1 (1 + 4 * n)))

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

100

101

102

103

104

105

106

107

108

109

110

111

import Data.Numbers.Primes (isPrime, primeFactors)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

factorizacionesH1 :: Integer -> [[Integer]]

factorizacionesH1 = aux primosH1

where

aux (x:xs) n

| x == n = [[n]]

| x > n = []

| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n

| otherwise = aux xs n

primosH1 :: [Integer]

primosH1 = [n | n <- tail numerosH,

divisoresH n == [1,n]]

-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,

-- take 15 numerosH == [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]

numerosH :: [Integer]

numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a

-- n. Por ejemplo,

-- divisoresH 117 == [1,9,13,117]

-- divisoresH 21 == [1,21]

divisoresH :: Integer -> [Integer]

divisoresH n = [x | x <- takeWhile (<=n) numerosH,

n `mod` x == 0]

-- 2ª solución

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

factorizacionesH2 :: Integer -> [[Integer]]

factorizacionesH2 = aux primosH2

where

aux (x:xs) n

| x == n = [[n]]

| x > n = []

| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n

| otherwise = aux xs n

primosH2 :: [Integer]

primosH2 = filter esPrimoH (tail numerosH)

where esPrimoH n = all noDivideAn [5,9..m]

where noDivideAn x = n `mod` x /= 0

m = ceiling (sqrt (fromIntegral n))

-- 3ª solución

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

-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo

-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)

-- (ver en https://bit.ly/3zq7h4e ).

factorizacionesH3 :: Integer -> [[Integer]]

factorizacionesH3 = aux primosH3

where

aux (x:xs) n

| x == n = [[n]]

| x > n = []

| n `mod` x == 0 = map (x:) (aux (x:xs) (n `div` x) ) ++ aux xs n

| otherwise = aux xs n

primosH3 :: [Integer]

primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]

where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs

where xs = primeFactors n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_factorizacionesH :: Positive Integer -> Bool

prop_factorizacionesH (Positive n) =

all (== factorizacionesH1 m)

[factorizacionesH2 m,

factorizacionesH3 m]

where m = 1 + 4 * n

-- La comprobación es

-- λ> quickCheck prop_factorizacionesH

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> factorizacionesH1 80109

-- [[9,9,989],[9,69,129]]

-- (42.77 secs, 14,899,787,640 bytes)

-- λ> factorizacionesH2 80109

-- [[9,9,989],[9,69,129]]

-- (0.26 secs, 156,051,104 bytes)

-- λ> factorizacionesH3 80109

-- [[9,9,989],[9,69,129]]

-- (0.35 secs, 1,118,236,536 bytes)

-- Propiedad de factorización

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

-- La propiedad es

prop_factorizable :: Positive Integer -> Bool

prop_factorizable (Positive n) =

not (null (factorizacionesH1 (1 + 4 * n)))

-- La comprobación es

-- λ> quickCheck prop_factorizable

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Referencias

Basado en el artículo Failure of unique factorization (A example of the failure of the fundamental theorem of arithmetic) de R.J. Lipton en el blog Gödel’s Lost Letter and P=NP.

Otras referencias

Wikipedia, Hilbert number.
E.W. Weisstein, Hilbert number en MathWorld.
N.J.A. Sloane, Sucesión A057948 en la OEIS.
N.J.A. Sloane, Sucesión A057949 en la OEIS.

Números primos de Hilbert

PorJosé A. Alonso 1-agosto-202223-julio-2022

Un número de Hilbert es un positivo de la forma 4n+1. Los primeros números de Hilbert son 1, 5, 9, 13, 17, 21, 25, 29, 33, 37, 41, 45, 49, 53, 57, 61, 65, 69, 73, 77, 81, 85, 89, 93, 97, …

Definir la sucesión

   primosH :: [Integer]

1	primosH :: [Integer]

tal que sus elementos son los primos de Hilbert. Por ejemplo,

   take 15 primosH     == [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73]
   primosH !! (3*10^4) == 313661

1 2	take 15 primosH == [5,9,13,17,21,29,33,37,41,49,53,57,61,69,73] primosH !! (3*10^4) == 313661

Soluciones

import Data.Numbers.Primes (isPrime, primeFactors) 
import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

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

primosH1 :: [Integer]
primosH1 = [n | n <- tail numerosH,
                divisoresH n == [1,n]]

-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,
--    take 15 numerosH  ==  [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]
numerosH :: [Integer]
numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a
-- n. Por ejemplo,
--   divisoresH 117  ==  [1,9,13,117]
--   divisoresH  21  ==  [1,21]
divisoresH :: Integer -> [Integer]
divisoresH n = [x | x <- takeWhile (<=n) numerosH,
                    n `mod` x == 0]

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

primosH2 :: [Integer]
primosH2 = filter esPrimoH (tail numerosH) 
  where esPrimoH n = all noDivideAn [5,9..m]
          where noDivideAn x = n `mod` x /= 0
                m            = ceiling (sqrt (fromIntegral n))

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

-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo 
-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)
-- (ver en https://bit.ly/3zq7h4e ).

primosH3 :: [Integer]
primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]
  where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs
          where xs = primeFactors n

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

-- La propiedad es
prop_primosH :: NonNegative Int -> Bool
prop_primosH (NonNegative n) =
  all (== primosH1 !! n)
      [primosH2 !! n,
       primosH3 !! n]

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

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

-- La comparación es
--    λ> primosH1 !! 2000
--    16957
--    (2.16 secs, 752,085,752 bytes)
--    λ> primosH2 !! 2000
--    16957
--    (0.03 secs, 19,771,008 bytes)
--    λ> primosH3 !! 2000
--    16957
--    (0.07 secs, 152,029,168 bytes)
--    
--    λ> primosH2 !! (3*10^4)
--    313661
--    (1.44 secs, 989,761,888 bytes)
--    λ> primosH3 !! (3*10^4)
--    313661
--    (2.06 secs, 6,554,068,992 bytes)

import Data.Numbers.Primes (isPrime, primeFactors)

import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

-- 1ª solución

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

primosH1 :: [Integer]

primosH1 = [n | n <- tail numerosH,

divisoresH n == [1,n]]

-- numerosH es la sucesión de los números de Hilbert. Por ejemplo,

-- take 15 numerosH == [1,5,9,13,17,21,25,29,33,37,41,45,49,53,57]

numerosH :: [Integer]

numerosH = [1,5..]

-- (divisoresH n) es la lista de los números de Hilbert que dividen a

-- n. Por ejemplo,

-- divisoresH 117 == [1,9,13,117]

-- divisoresH 21 == [1,21]

divisoresH :: Integer -> [Integer]

divisoresH n = [x | x <- takeWhile (<=n) numerosH,

n `mod` x == 0]

-- 2ª solución

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

primosH2 :: [Integer]

primosH2 = filter esPrimoH (tail numerosH)

where esPrimoH n = all noDivideAn [5,9..m]

where noDivideAn x = n `mod` x /= 0

m = ceiling (sqrt (fromIntegral n))

-- 3ª solución

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

-- Basada en la siguiente propiedad: Un primo de Hilbert es un primo

-- de la forma 4n + 1 o un semiprimo de la forma (4a + 3) × (4b + 3)

-- (ver en https://bit.ly/3zq7h4e ).

primosH3 :: [Integer]

primosH3 = [ n | n <- numerosH, isPrime n || semiPrimoH n ]

where semiPrimoH n = length xs == 2 && all (\x -> (x-3) `mod` 4 == 0) xs

where xs = primeFactors n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosH :: NonNegative Int -> Bool

prop_primosH (NonNegative n) =

all (== primosH1 !! n)

[primosH2 !! n,

primosH3 !! n]

-- La comprobación es

-- λ> quickCheck prop_primosH

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosH1 !! 2000

-- 16957

-- (2.16 secs, 752,085,752 bytes)

-- λ> primosH2 !! 2000

-- 16957

-- (0.03 secs, 19,771,008 bytes)

-- λ> primosH3 !! 2000

-- 16957

-- (0.07 secs, 152,029,168 bytes)

-- λ> primosH2 !! (3*10^4)

-- 313661

-- (1.44 secs, 989,761,888 bytes)

-- λ> primosH3 !! (3*10^4)

-- 313661

-- (2.06 secs, 6,554,068,992 bytes)

El código se encuentra en GitHub.

Sumas de dos primos

PorJosé A. Alonso 29-julio-202222-julio-2022

Definir la sucesión

   sumasDeDosPrimos :: [Integer]

1	sumasDeDosPrimos :: [Integer]

cuyos elementos son los números que se pueden escribir como suma de dos números primos. Por ejemplo,

   λ> take 23 sumasDeDosPrimos
   [4,5,6,7,8,9,10,12,13,14,15,16,18,19,20,21,22,24,25,26,28,30,31]
   λ> sumasDeDosPrimos !! (5*10^5)
   862878

λ> take 23 sumasDeDosPrimos

[4,5,6,7,8,9,10,12,13,14,15,16,18,19,20,21,22,24,25,26,28,30,31]

λ> sumasDeDosPrimos !! (5*10^5)

862878

Soluciones

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

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

sumasDeDosPrimos1 :: [Integer]
sumasDeDosPrimos1 =
  [n | n <- [1..], not (null (sumaDeDosPrimos1 n))]

-- (sumaDeDosPrimos1 n) es la lista de pares de primos cuya suma es
-- n. Por ejemplo,
--    sumaDeDosPrimos  9  ==  [(2,7),(7,2)]
--    sumaDeDosPrimos 16  ==  [(3,13),(5,11),(11,5),(13,3)]
--    sumaDeDosPrimos 17  ==  []
sumaDeDosPrimos1 :: Integer -> [(Integer,Integer)]
sumaDeDosPrimos1 n = 
  [(x,n-x) | x <- primosN, isPrime (n-x)]
  where primosN = takeWhile (< n) primes

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

sumasDeDosPrimos2 :: [Integer]
sumasDeDosPrimos2 =
  [n | n <- [1..], not (null (sumaDeDosPrimos2 n))]

-- (sumasDeDosPrimos2 n) es la lista de pares (x,y) de primos cuya suma
-- es n y tales que x <= y. Por ejemplo,
--    sumaDeDosPrimos2  9  ==  [(2,7)]
--    sumaDeDosPrimos2 16  ==  [(3,13),(5,11)]
--    sumaDeDosPrimos2 17  ==  []
sumaDeDosPrimos2 :: Integer -> [(Integer,Integer)]
sumaDeDosPrimos2 n = 
  [(x,n-x) | x <- primosN, isPrime (n-x)]
  where primosN = takeWhile (<= (n `div` 2)) primes

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

sumasDeDosPrimos3 :: [Integer]
sumasDeDosPrimos3 = filter esSumaDeDosPrimos3 [4..]

-- (esSumaDeDosPrimos3 n) se verifica si n es suma de dos primos. Por
-- ejemplo, 
--    esSumaDeDosPrimos3  9  ==  True
--    esSumaDeDosPrimos3 16  ==  True
--    esSumaDeDosPrimos3 17  ==  False
esSumaDeDosPrimos3 :: Integer -> Bool
esSumaDeDosPrimos3 n
  | odd n     = isPrime (n-2)
  | otherwise = any isPrime [n-x | x <- takeWhile (<= (n `div` 2)) primes]

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

-- Usando la conjetura de Goldbach que dice que "Todo número par mayor
-- que 2 puede escribirse como suma de dos números primos" .

sumasDeDosPrimos4 :: [Integer]
sumasDeDosPrimos4 = filter esSumaDeDosPrimos4 [4..]

-- (esSumaDeDosPrimos4 n) se verifica si n es suma de dos primos. Por
-- ejemplo, 
--    esSumaDeDosPrimos4  9  ==  True
--    esSumaDeDosPrimos4 16  ==  True
--    esSumaDeDosPrimos4 17  ==  False
esSumaDeDosPrimos4 :: Integer -> Bool
esSumaDeDosPrimos4 n = even n || isPrime (n-2)

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

-- La propiedad es
prop_sumasDeDosPrimos :: NonNegative Int -> Bool
prop_sumasDeDosPrimos (NonNegative n) =
  all (== sumasDeDosPrimos1 !! n)
      [sumasDeDosPrimos2 !! n,
       sumasDeDosPrimos3 !! n,
       sumasDeDosPrimos4 !! n]
  
-- La comprobación es
--    λ> quickCheck prop_sumasDeDosPrimos
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> sumasDeDosPrimos1 !! 5000
--    7994
--    (2.61 secs, 9,299,106,792 bytes)
--    λ> sumasDeDosPrimos2 !! 5000
--    7994
--    (1.48 secs, 5,190,651,760 bytes)
--    λ> sumasDeDosPrimos3 !! 5000
--    7994
--    (0.12 secs, 351,667,104 bytes)
--    λ> sumasDeDosPrimos4 !! 5000
--    7994
--    (0.04 secs, 63,464,320 bytes)
--
--    λ> sumasDeDosPrimos3 !! (5*10^4)
--    83674
--    (2.23 secs, 7,776,049,264 bytes)
--    λ> sumasDeDosPrimos4 !! (5*10^4)
--    83674
--    (0.34 secs, 1,183,604,984 bytes)

100

101

102

103

104

105

106

107

108

import Data.Numbers.Primes (isPrime, primes)

import Test.QuickCheck

-- 1ª solución

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

sumasDeDosPrimos1 :: [Integer]

sumasDeDosPrimos1 =

[n | n <- [1..], not (null (sumaDeDosPrimos1 n))]

-- (sumaDeDosPrimos1 n) es la lista de pares de primos cuya suma es

-- n. Por ejemplo,

-- sumaDeDosPrimos 9 == [(2,7),(7,2)]

-- sumaDeDosPrimos 16 == [(3,13),(5,11),(11,5),(13,3)]

-- sumaDeDosPrimos 17 == []

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

sumaDeDosPrimos1 n =

[(x,n-x) | x <- primosN, isPrime (n-x)]

where primosN = takeWhile (< n) primes

-- 2ª solución

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

sumasDeDosPrimos2 :: [Integer]

sumasDeDosPrimos2 =

[n | n <- [1..], not (null (sumaDeDosPrimos2 n))]

-- (sumasDeDosPrimos2 n) es la lista de pares (x,y) de primos cuya suma

-- es n y tales que x <= y. Por ejemplo,

-- sumaDeDosPrimos2 9 == [(2,7)]

-- sumaDeDosPrimos2 16 == [(3,13),(5,11)]

-- sumaDeDosPrimos2 17 == []

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

sumaDeDosPrimos2 n =

[(x,n-x) | x <- primosN, isPrime (n-x)]

where primosN = takeWhile (<= (n `div` 2)) primes

-- 3ª solución

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

sumasDeDosPrimos3 :: [Integer]

sumasDeDosPrimos3 = filter esSumaDeDosPrimos3 [4..]

-- (esSumaDeDosPrimos3 n) se verifica si n es suma de dos primos. Por

-- ejemplo,

-- esSumaDeDosPrimos3 9 == True

-- esSumaDeDosPrimos3 16 == True

-- esSumaDeDosPrimos3 17 == False

esSumaDeDosPrimos3 :: Integer -> Bool

esSumaDeDosPrimos3 n

| odd n = isPrime (n-2)

| otherwise = any isPrime [n-x | x <- takeWhile (<= (n `div` 2)) primes]

-- 4ª solución

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

-- Usando la conjetura de Goldbach que dice que "Todo número par mayor

-- que 2 puede escribirse como suma de dos números primos" .

sumasDeDosPrimos4 :: [Integer]

sumasDeDosPrimos4 = filter esSumaDeDosPrimos4 [4..]

-- (esSumaDeDosPrimos4 n) se verifica si n es suma de dos primos. Por

-- ejemplo,

-- esSumaDeDosPrimos4 9 == True

-- esSumaDeDosPrimos4 16 == True

-- esSumaDeDosPrimos4 17 == False

esSumaDeDosPrimos4 :: Integer -> Bool

esSumaDeDosPrimos4 n = even n || isPrime (n-2)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sumasDeDosPrimos :: NonNegative Int -> Bool

prop_sumasDeDosPrimos (NonNegative n) =

all (== sumasDeDosPrimos1 !! n)

[sumasDeDosPrimos2 !! n,

sumasDeDosPrimos3 !! n,

sumasDeDosPrimos4 !! n]

-- La comprobación es

-- λ> quickCheck prop_sumasDeDosPrimos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> sumasDeDosPrimos1 !! 5000

-- 7994

-- (2.61 secs, 9,299,106,792 bytes)

-- λ> sumasDeDosPrimos2 !! 5000

-- 7994

-- (1.48 secs, 5,190,651,760 bytes)

-- λ> sumasDeDosPrimos3 !! 5000

-- 7994

-- (0.12 secs, 351,667,104 bytes)

-- λ> sumasDeDosPrimos4 !! 5000

-- 7994

-- (0.04 secs, 63,464,320 bytes)

-- λ> sumasDeDosPrimos3 !! (5*10^4)

-- 83674

-- (2.23 secs, 7,776,049,264 bytes)

-- λ> sumasDeDosPrimos4 !! (5*10^4)

-- 83674

-- (0.34 secs, 1,183,604,984 bytes)

El código se encuentra en GitHub.

Referencia

N.J.A. Sloane, Sucesión A014091 en OEIS.

La sucesión de Thue-Morse

PorJosé A. Alonso 28-julio-202221-julio-2022

La serie de Thue-Morse comienza con el término [0] y sus siguientes términos se construyen añadiéndole al anterior su complementario. Los primeros términos de la serie son

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

[0]

[0,1]

[0,1,1,0]

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

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

De esta forma se va formando una sucesión

   0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,...

1	0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,...

que se conoce como la sucesión de Thue-Morse.

Definir la sucesión

   sucThueMorse :: [Int]

1	sucThueMorse :: [Int]

cuyos elementos son los de la sucesión de Thue-Morse. Por ejemplo,

   λ> take 30 sucThueMorse
   [0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0]
   λ> map (sucThueMorse4 !!) [1234567..1234596] 
   [1,1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,1,0,0,1,0]
   λ> map (sucThueMorse4 !!) [4000000..4000030] 
   [1,0,0,1,0,1,1,0,0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,1,0,0,1,0,1,1]

λ> take 30 sucThueMorse

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

λ> map (sucThueMorse4 !!) [1234567..1234596]

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

λ> map (sucThueMorse4 !!) [4000000..4000030]

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

Comprobar con QuickCheck que si s(n) representa el término n-ésimo de la sucesión de Thue-Morse, entonces

   s(2n)   = s(n)
   s(2n+1) = 1 - s(n)

1 2	s(2n) = s(n) s(2n+1) = 1 - s(n)

Soluciones

import Test.QuickCheck

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

sucThueMorse1 :: [Int]
sucThueMorse1 = map termSucThueMorse1 [0..]

-- (termSucThueMorse1 n) es el n-ésimo término de la sucesión de
-- Thue-Morse. Por ejemplo, 
--    termSucThueMorse1 0  ==  0
--    termSucThueMorse1 1  ==  1
--    termSucThueMorse1 2  ==  1
--    termSucThueMorse1 3  ==  0
--    termSucThueMorse1 4  ==  1
termSucThueMorse1 :: Int -> Int
termSucThueMorse1 0 = 0
termSucThueMorse1 n = 
  (serieThueMorse !! k) !! n
  where k = 1 + floor (logBase 2 (fromIntegral n))

-- serieThueMorse es la lista cuyos elementos son los términos de la
-- serie de Thue-Morse. Por ejemplo, 
--    λ> take 4 serieThueMorse3
--    [[0],[0,1],[0,1,1,0],[0,1,1,0,1,0,0,1]]
serieThueMorse :: [[Int]]
serieThueMorse = iterate paso [0]
  where paso xs = xs ++ map (1-) xs

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

sucThueMorse2 :: [Int]
sucThueMorse2 = 
  0 : intercala (map (1-) sucThueMorse2) (tail sucThueMorse2)

-- (intercala xs ys) es la lista obtenida intercalando los elementos de
-- las listas infinitas xs e ys. Por ejemplo, 
--    take 10 (intercala [1,5..] [2,4..])  ==  [1,2,5,4,9,6,13,8,17,10]
intercala :: [a] -> [a] -> [a]
intercala (x:xs) ys = x : intercala ys xs 

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

sucThueMorse3 :: [Int]
sucThueMorse3 = 0 : 1 : aux (tail sucThueMorse3) 
  where aux (x : xs) = x : (1 - x) : aux xs



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

sucThueMorse4 :: [Int]
sucThueMorse4 = 0 : aux [1]
  where aux xs = xs ++ aux (xs ++ map (1-) xs) 

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

-- La propiedad es
prop_termSucThueMorse :: NonNegative Int -> Bool
prop_termSucThueMorse (NonNegative n) =
  sucThueMorse1 !! (2*n)   == sn &&
  sucThueMorse1 !! (2*n+1) == 1 - sn 
  where sn = sucThueMorse1 !! n

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

-- 5ª solución
-- ===========

sucThueMorse5 :: [Int]
sucThueMorse5 = map termSucThueMorse5 [0..]

-- (termSucThueMorse5 n) es el n-ésimo término de la sucesión de
-- Thue-Morse. Por ejemplo, 
--    termSucThueMorse5 0  ==  0
--    termSucThueMorse5 1  ==  1
--    termSucThueMorse5 2  ==  1
--    termSucThueMorse5 3  ==  0
--    termSucThueMorse5 4  ==  1
termSucThueMorse5 :: Int -> Int
termSucThueMorse5 0 = 0
termSucThueMorse5 n 
  | even n    = termSucThueMorse5 (n `div` 2)
  | otherwise = 1 - termSucThueMorse5 (n `div` 2)

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

-- La propiedad es
prop_sucThueMorse :: NonNegative Int -> Bool
prop_sucThueMorse (NonNegative n) =
  all (== sucThueMorse1 !! n)
      [sucThueMorse2 !! n,
       sucThueMorse3 !! n,
       sucThueMorse4 !! n,
       sucThueMorse5 !! n]

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

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

-- La comparación es
--    λ> sucThueMorse1 !! (10^7)
--    0
--    (3.28 secs, 3,420,080,168 bytes)
--    λ> sucThueMorse2 !! (10^7)
--    0
--    (3.01 secs, 1,720,549,640 bytes)
--    λ> sucThueMorse3 !! (10^7)
--    0
--    (1.80 secs, 1,360,550,040 bytes)
--    λ> sucThueMorse4 !! (10^7)
--    0
--    (0.88 secs, 1,254,772,768 bytes)
--    λ> sucThueMorse5 !! (10^7)
--    0
--    (0.62 secs, 1,600,557,072 bytes)

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

import Test.QuickCheck

-- 1ª solución

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

sucThueMorse1 :: [Int]

sucThueMorse1 = map termSucThueMorse1 [0..]

-- (termSucThueMorse1 n) es el n-ésimo término de la sucesión de

-- Thue-Morse. Por ejemplo,

-- termSucThueMorse1 0 == 0

-- termSucThueMorse1 1 == 1

-- termSucThueMorse1 2 == 1

-- termSucThueMorse1 3 == 0

-- termSucThueMorse1 4 == 1

termSucThueMorse1 :: Int -> Int

termSucThueMorse1 0 = 0

termSucThueMorse1 n =

(serieThueMorse !! k) !! n

where k = 1 + floor (logBase 2 (fromIntegral n))

-- serieThueMorse es la lista cuyos elementos son los términos de la

-- serie de Thue-Morse. Por ejemplo,

-- λ> take 4 serieThueMorse3

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

serieThueMorse :: [[Int]]

serieThueMorse = iterate paso [0]

where paso xs = xs ++ map (1-) xs

-- 2ª solución

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

sucThueMorse2 :: [Int]

sucThueMorse2 =

0 : intercala (map (1-) sucThueMorse2) (tail sucThueMorse2)

-- (intercala xs ys) es la lista obtenida intercalando los elementos de

-- las listas infinitas xs e ys. Por ejemplo,

-- take 10 (intercala [1,5..] [2,4..]) == [1,2,5,4,9,6,13,8,17,10]

intercala :: [a] -> [a] -> [a]

intercala (x:xs) ys = x : intercala ys xs

-- 3ª solución

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

sucThueMorse3 :: [Int]

sucThueMorse3 = 0 : 1 : aux (tail sucThueMorse3)

where aux (x : xs) = x : (1 - x) : aux xs

-- 4ª solución

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

sucThueMorse4 :: [Int]

sucThueMorse4 = 0 : aux [1]

where aux xs = xs ++ aux (xs ++ map (1-) xs)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_termSucThueMorse :: NonNegative Int -> Bool

prop_termSucThueMorse (NonNegative n) =

sucThueMorse1 !! (2*n) == sn &&

sucThueMorse1 !! (2*n+1) == 1 - sn

where sn = sucThueMorse1 !! n

-- La comprobación es

-- λ> quickCheck prop_termSucThueMorse

-- +++ OK, passed 100 tests.

-- 5ª solución

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

sucThueMorse5 :: [Int]

sucThueMorse5 = map termSucThueMorse5 [0..]

-- (termSucThueMorse5 n) es el n-ésimo término de la sucesión de

-- Thue-Morse. Por ejemplo,

-- termSucThueMorse5 0 == 0

-- termSucThueMorse5 1 == 1

-- termSucThueMorse5 2 == 1

-- termSucThueMorse5 3 == 0

-- termSucThueMorse5 4 == 1

termSucThueMorse5 :: Int -> Int

termSucThueMorse5 0 = 0

termSucThueMorse5 n

| even n = termSucThueMorse5 (n `div` 2)

| otherwise = 1 - termSucThueMorse5 (n `div` 2)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sucThueMorse :: NonNegative Int -> Bool

prop_sucThueMorse (NonNegative n) =

all (== sucThueMorse1 !! n)

[sucThueMorse2 !! n,

sucThueMorse3 !! n,

sucThueMorse4 !! n,

sucThueMorse5 !! n]

-- La comprobación es

-- λ> quickCheck prop_sucThueMorse

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> sucThueMorse1 !! (10^7)

-- 0

-- (3.28 secs, 3,420,080,168 bytes)

-- λ> sucThueMorse2 !! (10^7)

-- 0

-- (3.01 secs, 1,720,549,640 bytes)

-- λ> sucThueMorse3 !! (10^7)

-- 0

-- (1.80 secs, 1,360,550,040 bytes)

-- λ> sucThueMorse4 !! (10^7)

-- 0

-- (0.88 secs, 1,254,772,768 bytes)

-- λ> sucThueMorse5 !! (10^7)

-- 0

-- (0.62 secs, 1,600,557,072 bytes)

El código se encuentra en GitHub.

Referencias

N.J.A. Sloane Sucesión A010060 en OEIS.
Programming Praxis Thue-Morse sequence.
Wikipedia Thue–Morse sequence

La serie de Thue-Morse

PorJosé A. Alonso 27-julio-202220-julio-2022

La serie de Thue-Morse comienza con el término [0] y sus siguientes términos se construyen añadiéndole al anterior su complementario (es decir, la lista obtenida cambiando el 0 por 1 y el 1 por 0). Los primeros términos de la serie son

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

[0]

[0,1]

[0,1,1,0]

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

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

Definir la lista

   serieThueMorse :: [[Int]]

1	serieThueMorse :: [[Int]]

tal que sus elementos son los términos de la serie de Thue-Morse. Por ejemplo,

   λ> take 4 serieThueMorse
   [[0],[0,1],[0,1,1,0],[0,1,1,0,1,0,0,1]]

1 2	λ> take 4 serieThueMorse [[0],[0,1],[0,1,1,0],[0,1,1,0,1,0,0,1]]

Soluciones

import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

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

serieThueMorse1 :: [[Int]]
serieThueMorse1 = map termSerieThueMorse [0..]

-- (termSerieThueMorse n) es el término n-ésimo de la serie de
-- Thue-Morse. Por ejemplo, 
--    termSerieThueMorse 1  ==  [0,1]
--    termSerieThueMorse 2  ==  [0,1,1,0]
--    termSerieThueMorse 3  ==  [0,1,1,0,1,0,0,1]
--    termSerieThueMorse 4  ==  [0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0]
termSerieThueMorse :: Int -> [Int]
termSerieThueMorse 0 = [0]
termSerieThueMorse n = xs ++ complementaria xs
  where xs = termSerieThueMorse (n-1)

-- (complementaria xs) es la complementaria de la lista xs (formada por
-- ceros y unos); es decir, la lista obtenida cambiando el 0 por 1 y el
-- 1 por 0. Por ejemplo, 
--    complementaria [1,0,0,1,1,0,1]  ==  [0,1,1,0,0,1,0]
complementaria :: [Int] -> [Int]
complementaria = map (1-)

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

serieThueMorse2 :: [[Int]]
serieThueMorse2 = [0] : map paso serieThueMorse2
  where paso xs = xs ++ complementaria xs

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

serieThueMorse3 :: [[Int]]
serieThueMorse3 = iterate paso [0]
  where paso xs = xs ++ complementaria xs

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

-- Observando que cada término de la serie de Thue-Morse se obtiene del
-- anterior sustituyendo los 1 por 1, 0 y los 0 por  0, 1. 

serieThueMorse4 :: [[Int]]
serieThueMorse4 = [0] : map (concatMap paso4) serieThueMorse4
  where paso4 x = [x,1-x]

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

-- La propiedad es
prop_serieThueMorse :: NonNegative Int -> Bool 
prop_serieThueMorse  (NonNegative n) =
  all (== serieThueMorse1 !! n)
      [serieThueMorse2 !! n,
       serieThueMorse3 !! n,
       serieThueMorse4 !! n]

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

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

-- La comparación es
--    λ> (serieThueMorse1 !! 23) !! (2^22)
--    1
--    (1.08 secs, 839,419,224 bytes)
--    λ> (serieThueMorse2 !! 23) !! (2^22)
--    1
--    (0.61 secs, 839,413,592 bytes)
--    λ> (serieThueMorse3 !! 23) !! (2^22)
--    1
--    (1.43 secs, 839,413,592 bytes)
--    λ> (serieThueMorse4 !! 23) !! (2^22)
--    1
--    (1.57 secs, 1,007,190,024 bytes)

import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

-- 1ª solución

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

serieThueMorse1 :: [[Int]]

serieThueMorse1 = map termSerieThueMorse [0..]

-- (termSerieThueMorse n) es el término n-ésimo de la serie de

-- Thue-Morse. Por ejemplo,

-- termSerieThueMorse 1 == [0,1]

-- termSerieThueMorse 2 == [0,1,1,0]

-- termSerieThueMorse 3 == [0,1,1,0,1,0,0,1]

-- termSerieThueMorse 4 == [0,1,1,0,1,0,0,1,1,0,0,1,0,1,1,0]

termSerieThueMorse :: Int -> [Int]

termSerieThueMorse 0 = [0]

termSerieThueMorse n = xs ++ complementaria xs

where xs = termSerieThueMorse (n-1)

-- (complementaria xs) es la complementaria de la lista xs (formada por

-- ceros y unos); es decir, la lista obtenida cambiando el 0 por 1 y el

-- 1 por 0. Por ejemplo,

-- complementaria [1,0,0,1,1,0,1] == [0,1,1,0,0,1,0]

complementaria :: [Int] -> [Int]

complementaria = map (1-)

-- 2ª solución

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

serieThueMorse2 :: [[Int]]

serieThueMorse2 = [0] : map paso serieThueMorse2

where paso xs = xs ++ complementaria xs

-- 3ª solución

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

serieThueMorse3 :: [[Int]]

serieThueMorse3 = iterate paso [0]

where paso xs = xs ++ complementaria xs

-- 4ª solución

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

-- Observando que cada término de la serie de Thue-Morse se obtiene del

-- anterior sustituyendo los 1 por 1, 0 y los 0 por 0, 1.

serieThueMorse4 :: [[Int]]

serieThueMorse4 = [0] : map (concatMap paso4) serieThueMorse4

where paso4 x = [x,1-x]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_serieThueMorse :: NonNegative Int -> Bool

prop_serieThueMorse (NonNegative n) =

all (== serieThueMorse1 !! n)

[serieThueMorse2 !! n,

serieThueMorse3 !! n,

serieThueMorse4 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=20}) prop_serieThueMorse

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> (serieThueMorse1 !! 23) !! (2^22)

-- 1

-- (1.08 secs, 839,419,224 bytes)

-- λ> (serieThueMorse2 !! 23) !! (2^22)

-- 1

-- (0.61 secs, 839,413,592 bytes)

-- λ> (serieThueMorse3 !! 23) !! (2^22)

-- 1

-- (1.43 secs, 839,413,592 bytes)

-- λ> (serieThueMorse4 !! 23) !! (2^22)

-- 1

-- (1.57 secs, 1,007,190,024 bytes)

El código se encuentra en [GitHub](https://github.com/jaalonso/Exercitium/blob/main/src/

Referencias

N.J.A. Sloane Sucesión A010060 en OEIS.
Programming Praxis Thue-Morse sequence.
Wikipedia Thue–Morse sequence

Números belgas

PorJosé A. Alonso 26-julio-202219-julio-2022

Un número n es k-belga si la sucesión cuyo primer elemento es k y cuyos elementos se obtienen sumando reiteradamente las cifras de n contiene a n.

El 18 es 0-belga, porque a partir del 0 vamos a ir sumando sucesivamente 1, 8, 1, 8, … hasta llegar o sobrepasar el 18: 0, 1, 9, 10, 18, … Como se alcanza el 18, resulta que el 18 es 0-belga.

El 19 no es 1-belga, porque a partir del 1 vamos a ir sucesivamente 1, 9, 1, 9, … hasta llegar o sobrepasar el 18: 0, 1, 10, 11, 20, 21, … Como no se alcanza el 19, resulta que el 19 no es 1-belga.

Definir la función

   esBelga :: Int -> Int -> Bool

1	esBelga :: Int -> Int -> Bool

tal que (esBelga k n) se verifica si n es k-belga. Por ejemplo,

   esBelga 0 18                              ==  True
   esBelga 1 19                              ==  False
   esBelga 0 2016                            ==  True
   [x | x <- [0..30], esBelga 7 x]           ==  [7,10,11,21,27,29]
   [x | x <- [0..30], esBelga 10 x]          ==  [10,11,20,21,22,24,26]
   length [n | n <- [1..10^6], esBelga 0 n]  ==  272049

esBelga 0 18 == True

esBelga 1 19 == False

esBelga 0 2016 == True

[x | x <- [0..30], esBelga 7 x] == [7,10,11,21,27,29]

[x | x <- [0..30], esBelga 10 x] == [10,11,20,21,22,24,26]

length [n | n <- [1..10^6], esBelga 0 n] == 272049

Comprobar con QuickCheck que para todo número entero positivo n, si k es el resto de n entre la suma de los dígitos de n, entonces n es k-belga.

Soluciones

import Data.Char (digitToInt)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

esBelga1 :: Int -> Int -> Bool
esBelga1 k n =
  n == head (dropWhile (<n) (scanl (+) k (cycle (digitos n))))

digitos :: Int -> [Int]
digitos n = map digitToInt (show n)

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

esBelga2 :: Int -> Int -> Bool
esBelga2 k n =
  k <= n && n == head (dropWhile (<n) (scanl (+) (k + q * s) ds))
  where ds = digitos n
        s  = sum ds
        q  = (n - k) `div` s

-- Equivalencia
-- ============

-- La propiedad es
prop_esBelga :: Positive Int -> Positive Int -> Bool
prop_esBelga (Positive k) (Positive n) = 
  esBelga1 k n == esBelga2 k n

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

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

-- La comparación es
--    λ> length [n | n <- [1..2*10^4], esBelga1 0 n]
--    6521
--    (6.27 secs, 6,508,102,192 bytes)
--    λ> length [n | n <- [1..2*10^4], esBelga2 0 n]
--    6521
--    (0.07 secs, 46,741,144 bytes)

-- Verificación de la propiedad
-- ============================

-- La propiedad es
prop_Belga :: Positive Int -> Bool
prop_Belga (Positive n) = 
  esBelga2 k n
  where k = n `mod` sum (digitos n)

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

import Data.Char (digitToInt)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

esBelga1 :: Int -> Int -> Bool

esBelga1 k n =

n == head (dropWhile (<n) (scanl (+) k (cycle (digitos n))))

digitos :: Int -> [Int]

digitos n = map digitToInt (show n)

-- 2ª solución

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

esBelga2 :: Int -> Int -> Bool

esBelga2 k n =

k <= n && n == head (dropWhile (<n) (scanl (+) (k + q * s) ds))

where ds = digitos n

s = sum ds

q = (n - k) `div` s

-- Equivalencia

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

-- La propiedad es

prop_esBelga :: Positive Int -> Positive Int -> Bool

prop_esBelga (Positive k) (Positive n) =

esBelga1 k n == esBelga2 k n

-- La comprobación es

-- λ> quickCheck prop_esBelga

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length [n | n <- [1..2*10^4], esBelga1 0 n]

-- 6521

-- (6.27 secs, 6,508,102,192 bytes)

-- λ> length [n | n <- [1..2*10^4], esBelga2 0 n]

-- 6521

-- (0.07 secs, 46,741,144 bytes)

-- Verificación de la propiedad

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

-- La propiedad es

prop_Belga :: Positive Int -> Bool

prop_Belga (Positive n) =

esBelga2 k n

where k = n `mod` sum (digitos n)

-- La comprobación es

-- λ> quickCheck prop_Belga

-- +++ OK, passed 100 tests.

Referencias

Basado en el artículo Números belgas del blog Números y hoja de cálculo de Antonio Roldán Martínez.

El código se encuentra en GitHub.

Huecos maximales entre primos

PorJosé A. Alonso 25-julio-202218-julio-2022

El hueco de un número primo p es la distancia entre p y primo siguiente de p. Por ejemplo, el hueco de 7 es 4 porque el primo siguiente de 7 es 11 y 4 = 11-7. Los huecos de los primeros números son

   Primo Hueco
    2    1
    3    2
    7    4
   11    2

Primo Hueco

2 1

3 2

7 4

11 2

El hueco de un número primo p es maximal si es mayor que huecos de todos los números menores que p. Por ejemplo, 4 es un hueco maximal de 7 ya que los huecos de los primos menores que 7 son 1 y 2 y ambos son menores que 4. La tabla de los primeros huecos maximales es

   Primo Hueco
     2    1
     3    2
     7    4
    23    6
    89    8
   113   14
   523   18
   887   20

Primo Hueco

2 1

3 2

7 4

23 6

89 8

113 14

523 18

887 20

Definir la sucesión

   primosYhuecosMaximales :: [(Integer,Integer)]

1	primosYhuecosMaximales :: [(Integer,Integer)]

cuyos elementos son los números primos con huecos maximales junto son sus huecos. Por ejemplo,

   λ> take 8 primosYhuecosMaximales
   [(2,1),(3,2),(7,4),(23,6),(89,8),(113,14),(523,18),(887,20)]
   λ> primosYhuecosMaximales !! 20
   (2010733,148)

λ> take 8 primosYhuecosMaximales

[(2,1),(3,2),(7,4),(23,6),(89,8),(113,14),(523,18),(887,20)]

λ> primosYhuecosMaximales !! 20

(2010733,148)

Soluciones

import Data.Numbers.Primes (primes)
import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

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

primosYhuecosMaximales1 :: [(Integer,Integer)]
primosYhuecosMaximales1 = 
  [(p,huecoPrimo p) | p <- primes, esMaximalHuecoPrimo p]

-- (siguientePrimo x) es el menor primo mayor que x. Por ejemplo,
--    siguientePrimo 7  ==  11
--    siguientePrimo 8  ==  11
siguientePrimo :: Integer -> Integer
siguientePrimo p =
  head (dropWhile (<= p) primes)

-- (huecoPrimo p) es la distancia del primo p hasta el siguiente
-- primo. Por ejemplo,
--    huecoPrimo 7  ==  4
huecoPrimo :: Integer -> Integer
huecoPrimo p = siguientePrimo p - p

-- (esMaximalHuecoPrimo p) se verifica si el hueco primo de p es
-- maximal. Por ejemplo,
--    esMaximalHuecoPrimo  7  ==  True
--    esMaximalHuecoPrimo 11  ==  False
esMaximalHuecoPrimo :: Integer -> Bool
esMaximalHuecoPrimo p =
  and [huecoPrimo n < h | n <- takeWhile (< p) primes]
  where h = huecoPrimo p

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

primosYhuecosMaximales2 :: [(Integer,Integer)]
primosYhuecosMaximales2 = aux primosYhuecos
  where aux ((x,y):ps) = (x,y) : aux (dropWhile (\(_,b) -> b <= y) ps)

-- primosYhuecos es la lista de los números primos junto son sus
-- huecos. Por ejemplo, 
--    λ> take 10 primosYhuecos
--    [(2,1),(3,2),(5,2),(7,4),(11,2),(13,4),(17,2),(19,4),(23,6),(29,2)]
primosYhuecos :: [(Integer,Integer)]
primosYhuecos =
  [(x,y-x) | (x,y) <- zip primes (tail primes)]

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

primosYhuecosMaximales3 :: [(Integer,Integer)]
primosYhuecosMaximales3 = aux 0 primes
  where aux n (x:y:zs) | y-x > n   = (x,y-x) : aux (y-x) (y:zs)
                       | otherwise = aux n (y:zs)

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

-- La propiedad es
prop_primosYhuecosMaximales :: NonNegative Int -> Bool
prop_primosYhuecosMaximales (NonNegative n) =
  all (== primosYhuecosMaximales1 !! n)
      [primosYhuecosMaximales2 !! n,
       primosYhuecosMaximales3 !! n]

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

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

-- La comparación es
--    λ> primosYhuecosMaximales1 !! 10
--    (9551,36)
--    (2.63 secs, 7,400,316,112 bytes)
--    λ> primosYhuecosMaximales2 !! 10
--    (9551,36)
--    (0.01 secs, 7,060,744 bytes)
--    λ> primosYhuecosMaximales3 !! 10
--    (9551,36)
--    (0.01 secs, 4,000,368 bytes)
--    
--    λ> primosYhuecosMaximales2 !! 22
--    (17051707,180)
--    (7.90 secs, 17,275,407,712 bytes)
--    λ> primosYhuecosMaximales3 !! 22
--    (17051707,180)
--    (3.78 secs, 8,808,779,096 bytes)

import Data.Numbers.Primes (primes)

import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

-- 1ª solución

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

primosYhuecosMaximales1 :: [(Integer,Integer)]

primosYhuecosMaximales1 =

[(p,huecoPrimo p) | p <- primes, esMaximalHuecoPrimo p]

-- (siguientePrimo x) es el menor primo mayor que x. Por ejemplo,

-- siguientePrimo 7 == 11

-- siguientePrimo 8 == 11

siguientePrimo :: Integer -> Integer

siguientePrimo p =

head (dropWhile (<= p) primes)

-- (huecoPrimo p) es la distancia del primo p hasta el siguiente

-- primo. Por ejemplo,

-- huecoPrimo 7 == 4

huecoPrimo :: Integer -> Integer

huecoPrimo p = siguientePrimo p - p

-- (esMaximalHuecoPrimo p) se verifica si el hueco primo de p es

-- maximal. Por ejemplo,

-- esMaximalHuecoPrimo 7 == True

-- esMaximalHuecoPrimo 11 == False

esMaximalHuecoPrimo :: Integer -> Bool

esMaximalHuecoPrimo p =

and [huecoPrimo n < h | n <- takeWhile (< p) primes]

where h = huecoPrimo p

-- 2ª solución

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

primosYhuecosMaximales2 :: [(Integer,Integer)]

primosYhuecosMaximales2 = aux primosYhuecos

where aux ((x,y):ps) = (x,y) : aux (dropWhile (\(_,b) -> b <= y) ps)

-- primosYhuecos es la lista de los números primos junto son sus

-- huecos. Por ejemplo,

-- λ> take 10 primosYhuecos

-- [(2,1),(3,2),(5,2),(7,4),(11,2),(13,4),(17,2),(19,4),(23,6),(29,2)]

primosYhuecos :: [(Integer,Integer)]

primosYhuecos =

[(x,y-x) | (x,y) <- zip primes (tail primes)]

-- 3ª solución

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

primosYhuecosMaximales3 :: [(Integer,Integer)]

primosYhuecosMaximales3 = aux 0 primes

where aux n (x:y:zs) | y-x > n = (x,y-x) : aux (y-x) (y:zs)

| otherwise = aux n (y:zs)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosYhuecosMaximales :: NonNegative Int -> Bool

prop_primosYhuecosMaximales (NonNegative n) =

all (== primosYhuecosMaximales1 !! n)

[primosYhuecosMaximales2 !! n,

primosYhuecosMaximales3 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=12}) prop_primosYhuecosMaximales

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosYhuecosMaximales1 !! 10

-- (9551,36)

-- (2.63 secs, 7,400,316,112 bytes)

-- λ> primosYhuecosMaximales2 !! 10

-- (9551,36)

-- (0.01 secs, 7,060,744 bytes)

-- λ> primosYhuecosMaximales3 !! 10

-- (9551,36)

-- (0.01 secs, 4,000,368 bytes)

-- λ> primosYhuecosMaximales2 !! 22

-- (17051707,180)

-- (7.90 secs, 17,275,407,712 bytes)

-- λ> primosYhuecosMaximales3 !! 22

-- (17051707,180)

-- (3.78 secs, 8,808,779,096 bytes)

Referencias

Basado en el ejercicio Maximal prime gaps de
Programming Praxis.

Otras referencias

C. Caldwell The gaps between primes.
J.K. Andersen Maximal prime gaps.
N.J.A. Sloane Sequence A002386 en OEIS.
N.J.A. Sloane Sequence A005250 en OEIS.
E.W. Weisstein Prime gaps en MathWorld.

El código se encuentra en GitHub.

La función indicatriz de Euler

PorJosé A. Alonso 22-julio-202217-julio-2022

La indicatriz de Euler (también función φ de Euler) es una función importante en teoría de números. Si n es un entero positivo, entonces φ(n) se define como el número de enteros positivos menores o iguales a n y coprimos con n. Por ejemplo, φ(36) = 12 ya que los números menores o iguales a 36 y coprimos con 36 son doce: 1, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, y 35.

Definir la función

   phi :: Integer -> Integer

1	phi :: Integer -> Integer

tal que (phi n) es igual a φ(n). Por ejemplo,

   phi 36                          ==  12
   map phi [10..20]                ==  [4,10,4,12,6,8,8,16,6,18,8]
   phi (3^10^5) `mod` (10^9)       ==  681333334
   length (show (phi (10^(10^5)))) == 100000

phi 36 == 12

map phi [10..20] == [4,10,4,12,6,8,8,16,6,18,8]

phi (3^10^5) `mod` (10^9) == 681333334

length (show (phi (10^(10^5)))) == 100000

Comprobar con QuickCheck que, para todo n > 0, φ(10^n) tiene n dígitos.

Soluciones

import Data.List (genericLength, group)
import Data.Numbers.Primes (primeFactors)
import Math.NumberTheory.ArithmeticFunctions (totient)
import Test.QuickCheck (Positive (Positive), quickCheck)

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

phi1 :: Integer -> Integer
phi1 n = genericLength [x | x <- [1..n], gcd x n == 1] 

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

phi2 :: Integer -> Integer
phi2 n = product [(p-1)*p^(e-1) | (p,e) <- factorizacion n] 
    
factorizacion :: Integer -> [(Integer,Integer)]
factorizacion n =
  [(head xs,genericLength xs) | xs <- group (primeFactors n)]

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

phi3 :: Integer -> Integer
phi3 n = 
  product [(x-1) * product xs | (x:xs) <- group (primeFactors n)]

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

phi4 :: Integer -> Integer
phi4 = totient 

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

-- La propiedad es
prop_phi :: Positive Integer -> Bool 
prop_phi (Positive n) =
  all (== phi1 n)
      [phi2 n,
       phi3 n,
       phi4 n]

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

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

-- La comparación es
--    λ> phi1 (2*10^6)
--    800000
--    (2.49 secs, 2,117,853,856 bytes)
--    λ> phi2 (2*10^6)
--    800000
--    (0.02 secs, 565,664 bytes)
--    
--    λ> length (show (phi2 (10^100000)))
--    100000
--    (2.80 secs, 5,110,043,208 bytes)
--    λ> length (show (phi3 (10^100000)))
--    100000
--    (4.81 secs, 7,249,353,896 bytes)
--    λ> length (show (phi4 (10^100000)))
--    100000
--    (0.78 secs, 1,467,573,768 bytes)

-- Verificación de la propiedad
-- ============================

-- La propiedad es
prop_phi2 :: Positive Integer -> Bool
prop_phi2 (Positive n) =
  genericLength (show (phi4 (10^n))) == n

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

import Data.List (genericLength, group)

import Data.Numbers.Primes (primeFactors)

import Math.NumberTheory.ArithmeticFunctions (totient)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

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

phi1 :: Integer -> Integer

phi1 n = genericLength [x | x <- [1..n], gcd x n == 1]

-- 2ª solución

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

phi2 :: Integer -> Integer

phi2 n = product [(p-1)*p^(e-1) | (p,e) <- factorizacion n]

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

factorizacion n =

[(head xs,genericLength xs) | xs <- group (primeFactors n)]

-- 3ª solución

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

phi3 :: Integer -> Integer

phi3 n =

product [(x-1) * product xs | (x:xs) <- group (primeFactors n)]

-- 4ª solución

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

phi4 :: Integer -> Integer

phi4 = totient

-- Comprobación de equivalencia

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

-- La propiedad es

prop_phi :: Positive Integer -> Bool

prop_phi (Positive n) =

all (== phi1 n)

[phi2 n,

phi3 n,

phi4 n]

-- La comprobación es

-- λ> quickCheck prop_phi

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> phi1 (2*10^6)

-- 800000

-- (2.49 secs, 2,117,853,856 bytes)

-- λ> phi2 (2*10^6)

-- 800000

-- (0.02 secs, 565,664 bytes)

-- λ> length (show (phi2 (10^100000)))

-- 100000

-- (2.80 secs, 5,110,043,208 bytes)

-- λ> length (show (phi3 (10^100000)))

-- 100000

-- (4.81 secs, 7,249,353,896 bytes)

-- λ> length (show (phi4 (10^100000)))

-- 100000

-- (0.78 secs, 1,467,573,768 bytes)

-- Verificación de la propiedad

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

-- La propiedad es

prop_phi2 :: Positive Integer -> Bool

prop_phi2 (Positive n) =

genericLength (show (phi4 (10^n))) == n

-- La comprobación es

-- λ> quickCheck prop_phi2

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Sucesión de suma de cuadrados de los dígitos

PorJosé A. Alonso 21-julio-202216-julio-2022

Definir la función

   sucSumaCuadradosDigitos :: Integer -> [Integer]

1	sucSumaCuadradosDigitos :: Integer -> [Integer]

tal que (sucSumaCuadradosDigitos n) es la sucesión cuyo primer término es n y los restantes se obtienen sumando los cuadrados de los dígitos de su término anterior. Por ejemplo,

   λ> take 20 (sucSumaCuadradosDigitos1 2000)
   [2000,4,16,37,58,89,145,42,20,4,16,37,58,89,145,42,20,4,16,37]
   λ> take 20 (sucSumaCuadradosDigitos 1976)
   [1976,167,86,100,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
   λ> sucSumaCuadradosDigitos 2000 !! (10^9)
   20

λ> take 20 (sucSumaCuadradosDigitos1 2000)

[2000,4,16,37,58,89,145,42,20,4,16,37,58,89,145,42,20,4,16,37]

λ> take 20 (sucSumaCuadradosDigitos 1976)

[1976,167,86,100,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]

λ> sucSumaCuadradosDigitos 2000 !! (10^9)

Soluciones

import Test.QuickCheck (Positive (Positive), NonNegative (NonNegative), quickCheck)

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

sucSumaCuadradosDigitos1 :: Integer -> [Integer]
sucSumaCuadradosDigitos1 n =
  n : [sumaCuadradosDigitos x | x <- sucSumaCuadradosDigitos1 n]

-- (sumaCuadradosDigitos n) es la suma de los cuadrados de los dígitos
-- de n. Por ejemplo, 
--    sumaCuadradosDigitos 2016  ==  41
sumaCuadradosDigitos :: Integer -> Integer
sumaCuadradosDigitos n = sum (map (^2) (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 [d] | d <- show n]

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

sucSumaCuadradosDigitos2 :: Integer -> [Integer]
sucSumaCuadradosDigitos2 = iterate sumaCuadradosDigitos

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

-- A partir de los cálculos con las definiciones anteriores, se observa
-- que para todo n (sucSumaCuadradosDigitos n) tiene una parte pura y
-- otra periódica. Por ejemplo, para n = 2016, 
--    λ> take 20 (sucSumaCuadradosDigitos 2016)
--    [2016,41,17,50,25,29,85,89,145,42,20,4,16,37,58,89,145,42,20,4]
-- la parte pura es
--    [2016,41,17,50,25,29,85]
-- y la parte periódica es
--    [89,145,42,20,4,16,37,58])

sucSumaCuadradosDigitos3 :: Integer -> [Integer]
sucSumaCuadradosDigitos3 n = xs ++ cycle ys
  where (xs,ys) = sucCompactaSumaCuadradosDigitos n

-- (sucCompactaSumaCuadradosDigitos n) es el par formado por la parte
-- pura y la periódica de (sucSumaCuadradosDigitos n). Por ejemplo, 
--    λ> sucCompactaSumaCuadradosDigitos 2016
--    ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58])
--    λ> sucCompactaSumaCuadradosDigitos 1976
--    ([1976,167,86,100],[1])
sucCompactaSumaCuadradosDigitos :: Integer -> ([Integer],[Integer])
sucCompactaSumaCuadradosDigitos = 
  partePuraPeriodica . sucSumaCuadradosDigitos1

-- (partePuraPeriodica xs) es el par formado por la parte pura y la
-- periódica de xs. Por ejemplo,
--    λ> partePuraPeriodica (sucSumaCuadradosDigitos 2016)
--    ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58])
--    λ> partePuraPeriodica (sucSumaCuadradosDigitos 1976)
--    ([1976,167,86,100],[1])
partePuraPeriodica :: [Integer] -> ([Integer],[Integer])
partePuraPeriodica = aux [] 
  where aux as (b:bs) | b `elem` as = span (/=b) (reverse as)
                      | otherwise = aux (b:as) bs

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

sucSumaCuadradosDigitos4 :: Integer -> [Integer]
sucSumaCuadradosDigitos4 1  = repeat 1
sucSumaCuadradosDigitos4 89 = cycle [89,145,42,20,4,16,37,58]
sucSumaCuadradosDigitos4 n  =
  n : sucSumaCuadradosDigitos4 (sumaCuadradosDigitos n)

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

-- La propiedad es
prop_sucSumaCuadradosDigitos :: Positive Integer -> NonNegative Int -> Bool
prop_sucSumaCuadradosDigitos (Positive n) (NonNegative k) =
  all (== sucSumaCuadradosDigitos1 n !! k)
      [sucSumaCuadradosDigitos2 n !! k,
       sucSumaCuadradosDigitos3 n !! k,
       sucSumaCuadradosDigitos4 n !! k]

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

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

-- La comparación es
--    λ> sucSumaCuadradosDigitos1 2000 !! (10^4)
--    20
--    (6.96 secs, 8,049,886,312 bytes)
--    λ> sucSumaCuadradosDigitos2 2000 !! (10^4)
--    20
--    (0.08 secs, 91,024,688 bytes)
--    λ> sucSumaCuadradosDigitos3 2000 !! (10^4)
--    20
--    (0.01 secs, 995,560 bytes)
--    λ> sucSumaCuadradosDigitos4 2000 !! (10^4)
--    20
--    (0.02 secs, 587,040 bytes)
--    
--    λ> sucSumaCuadradosDigitos2 2000 !! (3*10^5)
--    20
--    (1.96 secs, 2,715,501,416 bytes)
--    λ> sucSumaCuadradosDigitos3 2000 !! (3*10^5)
--    20
--    (0.02 secs, 995,872 bytes)
--    λ> sucSumaCuadradosDigitos4 2000 !! (3*10^5)
--    20
--    (0.02 secs, 587,352 bytes)
--    
--    λ> sucSumaCuadradosDigitos3 2000 !! (10^9)
--    20
--    (2.85 secs, 996,016 bytes)
--    λ> sucSumaCuadradosDigitos4 2000 !! (10^9)
--    20
--    (2.54 secs, 587,496 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

import Test.QuickCheck (Positive (Positive), NonNegative (NonNegative), quickCheck)

-- 1ª solución

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

sucSumaCuadradosDigitos1 :: Integer -> [Integer]

sucSumaCuadradosDigitos1 n =

n : [sumaCuadradosDigitos x | x <- sucSumaCuadradosDigitos1 n]

-- (sumaCuadradosDigitos n) es la suma de los cuadrados de los dígitos

-- de n. Por ejemplo,

-- sumaCuadradosDigitos 2016 == 41

sumaCuadradosDigitos :: Integer -> Integer

sumaCuadradosDigitos n = sum (map (^2) (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 [d] | d <- show n]

-- 2ª solución

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

sucSumaCuadradosDigitos2 :: Integer -> [Integer]

sucSumaCuadradosDigitos2 = iterate sumaCuadradosDigitos

-- 3ª solución

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

-- A partir de los cálculos con las definiciones anteriores, se observa

-- que para todo n (sucSumaCuadradosDigitos n) tiene una parte pura y

-- otra periódica. Por ejemplo, para n = 2016,

-- λ> take 20 (sucSumaCuadradosDigitos 2016)

-- [2016,41,17,50,25,29,85,89,145,42,20,4,16,37,58,89,145,42,20,4]

-- la parte pura es

-- [2016,41,17,50,25,29,85]

-- y la parte periódica es

-- [89,145,42,20,4,16,37,58])

sucSumaCuadradosDigitos3 :: Integer -> [Integer]

sucSumaCuadradosDigitos3 n = xs ++ cycle ys

where (xs,ys) = sucCompactaSumaCuadradosDigitos n

-- (sucCompactaSumaCuadradosDigitos n) es el par formado por la parte

-- pura y la periódica de (sucSumaCuadradosDigitos n). Por ejemplo,

-- λ> sucCompactaSumaCuadradosDigitos 2016

-- ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58])

-- λ> sucCompactaSumaCuadradosDigitos 1976

-- ([1976,167,86,100],[1])

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

sucCompactaSumaCuadradosDigitos =

partePuraPeriodica . sucSumaCuadradosDigitos1

-- (partePuraPeriodica xs) es el par formado por la parte pura y la

-- periódica de xs. Por ejemplo,

-- λ> partePuraPeriodica (sucSumaCuadradosDigitos 2016)

-- ([2016,41,17,50,25,29,85],[89,145,42,20,4,16,37,58])

-- λ> partePuraPeriodica (sucSumaCuadradosDigitos 1976)

-- ([1976,167,86,100],[1])

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

partePuraPeriodica = aux []

where aux as (b:bs) | b `elem` as = span (/=b) (reverse as)

| otherwise = aux (b:as) bs

-- 4ª solución

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

sucSumaCuadradosDigitos4 :: Integer -> [Integer]

sucSumaCuadradosDigitos4 1 = repeat 1

sucSumaCuadradosDigitos4 89 = cycle [89,145,42,20,4,16,37,58]

sucSumaCuadradosDigitos4 n =

n : sucSumaCuadradosDigitos4 (sumaCuadradosDigitos n)

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sucSumaCuadradosDigitos :: Positive Integer -> NonNegative Int -> Bool

prop_sucSumaCuadradosDigitos (Positive n) (NonNegative k) =

all (== sucSumaCuadradosDigitos1 n !! k)

[sucSumaCuadradosDigitos2 n !! k,

sucSumaCuadradosDigitos3 n !! k,

sucSumaCuadradosDigitos4 n !! k]

-- La comprobación es

-- λ> quickCheck prop_sucSumaCuadradosDigitos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> sucSumaCuadradosDigitos1 2000 !! (10^4)

-- 20

-- (6.96 secs, 8,049,886,312 bytes)

-- λ> sucSumaCuadradosDigitos2 2000 !! (10^4)

-- 20

-- (0.08 secs, 91,024,688 bytes)

-- λ> sucSumaCuadradosDigitos3 2000 !! (10^4)

-- 20

-- (0.01 secs, 995,560 bytes)

-- λ> sucSumaCuadradosDigitos4 2000 !! (10^4)

-- 20

-- (0.02 secs, 587,040 bytes)

-- λ> sucSumaCuadradosDigitos2 2000 !! (3*10^5)

-- 20

-- (1.96 secs, 2,715,501,416 bytes)

-- λ> sucSumaCuadradosDigitos3 2000 !! (3*10^5)

-- 20

-- (0.02 secs, 995,872 bytes)

-- λ> sucSumaCuadradosDigitos4 2000 !! (3*10^5)

-- 20

-- (0.02 secs, 587,352 bytes)

-- λ> sucSumaCuadradosDigitos3 2000 !! (10^9)

-- 20

-- (2.85 secs, 996,016 bytes)

-- λ> sucSumaCuadradosDigitos4 2000 !! (10^9)

-- 20

-- (2.54 secs, 587,496 bytes)

El código se encuentra en GitHub.

Potencias perfectas

PorJosé A. Alonso 20-julio-202213-julio-2022

Un número natural n es una potencia perfecta si existen dos números naturales m > 1 y k > 1 tales que n = m^k. Las primeras potencias perfectas son

   4 = 2², 8 = 2³, 9 = 3², 16 = 2⁴, 25 = 5², 27 = 3³, 32 = 2⁵, 
   36 = 6², 49 = 7², 64 = 2⁶, ...

1 2	4 = 2², 8 = 2³, 9 = 3², 16 = 2⁴, 25 = 5², 27 = 3³, 32 = 2⁵, 36 = 6², 49 = 7², 64 = 2⁶, ...

Definir la sucesión

   potenciasPerfectas :: [Integer]

1	potenciasPerfectas :: [Integer]

cuyos términos son las potencias perfectas. Por ejemplo,

   take 10 potenciasPerfectas  ==  [4,8,9,16,25,27,32,36,49,64]
   potenciasPerfectas !! 3000  ==  7778521

1 2	take 10 potenciasPerfectas == [4,8,9,16,25,27,32,36,49,64] potenciasPerfectas !! 3000 == 7778521

Definir el procedimiento

   grafica :: Int -> IO ()

1	grafica :: Int -> IO ()

tal que (grafica n) es la representación gráfica de las n primeras potencias perfectas. Por ejemplo, para (grafica 30) dibuja

Soluciones

import Data.List (group)
import Data.Numbers.Primes (primeFactors)
import Graphics.Gnuplot.Simple (Attribute (Key, PNG), plotList)
import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

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

potenciasPerfectas1 :: [Integer]
potenciasPerfectas1 = filter esPotenciaPerfecta [4..]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta 36  ==  True
--    esPotenciaPerfecta 72  ==  False
esPotenciaPerfecta :: Integer -> Bool
esPotenciaPerfecta = not . null. potenciasPerfectasDe 

-- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que 
-- x = a^b. Por ejemplo,
--    potenciasPerfectasDe 64  ==  [(2,6),(4,3),(8,2)]
--    potenciasPerfectasDe 72  ==  []
potenciasPerfectasDe :: Integer -> [(Integer,Integer)]
potenciasPerfectasDe n = 
  [(m,k) | m <- takeWhile (\x -> x*x <= n) [2..]
         , k <- takeWhile (\x -> m^x <= n) [2..]
         , m^k == n]

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

potenciasPerfectas2 :: [Integer]
potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x]

-- (esPotenciaPerfecta2 x) se verifica si x es una potencia perfecta. Por
-- ejemplo, 
--    esPotenciaPerfecta2 36  ==  True
--    esPotenciaPerfecta2 72  ==  False
esPotenciaPerfecta2 :: Integer -> Bool
esPotenciaPerfecta2 x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de l factorización prima
-- de x. Por ejemplos,
--    exponentes 36  ==  [2,2]
--    exponentes 72  ==  [3,2]
exponentes :: Integer -> [Int]
exponentes x = [length ys | ys <- group (primeFactors x)] 

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,
--    mcd [4,6,10]  ==  2
--    mcd [4,5,10]  ==  1
mcd :: [Int] -> Int
mcd = foldl1 gcd

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

potenciasPerfectas3 :: [Integer]
potenciasPerfectas3 = mezclaTodas potencias

-- potencias es la lista las listas de potencias de todos los números
-- mayores que 1 con exponentes mayores que 1. Por ejemplo,
--    λ> map (take 3) (take 4 potencias)
--    [[4,8,16],[9,27,81],[16,64,256],[25,125,625]]
potencias:: [[Integer]]
potencias = [[n^k | k <- [2..]] | n <- [2..]]

-- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las
-- listas ordenadas xss. Por ejemplo,
--    take 7 (mezclaTodas potencias)  ==  [4,8,9,16,25,27,32]
mezclaTodas :: Ord a => [[a]] -> [a]
mezclaTodas = foldr1 xmezcla
  where xmezcla (x:xs) ys = x : mezcla xs ys

-- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las
-- listas ordenadas xs e ys. Por ejemplo,
--    take 7 (mezcla [2,5..] [4,6..])  ==  [2,4,5,6,8,10,11]
mezcla :: Ord a => [a] -> [a] -> [a]
mezcla (x:xs) (y:ys) | x < y  = x : mezcla xs (y:ys)
                     | x == y = x : mezcla xs ys
                     | x > y  = y : mezcla (x:xs) ys

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

-- La propiedad es
prop_potenciasPerfectas :: NonNegative Int -> Bool
prop_potenciasPerfectas (NonNegative n) =
  all (== potenciasPerfectas1 !! n)
      [potenciasPerfectas2 !! n,
       potenciasPerfectas3 !! n]

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

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

-- La comparación es
--    λ> potenciasPerfectas1 !! 200
--    28224
--    (10.56 secs, 8,434,647,368 bytes)
--    λ> potenciasPerfectas2 !! 200
--    28224
--    (0.36 secs, 825,040,416 bytes)
--    λ> potenciasPerfectas3 !! 200
--    28224
--    (0.05 secs, 7,474,280 bytes)
--    
--    λ> potenciasPerfectas2 !! 500
--    191844
--    (4.16 secs, 9,899,367,112 bytes)
--    λ> potenciasPerfectas3 !! 500
--    191844
--    (0.09 secs, 51,275,464 bytes)

-- En lo que sigue se usa la 3ª solución
potenciasPerfectas :: [Integer]
potenciasPerfectas = potenciasPerfectas3

-- Representación gráfica
-- ======================

grafica :: Int -> IO ()
grafica n = 
  plotList [ Key Nothing
           , PNG "Potencias_perfectas.png"
           ]
           (take n potenciasPerfectas)

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

import Data.List (group)

import Data.Numbers.Primes (primeFactors)

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

import Test.QuickCheck (NonNegative (NonNegative), quickCheck)

-- 1ª solución

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

potenciasPerfectas1 :: [Integer]

potenciasPerfectas1 = filter esPotenciaPerfecta [4..]

-- (esPotenciaPerfecta x) se verifica si x es una potencia perfecta. Por

-- ejemplo,

-- esPotenciaPerfecta 36 == True

-- esPotenciaPerfecta 72 == False

esPotenciaPerfecta :: Integer -> Bool

esPotenciaPerfecta = not . null. potenciasPerfectasDe

-- (potenciasPerfectasDe x) es la lista de pares (a,b) tales que

-- x = a^b. Por ejemplo,

-- potenciasPerfectasDe 64 == [(2,6),(4,3),(8,2)]

-- potenciasPerfectasDe 72 == []

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

potenciasPerfectasDe n =

[(m,k) | m <- takeWhile (\x -> x*x <= n) [2..]

, k <- takeWhile (\x -> m^x <= n) [2..]

, m^k == n]

-- 2ª solución

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

potenciasPerfectas2 :: [Integer]

potenciasPerfectas2 = [x | x <- [4..], esPotenciaPerfecta2 x]

-- (esPotenciaPerfecta2 x) se verifica si x es una potencia perfecta. Por

-- ejemplo,

-- esPotenciaPerfecta2 36 == True

-- esPotenciaPerfecta2 72 == False

esPotenciaPerfecta2 :: Integer -> Bool

esPotenciaPerfecta2 x = mcd (exponentes x) > 1

-- (exponentes x) es la lista de los exponentes de l factorización prima

-- de x. Por ejemplos,

-- exponentes 36 == [2,2]

-- exponentes 72 == [3,2]

exponentes :: Integer -> [Int]

exponentes x = [length ys | ys <- group (primeFactors x)]

-- (mcd xs) es el máximo común divisor de la lista xs. Por ejemplo,

-- mcd [4,6,10] == 2

-- mcd [4,5,10] == 1

mcd :: [Int] -> Int

mcd = foldl1 gcd

-- 3ª solución

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

potenciasPerfectas3 :: [Integer]

potenciasPerfectas3 = mezclaTodas potencias

-- potencias es la lista las listas de potencias de todos los números

-- mayores que 1 con exponentes mayores que 1. Por ejemplo,

-- λ> map (take 3) (take 4 potencias)

-- [[4,8,16],[9,27,81],[16,64,256],[25,125,625]]

potencias:: [[Integer]]

potencias = [[n^k | k <- [2..]] | n <- [2..]]

-- (mezclaTodas xss) es la mezcla ordenada sin repeticiones de las

-- listas ordenadas xss. Por ejemplo,

-- take 7 (mezclaTodas potencias) == [4,8,9,16,25,27,32]

mezclaTodas :: Ord a => [[a]] -> [a]

mezclaTodas = foldr1 xmezcla

where xmezcla (x:xs) ys = x : mezcla xs ys

-- (mezcla xs ys) es la mezcla ordenada sin repeticiones de las

-- listas ordenadas xs e ys. Por ejemplo,

-- take 7 (mezcla [2,5..] [4,6..]) == [2,4,5,6,8,10,11]

mezcla :: Ord a => [a] -> [a] -> [a]

mezcla (x:xs) (y:ys) | x < y = x : mezcla xs (y:ys)

| x == y = x : mezcla xs ys

| x > y = y : mezcla (x:xs) ys

-- Comprobación de equivalencia

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

-- La propiedad es

prop_potenciasPerfectas :: NonNegative Int -> Bool

prop_potenciasPerfectas (NonNegative n) =

all (== potenciasPerfectas1 !! n)

[potenciasPerfectas2 !! n,

potenciasPerfectas3 !! n]

-- La comprobación es

-- λ> quickCheck prop_potenciasPerfectas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> potenciasPerfectas1 !! 200

-- 28224

-- (10.56 secs, 8,434,647,368 bytes)

-- λ> potenciasPerfectas2 !! 200

-- 28224

-- (0.36 secs, 825,040,416 bytes)

-- λ> potenciasPerfectas3 !! 200

-- 28224

-- (0.05 secs, 7,474,280 bytes)

-- λ> potenciasPerfectas2 !! 500

-- 191844

-- (4.16 secs, 9,899,367,112 bytes)

-- λ> potenciasPerfectas3 !! 500

-- 191844

-- (0.09 secs, 51,275,464 bytes)

-- En lo que sigue se usa la 3ª solución

potenciasPerfectas :: [Integer]

potenciasPerfectas = potenciasPerfectas3

-- Representación gráfica

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

grafica :: Int -> IO ()

grafica n =

plotList [ Key Nothing

, PNG "Potencias_perfectas.png"

]

(take n potenciasPerfectas)

El código se encuentra en GitHub.

Sumas alternas de factoriales

PorJosé A. Alonso 19-julio-202212-julio-2022

Las primeras sumas alternas de los factoriales son números primos; en efecto,

   3! - 2! + 1! = 5
   4! - 3! + 2! - 1! = 19
   5! - 4! + 3! - 2! + 1! = 101
   6! - 5! + 4! - 3! + 2! - 1! = 619
   7! - 6! + 5! - 4! + 3! - 2! + 1! = 4421
   8! - 7! + 6! - 5! + 4! - 3! + 2! - 1! = 35899

3! - 2! + 1! = 5

4! - 3! + 2! - 1! = 19

5! - 4! + 3! - 2! + 1! = 101

6! - 5! + 4! - 3! + 2! - 1! = 619

7! - 6! + 5! - 4! + 3! - 2! + 1! = 4421

8! - 7! + 6! - 5! + 4! - 3! + 2! - 1! = 35899

son primos, pero

   9! - 8! + 7! - 6! + 5! - 4! + 3! - 2! + 1! = 326981

1	9! - 8! + 7! - 6! + 5! - 4! + 3! - 2! + 1! = 326981

no es primo.

Definir las funciones

   sumaAlterna         :: Integer -> Integer
   sumasAlternas       :: [Integer]
   conSumaAlternaPrima :: [Integer]

sumaAlterna :: Integer -> Integer

sumasAlternas :: [Integer]

conSumaAlternaPrima :: [Integer]

tales que

(sumaAlterna n) es la suma alterna de los factoriales desde n hasta 1. Por ejemplo,

     sumaAlterna 3  ==  5
     sumaAlterna 4  ==  19
     sumaAlterna 5  ==  101
     sumaAlterna 6  ==  619
     sumaAlterna 7  ==  4421
     sumaAlterna 8  ==  35899
     sumaAlterna 9  ==  326981
     sumaAlterna (5*10^4) `mod` (10^6) == 577019

sumaAlterna 3 == 5

sumaAlterna 4 == 19

sumaAlterna 5 == 101

sumaAlterna 6 == 619

sumaAlterna 7 == 4421

sumaAlterna 8 == 35899

sumaAlterna 9 == 326981

sumaAlterna (5*10^4) `mod` (10^6) == 577019

sumasAlternas es la sucesión de las sumas alternas de factoriales. Por ejemplo,

     λ> take 10 sumasAlternas1
     [0,1,1,5,19,101,619,4421,35899,326981]

1 2	λ> take 10 sumasAlternas1 [0,1,1,5,19,101,619,4421,35899,326981]

conSumaAlternaPrima es la sucesión de los números cuya suma alterna de factoriales es prima. Por ejemplo,

     λ> take 8 conSumaAlternaPrima
     [3,4,5,6,7,8,10,15]

1 2	λ> take 8 conSumaAlternaPrima [3,4,5,6,7,8,10,15]

Soluciones

import Data.List (genericTake)
import Data.Numbers.Primes (isPrime)
import Test.QuickCheck

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

sumaAlterna1 :: Integer -> Integer
sumaAlterna1 1 = 1
sumaAlterna1 n = factorial n - sumaAlterna1 (n-1)

factorial :: Integer -> Integer
factorial n = product [1..n]

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

sumaAlterna2 :: Integer -> Integer
sumaAlterna2 n =
  sum (genericTake n (zipWith (*) signos (tail factoriales)))
  where
    signos | odd n     = concat (repeat [1,-1])
           | otherwise = concat (repeat [-1,1])

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

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

sumaAlterna3 :: Integer -> Integer
sumaAlterna3 n = 
  sum (genericTake n (zipWith (*) signos (tail factoriales)))
  where signos | odd n     = cycle [1,-1]
               | otherwise = cycle [-1,1]

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

sumaAlterna4 :: Integer -> Integer
sumaAlterna4 n =
  foldl (flip (-)) 0 (scanl1 (*) [1..n])

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

-- La propiedad es
prop_sumaAlterna :: Positive Integer -> Bool 
prop_sumaAlterna (Positive n) =
  all (== sumaAlterna1 n)
      [sumaAlterna2 n,
       sumaAlterna3 n,
       sumaAlterna4 n]

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

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

-- La comparación es
--    λ> sumaAlterna1 4000 `mod` (10^6)
--    577019
--    (6.21 secs, 16,154,113,192 bytes)
--    λ> sumaAlterna2 4000 `mod` (10^6)
--    577019
--    (0.01 secs, 24,844,664 bytes)
--    
--    λ> sumaAlterna2 (5*10^4) `mod` (10^6)
--    577019
--    (1.81 secs, 4,729,583,864 bytes)
--    λ> sumaAlterna3 (5*10^4) `mod` (10^6)
--    577019
--    (0.89 secs, 4,725,983,928 bytes)
--    λ> sumaAlterna4 (5*10^4) `mod` (10^6)
--    577019
--    (0.70 secs, 4,710,770,592 bytes)

-- En lo que sigue se usa la 3ª definición
sumaAlterna :: Integer -> Integer
sumaAlterna = sumaAlterna3

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

sumasAlternas1 :: [Integer]
sumasAlternas1 =
  map sumaAlterna [0..]

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

sumasAlternas2 :: [Integer]
sumasAlternas2 =
  0 : zipWith (-) (tail factoriales) sumasAlternas2

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

sumasAlternas3 :: [Integer]
sumasAlternas3 =
  scanl (flip (-)) 0 $ scanl1 (*) [1..]

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

-- La propiedad es
prop_sumasAlternas :: NonNegative Int -> Bool
prop_sumasAlternas (NonNegative n) =
  all (== sumasAlternas1 !! n)
      [sumasAlternas2 !! n,
       sumasAlternas3 !! n]
  
-- La comprobación es
--    λ> quickCheck prop_sumasAlternas
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> length (show (sumasAlternas1 !! (5*10^4)))
--    213237
--    (4.90 secs, 4,731,620,600 bytes)
--    λ> length (show (sumasAlternas2 !! (5*10^4)))
--    213237
--    (2.39 secs, 4,726,820,456 bytes)
--    λ> length (show (sumasAlternas3 !! (5*10^4)))
--    213237
--    (1.78 secs, 4,726,820,280 bytes)

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

conSumaAlternaPrima1 :: [Integer]
conSumaAlternaPrima1 =
  [n | n <- [0..], isPrime (sumaAlterna n)]

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

conSumaAlternaPrima2 :: [Integer]
conSumaAlternaPrima2 =
  [x | (x,y) <- zip [0..] sumasAlternas2, isPrime y]

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

conSumaAlternaPrima3 :: [Integer]
conSumaAlternaPrima3 =
  filter (isPrime . sumaAlterna) [0..]

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

-- La propiedad es
prop_conSumaAlternaPrima :: NonNegative Int -> Bool
prop_conSumaAlternaPrima (NonNegative n) =
  all (== conSumaAlternaPrima1 !! n)
      [conSumaAlternaPrima2 !! n,
       conSumaAlternaPrima3 !! n]

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=5}) prop_conSumaAlternaPrima
--    +++ 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

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

import Data.List (genericTake)

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck

-- 1ª definición de sumaAlterna

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

sumaAlterna1 :: Integer -> Integer

sumaAlterna1 1 = 1

sumaAlterna1 n = factorial n - sumaAlterna1 (n-1)

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª definición de sumaAlterna

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

sumaAlterna2 :: Integer -> Integer

sumaAlterna2 n =

sum (genericTake n (zipWith (*) signos (tail factoriales)))

where

signos | odd n = concat (repeat [1,-1])

| otherwise = concat (repeat [-1,1])

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

-- 3ª definición de sumaAlterna

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

sumaAlterna3 :: Integer -> Integer

sumaAlterna3 n =

sum (genericTake n (zipWith (*) signos (tail factoriales)))

where signos | odd n = cycle [1,-1]

| otherwise = cycle [-1,1]

-- 3ª definición de sumaAlterna

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

sumaAlterna4 :: Integer -> Integer

sumaAlterna4 n =

foldl (flip (-)) 0 (scanl1 (*) [1..n])

-- Comprobación de equivalencia de sumaAlterna

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

-- La propiedad es

prop_sumaAlterna :: Positive Integer -> Bool

prop_sumaAlterna (Positive n) =

all (== sumaAlterna1 n)

[sumaAlterna2 n,

sumaAlterna3 n,

sumaAlterna4 n]

-- La comprobación es

-- λ> quickCheck prop_sumaAlterna

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de sumaAlterna

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

-- La comparación es

-- λ> sumaAlterna1 4000 `mod` (10^6)

-- 577019

-- (6.21 secs, 16,154,113,192 bytes)

-- λ> sumaAlterna2 4000 `mod` (10^6)

-- 577019

-- (0.01 secs, 24,844,664 bytes)

-- λ> sumaAlterna2 (5*10^4) `mod` (10^6)

-- 577019

-- (1.81 secs, 4,729,583,864 bytes)

-- λ> sumaAlterna3 (5*10^4) `mod` (10^6)

-- 577019

-- (0.89 secs, 4,725,983,928 bytes)

-- λ> sumaAlterna4 (5*10^4) `mod` (10^6)

-- 577019

-- (0.70 secs, 4,710,770,592 bytes)

-- En lo que sigue se usa la 3ª definición

sumaAlterna :: Integer -> Integer

sumaAlterna = sumaAlterna3

-- 1ª definición de sumasAlternas

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

sumasAlternas1 :: [Integer]

sumasAlternas1 =

map sumaAlterna [0..]

-- 2ª definición de sumasAlternas

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

sumasAlternas2 :: [Integer]

sumasAlternas2 =

0 : zipWith (-) (tail factoriales) sumasAlternas2

-- 3ª definición de sumasAlternas

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

sumasAlternas3 :: [Integer]

sumasAlternas3 =

scanl (flip (-)) 0 $ scanl1 (*) [1..]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_sumasAlternas :: NonNegative Int -> Bool

prop_sumasAlternas (NonNegative n) =

all (== sumasAlternas1 !! n)

[sumasAlternas2 !! n,

sumasAlternas3 !! n]

-- La comprobación es

-- λ> quickCheck prop_sumasAlternas

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (show (sumasAlternas1 !! (5*10^4)))

-- 213237

-- (4.90 secs, 4,731,620,600 bytes)

-- λ> length (show (sumasAlternas2 !! (5*10^4)))

-- 213237

-- (2.39 secs, 4,726,820,456 bytes)

-- λ> length (show (sumasAlternas3 !! (5*10^4)))

-- 213237

-- (1.78 secs, 4,726,820,280 bytes)

-- 1ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima1 :: [Integer]

conSumaAlternaPrima1 =

[n | n <- [0..], isPrime (sumaAlterna n)]

-- 2ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima2 :: [Integer]

conSumaAlternaPrima2 =

[x | (x,y) <- zip [0..] sumasAlternas2, isPrime y]

-- 3ª definición de conSumaAlternaPrima

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

conSumaAlternaPrima3 :: [Integer]

conSumaAlternaPrima3 =

filter (isPrime . sumaAlterna) [0..]

-- Comprobación de equivalencia de conSumaAlternaPrima

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

-- La propiedad es

prop_conSumaAlternaPrima :: NonNegative Int -> Bool

prop_conSumaAlternaPrima (NonNegative n) =

all (== conSumaAlternaPrima1 !! n)

[conSumaAlternaPrima2 !! n,

conSumaAlternaPrima3 !! n]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

El código se encuentra en GitHub.

Primos con cubos

PorJosé A. Alonso 18-julio-202211-julio-2022

Un primo con cubo es un número primo p para el que existe algún entero positivo n tal que la expresión n^3 + n^2p es un cubo perfecto. Por ejemplo, 19 es un primo con cubo ya que 8^3 + 8^2×19 = 12^3.

Definir la sucesión

   primosConCubos :: [Integer]

1	primosConCubos :: [Integer]

tal que sus elementos son los primos con cubo. Por ejemplo,

   λ> take 6 primosConCubos
   [7,19,37,61,127,271]
   λ> length (takeWhile (< 1000000) primosConCubos)
   173

λ> take 6 primosConCubos

[7,19,37,61,127,271]

λ> length (takeWhile (< 1000000) primosConCubos)

173

Soluciones

import Data.Numbers.Primes (isPrime)
import Test.QuickCheck (NonNegative (NonNegative), maxSize, quickCheckWith, stdArgs) 

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

primosConCubos1 :: [Integer]
primosConCubos1 =
  [p | x <- [1..],
       n <- [1..x],
       (x^3 - n^3) `mod` (n^2) == 0,
       let p = (x^3 - n^3) `div` (n^2),
       isPrime p]

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

-- Para analizar la respuesta, en esta solución se calculan los pares
-- (p,n) tales que p es un primo con cubo y n es un número positivo tal
-- que n^3 + n^2p es un cubo

primosConCubos2' :: [(Integer,Integer)]
primosConCubos2' =
  [(p,n) | x <- [1..],
           n <- [1..x],
           (x^3 - n^3) `mod` (n^2) == 0,
           let p = (x^3 - n^3) `div` (n^2),
           isPrime p]

-- El cálculo es
--    λ> take 7 primosConCubos2
--    [(7,1),(19,8),(37,27),(61,64),(127,216),(271,729),(331,1000)]

-- Se observa que la sucesión de los segundos elementos [1,8,27,64,...]
-- es la de los cubos y que los primeros elementos se obtienen restando
-- los segundos elementos consecutivos; es decir,
--     7 =  8 -  1 = 2^3 - 1^3  
--    19 = 27 -  8 = 3^3 - 2^3
--    37 = 64 - 27 = 4^3 - 3^3
-- Continuando el patrón,
--     61 =  5^3 - 4^3 =  125 -   64
--     91 =  6^3 - 5^3 =  216 -  125
--    127 =  7^3 - 6^3 =  343 -  216
--    271 = 10^3 - 9^3 = 1000 -  729
--    331 = 11^3 -10^3 = 1331 - 1000
-- Por tanto, los primos con cubos son diferencias de dos cubos
-- consecutivos; es decir, coinciden con los números cubanos del
-- ejercicio anterior. A partir de la conjetura se obtienen las
-- siguientes definiciones

-- Basado en las anteriores observaciones se obtiene la siguiente
-- definición 
primosConCubos2 :: [Integer]
primosConCubos2 = 
  filter isPrime [(x+1)^3 - x^3 | x <- [1..]] 

-- 3ª definición
-- =============

primosConCubos3 :: [Integer]
primosConCubos3 =
  filter isPrime diferenciasCubosConsecutivos

diferenciasCubosConsecutivos :: [Integer]
diferenciasCubosConsecutivos =
  zipWith (-) (tail cubos) cubos 

cubos :: [Integer]
cubos = map (^3) [0..]

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

-- Simplificando la expresión (x+1)^3 - x^3 se obtiene 3*x^2+ 3*x + 1,
-- con lo que la 3ª definición se reduce a

primosConCubos4 :: [Integer]
primosConCubos4 =
  [p | x <- [1..],
       let p = 3*x^2+ 3*x + 1,
       isPrime p]

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

-- La propiedad es
prop_primosConCubos :: NonNegative Int -> Bool
prop_primosConCubos (NonNegative n) =
  all (== primosConCubos1 !! n)
      [primosConCubos2 !! n,
       primosConCubos3 !! n,
       primosConCubos4 !! n]
  
-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxSize=10}) prop_primosConCubos
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> primosConCubos1 !! 6
--    331
--    (1.15 secs, 1,105,612,336 bytes)
--    λ> primosConCubos1 !! 7
--    397
--    (1.96 secs, 1,909,369,592 bytes)
--    λ> primosConCubos2 !! 7
--    397
--    
--    (0.01 secs, 648,840 bytes)
--    λ> primosConCubos2 !! (10^3)
--    65580901
--    (0.53 secs, 1,726,837,688 bytes)
--    λ> primosConCubos3 !! (10^3)
--    65580901
--    (0.49 secs, 1,724,258,632 bytes)
--    λ> primosConCubos4 !! (10^3)
--    65580901
--    (0.47 secs, 1,724,833,992 bytes)

-- Demostración de la conjetura
-- ============================

-- Vamos a demostrar que los primos con cubos son diferencias de dos
-- cubos consecutivos.
--
-- Sea p un primo con cubo. Por la definición, existe un entero
-- positivo n tal que n³ + n²p es un cubo. 
--
-- Lema 1: Los números n y p son coprimos (es decir, mcd(n,p) = 1).
-- Dem.: En caso contrario, puesto que p es primo, existe un a tal que 
-- n = ap. Luego n³ + n²p = (a³+a²)p³ es un cubo y, por tanto,
-- a³+a² es un cubo lo que es imposible ya que el siguiente cubo de
-- a³ es a³+3a²+3a+1.
--
-- Lema 2: Los números n² y n+p son coprimos.
-- Dem.: Sea k = mcd(n^2,n+p). Por k divide n², luego k divide a n;
-- además, k divide a n+p y (usando el lema 1 y el ser p primo), se
-- tiene que k = 1.
--
-- Puesto que n³+n²p = n²(n+p) es un cubo, usando el lema 2, se tiene
-- que n² y n+p son cubos y, por serlo n², n también es un cubo. Es
-- decir, existen enteros positivos x e y tales que n = x³ y 
-- n+p = y³. Por tanto, p = y³-x³. Sea k = y-x. Se tiene que k = 1 ya
-- que 
--    p = y³-x³
--      = (n+k)³-n³
--      = 3k+3k²+k³
-- no es primo para k > 1.
--
-- Por consiguiente, p = (x+1)³-x³.

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

152

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck (NonNegative (NonNegative), maxSize, quickCheckWith, stdArgs)

-- 1ª solución

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

primosConCubos1 :: [Integer]

primosConCubos1 =

[p | x <- [1..],

n <- [1..x],

(x^3 - n^3) `mod` (n^2) == 0,

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

isPrime p]

-- 2ª solución

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

-- Para analizar la respuesta, en esta solución se calculan los pares

-- (p,n) tales que p es un primo con cubo y n es un número positivo tal

-- que n^3 + n^2p es un cubo

primosConCubos2' :: [(Integer,Integer)]

primosConCubos2' =

[(p,n) | x <- [1..],

n <- [1..x],

(x^3 - n^3) `mod` (n^2) == 0,

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

isPrime p]

-- El cálculo es

-- λ> take 7 primosConCubos2

-- [(7,1),(19,8),(37,27),(61,64),(127,216),(271,729),(331,1000)]

-- Se observa que la sucesión de los segundos elementos [1,8,27,64,...]

-- es la de los cubos y que los primeros elementos se obtienen restando

-- los segundos elementos consecutivos; es decir,

-- 7 = 8 - 1 = 2^3 - 1^3

-- 19 = 27 - 8 = 3^3 - 2^3

-- 37 = 64 - 27 = 4^3 - 3^3

-- Continuando el patrón,

-- 61 = 5^3 - 4^3 = 125 - 64

-- 91 = 6^3 - 5^3 = 216 - 125

-- 127 = 7^3 - 6^3 = 343 - 216

-- 271 = 10^3 - 9^3 = 1000 - 729

-- 331 = 11^3 -10^3 = 1331 - 1000

-- Por tanto, los primos con cubos son diferencias de dos cubos

-- consecutivos; es decir, coinciden con los números cubanos del

-- ejercicio anterior. A partir de la conjetura se obtienen las

-- siguientes definiciones

-- Basado en las anteriores observaciones se obtiene la siguiente

-- definición

primosConCubos2 :: [Integer]

primosConCubos2 =

filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

-- 3ª definición

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

primosConCubos3 :: [Integer]

primosConCubos3 =

filter isPrime diferenciasCubosConsecutivos

diferenciasCubosConsecutivos :: [Integer]

diferenciasCubosConsecutivos =

zipWith (-) (tail cubos) cubos

cubos :: [Integer]

cubos = map (^3) [0..]

-- 4ª solución

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

-- Simplificando la expresión (x+1)^3 - x^3 se obtiene 3*x^2+ 3*x + 1,

-- con lo que la 3ª definición se reduce a

primosConCubos4 :: [Integer]

primosConCubos4 =

[p | x <- [1..],

let p = 3*x^2+ 3*x + 1,

isPrime p]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_primosConCubos :: NonNegative Int -> Bool

prop_primosConCubos (NonNegative n) =

all (== primosConCubos1 !! n)

[primosConCubos2 !! n,

primosConCubos3 !! n,

primosConCubos4 !! n]

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxSize=10}) prop_primosConCubos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> primosConCubos1 !! 6

-- 331

-- (1.15 secs, 1,105,612,336 bytes)

-- λ> primosConCubos1 !! 7

-- 397

-- (1.96 secs, 1,909,369,592 bytes)

-- λ> primosConCubos2 !! 7

-- 397

-- (0.01 secs, 648,840 bytes)

-- λ> primosConCubos2 !! (10^3)

-- 65580901

-- (0.53 secs, 1,726,837,688 bytes)

-- λ> primosConCubos3 !! (10^3)

-- 65580901

-- (0.49 secs, 1,724,258,632 bytes)

-- λ> primosConCubos4 !! (10^3)

-- 65580901

-- (0.47 secs, 1,724,833,992 bytes)

-- Demostración de la conjetura

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

-- Vamos a demostrar que los primos con cubos son diferencias de dos

-- cubos consecutivos.

-- Sea p un primo con cubo. Por la definición, existe un entero

-- positivo n tal que n³ + n²p es un cubo.

-- Lema 1: Los números n y p son coprimos (es decir, mcd(n,p) = 1).

-- Dem.: En caso contrario, puesto que p es primo, existe un a tal que

-- n = ap. Luego n³ + n²p = (a³+a²)p³ es un cubo y, por tanto,

-- a³+a² es un cubo lo que es imposible ya que el siguiente cubo de

-- a³ es a³+3a²+3a+1.

-- Lema 2: Los números n² y n+p son coprimos.

-- Dem.: Sea k = mcd(n^2,n+p). Por k divide n², luego k divide a n;

-- además, k divide a n+p y (usando el lema 1 y el ser p primo), se

-- tiene que k = 1.

-- Puesto que n³+n²p = n²(n+p) es un cubo, usando el lema 2, se tiene

-- que n² y n+p son cubos y, por serlo n², n también es un cubo. Es

-- decir, existen enteros positivos x e y tales que n = x³ y

-- n+p = y³. Por tanto, p = y³-x³. Sea k = y-x. Se tiene que k = 1 ya

-- que

-- p = y³-x³

-- = (n+k)³-n³

-- = 3k+3k²+k³

-- no es primo para k > 1.

-- Por consiguiente, p = (x+1)³-x³.

El código se encuentra en GitHub.

Primos cubanos

PorJosé A. Alonso 15-julio-20229-julio-2022

Un primo cubano es un número primo que se puede escribir como diferencia de dos cubos consecutivos. Por ejemplo, el 61 es un primo cubano porque es primo y 61 = 5³-4³.

Definir la sucesión

   cubanos :: [Integer]

1	cubanos :: [Integer]

tal que sus elementos son los números cubanos. Por ejemplo,

   λ> take 15 cubanos
   [7,19,37,61,127,271,331,397,547,631,919,1657,1801,1951,2269]

1 2	λ> take 15 cubanos [7,19,37,61,127,271,331,397,547,631,919,1657,1801,1951,2269]

Soluciones

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

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

cubanos1 :: [Integer]
cubanos1 = filter isPrime (zipWith (-) (tail cubos) cubos) 

-- cubos es la lista de los cubos. Por ejemplo,
--    λ> take 10 cubos
--    [1,8,27,64,125,216,343,512,729,1000]
cubos :: [Integer]
cubos = map (^3) [1..]

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

cubanos2 :: [Integer]
cubanos2 = filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

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

cubanos3 :: [Integer]
cubanos3 = filter isPrime [3*x^2 + 3*x + 1 | x <- [1..]]

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

-- La propiedad es
prop_cubanos :: NonNegative Int -> Bool
prop_cubanos (NonNegative n) =
  all (== cubanos1 !! n)
      [cubanos2 !! n,
       cubanos3 !! n]

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

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

-- La comparación es
--    λ> cubanos1 !! 3000
--    795066361
--    (4.21 secs, 16,953,612,192 bytes)
--    λ> cubanos2 !! 3000
--    795066361
--    (4.27 secs, 16,962,597,288 bytes)
--    λ> cubanos3 !! 3000
--    795066361
--    (4.29 secs, 16,956,085,672 bytes)

import Data.Numbers.Primes (isPrime)

import Test.QuickCheck

-- 1ª solución

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

cubanos1 :: [Integer]

cubanos1 = filter isPrime (zipWith (-) (tail cubos) cubos)

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

-- λ> take 10 cubos

-- [1,8,27,64,125,216,343,512,729,1000]

cubos :: [Integer]

cubos = map (^3) [1..]

-- 2ª solución

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

cubanos2 :: [Integer]

cubanos2 = filter isPrime [(x+1)^3 - x^3 | x <- [1..]]

-- 3ª solución

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

cubanos3 :: [Integer]

cubanos3 = filter isPrime [3*x^2 + 3*x + 1 | x <- [1..]]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_cubanos :: NonNegative Int -> Bool

prop_cubanos (NonNegative n) =

all (== cubanos1 !! n)

[cubanos2 !! n,

cubanos3 !! n]

-- La comprobación es

-- λ> quickCheck prop_cubanos

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> cubanos1 !! 3000

-- 795066361

-- (4.21 secs, 16,953,612,192 bytes)

-- λ> cubanos2 !! 3000

-- 795066361

-- (4.27 secs, 16,962,597,288 bytes)

-- λ> cubanos3 !! 3000

-- 795066361

-- (4.29 secs, 16,956,085,672 bytes)

El código se encuentra en GitHub.

Clausura transitiva de una relación binaria

PorJosé A. Alonso 14-julio-20228-julio-2022

La clausura transitiva de una relación binaria R es la relación transitiva que contiene a R. Se puede calcular
usando la composición de relaciones. Veamos un ejemplo, en el que (R ∘ S) representa la composición de R y S: sea

   R = [(1,2),(2,5),(5,6)]

1	R = [(1,2),(2,5),(5,6)]

la relación R no es transitiva ya que (1,2) y (1,5) pertenecen a R pero (1,5) no pertenece; sea

   R1 = R ∪ (R ∘ R)
      = [(1,2),(2,5),(5,6),(1,5),(2,6)]

1 2	R1 = R ∪ (R ∘ R) = [(1,2),(2,5),(5,6),(1,5),(2,6)]

la relación R1 tampoco es transitiva ya que (1,2) y (2,6) pertenecen a R pero (1,6) no pertenece; sea

   R2 = R1 ∪ (R1 ∘ R1)
      = [(1,2),(2,5),(5,6),(1,5),(2,6),(1,6)]

1 2	R2 = R1 ∪ (R1 ∘ R1) = [(1,2),(2,5),(5,6),(1,5),(2,6),(1,6)]

La relación R2 es transitiva y contiene a R. Además, R2 es la clausura transitiva de R.

Definir la función

   clausuraTransitiva :: Ord a => [(a,a)] -> [(a,a)]

1	clausuraTransitiva :: Ord a => [(a,a)] -> [(a,a)]

tal que (clausuraTransitiva r) es la clausura transitiva de r; es decir, la menor relación transitiva que contiene a r. Por ejemplo,

   λ> clausuraTransitiva [(1,2),(2,5),(5,6)]
   [(1,2),(2,5),(5,6),(1,5),(2,6),(1,6)]
   λ> clausuraTransitiva [(1,2),(2,5),(5,6),(6,3)]
   [(1,2),(2,5),(5,6),(6,3),(1,5),(2,6),(5,3),(1,6),(2,3),(1,3)]
   λ> length (clausuraTransitiva [(n,n+1) | n <- [1..100]])
   5050

λ> clausuraTransitiva [(1,2),(2,5),(5,6)]

[(1,2),(2,5),(5,6),(1,5),(2,6),(1,6)]

λ> clausuraTransitiva [(1,2),(2,5),(5,6),(6,3)]

[(1,2),(2,5),(5,6),(6,3),(1,5),(2,6),(5,3),(1,6),(2,3),(1,3)]

λ> length (clausuraTransitiva [(n,n+1) | n <- [1..100]])

5050

Soluciones

import Data.List (union, nub, sort)
import Data.Maybe (mapMaybe)
import qualified Data.Map as M (Map, assocs, empty, insertWith, lookup, map)
import Test.QuickCheck (quickCheck)

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

clausuraTransitiva1 :: Ord a => [(a,a)] -> [(a,a)]  
clausuraTransitiva1 r
  | transitiva r = r
  | otherwise    = clausuraTransitiva1 r1
  where r1 = r `union` composicion r r

-- (transitiva r) se verifica si la relación r es transitiva. Por
-- ejemplo, 
--    transitiva [(1,1),(1,3),(3,1),(3,3),(5,5)]  ==  True
--    transitiva [(1,1),(1,3),(3,1),(5,5)]        ==  False
transitiva :: Ord a => [(a,a)] -> Bool
transitiva r = subconjunto (composicion r r) r

-- (composicion r s) es la composición de las relaciones binarias r y
-- s. Por ejemplo, 
--    λ> composicion [(1,2)] [(2,3),(2,4)]
--    [(1,3),(1,4)]
--    λ> composicion [(1,2),(5,2)] [(2,3),(2,4)]
--    [(1,3),(1,4),(5,3),(5,4)]
--    λ> composicion [(1,2),(1,4),(1,5)] [(2,3),(4,3)]
--    [(1,3)]
composicion :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion r s = nub [(x,y) | (x,u) <- r, (v,y) <- s, u == v] 

-- (subconjunto xs ys) se verifica si xs es un subconjunto de xs. Por
-- ejemplo, 
--    subconjunto [1,3] [3,1,5]  ==  True
--    subconjunto [3,1,5] [1,3]  ==  False
subconjunto :: Ord a => [a] -> [a] -> Bool
subconjunto xs ys = all (`elem` ys) xs

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

clausuraTransitiva2 :: Ord a => [(a,a)] -> [(a,a)]  
clausuraTransitiva2 r
  | r1 == r   = r
  | otherwise = clausuraTransitiva2 r1
  where r1 = r `union` composicion r r

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

-- La propiedad es
prop_clausuraTransitiva :: [(Int,Int)] -> Bool
prop_clausuraTransitiva r =
  all (== sort (clausuraTransitiva1 r))
      [sort (clausuraTransitiva2 r),
       sort (clausuraTransitiva3 r)]
  
-- La comprobación es
--    λ> quickCheck prop_clausuraTransitiva
--    +++ OK, passed 100 tests.

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

clausuraTransitiva3 :: Ord a => [(a,a)] -> [(a,a)]  
clausuraTransitiva3 r
  | transitiva3 r = r
  | otherwise     = clausuraTransitiva3 r1
  where r1 = r `union` composicion3 r r

transitiva3 :: Ord a => [(a,a)] -> Bool
transitiva3 r = subconjunto (composicion3 r r) r

composicion3 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]
composicion3 r s =
  relAlista (composicionRel (listaArel r) (listaArel s))

-- Una relación se puede representar por un diccionario donde las claves
-- son los elementos y los valores son las listas de los elementos con
-- los que se relaciona.
type Rel a = M.Map a [a]

-- (listaArel xys) es la relación correspondiente a la lista de pares
-- xys. Por ejemplo.
--    λ> listaArel [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]
--    fromList [(1,[1,2,3,6]),(2,[2,6]),(3,[3,6]),(6,[6])]
listaArel :: Ord a => [(a,a)] -> Rel a
listaArel []          = M.empty
listaArel ((x,y):xys) = M.insertWith (++) x [y] (listaArel xys)

-- (composicionRel r s) es la composición de las relaciones r y s. Por
-- ejemplo,
--    λ> r = listaArel [(1,2),(5,2)]
--    λ> s = listaArel [(2,3),(2,4)]
--    λ> composicionRel r s
--    fromList [(1,[3,4]),(5,[3,4])]
composicionRel :: Ord a => Rel a -> Rel a -> Rel a
composicionRel r s =
  M.map f r 
  where f xs = concat (mapMaybe (`M.lookup` s) xs)

-- (relAlista r) es la lista de pares correspondientes a la relación
-- r. Por ejemplo,
--    λ> relAlista (M.fromList [(1,[3,4]),(5,[3,4])])
--    [(1,3),(1,4),(5,3),(5,4)]
relAlista :: Ord a => Rel a -> [(a,a)]
relAlista r =
  nub [(x,y) | (x,ys) <- M.assocs r, y <- ys] 

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

-- La comparación es
--    λ> length (clausuraTransitiva1 [(n,n+1) | n <- [1..60]])
--    1830
--    (2.15 secs, 453,533,992 bytes)
--    λ> length (clausuraTransitiva2 [(n,n+1) | n <- [1..60]])
--    1830
--    (2.23 secs, 558,571,904 bytes)
--    λ> length (clausuraTransitiva3 [(n,n+1) | n <- [1..60]])
--    1830
--    (0.25 secs, 207,168,552 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

import Data.List (union, nub, sort)

import Data.Maybe (mapMaybe)

import qualified Data.Map as M (Map, assocs, empty, insertWith, lookup, map)

import Test.QuickCheck (quickCheck)

-- 1ª solución

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

clausuraTransitiva1 :: Ord a => [(a,a)] -> [(a,a)]

clausuraTransitiva1 r

| transitiva r = r

| otherwise = clausuraTransitiva1 r1

where r1 = r `union` composicion r r

-- (transitiva r) se verifica si la relación r es transitiva. Por

-- ejemplo,

-- transitiva [(1,1),(1,3),(3,1),(3,3),(5,5)] == True

-- transitiva [(1,1),(1,3),(3,1),(5,5)] == False

transitiva :: Ord a => [(a,a)] -> Bool

transitiva r = subconjunto (composicion r r) r

-- (composicion r s) es la composición de las relaciones binarias r y

-- s. Por ejemplo,

-- λ> composicion [(1,2)] [(2,3),(2,4)]

-- [(1,3),(1,4)]

-- λ> composicion [(1,2),(5,2)] [(2,3),(2,4)]

-- [(1,3),(1,4),(5,3),(5,4)]

-- λ> composicion [(1,2),(1,4),(1,5)] [(2,3),(4,3)]

-- [(1,3)]

composicion :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion r s = nub [(x,y) | (x,u) <- r, (v,y) <- s, u == v]

-- (subconjunto xs ys) se verifica si xs es un subconjunto de xs. Por

-- ejemplo,

-- subconjunto [1,3] [3,1,5] == True

-- subconjunto [3,1,5] [1,3] == False

subconjunto :: Ord a => [a] -> [a] -> Bool

subconjunto xs ys = all (`elem` ys) xs

-- 2ª solución

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

clausuraTransitiva2 :: Ord a => [(a,a)] -> [(a,a)]

clausuraTransitiva2 r

| r1 == r = r

| otherwise = clausuraTransitiva2 r1

where r1 = r `union` composicion r r

-- Comprobación de equivalencia

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

-- La propiedad es

prop_clausuraTransitiva :: [(Int,Int)] -> Bool

prop_clausuraTransitiva r =

all (== sort (clausuraTransitiva1 r))

[sort (clausuraTransitiva2 r),

sort (clausuraTransitiva3 r)]

-- La comprobación es

-- λ> quickCheck prop_clausuraTransitiva

-- +++ OK, passed 100 tests.

-- 3ª solución

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

clausuraTransitiva3 :: Ord a => [(a,a)] -> [(a,a)]

clausuraTransitiva3 r

| transitiva3 r = r

| otherwise = clausuraTransitiva3 r1

where r1 = r `union` composicion3 r r

transitiva3 :: Ord a => [(a,a)] -> Bool

transitiva3 r = subconjunto (composicion3 r r) r

composicion3 :: Ord a => [(a,a)] -> [(a,a)] -> [(a,a)]

composicion3 r s =

relAlista (composicionRel (listaArel r) (listaArel s))

-- Una relación se puede representar por un diccionario donde las claves

-- son los elementos y los valores son las listas de los elementos con

-- los que se relaciona.

type Rel a = M.Map a [a]

-- (listaArel xys) es la relación correspondiente a la lista de pares

-- xys. Por ejemplo.

-- λ> listaArel [(1,1),(1,2),(1,3),(1,6),(2,2),(2,6),(3,3),(3,6),(6,6)]

-- fromList [(1,[1,2,3,6]),(2,[2,6]),(3,[3,6]),(6,[6])]

listaArel :: Ord a => [(a,a)] -> Rel a

listaArel [] = M.empty

listaArel ((x,y):xys) = M.insertWith (++) x [y] (listaArel xys)

-- (composicionRel r s) es la composición de las relaciones r y s. Por

-- ejemplo,

-- λ> r = listaArel [(1,2),(5,2)]

-- λ> s = listaArel [(2,3),(2,4)]

-- λ> composicionRel r s

-- fromList [(1,[3,4]),(5,[3,4])]

composicionRel :: Ord a => Rel a -> Rel a -> Rel a

composicionRel r s =

M.map f r

where f xs = concat (mapMaybe (`M.lookup` s) xs)

-- (relAlista r) es la lista de pares correspondientes a la relación

-- r. Por ejemplo,

-- λ> relAlista (M.fromList [(1,[3,4]),(5,[3,4])])

-- [(1,3),(1,4),(5,3),(5,4)]

relAlista :: Ord a => Rel a -> [(a,a)]

relAlista r =

nub [(x,y) | (x,ys) <- M.assocs r, y <- ys]

-- Comparación de eficiencia

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

-- La comparación es

-- λ> length (clausuraTransitiva1 [(n,n+1) | n <- [1..60]])

-- 1830

-- (2.15 secs, 453,533,992 bytes)

-- λ> length (clausuraTransitiva2 [(n,n+1) | n <- [1..60]])

-- 1830

-- (2.23 secs, 558,571,904 bytes)

-- λ> length (clausuraTransitiva3 [(n,n+1) | n <- [1..60]])

-- 1830

-- (0.25 secs, 207,168,552 bytes)

El código se encuentra en GitHub.

Soluciones

Referencias

Soluciones

Referencias

Soluciones

Soluciones

Referencias

Soluciones

Soluciones

Referencia

Soluciones

Referencias

Soluciones

Referencias

Soluciones

Referencias

Soluciones

Referencias

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Soluciones

Entradas recientes

Buscar

Correo electrónico