My Project
Functions | Variables
cfGcdAlgExt.cc File Reference
#include "config.h"
#include <stdio.h>
#include <iostream>
#include "cf_assert.h"
#include "timing.h"
#include "templates/ftmpl_functions.h"
#include "cf_defs.h"
#include "canonicalform.h"
#include "cf_iter.h"
#include "cf_primes.h"
#include "cf_algorithm.h"
#include "cfGcdAlgExt.h"
#include "cfUnivarGcd.h"
#include "cf_map.h"
#include "cf_generator.h"
#include "facMul.h"
#include "cfNTLzzpEXGCD.h"
#include "NTLconvert.h"
#include "FLINTconvert.h"

Go to the source code of this file.

Functions

 TIMING_DEFINE_PRINT (alg_content_p) TIMING_DEFINE_PRINT(alg_content) TIMING_DEFINE_PRINT(alg_compress) TIMING_DEFINE_PRINT(alg_termination) TIMING_DEFINE_PRINT(alg_termination_p) TIMING_DEFINE_PRINT(alg_reconstruction) TIMING_DEFINE_PRINT(alg_newton_p) TIMING_DEFINE_PRINT(alg_recursion_p) TIMING_DEFINE_PRINT(alg_gcd_p) TIMING_DEFINE_PRINT(alg_euclid_p) static int myCompress(const CanonicalForm &F
 compressing two polynomials F and G, M is used for compressing, N to reverse the compression More...
 
 for (int i=0;i<=n;i++) degsf[i]
 
 if (topLevel)
 
 DELETE_ARRAY (degsg)
 
void tryInvert (const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
 
static CanonicalForm trycontent (const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
 
static CanonicalForm tryvcontent (const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
 
static CanonicalForm trycf_content (const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
 
static CanonicalForm tryNewtonInterp (const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x, const CanonicalForm &M, bool &fail)
 
void tryBrownGCD (const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M, CanonicalForm &result, bool &fail, bool topLevel)
 modular gcd over F_p[x]/(M) for not necessarily irreducible M. If a zero divisor is encountered fail is set to true. More...
 
static CanonicalForm myicontent (const CanonicalForm &f, const CanonicalForm &c)
 
static CanonicalForm myicontent (const CanonicalForm &f)
 
CanonicalForm QGCD (const CanonicalForm &F, const CanonicalForm &G)
 gcd over Q(a) More...
 
int * leadDeg (const CanonicalForm &f, int *degs)
 
bool isLess (int *a, int *b, int lower, int upper)
 
bool isEqual (int *a, int *b, int lower, int upper)
 
CanonicalForm firstLC (const CanonicalForm &f)
 

Variables

const CanonicalFormG
 
const CanonicalForm CFMapM
 
const CanonicalForm CFMap CFMapN
 
const CanonicalForm CFMap CFMap bool topLevel
 
int * degsf = NEW_ARRAY(int,n + 1)
 
int * degsg = NEW_ARRAY(int,n + 1)
 
int both_non_zero = 0
 
int f_zero = 0
 
int g_zero = 0
 
int both_zero = 0
 
int Flevel =F.level()
 
int Glevel =G.level()
 
 else
 
 return
 

Function Documentation

◆ DELETE_ARRAY()

DELETE_ARRAY ( degsg  )

◆ firstLC()

CanonicalForm firstLC ( const CanonicalForm f)

Definition at line 957 of file cfGcdAlgExt.cc.

958{ // returns the leading coefficient (LC) of level <= 1
959 CanonicalForm ret = f;
960 while(ret.level() > 1)
961 ret = LC(ret);
962 return ret;
963}
CanonicalForm LC(const CanonicalForm &f)
f
Definition: cfModGcd.cc:4083
factory's main class
Definition: canonicalform.h:86
int level() const
level() returns the level of CO.

◆ for()

for ( int  i = 0;i<=n;i++)

Definition at line 72 of file cfEzgcd.cc.

73 {
74 if (degsf[i] != 0 && degsg[i] != 0)
75 {
77 continue;
78 }
79 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
80 {
81 f_zero++;
82 continue;
83 }
84 if (degsg[i] == 0 && degsf[i] && i <= F.level())
85 {
86 g_zero++;
87 continue;
88 }
89 }
int * degsf
Definition: cfEzgcd.cc:59
int f_zero
Definition: cfEzgcd.cc:69
const CanonicalForm & G
Definition: cfEzgcd.cc:55
const CanonicalForm CFMap CFMap int & both_non_zero
Definition: cfEzgcd.cc:57
int i
Definition: cfEzgcd.cc:132
int g_zero
Definition: cfEzgcd.cc:70
int * degsg
Definition: cfEzgcd.cc:60

◆ if()

if ( topLevel  )

Definition at line 75 of file cfGcdAlgExt.cc.

76 {
77 for (int i= 1; i <= n; i++)
78 {
79 if (degsf[i] != 0 && degsg[i] != 0)
80 {
82 continue;
83 }
84 if (degsf[i] == 0 && degsg[i] != 0 && i <= Glevel)
85 {
86 f_zero++;
87 continue;
88 }
89 if (degsg[i] == 0 && degsf[i] && i <= Flevel)
90 {
91 g_zero++;
92 continue;
93 }
94 }
95
96 if (both_non_zero == 0)
97 {
100 return 0;
101 }
102
103 // map Variables which do not occur in both polynomials to higher levels
104 int k= 1;
105 int l= 1;
106 for (int i= 1; i <= n; i++)
107 {
108 if (degsf[i] != 0 && degsg[i] == 0 && i <= Flevel)
109 {
110 if (k + both_non_zero != i)
111 {
112 M.newpair (Variable (i), Variable (k + both_non_zero));
113 N.newpair (Variable (k + both_non_zero), Variable (i));
114 }
115 k++;
116 }
117 if (degsf[i] == 0 && degsg[i] != 0 && i <= Glevel)
118 {
119 if (l + g_zero + both_non_zero != i)
120 {
121 M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
122 N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
123 }
124 l++;
125 }
126 }
127
128 // sort Variables x_{i} in increasing order of
129 // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
130 int m= tmax (Flevel, Glevel);
131 int min_max_deg;
133 l= 0;
134 int i= 1;
135 while (k > 0)
136 {
137 if (degsf [i] != 0 && degsg [i] != 0)
138 min_max_deg= tmax (degsf[i], degsg[i]);
139 else
140 min_max_deg= 0;
141 while (min_max_deg == 0)
142 {
143 i++;
144 min_max_deg= tmax (degsf[i], degsg[i]);
145 if (degsf [i] != 0 && degsg [i] != 0)
146 min_max_deg= tmax (degsf[i], degsg[i]);
147 else
148 min_max_deg= 0;
149 }
150 for (int j= i + 1; j <= m; j++)
151 {
152 if (tmax (degsf[j],degsg[j]) <= min_max_deg && degsf[j] != 0 && degsg [j] != 0)
153 {
154 min_max_deg= tmax (degsf[j], degsg[j]);
155 l= j;
156 }
157 }
158 if (l != 0)
159 {
160 if (l != k)
161 {
162 M.newpair (Variable (l), Variable(k));
163 N.newpair (Variable (k), Variable(l));
164 degsf[l]= 0;
165 degsg[l]= 0;
166 l= 0;
167 }
168 else
169 {
170 degsf[l]= 0;
171 degsg[l]= 0;
172 l= 0;
173 }
174 }
175 else if (l == 0)
176 {
177 if (i != k)
178 {
179 M.newpair (Variable (i), Variable (k));
180 N.newpair (Variable (k), Variable (i));
181 degsf[i]= 0;
182 degsg[i]= 0;
183 }
184 else
185 {
186 degsf[i]= 0;
187 degsg[i]= 0;
188 }
189 i++;
190 }
191 k--;
192 }
193 }
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int k
Definition: cfEzgcd.cc:99
const CanonicalForm CFMap CFMap & N
Definition: cfGcdAlgExt.cc:56
int both_non_zero
Definition: cfGcdAlgExt.cc:68
int * degsf
Definition: cfGcdAlgExt.cc:59
int f_zero
Definition: cfGcdAlgExt.cc:69
int Flevel
Definition: cfGcdAlgExt.cc:72
int Glevel
Definition: cfGcdAlgExt.cc:73
const CanonicalForm CFMap & M
Definition: cfGcdAlgExt.cc:55
int g_zero
Definition: cfGcdAlgExt.cc:70
int * degsg
Definition: cfGcdAlgExt.cc:60
DELETE_ARRAY(degsg)
factory's class for variables
Definition: factory.h:134
int j
Definition: facHensel.cc:110
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)

◆ isEqual()

bool isEqual ( int *  a,
int *  b,
int  lower,
int  upper 
)

Definition at line 948 of file cfGcdAlgExt.cc.

949{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
950 for(int i=lower; i<=upper; i++)
951 if(a[i] != b[i])
952 return false;
953 return true;
954}
CanonicalForm b
Definition: cfModGcd.cc:4105

◆ isLess()

bool isLess ( int *  a,
int *  b,
int  lower,
int  upper 
)

Definition at line 937 of file cfGcdAlgExt.cc.

938{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
939 for(int i=upper; i>=lower; i--)
940 if(a[i] == b[i])
941 continue;
942 else
943 return a[i] < b[i];
944 return true;
945}

◆ leadDeg()

int * leadDeg ( const CanonicalForm f,
int *  degs 
)

Definition at line 920 of file cfGcdAlgExt.cc.

921{ // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
922 // if f is in a coeff domain, the zero pointer is returned
923 // 'a' should point to an array of sufficient size level(f)+1
924 if(f.inCoeffDomain())
925 return 0;
926 CanonicalForm tmp = f;
927 do
928 {
929 degs[tmp.level()] = tmp.degree();
930 tmp = LC(tmp);
931 }
932 while(!tmp.inCoeffDomain());
933 return degs;
934}
int degree() const
Returns -1 for the zero polynomial and 0 if CO is in a base domain.
bool inCoeffDomain() const

◆ myicontent() [1/2]

static CanonicalForm myicontent ( const CanonicalForm f)
static

Definition at line 707 of file cfGcdAlgExt.cc.

708{
709#if defined(HAVE_NTL) || defined(HAVE_FLINT)
710 return myicontent( f, 0 );
711#else
712 return 1;
713#endif
714}
static CanonicalForm myicontent(const CanonicalForm &f, const CanonicalForm &c)
Definition: cfGcdAlgExt.cc:657

◆ myicontent() [2/2]

static CanonicalForm myicontent ( const CanonicalForm f,
const CanonicalForm c 
)
static

Definition at line 657 of file cfGcdAlgExt.cc.

658{
659#if defined(HAVE_NTL) || defined(HAVE_FLINT)
660 if (f.isOne() || c.isOne())
661 return 1;
662 if ( f.inBaseDomain() && c.inBaseDomain())
663 {
664 if (c.isZero()) return abs(f);
665 return bgcd( f, c );
666 }
667 else if ( (f.inCoeffDomain() && c.inCoeffDomain()) ||
668 (f.inCoeffDomain() && c.inBaseDomain()) ||
669 (f.inBaseDomain() && c.inCoeffDomain()))
670 {
671 if (c.isZero()) return abs (f);
672#ifdef HAVE_FLINT
673 fmpz_poly_t FLINTf, FLINTc;
674 convertFacCF2Fmpz_poly_t (FLINTf, f);
675 convertFacCF2Fmpz_poly_t (FLINTc, c);
676 fmpz_poly_gcd (FLINTc, FLINTc, FLINTf);
678 if (f.inCoeffDomain())
679 result= convertFmpz_poly_t2FacCF (FLINTc, f.mvar());
680 else
681 result= convertFmpz_poly_t2FacCF (FLINTc, c.mvar());
682 fmpz_poly_clear (FLINTc);
683 fmpz_poly_clear (FLINTf);
684 return result;
685#else
686 ZZX NTLf= convertFacCF2NTLZZX (f);
687 ZZX NTLc= convertFacCF2NTLZZX (c);
688 NTLc= GCD (NTLc, NTLf);
689 if (f.inCoeffDomain())
690 return convertNTLZZX2CF(NTLc,f.mvar());
691 else
692 return convertNTLZZX2CF(NTLc,c.mvar());
693#endif
694 }
695 else
696 {
697 CanonicalForm g = c;
698 for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
699 g = myicontent( i.coeff(), g );
700 return g;
701 }
702#else
703 return 1;
704#endif
705}
CanonicalForm convertFmpz_poly_t2FacCF(const fmpz_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z to CanonicalForm
void convertFacCF2Fmpz_poly_t(fmpz_poly_t result, const CanonicalForm &f)
conversion of a factory univariate polynomial over Z to a fmpz_poly_t
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
Definition: NTLconvert.cc:691
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:285
CanonicalForm bgcd(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm bgcd ( const CanonicalForm & f, const CanonicalForm & g )
g
Definition: cfModGcd.cc:4092
CanonicalForm abs(const CanonicalForm &f)
inline CanonicalForm abs ( const CanonicalForm & f )
Definition: cf_algorithm.h:117
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
CF_NO_INLINE bool isZero() const
Definition: cf_inline.cc:379
bool inBaseDomain() const
CF_NO_INLINE bool isOne() const
CF_INLINE bool CanonicalForm::isOne, isZero () const.
Definition: cf_inline.cc:361
return result
Definition: facAbsBiFact.cc:75

◆ QGCD()

gcd over Q(a)

Definition at line 716 of file cfGcdAlgExt.cc.

717{ // f,g in Q(a)[x1,...,xn]
718 if(F.isZero())
719 {
720 if(G.isZero())
721 return G; // G is zero
722 if(G.inCoeffDomain())
723 return CanonicalForm(1);
724 CanonicalForm lcinv= 1/Lc (G);
725 return G*lcinv; // return monic G
726 }
727 if(G.isZero()) // F is non-zero
728 {
729 if(F.inCoeffDomain())
730 return CanonicalForm(1);
731 CanonicalForm lcinv= 1/Lc (F);
732 return F*lcinv; // return monic F
733 }
734 if(F.inCoeffDomain() || G.inCoeffDomain())
735 return CanonicalForm(1);
736 // here: both NOT inCoeffDomain
737 CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
738 int p, i;
739 int *bound, *other; // degree vectors
740 bool fail;
741 bool off_rational=!isOn(SW_RATIONAL);
742 On( SW_RATIONAL ); // needed by bCommonDen
743 f = F * bCommonDen(F);
744 g = G * bCommonDen(G);
746 CanonicalForm contf= myicontent (f);
747 CanonicalForm contg= myicontent (g);
748 f /= contf;
749 g /= contg;
750 CanonicalForm gcdcfcg= myicontent (contf, contg);
751 TIMING_END_AND_PRINT (alg_content, "time for content in alg gcd: ")
752 Variable a, b;
754 {
755 if(hasFirstAlgVar(g,b))
756 {
757 if(b.level() > a.level())
758 a = b;
759 }
760 }
761 else
762 {
763 if(!hasFirstAlgVar(g,a))// both not in extension
764 {
765 Off( SW_RATIONAL );
766 Off( SW_USE_QGCD );
767 tmp = gcdcfcg*gcd( f, g );
768 On( SW_USE_QGCD );
769 if (off_rational) Off(SW_RATIONAL);
770 return tmp;
771 }
772 }
773 // here: a is the biggest alg. var in f and g AND some of f,g is in extension
774 setReduce(a,false); // do not reduce expressions modulo mipo
775 tmp = getMipo(a);
776 M = tmp * bCommonDen(tmp);
777 // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
778 Off( SW_RATIONAL ); // needed by mod
779 // calculate upper bound for degree vector of gcd
780 int mv = f.level(); i = g.level();
781 if(i > mv)
782 mv = i;
783 // here: mv is level of the largest variable in f, g
784 bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
785 other = new int[mv+1];
786 for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
787 bound[i] = other[i] = 0;
788 bound = leadDeg(f,bound); // 'bound' is set the leading degree vector of f
789 other = leadDeg(g,other);
790 for(int i=1; i<=mv; i++) // entry at i=0 not visited
791 if(other[i] < bound[i])
792 bound[i] = other[i];
793 // now 'bound' is the smaller vector
794 cl = lc(M) * lc(f) * lc(g);
795 q = 1;
796 D = 0;
798 bool equal= false;
799 for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
800 {
801 p = cf_getBigPrime(i);
802 if( mod( cl, p ).isZero() ) // skip lc-bad primes
803 continue;
804 fail = false;
806 mipo = mapinto(M);
807 mipo /= mipo.lc();
808 // here: mipo is monic
809 TIMING_START (alg_gcd_p)
810 tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
811 TIMING_END_AND_PRINT (alg_gcd_p, "time for alg gcd mod p: ")
812 if( fail ) // mipo splits in char p
813 {
815 continue;
816 }
817 if( Dp.inCoeffDomain() ) // early termination
818 {
819 tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
821 if(fail)
822 continue;
823 setReduce(a,true);
824 if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
825 delete[] other;
826 delete[] bound;
827 return gcdcfcg;
828 }
830 // here: Dp NOT inCoeffDomain
831 for(int i=1; i<=mv; i++)
832 other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
833 other = leadDeg(Dp,other);
834
835 if(isEqual(bound, other, 1, mv)) // equal
836 {
837 chineseRemainder( D, q, mapinto(Dp), p, tmp, newq );
838 // tmp = Dp mod p
839 // tmp = D mod q
840 // newq = p*q
841 q = newq;
842 if( D != tmp )
843 D = tmp;
844 On( SW_RATIONAL );
845 TIMING_START (alg_reconstruction)
846 tmp = Farey( D, q ); // Farey
847 tmp *= bCommonDen (tmp);
848 TIMING_END_AND_PRINT (alg_reconstruction,
849 "time for rational reconstruction in alg gcd: ")
850 setReduce(a,true); // reduce expressions modulo mipo
851 On( SW_RATIONAL ); // needed by fdivides
852 if (test != tmp)
853 test= tmp;
854 else
855 equal= true; // modular image did not add any new information
856 TIMING_START (alg_termination)
857#ifdef HAVE_NTL
858#ifdef HAVE_FLINT
859 if (equal && tmp.isUnivariate() && f.isUnivariate() && g.isUnivariate()
860 && f.level() == tmp.level() && tmp.level() == g.level())
861 {
863 newtonDivrem (f, tmp, Q, R);
864 if (R.isZero())
865 {
866 newtonDivrem (g, tmp, Q, R);
867 if (R.isZero())
868 {
870 setReduce (a,true);
871 if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
872 TIMING_END_AND_PRINT (alg_termination,
873 "time for successful termination test in alg gcd: ")
874 delete[] other;
875 delete[] bound;
876 return tmp*gcdcfcg;
877 }
878 }
879 }
880 else
881#endif
882#endif
883 if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
884 {
885 Off( SW_RATIONAL );
886 setReduce(a,true);
887 if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
888 TIMING_END_AND_PRINT (alg_termination,
889 "time for successful termination test in alg gcd: ")
890 delete[] other;
891 delete[] bound;
892 return tmp*gcdcfcg;
893 }
894 TIMING_END_AND_PRINT (alg_termination,
895 "time for unsuccessful termination test in alg gcd: ")
896 Off( SW_RATIONAL );
897 setReduce(a,false); // do not reduce expressions modulo mipo
898 continue;
899 }
900 if( isLess(bound, other, 1, mv) ) // current prime unlucky
901 continue;
902 // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
903 q = p;
904 D = mapinto(Dp); // shortcut CRA // shortcut CRA
905 for(int i=1; i<=mv; i++) // tighten bound
906 bound[i] = other[i];
907 }
908 // hopefully, we never reach this point
909 setReduce(a,true);
910 delete[] other;
911 delete[] bound;
912 Off( SW_USE_QGCD );
913 D = gcdcfcg*gcd( f, g );
914 On( SW_USE_QGCD );
915 if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
916 return D;
917}
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm mapinto(const CanonicalForm &f)
CanonicalForm lc(const CanonicalForm &f)
CanonicalForm FACTORY_PUBLIC gcd(const CanonicalForm &, const CanonicalForm &)
Definition: cf_gcd.cc:685
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:571
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
bool hasFirstAlgVar(const CanonicalForm &f, Variable &a)
check if poly f contains an algebraic variable a
Definition: cf_ops.cc:679
CanonicalForm Lc(const CanonicalForm &f)
else
Definition: cfGcdAlgExt.cc:195
bool isLess(int *a, int *b, int lower, int upper)
Definition: cfGcdAlgExt.cc:937
bool isEqual(int *a, int *b, int lower, int upper)
Definition: cfGcdAlgExt.cc:948
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:72
if(topLevel)
Definition: cfGcdAlgExt.cc:75
return
Definition: cfGcdAlgExt.cc:218
const CanonicalForm & G
Definition: cfGcdAlgExt.cc:55
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
Definition: cfGcdAlgExt.cc:221
int * leadDeg(const CanonicalForm &f, int *degs)
Definition: cfGcdAlgExt.cc:920
int p
Definition: cfModGcd.cc:4080
return false
Definition: cfModGcd.cc:84
cl
Definition: cfModGcd.cc:4102
CanonicalForm gcdcfcg
Definition: cfModGcd.cc:4103
bool equal
Definition: cfModGcd.cc:4128
CanonicalForm test
Definition: cfModGcd.cc:4098
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
Definition: cf_chinese.cc:202
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_QGCD
set to 1 to use Encarnacion GCD over Q(a)
Definition: cf_defs.h:42
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:30
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
int cf_getNumBigPrimes()
Definition: cf_primes.cc:45
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
int level() const
Definition: factory.h:150
CanonicalForm mipo
Definition: facAlgExt.cc:57
TIMING_END_AND_PRINT(fac_alg_resultant, "time to compute resultant0: ")
TIMING_START(fac_alg_resultant)
CanonicalForm alg_content(const CanonicalForm &f, const CFList &as)
Definition: facAlgFunc.cc:42
void newtonDivrem(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Q, CanonicalForm &R)
division with remainder of univariate polynomials over Q and Q(a) using Newton inversion,...
Definition: facMul.cc:346
bool isZero(const CFArray &A)
checks if entries of A are zero
void setReduce(const Variable &alpha, bool reduce)
Definition: variable.cc:238
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
#define D(A)
Definition: gentable.cc:131
STATIC_VAR jList * Q
Definition: janet.cc:30
#define R
Definition: sirandom.c:27

◆ TIMING_DEFINE_PRINT()

TIMING_DEFINE_PRINT ( alg_content_p  ) const &

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

◆ tryBrownGCD()

void tryBrownGCD ( const CanonicalForm F,
const CanonicalForm G,
const CanonicalForm M,
CanonicalForm result,
bool &  fail,
bool  topLevel 
)

modular gcd over F_p[x]/(M) for not necessarily irreducible M. If a zero divisor is encountered fail is set to true.

Definition at line 372 of file cfGcdAlgExt.cc.

373{ // assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
374 // M is assumed to be monic
375 if(F.isZero())
376 {
377 if(G.isZero())
378 {
379 result = G; // G is zero
380 return;
381 }
382 if(G.inCoeffDomain())
383 {
384 tryInvert(G,M,result,fail);
385 if(fail)
386 return;
387 result = 1;
388 return;
389 }
390 // try to make G monic modulo M
391 CanonicalForm inv;
392 tryInvert(Lc(G),M,inv,fail);
393 if(fail)
394 return;
395 result = inv*G;
396 result= reduce (result, M);
397 return;
398 }
399 if(G.isZero()) // F is non-zero
400 {
401 if(F.inCoeffDomain())
402 {
403 tryInvert(F,M,result,fail);
404 if(fail)
405 return;
406 result = 1;
407 return;
408 }
409 // try to make F monic modulo M
410 CanonicalForm inv;
411 tryInvert(Lc(F),M,inv,fail);
412 if(fail)
413 return;
414 result = inv*F;
415 result= reduce (result, M);
416 return;
417 }
418 // here: F,G both nonzero
419 if(F.inCoeffDomain())
420 {
421 tryInvert(F,M,result,fail);
422 if(fail)
423 return;
424 result = 1;
425 return;
426 }
427 if(G.inCoeffDomain())
428 {
429 tryInvert(G,M,result,fail);
430 if(fail)
431 return;
432 result = 1;
433 return;
434 }
435 TIMING_START (alg_compress)
436 CFMap MM,NN;
437 int lev= myCompress (F, G, MM, NN, topLevel);
438 if (lev == 0)
439 {
440 result= 1;
441 return;
442 }
443 CanonicalForm f=MM(F);
444 CanonicalForm g=MM(G);
445 TIMING_END_AND_PRINT (alg_compress, "time to compress in alg gcd: ")
446 // here: f,g are compressed
447 // compute largest variable in f or g (least one is Variable(1))
448 int mv = f.level();
449 if(g.level() > mv)
450 mv = g.level();
451 // here: mv is level of the largest variable in f, g
452 Variable v1= Variable (1);
453#ifdef HAVE_NTL
454 Variable v= M.mvar();
455 int ch=getCharacteristic();
456 if (fac_NTL_char != ch)
457 {
458 fac_NTL_char= ch;
459 zz_p::init (ch);
460 }
461 zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
462 zz_pE::init (NTLMipo);
463 zz_pEX NTLResult;
464 zz_pEX NTLF;
465 zz_pEX NTLG;
466#endif
467 if(mv == 1) // f,g univariate
468 {
469 TIMING_START (alg_euclid_p)
470#ifdef HAVE_NTL
471 NTLF= convertFacCF2NTLzz_pEX (f, NTLMipo);
472 NTLG= convertFacCF2NTLzz_pEX (g, NTLMipo);
473 tryNTLGCD (NTLResult, NTLF, NTLG, fail);
474 if (fail)
475 return;
476 result= convertNTLzz_pEX2CF (NTLResult, f.mvar(), v);
477#else
478 tryEuclid(f,g,M,result,fail);
479 if(fail)
480 return;
481#endif
482 TIMING_END_AND_PRINT (alg_euclid_p, "time for euclidean alg mod p: ")
483 result= NN (reduce (result, M)); // do not forget to map back
484 return;
485 }
486 TIMING_START (alg_content_p)
487 // here: mv > 1
488 CanonicalForm cf = tryvcontent(f, Variable(2), M, fail); // cf is univariate poly in var(1)
489 if(fail)
490 return;
491 CanonicalForm cg = tryvcontent(g, Variable(2), M, fail);
492 if(fail)
493 return;
495#ifdef HAVE_NTL
496 NTLF= convertFacCF2NTLzz_pEX (cf, NTLMipo);
497 NTLG= convertFacCF2NTLzz_pEX (cg, NTLMipo);
498 tryNTLGCD (NTLResult, NTLF, NTLG, fail);
499 if (fail)
500 return;
501 c= convertNTLzz_pEX2CF (NTLResult, v1, v);
502#else
503 tryEuclid(cf,cg,M,c,fail);
504 if(fail)
505 return;
506#endif
507 // f /= cf
508 f.tryDiv (cf, M, fail);
509 if(fail)
510 return;
511 // g /= cg
512 g.tryDiv (cg, M, fail);
513 if(fail)
514 return;
515 TIMING_END_AND_PRINT (alg_content_p, "time for content in alg gcd mod p: ")
516 if(f.inCoeffDomain())
517 {
518 tryInvert(f,M,result,fail);
519 if(fail)
520 return;
521 result = NN(c);
522 return;
523 }
524 if(g.inCoeffDomain())
525 {
526 tryInvert(g,M,result,fail);
527 if(fail)
528 return;
529 result = NN(c);
530 return;
531 }
532 int *L = new int[mv+1]; // L is addressed by i from 2 to mv
533 int *N = new int[mv+1];
534 for(int i=2; i<=mv; i++)
535 L[i] = N[i] = 0;
536 L = leadDeg(f, L);
537 N = leadDeg(g, N);
538 CanonicalForm gamma;
539 TIMING_START (alg_euclid_p)
540#ifdef HAVE_NTL
541 NTLF= convertFacCF2NTLzz_pEX (firstLC (f), NTLMipo);
542 NTLG= convertFacCF2NTLzz_pEX (firstLC (g), NTLMipo);
543 tryNTLGCD (NTLResult, NTLF, NTLG, fail);
544 if (fail)
545 return;
546 gamma= convertNTLzz_pEX2CF (NTLResult, v1, v);
547#else
548 tryEuclid( firstLC(f), firstLC(g), M, gamma, fail );
549 if(fail)
550 return;
551#endif
552 TIMING_END_AND_PRINT (alg_euclid_p, "time for gcd of lcs in alg mod p: ")
553 for(int i=2; i<=mv; i++) // entries at i=0,1 not visited
554 if(N[i] < L[i])
555 L[i] = N[i];
556 // L is now upper bound for degrees of gcd
557 int *dg_im = new int[mv+1]; // for the degree vector of the image we don't need any entry at i=1
558 for(int i=2; i<=mv; i++)
559 dg_im[i] = 0; // initialize
560 CanonicalForm gamma_image, m=1;
561 CanonicalForm gm=0;
562 CanonicalForm g_image, alpha, gnew;
563 FFGenerator gen = FFGenerator();
564 Variable x= Variable (1);
565 bool divides= true;
566 for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next())
567 {
568 alpha = gen.item();
569 gamma_image = reduce(gamma(alpha, x),M); // plug in alpha for var(1)
570 if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
571 continue;
572 TIMING_START (alg_recursion_p)
573 tryBrownGCD( f(alpha, x), g(alpha, x), M, g_image, fail, false ); // recursive call with one var less
574 TIMING_END_AND_PRINT (alg_recursion_p,
575 "time for recursive calls in alg gcd mod p: ")
576 if(fail)
577 return;
578 g_image = reduce(g_image, M);
579 if(g_image.inCoeffDomain()) // early termination
580 {
581 tryInvert(g_image,M,result,fail);
582 if(fail)
583 return;
584 result = NN(c);
585 return;
586 }
587 for(int i=2; i<=mv; i++)
588 dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
589 dg_im = leadDeg(g_image, dg_im); // dg_im cannot be NIL-pointer
590 if(isEqual(dg_im, L, 2, mv))
591 {
592 CanonicalForm inv;
593 tryInvert (firstLC (g_image), M, inv, fail);
594 if (fail)
595 return;
596 g_image *= inv;
597 g_image *= gamma_image; // multiply by multiple of image lc(gcd)
598 g_image= reduce (g_image, M);
599 TIMING_START (alg_newton_p)
600 gnew= tryNewtonInterp (alpha, g_image, m, gm, x, M, fail);
601 TIMING_END_AND_PRINT (alg_newton_p,
602 "time for Newton interpolation in alg gcd mod p: ")
603 // gnew = gm mod m
604 // gnew = g_image mod var(1)-alpha
605 // mnew = m * (var(1)-alpha)
606 if(fail)
607 return;
608 m *= (x - alpha);
609 if((firstLC(gnew) == gamma) || (gnew == gm)) // gnew did not change
610 {
611 TIMING_START (alg_termination_p)
612 cf = tryvcontent(gnew, Variable(2), M, fail);
613 if(fail)
614 return;
615 divides = true;
616 g_image= gnew;
617 g_image.tryDiv (cf, M, fail);
618 if(fail)
619 return;
620 divides= tryFdivides (g_image,f, M, fail); // trial division (f)
621 if(fail)
622 return;
623 if(divides)
624 {
625 bool divides2= tryFdivides (g_image,g, M, fail); // trial division (g)
626 if(fail)
627 return;
628 if(divides2)
629 {
630 result = NN(reduce (c*g_image, M));
631 TIMING_END_AND_PRINT (alg_termination_p,
632 "time for successful termination test in alg gcd mod p: ")
633 return;
634 }
635 }
636 TIMING_END_AND_PRINT (alg_termination_p,
637 "time for unsuccessful termination test in alg gcd mod p: ")
638 }
639 gm = gnew;
640 continue;
641 }
642
643 if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky
644 continue;
645
646 // here: isLess(dg_im, L, 2, mv) --> all previous points were unlucky
647 m = CanonicalForm(1); // reset
648 gm = 0; // reset
649 for(int i=2; i<=mv; i++) // tighten bound
650 L[i] = dg_im[i];
651 }
652 // we are out of evaluation points
653 fail = true;
654}
zz_pEX convertFacCF2NTLzz_pEX(const CanonicalForm &f, const zz_pX &mipo)
Definition: NTLconvert.cc:1064
CanonicalForm convertNTLzz_pEX2CF(const zz_pEX &f, const Variable &x, const Variable &alpha)
Definition: NTLconvert.cc:1092
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:105
VAR long fac_NTL_char
Definition: NTLconvert.cc:46
CanonicalForm reduce(const CanonicalForm &f, const CanonicalForm &M)
polynomials in M.mvar() are considered coefficients M univariate monic polynomial the coefficients of...
Definition: cf_ops.cc:660
int level(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
static CanonicalForm tryNewtonInterp(const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x, const CanonicalForm &M, bool &fail)
Definition: cfGcdAlgExt.cc:357
const CanonicalForm CFMap CFMap bool topLevel
Definition: cfGcdAlgExt.cc:57
static CanonicalForm tryvcontent(const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
CanonicalForm firstLC(const CanonicalForm &f)
Definition: cfGcdAlgExt.cc:957
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:93
Variable x
Definition: cfModGcd.cc:4084
CanonicalForm cf
Definition: cfModGcd.cc:4085
CanonicalForm cg
Definition: cfModGcd.cc:4085
void tryNTLGCD(zz_pEX &x, const zz_pEX &a, const zz_pEX &b, bool &fail)
compute the GCD x of a and b, fail is set to true if a zero divisor is encountered
bool tryFdivides(const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
same as fdivides but handles zero divisors in Z_p[t]/(f)[x1,...,xn] for reducible f
class CFMap
Definition: cf_map.h:85
CanonicalForm & tryDiv(const CanonicalForm &, const CanonicalForm &, bool &)
same as divremt but handles zero divisors in case we are in Z_p[x]/(f) where f is not irreducible
generate all elements in F_p starting from 0
Definition: cf_generator.h:56
Variable alpha
Definition: facAbsBiFact.cc:51
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
ListNode * next
Definition: janet.h:31
void init()
Definition: lintree.cc:864

◆ trycf_content()

static CanonicalForm trycf_content ( const CanonicalForm f,
const CanonicalForm g,
const CanonicalForm M,
bool &  fail 
)
static

Definition at line 1068 of file cfGcdAlgExt.cc.

1069{ // as cf_content, but takes care of zero divisors
1070 if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
1071 {
1072 CFIterator i = f;
1073 CanonicalForm tmp = g, result;
1074 while ( i.hasTerms() && ! tmp.isOne() && ! fail )
1075 {
1076 tryBrownGCD( i.coeff(), tmp, M, result, fail );
1077 tmp = result;
1078 i++;
1079 }
1080 return result;
1081 }
1082 return abs( f );
1083}
void tryBrownGCD(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M, CanonicalForm &result, bool &fail, bool topLevel)
modular gcd over F_p[x]/(M) for not necessarily irreducible M. If a zero divisor is encountered fail ...
Definition: cfGcdAlgExt.cc:372
bool getReduce(const Variable &alpha)
Definition: variable.cc:232

◆ trycontent()

static CanonicalForm trycontent ( const CanonicalForm f,
const Variable x,
const CanonicalForm M,
bool &  fail 
)
static

Definition at line 1037 of file cfGcdAlgExt.cc.

1038{ // as 'content', but takes care of zero divisors
1039 ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
1040 Variable y = f.mvar();
1041 if ( y == x )
1042 return trycf_content( f, 0, M, fail );
1043 if ( y < x )
1044 return f;
1045 return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
1046}
CanonicalForm swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
static CanonicalForm trycontent(const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
static CanonicalForm trycf_content(const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
#define ASSERT(expression, message)
Definition: cf_assert.h:99
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53

◆ tryInvert()

void tryInvert ( const CanonicalForm F,
const CanonicalForm M,
CanonicalForm inv,
bool &  fail 
)

Definition at line 221 of file cfGcdAlgExt.cc.

222{ // F, M are required to be "univariate" polynomials in an algebraic variable
223 // we try to invert F modulo M
224 if(F.inBaseDomain())
225 {
226 if(F.isZero())
227 {
228 fail = true;
229 return;
230 }
231 inv = 1/F;
232 return;
233 }
235 Variable a = M.mvar();
236 Variable x = Variable(1);
237 if(!extgcd( replacevar( F, a, x ), replacevar( M, a, x ), inv, b ).isOne())
238 fail = true;
239 else
240 inv = replacevar( inv, x, a ); // change back to alg var
241}
CanonicalForm FACTORY_PUBLIC replacevar(const CanonicalForm &, const Variable &, const Variable &)
CanonicalForm replacevar ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
Definition: cf_ops.cc:271
CanonicalForm extgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
Definition: cfUnivarGcd.cc:174

◆ tryNewtonInterp()

static CanonicalForm tryNewtonInterp ( const CanonicalForm alpha,
const CanonicalForm u,
const CanonicalForm newtonPoly,
const CanonicalForm oldInterPoly,
const Variable x,
const CanonicalForm M,
bool &  fail 
)
inlinestatic

Definition at line 357 of file cfGcdAlgExt.cc.

360{
361 CanonicalForm interPoly;
362
363 CanonicalForm inv;
364 tryInvert (newtonPoly (alpha, x), M, inv, fail);
365 if (fail)
366 return 0;
367
368 interPoly= oldInterPoly+reduce ((u - oldInterPoly (alpha, x))*inv*newtonPoly, M);
369 return interPoly;
370}

◆ tryvcontent()

static CanonicalForm tryvcontent ( const CanonicalForm f,
const Variable x,
const CanonicalForm M,
bool &  fail 
)
static

Definition at line 1049 of file cfGcdAlgExt.cc.

1050{ // as vcontent, but takes care of zero divisors
1051 ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
1052 if ( f.mvar() <= x )
1053 return trycontent( f, x, M, fail );
1054 CFIterator i;
1055 CanonicalForm d = 0, e, ret;
1056 for ( i = f; i.hasTerms() && ! d.isOne() && ! fail; i++ )
1057 {
1058 e = tryvcontent( i.coeff(), x, M, fail );
1059 if(fail)
1060 break;
1061 tryBrownGCD( d, e, M, ret, fail );
1062 d = ret;
1063 }
1064 return d;
1065}

Variable Documentation

◆ both_non_zero

int both_non_zero = 0

Definition at line 68 of file cfGcdAlgExt.cc.

◆ both_zero

int both_zero = 0

Definition at line 71 of file cfGcdAlgExt.cc.

◆ degsf

degsf = NEW_ARRAY(int,n + 1)

Definition at line 59 of file cfGcdAlgExt.cc.

◆ degsg

degsg = NEW_ARRAY(int,n + 1)

Definition at line 60 of file cfGcdAlgExt.cc.

◆ else

else
Initial value:
{
for (int i= 1; i <= n; i++)
{
if (degsf[i] == 0 && degsg[i] == 0)
{
continue;
}
else
{
if (both_zero != 0)
{
M.newpair (Variable (i), Variable (i - both_zero));
N.newpair (Variable (i - both_zero), Variable (i));
}
}
}
}
int both_zero
Definition: cfGcdAlgExt.cc:71

Definition at line 194 of file cfGcdAlgExt.cc.

◆ f_zero

int f_zero = 0

Definition at line 69 of file cfGcdAlgExt.cc.

◆ Flevel

int Flevel =F.level()

Definition at line 72 of file cfGcdAlgExt.cc.

◆ G

Definition at line 55 of file cfGcdAlgExt.cc.

◆ g_zero

int g_zero = 0

Definition at line 70 of file cfGcdAlgExt.cc.

◆ Glevel

int Glevel =G.level()

Definition at line 73 of file cfGcdAlgExt.cc.

◆ M

Definition at line 55 of file cfGcdAlgExt.cc.

◆ N

Definition at line 56 of file cfGcdAlgExt.cc.

◆ return

return

Definition at line 218 of file cfGcdAlgExt.cc.

◆ topLevel

const CanonicalForm CFMap CFMap bool topLevel
Initial value:
{
int n= tmax (F.level(), G.level())

Definition at line 56 of file cfGcdAlgExt.cc.