bridge/class-master/depend/pari/src/basemath/nffactor.c
/* P,Q in Z[X,Y], T in Z[Y] irreducible. compute GCD in Q[Y]/(T)[X].
 *
 * M. Encarnacion "On a modular Algorithm for computing GCDs of polynomials
 * over number fields" (ISSAC'94).
 *
 * We procede as follows
 *  1:compute the gcd modulo primes discarding bad primes as they are detected.
 *  2:reconstruct the result via FpM_ratlift, stoping as soon as we get weird
 *    denominators.
 *  3:if FpM_ratlift succeeds, try the full division.
 * Suppose accuracy is insufficient to get the result right: FpM_ratlift will
 * rarely succeed, and even if it does the polynomial we get has sensible
 * coefficients, so the full division will not be too costly.
 *
 * If not NULL, den must be a multiple of the denominator of the gcd,
 * for example the discriminant of T.
 *
 * NOTE: if T is not irreducible, nfgcd may loop forever, esp. if gcd | T */
P,Q in Z[X,Y], T in Z[Y] irreducible. compute GCD in Q[Y]/(T)[X].
M. Encarnacion "On a modular Algorithm for computing GCDs of polynomials
over number fields" (ISSAC'94).

We procede as follows
1:compute the gcd modulo primes discarding bad primes as they are detected.
2:reconstruct the result via FpM_ratlift, stoping as soon as we get weird
denominators.
3:if FpM_ratlift succeeds, try the full division.
Suppose accuracy is insufficient to get the result right: FpM_ratlift will
rarely succeed, and even if it does the polynomial we get has sensible
coefficients, so the full division will not be too costly.

If not NULL, den must be a multiple of the denominator of the gcd,
for example the discriminant of T.

