My Project
Functions | Variables
cfModGcd.cc File Reference

This file implements the GCD of two polynomials over $ F_{p} $ , $ F_{p}(\alpha ) $, GF or Z based on Alg. More...

#include "config.h"
#include "cf_assert.h"
#include "debug.h"
#include "timing.h"
#include "canonicalform.h"
#include "cfGcdUtil.h"
#include "cf_map.h"
#include "cf_util.h"
#include "cf_irred.h"
#include "templates/ftmpl_functions.h"
#include "cf_random.h"
#include "cf_reval.h"
#include "facHensel.h"
#include "cf_iter.h"
#include "cfNewtonPolygon.h"
#include "cf_algorithm.h"
#include "cf_primes.h"
#include "cf_map_ext.h"
#include "NTLconvert.h"
#include "FLINTconvert.h"
#include "cfModGcd.h"

Go to the source code of this file.

Functions

 TIMING_DEFINE_PRINT (gcd_recursion) TIMING_DEFINE_PRINT(newton_interpolation) TIMING_DEFINE_PRINT(termination_test) TIMING_DEFINE_PRINT(ez_p_compress) TIMING_DEFINE_PRINT(ez_p_hensel_lift) TIMING_DEFINE_PRINT(ez_p_eval) TIMING_DEFINE_PRINT(ez_p_content) bool terminationTest(const CanonicalForm &F
 
 if (LCCand *abs(LC(coF))==abs(LC(F)))
 
int myCompress (const CanonicalForm &F, const CanonicalForm &G, CFMap &M, CFMap &N, bool topLevel)
 compressing two polynomials F and G, M is used for compressing, N to reverse the compression More...
 
static CanonicalForm uni_content (const CanonicalForm &F)
 compute the content of F, where F is considered as an element of $ R[x_{1}][x_{2},\ldots ,x_{n}] $ More...
 
CanonicalForm uni_content (const CanonicalForm &F, const Variable &x)
 
CanonicalForm extractContents (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &contentF, CanonicalForm &contentG, CanonicalForm &ppF, CanonicalForm &ppG, const int d)
 
static CanonicalForm uni_lcoeff (const CanonicalForm &F)
 compute the leading coefficient of F, where F is considered as an element of $ R[x_{1}][x_{2},\ldots ,x_{n}] $, order on $ Mon (x_{2},\ldots ,x_{n}) $ is dp. More...
 
static CanonicalForm newtonInterp (const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x)
 Newton interpolation - Incremental algorithm. Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n, computes the interpolation polynomial assuming that the polynomials in u are the results of evaluating the variabe x of the unknown polynomial at the alpha values. This incremental version receives only the values of alpha_n and u_n and the previous interpolation polynomial for points 1 <= i <= n-1, and computes the polynomial interpolating in all the points. newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1}) More...
 
static CanonicalForm randomElement (const CanonicalForm &F, const Variable &alpha, CFList &list, bool &fail)
 compute a random element a of $ F_{p}(\alpha ) $ , s.t. F(a) $ \neq 0 $ , F is a univariate polynomial, returns fail if there are no field elements left which have not been used before More...
 
static Variable chooseExtension (const Variable &alpha)
 
CanonicalForm modGCDFq (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, Variable &alpha, CFList &l, bool &topLevel)
 GCD of F and G over $ F_{p}(\alpha ) $ , l and topLevel are only used internally, output is monic based on Alg. 7.2. as described in "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn. More...
 
CanonicalForm modGCDFq (const CanonicalForm &F, const CanonicalForm &G, Variable &alpha, CFList &l, bool &topLevel)
 
static CanonicalForm GFRandomElement (const CanonicalForm &F, CFList &list, bool &fail)
 compute a random element a of GF, s.t. F(a) $ \neq 0 $ , F is a univariate polynomial, returns fail if there are no field elements left which have not been used before More...
 
CanonicalForm modGCDGF (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, CFList &l, bool &topLevel)
 GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn Usually this algorithm will be faster than modGCDFq since GF has faster field arithmetics, however it might fail if the input is large since the size of the base field is bounded by 2^16, output is monic. More...
 
CanonicalForm modGCDGF (const CanonicalForm &F, const CanonicalForm &G, CFList &l, bool &topLevel)
 
static CanonicalForm FpRandomElement (const CanonicalForm &F, CFList &list, bool &fail)
 
CanonicalForm modGCDFp (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l)
 
