快速多项式移位
Fast polynomial shift
我正在尝试实施 "F. The convolution method"(第 2.2 节):
来自 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))
但我只得到零,而正确的结果应该是:
[1, 10, 45, 112, 170, 172, 116, 52, 23]
拜托,有人知道我在这里做错了什么吗?
我还没有完全掌握算法,但是你有一些错误:
u
的权力从n
开始,到0
结束。为了使卷积起作用,您需要将其反转,因为您希望系数在卷积函数中按顺序排列。
v
多项式中的系数只取决于j
,而不是n-j
(你用i
)
- 只需要卷积的前
n+1
个元素(您不需要 n+1
...2n
. 的幂
- 最终的卷积(它不是真正的卷积,对吧?)是 "backwards",因为你的最终结果将从
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]
我不知道为什么这与您的预期不同,但我希望这能让您走上正轨。
既然你问了我的实现(这些有 f
"backwards"):
等式 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]
我正在尝试实施 "F. The convolution method"(第 2.2 节):
来自 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))
但我只得到零,而正确的结果应该是:
[1, 10, 45, 112, 170, 172, 116, 52, 23]
拜托,有人知道我在这里做错了什么吗?
我还没有完全掌握算法,但是你有一些错误:
u
的权力从n
开始,到0
结束。为了使卷积起作用,您需要将其反转,因为您希望系数在卷积函数中按顺序排列。v
多项式中的系数只取决于j
,而不是n-j
(你用i
)- 只需要卷积的前
n+1
个元素(您不需要n+1
...2n
. 的幂
- 最终的卷积(它不是真正的卷积,对吧?)是 "backwards",因为你的最终结果将从
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]
我不知道为什么这与您的预期不同,但我希望这能让您走上正轨。
既然你问了我的实现(这些有 f
"backwards"):
等式 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]