Changeset 597783 in git for factory/cf_gcd_smallp.cc


Ignore:
Timestamp:
May 2, 2012, 1:41:01 PM (12 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
e2666785b249262bbfdcdc9e4ca1a39a4054e994
Parents:
cb7827a842fd2e9878498f940e9f7735fdbe488e
git-author:
Martin Lee <martinlee84@web.de>2012-05-02 13:41:01+02:00
git-committer:
Martin Lee <martinlee84@web.de>2012-05-07 14:16:48+02:00
Message:
chg: compute cofactors of gcd and use them to test divisibility
chg: do not use bounds in gcd over Fq
File:
1 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_gcd_smallp.cc

    rcb7827 r597783  
    425425}
    426426
     427CanonicalForm
     428GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
     429                  CanonicalForm& coF, CanonicalForm& coG,
     430                  Variable & alpha, CFList& l, bool& topLevel);
     431
     432CanonicalForm
     433GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
     434                  Variable & alpha, CFList& l, bool& topLevel)
     435{
     436  CanonicalForm dummy1, dummy2;
     437  CanonicalForm result= GCD_Fp_extension (F, G, dummy1, dummy2, alpha, l,
     438                                          topLevel);
     439  return result;
     440}
     441
    427442/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
    428443/// l and topLevel are only used internally, output is monic
    … …  
    431446CanonicalForm
    432447GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
     448                  CanonicalForm& coF, CanonicalForm& coG,
    433449                  Variable & alpha, CFList& l, bool& topLevel)
    434450{
    435451  CanonicalForm A= F;
    436452  CanonicalForm B= G;
    437   if (F.isZero() && degree(G) > 0) return G/Lc(G);
    438   else if (G.isZero() && degree (F) > 0) return F/Lc(F);
    439   else if (F.isZero() && G.isZero()) return F.genOne();
    440   if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
    441   if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
    442   if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
    443   if (F == G) return F/Lc(F);
     453  if (F.isZero() && degree(G) > 0)
     454  {
     455    coF= 0;
     456    coG= Lc (G);
     457    return G/Lc(G);
     458  }
     459  else if (G.isZero() && degree (F) > 0)
     460  {
     461    coF= Lc (F);
     462    coG= 0;
     463    return F/Lc(F);
     464  }
     465  else if (F.isZero() && G.isZero())
     466  {
     467    coF= coG= 0;
     468    return F.genOne();
     469  }
     470  if (F.inBaseDomain() || G.inBaseDomain())
     471  {
     472    coF= F;
     473    coG= G;
     474    return F.genOne();
     475  }
     476  if (F.isUnivariate() && fdivides(F, G, coG))
     477  {
     478    coF= Lc (F);
     479    return F/Lc(F);
     480  }
     481  if (G.isUnivariate() && fdivides(G, F, coF))
     482  {
     483    coG= Lc (G);
     484    return G/Lc(G);
     485  }
     486  if (F == G)
     487  {
     488    coF= coG= Lc (F);
     489    return F/Lc(F);
     490  }
    444491
    445492  CFMap M,N;
    446493  int best_level= myCompress (A, B, M, N, topLevel);
    447494
    448   if (best_level == 0) return B.genOne();
     495  if (best_level == 0)
     496  {
     497    coF= F;
     498    coG= G;
     499    return B.genOne();
     500  }
    449501
    450502  A= M(A);
    … …  
    455507  //univariate case
    456508  if (A.isUnivariate() && B.isUnivariate())
    457     return N (gcd(A,B));
     509  {
     510    CanonicalForm result= gcd (A, B);
     511    coF= N (A/result);
     512    coG= N (B/result);
     513    return N (result);
     514  }
    458515
    459516  CanonicalForm cA, cB;    // content of A and B
    … …  
    484541  if (compressConvexDense)
    485542  {
     543    CanonicalForm bufcA= cA;
     544    CanonicalForm bufcB= cB;
    486545    cA= content (ppA, 1);
    487546    cB= content (ppB, 1);
    … …  
    489548    ppB /= cB;
    490549    gcdcAcB *= gcd (cA, cB);
     550    cA *= bufcA;
     551    cB *= bufcB;
    491552    if (ppA.isUnivariate() || ppB.isUnivariate())
    492553    {
    493554      if (ppA.level() == ppB.level())
    494         return N (gcd (ppA, ppB)*gcdcAcB);
     555      {
     556        CanonicalForm result= gcd (ppA, ppB);
     557        coF= N ((ppA/result)*(cA/gcdcAcB));
     558        coG= N ((ppB/result)*(cB/gcdcAcB));
     559        return N (result*gcdcAcB);
     560      }
    495561      else
     562      {
     563        coF= N (ppA*(cA/gcdcAcB));
     564        coG= N (ppB*(cB/gcdcAcB));
    496565        return N (gcdcAcB);
     566      }
    497567    }
    498568
    … …  
    511581    {
    512582      if (ppA.level() == ppB.level())
    513         return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
     583      {
     584        CanonicalForm result= gcd (ppA, ppB);
     585        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
     586        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
     587        return N (decompress (result, MM, V)*gcdcAcB);
     588      }
    514589      else
     590      {
     591        coF= N (decompress (ppA, MM, V));
     592        coG= N (decompress (ppB, MM, V));
    515593        return N (gcdcAcB);
     594      }
    516595    }
    517596  }
    … …  
    523602  lcB= uni_lcoeff (ppB);
    524603
    525   if (fdivides (lcA, lcB))
     604  /*if (fdivides (lcA, lcB))
    526605  {
    527606    if (fdivides (A, B))
    … …  
    532611    if (fdivides (B, A))
    533612      return G/Lc(G);
    534   }
     613  }*/
    535614
    536615  gcdlcAlcB= gcd (lcA, lcB);
    … …  
    539618
    540619  if (d == 0)
     620  {
     621    coF= N (ppA*(cA/gcdcAcB));
     622    coG= N (ppB*(cB/gcdcAcB));
    541623    return N(gcdcAcB);
     624  }
    542625
    543626  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
    … …  
    545628    d= d0;
    546629  if (d == 0)
     630  {
     631    coF= N (ppA*(cA/gcdcAcB));
     632    coG= N (ppB*(cB/gcdcAcB));
    547633    return N(gcdcAcB);
     634  }
    548635
    549636  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
    550   CanonicalForm newtonPoly;
     637  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
     638                coG_m, ppCoF, ppCoG;
    551639
    552640  newtonPoly= 1;
    553641  m= gcdlcAlcB;
    554642  G_m= 0;
     643  coF= 0;
     644  coG= 0;
    555645  H= 0;
    556646  bool fail= false;
    … …  
    560650  CanonicalForm prim_elem, im_prim_elem;
    561651  CFList source, dest;
    562   int bound= tmin (degree (ppA, 1), degree (ppB, 1));
    563   int count= 0;
     652  int bound1= degree (ppA, 1);
     653  int bound2= degree (ppB, 1);
    564654  do
    565655  {
    566     random_element= randomElement (m, V_buf, l, fail);
     656    random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
    567657    if (fail)
    568658    {
    … …  
    579669      {
    580670        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
    581         G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     671        G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
     672        coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
     673        coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
    582674        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
    583675                             source, dest);
    … …  
    586678        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
    587679                            source, dest);
     680        lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
     681        lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
    588682        for (CFListIterator i= l; i.hasItem(); i++)
    589683          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
    … …  
    605699      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
    606700      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     701      coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     702      coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
    607703      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
    608704                          source, dest);
    … …  
    611707      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
    612708                         source, dest);
     709      lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     710      lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
    613711
    614712      fail= false;
    615       random_element= randomElement (m, V_buf, l, fail );
     713      random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
    616714      DEBOUTLN (cerr, "fail= " << fail);
    617715      CFList list;
    618716      TIMING_START (gcd_recursion);
    619717      G_random_element=
    620       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
     718      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
     719                        coF_random_element, coG_random_element, V_buf,
    621720                        list, topLevel);
    622721      TIMING_END_AND_PRINT (gcd_recursion,
    … …  
    629728      TIMING_START (gcd_recursion);
    630729      G_random_element=
    631       GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf,
     730      GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x),
     731                        coF_random_element, coG_random_element, V_buf,
    632732                        list, topLevel);
    633733      TIMING_END_AND_PRINT (gcd_recursion,
    … …  
    643743
    644744    if (d0 == 0)
     745    {
     746      coF= N (ppA*(cA/gcdcAcB));
     747      coG= N (ppB*(cB/gcdcAcB));
    645748      return N(gcdcAcB);
     749    }
    646750    if (d0 >  d)
    647751    {
    … …  
    654758    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
    655759    * G_random_element;
     760    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
     761                        *coF_random_element;
     762    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
     763                        *coG_random_element;
    656764
    657765    if (!G_random_element.inCoeffDomain())
    … …  
    667775      G_m= 0;
    668776      d= d0;
    669       count= 0;
     777      coF_m= 0;
     778      coG_m= 0;
    670779    }
    671780
    672781    TIMING_START (newton_interpolation);
    673782    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
     783    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
     784    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
    674785    TIMING_END_AND_PRINT (newton_interpolation,
    675786                          "time for newton interpolation: ");
    676787
    677     count++;
    678788    //termination test
    679     if (uni_lcoeff (H) == gcdlcAlcB)
    680     {
    681       cH= uni_content (H);
     789    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
     790    {
     791      if (gcdlcAlcB.isOne())
     792        cH= 1;
     793      else
     794        cH= uni_content (H);
    682795      ppH= H/cH;
     796      CanonicalForm lcppH= gcdlcAlcB/cH;
     797      CanonicalForm ccoF= lcA/lcppH;
     798      ccoF /= Lc (ccoF);
     799      CanonicalForm ccoG= lcB/lcppH;
     800      ccoG /= Lc (ccoG);
     801      ppCoF= coF/ccoF;
     802      ppCoG= coG/ccoG;
    683803      if (inextension)
    684804      {
    685         CFList u, v;
    686         //maybe it's better to test if ppH is an element of F(\alpha) before
    687         //mapping down
    688         DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
    689         ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
    690         ppH /= Lc(ppH);
    691         DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
    692         if ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB)))
    693         {
     805        if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
     806             (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
     807             (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) ||
     808             (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
     809        {
     810          CFList u, v;
    694811          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
    695812          ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
     813          ppCoF= mapDown (ppCoF, prim_elem, im_prim_elem, alpha, u, v);
     814          ppCoF= mapDown (ppCoG, prim_elem, im_prim_elem, alpha, u, v);
    696815          ppH /= Lc(ppH);
    697816          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
    698817          if (compressConvexDense)
     818          {
    699819            ppH= decompress (ppH, MM, V);
     820            ppCoF= decompress (ppCoF, MM, V);
     821            ppCoG= decompress (ppCoG, MM, V);
     822          }
     823          coF= N ((cA/gcdcAcB)*ppCoF);
     824          coG= N ((cB/gcdcAcB)*ppCoG);
    700825          return N(gcdcAcB*ppH);
    701826        }
    702827      }
    703       else if ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB)))
     828      else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
     829                (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
     830                (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) ||
     831                (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
    704832      {
    705833        if (compressConvexDense)
     834        {
    706835          ppH= decompress (ppH, MM, V);
     836          ppCoF= decompress (ppCoF, MM, V);
     837          ppCoG= decompress (ppCoG, MM, V);
     838        }
     839        coF= N ((cA/gcdcAcB)*ppCoF);
     840        coG= N ((cB/gcdcAcB)*ppCoG);;
    707841        return N(gcdcAcB*ppH);
    708842      }
    … …  
    710844
    711845    G_m= H;
     846    coF_m= coF;
     847    coG_m= coG;
    712848    newtonPoly= newtonPoly*(x - random_element);
    713849    m= m*(x - random_element);
    … …  
    755891}
    756892
     893CanonicalForm
     894GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
     895        CanonicalForm& coF, CanonicalForm& coG,
     896        CFList& l, bool& topLevel);
     897
     898CanonicalForm
     899GCD_GF (const CanonicalForm& F, const CanonicalForm& G, CFList& l,
     900        bool& topLevel)
     901{
     902  CanonicalForm dummy1, dummy2;
     903  CanonicalForm result= GCD_GF (F, G, dummy1, dummy2, l, topLevel);
     904  return result;
     905}
     906
    757907/// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
    758908/// Computer Algebra" by Geddes, Czapor, Labahn
    … …  
    760910/// faster field arithmetics, however it might fail if the input is large since
    761911/// the size of the base field is bounded by 2^16, output is monic
    762 CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
     912CanonicalForm
     913GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
     914        CanonicalForm& coF, CanonicalForm& coG,
    763915        CFList& l, bool& topLevel)
    764916{
    765917  CanonicalForm A= F;
    766918  CanonicalForm B= G;
    767   if (F.isZero() && degree(G) > 0) return G/Lc(G);
    768   else if (G.isZero() && degree (F) > 0) return F/Lc(F);
    769   else if (F.isZero() && G.isZero()) return F.genOne();
    770   if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
    771   if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
    772   if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
    773   if (F == G) return F/Lc(F);
     919  if (F.isZero() && degree(G) > 0)
     920  {
     921    coF= 0;
     922    coG= Lc (G);
     923    return G/Lc(G);
     924  }
     925  else if (G.isZero() && degree (F) > 0)
     926  {
     927    coF= Lc (F);
     928    coG= 0;
     929    return F/Lc(F);
     930  }
     931  else if (F.isZero() && G.isZero())
     932  {
     933    coF= coG= 0;
     934    return F.genOne();
     935  }
     936  if (F.inBaseDomain() || G.inBaseDomain())
     937  {
     938    coF= F;
     939    coG= G;
     940    return F.genOne();
     941  }
     942  if (F.isUnivariate() && fdivides(F, G, coG))
     943  {
     944    coF= Lc (F);
     945    return F/Lc(F);
     946  }
     947  if (G.isUnivariate() && fdivides(G, F, coF))
     948  {
     949    coG= Lc (G);
     950    return G/Lc(G);
     951  }
     952  if (F == G)
     953  {
     954    coF= coG= Lc (F);
     955    return F/Lc(F);
     956  }
    774957
    775958  CFMap M,N;
    776959  int best_level= myCompress (A, B, M, N, topLevel);
    777960
    778   if (best_level == 0) return B.genOne();
     961  if (best_level == 0)
     962  {
     963    coF= F;
     964    coG= G;
     965    return B.genOne();
     966  }
    779967
    780968  A= M(A);
    … …  
    785973  //univariate case
    786974  if (A.isUnivariate() && B.isUnivariate())
    787     return N (gcd (A, B));
     975  {
     976    CanonicalForm result= gcd (A, B);
     977    coF= N (A/result);
     978    coG= N (B/result);
     979    return N (result);
     980  }
    788981
    789982  CanonicalForm cA, cB;    // content of A and B
    … …  
    8141007  if (compressConvexDense)
    8151008  {
     1009    CanonicalForm bufcA= cA;
     1010    CanonicalForm bufcB= cB;
    8161011    cA= content (ppA, 1);
    8171012    cB= content (ppB, 1);
    … …  
    8191014    ppB /= cB;
    8201015    gcdcAcB *= gcd (cA, cB);
     1016    cA *= bufcA;
     1017    cB *= bufcB;
    8211018    if (ppA.isUnivariate() || ppB.isUnivariate())
    8221019    {
    8231020      if (ppA.level() == ppB.level())
    824         return N (gcd (ppA, ppB)*gcdcAcB);
     1021      {
     1022        CanonicalForm result= gcd (ppA, ppB);
     1023        coF= N ((ppA/result)*(cA/gcdcAcB));
     1024        coG= N ((ppB/result)*(cB/gcdcAcB));
     1025        return N (result*gcdcAcB);
     1026      }
    8251027      else
     1028      {
     1029        coF= N (ppA*(cA/gcdcAcB));
     1030        coG= N (ppB*(cB/gcdcAcB));
    8261031        return N (gcdcAcB);
     1032      }
    8271033    }
    8281034
    … …  
    8411047    {
    8421048      if (ppA.level() == ppB.level())
    843         return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
     1049      {
     1050        CanonicalForm result= gcd (ppA, ppB);
     1051        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
     1052        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
     1053        return N (decompress (result, MM, V)*gcdcAcB);
     1054      }
    8441055      else
     1056      {
     1057        coF= N (decompress (ppA, MM, V));
     1058        coG= N (decompress (ppB, MM, V));
    8451059        return N (gcdcAcB);
     1060      }
    8461061    }
    8471062  }
    … …  
    8531068  lcB= uni_lcoeff (ppB);
    8541069
    855   if (fdivides (lcA, lcB))
    856   {
    857     if (fdivides (A, B))
     1070  /*if (fdivides (lcA, lcB))
     1071  {
     1072    if (fdivides (ppA, ppB, coG))
     1073    {
     1074      coF= 1;
     1075      if (compressConvexDense)
     1076        coG= decompress (coG, MM, V);
     1077      coG= N (coG*(cB/gcdcAcB));
    8581078      return F;
     1079    }
    8591080  }
    8601081  if (fdivides (lcB, lcA))
    8611082  {
    862     if (fdivides (B, A))
     1083    if (fdivides (ppB, ppA, coF))
     1084    {
     1085      coG= 1;
     1086      if (compressConvexDense)
     1087        coF= decompress (coF, MM, V);
     1088      coF= N (coF*(cA/gcdcAcB));
    8631089      return G;
    864   }
     1090    }
     1091  }*/
    8651092
    8661093  gcdlcAlcB= gcd (lcA, lcB);
    … …  
    8681095  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
    8691096  if (d == 0)
     1097  {
     1098    coF= N (ppA*(cA/gcdcAcB));
     1099    coG= N (ppB*(cB/gcdcAcB));
    8701100    return N(gcdcAcB);
     1101  }
    8711102  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
    8721103  if (d0 < d)
    8731104    d= d0;
    8741105  if (d == 0)
     1106  {
     1107    coF= N (ppA*(cA/gcdcAcB));
     1108    coG= N (ppB*(cB/gcdcAcB));
    8751109    return N(gcdcAcB);
     1110  }
    8761111
    8771112  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
    878   CanonicalForm newtonPoly;
     1113  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
     1114                coG_m, ppCoF, ppCoG;
    8791115
    8801116  newtonPoly= 1;
    8811117  m= gcdlcAlcB;
    8821118  G_m= 0;
     1119  coF= 0;
     1120  coG= 0;
    8831121  H= 0;
    8841122  bool fail= false;
    885   topLevel= false;
     1123  //topLevel= false;
    8861124  bool inextension= false;
    8871125  int p=-1;
    … …  
    8901128  int expon;
    8911129  char gf_name_buf= gf_name;
    892   int bound= tmin (degree (ppA, 1), degree (ppB, 1));
    893   int count= 0;
     1130  int bound1= degree (ppA, 1);
     1131  int bound2= degree (ppB, 1);
    8941132  do
    8951133  {
    896     random_element= GFRandomElement (m, l, fail);
     1134    random_element= GFRandomElement (m*lcA*lcB, l, fail);
    8971135    if (fail)
    8981136    {
    … …  
    9151153      G_m= GFMapUp (G_m, kk);
    9161154      newtonPoly= GFMapUp (newtonPoly, kk);
     1155      coF_m= GFMapUp (coF_m, kk);
     1156      coG_m= GFMapUp (coG_m, kk);
    9171157      ppA= GFMapUp (ppA, kk);
    9181158      ppB= GFMapUp (ppB, kk);
    9191159      gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
    920       random_element= GFRandomElement (m, l, fail);
     1160      lcA= GFMapUp (lcA, kk);
     1161      lcB= GFMapUp (lcB, kk);
     1162      random_element= GFRandomElement (m*lcA*lcB, l, fail);
    9211163      DEBOUTLN (cerr, "fail= " << fail);
    9221164      CFList list;
    9231165      TIMING_START (gcd_recursion);
    9241166      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
     1167                                coF_random_element, coG_random_element,
    9251168                                list, topLevel);
    9261169      TIMING_END_AND_PRINT (gcd_recursion,
    … …  
    9331176      TIMING_START (gcd_recursion);
    9341177      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
     1178                                coF_random_element, coG_random_element,
    9351179                                list, topLevel);
    9361180      TIMING_END_AND_PRINT (gcd_recursion,
    … …  
    9491193      if (inextension)
    9501194        setCharacteristic (p, k, gf_name_buf);
     1195      coF= N (ppA*(cA/gcdcAcB));
     1196      coG= N (ppB*(cB/gcdcAcB));
    9511197      return N(gcdcAcB);
    9521198    }
    … …  
    9611207    (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
    9621208     G_random_element;
     1209
     1210    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
     1211                        *coF_random_element;
     1212    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
     1213                        *coG_random_element;
     1214
    9631215    if (!G_random_element.inCoeffDomain())
    9641216      d0= totaldegree (G_random_element, Variable(2),
    … …  
    9731225      G_m= 0;
    9741226      d= d0;
    975       count= 0;
     1227      coF_m= 0;
     1228      coG_m= 0;
    9761229    }
    9771230
    9781231    TIMING_START (newton_interpolation);
    9791232    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
    980     TIMING_END_AND_PRINT (newton_interpolation, "time for newton interpolation: ");
    981 
    982     count++;
     1233    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
     1234    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
     1235    TIMING_END_AND_PRINT (newton_interpolation,
     1236                          "time for newton interpolation: ");
     1237
    9831238    //termination test
    984     if (uni_lcoeff (H) == gcdlcAlcB)
    985     {
    986       cH= uni_content (H);
     1239    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
     1240    {
     1241      if (gcdlcAlcB.isOne())
     1242        cH= 1;
     1243      else
     1244        cH= uni_content (H);
    9871245      ppH= H/cH;
     1246      CanonicalForm lcppH= gcdlcAlcB/cH;
     1247      CanonicalForm ccoF= lcA/lcppH;
     1248      ccoF /= Lc (ccoF);
     1249      CanonicalForm ccoG= lcB/lcppH;
     1250      ccoG /= Lc (ccoG);
     1251      ppCoF= coF/ccoF;
     1252      ppCoG= coG/ccoG;
    9881253      if (inextension)
    9891254      {
    990         if ((count == bound) || (fdivides(ppH, ppA) && fdivides(ppH, ppB)))
     1255        if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
     1256             (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
     1257             (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) ||
     1258             (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
    9911259        {
    9921260          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
    9931261          ppH= GFMapDown (ppH, k);
     1262          ppCoF= GFMapDown (ppCoF, k);
     1263          ppCoG= GFMapDown (ppCoG, k);
    9941264          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
    9951265          if (compressConvexDense)
     1266          {
    9961267            ppH= decompress (ppH, MM, V);
     1268            ppCoF= decompress (ppCoF, MM, V);
     1269            ppCoG= decompress (ppCoG, MM, V);
     1270          }
     1271          coF= N ((cA/gcdcAcB)*ppCoF);
     1272          coG= N ((cB/gcdcAcB)*ppCoG);
    9971273          setCharacteristic (p, k, gf_name_buf);
    9981274          return N(gcdcAcB*ppH);
    … …  
    10011277      else
    10021278      {
    1003         if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
     1279      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
     1280           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
     1281           (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) ||
     1282           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
    10041283        {
    10051284          if (compressConvexDense)
     1285          {
    10061286            ppH= decompress (ppH, MM, V);
     1287            ppCoF= decompress (ppCoF, MM, V);
     1288            ppCoG= decompress (ppCoG, MM, V);
     1289          }
     1290          coF= N ((cA/gcdcAcB)*ppCoF);
     1291          coG= N ((cB/gcdcAcB)*ppCoG);
    10071292          return N(gcdcAcB*ppH);
    10081293        }
    … …  
    10111296
    10121297    G_m= H;
     1298    coF_m= coF;
     1299    coG_m= coG;
    10131300    newtonPoly= newtonPoly*(x - random_element);
    10141301    m= m*(x - random_element);
    … …  
    10681355}
    10691356
    1070 CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
    1071                            bool& topLevel, CFList& l)
     1357CanonicalForm
     1358GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
     1359             CanonicalForm& coF, CanonicalForm& coG,
     1360             bool& topLevel, CFList& l);
     1361
     1362CanonicalForm
     1363GCD_small_p (const CanonicalForm& F, const CanonicalForm& G,
     1364             bool& topLevel, CFList& l)
     1365{
     1366  CanonicalForm dummy1, dummy2;
     1367  CanonicalForm result= GCD_small_p (F, G, dummy1, dummy2, topLevel, l);
     1368  return result;
     1369}
     1370
     1371CanonicalForm
     1372GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
     1373             CanonicalForm& coF, CanonicalForm& coG,
     1374             bool& topLevel, CFList& l)
    10721375{
    10731376  CanonicalForm A= F;
    10741377  CanonicalForm B= G;
    1075   if (F.isZero() && degree(G) > 0) return G/Lc(G);
    1076   else if (G.isZero() && degree (F) > 0) return F/Lc(F);
    1077   else if (F.isZero() && G.isZero()) return F.genOne();
    1078   if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
    1079   if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
    1080   if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
    1081   if (F == G) return F/Lc(F);
    1082 
     1378  if (F.isZero() && degree(G) > 0)
     1379  {
     1380    coF= 0;
     1381    coG= Lc (G);
     1382    return G/Lc(G);
     1383  }
     1384  else if (G.isZero() && degree (F) > 0)
     1385  {
     1386    coF= Lc (F);
     1387    coG= 0;
     1388    return F/Lc(F);
     1389  }
     1390  else if (F.isZero() && G.isZero())
     1391  {
     1392    coF= coG= 0;
     1393    return F.genOne();
     1394  }
     1395  if (F.inBaseDomain() || G.inBaseDomain())
     1396  {
     1397    coF= F;
     1398    coG= G;
     1399    return F.genOne();
     1400  }
     1401  if (F.isUnivariate() && fdivides(F, G, coG))
     1402  {
     1403    coF= Lc (F);
     1404    return F/Lc(F);
     1405  }
     1406  if (G.isUnivariate() && fdivides(G, F, coF))
     1407  {
     1408    coG= Lc (G);
     1409    return G/Lc(G);
     1410  }
     1411  if (F == G)
     1412  {
     1413    coF= coG= Lc (F);
     1414    return F/Lc(F);
     1415  }
    10831416  CFMap M,N;
    10841417  int best_level= myCompress (A, B, M, N, topLevel);
    10851418
    1086   if (best_level == 0) return B.genOne();
     1419  if (best_level == 0)
     1420  {
     1421    coF= F;
     1422    coG= G;
     1423    return B.genOne();
     1424  }
    10871425
    10881426  A= M(A);
    … …  
    10931431  //univariate case
    10941432  if (A.isUnivariate() && B.isUnivariate())
    1095     return N (gcd (A, B));
     1433  {
     1434    CanonicalForm result= gcd (A, B);
     1435    coF= N (A/result);
     1436    coG= N (B/result);
     1437    return N (result);
     1438  }
    10961439
    10971440  CanonicalForm cA, cB;    // content of A and B
    … …  
    11221465  if (compressConvexDense)
    11231466  {
     1467    CanonicalForm bufcA= cA;
     1468    CanonicalForm bufcB= cB;
    11241469    cA= content (ppA, 1);
    11251470    cB= content (ppB, 1);
    … …  
    11271472    ppB /= cB;
    11281473    gcdcAcB *= gcd (cA, cB);
     1474    cA *= bufcA;
     1475    cB *= bufcB;
    11291476    if (ppA.isUnivariate() || ppB.isUnivariate())
    11301477    {
    11311478      if (ppA.level() == ppB.level())
    1132         return N (gcd (ppA, ppB)*gcdcAcB);
     1479      {
     1480        CanonicalForm result= gcd (ppA, ppB);
     1481        coF= N ((ppA/result)*(cA/gcdcAcB));
     1482        coG= N ((ppB/result)*(cB/gcdcAcB));
     1483        return N (result*gcdcAcB);
     1484      }
    11331485      else
     1486      {
     1487        coF= N (ppA*(cA/gcdcAcB));
     1488        coG= N (ppB*(cB/gcdcAcB));
    11341489        return N (gcdcAcB);
     1490      }
    11351491    }
    11361492
    … …  
    11491505    {
    11501506      if (ppA.level() == ppB.level())
    1151         return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
     1507      {
     1508        CanonicalForm result= gcd (ppA, ppB);
     1509        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
     1510        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
     1511        return N (decompress (result, MM, V)*gcdcAcB);
     1512      }
    11521513      else
     1514      {
     1515        coF= N (decompress (ppA, MM, V));
     1516        coG= N (decompress (ppB, MM, V));
    11531517        return N (gcdcAcB);
     1518      }
    11541519    }
    11551520  }
    … …  
    11601525  lcB= uni_lcoeff (ppB);
    11611526
    1162   if (fdivides (lcA, lcB))
     1527  /*if (fdivides (lcA, lcB))
    11631528  {
    11641529    if (fdivides (A, B))
    … …  
    11691534    if (fdivides (B, A))
    11701535      return G/Lc(G);
    1171   }
     1536  }*/
    11721537
    11731538  gcdlcAlcB= gcd (lcA, lcB);
    … …  
    11771542
    11781543  if (d == 0)
     1544  {
     1545    coF= N (ppA*(cA/gcdcAcB));
     1546    coG= N (ppB*(cB/gcdcAcB));
    11791547    return N(gcdcAcB);
     1548  }
    11801549
    11811550  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
    … …  
    11851554
    11861555  if (d == 0)
     1556  {
     1557    coF= N (ppA*(cA/gcdcAcB));
     1558    coG= N (ppB*(cB/gcdcAcB));
    11871559    return N(gcdcAcB);
     1560  }
    11881561
    11891562  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
    1190   CanonicalForm newtonPoly= 1;
     1563  CanonicalForm newtonPoly, coF_random_element, coG_random_element,
     1564                coF_m, coG_m, ppCoF, ppCoG;
     1565
     1566  newtonPoly= 1;
    11911567  m= gcdlcAlcB;
    11921568  H= 0;
     1569  coF= 0;
     1570  coG= 0;
    11931571  G_m= 0;
    11941572  Variable alpha, V_buf;
    … …  
    11961574  bool inextension= false;
    11971575  bool inextensionextension= false;
    1198   topLevel= false;
     1576  bool topLevel2= false;
    11991577  CFList source, dest;
    1200   int bound= tmin (degree (ppA, 1), degree (ppB, 1));
    1201   int count= 0;
     1578  int bound1= degree (ppA, 1);
     1579  int bound2= degree (ppB, 1);
    12021580  do
    12031581  {
    12041582    if (inextension)
    1205       random_element= randomElement (m, V_buf, l, fail);
     1583      random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
    12061584    else
    1207       random_element= FpRandomElement (m, l, fail);
     1585      random_element= FpRandomElement (m*lcA*lcB, l, fail);
    12081586
    12091587    if (!fail && !inextension)
    … …  
    12121590      TIMING_START (gcd_recursion);
    12131591      G_random_element=
    1214       GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel,
    1215       list);
     1592      GCD_small_p (ppA (random_element,x), ppB (random_element,x),
     1593                   coF_random_element, coG_random_element, topLevel2,
     1594                   list);
    12161595      TIMING_END_AND_PRINT (gcd_recursion,
    12171596                            "time for recursive call: ");
    … …  
    12231602      TIMING_START (gcd_recursion);
    12241603      G_random_element=
    1225       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
    1226                         list, topLevel);
     1604      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
     1605                        coF_random_element, coG_random_element, alpha,
     1606                        list, topLevel2);
    12271607      TIMING_END_AND_PRINT (gcd_recursion,
    12281608                            "time for recursive call: ");
    … …  
    12411621        inextension= true;
    12421622        fail= false;
    1243         random_element= randomElement (m, alpha, l, fail);
     1623        random_element= randomElement (m*lcA*lcB, alpha, l, fail);
    12441624        deg++;
    12451625      } while (fail);
    … …  
    12481628      TIMING_START (gcd_recursion);
    12491629      G_random_element=
    1250       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
    1251                         list, topLevel);
     1630      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
     1631                        coF_random_element, coG_random_element, alpha,
     1632                        list, topLevel2);
    12521633      TIMING_END_AND_PRINT (gcd_recursion,
    12531634                            "time for recursive call: ");
    … …  
    12691650      {
    12701651        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
    1271         G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
     1652        G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
     1653        coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
     1654        coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
    12721655        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
    12731656                             source, dest);
    … …  
    12761659        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
    12771660                            dest);
     1661        lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
     1662        lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
    12781663        for (CFListIterator i= l; i.hasItem(); i++)
    12791664          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
    … …  
    12961681      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
    12971682      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     1683      coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     1684      coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
    12981685      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
    12991686                          source, dest);
    … …  
    13021689      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
    13031690                         source, dest);
     1691      lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     1692      lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
    13041693      fail= false;
    1305       random_element= randomElement (m, V_buf, l, fail );
     1694      random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
    13061695      DEBOUTLN (cerr, "fail= " << fail);
    13071696      CFList list;
    13081697      TIMING_START (gcd_recursion);
    13091698      G_random_element=
    1310       GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
    1311                         list, topLevel);
     1699      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
     1700                        coF_random_element, coG_random_element, V_buf,
     1701                        list, topLevel2);
    13121702      TIMING_END_AND_PRINT (gcd_recursion,
    13131703                            "time for recursive call: ");
    … …  
    13221712
    13231713    if (d0 == 0)
     1714    {
     1715      coF= N (ppA*(cA/gcdcAcB));
     1716      coG= N (ppB*(cB/gcdcAcB));
    13241717      return N(gcdcAcB);
     1718    }
     1719
    13251720    if (d0 >  d)
    13261721    {
    … …  
    13331728                       *G_random_element;
    13341729
     1730    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
     1731                        *coF_random_element;
     1732    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
     1733                        *coG_random_element;
    13351734
    13361735    if (!G_random_element.inCoeffDomain())
    … …  
    13461745      G_m= 0;
    13471746      d= d0;
    1348       count= 0;
     1747      coF_m= 0;
     1748      coG_m= 0;
    13491749    }
    13501750
    13511751    TIMING_START (newton_interpolation);
    13521752    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
     1753    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
     1754    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
    13531755    TIMING_END_AND_PRINT (newton_interpolation,
    13541756                          "time for newton_interpolation: ");
    13551757
    1356     count++;
    13571758    //termination test
    1358     if (uni_lcoeff (H) == gcdlcAlcB)
    1359     {
    1360       cH= uni_content (H);
     1759    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
     1760    {
     1761      if (gcdlcAlcB.isOne())
     1762        cH= 1;
     1763      else
     1764        cH= uni_content (H);
    13611765      ppH= H/cH;
    13621766      ppH /= Lc (ppH);
     1767      CanonicalForm lcppH= gcdlcAlcB/cH;
     1768      CanonicalForm ccoF= lcppH/Lc (lcppH);
     1769      CanonicalForm ccoG= lcppH/Lc (lcppH);
     1770      ppCoF= coF/ccoF;
     1771      ppCoG= coG/ccoG;
    13631772      DEBOUTLN (cerr, "ppH= " << ppH);
    1364       if ((count == bound) || (fdivides (ppH, ppA) && fdivides (ppH, ppB)))
     1773      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
     1774           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
     1775           (ppA == ppCoF*ppH && ppB == ppCoG*ppH)) ||
     1776           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
    13651777      {
    13661778        if (compressConvexDense)
     1779        {
    13671780          ppH= decompress (ppH, MM, V);
     1781          ppCoF= decompress (ppCoF, MM, V);
     1782          ppCoG= decompress (ppCoG, MM, V);
     1783        }
     1784        coF= N ((cA/gcdcAcB)*ppCoF);
     1785        coG= N ((cB/gcdcAcB)*ppCoG);
    13681786        return N(gcdcAcB*ppH);
    13691787      }
    … …  
    13711789
    13721790    G_m= H;
     1791    coF_m= coF;
     1792    coG_m= coG;
    13731793    newtonPoly= newtonPoly*(x - random_element);
    13741794    m= m*(x - random_element);
    … …  
    44524872      {
    44534873        cand = DD[2] / content( DD[2], Variable(1) );
    4454         gcdfound = fdivides( cand, G ) && fdivides ( cand, F );
     4874        gcdfound = fdivides( cand, G ) && fdivides ( cand, F ); //TODO use multiplication instead of fdivides
    44554875
    44564876        if (passToGF && gcdfound)
Note: See TracChangeset for help on using the changeset viewer.