My Project  UNKNOWN_GIT_VERSION
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 {
425  {
427  zz_p::init (getCharacteristic());
428  }
429  zz_pX NTLIrredpoly;
430  int i, m;
431  // extension of F_p needed
432  if (alpha.level() == 1)
433  {
434  i= 1;
435  m= 2;
436  } //extension of F_p(alpha)
437  if (alpha.level() != 1)
438  {
439  i= 4;
440  m= degree (getMipo (alpha));
441  }
442  BuildIrred (NTLIrredpoly, i*m);
443  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
444  return rootOf (newMipo);
445 }

◆ eval()

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

Definition at line 2126 of file cfModGcd.cc.

2128 {
2129  Aeval= A;
2130  Beval= B;
2131  int j= 1;
2132  for (CFListIterator i= L; i.hasItem(); i++, j++)
2133  {
2134  Aeval= Aeval (i.getItem(), j);
2135  Beval= Beval (i.getItem(), j);
2136  }
2137 }

◆ evaluate()

CFArray evaluate ( const CFArray A,
const CFList evalPoints 
)

Definition at line 1972 of file cfModGcd.cc.

1973 {
1974  CFArray result= A.size();
1975  CanonicalForm tmp;
1976  int k;
1977  for (int i= 0; i < A.size(); i++)
1978  {
1979  tmp= A[i];
1980  k= 1;
1981  for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
1982  tmp= tmp (j.getItem(), k);
1983  result[i]= tmp;
1984  }
1985  return result;
1986 }

◆ evaluateMonom()

CFArray evaluateMonom ( const CanonicalForm F,
const CFList evalPoints 
)

Definition at line 1933 of file cfModGcd.cc.

1934 {
1935  if (F.inCoeffDomain())
1936  {
1937  CFArray result= CFArray (1);
1938  result [0]= F;
1939  return result;
1940  }
1941  if (F.isUnivariate())
1942  {
1943  ASSERT (evalPoints.length() == 1,
1944  "expected an eval point with only one component");
1945  CFArray result= CFArray (size(F));
1946  int j= 0;
1948  for (CFIterator i= F; i.hasTerms(); i++, j++)
1949  result[j]= power (evalPoint, i.exp());
1950  return result;
1951  }
1952  int numMon= size (F);
1953  CFArray result= CFArray (numMon);
1954  int j= 0;
1957  buf.removeLast();
1958  CFArray recResult;
1959  CanonicalForm powEvalPoint;
1960  for (CFIterator i= F; i.hasTerms(); i++)
1961  {
1962  powEvalPoint= power (evalPoint, i.exp());
1963  recResult= evaluateMonom (i.coeff(), buf);
1964  for (int k= 0; k < recResult.size(); k++)
1965  result[j+k]= powEvalPoint*recResult[k];
1966  j += recResult.size();
1967  }
1968  return result;
1969 }

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

1994 {
1995  int k= tmax (F.level(), G.level()) - 1;
1996  Variable x= Variable (1);
1997  CFList result;
1998  FFRandom genFF;
1999  GFRandom genGF;
2000  int p= getCharacteristic ();
2001  double bound;
2002  if (alpha != Variable (1))
2003  {
2004  bound= pow ((double) p, (double) degree (getMipo(alpha)));
2005  bound= pow (bound, (double) k);
2006  }
2007  else if (GF)
2008  {
2009  bound= pow ((double) p, (double) getGFDegree());
2010  bound= pow ((double) bound, (double) k);
2011  }
2012  else
2013  bound= pow ((double) p, (double) k);
2014 
2015  CanonicalForm random;
2016  int j;
2017  bool zeroOneOccured= false;
2018  bool allEqual= false;
2020  do
2021  {
2022  random= 0;
2023  // possible overflow if list.length() does not fit into a int
2024  if (list.length() >= bound)
2025  {
2026  fail= true;
2027  break;
2028  }
2029  for (int i= 0; i < k; i++)
2030  {
2031  if (GF)
2032  {
2033  result.append (genGF.generate());
2034  random += result.getLast()*power (x, i);
2035  }
2036  else if (alpha.level() != 1)
2037  {
2038  AlgExtRandomF genAlgExt (alpha);
2039  result.append (genAlgExt.generate());
2040  random += result.getLast()*power (x, i);
2041  }
2042  else
2043  {
2044  result.append (genFF.generate());
2045  random += result.getLast()*power (x, i);
2046  }
2047  if (result.getLast().isOne() || result.getLast().isZero())
2048  zeroOneOccured= true;
2049  }
2050  if (find (list, random))
2051  {
2052  zeroOneOccured= false;
2053  allEqual= false;
2054  result= CFList();
2055  continue;
2056  }
2057  if (zeroOneOccured)
2058  {
2059  list.append (random);
2060  zeroOneOccured= false;
2061  allEqual= false;
2062  result= CFList();
2063  continue;
2064  }
2065  // no zero at this point
2066  if (k > 1)
2067  {
2068  allEqual= true;
2069  CFIterator iter= random;
2070  buf= iter.coeff();
2071  iter++;
2072  for (; iter.hasTerms(); iter++)
2073  if (buf != iter.coeff())
2074  allEqual= false;
2075  }
2076  if (allEqual)
2077  {
2078  list.append (random);
2079  allEqual= false;
2080  zeroOneOccured= false;
2081  result= CFList();
2082  continue;
2083  }
2084 
2085  Feval= F;
2086  Geval= G;
2087  CanonicalForm LCeval= LCF;
2088  j= 1;
2089  for (CFListIterator i= result; i.hasItem(); i++, j++)
2090  {
2091  Feval= Feval (i.getItem(), j);
2092  Geval= Geval (i.getItem(), j);
2093  LCeval= LCeval (i.getItem(), j);
2094  }
2095 
2096  if (LCeval.isZero())
2097  {
2098  if (!find (list, random))
2099  list.append (random);
2100  zeroOneOccured= false;
2101  allEqual= false;
2102  result= CFList();
2103  continue;
2104  }
2105 
2106  if (list.length() >= bound)
2107  {
2108  fail= true;
2109  break;
2110  }
2111  } while (find (list, random));
2112 
2113  return result;
2114 }

◆ 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 }

◆ for() [1/2]

for ( i,
0;i--   
)

Definition at line 4058 of file cfModGcd.cc.

4059  {
4060  if (degree (f, i) <= 0 || degree (g, i) <= 0)
4061  continue;
4062  else
4064  }

