# -*- coding: utf-8 -*-

#Copyright (c) 2013 Juntao ZH<zjuntor@gmail.com>

#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is
#furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in
#all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
#THE SOFTWARE.


#This program, as <numerical analysis> homework, is used to interpolate by being given some sample points, using Newtown and lagrange method. 






#@parm : x,y are the lists of sample point series
#@return : return all the first to (len(x)-1)th order of differential quotient of given sample points, or None for the length of x not #equals y. 

def diffQuotient(x,y):
    if len(x) != len(y):
        return None
    size = len(x)
    x = [xi*1.0 for xi in x]
    
    res = []

    for i in range(0,size-1):
        tmp = []
        for j in range(0,size-i-1):
            if i == 0:
                tmp.append([(y[j+1]-y[j])/(x[j+1]-x[j]),[j,j+1]])
            else:
                index = [subIndex for subIndex in res[i-1][j][1]]
                index.append(res[i-1][j+1][1][-1])
                tmp.append([(res[i-1][j+1][0]-res[i-1][j][0])/(x[res[i-1][j+1][1][-1]]-x[res[i-1][j][1][0]]),index])
        res.append(tmp)
                  
    return res


#@parm : x,y are the lists of sample point series
#@return : function of Newtown polynomial generated by given sample points. 

def newtownInterpolation(x,y):
    dq = diffQuotient(x,y)
    dq = [e[0] for e in dq]
    def func(v):
        res = y[0]
        for i in range(0,len(dq)):
            tmp = 1
            for j in dq[i][1][:-1]:
                tmp = tmp * (v-x[j])
            
            res = res + tmp * dq[i][0]
        return res
    return func



#@parm : x,y are the lists of sample point series
#@return : function of Lagrange polynomial generated by given sample points. 

def lagrangeInterpolation(x,y):
    size = len(x)
    
    def func(v):
        s = set(range(0,size))
        res = 0
        for i in range(0,size):
            tmp = 1
            rm = s - set([i])
            for j in rm:
                tmp = tmp * (v - x[j]) / (x[i] - x[j])
            
            res = res + tmp * y[i]      
        return res
    
    return func


