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
 
static const double log2exp = 1.442695041
 
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 422 of file cfModGcd.cc.

423 {
424  int i, m;
425  // extension of F_p needed
426  if (alpha.level() == 1)
427  {
428  i= 1;
429  m= 2;
430  } //extension of F_p(alpha)
431  if (alpha.level() != 1)
432  {
433  i= 4;
434  m= degree (getMipo (alpha));
435  }
436  #ifdef HAVE_FLINT
437  nmod_poly_t Irredpoly;
438  nmod_poly_init(Irredpoly,getCharacteristic());
439  nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, i*m+1);
440  CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
441  nmod_poly_clear(Irredpoly);
442  #else
444  {
447  }
448  zz_pX NTLIrredpoly;
449  BuildIrred (NTLIrredpoly, i*m);
450  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
451  #endif
452  return rootOf (newMipo);
453 }
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:4080
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:134
int level() const
Definition: factory.h:150
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 2187 of file cfModGcd.cc.

2189 {
2190  Aeval= A;
2191  Beval= B;
2192  int j= 1;
2193  for (CFListIterator i= L; i.hasItem(); i++, j++)
2194  {
2195  Aeval= Aeval (i.getItem(), j);
2196  Beval= Beval (i.getItem(), j);
2197  }
2198 }
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 2033 of file cfModGcd.cc.

2034 {
2035  CFArray result= A.size();
2036  CanonicalForm tmp;
2037  int k;
2038  for (int i= 0; i < A.size(); i++)
2039  {
2040  tmp= A[i];
2041  k= 1;
2042  for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2043  tmp= tmp (j.getItem(), k);
2044  result[i]= tmp;
2045  }
2046  return result;
2047 }
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 1994 of file cfModGcd.cc.

1995 {
1996  if (F.inCoeffDomain())
1997  {
1998  CFArray result= CFArray (1);
1999  result [0]= F;
2000  return result;
2001  }
2002  if (F.isUnivariate())
2003  {
2004  ASSERT (evalPoints.length() == 1,
2005  "expected an eval point with only one component");
2006  CFArray result= CFArray (size(F));
2007  int j= 0;
2009  for (CFIterator i= F; i.hasTerms(); i++, j++)
2010  result[j]= power (evalPoint, i.exp());
2011  return result;
2012  }
2013  int numMon= size (F);
2014  CFArray result= CFArray (numMon);
2015  int j= 0;
2018  buf.removeLast();
2019  CFArray recResult;
2020  CanonicalForm powEvalPoint;
2021  for (CFIterator i= F; i.hasTerms(); i++)
2022  {
2023  powEvalPoint= power (evalPoint, i.exp());
2024  recResult= evaluateMonom (i.coeff(), buf);
2025  for (int k= 0; k < recResult.size(); k++)
2026  result[j+k]= powEvalPoint*recResult[k];
2027  j += recResult.size();
2028  }
2029  return result;
2030 }
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:1994
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 2050 of file cfModGcd.cc.

