My Project
Loading...
Searching...
No Matches
cfGcdAlgExt.cc File Reference
#include "config.h"
#include <vector>
#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
 
 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)
 
static void leadDeg (const CanonicalForm &f, int degs[])
 
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.
 
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)
 
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()

Definition at line 956 of file cfGcdAlgExt.cc.

957{ // returns the leading coefficient (LC) of level <= 1
958 CanonicalForm ret = f;
959 while(ret.level() > 1)
960 ret = LC(ret);
961 return ret;
962}
CanonicalForm LC(const CanonicalForm &f)
FILE * f
Definition checklibs.c:9
factory's main class
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 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
STATIC_VAR TreeM * G
Definition janet.cc:31

◆ if()

if ( topLevel )

Definition at line 76 of file cfGcdAlgExt.cc.

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

◆ isEqual()

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

Definition at line 947 of file cfGcdAlgExt.cc.

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

◆ isLess()

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

Definition at line 936 of file cfGcdAlgExt.cc.

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

◆ leadDeg()

static void leadDeg ( const CanonicalForm & f,
int degs[] )
static

Definition at line 373 of file cfGcdAlgExt.cc.

374{ // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
375 // 'a' should point to an array of sufficient size level(f)+1
376 if(f.inCoeffDomain())
377 return;
378 CanonicalForm tmp = f;
379 do
380 {
381 degs[tmp.level()] = tmp.degree();
382 tmp = LC(tmp);
383 }
384 while(!tmp.inCoeffDomain());
385}
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 722 of file cfGcdAlgExt.cc.

723{
724#if defined(HAVE_NTL) || defined(HAVE_FLINT)
725 return myicontent( f, 0 );
726#else
727 return 1;
728#endif
729}
static CanonicalForm myicontent(const CanonicalForm &f, const CanonicalForm &c)

◆ myicontent() [2/2]

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

Definition at line 672 of file cfGcdAlgExt.cc.

