Haskell y Python

Factorizaciones de números de Hilbert

PorJosé A. Alonso 24-diciembre-202318-febrero-2024

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 divisible 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).
Read More «Factorizaciones de números de Hilbert»

El teorema de Navidad de Fermat

PorJosé A. Alonso 14-diciembre-20232-febrero-2024

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

Definir las funciones

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

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

primosImparesConRepresentacionUnica :: [Integer]

primos4nM1 :: [Integer]

tales que

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

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

representaciones 20 == [(2,4)]

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

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

representaciones 100000147984 == [(0,316228)]

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

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

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

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

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

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

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

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

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

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

Números de Pentanacci

PorJosé A. Alonso 9-diciembre-202322-diciembre-2023

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 las funciones

   pentanacci  :: Integer -> Integer
   pentanaccis :: [Integer]

1 2	pentanacci :: Integer -> Integer pentanaccis :: [Integer]

tales que

pentanacci n es el n-ésimo número de Pentanacci. Por ejemplo,

     λ> pentanacci 14
     3525
     λ> pentanacci (10^5) `mod` 10^30
     482929150584077921552549215816
     λ> length (show (pentanacci (10^5)))
     29357

λ> pentanacci 14

3525

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

482929150584077921552549215816

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

29357

pentanaccis es la lista de 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]

1 2	λ> take 15 pentanacci [0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]

Soluciones

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

Soluciones en Haskell

module Numeros_de_Pentanacci where

import Data.List (genericIndex, zipWith5)
import Test.Hspec (Spec, hspec, it, shouldBe)
import Test.QuickCheck (NonNegative (NonNegative), quickCheckWith, maxSize, stdArgs)

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

pentanacci1  :: Integer -> Integer
pentanacci1  0 = 0
pentanacci1  1 = 1
pentanacci1  2 = 1
pentanacci1  3 = 2
pentanacci1  4 = 4
pentanacci1  n = sum [pentanacci1 (n-k) | k <- [1..5]]

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

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

pentanaccis2 :: [Integer]
pentanaccis2 =
  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 pentanaccis2

pentanacci2  :: Integer -> Integer
pentanacci2 n = pentanaccis2 `genericIndex` n

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

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

pentanacci3  :: Integer -> Integer
pentanacci3 n = pentanaccis3 `genericIndex` n

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

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

pentanacci4  :: Integer -> Integer
pentanacci4 n = pentanaccis4 `genericIndex` n

-- Verificación
-- ============

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "ej1" $
    take 15 pentanaccis1 `shouldBe`
    [0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]
  it "ej2" $
    take 15 pentanaccis2 `shouldBe`
    [0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]
  it "ej3" $
    take 15 pentanaccis3 `shouldBe`
    [0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]
  it "ej4" $
    take 15 pentanaccis4 `shouldBe`
    [0,1,1,2,4,8,16,31,61,120,236,464,912,1793,3525]

-- La verificación es
--    λ> verifica
--
--    4 examples, 0 failures

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

-- La propiedad es
prop_pentanaccis :: NonNegative Int -> Bool
prop_pentanaccis (NonNegative n) =
  all (== pentanaccis1 !! n)
      [pentanaccis1 !! n,
       pentanaccis2 !! n,
       pentanaccis3 !! n]

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

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

-- La comparación es
--    λ> pentanacci1 25
--    5976577
--    (4.64 secs, 1,865,626,496 bytes)
--    λ> pentanacci2 25
--    5976577
--    (0.00 secs, 578,584 bytes)
--
--    λ> length (show (pentanacci2 (10^5)))
--    29357
--    (1.16 secs, 2,543,272,136 bytes)
--    λ> length (show (pentanacci3 (10^5)))
--    29357
--    (1.00 secs, 2,560,881,608 bytes)
--    λ> length (show (pentanacci4 (10^5)))
--    29357
--    (1.03 secs, 2,592,078,744 bytes)

100

101

102

103

104

105

106

107

108

109

110

111

112

113

module Numeros_de_Pentanacci where

import Data.List (genericIndex, zipWith5)