◆ for() [2/2]

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

Definition at line 4046 of file cfModGcd.cc.

4047  {
4048  if (degree (f, i) <= 0 || degree (g, i) <= 0)
4049  continue;
4050  else
4051  {
4052  minCommonDeg= tmin (degree (g, i), degree (f, i));
4053  break;
4054  }
4055  }

◆ FpRandomElement()

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

Definition at line 1159 of file cfModGcd.cc.

1160 {
1161  fail= false;
1162  Variable x= F.mvar();
1163  FFRandom genFF;
1164  CanonicalForm random;
1165  int p= getCharacteristic();
1166  int bound= p;
1167  do
1168  {
1169  if (list.length() == bound)
1170  {
1171  fail= true;
1172  break;
1173  }
1174  if (list.length() < 1)
1175  random= 0;
1176  else
1177  {
1178  random= genFF.generate();
1179  while (find (list, random))
1180  random= genFF.generate();
1181  }
1182  if (F (random, x) == 0)
1183  {
1184  list.append (random);
1185  continue;
1186  }
1187  } while (find (list, random));
1188  return random;
1189 }

◆ gaussianElimFp()

long gaussianElimFp ( CFMatrix M,
CFArray L 
)

Definition at line 1721 of file cfModGcd.cc.

1722 {
1723  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1724  CFMatrix *N;
1725  N= new CFMatrix (M.rows(), M.columns() + 1);
1726 
1727  for (int i= 1; i <= M.rows(); i++)
1728  for (int j= 1; j <= M.columns(); j++)
1729  (*N) (i, j)= M (i, j);
1730 
1731  int j= 1;
1732  for (int i= 0; i < L.size(); i++, j++)
1733  (*N) (j, M.columns() + 1)= L[i];
1734 #ifdef HAVE_FLINT
1735  nmod_mat_t FLINTN;
1736  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1737  long rk= nmod_mat_rref (FLINTN);
1738 
1739  delete N;
1740  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1741  nmod_mat_clear (FLINTN);
1742 #else
1743  int p= getCharacteristic ();
1744  if (fac_NTL_char != p)
1745  {
1746  fac_NTL_char= p;
1747  zz_p::init (p);
1748  }
1749  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1750  delete N;
1751  long rk= gauss (*NTLN);
1752 
1754  delete NTLN;
1755 #endif
1756 
1757  L= CFArray (M.rows());
1758  for (int i= 0; i < M.rows(); i++)
1759  L[i]= (*N) (i + 1, M.columns() + 1);
1760  M= (*N) (1, M.rows(), 1, M.columns());
1761  delete N;
1762  return rk;
1763 }

◆ gaussianElimFq()

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

Definition at line 1766 of file cfModGcd.cc.

1767 {
1768  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1769  CFMatrix *N;
1770  N= new CFMatrix (M.rows(), M.columns() + 1);
1771 
1772  for (int i= 1; i <= M.rows(); i++)
1773  for (int j= 1; j <= M.columns(); j++)
1774  (*N) (i, j)= M (i, j);
1775 
1776  int j= 1;
1777  for (int i= 0; i < L.size(); i++, j++)
1778  (*N) (j, M.columns() + 1)= L[i];
1779  int p= getCharacteristic ();
1780  if (fac_NTL_char != p)
1781  {
1782  fac_NTL_char= p;
1783  zz_p::init (p);
1784  }
1785  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1786  zz_pE::init (NTLMipo);
1787  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1788  long rk= gauss (*NTLN);
1789 
1790  delete N;
1792 
1793  delete NTLN;
1794 
1795  M= (*N) (1, M.rows(), 1, M.columns());
1796  L= CFArray (M.rows());
1797  for (int i= 0; i < M.rows(); i++)
1798  L[i]= (*N) (i + 1, M.columns() + 1);
1799 
1800  delete N;
1801  return rk;
1802 }

◆ getMonoms()

CFArray getMonoms ( const CanonicalForm F)

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

Parameters
[in]Fsome poly

Definition at line 1898 of file cfModGcd.cc.

1899 {
1900  if (F.inCoeffDomain())
1901  {
1902  CFArray result= CFArray (1);
1903  result [0]= 1;
1904  return result;
1905  }
1906  if (F.isUnivariate())
1907  {
1908  CFArray result= CFArray (size(F));
1909  int j= 0;
1910  for (CFIterator i= F; i.hasTerms(); i++, j++)
1911  result[j]= power (F.mvar(), i.exp());
1912  return result;
1913  }
1914  int numMon= size (F);
1915  CFArray result= CFArray (numMon);
1916  int j= 0;
1917  CFArray recResult;
1918  Variable x= F.mvar();
1919  CanonicalForm powX;
1920  for (CFIterator i= F; i.hasTerms(); i++)
1921  {
1922  powX= power (x, i.exp());
1923  recResult= getMonoms (i.coeff());
1924  for (int k= 0; k < recResult.size(); k++)
1925  result[j+k]= powX*recResult[k];
1926  j += recResult.size();
1927  }
1928  return result;
1929 }

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

808 {
809  fail= false;
810  Variable x= F.mvar();
811  GFRandom genGF;
812  CanonicalForm random;
813  int p= getCharacteristic();
814  int d= getGFDegree();
815  int bound= ipower (p, d);
816  do
817  {
818  if (list.length() == bound)
819  {
820  fail= true;
821  break;
822  }
823  if (list.length() < 1)
824  random= 0;
825  else
826  {
827  random= genGF.generate();
828  while (find (list, random))
829  random= genGF.generate();
830  }
831  if (F (random, x) == 0)
832  {
833  list.append (random);
834  continue;
835  }
836  } while (find (list, random));
837  return random;
838 }

◆ 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  }

◆ modGCDFp() [1/2]

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

Definition at line 1197 of file cfModGcd.cc.

1199 {
1200  CanonicalForm dummy1, dummy2;
1201  CanonicalForm result= modGCDFp (F, G, dummy1, dummy2, topLevel, l);
1202  return result;
1203 }

◆ modGCDFp() [2/2]

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

Definition at line 1206 of file cfModGcd.cc.

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

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

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

◆ modGCDFq() [2/2]

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

Definition at line 453 of file cfModGcd.cc.

455 {
456  CanonicalForm dummy1, dummy2;
457  CanonicalForm result= modGCDFq (F, G, dummy1, dummy2, alpha, l,
458  topLevel);
459  return result;
460 }

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

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

