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

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

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

Go to the source code of this file.

Functions

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

Variables

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

Detailed Description

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

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

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

Definition in file cfModGcd.cc.

Function Documentation

◆ chooseExtension()

static Variable chooseExtension ( const Variable alpha)
inlinestatic

Definition at line 420 of file cfModGcd.cc.

421{
422 int i, m;
423 // extension of F_p needed
424 if (alpha.level() == 1)
425 {
426 i= 1;
427 m= 2;
428 } //extension of F_p(alpha)
429 if (alpha.level() != 1)
430 {
431 i= 4;
432 m= degree (getMipo (alpha));
433 }
434 #ifdef HAVE_FLINT
440 #else
442 {
444 zz_p::init (getCharacteristic());
445 }
449 #endif
450 return rootOf (newMipo);
451}
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
VAR long fac_NTL_char
Definition NTLconvert.cc:46
int degree(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition cf_char.cc:70
int m
Definition cfEzgcd.cc:128
int i
Definition cfModGcd.cc:4078
GLOBAL_VAR flint_rand_t FLINTrandom
Definition cf_random.cc:25
factory's main class
factory's class for variables
Definition factory.h:127
int level() const
Definition factory.h:143
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 2185 of file cfModGcd.cc.

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

◆ evaluate()

CFArray evaluate ( const CFArray A,
const CFList evalPoints 
)

Definition at line 2031 of file cfModGcd.cc.

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

◆ evaluateMonom()

CFArray evaluateMonom ( const CanonicalForm F,
const CFList evalPoints 
)

Definition at line 1992 of file cfModGcd.cc.

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

◆ evaluationPoints()

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

Definition at line 2048 of file cfModGcd.cc.

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

◆ extractContents()

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

Definition at line 311 of file cfModGcd.cc.

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

◆ for() [1/2]

for ( i,
0;i--   
)

Definition at line 4117 of file cfModGcd.cc.

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

◆ for() [2/2]

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

Definition at line 4105 of file cfModGcd.cc.

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

◆ FpRandomElement()

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

Definition at line 1171 of file cfModGcd.cc.

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

◆ gaussianElimFp()

long gaussianElimFp ( CFMatrix M,
CFArray L 
)

Definition at line 1739 of file cfModGcd.cc.

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

◆ gaussianElimFq()

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

Definition at line 1784 of file cfModGcd.cc.

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

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

◆ GFRandomElement()

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

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

Definition at line 819 of file cfModGcd.cc.

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

◆ if() [1/2]

if ( i  = =0)

◆ if() [2/2]

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

Definition at line 71 of file cfModGcd.cc.

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

◆ modGCDFp() [1/2]

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

Definition at line 1212 of file cfModGcd.cc.

1214{
1217 return result;
1218}
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:1223

◆ modGCDFp() [2/2]

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

Definition at line 1223 of file cfModGcd.cc.

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

◆ modGCDFq() [1/2]

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

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

Definition at line 478 of file cfModGcd.cc.

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

◆ modGCDFq() [2/2]

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

Definition at line 462 of file cfModGcd.cc.

464{
467 topLevel);
468 return result;
469}

◆ modGCDGF() [1/2]

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

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

Definition at line 872 of file cfModGcd.cc.

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

860{
863 return result;
864}

◆ monicSparseInterpol()

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

Definition at line 2199 of file cfModGcd.cc.

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

◆ mult()

void mult ( CFList L1,
const CFList L2 
)

multiply two lists componentwise

Definition at line 2176 of file cfModGcd.cc.

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

◆ myCompress()

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

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

Definition at line 91 of file cfModGcd.cc.

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

◆ newtonInterp()

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

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

Definition at line 364 of file cfModGcd.cc.

367{
369
371 *newtonPoly;
372 return interPoly;
373}

◆ nonMonicSparseInterpol()

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

Definition at line 2483 of file cfModGcd.cc.

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

◆ randomElement()

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

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

Definition at line 379 of file cfModGcd.cc.

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

◆ readOffSolution() [1/2]

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

Definition at line 1710 of file cfModGcd.cc.

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

◆ readOffSolution() [2/2]

CFArray readOffSolution ( const CFMatrix M,
const long  rk 
)

M in row echolon form, rk rank of M.

Definition at line 1688 of file cfModGcd.cc.

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

◆ solveGeneralVandermonde()

CFArray solveGeneralVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1632 of file cfModGcd.cc.

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

◆ solveSystemFp()

