bridge/class-master/depend/pari/src/basemath/alglin2.c
/*******************************************************************/
/*                                                                 */
/*        CHARACTERISTIC POLYNOMIAL (BERKOWITZ'S ALGORITHM)        */
/*                                                                 */
/*******************************************************************/
CHARACTERISTIC POLYNOMIAL (BERKOWITZ'S ALGORITHM)
CHARACTERISTIC POLYNOMIAL ( BERKOWITZ ' S ALGORITHM )
GEN
carberkowitz(GEN x, long v)
{
  long lx, i, j, k, r;
  GEN V, S, C, Q;
  pari_sp av0, av;
  if ((V = easychar(x,v))) return V;
  lx = lg(x); av0 = avma;
  V = cgetg(lx+1, t_VEC);
  S = cgetg(lx+1, t_VEC);
  C = cgetg(lx+1, t_VEC);
  Q = cgetg(lx+1, t_VEC);
  av = avma;
  gel(C,1) = gen_m1;
  gel(V,1) = gen_m1;
  for (i=2;i<=lx; i++) gel(C,i) = gel(Q,i) = gel(S,i) = gel(V,i) = gen_0;
  gel(V,2) = gcoeff(x,1,1);
  for (r = 2; r < lx; r++)
  {
    pari_sp av2;
    GEN t;
    for (i = 1; i < r; i++) gel(S,i) = gcoeff(x,i,r);
    gel(C,2) = gcoeff(x,r,r);
    for (i = 1; i < r-1; i++)
    {
      av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
      for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
      gel(C,i+2) = gerepileupto(av2, t);
      for (j = 1; j < r; j++)
      {
        av2 = avma; t = gmul(gcoeff(x,j,1), gel(S,1));
        for (k = 2; k < r; k++) t = gadd(t, gmul(gcoeff(x,j,k), gel(S,k)));
        gel(Q,j) = gerepileupto(av2, t);
      }
      for (j = 1; j < r; j++) gel(S,j) = gel(Q,j);
    }
    av2 = avma; t = gmul(gcoeff(x,r,1), gel(S,1));
    for (j = 2; j < r; j++) t = gadd(t, gmul(gcoeff(x,r,j), gel(S,j)));
    gel(C,r+1) = gerepileupto(av2, t);
    if (gc_needed(av0,1))
    {
      if (DEBUGMEM>1) pari_warn(warnmem,"carberkowitz");
      gerepileall(av, 2, &C, &V);
    }
    for (i = 1; i <= r+1; i++)
    {
      av2 = avma; t = gmul(gel(C,i), gel(V,1));
      for (j = 2; j <= minss(r,i); j++)
        t = gadd(t, gmul(gel(C,i+1-j), gel(V,j)));
      gel(Q,i) = gerepileupto(av2, t);
    }
    for (i = 1; i <= r+1; i++) gel(V,i) = gel(Q,i);
  }
  V = RgV_to_RgX(vecreverse(V), v); 
  V = odd(lx)? gcopy(V): RgX_neg(V);
  return fix_pol(av0, V);
}
GEN carberkowitz ( GEN x , long v ) { long lx , i , j , k , r ; GEN V , S , C , Q ; pari_sp av0 , av ; if ( ( V = easychar ( x , v ) ) ) return V ; lx = lg ( x ) ; av0 = avma ; V = cgetg ( lx + 1 , t_VEC ) ; S = cgetg ( lx + 1 , t_VEC ) ; C = cgetg ( lx + 1 , t_VEC ) ; Q = cgetg ( lx + 1 , t_VEC ) ; av = avma ; gel ( C , 1 ) = gen_m1 ; gel ( V , 1 ) = gen_m1 ; for ( i = 2 ; i <= lx ; i ++ ) gel ( C , i ) = gel ( Q , i ) = gel ( S , i ) = gel ( V , i ) = gen_0 ; gel ( V , 2 ) = gcoeff ( x , 1 , 1 ) ; for ( r = 2 ; r < lx ; r ++ ) { pari_sp av2 ; GEN t ; for ( i = 1 ; i < r ; i ++ ) gel ( S , i ) = gcoeff ( x , i , r ) ; gel ( C , 2 ) = gcoeff ( x , r , r ) ; for ( i = 1 ; i < r - 1 ; i ++ ) { av2 = avma ; t = gmul ( gcoeff ( x , r , 1 ) , gel ( S , 1 ) ) ; for ( j = 2 ; j < r ; j ++ ) t = gadd ( t , gmul ( gcoeff ( x , r , j ) , gel ( S , j ) ) ) ; gel ( C , i + 2 ) = gerepileupto ( av2 , t ) ; for ( j = 1 ; j < r ; j ++ ) { av2 = avma ; t = gmul ( gcoeff ( x , j , 1 ) , gel ( S , 1 ) ) ; for ( k = 2 ; k < r ; k ++ ) t = gadd ( t , gmul ( gcoeff ( x , j , k ) , gel ( S , k ) ) ) ; gel ( Q , j ) = gerepileupto ( av2 , t ) ; } for ( j = 1 ; j < r ; j ++ ) gel ( S , j ) = gel ( Q , j ) ; } av2 = avma ; t = gmul ( gcoeff ( x , r , 1 ) , gel ( S , 1 ) ) ; for ( j = 2 ; j < r ; j ++ ) t = gadd ( t , gmul ( gcoeff ( x , r , j ) , gel ( S , j ) ) ) ; gel ( C , r + 1 ) = gerepileupto ( av2 , t ) ; if ( gc_needed ( av0 , 1 ) ) { if ( DEBUGMEM > 1 ) pari_warn ( warnmem , " " ) ; gerepileall ( av , 2 , & C , & V ) ; } for ( i = 1 ; i <= r + 1 ; i ++ ) { av2 = avma ; t = gmul ( gel ( C , i ) , gel ( V , 1 ) ) ; for ( j = 2 ; j <= minss ( r , i ) ; j ++ ) t = gadd ( t , gmul ( gel ( C , i + 1 - j ) , gel ( V , j ) ) ) ; gel ( Q , i ) = gerepileupto ( av2 , t ) ; } for ( i = 1 ; i <= r + 1 ; i ++ ) gel ( V , i ) = gel ( Q , i ) ; } V = RgV_to_RgX ( vecreverse ( V ) , v ) ; V = odd ( lx ) ? gcopy ( V ) : RgX_neg ( V ) ; return fix_pol ( av0 , V ) ; }
CHARACTERISTIC POLYNOMIAL (BERKOWITZ'S ALGORITHM)
CHARACTERISTIC POLYNOMIAL ( BERKOWITZ ' S ALGORITHM )
/* not gtopoly: fail if v > gvar(V) */