◆ modGCDGF() [2/2]

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

Definition at line 846 of file cfModGcd.cc.

848 {
849  CanonicalForm dummy1, dummy2;
850  CanonicalForm result= modGCDGF (F, G, dummy1, dummy2, l, topLevel);
851  return result;
852 }

◆ monicSparseInterpol()

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

Definition at line 2140 of file cfModGcd.cc.

2144 {
2145  CanonicalForm A= F;
2146  CanonicalForm B= G;
2147  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2148  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2149  else if (F.isZero() && G.isZero()) return F.genOne();
2150  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2151  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2152  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2153  if (F == G) return F/Lc(F);
2154 
2155  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2156  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2157 
2158  CFMap M,N;
2159  int best_level= myCompress (A, B, M, N, false);
2160 
2161  if (best_level == 0)
2162  return B.genOne();
2163 
2164  A= M(A);
2165  B= M(B);
2166 
2167  Variable x= Variable (1);
2168 
2169  //univariate case
2170  if (A.isUnivariate() && B.isUnivariate())
2171  return N (gcd (A, B));
2172 
2173  CanonicalForm skel= M(skeleton);
2174  CanonicalForm cA, cB; // content of A and B
2175  CanonicalForm ppA, ppB; // primitive part of A and B
2176  CanonicalForm gcdcAcB;
2177  cA = uni_content (A);
2178  cB = uni_content (B);
2179  gcdcAcB= gcd (cA, cB);
2180  ppA= A/cA;
2181  ppB= B/cB;
2182 
2183  CanonicalForm lcA, lcB; // leading coefficients of A and B
2184  CanonicalForm gcdlcAlcB;
2185  lcA= uni_lcoeff (ppA);
2186  lcB= uni_lcoeff (ppB);
2187 
2188  if (fdivides (lcA, lcB))
2189  {
2190  if (fdivides (A, B))
2191  return F/Lc(F);
2192  }
2193  if (fdivides (lcB, lcA))
2194  {
2195  if (fdivides (B, A))
2196  return G/Lc(G);
2197  }
2198 
2199  gcdlcAlcB= gcd (lcA, lcB);
2200  int skelSize= size (skel, skel.mvar());
2201 
2202  int j= 0;
2203  int biggestSize= 0;
2204 
2205  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2206  biggestSize= tmax (biggestSize, size (i.coeff()));
2207 
2208  CanonicalForm g, Aeval, Beval;
2209 
2211  bool evalFail= false;
2212  CFList list;
2213  bool GF= false;
2214  CanonicalForm LCA= LC (A);
2215  CanonicalForm tmp;
2216  CFArray gcds= CFArray (biggestSize);
2217  CFList * pEvalPoints= new CFList [biggestSize];
2218  Variable V_buf= alpha, V_buf4= alpha;
2219  CFList source, dest;
2220  CanonicalForm prim_elem, im_prim_elem;
2221  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2222  for (int i= 0; i < biggestSize; i++)
2223  {
2224  if (i == 0)
2225  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2226  list);
2227  else
2228  {
2229  mult (evalPoints, pEvalPoints [0]);
2230  eval (A, B, Aeval, Beval, evalPoints);
2231  }
2232 
2233  if (evalFail)
2234  {
2235  if (V_buf.level() != 1)
2236  {
2237  do
2238  {
2239  Variable V_buf3= V_buf;
2240  V_buf= chooseExtension (V_buf);
2241  source= CFList();
2242  dest= CFList();
2243 
2244  bool prim_fail= false;
2245  Variable V_buf2;
2246  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2247  if (V_buf4 == alpha && alpha.level() != 1)
2248  prim_elem_alpha= prim_elem;
2249 
2250  ASSERT (!prim_fail, "failure in integer factorizer");
2251  if (prim_fail)
2252  ; //ERROR
2253  else
2254  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2255 
2256  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2257  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2258 
2259  if (V_buf4 == alpha && alpha.level() != 1)
2260  im_prim_elem_alpha= im_prim_elem;
2261  else if (alpha.level() != 1)
2262  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2263  prim_elem, im_prim_elem, source, dest);
2264 
2265  for (CFListIterator j= list; j.hasItem(); j++)
2266  j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2267  im_prim_elem, source, dest);
2268  for (int k= 0; k < i; k++)
2269  {
2270  for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2271  j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2272  im_prim_elem, source, dest);
2273  gcds[k]= mapUp (gcds[k], V_buf4, V_buf, prim_elem, im_prim_elem,
2274  source, dest);
2275  }
2276 
2277  if (alpha.level() != 1)
2278  {
2279  A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2280  B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2281  }
2282  V_buf4= V_buf;
2283  evalFail= false;
2284  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2285  evalFail, list);
2286  } while (evalFail);
2287  }
2288  else
2289  {
2291  int deg= 2;
2292  bool initialized= false;
2293  do
2294  {
2295  mipo= randomIrredpoly (deg, x);
2296  if (initialized)
2297  setMipo (V_buf, mipo);
2298  else
2299  V_buf= rootOf (mipo);
2300  evalFail= false;
2301  initialized= true;
2302  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2303  evalFail, list);
2304  deg++;
2305  } while (evalFail);
2306  V_buf4= V_buf;
2307  }
2308  }
2309 
2310  g= gcd (Aeval, Beval);
2311  g /= Lc (g);
2312 
2313  if (degree (g) != degree (skel) || (skelSize != size (g)))
2314  {
2315  delete[] pEvalPoints;
2316  fail= true;
2317  if (alpha.level() != 1 && V_buf != alpha)
2318  prune1 (alpha);
2319  return 0;
2320  }
2321  CFIterator l= skel;
2322  for (CFIterator k= g; k.hasTerms(); k++, l++)
2323  {
2324  if (k.exp() != l.exp())
2325  {
2326  delete[] pEvalPoints;
2327  fail= true;
2328  if (alpha.level() != 1 && V_buf != alpha)
2329  prune1 (alpha);
2330  return 0;
2331  }
2332  }
2333  pEvalPoints[i]= evalPoints;
2334  gcds[i]= g;
2335 
2336  tmp= 0;
2337  int j= 0;
2338  for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2339  tmp += k.getItem()*power (x, j);
2340  list.append (tmp);
2341  }
2342 
2343  if (Monoms.size() == 0)
2344  Monoms= getMonoms (skel);
2345 
2346  coeffMonoms= new CFArray [skelSize];
2347  j= 0;
2348  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2349  coeffMonoms[j]= getMonoms (i.coeff());
2350 
2351  CFArray* pL= new CFArray [skelSize];
2352  CFArray* pM= new CFArray [skelSize];
2353  for (int i= 0; i < biggestSize; i++)
2354  {
2355  CFIterator l= gcds [i];
2356  evalPoints= pEvalPoints [i];
2357  for (int k= 0; k < skelSize; k++, l++)
2358  {
2359  if (i == 0)
2360  pL[k]= CFArray (biggestSize);
2361  pL[k] [i]= l.coeff();
2362 
2363  if (i == 0)
2364  pM[k]= evaluate (coeffMonoms [k], evalPoints);
2365  }
2366  }
2367 
2368  CFArray solution;
2369  CanonicalForm result= 0;
2370  int ind= 0;
2371  CFArray bufArray;
2372  CFMatrix Mat;
2373  for (int k= 0; k < skelSize; k++)
2374  {
2375  if (biggestSize != coeffMonoms[k].size())
2376  {
2377  bufArray= CFArray (coeffMonoms[k].size());
2378  for (int i= 0; i < coeffMonoms[k].size(); i++)
2379  bufArray [i]= pL[k] [i];
2380  solution= solveGeneralVandermonde (pM [k], bufArray);
2381  }
2382  else
2383  solution= solveGeneralVandermonde (pM [k], pL[k]);
2384 
2385  if (solution.size() == 0)
2386  {
2387  delete[] pEvalPoints;
2388  delete[] pM;
2389  delete[] pL;
2390  delete[] coeffMonoms;
2391  fail= true;
2392  if (alpha.level() != 1 && V_buf != alpha)
2393  prune1 (alpha);
2394  return 0;
2395  }
2396  for (int l= 0; l < solution.size(); l++)
2397  result += solution[l]*Monoms [ind + l];
2398  ind += solution.size();
2399  }
2400 
2401  delete[] pEvalPoints;
2402  delete[] pM;
2403  delete[] pL;
2404  delete[] coeffMonoms;
2405 
2406  if (alpha.level() != 1 && V_buf != alpha)
2407  {
2408  CFList u, v;
2409  result= mapDown (result, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2410  prune1 (alpha);
2411  }
2412 
2413  result= N(result);
2414  if (fdivides (result, F) && fdivides (result, G))
2415  return result;
2416  else
2417  {
2418  fail= true;
2419  return 0;
2420  }
2421 }