import Test.Hspec (Spec, hspec, it, shouldBe)

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

-- 1ª solución

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

pentanacci1 :: Integer -> Integer

pentanacci1 0 = 0

pentanacci1 1 = 1

pentanacci1 2 = 1

pentanacci1 3 = 2

pentanacci1 4 = 4

pentanacci1 n = sum [pentanacci1 (n-k) | k <- [1..5]]

pentanaccis1 :: [Integer]

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

-- 2ª solución

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

pentanaccis2 :: [Integer]

pentanaccis2 =

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 pentanaccis2

pentanacci2 :: Integer -> Integer

pentanacci2 n = pentanaccis2 `genericIndex` n

-- 3ª solución

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

pentanaccis3 :: [Integer]

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

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

pentanacci3 :: Integer -> Integer

pentanacci3 n = pentanaccis3 `genericIndex` n

-- 4ª solución

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

pentanaccis4 :: [Integer]

pentanaccis4 = 0: 1: 1: 2: 4: p pentanaccis4

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

pentanacci4 :: Integer -> Integer

pentanacci4 n = pentanaccis4 `genericIndex` n

-- Verificación

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

verifica :: IO ()

verifica = hspec spec

spec :: Spec

spec = do

it "ej1" $

take 15 pentanaccis1 `shouldBe`

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

it "ej2" $

take 15 pentanaccis2 `shouldBe`

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

it "ej3" $

take 15 pentanaccis3 `shouldBe`

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

it "ej4" $

take 15 pentanaccis4 `shouldBe`

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

-- La verificación es

-- λ> verifica

-- 4 examples, 0 failures

-- Comprobación de equivalencia

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

-- La propiedad es

prop_pentanaccis :: NonNegative Int -> Bool

prop_pentanaccis (NonNegative n) =

all (== pentanaccis1 !! n)

[pentanaccis1 !! n,

pentanaccis2 !! n,

pentanaccis3 !! n]

-- La comprobación es

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

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> pentanacci1 25

-- 5976577

-- (4.64 secs, 1,865,626,496 bytes)

-- λ> pentanacci2 25

-- 5976577

-- (0.00 secs, 578,584 bytes)

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

-- 29357

-- (1.16 secs, 2,543,272,136 bytes)

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

-- 29357

-- (1.00 secs, 2,560,881,608 bytes)

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

-- 29357

-- (1.03 secs, 2,592,078,744 bytes)

Soluciones en Python

from itertools import count, islice
from sys import set_int_max_str_digits
from timeit import Timer, default_timer
from typing import Iterator

set_int_max_str_digits(30000)

# 1ª solución
# ===========

def pentanacci1(n: int) -> int:
    if n == 0:
        return 0
    if n == 1:
        return 1
    if n == 2:
        return 1
    if n == 3:
        return 2
    if n == 4:
        return 4
    return sum((pentanacci1(n-k) for k in range(1, 6)))

def pentanaccis1() -> Iterator[int]:
    return (pentanacci1(n) for n in count())

# 2ª solución
# ===========

def pentanaccis2() -> Iterator[int]:
    seq = [0, 1, 1, 2, 4]
    while True:
        yield seq[0]
        seq.append(sum(seq))
        seq.pop(0)

def nth(i: Iterator[int], n: int) -> int:
    return list(islice(i, n, n+1))[0]

def pentanacci2(n: int) -> int:
    return nth(pentanaccis2(), n)

# Verificación
# ============

def test_pentanacci() -> None:
    assert pentanacci1(14) == 3525
    assert list(islice(pentanaccis1(), 15)) == \
        [0, 1, 1, 2, 4, 8, 16, 31, 61, 120, 236, 464, 912, 1793, 3525]
    assert pentanacci2(14) == 3525
    assert list(islice(pentanaccis2(), 15)) == \
        [0, 1, 1, 2, 4, 8, 16, 31, 61, 120, 236, 464, 912, 1793, 3525]
    print("Verificado")

# La verificación es
#    >>> test_pentanacci()
#    Verificado

# Comprobación de equivalencia
# ============================

