jmc_polynomial.py

# Symbolic integration of [(1 - x^a)^b * e^iwx] dx from -1 to 1
# see discussion on music-dsp mailing list
# "Re: [music-dsp] a family of simple polynomial windows and waveforms"

from math import factorial
from copy import deepcopy
import string

# binomial coefficient nCk
def binomial(n, k):
    try:
        result = factorial(n) // (factorial(n - k) * factorial(k))
    except ValueError:
        result = 0
    return result

#n = 4
#for k in range(0, n+1):
#    print n, "choose", k, "=", binomial(n, k)

class Polynomial: # polynomial in one variable
    def __init__(self, varName="x"):
        self.varName = varName
        self.coeffs = []
        
    def derivative(self):
        result = Polynomial(self.varName)
        result.coeffs = self.coeffs[1:]
        for k in range(0, len(result.coeffs)):
            result.coeffs[k] *= k+1
        return result

    def multiplied(self, m):
        result = deepcopy(self)
        for k in range(0, len(result.coeffs)):
            result.coeffs[k] *= m
        return result

    def evaluate(self, x):
        result = 0
        for k in range(0, len(self)):
            result += self[k] * (x**k)
        return result

    # sparse list implementation. returns 0 for empty cells
    def __setitem__(self, index, value):
        missing = index - len(self.coeffs) + 1
        if missing > 0:
            self.coeffs.extend([0] * missing)
        self.coeffs[index] = value
        
    def __getitem__(self, index):
        try: return self.coeffs[index]
        except IndexError: return 0

    def __len__(self):
        return len(self.coeffs)

    def __nonzero__(self):
        for c in self.coeffs:
            if c:
                return True
        return False

    def __str__(self):
        result = ""
        for k in range(len(self.coeffs)-1, -1, -1):
            if self.coeffs[k]:
                coeff = self.coeffs[k]

                if coeff > 0:
                    result += " + "
                else:
                    result += " - "

                if k == 0:
                    result += str(abs(coeff))
                else:
                    if abs(coeff) != 1:
                        result += str(abs(coeff))
                    if k == 1:
                        result += "x"
                    else:
                        result += "x^" + str(k)
        return result
    
    def __repr__(self):
        return self.__str__()

    

#p = Polynomial()
#p[0] = 0
#p[2] = 5
#print p


# expand a polynomial of the form (1 - x^a)^b using the binomial theorem
def expandJMC(a, b):
    p = Polynomial()
    for k in range(0, b+1): #[0,b]
        bCk = binomial(b, k)
        p[k*a] = bCk * (-1)**k
    return p

#print expandJMC(1, 3)
#print "~"
#p = expandJMC(2, 3)
#print p
#print p.derivative()


# special representation for a complex exponential of the form: (iw)^a e^(iw)
class ComplexExpIW:
    def __init__(self, varName="x", a=0):
        self.varName = varName
        self.a = a

    def antiderivative(self):
        return ComplexExpIW(self.varName, self.a-1)

    def __str__(self):
        if self.a == 0:
            result = "e^(iw"+self.varName+")"
        else:
            result = "(iw)^" + str(self.a) + "*e^(iw"+self.varName+")"
        return result

    def __repr__(self):
        return self.__str__()

eiw = ComplexExpIW()
#print eiw
#print eiw.antiderivative()

###########################################################################

# Simplify sum and difference of ComplexExpIW over [-1 1]
# Because our integrand is real, we expect the Fourier integral to be real
# hence the result of evaluating the antiderivative F as [F(1) - F(-1)]
# should result in ComplexExpIW collapsing into sin or cos.

# i^a
def i_a(a):
    powersOf_i = [ (1,0), (0,1), (-1,0), (0, -1) ] 
    return powersOf_i[a%4]

assert( i_a(-4) == (1, 0) )
assert( i_a(-3) == (0, 1) )
assert( i_a(-2) == (-1, 0) )
assert( i_a(-1) == (0, -1) )
assert( i_a(0) == (1, 0) )
assert( i_a(1) == (0, 1) )
assert( i_a(2) == (-1, 0) )
assert( i_a(3) == (0, -1) )
assert( i_a(4) == (1, 0) )
assert( i_a(5) == (0, 1) )     