◆ mult()

void mult ( CFList L1,
const CFList L2 
)

multiply two lists componentwise

Definition at line 2117 of file cfModGcd.cc.

2118 {
2119  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2120 
2121  CFListIterator j= L2;
2122  for (CFListIterator i= L1; i.hasItem(); i++, j++)
2123  i.getItem() *= j.getItem();
2124 }

◆ 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 }

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

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

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

1693 {
1694  CFArray result= CFArray (M.rows());
1695  CanonicalForm tmp1, tmp2, tmp3;
1696  int k;
1697  for (int i= M.rows(); i >= 1; i--)
1698  {
1699  tmp3= 0;
1700  tmp1= L[i - 1];
1701  k= 0;
1702  for (int j= M.columns(); j >= 1; j--, k++)
1703  {
1704  tmp2= M (i, j);
1705  if (j == i)
1706  break;
1707  else
1708  {
1709  if (k > partialSol.size() - 1)
1710  tmp3 += tmp2*result[j - 1];
1711  else
1712  tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
1713  }
1714  }
1715  result [i - 1]= (tmp1 - tmp3)/tmp2;
1716  }
1717  return result;
1718 }

◆ readOffSolution() [2/2]

CFArray readOffSolution ( const CFMatrix M,
const long  rk 
)

M in row echolon form, rk rank of M.

Definition at line 1670 of file cfModGcd.cc.

1671 {
1672  CFArray result= CFArray (rk);
1673  CanonicalForm tmp1, tmp2, tmp3;
1674  for (int i= rk; i >= 1; i--)
1675  {
1676  tmp3= 0;
1677  tmp1= M (i, M.columns());
1678  for (int j= M.columns() - 1; j >= 1; j--)
1679  {
1680  tmp2= M (i, j);
1681  if (j == i)
1682  break;
1683  else
1684  tmp3 += tmp2*result[j - 1];
1685  }
1686  result[i - 1]= (tmp1 - tmp3)/tmp2;
1687  }
1688  return result;
1689 }

◆ solveGeneralVandermonde()

CFArray solveGeneralVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1614 of file cfModGcd.cc.

1615 {
1616  int r= M.size();
1617  ASSERT (A.size() == r, "vector does not have right size");
1618  if (r == 1)
1619  {
1620  CFArray result= CFArray (1);
1621  result [0]= A[0] / M [0];
1622  return result;
1623  }
1624  // check solvability
1625  bool notDistinct= false;
1626  for (int i= 0; i < r - 1; i++)
1627  {
1628  for (int j= i + 1; j < r; j++)
1629  {
1630  if (M [i] == M [j])
1631  {
1632  notDistinct= true;
1633  break;
1634  }
1635  }
1636  }
1637  if (notDistinct)
1638  return CFArray();
1639 
1640  CanonicalForm master= 1;
1641  Variable x= Variable (1);
1642  for (int i= 0; i < r; i++)
1643  master *= x - M [i];
1644  master *= x;
1645  CFList Pj;
1646  CanonicalForm tmp;
1647  for (int i= 0; i < r; i++)
1648  {
1649  tmp= master/(x - M [i]);
1650  tmp /= tmp (M [i], 1);
1651  Pj.append (tmp);
1652  }
1653 
1654  CFArray result= CFArray (r);
1655 
1656  CFListIterator j= Pj;
1657  for (int i= 1; i <= r; i++, j++)
1658  {
1659  tmp= 0;
1660 
1661  for (int l= 1; l <= A.size(); l++)
1662  tmp += A[l - 1]*j.getItem()[l];
1663  result[i - 1]= tmp;
1664  }
1665  return result;
1666 }

◆ solveSystemFp()

CFArray solveSystemFp ( const CFMatrix M,
const CFArray L 
)

Definition at line 1805 of file cfModGcd.cc.