CanonicalForm modGCDFp (const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
 
CFArray solveVandermonde (const CFArray &M, const CFArray &A)
 
CFArray solveGeneralVandermonde (const CFArray &M, const CFArray &A)
 
CFArray readOffSolution (const CFMatrix &M, const long rk)
 M in row echolon form, rk rank of M. More...
 
CFArray readOffSolution (const CFMatrix &M, const CFArray &L, const CFArray &partialSol)
 
long gaussianElimFp (CFMatrix &M, CFArray &L)
 
long gaussianElimFq (CFMatrix &M, CFArray &L, const Variable &alpha)
 
CFArray solveSystemFp (const CFMatrix &M, const CFArray &L)
 
CFArray solveSystemFq (const CFMatrix &M, const CFArray &L, const Variable &alpha)
 
CFArray getMonoms (const CanonicalForm &F)
 extract monomials of F, parts in algebraic variable are considered coefficients More...
 
CFArray evaluateMonom (const CanonicalForm &F, const CFList &evalPoints)
 
CFArray evaluate (const CFArray &A, const CFList &evalPoints)
 
CFList evaluationPoints (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Feval, CanonicalForm &Geval, const CanonicalForm &LCF, const bool &GF, const Variable &alpha, bool &fail, CFList &list)
 
void mult (CFList &L1, const CFList &L2)
 multiply two lists componentwise More...
 
void eval (const CanonicalForm &A, const CanonicalForm &B, CanonicalForm &Aeval, CanonicalForm &Beval, const CFList &L)
 
CanonicalForm monicSparseInterpol (const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
 
CanonicalForm nonMonicSparseInterpol (const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
 
CanonicalForm sparseGCDFq (const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
 
CanonicalForm sparseGCDFp (const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
 
 TIMING_DEFINE_PRINT (modZ_termination) TIMING_DEFINE_PRINT(modZ_recursion) CanonicalForm modGCDZ(const CanonicalForm &FF
 modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer Algebra", Algorithm 7.1 More...
 
 for (i=tmax(f.level(), g.level());i > 0;i--)
 
 if (i==0) return gcdcfcg
 
 for (;i > 0;i--)
 
 while (true)
 

Variables

const CanonicalFormG
 
const CanonicalForm const CanonicalFormcoF
 
const CanonicalForm const CanonicalForm const CanonicalFormcoG
 
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalFormcand
 
return false
 
const CanonicalFormGG
 
int p
 
int i = cf_getNumBigPrimes() - 1
 
int dp_deg
 
int d_deg =-1
 
CanonicalForm cd (bCommonDen(FF)) = bCommonDen( GG )
 
 f =cd*FF
 
Variable x = Variable (1)
 
CanonicalForm cf = icontent (f)
 
CanonicalForm cg = icontent (g)
 
 g =cd*GG
 
CanonicalForm Dn
 
CanonicalForm test = 0
 
CanonicalForm lcf = f.lc()
 
CanonicalForm lcg = g.lc()
 
 cl = gcd (f.lc(),g.lc())
 
CanonicalForm gcdcfcg = gcd (cf, cg)
 
CanonicalForm fp
 
CanonicalForm gp
 
CanonicalForm b = 1
 
int minCommonDeg = 0
 
bool equal = false
 
CanonicalForm cof
 
CanonicalForm cog
 
CanonicalForm cofp
 
CanonicalForm cogp
 
CanonicalForm newCof
 
CanonicalForm newCog
 
CanonicalForm cofn
 
CanonicalForm cogn
 
CanonicalForm cDn
 
int maxNumVars = tmax (getNumVars (f), getNumVars (g))
 

Detailed Description

This file implements the GCD of two polynomials over $ F_{p} $ , $ F_{p}(\alpha ) $, GF or Z based on Alg.

7.1. and 7.2. as described in "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn via modular computations. And sparse modular variants as described in Zippel "Effective Polynomial Computation", deKleine, Monagan, Wittkopf "Algorithms for the non-monic case of the sparse modular GCD algorithm" and Javadi "A new solution to the normalization problem"

Author
Martin Lee
Date
22.10.2009
Copyright:
(c) by The SINGULAR Team, see LICENSE file

Definition in file cfModGcd.cc.

Function Documentation

◆ chooseExtension()

static Variable chooseExtension ( const Variable alpha)
inlinestatic

Definition at line 420 of file cfModGcd.cc.

421 {
422  int i, m;
423  // extension of F_p needed
424  if (alpha.level() == 1)
425  {
426  i= 1;
427  m= 2;
428  } //extension of F_p(alpha)
429  if (alpha.level() != 1)
430  {
431  i= 4;
432  m= degree (getMipo (alpha));
433  }
434  #ifdef HAVE_FLINT
435  nmod_poly_t Irredpoly;
436  nmod_poly_init(Irredpoly,getCharacteristic());
437  nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, i*m+1);
438  CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
439  nmod_poly_clear(Irredpoly);
440  #else
442  {
445  }
446  zz_pX NTLIrredpoly;
447  BuildIrred (NTLIrredpoly, i*m);
448  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
449  #endif
450  return rootOf (newMipo);
451 }
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
Definition: NTLconvert.cc:255
VAR long fac_NTL_char
Definition: NTLconvert.cc:46
int degree(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfModGcd.cc:4078
GLOBAL_VAR flint_rand_t FLINTrandom
Definition: cf_random.cc:25
factory's main class
Definition: canonicalform.h:86
factory's class for variables
Definition: factory.h:127
int level() const
Definition: factory.h:143
Variable alpha
Definition: facAbsBiFact.cc:51
nmod_poly_init(FLINTmipo, getCharacteristic())
nmod_poly_clear(FLINTmipo)
Variable FACTORY_PUBLIC rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables
Definition: variable.cc:162
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
void init()
Definition: lintree.cc:864

◆ eval()

void eval ( const CanonicalForm A,
const CanonicalForm B,
CanonicalForm Aeval,
CanonicalForm Beval,
const CFList L 
)

Definition at line 2185 of file cfModGcd.cc.

2187 {
2188  Aeval= A;
2189  Beval= B;
2190  int j= 1;
2191  for (CFListIterator i= L; i.hasItem(); i++, j++)
2192  {
2193  Aeval= Aeval (i.getItem(), j);
2194  Beval= Beval (i.getItem(), j);
2195  }
2196 }
b *CanonicalForm B
Definition: facBivar.cc:52
CFList *& Aeval
<[in] poly
Definition: facFactorize.h:31
int j
Definition: facHensel.cc:110
#define A
Definition: sirandom.c:24

◆ evaluate()

CFArray evaluate ( const CFArray A,
const CFList evalPoints 
)

Definition at line 2031 of file cfModGcd.cc.

2032 {
2033  CFArray result= A.size();
2034  CanonicalForm tmp;
2035  int k;
2036  for (int i= 0; i < A.size(); i++)
2037  {
2038  tmp= A[i];
2039  k= 1;
2040  for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2041  tmp= tmp (j.getItem(), k);
2042  result[i]= tmp;
2043  }
2044  return result;
2045 }
int k
Definition: cfEzgcd.cc:99
return result
Definition: facAbsBiFact.cc:75
CFList evalPoints(const CanonicalForm &F, CFList &eval, const Variable &alpha, CFList &list, const bool &GF, bool &fail)
evaluation point search for multivariate factorization, looks for a (F.level() - 1)-tuple such that t...

◆ evaluateMonom()

CFArray evaluateMonom ( const CanonicalForm F,
const CFList evalPoints 
)

Definition at line 1992 of file cfModGcd.cc.

1993 {
1994  if (F.inCoeffDomain())
1995  {
1996  CFArray result= CFArray (1);
1997  result [0]= F;
1998  return result;
1999  }
2000  if (F.isUnivariate())
2001  {
2002  ASSERT (evalPoints.length() == 1,
2003  "expected an eval point with only one component");
2004  CFArray result= CFArray (size(F));
2005  int j= 0;
2007  for (CFIterator i= F; i.hasTerms(); i++, j++)
2008  result[j]= power (evalPoint, i.exp());
2009  return result;
2010  }
2011  int numMon= size (F);
2012  CFArray result= CFArray (numMon);
2013  int j= 0;
2016  buf.removeLast();
2017  CFArray recResult;
2018  CanonicalForm powEvalPoint;
2019  for (CFIterator i= F; i.hasTerms(); i++)
2020  {
2021  powEvalPoint= power (evalPoint, i.exp());
2022  recResult= evaluateMonom (i.coeff(), buf);
2023  for (int k= 0; k < recResult.size(); k++)
2024  result[j+k]= powEvalPoint*recResult[k];
2025  j += recResult.size();
2026  }
2027  return result;
2028 }
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
Array< CanonicalForm > CFArray
CFArray evaluateMonom(const CanonicalForm &F, const CFList &evalPoints)
Definition: cfModGcd.cc:1992
static void evalPoint(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &FEval, CanonicalForm &GEval, CFGenerator &evalPoint)
#define ASSERT(expression, message)
Definition: cf_assert.h:99
int size() const
Definition: ftmpl_array.cc:92
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
bool inCoeffDomain() const
bool isUnivariate() const
int length() const
Definition: ftmpl_list.cc:273
T getLast() const
Definition: ftmpl_list.cc:309
int status int void * buf
Definition: si_signals.h:59

◆ evaluationPoints()

CFList evaluationPoints ( const CanonicalForm F,
const CanonicalForm G,
CanonicalForm Feval,
CanonicalForm Geval,
const CanonicalForm LCF,
const bool &  GF,
const Variable alpha,
bool &  fail,
CFList list 
)

Definition at line 2048 of file cfModGcd.cc.

2053 {
2054  int k= tmax (F.level(), G.level()) - 1;
2055  Variable x= Variable (1);
2056  CFList result;
2057  FFRandom genFF;
2058  GFRandom genGF;
2059  int p= getCharacteristic ();
2060  double bound;
2061  if (alpha != Variable (1))
2062  {
2063  bound= pow ((double) p, (double) degree (getMipo(alpha)));
2064  bound= pow (bound, (double) k);
2065  }
2066  else if (GF)
2067  {
2068  bound= pow ((double) p, (double) getGFDegree());
2069  bound= pow ((double) bound, (double) k);
2070  }
2071  else
2072  bound= pow ((double) p, (double) k);
2073 
2074  CanonicalForm random;
2075  int j;
2076  bool zeroOneOccured= false;
2077  bool allEqual= false;
2079  do
2080  {
2081  random= 0;
2082  // possible overflow if list.length() does not fit into a int
2083  if (list.length() >= bound)
2084  {
2085  fail= true;
2086  break;
2087  }
2088  for (int i= 0; i < k; i++)
2089  {
2090  if (GF)
2091  {
2092  result.append (genGF.generate());
2093  random += result.getLast()*power (x, i);
2094  }
2095  else if (alpha.level() != 1)
2096  {
2097  AlgExtRandomF genAlgExt (alpha);
2098  result.append (genAlgExt.generate());
2099  random += result.getLast()*power (x, i);
2100  }
2101  else
2102  {
2103  result.append (genFF.generate());
2104  random += result.getLast()*power (x, i);
2105  }
2106  if (result.getLast().isOne() || result.getLast().isZero())
2107  zeroOneOccured= true;
2108  }
2109  if (find (list, random))
2110  {
2111  zeroOneOccured= false;
2112  allEqual= false;
2113  result= CFList();
2114  continue;
2115  }
2116  if (zeroOneOccured)
2117  {
2118  list.append (random);
2119  zeroOneOccured= false;
2120  allEqual= false;
2121  result= CFList();
2122  continue;
2123  }
2124  // no zero at this point
2125  if (k > 1)
2126  {
2127  allEqual= true;
2128  CFIterator iter= random;
2129  buf= iter.coeff();
2130  iter++;
2131  for (; iter.hasTerms(); iter++)
2132  if (buf != iter.coeff())
2133  allEqual= false;
2134  }
2135  if (allEqual)
2136  {
2137  list.append (random);
2138  allEqual= false;
2139  zeroOneOccured= false;
2140  result= CFList();
2141  continue;
2142  }
2143 
2144  Feval= F;
2145  Geval= G;
2146  CanonicalForm LCeval= LCF;
2147  j= 1;
2148  for (CFListIterator i= result; i.hasItem(); i++, j++)
2149  {
2150  Feval= Feval (i.getItem(), j);
2151  Geval= Geval (i.getItem(), j);
2152  LCeval= LCeval (i.getItem(), j);
2153  }
2154 
2155  if (LCeval.isZero())
2156  {
2157  if (!find (list, random))
2158  list.append (random);
2159  zeroOneOccured= false;
2160  allEqual= false;
2161  result= CFList();
2162  continue;
2163  }
2164 
2165  if (list.length() >= bound)
2166  {
2167  fail= true;
2168  break;
2169  }
2170  } while (find (list, random));
2171 
2172  return result;
2173 }
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
int getGFDegree()
Definition: cf_char.cc:75
List< CanonicalForm > CFList
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
const CanonicalForm & G
Definition: cfModGcd.cc:66
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
generate random elements in F_p(alpha)
Definition: cf_random.h:70
CF_NO_INLINE bool isZero() const
int level() const
level() returns the level of CO.
generate random elements in F_p
Definition: cf_random.h:44
CanonicalForm generate() const
Definition: cf_random.cc:68
generate random elements in GF
Definition: cf_random.h:32
CanonicalForm generate() const
Definition: cf_random.cc:78
void append(const T &)
Definition: ftmpl_list.cc:256
CFFListIterator iter
Definition: facAbsBiFact.cc:53
CanonicalForm Feval
Definition: facAbsFact.cc:60
CanonicalForm LCF
Definition: facFactorize.cc:52
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
template bool find(const List< CanonicalForm > &, const CanonicalForm &)

◆ extractContents()

CanonicalForm extractContents ( const CanonicalForm F,
const CanonicalForm G,
CanonicalForm contentF,
CanonicalForm contentG,
CanonicalForm ppF,
CanonicalForm ppG,
const int  d 
)

Definition at line 311 of file cfModGcd.cc.

314 {
315  CanonicalForm uniContentF, uniContentG, gcdcFcG;
316  contentF= 1;
317  contentG= 1;
318  ppF= F;
319  ppG= G;
321  for (int i= 1; i <= d; i++)
322  {
323  uniContentF= uni_content (F, Variable (i));
324  uniContentG= uni_content (G, Variable (i));
325  gcdcFcG= gcd (uniContentF, uniContentG);
326  contentF *= uniContentF;
327  contentG *= uniContentG;
328  ppF /= uniContentF;
329  ppG /= uniContentG;
330  result *= gcdcFcG;
331  }
332  return result;
333 }
static CanonicalForm uni_content(const CanonicalForm &F)
compute the content of F, where F is considered as an element of
Definition: cfModGcd.cc:281
int gcd(int a, int b)
Definition: walkSupport.cc:836

◆ for() [1/2]

for ( i,
0;i--   
)

Definition at line 4117 of file cfModGcd.cc.

4118  {
4119  if (degree (f, i) <= 0 || degree (g, i) <= 0)
4120  continue;
4121  else
4123  }
f
Definition: cfModGcd.cc:4081
g
Definition: cfModGcd.cc:4090
int minCommonDeg
Definition: cfModGcd.cc:4104
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)

◆ for() [2/2]

for ( i  = tmax (f.level(), g.level()); i,
0;i--   
)

Definition at line 4105 of file cfModGcd.cc.

4106  {
4107  if (degree (f, i) <= 0 || degree (g, i) <= 0)
4108  continue;
4109  else
4110  {
4111  minCommonDeg= tmin (degree (g, i), degree (f, i));
4112  break;
4113  }
4114  }

◆ FpRandomElement()

static CanonicalForm FpRandomElement ( const CanonicalForm F,
CFList list,
bool &  fail 
)
inlinestatic

Definition at line 1171 of file cfModGcd.cc.

1172 {
1173  fail= false;
1174  Variable x= F.mvar();
1175  FFRandom genFF;
1176  CanonicalForm random;
1177  int p= getCharacteristic();
1178  int bound= p;
1179  do
1180  {
1181  if (list.length() == bound)
1182  {
1183  fail= true;
1184  break;
1185  }
1186  if (list.length() < 1)
1187  random= 0;
1188  else
1189  {
1190  random= genFF.generate();
1191  while (find (list, random))
1192  random= genFF.generate();
1193  }
1194  if (F (random, x) == 0)
1195  {
1196  list.append (random);
1197  continue;
1198  }
1199  } while (find (list, random));
1200  return random;
1201 }
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.

◆ gaussianElimFp()

long gaussianElimFp ( CFMatrix M,
CFArray L 
)

Definition at line 1739 of file cfModGcd.cc.

1740 {
1741  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1742  CFMatrix *N;
1743  N= new CFMatrix (M.rows(), M.columns() + 1);
1744 
1745  for (int i= 1; i <= M.rows(); i++)
1746  for (int j= 1; j <= M.columns(); j++)
1747  (*N) (i, j)= M (i, j);
1748 
1749  int j= 1;
1750  for (int i= 0; i < L.size(); i++, j++)
1751  (*N) (j, M.columns() + 1)= L[i];
1752 #ifdef HAVE_FLINT
1753  nmod_mat_t FLINTN;
1754  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1755  long rk= nmod_mat_rref (FLINTN);
1756 
1757  delete N;
1758  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1759  nmod_mat_clear (FLINTN);
1760 #else
1761  int p= getCharacteristic ();
1762  if (fac_NTL_char != p)
1763  {
1764  fac_NTL_char= p;
1765  zz_p::init (p);
1766  }
1767  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1768  delete N;
1769  long rk= gauss (*NTLN);
1770 
1772  delete NTLN;
1773 #endif
1774 
1775  L= CFArray (M.rows());
1776  for (int i= 0; i < M.rows(); i++)
1777  L[i]= (*N) (i + 1, M.columns() + 1);
1778  M= (*N) (1, M.rows(), 1, M.columns());
1779  delete N;
1780  return rk;
1781 }
void convertFacCFMatrix2nmod_mat_t(nmod_mat_t M, const CFMatrix &m)
conversion of a factory matrix over Z/p to a nmod_mat_t
CFMatrix * convertNmod_mat_t2FacCFMatrix(const nmod_mat_t m)
conversion of a FLINT matrix over Z/p to a factory matrix
CFMatrix * convertNTLmat_zz_p2FacCFMatrix(const mat_zz_p &m)
Definition: NTLconvert.cc:1183
mat_zz_p * convertFacCFMatrix2NTLmat_zz_p(const CFMatrix &m)
Definition: NTLconvert.cc:1167
Matrix< CanonicalForm > CFMatrix
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
#define M
Definition: sirandom.c:25

◆ gaussianElimFq()

long gaussianElimFq ( CFMatrix M,
CFArray L,
const Variable alpha 
)

Definition at line 1784 of file cfModGcd.cc.

1785 {
1786  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1787  CFMatrix *N;
1788  N= new CFMatrix (M.rows(), M.columns() + 1);
1789 
1790  for (int i= 1; i <= M.rows(); i++)
1791  for (int j= 1; j <= M.columns(); j++)
1792  (*N) (i, j)= M (i, j);
1793 
1794  int j= 1;
1795  for (int i= 0; i < L.size(); i++, j++)
1796  (*N) (j, M.columns() + 1)= L[i];
1797  #ifdef HAVE_FLINT
1798  // convert mipo
1799  nmod_poly_t mipo1;
1801  fq_nmod_ctx_t ctx;
1802  fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
1803  nmod_poly_clear(mipo1);
1804  // convert matrix
1805  fq_nmod_mat_t FLINTN;
1806  convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
1807  // rank
1808  long rk= fq_nmod_mat_rref (FLINTN,ctx);
1809  // clean up
1810  fq_nmod_mat_clear (FLINTN,ctx);
1811  fq_nmod_ctx_clear(ctx);
1812  #elif defined(HAVE_NTL)
1813  int p= getCharacteristic ();
1814  if (fac_NTL_char != p)
1815  {
1816  fac_NTL_char= p;
1817  zz_p::init (p);
1818  }
1819  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1820  zz_pE::init (NTLMipo);
1821  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1822  long rk= gauss (*NTLN);
1824  delete NTLN;
1825  #else
1826  factoryError("NTL/FLINT missing: gaussianElimFq");
1827  #endif
1828  delete N;
1829 
1830  M= (*N) (1, M.rows(), 1, M.columns());
1831  L= CFArray (M.rows());
1832  for (int i= 0; i < M.rows(); i++)
1833  L[i]= (*N) (i + 1, M.columns() + 1);
1834 
1835  delete N;
1836  return rk;
1837 }
void convertFacCFMatrix2Fq_nmod_mat_t(fq_nmod_mat_t M, const fq_nmod_ctx_t fq_con, const CFMatrix &m)
conversion of a factory matrix over F_q to a fq_nmod_mat_t
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:105
mat_zz_pE * convertFacCFMatrix2NTLmat_zz_pE(const CFMatrix &m)
Definition: NTLconvert.cc:1196
CFMatrix * convertNTLmat_zz_pE2FacCFMatrix(const mat_zz_pE &m, const Variable &alpha)
Definition: NTLconvert.cc:1212
VAR void(* factoryError)(const char *s)
Definition: cf_util.cc:80
fq_nmod_ctx_clear(fq_con)
fq_nmod_ctx_init_modulus(fq_con, FLINTmipo, "Z")
convertFacCF2nmod_poly_t(FLINTmipo, M)

◆ getMonoms()

CFArray getMonoms ( const CanonicalForm F)

extract monomials of F, parts in algebraic variable are considered coefficients

Parameters
[in]Fsome poly

Definition at line 1957 of file cfModGcd.cc.

1958 {
1959  if (F.inCoeffDomain())
1960  {
1961  CFArray result= CFArray (1);
1962  result [0]= 1;
1963  return result;
1964  }
1965  if (F.isUnivariate())
1966  {
1967  CFArray result= CFArray (size(F));
1968  int j= 0;
1969  for (CFIterator i= F; i.hasTerms(); i++, j++)
1970  result[j]= power (F.mvar(), i.exp());
1971  return result;
1972  }
1973  int numMon= size (F);
1974  CFArray result= CFArray (numMon);
1975  int j= 0;
1976  CFArray recResult;
1977  Variable x= F.mvar();
1978  CanonicalForm powX;
1979  for (CFIterator i= F; i.hasTerms(); i++)
1980  {
1981  powX= power (x, i.exp());
1982  recResult= getMonoms (i.coeff());
1983  for (int k= 0; k < recResult.size(); k++)
1984  result[j+k]= powX*recResult[k];
1985  j += recResult.size();
1986  }
1987  return result;
1988 }
CFArray getMonoms(const CanonicalForm &F)
extract monomials of F, parts in algebraic variable are considered coefficients
Definition: cfModGcd.cc:1957

◆ GFRandomElement()

static CanonicalForm GFRandomElement ( const CanonicalForm F,
CFList list,
bool &  fail 
)
inlinestatic

compute a random element a of GF, s.t. F(a) $ \neq 0 $ , F is a univariate polynomial, returns fail if there are no field elements left which have not been used before

Definition at line 819 of file cfModGcd.cc.

820 {
821  fail= false;
822  Variable x= F.mvar();
823  GFRandom genGF;
824  CanonicalForm random;
825  int p= getCharacteristic();
826  int d= getGFDegree();
827  int bound= ipower (p, d);
828  do
829  {
830  if (list.length() == bound)
831  {
832  fail= true;
833  break;
834  }
835  if (list.length() < 1)
836  random= 0;
837  else
838  {
839  random= genGF.generate();
840  while (find (list, random))
841  random= genGF.generate();
842  }
843  if (F (random, x) == 0)
844  {
845  list.append (random);
846  continue;
847  }
848  } while (find (list, random));
849  return random;
850 }
int ipower(int b, int m)
int ipower ( int b, int m )
Definition: cf_util.cc:27

◆ if() [1/2]

if ( i  = =0)

◆ if() [2/2]

if ( LCCand *  absLC(coF) = abs (LC (F)))

Definition at line 71 of file cfModGcd.cc.

72  {
73  if (LCCand*abs (LC (coG)) == abs (LC (G)))
74  {
75  if (abs (cand)*abs (coF) == abs (F))
76  {
77  if (abs (cand)*abs (coG) == abs (G))
78  return true;
79  }
80  return false;
81  }
82  return false;
83  }
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
CanonicalForm LC(const CanonicalForm &f)
const CanonicalForm const CanonicalForm const CanonicalForm & coG
Definition: cfModGcd.cc:67
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition: cfModGcd.cc:69
const CanonicalForm const CanonicalForm & coF
Definition: cfModGcd.cc:67

◆ modGCDFp() [1/2]

CanonicalForm modGCDFp ( const CanonicalForm F,
const CanonicalForm G,
bool &  topLevel,
CFList l 
)

Definition at line 1212 of file cfModGcd.cc.

1214 {
1215  CanonicalForm dummy1, dummy2;
1216  CanonicalForm result= modGCDFp (F, G, dummy1, dummy2, topLevel, l);
1217  return result;
1218 }
int l
Definition: cfEzgcd.cc:100
const CanonicalForm CFMap CFMap bool topLevel
Definition: cfGcdAlgExt.cc:57
CanonicalForm modGCDFp(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:1223

◆ modGCDFp() [2/2]

CanonicalForm modGCDFp ( const CanonicalForm F,
const CanonicalForm G,
CanonicalForm coF,
CanonicalForm coG,
bool &  topLevel,
CFList l 
)

Definition at line 1223 of file cfModGcd.cc.

1226 {
1227  CanonicalForm A= F;
1228  CanonicalForm B= G;
1229  if (F.isZero() && degree(G) > 0)
1230  {
1231  coF= 0;
1232  coG= Lc (G);
1233  return G/Lc(G);
1234  }
1235  else if (G.isZero() && degree (F) > 0)
1236  {
1237  coF= Lc (F);
1238  coG= 0;
1239  return F/Lc(F);
1240  }
1241  else if (F.isZero() && G.isZero())
1242  {
1243  coF= coG= 0;
1244  return F.genOne();
1245  }
1246  if (F.inBaseDomain() || G.inBaseDomain())
1247  {
1248  coF= F;
1249  coG= G;
1250  return F.genOne();
1251  }
1252  if (F.isUnivariate() && fdivides(F, G, coG))
1253  {
1254  coF= Lc (F);
1255  return F/Lc(F);
1256  }
1257  if (G.isUnivariate() && fdivides(G, F, coF))
1258  {
1259  coG= Lc (G);
1260  return G/Lc(G);
1261  }
1262  if (F == G)
1263  {
1264  coF= coG= Lc (F);
1265  return F/Lc(F);
1266  }
1267  CFMap M,N;
1268  int best_level= myCompress (A, B, M, N, topLevel);
1269 
1270  if (best_level == 0)
1271  {
1272  coF= F;
1273  coG= G;
1274  return B.genOne();
1275  }
1276 
1277  A= M(A);
1278  B= M(B);
1279 
1280  Variable x= Variable (1);
1281 
1282  //univariate case
1283  if (A.isUnivariate() && B.isUnivariate())
1284  {
1285  CanonicalForm result= gcd (A, B);
1286  coF= N (A/result);
1287  coG= N (B/result);
1288  return N (result);
1289  }
1290 
1291  CanonicalForm cA, cB; // content of A and B
1292  CanonicalForm ppA, ppB; // primitive part of A and B
1293  CanonicalForm gcdcAcB;
1294 
1295  cA = uni_content (A);
1296  cB = uni_content (B);
1297  gcdcAcB= gcd (cA, cB);
1298  ppA= A/cA;
1299  ppB= B/cB;
1300 
1301  CanonicalForm lcA, lcB; // leading coefficients of A and B
1302  CanonicalForm gcdlcAlcB;
1303  lcA= uni_lcoeff (ppA);
1304  lcB= uni_lcoeff (ppB);
1305 
1306  gcdlcAlcB= gcd (lcA, lcB);
1307 
1308  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
1309  int d0;
1310 
1311  if (d == 0)
1312  {
1313  coF= N (ppA*(cA/gcdcAcB));
1314  coG= N (ppB*(cB/gcdcAcB));
1315  return N(gcdcAcB);
1316  }
1317 
1318  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
1319 
1320  if (d0 < d)
1321  d= d0;
1322 
1323  if (d == 0)
1324  {
1325  coF= N (ppA*(cA/gcdcAcB));
1326  coG= N (ppB*(cB/gcdcAcB));
1327  return N(gcdcAcB);
1328  }
1329 
1330  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
1331  CanonicalForm newtonPoly, coF_random_element, coG_random_element,
1332  coF_m, coG_m, ppCoF, ppCoG;
1333 
1334  newtonPoly= 1;
1335  m= gcdlcAlcB;
1336  H= 0;
1337  coF= 0;
1338  coG= 0;
1339  G_m= 0;
1340  Variable alpha, V_buf, cleanUp;
1341  bool fail= false;
1342  bool inextension= false;
1343  topLevel= false;
1344  CFList source, dest;
1345  int bound1= degree (ppA, 1);
1346  int bound2= degree (ppB, 1);
1347  do
1348  {
1349  if (inextension)
1350  random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
1351  else
1352  random_element= FpRandomElement (m*lcA*lcB, l, fail);
1353 
1354  if (!fail && !inextension)
1355  {
1356  CFList list;
1357  TIMING_START (gcd_recursion);
1358  G_random_element=
1359  modGCDFp (ppA (random_element,x), ppB (random_element,x),
1360  coF_random_element, coG_random_element, topLevel,
1361  list);
1362  TIMING_END_AND_PRINT (gcd_recursion,
1363  "time for recursive call: ");
1364  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1365  }
1366  else if (!fail && inextension)
1367  {
1368  CFList list;
1369  TIMING_START (gcd_recursion);
1370  G_random_element=
1371  modGCDFq (ppA (random_element, x), ppB (random_element, x),
1372  coF_random_element, coG_random_element, V_buf,
1373  list, topLevel);
1374  TIMING_END_AND_PRINT (gcd_recursion,
1375  "time for recursive call: ");
1376  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1377  }
1378  else if (fail && !inextension)
1379  {
1380  source= CFList();
1381  dest= CFList();
1382  CFList list;
1384  int deg= 2;
1385  bool initialized= false;
1386  do
1387  {
1388  mipo= randomIrredpoly (deg, x);
1389  if (initialized)
1390  setMipo (alpha, mipo);
1391  else
1392  alpha= rootOf (mipo);
1393  inextension= true;
1394  initialized= true;
1395  fail= false;
1396  random_element= randomElement (m*lcA*lcB, alpha, l, fail);
1397  deg++;
1398  } while (fail);
1399  list= CFList();
1400  V_buf= alpha;
1401  cleanUp= alpha;
1402  TIMING_START (gcd_recursion);
1403  G_random_element=
1404  modGCDFq (ppA (random_element, x), ppB (random_element, x),
1405  coF_random_element, coG_random_element, alpha,
1406  list, topLevel);
1407  TIMING_END_AND_PRINT (gcd_recursion,
1408  "time for recursive call: ");
1409  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1410  }
1411  else if (fail && inextension)
1412  {
1413  source= CFList();
1414  dest= CFList();
1415 
1416  Variable V_buf3= V_buf;
1417  V_buf= chooseExtension (V_buf);
1418  bool prim_fail= false;
1419  Variable V_buf2;
1420  CanonicalForm prim_elem, im_prim_elem;
1421  prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
1422 
1423  if (V_buf3 != alpha)
1424  {
1425  m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
1426  G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
1427  coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
1428  coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
1429  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
1430  source, dest);
1431  ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
1432  ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
1433  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
1434  dest);
1435  lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
1436  lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
1437  for (CFListIterator i= l; i.hasItem(); i++)
1438  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
1439  source, dest);
1440  }
1441 
1442  ASSERT (!prim_fail, "failure in integer factorizer");
1443  if (prim_fail)
1444  ; //ERROR
1445  else
1446  im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
1447 
1448  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1449  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1450 
1451  for (CFListIterator i= l; i.hasItem(); i++)
1452  i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
1453  im_prim_elem, source, dest);
1454  m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1455  G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1456  coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1457  coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1458  newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
1459  source, dest);
1460  ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1461  ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1462  gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
1463  source, dest);
1464  lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1465  lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1466  fail= false;
1467  random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
1468  DEBOUTLN (cerr, "fail= " << fail);
1469  CFList list;
1470  TIMING_START (gcd_recursion);
1471  G_random_element=
1472  modGCDFq (ppA (random_element, x), ppB (random_element, x),
1473  coF_random_element, coG_random_element, V_buf,
1474  list, topLevel);
1475  TIMING_END_AND_PRINT (gcd_recursion,
1476  "time for recursive call: ");
1477  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1478  alpha= V_buf;
1479  }
1480 
1481  if (!G_random_element.inCoeffDomain())
1482  d0= totaldegree (G_random_element, Variable(2),
1483  Variable (G_random_element.level()));
1484  else
1485  d0= 0;
1486 
1487  if (d0 == 0)
1488  {
1489  if (inextension)
1490  prune (cleanUp);
1491  coF= N (ppA*(cA/gcdcAcB));
1492  coG= N (ppB*(cB/gcdcAcB));
1493  return N(gcdcAcB);
1494  }
1495 
1496  if (d0 > d)
1497  {
1498  if (!find (l, random_element))
1499  l.append (random_element);
1500  continue;
1501  }
1502 
1503  G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
1504  *G_random_element;
1505 
1506  coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1507  *coF_random_element;
1508  coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1509  *coG_random_element;
1510 
1511  if (!G_random_element.inCoeffDomain())
1512  d0= totaldegree (G_random_element, Variable(2),
1513  Variable (G_random_element.level()));
1514  else
1515  d0= 0;
1516 
1517  if (d0 < d)
1518  {
1519  m= gcdlcAlcB;
1520  newtonPoly= 1;
1521  G_m= 0;
1522  d= d0;
1523  coF_m= 0;
1524  coG_m= 0;
1525  }
1526 
1527  TIMING_START (newton_interpolation);
1528  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1529  coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1530  coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1531  TIMING_END_AND_PRINT (newton_interpolation,
1532  "time for newton_interpolation: ");
1533 
1534  //termination test
1535  if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1536  {
1537  TIMING_START (termination_test);
1538  if (gcdlcAlcB.isOne())
1539  cH= 1;
1540  else
1541  cH= uni_content (H);
1542  ppH= H/cH;
1543  ppH /= Lc (ppH);
1544  CanonicalForm lcppH= gcdlcAlcB/cH;
1545  CanonicalForm ccoF= lcppH/Lc (lcppH);
1546  CanonicalForm ccoG= lcppH/Lc (lcppH);
1547  ppCoF= coF/ccoF;
1548  ppCoG= coG/ccoG;
1549  DEBOUTLN (cerr, "ppH= " << ppH);
1550  if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1551  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1552  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1553  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1554  {
1555  if (inextension)
1556  prune (cleanUp);
1557  coF= N ((cA/gcdcAcB)*ppCoF);
1558  coG= N ((cB/gcdcAcB)*ppCoG);
1559  TIMING_END_AND_PRINT (termination_test,
1560  "time for successful termination Fp: ");
1561  return N(gcdcAcB*ppH);
1562  }
1563  TIMING_END_AND_PRINT (termination_test,
1564  "time for unsuccessful termination Fp: ");
1565  }
1566 
1567  G_m= H;
1568  coF_m= coF;
1569  coG_m= coG;
1570  newtonPoly= newtonPoly*(x - random_element);
1571  m= m*(x - random_element);
1572  if (!find (l, random_element))
1573  l.append (random_element);
1574  } while (1);
1575 }
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition: cf_ops.cc:523
CanonicalForm Lc(const CanonicalForm &f)
CanonicalForm modGCDFq(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, Variable &alpha, CFList &l, bool &topLevel)
GCD of F and G over , l and topLevel are only used internally, output is monic based on Alg....
Definition: cfModGcd.cc:478
int myCompress(const CanonicalForm &F, const CanonicalForm &G, CFMap &M, CFMap &N, bool topLevel)
compressing two polynomials F and G, M is used for compressing, N to reverse the compression
Definition: cfModGcd.cc:91
static Variable chooseExtension(const Variable &alpha)
Definition: cfModGcd.cc:420
static CanonicalForm uni_lcoeff(const CanonicalForm &F)
compute the leading coefficient of F, where F is considered as an element of , order on is dp.
Definition: cfModGcd.cc:339
static CanonicalForm randomElement(const CanonicalForm &F, const Variable &alpha, CFList &list, bool &fail)
compute a random element a of , s.t. F(a) , F is a univariate polynomial, returns fail if there are...
Definition: cfModGcd.cc:379
static CanonicalForm FpRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
Definition: cfModGcd.cc:1171
static CanonicalForm newtonInterp(const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x)
Newton interpolation - Incremental algorithm. Given a list of values alpha_i and a list of polynomial...
Definition: cfModGcd.cc:364
bool terminationTest(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &coF, const CanonicalForm &coG, const CanonicalForm &cand)
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
CanonicalForm randomIrredpoly(int i, const Variable &x)
computes a random monic irreducible univariate polynomial in x over Fp of degree i via NTL/FLINT
Definition: cf_irred.cc:26
CanonicalForm mapPrimElem(const CanonicalForm &primElem, const Variable &alpha, const Variable &beta)
compute the image of a primitive element of in . We assume .
Definition: cf_map_ext.cc:450
CanonicalForm primitiveElement(const Variable &alpha, Variable &beta, bool &fail)
determine a primitive element of , is a primitive element of a field which is isomorphic to
Definition: cf_map_ext.cc:342
static CanonicalForm mapDown(const CanonicalForm &F, const Variable &alpha, const CanonicalForm &G, CFList &source, CFList &dest)
the CanonicalForm G is the output of map_up, returns F considered as an element over ,...
Definition: cf_map_ext.cc:123
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
Definition: cf_map_ext.cc:70
class CFMap
Definition: cf_map.h:85
CanonicalForm genOne() const
CF_NO_INLINE bool isOne() const
bool inBaseDomain() const
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
CanonicalForm H
Definition: facAbsFact.cc:60
CanonicalForm mipo
Definition: facAlgExt.cc:57
TIMING_END_AND_PRINT(fac_alg_resultant, "time to compute resultant0: ")
TIMING_START(fac_alg_resultant)
void FACTORY_PUBLIC prune(Variable &alpha)
Definition: variable.cc:261
void setMipo(const Variable &alpha, const CanonicalForm &mipo)
Definition: variable.cc:219

