My Project
Loading...
Searching...
No Matches
Functions
cfGcdAlgExt.h File Reference

GCD over Q(a) More...

#include "canonicalform.h"
#include "variable.h"

Go to the source code of this file.

Functions

CanonicalForm QGCD (const CanonicalForm &, const CanonicalForm &)
 gcd over Q(a)
 
void tryInvert (const CanonicalForm &, const CanonicalForm &, CanonicalForm &, bool &)
 
void tryBrownGCD (const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M, CanonicalForm &result, bool &fail, bool topLevel=true)
 modular gcd over F_p[x]/(M) for not necessarily irreducible M. If a zero divisor is encountered fail is set to true.
 
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)
 

Detailed Description

GCD over Q(a)

ABSTRACT: Implementation of Encarnacion's GCD algorithm over number fields, see M.J. Encarnacion "Computing GCDs of polynomials over number fields", extended to the multivariate case.

See also
cfNTLzzpEXGCD.h

Definition in file cfGcdAlgExt.h.

Function Documentation

◆ firstLC()

Definition at line 955 of file cfGcdAlgExt.cc.

956{ // returns the leading coefficient (LC) of level <= 1
958 while(ret.level() > 1)
959 ret = LC(ret);
960 return ret;
961}
CanonicalForm LC(const CanonicalForm &f)
FILE * f
Definition checklibs.c:9
factory's main class

◆ isEqual()

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

Definition at line 946 of file cfGcdAlgExt.cc.

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

◆ isLess()

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

Definition at line 935 of file cfGcdAlgExt.cc.

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

◆ QGCD()

gcd over Q(a)

Definition at line 730 of file cfGcdAlgExt.cc.

