Finding an error in code for divide and conquer polynomial multiplication

Question

I'm trying to implement a simple divide and conquer algorithm for polynomial multiplication using Karatsuba's method, that is using that for p=a+b*x^k, q=c+d*x^k, with p*q=ac+(ad+bc)x^k+bd*x^(2k), ad+bc=ac+bd-(a-b)*(c-d) and recursively computing just ac,bd,(a-b)(c-d), where k is somewhere near half of the degree of p or q.

The following code crashes without an error message when used to square a polynomial (randomly generated integer coefficients between 0 and 9) of degree >= around 64, but seems to work for smaller degrees.

Addition: (The program doesn't terminate properly. Also, for degree 64 only approximately 3^(log(64))=729 recursive function calls should be used, so I think stack overflow from that can be ruled out if the code works as intended in that regard.)

The brute force function used alone seems to work fine also for large degrees.

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;
 }

I can add my code for the functions (polfirstpart,pollastpart, polneg, poladd,polshift,polfree) which are used, if required. They are nothing special. (And I tested them individually, they seemed to work).

Addition: The code for most of these functions:

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