673{
674#if defined(HAVE_NTL) || defined(HAVE_FLINT)
675 if (f.isOne() || c.isOne())
676 return 1;
677 if ( f.inBaseDomain() && c.inBaseDomain())
678 {
679 if (c.isZero()) return abs(f);
680 return bgcd( f, c );
681 }
682 else if ( (f.inCoeffDomain() && c.inCoeffDomain()) ||
683 (f.inCoeffDomain() && c.inBaseDomain()) ||
684 (f.inBaseDomain() && c.inCoeffDomain()))
685 {
686 if (c.isZero()) return abs (f);
687#ifdef HAVE_FLINT
688 fmpz_poly_t FLINTf, FLINTc;
689 convertFacCF2Fmpz_poly_t (FLINTf, f);
690 convertFacCF2Fmpz_poly_t (FLINTc, c);
691 fmpz_poly_gcd (FLINTc, FLINTc, FLINTf);
693 if (f.inCoeffDomain())
694 result= convertFmpz_poly_t2FacCF (FLINTc, f.mvar());
695 else
696 result= convertFmpz_poly_t2FacCF (FLINTc, c.mvar());
697 fmpz_poly_clear (FLINTc);
698 fmpz_poly_clear (FLINTf);
699 return result;
700#else
701 ZZX NTLf= convertFacCF2NTLZZX (f);
702 ZZX NTLc= convertFacCF2NTLZZX (c);
703 NTLc= GCD (NTLc, NTLf);
704 if (f.inCoeffDomain())
705 return convertNTLZZX2CF(NTLc,f.mvar());
706 else
707 return convertNTLZZX2CF(NTLc,c.mvar());
708#endif
709 }
710 else
711 {
712 CanonicalForm g = c;
713 for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
714 g = myicontent( i.coeff(), g );
715 return g;
716 }
717#else
718 return 1;
719#endif
720}
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
Rational abs(const Rational &a)
Definition GMPrat.cc:436
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
CanonicalForm bgcd(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm bgcd ( const CanonicalForm & f, const CanonicalForm & g )
g
Definition cfModGcd.cc:4098
class to iterate through CanonicalForm's
Definition cf_iter.h:44
CF_NO_INLINE bool isZero() const
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
CF_NO_INLINE bool isOne() const
bool inBaseDomain() const
return result

◆ QGCD()

gcd over Q(a)

Definition at line 731 of file cfGcdAlgExt.cc.

732{ // f,g in Q(a)[x1,...,xn]
733 if(F.isZero())
734 {
735 if(G.isZero())
736 return G; // G is zero
737 if(G.inCoeffDomain())
738 return CanonicalForm(1);
739 CanonicalForm lcinv= 1/Lc (G);
740 return G*lcinv; // return monic G
741 }
742 if(G.isZero()) // F is non-zero
743 {
744 if(F.inCoeffDomain())
745 return CanonicalForm(1);
746 CanonicalForm lcinv= 1/Lc (F);
747 return F*lcinv; // return monic F
748 }
749 if(F.inCoeffDomain() || G.inCoeffDomain())
750 return CanonicalForm(1);
751 // here: both NOT inCoeffDomain
752 CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
753 int p, i;
754 int *bound, *other; // degree vectors
755 bool fail;
756 bool off_rational=!isOn(SW_RATIONAL);
757 On( SW_RATIONAL ); // needed by bCommonDen
758 f = F * bCommonDen(F);
759 g = G * bCommonDen(G);
761 CanonicalForm contf= myicontent (f);
762 CanonicalForm contg= myicontent (g);
763 f /= contf;
764 g /= contg;
765 CanonicalForm gcdcfcg= myicontent (contf, contg);
766 TIMING_END_AND_PRINT (alg_content, "time for content in alg gcd: ")
767 Variable a, b;
769 {
770 if(hasFirstAlgVar(g,b))
771 {
772 if(b.level() > a.level())
773 a = b;
774 }
775 }
776 else
777 {
778 if(!hasFirstAlgVar(g,a))// both not in extension
779 {
780 Off( SW_RATIONAL );
781 Off( SW_USE_QGCD );
782 tmp = gcdcfcg*gcd( f, g );
783 On( SW_USE_QGCD );
784 if (off_rational) Off(SW_RATIONAL);
785 return tmp;
786 }
787 }
788 // here: a is the biggest alg. var in f and g AND some of f,g is in extension
789 setReduce(a,false); // do not reduce expressions modulo mipo
790 tmp = getMipo(a);
791 M = tmp * bCommonDen(tmp);
792 // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
793 Off( SW_RATIONAL ); // needed by mod
794 // calculate upper bound for degree vector of gcd
795 int mv = f.level(); i = g.level();
796 if(i > mv)
797 mv = i;
798 // here: mv is level of the largest variable in f, g
799 bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
800 other = new int[mv+1];
801 for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
802 bound[i] = other[i] = 0;
803 leadDeg(f,bound); // 'bound' is set the leading degree vector of f
804 leadDeg(g,other);
805 for(int i=1; i<=mv; i++) // entry at i=0 not visited
806 if(other[i] < bound[i])
807 bound[i] = other[i];
808 // now 'bound' is the smaller vector
809 cl = lc(M) * lc(f) * lc(g);
810 q = 1;
811 D = 0;
813 bool equal= false;
814 for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
815 {
816 p = cf_getBigPrime(i);
817 if( mod( cl, p ).isZero() ) // skip lc-bad primes
818 continue;
819 fail = false;
821 mipo = mapinto(M);
822 mipo /= mipo.lc();
823 // here: mipo is monic
824 TIMING_START (alg_gcd_p)
825 tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
826 TIMING_END_AND_PRINT (alg_gcd_p, "time for alg gcd mod p: ")
827 if( fail ) // mipo splits in char p
828 {
830 continue;
831 }
832 if( Dp.inCoeffDomain() ) // early termination
833 {
834 tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
836 if(fail)
837 continue;
838 setReduce(a,true);
839 if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
840 delete[] other;
841 delete[] bound;
842 return gcdcfcg;
843 }
845 // here: Dp NOT inCoeffDomain
846 for(int i=1; i<=mv; i++)
847 other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
848 leadDeg(Dp,other);
849
850 if(isEqual(bound, other, 1, mv)) // equal
851 {
852 chineseRemainder( D, q, mapinto(Dp), p, tmp, newq );
853 // tmp = Dp mod p
854 // tmp = D mod q
855 // newq = p*q
856 q = newq;
857 if( D != tmp )
858 D = tmp;
859 On( SW_RATIONAL );
860 TIMING_START (alg_reconstruction)
861 tmp = Farey( D, q ); // Farey
862 tmp *= bCommonDen (tmp);
863 TIMING_END_AND_PRINT (alg_reconstruction,
864 "time for rational reconstruction in alg gcd: ")
865 setReduce(a,true); // reduce expressions modulo mipo
866 On( SW_RATIONAL ); // needed by fdivides
867 if (test != tmp)
868 test= tmp;
869 else
870 equal= true; // modular image did not add any new information
871 TIMING_START (alg_termination)
872#ifdef HAVE_NTL
873#ifdef HAVE_FLINT
874 if (equal && tmp.isUnivariate() && f.isUnivariate() && g.isUnivariate()
875 && f.level() == tmp.level() && tmp.level() == g.level())
876 {
877 CanonicalForm Q, R;
878 newtonDivrem (f, tmp, Q, R);
879 if (R.isZero())
880 {
881 newtonDivrem (g, tmp, Q, R);
882 if (R.isZero())
883 {
884 Off (SW_RATIONAL);
885 setReduce (a,true);
886 if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
887 TIMING_END_AND_PRINT (alg_termination,
888 "time for successful termination test in alg gcd: ")
889 delete[] other;
890 delete[] bound;
891 return tmp*gcdcfcg;
892 }
893 }
894 }
895 else
896#endif
897#endif
898 if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
899 {
900 Off( SW_RATIONAL );
901 setReduce(a,true);
902 if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
903 TIMING_END_AND_PRINT (alg_termination,
904 "time for successful termination test in alg gcd: ")
905 delete[] other;
906 delete[] bound;
907 return tmp*gcdcfcg;
908 }
909 TIMING_END_AND_PRINT (alg_termination,
910 "time for unsuccessful termination test in alg gcd: ")
911 Off( SW_RATIONAL );
912 setReduce(a,false); // do not reduce expressions modulo mipo
913 continue;
914 }
915 if( isLess(bound, other, 1, mv) ) // current prime unlucky
916 continue;
917 // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
918 q = p;
919 D = mapinto(Dp); // shortcut CRA // shortcut CRA
920 for(int i=1; i<=mv; i++) // tighten bound
921 bound[i] = other[i];
922 }
923 // hopefully, we never reach this point
924 setReduce(a,true);
925 delete[] other;
926 delete[] bound;
927 Off( SW_USE_QGCD );
928 D = gcdcfcg*gcd( f, g );
929 On( SW_USE_QGCD );
930 if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
931 return D;
932}
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)
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
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)
for(int i=0;i<=n;i++) degsf[i]
Definition cfEzgcd.cc:72
else
bool isLess(int *a, int *b, int lower, int upper)
static void leadDeg(const CanonicalForm &f, int degs[])
bool isEqual(int *a, int *b, int lower, int upper)
return
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
int p
Definition cfModGcd.cc:4086
return false
Definition cfModGcd.cc:85
cl
Definition cfModGcd.cc:4108
CanonicalForm gcdcfcg
Definition cfModGcd.cc:4109
bool equal
Definition cfModGcd.cc:4134
CanonicalForm test
Definition cfModGcd.cc:4104
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 )
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:43
static const int SW_RATIONAL
set to 1 for computations over Q
Definition cf_defs.h:31
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
int level() const
Definition factory.h:143
CanonicalForm mipo
Definition facAlgExt.cc:57
CanonicalForm alg_content(const CanonicalForm &f, const CFList &as)
Definition facAlgFunc.cc:42
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
if(!FE_OPT_NO_SHELL_FLAG)
Definition fehelp.cc:1000
static number Farey(number, number, const coeffs)
Definition flintcf_Q.cc:495
#define D(A)
Definition gentable.cc:128
#define TIMING_START(t)
Definition timing.h:92
#define TIMING_END_AND_PRINT(t, msg)
Definition timing.h:94
int gcd(int a, int b)

