CODEWITHSUNDEEP
pythonalgorithmmathpolynomialstaylor-series
Ecir Hana
Ecir Hana

Reputation: 11508

Fast polynomial shift

I'm trying to implement "F. The convolution method" (section 2.2):

enter image description here

from Fast algorithms for Taylor shifts and certain difference equations (at the bottom, or here):

from math import factorial

def convolve(f, h):
    g = [0] * (len(f) + len(h) - 1)
    for hindex, hval in enumerate(h):
        for findex, fval in enumerate(f):
            g[hindex + findex] += fval * hval
    return g

def shift(f, a):
    n = len(f) - 1
    u = [factorial(i)*c for i, c in enumerate(f)]
    v = [factorial(n)*a**i//factorial(n-i) for i in range(n + 1)]
    g = convolve(u, v)
    g = [c//(factorial(n)*factorial(i)) for i, c in enumerate(g)]
    return g

f = [1, 2, 3, -4, 5, 6, -7, 8, 9]
print(shift(f, 1))

But I get only zeros, whereas the correct result should be:

[1, 10, 45, 112, 170, 172, 116, 52, 23]

Please, does anyone know what am I doing wrong here?

Upvotes: 4

Views: 774

Answers (2)

jedwards
jedwards

Reputation: 30240

Since you asked for my implementations (these have f "backwards"):

Equation 2:

from math import factorial
from collections import defaultdict

def binomial(n, k):
    try:
        binom = factorial(n) // factorial(k) // factorial(n - k)
    except ValueError:
        binom = 0
    return binom

f = [1, 2, 3, -4, 5, 6, -7, 8, 9][::-1]
k=0
n = len(f) - 1

g = defaultdict(int)
for k in range(n+1):
    for i in range(k, n+1):
        g[i-k] += binomial(i,k) * f[i]

print(g)
# defaultdict(<class 'int'>, {0: 23, 1: 52, 2: 116, 3: 172, 4: 170, 5: 112, 6: 45, 7: 10, 8: 1})

Equation in 2.2(F):

from math import factorial
from collections import defaultdict

def convolve(x, y):
    g = defaultdict(int)
    for (xi, xv) in x.items():
        for (yi, yv) in y.items():
            g[xi + yi] += xv * yv
    return g


def shift(f, a):
    n = len(f) - 1
    u = {n-i: factorial(i)*c for (i, c) in enumerate(f)}
    v = {j: factorial(n)*a**j//factorial(j) for j in range(n + 1)}
    uv = convolve(u, v)

    def g(k):
        ngk = uv[n-k]
        return ngk // factorial(n) // factorial(k)

    G = [g(k) for k in range(n+1)]
    return G

f = [1, 2, 3, -4, 5, 6, -7, 8, 9]

print(shift(f, 1))
# [23, 132, 396, 720, 840, 636, 301, 80, 9]

Upvotes: 1

kabanus
kabanus

Reputation: 25980

I still do not fully grasp the algorithm, but some mistakes you had:

  1. The powers of u start at n and end at 0. For your convolution to work you need to reverse it as you are expecting the coefficients to be in order in the convolution function.
  2. Coefficients in the v polynomial only depend on j, not n-j (you use i)
  3. Only the first n+1 elements of the convolution are needed (you do not need the powers of n+1...2n.
  4. The resulting convolution (it's not really a convolution is it?) is "backwards", as in you final result will start the calculation at i=0, so the power of x**(n-i=n).

Putting all of these together:

from math import factorial

def convolve(f, h):
    g = [0] * (len(f) + len(h) - 1)
    for hindex, hval in enumerate(h):
        for findex, fval in enumerate(f):
            g[hindex + findex] += fval * hval
    return g

def shift(f, a):
    n = len(f) - 1
    u = [factorial(i)*c for i, c in enumerate(f)][::-1]
    v = [factorial(n)*a**i//factorial(i) for i in range(n + 1)]
    g = convolve(u, v)
    g = [g[n-i]//(factorial(n)*factorial(i)) for i in range(n+1)][::-1]
    return g

f = [1, 2, 3, -4, 5, 6, -7, 8, 9]
print(shift(f, 1))

I get

[9, 80, 301, 636, 840, 720, 396, 132, 23]

I do not know why this is different than your expectation, but I hope this gets you on track.

Upvotes: 3

Related Questions