bridge/class-master/depend/pari/src/basemath/rootpol.c
/********************************************************************/
/**                                                                **/
/**              FACTORS THROUGH CIRCLE INTEGRATION                **/
/**                                                                **/
/********************************************************************/
/* l power of 2 */
FACTORS THROUGH CIRCLE INTEGRATION

l power of 2
FACTORS THROUGH CIRCLE INTEGRATION l power of 2
static void
fft(GEN Omega, GEN p, GEN f, long step, long l)
{
  pari_sp ltop;
  long i, l1, l2, l3, rapi, step4;
  GEN f1, f2, f3, f02, f13, g02, g13, ff;
  if (l == 2)
  {
    gel(f,0) = gadd(gel(p,0),gel(p,step));
    gel(f,1) = gsub(gel(p,0),gel(p,step)); return;
  }
  if (l == 4)
  {
    f1 = gadd(gel(p,0),   gel(p,step<<1));
    f2 = gsub(gel(p,0),   gel(p,step<<1));
    f3 = gadd(gel(p,step),gel(p,3*step));
    f02= gsub(gel(p,step),gel(p,3*step));
    f02 = mulcxI(f02);
    gel(f,0) = gadd(f1, f3);
    gel(f,1) = gadd(f2, f02);
    gel(f,2) = gsub(f1, f3);
    gel(f,3) = gsub(f2, f02); return;
  }
  ltop = avma;
  l1 = l>>2; l2 = 2*l1; l3 = l1+l2; step4 = step<<2;
  fft(Omega,p,          f,   step4,l1);
  fft(Omega,p+step,     f+l1,step4,l1);
  fft(Omega,p+(step<<1),f+l2,step4,l1);
  fft(Omega,p+3*step,   f+l3,step4,l1);
  ff = cgetg(l+1,t_VEC);
  for (i=0; i<l1; i++)
  {
    rapi = step*i;
    f1 = gmul(gel(Omega,rapi),    gel(f,i+l1));
    f2 = gmul(gel(Omega,rapi<<1), gel(f,i+l2));
    f3 = gmul(gel(Omega,3*rapi),  gel(f,i+l3));
    f02 = gadd(gel(f,i),f2);
    g02 = gsub(gel(f,i),f2);
    f13 = gadd(f1,f3);
    g13 = mulcxI(gsub(f1,f3));
    gel(ff,i+1)    = gadd(f02, f13);
    gel(ff,i+l1+1) = gadd(g02, g13);
    gel(ff,i+l2+1) = gsub(f02, f13);
    gel(ff,i+l3+1) = gsub(g02, g13);
  }
  ff = gerepilecopy(ltop,ff);
  for (i=0; i<l; i++) f[i] = ff[i+1];
}
static void fft ( GEN Omega , GEN p , GEN f , long step , long l ) { pari_sp ltop ; long i , l1 , l2 , l3 , rapi , step4 ; GEN f1 , f2 , f3 , f02 , f13 , g02 , g13 , ff ; if ( l == 2 ) { gel ( f , 0 ) = gadd ( gel ( p , 0 ) , gel ( p , step ) ) ; gel ( f , 1 ) = gsub ( gel ( p , 0 ) , gel ( p , step ) ) ; return ; } if ( l == 4 ) { f1 = gadd ( gel ( p , 0 ) , gel ( p , step << 1 ) ) ; f2 = gsub ( gel ( p , 0 ) , gel ( p , step << 1 ) ) ; f3 = gadd ( gel ( p , step ) , gel ( p , 3 * step ) ) ; f02 = gsub ( gel ( p , step ) , gel ( p , 3 * step ) ) ; f02 = mulcxI ( f02 ) ; gel ( f , 0 ) = gadd ( f1 , f3 ) ; gel ( f , 1 ) = gadd ( f2 , f02 ) ; gel ( f , 2 ) = gsub ( f1 , f3 ) ; gel ( f , 3 ) = gsub ( f2 , f02 ) ; return ; } ltop = avma ; l1 = l >> 2 ; l2 = 2 * l1 ; l3 = l1 + l2 ; step4 = step << 2 ; fft ( Omega , p , f , step4 , l1 ) ; fft ( Omega , p + step , f + l1 , step4 , l1 ) ; fft ( Omega , p + ( step << 1 ) , f + l2 , step4 , l1 ) ; fft ( Omega , p + 3 * step , f + l3 , step4 , l1 ) ; ff = cgetg ( l + 1 , t_VEC ) ; for ( i = 0 ; i < l1 ; i ++ ) { rapi = step * i ; f1 = gmul ( gel ( Omega , rapi ) , gel ( f , i + l1 ) ) ; f2 = gmul ( gel ( Omega , rapi << 1 ) , gel ( f , i + l2 ) ) ; f3 = gmul ( gel ( Omega , 3 * rapi ) , gel ( f , i + l3 ) ) ; f02 = gadd ( gel ( f , i ) , f2 ) ; g02 = gsub ( gel ( f , i ) , f2 ) ; f13 = gadd ( f1 , f3 ) ; g13 = mulcxI ( gsub ( f1 , f3 ) ) ; gel ( ff , i + 1 ) = gadd ( f02 , f13 ) ; gel ( ff , i + l1 + 1 ) = gadd ( g02 , g13 ) ; gel ( ff , i + l2 + 1 ) = gsub ( f02 , f13 ) ; gel ( ff , i + l3 + 1 ) = gsub ( g02 , g13 ) ; } ff = gerepilecopy ( ltop , ff ) ; for ( i = 0 ; i < l ; i ++ ) f [ i ] = ff [ i + 1 ] ; }
FACTORS THROUGH CIRCLE INTEGRATION
FACTORS THROUGH CIRCLE INTEGRATION