My Project
Loading...
Searching...
No Matches
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 <math.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 421 of file cfModGcd.cc.

422{
423 int i, m;
424 // extension of F_p needed
425 if (alpha.level() == 1)
426 {
427 i= 1;
428 m= 2;
429 } //extension of F_p(alpha)
430 if (alpha.level() != 1)
431 {
432 i= 4;
433 m= degree (getMipo (alpha));
434 }
435 #ifdef HAVE_FLINT
436 nmod_poly_t Irredpoly;
438 nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, i*m+1);
439 CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
440 nmod_poly_clear(Irredpoly);
441 #else
443 {
445 zz_p::init (getCharacteristic());
446 }
447 zz_pX NTLIrredpoly;
448 BuildIrred (NTLIrredpoly, i*m);
449 CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
450 #endif
451 return rootOf (newMipo);
452}
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)
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 cfEzgcd.cc:132
GLOBAL_VAR flint_rand_t FLINTrandom
Definition cf_random.cc:25
factory's main class
factory's class for variables
Definition factory.h:127
Variable alpha
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 2193 of file cfModGcd.cc.

2195{
2196 Aeval= A;
2197 Beval= B;
2198 int j= 1;
2199 for (CFListIterator i= L; i.hasItem(); i++, j++)
2200 {
2201 Aeval= Aeval (i.getItem(), j);
2202 Beval= Beval (i.getItem(), j);
2203 }
2204}
ListIterator< CanonicalForm > CFListIterator
b *CanonicalForm B
Definition facBivar.cc:52
CFList *& Aeval
<[in] poly
int j
Definition facHensel.cc:110
#define A
Definition sirandom.c:24

◆ evaluate()

CFArray evaluate ( const CFArray & A,
const CFList & evalPoints )

Definition at line 2039 of file cfModGcd.cc.

2040{
2041 CFArray result= A.size();
2042 CanonicalForm tmp;
2043 int k;
2044 for (int i= 0; i < A.size(); i++)
2045 {
2046 tmp= A[i];
2047 k= 1;
2048 for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2049 tmp= tmp (j.getItem(), k);
2050 result[i]= tmp;
2051 }
2052 return result;
2053}
Array< CanonicalForm > CFArray
int k
Definition cfEzgcd.cc:99
return result
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 2000 of file cfModGcd.cc.

2001{
2002 if (F.inCoeffDomain())
2003 {
2004 CFArray result= CFArray (1);
2005 result [0]= F;
2006 return result;
2007 }
2008 if (F.isUnivariate())
2009 {
2010 ASSERT (evalPoints.length() == 1,
2011 "expected an eval point with only one component");
2012 CFArray result= CFArray (size(F));
2013 int j= 0;
2015 for (CFIterator i= F; i.hasTerms(); i++, j++)
2016 result[j]= power (evalPoint, i.exp());
2017 return result;
2018 }
2019 int numMon= size (F);
2020 CFArray result= CFArray (numMon);
2021 int j= 0;
2024 buf.removeLast();
2025 CFArray recResult;
2026 CanonicalForm powEvalPoint;
2027 for (CFIterator i= F; i.hasTerms(); i++)
2028 {
2029 powEvalPoint= power (evalPoint, i.exp());
2030 recResult= evaluateMonom (i.coeff(), buf);
2031 for (int k= 0; k < recResult.size(); k++)
2032 result[j+k]= powEvalPoint*recResult[k];
2033 j += recResult.size();
2034 }
2035 return result;
2036}
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
List< CanonicalForm > CFList
CFArray evaluateMonom(const CanonicalForm &F, const CFList &evalPoints)
Definition cfModGcd.cc:2000
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
class to iterate through CanonicalForm's
Definition cf_iter.h:44
bool inCoeffDomain() const
bool isUnivariate() const
int status int void * buf
Definition si_signals.h:69

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

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

◆ extractContents()

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

Definition at line 312 of file cfModGcd.cc.

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

◆ for() [1/2]

for ( ; i,
0;i--  )

Definition at line 4125 of file cfModGcd.cc.

