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 <math.h>
#include "cf_assert.h"
#include "debug.h"
#include "timing.h"
#include "canonicalform.h"
#include "cfGcdUtil.h"
#include "cf_map.h"
#include "cf_util.h"
#include "cf_irred.h"
#include "templates/ftmpl_functions.h"
#include "cf_random.h"
#include "cf_reval.h"
#include "facHensel.h"
#include "cf_iter.h"
#include "cfNewtonPolygon.h"
#include "cf_algorithm.h"
#include "cf_primes.h"
#include "cf_map_ext.h"
#include "NTLconvert.h"
#include "FLINTconvert.h"
#include "cfModGcd.h"

Go to the source code of this file.

Functions

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

Variables

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

Detailed Description

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

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

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

Definition in file cfModGcd.cc.

Function Documentation

◆ chooseExtension()

static Variable chooseExtension ( const Variable alpha)
inlinestatic

Definition at line 421 of file cfModGcd.cc.

422{
423 int i, m;
424 // extension of F_p needed
425 if (alpha.level() == 1)
426 {
427 i= 1;
428 m= 2;
429 } //extension of F_p(alpha)
430 if (alpha.level() != 1)
431 {
432 i= 4;
433 m= degree (getMipo (alpha));
434 }
435 #ifdef HAVE_FLINT
441 #else
443 {
445 zz_p::init (getCharacteristic());
446 }
450 #endif
451 return rootOf (newMipo);
452}
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
VAR long fac_NTL_char
Definition NTLconvert.cc:46
int degree(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition cf_char.cc:70
int m
Definition cfEzgcd.cc:128
int i
Definition cfModGcd.cc:4079
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 2186 of file cfModGcd.cc.

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

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

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

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

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

◆ for() [1/2]

for ( i,
0;i--   
)

Definition at line 4118 of file cfModGcd.cc.

4119 {
4120 if (degree (f, i) <= 0 || degree (g, i) <= 0)
4121 continue;
4122 else
4124 }
f
Definition cfModGcd.cc:4082
g
Definition cfModGcd.cc:4091
int minCommonDeg
Definition cfModGcd.cc:4105
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)

◆ for() [2/2]

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

Definition at line 4106 of file cfModGcd.cc.

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

◆ FpRandomElement()

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

Definition at line 1172 of file cfModGcd.cc.

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

◆ gaussianElimFp()

long gaussianElimFp ( CFMatrix M,
CFArray L 
)

Definition at line 1740 of file cfModGcd.cc.

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

◆ gaussianElimFq()

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

Definition at line 1785 of file cfModGcd.cc.

1786{
1787 ASSERT (L.size() <= M.rows(), "dimension exceeded");
1788 CFMatrix *N;
1789 N= new CFMatrix (M.rows(), M.columns() + 1);
1790
1791 for (int i= 1; i <= M.rows(); i++)
1792 for (int j= 1; j <= M.columns(); j++)
1793 (*N) (i, j)= M (i, j);
1794
1795 int j= 1;
1796 for (int i= 0; i < L.size(); i++, j++)
1797 (*N) (j, M.columns() + 1)= L[i];
1798 #ifdef HAVE_FLINT
1799 // convert mipo
1805 // convert matrix
1808 // rank
1809 long rk= fq_nmod_mat_rref (FLINTN,ctx);
1810 // clean up
1813 #elif defined(HAVE_NTL)
1814 int p= getCharacteristic ();
1815 if (fac_NTL_char != p)
1816 {
1817 fac_NTL_char= p;
1818 zz_p::init (p);
1819 }
1821 zz_pE::init (NTLMipo);
1823 long rk= gauss (*NTLN);
1825 delete NTLN;
1826 #else
1827 factoryError("NTL/FLINT missing: gaussianElimFq");
1828 #endif
1829 delete N;
1830
1831 M= (*N) (1, M.rows(), 1, M.columns());
1832 L= CFArray (M.rows());
1833 for (int i= 0; i < M.rows(); i++)
1834 L[i]= (*N) (i + 1, M.columns() + 1);
1835
1836 delete N;
1837 return rk;
1838}
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 1958 of file cfModGcd.cc.

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

◆ GFRandomElement()

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

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

Definition at line 820 of file cfModGcd.cc.

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

◆ if() [1/2]

if ( i  = =0)

◆ if() [2/2]

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

Definition at line 72 of file cfModGcd.cc.

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

◆ modGCDFp() [1/2]

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

Definition at line 1213 of file cfModGcd.cc.

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

◆ modGCDFp() [2/2]

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

Definition at line 1224 of file cfModGcd.cc.