731{ // f,g in Q(a)[x1,...,xn]
732 if(F.isZero())
733 {
734 if(G.isZero())
735 return G; // G is zero
736 if(G.inCoeffDomain())
737 return CanonicalForm(1);
738 CanonicalForm lcinv= 1/Lc (G);
739 return G*lcinv; // return monic G
740 }
741 if(G.isZero()) // F is non-zero
742 {
743 if(F.inCoeffDomain())
744 return CanonicalForm(1);
745 CanonicalForm lcinv= 1/Lc (F);
746 return F*lcinv; // return monic F
747 }
748 if(F.inCoeffDomain() || G.inCoeffDomain())
749 return CanonicalForm(1);
750 // here: both NOT inCoeffDomain
751 CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
752 int p, i;
753 int *bound, *other; // degree vectors
754 bool fail;
756 On( SW_RATIONAL ); // needed by bCommonDen
757 f = F * bCommonDen(F);
758 g = G * bCommonDen(G);
762 f /= contf;
763 g /= contg;
765 TIMING_END_AND_PRINT (alg_content, "time for content in alg gcd: ")
766 Variable a, b;
768 {
769 if(hasFirstAlgVar(g,b))
770 {
771 if(b.level() > a.level())
772 a = b;
773 }
774 }
775 else
776 {
777 if(!hasFirstAlgVar(g,a))// both not in extension
778 {
779 Off( SW_RATIONAL );
780 Off( SW_USE_QGCD );
781 tmp = gcdcfcg*gcd( f, g );
782 On( SW_USE_QGCD );
784 return tmp;
785 }
786 }
787 // here: a is the biggest alg. var in f and g AND some of f,g is in extension
788 setReduce(a,false); // do not reduce expressions modulo mipo
789 tmp = getMipo(a);
790 M = tmp * bCommonDen(tmp);
791 // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
792 Off( SW_RATIONAL ); // needed by mod
793 // calculate upper bound for degree vector of gcd
794 int mv = f.level(); i = g.level();
795 if(i > mv)
796 mv = i;
797 // here: mv is level of the largest variable in f, g
798 bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
799 other = new int[mv+1];
800 for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
801 bound[i] = other[i] = 0;
802 leadDeg(f,bound); // 'bound' is set the leading degree vector of f
803 leadDeg(g,other);
804 for(int i=1; i<=mv; i++) // entry at i=0 not visited
805 if(other[i] < bound[i])
806 bound[i] = other[i];
807 // now 'bound' is the smaller vector
808 cl = lc(M) * lc(f) * lc(g);
809 q = 1;
810 D = 0;
812 bool equal= false;
813 for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
814 {
815 p = cf_getBigPrime(i);
816 if( mod( cl, p ).isZero() ) // skip lc-bad primes
817 continue;
818 fail = false;
820 mipo = mapinto(M);
821 mipo /= mipo.lc();
822 // here: mipo is monic
824 tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
825 TIMING_END_AND_PRINT (alg_gcd_p, "time for alg gcd mod p: ")
826 if( fail ) // mipo splits in char p
827 {
829 continue;
830 }
831 if( Dp.inCoeffDomain() ) // early termination
832 {
833 tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
835 if(fail)
836 continue;
837 setReduce(a,true);
839 delete[] other;
840 delete[] bound;
841 return gcdcfcg;
842 }
844 // here: Dp NOT inCoeffDomain
845 for(int i=1; i<=mv; i++)
846 other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
848
849 if(isEqual(bound, other, 1, mv)) // equal
850 {
851 chineseRemainder( D, q, mapinto(Dp), p, tmp, newq );
852 // tmp = Dp mod p
853 // tmp = D mod q
854 // newq = p*q
855 q = newq;
856 if( D != tmp )
857 D = tmp;
858 On( SW_RATIONAL );
860 tmp = Farey( D, q ); // Farey
861 tmp *= bCommonDen (tmp);
863 "time for rational reconstruction in alg gcd: ")
864 setReduce(a,true); // reduce expressions modulo mipo
865 On( SW_RATIONAL ); // needed by fdivides
866 if (test != tmp)
867 test= tmp;
868 else
869 equal= true; // modular image did not add any new information
871#ifdef HAVE_NTL
872#ifdef HAVE_FLINT
873 if (equal && tmp.isUnivariate() && f.isUnivariate() && g.isUnivariate()
874 && f.level() == tmp.level() && tmp.level() == g.level())
875 {
877 newtonDivrem (f, tmp, Q, R);
878 if (R.isZero())
879 {
880 newtonDivrem (g, tmp, Q, R);
881 if (R.isZero())
882 {
884 setReduce (a,true);
887 "time for successful termination test in alg gcd: ")
888 delete[] other;
889 delete[] bound;
891 }
892 }
893 }
894 else
895#endif
896#endif
897 if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
898 {
899 Off( SW_RATIONAL );
900 setReduce(a,true);
903 "time for successful termination test in alg gcd: ")
904 delete[] other;
905 delete[] bound;
907 }
910 Off( SW_RATIONAL );
911 setReduce(a,false); // do not reduce expressions modulo mipo
912 continue;
913 }
914 if( isLess(bound, other, 1, mv) ) // current prime unlucky
915 continue;
916 // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
917 q = p;
918 D = mapinto(Dp); // shortcut CRA // shortcut CRA
919 for(int i=1; i<=mv; i++) // tighten bound
920 bound[i] = other[i];
921 }
922 // hopefully, we never reach this point
923 setReduce(a,true);
924 delete[] other;
925 delete[] bound;
926 Off( SW_USE_QGCD );
927 D = gcdcfcg*gcd( f, g );
928 On( SW_USE_QGCD );
930 return D;
931}
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)
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)
static CanonicalForm myicontent(const CanonicalForm &f, const CanonicalForm &c)
if(topLevel)
return
const CanonicalForm & G
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
const CanonicalForm CFMap & M
int p
Definition cfModGcd.cc:4086
return false
Definition cfModGcd.cc:85
cl
Definition cfModGcd.cc:4108
g
Definition cfModGcd.cc:4098
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
CF_NO_INLINE bool isZero() const
bool inCoeffDomain() const
int level() const
level() returns the level of CO.
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
factory's class for variables
Definition factory.h:127
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
for(j=0;j< factors.length();j++)
Definition facHensel.cc:129
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:351
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
static number Farey(number, number, const coeffs)
Definition flintcf_Q.cc:494
#define D(A)
Definition gentable.cc:131
#define R
Definition sirandom.c:27
#define Q
Definition sirandom.c:26
#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)

◆ tryBrownGCD()

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

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

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

◆ 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,...
bool inBaseDomain() const