私は、カラツバ法を使用して多項式乗算の単純な分割統治アルゴリズムを実装しようとしています。つまり、、、を使用して、p=a+b*x^kまたはの次数の半分近くである、だけを再帰的に計算します。q=c+d*x^kp*q=ac+(ad+bc)x^k+bd*x^(2k)ad+bc=ac+bd-(a-b)*(c-d)ac,bd,(a-b)(c-d)kpq
次のコードは、次数が64以上の多項式（0〜9のランダムに生成された整数係数）を2乗するために使用すると、エラーメッセージなしでクラッシュしますが、それよりも小さい次数では機能するようです。
追加:(プログラムは正しく終了しません。また、次数64の場合、約3 ^（log（64））= 729の再帰関数呼び出しのみを使用する必要があるため、コードが機能すれば、そこからのスタックオーバーフローを除外できると思います。その点で意図されているように。）
単独で使用するブルートフォース関数は、大規模な場合でも正常に機能するようです。
struct poly {
int deg;
double* coeff;
};
poly stdpolmult(poly p,poly q) { // standard algorithm
poly r;
r.deg= p.deg+q.deg;
r.coeff = (double*) calloc (r.deg,sizeof(double));
int i,j;
for (i=0;i<=p.deg;i++)
for (j=0;j<=q.deg;j++)
r.coeff[i+j]=r.coeff[i+j]+p.coeff[i]*q.coeff[j];
return r;
}
poly fastpolmult(poly p,poly q) { // Divide & Conquer
if ((p.deg<=7)&&(q.deg<=7))
return stdpolmult(p,q); // brute force
poly a,b,c,d,u,v,x,y,w,z,s,r;
int k=p.deg/2;
if (p.deg<q.deg)
k=q.deg/2;
a=polfirstpart(p,k);
b=pollastpart(p,k);
c=polfirstpart(q,k);
d=pollastpart(q,k); /* let p=p_1+x^k*p_2, q=q_1+x^k+q_2, then
a= p_1,b=p_2,c=q_1,d=q_2 */
u = fastpolmult(a,c); // u =p_1*p_2
v = fastpolmult(b,d); // v =q_1*q_2
polneg(b); // b= -p_2
polneg(d); // d= -q_2
x=poladd(a,b); // x=p_1-p_2
y=poladd(c,d); // y=q_1-q_2
w=fastpolmult(x,y); // w=(p_1-p_2)*(q_1-q_2)
polneg(w); // w= -(p_1-p_2)*(q_1-q_2)
z=poladd(u,v); // z=p_1*p_2+q_1*q_2
s=poladd(z,w); // s=p_1*p_2+q_1*q_2-(p_1-p_2)*(q_1-q_2) = p_1*q_2+p_2*q_1
polfree(z); polfree(w); polfree(x); polfree(y);
x=polshift(s,k); // x=(p_1*q_2+p_2*q_1)*x^k
y=polshift(v,2*k); // y=q_1*q_2*x^(2k)
z=poladd(u,x); // z=p_1*p_2+(p_1*q_2+p_2*q_1)*x^k
r=poladd(z,y); // r = p_1*p_2+(p_1*q_2+p_2*q_1)*x^k +q_1*q_2*x^(2k) = p*q
polfree(x);polfree(y);polfree(z);
return r;
}
必要に応じて、使用する関数（polfirstpart、pollastpart、polneg、poladd、polshift、polfree）のコードを追加できます。それらは特別なものではありません。（そして私はそれらを個別にテストしました、それらは機能しているようでした）。
追加：これらの関数のほとんどのコード：
poly poladd(poly p,poly q) {
int i,n;
poly r;
if (p.deg>=q.deg)
r.deg=p.deg;
else
r.deg=q.deg;
r.coeff = (double*) calloc (r.deg+1,sizeof(double));
if (p.deg<=q.deg)
n=p.deg;
else
n=q.deg;
for (i=0;i<=n;i++)
r.coeff[i]=p.coeff[i]+q.coeff[i];
if (p.deg>q.deg)
for (i=n+1;i<=p.deg;i++)
r.coeff[i]=p.coeff[i];
else
for (i=n+1;i<=q.deg;i++)
r.coeff[i]=q.coeff[i];
return r;
}
poly polfirstpart(poly p, int k) { /* if p=a+x^k*b, take a */
poly r;
if (k<=0)
return zpol; // zero pol
if (k>p.deg)
return p;
r.coeff=p.coeff;
r.deg=k-1;
return r;
}
poly pollastpart(poly p,int k) { /* if p=a+x^k*b, take b */
poly r;
if (k<0)
return zpol;
if (k>p.deg)
return zpol;
r.coeff=(p.coeff)+k;
r.deg=p.deg-k;
return r;
}
poly polshift(poly p,int k) { /* x^k*p */
int i;
poly r;
r.deg=p.deg+k;
r.coeff= (double*) calloc(r.deg+1,sizeof(double));
for (i=0;i<=p.deg;i++)
r.coeff[k+i]=p.coeff[i];
return r;
}
void genzpol() { // called in main
zpol.deg=0;
zpol.coeff=(double*) calloc(1,sizeof(double));
}