My Project
Functions | Variables
cfModGcd.cc File Reference

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

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

Go to the source code of this file.

Functions

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

Variables

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

Detailed Description

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

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

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

Definition in file cfModGcd.cc.

Function Documentation

◆ chooseExtension()

static Variable chooseExtension ( const Variable alpha)
inlinestatic

Definition at line 422 of file cfModGcd.cc.

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

◆ eval()

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

Definition at line 2187 of file cfModGcd.cc.

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

◆ evaluate()

CFArray evaluate ( const CFArray A,
const CFList evalPoints 
)

Definition at line 2033 of file cfModGcd.cc.

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

◆ evaluateMonom()

CFArray evaluateMonom ( const CanonicalForm F,
const CFList evalPoints 
)

Definition at line 1994 of file cfModGcd.cc.

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

◆ evaluationPoints()

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

Definition at line 2050 of file cfModGcd.cc.

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

◆ extractContents()

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

Definition at line 313 of file cfModGcd.cc.

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

◆ for() [1/2]

for ( i,
0;i--   
)

Definition at line 4119 of file cfModGcd.cc.

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

◆ for() [2/2]

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

Definition at line 4107 of file cfModGcd.cc.

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

◆ FpRandomElement()

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

Definition at line 1173 of file cfModGcd.cc.

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

◆ gaussianElimFp()

long gaussianElimFp ( CFMatrix M,
CFArray L 
)

Definition at line 1741 of file cfModGcd.cc.

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

◆ gaussianElimFq()

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

Definition at line 1786 of file cfModGcd.cc.

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

◆ getMonoms()

CFArray getMonoms ( const CanonicalForm F)

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

Parameters
[in]Fsome poly

Definition at line 1959 of file cfModGcd.cc.

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

◆ GFRandomElement()

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

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

Definition at line 821 of file cfModGcd.cc.

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

◆ if() [1/2]

if ( i  = =0)

◆ if() [2/2]

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

Definition at line 71 of file cfModGcd.cc.

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

◆ modGCDFp() [1/2]

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

Definition at line 1214 of file cfModGcd.cc.

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

◆ modGCDFp() [2/2]

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

Definition at line 1225 of file cfModGcd.cc.

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

◆ modGCDFq() [1/2]

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

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

Definition at line 480 of file cfModGcd.cc.

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

◆ modGCDFq() [2/2]

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

Definition at line 464 of file cfModGcd.cc.

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

◆ modGCDGF() [1/2]

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

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

Definition at line 874 of file cfModGcd.cc.

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

◆ modGCDGF() [2/2]

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

Definition at line 860 of file cfModGcd.cc.

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

◆ monicSparseInterpol()

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

Definition at line 2201 of file cfModGcd.cc.

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

◆ mult()

void mult ( CFList L1,
const CFList L2 
)

multiply two lists componentwise

Definition at line 2178 of file cfModGcd.cc.

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

◆ myCompress()

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

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

Definition at line 93 of file cfModGcd.cc.