◆ modGCDFq() [1/2]

CanonicalForm modGCDFq ( const CanonicalForm F,
const CanonicalForm G,
CanonicalForm coF,
CanonicalForm coG,
Variable alpha,
CFList l,
bool &  topLevel 
)

GCD of F and G over $ F_{p}(\alpha ) $ , l and topLevel are only used internally, output is monic based on Alg. 7.2. as described in "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn.

Definition at line 478 of file cfModGcd.cc.

481 {
482  CanonicalForm A= F;
483  CanonicalForm B= G;
484  if (F.isZero() && degree(G) > 0)
485  {
486  coF= 0;
487  coG= Lc (G);
488  return G/Lc(G);
489  }
490  else if (G.isZero() && degree (F) > 0)
491  {
492  coF= Lc (F);
493  coG= 0;
494  return F/Lc(F);
495  }
496  else if (F.isZero() && G.isZero())
497  {
498  coF= coG= 0;
499  return F.genOne();
500  }
501  if (F.inBaseDomain() || G.inBaseDomain())
502  {
503  coF= F;
504  coG= G;
505  return F.genOne();
506  }
507  if (F.isUnivariate() && fdivides(F, G, coG))
508  {
509  coF= Lc (F);
510  return F/Lc(F);
511  }
512  if (G.isUnivariate() && fdivides(G, F, coF))
513  {
514  coG= Lc (G);
515  return G/Lc(G);
516  }
517  if (F == G)
518  {
519  coF= coG= Lc (F);
520  return F/Lc(F);
521  }
522 
523  CFMap M,N;
524  int best_level= myCompress (A, B, M, N, topLevel);
525 
526  if (best_level == 0)
527  {
528  coF= F;
529  coG= G;
530  return B.genOne();
531  }
532 
533  A= M(A);
534  B= M(B);
535 
536  Variable x= Variable(1);
537 
538  //univariate case
539  if (A.isUnivariate() && B.isUnivariate())
540  {
541  CanonicalForm result= gcd (A, B);
542  coF= N (A/result);
543  coG= N (B/result);
544  return N (result);
545  }
546 
547  CanonicalForm cA, cB; // content of A and B
548  CanonicalForm ppA, ppB; // primitive part of A and B
549  CanonicalForm gcdcAcB;
550 
551  cA = uni_content (A);
552  cB = uni_content (B);
553  gcdcAcB= gcd (cA, cB);
554  ppA= A/cA;
555  ppB= B/cB;
556 
557  CanonicalForm lcA, lcB; // leading coefficients of A and B
558  CanonicalForm gcdlcAlcB;
559 
560  lcA= uni_lcoeff (ppA);
561  lcB= uni_lcoeff (ppB);
562 
563  gcdlcAlcB= gcd (lcA, lcB);
564 
565  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
566 
567  if (d == 0)
568  {
569  coF= N (ppA*(cA/gcdcAcB));
570  coG= N (ppB*(cB/gcdcAcB));
571  return N(gcdcAcB);
572  }
573 
574  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
575  if (d0 < d)
576  d= d0;
577  if (d == 0)
578  {
579  coF= N (ppA*(cA/gcdcAcB));
580  coG= N (ppB*(cB/gcdcAcB));
581  return N(gcdcAcB);
582  }
583 
584  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
585  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
586  coG_m, ppCoF, ppCoG;
587 
588  newtonPoly= 1;
589  m= gcdlcAlcB;
590  G_m= 0;
591  coF= 0;
592  coG= 0;
593  H= 0;
594  bool fail= false;
595  topLevel= false;
596  bool inextension= false;
597  Variable V_buf= alpha, V_buf4= alpha;
598  CanonicalForm prim_elem, im_prim_elem;
599  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
600  CFList source, dest;
601  int bound1= degree (ppA, 1);
602  int bound2= degree (ppB, 1);
603  do
604  {
605  random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
606  if (fail)
607  {
608  source= CFList();
609  dest= CFList();
610 
611  Variable V_buf3= V_buf;
612  V_buf= chooseExtension (V_buf);
613  bool prim_fail= false;
614  Variable V_buf2;
615  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
616  if (V_buf4 == alpha)
617  prim_elem_alpha= prim_elem;
618 
619  if (V_buf3 != V_buf4)
620  {
621  m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
622  G_m= mapDown (G_m, prim_elem, im_prim_elem, V_buf4, source, dest);
623  coF_m= mapDown (coF_m, prim_elem, im_prim_elem, V_buf4, source, dest);
624  coG_m= mapDown (coG_m, prim_elem, im_prim_elem, V_buf4, source, dest);
625  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
626  source, dest);
627  ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
628  ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
629  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
630  source, dest);
631  lcA= mapDown (lcA, prim_elem, im_prim_elem, V_buf4, source, dest);
632  lcB= mapDown (lcB, prim_elem, im_prim_elem, V_buf4, source, dest);
633  for (CFListIterator i= l; i.hasItem(); i++)
634  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
635  source, dest);
636  }
637 
638  ASSERT (!prim_fail, "failure in integer factorizer");
639  if (prim_fail)
640  ; //ERROR
641  else
642  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
643 
644  if (V_buf4 == alpha)
645  im_prim_elem_alpha= im_prim_elem;
646  else
647  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
648  im_prim_elem, source, dest);
649  DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
650  DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
651  inextension= true;
652  for (CFListIterator i= l; i.hasItem(); i++)
653  i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
654  im_prim_elem, source, dest);
655  m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
656  G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
657  coF_m= mapUp (coF_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
658  coG_m= mapUp (coG_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
659  newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
660  source, dest);
661  ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
662  ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
663  gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
664  source, dest);
665  lcA= mapUp (lcA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
666  lcB= mapUp (lcB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
667 
668  fail= false;
669  random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
670  DEBOUTLN (cerr, "fail= " << fail);
671  CFList list;
672  TIMING_START (gcd_recursion);
673  G_random_element=
674  modGCDFq (ppA (random_element, x), ppB (random_element, x),
675  coF_random_element, coG_random_element, V_buf,
676  list, topLevel);
677  TIMING_END_AND_PRINT (gcd_recursion,
678  "time for recursive call: ");
679  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
680  V_buf4= V_buf;
681  }
682  else
683  {
684  CFList list;
685  TIMING_START (gcd_recursion);
686  G_random_element=
687  modGCDFq (ppA(random_element, x), ppB(random_element, x),
688  coF_random_element, coG_random_element, V_buf,
689  list, topLevel);
690  TIMING_END_AND_PRINT (gcd_recursion,
691  "time for recursive call: ");
692  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
693  }
694 
695  if (!G_random_element.inCoeffDomain())
696  d0= totaldegree (G_random_element, Variable(2),
697  Variable (G_random_element.level()));
698  else
699  d0= 0;
700 
701  if (d0 == 0)
702  {
703  if (inextension)
704  {
705  CFList u, v;
706  ppA= mapDown (ppA, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
707  ppB= mapDown (ppB, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
708  prune1 (alpha);
709  }
710  coF= N (ppA*(cA/gcdcAcB));
711  coG= N (ppB*(cB/gcdcAcB));
712  return N(gcdcAcB);
713  }
714  if (d0 > d)
715  {
716  if (!find (l, random_element))
717  l.append (random_element);
718  continue;
719  }
720 
721  G_random_element=
722  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
723  * G_random_element;
724  coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
725  *coF_random_element;
726  coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
727  *coG_random_element;
728 
729  if (!G_random_element.inCoeffDomain())
730  d0= totaldegree (G_random_element, Variable(2),
731  Variable (G_random_element.level()));
732  else
733  d0= 0;
734 
735  if (d0 < d)
736  {
737  m= gcdlcAlcB;
738  newtonPoly= 1;
739  G_m= 0;
740  d= d0;
741  coF_m= 0;
742  coG_m= 0;
743  }
744 
745  TIMING_START (newton_interpolation);
746  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
747  coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
748  coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
749  TIMING_END_AND_PRINT (newton_interpolation,
750  "time for newton interpolation: ");
751 
752  //termination test
753  if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
754  {
755  TIMING_START (termination_test);
756  if (gcdlcAlcB.isOne())
757  cH= 1;
758  else
759  cH= uni_content (H);
760  ppH= H/cH;
761  ppH /= Lc (ppH);
762  CanonicalForm lcppH= gcdlcAlcB/cH;
763  CanonicalForm ccoF= lcppH/Lc (lcppH);
764  CanonicalForm ccoG= lcppH/Lc (lcppH);
765  ppCoF= coF/ccoF;
766  ppCoG= coG/ccoG;
767  if (inextension)
768  {
769  if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
770  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
771  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
772  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
773  {
774  CFList u, v;
775  DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
776  ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
777  ppCoF= mapDown (ppCoF, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
778  ppCoG= mapDown (ppCoG, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
779  DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
780  coF= N ((cA/gcdcAcB)*ppCoF);
781  coG= N ((cB/gcdcAcB)*ppCoG);
782  TIMING_END_AND_PRINT (termination_test,
783  "time for successful termination test Fq: ");
784  prune1 (alpha);
785  return N(gcdcAcB*ppH);
786  }
787  }
788  else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
789  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
790  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
791  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
792  {
793  coF= N ((cA/gcdcAcB)*ppCoF);
794  coG= N ((cB/gcdcAcB)*ppCoG);
795  TIMING_END_AND_PRINT (termination_test,
796  "time for successful termination test Fq: ");
797  return N(gcdcAcB*ppH);
798  }
799  TIMING_END_AND_PRINT (termination_test,
800  "time for unsuccessful termination test Fq: ");
801  }
802 
803  G_m= H;
804  coF_m= coF;
805  coG_m= coG;
806  newtonPoly= newtonPoly*(x - random_element);
807  m= m*(x - random_element);
808  if (!find (l, random_element))
809  l.append (random_element);
810  } while (1);
811 }
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
void prune1(const Variable &alpha)
Definition: variable.cc:291

◆ modGCDFq() [2/2]

CanonicalForm modGCDFq ( const CanonicalForm F,
const CanonicalForm G,
Variable alpha,
CFList l,
bool &  topLevel 
)

Definition at line 462 of file cfModGcd.cc.

464 {
465  CanonicalForm dummy1, dummy2;
466  CanonicalForm result= modGCDFq (F, G, dummy1, dummy2, alpha, l,
467  topLevel);
468  return result;
469 }

◆ modGCDGF() [1/2]

CanonicalForm modGCDGF ( const CanonicalForm F,
const CanonicalForm G,
CanonicalForm coF,
CanonicalForm coG,
CFList l,
bool &  topLevel 
)

GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn Usually this algorithm will be faster than modGCDFq since GF has faster field arithmetics, however it might fail if the input is large since the size of the base field is bounded by 2^16, output is monic.

Definition at line 872 of file cfModGcd.cc.

875 {
876  CanonicalForm A= F;
877  CanonicalForm B= G;
878  if (F.isZero() && degree(G) > 0)
879  {
880  coF= 0;
881  coG= Lc (G);
882  return G/Lc(G);
883  }
884  else if (G.isZero() && degree (F) > 0)
885  {
886  coF= Lc (F);
887  coG= 0;
888  return F/Lc(F);
889  }
890  else if (F.isZero() && G.isZero())
891  {
892  coF= coG= 0;
893  return F.genOne();
894  }
895  if (F.inBaseDomain() || G.inBaseDomain())
896  {
897  coF= F;
898  coG= G;
899  return F.genOne();
900  }
901  if (F.isUnivariate() && fdivides(F, G, coG))
902  {
903  coF= Lc (F);
904  return F/Lc(F);
905  }
906  if (G.isUnivariate() && fdivides(G, F, coF))
907  {
908  coG= Lc (G);
909  return G/Lc(G);
910  }
911  if (F == G)
912  {
913  coF= coG= Lc (F);
914  return F/Lc(F);
915  }
916 
917  CFMap M,N;
918  int best_level= myCompress (A, B, M, N, topLevel);
919 
920  if (best_level == 0)
921  {
922  coF= F;
923  coG= G;
924  return B.genOne();
925  }
926 
927  A= M(A);
928  B= M(B);
929 
930  Variable x= Variable(1);
931 
932  //univariate case
933  if (A.isUnivariate() && B.isUnivariate())
934  {
935  CanonicalForm result= gcd (A, B);
936  coF= N (A/result);
937  coG= N (B/result);
938  return N (result);
939  }
940 
941  CanonicalForm cA, cB; // content of A and B
942  CanonicalForm ppA, ppB; // primitive part of A and B
943  CanonicalForm gcdcAcB;
944 
945  cA = uni_content (A);
946  cB = uni_content (B);
947  gcdcAcB= gcd (cA, cB);
948  ppA= A/cA;
949  ppB= B/cB;
950 
951  CanonicalForm lcA, lcB; // leading coefficients of A and B
952  CanonicalForm gcdlcAlcB;
953 
954  lcA= uni_lcoeff (ppA);
955  lcB= uni_lcoeff (ppB);
956 
957  gcdlcAlcB= gcd (lcA, lcB);
958 
959  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
960  if (d == 0)
961  {
962  coF= N (ppA*(cA/gcdcAcB));
963  coG= N (ppB*(cB/gcdcAcB));
964  return N(gcdcAcB);
965  }
966  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
967  if (d0 < d)
968  d= d0;
969  if (d == 0)
970  {
971  coF= N (ppA*(cA/gcdcAcB));
972  coG= N (ppB*(cB/gcdcAcB));
973  return N(gcdcAcB);
974  }
975 
976  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
977  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
978  coG_m, ppCoF, ppCoG;
979 
980  newtonPoly= 1;
981  m= gcdlcAlcB;
982  G_m= 0;
983  coF= 0;
984  coG= 0;
985  H= 0;
986  bool fail= false;
987  topLevel= false;
988  bool inextension= false;
989  int p=-1;
990  int k= getGFDegree();
991  int kk;
992  int expon;
993  char gf_name_buf= gf_name;
994  int bound1= degree (ppA, 1);
995  int bound2= degree (ppB, 1);
996  do
997  {
998  random_element= GFRandomElement (m*lcA*lcB, l, fail);
999  if (fail)
1000  {
1001  p= getCharacteristic();
1002  expon= 2;
1003  kk= getGFDegree();
1004  if (ipower (p, kk*expon) < (1 << 16))
1005  setCharacteristic (p, kk*(int)expon, 'b');
1006  else
1007  {
1008  expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
1009  ASSERT (expon >= 2, "not enough points in modGCDGF");
1010  setCharacteristic (p, (int)(kk*expon), 'b');
1011  }
1012  inextension= true;
1013  fail= false;
1014  for (CFListIterator i= l; i.hasItem(); i++)
1015  i.getItem()= GFMapUp (i.getItem(), kk);
1016  m= GFMapUp (m, kk);
1017  G_m= GFMapUp (G_m, kk);
1018  newtonPoly= GFMapUp (newtonPoly, kk);
1019  coF_m= GFMapUp (coF_m, kk);
1020  coG_m= GFMapUp (coG_m, kk);
1021  ppA= GFMapUp (ppA, kk);
1022  ppB= GFMapUp (ppB, kk);
1023  gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
1024  lcA= GFMapUp (lcA, kk);
1025  lcB= GFMapUp (lcB, kk);
1026  random_element= GFRandomElement (m*lcA*lcB, l, fail);
1027  DEBOUTLN (cerr, "fail= " << fail);
1028  CFList list;
1029  TIMING_START (gcd_recursion);
1030  G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
1031  coF_random_element, coG_random_element,
1032  list, topLevel);
1033  TIMING_END_AND_PRINT (gcd_recursion,
1034  "time for recursive call: ");
1035  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1036  }
1037  else
1038  {
1039  CFList list;
1040  TIMING_START (gcd_recursion);
1041  G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
1042  coF_random_element, coG_random_element,
1043  list, topLevel);
1044  TIMING_END_AND_PRINT (gcd_recursion,
1045  "time for recursive call: ");
1046  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1047  }
1048 
1049  if (!G_random_element.inCoeffDomain())
1050  d0= totaldegree (G_random_element, Variable(2),
1051  Variable (G_random_element.level()));
1052  else
1053  d0= 0;
1054 
1055  if (d0 == 0)
1056  {
1057  if (inextension)
1058  {
1059  ppA= GFMapDown (ppA, k);
1060  ppB= GFMapDown (ppB, k);
1061  setCharacteristic (p, k, gf_name_buf);
1062  }
1063  coF= N (ppA*(cA/gcdcAcB));
1064  coG= N (ppB*(cB/gcdcAcB));
1065  return N(gcdcAcB);
1066  }
1067  if (d0 > d)
1068  {
1069  if (!find (l, random_element))
1070  l.append (random_element);
1071  continue;
1072  }
1073 
1074  G_random_element=
1075  (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
1076  G_random_element;
1077 
1078  coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1079  *coF_random_element;
1080  coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1081  *coG_random_element;
1082 
1083  if (!G_random_element.inCoeffDomain())
1084  d0= totaldegree (G_random_element, Variable(2),
1085  Variable (G_random_element.level()));
1086  else
1087  d0= 0;
1088 
1089  if (d0 < d)
1090  {
1091  m= gcdlcAlcB;
1092  newtonPoly= 1;
1093  G_m= 0;
1094  d= d0;
1095  coF_m= 0;
1096  coG_m= 0;
1097  }
1098 
1099  TIMING_START (newton_interpolation);
1100  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1101  coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1102  coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1103  TIMING_END_AND_PRINT (newton_interpolation,
1104  "time for newton interpolation: ");
1105 
1106  //termination test
1107  if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1108  {
1109  TIMING_START (termination_test);
1110  if (gcdlcAlcB.isOne())
1111  cH= 1;
1112  else
1113  cH= uni_content (H);
1114  ppH= H/cH;
1115  ppH /= Lc (ppH);
1116  CanonicalForm lcppH= gcdlcAlcB/cH;
1117  CanonicalForm ccoF= lcppH/Lc (lcppH);
1118  CanonicalForm ccoG= lcppH/Lc (lcppH);
1119  ppCoF= coF/ccoF;
1120  ppCoG= coG/ccoG;
1121  if (inextension)
1122  {
1123  if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1124  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1125  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1126  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1127  {
1128  DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
1129  ppH= GFMapDown (ppH, k);
1130  ppCoF= GFMapDown (ppCoF, k);
1131  ppCoG= GFMapDown (ppCoG, k);
1132  DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
1133  coF= N ((cA/gcdcAcB)*ppCoF);
1134  coG= N ((cB/gcdcAcB)*ppCoG);
1135  setCharacteristic (p, k, gf_name_buf);
1136  TIMING_END_AND_PRINT (termination_test,
1137  "time for successful termination GF: ");
1138  return N(gcdcAcB*ppH);
1139  }
1140  }
1141  else
1142  {
1143  if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1144  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1145  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1146  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1147  {
1148  coF= N ((cA/gcdcAcB)*ppCoF);
1149  coG= N ((cB/gcdcAcB)*ppCoG);
1150  TIMING_END_AND_PRINT (termination_test,
1151  "time for successful termination GF: ");
1152  return N(gcdcAcB*ppH);
1153  }
1154  }
1155  TIMING_END_AND_PRINT (termination_test,
1156  "time for unsuccessful termination GF: ");
1157  }
1158 
1159  G_m= H;
1160  coF_m= coF;
1161  coG_m= coG;
1162  newtonPoly= newtonPoly*(x - random_element);
1163  m= m*(x - random_element);
1164  if (!find (l, random_element))
1165  l.append (random_element);
1166  } while (1);
1167 }
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
static CanonicalForm GFRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
compute a random element a of GF, s.t. F(a) , F is a univariate polynomial, returns fail if there ar...
Definition: cfModGcd.cc:819
CanonicalForm modGCDGF(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, CFList &l, bool &topLevel)
GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for Computer Algebra" by Gedde...
Definition: cfModGcd.cc:872
CanonicalForm GFMapDown(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition: cf_map_ext.cc:276
CanonicalForm GFMapUp(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition: cf_map_ext.cc:240
VAR char gf_name
Definition: gfops.cc:52
gmp_float log(const gmp_float &a)
Definition: mpr_complex.cc:343
const signed long floor(const ampf< Precision > &x)
Definition: amp.h:773

◆ modGCDGF() [2/2]

CanonicalForm modGCDGF ( const CanonicalForm F,
const CanonicalForm G,
CFList l,
bool &  topLevel 
)

Definition at line 858 of file cfModGcd.cc.

860 {
861  CanonicalForm dummy1, dummy2;
862  CanonicalForm result= modGCDGF (F, G, dummy1, dummy2, l, topLevel);
863  return result;
864 }

◆ monicSparseInterpol()

CanonicalForm monicSparseInterpol ( const CanonicalForm F,
const CanonicalForm G,
const CanonicalForm skeleton,
const Variable alpha,
bool &  fail,
CFArray *&  coeffMonoms,
CFArray Monoms 
)

Definition at line 2199 of file cfModGcd.cc.

2203 {
2204  CanonicalForm A= F;
2205  CanonicalForm B= G;
2206  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2207  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2208  else if (F.isZero() && G.isZero()) return F.genOne();
2209  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2210  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2211  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2212  if (F == G) return F/Lc(F);
2213 
2214  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2215  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2216 
2217  CFMap M,N;
2218  int best_level= myCompress (A, B, M, N, false);
2219 
2220  if (best_level == 0)
2221  return B.genOne();
2222 
2223  A= M(A);
2224  B= M(B);
2225 
2226  Variable x= Variable (1);
2227 
2228  //univariate case
2229  if (A.isUnivariate() && B.isUnivariate())
2230  return N (gcd (A, B));
2231 
2232  CanonicalForm skel= M(skeleton);
2233  CanonicalForm cA, cB; // content of A and B
2234  CanonicalForm ppA, ppB; // primitive part of A and B
2235  CanonicalForm gcdcAcB;
2236  cA = uni_content (A);
2237  cB = uni_content (B);
2238  gcdcAcB= gcd (cA, cB);
2239  ppA= A/cA;
2240  ppB= B/cB;
2241 
2242  CanonicalForm lcA, lcB; // leading coefficients of A and B
2243  CanonicalForm gcdlcAlcB;
2244  lcA= uni_lcoeff (ppA);
2245  lcB= uni_lcoeff (ppB);
2246 
2247  if (fdivides (lcA, lcB))
2248  {
2249  if (fdivides (A, B))
2250  return F/Lc(F);
2251  }
2252  if (fdivides (lcB, lcA))
2253  {
2254  if (fdivides (B, A))
2255  return G/Lc(G);
2256  }
2257 
2258  gcdlcAlcB= gcd (lcA, lcB);
2259  int skelSize= size (skel, skel.mvar());
2260 
2261  int j= 0;
2262  int biggestSize= 0;
2263 
2264  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2265  biggestSize= tmax (biggestSize, size (i.coeff()));
2266 
2267  CanonicalForm g, Aeval, Beval;
2268 
2270  bool evalFail= false;
2271  CFList list;
2272  bool GF= false;
2273  CanonicalForm LCA= LC (A);
2274  CanonicalForm tmp;
2275  CFArray gcds= CFArray (biggestSize);
2276  CFList * pEvalPoints= new CFList [biggestSize];
2277  Variable V_buf= alpha, V_buf4= alpha;
2278  CFList source, dest;
2279  CanonicalForm prim_elem, im_prim_elem;
2280  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2281  for (int i= 0; i < biggestSize; i++)
2282  {
2283  if (i == 0)
2284  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2285  list);
2286  else
2287  {
2288  mult (evalPoints, pEvalPoints [0]);
2289  eval (A, B, Aeval, Beval, evalPoints);
2290  }
2291 
2292  if (evalFail)
2293  {
2294  if (V_buf.level() != 1)
2295  {
2296  do
2297  {
2298  Variable V_buf3= V_buf;
2299  V_buf= chooseExtension (V_buf);
2300  source= CFList();
2301  dest= CFList();
2302 
2303  bool prim_fail= false;
2304  Variable V_buf2;
2305  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2306  if (V_buf4 == alpha && alpha.level() != 1)
2307  prim_elem_alpha= prim_elem;
2308 
2309  ASSERT (!prim_fail, "failure in integer factorizer");
2310  if (prim_fail)
2311  ; //ERROR
2312  else
2313  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2314 
2315  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2316  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2317 
2318  if (V_buf4 == alpha && alpha.level() != 1)
2319  im_prim_elem_alpha= im_prim_elem;
2320  else if (alpha.level() != 1)
2321  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2322  prim_elem, im_prim_elem, source, dest);
2323 
2324  for (CFListIterator j= list; j.hasItem(); j++)
2325  j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2326  im_prim_elem, source, dest);
2327  for (int k= 0; k < i; k++)
2328  {
2329  for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2330  j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2331  im_prim_elem, source, dest);
2332  gcds[k]= mapUp (gcds[k], V_buf4, V_buf, prim_elem, im_prim_elem,
2333  source, dest);
2334  }
2335 
2336  if (alpha.level() != 1)
2337  {
2338  A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2339  B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2340  }
2341  V_buf4= V_buf;
2342  evalFail= false;
2343  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2344  evalFail, list);
2345  } while (evalFail);
2346  }
2347  else
2348  {
2350  int deg= 2;
2351  bool initialized= false;
2352  do
2353  {
2354  mipo= randomIrredpoly (deg, x);
2355  if (initialized)
2356  setMipo (V_buf, mipo);
2357  else
2358  V_buf= rootOf (mipo);
2359  evalFail= false;
2360  initialized= true;
2361  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2362  evalFail, list);
2363  deg++;
2364  } while (evalFail);
2365  V_buf4= V_buf;
2366  }
2367  }
2368 
2369  g= gcd (Aeval, Beval);
2370  g /= Lc (g);
2371 
2372  if (degree (g) != degree (skel) || (skelSize != size (g)))
2373  {
2374  delete[] pEvalPoints;
2375  fail= true;
2376  if (alpha.level() != 1 && V_buf != alpha)
2377  prune1 (alpha);
2378  return 0;
2379  }
2380  CFIterator l= skel;
2381  for (CFIterator k= g; k.hasTerms(); k++, l++)
2382  {
2383  if (k.exp() != l.exp())
2384  {
2385  delete[] pEvalPoints;
2386  fail= true;
2387  if (alpha.level() != 1 && V_buf != alpha)
2388  prune1 (alpha);
2389  return 0;
2390  }
2391  }
2392  pEvalPoints[i]= evalPoints;
2393  gcds[i]= g;
2394 
2395  tmp= 0;
2396  int j= 0;
2397  for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2398  tmp += k.getItem()*power (x, j);
2399  list.append (tmp);
2400  }
2401 
2402  if (Monoms.size() == 0)
2403  Monoms= getMonoms (skel);
2404 
2405  coeffMonoms= new CFArray [skelSize];
2406  j= 0;
2407  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2408  coeffMonoms[j]= getMonoms (i.coeff());
2409 
2410  CFArray* pL= new CFArray [skelSize];
2411  CFArray* pM= new CFArray [skelSize];
2412  for (int i= 0; i < biggestSize; i++)
2413  {
2414  CFIterator l= gcds [i];
2415  evalPoints= pEvalPoints [i];
2416  for (int k= 0; k < skelSize; k++, l++)
2417  {
2418  if (i == 0)
2419  pL[k]= CFArray (biggestSize);
2420  pL[k] [i]= l.coeff();
2421 
2422  if (i == 0)
2423  pM[k]= evaluate (coeffMonoms [k], evalPoints);
2424  }
2425  }
2426 
2427  CFArray solution;
2428  CanonicalForm result= 0;
2429  int ind= 0;
2430  CFArray bufArray;
2431  CFMatrix Mat;
2432  for (int k= 0; k < skelSize; k++)
2433  {
2434  if (biggestSize != coeffMonoms[k].size())
2435  {
2436  bufArray= CFArray (coeffMonoms[k].size());
2437  for (int i= 0; i < coeffMonoms[k].size(); i++)
2438  bufArray [i]= pL[k] [i];
2439  solution= solveGeneralVandermonde (pM [k], bufArray);
2440  }
2441  else
2442  solution= solveGeneralVandermonde (pM [k], pL[k]);
2443 
2444  if (solution.size() == 0)
2445  {
2446  delete[] pEvalPoints;
2447  delete[] pM;
2448  delete[] pL;
2449  delete[] coeffMonoms;
2450  fail= true;
2451  if (alpha.level() != 1 && V_buf != alpha)
2452  prune1 (alpha);
2453  return 0;
2454  }
2455  for (int l= 0; l < solution.size(); l++)
2456  result += solution[l]*Monoms [ind + l];
2457  ind += solution.size();
2458  }
2459 
2460  delete[] pEvalPoints;
2461  delete[] pM;
2462  delete[] pL;
2463  delete[] coeffMonoms;
2464 
2465  if (alpha.level() != 1 && V_buf != alpha)
2466  {
2467  CFList u, v;
2468  result= mapDown (result, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2469  prune1 (alpha);
2470  }
2471 
2472  result= N(result);
2473  if (fdivides (result, F) && fdivides (result, G))
2474  return result;
2475  else
2476  {
2477  fail= true;
2478  return 0;
2479  }
2480 }
void mult(CFList &L1, const CFList &L2)
multiply two lists componentwise
Definition: cfModGcd.cc:2176
CFArray evaluate(const CFArray &A, const CFList &evalPoints)
Definition: cfModGcd.cc:2031
CFArray solveGeneralVandermonde(const CFArray &M, const CFArray &A)
Definition: cfModGcd.cc:1632
void eval(const CanonicalForm &A, const CanonicalForm &B, CanonicalForm &Aeval, CanonicalForm &Beval, const CFList &L)
Definition: cfModGcd.cc:2185
CFList evaluationPoints(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Feval, CanonicalForm &Geval, const CanonicalForm &LCF, const bool &GF, const Variable &alpha, bool &fail, CFList &list)
Definition: cfModGcd.cc:2048

