PFH: La semana en Exercitium (14 de octubre de 2022)

Esta semana he publicado en Exercitium las soluciones de los siguientes problemas:

1. Suma de múltiplos de 3 ó 5
2. Puntos en el círculo
3. Aproximación del número e
4. Aproximación al límite de sen(x)/x cuando x tiende a cero
5. Cálculo del número π mediante la fórmula de Leibniz

A continuación se muestran las soluciones.

1. Suma de múltiplos de 3 ó 5

Definir la función

   euler1 :: Integer -> Integer

1	euler1 :: Integer -> Integer

tal que euler1 n es la suma de todos los múltiplos de 3 ó 5 menores que n. Por ejemplo,

   euler1 10      == 23
   euler1 (10^2)  == 2318
   euler1 (10^3)  == 233168
   euler1 (10^4)  == 23331668
   euler1 (10^5)  == 2333316668
   euler1 (10^10) == 23333333331666666668
   euler1 (10^20) == 2333333333333333333316666666666666666668
   euler1 (10^30) == 233333333333333333333333333333166666666666666666666666666668

euler1 10 == 23

euler1 (10^2) == 2318

euler1 (10^3) == 233168

euler1 (10^4) == 23331668

euler1 (10^5) == 2333316668

euler1 (10^10) == 23333333331666666668

euler1 (10^20) == 2333333333333333333316666666666666666668

euler1 (10^30) == 233333333333333333333333333333166666666666666666666666666668

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

Soluciones en Haskell

import Data.List (nub, union)
import qualified Data.Set as S (fromAscList, union)
import Test.QuickCheck

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

euler1a :: Integer -> Integer
euler1a n =
  sum [x | x <- [1..n-1],
           multiplo x 3 || multiplo x 5]

-- (multiplo x y) se verifica si x es un múltiplo de y. Por ejemplo.
--    multiplo 12 3  ==  True
--    multiplo 14 3  ==  False
multiplo :: Integer -> Integer -> Bool
multiplo x y = mod x y == 0

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

euler1b :: Integer -> Integer
euler1b n =
  sum [x | x <- [1..n-1],
           gcd x 15 > 1]

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

euler1c :: Integer -> Integer
euler1c n =
  sum [3,6..n-1] + sum [5,10..n-1] - sum [15,30..n-1]

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

euler1d :: Integer -> Integer
euler1d n =
  sum (nub ([3,6..n-1] ++ [5,10..n-1]))

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

euler1e :: Integer -> Integer
euler1e n =
  sum ([3,6..n-1] `union` [5,10..n-1])

-- 6ª solución                                                        --
-- ===========

euler1f :: Integer -> Integer
euler1f n =
  sum (S.fromAscList [3,6..n-1] `S.union` S.fromAscList [5,10..n-1])

-- 7ª solución                                                      --
-- ===========

euler1g :: Integer -> Integer
euler1g n =
  suma 3 n + suma 5 n - suma 15 n

-- (suma d x) es la suma de los múltiplos de d menores que x. Por
-- ejemplo,
--    suma 3 20  ==  63
suma :: Integer -> Integer -> Integer
suma d x = (a+b)*n `div` 2
    where a = d
          b = d * ((x-1) `div` d)
          n = 1 + (b-a) `div` d

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

-- La propiedad es
prop_euler1 :: Positive Integer -> Bool
prop_euler1 (Positive n) =
  all (== euler1a n)
      [euler1b n,
       euler1c n,
       euler1d n,
       euler1e n,
       euler1f n,
       euler1g n]

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

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