2055 {
2056  int k= tmax (F.level(), G.level()) - 1;
2057  Variable x= Variable (1);
2058  CFList result;
2059  FFRandom genFF;
2060  GFRandom genGF;
2061  int p= getCharacteristic ();
2062  double bound;
2063  if (alpha != Variable (1))
2064  {
2065  bound= pow ((double) p, (double) degree (getMipo(alpha)));
2066  bound= pow (bound, (double) k);
2067  }
2068  else if (GF)
2069  {
2070  bound= pow ((double) p, (double) getGFDegree());
2071  bound= pow ((double) bound, (double) k);
2072  }
2073  else
2074  bound= pow ((double) p, (double) k);
2075 
2076  CanonicalForm random;
2077  int j;
2078  bool zeroOneOccured= false;
2079  bool allEqual= false;
2081  do
2082  {
2083  random= 0;
2084  // possible overflow if list.length() does not fit into a int
2085  if (list.length() >= bound)
2086  {
2087  fail= true;
2088  break;
2089  }
2090  for (int i= 0; i < k; i++)
2091  {
2092  if (GF)
2093  {
2094  result.append (genGF.generate());
2095  random += result.getLast()*power (x, i);
2096  }
2097  else if (alpha.level() != 1)
2098  {
2099  AlgExtRandomF genAlgExt (alpha);
2100  result.append (genAlgExt.generate());
2101  random += result.getLast()*power (x, i);
2102  }
2103  else
2104  {
2105  result.append (genFF.generate());
2106  random += result.getLast()*power (x, i);
2107  }
2108  if (result.getLast().isOne() || result.getLast().isZero())
2109  zeroOneOccured= true;
2110  }
2111  if (find (list, random))
2112  {
2113  zeroOneOccured= false;
2114  allEqual= false;
2115  result= CFList();
2116  continue;
2117  }
2118  if (zeroOneOccured)
2119  {
2120  list.append (random);
2121  zeroOneOccured= false;
2122  allEqual= false;
2123  result= CFList();
2124  continue;
2125  }
2126  // no zero at this point
2127  if (k > 1)
2128  {
2129  allEqual= true;
2130  CFIterator iter= random;
2131  buf= iter.coeff();
2132  iter++;
2133  for (; iter.hasTerms(); iter++)
2134  if (buf != iter.coeff())
2135  allEqual= false;
2136  }
2137  if (allEqual)
2138  {
2139  list.append (random);
2140  allEqual= false;
2141  zeroOneOccured= false;
2142  result= CFList();
2143  continue;
2144  }
2145 
2146  Feval= F;
2147  Geval= G;
2148  CanonicalForm LCeval= LCF;
2149  j= 1;
2150  for (CFListIterator i= result; i.hasItem(); i++, j++)
2151  {
2152  Feval= Feval (i.getItem(), j);
2153  Geval= Geval (i.getItem(), j);
2154  LCeval= LCeval (i.getItem(), j);
2155  }
2156 
2157  if (LCeval.isZero())
2158  {
2159  if (!find (list, random))
2160  list.append (random);
2161  zeroOneOccured= false;
2162  allEqual= false;
2163  result= CFList();
2164  continue;
2165  }
2166 
2167  if (list.length() >= bound)
2168  {
2169  fail= true;
2170  break;
2171  }
2172  } while (find (list, random));
2173 
2174  return result;
2175 }
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:4084
int p
Definition: cfModGcd.cc:4080
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 313 of file cfModGcd.cc.

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

◆ for() [1/2]

for ( i,
0;i--   
)

Definition at line 4119 of file cfModGcd.cc.

4120  {
4121  if (degree (f, i) <= 0 || degree (g, i) <= 0)
4122  continue;
4123  else
4125  }
f
Definition: cfModGcd.cc:4083
g
Definition: cfModGcd.cc:4092
int minCommonDeg
Definition: cfModGcd.cc:4106
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)

◆ for() [2/2]

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

Definition at line 4107 of file cfModGcd.cc.

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

◆ FpRandomElement()

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

Definition at line 1173 of file cfModGcd.cc.

1174 {
1175  fail= false;
1176  Variable x= F.mvar();
1177  FFRandom genFF;
1178  CanonicalForm random;
1179  int p= getCharacteristic();
1180  int bound= p;
1181  do
1182  {
1183  if (list.length() == bound)
1184  {
1185  fail= true;
1186  break;
1187  }
1188  if (list.length() < 1)
1189  random= 0;
1190  else
1191  {
1192  random= genFF.generate();
1193  while (find (list, random))
1194  random= genFF.generate();
1195  }
1196  if (F (random, x) == 0)
1197  {
1198  list.append (random);
1199  continue;
1200  }
1201  } while (find (list, random));
1202  return random;
1203 }
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 1741 of file cfModGcd.cc.

1742 {
1743  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1744  CFMatrix *N;
1745  N= new CFMatrix (M.rows(), M.columns() + 1);
1746 
1747  for (int i= 1; i <= M.rows(); i++)
1748  for (int j= 1; j <= M.columns(); j++)
1749  (*N) (i, j)= M (i, j);
1750 
1751  int j= 1;
1752  for (int i= 0; i < L.size(); i++, j++)
1753  (*N) (j, M.columns() + 1)= L[i];
1754 #ifdef HAVE_FLINT
1755  nmod_mat_t FLINTN;
1756  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1757  long rk= nmod_mat_rref (FLINTN);
1758 
1759  delete N;
1760  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1761  nmod_mat_clear (FLINTN);
1762 #else
1763  int p= getCharacteristic ();
1764  if (fac_NTL_char != p)
1765  {
1766  fac_NTL_char= p;
1767  zz_p::init (p);
1768  }
1769  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1770  delete N;
1771  long rk= gauss (*NTLN);
1772 
1774  delete NTLN;
1775 #endif
1776 
1777  L= CFArray (M.rows());
1778  for (int i= 0; i < M.rows(); i++)
1779  L[i]= (*N) (i + 1, M.columns() + 1);
1780  M= (*N) (1, M.rows(), 1, M.columns());
1781  delete N;
1782  return rk;
1783 }
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 1786 of file cfModGcd.cc.

