Changeset f7a4e9 in git


Ignore:
Timestamp:
Apr 14, 2012, 1:39:51 PM (12 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', '6e5adcba05493683b94648c659a729c189812c77')
Children:
c244d04150a9fd45331938d2c9fb6f89985acb05
Parents:
a9a6dcbcbc4f0be99516f14eae42b1776235ca9f
git-author:
Martin Lee <martinlee84@web.de>2012-04-14 13:39:51+02:00
git-committer:
Martin Lee <martinlee84@web.de>2012-05-07 14:14:17+02:00
Message:
chg: towards better EZGCD
Location:
factory
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_gcd.cc

    ra9a6dcb rf7a4e9  
    757757  else if (!fc_and_gc_Univariate)
    758758  {
    759     if (
    760     isOn(SW_USE_CHINREM_GCD)
    761     && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
    762     )
    763     {
    764     #if 0
    765       if ( p1 == fc.level() )
    766         fc = chinrem_gcd( fc, gc );
    767       else
    768       {
    769         fc = replacevar( fc, Variable(p1), Variable(mp) );
    770         gc = replacevar( gc, Variable(p1), Variable(mp) );
    771         fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
    772       }
    773     #else
    774       fc = chinrem_gcd( fc, gc);
    775     #endif
    776     }
    777     else if ( isOn( SW_USE_EZGCD ) )
    778     {
    779       if ( pe == 1 )
     759    if ( isOn( SW_USE_EZGCD ) )
     760    {
     761      fc= ezgcd (fc, gc);
     762      /*if ( pe == 1 )
    780763        fc = ezgcd( fc, gc );
    781764      else if ( pe > 0 )// no variable at position 1
    … …  
    791774        gc = swapvar( gc, Variable(pe), Variable(1) );
    792775        fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
     776      }*/
     777    }
     778    else if (
     779    isOn(SW_USE_CHINREM_GCD)
     780    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
     781    )
     782    {
     783    #if 0
     784      if ( p1 == fc.level() )
     785        fc = chinrem_gcd( fc, gc );
     786      else
     787      {
     788        fc = replacevar( fc, Variable(p1), Variable(mp) );
     789        gc = replacevar( gc, Variable(p1), Variable(mp) );
     790        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
    793791      }
     792    #else
     793      fc = chinrem_gcd( fc, gc);
     794    #endif
    794795    }
    795796    else

  • factory/fac_ezgcd.cc

    ra9a6dcb rf7a4e9  
    1616#include "fac_distrib.h"
    1717#include "templates/ftmpl_functions.h"
     18#include "cf_map.h"
     19#include "facHensel.h"
     20
     21static
     22int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
     23                    CFMap & N, int& both_non_zero)
     24{
     25  int n= tmax (F.level(), G.level());
     26  int * degsf= new int [n + 1];
     27  int * degsg= new int [n + 1];
     28
     29  for (int i = 0; i <= n; i++)
     30    degsf[i]= degsg[i]= 0;
     31
     32  degsf= degrees (F, degsf);
     33  degsg= degrees (G, degsg);
     34
     35  both_non_zero= 0;
     36  int f_zero= 0;
     37  int g_zero= 0;
     38
     39  for (int i= 1; i <= n; i++)
     40  {
     41    if (degsf[i] != 0 && degsg[i] != 0)
     42    {
     43      both_non_zero++;
     44      continue;
     45    }
     46    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
     47    {
     48      f_zero++;
     49      continue;
     50    }
     51    if (degsg[i] == 0 && degsf[i] && i <= F.level())
     52    {
     53      g_zero++;
     54      continue;
     55    }
     56  }
     57
     58  if (both_non_zero == 0)
     59  {
     60    delete [] degsf;
     61    delete [] degsg;
     62    return 0;
     63  }
     64
     65  // map Variables which do not occur in both polynomials to higher levels
     66  int k= 1;
     67  int l= 1;
     68  for (int i= 1; i <= n; i++)
     69  {
     70    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
     71    {
     72      if (k + both_non_zero != i)
     73      {
     74        M.newpair (Variable (i), Variable (k + both_non_zero));
     75        N.newpair (Variable (k + both_non_zero), Variable (i));
     76      }
     77      k++;
     78    }
     79    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
     80    {
     81      if (l + g_zero + both_non_zero != i)
     82      {
     83        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
     84        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
     85      }
     86      l++;
     87    }
     88  }
     89
     90  // sort Variables x_{i} in decreasing order of
     91  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
     92  int m= tmin (F.level(), G.level());
     93  int max_min_deg;
     94  k= both_non_zero;
     95  l= 0;
     96  int i= 1;
     97  while (k > 0)
     98  {
     99    max_min_deg= tmin (degsf[i], degsg[i]);
     100    while (max_min_deg == 0)
     101    {
     102      i++;
     103      max_min_deg= tmin (degsf[i], degsg[i]);
     104    }
     105    for (int j= i + 1; j <= m; j++)
     106    {
     107      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
     108          (tmin (degsf[j], degsg[j]) != 0))
     109      {
     110        max_min_deg= tmin (degsf[j], degsg[j]);
     111        l= j;
     112      }
     113    }
     114
     115    if (l != 0)
     116    {
     117      if (l != k)
     118      {
     119        M.newpair (Variable (l), Variable(k));
     120        N.newpair (Variable (k), Variable(l));
     121        degsf[l]= 0;
     122        degsg[l]= 0;
     123        l= 0;
     124      }
     125      else
     126      {
     127        degsf[l]= 0;
     128        degsg[l]= 0;
     129        l= 0;
     130      }
     131    }
     132    else if (l == 0)
     133    {
     134      if (i != k)
     135      {
     136        M.newpair (Variable (i), Variable (k));
     137        N.newpair (Variable (k), Variable (i));
     138        degsf[i]= 0;
     139        degsg[i]= 0;
     140      }
     141      else
     142      {
     143        degsf[i]= 0;
     144        degsg[i]= 0;
     145      }
     146      i++;
     147    }
     148    k--;
     149  }
     150
     151  delete [] degsf;
     152  delete [] degsg;
     153
     154  return both_non_zero;
     155}
     156
     157static inline
     158CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
     159                            const CFList& evaluation)
     160{
     161  CanonicalForm A= F;
     162  int k= 2;
     163  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
     164    A= A (Variable (k) + i.getItem(), k);
     165
     166  CanonicalForm buf= A;
     167  Feval= CFList();
     168  Feval.append (buf);
     169  for (k= evaluation.length() + 1; k > 2; k--)
     170  {
     171    buf= mod (buf, Variable (k));
     172    Feval.insert (buf);
     173  }
     174  return A;
     175}
     176
     177static inline
     178CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
     179{
     180  int l= evaluation.length() + 1;
     181  CanonicalForm result= F;
     182  CFListIterator j= evaluation;
     183  for (int i= 2; i < l + 1; i++, j++)
     184  {
     185    if (F.level() < i)
     186      continue;
     187    result= result (Variable (i) - j.getItem(), i);
     188  }
     189  return result;
     190}
     191
     192static inline
     193Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
     194                    CFMap & N, const Evaluation& A)
     195{
     196  int n= F.level();
     197  int * degsf= new int [n + 1];
     198
     199  for (int i = 0; i <= n; i++)
     200    degsf[i]= 0;
     201
     202  degsf= degrees (F, degsf);
     203
     204  Evaluation result= Evaluation (A.min(), A.max());
     205  ASSERT (A.min() == 2, "expected A.min() == 2");
     206  ASSERT (A.max() == n, "expected A.max() == n");
     207  int max_deg;
     208  int k= n;
     209  int l= 1;
     210  int i= 2;
     211  int pos= 2;
     212  while (k > 1)
     213  {
     214    max_deg= degsf [i];
     215    while (max_deg == 0)
     216    {
     217      i++;
     218      max_deg= degsf [i];
     219    }
     220    l= i;
     221    for (int j= i + 1; j <=  n; j++)
     222    {
     223      if (degsf[j] > max_deg)
     224      {
     225        max_deg= degsf[j];
     226        l= j;
     227      }
     228    }
     229
     230    if (l <= n)
     231    {
     232      if (l != pos)
     233      {
     234        result.setValue (pos, A [l]);
     235        M.newpair (Variable (l), Variable (pos));
     236        N.newpair (Variable (pos), Variable (l));
     237        degsf[l]= 0;
     238        l= 2;
     239        if (k == 2 && n == 3)
     240        {
     241          result.setValue (l, A [pos]);
     242          M.newpair (Variable (pos), Variable (l));
     243          N.newpair (Variable (l), Variable (pos));
     244          degsf[pos]= 0;
     245        }
     246      }
     247      else
     248      {
     249        result.setValue (l, A [l]);
     250        degsf [l]= 0;
     251      }
     252    }
     253    pos++;
     254    k--;
     255    l= 2;
     256  }
     257
     258  delete [] degsf;
     259
     260  return result;
     261}
     262
     263static inline
     264int Hensel (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
     265            const CFArray& LeadCoeffs )
     266{
     267  CFList factors;
     268  factors.append (G[1]);
     269  factors.append (G[2]);
     270
     271  CFMap NN, MM;
     272  Evaluation A= optimize4Lift (UU, MM, NN, AA);
     273
     274  CanonicalForm U= MM (UU);
     275  CFArray LCs= CFArray (1,2);
     276  LCs [1]= MM (LeadCoeffs [1]);
     277  LCs [2]= MM (LeadCoeffs [2]);
     278
     279  CFList evaluation;
     280  for (int i= A.min(); i <= A.max(); i++)
     281    evaluation.append (A [i]);
     282  CFList UEval;
     283  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
     284
     285  if (size (shiftedU)/getNumVars (U) > 500)
     286    return -1;
     287
     288  CFArray shiftedLCs= CFArray (2);
     289  CFList shiftedLCsEval1, shiftedLCsEval2;
     290  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
     291  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
     292  factors.insert (1);
     293  int liftBound= degree (UEval.getLast(), 2) + 1;
     294  CFArray Pi;
     295  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
     296  CFList diophant;
     297  CFArray lcs= CFArray (2);
     298  lcs [0]= shiftedLCsEval1.getFirst();
     299  lcs [1]= shiftedLCsEval2.getFirst();
     300  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
     301                        lcs, false);
     302
     303  for (CFListIterator i= factors; i.hasItem(); i++)
     304  {
     305    if (!fdivides (i.getItem(), UEval.getFirst()))
     306      return 0;
     307  }
     308
     309  int * liftBounds;
     310  bool noOneToOne= false;
     311  if (U.level() > 2)
     312  {
     313    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
     314    liftBounds[0]= liftBound;
     315    for (int i= 1; i < U.level() - 1; i++)
     316      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
     317    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
     318                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
     319                                  diophant, noOneToOne);
     320    delete [] liftBounds;
     321    if (noOneToOne)
     322      return 0;
     323  }
     324  G[1]= factors.getFirst();
     325  G[2]= factors.getLast();
     326  G[1]= myReverseShift (G[1], evaluation);
     327  G[2]= myReverseShift (G[2], evaluation);
     328  G[1]= NN (G[1]);
     329  G[2]= NN (G[2]);
     330  return 1;
     331}
    18332
    19333static void findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG );
    … …  
    37351ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal )
    38352{
    39     CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD;
     353    CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD, cand, result;
    40354    CFArray DD( 1, 2 ), lcDD( 1, 2 );
    41355    int degF, degG, delta, t;
    … …  
    45359    modpk bound;
    46360
     361    F= FF;
     362    G= GG;
     363
     364    CFMap M,N;
     365    int smallestDegLev;
     366    int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
     367
     368    if (best_level == 0) return G.genOne();
     369
     370    F= M (F);
     371    G= M (G);
     372
     373
    47374    DEBINCLEVEL( cerr, "ezgcd" );
    48375    DEBOUTLN( cerr, "FF = " << FF );
    49376    DEBOUTLN( cerr, "GG = " << GG );
    50     f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
     377    f = content( F, x ); g = content( G, x ); d = gcd( f, g );
    51378    DEBOUTLN( cerr, "f = " << f );
    52379    DEBOUTLN( cerr, "g = " << g );
    53     F = FF / f; G = GG / g;
     380    F /= f; G /= g;
    54381    if ( F.isUnivariate() && G.isUnivariate() )
    55382    {
    … …  
    57384        if(F.mvar()==G.mvar())
    58385          d*=gcd(F,G);
    59         return d;
     386        return N (d);
    60387    }
    61388    else  if ( gcd_test_one( F, G, false ) )
    62389    {
    63390        DEBDECLEVEL( cerr, "ezgcd" );
    64         return d;
     391        return N (d);
    65392    }
    66393    lcF = LC( F, x ); lcG = LC( G, x );
    … …  
    69396    degF = degree( F, x ); degG = degree( G, x );
    70397    t = tmax( F.level(), G.level() );
    71     bound = findBound( F, G, lcF, lcG, degF, degG );
     398    //bound = findBound( F, G, lcF, lcG, degF, degG );
    72399    if ( ! internal )
    73400        b = REvaluation( 2, t, IntRandom( 25 ) );
    74     while ( ! gcdfound ) {
     401    while ( ! gcdfound )
     402    {
    75403        /// ---> A2
    76404        DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
    … …  
    87415        {
    88416          DEBDECLEVEL( cerr, "ezgcd" );
    89           return d;
     417          return N (d);
    90418        }
    91419        /// ---> A4
    … …  
    98426            {
    99427                DEBDECLEVEL( cerr, "ezgcd" );
    100                 return d;
     428                return N (d);
    101429            }
    102430            if ( dd /*degree( Dbt )*/ == delta )
    … …  
    110438        DEBOUTLN( cerr, "now after A4, delta = " << delta );
    111439        /// ---> A5
    112         if ( degF <= degG && delta == degF && fdivides( F, G ) )
     440        if (delta == degF)
    113441        {
     442          if (degF <= degG  && fdivides (F, G))
     443          {
    114444            DEBDECLEVEL( cerr, "ezgcd" );
    115             return d*F;
    116         }
    117         if ( degG < degF && delta == degG && fdivides( G, F ) )
     445            return N (d*F);
     446          }
     447          else
     448            delta--;
     449        }
     450        else if (delta == degG)
    118451        {
     452          if (degG <= degF && fdivides( G, F ))
     453          {
    119454            DEBDECLEVEL( cerr, "ezgcd" );
    120             return d*G;
    121         }
    122         if ( delta != degF && delta != degG ) {
     455            return N (d*G);
     456          }
     457          else
     458            delta--;
     459        }
     460        if ( delta != degF && delta != degG )
     461        {
     462            /// ---> A6
    123463            bool B_is_F;
    124             /// ---> A6
    125             CanonicalForm xxx;
     464            CanonicalForm xxx1, xxx2;
     465            CanonicalForm buf;
     466            DD[1] = Fb / Db;
     467            buf= Gb/Db;
     468            xxx1 = gcd( DD[1], Db );
     469            xxx2 = gcd( buf, Db );
     470            if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
     471                (size (F) <= size (G)))
     472                 || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
     473            {
     474              B = F;
     475              DD[2] = Db;
     476              lcDD[1] = lcF;
     477              lcDD[2] = lcD;
     478              B_is_F = true;
     479            }
     480            else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
     481                     (size (G) < size (F)))
     482                     || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
     483            {
     484              DD[1] = buf;
     485              B = G;
     486              DD[2] = Db;
     487              lcDD[1] = lcG;
     488              lcDD[2] = lcD;
     489              B_is_F = false;
     490            }
     491            else
     492            {
     493              Off(SW_USE_EZGCD);
     494              result=gcd( F, G );
     495              On(SW_USE_EZGCD);
     496              return N (d*result);
     497            }
     498            /*CanonicalForm xxx;
    126499            //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
    127500            DD[1] = Fb / Db;
    … …  
    165538#endif
    166539              }
    167             }
     540            }*/
    168541            /// ---> A7
    169542            DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
    170543            DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
    171             bound = enlargeBound( B, DD[2], DD[1], bound );
     544            //bound = enlargeBound( B, DD[2], DD[1], bound );
    172545            DEBOUTLN( cerr, "(hensel) B    = " << B );
    173546            DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
    … …  
    175548            DEBOUTLN( cerr, "(hensel) DD   = " << DD );
    176549            DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
    177             gcdfound = Hensel( B, DD, lcDD, b, bound, x );
     550            gcdfound= Hensel (B*lcD, DD, b, lcDD);
    178551            DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
     552
     553            if (gcdfound == -1)
     554            {
     555              Off (SW_USE_EZGCD);
     556              result= gcd (F,G);
     557              On (SW_USE_EZGCD_P);
     558              return N (d*result);
     559            }
    179560
    180561            if (gcdfound)
    181562            {
    182               CanonicalForm xxx=content(DD[2],Variable(1));
    183               CanonicalForm cand=DD[2] / xxx; //content(DD[2],Variable(1));
    184 #if 0
    185               gcdfound= fdivides(cand,G) &&  fdivides(cand,F);
    186 #else
    187               if (B_is_F /*B==F*lcF*/)
    188               {
    189                 DEBOUTLN( cerr, "(test) G: "<<G<<" % gcd:"<<cand<<" -> " << G%cand );
    190                 gcdfound= fdivides(cand,G);
    191               }
    192               else
    193               {
    194                 DEBOUTLN( cerr, "(test) F: "<<F<<" % gcd:"<<cand<<" -> " << F%cand);
    195                 gcdfound= fdivides(cand,F);
    196               }
    197 #endif
     563              cand = DD[2] / content( DD[2], Variable(1) );
     564              gcdfound = fdivides( cand, G ) && fdivides ( cand, F );
    198565            }
    199566            /// ---> A8 (gcdfound)
    200567        }
    201         delta++;
     568        delta--;
    202569    }
    203570    /// ---> A9
    204571    DEBDECLEVEL( cerr, "ezgcd" );
    205     return d * DD[2] / content(DD[2],Variable(1));
     572    return N (d*cand);
    206573}
    207574
Note: See TracChangeset for help on using the changeset viewer.