4126 {
4127 if (degree (f, i) <= 0 || degree (g, i) <= 0)
4128 continue;
4129 else
4131 }
g
Definition cfModGcd.cc:4098
int minCommonDeg
Definition cfModGcd.cc:4112
FILE * f
Definition checklibs.c:9
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)

◆ for() [2/2]

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

Definition at line 4113 of file cfModGcd.cc.

4114 {
4115 if (degree (f, i) <= 0 || degree (g, i) <= 0)
4116 continue;
4117 else
4118 {
4119 minCommonDeg= tmin (degree (g, i), degree (f, i));
4120 break;
4121 }
4122 }

◆ FpRandomElement()

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

Definition at line 1172 of file cfModGcd.cc.

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

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

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

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

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

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

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

◆ modGCDFp() [1/2]

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

Definition at line 1213 of file cfModGcd.cc.

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

◆ modGCDFp() [2/2]

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

Definition at line 1224 of file cfModGcd.cc.

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

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

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

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

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

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

◆ monicSparseInterpol()

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

Definition at line 2207 of file cfModGcd.cc.

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

◆ mult()

void mult ( CFList & L1,
const CFList & L2 )

multiply two lists componentwise

Definition at line 2184 of file cfModGcd.cc.

2185{
2186 ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2187
2188 CFListIterator j= L2;
2189 for (CFListIterator i= L1; i.hasItem(); i++, j++)
2190 i.getItem() *= j.getItem();
2191}

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

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

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

◆ nonMonicSparseInterpol()

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

Definition at line 2491 of file cfModGcd.cc.

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

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

◆ readOffSolution() [1/2]

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

Definition at line 1711 of file cfModGcd.cc.

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

◆ readOffSolution() [2/2]

CFArray readOffSolution ( const CFMatrix & M,
const long rk )

M in row echolon form, rk rank of M.

Definition at line 1689 of file cfModGcd.cc.

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

◆ solveGeneralVandermonde()

CFArray solveGeneralVandermonde ( const CFArray & M,
const CFArray & A )

Definition at line 1633 of file cfModGcd.cc.

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

◆ solveSystemFp()

CFArray solveSystemFp ( const CFMatrix & M,
const CFArray & L )

Definition at line 1844 of file cfModGcd.cc.

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

◆ solveSystemFq()

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

Definition at line 1896 of file cfModGcd.cc.

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

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

◆ sparseGCDFp()

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

Definition at line 3570 of file cfModGcd.cc.

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

◆ sparseGCDFq()

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

Definition at line 3138 of file cfModGcd.cc.

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

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

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

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

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

◆ while()

while ( true )

Definition at line 4140 of file cfModGcd.cc.

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

Variable Documentation

◆ b

else L b = 1

Definition at line 4111 of file cfModGcd.cc.

◆ cand

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

Definition at line 69 of file cfModGcd.cc.

◆ cd

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

Definition at line 4097 of file cfModGcd.cc.

◆ cDn

Definition at line 4137 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4091 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4091 of file cfModGcd.cc.

◆ cl

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

Definition at line 4108 of file cfModGcd.cc.

◆ coF

Definition at line 68 of file cfModGcd.cc.

◆ cof

Definition at line 4137 of file cfModGcd.cc.

◆ cofn

Definition at line 4137 of file cfModGcd.cc.

◆ cofp

Definition at line 4137 of file cfModGcd.cc.

◆ coG

◆ cog

Definition at line 4137 of file cfModGcd.cc.

◆ cogn

Definition at line 4137 of file cfModGcd.cc.

◆ cogp

Definition at line 4137 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4086 of file cfModGcd.cc.

◆ Dn

Definition at line 4104 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4086 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4134 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4089 of file cfModGcd.cc.

◆ false

return false

Definition at line 85 of file cfModGcd.cc.

◆ fp

Definition at line 4110 of file cfModGcd.cc.

◆ G

Definition at line 67 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4098 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4109 of file cfModGcd.cc.

◆ GG

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

Definition at line 4083 of file cfModGcd.cc.

◆ gp

Definition at line 4110 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4086 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4105 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4105 of file cfModGcd.cc.

◆ maxNumVars

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

Definition at line 4138 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4112 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4137 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4137 of file cfModGcd.cc.

◆ p

int p

Definition at line 4086 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4104 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4090 of file cfModGcd.cc.