# L-7 MCS 507 Wed 6 Sep 2023 : intvalmotiv.py

"""
The example is taken from the paper of Stefano Taschini:
"Interval Arithmetic: Python Implementation and Applications"
Proceedings of the 7th Python in Science Conference (SciPy 2008).

While we get the final correct approximation when we use 36
decimal places using mpmath of SymPy, the choice of 36 is not obvious.
"""

import sympy as sp
import mpmath

def exact_versus_numeric(A, B):
    """
    Compares the exact to the numeric evaluation
    of a expression at (A, B).
    """
    F = lambda 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)
    Z = F(float(A), float(B))
    # to check the answer with SymPy :
    X, Y = sp.var('x, y')
    G = (sp.Rational(33375)/100 - X**2)*Y**6 \
        + X**2*(11*X**2*Y**2 - 121*Y**4 - 2) \
        + sp.Rational(55)/10*Y**8 \
        + sp.Rational(1)*X/(2*Y)
    print('evaluating', G, 'at', (A, B))
    E = sp.Subs(G, (X, Y), (A, B)).doit()
    E15 = E.evalf(15)
    print('numerical value :', Z)
    print('exact value :', E, '~', E15)
    print('error :', abs(E15 - Z))

def variable_precision(A, B):
    """
    The expression is evaluated with variable precision
    at (A, B).
    """
    # let us double the precision to 30
    mpmath.mp.dps = 30
    MP_F = lambda x, y: (mpmath.mpf('333.75') \
         - x**2)*y**6 \
         + x**2*(11*x**2*y**2 - 121*y**4 - 2) \
         + mpmath.mpf('5.5')*y**8 + x/(2*y)
    MP_A30 = mpmath.mpf(str(A))
    MP_B30 = mpmath.mpf(str(B))
    Z30 = MP_F(MP_A30, MP_B30)
    print('using 30 digits :', Z30)
    # adjusting the precision to 35
    mpmath.mp.dps = 35
    MP_A35 = mpmath.mpf(str(A))
    MP_B35 = mpmath.mpf(str(B))
    Z35 = MP_F(MP_A35, MP_B35)
    print('using 35 digits :', Z35)
    # adjusting the precision to 36
    mpmath.mp.dps = 36
    MP_A36 = mpmath.mpf(str(A))
    MP_B36 = mpmath.mpf(str(B))
    Z36 = MP_F(MP_A36, MP_B36)
    print('using 36 digits :', Z36)

def main():
    """
    Experimentation on a motivating example
    for the use of interval arithmetic.
    """
    (A, B) = (77617, 33096)
    exact_versus_numeric(A, B)
    variable_precision(A, B)

main()