Fast polynomial shift

Question

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

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?

Post the algorithm here to make it easier to help. As it stands, you are dividing small integers by big integers (and that's the reason for the 0s), so we have to understand what you expected `convolve(u,v)` to return. — kabanus, Jul 05 '18 at 19:34
Are you sure that's the correct output, from (2) in isaac97 (as well as my implementation of 2.2(F), I get `{0: 23, 1: 132, 2: 396, 3: 720, 4: 840, 5: 636, 6: 301, 7: 80, 8: 9})` — jedwards, Jul 05 '18 at 23:52
What @jedwards got is also what I got after fixing the code (see below,reversed). Double check. — kabanus, Jul 06 '18 at 04:36
@jedwards True, I mixed the order of x's powers up. Is it possible to see your implementation somewhere or does it look more or less like kabanus' answer below? — Ecir Hana, Jul 06 '18 at 06:35
@EcirHana I added it both, but yeah basically kabanus' (and using a `defaultdict` so I could wrap my head around all the index juggling). Since kabanus' answer does a great job of explaining the issues in your original code, mine is basically just the code you asked for. — jedwards, Jul 06 '18 at 10:06

score 3 · Accepted Answer · edited Jul 06 '18 at 08:18

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

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.
Coefficients in the v polynomial only depend on j, not n-j (you use i)
Only the first n+1 elements of the convolution are needed (you do not need the powers of n+1...2n.
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.

Looks like the OP had his list reversed (from what we were expecting at least). If `f = f[::-1]` your code should get their intended output (mine did). — jedwards, Jul 06 '18 at 06:09
@jedwards I agree, I think this is the correct version. (Note I reverse on purpose inside f). — kabanus, Jul 06 '18 at 06:24

score 1 · Answer 2 · answered Jul 06 '18 at 10:05

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]

Fast polynomial shift

2 Answers2