# La propiedad es
def test_pentanacci_equiv() -> bool:
    return list(islice(pentanaccis1(), 25)) == list(islice(pentanaccis2(), 25))

# La comprobación es
#    >>> test_pentanacci_equiv()
#    True

# Comparación de eficiencia
# =========================

# Comparación de eficiencia
# =========================

def tiempo(e: str) -> None:
    """Tiempo (en segundos) de evaluar la expresión e."""
    t = Timer(e, "", default_timer, globals()).timeit(1)
    print(f"{t:0.2f} segundos")

# La comparación es
#    >>> tiempo('pentanacci1(28)')
#    8.24 segundos
#    >>> tiempo('pentanacci2(28)')
#    0.00 segundos

from itertools import count, islice

from sys import set_int_max_str_digits

from timeit import Timer, default_timer

from typing import Iterator

set_int_max_str_digits(30000)

# 1ª solución

# ===========

def pentanacci1(n: int) -> int:

if n == 0:

return 0

if n == 1:

return 1

if n == 2:

return 1

if n == 3:

return 2

if n == 4:

return 4

return sum((pentanacci1(n-k) for k in range(1, 6)))

def pentanaccis1() -> Iterator[int]:

return (pentanacci1(n) for n in count())

# 2ª solución

# ===========

def pentanaccis2() -> Iterator[int]:

seq = [0, 1, 1, 2, 4]

while True:

yield seq[0]

seq.append(sum(seq))

seq.pop(0)

def nth(i: Iterator[int], n: int) -> int:

return list(islice(i, n, n+1))[0]

def pentanacci2(n: int) -> int:

return nth(pentanaccis2(), n)

# Verificación

# ============

def test_pentanacci() -> None:

assert pentanacci1(14) == 3525

assert list(islice(pentanaccis1(), 15)) == \

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

assert pentanacci2(14) == 3525

assert list(islice(pentanaccis2(), 15)) == \

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

print("Verificado")

# La verificación es

# >>> test_pentanacci()

# Verificado

# Comprobación de equivalencia

# ============================

# La propiedad es

def test_pentanacci_equiv() -> bool:

return list(islice(pentanaccis1(), 25)) == list(islice(pentanaccis2(), 25))

# La comprobación es

# >>> test_pentanacci_equiv()

# True

# Comparación de eficiencia

# =========================

# Comparación de eficiencia

# =========================

def tiempo(e: str) -> None:

"""Tiempo (en segundos) de evaluar la expresión e."""

t = Timer(e, "", default_timer, globals()).timeit(1)

print(f"{t:0.2f} segundos")

# La comparación es

# >>> tiempo('pentanacci1(28)')

# 8.24 segundos

# >>> tiempo('pentanacci2(28)')

# 0.00 segundos

Algoritmo de bajada para resolver un sistema triangular inferior

PorJosé A. Alonso 4-diciembre-202321-noviembre-2023

Un sistema de ecuaciones lineales $Ax = b$ es triangular inferior si todos los elementos de la matriz $A$ que están por encima de la diagonal principal son nulos; es decir, es de la forma
\begin{align}
&a_{1 1}x_1 &= b_1 \\
&a_{2 1}x_1 + a_{2 2}x_2 &= b_2 \\
&a_{3 1}x_1 + a_{3 2}x_2 + a_{3 3}x_3 &= b_3 \\
&… & \\
&a_{n 1}x_1 + a_{n 2}x_2 + a_{n 3}x_3 +…+ a_{n n}x_n &= b_n
\end{align}

El sistema es compatible si, y sólo si, el producto de los elementos de la diagonal principal es distinto de cero. En este caso, la solución se puede calcular mediante el algoritmo de bajada:
\begin{align}
x_1 &= \frac{b_1}{a_{1 1}} \\
x_2 &= \frac{b_2 – a_{2 1}x_1}{a_{2 2}} \\
x_3 &= \frac{b_3 – a_{3 1}x_1 – a_{3 2}x_2}{a_{3 3}} \\
… \\
x_n &= \frac{b_n – a_{n 1}x_1 – a_{n 2}x_2 – … – a_{n,n-1}x_{n-1}}{a_{n n}}
\end{align}