if T is not irreducible, nfgcd may loop forever, esp. if gcd | T
P Q in Z [ X Y ] T in Z [ Y ] irreducible . compute GCD in Q [ Y ] / ( T ) [ X ] . M . Encarnacion " On a modular Algorithm for computing GCDs of polynomials over number fields " ( ISSAC ' 94 ) . We procede as follows 1 : compute the gcd modulo primes discarding bad primes as they are detected . 2 : reconstruct the result via FpM_ratlift stoping as soon as we get weird denominators . 3 : if FpM_ratlift succeeds try the full division . Suppose accuracy is insufficient to get the result right : FpM_ratlift will rarely succeed and even if it does the polynomial we get has sensible coefficients so the full division will not be too costly . If not NULL den must be a multiple of the denominator of the gcd for example the discriminant of T . if T is not irreducible nfgcd may loop forever esp . if gcd | T
GEN
nfgcd_all(GEN P, GEN Q, GEN T, GEN den, GEN *Pnew)
{
  pari_sp btop, ltop = avma;
  GEN lP, lQ, M, dsol, R, bo, sol, mod = NULL, lden = NULL;
  long vP = varn(P), vT = varn(T), dT = degpol(T), dM = 0, dR;
  forprime_t S;
  if (!signe(P)) { if (Pnew) *Pnew = pol_0(vT); return gcopy(Q); }
  if (!signe(Q)) { if (Pnew) *Pnew = pol_1(vT);   return gcopy(P); }
  if ((lP = lead_simplify(P)) && (lQ = lead_simplify(Q)))
  {
    if (typ(lP) == t_INT && typ(lQ) == t_INT)
      lden = powiu(gcdii(lP, lQ), dT);
    else if (typ(lP) == t_INT)
      lden = gcdii(powiu(lP, dT), ZX_resultant(lQ, T));
    else if (typ(lQ) == t_INT)
      lden = gcdii(powiu(lQ, dT), ZX_resultant(lP, T));
    else
      lden = gcdii(ZX_resultant(lP, T), ZX_resultant(lQ, T));
    if (is_pm1(lden)) lden = NULL;
    den = mul_denom(den, lden);
  }
  init_modular_small(&S);
  btop = avma;
  for(;;)
  {
    ulong p = u_forprime_next(&S);
    GEN Tp;
    if (!p) pari_err_OVERFLOW("nfgcd [ran out of primes]");
    if (lden && !umodiu(lden, p)) continue;
    Tp = ZX_to_Flx(T,p);
    if (!Flx_is_squarefree(Tp, p)) continue;
    if ((R = FlxqX_safegcd(ZXX_to_FlxX(P,p,vT),
                           ZXX_to_FlxX(Q,p,vT),
                           Tp, p)) == NULL) continue;
    dR = degpol(R);
    if (dR == 0) { avma = ltop; if (Pnew) *Pnew = P; return pol_1(vP); }
    if (mod && dR > dM) continue; 
    R = FlxX_to_Flm(R, dT);
    if (!mod || dR < dM) { M = ZM_init_CRT(R, p); mod = utoipos(p); dM = dR; continue; }
    (void)ZM_incremental_CRT(&M,R, &mod,p);
    if (gc_needed(btop, 1))
    {
      if (DEBUGMEM>1) pari_warn(warnmem,"nfgcd");
      gerepileall(btop, 2, &M, &mod);
    }
    bo = sqrti(shifti(mod, -1));
    if ((sol = FpM_ratlift(M, mod, bo, bo, den)) == NULL) continue;
    sol = RgM_to_RgXX(sol,vP,vT);
    dsol = Q_primpart(sol);
    if (!ZXQX_dvd(Q, dsol, T)) continue;
    if (Pnew)
    {
      *Pnew = RgXQX_pseudodivrem(P, dsol, T, &R);
      if (signe(R)) continue;
    }
    else
    {
      if (!ZXQX_dvd(P, dsol, T)) continue;
    }
    gerepileall(ltop, Pnew? 2: 1, &dsol, Pnew);
    return dsol; 
  }
}
GEN nfgcd_all ( GEN P , GEN Q , GEN T , GEN den , GEN * Pnew ) { pari_sp btop , ltop = avma ; GEN lP , lQ , M , dsol , R , bo , sol , mod = NULL , lden = NULL ; long vP = varn ( P ) , vT = varn ( T ) , dT = degpol ( T ) , dM = 0 , dR ; forprime_t S ; if ( ! signe ( P ) ) { if ( Pnew ) * Pnew = pol_0 ( vT ) ; return gcopy ( Q ) ; } if ( ! signe ( Q ) ) { if ( Pnew ) * Pnew = pol_1 ( vT ) ; return gcopy ( P ) ; } if ( ( lP = lead_simplify ( P ) ) && ( lQ = lead_simplify ( Q ) ) ) { if ( typ ( lP ) == t_INT && typ ( lQ ) == t_INT ) lden = powiu ( gcdii ( lP , lQ ) , dT ) ; else if ( typ ( lP ) == t_INT ) lden = gcdii ( powiu ( lP , dT ) , ZX_resultant ( lQ , T ) ) ; else if ( typ ( lQ ) == t_INT ) lden = gcdii ( powiu ( lQ , dT ) , ZX_resultant ( lP , T ) ) ; else lden = gcdii ( ZX_resultant ( lP , T ) , ZX_resultant ( lQ , T ) ) ; if ( is_pm1 ( lden ) ) lden = NULL ; den = mul_denom ( den , lden ) ; } init_modular_small ( & S ) ; btop = avma ; for ( ; ; ) { ulong p = u_forprime_next ( & S ) ; GEN Tp ; if ( ! p ) pari_err_OVERFLOW ( " " ) ; if ( lden && ! umodiu ( lden , p ) ) continue ; Tp = ZX_to_Flx ( T , p ) ; if ( ! Flx_is_squarefree ( Tp , p ) ) continue ; if ( ( R = FlxqX_safegcd ( ZXX_to_FlxX ( P , p , vT ) , ZXX_to_FlxX ( Q , p , vT ) , Tp , p ) ) == NULL ) continue ; dR = degpol ( R ) ; if ( dR == 0 ) { avma = ltop ; if ( Pnew ) * Pnew = P ; return pol_1 ( vP ) ; } if ( mod && dR > dM ) continue ; R = FlxX_to_Flm ( R , dT ) ; if ( ! mod || dR < dM ) { M = ZM_init_CRT ( R , p ) ; mod = utoipos ( p ) ; dM = dR ; continue ; } ( void ) ZM_incremental_CRT ( & M , R , & mod , p ) ; if ( gc_needed ( btop , 1 ) ) { if ( DEBUGMEM > 1 ) pari_warn ( warnmem , " " ) ; gerepileall ( btop , 2 , & M , & mod ) ; } bo = sqrti ( shifti ( mod , -1 ) ) ; if ( ( sol = FpM_ratlift ( M , mod , bo , bo , den ) ) == NULL ) continue ; sol = RgM_to_RgXX ( sol , vP , vT ) ; dsol = Q_primpart ( sol ) ; if ( ! ZXQX_dvd ( Q , dsol , T ) ) continue ; if ( Pnew ) { * Pnew = RgXQX_pseudodivrem ( P , dsol , T , & R ) ; if ( signe ( R ) ) continue ; } else { if ( ! ZXQX_dvd ( P , dsol , T ) ) continue ; } gerepileall ( ltop , Pnew ? 2 : 1 , & dsol , Pnew ) ; return dsol ; } }
P,Q in Z[X,Y], T in Z[Y] irreducible.
P Q in Z [ X Y ] T in Z [ Y ] irreducible .
/* Compute denominators */ /*Discard primes dividing disc(T) or lc(PQ) */ /*Discard primes when modular gcd does not exist*/ /* p divides Res(P/gcd, Q/gcd). Discard. */ /* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */ /* I suspect it must be better to take amax > bmax*/ /* both remainders are 0 */