My Project
Functions | Variables
cfEzgcd.cc File Reference

This file implements the GCD of two multivariate polynomials over Q or F_q using EZ-GCD as described in "Algorithms for Computer Algebra" by Geddes, Czapor, Labahnn. More...

#include "config.h"
#include "timing.h"
#include "cf_assert.h"
#include "debug.h"
#include "cf_defs.h"
#include "canonicalform.h"
#include "cfEzgcd.h"
#include "cfModGcd.h"
#include "cf_util.h"
#include "cf_iter.h"
#include "cf_map_ext.h"
#include "cf_algorithm.h"
#include "cf_reval.h"
#include "cf_random.h"
#include "cf_primes.h"
#include "templates/ftmpl_functions.h"
#include "cf_map.h"
#include "facHensel.h"
#include "FLINTconvert.h"
#include "NTLconvert.h"

Go to the source code of this file.

Functions

 TIMING_DEFINE_PRINT (ez_eval) TIMING_DEFINE_PRINT(ez_compress) TIMING_DEFINE_PRINT(ez_hensel_lift) TIMING_DEFINE_PRINT(ez_content) TIMING_DEFINE_PRINT(ez_termination) static int compress4EZGCD(const CanonicalForm &F
 
 for (int i=0;i<=n;i++) degsf[i]
 
 if (both_non_zero==0)
 
 while (k > 0)
 
 DELETE_ARRAY (degsf)
 
 DELETE_ARRAY (degsg)
 
static CanonicalForm myShift2Zero (const CanonicalForm &F, CFList &Feval, const CFList &evaluation)
 
static CanonicalForm myReverseShift (const CanonicalForm &F, const CFList &evaluation)
 
static Evaluation optimize4Lift (const CanonicalForm &F, CFMap &M, CFMap &N, const Evaluation &A)
 
static int Hensel (const CanonicalForm &UU, CFArray &G, const Evaluation &AA, const CFArray &LeadCoeffs)
 
static bool findeval (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Fb, CanonicalForm &Gb, CanonicalForm &Db, REvaluation &b, int delta, int degF, int degG, int maxeval, int &count, int &k, int bound, int &l)
 
static void gcd_mon_rec (CanonicalForm G, CanonicalForm &cf, int *exp, int pl)
 
static CanonicalForm gcd_mon (CanonicalForm F, CanonicalForm G)
 
static CanonicalForm ezgcd (const CanonicalForm &FF, const CanonicalForm &GG, REvaluation &b, bool internal)
 real implementation of EZGCD over Z More...
 
CanonicalForm ezgcd (const CanonicalForm &FF, const CanonicalForm &GG)
 Extended Zassenhaus GCD over Z. In case things become too dense we switch to a modular algorithm. More...
 
CanonicalForm EZGCD_P (const CanonicalForm &FF, const CanonicalForm &GG)
 Extended Zassenhaus GCD for finite fields. In case things become too dense we switch to a modular algorithm. More...
 

Variables

static const double log2exp = 1.442695041
 
const CanonicalFormG
 
const CanonicalForm CFMapM
 
const CanonicalForm CFMap CFMapN
 
const CanonicalForm CFMap CFMap int & both_non_zero
 
int * degsf = NEW_ARRAY(int,n + 1)
 
int * degsg = NEW_ARRAY(int,n + 1)
 
int f_zero = 0
 
int g_zero = 0
 
int k = 1
 
int l = 1
 
int Flevel =F.level()
 
int Glevel =G.level()
 
int m = tmin (Flevel, Glevel)
 
int max_min_deg
 
int i = 1
 
STATIC_VAR int maxNumEval = 200
 
STATIC_VAR int sizePerVars1 = 500
 

Detailed Description

This file implements the GCD of two multivariate polynomials over Q or F_q using EZ-GCD as described in "Algorithms for Computer Algebra" by Geddes, Czapor, Labahnn.

Author
Martin Lee

Definition in file cfEzgcd.cc.

Function Documentation

◆ DELETE_ARRAY() [1/2]

DELETE_ARRAY ( degsf  )

◆ DELETE_ARRAY() [2/2]

DELETE_ARRAY ( degsg  )

◆ ezgcd() [1/2]

CanonicalForm ezgcd ( const CanonicalForm FF,
const CanonicalForm GG 
)

Extended Zassenhaus GCD over Z. In case things become too dense we switch to a modular algorithm.

Definition at line 861 of file cfEzgcd.cc.

862 {
863  REvaluation b;
864  return ezgcd( FF, GG, b, false );
865 }
static CanonicalForm ezgcd(const CanonicalForm &FF, const CanonicalForm &GG, REvaluation &b, bool internal)
real implementation of EZGCD over Z
Definition: cfEzgcd.cc:498
const CanonicalForm & GG
Definition: cfModGcd.cc:4076
CanonicalForm b
Definition: cfModGcd.cc:4103
class to generate random evaluation points
Definition: cf_reval.h:26

◆ ezgcd() [2/2]

static CanonicalForm ezgcd ( const CanonicalForm FF,
const CanonicalForm GG,
REvaluation b,
bool  internal 
)
static

real implementation of EZGCD over Z

—> A2

—> A3

—> A4

—> A5

—> A6

—> A7

—> A8 (gcdfound)

—> A9

Definition at line 498 of file cfEzgcd.cc.

500 {
501  bool isRat= isOn (SW_RATIONAL);
502 
503  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
504  int sizeF= size (FF);
505  int sizeG= size (GG);
506 
507 
508  if (sizeF==1)
509  {
510  Off(SW_USE_EZGCD);
512  On(SW_USE_EZGCD);
513  return result;
514  }
515  else if (sizeG==1)
516  {
517  Off(SW_USE_EZGCD);
519  On(SW_USE_EZGCD);
520  return result;
521  }
522  if (!isRat)
523  On (SW_RATIONAL);
524  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
525  {
526  Off(SW_USE_EZGCD);
527  CanonicalForm result=gcd( FF, GG );
528  On(SW_USE_EZGCD);
529  if (!isRat)
530  Off (SW_RATIONAL);
531  result /= icontent (result);
532  DEBDECLEVEL( cerr, "ezgcd" );
533  return result;
534  }
535 
536 
537  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
538  lcD, cand, contcand, result;
539  CFArray DD( 1, 2 ), lcDD( 1, 2 );
540  int degF, degG, delta, t, count, maxeval;
541  REvaluation bt;
542  int gcdfound = 0;
543  Variable x = Variable(1);
544  count= 0;
545  maxeval= 200;
546  int o, l;
547  o= 0;
548  l= 1;
549 
550  if (!isRat)
551  On (SW_RATIONAL);
552  F= FF*bCommonDen (FF);
553  G= GG*bCommonDen (GG);
554  if (!isRat)
555  Off (SW_RATIONAL);
556 
557  TIMING_START (ez_compress)
558  CFMap M,N;
559  int smallestDegLev;
560  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
561 
562  if (best_level == 0)
563  {
564  DEBDECLEVEL( cerr, "ezgcd" );
565  return G.genOne();
566  }
567 
568  F= M (F);
569  G= M (G);
570  TIMING_END_AND_PRINT (ez_compress, "time for compression in EZ: ")
571 
572  DEBINCLEVEL( cerr, "ezgcd" );
573  DEBOUTLN( cerr, "FF = " << FF );
574  DEBOUTLN( cerr, "GG = " << GG );
575  TIMING_START (ez_content)
576  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
577  DEBOUTLN( cerr, "f = " << f );
578  DEBOUTLN( cerr, "g = " << g );
579  F /= f; G /= g;
580  TIMING_END_AND_PRINT (ez_content, "time to extract content in EZ: ")
581  if ( F.isUnivariate() )
582  {
583  if ( G.isUnivariate() )
584  {
585  DEBDECLEVEL( cerr, "ezgcd" );
586  if(F.mvar()==G.mvar())
587  d*=gcd(F,G);
588  else
589  return N (d);
590  return N (d);
591  }
592  else
593  {
594  g= content (G,G.mvar());
595  return N(d*gcd(F,g));
596  }
597  }
598  if ( G.isUnivariate())
599  {
600  f= content (F,F.mvar());
601  return N(d*gcd(G,f));
602  }
603 
605  sizeF= size (F);
606  sizeG= size (G);
607 
608  if (!isRat)
609  On (SW_RATIONAL);
610  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
611  {
612  Off(SW_USE_EZGCD);
613  result=gcd( F, G );
614  On(SW_USE_EZGCD);
615  if (!isRat)
616  Off (SW_RATIONAL);
617  result /= icontent (result);
618  DEBDECLEVEL( cerr, "ezgcd" );
619  return N (d*result);
620  }
621 
622  int dummy= 0;
623  if ( gcd_test_one( F, G, false, dummy ) )
624  {
625  DEBDECLEVEL( cerr, "ezgcd" );
626  if (!isRat)
627  Off (SW_RATIONAL);
628  return N (d);
629  }
630  lcF = LC( F, x ); lcG = LC( G, x );
631  lcD = gcd( lcF, lcG );
632  delta = 0;
633  degF = degree( F, x ); degG = degree( G, x );
634  t = tmax( F.level(), G.level() );
635  if ( ! internal )
636  b = REvaluation( 2, t, IntRandom( 25 ) );
637  while ( ! gcdfound )
638  {
639  /// ---> A2
640  DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
641  DEBOUTLN( cerr, "F = " << F );
642  DEBOUTLN( cerr, "G = " << G );
643  TIMING_START (ez_eval)
644  if (!findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count,
645  o, 25, l))
646  {
647  Off(SW_USE_EZGCD);
648  result=gcd( F, G );
649  On(SW_USE_EZGCD);
650  if (!isRat)
651  Off (SW_RATIONAL);
652  DEBDECLEVEL( cerr, "ezgcd" );
653  result /= icontent (result);
654  return N (d*result);
655  }
656  TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ1: ")
657  DEBOUTLN( cerr, "found evaluation b = " << b );
658  DEBOUTLN( cerr, "F(b) = " << Fb );
659  DEBOUTLN( cerr, "G(b) = " << Gb );
660  DEBOUTLN( cerr, "D(b) = " << Db );
661  delta = degree( Db );
662  /// ---> A3
663  if (delta == degF)
664  {
665  if (degF <= degG && fdivides (F, G))
666  {
667  DEBDECLEVEL( cerr, "ezgcd" );
668  if (!isRat)
669  Off (SW_RATIONAL);
670  return N (d*F);
671  }
672  else
673  delta--;
674  }
675  else if (delta == degG)
676  {
677  if (degG <= degF && fdivides( G, F ))
678  {
679  DEBDECLEVEL( cerr, "ezgcd" );
680  if (!isRat)
681  Off (SW_RATIONAL);
682  return N (d*G);
683  }
684  else
685  delta--;
686  }
687  if ( delta == 0 )
688  {
689  DEBDECLEVEL( cerr, "ezgcd" );
690  if (!isRat)
691  Off (SW_RATIONAL);
692  return N (d);
693  }
694  /// ---> A4
695  //deltaold = delta;
696  while ( 1 )
697  {
698  bt = b;
699  TIMING_START (ez_eval)
700  if (!findeval( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count,
701  o, 25,l ))
702  {
703  Off(SW_USE_EZGCD);
704  result=gcd( F, G );
705  On(SW_USE_EZGCD);
706  if (!isRat)
707  Off (SW_RATIONAL);
708  DEBDECLEVEL( cerr, "ezgcd" );
709  result /= icontent (result);
710  return N (d*result);
711  }
712  TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ2: ")
713  int dd=degree( Dbt );
714  if ( dd /*degree( Dbt )*/ == 0 )
715  {
716  DEBDECLEVEL( cerr, "ezgcd" );
717  if (!isRat)
718  Off (SW_RATIONAL);
719  return N (d);
720  }
721  if ( dd /*degree( Dbt )*/ == delta )
722  break;
723  else if ( dd /*degree( Dbt )*/ < delta )
724  {
725  delta = dd /*degree( Dbt )*/;
726  b = bt;
727  Db = Dbt; Fb = Fbt; Gb = Gbt;
728  }
729  DEBOUTLN( cerr, "now after A4, delta = " << delta );
730  /// ---> A5
731  if (delta == degF)
732  {
733  if (degF <= degG && fdivides (F, G))
734  {
735  DEBDECLEVEL( cerr, "ezgcd" );
736  if (!isRat)
737  Off (SW_RATIONAL);
738  return N (d*F);
739  }
740  else
741  delta--;
742  }
743  else if (delta == degG)
744  {
745  if (degG <= degF && fdivides( G, F ))
746  {
747  DEBDECLEVEL( cerr, "ezgcd" );
748  if (!isRat)
749  Off (SW_RATIONAL);
750  return N (d*G);
751  }
752  else
753  delta--;
754  }
755  if ( delta == 0 )
756  {
757  DEBDECLEVEL( cerr, "ezgcd" );
758  if (!isRat)
759  Off (SW_RATIONAL);
760  return N (d);
761  }
762  }
763  if ( delta != degF && delta != degG )
764  {
765  /// ---> A6
766  bool B_is_F;
767  CanonicalForm xxx1, xxx2;
769  DD[1] = Fb / Db;
770  buf= Gb/Db;
771  xxx1 = gcd( DD[1], Db );
772  xxx2 = gcd( buf, Db );
773  if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
774  (size (F) <= size (G)))
775  || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
776  {
777  B = F;
778  DD[2] = Db;
779  lcDD[1] = lcF;
780  lcDD[2] = lcD;
781  B_is_F = true;
782  }
783  else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
784  (size (G) < size (F)))
785  || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
786  {
787  DD[1] = buf;
788  B = G;
789  DD[2] = Db;
790  lcDD[1] = lcG;
791  lcDD[2] = lcD;
792  B_is_F = false;
793  }
794  else
795  {
796  //special case
797  Off(SW_USE_EZGCD);
798  result=gcd( F, G );
799  On(SW_USE_EZGCD);
800  if (!isRat)
801  Off (SW_RATIONAL);
802  DEBDECLEVEL( cerr, "ezgcd" );
803  result /= icontent (result);
804  return N (d*result);
805  }
806  /// ---> A7
807  DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
808  DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
809  DEBOUTLN( cerr, "(hensel) B = " << B );
810  DEBOUTLN( cerr, "(hensel) lcB = " << LC( B, Variable(1) ) );
811  DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
812  DEBOUTLN( cerr, "(hensel) DD = " << DD );
813  DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
814  TIMING_START (ez_hensel_lift)
815  gcdfound= Hensel (B*lcD, DD, b, lcDD);
816  TIMING_END_AND_PRINT (ez_hensel_lift, "time to hensel lift in EZ: ")
817  DEBOUTLN( cerr, "(hensel finished) DD = " << DD );
818 
819  if (gcdfound == -1)
820  {
821  Off (SW_USE_EZGCD);
822  result= gcd (F,G);
823  On (SW_USE_EZGCD);
824  if (!isRat)
825  Off (SW_RATIONAL);
826  DEBDECLEVEL( cerr, "ezgcd" );
827  result /= icontent (result);
828  return N (d*result);
829  }
830 
831  if (gcdfound)
832  {
833  TIMING_START (ez_termination)
834  contcand= content (DD[2], Variable (1));
835  cand = DD[2] / contcand;
836  if (B_is_F)
837  gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
838  else
839  gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
840  TIMING_END_AND_PRINT (ez_termination,
841  "time for termination test in EZ: ")
842  }
843  /// ---> A8 (gcdfound)
844  }
845  delta--;
846  }
847  /// ---> A9
848  DEBDECLEVEL( cerr, "ezgcd" );
849  cand *= bCommonDen (cand);
850  if (!isRat)
851  Off (SW_RATIONAL);
852  cand /= icontent (cand);
853  return N (d*cand);
854 }
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm FACTORY_PUBLIC content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:603
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
int getNumVars(const CanonicalForm &f)
int getNumVars ( const CanonicalForm & f )
Definition: cf_ops.cc:314
CanonicalForm lc(const CanonicalForm &f)
CanonicalForm FACTORY_PUBLIC icontent(const CanonicalForm &f)
CanonicalForm icontent ( const CanonicalForm & f )
Definition: cf_gcd.cc:74
int degree(const CanonicalForm &f)
CanonicalForm LC(const CanonicalForm &f)
if(both_non_zero==0)
Definition: cfEzgcd.cc:91
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int l
Definition: cfEzgcd.cc:100
static CanonicalForm gcd_mon(CanonicalForm F, CanonicalForm G)
Definition: cfEzgcd.cc:470
const CanonicalForm & G
Definition: cfEzgcd.cc:55
static int Hensel(const CanonicalForm &UU, CFArray &G, const Evaluation &AA, const CFArray &LeadCoeffs)
Definition: cfEzgcd.cc:302
const CanonicalForm CFMap & M
Definition: cfEzgcd.cc:55
static bool findeval(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Fb, CanonicalForm &Gb, CanonicalForm &Db, REvaluation &b, int delta, int degF, int degG, int maxeval, int &count, int &k, int bound, int &l)
Definition: cfEzgcd.cc:386
bool gcd_test_one(const CanonicalForm &f, const CanonicalForm &g, bool swap, int &d)
Coprimality Check. f and g are assumed to have the same level. If swap is true, the main variables of...
Definition: cfGcdUtil.cc:25
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition: cfModGcd.cc:69
Variable x
Definition: cfModGcd.cc:4082
g
Definition: cfModGcd.cc:4090
int maxNumVars
Definition: cfModGcd.cc:4130
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 )
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:31
static const int SW_USE_EZGCD
set to 1 to use EZGCD over Z
Definition: cf_defs.h:35
FILE * f
Definition: checklibs.c:9
class CFMap
Definition: cf_map.h:85
factory's main class
Definition: canonicalform.h:86
CanonicalForm genOne() const
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
bool inCoeffDomain() const
int level() const
level() returns the level of CO.
bool isUnivariate() const
generate random integers
Definition: cf_random.h:56
factory's class for variables
Definition: factory.h:127
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int const CFList & evaluation
Definition: facAbsFact.cc:52
TIMING_END_AND_PRINT(fac_alg_resultant, "time to compute resultant0: ")
TIMING_START(fac_alg_resultant)
b *CanonicalForm B
Definition: facBivar.cc:52
bool found
Definition: facFactorize.cc:55
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
#define D(A)
Definition: gentable.cc:131
bool delta(X x, Y y, D d)
Definition: TestSuite.h:160
int status int void size_t count
Definition: si_signals.h:59
int status int void * buf
Definition: si_signals.h:59
int gcd(int a, int b)
Definition: walkSupport.cc:836

