Método de Herón para calcular la raíz cuadrada

José A. Alonso

29-10-2023

El método de Herón para calcular la raíz cuadrada de un número se basa en las siguientes propiedades:

Si \(y\) es una aproximación de la raíz cuadrada de \(x\), entonces \[\frac{y+\frac{x}{y}}{2}\] es una aproximación mejor.
El límite de la sucesión definida por \begin{align} x_{0} &= 1 \newline x_{n+1} &= \frac{x_n+\frac{x}{x_n}}{2} \end{align} es la raíz cuadrada de x.

Definir la función

   raiz :: Double -> Double

tal que raiz x es la raíz cuadrada de x calculada usando la propiedad anterior con una aproximación de 0.00001 y tomando como valor inicial 1. Por ejemplo,

   raiz 9  ==  3.000000001396984

Soluciones

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

Soluciones en Haskell

module Metodo_de_Heron_para_calcular_la_raiz_cuadrada where

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

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

raiz :: Double -> Double
raiz x = raizAux 1
  where raizAux y | aceptable y = y
                  | otherwise   = raizAux (mejora y)
        aceptable y = abs(y*y-x) < 0.00001
        mejora y    = 0.5*(y+x/y)

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

raiz2 :: Double -> Double
raiz2 x = until aceptable mejora 1
  where aceptable y = abs(y*y-x) < 0.00001
        mejora y    = 0.5*(y+x/y)

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

-- La propiedad es
prop_raiz :: Positive Double -> Bool
prop_raiz (Positive x) =
  raiz x ~= sqrt x &&
  raiz2 x ~= sqrt x
  where
    a ~= b = abs (a-b) < 0.001

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

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

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "e1" $
    raiz 9 `shouldBe`  3.000000001396984
  it "e2" $
    raiz2 9 `shouldBe`  3.000000001396984

-- La verificación es
--    λ> verifica
--
--    e1
--    e2
--
--    Finished in 0.0008 seconds
--    2 examples, 0 failures

Soluciones en Python

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

def raiz(x : float) -> float:
    def aceptable(y: float) -> bool:
        return abs(y*y-x) < 0.00001
    def mejora(y: float) -> float:
        return 0.5*(y+x/y)
    def raizAux(y: float) -> float:
        if aceptable(y):
            return y
        return raizAux(mejora(y))
    return raizAux(1)

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

def raiz2(x: float) -> float:
    def aceptable(y: float) -> bool:
        return abs(y*y-x) < 0.00001
    def mejora(y: float) -> float:
        return 0.5*(y+x/y)
    y = 1.0
    while not aceptable(y):
        y = mejora(y)
    return y

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

def test_raiz() -> None:
    assert raiz(9)  ==  3.000000001396984
    assert raiz2(9)  ==  3.000000001396984
    print("Verificado")

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