My Project
Loading...
Searching...
No Matches
Functions | Variables
cfModGcd.cc File Reference

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

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

Go to the source code of this file.

Functions

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

Variables

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

Detailed Description

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

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

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

Definition in file cfModGcd.cc.

Function Documentation

◆ chooseExtension()

static Variable chooseExtension ( const Variable alpha)
inlinestatic

Definition at line 420 of file cfModGcd.cc.

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

◆ eval()

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

Definition at line 2185 of file cfModGcd.cc.

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

◆ evaluate()

CFArray evaluate ( const CFArray A,
const CFList evalPoints 
)

Definition at line 2031 of file cfModGcd.cc.

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

◆ evaluateMonom()

CFArray evaluateMonom ( const CanonicalForm F,
const CFList evalPoints 
)

Definition at line 1992 of file cfModGcd.cc.

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

◆ evaluationPoints()

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

Definition at line 2048 of file cfModGcd.cc.

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

◆ extractContents()

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

Definition at line 311 of file cfModGcd.cc.

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

◆ for() [1/2]

for ( i,
0;i--   
)

Definition at line 4117 of file cfModGcd.cc.

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

◆ for() [2/2]

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

Definition at line 4105 of file cfModGcd.cc.

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

◆ FpRandomElement()

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

Definition at line 1171 of file cfModGcd.cc.

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

◆ gaussianElimFp()

long gaussianElimFp ( CFMatrix M,
CFArray L 
)

Definition at line 1739 of file cfModGcd.cc.

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

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

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

◆ GFRandomElement()

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

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

Definition at line 819 of file cfModGcd.cc.

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

◆ if() [1/2]

if ( i  = =0)

◆ if() [2/2]

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

Definition at line 71 of file cfModGcd.cc.

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

◆ modGCDFp() [1/2]

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

Definition at line 1212 of file cfModGcd.cc.

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

◆ modGCDFp() [2/2]

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

Definition at line 1223 of file cfModGcd.cc.

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

◆ modGCDFq() [1/2]

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

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

Definition at line 478 of file cfModGcd.cc.

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

◆ modGCDFq() [2/2]

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

Definition at line 462 of file cfModGcd.cc.

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

◆ modGCDGF() [1/2]

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

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

Definition at line 872 of file cfModGcd.cc.

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

◆ modGCDGF() [2/2]

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

Definition at line 858 of file cfModGcd.cc.

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

◆ monicSparseInterpol()

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

Definition at line 2199 of file cfModGcd.cc.

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

◆ mult()

void mult ( CFList L1,
const CFList L2 
)

multiply two lists componentwise

Definition at line 2176 of file cfModGcd.cc.

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

◆ myCompress()

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

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

Definition at line 91 of file cfModGcd.cc.

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

◆ newtonInterp()

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

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

Definition at line 364 of file cfModGcd.cc.

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

◆ nonMonicSparseInterpol()

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

Definition at line 2483 of file cfModGcd.cc.

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

◆ randomElement()

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

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

Definition at line 379 of file cfModGcd.cc.

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

◆ readOffSolution() [1/2]

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

Definition at line 1710 of file cfModGcd.cc.

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

◆ readOffSolution() [2/2]

CFArray readOffSolution ( const CFMatrix M,
const long  rk 
)

M in row echolon form, rk rank of M.

Definition at line 1688 of file cfModGcd.cc.

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

◆ solveGeneralVandermonde()

CFArray solveGeneralVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1632 of file cfModGcd.cc.

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

◆ solveSystemFp()

CFArray solveSystemFp ( const CFMatrix M,
const CFArray L 
)

Definition at line 1840 of file cfModGcd.cc.

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

◆ solveSystemFq()

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

Definition at line 1892 of file cfModGcd.cc.

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

◆ solveVandermonde()

CFArray solveVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1579 of file cfModGcd.cc.

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

◆ sparseGCDFp()

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

Definition at line 3562 of file cfModGcd.cc.

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

◆ sparseGCDFq()

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

Definition at line 3130 of file cfModGcd.cc.

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

◆ TIMING_DEFINE_PRINT() [1/2]

TIMING_DEFINE_PRINT ( gcd_recursion  ) const &

◆ TIMING_DEFINE_PRINT() [2/2]

TIMING_DEFINE_PRINT ( modZ_termination  ) const &

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

◆ uni_content() [1/2]

static CanonicalForm uni_content ( const CanonicalForm F)
inlinestatic

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

Definition at line 281 of file cfModGcd.cc.

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

◆ uni_content() [2/2]

CanonicalForm uni_content ( const CanonicalForm F,
const Variable x 
)

Definition at line 259 of file cfModGcd.cc.

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

◆ uni_lcoeff()

static CanonicalForm uni_lcoeff ( const CanonicalForm F)
inlinestatic

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

Definition at line 339 of file cfModGcd.cc.

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

◆ while()

while ( true  )

Definition at line 4132 of file cfModGcd.cc.

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

Variable Documentation

◆ b

else L b = 1

Definition at line 4103 of file cfModGcd.cc.

◆ cand

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

Definition at line 68 of file cfModGcd.cc.

◆ cd

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

Definition at line 4089 of file cfModGcd.cc.

◆ cDn

Definition at line 4129 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4083 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4083 of file cfModGcd.cc.

◆ cl

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

Definition at line 4100 of file cfModGcd.cc.

◆ coF

Definition at line 67 of file cfModGcd.cc.

◆ cof

Definition at line 4129 of file cfModGcd.cc.

◆ cofn

Definition at line 4129 of file cfModGcd.cc.

◆ cofp

Definition at line 4129 of file cfModGcd.cc.

◆ coG

Definition at line 67 of file cfModGcd.cc.

◆ cog

Definition at line 4129 of file cfModGcd.cc.

◆ cogn

Definition at line 4129 of file cfModGcd.cc.

◆ cogp

Definition at line 4129 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4078 of file cfModGcd.cc.

◆ Dn

Definition at line 4096 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4078 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4126 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4081 of file cfModGcd.cc.

◆ false

return false

Definition at line 84 of file cfModGcd.cc.

◆ fp

Definition at line 4102 of file cfModGcd.cc.

◆ G

Definition at line 66 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4090 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4101 of file cfModGcd.cc.

◆ GG

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

Definition at line 4075 of file cfModGcd.cc.

◆ gp

Definition at line 4102 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4078 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4097 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4097 of file cfModGcd.cc.

◆ maxNumVars

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

Definition at line 4130 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4104 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4129 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4129 of file cfModGcd.cc.

◆ p

int p

Definition at line 4078 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4096 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4082 of file cfModGcd.cc.