bridge/class-master/depend/pari/src/basemath/mftrace.c
/* If dimlim > 0, keep only the dimension <= dimlim eigenspaces.
 * See mfsplit for the meaning of flag. */
If dimlim > 0, keep only the dimension <= dimlim eigenspaces.
See mfsplit for the meaning of flag.
If dimlim > 0 keep only the dimension < = dimlim eigenspaces . See mfsplit for the meaning of flag .
static GEN
split_ii(GEN mf, long dimlim, long flag, long *pnewd)
{
  forprime_t iter;
  GEN CHI = MF_get_CHI(mf), empty = cgetg(1, t_VEC), mf0 = mf;
  GEN NF, POLCYC, todosp, Tpbigvec, simplesp;
  long N = MF_get_N(mf), k = MF_get_k(mf);
  long ord, FC, NEWT, dimsimple = 0, newd = -1;
  const long NBH = 5, vz = 1;
  ulong p;
  switch(MF_get_space(mf))
  {
    case mf_NEW: break;
    case mf_CUSP:
    case mf_FULL:
      if (k > 1) { mf0 = mfinittonew(mf); break; }
      newd = lg(MF_get_S(mf))-1 - mfolddim(N, k, CHI);
      break;
    default: pari_err_TYPE("mfsplit [space does not contain newspace]", mf);
      return NULL; 
  }
  if (newd < 0) newd = mf0? MF_get_dim(mf0): 0;
  *pnewd = newd;
  if (!newd) return mkvec2(cgetg(1, t_MAT), empty);
  NEWT = (k > 1 && MF_get_space(mf0) == mf_NEW);
  todosp = mkvec( split_starting_space(mf0) );
  simplesp = empty;
  FC = mfcharconductor(CHI);
  ord = mfcharorder_canon(CHI);
  if (ord == 1) NF = POLCYC = NULL;
  else
  {
    POLCYC = mfcharpol(CHI);
    NF = nfinit(POLCYC,DEFAULTPREC);
  }
  Tpbigvec = zerovec(NBH);
  u_forprime_init(&iter, 2, ULONG_MAX);
  while (dimsimple < newd && (p = u_forprime_next(&iter)))
  {
    GEN nextsp;
    long ind;
    if (N % (p*p) == 0 && N/p % FC == 0) continue; 
    vecpush(Tpbigvec, NEWT? mfnewmathecke_p(mf0,p): mfheckemat_p(mf0,p));
    nextsp = empty;
    for (ind = 1; ind < lg(todosp); ind++)
    {
      GEN tmp = gel(todosp, ind), fa, P, E, D, Tp, DTp;
      GEN A = gel(tmp, 1);
      GEN X = gel(tmp, 2);
      long lP, i;
      tmp = findbestsplit(NF, Tpbigvec, A, X, dimlim, vz);
      if (!tmp) continue; 
      Tp = gel(tmp, 1);
      fa = gel(tmp, 2);
      P = gel(fa, 1);
      E = gel(fa, 2); lP = lg(P);
      if (DEBUGLEVEL) err_printf("Exponents = %Ps\n", E);
      if (lP == 2)
      {
        GEN P1 = gel(P,1);
        long e1 = itos(gel(E,1)), d1 = degpol(P1);
        if (e1 * d1 == lg(Tp)-1)
        {
          if (e1 > 1) nextsp = vec_append(nextsp, mkvec2(A,X));
          else
          { 
            simplesp = vec_append(simplesp, mkvec3(A,Tp,P1));
            dimsimple += d1;
          }
          continue;
        }
      }
      DTp = Q_remove_denom(Tp, &D);
      for (i = 1; i < lP; i++)
      {
        GEN Ai, Xi, dXi, AAi, v, y, Pi = gel(P,i);
        Ai = RgX_RgM_eval(D? RgX_rescale(Pi,D): Pi, DTp);
        Ai = QabM_ker(Ai, POLCYC, ord);
        if (NF) Ai = nf_primpart(NF, Ai);
        AAi = RgM_mul(A, Ai);
        Xi = QabM_pseudoinv(Ai, POLCYC, ord, &v, &dXi);
        Xi = RgM_Rg_div(Xi, dXi);
        y = gel(v,1);
        if (isint1(gel(E,i)))
        {
          GEN Tpi = RgM_mul(Xi, RgM_mul(rowpermute(Tp,y), Ai));
          simplesp = vec_append(simplesp, mkvec3(AAi, Tpi, Pi));
          dimsimple += degpol(Pi);
        }
        else
        {
          Xi = RgM_mul(Xi, rowpermute(X,y));
          nextsp = vec_append(nextsp, mkvec2(AAi, Xi));
        }
      }
    }
    todosp = nextsp; if (lg(todosp) == 1) break;
  }
  if (DEBUGLEVEL) err_printf("end split, need to clean\n");
  return mfspclean(mf, mf0, NF, ord, sort_by_dim(simplesp), flag);
}
static GEN split_ii ( GEN mf , long dimlim , long flag , long * pnewd ) { forprime_t iter ; GEN CHI = MF_get_CHI ( mf ) , empty = cgetg ( 1 , t_VEC ) , mf0 = mf ; GEN NF , POLCYC , todosp , Tpbigvec , simplesp ; long N = MF_get_N ( mf ) , k = MF_get_k ( mf ) ; long ord , FC , NEWT , dimsimple = 0 , newd = -1 ; const long NBH = 5 , vz = 1 ; ulong p ; switch ( MF_get_space ( mf ) ) { case mf_NEW : break ; case mf_CUSP : case mf_FULL : if ( k > 1 ) { mf0 = mfinittonew ( mf ) ; break ; } newd = lg ( MF_get_S ( mf ) ) - 1 - mfolddim ( N , k , CHI ) ; break ; default : pari_err_TYPE ( " " , mf ) ; return NULL ; } if ( newd < 0 ) newd = mf0 ? MF_get_dim ( mf0 ) : 0 ; * pnewd = newd ; if ( ! newd ) return mkvec2 ( cgetg ( 1 , t_MAT ) , empty ) ; NEWT = ( k > 1 && MF_get_space ( mf0 ) == mf_NEW ) ; todosp = mkvec ( split_starting_space ( mf0 ) ) ; simplesp = empty ; FC = mfcharconductor ( CHI ) ; ord = mfcharorder_canon ( CHI ) ; if ( ord == 1 ) NF = POLCYC = NULL ; else { POLCYC = mfcharpol ( CHI ) ; NF = nfinit ( POLCYC , DEFAULTPREC ) ; } Tpbigvec = zerovec ( NBH ) ; u_forprime_init ( & iter , 2 , ULONG_MAX ) ; while ( dimsimple < newd && ( p = u_forprime_next ( & iter ) ) ) { GEN nextsp ; long ind ; if ( N % ( p * p ) == 0 && N / p % FC == 0 ) continue ; vecpush ( Tpbigvec , NEWT ? mfnewmathecke_p ( mf0 , p ) : mfheckemat_p ( mf0 , p ) ) ; nextsp = empty ; for ( ind = 1 ; ind < lg ( todosp ) ; ind ++ ) { GEN tmp = gel ( todosp , ind ) , fa , P , E , D , Tp , DTp ; GEN A = gel ( tmp , 1 ) ; GEN X = gel ( tmp , 2 ) ; long lP , i ; tmp = findbestsplit ( NF , Tpbigvec , A , X , dimlim , vz ) ; if ( ! tmp ) continue ; Tp = gel ( tmp , 1 ) ; fa = gel ( tmp , 2 ) ; P = gel ( fa , 1 ) ; E = gel ( fa , 2 ) ; lP = lg ( P ) ; if ( DEBUGLEVEL ) err_printf ( " \n " , E ) ; if ( lP == 2 ) { GEN P1 = gel ( P , 1 ) ; long e1 = itos ( gel ( E , 1 ) ) , d1 = degpol ( P1 ) ; if ( e1 * d1 == lg ( Tp ) - 1 ) { if ( e1 > 1 ) nextsp = vec_append ( nextsp , mkvec2 ( A , X ) ) ; else { simplesp = vec_append ( simplesp , mkvec3 ( A , Tp , P1 ) ) ; dimsimple += d1 ; } continue ; } } DTp = Q_remove_denom ( Tp , & D ) ; for ( i = 1 ; i < lP ; i ++ ) { GEN Ai , Xi , dXi , AAi , v , y , Pi = gel ( P , i ) ; Ai = RgX_RgM_eval ( D ? RgX_rescale ( Pi , D ) : Pi , DTp ) ; Ai = QabM_ker ( Ai , POLCYC , ord ) ; if ( NF ) Ai = nf_primpart ( NF , Ai ) ; AAi = RgM_mul ( A , Ai ) ; Xi = QabM_pseudoinv ( Ai , POLCYC , ord , & v , & dXi ) ; Xi = RgM_Rg_div ( Xi , dXi ) ; y = gel ( v , 1 ) ; if ( isint1 ( gel ( E , i ) ) ) { GEN Tpi = RgM_mul ( Xi , RgM_mul ( rowpermute ( Tp , y ) , Ai ) ) ; simplesp = vec_append ( simplesp , mkvec3 ( AAi , Tpi , Pi ) ) ; dimsimple += degpol ( Pi ) ; } else { Xi = RgM_mul ( Xi , rowpermute ( X , y ) ) ; nextsp = vec_append ( nextsp , mkvec2 ( AAi , Xi ) ) ; } } } todosp = nextsp ; if ( lg ( todosp ) == 1 ) break ; } if ( DEBUGLEVEL ) err_printf ( " \n " ) ; return mfspclean ( mf , mf0 , NF , ord , sort_by_dim ( simplesp ) , flag ) ; }
If dimlim > 0, keep only the dimension <= dimlim eigenspaces.
If dimlim > 0 keep only the dimension < = dimlim eigenspaces .
/*LCOV_EXCL_LINE*/ /* T_p = 0 in this case */ /* nothing there */ /* lP > 1 */ /* simple module */ /* Found splitting */ /* gives section, works on nonsquare matrices */