◆ 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 387 of file cfGcdAlgExt.cc.

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

◆ trycf_content()

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

Definition at line 1067 of file cfGcdAlgExt.cc.

1068{ // as cf_content, but takes care of zero divisors
1069 if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
1070 {
1071 CFIterator i = f;
1072 CanonicalForm tmp = g, result;
1073 while ( i.hasTerms() && ! tmp.isOne() && ! fail )
1074 {
1075 tryBrownGCD( i.coeff(), tmp, M, result, fail );
1076 tmp = result;
1077 i++;
1078 }
1079 return result;
1080 }
1081 return abs( f );
1082}
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 ...
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 1036 of file cfGcdAlgExt.cc.

1037{ // as 'content', but takes care of zero divisors
1038 ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
1039 Variable y = f.mvar();
1040 if ( y == x )
1041 return trycf_content( f, 0, M, fail );
1042 if ( y < x )
1043 return f;
1044 return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
1045}
CanonicalForm FACTORY_PUBLIC 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 222 of file cfGcdAlgExt.cc.

223{ // F, M are required to be "univariate" polynomials in an algebraic variable
224 // we try to invert F modulo M
225 if(F.inBaseDomain())
226 {
227 if(F.isZero())
228 {
229 fail = true;
230 return;
231 }
232 inv = 1/F;
233 return;
234 }
236 Variable a = M.mvar();
237 Variable x = Variable(1);
238 if(!extgcd( replacevar( F, a, x ), replacevar( M, a, x ), inv, b ).isOne())
239 fail = true;
240 else
241 inv = replacevar( inv, x, a ); // change back to alg var
242}
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,...

