Нахождение ошибки в коде для умножения полинома «разделяй и властвуй»

Я пытаюсь реализовать простой алгоритм «разделяй и властвуй» для полиномиального умножения, используя метод Карацубы, который использует его для p=a+b*x^k, q=c+d*x^k, с p*q=ac+(ad+bc)x^k+bd*x^(2k), ad+bc=ac+bd-(a-b)*(c-d) и рекурсивно вычисляя просто ac,bd,(a-b)(c-d), где k где-то около половины степени p или же q,

В следующем коде происходит сбой без сообщения об ошибке, когда он используется для возведения в квадрат многочлена (случайным образом сгенерированных целочисленных коэффициентов от 0 до 9) степени> = около 64, но, похоже, работает для меньших степеней.

Дополнение: (Программа не завершается должным образом. Кроме того, для степени 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));
}