1227{
1228 CanonicalForm A= F;
1229 CanonicalForm B= G;
1230 if (F.isZero() && degree(G) > 0)
1231 {
1232 coF= 0;
1233 coG= Lc (G);
1234 return G/Lc(G);
1235 }
1236 else if (G.isZero() && degree (F) > 0)
1237 {
1238 coF= Lc (F);
1239 coG= 0;
1240 return F/Lc(F);
1241 }
1242 else if (F.isZero() && G.isZero())
1243 {
1244 coF= coG= 0;
1245 return F.genOne();
1246 }
1247 if (F.inBaseDomain() || G.inBaseDomain())
1248 {
1249 coF= F;
1250 coG= G;
1251 return F.genOne();
1252 }
1253 if (F.isUnivariate() && fdivides(F, G, coG))
1254 {
1255 coF= Lc (F);
1256 return F/Lc(F);
1257 }
1258 if (G.isUnivariate() && fdivides(G, F, coF))
1259 {
1260 coG= Lc (G);
1261 return G/Lc(G);
1262 }
1263 if (F == G)
1264 {
1265 coF= coG= Lc (F);
1266 return F/Lc(F);
1267 }
1268 CFMap M,N;
1269 int best_level= myCompress (A, B, M, N, topLevel);
1270
1271 if (best_level == 0)
1272 {
1273 coF= F;
1274 coG= G;
1275 return B.genOne();
1276 }
1277
1278 A= M(A);
1279 B= M(B);
1280
1281 Variable x= Variable (1);
1282
1283 //univariate case
1284 if (A.isUnivariate() && B.isUnivariate())
1285 {
1287 coF= N (A/result);
1288 coG= N (B/result);
1289 return N (result);
1290 }
1291
1292 CanonicalForm cA, cB; // content of A and B
1293 CanonicalForm ppA, ppB; // primitive part of A and B
1295
1296 cA = uni_content (A);
1297 cB = uni_content (B);
1298 gcdcAcB= gcd (cA, cB);
1299 ppA= A/cA;
1300 ppB= B/cB;
1301
1302 CanonicalForm lcA, lcB; // leading coefficients of A and B
1304 lcA= uni_lcoeff (ppA);
1305 lcB= uni_lcoeff (ppB);
1306
1307 gcdlcAlcB= gcd (lcA, lcB);
1308
1309 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
1310 int d0;
1311
1312 if (d == 0)
1313 {
1314 coF= N (ppA*(cA/gcdcAcB));
1315 coG= N (ppB*(cB/gcdcAcB));
1316 return N(gcdcAcB);
1317 }
1318
1319 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
1320
1321 if (d0 < d)
1322 d= d0;
1323
1324 if (d == 0)
1325 {
1326 coF= N (ppA*(cA/gcdcAcB));
1327 coG= N (ppB*(cB/gcdcAcB));
1328 return N(gcdcAcB);
1329 }
1330
1334
1335 newtonPoly= 1;
1336 m= gcdlcAlcB;
1337 H= 0;
1338 coF= 0;
1339 coG= 0;
1340 G_m= 0;
1342 bool fail= false;
1343 bool inextension= false;
1344 topLevel= false;
1345 CFList source, dest;
1346 int bound1= degree (ppA, 1);
1347 int bound2= degree (ppB, 1);
1348 do
1349 {
1350 if (inextension)
1352 else
1354
1355 if (!fail && !inextension)
1356 {
1357 CFList list;
1362 list);
1364 "time for recursive call: ");
1365 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1366 }
1367 else if (!fail && inextension)
1368 {
1369 CFList list;
1374 list, topLevel);
1376 "time for recursive call: ");
1377 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1378 }
1379 else if (fail && !inextension)
1380 {
1381 source= CFList();
1382 dest= CFList();
1383 CFList list;
1385 int deg= 2;
1386 bool initialized= false;
1387 do
1388 {
1389 mipo= randomIrredpoly (deg, x);
1390 if (initialized)
1391 setMipo (alpha, mipo);
1392 else
1393 alpha= rootOf (mipo);
1394 inextension= true;
1395 initialized= true;
1396 fail= false;
1398 deg++;
1399 } while (fail);
1400 list= CFList();
1401 V_buf= alpha;
1402 cleanUp= alpha;
1407 list, topLevel);
1409 "time for recursive call: ");
1410 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1411 }
1412 else if (fail && inextension)
1413 {
1414 source= CFList();
1415 dest= CFList();
1416
1419 bool prim_fail= false;
1423
1424 if (V_buf3 != alpha)
1425 {
1431 source, dest);
1435 dest);
1438 for (CFListIterator i= l; i.hasItem(); i++)
1439 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
1440 source, dest);
1441 }
1442
1443 ASSERT (!prim_fail, "failure in integer factorizer");
1444 if (prim_fail)
1445 ; //ERROR
1446 else
1448
1449 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1450 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1451
1452 for (CFListIterator i= l; i.hasItem(); i++)
1453 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
1454 im_prim_elem, source, dest);
1460 source, dest);
1464 source, dest);
1467 fail= false;
1469 DEBOUTLN (cerr, "fail= " << fail);
1470 CFList list;
1475 list, topLevel);
1477 "time for recursive call: ");
1478 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1479 alpha= V_buf;
1480 }
1481
1482 if (!G_random_element.inCoeffDomain())
1484 Variable (G_random_element.level()));
1485 else
1486 d0= 0;
1487
1488 if (d0 == 0)
1489 {
1490 if (inextension)
1491 prune (cleanUp);
1492 coF= N (ppA*(cA/gcdcAcB));
1493 coG= N (ppB*(cB/gcdcAcB));
1494 return N(gcdcAcB);
1495 }
1496
1497 if (d0 > d)
1498 {
1499 if (!find (l, random_element))
1501 continue;
1502 }
1503
1506
1511
1512 if (!G_random_element.inCoeffDomain())
1514 Variable (G_random_element.level()));
1515 else
1516 d0= 0;
1517
1518 if (d0 < d)
1519 {
1520 m= gcdlcAlcB;
1521 newtonPoly= 1;
1522 G_m= 0;
1523 d= d0;
1524 coF_m= 0;
1525 coG_m= 0;
1526 }
1527
1533 "time for newton_interpolation: ");
1534
1535 //termination test
1536 if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1537 {
1539 if (gcdlcAlcB.isOne())
1540 cH= 1;
1541 else
1542 cH= uni_content (H);
1543 ppH= H/cH;
1544 ppH /= Lc (ppH);
1548 ppCoF= coF/ccoF;
1549 ppCoG= coG/ccoG;
1550 DEBOUTLN (cerr, "ppH= " << ppH);
1551 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1552 (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1554 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1555 {
1556 if (inextension)
1557 prune (cleanUp);
1558 coF= N ((cA/gcdcAcB)*ppCoF);
1559 coG= N ((cB/gcdcAcB)*ppCoG);
1561 "time for successful termination Fp: ");
1562 return N(gcdcAcB*ppH);
1563 }
1565 "time for unsuccessful termination Fp: ");
1566 }
1567
1568 G_m= H;
1569 coF_m= coF;
1570 coG_m= coG;
1572 m= m*(x - random_element);
1573 if (!find (l, random_element))
1575 } while (1);
1576}
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition cf_ops.cc:523
CanonicalForm Lc(const CanonicalForm &f)
CanonicalForm modGCDFq(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, Variable &alpha, CFList &l, bool &topLevel)
GCD of F and G over , l and topLevel are only used internally, output is monic based on Alg....
Definition cfModGcd.cc:479
int myCompress(const CanonicalForm &F, const CanonicalForm &G, CFMap &M, CFMap &N, bool topLevel)
compressing two polynomials F and G, M is used for compressing, N to reverse the compression
Definition cfModGcd.cc:92
static Variable chooseExtension(const Variable &alpha)
Definition cfModGcd.cc:421
static CanonicalForm uni_lcoeff(const CanonicalForm &F)
compute the leading coefficient of F, where F is considered as an element of , order on is dp.
Definition cfModGcd.cc:340
static CanonicalForm randomElement(const CanonicalForm &F, const Variable &alpha, CFList &list, bool &fail)
compute a random element a of , s.t. F(a) , F is a univariate polynomial, returns fail if there are...
Definition cfModGcd.cc:380
static CanonicalForm FpRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
Definition cfModGcd.cc:1172
static CanonicalForm newtonInterp(const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x)
Newton interpolation - Incremental algorithm. Given a list of values alpha_i and a list of polynomial...
Definition cfModGcd.cc:365
bool terminationTest(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &coF, const CanonicalForm &coG, const CanonicalForm &cand)
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
CanonicalForm randomIrredpoly(int i, const Variable &x)
computes a random monic irreducible univariate polynomial in x over Fp of degree i via NTL/FLINT
Definition cf_irred.cc:26
CanonicalForm mapPrimElem(const CanonicalForm &primElem, const Variable &alpha, const Variable &beta)
compute the image of a primitive element of in . We assume .
CanonicalForm primitiveElement(const Variable &alpha, Variable &beta, bool &fail)
determine a primitive element of , is a primitive element of a field which is isomorphic to
static CanonicalForm mapDown(const CanonicalForm &F, const Variable &alpha, const CanonicalForm &G, CFList &source, CFList &dest)
the CanonicalForm G is the output of map_up, returns F considered as an element over ,...
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
Definition cf_map_ext.cc:71
class CFMap
Definition cf_map.h:85
CanonicalForm genOne() const
CF_NO_INLINE bool 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 479 of file cfModGcd.cc.

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

