 |
My Project
UNKNOWN_GIT_VERSION
|
#include "config.h"
#include "cf_assert.h"
#include "debug.h"
#include "timing.h"
#include "canonicalform.h"
#include "cfGcdUtil.h"
#include "cf_map.h"
#include "cf_util.h"
#include "cf_irred.h"
#include "templates/ftmpl_functions.h"
#include "cf_random.h"
#include "cf_reval.h"
#include "facHensel.h"
#include "cf_iter.h"
#include "cfNewtonPolygon.h"
#include "cf_algorithm.h"
#include "cf_primes.h"
#include "cf_map_ext.h"
#include "NTLconvert.h"
#include "FLINTconvert.h"
#include "cfModGcd.h"
Go to the source code of this file.
|
| TIMING_DEFINE_PRINT (gcd_recursion) TIMING_DEFINE_PRINT(newton_interpolation) TIMING_DEFINE_PRINT(termination_test) TIMING_DEFINE_PRINT(ez_p_compress) TIMING_DEFINE_PRINT(ez_p_hensel_lift) TIMING_DEFINE_PRINT(ez_p_eval) TIMING_DEFINE_PRINT(ez_p_content) bool terminationTest(const CanonicalForm &F |
|
| if (LCCand *abs(LC(coF))==abs(LC(F))) |
|
int | myCompress (const CanonicalForm &F, const CanonicalForm &G, CFMap &M, CFMap &N, bool topLevel) |
| compressing two polynomials F and G, M is used for compressing, N to reverse the compression More...
|
|
static CanonicalForm | uni_content (const CanonicalForm &F) |
| compute the content of F, where F is considered as an element of More...
|
|
CanonicalForm | uni_content (const CanonicalForm &F, const Variable &x) |
|
CanonicalForm | extractContents (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &contentF, CanonicalForm &contentG, CanonicalForm &ppF, CanonicalForm &ppG, const int d) |
|
static CanonicalForm | uni_lcoeff (const CanonicalForm &F) |
| compute the leading coefficient of F, where F is considered as an element of , order on is dp. More...
|
|
static CanonicalForm | newtonInterp (const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x) |
| Newton interpolation - Incremental algorithm. Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n, computes the interpolation polynomial assuming that the polynomials in u are the results of evaluating the variabe x of the unknown polynomial at the alpha values. This incremental version receives only the values of alpha_n and u_n and the previous interpolation polynomial for points 1 <= i <= n-1, and computes the polynomial interpolating in all the points. newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1}) More...
|
|
static CanonicalForm | randomElement (const CanonicalForm &F, const Variable &alpha, CFList &list, bool &fail) |
| compute a random element a of , s.t. F(a) , F is a univariate polynomial, returns fail if there are no field elements left which have not been used before More...
|
|
static Variable | chooseExtension (const Variable &alpha) |
|
CanonicalForm | modGCDFq (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, Variable &alpha, CFList &l, bool &topLevel) |
| GCD of F and G over , l and topLevel are only used internally, output is monic based on Alg. 7.2. as described in "Algorithms for
Computer Algebra" by Geddes, Czapor, Labahn. More...
|
|
CanonicalForm | modGCDFq (const CanonicalForm &F, const CanonicalForm &G, Variable &alpha, CFList &l, bool &topLevel) |
|
static CanonicalForm | GFRandomElement (const CanonicalForm &F, CFList &list, bool &fail) |
| compute a random element a of GF, s.t. F(a) , F is a univariate polynomial, returns fail if there are no field elements left which have not been used before More...
|
|
CanonicalForm | modGCDGF (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, CFList &l, bool &topLevel) |
| GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
Computer Algebra" by Geddes, Czapor, Labahn Usually this algorithm will be faster than modGCDFq since GF has faster field arithmetics, however it might fail if the input is large since the size of the base field is bounded by 2^16, output is monic. More...
|
|
CanonicalForm | modGCDGF (const CanonicalForm &F, const CanonicalForm &G, CFList &l, bool &topLevel) |
|
static CanonicalForm | FpRandomElement (const CanonicalForm &F, CFList &list, bool &fail) |
|
CanonicalForm | modGCDFp (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l) |
|
CanonicalForm | modGCDFp (const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l) |
|
CFArray | solveVandermonde (const CFArray &M, const CFArray &A) |
|
CFArray | solveGeneralVandermonde (const CFArray &M, const CFArray &A) |
|
CFArray | readOffSolution (const CFMatrix &M, const long rk) |
| M in row echolon form, rk rank of M. More...
|
|
CFArray | readOffSolution (const CFMatrix &M, const CFArray &L, const CFArray &partialSol) |
|
long | gaussianElimFp (CFMatrix &M, CFArray &L) |
|
long | gaussianElimFq (CFMatrix &M, CFArray &L, const Variable &alpha) |
|
CFArray | solveSystemFp (const CFMatrix &M, const CFArray &L) |
|
CFArray | solveSystemFq (const CFMatrix &M, const CFArray &L, const Variable &alpha) |
|
CFArray | getMonoms (const CanonicalForm &F) |
| extract monomials of F, parts in algebraic variable are considered coefficients More...
|
|
CFArray | evaluateMonom (const CanonicalForm &F, const CFList &evalPoints) |
|
CFArray | evaluate (const CFArray &A, const CFList &evalPoints) |
|
CFList | evaluationPoints (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Feval, CanonicalForm &Geval, const CanonicalForm &LCF, const bool &GF, const Variable &alpha, bool &fail, CFList &list) |
|
void | mult (CFList &L1, const CFList &L2) |
| multiply two lists componentwise More...
|
|
void | eval (const CanonicalForm &A, const CanonicalForm &B, CanonicalForm &Aeval, CanonicalForm &Beval, const CFList &L) |
|
CanonicalForm | monicSparseInterpol (const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms) |
|
CanonicalForm | nonMonicSparseInterpol (const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms) |
|
CanonicalForm | sparseGCDFq (const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel) |
|
CanonicalForm | sparseGCDFp (const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l) |
|
| TIMING_DEFINE_PRINT (modZ_termination) TIMING_DEFINE_PRINT(modZ_recursion) CanonicalForm modGCDZ(const CanonicalForm &FF |
| modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer
Algebra", Algorithm 7.1 More...
|
|
| for (i=tmax(f.level(), g.level());i > 0;i--) |
|
| if (i==0) return gcdcfcg |
|
| for (;i > 0;i--) |
|
| while (true) |
|
This file implements the GCD of two polynomials over
,
, 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.
◆ chooseExtension()
◆ eval()
◆ evaluate()
Definition at line 1972 of file cfModGcd.cc.
1977 for (
int i= 0;
i <
A.size();
i++)
1982 tmp= tmp (
j.getItem(),
k);
◆ evaluateMonom()
Definition at line 1933 of file cfModGcd.cc.
1944 "expected an eval point with only one component");
1952 int numMon=
size (F);
1964 for (
int k= 0;
k < recResult.
size();
k++)
1966 j += recResult.
size();
◆ evaluationPoints()
Definition at line 1989 of file cfModGcd.cc.
2017 bool zeroOneOccured=
false;
2018 bool allEqual=
false;
2029 for (
int i= 0;
i <
k;
i++)
2039 result.append (genAlgExt.generate());
2047 if (
result.getLast().isOne() ||
result.getLast().isZero())
2048 zeroOneOccured=
true;
2050 if (
find (list, random))
2052 zeroOneOccured=
false;
2060 zeroOneOccured=
false;
2080 zeroOneOccured=
false;
2092 Geval= Geval (
i.getItem(),
j);
2093 LCeval= LCeval (
i.getItem(),
j);
2098 if (!
find (list, random))
2100 zeroOneOccured=
false;
2111 }
while (
find (list, random));
◆ extractContents()
Definition at line 313 of file cfModGcd.cc.
323 for (
int i= 1;
i <= d;
i++)
327 gcdcFcG=
gcd (uniContentF, uniContentG);
328 contentF *= uniContentF;
329 contentG *= uniContentG;
◆ for() [1/2]
◆ for() [2/2]
◆ FpRandomElement()
Definition at line 1159 of file cfModGcd.cc.
1179 while (
find (list, random))
1182 if (F (random,
x) == 0)
1187 }
while (
find (list, random));
◆ gaussianElimFp()
Definition at line 1721 of file cfModGcd.cc.
1723 ASSERT (L.
size() <=
M.rows(),
"dimension exceeded");
1727 for (
int i= 1;
i <=
M.rows();
i++)
1728 for (
int j= 1;
j <=
M.columns();
j++)
1732 for (
int i= 0;
i < L.
size();
i++,
j++)
1733 (*
N) (
j,
M.columns() + 1)= L[
i];
1737 long rk= nmod_mat_rref (FLINTN);
1741 nmod_mat_clear (FLINTN);
1751 long rk= gauss (*NTLN);
1758 for (
int i= 0;
i <
M.rows();
i++)
1759 L[
i]= (*
N) (
i + 1,
M.columns() + 1);
1760 M= (*N) (1,
M.rows(), 1,
M.columns());
◆ gaussianElimFq()
Definition at line 1766 of file cfModGcd.cc.
1768 ASSERT (L.
size() <=
M.rows(),
"dimension exceeded");
1772 for (
int i= 1;
i <=
M.rows();
i++)
1773 for (
int j= 1;
j <=
M.columns();
j++)
1777 for (
int i= 0;
i < L.
size();
i++,
j++)
1778 (*
N) (
j,
M.columns() + 1)= L[
i];
1786 zz_pE::init (NTLMipo);
1788 long rk= gauss (*NTLN);
1795 M= (*N) (1,
M.rows(), 1,
M.columns());
1797 for (
int i= 0;
i <
M.rows();
i++)
1798 L[
i]= (*
N) (
i + 1,
M.columns() + 1);
◆ getMonoms()
extract monomials of F, parts in algebraic variable are considered coefficients
Definition at line 1898 of file cfModGcd.cc.
1914 int numMon=
size (F);
1924 for (
int k= 0;
k < recResult.
size();
k++)
1926 j += recResult.
size();
◆ GFRandomElement()
compute a random element a of GF, s.t. F(a)
, F is a univariate polynomial, returns fail if there are no field elements left which have not been used before
Definition at line 807 of file cfModGcd.cc.
828 while (
find (list, random))
831 if (F (random,
x) == 0)
836 }
while (
find (list, random));
◆ if() [1/2]
◆ if() [2/2]
if |
( |
LCCand * |
absLC(coF) = = abs (LC (F)) | ) |
|
◆ modGCDFp() [1/2]
◆ modGCDFp() [2/2]
Definition at line 1206 of file cfModGcd.cc.
1253 if (best_level == 0)
1280 gcdcAcB=
gcd (cA, cB);
1289 gcdlcAlcB=
gcd (lcA, lcB);
1296 coF=
N (ppA*(cA/gcdcAcB));
1297 coG=
N (ppB*(cB/gcdcAcB));
1308 coF=
N (ppA*(cA/gcdcAcB));
1309 coG=
N (ppB*(cB/gcdcAcB));
1314 CanonicalForm newtonPoly, coF_random_element, coG_random_element,
1315 coF_m, coG_m, ppCoF, ppCoG;
1325 bool inextension=
false;
1328 int bound1=
degree (ppA, 1);
1329 int bound2=
degree (ppB, 1);
1337 if (!fail && !inextension)
1342 modGCDFp (ppA (random_element,
x), ppB (random_element,
x),
1343 coF_random_element, coG_random_element,
topLevel,
1346 "time for recursive call: ");
1347 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
1349 else if (!fail && inextension)
1354 modGCDFq (ppA (random_element,
x), ppB (random_element,
x),
1355 coF_random_element, coG_random_element, V_buf,
1358 "time for recursive call: ");
1359 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
1361 else if (fail && !inextension)
1368 bool initialized=
false;
1387 modGCDFq (ppA (random_element,
x), ppB (random_element,
x),
1388 coF_random_element, coG_random_element,
alpha,
1391 "time for recursive call: ");
1392 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
1394 else if (fail && inextension)
1401 bool prim_fail=
false;
1406 if (V_buf3 !=
alpha)
1409 G_m=
mapDown (G_m, prim_elem, im_prim_elem,
alpha, source, dest);
1410 coF_m=
mapDown (coF_m, prim_elem, im_prim_elem,
alpha, source, dest);
1411 coG_m=
mapDown (coG_m, prim_elem, im_prim_elem,
alpha, source, dest);
1412 newtonPoly=
mapDown (newtonPoly, prim_elem, im_prim_elem,
alpha,
1414 ppA=
mapDown (ppA, prim_elem, im_prim_elem,
alpha, source, dest);
1415 ppB=
mapDown (ppB, prim_elem, im_prim_elem,
alpha, source, dest);
1416 gcdlcAlcB=
mapDown (gcdlcAlcB, prim_elem, im_prim_elem,
alpha, source,
1418 lcA=
mapDown (lcA, prim_elem, im_prim_elem,
alpha, source, dest);
1419 lcB=
mapDown (lcB, prim_elem, im_prim_elem,
alpha, source, dest);
1421 i.getItem()=
mapDown (
i.getItem(), prim_elem, im_prim_elem,
alpha,
1425 ASSERT (!prim_fail,
"failure in integer factorizer");
1435 i.getItem()=
mapUp (
i.getItem(),
alpha, V_buf, prim_elem,
1436 im_prim_elem, source, dest);
1437 m=
mapUp (
m,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1438 G_m=
mapUp (G_m,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1439 coF_m=
mapUp (coF_m,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1440 coG_m=
mapUp (coG_m,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1441 newtonPoly=
mapUp (newtonPoly,
alpha, V_buf, prim_elem, im_prim_elem,
1443 ppA=
mapUp (ppA,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1444 ppB=
mapUp (ppB,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1445 gcdlcAlcB=
mapUp (gcdlcAlcB,
alpha, V_buf, prim_elem, im_prim_elem,
1447 lcA=
mapUp (lcA,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1448 lcB=
mapUp (lcB,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1455 modGCDFq (ppA (random_element,
x), ppB (random_element,
x),
1456 coF_random_element, coG_random_element, V_buf,
1459 "time for recursive call: ");
1460 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
1474 coF=
N (ppA*(cA/gcdcAcB));
1475 coG=
N (ppB*(cB/gcdcAcB));
1481 if (!
find (
l, random_element))
1482 l.append (random_element);
1486 G_random_element= (gcdlcAlcB(random_element,
x)/
uni_lcoeff(G_random_element))
1489 coF_random_element= (lcA(random_element,
x)/
uni_lcoeff(coF_random_element))
1490 *coF_random_element;
1491 coG_random_element= (lcB(random_element,
x)/
uni_lcoeff(coG_random_element))
1492 *coG_random_element;
1511 H=
newtonInterp (random_element, G_random_element, newtonPoly, G_m,
x);
1512 coF=
newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,
x);
1513 coG=
newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,
x);
1515 "time for newton_interpolation: ");
1521 if (gcdlcAlcB.
isOne())
1540 coF=
N ((cA/gcdcAcB)*ppCoF);
1541 coG=
N ((cB/gcdcAcB)*ppCoG);
1543 "time for successful termination Fp: ");
1544 return N(gcdcAcB*ppH);
1547 "time for unsuccessful termination Fp: ");
1553 newtonPoly= newtonPoly*(
x - random_element);
1554 m=
m*(
x - random_element);
1555 if (!
find (
l, random_element))
1556 l.append (random_element);
◆ modGCDFq() [1/2]
GCD of F and G over
, 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 467 of file cfModGcd.cc.
542 gcdcAcB=
gcd (cA, cB);
552 gcdlcAlcB=
gcd (lcA, lcB);
558 coF=
N (ppA*(cA/gcdcAcB));
559 coG=
N (ppB*(cB/gcdcAcB));
568 coF=
N (ppA*(cA/gcdcAcB));
569 coG=
N (ppB*(cB/gcdcAcB));
574 CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
585 bool inextension=
false;
590 int bound1=
degree (ppA, 1);
591 int bound2=
degree (ppB, 1);
602 bool prim_fail=
false;
606 prim_elem_alpha= prim_elem;
608 if (V_buf3 != V_buf4)
610 m=
mapDown (
m, prim_elem, im_prim_elem, V_buf4, source, dest);
611 G_m=
mapDown (G_m, prim_elem, im_prim_elem, V_buf4, source, dest);
612 coF_m=
mapDown (coF_m, prim_elem, im_prim_elem, V_buf4, source, dest);
613 coG_m=
mapDown (coG_m, prim_elem, im_prim_elem, V_buf4, source, dest);
614 newtonPoly=
mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
616 ppA=
mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
617 ppB=
mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
618 gcdlcAlcB=
mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
620 lcA=
mapDown (lcA, prim_elem, im_prim_elem, V_buf4, source, dest);
621 lcB=
mapDown (lcB, prim_elem, im_prim_elem, V_buf4, source, dest);
623 i.getItem()=
mapDown (
i.getItem(), prim_elem, im_prim_elem, V_buf4,
627 ASSERT (!prim_fail,
"failure in integer factorizer");
631 im_prim_elem=
mapPrimElem (prim_elem, V_buf4, V_buf);
634 im_prim_elem_alpha= im_prim_elem;
636 im_prim_elem_alpha=
mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
637 im_prim_elem, source, dest);
642 i.getItem()=
mapUp (
i.getItem(), V_buf4, V_buf, prim_elem,
643 im_prim_elem, source, dest);
644 m=
mapUp (
m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
645 G_m=
mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
646 coF_m=
mapUp (coF_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
647 coG_m=
mapUp (coG_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
648 newtonPoly=
mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
650 ppA=
mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
651 ppB=
mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
652 gcdlcAlcB=
mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
654 lcA=
mapUp (lcA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
655 lcB=
mapUp (lcB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
663 modGCDFq (ppA (random_element,
x), ppB (random_element,
x),
664 coF_random_element, coG_random_element, V_buf,
667 "time for recursive call: ");
668 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
676 modGCDFq (ppA(random_element,
x), ppB(random_element,
x),
677 coF_random_element, coG_random_element, V_buf,
680 "time for recursive call: ");
681 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
695 ppA=
mapDown (ppA, prim_elem_alpha, im_prim_elem_alpha,
alpha, u,
v);
696 ppB=
mapDown (ppB, prim_elem_alpha, im_prim_elem_alpha,
alpha, u,
v);
699 coF=
N (ppA*(cA/gcdcAcB));
700 coG=
N (ppB*(cB/gcdcAcB));
705 if (!
find (
l, random_element))
706 l.append (random_element);
711 (gcdlcAlcB(random_element,
x)/
uni_lcoeff (G_random_element))
713 coF_random_element= (lcA(random_element,
x)/
uni_lcoeff(coF_random_element))
715 coG_random_element= (lcB(random_element,
x)/
uni_lcoeff(coG_random_element))
735 H=
newtonInterp (random_element, G_random_element, newtonPoly, G_m,
x);
736 coF=
newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,
x);
737 coG=
newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,
x);
739 "time for newton interpolation: ");
745 if (gcdlcAlcB.
isOne())
764 DEBOUTLN (cerr,
"ppH before mapDown= " << ppH);
765 ppH=
mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha,
alpha, u,
v);
766 ppCoF=
mapDown (ppCoF, prim_elem_alpha, im_prim_elem_alpha,
alpha, u,
v);
767 ppCoG=
mapDown (ppCoG, prim_elem_alpha, im_prim_elem_alpha,
alpha, u,
v);
768 DEBOUTLN (cerr,
"ppH after mapDown= " << ppH);
769 coF=
N ((cA/gcdcAcB)*ppCoF);
770 coG=
N ((cB/gcdcAcB)*ppCoG);
772 "time for successful termination test Fq: ");
774 return N(gcdcAcB*ppH);
777 else if (((
degree (ppCoF,1)+
degree (ppH,1) == bound1) &&
782 coF=
N ((cA/gcdcAcB)*ppCoF);
783 coG=
N ((cB/gcdcAcB)*ppCoG);
785 "time for successful termination test Fq: ");
786 return N(gcdcAcB*ppH);
789 "time for unsuccessful termination test Fq: ");
795 newtonPoly= newtonPoly*(
x - random_element);
796 m=
m*(
x - random_element);
797 if (!
find (
l, random_element))
798 l.append (random_element);
◆ modGCDFq() [2/2]
◆ modGCDGF() [1/2]
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 860 of file cfModGcd.cc.
935 gcdcAcB=
gcd (cA, cB);
945 gcdlcAlcB=
gcd (lcA, lcB);
950 coF=
N (ppA*(cA/gcdcAcB));
951 coG=
N (ppB*(cB/gcdcAcB));
959 coF=
N (ppA*(cA/gcdcAcB));
960 coG=
N (ppB*(cB/gcdcAcB));
965 CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
976 bool inextension=
false;
982 int bound1=
degree (ppA, 1);
983 int bound2=
degree (ppB, 1);
992 if (
ipower (
p, kk*expon) < (1 << 16))
996 expon= (int)
floor((
log ((
double)((1<<16) - 1)))/(
log((
double)
p)*kk));
997 ASSERT (expon >= 2,
"not enough points in modGCDGF");
1006 newtonPoly=
GFMapUp (newtonPoly, kk);
1011 gcdlcAlcB=
GFMapUp (gcdlcAlcB, kk);
1018 G_random_element=
modGCDGF (ppA(random_element,
x), ppB(random_element,
x),
1019 coF_random_element, coG_random_element,
1022 "time for recursive call: ");
1023 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
1029 G_random_element=
modGCDGF (ppA(random_element,
x), ppB(random_element,
x),
1030 coF_random_element, coG_random_element,
1033 "time for recursive call: ");
1034 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
1051 coF=
N (ppA*(cA/gcdcAcB));
1052 coG=
N (ppB*(cB/gcdcAcB));
1057 if (!
find (
l, random_element))
1058 l.append (random_element);
1063 (gcdlcAlcB(random_element,
x)/
uni_lcoeff(G_random_element)) *
1066 coF_random_element= (lcA(random_element,
x)/
uni_lcoeff(coF_random_element))
1067 *coF_random_element;
1068 coG_random_element= (lcB(random_element,
x)/
uni_lcoeff(coG_random_element))
1069 *coG_random_element;
1088 H=
newtonInterp (random_element, G_random_element, newtonPoly, G_m,
x);
1089 coF=
newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,
x);
1090 coG=
newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,
x);
1092 "time for newton interpolation: ");
1098 if (gcdlcAlcB.
isOne())
1116 DEBOUTLN (cerr,
"ppH before mapDown= " << ppH);
1120 DEBOUTLN (cerr,
"ppH after mapDown= " << ppH);
1121 coF=
N ((cA/gcdcAcB)*ppCoF);
1122 coG=
N ((cB/gcdcAcB)*ppCoG);
1125 "time for successful termination GF: ");
1126 return N(gcdcAcB*ppH);
1136 coF=
N ((cA/gcdcAcB)*ppCoF);
1137 coG=
N ((cB/gcdcAcB)*ppCoG);
1139 "time for successful termination GF: ");
1140 return N(gcdcAcB*ppH);
1144 "time for unsuccessful termination GF: ");
1150 newtonPoly= newtonPoly*(
x - random_element);
1151 m=
m*(
x - random_element);
1152 if (!
find (
l, random_element))
1153 l.append (random_element);
◆ modGCDGF() [2/2]
◆ monicSparseInterpol()
Definition at line 2140 of file cfModGcd.cc.
2153 if (F ==
G)
return F/
Lc(F);
2161 if (best_level == 0)
2179 gcdcAcB=
gcd (cA, cB);
2199 gcdlcAlcB=
gcd (lcA, lcB);
2200 int skelSize=
size (skel, skel.
mvar());
2206 biggestSize=
tmax (biggestSize,
size (
i.coeff()));
2211 bool evalFail=
false;
2222 for (
int i= 0;
i < biggestSize;
i++)
2235 if (V_buf.
level() != 1)
2244 bool prim_fail=
false;
2248 prim_elem_alpha= prim_elem;
2250 ASSERT (!prim_fail,
"failure in integer factorizer");
2254 im_prim_elem=
mapPrimElem (prim_elem, V_buf4, V_buf);
2260 im_prim_elem_alpha= im_prim_elem;
2262 im_prim_elem_alpha=
mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2263 prim_elem, im_prim_elem, source, dest);
2266 j.getItem()=
mapUp (
j.getItem(), V_buf4, V_buf, prim_elem,
2267 im_prim_elem, source, dest);
2268 for (
int k= 0;
k <
i;
k++)
2271 j.getItem()=
mapUp (
j.getItem(), V_buf4, V_buf, prim_elem,
2272 im_prim_elem, source, dest);
2273 gcds[
k]=
mapUp (gcds[
k], V_buf4, V_buf, prim_elem, im_prim_elem,
2279 A=
mapUp (
A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2280 B=
mapUp (
B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2292 bool initialized=
false;
2315 delete[] pEvalPoints;
2324 if (
k.exp() !=
l.exp())
2326 delete[] pEvalPoints;
2343 if (Monoms.
size() == 0)
2346 coeffMonoms=
new CFArray [skelSize];
2353 for (
int i= 0;
i < biggestSize;
i++)
2357 for (
int k= 0;
k < skelSize;
k++,
l++)
2361 pL[
k] [
i]=
l.coeff();
2373 for (
int k= 0;
k < skelSize;
k++)
2375 if (biggestSize != coeffMonoms[
k].
size())
2378 for (
int i= 0;
i < coeffMonoms[
k].
size();
i++)
2379 bufArray [
i]= pL[
k] [
i];
2385 if (solution.
size() == 0)
2387 delete[] pEvalPoints;
2390 delete[] coeffMonoms;
2396 for (
int l= 0;
l < solution.
size();
l++)
2397 result += solution[
l]*Monoms [ind +
l];
2398 ind += solution.
size();
2401 delete[] pEvalPoints;
2404 delete[] coeffMonoms;
◆ mult()
multiply two lists componentwise
Definition at line 2117 of file cfModGcd.cc.
2123 i.getItem() *=
j.getItem();
◆ myCompress()
compressing two polynomials F and G, M is used for compressing, N to reverse the compression
Definition at line 93 of file cfModGcd.cc.
100 for (
int i = n;
i >= 0;
i--)
113 for (
int i= 1;
i <= n;
i++)
142 for (
int i= 1;
i <= n;
i++)
177 while (min_max_deg == 0)
185 for (
int j=
i + 1;
j <=
m;
j++)
233 for (
int i= 1;
i <= n;
i++)
◆ newtonInterp()
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 366 of file cfModGcd.cc.
372 interPoly= oldInterPoly + ((u - oldInterPoly(
alpha,
x))/newtonPoly(
alpha,
x))
◆ nonMonicSparseInterpol()
Definition at line 2424 of file cfModGcd.cc.
2437 if (F ==
G)
return F/
Lc(F);
2445 if (best_level == 0)
2464 gcdcAcB=
gcd (cA, cB);
2484 gcdlcAlcB=
gcd (lcA, lcB);
2485 int skelSize=
size (skel, skel.
mvar());
2493 bufSize=
size (
i.coeff());
2494 biggestSize=
tmax (biggestSize, bufSize);
2495 numberUni += bufSize;
2498 numberUni= (int)
ceil ( (
double) numberUni / (skelSize - 1));
2499 biggestSize=
tmax (biggestSize , numberUni);
2501 numberUni= biggestSize +
size (
LC(skel)) - 1;
2502 int biggestSize2=
tmax (numberUni, biggestSize);
2508 bool evalFail=
false;
2519 for (
int i= 0;
i < biggestSize;
i++)
2531 if (V_buf.
level() != 1)
2540 bool prim_fail=
false;
2544 prim_elem_alpha= prim_elem;
2546 ASSERT (!prim_fail,
"failure in integer factorizer");
2550 im_prim_elem=
mapPrimElem (prim_elem, V_buf4, V_buf);
2556 im_prim_elem_alpha= im_prim_elem;
2558 im_prim_elem_alpha=
mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2559 prim_elem, im_prim_elem, source, dest);
2562 i.getItem()=
mapUp (
i.getItem(), V_buf4, V_buf, prim_elem,
2563 im_prim_elem, source, dest);
2566 A=
mapUp (
A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2567 B=
mapUp (
B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2579 bool initialized=
false;
2608 delete[] pEvalPoints;
2617 if (
k.exp() !=
l.exp())
2619 delete[] pEvalPoints;
2636 if (Monoms.
size() == 0)
2639 coeffMonoms=
new CFArray [skelSize];
2645 int minimalColumnsIndex;
2647 minimalColumnsIndex= 1;
2649 minimalColumnsIndex= 0;
2650 int minimalColumns=-1;
2655 for (
int i= 0;
i < skelSize;
i++)
2659 minimalColumns= coeffMonoms[
i].
size();
2662 if (minimalColumns > coeffMonoms[
i].
size())
2664 minimalColumns= coeffMonoms[
i].
size();
2665 minimalColumnsIndex=
i;
2670 pMat[0]=
CFMatrix (biggestSize, coeffMonoms[0].
size());
2671 pMat[1]=
CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].
size());
2673 for (
int i= 0;
i < biggestSize;
i++)
2677 for (
int k= 0;
k < skelSize;
k++,
l++)
2681 pL[
k] [
i]=
l.coeff();
2683 if (
i == 0 &&
k != 0 &&
k != minimalColumnsIndex)
2685 else if (
i == 0 && (
k == 0 ||
k == minimalColumnsIndex))
2687 if (
i > 0 && (
k == 0 ||
k == minimalColumnsIndex))
2692 if (pMat[
k].rows() >=
i + 1)
2694 for (
int ind= 1; ind < coeffEval.
size() + 1; ind++)
2695 pMat[
k] (
i + 1, ind)= coeffEval[ind - 1];
2698 if (
k == minimalColumnsIndex)
2700 if (pMat[1].rows() >=
i + 1)
2702 for (
int ind= 1; ind < coeffEval.
size() + 1; ind++)
2703 pMat[1] (
i + 1, ind)= coeffEval[ind - 1];
2712 int matRows, matColumns;
2713 matRows= pMat[1].
rows();
2714 matColumns= pMat[0].
columns() - 1;
2715 matColumns += pMat[1].
columns();
2717 Mat=
CFMatrix (matRows, matColumns);
2718 for (
int i= 1;
i <= matRows;
i++)
2720 Mat (
i,
j)= pMat[1] (
i,
j);
2722 for (
int j= pMat[1].columns() + 1;
j <= matColumns;
2725 for (
int i= 1;
i <= matRows;
i++)
2726 Mat (
i,
j)= (-pMat [0] (
i, ind + 1))*pL[minimalColumnsIndex] [
i - 1];
2730 for (
int i= 0;
i < pMat[0].
rows();
i++)
2731 firstColumn [
i]= pMat[0] (
i + 1, 1);
2733 CFArray bufArray= pL[minimalColumnsIndex];
2735 for (
int i= 0;
i < biggestSize;
i++)
2736 pL[minimalColumnsIndex] [
i] *= firstColumn[
i];
2741 if (V_buf.
level() != 1)
2743 pL[minimalColumnsIndex], V_buf);
2746 pL[minimalColumnsIndex]);
2748 if (solution.
size() == 0)
2753 pL[minimalColumnsIndex]= bufArray;
2756 if (biggestSize != biggestSize2)
2758 bufpEvalPoints= pEvalPoints;
2759 pEvalPoints=
new CFList [biggestSize2];
2762 for (
int i= 0;
i < biggestSize;
i++)
2764 pEvalPoints[
i]= bufpEvalPoints [
i];
2765 gcds[
i]= bufGcds[
i];
2767 for (
int i= 0;
i < biggestSize2 - biggestSize;
i++)
2775 delete[] pEvalPoints;
2778 delete[] coeffMonoms;
2780 if (bufpEvalPoints !=
NULL)
2781 delete [] bufpEvalPoints;
2790 if (
k.exp() !=
l.exp())
2792 delete[] pEvalPoints;
2795 delete[] coeffMonoms;
2797 if (bufpEvalPoints !=
NULL)
2798 delete [] bufpEvalPoints;
2806 gcds[
i + biggestSize]=
g;
2809 for (
int i= 0;
i < biggestSize;
i++)
2813 for (
int k= 1;
k < skelSize;
k++,
l++)
2816 pMat[
k]=
CFMatrix (biggestSize2,coeffMonoms[
k].
size()+biggestSize2-1);
2817 if (
k == minimalColumnsIndex)
2820 if (pMat[
k].rows() >=
i + 1)
2822 for (
int ind= 1; ind < coeffEval.
size() + 1; ind++)
2823 pMat[
k] (
i + 1, ind)= coeffEval[ind - 1];
2828 pMat[0]=
CFMatrix (biggestSize2, coeffMonoms[0].
size() + biggestSize2 - 1);
2829 for (
int j= 1;
j <= Mat.
rows();
j++)
2831 pMat [0] (
j,
k)= Mat (
j,
k);
2833 for (
int j= 1;
j <= Mat.
rows();
j++)
2835 pMat [minimalColumnsIndex] (
j,
k)= Mat (
j,
k);
2837 for (
int i= 0;
i < skelSize;
i++)
2841 for (
int j= 0;
j < bufArray.
size();
j++)
2842 pL[
i] [
j]= bufArray [
j];
2845 for (
int i= 0;
i < biggestSize2 - biggestSize;
i++)
2849 for (
int k= 0;
k < skelSize;
k++,
l++)
2851 pL[
k] [
i + biggestSize]=
l.coeff();
2853 if (pMat[
k].rows() >=
i + biggestSize + 1)
2855 for (
int ind= 1; ind < coeffEval.
size() + 1; ind++)
2856 pMat[
k] (
i + biggestSize + 1, ind)= coeffEval[ind - 1];
2861 for (
int i= 0;
i < skelSize;
i++)
2863 if (pL[
i].
size() > 1)
2865 for (
int j= 2;
j <= pMat[
i].
rows();
j++)
2866 pMat[
i] (
j, coeffMonoms[
i].
size() +
j - 1)=
2871 matColumns= biggestSize2 - 1;
2873 for (
int i= 0;
i < skelSize;
i++)
2875 if (V_buf.
level() == 1)
2880 if (pMat[
i] (coeffMonoms[
i].
size(), coeffMonoms[
i].
size()) == 0)
2882 delete[] pEvalPoints;
2885 delete[] coeffMonoms;
2887 if (bufpEvalPoints !=
NULL)
2888 delete [] bufpEvalPoints;
2894 matRows += pMat[
i].
rows() - coeffMonoms[
i].
size();
2898 Mat=
CFMatrix (matRows, matColumns);
2902 for (
int i= 0;
i < skelSize;
i++)
2904 if (coeffMonoms[
i].
size() + 1 >= pMat[
i].rows() || coeffMonoms[
i].
size() + 1 >= pMat[
i].columns())
2906 delete[] pEvalPoints;
2909 delete[] coeffMonoms;
2911 if (bufpEvalPoints !=
NULL)
2912 delete [] bufpEvalPoints;
2918 bufMat= pMat[
i] (coeffMonoms[
i].
size() + 1, pMat[
i].
rows(),
2921 for (
int j= 1;
j <= bufMat.
rows();
j++)
2923 Mat (
j + ind,
k)= bufMat(
j,
k);
2924 bufArray2= coeffMonoms[
i].
size();
2925 for (
int j= 1;
j <= pMat[
i].
rows();
j++)
2927 if (
j > coeffMonoms[
i].
size())
2928 bufArray [
j-coeffMonoms[
i].
size() + ind - 1]= pL[
i] [
j - 1];
2930 bufArray2 [
j - 1]= pL[
i] [
j - 1];
2933 ind += bufMat.
rows();
2937 if (V_buf.
level() != 1)
2942 if (solution.
size() == 0)
2944 delete[] pEvalPoints;
2947 delete[] coeffMonoms;
2949 if (bufpEvalPoints !=
NULL)
2950 delete [] bufpEvalPoints;
2960 for (
int i= 0;
i < skelSize;
i++)
2963 for (
int i= 0;
i < bufSolution.
size();
i++, ind++)
2964 result += Monoms [ind]*bufSolution[
i];
2973 delete[] pEvalPoints;
2976 delete[] coeffMonoms;
2979 if (bufpEvalPoints !=
NULL)
2980 delete [] bufpEvalPoints;
2991 int ind2= 0, ind3= 2;
2993 for (
int l= 0;
l < minimalColumnsIndex;
l++)
2994 ind += coeffMonoms[
l].
size();
2995 for (
int l= solution.
size() - pMat[0].columns() + 1;
l < solution.
size();
2996 l++,
ind2++, ind3++)
2999 for (
int i= 0;
i < pMat[0].
rows();
i++)
3000 firstColumn[
i] += solution[
l]*pMat[0] (
i + 1, ind3);
3002 for (
int l= 0;
l < solution.
size() + 1 - pMat[0].columns();
l++)
3003 result += solution[
l]*Monoms [ind +
l];
3004 ind= coeffMonoms[0].
size();
3005 for (
int k= 1;
k < skelSize;
k++)
3007 if (
k == minimalColumnsIndex)
3009 ind += coeffMonoms[
k].
size();
3012 if (
k != minimalColumnsIndex)
3014 for (
int i= 0;
i < biggestSize;
i++)
3015 pL[
k] [
i] *= firstColumn [
i];
3018 if (biggestSize != coeffMonoms[
k].
size() &&
k != minimalColumnsIndex)
3021 for (
int i= 0;
i < bufArray.
size();
i++)
3022 bufArray [
i]= pL[
k] [
i];
3028 if (solution.
size() == 0)
3030 delete[] pEvalPoints;
3033 delete[] coeffMonoms;
3040 if (
k != minimalColumnsIndex)
3042 for (
int l= 0;
l < solution.
size();
l++)
3043 result += solution[
l]*Monoms [ind +
l];
3044 ind += solution.
size();
3048 delete[] pEvalPoints;
3052 delete[] coeffMonoms;
◆ randomElement()
compute a random element a of
, s.t. F(a)
, F is a univariate polynomial, returns fail if there are no field elements left which have not been used before
Definition at line 381 of file cfModGcd.cc.
392 double bound=
pow ((
double)
p, (
double) d);
403 while (
find (list, random))
408 random= genAlgExt.generate();
409 while (
find (list, random))
410 random= genAlgExt.generate();
412 if (F (random,
x) == 0)
417 }
while (
find (list, random));
◆ readOffSolution() [1/2]
Definition at line 1692 of file cfModGcd.cc.
1697 for (
int i=
M.rows();
i >= 1;
i--)
1702 for (
int j=
M.columns();
j >= 1;
j--,
k++)
1709 if (
k > partialSol.
size() - 1)
1712 tmp3 +=
tmp2*partialSol[partialSol.
size() -
k - 1];
◆ readOffSolution() [2/2]
M in row echolon form, rk rank of M.
Definition at line 1670 of file cfModGcd.cc.
1674 for (
int i= rk;
i >= 1;
i--)
1678 for (
int j=
M.columns() - 1;
j >= 1;
j--)
◆ solveGeneralVandermonde()
Definition at line 1614 of file cfModGcd.cc.
1617 ASSERT (
A.size() == r,
"vector does not have right size");
1625 bool notDistinct=
false;
1626 for (
int i= 0;
i < r - 1;
i++)
1628 for (
int j=
i + 1;
j < r;
j++)
1642 for (
int i= 0;
i < r;
i++)
1643 master *=
x -
M [
i];
1647 for (
int i= 0;
i < r;
i++)
1649 tmp= master/(
x -
M [
i]);
1650 tmp /= tmp (
M [
i], 1);
1657 for (
int i= 1;
i <= r;
i++,
j++)
1661 for (
int l= 1;
l <=
A.size();
l++)
1662 tmp +=
A[
l - 1]*
j.getItem()[
l];
◆ solveSystemFp()
Definition at line 1805 of file cfModGcd.cc.
1807 ASSERT (L.
size() <=
M.rows(),
"dimension exceeded");
1811 for (
int i= 1;
i <=
M.rows();
i++)
1812 for (
int j= 1;
j <=
M.columns();
j++)
1816 for (
int i= 0;
i < L.
size();
i++,
j++)
1817 (*
N) (
j,
M.columns() + 1)= L[
i];
1822 long rk= nmod_mat_rref (FLINTN);
1831 long rk= gauss (*NTLN);
1834 if (rk !=
M.columns())
1837 nmod_mat_clear (FLINTN);
1845 nmod_mat_clear (FLINTN);
◆ solveSystemFq()
Definition at line 1857 of file cfModGcd.cc.
1859 ASSERT (L.
size() <=
M.rows(),
"dimension exceeded");
1863 for (
int i= 1;
i <=
M.rows();
i++)
1864 for (
int j= 1;
j <=
M.columns();
j++)
1867 for (
int i= 0;
i < L.
size();
i++,
j++)
1868 (*
N) (
j,
M.columns() + 1)= L[
i];
1876 zz_pE::init (NTLMipo);
1878 long rk= gauss (*NTLN);
1881 if (rk !=
M.columns())
◆ solveVandermonde()
Definition at line 1561 of file cfModGcd.cc.
1564 ASSERT (
A.size() == r,
"vector does not have right size");
1573 bool notDistinct=
false;
1574 for (
int i= 0;
i < r - 1;
i++)
1576 for (
int j=
i + 1;
j < r;
j++)
1590 for (
int i= 0;
i < r;
i++)
1591 master *=
x -
M [
i];
1594 for (
int i= 0;
i < r;
i++)
1596 tmp= master/(
x -
M [
i]);
1597 tmp /= tmp (
M [
i], 1);
1603 for (
int i= 1;
i <= r;
i++,
j++)
1606 for (
int l= 0;
l <
A.size();
l++)
1607 tmp +=
A[
l]*
j.getItem()[
l];
◆ sparseGCDFp()
Definition at line 3503 of file cfModGcd.cc.
3514 if (F ==
G)
return F/
Lc(F);
3519 if (best_level == 0)
return B.
genOne();
3536 gcdcAcB=
gcd (cA, cB);
3556 gcdlcAlcB=
gcd (lcA, lcB);
3572 CanonicalForm m, random_element, G_m, G_random_element,
H, cH, ppH, skeleton;
3579 bool inextension=
false;
3589 if (random_element == 0 && !fail)
3591 if (!
find (
l, random_element))
3592 l.append (random_element);
3596 if (!fail && !inextension)
3604 "time for recursive call: ");
3605 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3607 else if (!fail && inextension)
3615 "time for recursive call: ");
3616 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3618 else if (fail && !inextension)
3625 bool initialized=
false;
3647 "time for recursive call: ");
3648 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3650 else if (fail && inextension)
3657 bool prim_fail=
false;
3662 if (V_buf3 !=
alpha)
3665 G_m=
mapDown (
m, prim_elem, im_prim_elem,
alpha, source, dest);
3666 newtonPoly=
mapDown (newtonPoly, prim_elem, im_prim_elem,
alpha, source,
3668 ppA=
mapDown (ppA, prim_elem, im_prim_elem,
alpha, source, dest);
3669 ppB=
mapDown (ppB, prim_elem, im_prim_elem,
alpha, source, dest);
3670 gcdlcAlcB=
mapDown (gcdlcAlcB, prim_elem, im_prim_elem,
alpha, source,
3673 i.getItem()=
mapDown (
i.getItem(), prim_elem, im_prim_elem,
alpha,
3677 ASSERT (!prim_fail,
"failure in integer factorizer");
3687 i.getItem()=
mapUp (
i.getItem(),
alpha, V_buf, prim_elem,
3688 im_prim_elem, source, dest);
3689 m=
mapUp (
m,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3690 G_m=
mapUp (G_m,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3691 newtonPoly=
mapUp (newtonPoly,
alpha, V_buf, prim_elem, im_prim_elem,
3693 ppA=
mapUp (ppA,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3694 ppB=
mapUp (ppB,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3695 gcdlcAlcB=
mapUp (gcdlcAlcB,
alpha, V_buf, prim_elem, im_prim_elem,
3703 sparseGCDFq (ppA (random_element,
x), ppB (random_element,
x), V_buf,
3706 "time for recursive call: ");
3707 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3725 if (!
find (
l, random_element))
3726 l.append (random_element);
3731 (gcdlcAlcB(random_element,
x)/
uni_lcoeff (G_random_element))
3734 skeleton= G_random_element;
3750 H=
newtonInterp (random_element, G_random_element, newtonPoly, G_m,
x);
3763 return N(gcdcAcB*ppH);
3767 newtonPoly= newtonPoly*(
x - random_element);
3768 m=
m*(
x - random_element);
3769 if (!
find (
l, random_element))
3770 l.append (random_element);
3783 if (random_element == 0 && !fail)
3785 if (!
find (
l, random_element))
3786 l.append (random_element);
3790 bool sparseFail=
false;
3791 if (!fail && !inextension)
3795 if (
LC (skeleton).inCoeffDomain())
3798 skeleton,
x, sparseFail, coeffMonoms,
3803 skeleton,
x, sparseFail,
3804 coeffMonoms, Monoms);
3806 "time for recursive call: ");
3807 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3809 else if (!fail && inextension)
3813 if (
LC (skeleton).inCoeffDomain())
3816 skeleton,
alpha, sparseFail, coeffMonoms,
3821 skeleton,
alpha, sparseFail, coeffMonoms,
3824 "time for recursive call: ");
3825 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3827 else if (fail && !inextension)
3834 bool initialized=
false;
3852 if (
LC (skeleton).inCoeffDomain())
3855 skeleton,
alpha, sparseFail, coeffMonoms,
3860 skeleton,
alpha, sparseFail, coeffMonoms,
3863 "time for recursive call: ");
3864 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3866 else if (fail && inextension)
3873 bool prim_fail=
false;
3878 if (V_buf3 !=
alpha)
3881 G_m=
mapDown (
m, prim_elem, im_prim_elem,
alpha, source, dest);
3882 newtonPoly=
mapDown (newtonPoly, prim_elem, im_prim_elem,
alpha,
3884 ppA=
mapDown (ppA, prim_elem, im_prim_elem,
alpha, source, dest);
3885 ppB=
mapDown (ppB, prim_elem, im_prim_elem,
alpha, source, dest);
3886 gcdlcAlcB=
mapDown (gcdlcAlcB, prim_elem, im_prim_elem,
alpha,
3889 i.getItem()=
mapDown (
i.getItem(), prim_elem, im_prim_elem,
alpha,
3893 ASSERT (!prim_fail,
"failure in integer factorizer");
3903 i.getItem()=
mapUp (
i.getItem(),
alpha, V_buf, prim_elem,
3904 im_prim_elem, source, dest);
3905 m=
mapUp (
m,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3906 G_m=
mapUp (G_m,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3907 newtonPoly=
mapUp (newtonPoly,
alpha, V_buf, prim_elem, im_prim_elem,
3909 ppA=
mapUp (ppA,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3910 ppB=
mapUp (ppB,
alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3911 gcdlcAlcB=
mapUp (gcdlcAlcB,
alpha, V_buf, prim_elem, im_prim_elem,
3918 if (
LC (skeleton).inCoeffDomain())
3921 skeleton, V_buf, sparseFail, coeffMonoms,
3926 skeleton, V_buf, sparseFail,
3927 coeffMonoms, Monoms);
3929 "time for recursive call: ");
3930 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3951 if (!
find (
l, random_element))
3952 l.append (random_element);
3957 (gcdlcAlcB(random_element,
x)/
uni_lcoeff (G_random_element))
3975 H=
newtonInterp (random_element, G_random_element, newtonPoly, G_m,
x);
3977 "time for newton interpolation: ");
3990 return N(gcdcAcB*ppH);
3995 newtonPoly= newtonPoly*(
x - random_element);
3996 m=
m*(
x - random_element);
3997 if (!
find (
l, random_element))
3998 l.append (random_element);
◆ sparseGCDFq()
Definition at line 3071 of file cfModGcd.cc.
3082 if (F ==
G)
return F/
Lc(F);
3087 if (best_level == 0)
return B.
genOne();
3104 gcdcAcB=
gcd (cA, cB);
3124 gcdlcAlcB=
gcd (lcA, lcB);
3139 CanonicalForm m, random_element, G_m, G_random_element,
H, cH, ppH, skeleton;
3146 bool inextension=
false;
3154 if (random_element == 0 && !fail)
3156 if (!
find (
l, random_element))
3157 l.append (random_element);
3167 bool prim_fail=
false;
3170 if (V_buf4 ==
alpha)
3171 prim_elem_alpha= prim_elem;
3173 if (V_buf3 != V_buf4)
3175 m=
mapDown (
m, prim_elem, im_prim_elem, V_buf4, source, dest);
3176 G_m=
mapDown (
m, prim_elem, im_prim_elem, V_buf4, source, dest);
3177 newtonPoly=
mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3179 ppA=
mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3180 ppB=
mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3181 gcdlcAlcB=
mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4, source,
3184 i.getItem()=
mapDown (
i.getItem(), prim_elem, im_prim_elem, V_buf4,
3188 ASSERT (!prim_fail,
"failure in integer factorizer");
3192 im_prim_elem=
mapPrimElem (prim_elem, V_buf4, V_buf);
3194 if (V_buf4 ==
alpha)
3195 im_prim_elem_alpha= im_prim_elem;
3197 im_prim_elem_alpha=
mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
3198 im_prim_elem, source, dest);
3204 i.getItem()=
mapUp (
i.getItem(), V_buf4, V_buf, prim_elem,
3205 im_prim_elem, source, dest);
3206 m=
mapUp (
m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3207 G_m=
mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3208 newtonPoly=
mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3210 ppA=
mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3211 ppB=
mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3212 gcdlcAlcB=
mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3221 sparseGCDFq (ppA (random_element,
x), ppB (random_element,
x), V_buf,
3224 "time for recursive call: ");
3225 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3233 sparseGCDFq (ppA(random_element,
x), ppB(random_element,
x), V_buf,
3236 "time for recursive call: ");
3237 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3254 if (!
find (
l, random_element))
3255 l.append (random_element);
3260 (gcdlcAlcB(random_element,
x)/
uni_lcoeff (G_random_element))
3263 skeleton= G_random_element;
3278 H=
newtonInterp (random_element, G_random_element, newtonPoly, G_m,
x);
3290 DEBOUTLN (cerr,
"ppH before mapDown= " << ppH);
3291 ppH=
mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha,
alpha, u,
v);
3293 DEBOUTLN (cerr,
"ppH after mapDown= " << ppH);
3295 return N(gcdcAcB*ppH);
3299 return N(gcdcAcB*ppH);
3302 newtonPoly= newtonPoly*(
x - random_element);
3303 m=
m*(
x - random_element);
3304 if (!
find (
l, random_element))
3305 l.append (random_element);
3314 if (random_element == 0 && !fail)
3316 if (!
find (
l, random_element))
3317 l.append (random_element);
3327 bool prim_fail=
false;
3330 if (V_buf4 ==
alpha)
3331 prim_elem_alpha= prim_elem;
3333 if (V_buf3 != V_buf4)
3335 m=
mapDown (
m, prim_elem, im_prim_elem, V_buf4, source, dest);
3336 G_m=
mapDown (
m, prim_elem, im_prim_elem, V_buf4, source, dest);
3337 newtonPoly=
mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3339 ppA=
mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3340 ppB=
mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3341 gcdlcAlcB=
mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
3344 i.getItem()=
mapDown (
i.getItem(), prim_elem, im_prim_elem, V_buf4,
3348 ASSERT (!prim_fail,
"failure in integer factorizer");
3352 im_prim_elem=
mapPrimElem (prim_elem, V_buf4, V_buf);
3354 if (V_buf4 ==
alpha)
3355 im_prim_elem_alpha= im_prim_elem;
3357 im_prim_elem_alpha=
mapUp (im_prim_elem_alpha, V_buf4, V_buf,
3358 prim_elem, im_prim_elem, source, dest);
3364 i.getItem()=
mapUp (
i.getItem(), V_buf4, V_buf, prim_elem,
3365 im_prim_elem, source, dest);
3366 m=
mapUp (
m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3367 G_m=
mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3368 newtonPoly=
mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3370 ppA=
mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3371 ppB=
mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3373 gcdlcAlcB=
mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3385 bool sparseFail=
false;
3386 if (
LC (skeleton).inCoeffDomain())
3389 skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3393 skeleton, V_buf, sparseFail, coeffMonoms,
3399 "time for recursive call: ");
3400 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3406 bool sparseFail=
false;
3407 if (
LC (skeleton).inCoeffDomain())
3410 skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3414 skeleton, V_buf, sparseFail, coeffMonoms,
3420 "time for recursive call: ");
3421 DEBOUTLN (cerr,
"G_random_element= " << G_random_element);
3438 if (!
find (
l, random_element))
3439 l.append (random_element);
3444 (gcdlcAlcB(random_element,
x)/
uni_lcoeff (G_random_element))
3462 H=
newtonInterp (random_element, G_random_element, newtonPoly, G_m,
x);
3464 "time for newton interpolation: ");
3478 DEBOUTLN (cerr,
"ppH before mapDown= " << ppH);
3479 ppH=
mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha,
alpha, u,
v);
3481 DEBOUTLN (cerr,
"ppH after mapDown= " << ppH);
3483 return N(gcdcAcB*ppH);
3488 return N(gcdcAcB*ppH);
3493 newtonPoly= newtonPoly*(
x - random_element);
3494 m=
m*(
x - random_element);
3495 if (!
find (
l, random_element))
3496 l.append (random_element);
◆ 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]
compute the content of F, where F is considered as an element of
Definition at line 283 of file cfModGcd.cc.
300 for (;
i.hasTerms();
i++)
◆ uni_content() [2/2]
◆ uni_lcoeff()
compute the leading coefficient of F, where F is considered as an element of
, order on
is dp.
Definition at line 341 of file cfModGcd.cc.
◆ while()
Definition at line 4073 of file cfModGcd.cc.
4102 "time for gcd mod p in modular gcd: ");
4181 "time for successful termination in modular gcd: ");
4186 "time for unsuccessful termination in modular gcd: ");
◆ cand
◆ cd
◆ cDn
◆ cf
◆ cg
◆ cl
◆ coF
◆ cof
◆ cofn
◆ cofp
◆ coG
◆ cog
◆ cogn
◆ cogp
◆ d_deg
◆ Dn
◆ dp_deg
◆ equal
◆ false
◆ fp
◆ gcdcfcg
◆ GG
◆ gp
◆ lcf
◆ lcg
◆ log2exp
const double log2exp = 1.442695041 |
|
static |
◆ maxNumVars
◆ minCommonDeg
◆ newCof
◆ newCog
◆ test
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...
CFArray evaluate(const CFArray &A, const CFList &evalPoints)
generate random elements in F_p
CFMatrix * convertNTLmat_zz_pE2FacCFMatrix(const mat_zz_pE &m, const Variable &alpha)
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
bool terminationTest(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &coF, const CanonicalForm &coG, const CanonicalForm &cand)
mat_zz_pE * convertFacCFMatrix2NTLmat_zz_pE(const CFMatrix &m)
class to iterate through CanonicalForm's
#define DEBOUTLN(stream, objects)
static CanonicalForm uni_content(const CanonicalForm &F)
compute the content of F, where F is considered as an element of
void setMipo(const Variable &alpha, const CanonicalForm &mipo)
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...
static CanonicalForm FpRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
CFArray solveSystemFp(const CFMatrix &M, const CFArray &L)
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...
generate random elements in F_p(alpha)
CanonicalForm generate() const
CFMatrix * convertNTLmat_zz_p2FacCFMatrix(const mat_zz_p &m)
Variable rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables
CFArray readOffSolution(const CFMatrix &M, const long rk)
M in row echolon form, rk rank of M.
const signed long floor(const ampf< Precision > &x)
const CanonicalForm CFMap CFMap & N
CanonicalForm randomIrredpoly(int i, const Variable &x)
computes a random monic irreducible univariate polynomial in x over Fp of degree i via NTL
int cf_getBigPrime(int i)
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
CFArray solveGeneralVandermonde(const CFArray &M, const CFArray &A)
const signed long ceil(const ampf< Precision > &x)
CanonicalForm GFMapDown(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
CFArray getMonoms(const CanonicalForm &F)
extract monomials of F, parts in algebraic variable are considered coefficients
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
#define ASSERT(expression, message)
Rational abs(const Rational &a)
int status int void * buf
TIMING_START(fac_alg_resultant)
CFList *& Aeval
<[in] poly
static void evalPoint(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &FEval, CanonicalForm &GEval, CFGenerator &evalPoint)
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...
void mult(CFList &L1, const CFList &L2)
multiply two lists componentwise
const CanonicalForm const CanonicalForm & coF
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....
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.
void prune(Variable &alpha)
const CanonicalForm const CanonicalForm const CanonicalForm & coG
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
void 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,...
static const int SW_USE_CHINREM_GCD
set to 1 to use modular gcd over Z
mat_zz_p * convertFacCFMatrix2NTLmat_zz_p(const CFMatrix &m)
generate random elements in GF
int ipower(int b, int m)
int ipower ( int b, int m )
CanonicalForm monicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
const CanonicalForm CFMap CFMap int & both_non_zero
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
CanonicalForm modGCDFp(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l)
TIMING_END_AND_PRINT(fac_alg_resultant, "time to compute resultant0: ")
factory's class for variables
gmp_float log(const gmp_float &a)
static CanonicalForm bound(const CFMatrix &M)
CanonicalForm nonMonicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
const CanonicalForm CFMap CFMap bool topLevel
long gaussianElimFp(CFMatrix &M, CFArray &L)
CanonicalForm mapPrimElem(const CanonicalForm &primElem, const Variable &alpha, const Variable &beta)
compute the image of a primitive element of in . We assume .
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 eval(const CanonicalForm &A, const CanonicalForm &B, CanonicalForm &Aeval, CanonicalForm &Beval, const CFList &L)
Rational pow(const Rational &a, int e)
CanonicalForm generate() const
CanonicalForm balance_p(const CanonicalForm &f, const CanonicalForm &q, const CanonicalForm &qh)
same as balance_p ( const CanonicalForm & f, const CanonicalForm & q ) but qh= q/2 is provided,...
template bool find(const List< CanonicalForm > &, const CanonicalForm &)
const Variable & v
< [in] a sqrfree bivariate poly
CFMatrix * convertNmod_mat_t2FacCFMatrix(const nmod_mat_t m)
conversion of a FLINT matrix over Z/p to a factory matrix
CFArray evaluateMonom(const CanonicalForm &F, const CFList &evalPoints)
long gaussianElimFq(CFMatrix &M, CFArray &L, const Variable &alpha)
CanonicalForm GFMapUp(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
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
void convertFacCFMatrix2nmod_mat_t(nmod_mat_t M, const CFMatrix &m)
conversion of a factory matrix over Z/p to a nmod_mat_t
void prune1(const Variable &alpha)
CanonicalForm sparseGCDFq(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
CFArray solveSystemFq(const CFMatrix &M, const CFArray &L, const Variable &alpha)
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...
static Variable chooseExtension(const Variable &alpha)
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 ,...
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)