1806 {
1807  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1808  CFMatrix *N;
1809  N= new CFMatrix (M.rows(), M.columns() + 1);
1810 
1811  for (int i= 1; i <= M.rows(); i++)
1812  for (int j= 1; j <= M.columns(); j++)
1813  (*N) (i, j)= M (i, j);
1814 
1815  int j= 1;
1816  for (int i= 0; i < L.size(); i++, j++)
1817  (*N) (j, M.columns() + 1)= L[i];
1818 
1819 #ifdef HAVE_FLINT
1820  nmod_mat_t FLINTN;
1821  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1822  long rk= nmod_mat_rref (FLINTN);
1823 #else
1824  int p= getCharacteristic ();
1825  if (fac_NTL_char != p)
1826  {
1827  fac_NTL_char= p;
1828  zz_p::init (p);
1829  }
1830  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1831  long rk= gauss (*NTLN);
1832 #endif
1833  delete N;
1834  if (rk != M.columns())
1835  {
1836 #ifdef HAVE_FLINT
1837  nmod_mat_clear (FLINTN);
1838 #else
1839  delete NTLN;
1840 #endif
1841  return CFArray();
1842  }
1843 #ifdef HAVE_FLINT
1844  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1845  nmod_mat_clear (FLINTN);
1846 #else
1848  delete NTLN;
1849 #endif
1850  CFArray A= readOffSolution (*N, rk);
1851 
1852  delete N;
1853  return A;
1854 }

◆ solveSystemFq()

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

Definition at line 1857 of file cfModGcd.cc.

1858 {
1859  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1860  CFMatrix *N;
1861  N= new CFMatrix (M.rows(), M.columns() + 1);
1862 
1863  for (int i= 1; i <= M.rows(); i++)
1864  for (int j= 1; j <= M.columns(); j++)
1865  (*N) (i, j)= M (i, j);
1866  int j= 1;
1867  for (int i= 0; i < L.size(); i++, j++)
1868  (*N) (j, M.columns() + 1)= L[i];
1869  int p= getCharacteristic ();
1870  if (fac_NTL_char != p)
1871  {
1872  fac_NTL_char= p;
1873  zz_p::init (p);
1874  }
1875  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1876  zz_pE::init (NTLMipo);
1877  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1878  long rk= gauss (*NTLN);
1879 
1880  delete N;
1881  if (rk != M.columns())
1882  {
1883  delete NTLN;
1884  return CFArray();
1885  }
1887 
1888  delete NTLN;
1889 
1890  CFArray A= readOffSolution (*N, rk);
1891 
1892  delete N;
1893  return A;
1894 }

◆ solveVandermonde()

CFArray solveVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1561 of file cfModGcd.cc.

1562 {
1563  int r= M.size();
1564  ASSERT (A.size() == r, "vector does not have right size");
1565 
1566  if (r == 1)
1567  {
1568  CFArray result= CFArray (1);
1569  result [0]= A [0] / M [0];
1570  return result;
1571  }
1572  // check solvability
1573  bool notDistinct= false;
1574  for (int i= 0; i < r - 1; i++)
1575  {
1576  for (int j= i + 1; j < r; j++)
1577  {
1578  if (M [i] == M [j])
1579  {
1580  notDistinct= true;
1581  break;
1582  }
1583  }
1584  }
1585  if (notDistinct)
1586  return CFArray();
1587 
1588  CanonicalForm master= 1;
1589  Variable x= Variable (1);
1590  for (int i= 0; i < r; i++)
1591  master *= x - M [i];
1592  CFList Pj;
1593  CanonicalForm tmp;
1594  for (int i= 0; i < r; i++)
1595  {
1596  tmp= master/(x - M [i]);
1597  tmp /= tmp (M [i], 1);
1598  Pj.append (tmp);
1599  }
1600  CFArray result= CFArray (r);
1601 
1602  CFListIterator j= Pj;
1603  for (int i= 1; i <= r; i++, j++)
1604  {
1605  tmp= 0;
1606  for (int l= 0; l < A.size(); l++)
1607  tmp += A[l]*j.getItem()[l];
1608  result[i - 1]= tmp;
1609  }
1610  return result;
1611 }

◆ sparseGCDFp()

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

Definition at line 3503 of file cfModGcd.cc.

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

◆ sparseGCDFq()

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

Definition at line 3071 of file cfModGcd.cc.

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

◆ 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 }

◆ 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 }

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