class EiwSinusoid:
    def __init__(self, coeff, a, f):
        self.coeff = coeff
        self.a = a
        self.f = f

    def __str__(self):
        return str(self.coeff) + " w^" + str(self.a) + " " + f + "(w)"

# (iw)^a (e^iw + e^-iw)
def simplifyEiwSum(eiw):
    a = eiw.a
    # (e^iw + e^-iw) = 2 cos(w)
    # i^a
    ia = i_a(a)
    assert( ia[1] == 0 ) # imaginary part should be zero
    coeff = ia[0] * 2
    return EiwSinusoid(coeff, a, "cos")
    
# (iw)^a (e^iw - e^-iw)
def simplifyEiwDifference(eiw):
    a = eiw.a
    # (e^iw - e^-iw) = 2i cos(w)
    # i^a
    ia = i_a(a)
    assert( ia[0] == 0 ) # real part should be zero
    if ia[1] == 1:
        # i * i = -1
        coeff = -2
    else:
        assert(ia[1] == -1)
        # -i * i = 1
        coeff = 2    
    return EiwSinusoid(coeff, a, "sin")

###########################################################################


def integrateJMC(a, b, verbose=False):
    print "-----------"

    # (1 - x^a)^b

    print "a = " + str(a)
    print "b = " + str(b)
    print "integrating: (1 - x^" + str(a) + ")^" + str(b) + " * e^iwx dx from -1 to 1"

    if verbose:
        print "-----------"

    # generate the anti-derivative terms using integration by parts
    # this follows the tabular method described on the last page here:
    # http://www.willamette.edu/~mjaneba/courses/ma256/integrationbyparts.pdf
    p = expandJMC(a, b)
    eiw = ComplexExpIW().antiderivative()
    sign = +1
    antiderivativeTerms = []
    while p:
        t = (p.multiplied(sign), eiw)
        antiderivativeTerms.append(t)
        if verbose:
            print str(t)
        p = p.derivative()
        eiw = eiw.antiderivative()
        sign = -sign

    # evaluate antiderivative over [-1, 1]

    if verbose:
        print "-----------"

    simplifiedTerms = []

    for term in antiderivativeTerms:
        # polynomial evaluated at upper and lower bounds
        p_top = term[0].evaluate(1)
        p_bottom = term[0].evaluate(-1)

        eiw = term[1]

        #   (p_top * k e^iw - p_bottom * e^-iw)
        # p_top == p_bottom:
        #   p_top k (e^iw - e^-iw)  // the (e-e) bit is 2i sin(w). k had better have an imaginary part to cancel it
        # p_top == -p_bottom:
        #   p_top k (e^iw + e^-iw)  // the (e+e) bit is 2 cos(w). k had better not have an imaginary part

        if p_top == 0:
            assert(p_bottom == 0)
            # polynomials cancel
        else:
            #print p_top, p_bottom

            if p_top == p_bottom:
                #print str(p_top) + " (" + str(eiw) + " - " + str(eiw) + ")"
                sinusoid = simplifyEiwDifference(eiw)

                # result is of the form: k * w^a *{sin,cos}(w)
                k = p_top*sinusoid.coeff
                       
            else:
                assert(p_bottom == -p_top)
                #print str(p_top) + " (" + str(eiw) + " + " + str(eiw) + ")"
                sinusoid = simplifyEiwSum(eiw)

                # result is of the form: k * w^a *{sin,cos}(w)
                k = p_top*sinusoid.coeff

            expr = string.join([str(k), "(w^" + str(sinusoid.a)+")", sinusoid.f + "(w)"], "*")
            if verbose:
                print expr
            simplifiedTerms.append(expr)

    if verbose:
        print "-----------"

    print "result:"
    print string.join(simplifiedTerms, " + ")


#integrateJMC(2,1)
#integrateJMC(2,3)

for a in range(2, 8, 2):
    for b in range(1, 5, 1):
        integrateJMC(a,b, True)

