Changeset e368746 in git


Ignore:
Timestamp:
Jun 6, 2011, 3:33:48 PM (12 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
Children:
c227a2fe328d29d1cfbd56bc8d8824aed8aa2d38
Parents:
1bebbfe6863848ca47716f2e6bab1742430df913
Message:
added new functions for multivariate Hensel lifting
minor changes of existing Hensel lift functions


git-svn-id: file:///usr/local/Singular/svn/trunk@14259 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
factory
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_gcd_smallp.cc

    r1bebbf re368746  
    38623862    return !noChange;
    38633863  int * liftBounds;
     3864  noChange= false;
    38643865  if (U.level() > 2)
    38653866  {
    … …  
    38693870      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
    38703871    factors= henselLift2 (UEval, factors, liftBounds, U.level() - 1, false,
    3871                           shiftedLCsEval1, shiftedLCsEval2, Pi, diophant);
     3872                          shiftedLCsEval1, shiftedLCsEval2, Pi, diophant, noChange);
    38723873   delete [] liftBounds;
     3874   if (noChange)
     3875     return !noChange;
    38733876  }
    38743877  G[1]= factors.getFirst();

  • factory/facHensel.cc

    r1bebbf re368746  
    11861186  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
    11871187  CFList splitA= split (A, m, x);
    1188   CFList splitB= split (B, m, x);
    11891188  if (splitA.length() == 3)
    11901189    splitA.insert (0);
    … …  
    13091308  CanonicalForm A= mod (F, M);
    13101309  CanonicalForm B= mod (G, M);
     1310
     1311  if (B.inCoeffDomain())
     1312  {
     1313    divrem (A, B, Q, R);
     1314    return;
     1315  }
     1316  if (A.inCoeffDomain() && !B.inCoeffDomain())
     1317  {
     1318    Q= 0;
     1319    R= A;
     1320    return;
     1321  }
     1322
     1323  if (B.level() < A.level())
     1324  {
     1325    divrem (A, B, Q, R);
     1326    return;
     1327  }
     1328  if (A.level() > B.level())
     1329  {
     1330    R= A;
     1331    Q= 0;
     1332    return;
     1333  }
     1334  if (B.level() == 1 && B.isUnivariate())
     1335  {
     1336    divrem (A, B, Q, R);
     1337    return;
     1338  }
     1339  if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
     1340  {
     1341    Q= 0;
     1342    R= A;
     1343    return;
     1344  }
     1345
    13111346  Variable x= Variable (1);
    13121347  int degB= degree (B, x);
    … …  
    15291564      if (k != j - k + 1)
    15301565      {
    1531         if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     1566        if ((one.hasTerms() && one.exp() == j - k + 1) && (two.hasTerms() && two.exp() == j - k + 1))
    15321567        {
    15331568          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
    … …  
    15361571          two++;
    15371572        }
    1538         else if (one.exp() == j - k + 1)
     1573        else if (one.hasTerms() && one.exp() == j - k + 1)
    15391574        {
    15401575          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) -
    … …  
    15421577          one++;
    15431578        }
    1544         else if (two.exp() == j - k + 1)
     1579        else if (two.hasTerms() && two.exp() == j - k + 1)
    15451580        {
    15461581          tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) -
    … …  
    15941629    one= bufFactors [l + 1];
    15951630    two= Pi [l - 1];
    1596     if (two.exp() == j + 1)
     1631    if (two.hasTerms() && two.exp() == j + 1)
    15971632    {
    15981633      if (degBuf > 0 && degPi > 0)
    … …  
    16131648        if (k != j - k + 1)
    16141649        {
    1615           if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     1650          if ((one.hasTerms() && one.exp() == j - k + 1) &&
     1651              (two.hasTerms() && two.exp() == j - k + 1))
    16161652          {
    16171653            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
    … …  
    16201656            two++;
    16211657          }
    1622           else if (one.exp() == j - k + 1)
     1658          else if (one.hasTerms() && one.exp() == j - k + 1)
    16231659          {
    16241660            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) -
    … …  
    16261662            one++;
    16271663          }
    1628           else if (two.exp() == j - k + 1)
     1664          else if (two.hasTerms() && two.exp() == j - k + 1)
    16291665          {
    16301666            tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) -
    … …  
    16441680void
    16451681henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
    1646               CFList& diophant, CFMatrix& M)
    1647 {
    1648   sortList (factors, Variable (1));
     1682              CFList& diophant, CFMatrix& M, bool sort)
     1683{
     1684  if (sort)
     1685    sortList (factors, Variable (1));
    16491686  Pi= CFArray (factors.length() - 1);
    16501687  CFListIterator j= factors;
    … …  
    20182055      if (k != j - k + 1)
    20192056      {
    2020         if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     2057        if ((one.hasTerms() && one.exp() == j - k + 1) &&
     2058            (two.hasTerms() && two.exp() == j - k + 1))
    20212059        {
    20222060          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
    … …  
    20262064          two++;
    20272065        }
    2028         else if (one.exp() == j - k + 1)
     2066        else if (one.hasTerms() && one.exp() == j - k + 1)
    20292067        {
    20302068          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
    … …  
    20322070          one++;
    20332071        }
    2034         else if (two.exp() == j - k + 1)
     2072        else if (two.hasTerms() && two.exp() == j - k + 1)
    20352073        {
    20362074          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
    … …  
    20842122    one= bufFactors [l + 1];
    20852123    two= Pi [l - 1];
    2086     if (two.exp() == j + 1)
     2124    if (two.hasTerms() && two.exp() == j + 1)
    20872125    {
    20882126      if (degBuf > 0 && degPi > 0)
    … …  
    21032141        if (k != j - k + 1)
    21042142        {
    2105           if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     2143          if ((one.hasTerms() && one.exp() == j - k + 1) &&
     2144              (two.hasTerms() && two.exp() == j - k + 1))
    21062145          {
    21072146            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
    … …  
    21112150            two++;
    21122151          }
    2113           else if (one.exp() == j - k + 1)
     2152          else if (one.hasTerms() && one.exp() == j - k + 1)
    21142153          {
    21152154            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
    … …  
    21172156            one++;
    21182157          }
    2119           else if (two.exp() == j - k + 1)
     2158          else if (two.hasTerms() && two.exp() == j - k + 1)
    21202159          {
    21212160            tmp[l] += mulMod (bufFactors[l + 1] [k],
    … …  
    22452284CFList
    22462285henselLift (const CFList& eval, const CFList& factors, const int* l, const int
    2247             lLength)
     2286            lLength, bool sort)
    22482287{
    22492288  CFList diophant;
    22502289  CFList buf= factors;
    22512290  buf.insert (LC (eval.getFirst(), 1));
    2252   sortList (buf, Variable (1));
     2291  if (sort)
     2292    sortList (buf, Variable (1));
    22532293  CFArray Pi;
    22542294  CFMatrix M= CFMatrix (l[1], factors.length());
    … …  
    23452385    for (k= 1; k <= (int) ceil (j/2.0); k++)
    23462386    {
    2347       if (k != j - k + 1)
     2387      if (one.hasTerms() && two.hasTerms())
    23482388      {
    2349         if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     2389        if (k != j - k + 1)
     2390        {
     2391          if ((one.hasTerms() && one.exp() == j - k + 1) && +
     2392              (two.hasTerms() && two.exp() == j - k + 1))
    23502393        {
    23512394          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
    … …  
    23542397          two++;
    23552398        }
    2356         else if (one.exp() == j - k + 1)
     2399        else if (one.hasTerms() && one.exp() == j - k + 1)
    23572400        {
    23582401          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
    … …  
    23602403          one++;
    23612404        }
    2362         else if (two.exp() == j - k + 1)
     2405        else if (two.hasTerms() && two.exp() == j - k + 1)
    23632406        {
    23642407          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
    … …  
    23692412      else
    23702413        tmp[0] += M (k + 1, 1);
     2414      }
    23712415    }
    23722416  }
    … …  
    23952439
    23962440  Pi [0] += tmp[0]*xToJ*F.mvar();
    2397 
    2398   /*// update Pi [l]
    2399   int degPi, degBuf;
    2400   for (int l= 1; l < factors.length() - 1; l++)
    2401   {
    2402     degPi= degree (Pi [l - 1], x);
    2403     degBuf= degree (bufFactors[l + 1], x);
    2404     if (degPi > 0 && degBuf > 0)
    2405       M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
    2406     if (j == 1)
    2407     {
    2408       if (degPi > 0 && degBuf > 0)
    2409         Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
    2410                   bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
    2411                   M (1, l + 1));
    2412       else if (degPi > 0)
    2413         Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
    2414       else if (degBuf > 0)
    2415         Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
    2416     }
    2417     else
    2418     {
    2419       if (degPi > 0 && degBuf > 0)
    2420       {
    2421         uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
    2422         uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
    2423       }
    2424       else if (degPi > 0)
    2425         uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
    2426       else if (degBuf > 0)
    2427       {
    2428         uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
    2429         uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
    2430       }
    2431       Pi[l] += xToJ*uIZeroJ;
    2432     }
    2433     one= bufFactors [l + 1];
    2434     two= Pi [l - 1];
    2435     if (two.exp() == j + 1)
    2436     {
    2437       if (degBuf > 0 && degPi > 0)
    2438       {
    2439           tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
    2440           two++;
    2441       }
    2442       else if (degPi > 0)
    2443       {
    2444           tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
    2445           two++;
    2446       }
    2447     }
    2448     if (degBuf > 0 && degPi > 0)
    2449     {
    2450       for (k= 1; k <= (int) ceil (j/2.0); k++)
    2451       {
    2452         if (k != j - k + 1)
    2453         {
    2454           if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
    2455           {
    2456             tmp[l] += mulNTL ((bufFactors[l + 1][k]+one.coeff()),(Pi[l-1] [k] +
    2457                         two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
    2458             one++;
    2459             two++;
    2460           }
    2461           else if (one.exp() == j - k + 1)
    2462           {
    2463             tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()),Pi[l - 1] [k]) -
    2464                         M (k + 1, l + 1);
    2465             one++;
    2466           }
    2467           else if (two.exp() == j - k + 1)
    2468           {
    2469             tmp[l] += mulNTL (bufFactors[l+1][k],(Pi[l-1] [k] + two.coeff())) -
    2470                       M (k + 1, l + 1);
    2471             two++;
    2472           }
    2473         }
    2474         else
    2475           tmp[l] += M (k + 1, l + 1);
    2476       }
    2477     }
    2478     Pi[l] += tmp[l]*xToJ*F.mvar();
    2479   }*/
    24802441  return;
    24812442}
    … …  
    25432504CFList
    25442505diophantine (const CFList& recResult, const CFList& factors,
    2545              const CFList& products, const CFList& M, const CanonicalForm& E)
     2506             const CFList& products, const CFList& M, const CanonicalForm& E,
     2507             bool& bad)
    25462508{
    25472509  if (M.isEmpty())
    … …  
    25742536  CanonicalForm bufE= mod (E, y);
    25752537  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
    2576                                       bufE);
     2538                                      bufE, bad);
     2539
     2540  if (bad)
     2541    return CFList();
    25772542
    25782543  CanonicalForm e= E;
    25792544  CFListIterator j= products;
    25802545  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
    2581     e -= mulMod (i.getItem(), j.getItem(), M);
     2546    e -= i.getItem()*j.getItem();
    25822547
    25832548  CFList result= recDiophantine;
    … …  
    25942559      CFListIterator k= result;
    25952560      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
    2596                                    coeffE);
     2561                                   coeffE, bad);
     2562      if (bad)
     2563        return CFList();
    25972564      CFListIterator l= products;
    25982565      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
    25992566      {
    26002567        k.getItem() += j.getItem()*power (y, i);
    2601         e -= mulMod (j.getItem()*power (y, i), l.getItem(), M);
     2568        e -= j.getItem()*power (y, i)*l.getItem();
    26022569      }
    26032570    }
    … …  
    26052572      break;
    26062573  }
     2574  if (!e.isZero())
     2575  {
     2576    bad= true;
     2577    return CFList();
     2578  }
    26072579  return result;
    26082580}
    26092581
    2610 // non monic case
    2611 static inline
    2612 CFList
    2613 multiRecDiophantine2 (const CanonicalForm& F, const CFList& factors,
    2614                      const CFList& recResult, const CFList& M, const int d)
    2615 {
    2616   Variable y= F.mvar();
    2617   CFList result;
    2618   CFListIterator i;
    2619   CanonicalForm e= 1;
    2620   CFListIterator j= factors;
    2621   CFList p;
    2622   CFArray bufFactors= CFArray (factors.length());
    2623   CanonicalForm yToD= power (y, d);
    2624   int k= 0;
    2625   for (CFListIterator i= factors; i.hasItem(); i++, k++)
    2626     bufFactors [k]= i.getItem();
    2627   CanonicalForm b;
    2628   CFList buf= M;
    2629   buf.removeLast();
    2630   buf.append (yToD);
    2631   for (k= 0; k < factors.length(); k++) //TODO compute b's faster
    2632   {
    2633     b= 1;
    2634     if (fdivides (bufFactors[k], F))
    2635       b= F/bufFactors[k];
    2636     else
    2637     {
    2638       for (int l= 0; l < factors.length(); l++)
    2639       {
    2640         if (l == k)
    2641           continue;
    2642         else
    2643         {
    2644           b= mulMod (b, bufFactors[l], buf);
    2645         }
    2646       }
    2647     }
    2648     p.append (b);
    2649   }
    2650   j= p;
    2651   for (CFListIterator i= recResult; i.hasItem(); i++, j++)
    2652     e -= mulMod (i.getItem(), j.getItem(), M);
    2653   if (e.isZero())
    2654     return recResult;
    2655   CanonicalForm coeffE;
    2656   result= recResult;
    2657   CanonicalForm g;
    2658   buf.removeLast();
    2659   for (int i= 1; i < d; i++)
    2660   {
    2661     if (degree (e, y) > 0)
    2662       coeffE= e[i];
    2663     else
    2664       coeffE= 0;
    2665     if (!coeffE.isZero())
    2666     {
    2667       CFListIterator k= result;
    2668       CFListIterator l= p;
    2669       j= recResult;
    2670       int ii= 0;
    2671       CanonicalForm dummy;
    2672       CFList recDiophantine;
    2673       CFList buf2, buf3;
    2674       buf2= factors;
    2675       buf3= p;
    2676       for (CFListIterator iii= buf2; iii.hasItem(); iii++)
    2677         iii.getItem()= mod (iii.getItem(), y);
    2678       for (CFListIterator iii= buf3; iii.hasItem(); iii++)
    2679         iii.getItem()= mod (iii.getItem(), y);
    2680       recDiophantine= diophantine (recResult, buf2, buf3, buf, coeffE);
    2681       CFListIterator iter= recDiophantine;
    2682       for (; j.hasItem(); j++, k++, l++, ii++, iter++)
    2683       {
    2684         k.getItem() += iter.getItem()*power (y, i);
    2685         e -= mulMod (iter.getItem()*power (y, i), l.getItem(), M);
    2686       }
    2687     }
    2688     if (e.isZero())
    2689       break;
    2690   }
    2691 
    2692 #ifdef DEBUGOUTPUT
    2693     CanonicalForm test= 0;
    2694     j= p;
    2695     for (CFListIterator i= result; i.hasItem(); i++, j++)
    2696       test += mod (i.getItem()*j.getItem(), power (y, d));
    2697     DEBOUTLN (cerr, "test= " << test);
    2698 #endif
    2699   return result;
    2700 }
    2701 
    2702 // so far only tested for two! factor Hensel lifting
    27032582static inline
    27042583void
    27052584henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
    27062585            const CFList& diophant, CFMatrix& M, CFArray& Pi,
    2707             const CFList& products, int j, const CFList& MOD)
     2586            const CFList& products, int j, const CFList& MOD, bool& noOneToOne)
    27082587{
    27092588  CanonicalForm E;
    … …  
    27352614
    27362615  // actual lifting
    2737   CFList diophantine2= diophantine (diophant, factors, products, MOD, E);
     2616  CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
     2617                                    noOneToOne);
     2618
     2619  if (noOneToOne)
     2620    return;
    27382621
    27392622  int k= 0;
    … …  
    27802663      if (k != j - k + 1)
    27812664      {
    2782         if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     2665        if ((one.hasTerms() && one.exp() == j - k + 1) &&
     2666            (two.hasTerms() && two.exp() == j - k + 1))
    27832667        {
    27842668          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
    … …  
    27882672          two++;
    27892673        }
    2790         else if (one.exp() == j - k + 1)
     2674        else if (one.hasTerms() && one.exp() == j - k + 1)
    27912675        {
    27922676          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
    … …  
    27942678          one++;
    27952679        }
    2796         else if (two.exp() == j - k + 1)
     2680        else if (two.hasTerms() && two.exp() == j - k + 1)
    27972681        {
    27982682          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
    … …  
    28382722    degBuf= degree (bufFactors[l + 1], x);
    28392723    if (degPi > 0 && degBuf > 0)
     2724    {
    28402725      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
    2841     if (j == 1)
    2842     {
    2843       if (degPi > 0 && degBuf > 0)
    2844         Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
    2845                   (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
    2846                   M (1, l + 1));
    2847       else if (degPi > 0)
    2848         Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
    2849       else if (degBuf > 0)
    2850         Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
    2851     }
    2852     else
    2853     {
    2854       if (degPi > 0 && degBuf > 0)
    2855       {
    2856         uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
    2857         uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
    2858       }
    2859       else if (degPi > 0)
    2860         uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
    2861       else if (degBuf > 0)
    2862       {
    2863         uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
    2864         uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
    2865       }
    2866       Pi[l] += xToJ*uIZeroJ;
    2867     }
     2726      if (j + 2 <= M.rows())
     2727        M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
     2728                                  MOD);
     2729    }
     2730
     2731    if (degPi > 0 && degBuf > 0)
     2732      uIZeroJ= mulMod (Pi[l -1] [0], buf[l + 1], MOD) +
     2733               mulMod (uIZeroJ, bufFactors[l+1] [0], MOD);
     2734    else if (degPi > 0)
     2735      uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD);
     2736    else if (degBuf > 0)
     2737      uIZeroJ= mulMod (Pi[l - 1], buf[1], MOD);
     2738    else
     2739      uIZeroJ= 0;
     2740
     2741    Pi [l] += xToJ*uIZeroJ;
     2742
    28682743    one= bufFactors [l + 1];
    28692744    two= Pi [l - 1];
    2870     if (two.exp() == j + 1)
    2871     {
    2872       if (degBuf > 0 && degPi > 0)
    2873       {
    2874           tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
    2875           two++;
    2876       }
    2877       else if (degPi > 0)
    2878       {
    2879           tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
    2880           two++;
    2881       }
    2882     }
    28832745    if (degBuf > 0 && degPi > 0)
    28842746    {
     2747      while (one.hasTerms() && one.exp() > j) one++;
     2748      while (two.hasTerms() && two.exp() > j) two++;
    28852749      for (k= 1; k <= (int) ceil (j/2.0); k++)
    28862750      {
    28872751        if (k != j - k + 1)
    28882752        {
    2889           if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     2753          if ((one.hasTerms() && one.exp() == j - k + 1) &&
     2754              (two.hasTerms() && two.exp() == j - k + 1))
    28902755          {
    28912756            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
    … …  
    28952760            two++;
    28962761          }
    2897           else if (one.exp() == j - k + 1)
     2762          else if (one.hasTerms() && one.exp() == j - k + 1)
    28982763          {
    28992764            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
    … …  
    29012766            one++;
    29022767          }
    2903           else if (two.exp() == j - k + 1)
     2768          else if (two.hasTerms() && two.exp() == j - k + 1)
    29042769          {
    2905             tmp[l] += mulMod (bufFactors[l + 1] [k],
     2770            tmp[l] += mulMod (bufFactors[l + 1] [k], 
    29062771                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
    29072772            two++;
    2908           }
     2773           }
    29092774        }
    29102775        else
    … …  
    29122777      }
    29132778    }
     2779
     2780    if (degPi >= j + 1 && degBuf >= j + 1)
     2781    {
     2782      if (j + 2 <= M.rows())
     2783        tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
     2784                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
     2785                            , MOD) - M(1,l+1) - M (j + 2,l+1);
     2786    }
     2787    else if (degPi >= j + 1)
     2788    {
     2789      if (degBuf > 0)
     2790        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
     2791      else
     2792        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
     2793    }
     2794    else if (degBuf >= j + 1)
     2795    {
     2796      if (degPi > 0)
     2797        tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
     2798      else
     2799        tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
     2800    }
     2801
    29142802    Pi[l] += tmp[l]*xToJ*F.mvar();
    29152803  }
    2916 
    29172804  return;
    29182805}
    … …  
    29322819}
    29332820
    2934 // so far only for two factor Hensel lifting
    29352821CFList
    29362822henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
    29372823              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
    2938               const CFList& LCs2)
     2824              const CFList& LCs2, bool& bad)
    29392825{
    29402826  CFList buf= factors;
    … …  
    29772863  {
    29782864    if (degree (bufFactors[i], y) > 0)
    2979       products.append (M (1, 1)/bufFactors[i] [0]);
    2980     else
    2981       products.append (M (1, 1)/bufFactors[i]);
     2865      products.append (eval.getFirst()/bufFactors[i] [0]);
     2866    else
     2867      products.append (eval.getFirst()/bufFactors[i]);
    29822868  }
    29832869
    29842870  for (int d= 1; d < l[1]; d++)
    2985     henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
     2871  {
     2872    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
     2873    if (bad)
     2874      return CFList();
     2875  }
    29862876  CFList result;
    29872877  for (k= 0; k < factors.length(); k++)
    … …  
    29942884henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
    29952885            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
    2996             lNew, const CFList& LCs1, const CFList& LCs2)
     2886            lNew, const CFList& LCs1, const CFList& LCs2, bool& bad)
    29972887{
    29982888  int k= 0;
    … …  
    30202910    if (degree (bufFactors[i], y) > 0)
    30212911    {
    3022       ASSERT (fdivides (bufFactors[i][0], M (1, 1)), "expected exact division");
    3023       products.append (M (1, 1)/bufFactors[i] [0]);
    3024     }
    3025     else
    3026     {
    3027       ASSERT (fdivides (bufFactors[i], M (1, 1)), "expected exact division");
    3028       products.append (M (1, 1)/bufFactors[i]);
     2912      if (!fdivides (bufFactors[i] [0], F.getFirst()))
     2913      {
     2914        bad= true;
     2915        return CFList();
     2916      }
     2917      products.append (F.getFirst()/bufFactors[i] [0]);
     2918    }
     2919    else
     2920    {
     2921      if (!fdivides (bufFactors[i], F.getFirst()))
     2922      {
     2923        bad= true;
     2924        return CFList();
     2925      }
     2926      products.append (F.getFirst()/bufFactors[i]);
    30292927    }
    30302928  }
    30312929
    30322930  for (int d= 1; d < lNew; d++)
    3033     henselStep2 (F.getLast(),buf,bufFactors, diophant, M, Pi, products, d, MOD);
     2931  {
     2932    henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
     2933    if (bad)
     2934      return CFList();
     2935  }
    30342936
    30352937  CFList result;
    … …  
    30392941}
    30402942
    3041 // so far only for two! factor Hensel lifting
    30422943CFList
    30432944henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
    30442945             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
    3045              const CFArray& Pi, const CFList& diophant)
     2946             const CFArray& Pi, const CFList& diophant, bool& bad)
    30462947{
    30472948  CFList bufDiophant= diophant;
    … …  
    30512952  CFArray bufPi= Pi;
    30522953  CFMatrix M= CFMatrix (l[1], factors.length());
    3053   CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2);
     2954  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2,
     2955                               bad);
     2956  if (bad)
     2957    return CFList();
     2958
    30542959  if (eval.length() == 2)
    30552960    return result;
    … …  
    30772982    M= CFMatrix (l[i], factors.length());
    30782983    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
    3079                          l[i], bufLCs1, bufLCs2);
     2984                         l[i], bufLCs1, bufLCs2, bad);
     2985    if (bad)
     2986      return CFList();
    30802987    MOD.append (power (Variable (i + 2), l[i]));
    30812988    bufEval.removeFirst();
    … …  
    30862993}
    30872994
     2995CFList
     2996nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
     2997                      CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
     2998                      int bivarLiftBound, bool& bad)
     2999{
     3000  CFList bufFactors2= factors;
     3001
     3002  Variable y= Variable (2);
     3003  for (CFListIterator i= bufFactors2; i.hasItem(); i++)
     3004    i.getItem()= mod (i.getItem(), y);
     3005
     3006  CanonicalForm bufF= F;
     3007  bufF= mod (bufF, y);
     3008  bufF= mod (bufF, Variable (3));
     3009
     3010  diophant= diophantine (bufF, bufFactors2);
     3011
     3012  CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
     3013
     3014  Pi= CFArray (bufFactors2.length() - 1);
     3015
     3016  CFArray bufFactors= CFArray (bufFactors2.length());
     3017  CFListIterator j= LCs;
     3018  int i= 0;
     3019  for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
     3020    bufFactors[i]= replaceLC (k.getItem(), j.getItem());
     3021
     3022  //initialise Pi
     3023  Variable v= Variable (3);
     3024  CanonicalForm yToL= power (y, bivarLiftBound);
     3025  if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
     3026  {
     3027    M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
     3028    Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
     3029                       mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
     3030  }
     3031  else if (degree (bufFactors[0], v) > 0)
     3032  {
     3033    M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
     3034    Pi [0]=  M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
     3035  }
     3036  else if (degree (bufFactors[1], v) > 0)
     3037  {
     3038    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
     3039    Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
     3040  }
     3041  else
     3042  {
     3043    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
     3044    Pi [0]= M (1,1);
     3045  }
     3046
     3047  for (i= 1; i < Pi.size(); i++)
     3048  {
     3049    if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
     3050    {
     3051      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
     3052      Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
     3053                       mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
     3054    }
     3055    else if (degree (Pi[i-1], v) > 0)
     3056    {
     3057      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
     3058      Pi [i]=  M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
     3059    }
     3060    else if (degree (bufFactors[i+1], v) > 0)
     3061    {
     3062      M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
     3063      Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
     3064    }
     3065    else
     3066    {
     3067      M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
     3068      Pi [i]= M (1,i+1);
     3069    }
     3070  }
     3071
     3072  CFList products;
     3073  bufF= mod (F, Variable (3));
     3074  for (CFListIterator k= factors; k.hasItem(); k++)
     3075    products.append (bufF/k.getItem());
     3076
     3077  CFList MOD= CFList (power (v, liftBound));
     3078  MOD.insert (yToL);
     3079  for (int d= 1; d < liftBound; d++)
     3080  {
     3081    henselStep2 (F, factors, bufFactors, diophant, M, Pi, products, d, MOD, bad);
     3082    if (bad)
     3083      return CFList();
     3084  }
     3085
     3086  CFList result;
     3087  for (i= 0; i < factors.length(); i++)
     3088    result.append (bufFactors[i]);
     3089  return result;
     3090}
     3091
     3092CFList
     3093nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
     3094                    CFList& diophant, CFArray& Pi, CFMatrix& M, const int lOld,
     3095                    int& lNew, const CFList& MOD, bool& noOneToOne
     3096                   )
     3097{
     3098
     3099  int k= 0;
     3100  CFArray bufFactors= CFArray (factors.length());
     3101  CFListIterator j= LCs;
     3102  for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
     3103    bufFactors [k]= replaceLC (i.getItem(), j.getItem());
     3104
     3105  Variable y= F.getLast().mvar();
     3106  Variable x= F.getFirst().mvar();
     3107  CanonicalForm xToLOld= power (x, lOld);
     3108
     3109  Pi [0]= mod (Pi[0], xToLOld);
     3110  M (1, 1)= Pi [0];
     3111
     3112  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
     3113    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
     3114                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
     3115  else if (degree (bufFactors[0], y) > 0)
     3116    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
     3117  else if (degree (bufFactors[1], y) > 0)
     3118    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
     3119
     3120  for (int i= 1; i < Pi.size(); i++)
     3121  {
     3122    Pi [i]= mod (Pi [i], xToLOld);
     3123    M (1, i + 1)= Pi [i];
     3124
     3125    if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
     3126      Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
     3127                 mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
     3128    else if (degree (Pi[i-1], y) > 0)
     3129      Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
     3130    else if (degree (bufFactors[i+1], y) > 0)
     3131      Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
     3132  }
     3133
     3134  CFList products;
     3135  CanonicalForm bufF= F.getFirst();
     3136
     3137  for (int i= 0; i < bufFactors.size(); i++)
     3138  {
     3139    if (degree (bufFactors[i], y) > 0)
     3140    {
     3141      if (!fdivides (bufFactors[i] [0], bufF))
     3142      {
     3143        noOneToOne= true;
     3144        return factors;
     3145      }
     3146      products.append (bufF/bufFactors[i] [0]);
     3147    }
     3148    else
     3149    {
     3150      if (!fdivides (bufFactors[i], bufF))
     3151      {
     3152        noOneToOne= true;
     3153        return factors;
     3154      }
     3155      products.append (bufF/bufFactors[i]);
     3156    }
     3157  }
     3158
     3159  for (int d= 1; d < lNew; d++)
     3160  {
     3161    henselStep2 (F.getLast(), factors, bufFactors, diophant, M, Pi, products, d,
     3162                 MOD, noOneToOne);
     3163    if (noOneToOne)
     3164      return CFList();
     3165  }
     3166
     3167  CFList result;
     3168  for (k= 0; k < factors.length(); k++)
     3169    result.append (bufFactors[k]);
     3170  return result;
     3171}
     3172
     3173CFList
     3174nonMonicHenselLift (const CFList& eval, const CFList& factors,
     3175                    CFList* const& LCs, CFList& diophant, CFArray& Pi,
     3176                    int* liftBound, int length, bool& noOneToOne
     3177                   )
     3178{
     3179  CFList bufDiophant= diophant;
     3180  CFList buf= factors;
     3181  CFArray bufPi= Pi;
     3182  CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
     3183  int k= 0;
     3184
     3185  CFList result=
     3186  nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
     3187                        liftBound[1], liftBound[0], noOneToOne);
     3188
     3189  if (noOneToOne)
     3190    return CFList();
     3191
     3192  if (eval.length() == 1)
     3193    return result;
     3194
     3195  k++;
     3196  CFList MOD;
     3197  for (int i= 0; i < 2; i++)
     3198    MOD.append (power (Variable (i + 2), liftBound[i]));
     3199
     3200  CFListIterator j= eval;
     3201  CFList bufEval;
     3202  bufEval.append (j.getItem());
     3203  j++;
     3204
     3205  for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
     3206  {
     3207    bufEval.append (j.getItem());
     3208    M= CFMatrix (liftBound[i], factors.length() - 1);
     3209    result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
     3210                                liftBound[i-1], liftBound[i], MOD, noOneToOne);
     3211    if (noOneToOne)
     3212      return result;
     3213    MOD.append (power (Variable (i + 2), liftBound[i]));
     3214    bufEval.removeFirst();
     3215  }
     3216
     3217  return result;
     3218}
     3219
    30883220#endif
    30893221/* HAVE_NTL */

  • factory/facHensel.h

    r1bebbf re368746  
    243243              CFArray& Pi,            ///< [in,out] stores intermediate results
    244244              CFList& diophant,       ///< [in,out] result of diophantine()
    245               CFMatrix& M             ///< [in,out] stores intermediate results
     245              CFMatrix& M,            ///< [in,out] stores intermediate results
     246              bool sort= true         ///< [in] sort factors by degree in
     247                                      ///< Variable(1)
    246248             );
    247249
    … …  
    383385                                   ///< leading coefficient
    384386            const int* l,          ///< [in] lifting bounds
    385             const int lLength     ///< [in] length of l
     387            const int lLength,     ///< [in] length of l
     388            bool sort= true        ///< [in] sort factors by degree in
     389                                   ///< Variable(1)
    386390           );
    387391
    … …  
    422426                                   ///< factor
    423427             const CFArray& Pi,    ///< [in] intermediate result
    424              const CFList& diophant///< [in] result of diophantine
     428             const CFList& diophant,///< [in] result of diophantine
     429             bool& noOneToOne      ///< [in,out] check for one to one
     430                                   ///< correspondence
    425431            );
    426432
     433/// Hensel lifting of non monic factors, needs correct leading coefficients of
     434/// factors and a one to one corresponds between bivariate and multivariate
     435/// factors to succeed
     436///
     437/// @return @a nonMonicHenselLift returns a list of lifted factors
     438/// such that prod (factors) == eval.getLast() if there is a one to one
     439/// correspondence
     440CFList
     441nonMonicHenselLift (const CFList& eval,    ///< [in] a list of polys the last
     442                                           ///< element is a compressed
     443                                           ///< multivariate poly, last but one
     444                                           ///< element equals the last elements
     445                                           ///< modulo its main variable and so
     446                                           ///< on. The first element is a
     447                                           ///< compressed poly in 3 variables
     448                    const CFList& factors, ///< [in] a list of bivariate factors
     449                    CFList* const& LCs,    ///< [in] leading coefficients,
     450                                           ///< evaluated in the same way as
     451                                           ///< eval
     452                    CFList& diophant,      ///< [in, out] solution of univariate
     453                                           ///< diophantine equation
     454                    CFArray& Pi,           ///< [in, out] buffer intermediate
     455                                           ///< results
     456                    int* liftBound,        ///< [in] lifting bounds
     457                    int length,            ///< [in] length of lifting bounds
     458                    bool& noOneToOne       ///< [in, out] check for one to one
     459                                           ///< correspondence
     460                   );
    427461#endif
    428462/* FAC_HENSEL_H */
Note: See TracChangeset for help on using the changeset viewer.