Rump's example

rump1988.py

from decimal import *

# Siegfried M. Rump (1988), 
# “Algorithms for Verified Inclusions: Theory and Practice,”,
# in R. E. Moore, Reliability in Computing: The Role of Interval Methods in Scientific Computing, 
# Boston: Academic Press, pp. 109-126., doi:10.15480/882.316

def f(a, b):
    return (333.75*(b**6)) + (a**2)*(11*(a**2)*(b**2) - (b**6) - (121*(b**4)) - 2) + (5.5*(b**8)) + (a / (2*b))

def f2(x, y):
    return (333.75 - x**2) * (y**6) + (x**2) * (11*(x**2)*(y**2) - 121*(y**4) - 2) + 5.5*(y**8) + (x/(2*y))

def f_dec(a_r, b_r):
    a = Decimal(a_r)
    b = Decimal(b_r)
    return (Decimal('333.75')*(b**Decimal('6'))) \
        + (a**Decimal('2'))*(Decimal('11')*(a**Decimal('2'))*(b**Decimal('2')) \
            - (b**Decimal('6')) - (Decimal('121')*(b**Decimal('4'))) - Decimal('2')) \
        + (Decimal('5.5')*(b**Decimal('8'))) \
        + (a / (Decimal('2')*b))

def f_dec2(x_r, y_r):
    x = Decimal(x_r)
    y = Decimal(y_r)
    return (Decimal('333.75') - x**Decimal('2')) * (y**Decimal('6')) + (x**Decimal('2')) * (Decimal('11')*(x**Decimal('2'))*(y**Decimal('2')) - Decimal('121')*(y**Decimal('4')) - Decimal('2')) + Decimal('5.5')*(y**Decimal('8')) + (x/(Decimal('2')*y))


if __name__ == '__main__':
    a = 77617.0
    b = 33096.0
    print(f(a, b))
    print(f2(a, b))

    for p in range(8, 40):
        getcontext().prec = p
        print(f"f{p}=", f_dec('77617.0', '33096.0'))

    for p in range(8, 40):
        getcontext().prec = p
        print(f"f{p}=", f_dec2('77617.0', '33096.0'))

rump1988.rb

require 'bigdecimal'

# Siegfried M. Rump (1988), 
# “Algorithms for Verified Inclusions: Theory and Practice,”,
# in R. E. Moore, Reliability in Computing: The Role of Interval Methods in Scientific Computing, 
# Boston: Academic Press, pp. 109-126., doi:10.15480/882.316
def f(a, b)
    (333.75*(b**6)) + (a**2)*(11*(a**2)*(b**2) - (b**6) - (121*(b**4)) - 2) + (5.5*(b**8)) + (a / (2*b))
end

# Loh, E., Walster, G.W. 
# Rump's Example Revisited.
# Reliable Computing 8, 245–248 (2002).
# https://doi.org/10.1023/A:1015569431383
def f2(x, y)
    (333.75 - x**2) * (y**6) + (x**2) * (11*(x**2)*(y**2) - 121*(y**4) - 2) + 5.5*(y**8) + (x/(2*y))
end

if $0 == __FILE__
    p f(77617.0, 33096.0)       # => -1.1805916207174113e+21
    p f2(77617.0, 33096.0)      # => 1.1726039400531787
    p f(BigDecimal('77617.0'), BigDecimal('33096.0'))       # => -0.827396059946821368e0
    p f2(BigDecimal('77617.0'), BigDecimal('33096.0'))      # => -0.827396059946821368e0
end