我正在尝试编写一个算法,在给定
a(0),..., a(n-1)
的值的情况下找到 n, x_1, ..., x_n, a(n)
,这样:
a(n)*p^n + a(n-1)*p^(n-1) + ... + a(1)*p + a(0) = a(n)(p-x_1)(p-x_2)...(p-x_n)
对于所有真实的p。
乘以 a(n)(p-x_1)(p-x_2) 之后,我想到使用 Viete 公式来查找系数。
但事实证明,把代码写下来并不像我想象的那么明显。
我只想在代码中使用基础知识 - 即循环、if-s 加法和乘法 - 没有现成/复杂的函数。
首先,我想强调的是,我只需要一个伪代码,我不关心为根和系数定义数组。这就是为什么我只写a(n), xn。哦,我希望如果我从 i=1 而不是 i=0 开始索引以便与数学符号保持同步,这不会太打扰您。为了从 i=0 开始,我必须重新枚举根并引入更多括号。
这就是我到目前为止所想到的:
a(n-1)=0;
for(i=1; i <= n; i++){
a(n-1) = a(n-1) + x_i;
}
a(n-1) = -a(n)*a(n-1);
a(n-2)=0;
for(i=1; i <= n; i++){
for(j=i; j <= n; j++){
a(n-2) = a(n-2)+ x_i * x_j;
}
}
a(n-2) = -a(n)*a(n-2);
a(n-3)=0;
for(i=1; i <= n; i++){
for(j=i; j <= n; j++){
for(k=j; k <= n; k++){
a(n-3) = a(n-3)+ x_i * x_j * x_k;
}
}
}
a(n-3) = a(n)*a(n-3);
...
a(0)=1;
for(i=1; i<=n; i++){
a(0) = a(0) * x_i;
}
if(n%2 == 0) a(0) = a(n) * a(0);
else a(0) = -a(n) * a(0);
如你所见,看起来不太好。
我想将所有这些循环链接到一个循环中,因为如果没有,我无法编写完整的代码,我无法填补固定 j 的 a(0) 和 a(n-j) 之间的空白。
你能帮我吗?
这就是我所拥有的,基于 Nico Schertler 的回答:
for(i=1; i<=n; i++)
{a(i)=1;
for(j=1; j <= n; j++)
{b(i)= clone( a(i) );
a(i) = a(i-1);
b(i) = x_j * b(i);
c(i) = a(i) - b(i);
}
}
如果我们这样写,会一样吗
for(i=1; i<=n; i++)
{a(i)=1; b(i)=1;
for(j=1; j <= n; j++)
{t = a(i) ;
a(i) = a(i-1);
b(i) = x_j * t;
c(i) = a(i) - b(i);
}
}
(这就是我们如何通过将 a[i] 的值保留在某个变量 t 中来交换数组的两个元素)。
您可以增量创建多项式。
从
p = 1
开始。 IE。 a(0) = 1
。
为了添加根,您必须将当前多项式乘以
x - x_i
。这是:
p * (x - x_i) = p * x - p * x_i
所以需要支持三种操作:
这很简单。只需将所有系数向左移动一位即可。即
a(i ) := a(i - 1)
a(i - 1) := a(i - 2)
...
a(1 ) := a(0)
a(0 ) := 0
这同样简单。将每个系数相乘:
a(i ) *= s
a(i - 1) *= s
...
只需减去各自的系数即可:
c(i ) = a(i ) - b(i )
c(i - 1) = a(i - 1) - b(i - 1)
...
逐根添加根。首先,克隆当前的多项式。然后,进行上述操作:
p := 1
for each root r
p' = clone(p)
multiply p with x
multiply p' with r
p := p - p'
next
C# 中用于此目的的静态函数。 x^4-11x^3+44x^2-76x+48 的根是 {2,2,3,4} 并给定参数
roots = new Complex[4] {2, 2, 3, 4}
此函数返回 [48,-76,44,-11,1]
public static double[] FromRoots(Complex[] roots)
{
int N = roots.Length;
Complex[] coefs = new Complex[N + 1];
coefs[0] = -roots[0];
coefs[1] = 1.0;
for (int k = 2; k <= N; k++)
{
coefs[k] = 1.0;
for (int i = k - 2; i >= 0; i--)
{
coefs[i + 1] = coefs[i] - roots[k - 1] * coefs[i + 1];
}
coefs[0] *= -roots[k - 1];
if (Math.IEEERemainder(k, 2) == 1)
coefs[k] = -coefs[k];
}
double[] realCoefs = new double[N + 1];
for (int i = 0; i < N + 1; i++)
realCoefs[i] = coefs[i].Real; // Not sure about this part!
return realCoefs;
}
我想向 PDL 添加
polyfromroots
,因此需要一个 C 版本。
Octave 中的算法(撰写本文时的 GitHub 版本,带有我的评论):
y = zeros (1, n+1);
y(1) = 1;
for j = 1:n
# this invisibly creates temporary copy before assign
y(2:(j+1)) -= v(j) .* y(1:j);
endfor
翻译为C99:
/* pc must point to a block of degree+1 complex doubles */
void polyfromroots(complex double root[], int degree, complex double coeff[])
{
int j, k;
coeff[0] = 1.0;
for (j = 1; j <= degree; j++) coeff[j] = 0.0;
for (j = 0; j < degree; j++)
for (k = j; k >= 0; k--) /* count down to use data before we mutate */
coeff[k+1] -= root[j] * coeff[k];
}