Definir la función

   bajada :: Matrix Double -> Matrix Double -> Matrix Double

1	bajada :: Matrix Double -> Matrix Double -> Matrix Double

tal que bajada a b es la solución, mediante el algoritmo de bajada, del sistema compatible triangular superior ax = b. Por ejemplo,

   λ> let a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]
   λ> let b = fromLists [[3],[6.5],[10]]
   λ> bajada a b
   ( 1.5 )
   ( 2.0 )
   ( 0.0 )

λ> let a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]

λ> let b = fromLists [[3],[6.5],[10]]

λ> bajada a b

( 1.5 )

( 2.0 )

( 0.0 )

Es decir, la solución del sistema
\begin{align}
2x &= 3 \\
3x + y &= 6.5 \\
4x + 2y + 5z &= 10
\end{align}
es $x=1.5$, $y=2$ y $z=0$.

Soluciones

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

Soluciones en Haskell

module Algoritmo_de_bajada where

import Data.Matrix (Matrix, (!), fromLists, nrows, toLists)
import Test.Hspec (Spec, hspec, it, shouldBe)

bajada :: Matrix Double -> Matrix Double -> Matrix Double
bajada a b = fromLists [[x i] | i <- [1..m]]
  where m   = nrows a
        x k = (b!(k,1) - sum [a!(k,j) * x j | j <- [1..k-1]]) / a!(k,k)

-- Verificación
-- ============

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "ej1" $
    toLists (bajada a b) `shouldBe` [[1.5],[2.0],[0.0]]
    where
      a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]
      b = fromLists [[3],[6.5],[10]]

-- La verificación es
--    λ> verifica
--
--    Finished in 0.0007 seconds
--    1 example, 0 failures

module Algoritmo_de_bajada where

import Data.Matrix (Matrix, (!), fromLists, nrows, toLists)

import Test.Hspec (Spec, hspec, it, shouldBe)

bajada :: Matrix Double -> Matrix Double -> Matrix Double

bajada a b = fromLists [[x i] | i <- [1..m]]

where m = nrows a

x k = (b!(k,1) - sum [a!(k,j) * x j | j <- [1..k-1]]) / a!(k,k)

-- Verificación

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

verifica :: IO ()

verifica = hspec spec

spec :: Spec

spec = do

it "ej1" $

toLists (bajada a b) `shouldBe` [[1.5],[2.0],[0.0]]

where

a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]

b = fromLists [[3],[6.5],[10]]

-- La verificación es

-- λ> verifica

-- Finished in 0.0007 seconds

-- 1 example, 0 failures

Soluciones en Python

def bajada(a: list[list[float]], b: list[list[float]]) -> list[list[float]]:
    n = len(a)
    def x(k: int) -> float:
        return (b[k][0] - sum((a[k][j] * x(j) for j in range(0, k)))) / a[k][k]
    return [[x(i)] for i in range(0, n)]

# Verificación
# ============

def test_bajada() -> None:
    assert bajada([[2,0,0],[3,1,0],[4,2,5.0]], [[3],[6.5],[10]]) == \
        [[1.5], [2.0], [0.0]]
    print("Verificado")

# La verificación es
#    >>> test_bajada()
#    Verificado

def bajada(a: list[list[float]], b: list[list[float]]) -> list[list[float]]:

n = len(a)

def x(k: int) -> float:

return (b[k][0] - sum((a[k][j] * x(j) for j in range(0, k)))) / a[k][k]

return [[x(i)] for i in range(0, n)]

# Verificación

# ============

def test_bajada() -> None:

assert bajada([[2,0,0],[3,1,0],[4,2,5.0]], [[3],[6.5],[10]]) == \

[[1.5], [2.0], [0.0]]

print("Verificado")

# La verificación es

# >>> test_bajada()

# Verificado

Integración por el método de los rectángulos

PorJosé A. Alonso 29-noviembre-202320-noviembre-2023