◆ 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 358 of file cfGcdAlgExt.cc.

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

◆ tryvcontent()

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

Definition at line 1048 of file cfGcdAlgExt.cc.

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

Variable Documentation

◆ both_non_zero

int both_non_zero = 0

Definition at line 69 of file cfGcdAlgExt.cc.

◆ both_zero

int both_zero = 0

Definition at line 72 of file cfGcdAlgExt.cc.

◆ degsf

degsf = NEW_ARRAY(int,n + 1)

Definition at line 60 of file cfGcdAlgExt.cc.

◆ degsg

degsg = NEW_ARRAY(int,n + 1)

Definition at line 61 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 at line 195 of file cfGcdAlgExt.cc.

196 {
197 //arrange Variables such that no gaps occur
198 for (int i= 1; i <= n; i++)
199 {
200 if (degsf[i] == 0 && degsg[i] == 0)
201 {
202 both_zero++;
203 continue;
204 }
205 else
206 {
207 if (both_zero != 0)
208 {
209 M.newpair (Variable (i), Variable (i - both_zero));
210 N.newpair (Variable (i - both_zero), Variable (i));
211 }
212 }
213 }
214 }

◆ f_zero

int f_zero = 0

Definition at line 70 of file cfGcdAlgExt.cc.

◆ Flevel

int Flevel =F.level()

Definition at line 73 of file cfGcdAlgExt.cc.

◆ G

Definition at line 56 of file cfGcdAlgExt.cc.

◆ g_zero

int g_zero = 0

Definition at line 71 of file cfGcdAlgExt.cc.

◆ Glevel

int Glevel =G.level()

Definition at line 74 of file cfGcdAlgExt.cc.

◆ M

Definition at line 56 of file cfGcdAlgExt.cc.

◆ N

Definition at line 57 of file cfGcdAlgExt.cc.

◆ return

return

Definition at line 219 of file cfGcdAlgExt.cc.

◆ topLevel

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

Definition at line 57 of file cfGcdAlgExt.cc.