"""
Sample output:

-----------
a = 2
b = 1
integrating: (1 - x^2)^1 * e^iwx dx from -1 to 1
-----------
( - x^2 + 1, (iw)^-1*e^(iwx))
( + 2x, (iw)^-2*e^(iwx))
( - 2, (iw)^-3*e^(iwx))
-----------
-4*(w^-2)*cos(w)
4*(w^-3)*sin(w)
-----------
result:
-4*(w^-2)*cos(w) + 4*(w^-3)*sin(w)
-----------
a = 2
b = 2
integrating: (1 - x^2)^2 * e^iwx dx from -1 to 1
-----------
( + x^4 - 2x^2 + 1, (iw)^-1*e^(iwx))
( - 4x^3 + 4x, (iw)^-2*e^(iwx))
( + 12x^2 - 4, (iw)^-3*e^(iwx))
( - 24x, (iw)^-4*e^(iwx))
( + 24, (iw)^-5*e^(iwx))
-----------
-16*(w^-3)*sin(w)
-48*(w^-4)*cos(w)
48*(w^-5)*sin(w)
-----------
result:
-16*(w^-3)*sin(w) + -48*(w^-4)*cos(w) + 48*(w^-5)*sin(w)
-----------
a = 2
b = 3
integrating: (1 - x^2)^3 * e^iwx dx from -1 to 1
-----------
( - x^6 + 3x^4 - 3x^2 + 1, (iw)^-1*e^(iwx))
( + 6x^5 - 12x^3 + 6x, (iw)^-2*e^(iwx))
( - 30x^4 + 36x^2 - 6, (iw)^-3*e^(iwx))
( + 120x^3 - 72x, (iw)^-4*e^(iwx))
( - 360x^2 + 72, (iw)^-5*e^(iwx))
( + 720x, (iw)^-6*e^(iwx))
( - 720, (iw)^-7*e^(iwx))
-----------
96*(w^-4)*cos(w)
-576*(w^-5)*sin(w)
-1440*(w^-6)*cos(w)
1440*(w^-7)*sin(w)
-----------
result:
96*(w^-4)*cos(w) + -576*(w^-5)*sin(w) + -1440*(w^-6)*cos(w) + 1440*(w^-7)*sin(w)
-----------
a = 2
b = 4
integrating: (1 - x^2)^4 * e^iwx dx from -1 to 1
-----------
( + x^8 - 4x^6 + 6x^4 - 4x^2 + 1, (iw)^-1*e^(iwx))
( - 8x^7 + 24x^5 - 24x^3 + 8x, (iw)^-2*e^(iwx))
( + 56x^6 - 120x^4 + 72x^2 - 8, (iw)^-3*e^(iwx))
( - 336x^5 + 480x^3 - 144x, (iw)^-4*e^(iwx))
( + 1680x^4 - 1440x^2 + 144, (iw)^-5*e^(iwx))
( - 6720x^3 + 2880x, (iw)^-6*e^(iwx))
( + 20160x^2 - 2880, (iw)^-7*e^(iwx))
( - 40320x, (iw)^-8*e^(iwx))
( + 40320, (iw)^-9*e^(iwx))
-----------
768*(w^-5)*sin(w)
7680*(w^-6)*cos(w)
-34560*(w^-7)*sin(w)
-80640*(w^-8)*cos(w)
80640*(w^-9)*sin(w)
-----------
result:
768*(w^-5)*sin(w) + 7680*(w^-6)*cos(w) + -34560*(w^-7)*sin(w) + -80640*(w^-8)*cos(w) + 80640*(w^-9)*sin(w)
-----------
a = 4
b = 1
integrating: (1 - x^4)^1 * e^iwx dx from -1 to 1
-----------
( - x^4 + 1, (iw)^-1*e^(iwx))
( + 4x^3, (iw)^-2*e^(iwx))
( - 12x^2, (iw)^-3*e^(iwx))
( + 24x, (iw)^-4*e^(iwx))
( - 24, (iw)^-5*e^(iwx))
-----------
-8*(w^-2)*cos(w)
24*(w^-3)*sin(w)
48*(w^-4)*cos(w)
-48*(w^-5)*sin(w)
-----------
result:
-8*(w^-2)*cos(w) + 24*(w^-3)*sin(w) + 48*(w^-4)*cos(w) + -48*(w^-5)*sin(w)
-----------
a = 4
b = 2
integrating: (1 - x^4)^2 * e^iwx dx from -1 to 1
-----------
( + x^8 - 2x^4 + 1, (iw)^-1*e^(iwx))
( - 8x^7 + 8x^3, (iw)^-2*e^(iwx))
( + 56x^6 - 24x^2, (iw)^-3*e^(iwx))
( - 336x^5 + 48x, (iw)^-4*e^(iwx))
( + 1680x^4 - 48, (iw)^-5*e^(iwx))
( - 6720x^3, (iw)^-6*e^(iwx))
( + 20160x^2, (iw)^-7*e^(iwx))
( - 40320x, (iw)^-8*e^(iwx))
( + 40320, (iw)^-9*e^(iwx))
-----------
-64*(w^-3)*sin(w)
-576*(w^-4)*cos(w)
3264*(w^-5)*sin(w)
13440*(w^-6)*cos(w)
-40320*(w^-7)*sin(w)
-80640*(w^-8)*cos(w)
80640*(w^-9)*sin(w)
-----------
result:
-64*(w^-3)*sin(w) + -576*(w^-4)*cos(w) + 3264*(w^-5)*sin(w) + 13440*(w^-6)*cos(w) + -40320*(w^-7)*sin(w) + -80640*(w^-8)*cos(w) + 80640*(w^-9)*sin(w)
-----------
a = 4
b = 3
integrating: (1 - x^4)^3 * e^iwx dx from -1 to 1
-----------
( - x^12 + 3x^8 - 3x^4 + 1, (iw)^-1*e^(iwx))
( + 12x^11 - 24x^7 + 12x^3, (iw)^-2*e^(iwx))
( - 132x^10 + 168x^6 - 36x^2, (iw)^-3*e^(iwx))
( + 1320x^9 - 1008x^5 + 72x, (iw)^-4*e^(iwx))
( - 11880x^8 + 5040x^4 - 72, (iw)^-5*e^(iwx))
( + 95040x^7 - 20160x^3, (iw)^-6*e^(iwx))
( - 665280x^6 + 60480x^2, (iw)^-7*e^(iwx))
( + 3991680x^5 - 120960x, (iw)^-8*e^(iwx))
( - 19958400x^4 + 120960, (iw)^-9*e^(iwx))
( + 79833600x^3, (iw)^-10*e^(iwx))
( - 239500800x^2, (iw)^-11*e^(iwx))
( + 479001600x, (iw)^-12*e^(iwx))
( - 479001600, (iw)^-13*e^(iwx))
-----------
768*(w^-4)*cos(w)
-13824*(w^-5)*sin(w)
-149760*(w^-6)*cos(w)
1209600*(w^-7)*sin(w)
7741440*(w^-8)*cos(w)
-39674880*(w^-9)*sin(w)
-159667200*(w^-10)*cos(w)
479001600*(w^-11)*sin(w)
958003200*(w^-12)*cos(w)
-958003200*(w^-13)*sin(w)
-----------
result:
768*(w^-4)*cos(w) + -13824*(w^-5)*sin(w) + -149760*(w^-6)*cos(w) + 1209600*(w^-7)*sin(w) + 7741440*(w^-8)*cos(w) + -39674880*(w^-9)*sin(w) + -159667200*(w^-10)*cos(w) + 479001600*(w^-11)*sin(w) + 958003200*(w^-12)*cos(w) + -958003200*(w^-13)*sin(w)
-----------
a = 4
b = 4
integrating: (1 - x^4)^4 * e^iwx dx from -1 to 1
-----------
( + x^16 - 4x^12 + 6x^8 - 4x^4 + 1, (iw)^-1*e^(iwx))
( - 16x^15 + 48x^11 - 48x^7 + 16x^3, (iw)^-2*e^(iwx))
( + 240x^14 - 528x^10 + 336x^6 - 48x^2, (iw)^-3*e^(iwx))
( - 3360x^13 + 5280x^9 - 2016x^5 + 96x, (iw)^-4*e^(iwx))
( + 43680x^12 - 47520x^8 + 10080x^4 - 96, (iw)^-5*e^(iwx))
( - 524160x^11 + 380160x^7 - 40320x^3, (iw)^-6*e^(iwx))
( + 5765760x^10 - 2661120x^6 + 120960x^2, (iw)^-7*e^(iwx))
( - 57657600x^9 + 15966720x^5 - 241920x, (iw)^-8*e^(iwx))
( + 518918400x^8 - 79833600x^4 + 241920, (iw)^-9*e^(iwx))
( - 4151347200x^7 + 319334400x^3, (iw)^-10*e^(iwx))
( + 29059430400x^6 - 958003200x^2, (iw)^-11*e^(iwx))
( - 174356582400x^5 + 1916006400x, (iw)^-12*e^(iwx))
( + 871782912000x^4 - 1916006400, (iw)^-13*e^(iwx))
( - 3487131648000x^3, (iw)^-14*e^(iwx))
( + 10461394944000x^2, (iw)^-15*e^(iwx))
( - 20922789888000x, (iw)^-16*e^(iwx))
( + 20922789888000, (iw)^-17*e^(iwx))
-----------
12288*(w^-5)*sin(w)
368640*(w^-6)*cos(w)
-6451200*(w^-7)*sin(w)
-83865600*(w^-8)*cos(w)
878653440*(w^-9)*sin(w)
7664025600*(w^-10)*cos(w)
-56202854400*(w^-11)*sin(w)
-344881152000*(w^-12)*cos(w)
1739733811200*(w^-13)*sin(w)
6974263296000*(w^-14)*cos(w)
-20922789888000*(w^-15)*sin(w)
-41845579776000*(w^-16)*cos(w)
41845579776000*(w^-17)*sin(w)
-----------
result:
12288*(w^-5)*sin(w) + 368640*(w^-6)*cos(w) + -6451200*(w^-7)*sin(w) + -83865600*(w^-8)*cos(w) + 878653440*(w^-9)*sin(w) + 7664025600*(w^-10)*cos(w) + -56202854400*(w^-11)*sin(w) + -344881152000*(w^-12)*cos(w) + 1739733811200*(w^-13)*sin(w) + 6974263296000*(w^-14)*cos(w) + -20922789888000*(w^-15)*sin(w) + -41845579776000*(w^-16)*cos(w) + 41845579776000*(w^-17)*sin(w)
-----------
a = 6
b = 1
integrating: (1 - x^6)^1 * e^iwx dx from -1 to 1
-----------
( - x^6 + 1, (iw)^-1*e^(iwx))
( + 6x^5, (iw)^-2*e^(iwx))
( - 30x^4, (iw)^-3*e^(iwx))
( + 120x^3, (iw)^-4*e^(iwx))
( - 360x^2, (iw)^-5*e^(iwx))
( + 720x, (iw)^-6*e^(iwx))
( - 720, (iw)^-7*e^(iwx))
-----------
-12*(w^-2)*cos(w)
60*(w^-3)*sin(w)
240*(w^-4)*cos(w)
-720*(w^-5)*sin(w)
-1440*(w^-6)*cos(w)
1440*(w^-7)*sin(w)
-----------
result:
-12*(w^-2)*cos(w) + 60*(w^-3)*sin(w) + 240*(w^-4)*cos(w) + -720*(w^-5)*sin(w) + -1440*(w^-6)*cos(w) + 1440*(w^-7)*sin(w)
-----------
a = 6
b = 2
integrating: (1 - x^6)^2 * e^iwx dx from -1 to 1
-----------
( + x^12 - 2x^6 + 1, (iw)^-1*e^(iwx))
( - 12x^11 + 12x^5, (iw)^-2*e^(iwx))
( + 132x^10 - 60x^4, (iw)^-3*e^(iwx))
( - 1320x^9 + 240x^3, (iw)^-4*e^(iwx))
( + 11880x^8 - 720x^2, (iw)^-5*e^(iwx))
( - 95040x^7 + 1440x, (iw)^-6*e^(iwx))
( + 665280x^6 - 1440, (iw)^-7*e^(iwx))
( - 3991680x^5, (iw)^-8*e^(iwx))
( + 19958400x^4, (iw)^-9*e^(iwx))
( - 79833600x^3, (iw)^-10*e^(iwx))
( + 239500800x^2, (iw)^-11*e^(iwx))
( - 479001600x, (iw)^-12*e^(iwx))
( + 479001600, (iw)^-13*e^(iwx))
-----------
-144*(w^-3)*sin(w)
-2160*(w^-4)*cos(w)
22320*(w^-5)*sin(w)
187200*(w^-6)*cos(w)
-1327680*(w^-7)*sin(w)
-7983360*(w^-8)*cos(w)
39916800*(w^-9)*sin(w)
159667200*(w^-10)*cos(w)
-479001600*(w^-11)*sin(w)
-958003200*(w^-12)*cos(w)
958003200*(w^-13)*sin(w)
-----------
result:
-144*(w^-3)*sin(w) + -2160*(w^-4)*cos(w) + 22320*(w^-5)*sin(w) + 187200*(w^-6)*cos(w) + -1327680*(w^-7)*sin(w) + -7983360*(w^-8)*cos(w) + 39916800*(w^-9)*sin(w) + 159667200*(w^-10)*cos(w) + -479001600*(w^-11)*sin(w) + -958003200*(w^-12)*cos(w) + 958003200*(w^-13)*sin(w)
-----------
a = 6
b = 3
integrating: (1 - x^6)^3 * e^iwx dx from -1 to 1
-----------
( - x^18 + 3x^12 - 3x^6 + 1, (iw)^-1*e^(iwx))
( + 18x^17 - 36x^11 + 18x^5, (iw)^-2*e^(iwx))
( - 306x^16 + 396x^10 - 90x^4, (iw)^-3*e^(iwx))
( + 4896x^15 - 3960x^9 + 360x^3, (iw)^-4*e^(iwx))
( - 73440x^14 + 35640x^8 - 1080x^2, (iw)^-5*e^(iwx))
( + 1028160x^13 - 285120x^7 + 2160x, (iw)^-6*e^(iwx))
( - 13366080x^12 + 1995840x^6 - 2160, (iw)^-7*e^(iwx))
( + 160392960x^11 - 11975040x^5, (iw)^-8*e^(iwx))
( - 1764322560x^10 + 59875200x^4, (iw)^-9*e^(iwx))
( + 17643225600x^9 - 239500800x^3, (iw)^-10*e^(iwx))
( - 158789030400x^8 + 718502400x^2, (iw)^-11*e^(iwx))
( + 1270312243200x^7 - 1437004800x, (iw)^-12*e^(iwx))
( - 8892185702400x^6 + 1437004800, (iw)^-13*e^(iwx))
( + 53353114214400x^5, (iw)^-14*e^(iwx))
( - 266765571072000x^4, (iw)^-15*e^(iwx))
( + 1067062284288000x^3, (iw)^-16*e^(iwx))
( - 3201186852864000x^2, (iw)^-17*e^(iwx))
( + 6402373705728000x, (iw)^-18*e^(iwx))
( - 6402373705728000, (iw)^-19*e^(iwx))
-----------
2592*(w^-4)*cos(w)
-77760*(w^-5)*sin(w)
-1490400*(w^-6)*cos(w)
22744800*(w^-7)*sin(w)
296835840*(w^-8)*cos(w)
-3408894720*(w^-9)*sin(w)
-34807449600*(w^-10)*cos(w)
316141056000*(w^-11)*sin(w)
2537750476800*(w^-12)*cos(w)
-17781497395200*(w^-13)*sin(w)
-106706228428800*(w^-14)*cos(w)
533531142144000*(w^-15)*sin(w)
2134124568576000*(w^-16)*cos(w)
-6402373705728000*(w^-17)*sin(w)
-12804747411456000*(w^-18)*cos(w)
12804747411456000*(w^-19)*sin(w)
-----------
result:
2592*(w^-4)*cos(w) + -77760*(w^-5)*sin(w) + -1490400*(w^-6)*cos(w) + 22744800*(w^-7)*sin(w) + 296835840*(w^-8)*cos(w) + -3408894720*(w^-9)*sin(w) + -34807449600*(w^-10)*cos(w) + 316141056000*(w^-11)*sin(w) + 2537750476800*(w^-12)*cos(w) + -17781497395200*(w^-13)*sin(w) + -106706228428800*(w^-14)*cos(w) + 533531142144000*(w^-15)*sin(w) + 2134124568576000*(w^-16)*cos(w) + -6402373705728000*(w^-17)*sin(w) + -12804747411456000*(w^-18)*cos(w) + 12804747411456000*(w^-19)*sin(w)
-----------
a = 6
b = 4
integrating: (1 - x^6)^4 * e^iwx dx from -1 to 1
-----------
( + x^24 - 4x^18 + 6x^12 - 4x^6 + 1, (iw)^-1*e^(iwx))
( - 24x^23 + 72x^17 - 72x^11 + 24x^5, (iw)^-2*e^(iwx))
( + 552x^22 - 1224x^16 + 792x^10 - 120x^4, (iw)^-3*e^(iwx))
( - 12144x^21 + 19584x^15 - 7920x^9 + 480x^3, (iw)^-4*e^(iwx))
( + 255024x^20 - 293760x^14 + 71280x^8 - 1440x^2, (iw)^-5*e^(iwx))
( - 5100480x^19 + 4112640x^13 - 570240x^7 + 2880x, (iw)^-6*e^(iwx))
( + 96909120x^18 - 53464320x^12 + 3991680x^6 - 2880, (iw)^-7*e^(iwx))
( - 1744364160x^17 + 641571840x^11 - 23950080x^5, (iw)^-8*e^(iwx))
( + 29654190720x^16 - 7057290240x^10 + 119750400x^4, (iw)^-9*e^(iwx))
( - 474467051520x^15 + 70572902400x^9 - 479001600x^3, (iw)^-10*e^(iwx))
( + 7117005772800x^14 - 635156121600x^8 + 1437004800x^2, (iw)^-11*e^(iwx))
( - 99638080819200x^13 + 5081248972800x^7 - 2874009600x, (iw)^-12*e^(iwx))
( + 1295295050649600x^12 - 35568742809600x^6 + 2874009600, (iw)^-13*e^(iwx))
( - 15543540607795200x^11 + 213412456857600x^5, (iw)^-14*e^(iwx))
( + 170978946685747200x^10 - 1067062284288000x^4, (iw)^-15*e^(iwx))
( - 1709789466857472000x^9 + 4268249137152000x^3, (iw)^-16*e^(iwx))
( + 15388105201717248000x^8 - 12804747411456000x^2, (iw)^-17*e^(iwx))
( - 123104841613737984000x^7 + 25609494822912000x, (iw)^-18*e^(iwx))
( + 861733891296165888000x^6 - 25609494822912000, (iw)^-19*e^(iwx))
( - 5170403347776995328000x^5, (iw)^-20*e^(iwx))
( + 25852016738884976640000x^4, (iw)^-21*e^(iwx))
( - 103408066955539906560000x^3, (iw)^-22*e^(iwx))
( + 310224200866619719680000x^2, (iw)^-23*e^(iwx))
( - 620448401733239439360000x, (iw)^-24*e^(iwx))
( + 620448401733239439360000, (iw)^-25*e^(iwx))
-----------
62208*(w^-5)*sin(w)
3110400*(w^-6)*cos(w)
-94867200*(w^-7)*sin(w)
-2253484800*(w^-8)*cos(w)
45433301760*(w^-9)*sin(w)
808746301440*(w^-10)*cos(w)
-12966573312000*(w^-11)*sin(w)
-189119411712000*(w^-12)*cos(w)
2519458363699200*(w^-13)*sin(w)
30660256301875200*(w^-14)*cos(w)
-339823768802918400*(w^-15)*sin(w)
-3411042435440640000*(w^-16)*cos(w)
30750600908611584000*(w^-17)*sin(w)
246158464237830144000*(w^-18)*cos(w)
-1723416563602685952000*(w^-19)*sin(w)
-10340806695553990656000*(w^-20)*cos(w)
51704033477769953280000*(w^-21)*sin(w)
206816133911079813120000*(w^-22)*cos(w)
-620448401733239439360000*(w^-23)*sin(w)
-1240896803466478878720000*(w^-24)*cos(w)
1240896803466478878720000*(w^-25)*sin(w)
-----------
result:
62208*(w^-5)*sin(w) + 3110400*(w^-6)*cos(w) + -94867200*(w^-7)*sin(w) + -2253484800*(w^-8)*cos(w) + 45433301760*(w^-9)*sin(w) + 808746301440*(w^-10)*cos(w) + -12966573312000*(w^-11)*sin(w) + -189119411712000*(w^-12)*cos(w) + 2519458363699200*(w^-13)*sin(w) + 30660256301875200*(w^-14)*cos(w) + -339823768802918400*(w^-15)*sin(w) + -3411042435440640000*(w^-16)*cos(w) + 30750600908611584000*(w^-17)*sin(w) + 246158464237830144000*(w^-18)*cos(w) + -1723416563602685952000*(w^-19)*sin(w) + -10340806695553990656000*(w^-20)*cos(w) + 51704033477769953280000*(w^-21)*sin(w) + 206816133911079813120000*(w^-22)*cos(w) + -620448401733239439360000*(w^-23)*sin(w) + -1240896803466478878720000*(w^-24)*cos(w) + 1240896803466478878720000*(w^-25)*sin(w)

"""