4074  {
4075  p = cf_getBigPrime( i );
4076  i--;
4077  while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
4078  {
4079  p = cf_getBigPrime( i );
4080  i--;
4081  }
4082  //printf("try p=%d\n",p);
4083  setCharacteristic( p );
4084  fp= mapinto (f);
4085  gp= mapinto (g);
4086  TIMING_START (modZ_recursion)
4087 #ifdef HAVE_NTL
4088  if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
4089  Dp = modGCDFp (fp, gp, cofp, cogp);
4090  else
4091  {
4092  Dp= gcd_poly (fp, gp);
4093  cofp= fp/Dp;
4094  cogp= gp/Dp;
4095  }
4096 #else
4097  Dp= gcd_poly (fp, gp);
4098  cofp= fp/Dp;
4099  cogp= gp/Dp;
4100 #endif
4101  TIMING_END_AND_PRINT (modZ_recursion,
4102  "time for gcd mod p in modular gcd: ");
4103  Dp /=Dp.lc();
4104  Dp *= mapinto (cl);
4105  cofp /= lc (cofp);
4106  cofp *= mapinto (lcf);
4107  cogp /= lc (cogp);
4108  cogp *= mapinto (lcg);
4109  setCharacteristic( 0 );
4110  dp_deg=totaldegree(Dp);
4111  if ( dp_deg == 0 )
4112  {
4113  //printf(" -> 1\n");
4114  return CanonicalForm(gcdcfcg);
4115  }
4116  if ( q.isZero() )
4117  {
4118  D = mapinto( Dp );
4119  cof= mapinto (cofp);
4120  cog= mapinto (cogp);
4121  d_deg=dp_deg;
4122  q = p;
4123  Dn= balance_p (D, p);
4124  cofn= balance_p (cof, p);
4125  cogn= balance_p (cog, p);
4126  }
4127  else
4128  {
4129  if ( dp_deg == d_deg )
4130  {
4131  chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
4132  chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
4133  chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
4134  cof= newCof;
4135  cog= newCog;
4136  newqh= newq/2;
4137  Dn= balance_p (newD, newq, newqh);
4138  cofn= balance_p (newCof, newq, newqh);
4139  cogn= balance_p (newCog, newq, newqh);
4140  if (test != Dn) //balance_p (newD, newq))
4141  test= balance_p (newD, newq);
4142  else
4143  equal= true;
4144  q = newq;
4145  D = newD;
4146  }
4147  else if ( dp_deg < d_deg )
4148  {
4149  //printf(" were all bad, try more\n");
4150  // all previous p's are bad primes
4151  q = p;
4152  D = mapinto( Dp );
4153  cof= mapinto (cof);
4154  cog= mapinto (cog);
4155  d_deg=dp_deg;
4156  test= 0;
4157  equal= false;
4158  Dn= balance_p (D, p);
4159  cofn= balance_p (cof, p);
4160  cogn= balance_p (cog, p);
4161  }
4162  else
4163  {
4164  //printf(" was bad, try more\n");
4165  }
4166  //else dp_deg > d_deg: bad prime
4167  }
4168  if ( i >= 0 )
4169  {
4170  cDn= icontent (Dn);
4171  Dn /= cDn;
4172  cofn /= cl/cDn;
4173  //cofn /= icontent (cofn);
4174  cogn /= cl/cDn;
4175  //cogn /= icontent (cogn);
4176  TIMING_START (modZ_termination);
4177  if ((terminationTest (f,g, cofn, cogn, Dn)) ||
4178  ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
4179  {
4180  TIMING_END_AND_PRINT (modZ_termination,
4181  "time for successful termination in modular gcd: ");
4182  //printf(" -> success\n");
4183  return Dn*gcdcfcg;
4184  }
4185  TIMING_END_AND_PRINT (modZ_termination,
4186  "time for unsuccessful termination in modular gcd: ");
4187  equal= false;
4188  //else: try more primes
4189  }
4190  else
4191  { // try other method
4192  //printf("try other gcd\n");
4194  D=gcd_poly( f, g );
4196  return D*gcdcfcg;
4197  }
4198  }

Variable Documentation

◆ b

else L b = 1

Definition at line 4044 of file cfModGcd.cc.

◆ cand

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

Definition at line 68 of file cfModGcd.cc.

◆ cd

cd = bCommonDen( GG )

Definition at line 4030 of file cfModGcd.cc.

◆ cDn

Definition at line 4070 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4024 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4024 of file cfModGcd.cc.

◆ cl

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

Definition at line 4041 of file cfModGcd.cc.

◆ coF

Definition at line 67 of file cfModGcd.cc.

◆ cof

Definition at line 4070 of file cfModGcd.cc.

◆ cofn

Definition at line 4070 of file cfModGcd.cc.

◆ cofp

Definition at line 4070 of file cfModGcd.cc.

◆ coG

Definition at line 67 of file cfModGcd.cc.

◆ cog

Definition at line 4070 of file cfModGcd.cc.

◆ cogn

Definition at line 4070 of file cfModGcd.cc.

◆ cogp

Definition at line 4070 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4019 of file cfModGcd.cc.

◆ Dn

Definition at line 4037 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4019 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4067 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4022 of file cfModGcd.cc.

◆ false

return false

Definition at line 84 of file cfModGcd.cc.

◆ fp

Definition at line 4043 of file cfModGcd.cc.

◆ G

Definition at line 66 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4031 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4042 of file cfModGcd.cc.

◆ GG

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

Definition at line 4016 of file cfModGcd.cc.

◆ gp

Definition at line 4043 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4019 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4038 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4038 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 4071 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4045 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4070 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4070 of file cfModGcd.cc.

◆ p

int p

Definition at line 4019 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4037 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4023 of file cfModGcd.cc.

test
CanonicalForm test
Definition: cfModGcd.cc:4037
Matrix
Definition: ftmpl_matrix.h:29
fac_NTL_char
long fac_NTL_char
Definition: NTLconvert.cc:41
newtonInterp
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
evaluate
CFArray evaluate(const CFArray &A, const CFList &evalPoints)
Definition: cfModGcd.cc:1972
cog
CanonicalForm cog
Definition: cfModGcd.cc:4070
FFRandom
generate random elements in F_p
Definition: cf_random.h:44
convertNTLmat_zz_pE2FacCFMatrix
CFMatrix * convertNTLmat_zz_pE2FacCFMatrix(const mat_zz_pE &m, const Variable &alpha)
Definition: NTLconvert.cc:1211
convertFacCF2NTLzzpX
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:100
j
int j
Definition: facHensel.cc:105
k
int k
Definition: cfEzgcd.cc:92
equal
bool equal
Definition: cfModGcd.cc:4067
terminationTest
bool terminationTest(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &coF, const CanonicalForm &coG, const CanonicalForm &cand)
icontent
CanonicalForm icontent(const CanonicalForm &f)
CanonicalForm icontent ( const CanonicalForm & f )
Definition: cf_gcd.cc:71
convertFacCFMatrix2NTLmat_zz_pE
mat_zz_pE * convertFacCFMatrix2NTLmat_zz_pE(const CFMatrix &m)
Definition: NTLconvert.cc:1195
CFIterator
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
x
Variable x
Definition: cfModGcd.cc:4023
DEBOUTLN
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
dp_deg
int dp_deg
Definition: cfModGcd.cc:4019
uni_content
static CanonicalForm uni_content(const CanonicalForm &F)
compute the content of F, where F is considered as an element of
Definition: cfModGcd.cc:283
result
return result
Definition: facAbsBiFact.cc:76
setMipo
void setMipo(const Variable &alpha, const CanonicalForm &mipo)
Definition: variable.cc:219
CFArray
Array< CanonicalForm > CFArray
Definition: canonicalform.h:390
i
int i
Definition: cfModGcd.cc:4019
CanonicalForm::inBaseDomain
bool inBaseDomain() const
Definition: canonicalform.cc:101
GFRandomElement
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:807
degrees
int * degrees(const CanonicalForm &f, int *degs=0)
int * degrees ( const CanonicalForm & f, int * degs )
Definition: cf_ops.cc:493
FpRandomElement
static CanonicalForm FpRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
Definition: cfModGcd.cc:1159
solveSystemFp
CFArray solveSystemFp(const CFMatrix &M, const CFArray &L)
Definition: cfModGcd.cc:1805
DELETE_ARRAY
#define DELETE_ARRAY(P)
Definition: cf_defs.h:49
cof
CanonicalForm cof
Definition: cfModGcd.cc:4070
Matrix::rows
int rows() const
Definition: ftmpl_matrix.h:43
CFMap
class CFMap
Definition: cf_map.h:85
lcf
CanonicalForm lcf
Definition: cfModGcd.cc:4038
power
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Definition: canonicalform.cc:1837
randomElement
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
AlgExtRandomF
generate random elements in F_p(alpha)
Definition: cf_random.h:70
GFRandom::generate
CanonicalForm generate() const
Definition: cf_random.cc:66
g
g
Definition: cfModGcd.cc:4031
convertNTLmat_zz_p2FacCFMatrix
CFMatrix * convertNTLmat_zz_p2FacCFMatrix(const mat_zz_p &m)
Definition: NTLconvert.cc:1182
mod
CF_NO_INLINE CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:564
rootOf
Variable rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables
Definition: variable.cc:162
readOffSolution
CFArray readOffSolution(const CFMatrix &M, const long rk)
M in row echolon form, rk rank of M.
Definition: cfModGcd.cc:1670
amp::floor
const signed long floor(const ampf< Precision > &x)
Definition: amp.h:774
CFMatrix
Matrix< CanonicalForm > CFMatrix
Definition: canonicalform.h:391
N
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:49
randomIrredpoly
CanonicalForm randomIrredpoly(int i, const Variable &x)
computes a random monic irreducible univariate polynomial in x over Fp of degree i via NTL
Definition: cf_irred.cc:42
iter
CFFListIterator iter
Definition: facAbsBiFact.cc:54
cf_getBigPrime
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
gf_name
char gf_name
Definition: gfops.cc:52
primitiveElement
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:310
getCharacteristic
int getCharacteristic()
Definition: cf_char.cc:51
b
CanonicalForm b
Definition: cfModGcd.cc:4044
ind2
long ind2(long arg)
Definition: kutil.cc:4060
cl
cl
Definition: cfModGcd.cc:4041
solveGeneralVandermonde
CFArray solveGeneralVandermonde(const CFArray &M, const CFArray &A)
Definition: cfModGcd.cc:1614
CanonicalForm
factory's main class
Definition: canonicalform.h:83
amp::ceil
const signed long ceil(const ampf< Precision > &x)
Definition: amp.h:789
degsg
int * degsg
Definition: cfEzgcd.cc:53
GFMapDown
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:243
cogp
CanonicalForm cogp
Definition: cfModGcd.cc:4070
getMonoms
CFArray getMonoms(const CanonicalForm &F)
extract monomials of F, parts in algebraic variable are considered coefficients
Definition: cfModGcd.cc:1898
CanonicalForm::isOne
CF_NO_INLINE bool isOne() const
CF_INLINE bool CanonicalForm::isOne, isZero () const.
Definition: cf_inline.cc:354
CanonicalForm::genOne
CanonicalForm genOne() const
Definition: canonicalform.cc:1822
gcdcfcg
CanonicalForm gcdcfcg
Definition: cfModGcd.cc:4042
mapUp
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
Definition: cf_map_ext.cc:67
Lc
CanonicalForm Lc(const CanonicalForm &f)
Definition: canonicalform.h:300
swapvar
CanonicalForm swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
Array
Definition: ftmpl_array.h:17
getMipo
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
fp
CanonicalForm fp
Definition: cfModGcd.cc:4043
ASSERT
#define ASSERT(expression, message)
Definition: cf_assert.h:99
abs
Rational abs(const Rational &a)
Definition: GMPrat.cc:439
M
#define M
Definition: sirandom.c:24
buf
int status int void * buf
Definition: si_signals.h:59
content
CanonicalForm content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:175
Array::size
int size() const
Definition: ftmpl_array.cc:92
TIMING_START
TIMING_START(fac_alg_resultant)
lc
CanonicalForm lc(const CanonicalForm &f)
Definition: canonicalform.h:297
Aeval
CFList *& Aeval
<[in] poly
Definition: facFactorize.h:31
evalPoint
static void evalPoint(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &FEval, CanonicalForm &GEval, CFGenerator &evalPoint)
Definition: cfModResultant.cc:311
evalPoints
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...
Definition: facFqFactorize.cc:740
alpha
Variable alpha
Definition: facAbsBiFact.cc:52
f
f
Definition: cfModGcd.cc:4022
lcg
CanonicalForm lcg
Definition: cfModGcd.cc:4038
mult
void mult(CFList &L1, const CFList &L2)
multiply two lists componentwise
Definition: cfModGcd.cc:2117
Matrix::columns
int columns() const
Definition: ftmpl_matrix.h:44
D
#define D(A)
Definition: gentable.cc:131
gp
CanonicalForm gp
Definition: cfModGcd.cc:4043
setCharacteristic
void setCharacteristic(int c)
Definition: cf_char.cc:23
coF
const CanonicalForm const CanonicalForm & coF
Definition: cfModGcd.cc:67
modGCDFq
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:467
uni_lcoeff
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
size
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
prune
void prune(Variable &alpha)
Definition: variable.cc:261
coG
const CanonicalForm const CanonicalForm const CanonicalForm & coG
Definition: cfModGcd.cc:67
tmin
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
chineseRemainder
void 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:52
SW_USE_CHINREM_GCD
static const int SW_USE_CHINREM_GCD
set to 1 to use modular gcd over Z
Definition: cf_defs.h:38
f_zero
int f_zero
Definition: cfEzgcd.cc:62
cogn
CanonicalForm cogn
Definition: cfModGcd.cc:4070
both_zero
int both_zero
Definition: cfGcdAlgExt.cc:71
Variable::level
int level() const
Definition: factory.h:134
Feval
CanonicalForm Feval
Definition: facAbsFact.cc:64
tmp1
CFList tmp1
Definition: facFqBivar.cc:70
convertFacCFMatrix2NTLmat_zz_p
mat_zz_p * convertFacCFMatrix2NTLmat_zz_p(const CFMatrix &m)
Definition: NTLconvert.cc:1166
GFRandom
generate random elements in GF
Definition: cf_random.h:32
ipower
int ipower(int b, int m)
int ipower ( int b, int m )
Definition: cf_util.cc:25
g_zero
int g_zero
Definition: cfEzgcd.cc:63
monicSparseInterpol
CanonicalForm monicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition: cfModGcd.cc:2140
newCof
CanonicalForm newCof
Definition: cfModGcd.cc:4070
fdivides
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
Definition: cf_algorithm.cc:338
Off
void Off(int sw)
switches
Definition: canonicalform.cc:1905
Dn
CanonicalForm Dn
Definition: cfModGcd.cc:4037
both_non_zero
const CanonicalForm CFMap CFMap int & both_non_zero
Definition: cfEzgcd.cc:50
getGFDegree
int getGFDegree()
Definition: cf_char.cc:56
tmax
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
cand
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition: cfModGcd.cc:69
modGCDFp
CanonicalForm modGCDFp(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:1206
CanonicalForm::inCoeffDomain
bool inCoeffDomain() const
Definition: canonicalform.cc:119
List::length
int length() const
Definition: ftmpl_list.cc:273
TIMING_END_AND_PRINT
TIMING_END_AND_PRINT(fac_alg_resultant, "time to compute resultant0: ")
degsf
int * degsf
Definition: cfEzgcd.cc:52
Variable
factory's class for variables
Definition: factory.h:118
log
gmp_float log(const gmp_float &a)
Definition: mpr_complex.cc:344
bound
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
nonMonicSparseInterpol
CanonicalForm nonMonicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition: cfModGcd.cc:2424
B
b *CanonicalForm B
Definition: facBivar.cc:52
topLevel
const CanonicalForm CFMap CFMap bool topLevel
Definition: cfGcdAlgExt.cc:57
gaussianElimFp
long gaussianElimFp(CFMatrix &M, CFArray &L)
Definition: cfModGcd.cc:1721
gcd_poly
CanonicalForm gcd_poly(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
Definition: cf_gcd.cc:91
CanonicalForm::mvar
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
Definition: canonicalform.cc:560
mapPrimElem
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:377
LC
CanonicalForm LC(const CanonicalForm &f)
Definition: canonicalform.h:303
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: cfModGcd.cc:1989
mapinto
CanonicalForm mapinto(const CanonicalForm &f)
Definition: canonicalform.h:348
d_deg
int d_deg
Definition: cfModGcd.cc:4019
eval
void eval(const CanonicalForm &A, const CanonicalForm &B, CanonicalForm &Aeval, CanonicalForm &Beval, const CFList &L)
Definition: cfModGcd.cc:2126
pow
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:414
FFRandom::generate
CanonicalForm generate() const
Definition: cf_random.cc:56
balance_p
CanonicalForm balance_p(const CanonicalForm &f, const CanonicalForm &q, const CanonicalForm &qh)
same as balance_p ( const CanonicalForm & f, const CanonicalForm & q ) but qh= q/2 is provided,...
Definition: cfGcdUtil.cc:227
m
int m
Definition: cfEzgcd.cc:121
find
template bool find(const List< CanonicalForm > &, const CanonicalForm &)
NULL
#define NULL
Definition: omList.c:10
l
int l
Definition: cfEzgcd.cc:93
maxNumVars
int maxNumVars
Definition: cfModGcd.cc:4071
CanonicalForm::isUnivariate
bool isUnivariate() const
Definition: canonicalform.cc:152
gcd
int gcd(int a, int b)
Definition: walkSupport.cc:836
newCog
CanonicalForm newCog
Definition: cfModGcd.cc:4070
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
p
int p
Definition: cfModGcd.cc:4019
mipo
CanonicalForm mipo
Definition: facAlgExt.cc:57
List< CanonicalForm >
convertNmod_mat_t2FacCFMatrix
CFMatrix * convertNmod_mat_t2FacCFMatrix(const nmod_mat_t m)
conversion of a FLINT matrix over Z/p to a factory matrix
Definition: FLINTconvert.cc:489
evaluateMonom
CFArray evaluateMonom(const CanonicalForm &F, const CFList &evalPoints)
Definition: cfModGcd.cc:1933
CFList
List< CanonicalForm > CFList
Definition: canonicalform.h:388
H
CanonicalForm H
Definition: facAbsFact.cc:64
CanonicalForm::isZero
CF_NO_INLINE bool isZero() const
Definition: cf_inline.cc:372
tmp2
CFList tmp2
Definition: facFqBivar.cc:70
gaussianElimFq
long gaussianElimFq(CFMatrix &M, CFArray &L, const Variable &alpha)
Definition: cfModGcd.cc:1766
GFMapUp
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:207
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: cfModGcd.cc:93
CanonicalForm::level
int level() const
level() returns the level of CO.
Definition: canonicalform.cc:543
convertFacCFMatrix2nmod_mat_t
void convertFacCFMatrix2nmod_mat_t(nmod_mat_t M, const CFMatrix &m)
conversion of a factory matrix over Z/p to a nmod_mat_t
Definition: FLINTconvert.cc:471
List::append
void append(const T &)
Definition: ftmpl_list.cc:256
A
#define A
Definition: sirandom.c:23
prune1
void prune1(const Variable &alpha)
Definition: variable.cc:291
degree
int degree(const CanonicalForm &f)
Definition: canonicalform.h:309
minCommonDeg
int minCommonDeg
Definition: cfModGcd.cc:4045
LCF
CanonicalForm LCF
Definition: facFactorize.cc:53
G
const CanonicalForm & G
Definition: cfModGcd.cc:66
On
void On(int sw)
switches
Definition: canonicalform.cc:1898
sparseGCDFq
CanonicalForm sparseGCDFq(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
Definition: cfModGcd.cc:3071
convertNTLzzpX2CF
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
Definition: NTLconvert.cc:248
solveSystemFq
CFArray solveSystemFq(const CFMatrix &M, const CFArray &L, const Variable &alpha)
Definition: cfModGcd.cc:1857
List::getLast
T getLast() const
Definition: ftmpl_list.cc:309
totaldegree
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition: cf_ops.cc:523
modGCDGF
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:860
cofp
CanonicalForm cofp
Definition: cfModGcd.cc:4070
chooseExtension
static Variable chooseExtension(const Variable &alpha)
Definition: cfModGcd.cc:422
NEW_ARRAY
#define NEW_ARRAY(T, N)
Definition: cf_defs.h:48
mapDown
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:90
ListIterator
Definition: ftmpl_list.h:87
cofn
CanonicalForm cofn
Definition: cfModGcd.cc:4070
cDn
CanonicalForm cDn
Definition: cfModGcd.cc:4070
sparseGCDFp
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:3503