◆ mult()

void mult ( CFList L1,
const CFList L2 
)

multiply two lists componentwise

Definition at line 2176 of file cfModGcd.cc.

2177 {
2178  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2179 
2180  CFListIterator j= L2;
2181  for (CFListIterator i= L1; i.hasItem(); i++, j++)
2182  i.getItem() *= j.getItem();
2183 }

◆ myCompress()

int myCompress ( const CanonicalForm F,
const CanonicalForm G,
CFMap M,
CFMap N,
bool  topLevel 
)

compressing two polynomials F and G, M is used for compressing, N to reverse the compression

Definition at line 91 of file cfModGcd.cc.

93 {
94  int n= tmax (F.level(), G.level());
95  int * degsf= NEW_ARRAY(int,n + 1);
96  int * degsg= NEW_ARRAY(int,n + 1);
97 
98  for (int i = n; i >= 0; i--)
99  degsf[i]= degsg[i]= 0;
100 
101  degsf= degrees (F, degsf);
102  degsg= degrees (G, degsg);
103 
104  int both_non_zero= 0;
105  int f_zero= 0;
106  int g_zero= 0;
107  int both_zero= 0;
108 
109  if (topLevel)
110  {
111  for (int i= 1; i <= n; i++)
112  {
113  if (degsf[i] != 0 && degsg[i] != 0)
114  {
115  both_non_zero++;
116  continue;
117  }
118  if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
119  {
120  f_zero++;
121  continue;
122  }
123  if (degsg[i] == 0 && degsf[i] && i <= F.level())
124  {
125  g_zero++;
126  continue;
127  }
128  }
129 
130  if (both_non_zero == 0)
131  {
134  return 0;
135  }
136 
137  // map Variables which do not occur in both polynomials to higher levels
138  int k= 1;
139  int l= 1;
140  for (int i= 1; i <= n; i++)
141  {
142  if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
143  {
144  if (k + both_non_zero != i)
145  {
146  M.newpair (Variable (i), Variable (k + both_non_zero));
147  N.newpair (Variable (k + both_non_zero), Variable (i));
148  }
149  k++;
150  }
151  if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
152  {
153  if (l + g_zero + both_non_zero != i)
154  {
155  M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
156  N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
157  }
158  l++;
159  }
160  }
161 
162  // sort Variables x_{i} in increasing order of
163  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
164  int m= tmax (F.level(), G.level());
165  int min_max_deg;
166  k= both_non_zero;
167  l= 0;
168  int i= 1;
169  while (k > 0)
170  {
171  if (degsf [i] != 0 && degsg [i] != 0)
172  min_max_deg= tmax (degsf[i], degsg[i]);
173  else
174  min_max_deg= 0;
175  while (min_max_deg == 0)
176  {
177  i++;
178  if (degsf [i] != 0 && degsg [i] != 0)
179  min_max_deg= tmax (degsf[i], degsg[i]);
180  else
181  min_max_deg= 0;
182  }
183  for (int j= i + 1; j <= m; j++)
184  {
185  if (degsf[j] != 0 && degsg [j] != 0 &&
186  tmax (degsf[j],degsg[j]) <= min_max_deg)
187  {
188  min_max_deg= tmax (degsf[j], degsg[j]);
189  l= j;
190  }
191  }
192  if (l != 0)
193  {
194  if (l != k)
195  {
196  M.newpair (Variable (l), Variable(k));
197  N.newpair (Variable (k), Variable(l));
198  degsf[l]= 0;
199  degsg[l]= 0;
200  l= 0;
201  }
202  else
203  {
204  degsf[l]= 0;
205  degsg[l]= 0;
206  l= 0;
207  }
208  }
209  else if (l == 0)
210  {
211  if (i != k)
212  {
213  M.newpair (Variable (i), Variable (k));
214  N.newpair (Variable (k), Variable (i));
215  degsf[i]= 0;
216  degsg[i]= 0;
217  }
218  else
219  {
220  degsf[i]= 0;
221  degsg[i]= 0;
222  }
223  i++;
224  }
225  k--;
226  }
227  }
228  else
229  {
230  //arrange Variables such that no gaps occur
231  for (int i= 1; i <= n; i++)
232  {
233  if (degsf[i] == 0 && degsg[i] == 0)
234  {
235  both_zero++;
236  continue;
237  }
238  else
239  {
240  if (both_zero != 0)
241  {
242  M.newpair (Variable (i), Variable (i - both_zero));
243  N.newpair (Variable (i - both_zero), Variable (i));
244  }
245  }
246  }
247  }
248 
251 
252  return 1;
253 }
int * degrees(const CanonicalForm &f, int *degs=0)
int * degrees ( const CanonicalForm & f, int * degs )
Definition: cf_ops.cc:493
int * degsf
Definition: cfEzgcd.cc:59
int f_zero
Definition: cfEzgcd.cc:69
const CanonicalForm CFMap CFMap int & both_non_zero
Definition: cfEzgcd.cc:57
int g_zero
Definition: cfEzgcd.cc:70
int * degsg
Definition: cfEzgcd.cc:60
int both_zero
Definition: cfGcdAlgExt.cc:71
#define DELETE_ARRAY(P)
Definition: cf_defs.h:65
#define NEW_ARRAY(T, N)
Definition: cf_defs.h:64