La integral definida de una función f entre los límites a y b puede calcularse mediante la regla del rectángulo (ver en http://bit.ly/1FDhZ1z) usando la fórmula
\[ h(f(a+\frac{h}{2}) + f(a+h+\frac{h}{2}) + f(a+2h+\frac{h}{2}) + … + f(a+nh+\frac{h}{2}))\]
con $a+nh+\dfrac{h}{2} \leq b < a+(n+1)h+\dfrac{h}{2}$ y usando valores pequeños para $h$.

Definir la función

   integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

1	integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

tal que integral a b f h es el valor de dicha expresión. Por ejemplo, el cálculo de la integral de (f(x) = x^3) entre 0 y 1, con paso 0.01, es

   integral 0 1 (^3) 0.01  ==  0.24998750000000042

1	integral 0 1 (^3) 0.01 == 0.24998750000000042

Otros ejemplos son

   integral 0 1 (^4) 0.01                        ==  0.19998333362500048
   integral 0 1 (\x -> 3*x^2 + 4*x^3) 0.01       ==  1.9999250000000026
   log 2 - integral 1 2 (\x -> 1/x) 0.01         ==  3.124931644782336e-6
   pi - 4 * integral 0 1 (\x -> 1/(x^2+1)) 0.01  ==  -8.333333331389525e-6

integral 0 1 (^4) 0.01 == 0.19998333362500048

integral 0 1 (\x -> 3*x^2 + 4*x^3) 0.01 == 1.9999250000000026

log 2 - integral 1 2 (\x -> 1/x) 0.01 == 3.124931644782336e-6

pi - 4 * integral 0 1 (\x -> 1/(x^2+1)) 0.01 == -8.333333331389525e-6

Soluciones

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

Soluciones en Haskell

module Integracion_por_rectangulos where

import Test.Hspec (Spec, hspec, it, shouldBe, shouldSatisfy)

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

integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral a b f h = h * suma (a+h/2) b (+h) f

-- (suma a b s f) es l valor de
--    f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))
-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,
--    suma 2 5 (1+) (^3)  ==  224
suma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a
suma a b s f = sum [f x | x <- sucesion a b s]

-- (sucesion x y s) es la lista
--    [a, s(a), s(s(a), ..., s(...(s(a))...)]
-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,
--    sucesion 3 20 (+2)  ==  [3,5,7,9,11,13,15,17,19]
sucesion :: Ord a => a -> a -> (a -> a) -> [a]
sucesion a b s = takeWhile (<=b) (iterate s a)

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

integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral2 a b f h
  | a+h/2 > b = 0
  | otherwise = h * f (a+h/2) + integral2 (a+h) b f h

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

integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral3 a b f h = aux a where
  aux x | x+h/2 > b = 0
        | otherwise = h * f (x+h/2) + aux (x+h)

-- Comparación de eficiencia
--    λ> integral 0 10 (^3) 0.00001
--    2499.9999998811422
--    (4.62 secs, 1084774336 bytes)
--    λ> integral2 0 10 (^3) 0.00001
--    2499.999999881125
--    (7.90 secs, 1833360768 bytes)
--    λ> integral3 0 10 (^3) 0.00001
--    2499.999999881125
--    (7.27 secs, 1686056080 bytes)

-- Verificación
-- ============

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "e1" $
    integral' 0 1 (^(3::Int)) 0.01 `shouldBe` 0.24998750000000042
  it "e2" $
    integral' 0 1 (^(4::Int)) 0.01 `shouldBe` 0.19998333362500048
  it "e3" $
    integral' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldBe` 1.9999250000000026
  it "e4" $
    log 2 - integral' 1 2 (1 /) 0.01 `shouldBe` 3.124931644782336e-6
  it "e5" $
    pi - 4 * integral' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldBe` -8.333333331389525e-6
  it "e1b" $
    integral2' 0 1 (^(3::Int)) 0.01 `shouldSatisfy` (~= 0.24998750000000042)
  it "e2b" $
    integral2' 0 1 (^(4::Int)) 0.01 `shouldSatisfy` (~= 0.19998333362500048)
  it "e3b" $
    integral2' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldSatisfy` (~= 1.9999250000000026)
  it "e4b" $
    log 2 - integral2' 1 2 (1 /) 0.01 `shouldSatisfy` (~= 3.124931644782336e-6)
  it "e5b" $
    pi - 4 * integral2' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldSatisfy` (~= (-8.333333331389525e-6))
  it "e1c" $
    integral3' 0 1 (^(3::Int)) 0.01 `shouldSatisfy` (~= 0.24998750000000042)
  it "e2c" $
    integral3' 0 1 (^(4::Int)) 0.01 `shouldSatisfy` (~= 0.19998333362500048)
  it "e3c" $
    integral3' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldSatisfy` (~= 1.9999250000000026)
  it "e4c" $
    log 2 - integral3' 1 2 (1 /) 0.01 `shouldSatisfy` (~= 3.124931644782336e-6)
  it "e5c" $
    pi - 4 * integral3' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldSatisfy` (~= (-8.333333331389525e-6))
  where
    integral', integral2', integral3' :: Double -> Double -> (Double -> Double) -> Double -> Double
    integral'  = integral
    integral2' = integral2
    integral3' = integral3
    a ~= b = abs (a - b) < 0.00001

-- La verificación es
--    λ> verifica
--
--    Finished in 0.0058 seconds
--    15 examples, 0 failures

100

101

module Integracion_por_rectangulos where

import Test.Hspec (Spec, hspec, it, shouldBe, shouldSatisfy)

-- 1ª solución

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

integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral a b f h = h * suma (a+h/2) b (+h) f

-- (suma a b s f) es l valor de

-- f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))