◆ modGCDFq() [2/2]

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

Definition at line 463 of file cfModGcd.cc.

465{
468 topLevel);
469 return result;
470}

◆ modGCDGF() [1/2]

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

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

Definition at line 873 of file cfModGcd.cc.

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

◆ modGCDGF() [2/2]

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

Definition at line 859 of file cfModGcd.cc.

861{
864 return result;
865}

◆ monicSparseInterpol()

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

Definition at line 2200 of file cfModGcd.cc.

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

◆ mult()

void mult ( CFList L1,
const CFList L2 
)

multiply two lists componentwise

Definition at line 2177 of file cfModGcd.cc.

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

◆ myCompress()

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

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

Definition at line 92 of file cfModGcd.cc.

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

◆ newtonInterp()

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

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

Definition at line 365 of file cfModGcd.cc.

368{
370
372 *newtonPoly;
373 return interPoly;
374}

◆ nonMonicSparseInterpol()

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

Definition at line 2484 of file cfModGcd.cc.

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

◆ randomElement()

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

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

Definition at line 380 of file cfModGcd.cc.

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

◆ readOffSolution() [1/2]

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

Definition at line 1711 of file cfModGcd.cc.

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

◆ readOffSolution() [2/2]

CFArray readOffSolution ( const CFMatrix M,
const long  rk 
)

