source: git/factory/facFqBivar.cc @ dd8047

fieker-DuValspielwiese
Last change on this file since dd8047 was dd8047, checked in by Martin Lee <martinlee84@…>, 13 years ago
fix in bivariate factorization git-svn-id: file:///usr/local/Singular/svn/trunk@14095 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 38.0 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facFqBivar.cc
5 *
6 * This file provides functions for factorizing a bivariate polynomial over
7 * \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF, based on "Modern Computer
8 * Algebra, Chapter 15" by J. von zur Gathen & J. Gerhard and "Factoring
9 * multivariate polynomials over a finite field" by L. Bernardin.
10 *
11 * ABSTRACT: In contrast to biFactorizer() in facFqFactorice.cc we evaluate and
12 * factorize the polynomial in both variables. So far factor recombination is
13 * done naive!
14 *
15 * @author Martin Lee
16 *
17 * @internal @version \$Id$
18 *
19 **/
20/*****************************************************************************/
21
22#include <config.h>
23
24#include "assert.h"
25#include "debug.h"
26#include "timing.h"
27
28#include "canonicalform.h"
29#include "cf_defs.h"
30#include "cf_map_ext.h"
31#include "cf_random.h"
32#include "facHensel.h"
33#include "cf_map.h"
34#include "cf_gcd_smallp.h"
35#include "facFqBivar.h"
36
37
38#ifdef HAVE_NTL
39#include "NTLconvert.h"
40
41TIMING_DEFINE_PRINT(fac_uni_factorizer);
42TIMING_DEFINE_PRINT(fac_hensel_lift);
43TIMING_DEFINE_PRINT(fac_factor_recombination);
44
45CanonicalForm prodMod0 (const CFList& L, const CanonicalForm& M)
46{
47  if (L.isEmpty())
48    return 1;
49  else if (L.length() == 1)
50    return mod (L.getFirst()(0, 1) , M);
51  else if (L.length() == 2)
52    return mod (L.getFirst()(0, 1)*L.getLast()(0, 1), M);
53  else
54  {
55    int l= L.length()/2;
56    CFListIterator i= L;
57    CFList tmp1, tmp2;
58    CanonicalForm buf1, buf2;
59    for (int j= 1; j <= l; j++, i++)
60      tmp1.append (i.getItem());
61    tmp2= Difference (L, tmp1);
62    buf1= prodMod0 (tmp1, M);
63    buf2= prodMod0 (tmp2, M);
64    return mod (buf1*buf2, M);
65  }
66}
67
68CanonicalForm evalPoint (const CanonicalForm& F, CanonicalForm & eval,
69                         const Variable& alpha, CFList& list, const bool& GF,
70                         bool& fail)
71{
72  fail= false;
73  Variable x= Variable(2);
74  FFRandom genFF;
75  GFRandom genGF;
76  CanonicalForm random, mipo;
77  double bound;
78  int p= getCharacteristic ();
79  if (alpha.level() != 1)
80  {
81    mipo= getMipo (alpha);
82    int d= degree (mipo);
83    bound= ipower (p, d);
84  }
85  else if (GF)
86  {
87    int d= getGFDegree();
88    bound= ipower (p, d);
89  }
90  else
91    bound= p;
92
93  do
94  {
95    if (list.length() >= bound)
96    {
97      fail= true;
98      break;
99    }
100    if (list.isEmpty())
101      random= 0;
102    else if (GF)
103      random= genGF.generate();
104    else if (list.length() < p || alpha == x)
105      random= genFF.generate();
106    else if (alpha != x && list.length() >= p)
107    {
108      AlgExtRandomF genAlgExt (alpha);
109      random= genAlgExt.generate();
110    }
111    if (find (list, random)) continue;
112    eval= F (random, x);
113    if (degree (eval) != degree (F, Variable (1)))
114    { //leading coeff vanishes
115      if (!find (list, random))
116        list.append (random);
117      continue;
118    }
119    if (degree (gcd (deriv (eval, eval.mvar()), eval), eval.mvar()) > 0)
120    { //evaluated polynomial is not squarefree
121      if (!find (list, random))
122        list.append (random);
123      continue;
124    }
125  } while (find (list, random));
126
127  return random;
128}
129
130CFList
131uniFactorizer (const CanonicalForm& A, const Variable& alpha, const bool& GF)
132{
133  Variable x= A.mvar();
134  ASSERT (A.isUnivariate(), "univariate polynomial expected");
135  CFFList factorsA;
136  ZZ p= to_ZZ (getCharacteristic());
137  ZZ_p::init (p);
138  if (GF)
139  {
140    Variable beta= rootOf (gf_mipo);
141    int k= getGFDegree();
142    char cGFName= gf_name;
143    setCharacteristic (getCharacteristic());
144    CanonicalForm buf= GF2FalphaRep (A, beta);
145    if (getCharacteristic() > 2)
146    {
147      ZZ_pX NTLMipo= convertFacCF2NTLZZpX (gf_mipo);
148      ZZ_pE::init (NTLMipo);
149      ZZ_pEX NTLA= convertFacCF2NTLZZ_pEX (buf, NTLMipo);
150      MakeMonic (NTLA);
151      vec_pair_ZZ_pEX_long NTLFactorsA= CanZass (NTLA);
152      ZZ_pE multi= to_ZZ_pE (1);
153      factorsA= convertNTLvec_pair_ZZpEX_long2FacCFFList (NTLFactorsA, multi,
154                                                         x, beta);
155    }
156    else
157    {
158      GF2X NTLMipo= convertFacCF2NTLGF2X (gf_mipo);
159      GF2E::init (NTLMipo);
160      GF2EX NTLA= convertFacCF2NTLGF2EX (buf, NTLMipo);
161      MakeMonic (NTLA);
162      vec_pair_GF2EX_long NTLFactorsA= CanZass (NTLA);
163      GF2E multi= to_GF2E (1);
164      factorsA= convertNTLvec_pair_GF2EX_long2FacCFFList (NTLFactorsA, multi,
165                                                           x, beta);
166    }
167    setCharacteristic (getCharacteristic(), k, cGFName);
168    for (CFFListIterator i= factorsA; i.hasItem(); i++)
169    {
170      buf= i.getItem().factor();
171      buf= Falpha2GFRep (buf);
172      i.getItem()= CFFactor (buf, i.getItem().exp());
173    }
174  }
175  else if (alpha.level() != 1)
176  {
177    if (getCharacteristic() > 2)
178    {
179      ZZ_pX NTLMipo= convertFacCF2NTLZZpX (getMipo (alpha));
180      ZZ_pE::init (NTLMipo);
181      ZZ_pEX NTLA= convertFacCF2NTLZZ_pEX (A, NTLMipo);
182      MakeMonic (NTLA);
183      vec_pair_ZZ_pEX_long NTLFactorsA= CanZass (NTLA);
184      ZZ_pE multi= to_ZZ_pE (1);
185      factorsA= convertNTLvec_pair_ZZpEX_long2FacCFFList (NTLFactorsA, multi,
186                                                           x, alpha);
187    }
188    else
189    {
190      GF2X NTLMipo= convertFacCF2NTLGF2X (getMipo (alpha));
191      GF2E::init (NTLMipo);
192      GF2EX NTLA= convertFacCF2NTLGF2EX (A, NTLMipo);
193      MakeMonic (NTLA);
194      vec_pair_GF2EX_long NTLFactorsA= CanZass (NTLA);
195      GF2E multi= to_GF2E (1);
196      factorsA= convertNTLvec_pair_GF2EX_long2FacCFFList (NTLFactorsA, multi,
197                                                           x, alpha);
198    }
199  }
200  else
201  {
202    if (getCharacteristic() > 2)
203    {
204      ZZ_pX NTLA= convertFacCF2NTLZZpX (A);
205      MakeMonic (NTLA);
206      vec_pair_ZZ_pX_long NTLFactorsA= CanZass (NTLA);
207      ZZ_p multi= to_ZZ_p (1);
208      factorsA= convertNTLvec_pair_ZZpX_long2FacCFFList (NTLFactorsA, multi,
209                                                          x);
210    }
211    else
212    {
213      GF2X NTLA= convertFacCF2NTLGF2X (A);
214      vec_pair_GF2X_long NTLFactorsA= CanZass (NTLA);
215      GF2 multi= to_GF2 (1);
216      factorsA= convertNTLvec_pair_GF2X_long2FacCFFList (NTLFactorsA, multi,
217                                                          x);
218    }
219  }
220  CFList uniFactors;
221  for (CFFListIterator i= factorsA; i.hasItem(); i++)
222    uniFactors.append (i.getItem().factor());
223  return uniFactors;
224}
225
226/// naive factor recombination as decribed in "Factoring
227/// multivariate polynomials over a finite field" by L Bernardin.
228CFList
229extFactorRecombination (const CFList& factors, const CanonicalForm& F,
230                        const CanonicalForm& M, const ExtensionInfo& info,
231                        const DegreePattern& degs, const CanonicalForm& eval)
232{
233  if (factors.length() == 0)
234    return CFList();
235
236  Variable alpha= info.getAlpha();
237  Variable beta= info.getBeta();
238  CanonicalForm gamma= info.getGamma();
239  CanonicalForm delta= info.getDelta();
240  int k= info.getGFDegree();
241
242  Variable y= F.mvar();
243  CFList source, dest;
244  if (degs.getLength() <= 1 || factors.length() == 1)
245    return CFList(mapDown (F(y-eval, y), info, source, dest));
246
247  DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, M) == F " <<
248            (LC (F, 1)*prodMod (factors, M) == F));
249  int degMipoBeta;
250  if (!k && beta.level() == 1)
251    degMipoBeta= 1;
252  else if (!k && beta.level() != 1)
253    degMipoBeta= degree (getMipo (beta));
254
255  CFList T, S, Diff;
256  T= factors;
257
258  int s= 1;
259  CFList result;
260  CanonicalForm buf, buf2;
261  if (beta != Variable (1))
262    buf= mapDown (F, gamma, delta, alpha, source, dest);
263  else
264    buf= F;
265
266  CanonicalForm g, LCBuf= LC (buf, Variable (1));
267  int v [T.length()];
268  for (int i= 0; i < T.length(); i++)
269    v[i]= 0;
270
271  CFArray TT;
272  DegreePattern bufDegs1, bufDegs2;
273  bufDegs1= degs;
274  int subsetDeg;
275  TT= copy (factors);
276  bool nosubset= false;
277  bool recombination= false;
278  bool trueFactor= false;
279  while (T.length() >= 2*s)
280  {
281    while (nosubset == false)
282    {
283      if (T.length() == s)
284      {
285        if (recombination)
286        {
287          T.insert (LCBuf);
288          g= prodMod (T, M);
289          T.removeFirst();
290          g /= content(g);
291          g= g (y - eval, y);
292          g /= Lc (g);
293          appendTestMapDown (result, g, info, source, dest);
294          return result;
295        }
296        else
297        {
298          appendMapDown (result, F (y - eval, y), info, source, dest);
299          return result;
300        }
301      }
302      S= subset (v, s, TT, nosubset);
303      if (nosubset) break;
304      subsetDeg= subsetDegree (S);
305      // skip those combinations that are not possible
306      if (!degs.find (subsetDeg))
307        continue;
308      else
309      {
310        g= prodMod0 (S, M);
311        g= mod (g*LCBuf, M);
312        g /= content (g);
313        if (fdivides (LC (g), LCBuf))
314        {
315          S.insert (LCBuf);
316          g= prodMod (S, M);
317          S.removeFirst();
318          g /= content (g, Variable (1));
319          if (fdivides (g, buf))
320          {
321            buf2= g (y - eval, y);
322            buf2 /= Lc (buf2);
323
324            if (!k && beta == Variable (1))
325            {
326              if (degree (buf2, alpha) < degMipoBeta)
327              {
328                buf /= g;
329                LCBuf= LC (buf, Variable (1));
330                recombination= true;
331                appendTestMapDown (result, buf2, info, source, dest);
332                trueFactor= true;
333              }
334            }
335            else
336            {
337              if (!isInExtension (buf2, delta, k))
338              {
339                buf /= g;
340                LCBuf= LC (buf, Variable (1));
341                recombination= true;
342                appendTestMapDown (result, buf2, info, source, dest);
343                trueFactor= true;
344              }
345            }
346            if (trueFactor)
347            {
348              T= Difference (T, S);
349              // compute new possible degree pattern
350              bufDegs2= DegreePattern (T);
351              bufDegs1.intersect (bufDegs2);
352              bufDegs1.refine ();
353              if (T.length() < 2*s || T.length() == s ||
354                  bufDegs1.getLength() == 1)
355              {
356                if (recombination)
357                {
358                  appendTestMapDown (result, buf (y - eval, y), info, source,
359                                     dest);
360                  return result;
361                }
362                else
363                {
364                  appendMapDown (result, F (y - eval, y), info, source, dest);
365                  return result;
366                }
367              }
368              trueFactor= false;
369              TT= copy (T);
370              indexUpdate (v, s, T.length(), nosubset);
371              if (nosubset) break;
372            }
373          }
374        }
375      }
376    }
377    s++;
378    if (T.length() < 2*s || T.length() == s)
379    {
380      if (recombination)
381      {
382        appendTestMapDown (result, buf (y - eval, y), info, source, dest);
383        return result;
384      }
385      else
386      {
387        appendMapDown (result, F (y - eval, y), info, source, dest);
388        return result;
389      }
390    }
391    int v [T.length()];
392    for (int i= 0; i < T.length(); i++)
393      v[i]= 0;
394    nosubset= false;
395  }
396  if (T.length() < 2*s)
397    appendMapDown (result, F (y - eval, y), info, source, dest);
398
399  return result;
400}
401
402/// naive factor recombination as decribed in "Factoring
403/// multivariate polynomials over a finite field" by L Bernardin.
404CFList
405factorRecombination (const CFList& factors, const CanonicalForm& F,
406                     const CanonicalForm& M, const DegreePattern& degs)
407{
408  if (factors.length() == 0)
409    return CFList ();
410  if (degs.getLength() <= 1 || factors.length() == 1)
411    return CFList(F);
412  DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, M) == F " <<
413            (LC (F, 1)*prodMod (factors, M) == F));
414  CFList T, S;
415
416  T= factors;
417  int s= 1;
418  CFList result;
419  CanonicalForm LCBuf= LC (F, Variable (1));
420  CanonicalForm g, buf= F;
421  int v [T.length()];
422  for (int i= 0; i < T.length(); i++)
423    v[i]= 0;
424  bool nosubset= false;
425  CFArray TT;
426  DegreePattern bufDegs1, bufDegs2;
427  bufDegs1= degs;
428  int subsetDeg;
429  TT= copy (factors);
430  bool recombination= false;
431  while (T.length() >= 2*s)
432  {
433    while (nosubset == false)
434    {
435      if (T.length() == s)
436      {
437        if (recombination)
438        {
439          T.insert (LCBuf);
440          g= prodMod (T, M);
441          T.removeFirst();
442          result.append (g/content (g, Variable (1)));
443          return result;
444        }
445        else
446          return CFList (F);
447      }
448      S= subset (v, s, TT, nosubset);
449      if (nosubset) break;
450      subsetDeg= subsetDegree (S);
451      // skip those combinations that are not possible
452      if (!degs.find (subsetDeg))
453        continue;
454      else
455      {
456        g= prodMod0 (S, M);
457        g= mod (g*LCBuf, M);
458        g /= content (g);
459        if (fdivides (LC(g), LCBuf))
460        {
461          S.insert (LCBuf);
462          g= prodMod (S, M);
463          S.removeFirst();
464          g /= content (g, Variable (1));
465          if (fdivides (g, buf))
466          {
467            recombination= true;
468            result.append (g);
469            buf /= g;
470            LCBuf= LC (buf, Variable(1));
471            T= Difference (T, S);
472
473            // compute new possible degree pattern
474            bufDegs2= DegreePattern (T);
475            bufDegs1.intersect (bufDegs2);
476            bufDegs1.refine ();
477            if (T.length() < 2*s || T.length() == s ||
478                bufDegs1.getLength() == 1)
479            {
480              if (recombination)
481              {
482                result.append (buf);
483                return result;
484              }
485              else
486                return CFList (F);
487            }
488            TT= copy (T);
489            indexUpdate (v, s, T.length(), nosubset);
490            if (nosubset) break;
491          }
492        }
493      }
494    }
495    s++;
496    if (T.length() < 2*s || T.length() == s)
497    {
498      if (recombination)
499      {
500        result.append (buf);
501        return result;
502      }
503      else
504        return CFList (F);
505    }
506    int v [T.length()];
507    for (int i= 0; i < T.length(); i++)
508      v[i]= 0;
509    nosubset= false;
510  }
511  if (T.length() < 2*s)
512    result.append (F);
513
514  return result;
515}
516
517Variable chooseExtension (const CanonicalForm & A, const Variable & alpha)
518{
519  int p= getCharacteristic();
520  ZZ NTLp= to_ZZ (p);
521  ZZ_p::init (NTLp);
522  ZZ_pX NTLirredpoly;
523  int d= degree (A);
524  int m= 1;
525  if (alpha != Variable (1))
526    m= degree (getMipo (alpha));
527  int i= (int) floor((double) d/(double) m) - 1;
528  if (i < 2)
529    i= 2;
530  if (i > 8)
531    i= i/4;
532  BuildIrred (NTLirredpoly, i*m);
533  Variable x= A.mvar();
534  CanonicalForm newMipo= convertNTLZZpX2CF (NTLirredpoly, x);
535  return rootOf (newMipo);
536}
537
538CFList
539earlyFactorDetection (CanonicalForm& F, CFList& factors,int& adaptedLiftBound,
540                      DegreePattern& degs, bool& success, int deg)
541{
542  DegreePattern bufDegs1= degs;
543  DegreePattern bufDegs2;
544  CFList result;
545  CFList T= factors;
546  CanonicalForm buf= F;
547  CanonicalForm LCBuf= LC (buf, Variable (1));
548  CanonicalForm g;
549  CanonicalForm M= power (F.mvar(), deg);
550  adaptedLiftBound= 0;
551  int d;
552  if (degree (LCBuf) == degree (F))
553    d= degree (F);
554  else
555    d= degree (F) + degree (LCBuf);
556  for (CFListIterator i= factors; i.hasItem(); i++)
557  {
558    if (!bufDegs1.find (degree (i.getItem(), 1)))
559      continue;
560    else
561    {
562      g= i.getItem() (0, 1);
563      g *= LCBuf;
564      g= mod (g, M);
565      if (fdivides (LC (g), LCBuf))
566      {
567        g= mulMod2 (i.getItem(), LCBuf, M);
568        g /= content (g, Variable (1));
569        if (fdivides (g, buf))
570        {
571          result.append (g);
572          buf /= g;
573          d -= degree (g) + degree (LC (g, Variable (1)));
574          LCBuf= LC (buf, Variable (1));
575          T= Difference (T, CFList (i.getItem()));
576
577          // compute new possible degree pattern
578          bufDegs2= DegreePattern (T);
579          bufDegs1.intersect (bufDegs2);
580          bufDegs1.refine ();
581          if (bufDegs1.getLength() <= 1)
582          {
583            result.append (buf);
584            break;
585          }
586        }
587      }
588    }
589  }
590  adaptedLiftBound= d + 1;
591  if (d < deg)
592  {
593    factors= T;
594    degs= bufDegs1;
595    F= buf;
596    success= true;
597  }
598  if (bufDegs1.getLength() <= 1)
599    degs= bufDegs1;
600  return result;
601}
602
603CFList
604extEarlyFactorDetection (CanonicalForm& F, CFList& factors,
605                         int& adaptedLiftBound, DegreePattern& degs,
606                         bool& success, const ExtensionInfo& info,
607                         const CanonicalForm& eval, int deg)
608{
609  Variable alpha= info.getAlpha();
610  Variable beta= info.getBeta();
611  CanonicalForm gamma= info.getGamma();
612  CanonicalForm delta= info.getDelta();
613  int k= info.getGFDegree();
614  DegreePattern bufDegs1= degs, bufDegs2;
615  CFList result;
616  CFList T= factors;
617  Variable y= F.mvar();
618  CanonicalForm buf= F, LCBuf= LC (buf, Variable (1)), g, buf2;
619  CanonicalForm M= power (y, deg);
620  adaptedLiftBound= 0;
621  bool trueFactor= false;
622  int d;
623  if (degree (F) == degree (LCBuf))
624    d= degree (F);
625  else
626    d= degree (F) + degree (LCBuf);
627  CFList source, dest;
628  int degMipoBeta;
629  if (!k && beta.level() == 1)
630    degMipoBeta= 1;
631  else if (!k && beta.level() != 1)
632    degMipoBeta= degree (getMipo (beta));
633  for (CFListIterator i= factors; i.hasItem(); i++)
634  {
635    if (!bufDegs1.find (degree (i.getItem(), 1)))
636      continue;
637    else
638    {
639      g= i.getItem() (0, 1);
640      g *= LCBuf;
641      g= mod (g, M);
642      if (fdivides (LC (g), LCBuf))
643      {
644        g= mulMod2 (i.getItem(), LCBuf, M);
645        g /= content (g, Variable (1));
646        if (fdivides (g, buf))
647        {
648          buf2= g (y - eval, y);
649          buf2 /= Lc (buf2);
650
651          if (!k && beta == Variable (1))
652          {
653            if (degree (buf2, alpha) < degMipoBeta)
654            {
655              appendTestMapDown (result, buf2, info, source, dest);
656              buf /= g;
657              d -= degree (g) + degree (LC (g, Variable (1)));
658              LCBuf= LC (buf, Variable (1));
659              trueFactor= true;
660            }
661          }
662          else
663          {
664            if (!isInExtension (buf2, delta, k))
665            {
666              appendTestMapDown (result, buf2, info, source, dest);
667              buf /= g;
668              d -= degree (g) + degree (LC (g, Variable (1)));
669              LCBuf= LC (buf, Variable (1));
670              trueFactor= true;
671            }
672          }
673          if (trueFactor)
674          {
675            T= Difference (T, CFList (i.getItem()));
676
677            // compute new possible degree pattern
678            bufDegs2= DegreePattern (T);
679            bufDegs1.intersect (bufDegs2);
680            bufDegs1.refine ();
681            trueFactor= false;
682            if (bufDegs1.getLength() <= 1)
683            {
684              buf= buf (y - eval, y);
685              buf /= Lc (buf);
686              appendMapDown (result, buf, info, source, dest);
687              break;
688            }
689          }
690        }
691      }
692    }
693  }
694  adaptedLiftBound= d + 1;
695  if (adaptedLiftBound < deg)
696  {
697    success= true;
698    factors= T;
699    degs= bufDegs1;
700    F= buf;
701  }
702  if (bufDegs1.getLength() <= 1)
703    degs= bufDegs1;
704
705  return result;
706}
707
708CFList
709henselLiftAndEarly (CanonicalForm& A, bool& earlySuccess, CFList&
710                    earlyFactors, DegreePattern& degs, int& liftBound,
711                    const CFList& uniFactors, const ExtensionInfo& info,
712                    const CanonicalForm& eval)
713{
714  Variable alpha= info.getAlpha();
715  Variable beta= info.getBeta();
716  CanonicalForm gamma= info.getGamma();
717  CanonicalForm delta= info.getDelta();
718  int k= info.getGFDegree();
719  bool extension= info.isInExtension();
720
721  Variable x= Variable (1);
722  Variable y= Variable (2);
723  CFArray Pi;
724  CFList diophant;
725  CFList bufUniFactors= uniFactors;
726  bufUniFactors.insert (LC (A, x));
727  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
728  earlySuccess= false;
729  int newLiftBound= 0;
730  int smallFactorDeg= 11; //this is a tunable parameter
731  if (smallFactorDeg >= liftBound)
732    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
733  else if (smallFactorDeg >= degree (A, y) + 1)
734  {
735    henselLift12 (A, bufUniFactors, degree (A, y) + 1, Pi, diophant, M);
736    if (!extension)
737      earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
738                     degs, earlySuccess, degree (A, y) + 1);
739    else
740      earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
741                     newLiftBound, degs, earlySuccess, info, eval,
742                     degree (A, y) + 1);
743    if (degs.getLength() > 1 && !earlySuccess)
744    {
745      if (newLiftBound > degree (A, y) + 1)
746      {
747        liftBound= newLiftBound;
748        bufUniFactors.insert (LC(A, x));
749        henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
750                            Pi, diophant, M);
751      }
752    }
753    else if (earlySuccess)
754      liftBound= newLiftBound;
755  }
756  else if (smallFactorDeg < degree (A, y) + 1)
757  {
758    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
759    if (!extension)
760      earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
761                                           degs, earlySuccess,
762                                           smallFactorDeg);
763    else
764      earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
765                     newLiftBound, degs, earlySuccess,
766                     info, eval, smallFactorDeg);
767    if (degs.getLength() > 1 && !earlySuccess)
768    {
769      bufUniFactors.insert (LC (A, x));
770      henselLiftResume12 (A, bufUniFactors, smallFactorDeg, degree (A, y)
771                          + 1, Pi, diophant, M);
772      if (!extension)
773        earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
774                       degs, earlySuccess, degree (A, y) + 1);
775      else
776        earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
777                       newLiftBound, degs, earlySuccess,
778                       info, eval, degree(A,y) + 1);
779      if (degs.getLength() > 1 && !earlySuccess)
780      {
781        if (newLiftBound > degree (A, y) + 1)
782        {
783          if (newLiftBound < newLiftBound)
784            liftBound= newLiftBound;
785          bufUniFactors.insert (LC(A, x));
786          henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
787                              Pi, diophant, M);
788        }
789      }
790      else if (earlySuccess)
791        liftBound= newLiftBound;
792    }
793    else if (earlySuccess)
794      liftBound= newLiftBound;
795  }
796  if (newLiftBound > 0)
797    liftBound= newLiftBound;
798  return bufUniFactors;
799}
800
801CFList
802extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info);
803
804/// bivariate factorization over finite fields as decribed in "Factoring
805/// multivariate polynomials over a finite field" by L Bernardin.
806CFList
807biFactorize (const CanonicalForm& F, const ExtensionInfo& info)
808{
809  if (F.inCoeffDomain())
810    return CFList(F);
811
812  CanonicalForm A= F;
813  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
814
815  Variable alpha= info.getAlpha();
816  Variable beta= info.getBeta();
817  CanonicalForm gamma= info.getGamma();
818  CanonicalForm delta= info.getDelta();
819  int k= info.getGFDegree();
820  bool extension= info.isInExtension();
821  if (A.isUnivariate())
822  {
823    if (extension == false)
824      return uniFactorizer (F, alpha, GF);
825    else
826    {
827      CFList source, dest;
828      A= mapDown (A, info, source, dest);
829      return uniFactorizer (A, beta, GF);
830    }
831  }
832
833  CFMap N;
834  A= compress (A, N);
835  Variable y= A.mvar();
836
837  if (y.level() > 2) return CFList (F);
838  Variable x= Variable (1);
839
840  //remove and factorize content
841  CanonicalForm contentAx= content (A, x);
842  CanonicalForm contentAy= content (A);
843
844  A= A/(contentAx*contentAy);
845  CFList contentAxFactors, contentAyFactors;
846
847  if (!extension)
848  {
849    contentAxFactors= uniFactorizer (contentAx, alpha, GF);
850    contentAyFactors= uniFactorizer (contentAy, alpha, GF);
851  }
852
853  //trivial case
854  CFList factors;
855  if (A.inCoeffDomain())
856  {
857    append (factors, contentAxFactors);
858    append (factors, contentAyFactors);
859    decompress (factors, N);
860    return factors;
861  }
862  else if (A.isUnivariate())
863  {
864    factors= uniFactorizer (A, alpha, GF);
865    append (factors, contentAxFactors);
866    append (factors, contentAyFactors);
867    decompress (factors, N);
868    return factors;
869  }
870
871  // check derivatives
872  bool derivXZero= false;
873  CanonicalForm derivX= deriv (A, x);
874  CanonicalForm gcdDerivX;
875  if (derivX.isZero())
876    derivXZero= true;
877  else
878  {
879    gcdDerivX= gcd (A, derivX);
880    if (degree (gcdDerivX) > 0)
881    {
882      CanonicalForm g= A/gcdDerivX;
883      CFList factorsG=
884      Union (biFactorize (g, info), biFactorize (gcdDerivX, info));
885      append (factorsG, contentAxFactors);
886      append (factorsG, contentAyFactors);
887      decompress (factorsG, N);
888      normalize (factors);
889      return factorsG;
890    }
891  }
892  bool derivYZero= false;
893  CanonicalForm derivY= deriv (A, y);
894  CanonicalForm gcdDerivY;
895  if (derivY.isZero())
896    derivYZero= true;
897  else
898  {
899    gcdDerivY= gcd (A, derivY);
900    if (degree (gcdDerivY) > 0)
901    {
902      CanonicalForm g= A/gcdDerivY;
903      CFList factorsG=
904      Union (biFactorize (g, info), biFactorize (gcdDerivY, info));
905      append (factorsG, contentAxFactors);
906      append (factorsG, contentAyFactors);
907      decompress (factorsG, N);
908      normalize (factors);
909      return factorsG;
910    }
911  }
912  //main variable is chosen s.t. the degree in x is minimal
913  bool swap= false;
914  if ((degree (A) > degree (A, x)) || derivXZero)
915  {
916    if (!derivYZero)
917    {
918      A= swapvar (A, y, x);
919      swap= derivXZero;
920      derivXZero= derivYZero;
921      derivYZero= swap;
922      swap= true;
923    }
924  }
925
926  bool fail;
927  CanonicalForm Aeval, evaluation, bufAeval, bufEvaluation, buf;
928  CFList uniFactors, list, bufUniFactors;
929  DegreePattern degs;
930  DegreePattern bufDegs;
931
932  bool fail2= false;
933  CanonicalForm Aeval2, evaluation2, bufAeval2, bufEvaluation2;
934  CFList bufUniFactors2, list2, uniFactors2;
935  DegreePattern degs2;
936  DegreePattern bufDegs2;
937  bool swap2= false;
938
939  // several univariate factorizations to obtain more information about the
940  // degree pattern therefore usually less combinations have to be tried during
941  // the recombination process
942  int factorNums= 5;
943  double logarithm= (double) ilog2 (totaldegree (A));
944  logarithm /= log2exp;
945  logarithm= ceil (logarithm);
946  if (factorNums < (int) logarithm)
947    factorNums= (int) logarithm;
948  int subCheck1= substituteCheck (A, x);
949  int subCheck2= substituteCheck (A, y);
950  for (int i= 0; i < factorNums; i++)
951  {
952    bufAeval= A;
953    bufEvaluation= evalPoint (A, bufAeval, alpha, list, GF, fail);
954    if (!derivXZero && !fail2)
955    {
956      buf= swapvar (A, x, y);
957      bufAeval2= buf;
958      bufEvaluation2= evalPoint (buf, bufAeval2, alpha, list2, GF, fail2);
959    }
960    // first try to change main variable if there is no valid evaluation point
961    if (fail && (i == 0))
962    {
963      if (!derivXZero)
964        bufEvaluation= evalPoint (buf, bufAeval, alpha, list, GF, fail);
965      else
966        fail= true;
967
968      if (!fail)
969      {
970        int dummy= subCheck2;
971        subCheck2= subCheck1;
972        subCheck1= dummy;
973        A= buf;
974        swap2= true;
975      }
976    }
977
978    // if there is no valid evaluation point pass to a field extension
979    if (fail && (i == 0))
980    {
981      factors= extBiFactorize (A, info);
982      appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
983                            swap, swap2, N);
984      normalize (factors);
985      return factors;
986    }
987
988    // there is at least one valid evaluation point
989    // but we do not compute more univariate factorization over an extension
990    if (fail && (i != 0))
991      break;
992
993    // univariate factorization
994    TIMING_START (fac_uni_factorizer);
995    bufUniFactors= uniFactorizer (bufAeval, alpha, GF);
996    TIMING_END_AND_PRINT (fac_uni_factorizer,
997                          "time for univariate factorization: ");
998    DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " <<
999              prod (bufUniFactors)*Lc (bufAeval) == bufAeval);
1000
1001    if (!derivXZero && !fail2)
1002    {
1003      TIMING_START (fac_uni_factorizer);
1004      bufUniFactors2= uniFactorizer (bufAeval2, alpha, GF);
1005      TIMING_END_AND_PRINT (fac_uni_factorizer,
1006                            "time for univariate factorization in y: ");
1007      DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " <<
1008                prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2);
1009    }
1010
1011    if (bufUniFactors.length() == 1 ||
1012        (!fail2 && !derivXZero && (bufUniFactors2.length() == 1)))
1013    {
1014      if (extension)
1015      {
1016        CFList source, dest;
1017        ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k,
1018                             info.getGFName(), info.isInExtension());
1019        appendMapDown (factors, A, info2, source, dest);
1020      }
1021      else
1022        factors.append (A);
1023
1024      appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1025                            swap, swap2, N);
1026
1027      normalize (factors);
1028      return factors;
1029    }
1030
1031    if (i == 0)
1032    {
1033      if (subCheck1 > 0)
1034      {
1035        int subCheck= substituteCheck (bufUniFactors);
1036
1037        if (subCheck > 1 && (subCheck1%subCheck == 0))
1038        {
1039          CanonicalForm bufA= A;
1040          subst (bufA, bufA, subCheck, x);
1041          factors= biFactorize (bufA, info);
1042          reverseSubst (factors, subCheck, x);
1043          appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1044                                swap, swap2, N);
1045          normalize (factors);
1046          return factors;
1047        }
1048      }
1049
1050      if (!derivXZero && !fail2 && subCheck2 > 0)
1051      {
1052        int subCheck= substituteCheck (bufUniFactors2);
1053
1054        if (subCheck > 1 && (subCheck2%subCheck == 0))
1055        {
1056          CanonicalForm bufA= A;
1057          subst (bufA, bufA, subCheck, y);
1058          factors= biFactorize (bufA, info);
1059          reverseSubst (factors, subCheck, y);
1060          appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1061                                swap, swap2, N);
1062          normalize (factors);
1063          return factors;
1064        }
1065      }
1066    }
1067
1068    // degree analysis
1069    bufDegs = DegreePattern (bufUniFactors);
1070    if (!derivXZero && !fail2)
1071      bufDegs2= DegreePattern (bufUniFactors2);
1072
1073    if (i == 0)
1074    {
1075      Aeval= bufAeval;
1076      evaluation= bufEvaluation;
1077      uniFactors= bufUniFactors;
1078      degs= bufDegs;
1079      if (!derivXZero && !fail2)
1080      {
1081        Aeval2= bufAeval2;
1082        evaluation2= bufEvaluation2;
1083        uniFactors2= bufUniFactors2;
1084        degs2= bufDegs2;
1085      }
1086    }
1087    else
1088    {
1089      degs.intersect (bufDegs);
1090      if (!derivXZero && !fail2)
1091      {
1092        degs2.intersect (bufDegs2);
1093        if (bufUniFactors2.length() < uniFactors2.length())
1094        {
1095          uniFactors2= bufUniFactors2;
1096          Aeval2= bufAeval2;
1097          evaluation2= bufEvaluation2;
1098        }
1099      }
1100      if (bufUniFactors.length() < uniFactors.length())
1101      {
1102        uniFactors= bufUniFactors;
1103        Aeval= bufAeval;
1104        evaluation= bufEvaluation;
1105      }
1106    }
1107    list.append (bufEvaluation);
1108    if (!derivXZero && !fail2)
1109      list2.append (bufEvaluation2);
1110  }
1111
1112  if (!derivXZero && !fail2)
1113  {
1114    if (uniFactors.length() > uniFactors2.length() ||
1115        (uniFactors.length() == uniFactors2.length()
1116         && degs.getLength() > degs2.getLength()))
1117    {
1118      degs= degs2;
1119      uniFactors= uniFactors2;
1120      evaluation= evaluation2;
1121      Aeval= Aeval2;
1122      A= buf;
1123      swap2= true;
1124    }
1125  }
1126
1127  if (degs.getLength() == 1) // A is irreducible
1128  {
1129    if (extension)
1130    {
1131      CFList source, dest;
1132      ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k,
1133                           info.getGFName(), info.isInExtension());
1134      appendMapDown (factors, A, info2, source, dest);
1135    }
1136    else
1137      factors.append (A);
1138    appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1139                            swap, swap2, N);
1140    normalize (factors);
1141    return factors;
1142  }
1143
1144  A= A (y + evaluation, y);
1145
1146  int liftBound;
1147  if (degree (A, y) == degree (LC (A, x)))
1148    liftBound= degree (LC (A, x)) + 1;
1149  else
1150    liftBound= degree (A, y) + 1 + degree (LC(A, x));
1151
1152  DEBOUTLN (cerr, "uniFactors= " << uniFactors);
1153  bool earlySuccess= false;
1154  CFList earlyFactors;
1155  TIMING_START (fac_hensel_lift);
1156  uniFactors= henselLiftAndEarly
1157               (A, earlySuccess, earlyFactors, degs, liftBound,
1158                uniFactors, info, evaluation);
1159  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1160  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
1161
1162  CanonicalForm MODl= power (y, liftBound);
1163  TIMING_START (fac_factor_recombination);
1164  if (!extension)
1165    factors= factorRecombination (uniFactors, A, MODl, degs);
1166  else
1167    factors= extFactorRecombination (uniFactors, A, MODl, info, degs,
1168                                     evaluation);
1169  TIMING_END_AND_PRINT (fac_factor_recombination,
1170                        "time for factor recombination: ");
1171  if (earlySuccess)
1172    factors= Union (earlyFactors, factors);
1173  else if (!earlySuccess && degs.getLength() == 1)
1174    factors= earlyFactors;
1175
1176  if (!extension)
1177  {
1178    for (CFListIterator i= factors; i.hasItem(); i++)
1179      i.getItem()= i.getItem() (y - evaluation, y);
1180  }
1181
1182  appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
1183                        swap, swap2, N);
1184  normalize (factors);
1185
1186  return factors;
1187}
1188
1189CFList
1190extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1191{
1192
1193  CanonicalForm A= F;
1194  Variable alpha= info.getAlpha();
1195  Variable beta= info.getBeta();
1196  int k= info.getGFDegree();
1197  char cGFName= info.getGFName();
1198  CanonicalForm delta= info.getDelta();
1199
1200  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1201  Variable x= Variable (1);
1202  CFList factors;
1203  if (!GF && alpha == x)  // we are in F_p
1204  {
1205    bool extension= true;
1206    int p= getCharacteristic();
1207    if (p*p < (1<<16)) // pass to GF if possible
1208    {
1209      setCharacteristic (getCharacteristic(), 2, 'Z');
1210      A= A.mapinto();
1211      ExtensionInfo info= ExtensionInfo (extension);
1212      factors= biFactorize (A, info);
1213
1214      Variable vBuf= rootOf (gf_mipo);
1215      setCharacteristic (getCharacteristic());
1216      for (CFListIterator j= factors; j.hasItem(); j++)
1217        j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
1218    }
1219    else // not able to pass to GF, pass to F_p(\alpha)
1220    {
1221      CanonicalForm mipo= randomIrredpoly (3, Variable (1));
1222      Variable v= rootOf (mipo);
1223      ExtensionInfo info= ExtensionInfo (v);
1224      factors= biFactorize (A, info);
1225    }
1226    return factors;
1227  }
1228  else if (!GF && (alpha != x)) // we are in F_p(\alpha)
1229  {
1230    if (k == 1) // need factorization over F_p
1231    {
1232      int extDeg= degree (getMipo (alpha));
1233      extDeg++;
1234      CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
1235      Variable v= rootOf (mipo);
1236      ExtensionInfo info= ExtensionInfo (v);
1237      factors= biFactorize (A, info);
1238    }
1239    else
1240    {
1241      if (beta == Variable (1))
1242      {
1243        Variable v= chooseExtension (A, alpha);
1244        CanonicalForm primElem, imPrimElem;
1245        bool primFail= false;
1246        Variable vBuf;
1247        primElem= primitiveElement (alpha, vBuf, primFail);
1248        ASSERT (!primFail, "failure in integer factorizer");
1249        if (primFail)
1250          ; //ERROR
1251        else
1252          imPrimElem= mapPrimElem (primElem, alpha, v);
1253
1254        CFList source, dest;
1255        CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem,
1256                                   source, dest);
1257        ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
1258        factors= biFactorize (bufA, info);
1259      }
1260      else
1261      {
1262        Variable v= chooseExtension (A, alpha);
1263        CanonicalForm primElem, imPrimElem;
1264        bool primFail= false;
1265        Variable vBuf;
1266        ASSERT (!primFail, "failure in integer factorizer");
1267        if (primFail)
1268          ; //ERROR
1269        else
1270          imPrimElem= mapPrimElem (delta, beta, v);
1271
1272        CFList source, dest;
1273        CanonicalForm bufA= mapDown (A, info, source, dest);
1274        source= CFList();
1275        dest= CFList();
1276        bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
1277        ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
1278        factors= biFactorize (bufA, info);
1279      }
1280    }
1281    return factors;
1282  }
1283  else // we are in GF (p^k)
1284  {
1285    int p= getCharacteristic();
1286    int extensionDeg= getGFDegree();
1287    bool extension= true;
1288    if (k == 1) // need factorization over F_p
1289    {
1290      extensionDeg++;
1291      if (ipower (p, extensionDeg) < (1<<16))
1292      // pass to GF(p^k+1)
1293      {
1294        setCharacteristic (p);
1295        Variable vBuf= rootOf (gf_mipo);
1296        A= GF2FalphaRep (A, vBuf);
1297        setCharacteristic (p, extensionDeg, 'Z');
1298        ExtensionInfo info= ExtensionInfo (extension);
1299        factors= biFactorize (A.mapinto(), info);
1300      }
1301      else // not able to pass to another GF, pass to F_p(\alpha)
1302      {
1303        setCharacteristic (p);
1304        Variable vBuf= rootOf (gf_mipo);
1305        A= GF2FalphaRep (A, vBuf);
1306        Variable v= chooseExtension (A, beta);
1307        ExtensionInfo info= ExtensionInfo (v, extension);
1308        factors= biFactorize (A, info);
1309      }
1310    }
1311    else // need factorization over GF (p^k)
1312    {
1313      if (ipower (p, 2*extensionDeg) < (1<<16))
1314      // pass to GF (p^2k)
1315      {
1316        setCharacteristic (p, 2*extensionDeg, 'Z');
1317        ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
1318        factors= biFactorize (GFMapUp (A, extensionDeg), info);
1319        setCharacteristic (p, extensionDeg, cGFName);
1320      }
1321      else // not able to pass to GF (p^2k), pass to F_p (\alpha)
1322      {
1323        setCharacteristic (p);
1324        Variable v1= rootOf (gf_mipo);
1325        A= GF2FalphaRep (A, v1);
1326        Variable v2= chooseExtension (A, v1);
1327        CanonicalForm primElem, imPrimElem;
1328        bool primFail= false;
1329        Variable vBuf;
1330        primElem= primitiveElement (v1, vBuf, primFail);
1331        ASSERT (!primFail, "failure in integer factorizer");
1332        if (primFail)
1333          ; //ERROR
1334        else
1335          imPrimElem= mapPrimElem (primElem, v1, v2);
1336
1337        CFList source, dest;
1338        CanonicalForm bufA= mapUp (A, v1, v2, primElem, imPrimElem,
1339                                     source, dest);
1340        ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
1341        factors= biFactorize (bufA, info);
1342        setCharacteristic (p, k, cGFName);
1343        for (CFListIterator i= factors; i.hasItem(); i++)
1344          i.getItem()= Falpha2GFRep (i.getItem());
1345      }
1346    }
1347    return factors;
1348  }
1349}
1350
1351#endif
1352/* HAVE_NTL */
1353
1354
Note: See TracBrowser for help on using the repository browser.