-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,

-- suma 2 5 (1+) (^3) == 224

suma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a

suma a b s f = sum [f x | x <- sucesion a b s]

-- (sucesion x y s) es la lista

-- [a, s(a), s(s(a), ..., s(...(s(a))...)]

-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,

-- sucesion 3 20 (+2) == [3,5,7,9,11,13,15,17,19]

sucesion :: Ord a => a -> a -> (a -> a) -> [a]

sucesion a b s = takeWhile (<=b) (iterate s a)

-- 2ª solución

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

integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral2 a b f h

| a+h/2 > b = 0

| otherwise = h * f (a+h/2) + integral2 (a+h) b f h

-- 3ª solución

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

integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral3 a b f h = aux a where

aux x | x+h/2 > b = 0

| otherwise = h * f (x+h/2) + aux (x+h)

-- Comparación de eficiencia

-- λ> integral 0 10 (^3) 0.00001

-- 2499.9999998811422

-- (4.62 secs, 1084774336 bytes)

-- λ> integral2 0 10 (^3) 0.00001

-- 2499.999999881125

-- (7.90 secs, 1833360768 bytes)

-- λ> integral3 0 10 (^3) 0.00001

-- 2499.999999881125

-- (7.27 secs, 1686056080 bytes)

-- Verificación

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

verifica :: IO ()

verifica = hspec spec

spec :: Spec

spec = do

it "e1" $

integral' 0 1 (^(3::Int)) 0.01 `shouldBe` 0.24998750000000042

it "e2" $

integral' 0 1 (^(4::Int)) 0.01 `shouldBe` 0.19998333362500048

it "e3" $

integral' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldBe` 1.9999250000000026

it "e4" $

log 2 - integral' 1 2 (1 /) 0.01 `shouldBe` 3.124931644782336e-6

it "e5" $

pi - 4 * integral' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldBe` -8.333333331389525e-6

it "e1b" $

integral2' 0 1 (^(3::Int)) 0.01 `shouldSatisfy` (~= 0.24998750000000042)

it "e2b" $

integral2' 0 1 (^(4::Int)) 0.01 `shouldSatisfy` (~= 0.19998333362500048)

it "e3b" $

integral2' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldSatisfy` (~= 1.9999250000000026)

it "e4b" $

log 2 - integral2' 1 2 (1 /) 0.01 `shouldSatisfy` (~= 3.124931644782336e-6)

it "e5b" $

pi - 4 * integral2' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldSatisfy` (~= (-8.333333331389525e-6))

it "e1c" $

integral3' 0 1 (^(3::Int)) 0.01 `shouldSatisfy` (~= 0.24998750000000042)

it "e2c" $

integral3' 0 1 (^(4::Int)) 0.01 `shouldSatisfy` (~= 0.19998333362500048)

it "e3c" $

integral3' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldSatisfy` (~= 1.9999250000000026)

it "e4c" $

log 2 - integral3' 1 2 (1 /) 0.01 `shouldSatisfy` (~= 3.124931644782336e-6)

it "e5c" $

pi - 4 * integral3' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldSatisfy` (~= (-8.333333331389525e-6))

where

integral', integral2', integral3' :: Double -> Double -> (Double -> Double) -> Double -> Double

integral' = integral

integral2' = integral2

integral3' = integral3

a ~= b = abs (a - b) < 0.00001

-- La verificación es

-- λ> verifica

-- Finished in 0.0058 seconds

-- 15 examples, 0 failures

Soluciones en Python

from math import log, pi
from typing import Callable

# 1ª solución
# ===========

# sucesion(x, y, s) es la lista
#    [a, s(a), s(s(a), ..., s(...(s(a))...)]
# hasta que s(s(...(s(a))...)) > b. Por ejemplo,
#    sucesion(3, 20, lambda x : x+2)  ==  [3,5,7,9,11,13,15,17,19]
def sucesion(a: float, b: float, s: Callable[[float], float]) -> list[float]:
    xs = []
    while a <= b:
        xs.append(a)
        a = s(a)
    return xs

# suma(a, b, s, f) es el valor de
#    f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))
# hasta que s(s(...(s(a))...)) > b. Por ejemplo,
#    suma(2, 5, lambda x: x+1, lambda x: x**3)  ==  224
def suma(a: float,
         b: float,
         s: Callable[[float], float],
         f: Callable[[float], float]) -> float:
    return sum(f(x) for x in sucesion(a, b, s))

def integral(a: float,
             b: float,
             f: Callable[[float], float],
             h: float) -> float:
    return h * suma(a+h/2, b, lambda x: x+h, f)

# 2ª solución
# ===========

def integral2(a: float,
              b: float,
              f: Callable[[float], float],
              h: float) -> float:
    if a+h/2 > b:
        return 0
    return h * f(a+h/2) + integral2(a+h, b, f, h)

# 3ª solución
# ===========

def integral3(a: float,
              b: float,
              f: Callable[[float], float],
              h: float) -> float:
    def aux(x: float) -> float:
        if x+h/2 > b:
            return 0
        return h * f(x+h/2) + aux(x+h)
    return aux(a)

# Verificación
# ============

def test_integral() -> None:
    def aproximado(a: float, b: float) -> bool:
        return abs(a - b) < 0.00001
    assert integral(0, 1, lambda x : x**3, 0.01) == 0.24998750000000042
    assert integral(0, 1, lambda x : x**4, 0.01) == 0.19998333362500054
    assert integral(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01) == 1.9999250000000026
    assert log(2) - integral(1, 2, lambda x : 1/x, 0.01) == 3.124931644782336e-6
    assert pi - 4 * integral(0, 1, lambda x : 1/(x**2+1), 0.01) == -8.333333331389525e-6
    assert aproximado(integral2(0, 1, lambda x : x**3, 0.01),
                      0.24998750000000042)
    assert aproximado(integral2(0, 1, lambda x : x**4, 0.01),
                      0.19998333362500054)
    assert aproximado(integral2(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01),
                      1.9999250000000026)
    assert aproximado(log(2) - integral2(1, 2, lambda x : 1/x, 0.01),
                      3.124931644782336e-6)
    assert aproximado(pi - 4 * integral2(0, 1, lambda x : 1/(x**2+1), 0.01),
                      -8.333333331389525e-6)
    assert aproximado(integral3(0, 1, lambda x : x**3, 0.01),
                      0.24998750000000042)
    assert aproximado(integral3(0, 1, lambda x : x**4, 0.01),
                      0.19998333362500054)
    assert aproximado(integral3(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01),
                      1.9999250000000026)
    assert aproximado(log(2) - integral3(1, 2, lambda x : 1/x, 0.01),
                      3.124931644782336e-6)
    assert aproximado(pi - 4 * integral3(0, 1, lambda x : 1/(x**2+1), 0.01),
                      -8.333333331389525e-6)
    print("Verificado")

# La verificación es
#    >>> test_integral()
#    Verificado

from math import log, pi

from typing import Callable

# 1ª solución

# ===========

# sucesion(x, y, s) es la lista

# [a, s(a), s(s(a), ..., s(...(s(a))...)]

# hasta que s(s(...(s(a))...)) > b. Por ejemplo,

# sucesion(3, 20, lambda x : x+2) == [3,5,7,9,11,13,15,17,19]

def sucesion(a: float, b: float, s: Callable[[float], float]) -> list[float]:

xs = []

while a <= b:

xs.append(a)

a = s(a)

return xs

# suma(a, b, s, f) es el valor de

# f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))

# hasta que s(s(...(s(a))...)) > b. Por ejemplo,

# suma(2, 5, lambda x: x+1, lambda x: x**3) == 224

def suma(a: float,

b: float,

s: Callable[[float], float],

f: Callable[[float], float]) -> float:

return sum(f(x) for x in sucesion(a, b, s))

def integral(a: float,

b: float,

f: Callable[[float], float],

h: float) -> float:

return h * suma(a+h/2, b, lambda x: x+h, f)

# 2ª solución

# ===========

def integral2(a: float,

b: float,

f: Callable[[float], float],

h: float) -> float:

if a+h/2 > b:

return 0

return h * f(a+h/2) + integral2(a+h, b, f, h)

# 3ª solución

# ===========

def integral3(a: float,

b: float,

f: Callable[[float], float],

h: float) -> float:

def aux(x: float) -> float:

if x+h/2 > b:

return 0

return h * f(x+h/2) + aux(x+h)

return aux(a)

# Verificación

# ============

def test_integral() -> None:

def aproximado(a: float, b: float) -> bool:

return abs(a - b) < 0.00001

assert integral(0, 1, lambda x : x**3, 0.01) == 0.24998750000000042

assert integral(0, 1, lambda x : x**4, 0.01) == 0.19998333362500054

assert integral(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01) == 1.9999250000000026

assert log(2) - integral(1, 2, lambda x : 1/x, 0.01) == 3.124931644782336e-6

assert pi - 4 * integral(0, 1, lambda x : 1/(x**2+1), 0.01) == -8.333333331389525e-6

assert aproximado(integral2(0, 1, lambda x : x**3, 0.01),

0.24998750000000042)

assert aproximado(integral2(0, 1, lambda x : x**4, 0.01),

0.19998333362500054)

assert aproximado(integral2(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01),

1.9999250000000026)

assert aproximado(log(2) - integral2(1, 2, lambda x : 1/x, 0.01),

3.124931644782336e-6)

assert aproximado(pi - 4 * integral2(0, 1, lambda x : 1/(x**2+1), 0.01),

-8.333333331389525e-6)

assert aproximado(integral3(0, 1, lambda x : x**3, 0.01),

0.24998750000000042)

assert aproximado(integral3(0, 1, lambda x : x**4, 0.01),

0.19998333362500054)

assert aproximado(integral3(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01),

1.9999250000000026)

assert aproximado(log(2) - integral3(1, 2, lambda x : 1/x, 0.01),

3.124931644782336e-6)

assert aproximado(pi - 4 * integral3(0, 1, lambda x : 1/(x**2+1), 0.01),

-8.333333331389525e-6)

print("Verificado")

# La verificación es

# >>> test_integral()

# Verificado

Factorizaciones de números de Hilbert

El teorema de Navidad de Fermat

Números de Pentanacci

Algoritmo de bajada para resolver un sistema triangular inferior

Integración por el método de los rectángulos

Entradas recientes

Buscar

Correo electrónico