◆ EZGCD_P()

CanonicalForm EZGCD_P ( const CanonicalForm FF,
const CanonicalForm GG 
)

Extended Zassenhaus GCD for finite fields. In case things become too dense we switch to a modular algorithm.

Definition at line 876 of file cfEzgcd.cc.

877 {
878  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
879  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
880  else if (FF.isZero() && GG.isZero()) return FF.genOne();
881  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
882  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
883  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
884  if (FF == GG) return FF/Lc(FF);
885 
886  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
887  Variable a, oldA;
888  int sizeF= size (FF);
889  int sizeG= size (GG);
890 
891  if (sizeF==1)
892  {
893  return gcd_mon( FF, GG );
894  }
895  else if (sizeG==1)
896  {
897  return gcd_mon( GG, FF );
898  }
899 
900  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
901  {
902  if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
903  return modGCDFq (FF, GG, a);
905  return modGCDGF (FF, GG);
906  else
907  return modGCDFp (FF, GG);
908  }
909 
910  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
911  lcD;
912  CFArray DD( 1, 2 ), lcDD( 1, 2 );
913  int degF, degG, delta, count;
914  int maxeval;
915  int ch=getCharacteristic();
916  maxeval= tmin((ch/
917  (int)(ilog2(ch)*log2exp))*2, maxNumEval);
918  count= 0; // number of eval. used
919  REvaluation b, bt;
920  int gcdfound = 0;
921  Variable x = Variable(1);
922 
923  F= FF;
924  G= GG;
925 
926  CFMap M,N;
927  int smallestDegLev;
928  TIMING_DEFINE(ez_p_compress);
929  TIMING_START (ez_p_compress);
930  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
931 
932  if (best_level == 0) return G.genOne();
933 
934  F= M (F);
935  G= M (G);
936  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
937 
938  TIMING_DEFINE (ez_p_content)
939  TIMING_START (ez_p_content)
940  f = content( F, x ); g = content( G, x );
941  d = gcd( f, g );
942  F /= f; G /= g;
943  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
944 
945  if( F.isUnivariate() && G.isUnivariate() )
946  {
947  if( F.mvar() == G.mvar() )
948  d *= gcd( F, G );
949  else
950  return N (d);
951  return N (d);
952  }
953  if ( F.isUnivariate())
954  {
955  g= content (G,G.mvar());
956  return N(d*gcd(F,g));
957  }
958  if ( G.isUnivariate())
959  {
960  f= content (F,F.mvar());
961  return N(d*gcd(G,f));
962  }
963 
965  sizeF= size (F);
966  sizeG= size (G);
967 
968  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
969  {
970  if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
971  return N (d*modGCDFq (F, G, a));
973  return N (d*modGCDGF (F, G));
974  else
975  return N (d*modGCDFp (F, G));
976  }
977 
978  int dummy= 0;
979  if( gcd_test_one( F, G, false, dummy ) )
980  {
981  return N (d);
982  }
983 
984  bool passToGF= false;
985  bool extOfExt= false;
986  int p= getCharacteristic();
987  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
988  int k= 1;
989  CanonicalForm primElem, imPrimElem;
990  CFList source, dest;
991  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
992  {
993  if (p == 2)
994  setCharacteristic (2, 12, 'Z');
995  else if (p == 3)
996  setCharacteristic (3, 4, 'Z');
997  else if (p == 5 || p == 7)
998  setCharacteristic (p, 3, 'Z');
999  else
1000  setCharacteristic (p, 2, 'Z');
1001  passToGF= true;
1002  F= F.mapinto();
1003  G= G.mapinto();
1004  maxeval= 2*ipower (p, getGFDegree());
1005  }
1006  else if (CFFactory::gettype() == GaloisFieldDomain &&
1007  ipower (p , getGFDegree()) < 50)
1008  {
1009  k= getGFDegree();
1010  if (ipower (p, 2*k) > 50)
1011  setCharacteristic (p, 2*k, gf_name);
1012  else
1013  setCharacteristic (p, 3*k, gf_name);
1014  F= GFMapUp (F, k);
1015  G= GFMapUp (G, k);
1016  maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
1017  }
1018  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
1019  {
1020  int d= degree (getMipo (a));
1021  oldA= a;
1022  Variable v2;
1023  if (p == 2 && d < 6)
1024  {
1025  bool primFail= false;
1026  Variable vBuf;
1027  primElem= primitiveElement (a, vBuf, primFail);
1028  ASSERT (!primFail, "failure in integer factorizer");
1029  if (d < 3)
1030  {
1031  #ifdef HAVE_FLINT
1032  nmod_poly_t Irredpoly;
1033  nmod_poly_init(Irredpoly,p);
1034  nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 3*d+1);
1035  CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
1036  nmod_poly_clear(Irredpoly);
1037  #elif defined(HAVE_NTL)
1038  if (fac_NTL_char != p)
1039  {
1040  fac_NTL_char= p;
1041  zz_p::init (p);
1042  }
1043  zz_pX NTLIrredpoly;
1044  BuildIrred (NTLIrredpoly, d*3);
1045  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1046  #else
1047  factoryError("NTL/FLINT missing: EZGCD_P");
1048  #endif
1049  v2= rootOf (newMipo);
1050  }
1051  else
1052  {
1053  #ifdef HAVE_FLINT
1054  nmod_poly_t Irredpoly;
1055  nmod_poly_init(Irredpoly,p);
1056  nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 2*d+1);
1057  CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
1058  nmod_poly_clear(Irredpoly);
1059  #elif defined(HAVE_NTL)
1060  if (fac_NTL_char != p)
1061  {
1062  fac_NTL_char= p;
1063  zz_p::init (p);
1064  }
1065  zz_pX NTLIrredpoly;
1066  BuildIrred (NTLIrredpoly, d*2);
1067  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1068  #else
1069  factoryError("NTL/FLINT missing: EZGCD_P");
1070  #endif
1071  v2= rootOf (newMipo);
1072  }
1073  imPrimElem= mapPrimElem (primElem, a, v2);
1074  extOfExt= true;
1075  }
1076  else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
1077  {
1078  bool primFail= false;
1079  Variable vBuf;
1080  primElem= primitiveElement (a, vBuf, primFail);
1081  ASSERT (!primFail, "failure in integer factorizer");
1082  #ifdef HAVE_FLINT
1083  nmod_poly_t Irredpoly;
1084  nmod_poly_init(Irredpoly,p);
1085  nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 2*d+1);
1086  CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
1087  nmod_poly_clear(Irredpoly);
1088  #elif defined(HAVE_NTL)
1089  if (fac_NTL_char != p)
1090  {
1091  fac_NTL_char= p;
1092  zz_p::init (p);
1093  }
1094  zz_pX NTLIrredpoly;
1095  BuildIrred (NTLIrredpoly, d*2);
1096  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1097  #else
1098  factoryError("NTL/FLINT missing: EZGCD_P");
1099  #endif
1100  v2= rootOf (newMipo);
1101  imPrimElem= mapPrimElem (primElem, a, v2);
1102  extOfExt= true;
1103  }
1104  if (extOfExt)
1105  {
1106  maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
1107  F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
1108  G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
1109  a= v2;
1110  }
1111  }
1112 
1113  lcF = LC( F, x ); lcG = LC( G, x );
1114  lcD = gcd( lcF, lcG );
1115 
1116  delta = 0;
1117  degF = degree( F, x ); degG = degree( G, x );
1118 
1119  if (algExtension)
1120  b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
1121  else
1122  { // both not in extension given by algebraic variable
1124  b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
1125  else
1126  b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
1127  }
1128 
1129  CanonicalForm cand, contcand;
1131  int o, t;
1132  o= 0;
1133  t= 1;
1134  int goodPointCount= 0;
1135  TIMING_DEFINE(ez_p_eval);
1136  while( !gcdfound )
1137  {
1138  TIMING_START (ez_p_eval);
1139  if( !findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
1140  maxeval/maxNumVars, t ))
1141  { // too many eval. used --> try another method
1142  Off (SW_USE_EZGCD_P);
1143  result= gcd (F,G);
1144  On (SW_USE_EZGCD_P);
1145  if (passToGF)
1146  {
1148  setCharacteristic (p);
1151  prune (alpha);
1152  }
1153  if (k > 1)
1154  {
1155  result= GFMapDown (result, k);
1157  }
1158  if (extOfExt)
1159  {
1160  result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1161  prune1 (oldA);
1162  }
1163  return N (d*result);
1164  }
1165  TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
1166  delta = degree( Db );
1167  if (delta == degF)
1168  {
1169  if (degF <= degG && fdivides (F, G))
1170  {
1171  if (passToGF)
1172  {
1174  setCharacteristic (p);
1176  F= GF2FalphaRep (F, alpha);
1177  prune (alpha);
1178  }
1179  if (k > 1)
1180  {
1181  F= GFMapDown (F, k);
1183  }
1184  if (extOfExt)
1185  {
1186  F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1187  prune1 (oldA);
1188  }
1189  return N (d*F);
1190  }
1191  else
1192  delta--;
1193  }
1194  else if (delta == degG)
1195  {
1196  if (degG <= degF && fdivides (G, F))
1197  {
1198  if (passToGF)
1199  {
1201  setCharacteristic (p);
1203  G= GF2FalphaRep (G, alpha);
1204  prune (alpha);
1205  }
1206  if (k > 1)
1207  {
1208  G= GFMapDown (G, k);
1210  }
1211  if (extOfExt)
1212  {
1213  G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1214  prune1 (oldA);
1215  }
1216  return N (d*G);
1217  }
1218  else
1219  delta--;
1220  }
1221  if( delta == 0 )
1222  {
1223  if (passToGF)
1224  setCharacteristic (p);
1225  if (k > 1)
1227  return N (d);
1228  }
1229  while( true )
1230  {
1231  bt = b;
1232  TIMING_START (ez_p_eval);
1233  if( !findeval(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
1234  maxeval/maxNumVars, t ))
1235  { // too many eval. used --> try another method
1236  Off (SW_USE_EZGCD_P);
1237  result= gcd (F,G);
1238  On (SW_USE_EZGCD_P);
1239  if (passToGF)
1240  {
1242  setCharacteristic (p);
1245  prune (alpha);
1246  }
1247  if (k > 1)
1248  {
1249  result= GFMapDown (result, k);
1251  }
1252  if (extOfExt)
1253  {
1254  result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1255  prune1 (oldA);
1256  }
1257  return N (d*result);
1258  }
1259  TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
1260  int dd = degree( Dbt );
1261  if( dd == 0 )
1262  {
1263  if (passToGF)
1264  setCharacteristic (p);
1265  if (k > 1)
1267  return N (d);
1268  }
1269  if( dd == delta )
1270  {
1271  goodPointCount++;
1272  if (goodPointCount == 5)
1273  break;
1274  }
1275  if( dd < delta )
1276  {
1277  goodPointCount= 0;
1278  delta = dd;
1279  b = bt;
1280  Db = Dbt; Fb = Fbt; Gb = Gbt;
1281  }
1282  if (delta == degF)
1283  {
1284  if (degF <= degG && fdivides (F, G))
1285  {
1286  if (passToGF)
1287  {
1289  setCharacteristic (p);
1291  F= GF2FalphaRep (F, alpha);
1292  prune (alpha);
1293  }
1294  if (k > 1)
1295  {
1296  F= GFMapDown (F, k);
1298  }
1299  if (extOfExt)
1300  {
1301  F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1302  prune1 (oldA);
1303  }
1304  return N (d*F);
1305  }
1306  else
1307  delta--;
1308  }
1309  else if (delta == degG)
1310  {
1311  if (degG <= degF && fdivides (G, F))
1312  {
1313  if (passToGF)
1314  {
1316  setCharacteristic (p);
1318  G= GF2FalphaRep (G, alpha);
1319  prune (alpha);
1320  }
1321  if (k > 1)
1322  {
1323  G= GFMapDown (G, k);
1325  }
1326  if (extOfExt)
1327  {
1328  G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1329  prune1 (oldA);
1330  }
1331  return N (d*G);
1332  }
1333  else
1334  delta--;
1335  }
1336  if( delta == 0 )
1337  {
1338  if (passToGF)
1339  setCharacteristic (p);
1340  if (k > 1)
1342  return N (d);
1343  }
1344  }
1345  if( delta != degF && delta != degG )
1346  {
1347  bool B_is_F;
1348  CanonicalForm xxx1, xxx2;
1350  DD[1] = Fb / Db;
1351  buf= Gb/Db;
1352  xxx1 = gcd( DD[1], Db );
1353  xxx2 = gcd( buf, Db );
1354  if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1355  (size (F) <= size (G)))
1356  || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
1357  {
1358  B = F;
1359  DD[2] = Db;
1360  lcDD[1] = lcF;
1361  lcDD[2] = lcD;
1362  B_is_F = true;
1363  }
1364  else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1365  (size (G) < size (F)))
1366  || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
1367  {
1368  DD[1] = buf;
1369  B = G;
1370  DD[2] = Db;
1371  lcDD[1] = lcG;
1372  lcDD[2] = lcD;
1373  B_is_F = false;
1374  }
1375  else // special case handling
1376  {
1377  Off (SW_USE_EZGCD_P);
1378  result= gcd (F,G);
1379  On (SW_USE_EZGCD_P);
1380  if (passToGF)
1381  {
1383  setCharacteristic (p);
1386  prune (alpha);
1387  }
1388  if (k > 1)
1389  {
1390  result= GFMapDown (result, k);
1392  }
1393  if (extOfExt)
1394  {
1395  result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1396  prune1 (oldA);
1397  }
1398  return N (d*result);
1399  }
1400  DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
1401  DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
1402 
1403  if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
1404  {
1405  if (algExtension)
1406  {
1407  result= modGCDFq (F, G, a);
1408  if (extOfExt)
1409  {
1410  result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1411  prune1 (oldA);
1412  }
1413  return N (d*result);
1414  }
1416  {
1417  result= modGCDGF (F, G);
1418  if (passToGF)
1419  {
1421  setCharacteristic (p);
1424  prune (alpha);
1425  }
1426  if (k > 1)
1427  {
1428  result= GFMapDown (result, k);
1430  }
1431  return N (d*result);
1432  }
1433  else
1434  return N (d*modGCDFp (F,G));
1435  }
1436 
1437  TIMING_DEFINE(ez_p_hensel_lift);
1438  TIMING_START (ez_p_hensel_lift);
1439  gcdfound= Hensel (B*lcD, DD, b, lcDD);
1440  TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
1441 
1442  if (gcdfound == -1) //things became dense
1443  {
1444  if (algExtension)
1445  {
1446  result= modGCDFq (F, G, a);
1447  if (extOfExt)
1448  {
1449  result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1450  prune1 (oldA);
1451  }
1452  return N (d*result);
1453  }
1455  {
1456  result= modGCDGF (F, G);
1457  if (passToGF)
1458  {
1460  setCharacteristic (p);
1463  prune (alpha);
1464  }
1465  if (k > 1)
1466  {
1467  result= GFMapDown (result, k);
1469  }
1470  return N (d*result);
1471  }
1472  else
1473  {
1474  if (p >= cf_getBigPrime(0))
1475  return N (d*sparseGCDFp (F,G));
1476  else
1477  return N (d*modGCDFp (F,G));
1478  }
1479  }
1480 
1481  if (gcdfound == 1)
1482  {
1483  TIMING_DEFINE(termination_test);
1484  TIMING_START (termination_test);
1485  contcand= content (DD[2], Variable (1));
1486  cand = DD[2] / contcand;
1487  if (B_is_F)
1488  gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
1489  else
1490  gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
1491  TIMING_END_AND_PRINT (termination_test,
1492  "time for termination test EZ_P: ");
1493 
1494  if (passToGF && gcdfound)
1495  {
1497  setCharacteristic (p);
1500  prune (alpha);
1501  }
1502  if (k > 1 && gcdfound)
1503  {
1504  cand= GFMapDown (cand, k);
1506  }
1507  if (extOfExt && gcdfound)
1508  {
1509  cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
1510  prune1 (oldA);
1511  }
1512  }
1513  }
1514  delta--;
1515  goodPointCount= 0;
1516  }
1517  return N (d*cand);
1518 }
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
Definition: NTLconvert.cc:255
VAR long fac_NTL_char
Definition: NTLconvert.cc:46
int ilog2(const CanonicalForm &a)
int getGFDegree()
Definition: cf_char.cc:75
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)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
STATIC_VAR int maxNumEval
Definition: cfEzgcd.cc:871
static const double log2exp
Definition: cfEzgcd.cc:45
STATIC_VAR int sizePerVars1
Definition: cfEzgcd.cc:872
int k
Definition: cfEzgcd.cc:99
CanonicalForm modGCDFq(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, Variable &alpha, CFList &l, bool &topLevel)
GCD of F and G over , l and topLevel are only used internally, output is monic based on Alg....
Definition: cfModGcd.cc:478
int p
Definition: cfModGcd.cc:4078
CanonicalForm modGCDFp(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:1223
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:3562
CanonicalForm modGCDGF(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, CFList &l, bool &topLevel)
GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for Computer Algebra" by Gedde...
Definition: cfModGcd.cc:872
#define ASSERT(expression, message)
Definition: cf_assert.h:99
static const int SW_USE_EZGCD_P
set to 1 to use EZGCD over F_q
Definition: cf_defs.h:37
#define GaloisFieldDomain
Definition: cf_defs.h:18
CanonicalForm mapPrimElem(const CanonicalForm &primElem, const Variable &alpha, const Variable &beta)
compute the image of a primitive element of in . We assume .
Definition: cf_map_ext.cc:450
CanonicalForm GFMapDown(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition: cf_map_ext.cc:276
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
Definition: cf_map_ext.cc:342
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 ,...
Definition: cf_map_ext.cc:123
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
Definition: cf_map_ext.cc:70
CanonicalForm GFMapUp(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition: cf_map_ext.cc:240
CanonicalForm GF2FalphaRep(const CanonicalForm &F, const Variable &alpha)
changes representation by primitive element to representation by residue classes modulo a Conway poly...
Definition: cf_map_ext.cc:195
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
GLOBAL_VAR flint_rand_t FLINTrandom
Definition: cf_random.cc:25
VAR void(* factoryError)(const char *s)
Definition: cf_util.cc:80
int ipower(int b, int m)
int ipower ( int b, int m )
Definition: cf_util.cc:27
generate random elements in F_p(alpha)
Definition: cf_random.h:70
static int gettype()
Definition: cf_factory.h:28
CF_NO_INLINE bool isZero() const
bool inBaseDomain() const
CanonicalForm mapinto() const
generate random elements in F_p
Definition: cf_random.h:44
generate random elements in GF
Definition: cf_random.h:32
Variable alpha
Definition: facAbsBiFact.cc:51
CanonicalForm mipo
Definition: facAlgExt.cc:57
nmod_poly_init(FLINTmipo, getCharacteristic())
nmod_poly_clear(FLINTmipo)
Variable FACTORY_PUBLIC rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables
Definition: variable.cc:162
void prune1(const Variable &alpha)
Definition: variable.cc:291
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
void FACTORY_PUBLIC prune(Variable &alpha)
Definition: variable.cc:261
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
VAR char gf_name
Definition: gfops.cc:52
INST_VAR CanonicalForm gf_mipo
Definition: gfops.cc:56
void init()
Definition: lintree.cc:864
#define TIMING_DEFINE(t)
Definition: timing.h:91

◆ findeval()

static bool findeval ( const CanonicalForm F,
const CanonicalForm G,
CanonicalForm Fb,
CanonicalForm Gb,
CanonicalForm Db,
REvaluation b,
int  delta,
int  degF,
int  degG,
int  maxeval,
int &  count,
int &  k,
int  bound,
int &  l 
)
static

Definition at line 386 of file cfEzgcd.cc.

390 {
391  if( count == 0 && delta != 0)
392  {
393  if( count++ > maxeval )
394  return false;
395  }
396  if (count > 0)
397  {
398  b.nextpoint(k);
399  if (k == 0)
400  k++;
401  l++;
402  if (l > bound)
403  {
404  l= 1;
405  k++;
406  if (k > tmax (F.level(), G.level()) - 1)
407  return false;
408  b.nextpoint (k);
409  }
410  if (count++ > maxeval)
411  return false;
412  }
413  while( true )
414  {
415  Fb = b( F );
416  if( degree( Fb, 1 ) == degF )
417  {
418  Gb = b( G );
419  if( degree( Gb, 1 ) == degG )
420  {
421  Db = gcd( Fb, Gb );
422  if( delta > 0 )
423  {
424  if( degree( Db, 1 ) <= delta )
425  return true;
426  }
427  else
428  {
429  k++;
430  return true;
431  }
432  }
433  }
434  if (k == 0)
435  k++;
436  b.nextpoint(k);
437  l++;
438  if (l > bound)
439  {
440  l= 1;
441  k++;
442  if (k > tmax (F.level(), G.level()) - 1)
443  return false;
444  b.nextpoint (k);
445  }
446  if( count++ > maxeval )
447  return false;
448  }
449 }
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460

◆ 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  {
76  both_non_zero++;
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

◆ gcd_mon()

static CanonicalForm gcd_mon ( CanonicalForm  F,
CanonicalForm  G 
)
static

Definition at line 470 of file cfEzgcd.cc.

471 {
472  // assume: size(F)==1
473  CanonicalForm cf=F;
474  int ll=tmax(F.level(),G.level());
475  int *exp=NEW_ARRAY(int,ll+1);
476  for(int i=ll;i>=0;i--) exp[i]=0;
477  CanonicalForm c=F;
478  while(!c.inCoeffDomain())
479  {
480  exp[c.level()]=c.degree();
481  c=c.LC();
482  cf=c;
483  }
484  gcd_mon_rec(G,cf,exp,G.level()+1);
486  for(int i=0;i<=ll;i++)
487  {
488  if (exp[i]>0) res*=power(Variable(i),exp[i]);
489  }
490  DELETE_ARRAY(exp);
491  return res;
492 }
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
static void gcd_mon_rec(CanonicalForm G, CanonicalForm &cf, int *exp, int pl)
Definition: cfEzgcd.cc:452
DELETE_ARRAY(degsf)
CanonicalForm cf
Definition: cfModGcd.cc:4083
#define NEW_ARRAY(T, N)
Definition: cf_defs.h:64
CanonicalForm LC() const
int degree() const
Returns -1 for the zero polynomial and 0 if CO is in a base domain.
CanonicalForm res
Definition: facAbsFact.cc:60
gmp_float exp(const gmp_float &a)
Definition: mpr_complex.cc:357

◆ gcd_mon_rec()

static void gcd_mon_rec ( CanonicalForm  G,
CanonicalForm cf,
int *  exp,
int  pl 
)
static

Definition at line 452 of file cfEzgcd.cc.

453 { // prevoius level: pl
454  if (G.inCoeffDomain())
455  {
456  for(int i=pl-1;i>0;i--) exp[i]=0;
457  cf=gcd(cf,G);
458  return;
459  }
460  int l=G.level();
461  for(int i=pl-1;i>l;i--) exp[i]=0;
462  for(CFIterator i=G; i.hasTerms(); i++)
463  {
464  if (i.exp()<exp[l]) exp[l]=i.exp();
465  gcd_mon_rec(i.coeff(),cf,exp,l);
466  }
467 }
class to iterate through CanonicalForm's
Definition: cf_iter.h:44

◆ Hensel()

static int Hensel ( const CanonicalForm UU,
CFArray G,
const Evaluation AA,
const CFArray LeadCoeffs 
)
inlinestatic

Definition at line 302 of file cfEzgcd.cc.

304 {
305  CFList factors;
306  factors.append (G[1]);
307  factors.append (G[2]);
308 
309  CFMap NN, MM;
310  Evaluation A= optimize4Lift (UU, MM, NN, AA);
311 
312  CanonicalForm U= MM (UU);
313  CFArray LCs= CFArray (1,2);
314  LCs [1]= MM (LeadCoeffs [1]);
315  LCs [2]= MM (LeadCoeffs [2]);
316 
318  long termEstimate= size (U);
319  int ch=getCharacteristic();
320  for (int i= A.min(); i <= A.max(); i++)
321  {
322  if (!A[i].isZero() &&
323  ((ch > degree (U,i)) || ch == 0))
324  {
325  termEstimate *= degree (U,i)*2;
326  termEstimate /= 3;
327  }
328  evaluation.append (A [i]);
329  }
330  if (termEstimate/getNumVars(U) > 500)
331  return -1;
332  CFList UEval;
333  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
334 
335  if (size (shiftedU)/getNumVars (U) > 500)
336  return -1;
337 
338  CFArray shiftedLCs= CFArray (2);
339  CFList shiftedLCsEval1, shiftedLCsEval2;
340  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
341  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
342  factors.insert (1);
343  int liftBound= degree (UEval.getLast(), 2) + 1;
344  CFArray Pi;
345  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
346  CFList diophant;
347  CFArray lcs= CFArray (2);
348  lcs [0]= shiftedLCsEval1.getFirst();
349  lcs [1]= shiftedLCsEval2.getFirst();
350  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
351  lcs, false);
352 
353  for (CFListIterator i= factors; i.hasItem(); i++)
354  {
355  if (!fdivides (i.getItem(), UEval.getFirst()))
356  return 0;
357  }
358 
359  int * liftBounds;
360  bool noOneToOne= false;
361  if (U.level() > 2)
362  {
363  liftBounds= NEW_ARRAY(int,U.level() - 1); /* index: 0.. U.level()-2 */
364  liftBounds[0]= liftBound;
365  for (int i= 1; i < U.level() - 1; i++)
366  liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
367  factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
368  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
369  diophant, noOneToOne);
370  DELETE_ARRAY(liftBounds);
371  if (noOneToOne)
372  return 0;
373  }
374  G[1]= factors.getFirst();
375  G[2]= factors.getLast();
376  G[1]= myReverseShift (G[1], evaluation);
377  G[2]= myReverseShift (G[2], evaluation);
378  G[1]= NN (G[1]);
379  G[2]= NN (G[2]);
380  return 1;
381 }
Array< CanonicalForm > CFArray
Matrix< CanonicalForm > CFMatrix
static CanonicalForm myShift2Zero(const CanonicalForm &F, CFList &Feval, const CFList &evaluation)
Definition: cfEzgcd.cc:195
static Evaluation optimize4Lift(const CanonicalForm &F, CFMap &M, CFMap &N, const Evaluation &A)
Definition: cfEzgcd.cc:231
static CanonicalForm myReverseShift(const CanonicalForm &F, const CFList &evaluation)
Definition: cfEzgcd.cc:215
class to evaluate a polynomial at points
Definition: cf_eval.h:32
T getFirst() const
Definition: ftmpl_list.cc:279
int length() const
Definition: ftmpl_list.cc:273
void append(const T &)
Definition: ftmpl_list.cc:256
T getLast() const
Definition: ftmpl_list.cc:309
void insert(const T &)
Definition: ftmpl_list.cc:193
CFList nonMonicHenselLift2(const CFList &F, const CFList &factors, const CFList &MOD, CFList &diophant, CFArray &Pi, CFMatrix &M, int lOld, int &lNew, const CFList &LCs1, const CFList &LCs2, bool &bad)
Definition: facHensel.cc:2632
void nonMonicHenselLift12(const CanonicalForm &F, CFList &factors, int l, CFArray &Pi, CFList &diophant, CFMatrix &M, const CFArray &LCs, bool sort)
Hensel lifting from univariate to bivariate, factors need not to be monic.
Definition: facHensel.cc:2154
bool isZero(const CFArray &A)
checks if entries of A are zero
#define A
Definition: sirandom.c:24

◆ if()

if ( both_non_zero  = = 0)

Definition at line 91 of file cfEzgcd.cc.

92  {
95  return 0;
96  }

◆ myReverseShift()

static CanonicalForm myReverseShift ( const CanonicalForm F,
const CFList evaluation 
)
inlinestatic

Definition at line 215 of file cfEzgcd.cc.

216 {
217  int l= evaluation.length() + 1;
220  int Flevel=F.level();
221  for (int i= 2; i < l + 1; i++, j++)
222  {
223  if (Flevel < i)
224  continue;
225  result= result (Variable (i) - j.getItem(), i);
226  }
227  return result;
228 }
int Flevel
Definition: cfEzgcd.cc:101
int j
Definition: facHensel.cc:110

◆ myShift2Zero()

static CanonicalForm myShift2Zero ( const CanonicalForm F,
CFList Feval,
const CFList evaluation 
)
inlinestatic

Definition at line 195 of file cfEzgcd.cc.

197 {
198  CanonicalForm A= F;
199  int k= 2;
200  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
201  A= A (Variable (k) + i.getItem(), k);
202 
204  Feval= CFList();
205  Feval.append (buf);
206  for (k= evaluation.length() + 1; k > 2; k--)
207  {
208  buf= mod (buf, Variable (k));
209  Feval.insert (buf);
210  }
211  return A;
212 }
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
List< CanonicalForm > CFList
CanonicalForm Feval
Definition: facAbsFact.cc:60

◆ optimize4Lift()

static Evaluation optimize4Lift ( const CanonicalForm F,
CFMap M,
CFMap N,
const Evaluation A 
)
inlinestatic

Definition at line 231 of file cfEzgcd.cc.

233 {
234  int n= F.level();
235  int * degsf= NEW_ARRAY(int,n + 1);
236 
237  for (int i = n; i >= 0; i--)
238  degsf[i]= 0;
239 
240  degsf= degrees (F, degsf);
241 
242  Evaluation result= Evaluation (A.min(), A.max());
243  ASSERT (A.min() == 2, "expected A.min() == 2");
244  int max_deg;
245  int k= n;
246  int l= 1;
247  int i= 2;
248  int pos= 2;
249  while (k > 1)
250  {
251  max_deg= degsf [i]; // i is always 2 here, n>=2
252  while ((i<n) &&(max_deg == 0))
253  {
254  i++;
255  max_deg= degsf [i];
256  }
257  l= i;
258  for (int j= i + 1; j <= n; j++)
259  {
260  if (degsf[j] > max_deg)
261  {
262  max_deg= degsf[j];
263  l= j;
264  }
265  }
266 
267  if (l <= n)
268  {
269  if (l != pos)
270  {
271  result.setValue (pos, A [l]);
272  M.newpair (Variable (l), Variable (pos));
273  N.newpair (Variable (pos), Variable (l));
274  degsf[l]= 0;
275  l= 2;
276  if (k == 2 && n == 3)
277  {
278  result.setValue (l, A [pos]);
279  M.newpair (Variable (pos), Variable (l));
280  N.newpair (Variable (l), Variable (pos));
281  degsf[pos]= 0;
282  }
283  }
284  else
285  {
286  result.setValue (l, A [l]);
287  degsf [l]= 0;
288  }
289  }
290  pos++;
291  k--;
292  l= 2;
293  }
294 
296 
297  return result;
298 }
int * degrees(const CanonicalForm &f, int *degs=0)
int * degrees ( const CanonicalForm & f, int * degs )
Definition: cf_ops.cc:493

◆ TIMING_DEFINE_PRINT()

TIMING_DEFINE_PRINT ( ez_eval  ) const &

◆ while()

while ( k  ,
 
)

Definition at line 133 of file cfEzgcd.cc.

134  {
135  max_min_deg= tmin (degsf[i], degsg[i]);
136  while (max_min_deg == 0)
137  {
138  i++;
139  max_min_deg= tmin (degsf[i], degsg[i]);
140  }
141  for (int j= i + 1; j <= m; j++)
142  {
143  if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
144  (tmin (degsf[j], degsg[j]) != 0))
145  {
146  max_min_deg= tmin (degsf[j], degsg[j]);
147  l= j;
148  }
149  }
150 
151  if (l != 0)
152  {
153  if (l != k)
154  {
155  M.newpair (Variable (l), Variable(k));
156  N.newpair (Variable (k), Variable(l));
157  degsf[l]= 0;
158  degsg[l]= 0;
159  l= 0;
160  }
161  else
162  {
163  degsf[l]= 0;
164  degsg[l]= 0;
165  l= 0;
166  }
167  }
168  else if (l == 0)
169  {
170  if (i != k)
171  {
172  M.newpair (Variable (i), Variable (k));
173  N.newpair (Variable (k), Variable (i));
174  degsf[i]= 0;
175  degsg[i]= 0;
176  }
177  else
178  {
179  degsf[i]= 0;
180  degsg[i]= 0;
181  }
182  i++;
183  }
184  k--;
185  }
int m
Definition: cfEzgcd.cc:128
int max_min_deg
Definition: cfEzgcd.cc:129

Variable Documentation

◆ both_non_zero

return both_non_zero
Initial value:
{
int n= tmax (F.level(), G.level())

Definition at line 56 of file cfEzgcd.cc.

◆ degsf

degsf = NEW_ARRAY(int,n + 1)

Definition at line 59 of file cfEzgcd.cc.

◆ degsg

degsg = NEW_ARRAY(int,n + 1)

Definition at line 60 of file cfEzgcd.cc.

◆ f_zero

int f_zero = 0

Definition at line 69 of file cfEzgcd.cc.

◆ Flevel

int Flevel =F.level()

Definition at line 101 of file cfEzgcd.cc.

◆ G

Definition at line 55 of file cfEzgcd.cc.

◆ g_zero

int g_zero = 0

Definition at line 70 of file cfEzgcd.cc.

◆ Glevel

int Glevel =G.level()

Definition at line 102 of file cfEzgcd.cc.

◆ i

int i = 1

Definition at line 132 of file cfEzgcd.cc.

◆ k

k = 1

Definition at line 99 of file cfEzgcd.cc.

◆ l

l = 1

Definition at line 100 of file cfEzgcd.cc.

◆ log2exp

const double log2exp = 1.442695041
static

Definition at line 45 of file cfEzgcd.cc.

◆ M

Definition at line 55 of file cfEzgcd.cc.

◆ m

int m = tmin (Flevel, Glevel)

Definition at line 128 of file cfEzgcd.cc.

◆ max_min_deg

int max_min_deg

Definition at line 129 of file cfEzgcd.cc.

◆ maxNumEval

STATIC_VAR int maxNumEval = 200

Definition at line 871 of file cfEzgcd.cc.

◆ N

Definition at line 56 of file cfEzgcd.cc.

◆ sizePerVars1

STATIC_VAR int sizePerVars1 = 500

Definition at line 872 of file cfEzgcd.cc.