◆ newtonInterp()

static CanonicalForm newtonInterp ( const CanonicalForm alpha,
const CanonicalForm u,
const CanonicalForm newtonPoly,
const CanonicalForm oldInterPoly,
const Variable x 
)
inlinestatic

Newton interpolation - Incremental algorithm. Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n, computes the interpolation polynomial assuming that the polynomials in u are the results of evaluating the variabe x of the unknown polynomial at the alpha values. This incremental version receives only the values of alpha_n and u_n and the previous interpolation polynomial for points 1 <= i <= n-1, and computes the polynomial interpolating in all the points. newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})

Definition at line 364 of file cfModGcd.cc.

367 {
368  CanonicalForm interPoly;
369 
370  interPoly= oldInterPoly + ((u - oldInterPoly(alpha, x))/newtonPoly(alpha, x))
371  *newtonPoly;
372  return interPoly;
373 }

◆ nonMonicSparseInterpol()

CanonicalForm nonMonicSparseInterpol ( const CanonicalForm F,
const CanonicalForm G,
const CanonicalForm skeleton,
const Variable alpha,
bool &  fail,
CFArray *&  coeffMonoms,
CFArray Monoms 
)

Definition at line 2483 of file cfModGcd.cc.

2487 {
2488  CanonicalForm A= F;
2489  CanonicalForm B= G;
2490  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2491  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2492  else if (F.isZero() && G.isZero()) return F.genOne();
2493  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2494  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2495  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2496  if (F == G) return F/Lc(F);
2497 
2498  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2499  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2500 
2501  CFMap M,N;
2502  int best_level= myCompress (A, B, M, N, false);
2503 
2504  if (best_level == 0)
2505  return B.genOne();
2506 
2507  A= M(A);
2508  B= M(B);
2509 
2510  Variable x= Variable (1);
2511 
2512  //univariate case
2513  if (A.isUnivariate() && B.isUnivariate())
2514  return N (gcd (A, B));
2515 
2516  CanonicalForm skel= M(skeleton);
2517 
2518  CanonicalForm cA, cB; // content of A and B
2519  CanonicalForm ppA, ppB; // primitive part of A and B
2520  CanonicalForm gcdcAcB;
2521  cA = uni_content (A);
2522  cB = uni_content (B);
2523  gcdcAcB= gcd (cA, cB);
2524  ppA= A/cA;
2525  ppB= B/cB;
2526 
2527  CanonicalForm lcA, lcB; // leading coefficients of A and B
2528  CanonicalForm gcdlcAlcB;
2529  lcA= uni_lcoeff (ppA);
2530  lcB= uni_lcoeff (ppB);
2531 
2532  if (fdivides (lcA, lcB))
2533  {
2534  if (fdivides (A, B))
2535  return F/Lc(F);
2536  }
2537  if (fdivides (lcB, lcA))
2538  {
2539  if (fdivides (B, A))
2540  return G/Lc(G);
2541  }
2542 
2543  gcdlcAlcB= gcd (lcA, lcB);
2544  int skelSize= size (skel, skel.mvar());
2545 
2546  int j= 0;
2547  int biggestSize= 0;
2548  int bufSize;
2549  int numberUni= 0;
2550  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2551  {
2552  bufSize= size (i.coeff());
2553  biggestSize= tmax (biggestSize, bufSize);
2554  numberUni += bufSize;
2555  }
2556  numberUni--;
2557  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2558  biggestSize= tmax (biggestSize , numberUni);
2559 
2560  numberUni= biggestSize + size (LC(skel)) - 1;
2561  int biggestSize2= tmax (numberUni, biggestSize);
2562 
2563  CanonicalForm g, Aeval, Beval;
2564 
2566  CFArray coeffEval;
2567  bool evalFail= false;
2568  CFList list;
2569  bool GF= false;
2570  CanonicalForm LCA= LC (A);
2571  CanonicalForm tmp;
2572  CFArray gcds= CFArray (biggestSize);
2573  CFList * pEvalPoints= new CFList [biggestSize];
2574  Variable V_buf= alpha, V_buf4= alpha;
2575  CFList source, dest;
2576  CanonicalForm prim_elem, im_prim_elem;
2577  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2578  for (int i= 0; i < biggestSize; i++)
2579  {
2580  if (i == 0)
2581  {
2582  if (getCharacteristic() > 3)
2583  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2584  evalFail, list);
2585  else
2586  evalFail= true;
2587 
2588  if (evalFail)
2589  {
2590  if (V_buf.level() != 1)
2591  {
2592  do
2593  {
2594  Variable V_buf3= V_buf;
2595  V_buf= chooseExtension (V_buf);
2596  source= CFList();
2597  dest= CFList();
2598 
2599  bool prim_fail= false;
2600  Variable V_buf2;
2601  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2602  if (V_buf4 == alpha && alpha.level() != 1)
2603  prim_elem_alpha= prim_elem;
2604 
2605  ASSERT (!prim_fail, "failure in integer factorizer");
2606  if (prim_fail)
2607  ; //ERROR
2608  else
2609  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2610 
2611  DEBOUTLN (cerr, "getMipo (V_buf)= " << getMipo (V_buf));
2612  DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
2613 
2614  if (V_buf4 == alpha && alpha.level() != 1)
2615  im_prim_elem_alpha= im_prim_elem;
2616  else if (alpha.level() != 1)
2617  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2618  prim_elem, im_prim_elem, source, dest);
2619 
2620  for (CFListIterator i= list; i.hasItem(); i++)
2621  i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
2622  im_prim_elem, source, dest);
2623  if (alpha.level() != 1)
2624  {
2625  A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2626  B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2627  }
2628  evalFail= false;
2629  V_buf4= V_buf;
2630  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2631  evalFail, list);
2632  } while (evalFail);
2633  }
2634  else
2635  {
2637  int deg= 2;
2638  bool initialized= false;
2639  do
2640  {
2641  mipo= randomIrredpoly (deg, x);
2642  if (initialized)
2643  setMipo (V_buf, mipo);
2644  else
2645  V_buf= rootOf (mipo);
2646  evalFail= false;
2647  initialized= true;
2648  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2649  evalFail, list);
2650  deg++;
2651  } while (evalFail);
2652  V_buf4= V_buf;
2653  }
2654  }
2655  }
2656  else
2657  {
2658  mult (evalPoints, pEvalPoints[0]);
2659  eval (A, B, Aeval, Beval, evalPoints);
2660  }
2661 
2662  g= gcd (Aeval, Beval);
2663  g /= Lc (g);
2664 
2665  if (degree (g) != degree (skel) || (skelSize != size (g)))
2666  {
2667  delete[] pEvalPoints;
2668  fail= true;
2669  if (alpha.level() != 1 && V_buf != alpha)
2670  prune1 (alpha);
2671  return 0;
2672  }
2673  CFIterator l= skel;
2674  for (CFIterator k= g; k.hasTerms(); k++, l++)
2675  {
2676  if (k.exp() != l.exp())
2677  {
2678  delete[] pEvalPoints;
2679  fail= true;
2680  if (alpha.level() != 1 && V_buf != alpha)
2681  prune1 (alpha);
2682  return 0;
2683  }
2684  }
2685  pEvalPoints[i]= evalPoints;
2686  gcds[i]= g;
2687 
2688  tmp= 0;
2689  int j= 0;
2690  for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2691  tmp += k.getItem()*power (x, j);
2692  list.append (tmp);
2693  }
2694 
2695  if (Monoms.size() == 0)
2696  Monoms= getMonoms (skel);
2697 
2698  coeffMonoms= new CFArray [skelSize];
2699 
2700  j= 0;
2701  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2702  coeffMonoms[j]= getMonoms (i.coeff());
2703 
2704  int minimalColumnsIndex;
2705  if (skelSize > 1)
2706  minimalColumnsIndex= 1;
2707  else
2708  minimalColumnsIndex= 0;
2709  int minimalColumns=-1;
2710 
2711  CFArray* pM= new CFArray [skelSize];
2712  CFMatrix Mat;
2713  // find the Matrix with minimal number of columns
2714  for (int i= 0; i < skelSize; i++)
2715  {
2716  pM[i]= CFArray (coeffMonoms[i].size());
2717  if (i == 1)
2718  minimalColumns= coeffMonoms[i].size();
2719  if (i > 1)
2720  {
2721  if (minimalColumns > coeffMonoms[i].size())
2722  {
2723  minimalColumns= coeffMonoms[i].size();
2724  minimalColumnsIndex= i;
2725  }
2726  }
2727  }
2728  CFMatrix* pMat= new CFMatrix [2];
2729  pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2730  pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2731  CFArray* pL= new CFArray [skelSize];
2732  for (int i= 0; i < biggestSize; i++)
2733  {
2734  CFIterator l= gcds [i];
2735  evalPoints= pEvalPoints [i];
2736  for (int k= 0; k < skelSize; k++, l++)
2737  {
2738  if (i == 0)
2739  pL[k]= CFArray (biggestSize);
2740  pL[k] [i]= l.coeff();
2741 
2742  if (i == 0 && k != 0 && k != minimalColumnsIndex)
2743  pM[k]= evaluate (coeffMonoms[k], evalPoints);
2744  else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
2745  coeffEval= evaluate (coeffMonoms[k], evalPoints);
2746  if (i > 0 && (k == 0 || k == minimalColumnsIndex))
2747  coeffEval= evaluate (coeffMonoms[k], evalPoints);
2748 
2749  if (k == 0)
2750  {
2751  if (pMat[k].rows() >= i + 1)
2752  {
2753  for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2754  pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2755  }
2756  }
2757  if (k == minimalColumnsIndex)
2758  {
2759  if (pMat[1].rows() >= i + 1)
2760  {
2761  for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2762  pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2763  }
2764  }
2765  }
2766  }
2767 
2768  CFArray solution;
2769  CanonicalForm result= 0;
2770  int ind= 1;
2771  int matRows, matColumns;
2772  matRows= pMat[1].rows();
2773  matColumns= pMat[0].columns() - 1;
2774  matColumns += pMat[1].columns();
2775 
2776  Mat= CFMatrix (matRows, matColumns);
2777  for (int i= 1; i <= matRows; i++)
2778  for (int j= 1; j <= pMat[1].columns(); j++)
2779  Mat (i, j)= pMat[1] (i, j);
2780 
2781  for (int j= pMat[1].columns() + 1; j <= matColumns;
2782  j++, ind++)
2783  {
2784  for (int i= 1; i <= matRows; i++)
2785  Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2786  }
2787 
2788  CFArray firstColumn= CFArray (pMat[0].rows());
2789  for (int i= 0; i < pMat[0].rows(); i++)
2790  firstColumn [i]= pMat[0] (i + 1, 1);
2791 
2792  CFArray bufArray= pL[minimalColumnsIndex];
2793 
2794  for (int i= 0; i < biggestSize; i++)
2795  pL[minimalColumnsIndex] [i] *= firstColumn[i];
2796 
2797  CFMatrix bufMat= pMat[1];
2798  pMat[1]= Mat;
2799 
2800  if (V_buf.level() != 1)
2801  solution= solveSystemFq (pMat[1],
2802  pL[minimalColumnsIndex], V_buf);
2803  else
2804  solution= solveSystemFp (pMat[1],
2805  pL[minimalColumnsIndex]);
2806 
2807  if (solution.size() == 0)
2808  { //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2809  CFMatrix bufMat0= pMat[0];
2810  delete [] pMat;
2811  pMat= new CFMatrix [skelSize];
2812  pL[minimalColumnsIndex]= bufArray;
2813  CFList* bufpEvalPoints= NULL;
2814  CFArray bufGcds;
2815  if (biggestSize != biggestSize2)
2816  {
2817  bufpEvalPoints= pEvalPoints;
2818  pEvalPoints= new CFList [biggestSize2];
2819  bufGcds= gcds;
2820  gcds= CFArray (biggestSize2);
2821  for (int i= 0; i < biggestSize; i++)
2822  {
2823  pEvalPoints[i]= bufpEvalPoints [i];
2824  gcds[i]= bufGcds[i];
2825  }
2826  for (int i= 0; i < biggestSize2 - biggestSize; i++)
2827  {
2828  mult (evalPoints, pEvalPoints[0]);
2829  eval (A, B, Aeval, Beval, evalPoints);
2830  g= gcd (Aeval, Beval);
2831  g /= Lc (g);
2832  if (degree (g) != degree (skel) || (skelSize != size (g)))
2833  {
2834  delete[] pEvalPoints;
2835  delete[] pMat;
2836  delete[] pL;
2837  delete[] coeffMonoms;
2838  delete[] pM;
2839  if (bufpEvalPoints != NULL)
2840  delete [] bufpEvalPoints;
2841  fail= true;
2842  if (alpha.level() != 1 && V_buf != alpha)
2843  prune1 (alpha);
2844  return 0;
2845  }
2846  CFIterator l= skel;
2847  for (CFIterator k= g; k.hasTerms(); k++, l++)
2848  {
2849  if (k.exp() != l.exp())
2850  {
2851  delete[] pEvalPoints;
2852  delete[] pMat;
2853  delete[] pL;
2854  delete[] coeffMonoms;
2855  delete[] pM;
2856  if (bufpEvalPoints != NULL)
2857  delete [] bufpEvalPoints;
2858  fail= true;
2859  if (alpha.level() != 1 && V_buf != alpha)
2860  prune1 (alpha);
2861  return 0;
2862  }
2863  }
2864  pEvalPoints[i + biggestSize]= evalPoints;
2865  gcds[i + biggestSize]= g;
2866  }
2867  }
2868  for (int i= 0; i < biggestSize; i++)
2869  {
2870  CFIterator l= gcds [i];
2871  evalPoints= pEvalPoints [i];
2872  for (int k= 1; k < skelSize; k++, l++)
2873  {
2874  if (i == 0)
2875  pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
2876  if (k == minimalColumnsIndex)
2877  continue;
2878  coeffEval= evaluate (coeffMonoms[k], evalPoints);
2879  if (pMat[k].rows() >= i + 1)
2880  {
2881  for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2882  pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2883  }
2884  }
2885  }
2886  Mat= bufMat0;
2887  pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
2888  for (int j= 1; j <= Mat.rows(); j++)
2889  for (int k= 1; k <= Mat.columns(); k++)
2890  pMat [0] (j,k)= Mat (j,k);
2891  Mat= bufMat;
2892  for (int j= 1; j <= Mat.rows(); j++)
2893  for (int k= 1; k <= Mat.columns(); k++)
2894  pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
2895  // write old matrix entries into new matrices
2896  for (int i= 0; i < skelSize; i++)
2897  {
2898  bufArray= pL[i];
2899  pL[i]= CFArray (biggestSize2);
2900  for (int j= 0; j < bufArray.size(); j++)
2901  pL[i] [j]= bufArray [j];
2902  }
2903  //write old vector entries into new and add new entries to old matrices
2904  for (int i= 0; i < biggestSize2 - biggestSize; i++)
2905  {
2906  CFIterator l= gcds [i + biggestSize];
2907  evalPoints= pEvalPoints [i + biggestSize];
2908  for (int k= 0; k < skelSize; k++, l++)
2909  {
2910  pL[k] [i + biggestSize]= l.coeff();
2911  coeffEval= evaluate (coeffMonoms[k], evalPoints);
2912  if (pMat[k].rows() >= i + biggestSize + 1)
2913  {
2914  for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2915  pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
2916  }
2917  }
2918  }
2919  // begin new
2920  for (int i= 0; i < skelSize; i++)
2921  {
2922  if (pL[i].size() > 1)
2923  {
2924  for (int j= 2; j <= pMat[i].rows(); j++)
2925  pMat[i] (j, coeffMonoms[i].size() + j - 1)=
2926  -pL[i] [j - 1];
2927  }
2928  }
2929 
2930  matColumns= biggestSize2 - 1;
2931  matRows= 0;
2932  for (int i= 0; i < skelSize; i++)
2933  {
2934  if (V_buf.level() == 1)
2935  (void) gaussianElimFp (pMat[i], pL[i]);
2936  else
2937  (void) gaussianElimFq (pMat[i], pL[i], V_buf);
2938 
2939  if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
2940  {
2941  delete[] pEvalPoints;
2942  delete[] pMat;
2943  delete[] pL;
2944  delete[] coeffMonoms;
2945  delete[] pM;
2946  if (bufpEvalPoints != NULL)
2947  delete [] bufpEvalPoints;
2948  fail= true;
2949  if (alpha.level() != 1 && V_buf != alpha)
2950  prune1 (alpha);
2951  return 0;
2952  }
2953  matRows += pMat[i].rows() - coeffMonoms[i].size();
2954  }
2955 
2956  CFMatrix bufMat;
2957  Mat= CFMatrix (matRows, matColumns);
2958  ind= 0;
2959  bufArray= CFArray (matRows);
2960  CFArray bufArray2;
2961  for (int i= 0; i < skelSize; i++)
2962  {
2963  if (coeffMonoms[i].size() + 1 >= pMat[i].rows() || coeffMonoms[i].size() + 1 >= pMat[i].columns())
2964  {
2965  delete[] pEvalPoints;
2966  delete[] pMat;
2967  delete[] pL;
2968  delete[] coeffMonoms;
2969  delete[] pM;
2970  if (bufpEvalPoints != NULL)
2971  delete [] bufpEvalPoints;
2972  fail= true;
2973  if (alpha.level() != 1 && V_buf != alpha)
2974  prune1 (alpha);
2975  return 0;
2976  }
2977  bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
2978  coeffMonoms[i].size() + 1, pMat[i].columns());
2979 
2980  for (int j= 1; j <= bufMat.rows(); j++)
2981  for (int k= 1; k <= bufMat.columns(); k++)
2982  Mat (j + ind, k)= bufMat(j, k);
2983  bufArray2= coeffMonoms[i].size();
2984  for (int j= 1; j <= pMat[i].rows(); j++)
2985  {
2986  if (j > coeffMonoms[i].size())
2987  bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
2988  else
2989  bufArray2 [j - 1]= pL[i] [j - 1];
2990  }
2991  pL[i]= bufArray2;
2992  ind += bufMat.rows();
2993  pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
2994  }
2995 
2996  if (V_buf.level() != 1)
2997  solution= solveSystemFq (Mat, bufArray, V_buf);
2998  else
2999  solution= solveSystemFp (Mat, bufArray);
3000 
3001  if (solution.size() == 0)
3002  {
3003  delete[] pEvalPoints;
3004  delete[] pMat;
3005  delete[] pL;
3006  delete[] coeffMonoms;
3007  delete[] pM;
3008  if (bufpEvalPoints != NULL)
3009  delete [] bufpEvalPoints;
3010  fail= true;
3011  if (alpha.level() != 1 && V_buf != alpha)
3012  prune1 (alpha);
3013  return 0;
3014  }
3015 
3016  ind= 0;
3017  result= 0;
3018  CFArray bufSolution;
3019  for (int i= 0; i < skelSize; i++)
3020  {
3021  bufSolution= readOffSolution (pMat[i], pL[i], solution);
3022  for (int i= 0; i < bufSolution.size(); i++, ind++)
3023  result += Monoms [ind]*bufSolution[i];
3024  }
3025  if (alpha.level() != 1 && V_buf != alpha)
3026  {
3027  CFList u, v;
3028  result= mapDown (result,prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3029  prune1 (alpha);
3030  }
3031  result= N(result);
3032  delete[] pEvalPoints;
3033  delete[] pMat;
3034  delete[] pL;
3035  delete[] coeffMonoms;
3036  delete[] pM;
3037 
3038  if (bufpEvalPoints != NULL)
3039  delete [] bufpEvalPoints;
3040  if (fdivides (result, F) && fdivides (result, G))
3041  return result;
3042  else
3043  {
3044  fail= true;
3045  return 0;
3046  }
3047  } // end of deKleine, Monagan & Wittkopf
3048 
3049  result += Monoms[0];
3050  int ind2= 0, ind3= 2;
3051  ind= 0;
3052  for (int l= 0; l < minimalColumnsIndex; l++)
3053  ind += coeffMonoms[l].size();
3054  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
3055  l++, ind2++, ind3++)
3056  {
3057  result += solution[l]*Monoms [1 + ind2];
3058  for (int i= 0; i < pMat[0].rows(); i++)
3059  firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
3060  }
3061  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
3062  result += solution[l]*Monoms [ind + l];
3063  ind= coeffMonoms[0].size();
3064  for (int k= 1; k < skelSize; k++)
3065  {
3066  if (k == minimalColumnsIndex)
3067  {
3068  ind += coeffMonoms[k].size();
3069  continue;
3070  }
3071  if (k != minimalColumnsIndex)
3072  {
3073  for (int i= 0; i < biggestSize; i++)
3074  pL[k] [i] *= firstColumn [i];
3075  }
3076 
3077  if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
3078  {
3079  bufArray= CFArray (coeffMonoms[k].size());
3080  for (int i= 0; i < bufArray.size(); i++)
3081  bufArray [i]= pL[k] [i];
3082  solution= solveGeneralVandermonde (pM[k], bufArray);
3083  }
3084  else
3085  solution= solveGeneralVandermonde (pM[k], pL[k]);
3086 
3087  if (solution.size() == 0)
3088  {
3089  delete[] pEvalPoints;
3090  delete[] pMat;
3091  delete[] pL;
3092  delete[] coeffMonoms;
3093  delete[] pM;
3094  fail= true;
3095  if (alpha.level() != 1 && V_buf != alpha)
3096  prune1 (alpha);
3097  return 0;
3098  }
3099  if (k != minimalColumnsIndex)
3100  {
3101  for (int l= 0; l < solution.size(); l++)
3102  result += solution[l]*Monoms [ind + l];
3103  ind += solution.size();
3104  }
3105  }
3106 
3107  delete[] pEvalPoints;
3108  delete[] pMat;
3109  delete[] pL;
3110  delete[] pM;
3111  delete[] coeffMonoms;
3112 
3113  if (alpha.level() != 1 && V_buf != alpha)
3114  {
3115  CFList u, v;
3116  result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3117  prune1 (alpha);
3118  }
3119  result= N(result);
3120 
3121  if (fdivides (result, F) && fdivides (result, G))
3122  return result;
3123  else
3124  {
3125  fail= true;
3126  return 0;
3127  }
3128 }
CFArray readOffSolution(const CFMatrix &M, const long rk)
M in row echolon form, rk rank of M.
Definition: cfModGcd.cc:1688
long gaussianElimFq(CFMatrix &M, CFArray &L, const Variable &alpha)
Definition: cfModGcd.cc:1784
long gaussianElimFp(CFMatrix &M, CFArray &L)
Definition: cfModGcd.cc:1739
CFArray solveSystemFp(const CFMatrix &M, const CFArray &L)
Definition: cfModGcd.cc:1840
CFArray solveSystemFq(const CFMatrix &M, const CFArray &L, const Variable &alpha)
Definition: cfModGcd.cc:1892
int rows() const
Definition: ftmpl_matrix.h:43
int columns() const
Definition: ftmpl_matrix.h:44
const signed long ceil(const ampf< Precision > &x)
Definition: amp.h:788
#define NULL
Definition: omList.c:12