M in row echolon form, rk rank of M.

Definition at line 1689 of file cfModGcd.cc.

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

◆ solveGeneralVandermonde()

CFArray solveGeneralVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1633 of file cfModGcd.cc.

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

◆ solveSystemFp()

CFArray solveSystemFp ( const CFMatrix M,
const CFArray L 
)

Definition at line 1841 of file cfModGcd.cc.

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

◆ solveSystemFq()

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

Definition at line 1893 of file cfModGcd.cc.

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

◆ solveVandermonde()

CFArray solveVandermonde ( const CFArray M,
const CFArray A 
)

Definition at line 1580 of file cfModGcd.cc.

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

◆ sparseGCDFp()

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

Definition at line 3563 of file cfModGcd.cc.

3565{
3566 CanonicalForm A= F;
3567 CanonicalForm B= G;
3568 if (F.isZero() && degree(G) > 0) return G/Lc(G);
3569 else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3570 else if (F.isZero() && G.isZero()) return F.genOne();
3571 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3572 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3573 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3574 if (F == G) return F/Lc(F);
3575
3576 CFMap M,N;
3577 int best_level= myCompress (A, B, M, N, topLevel);
3578
3579 if (best_level == 0) return B.genOne();
3580
3581 A= M(A);
3582 B= M(B);
3583
3584 Variable x= Variable (1);
3585
3586 //univariate case
3587 if (A.isUnivariate() && B.isUnivariate())
3588 return N (gcd (A, B));
3589
3590 CanonicalForm cA, cB; // content of A and B
3591 CanonicalForm ppA, ppB; // primitive part of A and B
3593
3594 cA = uni_content (A);
3595 cB = uni_content (B);
3596 gcdcAcB= gcd (cA, cB);
3597 ppA= A/cA;
3598 ppB= B/cB;
3599
3600 CanonicalForm lcA, lcB; // leading coefficients of A and B
3602 lcA= uni_lcoeff (ppA);
3603 lcB= uni_lcoeff (ppB);
3604
3605 if (fdivides (lcA, lcB))
3606 {
3607 if (fdivides (A, B))
3608 return F/Lc(F);
3609 }
3610 if (fdivides (lcB, lcA))
3611 {
3612 if (fdivides (B, A))
3613 return G/Lc(G);
3614 }
3615
3616 gcdlcAlcB= gcd (lcA, lcB);
3617
3618 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3619 int d0;
3620
3621 if (d == 0)
3622 return N(gcdcAcB);
3623
3624 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3625
3626 if (d0 < d)
3627 d= d0;
3628
3629 if (d == 0)
3630 return N(gcdcAcB);
3631
3634 m= gcdlcAlcB;
3635 G_m= 0;
3636 H= 0;
3637 bool fail= false;
3638 topLevel= false;
3639 bool inextension= false;
3642 CFList source, dest;
3643 do //first do
3644 {
3645 if (inextension)
3647 else
3649 if (random_element == 0 && !fail)
3650 {
3651 if (!find (l, random_element))
3653 continue;
3654 }
3655
3656 if (!fail && !inextension)
3657 {
3658 CFList list;
3662 list);
3664 "time for recursive call: ");
3665 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3666 }
3667 else if (!fail && inextension)
3668 {
3669 CFList list;
3673 list, topLevel);
3675 "time for recursive call: ");
3676 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3677 }
3678 else if (fail && !inextension)
3679 {
3680 source= CFList();
3681 dest= CFList();
3682 CFList list;
3684 int deg= 2;
3685 bool initialized= false;
3686 do
3687 {
3688 mipo= randomIrredpoly (deg, x);
3689 if (initialized)
3690 setMipo (alpha, mipo);
3691 else
3692 alpha= rootOf (mipo);
3693 inextension= true;
3694 fail= false;
3695 initialized= true;
3697 deg++;
3698 } while (fail);
3699 cleanUp= alpha;
3700 V_buf= alpha;
3701 list= CFList();
3705 list, topLevel);
3707 "time for recursive call: ");
3708 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3709 }
3710 else if (fail && inextension)
3711 {
3712 source= CFList();
3713 dest= CFList();
3714
3717 bool prim_fail= false;
3721
3722 if (V_buf3 != alpha)
3723 {
3727 dest);
3731 dest);
3732 for (CFListIterator i= l; i.hasItem(); i++)
3733 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3734 source, dest);
3735 }
3736
3737 ASSERT (!prim_fail, "failure in integer factorizer");
3738 if (prim_fail)
3739 ; //ERROR
3740 else
3742
3743 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3744 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3745
3746 for (CFListIterator i= l; i.hasItem(); i++)
3747 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3748 im_prim_elem, source, dest);
3752 source, dest);
3756 source, dest);
3757 fail= false;
3759 DEBOUTLN (cerr, "fail= " << fail);
3760 CFList list;
3764 list, topLevel);
3766 "time for recursive call: ");
3767 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3768 alpha= V_buf;
3769 }
3770
3771 if (!G_random_element.inCoeffDomain())
3773 Variable (G_random_element.level()));
3774 else
3775 d0= 0;
3776
3777 if (d0 == 0)
3778 {
3779 if (inextension)
3780 prune (cleanUp);
3781 return N(gcdcAcB);
3782 }
3783 if (d0 > d)
3784 {
3785 if (!find (l, random_element))
3787 continue;
3788 }
3789
3793
3795
3796 if (!G_random_element.inCoeffDomain())
3798 Variable (G_random_element.level()));
3799 else
3800 d0= 0;
3801
3802 if (d0 < d)
3803 {
3804 m= gcdlcAlcB;
3805 newtonPoly= 1;
3806 G_m= 0;
3807 d= d0;
3808 }
3809
3811
3812 if (uni_lcoeff (H) == gcdlcAlcB)
3813 {
3814 cH= uni_content (H);
3815 ppH= H/cH;
3816 ppH /= Lc (ppH);
3817 DEBOUTLN (cerr, "ppH= " << ppH);
3818
3819 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3820 {
3821 if (inextension)
3822 prune (cleanUp);
3823 return N(gcdcAcB*ppH);
3824 }
3825 }
3826 G_m= H;
3828 m= m*(x - random_element);
3829 if (!find (l, random_element))
3831
3832 if ((getCharacteristic() > 3 && size (skeleton) < 200))
3833 {
3836
3837 do //second do
3838 {
3839 if (inextension)
3841 else
3843 if (random_element == 0 && !fail)
3844 {
3845 if (!find (l, random_element))
3847 continue;
3848 }
3849
3850 bool sparseFail= false;
3851 if (!fail && !inextension)
3852 {
3853 CFList list;
3855 if (LC (skeleton).inCoeffDomain())
3859 Monoms);
3860 else
3866 "time for recursive call: ");
3867 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3868 }
3869 else if (!fail && inextension)
3870 {
3871 CFList list;
3873 if (LC (skeleton).inCoeffDomain())
3877 Monoms);
3878 else
3882 Monoms);
3884 "time for recursive call: ");
3885 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3886 }
3887 else if (fail && !inextension)
3888 {
3889 source= CFList();
3890 dest= CFList();
3891 CFList list;
3893 int deg= 2;
3894 bool initialized= false;
3895 do
3896 {
3897 mipo= randomIrredpoly (deg, x);
3898 if (initialized)
3899 setMipo (alpha, mipo);
3900 else
3901 alpha= rootOf (mipo);
3902 inextension= true;
3903 fail= false;
3904 initialized= true;
3906 deg++;
3907 } while (fail);
3908 cleanUp= alpha;
3909 V_buf= alpha;
3910 list= CFList();
3912 if (LC (skeleton).inCoeffDomain())
3916 Monoms);
3917 else
3921 Monoms);
3923 "time for recursive call: ");
3924 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3925 }
3926 else if (fail && inextension)
3927 {
3928 source= CFList();
3929 dest= CFList();
3930
3933 bool prim_fail= false;
3937
3938 if (V_buf3 != alpha)
3939 {
3943 source, dest);
3947 source, dest);
3948 for (CFListIterator i= l; i.hasItem(); i++)
3949 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3950 source, dest);
3951 }
3952
3953 ASSERT (!prim_fail, "failure in integer factorizer");
3954 if (prim_fail)
3955 ; //ERROR
3956 else
3958
3959 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3960 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3961
3962 for (CFListIterator i= l; i.hasItem(); i++)
3963 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3964 im_prim_elem, source, dest);
3968 source, dest);
3972 source, dest);
3973 fail= false;
3975 DEBOUTLN (cerr, "fail= " << fail);
3976 CFList list;
3978 if (LC (skeleton).inCoeffDomain())
3982 Monoms);
3983 else
3989 "time for recursive call: ");
3990 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3991 alpha= V_buf;
3992 }
3993
3994 if (sparseFail)
3995 break;
3996
3997 if (!G_random_element.inCoeffDomain())
3999 Variable (G_random_element.level()));
4000 else
4001 d0= 0;
4002
4003 if (d0 == 0)
4004 {
4005 if (inextension)
4006 prune (cleanUp);
4007 return N(gcdcAcB);
4008 }
4009 if (d0 > d)
4010 {
4011 if (!find (l, random_element))
4013 continue;
4014 }
4015
4019
4020 if (!G_random_element.inCoeffDomain())
4022 Variable (G_random_element.level()));
4023 else
4024 d0= 0;
4025
4026 if (d0 < d)
4027 {
4028 m= gcdlcAlcB;
4029 newtonPoly= 1;
4030 G_m= 0;
4031 d= d0;
4032 }
4033
4037 "time for newton interpolation: ");
4038
4039 //termination test
4040 if (uni_lcoeff (H) == gcdlcAlcB)
4041 {
4042 cH= uni_content (H);
4043 ppH= H/cH;
4044 ppH /= Lc (ppH);
4045 DEBOUTLN (cerr, "ppH= " << ppH);
4046 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4047 {
4048 if (inextension)
4049 prune (cleanUp);
4050 return N(gcdcAcB*ppH);
4051 }
4052 }
4053
4054 G_m= H;
4056 m= m*(x - random_element);
4057 if (!find (l, random_element))
4059
4060 } while (1); //end of second do
4061 }
4062 else
4063 {
4064 if (inextension)
4065 prune (cleanUp);
4066 return N(gcdcAcB*modGCDFp (ppA, ppB));
4067 }
4068 } while (1); //end of first do
4069}
CanonicalForm sparseGCDFq(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
Definition cfModGcd.cc:3131
CanonicalForm monicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition cfModGcd.cc:2200
CanonicalForm nonMonicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition cfModGcd.cc:2484
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
Definition cfModGcd.cc:3563

◆ sparseGCDFq()

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

Definition at line 3131 of file cfModGcd.cc.

3133{
3134 CanonicalForm A= F;
3135 CanonicalForm B= G;
3136 if (F.isZero() && degree(G) > 0) return G/Lc(G);
3137 else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3138 else if (F.isZero() && G.isZero()) return F.genOne();
3139 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3140 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3141 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3142 if (F == G) return F/Lc(F);
3143
3144 CFMap M,N;
3145 int best_level= myCompress (A, B, M, N, topLevel);
3146
3147 if (best_level == 0) return B.genOne();
3148
3149 A= M(A);
3150 B= M(B);
3151
3152 Variable x= Variable (1);
3153
3154 //univariate case
3155 if (A.isUnivariate() && B.isUnivariate())
3156 return N (gcd (A, B));
3157
3158 CanonicalForm cA, cB; // content of A and B
3159 CanonicalForm ppA, ppB; // primitive part of A and B
3161
3162 cA = uni_content (A);
3163 cB = uni_content (B);
3164 gcdcAcB= gcd (cA, cB);
3165 ppA= A/cA;
3166 ppB= B/cB;
3167
3168 CanonicalForm lcA, lcB; // leading coefficients of A and B
3170 lcA= uni_lcoeff (ppA);
3171 lcB= uni_lcoeff (ppB);
3172
3173 if (fdivides (lcA, lcB))
3174 {
3175 if (fdivides (A, B))
3176 return F/Lc(F);
3177 }
3178 if (fdivides (lcB, lcA))
3179 {
3180 if (fdivides (B, A))
3181 return G/Lc(G);
3182 }
3183
3184 gcdlcAlcB= gcd (lcA, lcB);
3185
3186 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3187 int d0;
3188
3189 if (d == 0)
3190 return N(gcdcAcB);
3191 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3192
3193 if (d0 < d)
3194 d= d0;
3195
3196 if (d == 0)
3197 return N(gcdcAcB);
3198
3201 m= gcdlcAlcB;
3202 G_m= 0;
3203 H= 0;
3204 bool fail= false;
3205 topLevel= false;
3206 bool inextension= false;
3210 CFList source, dest;
3211 do // first do
3212 {
3214 if (random_element == 0 && !fail)
3215 {
3216 if (!find (l, random_element))
3218 continue;
3219 }
3220 if (fail)
3221 {
3222 source= CFList();
3223 dest= CFList();
3224
3227 bool prim_fail= false;
3230 if (V_buf4 == alpha)
3232
3233 if (V_buf3 != V_buf4)
3234 {
3238 source, dest);
3242 dest);
3243 for (CFListIterator i= l; i.hasItem(); i++)
3244 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3245 source, dest);
3246 }
3247
3248 ASSERT (!prim_fail, "failure in integer factorizer");
3249 if (prim_fail)
3250 ; //ERROR
3251 else
3253
3254 if (V_buf4 == alpha)
3256 else
3258 im_prim_elem, source, dest);
3259
3260 DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3261 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3262 inextension= true;
3263 for (CFListIterator i= l; i.hasItem(); i++)
3264 i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3265 im_prim_elem, source, dest);
3269 source, dest);
3273 source, dest);
3274
3275 fail= false;
3277 DEBOUTLN (cerr, "fail= " << fail);
3278 CFList list;
3282 list, topLevel);
3284 "time for recursive call: ");
3285 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3286 V_buf4= V_buf;
3287 }
3288 else
3289 {
3290 CFList list;
3294 list, topLevel);
3296 "time for recursive call: ");
3297 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3298 }
3299
3300 if (!G_random_element.inCoeffDomain())
3302 Variable (G_random_element.level()));
3303 else
3304 d0= 0;
3305
3306 if (d0 == 0)
3307 {
3308 if (inextension)
3309 prune1 (alpha);
3310 return N(gcdcAcB);
3311 }
3312 if (d0 > d)
3313 {
3314 if (!find (l, random_element))
3316 continue;
3317 }
3318
3322
3324 if (!G_random_element.inCoeffDomain())
3326 Variable (G_random_element.level()));
3327 else
3328 d0= 0;
3329
3330 if (d0 < d)
3331 {
3332 m= gcdlcAlcB;
3333 newtonPoly= 1;
3334 G_m= 0;
3335 d= d0;
3336 }
3337
3339 if (uni_lcoeff (H) == gcdlcAlcB)
3340 {
3341 cH= uni_content (H);
3342 ppH= H/cH;
3343 if (inextension)
3344 {
3345 CFList u, v;
3346 //maybe it's better to test if ppH is an element of F(\alpha) before
3347 //mapping down
3348 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3349 {
3350 DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3352 ppH /= Lc(ppH);
3353 DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3354 prune1 (alpha);
3355 return N(gcdcAcB*ppH);
3356 }
3357 }
3358 else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3359 return N(gcdcAcB*ppH);
3360 }
3361 G_m= H;
3363 m= m*(x - random_element);
3364 if (!find (l, random_element))
3366
3367 if (getCharacteristic () > 3 && size (skeleton) < 100)
3368 {
3371 do //second do
3372 {
3374 if (random_element == 0 && !fail)
3375 {
3376 if (!find (l, random_element))
3378 continue;
3379 }
3380 if (fail)
3381 {
3382 source= CFList();
3383 dest= CFList();
3384
3387 bool prim_fail= false;
3390 if (V_buf4 == alpha)
3392
3393 if (V_buf3 != V_buf4)
3394 {
3398 source, dest);
3402 source, dest);
3403 for (CFListIterator i= l; i.hasItem(); i++)
3404 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3405 source, dest);
3406 }
3407
3408 ASSERT (!prim_fail, "failure in integer factorizer");
3409 if (prim_fail)
3410 ; //ERROR
3411 else
3413
3414 if (V_buf4 == alpha)
3416 else
3418 prim_elem, im_prim_elem, source, dest);
3419
3420 DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3421 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3422 inextension= true;
3423 for (CFListIterator i= l; i.hasItem(); i++)
3424 i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3425 im_prim_elem, source, dest);
3429 source, dest);
3432
3434 source, dest);
3435
3436 fail= false;
3438 DEBOUTLN (cerr, "fail= " << fail);
3439 CFList list;
3441
3442 V_buf4= V_buf;
3443
3444 //sparseInterpolation
3445 bool sparseFail= false;
3446 if (LC (skeleton).inCoeffDomain())
3450 else
3454 Monoms);
3455 if (sparseFail)
3456 break;
3457
3459 "time for recursive call: ");
3460 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3461 }
3462 else
3463 {
3464 CFList list;
3466 bool sparseFail= false;
3467 if (LC (skeleton).inCoeffDomain())
3471 else
3475 Monoms);
3476 if (sparseFail)
3477 break;
3478
3480 "time for recursive call: ");
3481 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3482 }
3483
3484 if (!G_random_element.inCoeffDomain())
3486 Variable (G_random_element.level()));
3487 else
3488 d0= 0;
3489
3490 if (d0 == 0)
3491 {
3492 if (inextension)
3493 prune1 (alpha);
3494 return N(gcdcAcB);
3495 }
3496 if (d0 > d)
3497 {
3498 if (!find (l, random_element))
3500 continue;
3501 }
3502
3506
3507 if (!G_random_element.inCoeffDomain())
3509 Variable (G_random_element.level()));
3510 else
3511 d0= 0;
3512
3513 if (d0 < d)
3514 {
3515 m= gcdlcAlcB;
3516 newtonPoly= 1;
3517 G_m= 0;
3518 d= d0;
3519 }
3520
3524 "time for newton interpolation: ");
3525
3526 //termination test
3527 if (uni_lcoeff (H) == gcdlcAlcB)
3528 {
3529 cH= uni_content (H);
3530 ppH= H/cH;
3531 if (inextension)
3532 {
3533 CFList u, v;
3534 //maybe it's better to test if ppH is an element of F(\alpha) before
3535 //mapping down
3536 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3537 {
3538 DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3540 ppH /= Lc(ppH);
3541 DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3542 prune1 (alpha);
3543 return N(gcdcAcB*ppH);
3544 }
3545 }
3546 else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3547 {
3548 return N(gcdcAcB*ppH);
3549 }
3550 }
3551
3552 G_m= H;
3554 m= m*(x - random_element);
3555 if (!find (l, random_element))
3557
3558 } while (1);
3559 }
3560 } while (1);
3561}

