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
 
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}] $
 
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.
 
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})
 
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
 
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.
 
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
 
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.
 
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.
 
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
 
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
 
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
 
 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:4077
GLOBAL_VAR flint_rand_t FLINTrandom
Definition: cf_random.cc:25
factory's main class
Definition: canonicalform.h:86
factory's class for variables
Definition: factory.h:127
int level() const
Definition: factory.h:143
Variable alpha
Definition: facAbsBiFact.cc:51
nmod_poly_init(FLINTmipo, getCharacteristic())
nmod_poly_clear(FLINTmipo)
Variable FACTORY_PUBLIC rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables
Definition: variable.cc:162
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207

◆ eval()

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

Definition at line 2184 of file cfModGcd.cc.

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

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

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

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

4117 {
4118 if (degree (f, i) <= 0 || degree (g, i) <= 0)
4119 continue;
4120 else
4122 }
f
Definition: cfModGcd.cc:4080
g
Definition: cfModGcd.cc:4089
int minCommonDeg
Definition: cfModGcd.cc:4103
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)

◆ for() [2/2]

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

Definition at line 4104 of file cfModGcd.cc.

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

◆ 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
1829 M= (*N) (1, M.rows(), 1, M.columns());
1830 L= CFArray (M.rows());
1831 for (int i= 0; i < M.rows(); i++)
1832 L[i]= (*N) (i + 1, M.columns() + 1);
1833
1834 delete N;
1835 return rk;
1836}
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 1956 of file cfModGcd.cc.

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

◆ 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
void FACTORY_PUBLIC prune(Variable &alpha)
Definition: variable.cc:261
void setMipo(const Variable &alpha, const CanonicalForm &mipo)
Definition: variable.cc:219
#define TIMING_START(t)
Definition: timing.h:92
#define TIMING_END_AND_PRINT(t, msg)
Definition: timing.h:94

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

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

◆ mult()

void mult ( CFList L1,
const CFList L2 
)

multiply two lists componentwise

Definition at line 2175 of file cfModGcd.cc.

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

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

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

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

◆ solveSystemFq()

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

Definition at line 1891 of file cfModGcd.cc.

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

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

◆ sparseGCDFq()

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

Definition at line 3129 of file cfModGcd.cc.

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

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

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

◆ cDn

Definition at line 4128 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4082 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4082 of file cfModGcd.cc.

◆ cl

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

Definition at line 4099 of file cfModGcd.cc.

◆ coF

Definition at line 67 of file cfModGcd.cc.

◆ cof

Definition at line 4128 of file cfModGcd.cc.

◆ cofn

Definition at line 4128 of file cfModGcd.cc.

◆ cofp

Definition at line 4128 of file cfModGcd.cc.

◆ coG

Definition at line 67 of file cfModGcd.cc.

◆ cog

Definition at line 4128 of file cfModGcd.cc.

◆ cogn

Definition at line 4128 of file cfModGcd.cc.

◆ cogp

Definition at line 4128 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4077 of file cfModGcd.cc.

◆ Dn

Definition at line 4095 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4077 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4125 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4080 of file cfModGcd.cc.

◆ false

return false

Definition at line 84 of file cfModGcd.cc.

◆ fp

Definition at line 4101 of file cfModGcd.cc.

◆ G

Definition at line 66 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4089 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4100 of file cfModGcd.cc.

◆ GG

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

Definition at line 4074 of file cfModGcd.cc.

◆ gp

Definition at line 4101 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4077 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4096 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4096 of file cfModGcd.cc.

◆ maxNumVars

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

Definition at line 4129 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4103 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4128 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4128 of file cfModGcd.cc.

◆ p

int p

Definition at line 4077 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4095 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4081 of file cfModGcd.cc.