95{
96 int n= tmax (F.level(), G.level());
97 int * degsf= NEW_ARRAY(int,n + 1);
98 int * degsg= NEW_ARRAY(int,n + 1);
99
100 for (int i = n; i >= 0; i--)
101 degsf[i]= degsg[i]= 0;
102
103 degsf= degrees (F, degsf);
104 degsg= degrees (G, degsg);
105
106 int both_non_zero= 0;
107 int f_zero= 0;
108 int g_zero= 0;
109 int both_zero= 0;
110
111 if (topLevel)
112 {
113 for (int i= 1; i <= n; i++)
114 {
115 if (degsf[i] != 0 && degsg[i] != 0)
116 {
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;
169 l= 0;
170 int i= 1;
171 while (k > 0)
172 {
173 if (degsf [i] != 0 && degsg [i] != 0)
174 min_max_deg= tmax (degsf[i], degsg[i]);
175 else
176 min_max_deg= 0;
177 while (min_max_deg == 0)
178 {
179 i++;
180 if (degsf [i] != 0 && degsg [i] != 0)
181 min_max_deg= tmax (degsf[i], degsg[i]);
182 else
183 min_max_deg= 0;
184 }
185 for (int j= i + 1; j <= m; j++)
186 {
187 if (degsf[j] != 0 && degsg [j] != 0 &&
188 tmax (degsf[j],degsg[j]) <= min_max_deg)
189 {
190 min_max_deg= tmax (degsf[j], degsg[j]);
191 l= j;
192 }
193 }
194 if (l != 0)
195 {
196 if (l != k)
197 {
198 M.newpair (Variable (l), Variable(k));
199 N.newpair (Variable (k), Variable(l));
200 degsf[l]= 0;
201 degsg[l]= 0;
202 l= 0;
203 }
204 else
205 {
206 degsf[l]= 0;
207 degsg[l]= 0;
208 l= 0;
209 }
210 }
211 else if (l == 0)
212 {
213 if (i != k)
214 {
215 M.newpair (Variable (i), Variable (k));
216 N.newpair (Variable (k), Variable (i));
217 degsf[i]= 0;
218 degsg[i]= 0;
219 }
220 else
221 {
222 degsf[i]= 0;
223 degsg[i]= 0;
224 }
225 i++;
226 }
227 k--;
228 }
229 }
230 else
231 {
232 //arrange Variables such that no gaps occur
233 for (int i= 1; i <= n; i++)
234 {
235 if (degsf[i] == 0 && degsg[i] == 0)
236 {
237 both_zero++;
238 continue;
239 }
240 else
241 {
242 if (both_zero != 0)
243 {
244 M.newpair (Variable (i), Variable (i - both_zero));
245 N.newpair (Variable (i - both_zero), Variable (i));
246 }
247 }
248 }
249 }
250
253
254 return 1;
255}
int * degrees(const CanonicalForm &f, int *degs=0)
int * degrees ( const CanonicalForm & f, int * degs )
Definition: cf_ops.cc:493
int * degsf
Definition: cfEzgcd.cc:59
int f_zero
Definition: cfEzgcd.cc:69
const CanonicalForm CFMap CFMap int & both_non_zero
Definition: cfEzgcd.cc:57
int g_zero
Definition: cfEzgcd.cc:70
int * degsg
Definition: cfEzgcd.cc:60
int both_zero
Definition: cfGcdAlgExt.cc:71
#define DELETE_ARRAY(P)
Definition: cf_defs.h:64
#define NEW_ARRAY(T, N)
Definition: cf_defs.h:63

◆ newtonInterp()

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

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

Definition at line 366 of file cfModGcd.cc.

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

◆ nonMonicSparseInterpol()

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

Definition at line 2485 of file cfModGcd.cc.

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

◆ randomElement()

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

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

Definition at line 381 of file cfModGcd.cc.

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

◆ readOffSolution() [1/2]

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

Definition at line 1712 of file cfModGcd.cc.

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

◆ readOffSolution() [2/2]

CFArray readOffSolution ( const CFMatrix M,
const long  rk 
)

M in row echolon form, rk rank of M.

Definition at line 1690 of file cfModGcd.cc.

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

◆ solveGeneralVandermonde()

CFArray solveGeneralVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1634 of file cfModGcd.cc.

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

◆ solveSystemFp()

CFArray solveSystemFp ( const CFMatrix M,
const CFArray L 
)

Definition at line 1842 of file cfModGcd.cc.

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

◆ solveSystemFq()

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

Definition at line 1894 of file cfModGcd.cc.

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

◆ solveVandermonde()

CFArray solveVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1581 of file cfModGcd.cc.

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

◆ sparseGCDFp()

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

Definition at line 3564 of file cfModGcd.cc.

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

◆ sparseGCDFq()

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

Definition at line 3132 of file cfModGcd.cc.

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

◆ TIMING_DEFINE_PRINT() [1/2]

TIMING_DEFINE_PRINT ( gcd_recursion  ) const &

◆ TIMING_DEFINE_PRINT() [2/2]

TIMING_DEFINE_PRINT ( modZ_termination  ) const &

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

◆ uni_content() [1/2]

static CanonicalForm uni_content ( const CanonicalForm F)
inlinestatic

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

Definition at line 283 of file cfModGcd.cc.

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

◆ uni_content() [2/2]

CanonicalForm uni_content ( const CanonicalForm F,
const Variable x 
)

Definition at line 261 of file cfModGcd.cc.

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

◆ uni_lcoeff()

static CanonicalForm uni_lcoeff ( const CanonicalForm F)
inlinestatic

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

Definition at line 341 of file cfModGcd.cc.

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

◆ while()

while ( true  )

Definition at line 4134 of file cfModGcd.cc.

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

Variable Documentation

◆ b

else L b = 1

Definition at line 4105 of file cfModGcd.cc.

◆ cand

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

Definition at line 68 of file cfModGcd.cc.

◆ cd

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

Definition at line 4091 of file cfModGcd.cc.

◆ cDn

Definition at line 4131 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4085 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4085 of file cfModGcd.cc.

◆ cl

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

Definition at line 4102 of file cfModGcd.cc.

◆ coF

Definition at line 67 of file cfModGcd.cc.

◆ cof

Definition at line 4131 of file cfModGcd.cc.

◆ cofn

Definition at line 4131 of file cfModGcd.cc.

◆ cofp

Definition at line 4131 of file cfModGcd.cc.

◆ coG

Definition at line 67 of file cfModGcd.cc.

◆ cog

Definition at line 4131 of file cfModGcd.cc.

◆ cogn

Definition at line 4131 of file cfModGcd.cc.

◆ cogp

Definition at line 4131 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4080 of file cfModGcd.cc.

◆ Dn

Definition at line 4098 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4080 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4128 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4083 of file cfModGcd.cc.

◆ false

return false

Definition at line 84 of file cfModGcd.cc.

◆ fp

Definition at line 4104 of file cfModGcd.cc.

◆ G

Definition at line 66 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4092 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4103 of file cfModGcd.cc.

◆ GG

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

Definition at line 4077 of file cfModGcd.cc.

◆ gp

Definition at line 4104 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4080 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4099 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4099 of file cfModGcd.cc.

◆ log2exp

const double log2exp = 1.442695041
static

Definition at line 89 of file cfModGcd.cc.

◆ maxNumVars

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

Definition at line 4132 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4106 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4131 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4131 of file cfModGcd.cc.

◆ p

int p

Definition at line 4080 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4098 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4084 of file cfModGcd.cc.