#----------------------------------------------------------------------
# Python Polynomials
#
# Copyright (C) December 2018 -- Dr. William T. Verts
#
# Polynomials are represented by lists of numbers, where the index
# into the list is the exponent of X.  For example, [1,2,3] means
# 3X^2 + 2X^1 + 1X^0 = 3X^2 + 2X + 1
#
# For more information and more functions, see pages 272-275 in the
# Companion.
#----------------------------------------------------------------------

#----------------------------------------------------------------------
# polyNormalize changes floats very close to an integer value into ints
# and also eliminates high-order terms with zero coefficients.
#----------------------------------------------------------------------

def polyNormalize (P):
    Q = [0] * max(len(P),1)
    
    for I in range(len(P)):
        if abs(P[I] - float(int(P[I]))) < 1.0E-14:
            Q[I] = int(P[I])
        else:
            Q[I] = P[I]
            
    while (len(Q) > 1) and (Q[-1] == 0): del Q[-1]
    return Q

#----------------------------------------------------------------------
# polyAdd and polySubtract add corresponding terms from two or three
# polynomials (either P0+P1+P2 or P0-P1-P2).  P1 or P2 default to empty
# lists to indicate no polynomial in those positions unless overridden.
#----------------------------------------------------------------------

def polyAdd (P0,P1=[],P2=[]):
    Q = [0] * max(len(P0),len(P1),len(P2))
    for I in range(len(P0)): Q[I] = Q[I] + P0[I]
    for I in range(len(P1)): Q[I] = Q[I] + P1[I]
    for I in range(len(P2)): Q[I] = Q[I] + P2[I]
    return polyNormalize(Q)
    
def polySubtract (P0,P1=[],P2=[]):
    Q = [0] * max(len(P0),len(P1),len(P2))
    for I in range(len(P0)): Q[I] = Q[I] + P0[I]
    for I in range(len(P1)): Q[I] = Q[I] - P1[I]
    for I in range(len(P2)): Q[I] = Q[I] - P2[I]
    return polyNormalize(Q)
    
#----------------------------------------------------------------------
# polyMultiply has to generate every inner product and add it to the
# correct slot.  For example, 2X^3 * 4X^2 = (2*4)X^(3+2) = 8X^5 
#----------------------------------------------------------------------

def polyMultiply (P0,P1):
    Q = [0] * (len(P0) + len(P1))
    for I in range(len(P0)):
        for J in range(len(P1)):
            Q[I+J] = Q[I+J] + P0[I] * P1[J]
    return polyNormalize(Q)
    
#----------------------------------------------------------------------
# polyDifferentiate computes the calculus derivative of a polynomial
# with respect to X.
#----------------------------------------------------------------------

def polyDifferentiate (P0):
    Q = [0] * (len(P0) - 1)
    for I in range(1,len(P0)): Q[I-1] = I * P0[I]
    return polyNormalize(Q)
    
#----------------------------------------------------------------------
# Take a list indicating a polynomial and turn it into a printable
# string that looks close to how we would write it in math notation.
#----------------------------------------------------------------------

def polyStr(P0):
    S = ""
    for Exponent in range(len(P0)-1,-1,-1):
        Coefficient = P0[Exponent]
        if (Coefficient != 0):                  # This was the line I forgot in class
            if (Coefficient < 0):
                if (S != ""):
                    S = S + " - "
                else:
                    S = S + "-"
            else:
                if (S != ""):
                    S = S + " + "
                 
            if (Exponent == 0) or (abs(Coefficient) != 1):
                S = S + str(abs(Coefficient))
            
            if (Exponent >= 1):
                S = S + "X"
            
            if (Exponent > 1):
                S = S + "^" + str(Exponent)
          
    if (S == ""): S = "0"  
    return S    
    
#----------------------------------------------------------------------
# Test all the functions.
#----------------------------------------------------------------------
        
P0 = [-1.5, 2, 0, 4]
P1 = [3.5, -2, 3]
print "P0 = ", P0, " = ", polyStr(P0)
print "P1 = ", P1, " = ", polyStr(P1)
P2 = polyAdd(P0,P1)
print "P0 + P1 = ", P2, " = ", polyStr(P2)
P3 = polySubtract(P0,P1)
print "P0 - P1 = ", P3, " = ", polyStr(P3)
P4 = polyMultiply(P0,P1)
print "P0 * P1 = ", P4, " = ", polyStr(P4)
P5 = polyDifferentiate(P0)
print "d/dx(P0) = ", P5, " = ", polyStr(P5)
P6 = polySubtract(P0,P0)
print "P0 - P0 = ", P6, " = ", polyStr(P6)
P7 = polyAdd(P0,P0,P1)
print "P0 + P0 + P1 = ", P7, " = ", polyStr(P7)
P8 = polyDifferentiate(polyMultiply(P0,polyAdd(P0,P1)))
print "d/dx(P0 * (P0 + P1)) = ", P8, " = ", polyStr(P8)