CFArray solveSystemFp ( const CFMatrix M,
const CFArray L 
)

Definition at line 1840 of file cfModGcd.cc.

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

◆ solveSystemFq()

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

Definition at line 1892 of file cfModGcd.cc.

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

◆ solveVandermonde()

CFArray solveVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1579 of file cfModGcd.cc.

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

◆ sparseGCDFp()

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

Definition at line 3562 of file cfModGcd.cc.

3564{
3565 CanonicalForm A= F;
3566 CanonicalForm B= G;
3567 if (F.isZero() && degree(G) > 0) return G/Lc(G);
3568 else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3569 else if (F.isZero() && G.isZero()) return F.genOne();
3570 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3571 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3572 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3573 if (F == G) return F/Lc(F);
3574
3575 CFMap M,N;
3576 int best_level= myCompress (A, B, M, N, topLevel);
3577
3578 if (best_level == 0) return B.genOne();
3579
3580 A= M(A);
3581 B= M(B);
3582
3583 Variable x= Variable (1);
3584
3585 //univariate case
3586 if (A.isUnivariate() && B.isUnivariate())
3587 return N (gcd (A, B));
3588
3589 CanonicalForm cA, cB; // content of A and B
3590 CanonicalForm ppA, ppB; // primitive part of A and B
3592
3593 cA = uni_content (A);
3594 cB = uni_content (B);
3595 gcdcAcB= gcd (cA, cB);
3596 ppA= A/cA;
3597 ppB= B/cB;
3598
3599 CanonicalForm lcA, lcB; // leading coefficients of A and B
3601 lcA= uni_lcoeff (ppA);
3602 lcB= uni_lcoeff (ppB);
3603
3604 if (fdivides (lcA, lcB))
3605 {
3606 if (fdivides (A, B))
3607 return F/Lc(F);
3608 }
3609 if (fdivides (lcB, lcA))
3610 {
3611 if (fdivides (B, A))
3612 return G/Lc(G);
3613 }
3614
3615 gcdlcAlcB= gcd (lcA, lcB);
3616
3617 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3618 int d0;
3619
3620 if (d == 0)
3621 return N(gcdcAcB);
3622
3623 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3624
3625 if (d0 < d)
3626 d= d0;
3627
3628 if (d == 0)
3629 return N(gcdcAcB);
3630
3633 m= gcdlcAlcB;
3634 G_m= 0;
3635 H= 0;
3636 bool fail= false;
3637 topLevel= false;
3638 bool inextension= false;
3641 CFList source, dest;
3642 do //first do
3643 {
3644 if (inextension)
3646 else
3648 if (random_element == 0 && !fail)
3649 {
3650 if (!find (l, random_element))
3652 continue;
3653 }
3654
3655 if (!fail && !inextension)
3656 {
3657 CFList list;
3661 list);
3663 "time for recursive call: ");
3664 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3665 }
3666 else if (!fail && inextension)
3667 {
3668 CFList list;
3672 list, topLevel);
3674 "time for recursive call: ");
3675 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3676 }
3677 else if (fail && !inextension)
3678 {
3679 source= CFList();
3680 dest= CFList();
3681 CFList list;
3683 int deg= 2;
3684 bool initialized= false;
3685 do
3686 {
3687 mipo= randomIrredpoly (deg, x);
3688 if (initialized)
3689 setMipo (alpha, mipo);
3690 else
3691 alpha= rootOf (mipo);
3692 inextension= true;
3693 fail= false;
3694 initialized= true;
3696 deg++;
3697 } while (fail);
3698 cleanUp= alpha;
3699 V_buf= alpha;
3700 list= CFList();
3704 list, topLevel);
3706 "time for recursive call: ");
3707 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3708 }
3709 else if (fail && inextension)
3710 {
3711 source= CFList();
3712 dest= CFList();
3713
3716 bool prim_fail= false;
3720
3721 if (V_buf3 != alpha)
3722 {
3726 dest);
3730 dest);
3731 for (CFListIterator i= l; i.hasItem(); i++)
3732 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3733 source, dest);
3734 }
3735
3736 ASSERT (!prim_fail, "failure in integer factorizer");
3737 if (prim_fail)
3738 ; //ERROR
3739 else
3741
3742 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3743 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3744
3745 for (CFListIterator i= l; i.hasItem(); i++)
3746 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3747 im_prim_elem, source, dest);
3751 source, dest);
3755 source, dest);
3756 fail= false;
3758 DEBOUTLN (cerr, "fail= " << fail);
3759 CFList list;
3763 list, topLevel);
3765 "time for recursive call: ");
3766 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3767 alpha= V_buf;
3768 }
3769
3770 if (!G_random_element.inCoeffDomain())
3772 Variable (G_random_element.level()));
3773 else
3774 d0= 0;
3775
3776 if (d0 == 0)
3777 {
3778 if (inextension)
3779 prune (cleanUp);
3780 return N(gcdcAcB);
3781 }
3782 if (d0 > d)
3783 {
3784 if (!find (l, random_element))
3786 continue;
3787 }
3788
3792
3794
3795 if (!G_random_element.inCoeffDomain())
3797 Variable (G_random_element.level()));
3798 else
3799 d0= 0;
3800
3801 if (d0 < d)
3802 {
3803 m= gcdlcAlcB;
3804 newtonPoly= 1;
3805 G_m= 0;
3806 d= d0;
3807 }
3808
3810
3811 if (uni_lcoeff (H) == gcdlcAlcB)
3812 {
3813 cH= uni_content (H);
3814 ppH= H/cH;
3815 ppH /= Lc (ppH);
3816 DEBOUTLN (cerr, "ppH= " << ppH);
3817
3818 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3819 {
3820 if (inextension)
3821 prune (cleanUp);
3822 return N(gcdcAcB*ppH);
3823 }
3824 }
3825 G_m= H;
3827 m= m*(x - random_element);
3828 if (!find (l, random_element))
3830
3831 if ((getCharacteristic() > 3 && size (skeleton) < 200))
3832 {
3835
3836 do //second do
3837 {
3838 if (inextension)
3840 else
3842 if (random_element == 0 && !fail)
3843 {
3844 if (!find (l, random_element))
3846 continue;
3847 }
3848
3849 bool sparseFail= false;
3850 if (!fail && !inextension)
3851 {
3852 CFList list;
3854 if (LC (skeleton).inCoeffDomain())
3858 Monoms);
3859 else
3865 "time for recursive call: ");
3866 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3867 }
3868 else if (!fail && inextension)
3869 {
3870 CFList list;
3872 if (LC (skeleton).inCoeffDomain())
3876 Monoms);
3877 else
3881 Monoms);
3883 "time for recursive call: ");
3884 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3885 }
3886 else if (fail && !inextension)
3887 {
3888 source= CFList();
3889 dest= CFList();
3890 CFList list;
3892 int deg= 2;
3893 bool initialized= false;
3894 do
3895 {
3896 mipo= randomIrredpoly (deg, x);
3897 if (initialized)
3898 setMipo (alpha, mipo);
3899 else
3900 alpha= rootOf (mipo);
3901 inextension= true;
3902 fail= false;
3903 initialized= true;
3905 deg++;
3906 } while (fail);
3907 cleanUp= alpha;
3908 V_buf= alpha;
3909 list= CFList();
3911 if (LC (skeleton).inCoeffDomain())
3915 Monoms);
3916 else
3920 Monoms);
3922 "time for recursive call: ");
3923 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3924 }
3925 else if (fail && inextension)
3926 {
3927 source= CFList();
3928 dest= CFList();
3929
3932 bool prim_fail= false;
3936
3937 if (V_buf3 != alpha)
3938 {
3942 source, dest);
3946 source, dest);
3947 for (CFListIterator i= l; i.hasItem(); i++)
3948 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3949 source, dest);
3950 }
3951
3952 ASSERT (!prim_fail, "failure in integer factorizer");
3953 if (prim_fail)
3954 ; //ERROR
3955 else
3957
3958 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3959 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3960
3961 for (CFListIterator i= l; i.hasItem(); i++)
3962 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3963 im_prim_elem, source, dest);
3967 source, dest);
3971 source, dest);
3972 fail= false;
3974 DEBOUTLN (cerr, "fail= " << fail);
3975 CFList list;
3977 if (LC (skeleton).inCoeffDomain())
3981 Monoms);
3982 else
3988 "time for recursive call: ");
3989 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3990 alpha= V_buf;
3991 }
3992
3993 if (sparseFail)
3994 break;
3995
3996 if (!G_random_element.inCoeffDomain())
3998 Variable (G_random_element.level()));
3999 else
4000 d0= 0;
4001
4002 if (d0 == 0)
4003 {
4004 if (inextension)
4005 prune (cleanUp);
4006 return N(gcdcAcB);
4007 }
4008 if (d0 > d)
4009 {
4010 if (!find (l, random_element))
4012 continue;
4013 }
4014
4018
4019 if (!G_random_element.inCoeffDomain())
4021 Variable (G_random_element.level()));
4022 else
4023 d0= 0;
4024
4025 if (d0 < d)
4026 {
4027 m= gcdlcAlcB;
4028 newtonPoly= 1;
4029 G_m= 0;
4030 d= d0;
4031 }
4032
4036 "time for newton interpolation: ");
4037
4038 //termination test
4039 if (uni_lcoeff (H) == gcdlcAlcB)
4040 {
4041 cH= uni_content (H);
4042 ppH= H/cH;
4043 ppH /= Lc (ppH);
4044 DEBOUTLN (cerr, "ppH= " << ppH);
4045 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4046 {
4047 if (inextension)
4048 prune (cleanUp);
4049 return N(gcdcAcB*ppH);
4050 }
4051 }
4052
4053 G_m= H;
4055 m= m*(x - random_element);
4056 if (!find (l, random_element))
4058
4059 } while (1); //end of second do
4060 }
4061 else
4062 {
4063 if (inextension)
4064 prune (cleanUp);
4065 return N(gcdcAcB*modGCDFp (ppA, ppB));
4066 }
4067 } while (1); //end of first do
4068}
CanonicalForm sparseGCDFq(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
Definition cfModGcd.cc:3130
CanonicalForm monicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition cfModGcd.cc:2199
CanonicalForm nonMonicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition cfModGcd.cc:2483
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
Definition cfModGcd.cc:3562

◆ sparseGCDFq()

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

Definition at line 3130 of file cfModGcd.cc.

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

◆ TIMING_DEFINE_PRINT() [1/2]

TIMING_DEFINE_PRINT ( gcd_recursion  ) const &

◆ TIMING_DEFINE_PRINT() [2/2]

TIMING_DEFINE_PRINT ( modZ_termination  ) const &

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

◆ uni_content() [1/2]

static CanonicalForm uni_content ( const CanonicalForm F)
inlinestatic

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

Definition at line 281 of file cfModGcd.cc.

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

◆ uni_content() [2/2]

CanonicalForm uni_content ( const CanonicalForm F,
const Variable x 
)

Definition at line 259 of file cfModGcd.cc.

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

◆ uni_lcoeff()

static CanonicalForm uni_lcoeff ( const CanonicalForm F)
inlinestatic

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

Definition at line 339 of file cfModGcd.cc.

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

◆ while()

while ( true  )

Definition at line 4132 of file cfModGcd.cc.

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

Variable Documentation

◆ b

else L b = 1

Definition at line 4103 of file cfModGcd.cc.

◆ cand

Initial value:

Definition at line 68 of file cfModGcd.cc.

◆ cd

cd ( bCommonDen(FF ) = bCommonDen( GG )

Definition at line 4089 of file cfModGcd.cc.

◆ cDn

Definition at line 4129 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4083 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4083 of file cfModGcd.cc.

◆ cl

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

Definition at line 4100 of file cfModGcd.cc.

◆ coF

Definition at line 67 of file cfModGcd.cc.

◆ cof

Definition at line 4129 of file cfModGcd.cc.

◆ cofn

Definition at line 4129 of file cfModGcd.cc.

◆ cofp

Definition at line 4129 of file cfModGcd.cc.

◆ coG

Definition at line 67 of file cfModGcd.cc.

◆ cog

Definition at line 4129 of file cfModGcd.cc.

◆ cogn

Definition at line 4129 of file cfModGcd.cc.

◆ cogp

Definition at line 4129 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4078 of file cfModGcd.cc.

◆ Dn

Definition at line 4096 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4078 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4126 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4081 of file cfModGcd.cc.

◆ false

return false

Definition at line 84 of file cfModGcd.cc.

◆ fp

Definition at line 4102 of file cfModGcd.cc.

◆ G

Definition at line 66 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4090 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4101 of file cfModGcd.cc.

◆ GG

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

Definition at line 4075 of file cfModGcd.cc.

◆ gp

Definition at line 4102 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4078 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4097 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4097 of file cfModGcd.cc.

◆ maxNumVars

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

Definition at line 4130 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4104 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4129 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4129 of file cfModGcd.cc.

◆ p

int p

Definition at line 4078 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4096 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4082 of file cfModGcd.cc.