◆ TIMING_DEFINE_PRINT() [1/2]

TIMING_DEFINE_PRINT ( gcd_recursion  ) const &

◆ TIMING_DEFINE_PRINT() [2/2]

TIMING_DEFINE_PRINT ( modZ_termination  ) const &

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

◆ uni_content() [1/2]

static CanonicalForm uni_content ( const CanonicalForm F)
inlinestatic

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

Definition at line 282 of file cfModGcd.cc.

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

◆ uni_content() [2/2]

CanonicalForm uni_content ( const CanonicalForm F,
const Variable x 
)

Definition at line 260 of file cfModGcd.cc.

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

◆ uni_lcoeff()

static CanonicalForm uni_lcoeff ( const CanonicalForm F)
inlinestatic

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

Definition at line 340 of file cfModGcd.cc.

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

◆ while()

while ( true  )

Definition at line 4133 of file cfModGcd.cc.

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

◆ cand

Initial value:

Definition at line 69 of file cfModGcd.cc.

◆ cd

cd ( bCommonDen(FF ) = bCommonDen( GG )

Definition at line 4090 of file cfModGcd.cc.

◆ cDn

Definition at line 4130 of file cfModGcd.cc.

◆ cf

cf = icontent (f)

Definition at line 4084 of file cfModGcd.cc.

◆ cg

cg = icontent (g)

Definition at line 4084 of file cfModGcd.cc.

◆ cl

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

Definition at line 4101 of file cfModGcd.cc.

◆ coF

Definition at line 68 of file cfModGcd.cc.

◆ cof

Definition at line 4130 of file cfModGcd.cc.

◆ cofn

Definition at line 4130 of file cfModGcd.cc.

◆ cofp

Definition at line 4130 of file cfModGcd.cc.

◆ coG

Definition at line 68 of file cfModGcd.cc.

◆ cog

Definition at line 4130 of file cfModGcd.cc.

◆ cogn

Definition at line 4130 of file cfModGcd.cc.

◆ cogp

Definition at line 4130 of file cfModGcd.cc.

◆ d_deg

int d_deg =-1

Definition at line 4079 of file cfModGcd.cc.

◆ Dn

Definition at line 4097 of file cfModGcd.cc.

◆ dp_deg

int dp_deg

Definition at line 4079 of file cfModGcd.cc.

◆ equal

bool equal = false

Definition at line 4127 of file cfModGcd.cc.

◆ f

f =cd*FF

Definition at line 4082 of file cfModGcd.cc.

◆ false

return false

Definition at line 85 of file cfModGcd.cc.

◆ fp

Definition at line 4103 of file cfModGcd.cc.

◆ G

Definition at line 67 of file cfModGcd.cc.

◆ g

g =cd*GG

Definition at line 4091 of file cfModGcd.cc.

◆ gcdcfcg

CanonicalForm gcdcfcg = gcd (cf, cg)

Definition at line 4102 of file cfModGcd.cc.

◆ GG

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

Definition at line 4076 of file cfModGcd.cc.

◆ gp

Definition at line 4103 of file cfModGcd.cc.

◆ i

i = cf_getNumBigPrimes() - 1

Definition at line 4079 of file cfModGcd.cc.

◆ lcf

lcf = f.lc()

Definition at line 4098 of file cfModGcd.cc.

◆ lcg

lcg = g.lc()

Definition at line 4098 of file cfModGcd.cc.

◆ maxNumVars

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

Definition at line 4131 of file cfModGcd.cc.

◆ minCommonDeg

int minCommonDeg = 0

Definition at line 4105 of file cfModGcd.cc.

◆ newCof

CanonicalForm newCof

Definition at line 4130 of file cfModGcd.cc.

◆ newCog

CanonicalForm newCog

Definition at line 4130 of file cfModGcd.cc.

◆ p

int p

Definition at line 4079 of file cfModGcd.cc.

◆ test

CanonicalForm test = 0

Definition at line 4097 of file cfModGcd.cc.

◆ x

Variable x = Variable (1)

Definition at line 4083 of file cfModGcd.cc.