# -*- coding: utf-8 -*-
# cos_quadratic_interpolation.py
# 
# 2014.10.27
#
# 2次ラグランジュ補間とテイラー展開の比較
#
# ラグランジュ補間：区間内で満遍なくそこそこ近似
# テイラー展開：　　展開点の近くで精度よく近似

import math

#------------------------------#
#  polynomial operations       #
#------------------------------#

#   f = [f[0], f[1], ..., f[d]]
# represents the polynomial
#   f(x) = c[d] x^d + ... + c[1] * x + c[0]
# zero polynomial = [0]
# degree of f = len(f) - 1

# weather f is zero polynomial or not
def polzeroq(f):
    return (f == [0])

# degree of f
def deg(f):
    return len(f) - 1

# set the degree
# ( necessary after addition, subtraction and so on )
def cutzero(f):
    while f[-1] == 0:
        if f == [0]:
            break
        del f[-1]
    return

# addition
def poladd(f, g):
    if polzeroq(f):
        return g
    if polzeroq(g):
        return f
    m = deg(f)
    n = deg(g)
    c = []
    if m >= n:
        for i in range(n + 1):
            c.append(f[i] + g[i])
        for i in range(n + 1, m + 1):
            c.append(f[i])
        cutzero(c)
    else:
        for i in range(m + 1):
            c.append(f[i] + g[i])
        for i in range(m + 1, n + 1):
            c.append(g[i])
    return c

# scalar multiple
def polscmul(a, f):
    if polzeroq(f):
        return [0]
    c = []
    for i in range(len(f)):
        c.append(a * f[i])
    return c

# multiplication
def polmul(f, g):
    if polzeroq(f) or polzeroq(g):
        return [0]
    m = deg(f)
    n = deg(g)
    c = [0] * (m + n + 1)
    for i in range(m + 1):
        a = f[i]
        for j in range(n + 1):
            c[i + j] = c[i + j] + (a * g[j])
    return c

# evaluate f(x) at x = a
def poleval(f, a):
    w = 0
    d = deg(f)
    for i in range(d, -1, -1):
        w = a * w + f[i]
    return w

# string expressing polynomial f with real coefficients
def realpolstr(f):
    if f == [0]:
        return '0'
    s = ''
    d = deg(f)
    for i in range(d, -1, -1):
        c = f[i]
        if c != 0:
            if c > 0:
                if i < d:
                    s = s + ' + '
            else:
                s = s + ' - '
            c = abs(c)
            if c != 1 or i == 0:
                s = s + ('%f' % c)
                if i > 0:
                    s = s + ' * '
            if i > 1:
                s = s + ('x^%d' % i)
            if i == 1:
                s = s + 'x'
    return s

# output a polynomial with real coefficients
def realpolwrite(f):
    print(realpolstr(f), end = '')
    return

# + line feed
def realpolwriteln(f):
    print(realpolstr(f))
    return

#------------------------------#
#  main routine                #
#------------------------------#

#x0, x1, x2 = - math.pi / 8, 0.0, math.pi / 8
x0, x1, x2 = - 0.5, 0.0, 0.5
y0 = math.cos(x0)
y1 = math.cos(x1)
y2 = math.cos(x2)
n0 = polscmul(1.0 / ((x0 - x1) * (x0 - x2)), polmul([-x1, 1], [-x2, 1]))
n1 = polscmul(1.0 / ((x1 - x0) * (x1 - x2)), polmul([-x0, 1], [-x2, 1]))
n2 = polscmul(1.0 / ((x2 - x0) * (x2 - x1)), polmul([-x0, 1], [-x1, 1]))
p = poladd(poladd(polscmul(y0, n0), polscmul(y1, n1)), polscmul(y2, n2))

print('ラグランジュ補間とテイラー展開の比較')
print()
print('cos(x) の２次のラグランジュ補間：')
print('p(x) =', end = '')
realpolwriteln(p)
print()

f = [1, 0, - 0.5]
print('cos(x) の２次までのテイラー展開：')
print('f(x) =', end = '')
realpolwriteln(f)
print()

x = 0
while x <= 1.01:
    y = math.cos(x)
    px = poleval(p, x)
    fx = poleval(f, x)
    print('x = %1.2lf' % x)
    print('cos(x) = %12.10lf' % y)
    print('  p(x) = %12.10lf,' % px, end = '')
    print(' ラグランジュ補間の誤差 = %12.10lf' % abs(y - px))
    print('  f(x) = %12.10lf,' % fx, end = '')
    print(' テイラー展開の誤差     = %12.10lf' % abs(y - fx))
    print()
    x += 0.05