◆ randomElement()

static CanonicalForm randomElement ( const CanonicalForm F,
const Variable alpha,
CFList list,
bool &  fail 
)
inlinestatic

compute a random element a of $ F_{p}(\alpha ) $ , s.t. F(a) $ \neq 0 $ , F is a univariate polynomial, returns fail if there are no field elements left which have not been used before

Definition at line 379 of file cfModGcd.cc.

381 {
382  fail= false;
383  Variable x= F.mvar();
384  AlgExtRandomF genAlgExt (alpha);
385  FFRandom genFF;
386  CanonicalForm random, mipo;
387  mipo= getMipo (alpha);
388  int p= getCharacteristic ();
389  int d= degree (mipo);
390  double bound= pow ((double) p, (double) d);
391  do
392  {
393  if (list.length() == bound)
394  {
395  fail= true;
396  break;
397  }
398  if (list.length() < p)
399  {
400  random= genFF.generate();
401  while (find (list, random))
402  random= genFF.generate();
403  }
404  else
405  {
406  random= genAlgExt.generate();
407  while (find (list, random))
408  random= genAlgExt.generate();
409  }
410  if (F (random, x) == 0)
411  {
412  list.append (random);
413  continue;
414  }
415  } while (find (list, random));
416  return random;
417 }

◆ readOffSolution() [1/2]

CFArray readOffSolution ( const CFMatrix M,
const CFArray L,
const CFArray partialSol 
)

Definition at line 1710 of file cfModGcd.cc.

1711 {
1712  CFArray result= CFArray (M.rows());
1713  CanonicalForm tmp1, tmp2, tmp3;
1714  int k;
1715  for (int i= M.rows(); i >= 1; i--)
1716  {
1717  tmp3= 0;
1718  tmp1= L[i - 1];
1719  k= 0;
1720  for (int j= M.columns(); j >= 1; j--, k++)
1721  {
1722  tmp2= M (i, j);
1723  if (j == i)
1724  break;
1725  else
1726  {
1727  if (k > partialSol.size() - 1)
1728  tmp3 += tmp2*result[j - 1];
1729  else
1730  tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
1731  }
1732  }
1733  result [i - 1]= (tmp1 - tmp3)/tmp2;
1734  }
1735  return result;
1736 }
CFList tmp1
Definition: facFqBivar.cc:72
CFList tmp2
Definition: facFqBivar.cc:72

◆ readOffSolution() [2/2]

CFArray readOffSolution ( const CFMatrix M,
const long  rk 
)

M in row echolon form, rk rank of M.

Definition at line 1688 of file cfModGcd.cc.

1689 {
1690  CFArray result= CFArray (rk);
1691  CanonicalForm tmp1, tmp2, tmp3;
1692  for (int i= rk; i >= 1; i--)
1693  {
1694  tmp3= 0;
1695  tmp1= M (i, M.columns());
1696  for (int j= M.columns() - 1; j >= 1; j--)
1697  {
1698  tmp2= M (i, j);
1699  if (j == i)
1700  break;
1701  else
1702  tmp3 += tmp2*result[j - 1];
1703  }
1704  result[i - 1]= (tmp1 - tmp3)/tmp2;
1705  }
1706  return result;
1707 }

◆ solveGeneralVandermonde()

CFArray solveGeneralVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1632 of file cfModGcd.cc.

1633 {
1634  int r= M.size();
1635  ASSERT (A.size() == r, "vector does not have right size");
1636  if (r == 1)
1637  {
1638  CFArray result= CFArray (1);
1639  result [0]= A[0] / M [0];
1640  return result;
1641  }
1642  // check solvability
1643  bool notDistinct= false;
1644  for (int i= 0; i < r - 1; i++)
1645  {
1646  for (int j= i + 1; j < r; j++)
1647  {
1648  if (M [i] == M [j])
1649  {
1650  notDistinct= true;
1651  break;
1652  }
1653  }
1654  }
1655  if (notDistinct)
1656  return CFArray();
1657 
1658  CanonicalForm master= 1;
1659  Variable x= Variable (1);
1660  for (int i= 0; i < r; i++)
1661  master *= x - M [i];
1662  master *= x;
1663  CFList Pj;
1664  CanonicalForm tmp;
1665  for (int i= 0; i < r; i++)
1666  {
1667  tmp= master/(x - M [i]);
1668  tmp /= tmp (M [i], 1);
1669  Pj.append (tmp);
1670  }
1671 
1672  CFArray result= CFArray (r);
1673 
1674  CFListIterator j= Pj;
1675  for (int i= 1; i <= r; i++, j++)
1676  {
1677  tmp= 0;
1678 
1679  for (int l= 1; l <= A.size(); l++)
1680  tmp += A[l - 1]*j.getItem()[l];
1681  result[i - 1]= tmp;
1682  }
1683  return result;
1684 }

◆ solveSystemFp()

CFArray solveSystemFp ( const CFMatrix M,
const CFArray L 
)

Definition at line 1840 of file cfModGcd.cc.

1841 {
1842  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1843  CFMatrix *N;
1844  N= new CFMatrix (M.rows(), M.columns() + 1);
1845 
1846  for (int i= 1; i <= M.rows(); i++)
1847  for (int j= 1; j <= M.columns(); j++)
1848  (*N) (i, j)= M (i, j);
1849 
1850  int j= 1;
1851  for (int i= 0; i < L.size(); i++, j++)
1852  (*N) (j, M.columns() + 1)= L[i];
1853 
1854 #ifdef HAVE_FLINT
1855  nmod_mat_t FLINTN;
1856  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1857  long rk= nmod_mat_rref (FLINTN);
1858 #else
1859  int p= getCharacteristic ();
1860  if (fac_NTL_char != p)
1861  {
1862  fac_NTL_char= p;
1863  zz_p::init (p);
1864  }
1865  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1866  long rk= gauss (*NTLN);
1867 #endif
1868  delete N;
1869  if (rk != M.columns())
1870  {
1871 #ifdef HAVE_FLINT
1872  nmod_mat_clear (FLINTN);
1873 #else
1874  delete NTLN;
1875 #endif
1876  return CFArray();
1877  }
1878 #ifdef HAVE_FLINT
1879  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1880  nmod_mat_clear (FLINTN);
1881 #else
1883  delete NTLN;
1884 #endif
1885  CFArray A= readOffSolution (*N, rk);
1886 
1887  delete N;
1888  return A;
1889 }

◆ solveSystemFq()

CFArray solveSystemFq ( const CFMatrix M,
const CFArray L,
const Variable alpha 
)

Definition at line 1892 of file cfModGcd.cc.

1893 {
1894  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1895  CFMatrix *N;
1896  N= new CFMatrix (M.rows(), M.columns() + 1);
1897 
1898  for (int i= 1; i <= M.rows(); i++)
1899  for (int j= 1; j <= M.columns(); j++)
1900  (*N) (i, j)= M (i, j);
1901  int j= 1;
1902  for (int i= 0; i < L.size(); i++, j++)
1903  (*N) (j, M.columns() + 1)= L[i];
1904  #ifdef HAVE_FLINT
1905  // convert mipo
1906  nmod_poly_t mipo1;
1908  fq_nmod_ctx_t ctx;
1909  fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
1910  nmod_poly_clear(mipo1);
1911  // convert matrix
1912  fq_nmod_mat_t FLINTN;
1913  convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
1914  // rank
1915  long rk= fq_nmod_mat_rref (FLINTN,ctx);
1916  #elif defined(HAVE_NTL)
1917  int p= getCharacteristic ();
1918  if (fac_NTL_char != p)
1919  {
1920  fac_NTL_char= p;
1921  zz_p::init (p);
1922  }
1923  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1924  zz_pE::init (NTLMipo);
1925  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1926  long rk= gauss (*NTLN);
1927  #else
1928  factoryError("NTL/FLINT missing: solveSystemFq");
1929  #endif
1930 
1931  delete N;
1932  if (rk != M.columns())
1933  {
1934  #if defined(HAVE_NTL) && !defined(HAVE_FLINT)
1935  delete NTLN;
1936  #endif
1937  return CFArray();
1938  }
1939  #ifdef HAVE_FLINT
1940  // convert and clean up
1942  fq_nmod_mat_clear (FLINTN,ctx);
1943  fq_nmod_ctx_clear(ctx);
1944  #elif defined(HAVE_NTL)
1946  delete NTLN;
1947  #endif
1948 
1949  CFArray A= readOffSolution (*N, rk);
1950 
1951  delete N;
1952  return A;
1953 }
CFMatrix * convertFq_nmod_mat_t2FacCFMatrix(const fq_nmod_mat_t m, const fq_nmod_ctx_t &fq_con, const Variable &alpha)
conversion of a FLINT matrix over F_q to a factory matrix

◆ solveVandermonde()

CFArray solveVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1579 of file cfModGcd.cc.

1580 {
1581  int r= M.size();
1582  ASSERT (A.size() == r, "vector does not have right size");
1583 
1584  if (r == 1)
1585  {
1586  CFArray result= CFArray (1);
1587  result [0]= A [0] / M [0];
1588  return result;
1589  }
1590  // check solvability
1591  bool notDistinct= false;
1592  for (int i= 0; i < r - 1; i++)
1593  {
1594  for (int j= i + 1; j < r; j++)
1595  {
1596  if (M [i] == M [j])
1597  {
1598  notDistinct= true;
1599  break;
1600  }
1601  }
1602  }
1603  if (notDistinct)
1604  return CFArray();
1605 
1606  CanonicalForm master= 1;
1607  Variable x= Variable (1);
1608  for (int i= 0; i < r; i++)
1609  master *= x - M [i];
1610  CFList Pj;
1611  CanonicalForm tmp;
1612  for (int i= 0; i < r; i++)
1613  {
1614  tmp= master/(x - M [i]);
1615  tmp /= tmp (M [i], 1);
1616  Pj.append (tmp);
1617  }
1618  CFArray result= CFArray (r);
1619 
1620  CFListIterator j= Pj;
1621  for (int i= 1; i <= r; i++, j++)
1622  {
1623  tmp= 0;
1624  for (int l= 0; l < A.size(); l++)
1625  tmp += A[l]*j.getItem()[l];
1626  result[i - 1]= tmp;
1627  }
1628  return result;
1629 }