1787 {
1788  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1789  CFMatrix *N;
1790  N= new CFMatrix (M.rows(), M.columns() + 1);
1791 
1792  for (int i= 1; i <= M.rows(); i++)
1793  for (int j= 1; j <= M.columns(); j++)
1794  (*N) (i, j)= M (i, j);
1795 
1796  int j= 1;
1797  for (int i= 0; i < L.size(); i++, j++)
1798  (*N) (j, M.columns() + 1)= L[i];
1799  #ifdef HAVE_FLINT
1800  // convert mipo
1801  nmod_poly_t mipo1;
1803  fq_nmod_ctx_t ctx;
1804  fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
1805  nmod_poly_clear(mipo1);
1806  // convert matrix
1807  fq_nmod_mat_t FLINTN;
1808  convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
1809  // rank
1810  long rk= fq_nmod_mat_rref (FLINTN,ctx);
1811  // clean up
1812  fq_nmod_mat_clear (FLINTN,ctx);
1813  fq_nmod_ctx_clear(ctx);
1814  #elif defined(HAVE_NTL)
1815  int p= getCharacteristic ();
1816  if (fac_NTL_char != p)
1817  {
1818  fac_NTL_char= p;
1819  zz_p::init (p);
1820  }
1821  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1822  zz_pE::init (NTLMipo);
1823  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1824  long rk= gauss (*NTLN);
1826  delete NTLN;
1827  #else
1828  factoryError("NTL/FLINT missing: gaussianElimFq");
1829  #endif
1830  delete N;
1831 
1832  M= (*N) (1, M.rows(), 1, M.columns());
1833  L= CFArray (M.rows());
1834  for (int i= 0; i < M.rows(); i++)
1835  L[i]= (*N) (i + 1, M.columns() + 1);
1836 
1837  delete N;
1838  return rk;
1839 }
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 1959 of file cfModGcd.cc.

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

◆ 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 821 of file cfModGcd.cc.

822 {
823  fail= false;
824  Variable x= F.mvar();
825  GFRandom genGF;
826  CanonicalForm random;
827  int p= getCharacteristic();
828  int d= getGFDegree();
829  int bound= ipower (p, d);
830  do
831  {
832  if (list.length() == bound)
833  {
834  fail= true;
835  break;
836  }
837  if (list.length() < 1)
838  random= 0;
839  else
840  {
841  random= genGF.generate();
842  while (find (list, random))
843  random= genGF.generate();
844  }
845  if (F (random, x) == 0)
846  {
847  list.append (random);
848  continue;
849  }
850  } while (find (list, random));
851  return random;
852 }
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 1214 of file cfModGcd.cc.

1216 {
1217  CanonicalForm dummy1, dummy2;
1218  CanonicalForm result= modGCDFp (F, G, dummy1, dummy2, topLevel, l);
1219  return result;
1220 }
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:1225

◆ modGCDFp() [2/2]

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

Definition at line 1225 of file cfModGcd.cc.

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

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

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

◆ 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 874 of file cfModGcd.cc.

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

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

◆ monicSparseInterpol()

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

Definition at line 2201 of file cfModGcd.cc.

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

◆ mult()

void mult ( CFList L1,
const CFList L2 
)

multiply two lists componentwise

Definition at line 2178 of file cfModGcd.cc.

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

◆ 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 93 of file cfModGcd.cc.

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

◆ 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 366 of file cfModGcd.cc.

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

◆ nonMonicSparseInterpol()

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

Definition at line 2485 of file cfModGcd.cc.

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

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

◆ readOffSolution() [1/2]

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

Definition at line 1712 of file cfModGcd.cc.