-- La comparación es
--    λ> euler1a (5*10^4)
--    583291668
--    (0.05 secs, 21,895,296 bytes)
--    λ> euler1b (5*10^4)
--    583291668
--    (0.05 secs, 26,055,096 bytes)
--    λ> euler1c (5*10^4)
--    583291668
--    (0.01 secs, 5,586,072 bytes)
--    λ> euler1d (5*10^4)
--    583291668
--    (2.83 secs, 7,922,304 bytes)
--    λ> euler1e (5*10^4)
--    583291668
--    (4.56 secs, 12,787,705,248 bytes)
--    λ> euler1f (5*10^4)
--    583291668
--    (0.01 secs, 8,168,584 bytes)
--    λ> euler1g (5*10^4)
--    583291668
--    (0.02 secs, 557,488 bytes)
--
--    λ> euler1a (3*10^6)
--    2099998500000
--    (2.72 secs, 1,282,255,816 bytes)
--    λ> euler1b (3*10^6)
--    2099998500000
--    (2.06 secs, 1,531,855,776 bytes)
--    λ> euler1c (3*10^6)
--    2099998500000
--    (0.38 secs, 305,127,480 bytes)
--    λ> euler1f (3*10^6)
--    2099998500000
--    (0.54 secs, 457,358,232 bytes)
--    λ> euler1g (3*10^6)
--    2099998500000
--    (0.01 secs, 560,472 bytes)
--
--    λ> euler1c (10^7)
--    23333331666668
--    (1.20 secs, 1,015,920,024 bytes)
--    λ> euler1f (10^7)
--    23333331666668
--    (2.00 secs, 1,523,225,648 bytes)
--    λ> euler1g (10^7)
--    23333331666668
--    (0.01 secs, 561,200 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

127

128

129

130

131

132

133

134

135

136

137

138

139

import Data.List (nub, union)

import qualified Data.Set as S (fromAscList, union)

import Test.QuickCheck

-- 1ª solución

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

euler1a :: Integer -> Integer

euler1a n =

sum [x | x <- [1..n-1],

multiplo x 3 || multiplo x 5]

-- (multiplo x y) se verifica si x es un múltiplo de y. Por ejemplo.

-- multiplo 12 3 == True

-- multiplo 14 3 == False

multiplo :: Integer -> Integer -> Bool

multiplo x y = mod x y == 0

-- 2ª solución --

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

euler1b :: Integer -> Integer

euler1b n =

sum [x | x <- [1..n-1],

gcd x 15 > 1]

-- 3ª solución --

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

euler1c :: Integer -> Integer

euler1c n =

sum [3,6..n-1] + sum [5,10..n-1] - sum [15,30..n-1]

-- 4ª solución --

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

euler1d :: Integer -> Integer

euler1d n =

sum (nub ([3,6..n-1] ++ [5,10..n-1]))

-- 5ª solución --

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

euler1e :: Integer -> Integer

euler1e n =

sum ([3,6..n-1] `union` [5,10..n-1])

-- 6ª solución --

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

euler1f :: Integer -> Integer

euler1f n =

sum (S.fromAscList [3,6..n-1] `S.union` S.fromAscList [5,10..n-1])

-- 7ª solución --

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

euler1g :: Integer -> Integer

euler1g n =

suma 3 n + suma 5 n - suma 15 n

-- (suma d x) es la suma de los múltiplos de d menores que x. Por

-- ejemplo,

-- suma 3 20 == 63

suma :: Integer -> Integer -> Integer

suma d x = (a+b)*n `div` 2

where a = d

b = d * ((x-1) `div` d)

n = 1 + (b-a) `div` d

-- Comprobación de equivalencia

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

-- La propiedad es

prop_euler1 :: Positive Integer -> Bool

prop_euler1 (Positive n) =

all (== euler1a n)

[euler1b n,

euler1c n,

euler1d n,

euler1e n,

euler1f n,

euler1g n]

-- La comprobación es

-- λ> quickCheck prop_euler1

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> euler1a (5*10^4)

-- 583291668

-- (0.05 secs, 21,895,296 bytes)

-- λ> euler1b (5*10^4)

-- 583291668

-- (0.05 secs, 26,055,096 bytes)

-- λ> euler1c (5*10^4)

-- 583291668

-- (0.01 secs, 5,586,072 bytes)

-- λ> euler1d (5*10^4)

-- 583291668

-- (2.83 secs, 7,922,304 bytes)

-- λ> euler1e (5*10^4)

-- 583291668

-- (4.56 secs, 12,787,705,248 bytes)

-- λ> euler1f (5*10^4)

-- 583291668

-- (0.01 secs, 8,168,584 bytes)

-- λ> euler1g (5*10^4)

-- 583291668

-- (0.02 secs, 557,488 bytes)

-- λ> euler1a (3*10^6)

-- 2099998500000

-- (2.72 secs, 1,282,255,816 bytes)

-- λ> euler1b (3*10^6)

-- 2099998500000

-- (2.06 secs, 1,531,855,776 bytes)

-- λ> euler1c (3*10^6)

-- 2099998500000

-- (0.38 secs, 305,127,480 bytes)

-- λ> euler1f (3*10^6)

-- 2099998500000

-- (0.54 secs, 457,358,232 bytes)

-- λ> euler1g (3*10^6)

-- 2099998500000

-- (0.01 secs, 560,472 bytes)

-- λ> euler1c (10^7)

-- 23333331666668

-- (1.20 secs, 1,015,920,024 bytes)

-- λ> euler1f (10^7)

-- 23333331666668

-- (2.00 secs, 1,523,225,648 bytes)

-- λ> euler1g (10^7)

-- 23333331666668

-- (0.01 secs, 561,200 bytes)

El código se encuentra en GitHub.

Soluciones en Python

from math import gcd
from timeit import Timer, default_timer

from hypothesis import given
from hypothesis import strategies as st

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

# multiplo(x, y) se verifica si x es un múltiplo de y. Por ejemplo.
#    multiplo(12, 3)  ==  True
#    multiplo(14, 3)  ==  False
def multiplo(x: int, y: int) -> int:
    return x % y == 0

def euler1a(n: int) -> int:
    return sum(x for x in range(1, n)
               if (multiplo(x, 3) or multiplo(x, 5)))

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

def euler1b(n: int) -> int:
    return sum(x for x in range(1, n)
               if gcd(x, 15) > 1)

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

def euler1c(n: int) -> int:
    return sum(range(3, n, 3)) + \
           sum(range(5, n, 5)) - \
           sum(range(15, n, 15))

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

def euler1d(n: int) -> int:
    return sum(set(list(range(3, n, 3)) + list(range(5, n, 5))))

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

def euler1e(n: int) -> int:
    return sum(set(list(range(3, n, 3))) | set(list(range(5, n, 5))))

# 6ª solución                                                        --
# ===========

# suma(d, x) es la suma de los múltiplos de d menores que x. Por
# ejemplo,
#    suma(3, 20)  ==  63
def suma(d: int, x: int) -> int:
    a = d
    b = d * ((x - 1) // d)
    n = 1 + (b - a) // d
    return (a + b) * n // 2

def euler1f(n: int) -> int:
    return suma(3, n) + suma(5, n) - suma(15, n)

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

# La propiedad es
@given(st.integers(min_value=1, max_value=1000))
def test_euler1(n: int) -> None:
    r = euler1a(n)
    assert euler1b(n) == r
    assert euler1c(n) == r
    assert euler1d(n) == r
    assert euler1e(n) == r

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

# 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('euler1a(10**7)')
#    1.49 segundos
#    >>> tiempo('euler1b(10**7)')
#    0.93 segundos
#    >>> tiempo('euler1c(10**7)')
#    0.07 segundos
#    >>> tiempo('euler1d(10**7)')
#    0.42 segundos
#    >>> tiempo('euler1e(10**7)')
#    0.69 segundos
#    >>> tiempo('euler1f(10**7)')
#    0.00 segundos
#
#    >>> tiempo('euler1c(10**8)')
#    0.72 segundos
#    >>> tiempo('euler1f(10**8)')
#    0.00 segundos

100

101

102

103

from math import gcd

from timeit import Timer, default_timer

from hypothesis import given

from hypothesis import strategies as st

# 1ª solución

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

# multiplo(x, y) se verifica si x es un múltiplo de y. Por ejemplo.

# multiplo(12, 3) == True

# multiplo(14, 3) == False

def multiplo(x: int, y: int) -> int:

return x % y == 0

def euler1a(n: int) -> int:

return sum(x for x in range(1, n)

if (multiplo(x, 3) or multiplo(x, 5)))

# 2ª solución --

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

def euler1b(n: int) -> int:

return sum(x for x in range(1, n)

if gcd(x, 15) > 1)

# 3ª solución --

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

def euler1c(n: int) -> int:

return sum(range(3, n, 3)) + \

sum(range(5, n, 5)) - \

sum(range(15, n, 15))

# 4ª solución --

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

def euler1d(n: int) -> int:

return sum(set(list(range(3, n, 3)) + list(range(5, n, 5))))

# 5ª solución --

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

def euler1e(n: int) -> int:

return sum(set(list(range(3, n, 3))) | set(list(range(5, n, 5))))

# 6ª solución --

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

# suma(d, x) es la suma de los múltiplos de d menores que x. Por

# ejemplo,

# suma(3, 20) == 63

def suma(d: int, x: int) -> int:

a = d

b = d * ((x - 1) // d)

n = 1 + (b - a) // d

return (a + b) * n // 2

def euler1f(n: int) -> int:

return suma(3, n) + suma(5, n) - suma(15, n)

# Comprobación de equivalencia

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

# La propiedad es

@given(st.integers(min_value=1, max_value=1000))

def test_euler1(n: int) -> None:

r = euler1a(n)

assert euler1b(n) == r

assert euler1c(n) == r

assert euler1d(n) == r

assert euler1e(n) == r

# La comprobación es

# src> poetry run pytest -q suma_de_multiplos_de_3_o_5.py

# 1 passed in 0.16s

# 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('euler1a(10**7)')

# 1.49 segundos

# >>> tiempo('euler1b(10**7)')

# 0.93 segundos

# >>> tiempo('euler1c(10**7)')

# 0.07 segundos

# >>> tiempo('euler1d(10**7)')

# 0.42 segundos

# >>> tiempo('euler1e(10**7)')

# 0.69 segundos

# >>> tiempo('euler1f(10**7)')

# 0.00 segundos

# >>> tiempo('euler1c(10**8)')

# 0.72 segundos

# >>> tiempo('euler1f(10**8)')

# 0.00 segundos

El código se encuentra en GitHub.

2. Puntos en el círculo

En el círculo de radio 2 hay 6 puntos cuyas coordenadas son puntos naturales:

   (0,0),(0,1),(0,2),(1,0),(1,1),(2,0)

1	(0,0),(0,1),(0,2),(1,0),(1,1),(2,0)

y en de radio 3 hay 11:

   (0,0),(0,1),(0,2),(0,3),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2),(3,0)

1	(0,0),(0,1),(0,2),(0,3),(1,0),(1,1),(1,2),(2,0),(2,1),(2,2),(3,0)

Definir la función

   circulo :: Int -> Int

1	circulo :: Int -> Int

tal que circulo n es el la cantidad de pares de números naturales (x,y) que se encuentran en el círculo de radio n. Por ejemplo,

   circulo 1    ==  3
   circulo 2    ==  6
   circulo 3    ==  11
   circulo 4    ==  17
   circulo 100  ==  7955

circulo 1 == 3

circulo 2 == 6

circulo 3 == 11

circulo 4 == 17

circulo 100 == 7955

Soluciones en Haskell

import Test.QuickCheck

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

circulo1 :: Int -> Int
circulo1 n = length (enCirculo1 n)

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

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

circulo2 :: Int -> Int
circulo2 0 = 1
circulo2 n =
  2 * length (enSemiCirculo n) + ceiling(fromIntegral n / sqrt 2)

enSemiCirculo :: Int -> [(Int, Int)]
enSemiCirculo n =
  [(x,y) | x <- [0..floor (sqrt (fromIntegral (n * n)))],
           y <- [x+1..truncate (sqrt (fromIntegral (n*n - x*x)))]]

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

-- La propiedad es
prop_circulo :: Positive Int -> Bool
prop_circulo (Positive n) =
  circulo1 n == circulo2 n

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

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

-- La comparación es
--    λ> circulo1 (2*10^3)
--    3143587
--    (3.58 secs, 1,744,162,600 bytes)
--    λ> circulo2 (2*10^3)
--    3143587
--    (0.41 secs, 266,374,208 bytes)

import Test.QuickCheck

-- 1ª solución

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

circulo1 :: Int -> Int

circulo1 n = length (enCirculo1 n)

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

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

y <- [0..n],

x*x+y*y <= n*n]

-- 2ª solución

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

circulo2 :: Int -> Int

circulo2 0 = 1

circulo2 n =

2 * length (enSemiCirculo n) + ceiling(fromIntegral n / sqrt 2)

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

enSemiCirculo n =

[(x,y) | x <- [0..floor (sqrt (fromIntegral (n * n)))],

y <- [x+1..truncate (sqrt (fromIntegral (n*n - x*x)))]]

-- Comprobación de equivalencia

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

-- La propiedad es

prop_circulo :: Positive Int -> Bool

prop_circulo (Positive n) =

circulo1 n == circulo2 n

-- La comprobación es

-- λ> quickCheck prop_circulo

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> circulo1 (2*10^3)

-- 3143587

-- (3.58 secs, 1,744,162,600 bytes)

-- λ> circulo2 (2*10^3)

-- 3143587

-- (0.41 secs, 266,374,208 bytes)

El código se encuentra en GitHub.

Soluciones en Python

from math import sqrt, trunc, ceil
from timeit import Timer, default_timer

from hypothesis import given
from hypothesis import strategies as st

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

def circulo1(n: int) -> int:
    return len([(x, y)
                for x in range(0, n + 1)
                for y in range(0, n + 1)
                if x * x + y * y <= n * n])

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

def enSemiCirculo(n: int) -> list[tuple[int, int]]:
    return [(x, y)
            for x in range(0, ceil(sqrt(n**2)) + 1)
            for y in range(x+1, trunc(sqrt(n**2 - x**2)) + 1)]

def circulo2(n: int) -> int:
    if n == 0:
        return 1
    return (2 * len(enSemiCirculo(n)) + ceil(n / sqrt(2)))

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

def circulo3(n: int) -> int:
    r = 0
    for x in range(0, n + 1):
        for y in range(0, n + 1):
            if x**2 + y**2 <= n**2:
                r = r + 1
    return r

# 4ª solución
# ===========

def circulo4(n: int) -> int:
    r = 0
    for x in range(0, ceil(sqrt(n**2)) + 1):
        for y in range(x + 1, trunc(sqrt(n**2 - x**2)) + 1):
            if x**2 + y**2 <= n**2:
                r = r + 1
    return 2 * r + ceil(n / sqrt(2))

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

# La propiedad es
@given(st.integers(min_value=1, max_value=100))
def test_circulo(n: int) -> None:
    r = circulo1(n)
    assert circulo2(n) == r
    assert circulo3(n) == r
    assert circulo4(n) == r

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

# 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('circulo1(2000)')
#    0.71 segundos
#    >>> tiempo('circulo2(2000)')
#    0.76 segundos
#    >>> tiempo('circulo3(2000)')
#    2.63 segundos
#    >>> tiempo('circulo4(2000)')
#    1.06 segundos

from math import sqrt, trunc, ceil

from timeit import Timer, default_timer

from hypothesis import given

from hypothesis import strategies as st

# 1ª solución

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

def circulo1(n: int) -> int:

return len([(x, y)

for x in range(0, n + 1)

for y in range(0, n + 1)

if x * x + y * y <= n * n])

# 2ª solución

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

def enSemiCirculo(n: int) -> list[tuple[int, int]]:

return [(x, y)

for x in range(0, ceil(sqrt(n**2)) + 1)

for y in range(x+1, trunc(sqrt(n**2 - x**2)) + 1)]

def circulo2(n: int) -> int:

if n == 0:

return 1

return (2 * len(enSemiCirculo(n)) + ceil(n / sqrt(2)))

# 3ª solución

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

def circulo3(n: int) -> int:

r = 0

for x in range(0, n + 1):

for y in range(0, n + 1):

if x**2 + y**2 <= n**2:

r = r + 1

return r

# 4ª solución

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

def circulo4(n: int) -> int:

r = 0

for x in range(0, ceil(sqrt(n**2)) + 1):

for y in range(x + 1, trunc(sqrt(n**2 - x**2)) + 1):

if x**2 + y**2 <= n**2:

r = r + 1

return 2 * r + ceil(n / sqrt(2))

# Comprobación de equivalencia

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

# La propiedad es

@given(st.integers(min_value=1, max_value=100))

def test_circulo(n: int) -> None:

r = circulo1(n)

assert circulo2(n) == r

assert circulo3(n) == r

assert circulo4(n) == r

# La comprobación es

# src> poetry run pytest -q puntos_dentro_del_circulo.py

# 1 passed in 0.60s

# 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('circulo1(2000)')

# 0.71 segundos

# >>> tiempo('circulo2(2000)')

# 0.76 segundos

# >>> tiempo('circulo3(2000)')

# 2.63 segundos

# >>> tiempo('circulo4(2000)')

# 1.06 segundos

El código se encuentra en GitHub.

3. Aproximación del número e

El número e se define como el límite de la sucesión

$\left(1+\dfrac{1}{n}\right)^n$ .

Definir las funciones

   aproxE      :: Int -> [Double]
   errorAproxE :: Double -> Int

1 2	aproxE :: Int -> [Double] errorAproxE :: Double -> Int

tales que

aproxE k es la lista de los k primeros términos de la sucesión. Por ejemplo,

     aproxE 4 == [2.0,2.25,2.37037037037037,2.44140625]
     last (aproxE (7*10^7))  ==  2.7182818287372563

1 2	aproxE 4 == [2.0,2.25,2.37037037037037,2.44140625] last (aproxE (7*10^7)) == 2.7182818287372563

errorE x es el menor número de términos de la sucesión necesarios para obtener su límite con un error menor que x. Por ejemplo,

     errorAproxE 0.1    ==  13
     errorAproxE 0.01   ==  135
     errorAproxE 0.001  ==  1359

errorAproxE 0.1 == 13

errorAproxE 0.01 == 135

errorAproxE 0.001 == 1359

Soluciones en Haskell

import Test.QuickCheck

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

aproxE1 :: Int -> [Double]
aproxE1 k = [(1+1/n)**n | n <- [1..k']]
  where k' = fromIntegral k

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

aproxE2 :: Int -> [Double]
aproxE2 0 = []
aproxE2 n = aproxE2 (n-1) ++ [(1+1/n')**n']
  where n' = fromIntegral n

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

aproxE3 :: Int -> [Double]
aproxE3 = reverse . aux . fromIntegral
  where aux 0 = []
        aux n = (1+1/n)**n : aux (n-1)

-- 4ª definición de aproxE
-- =======================

aproxE4 :: Int -> [Double]
aproxE4 k = aux [] (fromIntegral k)
  where aux xs 0 = xs
        aux xs n = aux ((1+1/n)**n : xs) (n-1)

-- 5ª definición de aproxE
-- =======================

aproxE5 :: Int -> [Double]
aproxE5 k = map (\ n -> (1+1/n)**n) [1..k']
  where k' = fromIntegral k

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

-- La propiedad es
prop_aproxE :: Positive Int -> Bool
prop_aproxE (Positive k) =
  all (== aproxE1 k)
      [aproxE2 k,
       aproxE3 k,
       aproxE4 k,
       aproxE5 k]

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

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

-- La comparación es
--    λ> last (aproxE1 (2*10^4))
--    2.718213874533619
--    (0.04 secs, 5,368,968 bytes)
--    λ> last (aproxE2 (2*10^4))
--    2.718213874533619
--    (5.93 secs, 17,514,767,104 bytes)
--    λ> last (aproxE3 (2*10^4))
--    2.718213874533619
--    (0.05 secs, 9,529,336 bytes)
--    λ> last (aproxE4 (2*10^4))
--    2.718213874533619
--    (0.05 secs, 9,529,184 bytes)
--    λ> last (aproxE5 (2*10^4))
--    2.718213874533619
--    (0.01 secs, 4,888,960 bytes)
--
--    λ> last (aproxE1 (2*10^6))
--    2.7182811492688552
--    (0.54 secs, 480,570,120 bytes)
--    λ> last (aproxE3 (2*10^6))
--    2.7182811492688552
--    (2.07 secs, 896,570,280 bytes)
--    λ> last (aproxE4 (2*10^6))
--    2.7182811492688552
--    (2.18 secs, 896,570,336 bytes)
--    λ> last (aproxE5 (2*10^6))
--    2.7182811492688552
--    (0.09 secs, 432,570,112 bytes)

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

errorAproxE1 :: Double -> Int
errorAproxE1 x =
  round (head [n | n <- [1..], abs (exp 1 - (1+1/n)**n) < x])

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

errorAproxE2 :: Double -> Int
errorAproxE2 x = aux 1
  where aux n | abs (exp 1 - (1+1/n)**n) < x = round n
              | otherwise                    = aux (n+1)

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

errorAproxE3 :: Double -> Int
errorAproxE3 x =
  round (head (dropWhile (\ n -> abs (exp 1 - (1+1/n)**n) >= x) [1..]))

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

-- La propiedad es
prop_errorAproxE :: Positive Double -> Bool
prop_errorAproxE (Positive x) =
  all (== errorAproxE1 x)
      [errorAproxE2 x,
       errorAproxE3 x]

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

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

-- La comparación es
--    λ> errorAproxE1 0.000001
--    1358611
--    (1.70 secs, 674,424,552 bytes)
--    λ> errorAproxE2 0.000001
--    1358611
--    (1.79 secs, 739,637,704 bytes)
--    λ> errorAproxE3 0.000001
--    1358611
--    (1.20 secs, 609,211,144 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

127

128

129

130

131

132

133

134

135

136

137

138

import Test.QuickCheck

-- 1ª definición de aproxE

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

aproxE1 :: Int -> [Double]

aproxE1 k = [(1+1/n)**n | n <- [1..k']]

where k' = fromIntegral k

-- 2ª definición de aproxE

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

aproxE2 :: Int -> [Double]

aproxE2 0 = []

aproxE2 n = aproxE2 (n-1) ++ [(1+1/n')**n']

where n' = fromIntegral n

-- 3ª definición de aproxE

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

aproxE3 :: Int -> [Double]

aproxE3 = reverse . aux . fromIntegral

where aux 0 = []

aux n = (1+1/n)**n : aux (n-1)

-- 4ª definición de aproxE

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

aproxE4 :: Int -> [Double]

aproxE4 k = aux [] (fromIntegral k)

where aux xs 0 = xs

aux xs n = aux ((1+1/n)**n : xs) (n-1)

-- 5ª definición de aproxE

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

aproxE5 :: Int -> [Double]

aproxE5 k = map (\ n -> (1+1/n)**n) [1..k']

where k' = fromIntegral k

-- Comprobación de equivalencia de aproxE

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

-- La propiedad es

prop_aproxE :: Positive Int -> Bool

prop_aproxE (Positive k) =

all (== aproxE1 k)

[aproxE2 k,

aproxE3 k,

aproxE4 k,

aproxE5 k]

-- La comprobación es

-- λ> quickCheck prop_aproxE

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de aproxE

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

-- La comparación es

-- λ> last (aproxE1 (2*10^4))

-- 2.718213874533619

-- (0.04 secs, 5,368,968 bytes)

-- λ> last (aproxE2 (2*10^4))

-- 2.718213874533619

-- (5.93 secs, 17,514,767,104 bytes)

-- λ> last (aproxE3 (2*10^4))

-- 2.718213874533619

-- (0.05 secs, 9,529,336 bytes)

-- λ> last (aproxE4 (2*10^4))

-- 2.718213874533619

-- (0.05 secs, 9,529,184 bytes)

-- λ> last (aproxE5 (2*10^4))

-- 2.718213874533619

-- (0.01 secs, 4,888,960 bytes)

-- λ> last (aproxE1 (2*10^6))

-- 2.7182811492688552

-- (0.54 secs, 480,570,120 bytes)

-- λ> last (aproxE3 (2*10^6))

-- 2.7182811492688552

-- (2.07 secs, 896,570,280 bytes)

-- λ> last (aproxE4 (2*10^6))

-- 2.7182811492688552

-- (2.18 secs, 896,570,336 bytes)

-- λ> last (aproxE5 (2*10^6))

-- 2.7182811492688552

-- (0.09 secs, 432,570,112 bytes)

-- 1ª definición de errorAproxE

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

errorAproxE1 :: Double -> Int

errorAproxE1 x =

round (head [n | n <- [1..], abs (exp 1 - (1+1/n)**n) < x])

-- 2ª definición de errorAproxE

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

errorAproxE2 :: Double -> Int

errorAproxE2 x = aux 1

where aux n | abs (exp 1 - (1+1/n)**n) < x = round n

| otherwise = aux (n+1)

-- 3ª definición de errorAproxE

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

errorAproxE3 :: Double -> Int

errorAproxE3 x =

round (head (dropWhile (\ n -> abs (exp 1 - (1+1/n)**n) >= x) [1..]))

-- Comprobación de equivalencia de errorAproxE

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

-- La propiedad es

prop_errorAproxE :: Positive Double -> Bool

prop_errorAproxE (Positive x) =

all (== errorAproxE1 x)

[errorAproxE2 x,

errorAproxE3 x]

-- La comprobación es

-- λ> quickCheck prop_errorAproxE

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de errorAproxE

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

-- La comparación es

-- λ> errorAproxE1 0.000001

-- 1358611

-- (1.70 secs, 674,424,552 bytes)

-- λ> errorAproxE2 0.000001

-- 1358611

-- (1.79 secs, 739,637,704 bytes)

-- λ> errorAproxE3 0.000001

-- 1358611

-- (1.20 secs, 609,211,144 bytes)

El código se encuentra en GitHub.

Soluciones en Python

from itertools import dropwhile, islice
from math import e
from sys import setrecursionlimit
from timeit import Timer, default_timer
from typing import Iterator

from hypothesis import given
from hypothesis import strategies as st

setrecursionlimit(10**6)

# 1ª definición de aproxE
# =======================

def aproxE1(k: int) -> list[float]:
    return [(1 + 1/n)**n for n in range(1, k + 1)]

# 2ª definición de aproxE
# =======================

def aproxE2(n: int) -> list[float]:
    if n == 0:
        return []
    return aproxE2(n - 1) + [(1 + 1/n)**n]

# 3ª definición de aproxE
# =======================

def aproxE3(n: int) -> list[float]:
    def aux(n: int) -> list[float]:
        if n == 0:
            return []
        return [(1 + 1/n)**n] + aux(n - 1)

    return list(reversed(aux(n)))

# 4ª definición de aproxE
# =======================

def aproxE4(n: int) -> list[float]:
    def aux(xs: list[float], n: int) -> list[float]:
        if n == 0:
            return xs
        return aux([(1 + 1/n)**n] + xs, n - 1)

    return aux([], n)

# 5ª definición de aproxE
# =======================

def aproxE5(n: int) -> list[float]:
    return list(map((lambda k: (1+1/k)**k), range(1, n+1)))

# 6ª definición de aproxE
# =======================

def aproxE6(n: int) -> list[float]:
    r = []
    for k in range(1, n+1):
        r.append((1+1/k)**k)
    return r

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

# La propiedad es
@given(st.integers(min_value=1, max_value=100))
def test_aproxE(n: int) -> None:
    r = aproxE1(n)
    assert aproxE2(n) == r
    assert aproxE3(n) == r
    assert aproxE4(n) == r
    assert aproxE5(n) == r
    assert aproxE6(n) == r

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

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

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

# La comparación es
#    >>> tiempo('aproxE1(20000)')
#    0.00 segundos
#    >>> tiempo('aproxE2(20000)')
#    0.43 segundos
#    >>> tiempo('aproxE3(20000)')
#    0.60 segundos
#    >>> tiempo('aproxE4(20000)')
#    1.23 segundos
#    >>> tiempo('aproxE5(20000)')
#    0.00 segundos
#    >>> tiempo('aproxE6(20000)')
#    0.00 segundos

#    >>> tiempo('aproxE1(10**7)')
#    1.18 segundos
#    >>> tiempo('aproxE5(10**7)')
#    1.48 segundos
#    >>> tiempo('aproxE6(10**7)')
#    1.43 segundos

# 1ª definición de errorAproxE
# ============================

# naturales es el generador de los números naturales positivos, Por
# ejemplo,
#    >>> list(islice(naturales(), 5))
#    [1, 2, 3, 4, 5]
def naturales() -> Iterator[int]:
    i = 1
    while True:
        yield i
        i += 1

def errorAproxE1(x: float) -> int:
    return list(islice((n for n in naturales()
                        if abs(e - (1 + 1/n)**n) < x), 1))[0]

# # 2ª definición de errorAproxE
# # ============================

def errorAproxE2(x: float) -> int:
    def aux(n: int) -> int:
        if abs(e - (1 + 1/n)**n) < x:
            return n
        return aux(n + 1)

    return aux(1)

# 3ª definición de errorAproxE
# ============================

def errorAproxE3(x: float) -> int:
    return list(islice(dropwhile(lambda n: abs(e - (1 + 1/n)**n) >= x,
                                 naturales()),
                       1))[0]


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

@given(st.integers(min_value=1, max_value=100))
def test_errorAproxE(n: int) -> None:
    r = errorAproxE1(n)
    assert errorAproxE2(n) == r
    assert errorAproxE3(n) == r

# La comprobación es
#    src> poetry run pytest -q aproximacion_del_numero_e.py
#    2 passed in 0.60s

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

# La comparación es
#    >>> tiempo('errorAproxE1(0.0001)')
#    0.00 segundos
#    >>> tiempo('errorAproxE2(0.0001)')
#    0.00 segundos
#    >>> tiempo('errorAproxE3(0.0001)')
#    0.00 segundos
#
#    >>> tiempo('errorAproxE1(0.0000001)')
#    2.48 segundos
#    >>> tiempo('errorAproxE3(0.0000001)')
#    2.61 segundos

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

169

170

171

172

173

from itertools import dropwhile, islice

from math import e

from sys import setrecursionlimit

from timeit import Timer, default_timer

from typing import Iterator

from hypothesis import given

from hypothesis import strategies as st

setrecursionlimit(10**6)

# 1ª definición de aproxE

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

def aproxE1(k: int) -> list[float]:

return [(1 + 1/n)**n for n in range(1, k + 1)]

# 2ª definición de aproxE

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

def aproxE2(n: int) -> list[float]:

if n == 0:

return []

return aproxE2(n - 1) + [(1 + 1/n)**n]

# 3ª definición de aproxE

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

def aproxE3(n: int) -> list[float]:

def aux(n: int) -> list[float]:

if n == 0:

return []

return [(1 + 1/n)**n] + aux(n - 1)

return list(reversed(aux(n)))

# 4ª definición de aproxE

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

def aproxE4(n: int) -> list[float]:

def aux(xs: list[float], n: int) -> list[float]:

if n == 0:

return xs

return aux([(1 + 1/n)**n] + xs, n - 1)

return aux([], n)

# 5ª definición de aproxE

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

def aproxE5(n: int) -> list[float]:

return list(map((lambda k: (1+1/k)**k), range(1, n+1)))

# 6ª definición de aproxE

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

def aproxE6(n: int) -> list[float]:

r = []

for k in range(1, n+1):

r.append((1+1/k)**k)

return r

# Comprobación de equivalencia de aproxE

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

# La propiedad es

@given(st.integers(min_value=1, max_value=100))

def test_aproxE(n: int) -> None:

r = aproxE1(n)

assert aproxE2(n) == r

assert aproxE3(n) == r

assert aproxE4(n) == r

assert aproxE5(n) == r

assert aproxE6(n) == r

# La comprobación es

# src> poetry run pytest -q aproximacion_del_numero_e.py

# 1 passed in 0.60s

# Comparación de eficiencia de aproxE

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

def tiempo(ex: str) -> None:

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

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

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

# La comparación es

# >>> tiempo('aproxE1(20000)')

# 0.00 segundos

# >>> tiempo('aproxE2(20000)')

# 0.43 segundos

# >>> tiempo('aproxE3(20000)')

# 0.60 segundos

# >>> tiempo('aproxE4(20000)')

# 1.23 segundos

# >>> tiempo('aproxE5(20000)')

# 0.00 segundos

# >>> tiempo('aproxE6(20000)')

# 0.00 segundos

# >>> tiempo('aproxE1(10**7)')

# 1.18 segundos

# >>> tiempo('aproxE5(10**7)')

# 1.48 segundos

# >>> tiempo('aproxE6(10**7)')

# 1.43 segundos

# 1ª definición de errorAproxE

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

# naturales es el generador de los números naturales positivos, Por

# ejemplo,

# >>> list(islice(naturales(), 5))

# [1, 2, 3, 4, 5]

def naturales() -> Iterator[int]:

i = 1

while True:

yield i

i += 1

def errorAproxE1(x: float) -> int:

return list(islice((n for n in naturales()

if abs(e - (1 + 1/n)**n) < x), 1))[0]

# # 2ª definición de errorAproxE

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

def errorAproxE2(x: float) -> int:

def aux(n: int) -> int:

if abs(e - (1 + 1/n)**n) < x:

return n

return aux(n + 1)

return aux(1)

# 3ª definición de errorAproxE

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

def errorAproxE3(x: float) -> int:

return list(islice(dropwhile(lambda n: abs(e - (1 + 1/n)**n) >= x,

naturales()),

1))[0]

# Comprobación de equivalencia de errorAproxE

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

@given(st.integers(min_value=1, max_value=100))

def test_errorAproxE(n: int) -> None:

r = errorAproxE1(n)

assert errorAproxE2(n) == r

assert errorAproxE3(n) == r

# La comprobación es

# src> poetry run pytest -q aproximacion_del_numero_e.py

# 2 passed in 0.60s

# Comparación de eficiencia de aproxE

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

# La comparación es

# >>> tiempo('errorAproxE1(0.0001)')

# 0.00 segundos

# >>> tiempo('errorAproxE2(0.0001)')

# 0.00 segundos

# >>> tiempo('errorAproxE3(0.0001)')

# 0.00 segundos

# >>> tiempo('errorAproxE1(0.0000001)')

# 2.48 segundos

# >>> tiempo('errorAproxE3(0.0000001)')

# 2.61 segundos

El código se encuentra en GitHub.

4. Aproximación al límite de sen(x)/x cuando x tiende a cero

El limite de sen(x)/x, cuando x tiende a cero, se puede calcular como el límite de la sucesión sen(1/n)/(1/n), cuando n tiende a infinito.

Definir las funciones

   aproxLimSeno :: Int -> [Double]
   errorLimSeno :: Double -> Int

1 2	aproxLimSeno :: Int -> [Double] errorLimSeno :: Double -> Int

tales que

aproxLimSeno n es la lista de los n primeros términos de la sucesión sen(1/m)/(1/m). Por ejemplo,

     aproxLimSeno 1 == [0.8414709848078965]
     aproxLimSeno 2 == [0.8414709848078965,0.958851077208406]

1 2	aproxLimSeno 1 == [0.8414709848078965] aproxLimSeno 2 == [0.8414709848078965,0.958851077208406]

errorLimSeno x es el menor número de términos de la sucesión sen(1/m)/(1/m) necesarios para obtener su límite con un error menor que x. Por ejemplo,

     errorLimSeno 0.1     ==   2
     errorLimSeno 0.01    ==   5
     errorLimSeno 0.001   ==  13
     errorLimSeno 0.0001  ==  41

errorLimSeno 0.1 == 2

errorLimSeno 0.01 == 5

errorLimSeno 0.001 == 13

errorLimSeno 0.0001 == 41

Soluciones en Haskell

import Test.QuickCheck

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

aproxLimSeno1 :: Int -> [Double]
aproxLimSeno1 k = [sin(1/n)/(1/n) | n <- [1..k']]
  where k' = fromIntegral k

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

aproxLimSeno2 :: Int -> [Double]
aproxLimSeno2 0 = []
aproxLimSeno2 n = aproxLimSeno2 (n-1) ++ [sin(1/n')/(1/n')]
  where n' = fromIntegral n

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

aproxLimSeno3 :: Int -> [Double]
aproxLimSeno3 = reverse . aux . fromIntegral
  where aux 0 = []
        aux n = sin(1/n)/(1/n) : aux (n-1)

-- 4ª definición de aproxLimSeno
-- =============================

aproxLimSeno4 :: Int -> [Double]
aproxLimSeno4 k = aux [] (fromIntegral k)
  where aux xs 0 = xs
        aux xs n = aux (sin(1/n)/(1/n) : xs) (n-1)

-- 5ª definición de aproxLimSeno
-- =============================

aproxLimSeno5 :: Int -> [Double]
aproxLimSeno5 k =  map (\ n -> sin(1/n)/(1/n)) [1..k']
  where k' = fromIntegral k

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

-- La propiedad es
prop_aproxLimSeno :: Positive Int -> Bool
prop_aproxLimSeno (Positive k) =
  all (== aproxLimSeno1 k)
      [aproxLimSeno2 k,
       aproxLimSeno3 k,
       aproxLimSeno4 k,
       aproxLimSeno5 k]

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

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

-- La comparación es
--    λ> last (aproxLimSeno1 (2*10^4))
--    0.9999999995833334
--    (0.01 secs, 5,415,816 bytes)
--    λ> last (aproxLimSeno2 (2*10^4))
--    0.9999999995833334
--    (4.48 secs, 17,514,768,064 bytes)
--    λ> last (aproxLimSeno3 (2*10^4))
--    0.9999999995833334
--    (0.02 secs, 9,530,120 bytes)
--    λ> last (aproxLimSeno4 (2*10^4))
--    0.9999999995833334
--    (0.02 secs, 9,529,968 bytes)
--    λ> last (aproxLimSeno5 (2*10^4))
--    0.9999999995833334
--    (0.01 secs, 4,889,720 bytes)
--
--    λ> last (aproxLimSeno1 (2*10^6))
--    0.9999999999999583
--    (0.46 secs, 480,569,808 bytes)
--    λ> last (aproxLimSeno3 (2*10^6))
--    0.9999999999999583
--    (1.96 secs, 896,569,992 bytes)
--    λ> last (aproxLimSeno4 (2*10^6))
--    0.9999999999999583
--    (1.93 secs, 896,570,048 bytes)
--    λ> last (aproxLimSeno5 (2*10^6))
--    0.9999999999999583
--    (0.05 secs, 432,569,800 bytes)
--
--    λ> last (aproxLimSeno1 (10^7))
--    0.9999999999999983
--    (2.26 secs, 2,400,569,760 bytes)
--    λ> last (aproxLimSeno5 (10^7))
--    0.9999999999999983
--    (0.24 secs, 2,160,569,752 bytes)

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

errorLimSeno1 :: Double -> Int
errorLimSeno1 x =
  round (head [m | m <- [1..], abs (1 - sin(1/m)/(1/m)) < x])

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

errorLimSeno2 :: Double -> Int
errorLimSeno2 x = aux 1
  where aux n | abs (1 - sin(1/n)/(1/n)) < x = round n
              | otherwise                    = aux (n+1)

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

errorLimSeno3 :: Double -> Int
errorLimSeno3 x =
  round (head (dropWhile (\ n -> abs (1 - sin(1/n)/(1/n)) >= x) [1..]))

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

-- La propiedad es
prop_errorLimSeno :: Positive Double -> Bool
prop_errorLimSeno (Positive x) =
  all (== errorLimSeno1 x)
      [errorLimSeno2 x,
       errorLimSeno3 x]

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

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

-- La comparación es
--    λ> errorLimSeno1 (10**(-12))
--    408230
--    (0.41 secs, 206,300,808 bytes)
--    λ> errorLimSeno2 (10**(-12))
--    408230
--    (0.46 secs, 225,895,672 bytes)
--    λ> errorLimSeno3 (10**(-12))
--    408230
--    (0.37 secs, 186,705,688 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

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

import Test.QuickCheck

-- 1ª definición de aproxLimSeno

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

aproxLimSeno1 :: Int -> [Double]

aproxLimSeno1 k = [sin(1/n)/(1/n) | n <- [1..k']]

where k' = fromIntegral k

-- 2ª definición de aproxLimSeno

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

aproxLimSeno2 :: Int -> [Double]

aproxLimSeno2 0 = []

aproxLimSeno2 n = aproxLimSeno2 (n-1) ++ [sin(1/n')/(1/n')]

where n' = fromIntegral n

-- 3ª definición de aproxLimSeno

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

aproxLimSeno3 :: Int -> [Double]

aproxLimSeno3 = reverse . aux . fromIntegral

where aux 0 = []

aux n = sin(1/n)/(1/n) : aux (n-1)

-- 4ª definición de aproxLimSeno

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

aproxLimSeno4 :: Int -> [Double]

aproxLimSeno4 k = aux [] (fromIntegral k)

where aux xs 0 = xs

aux xs n = aux (sin(1/n)/(1/n) : xs) (n-1)

-- 5ª definición de aproxLimSeno

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

aproxLimSeno5 :: Int -> [Double]

aproxLimSeno5 k = map (\ n -> sin(1/n)/(1/n)) [1..k']

where k' = fromIntegral k

-- Comprobación de equivalencia de aproxLimSeno

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

-- La propiedad es

prop_aproxLimSeno :: Positive Int -> Bool

prop_aproxLimSeno (Positive k) =

all (== aproxLimSeno1 k)

[aproxLimSeno2 k,

aproxLimSeno3 k,

aproxLimSeno4 k,

aproxLimSeno5 k]

-- La comprobación es

-- λ> quickCheck prop_aproxLimSeno

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de aproxLimSeno

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

-- La comparación es

-- λ> last (aproxLimSeno1 (2*10^4))

-- 0.9999999995833334

-- (0.01 secs, 5,415,816 bytes)

-- λ> last (aproxLimSeno2 (2*10^4))

-- 0.9999999995833334

-- (4.48 secs, 17,514,768,064 bytes)

-- λ> last (aproxLimSeno3 (2*10^4))

-- 0.9999999995833334

-- (0.02 secs, 9,530,120 bytes)

-- λ> last (aproxLimSeno4 (2*10^4))

-- 0.9999999995833334

-- (0.02 secs, 9,529,968 bytes)

-- λ> last (aproxLimSeno5 (2*10^4))

-- 0.9999999995833334

-- (0.01 secs, 4,889,720 bytes)

-- λ> last (aproxLimSeno1 (2*10^6))

-- 0.9999999999999583

-- (0.46 secs, 480,569,808 bytes)

-- λ> last (aproxLimSeno3 (2*10^6))

-- 0.9999999999999583

-- (1.96 secs, 896,569,992 bytes)

-- λ> last (aproxLimSeno4 (2*10^6))

-- 0.9999999999999583

-- (1.93 secs, 896,570,048 bytes)

-- λ> last (aproxLimSeno5 (2*10^6))

-- 0.9999999999999583

-- (0.05 secs, 432,569,800 bytes)

-- λ> last (aproxLimSeno1 (10^7))

-- 0.9999999999999983

-- (2.26 secs, 2,400,569,760 bytes)

-- λ> last (aproxLimSeno5 (10^7))

-- 0.9999999999999983

-- (0.24 secs, 2,160,569,752 bytes)

-- 1ª definición de errorLimSeno

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

errorLimSeno1 :: Double -> Int

errorLimSeno1 x =

round (head [m | m <- [1..], abs (1 - sin(1/m)/(1/m)) < x])

-- 2ª definición de errorLimSeno

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

errorLimSeno2 :: Double -> Int

errorLimSeno2 x = aux 1

where aux n | abs (1 - sin(1/n)/(1/n)) < x = round n

| otherwise = aux (n+1)

-- 3ª definición de errorLimSeno

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

errorLimSeno3 :: Double -> Int

errorLimSeno3 x =

round (head (dropWhile (\ n -> abs (1 - sin(1/n)/(1/n)) >= x) [1..]))

-- Comprobación de equivalencia de errorLimSeno

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

-- La propiedad es

prop_errorLimSeno :: Positive Double -> Bool

prop_errorLimSeno (Positive x) =

all (== errorLimSeno1 x)

[errorLimSeno2 x,

errorLimSeno3 x]

-- La comprobación es

-- λ> quickCheck prop_errorLimSeno

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de errorLimSeno

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

-- La comparación es

-- λ> errorLimSeno1 (10**(-12))

-- 408230

-- (0.41 secs, 206,300,808 bytes)

-- λ> errorLimSeno2 (10**(-12))

-- 408230

-- (0.46 secs, 225,895,672 bytes)

-- λ> errorLimSeno3 (10**(-12))

-- 408230

-- (0.37 secs, 186,705,688 bytes)

El código se encuentra en GitHub.

Soluciones en Python

from itertools import dropwhile, islice
from math import sin
from sys import setrecursionlimit
from timeit import Timer, default_timer
from typing import Iterator

from hypothesis import given
from hypothesis import strategies as st

setrecursionlimit(10**6)

# 1ª definición de aproxLimSeno
# =============================

def aproxLimSeno1(k: int) -> list[float]:
    return [sin(1/n)/(1/n) for n in range(1, k + 1)]

# 2ª definición de aproxLimSeno
# =============================

def aproxLimSeno2(n: int) -> list[float]:
    if n == 0:
        return []
    return aproxLimSeno2(n - 1) + [sin(1/n)/(1/n)]

# 3ª definición de aproxLimSeno
# =============================

def aproxLimSeno3(n: int) -> list[float]:
    def aux(n: int) -> list[float]:
        if n == 0:
            return []
        return [sin(1/n)/(1/n)] + aux(n - 1)

    return list(reversed(aux(n)))

# 4ª definición de aproxLimSeno
# =============================

def aproxLimSeno4(n: int) -> list[float]:
    def aux(xs: list[float], n: int) -> list[float]:
        if n == 0:
            return xs
        return aux([sin(1/n)/(1/n)] + xs, n - 1)

    return aux([], n)

# 5ª definición de aproxLimSeno
# =============================

def aproxLimSeno5(n: int) -> list[float]:
    return list(map((lambda k: sin(1/k)/(1/k)), range(1, n+1)))

# 6ª definición de aproxLimSeno
# =============================

def aproxLimSeno6(n: int) -> list[float]:
    r = []
    for k in range(1, n+1):
        r.append(sin(1/k)/(1/k))
    return r

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

# La propiedad es
@given(st.integers(min_value=1, max_value=100))
def test_aproxLimSeno(n: int) -> None:
    r = aproxLimSeno1(n)
    assert aproxLimSeno2(n) == r
    assert aproxLimSeno3(n) == r
    assert aproxLimSeno4(n) == r
    assert aproxLimSeno5(n) == r
    assert aproxLimSeno6(n) == r

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

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

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('aproxLimSeno1(3*10**5)')
#    0.03 segundos
#    >>> tiempo('aproxLimSeno2(3*10**5)')
#    Process Python violación de segmento (core dumped)
#    >>> tiempo('aproxLimSeno3(3*10**5)')
#    Process Python violación de segmento (core dumped)
#    >>> tiempo('aproxLimSeno4(3*10**5)')
#    Process Python violación de segmento (core dumped)
#    >>> tiempo('aproxLimSeno5(3*10**5)')
#    0.04 segundos
#    >>> tiempo('aproxLimSeno6(3*10**5)')
#    0.07 segundos
#
#    >>> tiempo('aproxLimSeno1(10**7)')
#    1.29 segundos
#    >>> tiempo('aproxLimSeno5(10**7)')
#    1.40 segundos
#    >>> tiempo('aproxLimSeno6(10**7)')
#    1.45 segundos

# 1ª definición de errorLimSeno
# ============================

# naturales es el generador de los números naturales positivos, Por
# ejemplo,
#    >>> list(islice(naturales(), 5))
#    [1, 2, 3, 4, 5]
def naturales() -> Iterator[int]:
    i = 1
    while True:
        yield i
        i += 1

def errorLimSeno1(x: float) -> int:
    return list(islice((n for n in naturales()
                        if abs(1 - sin(1/n)/(1/n)) < x), 1))[0]

# # 2ª definición de errorLimSeno
# # ============================

def errorLimSeno2(x: float) -> int:
    def aux(n: int) -> int:
        if abs(1 - sin(1/n)/(1/n)) < x:
            return n
        return aux(n + 1)

    return aux(1)

# 3ª definición de errorLimSeno
# ============================

def errorLimSeno3(x: float) -> int:
    return list(islice(dropwhile(lambda n: abs(1 - sin(1/n)/(1/n)) >= x,
                                 naturales()),
                       1))[0]


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

@given(st.integers(min_value=1, max_value=100))
def test_errorLimSeno(n: int) -> None:
    r = errorLimSeno1(n)
    assert errorLimSeno2(n) == r
    assert errorLimSeno3(n) == r

# La comprobación es
#    src> poetry run pytest -q limite_del_seno.py
#    2 passed in 0.60s

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

# La comparación es
#    >>> tiempo('errorLimSeno1(10**(-12))')
#    0.07 segundos
#    >>> tiempo('errorLimSeno2(10**(-12))')
#    Process Python violación de segmento (core dumped)
#    >>> tiempo('errorLimSeno3(10**(-12))')
#    0.10 segundos

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

from itertools import dropwhile, islice

from math import sin

from sys import setrecursionlimit

from timeit import Timer, default_timer

from typing import Iterator

from hypothesis import given

from hypothesis import strategies as st

setrecursionlimit(10**6)

# 1ª definición de aproxLimSeno

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

def aproxLimSeno1(k: int) -> list[float]:

return [sin(1/n)/(1/n) for n in range(1, k + 1)]

# 2ª definición de aproxLimSeno

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

def aproxLimSeno2(n: int) -> list[float]:

if n == 0:

return []

return aproxLimSeno2(n - 1) + [sin(1/n)/(1/n)]

# 3ª definición de aproxLimSeno

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

def aproxLimSeno3(n: int) -> list[float]:

def aux(n: int) -> list[float]:

if n == 0:

return []

return [sin(1/n)/(1/n)] + aux(n - 1)

return list(reversed(aux(n)))

# 4ª definición de aproxLimSeno

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

def aproxLimSeno4(n: int) -> list[float]:

def aux(xs: list[float], n: int) -> list[float]:

if n == 0:

return xs

return aux([sin(1/n)/(1/n)] + xs, n - 1)

return aux([], n)

# 5ª definición de aproxLimSeno

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

def aproxLimSeno5(n: int) -> list[float]:

return list(map((lambda k: sin(1/k)/(1/k)), range(1, n+1)))

# 6ª definición de aproxLimSeno

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

def aproxLimSeno6(n: int) -> list[float]:

r = []

for k in range(1, n+1):

r.append(sin(1/k)/(1/k))

return r

# Comprobación de equivalencia de aproxLimSeno

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

# La propiedad es

@given(st.integers(min_value=1, max_value=100))

def test_aproxLimSeno(n: int) -> None:

r = aproxLimSeno1(n)

assert aproxLimSeno2(n) == r

assert aproxLimSeno3(n) == r

assert aproxLimSeno4(n) == r

assert aproxLimSeno5(n) == r

assert aproxLimSeno6(n) == r

# La comprobación es

# src> poetry run pytest -q limite_del_seno.py

# 1 passed in 0.60s

# Comparación de eficiencia de aproxLimSeno

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

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('aproxLimSeno1(3*10**5)')

# 0.03 segundos

# >>> tiempo('aproxLimSeno2(3*10**5)')

# Process Python violación de segmento (core dumped)

# >>> tiempo('aproxLimSeno3(3*10**5)')

# Process Python violación de segmento (core dumped)

# >>> tiempo('aproxLimSeno4(3*10**5)')

# Process Python violación de segmento (core dumped)

# >>> tiempo('aproxLimSeno5(3*10**5)')

# 0.04 segundos

# >>> tiempo('aproxLimSeno6(3*10**5)')

# 0.07 segundos

# >>> tiempo('aproxLimSeno1(10**7)')

# 1.29 segundos

# >>> tiempo('aproxLimSeno5(10**7)')

# 1.40 segundos

# >>> tiempo('aproxLimSeno6(10**7)')

# 1.45 segundos

# 1ª definición de errorLimSeno

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

# naturales es el generador de los números naturales positivos, Por

# ejemplo,

# >>> list(islice(naturales(), 5))

# [1, 2, 3, 4, 5]

def naturales() -> Iterator[int]:

i = 1

while True:

yield i

i += 1

def errorLimSeno1(x: float) -> int:

return list(islice((n for n in naturales()

if abs(1 - sin(1/n)/(1/n)) < x), 1))[0]

# # 2ª definición de errorLimSeno

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

def errorLimSeno2(x: float) -> int:

def aux(n: int) -> int:

if abs(1 - sin(1/n)/(1/n)) < x:

return n

return aux(n + 1)

return aux(1)

# 3ª definición de errorLimSeno

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

def errorLimSeno3(x: float) -> int:

return list(islice(dropwhile(lambda n: abs(1 - sin(1/n)/(1/n)) >= x,

naturales()),

1))[0]

# Comprobación de equivalencia de errorLimSeno

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

@given(st.integers(min_value=1, max_value=100))

def test_errorLimSeno(n: int) -> None:

r = errorLimSeno1(n)

assert errorLimSeno2(n) == r

assert errorLimSeno3(n) == r

# La comprobación es

# src> poetry run pytest -q limite_del_seno.py

# 2 passed in 0.60s

# Comparación de eficiencia de errorLimSeno

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

# La comparación es

# >>> tiempo('errorLimSeno1(10**(-12))')

# 0.07 segundos

# >>> tiempo('errorLimSeno2(10**(-12))')

# Process Python violación de segmento (core dumped)

# >>> tiempo('errorLimSeno3(10**(-12))')

# 0.10 segundos

El código se encuentra en GitHub.

5. Cálculo del número π mediante la fórmula de Leibniz

El número π puede calcularse con la fórmula de Leibniz

   π/4 = 1 - 1/3 + 1/5 - 1/7 + ...+ (-1)**n/(2*n+1) + ...

1	π/4 = 1 - 1/3 + 1/5 - 1/7 + ...+ (-1)*n/(2n+1) + ...

Definir las funciones

   calculaPi :: Int -> Double
   errorPi   :: Double -> Int

1 2	calculaPi :: Int -> Double errorPi :: Double -> Int

tales que

calculaPi n es la aproximación del número π calculada mediante la expresión

     4*(1 - 1/3 + 1/5 - 1/7 + ...+ (-1)**n/(2*n+1))

1	4(1 - 1/3 + 1/5 - 1/7 + ...+ (-1)n/(2n+1))

Por ejemplo,

     calculaPi 3    ==  2.8952380952380956
     calculaPi 300  ==  3.1449149035588526

1 2	calculaPi 3 == 2.8952380952380956 calculaPi 300 == 3.1449149035588526

errorPi x es el menor número de términos de la serie necesarios para obtener pi con un error menor que x. Por ejemplo,

     errorPi 0.1    ==    9
     errorPi 0.01   ==   99
     errorPi 0.001  ==  999

errorPi 0.1 == 9

errorPi 0.01 == 99

errorPi 0.001 == 999

Soluciones en Haskell

import Test.QuickCheck

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

calculaPi1 :: Int -> Double
calculaPi1 k = 4 * sum [(-1)**n/(2*n+1) | n <- [0..k']]
  where k' = fromIntegral k

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

calculaPi2 :: Int -> Double
calculaPi2 0 = 4
calculaPi2 n = calculaPi2 (n-1) + 4*(-1)**n'/(2*n'+1)
  where n' = fromIntegral n

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

calculaPi3 :: Int -> Double
calculaPi3 = aux . fromIntegral
  where aux 0 = 4
        aux n = 4*(-1)**n/(2*n+1) + aux (n-1)

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

-- La propiedad es
prop_calculaPi :: Positive Int -> Bool
prop_calculaPi (Positive k) =
  all (== calculaPi1 k)
      [calculaPi2 k,
       calculaPi3 k]

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

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

-- La comparación es
--    λ> calculaPi1 (10^6)
--    3.1415936535887745
--    (1.31 secs, 609,797,408 bytes)
--    λ> calculaPi2 (10^6)
--    3.1415936535887745
--    (1.68 secs, 723,032,272 bytes)
--    λ> calculaPi3 (10^6)
--    3.1415936535887745
--    (2.22 secs, 1,099,032,608 bytes)


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

errorPi1 :: Double -> Int
errorPi1 x =
  head [n | n <- [1..]
          , abs (pi - calculaPi1 n) < x]

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

errorPi2 :: Double -> Int
errorPi2 x = aux 1
  where aux n | abs (pi - calculaPi1 n) < x = n
              | otherwise                   = aux (n+1)

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

-- La propiedad es
prop_errorPi :: Positive Double -> Bool
prop_errorPi (Positive x) =
  errorPi1 x == errorPi2 x

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

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

-- La comparación es
--    λ> errorPi1 0.0005
--    1999
--    (1.88 secs, 1,189,226,384 bytes)
--    λ> errorPi2 0.0005
--    1999
--    (1.87 secs, 1,213,341,096 bytes)

import Test.QuickCheck

-- 1ª definición de calculaPi

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

calculaPi1 :: Int -> Double

calculaPi1 k = 4 * sum [(-1)**n/(2*n+1) | n <- [0..k']]

where k' = fromIntegral k

-- 2ª definición de calculaPi

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

calculaPi2 :: Int -> Double

calculaPi2 0 = 4

calculaPi2 n = calculaPi2 (n-1) + 4*(-1)**n'/(2*n'+1)

where n' = fromIntegral n

-- 3ª definición de calculaPi

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

calculaPi3 :: Int -> Double

calculaPi3 = aux . fromIntegral

where aux 0 = 4

aux n = 4*(-1)**n/(2*n+1) + aux (n-1)

-- Comprobación de equivalencia de calculaPi

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

-- La propiedad es

prop_calculaPi :: Positive Int -> Bool

prop_calculaPi (Positive k) =

all (== calculaPi1 k)

[calculaPi2 k,

calculaPi3 k]

-- La comprobación es

-- λ> quickCheck prop_calculaPi

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de calculaPi

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

-- La comparación es

-- λ> calculaPi1 (10^6)

-- 3.1415936535887745

-- (1.31 secs, 609,797,408 bytes)

-- λ> calculaPi2 (10^6)

-- 3.1415936535887745

-- (1.68 secs, 723,032,272 bytes)

-- λ> calculaPi3 (10^6)

-- 3.1415936535887745

-- (2.22 secs, 1,099,032,608 bytes)

-- 1ª definición de errorPi

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

errorPi1 :: Double -> Int

errorPi1 x =

head [n | n <- [1..]

, abs (pi - calculaPi1 n) < x]

-- 2ª definición de errorPi

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

errorPi2 :: Double -> Int

errorPi2 x = aux 1

where aux n | abs (pi - calculaPi1 n) < x = n

| otherwise = aux (n+1)

-- Comprobación de equivalencia de errorPi

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

-- La propiedad es

prop_errorPi :: Positive Double -> Bool

prop_errorPi (Positive x) =

errorPi1 x == errorPi2 x

-- La comprobación es

-- λ> quickCheck prop_errorPi

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia de errorPi

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

-- La comparación es

-- λ> errorPi1 0.0005

-- 1999

-- (1.88 secs, 1,189,226,384 bytes)

-- λ> errorPi2 0.0005

-- 1999

-- (1.87 secs, 1,213,341,096 bytes)

El código se encuentra en GitHub.

Soluciones en Python

from itertools import islice
from math import pi
from sys import setrecursionlimit
from timeit import Timer, default_timer
from typing import Iterator

from hypothesis import given
from hypothesis import strategies as st

setrecursionlimit(10**6)

# 1ª definición de calculaPi
# ==========================

def calculaPi1(k: int) -> float:
    return 4 * sum(((-1)**n/(2*n+1) for n in range(0, k+1)))

# 2ª definición de calculaPi
# ==========================

def calculaPi2(n: int) -> float:
    if n == 0:
        return 4
    return calculaPi2(n-1) + 4*(-1)**n/(2*n+1)

# 3ª definición de calculaPi
# ==========================

def calculaPi3(n: int) -> float:
    r = 1
    for k in range(1, n+1):
        r = r + (-1)**k/(2*k+1)
    return 4 * r

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

# La propiedad es
@given(st.integers(min_value=1, max_value=100))
def test_calculaPi(n: int) -> None:
    r = calculaPi1(n)
    assert calculaPi2(n) == r
    assert calculaPi3(n) == r

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

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

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('calculaPi1(10**6)')
#    0.37 segundos
#    >>> tiempo('calculaPi2(10**6)')
#    Process Python violación de segmento (core dumped)
#    >>> tiempo('calculaPi3(10**6)')
#    0.39 segundos

# 1ª definición de errorPi
# ========================

# naturales es el generador de los números naturales positivos, Por
# ejemplo,
#    >>> list(islice(naturales(), 5))
#    [1, 2, 3, 4, 5]
def naturales() -> Iterator[int]:
    i = 1
    while True:
        yield i
        i += 1

def errorPi1(x: float) -> int:
    return list(islice((n for n in naturales()
                        if abs(pi - calculaPi1(n)) < x), 1))[0]

# 2ª definición de errorPi
# ========================

def errorPi2(x: float) -> int:
    def aux(n: int) -> int:
        if abs(pi - calculaPi1(n)) < x:
            return n
        return aux(n + 1)

    return aux(1)

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

@given(st.integers(min_value=1, max_value=100))
def test_errorPi(n: int) -> None:
    assert errorPi1(n) == errorPi2(n)

# La comprobación es
#    src> poetry run pytest -q calculo_de_pi_mediante_la_formula_de_Leibniz.py
#    2 passed in 0.60s

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

# La comparación es
#    >>> tiempo('errorPi1(0.0005)')
#    0.63 segundos
#    >>> tiempo('errorPi2(0.0005)')
#    0.58 segundos

100

101

102

103

104

105

106

107

108

109

110

111

from itertools import islice

from math import pi

from sys import setrecursionlimit

from timeit import Timer, default_timer

from typing import Iterator

from hypothesis import given

from hypothesis import strategies as st

setrecursionlimit(10**6)

# 1ª definición de calculaPi

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

def calculaPi1(k: int) -> float:

return 4 * sum(((-1)**n/(2*n+1) for n in range(0, k+1)))

# 2ª definición de calculaPi

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

def calculaPi2(n: int) -> float:

if n == 0:

return 4

return calculaPi2(n-1) + 4*(-1)**n/(2*n+1)

# 3ª definición de calculaPi

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

def calculaPi3(n: int) -> float:

r = 1

for k in range(1, n+1):

r = r + (-1)**k/(2*k+1)

return 4 * r

# Comprobación de equivalencia de calculaPi

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

# La propiedad es

@given(st.integers(min_value=1, max_value=100))

def test_calculaPi(n: int) -> None:

r = calculaPi1(n)

assert calculaPi2(n) == r

assert calculaPi3(n) == r

# La comprobación es

# src> poetry run pytest -q calculo_de_pi_mediante_la_formula_de_Leibniz.py

# 1 passed in 0.14s

# Comparación de eficiencia de calculaPi

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

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('calculaPi1(10**6)')

# 0.37 segundos

# >>> tiempo('calculaPi2(10**6)')

# Process Python violación de segmento (core dumped)

# >>> tiempo('calculaPi3(10**6)')

# 0.39 segundos

# 1ª definición de errorPi

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

# naturales es el generador de los números naturales positivos, Por

# ejemplo,

# >>> list(islice(naturales(), 5))

# [1, 2, 3, 4, 5]

def naturales() -> Iterator[int]:

i = 1

while True:

yield i

i += 1

def errorPi1(x: float) -> int:

return list(islice((n for n in naturales()

if abs(pi - calculaPi1(n)) < x), 1))[0]

# 2ª definición de errorPi

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

def errorPi2(x: float) -> int:

def aux(n: int) -> int:

if abs(pi - calculaPi1(n)) < x:

return n

return aux(n + 1)

return aux(1)

# Comprobación de equivalencia de errorPi

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

@given(st.integers(min_value=1, max_value=100))

def test_errorPi(n: int) -> None:

assert errorPi1(n) == errorPi2(n)

# La comprobación es

# src> poetry run pytest -q calculo_de_pi_mediante_la_formula_de_Leibniz.py

# 2 passed in 0.60s

# Comparación de eficiencia de errorPi

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

# La comparación es

# >>> tiempo('errorPi1(0.0005)')

# 0.63 segundos

# >>> tiempo('errorPi2(0.0005)')

# 0.58 segundos

El código se encuentra en GitHub.

1. Suma de múltiplos de 3 ó 5

2. Puntos en el círculo

3. Aproximación del número e

4. Aproximación al límite de sen(x)/x cuando x tiende a cero

5. Cálculo del número π mediante la fórmula de Leibniz

Se puede imprimir o compartir con