◆ sparseGCDFp()

CanonicalForm sparseGCDFp ( const CanonicalForm F,
const CanonicalForm G,
bool &  topLevel,
CFList l 
)

Definition at line 3562 of file cfModGcd.cc.

3564 {
3565  CanonicalForm A= F;
3566  CanonicalForm B= G;
3567  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3568  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3569  else if (F.isZero() && G.isZero()) return F.genOne();
3570  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3571  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3572  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3573  if (F == G) return F/Lc(F);
3574 
3575  CFMap M,N;
3576  int best_level= myCompress (A, B, M, N, topLevel);
3577 
3578  if (best_level == 0) return B.genOne();
3579 
3580  A= M(A);
3581  B= M(B);
3582 
3583  Variable x= Variable (1);
3584 
3585  //univariate case
3586  if (A.isUnivariate() && B.isUnivariate())
3587  return N (gcd (A, B));
3588 
3589  CanonicalForm cA, cB; // content of A and B
3590  CanonicalForm ppA, ppB; // primitive part of A and B
3591  CanonicalForm gcdcAcB;
3592 
3593  cA = uni_content (A);
3594  cB = uni_content (B);
3595  gcdcAcB= gcd (cA, cB);
3596  ppA= A/cA;
3597  ppB= B/cB;
3598 
3599  CanonicalForm lcA, lcB; // leading coefficients of A and B
3600  CanonicalForm gcdlcAlcB;
3601  lcA= uni_lcoeff (ppA);
3602  lcB= uni_lcoeff (ppB);
3603 
3604  if (fdivides (lcA, lcB))
3605  {
3606  if (fdivides (A, B))
3607  return F/Lc(F);
3608  }
3609  if (fdivides (lcB, lcA))
3610  {
3611  if (fdivides (B, A))
3612  return G/Lc(G);
3613  }
3614 
3615  gcdlcAlcB= gcd (lcA, lcB);
3616 
3617  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3618  int d0;
3619 
3620  if (d == 0)
3621  return N(gcdcAcB);
3622 
3623  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3624 
3625  if (d0 < d)
3626  d= d0;
3627 
3628  if (d == 0)
3629  return N(gcdcAcB);
3630 
3631  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3632  CanonicalForm newtonPoly= 1;
3633  m= gcdlcAlcB;
3634  G_m= 0;
3635  H= 0;
3636  bool fail= false;
3637  topLevel= false;
3638  bool inextension= false;
3639  Variable V_buf, alpha, cleanUp;
3640  CanonicalForm prim_elem, im_prim_elem;
3641  CFList source, dest;
3642  do //first do
3643  {
3644  if (inextension)
3645  random_element= randomElement (m, V_buf, l, fail);
3646  else
3647  random_element= FpRandomElement (m, l, fail);
3648  if (random_element == 0 && !fail)
3649  {
3650  if (!find (l, random_element))
3651  l.append (random_element);
3652  continue;
3653  }
3654 
3655  if (!fail && !inextension)
3656  {
3657  CFList list;
3658  TIMING_START (gcd_recursion);
3659  G_random_element=
3660  sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3661  list);
3662  TIMING_END_AND_PRINT (gcd_recursion,
3663  "time for recursive call: ");
3664  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3665  }
3666  else if (!fail && inextension)
3667  {
3668  CFList list;
3669  TIMING_START (gcd_recursion);
3670  G_random_element=
3671  sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3672  list, topLevel);
3673  TIMING_END_AND_PRINT (gcd_recursion,
3674  "time for recursive call: ");
3675  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3676  }
3677  else if (fail && !inextension)
3678  {
3679  source= CFList();
3680  dest= CFList();
3681  CFList list;
3683  int deg= 2;
3684  bool initialized= false;
3685  do
3686  {
3687  mipo= randomIrredpoly (deg, x);
3688  if (initialized)
3689  setMipo (alpha, mipo);
3690  else
3691  alpha= rootOf (mipo);
3692  inextension= true;
3693  fail= false;
3694  initialized= true;
3695  random_element= randomElement (m, alpha, l, fail);
3696  deg++;
3697  } while (fail);
3698  cleanUp= alpha;
3699  V_buf= alpha;
3700  list= CFList();
3701  TIMING_START (gcd_recursion);
3702  G_random_element=
3703  sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3704  list, topLevel);
3705  TIMING_END_AND_PRINT (gcd_recursion,
3706  "time for recursive call: ");
3707  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3708  }
3709  else if (fail && inextension)
3710  {
3711  source= CFList();
3712  dest= CFList();
3713 
3714  Variable V_buf3= V_buf;
3715  V_buf= chooseExtension (V_buf);
3716  bool prim_fail= false;
3717  Variable V_buf2;
3718  CanonicalForm prim_elem, im_prim_elem;
3719  prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3720 
3721  if (V_buf3 != alpha)
3722  {
3723  m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3724  G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3725  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3726  dest);
3727  ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3728  ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3729  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3730  dest);
3731  for (CFListIterator i= l; i.hasItem(); i++)
3732  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3733  source, dest);
3734  }
3735 
3736  ASSERT (!prim_fail, "failure in integer factorizer");
3737  if (prim_fail)
3738  ; //ERROR
3739  else
3740  im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3741 
3742  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3743  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3744 
3745  for (CFListIterator i= l; i.hasItem(); i++)
3746  i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3747  im_prim_elem, source, dest);
3748  m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3749  G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3750  newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3751  source, dest);
3752  ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3753  ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3754  gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3755  source, dest);
3756  fail= false;
3757  random_element= randomElement (m, V_buf, l, fail );
3758  DEBOUTLN (cerr, "fail= " << fail);
3759  CFList list;
3760  TIMING_START (gcd_recursion);
3761  G_random_element=
3762  sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3763  list, topLevel);
3764  TIMING_END_AND_PRINT (gcd_recursion,
3765  "time for recursive call: ");
3766  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3767  alpha= V_buf;
3768  }
3769 
3770  if (!G_random_element.inCoeffDomain())
3771  d0= totaldegree (G_random_element, Variable(2),
3772  Variable (G_random_element.level()));
3773  else
3774  d0= 0;
3775 
3776  if (d0 == 0)
3777  {
3778  if (inextension)
3779  prune (cleanUp);
3780  return N(gcdcAcB);
3781  }
3782  if (d0 > d)
3783  {
3784  if (!find (l, random_element))
3785  l.append (random_element);
3786  continue;
3787  }
3788 
3789  G_random_element=
3790  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3791  * G_random_element;
3792 
3793  skeleton= G_random_element;
3794 
3795  if (!G_random_element.inCoeffDomain())
3796  d0= totaldegree (G_random_element, Variable(2),
3797  Variable (G_random_element.level()));
3798  else
3799  d0= 0;
3800 
3801  if (d0 < d)
3802  {
3803  m= gcdlcAlcB;
3804  newtonPoly= 1;
3805  G_m= 0;
3806  d= d0;
3807  }
3808 
3809  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3810 
3811  if (uni_lcoeff (H) == gcdlcAlcB)
3812  {
3813  cH= uni_content (H);
3814  ppH= H/cH;
3815  ppH /= Lc (ppH);
3816  DEBOUTLN (cerr, "ppH= " << ppH);
3817 
3818  if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3819  {
3820  if (inextension)
3821  prune (cleanUp);
3822  return N(gcdcAcB*ppH);
3823  }
3824  }
3825  G_m= H;
3826  newtonPoly= newtonPoly*(x - random_element);
3827  m= m*(x - random_element);
3828  if (!find (l, random_element))
3829  l.append (random_element);
3830 
3831  if ((getCharacteristic() > 3 && size (skeleton) < 200))
3832  {
3833  CFArray Monoms;
3834  CFArray* coeffMonoms;
3835 
3836  do //second do
3837  {
3838  if (inextension)
3839  random_element= randomElement (m, alpha, l, fail);
3840  else
3841  random_element= FpRandomElement (m, l, fail);
3842  if (random_element == 0 && !fail)
3843  {
3844  if (!find (l, random_element))
3845  l.append (random_element);
3846  continue;
3847  }
3848 
3849  bool sparseFail= false;
3850  if (!fail && !inextension)
3851  {
3852  CFList list;
3853  TIMING_START (gcd_recursion);
3854  if (LC (skeleton).inCoeffDomain())
3855  G_random_element=
3856  monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3857  skeleton, x, sparseFail, coeffMonoms,
3858  Monoms);
3859  else
3860  G_random_element=
3861  nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3862  skeleton, x, sparseFail,
3863  coeffMonoms, Monoms);
3864  TIMING_END_AND_PRINT (gcd_recursion,
3865  "time for recursive call: ");
3866  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3867  }
3868  else if (!fail && inextension)
3869  {
3870  CFList list;
3871  TIMING_START (gcd_recursion);
3872  if (LC (skeleton).inCoeffDomain())
3873  G_random_element=
3874  monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3875  skeleton, alpha, sparseFail, coeffMonoms,
3876  Monoms);
3877  else
3878  G_random_element=
3879  nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3880  skeleton, alpha, sparseFail, coeffMonoms,
3881  Monoms);
3882  TIMING_END_AND_PRINT (gcd_recursion,
3883  "time for recursive call: ");
3884  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3885  }
3886  else if (fail && !inextension)
3887  {
3888  source= CFList();
3889  dest= CFList();
3890  CFList list;
3892  int deg= 2;
3893  bool initialized= false;
3894  do
3895  {
3896  mipo= randomIrredpoly (deg, x);
3897  if (initialized)
3898  setMipo (alpha, mipo);
3899  else
3900  alpha= rootOf (mipo);
3901  inextension= true;
3902  fail= false;
3903  initialized= true;
3904  random_element= randomElement (m, alpha, l, fail);
3905  deg++;
3906  } while (fail);
3907  cleanUp= alpha;
3908  V_buf= alpha;
3909  list= CFList();
3910  TIMING_START (gcd_recursion);
3911  if (LC (skeleton).inCoeffDomain())
3912  G_random_element=
3913  monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3914  skeleton, alpha, sparseFail, coeffMonoms,
3915  Monoms);
3916  else
3917  G_random_element=
3918  nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3919  skeleton, alpha, sparseFail, coeffMonoms,
3920  Monoms);
3921  TIMING_END_AND_PRINT (gcd_recursion,
3922  "time for recursive call: ");
3923  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3924  }
3925  else if (fail && inextension)
3926  {
3927  source= CFList();
3928  dest= CFList();
3929 
3930  Variable V_buf3= V_buf;
3931  V_buf= chooseExtension (V_buf);
3932  bool prim_fail= false;
3933  Variable V_buf2;
3934  CanonicalForm prim_elem, im_prim_elem;
3935  prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3936 
3937  if (V_buf3 != alpha)
3938  {
3939  m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3940  G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3941  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3942  source, dest);
3943  ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3944  ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3945  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3946  source, dest);
3947  for (CFListIterator i= l; i.hasItem(); i++)
3948  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3949  source, dest);
3950  }
3951 
3952  ASSERT (!prim_fail, "failure in integer factorizer");
3953  if (prim_fail)
3954  ; //ERROR
3955  else
3956  im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3957 
3958  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3959  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3960 
3961  for (CFListIterator i= l; i.hasItem(); i++)
3962  i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3963  im_prim_elem, source, dest);
3964  m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3965  G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3966  newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3967  source, dest);
3968  ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3969  ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3970  gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3971  source, dest);
3972  fail= false;
3973  random_element= randomElement (m, V_buf, l, fail );
3974  DEBOUTLN (cerr, "fail= " << fail);
3975  CFList list;
3976  TIMING_START (gcd_recursion);
3977  if (LC (skeleton).inCoeffDomain())
3978  G_random_element=
3979  monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
3980  skeleton, V_buf, sparseFail, coeffMonoms,
3981  Monoms);
3982  else
3983  G_random_element=
3984  nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
3985  skeleton, V_buf, sparseFail,
3986  coeffMonoms, Monoms);
3987  TIMING_END_AND_PRINT (gcd_recursion,
3988  "time for recursive call: ");
3989  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3990  alpha= V_buf;
3991  }
3992 
3993  if (sparseFail)
3994  break;
3995 
3996  if (!G_random_element.inCoeffDomain())
3997  d0= totaldegree (G_random_element, Variable(2),
3998  Variable (G_random_element.level()));
3999  else
4000  d0= 0;
4001 
4002  if (d0 == 0)
4003  {
4004  if (inextension)
4005  prune (cleanUp);
4006  return N(gcdcAcB);
4007  }
4008  if (d0 > d)
4009  {
4010  if (!find (l, random_element))
4011  l.append (random_element);
4012  continue;
4013  }
4014 
4015  G_random_element=
4016  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
4017  * G_random_element;
4018 
4019  if (!G_random_element.inCoeffDomain())
4020  d0= totaldegree (G_random_element, Variable(2),
4021  Variable (G_random_element.level()));
4022  else
4023  d0= 0;
4024 
4025  if (d0 < d)
4026  {
4027  m= gcdlcAlcB;
4028  newtonPoly= 1;
4029  G_m= 0;
4030  d= d0;
4031  }
4032 
4033  TIMING_START (newton_interpolation);
4034  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
4035  TIMING_END_AND_PRINT (newton_interpolation,
4036  "time for newton interpolation: ");
4037 
4038  //termination test
4039  if (uni_lcoeff (H) == gcdlcAlcB)
4040  {
4041  cH= uni_content (H);
4042  ppH= H/cH;
4043  ppH /= Lc (ppH);
4044  DEBOUTLN (cerr, "ppH= " << ppH);
4045  if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4046  {
4047  if (inextension)
4048  prune (cleanUp);
4049  return N(gcdcAcB*ppH);
4050  }
4051  }
4052 
4053  G_m= H;
4054  newtonPoly= newtonPoly*(x - random_element);
4055  m= m*(x - random_element);
4056  if (!find (l, random_element))
4057  l.append (random_element);
4058 
4059  } while (1); //end of second do
4060  }
4061  else
4062  {
4063  if (inextension)
4064  prune (cleanUp);
4065  return N(gcdcAcB*modGCDFp (ppA, ppB));
4066  }
4067  } while (1); //end of first do
4068 }
CanonicalForm sparseGCDFq(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
Definition: cfModGcd.cc:3130
CanonicalForm monicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition: cfModGcd.cc:2199
CanonicalForm nonMonicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition: cfModGcd.cc:2483
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:3562

◆ sparseGCDFq()

CanonicalForm sparseGCDFq ( const CanonicalForm F,
const CanonicalForm G,
const Variable alpha,
CFList l,
bool &  topLevel 
)

Definition at line 3130 of file cfModGcd.cc.

3132 {
3133  CanonicalForm A= F;
3134  CanonicalForm B= G;
3135  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3136  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3137  else if (F.isZero() && G.isZero()) return F.genOne();
3138  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3139  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3140  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3141  if (F == G) return F/Lc(F);
3142 
3143  CFMap M,N;
3144  int best_level= myCompress (A, B, M, N, topLevel);
3145 
3146  if (best_level == 0) return B.genOne();
3147 
3148  A= M(A);
3149  B= M(B);
3150 
3151  Variable x= Variable (1);
3152 
3153  //univariate case
3154  if (A.isUnivariate() && B.isUnivariate())
3155  return N (gcd (A, B));
3156 
3157  CanonicalForm cA, cB; // content of A and B
3158  CanonicalForm ppA, ppB; // primitive part of A and B
3159  CanonicalForm gcdcAcB;
3160 
3161  cA = uni_content (A);
3162  cB = uni_content (B);
3163  gcdcAcB= gcd (cA, cB);
3164  ppA= A/cA;
3165  ppB= B/cB;
3166 
3167  CanonicalForm lcA, lcB; // leading coefficients of A and B
3168  CanonicalForm gcdlcAlcB;
3169  lcA= uni_lcoeff (ppA);
3170  lcB= uni_lcoeff (ppB);
3171 
3172  if (fdivides (lcA, lcB))
3173  {
3174  if (fdivides (A, B))
3175  return F/Lc(F);
3176  }
3177  if (fdivides (lcB, lcA))
3178  {
3179  if (fdivides (B, A))
3180  return G/Lc(G);
3181  }
3182 
3183  gcdlcAlcB= gcd (lcA, lcB);
3184 
3185  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3186  int d0;
3187 
3188  if (d == 0)
3189  return N(gcdcAcB);
3190  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3191 
3192  if (d0 < d)
3193  d= d0;
3194 
3195  if (d == 0)
3196  return N(gcdcAcB);
3197 
3198  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3199  CanonicalForm newtonPoly= 1;
3200  m= gcdlcAlcB;
3201  G_m= 0;
3202  H= 0;
3203  bool fail= false;
3204  topLevel= false;
3205  bool inextension= false;
3206  Variable V_buf= alpha, V_buf4= alpha;
3207  CanonicalForm prim_elem, im_prim_elem;
3208  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
3209  CFList source, dest;
3210  do // first do
3211  {
3212  random_element= randomElement (m, V_buf, l, fail);
3213  if (random_element == 0 && !fail)
3214  {
3215  if (!find (l, random_element))
3216  l.append (random_element);
3217  continue;
3218  }
3219  if (fail)
3220  {
3221  source= CFList();
3222  dest= CFList();
3223 
3224  Variable V_buf3= V_buf;
3225  V_buf= chooseExtension (V_buf);
3226  bool prim_fail= false;
3227  Variable V_buf2;
3228  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3229  if (V_buf4 == alpha)
3230  prim_elem_alpha= prim_elem;
3231 
3232  if (V_buf3 != V_buf4)
3233  {
3234  m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3235  G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3236  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3237  source, dest);
3238  ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3239  ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3240  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4, source,
3241  dest);
3242  for (CFListIterator i= l; i.hasItem(); i++)
3243  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3244  source, dest);
3245  }
3246 
3247  ASSERT (!prim_fail, "failure in integer factorizer");
3248  if (prim_fail)
3249  ; //ERROR
3250  else
3251  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3252 
3253  if (V_buf4 == alpha)
3254  im_prim_elem_alpha= im_prim_elem;
3255  else
3256  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
3257  im_prim_elem, source, dest);
3258 
3259  DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3260  DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3261  inextension= true;
3262  for (CFListIterator i= l; i.hasItem(); i++)
3263  i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3264  im_prim_elem, source, dest);
3265  m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3266  G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3267  newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3268  source, dest);
3269  ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3270  ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3271  gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3272  source, dest);
3273 
3274  fail= false;
3275  random_element= randomElement (m, V_buf, l, fail );
3276  DEBOUTLN (cerr, "fail= " << fail);
3277  CFList list;
3278  TIMING_START (gcd_recursion);
3279  G_random_element=
3280  sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3281  list, topLevel);
3282  TIMING_END_AND_PRINT (gcd_recursion,
3283  "time for recursive call: ");
3284  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3285  V_buf4= V_buf;
3286  }
3287  else
3288  {
3289  CFList list;
3290  TIMING_START (gcd_recursion);
3291  G_random_element=
3292  sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3293  list, topLevel);
3294  TIMING_END_AND_PRINT (gcd_recursion,
3295  "time for recursive call: ");
3296  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3297  }
3298 
3299  if (!G_random_element.inCoeffDomain())
3300  d0= totaldegree (G_random_element, Variable(2),
3301  Variable (G_random_element.level()));
3302  else
3303  d0= 0;
3304 
3305  if (d0 == 0)
3306  {
3307  if (inextension)
3308  prune1 (alpha);
3309  return N(gcdcAcB);
3310  }
3311  if (d0 > d)
3312  {
3313  if (!find (l, random_element))
3314  l.append (random_element);
3315  continue;
3316  }
3317 
3318  G_random_element=
3319  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3320  * G_random_element;
3321 
3322  skeleton= G_random_element;
3323  if (!G_random_element.inCoeffDomain())
3324  d0= totaldegree (G_random_element, Variable(2),
3325  Variable (G_random_element.level()));
3326  else
3327  d0= 0;
3328 
3329  if (d0 < d)
3330  {
3331  m= gcdlcAlcB;
3332  newtonPoly= 1;
3333  G_m= 0;
3334  d= d0;
3335  }
3336 
3337  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3338  if (uni_lcoeff (H) == gcdlcAlcB)
3339  {
3340  cH= uni_content (H);
3341  ppH= H/cH;
3342  if (inextension)
3343  {
3344  CFList u, v;
3345  //maybe it's better to test if ppH is an element of F(\alpha) before
3346  //mapping down
3347  if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3348  {
3349  DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3350  ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3351  ppH /= Lc(ppH);
3352  DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3353  prune1 (alpha);
3354  return N(gcdcAcB*ppH);
3355  }
3356  }
3357  else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3358  return N(gcdcAcB*ppH);
3359  }
3360  G_m= H;
3361  newtonPoly= newtonPoly*(x - random_element);
3362  m= m*(x - random_element);
3363  if (!find (l, random_element))
3364  l.append (random_element);
3365 
3366  if (getCharacteristic () > 3 && size (skeleton) < 100)
3367  {
3368  CFArray Monoms;
3369  CFArray *coeffMonoms;
3370  do //second do
3371  {
3372  random_element= randomElement (m, V_buf, l, fail);
3373  if (random_element == 0 && !fail)
3374  {
3375  if (!find (l, random_element))
3376  l.append (random_element);
3377  continue;
3378  }
3379  if (fail)
3380  {
3381  source= CFList();
3382  dest= CFList();
3383 
3384  Variable V_buf3= V_buf;
3385  V_buf= chooseExtension (V_buf);
3386  bool prim_fail= false;
3387  Variable V_buf2;
3388  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3389  if (V_buf4 == alpha)
3390  prim_elem_alpha= prim_elem;
3391 
3392  if (V_buf3 != V_buf4)
3393  {
3394  m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3395  G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3396  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3397  source, dest);
3398  ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3399  ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3400  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
3401  source, dest);
3402  for (CFListIterator i= l; i.hasItem(); i++)
3403  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3404  source, dest);
3405  }
3406 
3407  ASSERT (!prim_fail, "failure in integer factorizer");
3408  if (prim_fail)
3409  ; //ERROR
3410  else
3411  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3412 
3413  if (V_buf4 == alpha)
3414  im_prim_elem_alpha= im_prim_elem;
3415  else
3416  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
3417  prim_elem, im_prim_elem, source, dest);
3418 
3419  DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3420  DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3421  inextension= true;
3422  for (CFListIterator i= l; i.hasItem(); i++)
3423  i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3424  im_prim_elem, source, dest);
3425  m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3426  G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3427  newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3428  source, dest);
3429  ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3430  ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3431 
3432  gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3433  source, dest);
3434 
3435  fail= false;
3436  random_element= randomElement (m, V_buf, l, fail);
3437  DEBOUTLN (cerr, "fail= " << fail);
3438  CFList list;
3439  TIMING_START (gcd_recursion);
3440 
3441  V_buf4= V_buf;
3442 
3443  //sparseInterpolation
3444  bool sparseFail= false;
3445  if (LC (skeleton).inCoeffDomain())
3446  G_random_element=
3447  monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3448  skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3449  else
3450  G_random_element=
3451  nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3452  skeleton, V_buf, sparseFail, coeffMonoms,
3453  Monoms);
3454  if (sparseFail)
3455  break;
3456 
3457  TIMING_END_AND_PRINT (gcd_recursion,
3458  "time for recursive call: ");
3459  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3460  }
3461  else
3462  {
3463  CFList list;
3464  TIMING_START (gcd_recursion);
3465  bool sparseFail= false;
3466  if (LC (skeleton).inCoeffDomain())
3467  G_random_element=
3468  monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3469  skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3470  else
3471  G_random_element=
3472  nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3473  skeleton, V_buf, sparseFail, coeffMonoms,
3474  Monoms);
3475  if (sparseFail)
3476  break;
3477 
3478  TIMING_END_AND_PRINT (gcd_recursion,
3479  "time for recursive call: ");
3480  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3481  }
3482 
3483  if (!G_random_element.inCoeffDomain())
3484  d0= totaldegree (G_random_element, Variable(2),
3485  Variable (G_random_element.level()));
3486  else
3487  d0= 0;
3488 
3489  if (d0 == 0)
3490  {
3491  if (inextension)
3492  prune1 (alpha);
3493  return N(gcdcAcB);
3494  }
3495  if (d0 > d)
3496  {
3497  if (!find (l, random_element))
3498  l.append (random_element);
3499  continue;
3500  }
3501 
3502  G_random_element=
3503  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3504  * G_random_element;
3505 
3506  if (!G_random_element.inCoeffDomain())
3507  d0= totaldegree (G_random_element, Variable(2),
3508  Variable (G_random_element.level()));
3509  else
3510  d0= 0;
3511 
3512  if (d0 < d)
3513  {
3514  m= gcdlcAlcB;
3515  newtonPoly= 1;
3516  G_m= 0;
3517  d= d0;
3518  }
3519 
3520  TIMING_START (newton_interpolation);
3521  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3522  TIMING_END_AND_PRINT (newton_interpolation,
3523  "time for newton interpolation: ");
3524 
3525  //termination test
3526  if (uni_lcoeff (H) == gcdlcAlcB)
3527  {
3528  cH= uni_content (H);
3529  ppH= H/cH;
3530  if (inextension)
3531  {
3532  CFList u, v;
3533  //maybe it's better to test if ppH is an element of F(\alpha) before
3534  //mapping down
3535  if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3536  {
3537  DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3538  ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3539  ppH /= Lc(ppH);
3540  DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3541  prune1 (alpha);
3542  return N(gcdcAcB*ppH);
3543  }
3544  }
3545  else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3546  {
3547  return N(gcdcAcB*ppH);
3548  }
3549  }
3550 
3551  G_m= H;
3552  newtonPoly= newtonPoly*(x - random_element);
3553  m= m*(x - random_element);
3554  if (!find (l, random_element))
3555  l.append (random_element);
3556 
3557  } while (1);
3558  }
3559  } while (1);
3560 }