1713 {
1714  CFArray result= CFArray (M.rows());
1715  CanonicalForm tmp1, tmp2, tmp3;
1716  int k;
1717  for (int i= M.rows(); i >= 1; i--)
1718  {
1719  tmp3= 0;
1720  tmp1= L[i - 1];
1721  k= 0;
1722  for (int j= M.columns(); j >= 1; j--, k++)
1723  {
1724  tmp2= M (i, j);
1725  if (j == i)
1726  break;
1727  else
1728  {
1729  if (k > partialSol.size() - 1)
1730  tmp3 += tmp2*result[j - 1];
1731  else
1732  tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
1733  }
1734  }
1735  result [i - 1]= (tmp1 - tmp3)/tmp2;
1736  }
1737  return result;
1738 }
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 1690 of file cfModGcd.cc.

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

◆ solveGeneralVandermonde()

CFArray solveGeneralVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1634 of file cfModGcd.cc.

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

◆ solveSystemFp()

CFArray solveSystemFp ( const CFMatrix M,
const CFArray L 
)

Definition at line 1842 of file cfModGcd.cc.

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

◆ solveSystemFq()

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

Definition at line 1894 of file cfModGcd.cc.

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

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

◆ sparseGCDFp()

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

Definition at line 3564 of file cfModGcd.cc.

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

◆ sparseGCDFq()

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

Definition at line 3132 of file cfModGcd.cc.

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

◆ 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 283 of file cfModGcd.cc.

284 {
285  if (F.inBaseDomain())
286  return F.genOne();
287  if (F.level() == 1 && F.isUnivariate())
288  return F;
289  if (F.level() != 1 && F.isUnivariate())
290  return F.genOne();
291  if (degree (F,1) == 0) return F.genOne();
292 
293  int l= F.level();
294  if (l == 2)
295  return content(F);
296  else
297  {
298  CanonicalForm pol, c = 0;
299  CFIterator i = F;
300  for (; i.hasTerms(); i++)
301  {
302  pol= i.coeff();
303  pol= uni_content (pol);
304  c= gcd (c, pol);
305  if (c.isOne())
306  return c;
307  }
308  return c;
309  }
310 }
CanonicalForm 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 261 of file cfModGcd.cc.

262 {
263  if (F.inCoeffDomain())
264  return F.genOne();
265  if (F.level() == x.level() && F.isUnivariate())
266  return F;
267  if (F.level() != x.level() && F.isUnivariate())
268  return F.genOne();
269 
270  if (x.level() != 1)
271  {
272  CanonicalForm f= swapvar (F, x, Variable (1));
274  return swapvar (result, x, Variable (1));
275  }
276  else
277  return uni_content (F);
278 }
CanonicalForm 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 341 of file cfModGcd.cc.

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

◆ while()

while ( true  )

Definition at line 4134 of file cfModGcd.cc.

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

◆ cDn

Definition at line 4131 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4085 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4085 of file cfModGcd.cc.

◆ cl

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

Definition at line 4102 of file cfModGcd.cc.

◆ coF

Definition at line 67 of file cfModGcd.cc.

◆ cof

Definition at line 4131 of file cfModGcd.cc.

◆ cofn

Definition at line 4131 of file cfModGcd.cc.

◆ cofp

Definition at line 4131 of file cfModGcd.cc.

◆ coG

Definition at line 67 of file cfModGcd.cc.

◆ cog

Definition at line 4131 of file cfModGcd.cc.

◆ cogn

Definition at line 4131 of file cfModGcd.cc.

◆ cogp

Definition at line 4131 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4080 of file cfModGcd.cc.

◆ Dn

Definition at line 4098 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4080 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4128 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4083 of file cfModGcd.cc.

◆ false

return false

Definition at line 84 of file cfModGcd.cc.

◆ fp

Definition at line 4104 of file cfModGcd.cc.

◆ G

Definition at line 66 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4092 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4103 of file cfModGcd.cc.

◆ GG

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

Definition at line 4077 of file cfModGcd.cc.

◆ gp

Definition at line 4104 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4080 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4099 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4099 of file cfModGcd.cc.

◆ log2exp

const double log2exp = 1.442695041
static

Definition at line 89 of file cfModGcd.cc.

◆ maxNumVars

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

Definition at line 4132 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4106 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4131 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4131 of file cfModGcd.cc.

◆ p

int p

Definition at line 4080 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4098 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4084 of file cfModGcd.cc.