Changeset 9ff686 in git for factory/cf_gcd_smallp.cc


Ignore:
Timestamp:
Jun 14, 2011, 11:29:25 AM (13 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
d3b70ae172ff67aecbfd2f47077e21aa319554fd
Parents:
3ffd6d2e233ebf7267b3b68dc7d6d7d506305516
Message:
changes to EZGCD_P
changes to passing to extensions


git-svn-id: file:///usr/local/Singular/svn/trunk@14275 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_gcd_smallp.cc

    r3ffd6d r9ff686  
    5252/// compressing two polynomials F and G, M is used for compressing,
    5353/// N to reverse the compression
    54 static inline
    5554int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
    5655                CFMap & N, bool topLevel)
    … …  
    126125    // sort Variables x_{i} in increasing order of
    127126    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
    128     int m= tmin (F.level(), G.level());
    129     int max_min_deg;
     127    int m= tmax (F.level(), G.level());
     128    int min_max_deg;
    130129    k= both_non_zero;
    131130    l= 0;
    … …  
    133132    while (k > 0)
    134133    {
    135       max_min_deg= tmin (degsf[i], degsg[i]);
    136       while (max_min_deg == 0)
     134      if (degsf [i] != 0 && degsg [i] != 0)
     135        min_max_deg= tmax (degsf[i], degsg[i]);
     136      else
     137        min_max_deg= 0;
     138      while (min_max_deg == 0)
    137139      {
    138140        i++;
    139         max_min_deg= tmin (degsf[i], degsg[i]);
     141        if (degsf [i] != 0 && degsg [i] != 0)
     142          min_max_deg= tmax (degsf[i], degsg[i]);
     143        else
     144          min_max_deg= 0;
    140145      }
    141146      for (int j= i + 1; j <=  m; j++)
    142147      {
    143         if (tmin (degsf[j],degsg[j]) >= max_min_deg)
    144         {
    145           max_min_deg= tmin (degsf[j], degsg[j]);
     148        if (degsf[j] != 0 && degsg [j] != 0 &&
     149            tmax (degsf[j],degsg[j]) <= min_max_deg)
     150        {
     151          min_max_deg= tmax (degsf[j], degsg[j]);
    146152          l= j;
    147153        }
    … …  
    209215  return 1;
    210216}
    211 
    212217
    213218int
    … …  
    219224  if (degree (F, x) <= 1 || degree (G, x) <= 1)
    220225    return 0;
    221   CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive
     226  CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping is pretty expensive
    222227  CanonicalForm g= swapvar (G, G.mvar(), x);
    223228  int sizef= 0;
    … …  
    437442/// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})
    438443static inline CanonicalForm
    439 newtonInterp(const CanonicalForm alpha, const CanonicalForm u, const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly, const Variable & x)
     444newtonInterp(const CanonicalForm alpha, const CanonicalForm u,
     445             const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly,
     446             const Variable & x)
    440447{
    441448  CanonicalForm interPoly;
    442449
    443   interPoly = oldInterPoly + ((u - oldInterPoly(alpha, x)) / newtonPoly(alpha, x)) * newtonPoly;
     450  interPoly= oldInterPoly + ((u - oldInterPoly(alpha, x))/newtonPoly(alpha, x))
     451             *newtonPoly;
    444452  return interPoly;
    445453}
    … …  
    487495  } while (find (list, random));
    488496  return random;
     497}
     498
     499static inline
     500Variable chooseExtension (const Variable & alpha)
     501{
     502  zz_p::init (getCharacteristic());
     503  zz_pX NTLIrredpoly;
     504  int i, m;
     505  // extension of F_p needed
     506  if (alpha.level() == 1)
     507  {
     508    i= 1;
     509    m= 2;
     510  } //extension of F_p(alpha)
     511  if (alpha.level() != 1)
     512  {
     513    i= 4;
     514    m= degree (getMipo (alpha));
     515  }
     516  BuildIrred (NTLIrredpoly, i*m);
     517  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
     518  return rootOf (newMipo);
    489519}
    490520
    … …  
    625655      source= CFList();
    626656      dest= CFList();
    627       int num_vars= tmin (getNumVars(A), getNumVars(B));
    628       num_vars--;
    629 
    630       choose_extension (d, num_vars, V_buf);
     657
     658      Variable V_buf3= V_buf;
     659      V_buf= chooseExtension (V_buf);
    631660      bool prim_fail= false;
    632661      Variable V_buf2;
    633662      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     663
     664      if (V_buf3 != alpha)
     665      {
     666        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     667        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     668        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
     669                             source, dest);
     670        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
     671        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
     672        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
     673                            source, dest);
     674        for (CFListIterator i= l; i.hasItem(); i++)
     675          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
     676                                source, dest);
     677      }
    634678
    635679      ASSERT (!prim_fail, "failure in integer factorizer");
    … …  
    909953    if (fail)
    910954    {
    911       int num_vars= tmin (getNumVars(A), getNumVars(B));
    912       num_vars--;
    913955      p= getCharacteristic();
    914       expon= (int) floor ((log ((double) ipower (d + 1, num_vars))/log ((double) p)));
    915       if (expon < 2)
    916         expon= 2;
     956      expon= 2;
    917957      kk= getGFDegree();
    918958      if (ipower (p, kk*expon) < (1 << 16))
    … …  
    12021242  {
    12031243    if (inextension)
    1204       random_element= randomElement (m, alpha, l, fail);
     1244      random_element= randomElement (m, V_buf, l, fail);
    12051245    else
    12061246      random_element= FpRandomElement (m, l, fail);
    … …  
    12441284      } while (fail);
    12451285      list= CFList();
     1286      V_buf= alpha;
    12461287      TIMING_START (gcd_recursion);
    12471288      G_random_element=
    … …  
    12561297      source= CFList();
    12571298      dest= CFList();
    1258       int num_vars= tmin (getNumVars(A), getNumVars(B));
    1259       num_vars--;
    1260       V_buf= alpha;
    1261       choose_extension (d, num_vars, V_buf);
     1299
     1300      Variable V_buf3= V_buf;
     1301      V_buf= chooseExtension (V_buf);
    12621302      bool prim_fail= false;
    12631303      Variable V_buf2;
    12641304      CanonicalForm prim_elem, im_prim_elem;
    12651305      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     1306
     1307      if (V_buf3 != alpha)
     1308      {
     1309        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     1310        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     1311        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
     1312                             source, dest);
     1313        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
     1314        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
     1315        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
     1316                            dest);
     1317        for (CFListIterator i= l; i.hasItem(); i++)
     1318          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
     1319                                source, dest);
     1320      }
    12661321
    12671322      ASSERT (!prim_fail, "failure in integer factorizer");
    … …  
    12971352                            "time for recursive call: ");
    12981353      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
    1299       alpha= V_buf;
    13001354    }
    13011355
    … …  
    19832037    if (evalFail)
    19842038    {
    1985       if (V_buf != Variable (1))
     2039      if (V_buf.level() != 1)
    19862040      {
    19872041        do
    19882042        {
    1989           int num_vars= tmin (getNumVars(A), getNumVars(B));
    1990           int d= totaldegree (A, Variable(2), Variable (A.level()));
    1991           d= tmin (d, totaldegree (B, Variable(2), Variable (B.level())));
    1992           Variable V_buf2= V_buf;
    1993           choose_extension (d, num_vars, V_buf2);
     2043          Variable V_buf2= chooseExtension (V_buf);
    19942044          source= CFList();
    19952045          dest= CFList();
    … …  
    20202070          }
    20212071
    2022           if (alpha != Variable (1))
     2072          if (alpha.level() != 1)
    20232073          {
    20242074            A= mapUp (A, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
    … …  
    21342184  delete[] pL;
    21352185
    2136   if (alpha != Variable (1) && V_buf != alpha)
     2186  if (alpha.level() != 1 && V_buf != alpha)
    21372187  {
    21382188    CFList u, v;
    … …  
    22572307      if (evalFail)
    22582308      {
    2259         if (V_buf != Variable (1))
     2309        if (V_buf.level() != 1)
    22602310        {
    22612311          do
    22622312          {
    2263             int num_vars= tmin (getNumVars(A), getNumVars(B));
    2264             int d= totaldegree (A, Variable(2), Variable (A.level()));
    2265             d= tmin (d, totaldegree (B, Variable(2), Variable (B.level())));
    2266             Variable V_buf2= V_buf;
    2267             choose_extension (d, num_vars, V_buf2);
     2313            Variable V_buf2= chooseExtension (V_buf);
    22682314            source= CFList();
    22692315            dest= CFList();
    … …  
    24462492  pMat[1]= Mat;
    24472493
    2448   if (V_buf != x)
     2494  if (V_buf.level() != 1)
    24492495    solution= solveSystemFq (pMat[1],
    24502496                             pL[minimalColumnsIndex], V_buf);
    … …  
    25772623    for (int i= 0; i < skelSize; i++)
    25782624    {
    2579       if (V_buf == x)
     2625      if (V_buf.level() == 1)
    25802626        rk= gaussianElimFp (pMat[i], pL[i]);
    25812627      else
    … …  
    26232669    }
    26242670
    2625     if (V_buf != x)
     2671    if (V_buf.level() != 1)
    26262672      solution= solveSystemFq (Mat, bufArray, V_buf);
    26272673    else
    … …  
    26502696        result += Monoms [ind]*bufSolution[i];
    26512697    }
    2652     if (alpha != Variable (1) && V_buf != alpha)
     2698    if (alpha.level() != 1 && V_buf != alpha)
    26532699    {
    26542700      CFList u, v;
    … …  
    27412787  delete[] pM;
    27422788
    2743   if (alpha != Variable (1) && V_buf != alpha)
     2789  if (alpha.level() != 1 && V_buf != alpha)
    27442790  {
    27452791    CFList u, v;
    … …  
    28742920      source= CFList();
    28752921      dest= CFList();
    2876       int num_vars= tmin (getNumVars(A), getNumVars(B));
    2877       num_vars--;
    2878 
    2879       choose_extension (d, num_vars, V_buf);
     2922
     2923      Variable V_buf3= V_buf;
     2924      V_buf= chooseExtension (V_buf);
    28802925      bool prim_fail= false;
    28812926      Variable V_buf2;
    28822927      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     2928
     2929      if (V_buf3 != alpha)
     2930      {
     2931        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     2932        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     2933        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
     2934                             source, dest);
     2935        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
     2936        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
     2937        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
     2938                            dest);
     2939        for (CFListIterator i= l; i.hasItem(); i++)
     2940          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
     2941                                source, dest);
     2942      }
    28832943
    28842944      ASSERT (!prim_fail, "failure in integer factorizer");
    … …  
    30153075          source= CFList();
    30163076          dest= CFList();
    3017           int num_vars= tmin (getNumVars(A), getNumVars(B));
    3018           num_vars--;
    3019 
    3020           choose_extension (d, num_vars, V_buf);
     3077
     3078          Variable V_buf3= V_buf;
     3079          V_buf= chooseExtension (V_buf);
    30213080          bool prim_fail= false;
    30223081          Variable V_buf2;
    30233082          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     3083
     3084          if (V_buf3 != alpha)
     3085          {
     3086            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     3087            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     3088            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
     3089                                 source, dest);
     3090            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
     3091            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
     3092            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
     3093                                source, dest);
     3094            for (CFListIterator i= l; i.hasItem(); i++)
     3095              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
     3096                                    source, dest);
     3097          }
    30243098
    30253099          ASSERT (!prim_fail, "failure in integer factorizer");
    … …  
    32793353  {
    32803354    if (inextension)
    3281       random_element= randomElement (m, alpha, l, fail);
     3355      random_element= randomElement (m, V_buf, l, fail);
    32823356    else
    32833357      random_element= FpRandomElement (m, l, fail);
    … …  
    33273401        deg++;
    33283402      } while (fail);
     3403      V_buf= alpha;
    33293404      list= CFList();
    33303405      TIMING_START (gcd_recursion);
    … …  
    33403415      source= CFList();
    33413416      dest= CFList();
    3342       int num_vars= tmin (getNumVars(A), getNumVars(B));
    3343       num_vars--;
    3344       V_buf= alpha;
    3345       choose_extension (d, num_vars, V_buf);
     3417
     3418      Variable V_buf3= V_buf;
     3419      V_buf= chooseExtension (V_buf);
    33463420      bool prim_fail= false;
    33473421      Variable V_buf2;
    33483422      CanonicalForm prim_elem, im_prim_elem;
    33493423      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     3424
     3425      if (V_buf3 != alpha)
     3426      {
     3427        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     3428        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     3429        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
     3430                             dest);
     3431        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
     3432        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
     3433        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
     3434                            dest);
     3435        for (CFListIterator i= l; i.hasItem(); i++)
     3436          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
     3437                                source, dest);
     3438      }
    33503439
    33513440      ASSERT (!prim_fail, "failure in integer factorizer");
    … …  
    33813470                            "time for recursive call: ");
    33823471      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
    3383       alpha= V_buf;
    33843472    }
    33853473
    … …  
    35123600            deg++;
    35133601          } while (fail);
     3602          V_buf= alpha;
    35143603          list= CFList();
    35153604          TIMING_START (gcd_recursion);
    … …  
    35323621          source= CFList();
    35333622          dest= CFList();
    3534           int num_vars= tmin (getNumVars(A), getNumVars(B));
    3535           num_vars--;
    3536           V_buf= alpha;
    3537           choose_extension (d, num_vars, V_buf);
     3623
     3624          Variable V_buf3= V_buf;
     3625          V_buf= chooseExtension (V_buf);
    35383626          bool prim_fail= false;
    35393627          Variable V_buf2;
    35403628          CanonicalForm prim_elem, im_prim_elem;
    35413629          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     3630
     3631          if (V_buf3 != alpha)
     3632          {
     3633            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     3634            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     3635            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
     3636                                 source, dest);
     3637            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
     3638            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
     3639            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
     3640                                source, dest);
     3641            for (CFListIterator i= l; i.hasItem(); i++)
     3642              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
     3643                                    source, dest);
     3644          }
    35423645
    35433646          ASSERT (!prim_fail, "failure in integer factorizer");
    … …  
    35803683                                "time for recursive call: ");
    35813684          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
    3582           alpha= V_buf;
    35833685        }
    35843686
    … …  
    38223924
    38233925static inline
    3824 bool Hensel_P (const CanonicalForm & U, CFArray & G, const Evaluation & A,
    3825                 const Variable & x, const CFArray& LCs )
     3926Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
     3927                    CFMap & N, const Evaluation& A)
     3928{
     3929  int n= F.level();
     3930  int * degsf= new int [n + 1];
     3931
     3932  for (int i = 0; i <= n; i++)
     3933    degsf[i]= 0;
     3934
     3935  degsf= degrees (F, degsf);
     3936
     3937  Evaluation result= Evaluation (A.min(), A.max());
     3938  ASSERT (A.min() == 2, "expected A.min() == 2");
     3939  ASSERT (A.max() == n, "expected A.max() == n");
     3940  int max_deg;
     3941  int k= n;
     3942  int l= 1;
     3943  int i= 2;
     3944  int pos= 2;
     3945  while (k > 1)
     3946  {
     3947    max_deg= degsf [i];
     3948    while (max_deg == 0)
     3949    {
     3950      i++;
     3951      max_deg= degsf [i];
     3952    }
     3953    l= i;
     3954    for (int j= i + 1; j <=  n; j++)
     3955    {
     3956      if (degsf[j] > max_deg)
     3957      {
     3958        max_deg= degsf[j];
     3959        l= j;
     3960      }
     3961    }
     3962
     3963    if (l <= n)
     3964    {
     3965      if (l != pos)
     3966      {
     3967        result.setValue (pos, A [l]);
     3968        M.newpair (Variable (l), Variable (pos));
     3969        N.newpair (Variable (pos), Variable (l));
     3970        degsf[l]= 0;
     3971        l= 2;
     3972        if (k == 2 && n == 3)
     3973        {
     3974          result.setValue (l, A [pos]);
     3975          M.newpair (Variable (pos), Variable (l));
     3976          N.newpair (Variable (l), Variable (pos));
     3977          degsf[pos]= 0;
     3978        }
     3979      }
     3980      else
     3981      {
     3982        result.setValue (l, A [l]);
     3983        degsf [l]= 0;
     3984      }
     3985    }
     3986    pos++;
     3987    k--;
     3988    l= 2;
     3989  }
     3990
     3991  delete [] degsf;
     3992
     3993  return result;
     3994}
     3995
     3996static inline
     3997int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
     3998                const Variable & x, const CFArray& LeadCoeffs )
    38263999{
    38274000  CFList factors;
    38284001  factors.append (G[1]);
    38294002  factors.append (G[2]);
     4003
     4004  CFMap NN, MM;
     4005  Evaluation A= optimize4Lift (UU, MM, NN, AA);
     4006
     4007  CanonicalForm U= MM (UU);
     4008  CFArray LCs= CFArray (1,2);
     4009  LCs [1]= MM (LeadCoeffs [1]);
     4010  LCs [2]= MM (LeadCoeffs [2]);
     4011
    38304012  CFList evaluation;
    38314013  for (int i= A.min(); i <= A.max(); i++)
    … …  
    38334015  CFList UEval;
    38344016  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
     4017
     4018  if (size (shiftedU)/getNumVars (U) > 500)
     4019    return -1;
     4020
    38354021  CFArray shiftedLCs= CFArray (2);
    38364022  CFList shiftedLCsEval1, shiftedLCsEval2;
    38374023  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
    38384024  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
    3839   CanonicalForm LcUEval= LC (UEval.getFirst(), x);
    38404025  factors.insert (1);
    38414026  int liftBound= degree (UEval.getLast(), 2) + 1;
    … …  
    38494034                 lcs, false);
    38504035
    3851   bool noChange= true;
    38524036  for (CFListIterator i= factors; i.hasItem(); i++)
    38534037  {
    3854     if (degree (i.getItem(), 2) != 0)
    3855       noChange= false;
    3856   }
    3857   if (noChange)
    3858     return !noChange;
     4038    if (!fdivides (i.getItem(), UEval.getFirst()))
     4039      return 0;
     4040  }
     4041
    38594042  int * liftBounds;
    3860   noChange= false;
     4043  bool noOneToOne;
    38614044  if (U.level() > 2)
    38624045  {
    … …  
    38664049      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
    38674050    factors= henselLift2 (UEval, factors, liftBounds, U.level() - 1, false,
    3868                           shiftedLCsEval1, shiftedLCsEval2, Pi, diophant, noChange);
    3869    delete [] liftBounds;
    3870    if (noChange)
    3871      return !noChange;
     4051                          shiftedLCsEval1, shiftedLCsEval2, Pi, diophant,
     4052                          noOneToOne);
     4053    delete [] liftBounds;
     4054    if (noOneToOne)
     4055      return 0;
    38724056  }
    38734057  G[1]= factors.getFirst();
    … …  
    38754059  G[1]= myReverseShift (G[1], evaluation);
    38764060  G[2]= myReverseShift (G[2], evaluation);
    3877   return true;
     4061  G[1]= NN (G[1]);
     4062  G[2]= NN (G[2]);
     4063  return 1;
    38784064}
    38794065
    … …  
    38814067bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
    38824068                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
    3883                  FFREvaluation & b, int delta, int degF, int degG, int maxeval,
    3884                  int & count)
     4069                 REvaluation & b, int delta, int degF, int degG, int maxeval,
     4070                 int & count, int& k, int bound, int& l)
    38854071{
    38864072  if( count == 0 && delta != 0)
    … …  
    38894075      return false;
    38904076  }
     4077  CanonicalForm evalPoint;
    38914078  if (count > 0)
    38924079  {
    3893     b.nextpoint();
     4080    b.nextpoint(k);
     4081    l++;
     4082    if (l > bound)
     4083    {
     4084      l= 1;
     4085      k++;
     4086      if (k > tmax (F.level(), G.level()) - 1)
     4087        return false;
     4088      b.nextpoint (k);
     4089    }
    38944090    if (count++ > maxeval)
    38954091      return false;
    … …  
    39134109      }
    39144110    }
    3915     b.nextpoint();
     4111    if (k == 0)
     4112      k++;
     4113    b.nextpoint(k);
     4114    l++;
     4115    if (l > bound)
     4116    {
     4117      l= 1;
     4118      k++;
     4119      if (k > tmax (F.level(), G.level()) - 1)
     4120        return false;
     4121      b.nextpoint (k);
     4122    }
    39164123    if( count++ > maxeval )
    39174124      return false;
    … …  
    39224129static int maxNumEval= 200;
    39234130static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
    3924 static int sizePerVars2= 290; //try psr gcd if size/#variables is less
    39254131
    39264132/// Extended Zassenhaus GCD for finite fields
    … …  
    39434149                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
    39444150  count= 0; // number of eval. used
    3945   FFREvaluation b, bt;
    3946   bool gcdfound = false;
     4151  REvaluation b, bt;
     4152  int gcdfound = 0;
    39474153  Variable x = Variable(1);
    39484154
    … …  
    39714177
    39724178  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
    3973   Variable a,bv;
     4179  Variable a, oldA;
    39744180  int sizeF= size (F);
    39754181  int sizeG= size (G);
    … …  
    39904196  }
    39914197
     4198  bool passToGF= false;
     4199  bool extOfExt= false;
     4200  int p= getCharacteristic();
     4201  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
     4202  int k= 1;
     4203  CanonicalForm primElem, imPrimElem;
     4204  CFList source, dest;
     4205  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
     4206  {
     4207    if (p == 2)
     4208      setCharacteristic (2, 6, 'Z');
     4209    else if (p == 3)
     4210      setCharacteristic (3, 4, 'Z');
     4211    else if (p == 5 || p == 7)
     4212      setCharacteristic (p, 3, 'Z');
     4213    else
     4214      setCharacteristic (p, 2, 'Z');
     4215    passToGF= true;
     4216    F= F.mapinto();
     4217    G= G.mapinto();
     4218    maxeval= 2*ipower (p, getGFDegree());
     4219  }
     4220  else if (CFFactory::gettype() == GaloisFieldDomain &&
     4221           ipower (p , getGFDegree()) < 50)
     4222  {
     4223    k= getGFDegree();
     4224    if (ipower (p, 2*k) > 50)
     4225      setCharacteristic (p, 2*k, gf_name);
     4226    else
     4227      setCharacteristic (p, 3*k, gf_name);
     4228    F= GFMapUp (F, k);
     4229    G= GFMapUp (G, k);
     4230    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
     4231  }
     4232  else if (p < 50 && algExtension && !CFFactory::gettype() == GaloisFieldDomain)
     4233  {
     4234    int d= degree (getMipo (a));
     4235    oldA= a;
     4236    Variable v2;
     4237    if (p == 2 && d < 6)
     4238    {
     4239      zz_p::init (p);
     4240      bool primFail= false;
     4241      Variable vBuf;
     4242      primElem= primitiveElement (a, vBuf, primFail);
     4243      ASSERT (!primFail, "failure in integer factorizer");
     4244      if (d < 3)
     4245      {
     4246        zz_pX NTLIrredpoly;
     4247        BuildIrred (NTLIrredpoly, d*3);
     4248        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
     4249        v2= rootOf (newMipo);
     4250      }
     4251      else
     4252      {
     4253        zz_pX NTLIrredpoly;
     4254        BuildIrred (NTLIrredpoly, d*2);
     4255        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
     4256        v2= rootOf (newMipo);
     4257      }
     4258      imPrimElem= mapPrimElem (primElem, a, v2);
     4259      extOfExt= true;
     4260    }
     4261    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
     4262    {
     4263      zz_p::init (p);
     4264      bool primFail= false;
     4265      Variable vBuf;
     4266      primElem= primitiveElement (a, vBuf, primFail);
     4267      ASSERT (!primFail, "failure in integer factorizer");
     4268      zz_pX NTLIrredpoly;
     4269      BuildIrred (NTLIrredpoly, d*2);
     4270      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
     4271      v2= rootOf (newMipo);
     4272      imPrimElem= mapPrimElem (primElem, a, v2);
     4273      extOfExt= true;
     4274    }
     4275    if (extOfExt)
     4276    {
     4277      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
     4278      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
     4279      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
     4280      a= v2;
     4281    }
     4282  }
     4283
    39924284  lcF = LC( F, x ); lcG = LC( G, x );
    39934285  lcD = gcd( lcF, lcG );
    … …  
    39954287  delta = 0;
    39964288  degF = degree( F, x ); degG = degree( G, x );
    3997   if(hasFirstAlgVar(F,a))
    3998   {
    3999     if(hasFirstAlgVar(G,bv))
    4000     {
    4001       if(bv.level() > a.level())
    4002         a = bv;
    4003     }
    4004     b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
    4005   }
    4006   else // F not in extension
    4007   {
    4008     if(hasFirstAlgVar(G,a))
    4009       b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
     4289
     4290  if(hasFirstAlgVar(G,a))
     4291    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
     4292  else
     4293  { // both not in extension given by algebraic variable
     4294    if (CFFactory::gettype() != GaloisFieldDomain)
     4295      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
    40104296    else
    4011     { // both not in extension given by algebraic variable
    4012       if (CFFactory::gettype() != GaloisFieldDomain)
    4013         b = FFREvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
    4014       else
    4015         b = FFREvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
    4016     }
    4017   }
     4297      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
     4298  }
     4299
    40184300  CanonicalForm cand;
     4301  CanonicalForm result;
     4302  int o, t;
     4303  o= 0;
     4304  t= 1;
     4305  int goodPointCount= 0;
    40194306  while( !gcdfound )
    40204307  {
    4021     if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count ))
     4308    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
     4309         maxeval/maxNumVars, t ))
    40224310    { // too many eval. used --> try another method
    4023       if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2)
    4024       {
    4025         CanonicalForm dummy1, dummy2;
    4026         Variable y= Variable (tmax (F.level(), G.level()));
    4027         Variable z= Variable (smallestDegLev);
    4028         dummy1= swapvar (F, x, y);
    4029         dummy2= swapvar (G, x, y);
    4030         if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G)))
    4031         {
    4032           Off (SW_USE_EZGCD_P);
    4033           CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y);
    4034           On (SW_USE_EZGCD_P);
    4035           return N (d*result);
    4036         }
    4037       }
    4038       if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
    4039         return N (d*sparseGCDFq (F, G, a));
    4040       if (CFFactory::gettype() == GaloisFieldDomain)
    4041         return N (d*GCD_GF (F, G));
    4042       else
    4043         return N (d*sparseGCDFp (F,G));
     4311      Off (SW_USE_EZGCD_P);
     4312      result= gcd (F,G);
     4313      On (SW_USE_EZGCD_P);
     4314      if (passToGF)
     4315      {
     4316        Variable alpha= rootOf (gf_mipo);
     4317        setCharacteristic (p);
     4318        result= GF2FalphaRep (result, alpha);
     4319      }
     4320      if (k > 1)
     4321      {
     4322        result= GFMapDown (result, k);
     4323        setCharacteristic (p, k, gf_name);
     4324      }
     4325      if (extOfExt)
     4326        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     4327      return N (d*result);
    40444328    }
    40454329    delta = degree( Db );
    40464330    if( delta == 0 )
     4331    {
     4332      if (passToGF)
     4333        setCharacteristic (p);
     4334      if (k > 1)
     4335        setCharacteristic (p, k, gf_name);
    40474336      return N (d);
     4337    }
    40484338    while( true )
    40494339    {
    40504340      bt = b;
    4051       if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count ))
     4341      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
     4342           maxeval/maxNumVars, t ))
    40524343      { // too many eval. used --> try another method
    4053         if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2)
    4054         {
    4055           CanonicalForm dummy1, dummy2;
    4056           Variable y= Variable (tmax (F.level(), G.level()));
    4057           Variable z= Variable (smallestDegLev);
    4058           dummy1= swapvar (F, x, y);
    4059           dummy2= swapvar (G, x, y);
    4060           if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G)))
    4061           {
    4062             Off (SW_USE_EZGCD_P);
    4063             CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y);
    4064             On (SW_USE_EZGCD_P);
    4065             return N (d*result);
    4066           }
    4067         }
    4068         if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
    4069           return N (d*sparseGCDFq (F, G, a));
    4070         if (CFFactory::gettype() == GaloisFieldDomain)
    4071           return N (d*GCD_GF (F, G));
    4072         else
    4073           return N (d*sparseGCDFp (F,G));
     4344        Off (SW_USE_EZGCD_P);
     4345        result= gcd (F,G);
     4346        On (SW_USE_EZGCD_P);
     4347        if (passToGF)
     4348        {
     4349          Variable alpha= rootOf (gf_mipo);
     4350          setCharacteristic (p);
     4351          result= GF2FalphaRep (result, alpha);
     4352        }
     4353        if (k > 1)
     4354        {
     4355          result= GFMapDown (result, k);
     4356          setCharacteristic (p, k, gf_name);
     4357        }
     4358        if (extOfExt)
     4359          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     4360        return N (d*result);
    40744361      }
    40754362      int dd = degree( Dbt );
    40764363      if( dd == 0 )
     4364      {
     4365        if (passToGF)
     4366          setCharacteristic (p);
     4367        if (k > 1)
     4368          setCharacteristic (p, k, gf_name);
    40774369        return N (d);
     4370      }
    40784371      if( dd == delta )
    4079         break;
     4372      {
     4373        goodPointCount++;
     4374        if (goodPointCount == 5)
     4375          break;
     4376      }
    40804377      if( dd < delta )
    40814378      {
     4379        goodPointCount= 0;
    40824380        delta = dd;
    40834381        b = bt;
    40844382        Db = Dbt; Fb = Fbt; Gb = Gbt;
    40854383      }
    4086     }
    4087     if( degF <= degG && delta == degF && fdivides( F, G ) )
    4088       return N (d*F);
    4089     if( degG < degF && delta == degG && fdivides( G, F ) )
    4090       return N (d*G);
     4384      if (delta == degF)
     4385      {
     4386        if (degF <= degG && fdivides (F, G))
     4387        {
     4388          if (passToGF)
     4389          {
     4390            Variable alpha= rootOf (gf_mipo);
     4391            setCharacteristic (p);
     4392            F= GF2FalphaRep (F, alpha);
     4393          }
     4394          if (k > 1)
     4395          {
     4396            F= GFMapDown (F, k);
     4397            setCharacteristic (p, k, gf_name);
     4398          }
     4399          if (extOfExt)
     4400            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
     4401          return N (d*F);
     4402        }
     4403        else
     4404          delta--;
     4405      }
     4406      else if (delta == degG)
     4407      {
     4408        if (degG <= degF && fdivides (G, F))
     4409        {
     4410          if (passToGF)
     4411          {
     4412            Variable alpha= rootOf (gf_mipo);
     4413            setCharacteristic (p);
     4414            G= GF2FalphaRep (G, alpha);
     4415          }
     4416          if (k > 1)
     4417          {
     4418            G= GFMapDown (G, k);
     4419            setCharacteristic (p, k, gf_name);
     4420          }
     4421          if (extOfExt)
     4422            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
     4423          return N (d*G);
     4424        }
     4425        else
     4426          delta--;
     4427      }
     4428    }
    40914429    if( delta != degF && delta != degG )
    40924430    {
    … …  
    41214459      else // special case handling
    41224460      {
    4123         if (sizeF/maxNumVars < sizePerVars2 && sizeG/maxNumVars < sizePerVars2)
    4124         {
    4125           CanonicalForm dummy1, dummy2;
    4126           Variable y= Variable (tmax (F.level(), G.level()));
    4127           Variable z= Variable (smallestDegLev);
    4128           dummy1= swapvar (F, x, y);
    4129           dummy2= swapvar (G, x, y);
    4130           if (size (psr (dummy1, dummy2, y)) <= tmin (size (F), size (G)))
    4131           {
    4132             Off (SW_USE_EZGCD_P);
    4133             CanonicalForm result= swapvar (gcd (dummy1,dummy2), x, y);
    4134             On (SW_USE_EZGCD_P);
    4135             return N (d*result);
    4136           }
    4137         }
    4138         if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
    4139           return N (d*sparseGCDFq (F, G, a));
    4140         if (CFFactory::gettype() == GaloisFieldDomain)
    4141           return N (d*GCD_GF (F, G));
    4142         else
    4143           return N (d*sparseGCDFp (F,G));
     4461        Off (SW_USE_EZGCD_P);
     4462        result= gcd (F,G);
     4463        On (SW_USE_EZGCD_P);
     4464        if (passToGF)
     4465        {
     4466          Variable alpha= rootOf (gf_mipo);
     4467          setCharacteristic (p);
     4468          result= GF2FalphaRep (result, alpha);
     4469        }
     4470        if (k > 1)
     4471        {
     4472          result= GFMapDown (result, k);
     4473          setCharacteristic (p, k, gf_name);
     4474        }
     4475        if (extOfExt)
     4476          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     4477        return N (d*result);
    41444478      }
    41454479      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
    41464480      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
    41474481
     4482      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
     4483      {
     4484        if (algExtension)
     4485        {
     4486          result= GCD_Fp_extension (F, G, a);
     4487          if (extOfExt)
     4488            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     4489          return N (d*result);
     4490        }
     4491        if (CFFactory::gettype() == GaloisFieldDomain)
     4492        {
     4493          result= GCD_GF (F, G);
     4494          if (passToGF)
     4495          {
     4496            Variable alpha= rootOf (gf_mipo);
     4497            setCharacteristic (p);
     4498            result= GF2FalphaRep (result, alpha);
     4499          }
     4500          if (k > 1)
     4501          {
     4502            result= GFMapDown (result, k);
     4503            setCharacteristic (p, k, gf_name);
     4504          }
     4505          return N (d*result);
     4506        }
     4507        else
     4508          return N (d*GCD_small_p (F,G));
     4509      }
     4510
    41484511      gcdfound= Hensel_P (B*lcD, DD, b, x, lcDD);
    41494512
    4150       if( gcdfound )
     4513      if (gcdfound == -1)
     4514      {
     4515        Off (SW_USE_EZGCD_P);
     4516        result= gcd (F,G);
     4517        On (SW_USE_EZGCD_P);
     4518        if (passToGF)
     4519        {
     4520          Variable alpha= rootOf (gf_mipo);
     4521          setCharacteristic (p);
     4522          result= GF2FalphaRep (result, alpha);
     4523        }
     4524        if (k > 1)
     4525        {
     4526          result= GFMapDown (result, k);
     4527          setCharacteristic (p, k, gf_name);
     4528        }
     4529        if (extOfExt)
     4530          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     4531        return N (d*result);
     4532      }
     4533
     4534      if (gcdfound == 1)
    41514535      {
    41524536        cand = DD[2] / content( DD[2], Variable(1) );
    4153         if( B_is_F )
    4154           gcdfound = fdivides( cand, G );
    4155         else
    4156           gcdfound = fdivides( cand, F );
    4157       }
    4158     }
    4159     delta= 0;
     4537        gcdfound = fdivides( cand, G ) && fdivides ( cand, F );
     4538
     4539        if (passToGF && gcdfound)
     4540        {
     4541          Variable alpha= rootOf (gf_mipo);
     4542          setCharacteristic (p);
     4543          cand= GF2FalphaRep (cand, alpha);
     4544        }
     4545        if (k > 1 && gcdfound)
     4546        {
     4547          cand= GFMapDown (cand, k);
     4548          setCharacteristic (p, k, gf_name);
     4549        }
     4550        if (extOfExt && gcdfound)
     4551          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
     4552      }
     4553    }
     4554    delta--;
     4555    goodPointCount= 0;
    41604556  }
    41614557  return N (d*cand);
Note: See TracChangeset for help on using the changeset viewer.