◆ TIMING_DEFINE_PRINT() [1/2]

TIMING_DEFINE_PRINT ( gcd_recursion  ) const &

◆ TIMING_DEFINE_PRINT() [2/2]

TIMING_DEFINE_PRINT ( modZ_termination  ) const &

modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer Algebra", Algorithm 7.1

◆ uni_content() [1/2]

static CanonicalForm uni_content ( const CanonicalForm F)
inlinestatic

compute the content of F, where F is considered as an element of $ R[x_{1}][x_{2},\ldots ,x_{n}] $

Definition at line 281 of file cfModGcd.cc.

282 {
283  if (F.inBaseDomain())
284  return F.genOne();
285  if (F.level() == 1 && F.isUnivariate())
286  return F;
287  if (F.level() != 1 && F.isUnivariate())
288  return F.genOne();
289  if (degree (F,1) == 0) return F.genOne();
290 
291  int l= F.level();
292  if (l == 2)
293  return content(F);
294  else
295  {
296  CanonicalForm pol, c = 0;
297  CFIterator i = F;
298  for (; i.hasTerms(); i++)
299  {
300  pol= i.coeff();
301  pol= uni_content (pol);
302  c= gcd (c, pol);
303  if (c.isOne())
304  return c;
305  }
306  return c;
307  }
308 }
CanonicalForm FACTORY_PUBLIC content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:603

◆ uni_content() [2/2]

CanonicalForm uni_content ( const CanonicalForm F,
const Variable x 
)

Definition at line 259 of file cfModGcd.cc.

260 {
261  if (F.inCoeffDomain())
262  return F.genOne();
263  if (F.level() == x.level() && F.isUnivariate())
264  return F;
265  if (F.level() != x.level() && F.isUnivariate())
266  return F.genOne();
267 
268  if (x.level() != 1)
269  {
270  CanonicalForm f= swapvar (F, x, Variable (1));
272  return swapvar (result, x, Variable (1));
273  }
274  else
275  return uni_content (F);
276 }
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168

◆ uni_lcoeff()

static CanonicalForm uni_lcoeff ( const CanonicalForm F)
inlinestatic

compute the leading coefficient of F, where F is considered as an element of $ R[x_{1}][x_{2},\ldots ,x_{n}] $, order on $ Mon (x_{2},\ldots ,x_{n}) $ is dp.

Definition at line 339 of file cfModGcd.cc.

340 {
341  if (F.level() > 1)
342  {
343  Variable x= Variable (2);
344  int deg= totaldegree (F, x, F.mvar());
345  for (CFIterator i= F; i.hasTerms(); i++)
346  {
347  if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
348  return uni_lcoeff (i.coeff());
349  }
350  }
351  return F;
352 }

◆ while()

while ( true  )

Definition at line 4132 of file cfModGcd.cc.

4133  {
4134  p = cf_getBigPrime( i );
4135  i--;
4136  while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
4137  {
4138  p = cf_getBigPrime( i );
4139  i--;
4140  }
4141  //printf("try p=%d\n",p);
4142  setCharacteristic( p );
4143  fp= mapinto (f);
4144  gp= mapinto (g);
4145  TIMING_START (modZ_recursion)
4146 #if defined(HAVE_NTL) || defined(HAVE_FLINT)
4147  if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
4148  Dp = modGCDFp (fp, gp, cofp, cogp);
4149  else
4150  {
4151  Dp= gcd_poly (fp, gp);
4152  cofp= fp/Dp;
4153  cogp= gp/Dp;
4154  }
4155 #else
4156  Dp= gcd_poly (fp, gp);
4157  cofp= fp/Dp;
4158  cogp= gp/Dp;
4159 #endif
4160  TIMING_END_AND_PRINT (modZ_recursion,
4161  "time for gcd mod p in modular gcd: ");
4162  Dp /=Dp.lc();
4163  Dp *= mapinto (cl);
4164  cofp /= lc (cofp);
4165  cofp *= mapinto (lcf);
4166  cogp /= lc (cogp);
4167  cogp *= mapinto (lcg);
4168  setCharacteristic( 0 );
4169  dp_deg=totaldegree(Dp);
4170  if ( dp_deg == 0 )
4171  {
4172  //printf(" -> 1\n");
4173  return CanonicalForm(gcdcfcg);
4174  }
4175  if ( q.isZero() )
4176  {
4177  D = mapinto( Dp );
4178  cof= mapinto (cofp);
4179  cog= mapinto (cogp);
4180  d_deg=dp_deg;
4181  q = p;
4182  Dn= balance_p (D, p);
4183  cofn= balance_p (cof, p);
4184  cogn= balance_p (cog, p);
4185  }
4186  else
4187  {
4188  if ( dp_deg == d_deg )
4189  {
4190  chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
4191  chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
4192  chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
4193  cof= newCof;
4194  cog= newCog;
4195  newqh= newq/2;
4196  Dn= balance_p (newD, newq, newqh);
4197  cofn= balance_p (newCof, newq, newqh);
4198  cogn= balance_p (newCog, newq, newqh);
4199  if (test != Dn) //balance_p (newD, newq))
4200  test= balance_p (newD, newq);
4201  else
4202  equal= true;
4203  q = newq;
4204  D = newD;
4205  }
4206  else if ( dp_deg < d_deg )
4207  {
4208  //printf(" were all bad, try more\n");
4209  // all previous p's are bad primes
4210  q = p;
4211  D = mapinto( Dp );
4212  cof= mapinto (cof);
4213  cog= mapinto (cog);
4214  d_deg=dp_deg;
4215  test= 0;
4216  equal= false;
4217  Dn= balance_p (D, p);
4218  cofn= balance_p (cof, p);
4219  cogn= balance_p (cog, p);
4220  }
4221  else
4222  {
4223  //printf(" was bad, try more\n");
4224  }
4225  //else dp_deg > d_deg: bad prime
4226  }
4227  if ( i >= 0 )
4228  {
4229  cDn= icontent (Dn);
4230  Dn /= cDn;
4231  cofn /= cl/cDn;
4232  //cofn /= icontent (cofn);
4233  cogn /= cl/cDn;
4234  //cogn /= icontent (cogn);
4235  TIMING_START (modZ_termination);
4236  if ((terminationTest (f,g, cofn, cogn, Dn)) ||
4237  ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
4238  {
4239  TIMING_END_AND_PRINT (modZ_termination,
4240  "time for successful termination in modular gcd: ");
4241  //printf(" -> success\n");
4242  return Dn*gcdcfcg;
4243  }
4244  TIMING_END_AND_PRINT (modZ_termination,
4245  "time for unsuccessful termination in modular gcd: ");
4246  equal= false;
4247  //else: try more primes
4248  }
4249  else
4250  { // try other method
4251  //printf("try other gcd\n");
4253  D=gcd_poly( f, g );
4255  return D*gcdcfcg;
4256  }
4257  }
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm mapinto(const CanonicalForm &f)
CanonicalForm lc(const CanonicalForm &f)
CanonicalForm FACTORY_PUBLIC icontent(const CanonicalForm &f)
CanonicalForm icontent ( const CanonicalForm & f )
Definition: cf_gcd.cc:74
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
CanonicalForm FACTORY_PUBLIC gcd_poly(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
Definition: cf_gcd.cc:492
CanonicalForm cofp
Definition: cfModGcd.cc:4129
CanonicalForm lcg
Definition: cfModGcd.cc:4097
int dp_deg
Definition: cfModGcd.cc:4078
CanonicalForm newCog
Definition: cfModGcd.cc:4129
CanonicalForm cogn
Definition: cfModGcd.cc:4129
cl
Definition: cfModGcd.cc:4100
CanonicalForm lcf
Definition: cfModGcd.cc:4097
int maxNumVars
Definition: cfModGcd.cc:4130
CanonicalForm fp
Definition: cfModGcd.cc:4102
CanonicalForm cogp
Definition: cfModGcd.cc:4129
CanonicalForm cofn
Definition: cfModGcd.cc:4129
CanonicalForm cof
Definition: cfModGcd.cc:4129
CanonicalForm cog
Definition: cfModGcd.cc:4129
CanonicalForm gcdcfcg
Definition: cfModGcd.cc:4101
CanonicalForm Dn
Definition: cfModGcd.cc:4096
CanonicalForm newCof
Definition: cfModGcd.cc:4129
CanonicalForm gp
Definition: cfModGcd.cc:4102
bool equal
Definition: cfModGcd.cc:4126
CanonicalForm test
Definition: cfModGcd.cc:4096
CanonicalForm b
Definition: cfModGcd.cc:4103
CanonicalForm cDn
Definition: cfModGcd.cc:4129
int d_deg
Definition: cfModGcd.cc:4078
void FACTORY_PUBLIC chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition: cf_chinese.cc:57
static const int SW_USE_CHINREM_GCD
set to 1 to use modular gcd over Z
Definition: cf_defs.h:41
static CanonicalForm balance_p(const CanonicalForm &f, const CanonicalForm &q, const CanonicalForm &qh)
Definition: cf_gcd.cc:149
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
#define D(A)
Definition: gentable.cc:131

Variable Documentation

◆ b

else L b = 1

Definition at line 4103 of file cfModGcd.cc.

◆ cand

Initial value:
{
CanonicalForm LCCand= abs (LC (cand))

Definition at line 68 of file cfModGcd.cc.

◆ cd

cd ( bCommonDen(FF)  ) = bCommonDen( GG )

Definition at line 4089 of file cfModGcd.cc.

◆ cDn

Definition at line 4129 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4083 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4083 of file cfModGcd.cc.

◆ cl

cl = gcd (f.lc(),g.lc())

Definition at line 4100 of file cfModGcd.cc.

◆ coF

Definition at line 67 of file cfModGcd.cc.

◆ cof

Definition at line 4129 of file cfModGcd.cc.

◆ cofn

Definition at line 4129 of file cfModGcd.cc.

◆ cofp

Definition at line 4129 of file cfModGcd.cc.

◆ coG

Definition at line 67 of file cfModGcd.cc.

◆ cog

Definition at line 4129 of file cfModGcd.cc.

◆ cogn

Definition at line 4129 of file cfModGcd.cc.

◆ cogp

Definition at line 4129 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4078 of file cfModGcd.cc.

◆ Dn

Definition at line 4096 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4078 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4126 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4081 of file cfModGcd.cc.

◆ false

return false

Definition at line 84 of file cfModGcd.cc.

◆ fp

Definition at line 4102 of file cfModGcd.cc.

◆ G

Definition at line 66 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4090 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4101 of file cfModGcd.cc.

◆ GG

Initial value:
{
CanonicalForm f, g, cl, q(0), Dp, newD, D, newq, newqh

Definition at line 4075 of file cfModGcd.cc.

◆ gp

Definition at line 4102 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4078 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4097 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4097 of file cfModGcd.cc.

◆ maxNumVars

int maxNumVars = tmax (getNumVars (f), getNumVars (g))

Definition at line 4130 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4104 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4129 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4129 of file cfModGcd.cc.

◆ p

int p

Definition at line 4078 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4096 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4082 of file cfModGcd.cc.