Быстрый полиномиальный сдвиг - PullRequest
0 голосов
/ 05 июля 2018

Я пытаюсь реализовать "F. Метод свертки" (раздел 2.2):

enter image description 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))

Но я получаю только нули, тогда как правильный результат должен быть:

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

Пожалуйста, кто-нибудь знает, что я здесь не так делаю?

Ответы [ 2 ]

0 голосов
/ 06 июля 2018

Поскольку вы просили мои реализации (у них f "назад"):

Уравнение 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})

Уравнение в 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]
0 голосов
/ 05 июля 2018

Я до сих пор не до конца разбираюсь в алгоритме, но некоторые ошибки у вас были:

  1. Полномочия u начинаются с n и заканчиваются 0. Чтобы ваша свертка работала, вам нужно повернуть ее вспять, поскольку вы ожидаете, что коэффициенты будут в порядке в функции свертки.
  2. Коэффициенты в полиноме v зависят только от j, а не n-j (вы используете i)
  3. Необходимы только первые n+1 элементы свертки (вам не нужны полномочия n+1 ... 2n.
  4. Результирующая свертка (разве это не свертка, не так ли?) Является «задом наперед», так как в вашем конечном результате вы начнете вычисление с i=0, поэтому мощность x**(n-i=n).

Собираем все это вместе:

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))

Я получаю

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

Я не знаю, почему это отличается от ваших ожиданий, но я надеюсь, что это поможет вам.

...