#!/bin/env python
# -*- coding: utf-8 -*-
# cos_quadratic_interpolation.py
# 
# 2014.10.27

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),
    return

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

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

x0 = - math.pi / 8
x1 = 0
x2 = math.pi / 8
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 "quadratic Lagrange interpolation of cos(x) at -pi/8, 0, pi/8:"
print "p(x) =",
realpolwriteln(p)
print

f = [1, 0, - 0.5]
print "Taylor expansion at x = 0 up to the 2nd term:"
print "f(x) =",
realpolwriteln(f)
print

x = 0
while x <= x2:
    y = math.cos(x)
    px = poleval(p, x)
    fx = poleval(f, x)
    print "x =", x
    print "cos(x) =", y
    print "p(x) =", px,
    print ",  error =", abs(y - px)
    print "f(x) =", fx,
    print ",  error =", abs(y - fx)
    print
    x += 0.05