source: git/factory/facFqBivar.cc @ 6720a1

spielwiese
Last change on this file since 6720a1 was 6720a1, checked in by Hans Schoenemann <hannes@…>, 14 years ago
syntax fix for log2exp, part 1 git-svn-id: file:///usr/local/Singular/svn/trunk@12930 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 36.7 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
45const double log2exp= 1.442695041;
46
47CanonicalForm prodMod0 (const CFList& L, const CanonicalForm& M) 
48{
49  if (L.isEmpty())
50    return 1;
51  else if (L.length() == 1)
52    return mod (L.getFirst()(0, 1) , M);
53  else if (L.length() == 2)
54    return mod (L.getFirst()(0, 1)*L.getLast()(0, 1), M);
55  else
56  {
57    int l= L.length()/2;
58    CFListIterator i= L;
59    CFList tmp1, tmp2;
60    CanonicalForm buf1, buf2;
61    for (int j= 1; j <= l; j++, i++) 
62      tmp1.append (i.getItem());
63    tmp2= Difference (L, tmp1);
64    buf1= prodMod0 (tmp1, M);
65    buf2= prodMod0 (tmp2, M);
66    return mod (buf1*buf2, M);
67  }
68}
69
70CanonicalForm evalPoint (const CanonicalForm& F, CanonicalForm & eval, 
71                         const Variable& alpha, CFList& list, const bool& GF,
72                         bool& fail) 
73{
74  fail= false;
75  Variable x= Variable(2);
76  FFRandom genFF;
77  GFRandom genGF;
78  CanonicalForm random, mipo;
79  double bound;
80  int p= getCharacteristic ();
81  if (alpha != Variable(1)) 
82  {
83    mipo= getMipo (alpha);
84    int d= degree (mipo);
85    bound= ipower (p, d);
86  }
87  else if (GF) 
88  {
89    int d= getGFDegree();
90    bound= ipower (p, d);
91  }
92  else
93    bound= p;
94
95  do 
96  {
97    if (list.length() >= bound)
98    {
99      fail= true;
100      break;
101    }
102    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 != Variable(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 == Variable(1))
251    degMipoBeta= 1;
252  else if (!k && beta != Variable(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  while (T.length() >= 2*s) 
279  {
280    while (nosubset == false) 
281    {
282      if (T.length() == s) 
283      {
284        if (recombination) 
285        {
286          T.insert (LCBuf);
287          g= prodMod (T, M);
288          T.removeFirst();
289          g /= content(g);
290          g= g (y - eval, y);
291          g /= Lc (g);
292          appendTestMapDown (result, g, info, source, dest);
293          return result;
294        }
295        else
296        {
297          appendMapDown (result, F (y - eval, y), info, source, dest);
298          return result;
299        }
300      }
301      S= subset (v, s, TT, nosubset);
302      if (nosubset) break;
303      subsetDeg= subsetDegree (S);
304      // skip those combinations that are not possible
305      if (!degs.find (subsetDeg)) 
306        continue;
307      else 
308      { 
309        g= prodMod0 (S, M);
310        g= mod (g*LCBuf, M);
311        g /= content (g);
312        if (fdivides (LC (g), LCBuf))
313        {
314          S.insert (LCBuf);
315          g= prodMod (S, M);
316          S.removeFirst();
317          g /= content (g, Variable (1));
318          if (fdivides (g, buf)) 
319          {
320            buf2= g (y - eval, y);
321            appendTestMapDown (result, buf2, info, source, dest); 
322            if (!k && beta == Variable (1)) 
323            {
324              if (degree (buf2, alpha) < degMipoBeta)
325              { 
326                buf /= g;
327                LCBuf= LC (buf, Variable (1));
328                recombination= true;
329              }
330            }
331            else 
332            {
333              if (!isInExtension (buf2, delta, k))
334              {
335                buf /= g;
336                LCBuf= LC (buf, Variable (1));
337                recombination= true;
338              }
339            }
340            T= Difference (T, S);
341            // compute new possible degree pattern
342            bufDegs2= DegreePattern (T);
343            bufDegs1.intersect (bufDegs2);
344            bufDegs1.refine ();
345            if (T.length() < 2*s || T.length() == s || 
346                bufDegs1.getLength() == 1) 
347            {
348              if (recombination)
349              {
350                appendTestMapDown (result, buf (y - eval, y), info, source, 
351                                   dest);
352                return result; 
353              }
354              else
355              {
356                appendMapDown (result, F (y - eval, y), info, source, dest);
357                return result;               
358              }
359            }
360            TT= copy (T);
361            indexUpdate (v, s, T.length(), nosubset);
362            if (nosubset) break;       
363          }
364        }
365      }
366    }
367    s++;
368    if (T.length() < 2*s || T.length() == s) 
369    {
370      if (recombination)
371      {
372        appendTestMapDown (result, buf (y - eval, y), info, source, dest);
373        return result;
374      }
375      else
376      {
377        appendMapDown (result, F (y - eval, y), info, source, dest);
378        return result;
379      }
380    }
381    int v [T.length()];
382    for (int i= 0; i < T.length(); i++)
383      v[i]= 0;
384    nosubset= false;
385  }
386  if (T.length() < 2*s) 
387    appendMapDown (result, F (y - eval, y), info, source, dest);
388
389  return result; 
390}     
391
392/// naive factor recombination as decribed in "Factoring
393/// multivariate polynomials over a finite field" by L Bernardin.
394CFList
395factorRecombination (const CFList& factors, const CanonicalForm& F, 
396                     const CanonicalForm& M, const DegreePattern& degs)
397{
398  if (factors.length() == 0)
399    return CFList ();
400  if (degs.getLength() <= 1 || factors.length() == 1) 
401    return CFList(F);
402  DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, M) == F " << 
403            (LC (F, 1)*prodMod (factors, M) == F));
404  CFList T, S;
405
406  T= factors;
407  int s= 1;
408  CFList result;
409  CanonicalForm LCBuf= LC (F, Variable (1));
410  CanonicalForm g, buf= F;
411  int v [T.length()];
412  for (int i= 0; i < T.length(); i++)
413    v[i]= 0;
414  bool nosubset= false;
415  CFArray TT;
416  DegreePattern bufDegs1, bufDegs2;
417  bufDegs1= degs;
418  int subsetDeg;
419  TT= copy (factors);
420  bool recombination= false;
421  while (T.length() >= 2*s) 
422  {
423    while (nosubset == false) 
424    {
425      if (T.length() == s) 
426      {
427        if (recombination)
428        {
429          T.insert (LCBuf);
430          g= prodMod (T, M);
431          T.removeFirst();
432          result.append (g/content (g, Variable (1)));
433          return result;
434        }
435        else
436          return CFList (F);
437      }
438      S= subset (v, s, TT, nosubset);
439      if (nosubset) break;
440      subsetDeg= subsetDegree (S);
441      // skip those combinations that are not possible
442      if (!degs.find (subsetDeg)) 
443        continue;
444      else 
445      { 
446        g= prodMod0 (S, M);
447        g= mod (g*LCBuf, M);
448        g /= content (g);
449        if (fdivides (LC(g), LCBuf)) 
450        {
451          S.insert (LCBuf);
452          g= prodMod (S, M);
453          S.removeFirst();
454          g /= content (g, Variable (1));
455          if (fdivides (g, buf))
456          {
457            recombination= true;
458            result.append (g);
459            buf /= g;
460            LCBuf= LC (buf, Variable(1));
461            T= Difference (T, S);
462
463            // compute new possible degree pattern
464            bufDegs2= DegreePattern (T);
465            bufDegs1.intersect (bufDegs2);
466            bufDegs1.refine ();
467            if (T.length() < 2*s || T.length() == s || 
468                bufDegs1.getLength() == 1) 
469            {
470              if (recombination)
471              {
472                result.append (buf);
473                return result;
474              }
475              else
476                return CFList (F); 
477            }
478            TT= copy (T);
479            indexUpdate (v, s, T.length(), nosubset);
480            if (nosubset) break; 
481          }       
482        }
483      }
484    }
485    s++;
486    if (T.length() < 2*s || T.length() == s) 
487    {
488      if (recombination)
489      {
490        result.append (buf);
491        return result;
492      }
493      else
494        return CFList (F);
495    }
496    int v [T.length()];
497    for (int i= 0; i < T.length(); i++)
498      v[i]= 0;
499    nosubset= false;
500  }
501  if (T.length() < 2*s) 
502    result.append (F);
503
504  return result; 
505} 
506
507Variable chooseExtension (const CanonicalForm & A, const Variable & alpha) 
508{
509  int p= getCharacteristic();
510  ZZ NTLp= to_ZZ (p);
511  ZZ_p::init (NTLp);
512  ZZ_pX NTLirredpoly; 
513  int d= degree (A); 
514  int m= 1;
515  if (alpha != Variable (1)) 
516    m= degree (getMipo (alpha));   
517  int i= (int) floor((double) d/(double) m) - 1;
518  if (i < 2) 
519    i= 2;
520  if (i > 8)
521    i= i/4;
522  BuildIrred (NTLirredpoly, i*m);
523  Variable x= A.mvar();
524  CanonicalForm newMipo= convertNTLZZpX2CF (NTLirredpoly, x);
525  return rootOf (newMipo);
526}
527
528CFList
529earlyFactorDetection (CanonicalForm& F, CFList& factors,int& adaptedLiftBound, 
530                      DegreePattern& degs, bool& success, int deg)
531{
532  DegreePattern bufDegs1= degs;
533  DegreePattern bufDegs2;
534  CFList result;
535  CFList T= factors;
536  CanonicalForm buf= F;
537  CanonicalForm LCBuf= LC (buf, Variable (1));
538  CanonicalForm g;
539  CanonicalForm M= power (F.mvar(), deg);
540  adaptedLiftBound= 0;
541  int d;
542  if (degree (LCBuf) == degree (F))
543    d= degree (F);
544  else
545    d= degree (F) + degree (LCBuf);
546  for (CFListIterator i= factors; i.hasItem(); i++) 
547  {
548    if (!bufDegs1.find (degree (i.getItem(), 1)))
549      continue;
550    else
551    {
552      g= i.getItem() (0, 1);
553      g *= LCBuf;
554      g= mod (g, M);
555      if (fdivides (LC (g), LCBuf))
556      {
557        g= mulMod2 (i.getItem(), LCBuf, M);
558        g /= content (g, Variable (1));
559        if (fdivides (g, buf))
560        {
561          result.append (g);
562          buf /= g;
563          d -= degree (g) + degree (LC (g, Variable (1)));
564          LCBuf= LC (buf, Variable (1)); 
565          T= Difference (T, CFList (i.getItem()));
566
567          // compute new possible degree pattern
568          bufDegs2= DegreePattern (T);
569          bufDegs1.intersect (bufDegs2);
570          bufDegs1.refine (); 
571          if (bufDegs1.getLength() <= 1)
572          {
573            result.append (buf);
574            break;
575          }
576        }
577      }
578    }
579  }
580  adaptedLiftBound= d + 1;
581  if (d < deg)
582  {
583    factors= T;
584    degs= bufDegs1;
585    F= buf;
586    success= true;
587  }
588  if (bufDegs1.getLength() <= 1)
589    degs= bufDegs1;
590  return result;
591}
592
593CFList
594extEarlyFactorDetection (CanonicalForm& F, CFList& factors, 
595                         int& adaptedLiftBound, DegreePattern& degs, 
596                         bool& success, const ExtensionInfo& info,
597                         const CanonicalForm& eval, int deg)
598{
599  Variable alpha= info.getAlpha();
600  Variable beta= info.getBeta();
601  CanonicalForm gamma= info.getGamma();
602  CanonicalForm delta= info.getDelta();
603  int k= info.getGFDegree();
604  DegreePattern bufDegs1= degs, bufDegs2;
605  CFList result; 
606  CFList T= factors;
607  Variable y= F.mvar();
608  CanonicalForm buf= F, LCBuf= LC (buf, Variable (1)), g, buf2;
609  CanonicalForm M= power (y, deg);
610  adaptedLiftBound= 0;
611  int d;
612  if (degree (F) == degree (LCBuf))
613    d= degree (F);
614  else
615    d= degree (F) + degree (LCBuf);
616  CFList source, dest;
617  int degMipoBeta;
618  if (!k && beta == Variable(1))
619    degMipoBeta= 1;
620  else if (!k && beta != Variable(1))
621    degMipoBeta= degree (getMipo (beta));
622  for (CFListIterator i= factors; i.hasItem(); i++) 
623  {
624    if (!bufDegs1.find (degree (i.getItem(), 1)))
625      continue;
626    else
627    {
628      g= i.getItem() (0, 1);
629      g *= LCBuf;
630      g= mod (g, M);
631      if (fdivides (LC (g), LCBuf))
632      {
633        g= mulMod2 (i.getItem(), LCBuf, M);
634        g /= content (g, Variable (1));
635        if (fdivides (g, buf))
636        {
637          buf2= g (y - eval, y);
638          buf2 /= Lc (buf2);
639           
640          if (!k && beta == Variable (1)) 
641          {
642            if (degree (buf2, alpha) < degMipoBeta)
643            {
644              appendTestMapDown (result, buf2, info, source, dest);
645              buf /= g;
646              d -= degree (g) + degree (LC (g, Variable (1)));
647              LCBuf= LC (buf, Variable (1));           
648            }
649          }
650          else 
651          {
652            if (!isInExtension (buf2, delta, k))
653            {
654              appendTestMapDown (result, buf2, info, source, dest);
655              buf /= g;
656              d -= degree (g) + degree (LC (g, Variable (1)));
657              LCBuf= LC (buf, Variable (1));
658            }
659          } 
660          T= Difference (T, CFList (i.getItem()));
661
662          // compute new possible degree pattern
663          bufDegs2= DegreePattern (T);
664          bufDegs1.intersect (bufDegs2);
665          bufDegs1.refine ();
666          if (bufDegs1.getLength() <= 1)
667          {
668            buf= buf (y - eval, y);
669            buf /= Lc (buf);
670            appendMapDown (result, buf, info, source, dest);
671            break;
672          }
673        }
674      }
675    }
676  }
677  adaptedLiftBound= d + 1;
678  if (adaptedLiftBound < deg)
679  {
680    success= true;
681    factors= T;
682    degs= bufDegs1;
683    F= buf;
684  }
685  if (bufDegs1.getLength() <= 1)
686    degs= bufDegs1;
687
688  return result;
689}
690
691CFList
692henselLiftAndEarly (CanonicalForm& A, bool& earlySuccess, CFList& 
693                    earlyFactors, DegreePattern& degs, int& liftBound,
694                    const CFList& uniFactors, const ExtensionInfo& info, 
695                    const CanonicalForm& eval)
696{
697  Variable alpha= info.getAlpha();
698  Variable beta= info.getBeta();
699  CanonicalForm gamma= info.getGamma();
700  CanonicalForm delta= info.getDelta();
701  int k= info.getGFDegree();
702  bool extension= info.isInExtension();
703
704  Variable x= Variable (1);
705  Variable y= Variable (2);
706  CFArray Pi;
707  CFList diophant;
708  CFList bufUniFactors= uniFactors;
709  bufUniFactors.insert (LC (A, x));
710  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
711  earlySuccess= false;
712  int newLiftBound= 0;
713  int smallFactorDeg= 11; //this is a tunable parameter
714  if (smallFactorDeg >= liftBound)
715    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
716  else if (smallFactorDeg >= degree (A, y) + 1)
717  {
718    henselLift12 (A, bufUniFactors, degree (A, y) + 1, Pi, diophant, M);
719    if (!extension)
720      earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
721                     degs, earlySuccess, degree (A, y) + 1);
722    else
723      earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
724                     newLiftBound, degs, earlySuccess, info, eval, 
725                     degree (A, y) + 1);
726    if (degs.getLength() > 1 && !earlySuccess)
727    {
728      if (newLiftBound > degree (A, y) + 1)
729      {
730        liftBound= newLiftBound;
731        bufUniFactors.insert (LC(A, x));
732        henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound, 
733                            Pi, diophant, M);
734      }
735    }
736    else if (earlySuccess)
737      liftBound= newLiftBound;
738  }
739  else if (smallFactorDeg < degree (A, y) + 1)
740  { 
741    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
742    if (!extension)
743      earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
744                                           degs, earlySuccess,
745                                           smallFactorDeg);
746    else
747      earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
748                     newLiftBound, degs, earlySuccess,
749                     info, eval, smallFactorDeg);
750    if (degs.getLength() > 1 && !earlySuccess)
751    {
752      bufUniFactors.insert (LC (A, x));
753      henselLiftResume12 (A, bufUniFactors, smallFactorDeg, degree (A, y) 
754                          + 1, Pi, diophant, M);
755      if (!extension)
756        earlyFactors= earlyFactorDetection (A, bufUniFactors, newLiftBound,
757                       degs, earlySuccess, degree (A, y) + 1);
758      else
759        earlyFactors= extEarlyFactorDetection (A, bufUniFactors,
760                       newLiftBound, degs, earlySuccess,
761                       info, eval, degree(A,y) + 1);
762      if (degs.getLength() > 1 && !earlySuccess)
763      {
764        if (newLiftBound > degree (A, y) + 1)
765        {
766          if (newLiftBound < newLiftBound)
767            liftBound= newLiftBound;
768          bufUniFactors.insert (LC(A, x));
769          henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
770                              Pi, diophant, M);
771        }
772      }
773      else if (earlySuccess)
774        liftBound= newLiftBound;
775    }
776    else if (earlySuccess)
777      liftBound= newLiftBound;
778  }
779  if (newLiftBound > 0)
780    liftBound= newLiftBound;
781  return bufUniFactors;
782}
783
784CFList
785extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info);
786
787/// bivariate factorization over finite fields as decribed in "Factoring
788/// multivariate polynomials over a finite field" by L Bernardin.
789CFList
790biFactorize (const CanonicalForm& F, const ExtensionInfo& info)
791{
792  if (F.inCoeffDomain())
793    return CFList(F);
794 
795  CanonicalForm A= F;
796  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
797
798  Variable alpha= info.getAlpha();
799  Variable beta= info.getBeta();
800  CanonicalForm gamma= info.getGamma();
801  CanonicalForm delta= info.getDelta();
802  int k= info.getGFDegree();
803  bool extension= info.isInExtension();
804  if (A.isUnivariate()) 
805  {
806    if (extension == false) 
807      return uniFactorizer (F, alpha, GF);
808    else 
809    {
810      CFList source, dest;
811      A= mapDown (A, info, source, dest);
812      return uniFactorizer (A, beta, GF);
813    }
814  }
815
816  CFMap N;
817  A= compress (A, N);
818  Variable y= A.mvar();
819
820  if (y.level() > 2) return CFList (F);
821  Variable x= Variable (1);
822
823  //remove and factorize content
824  CanonicalForm contentAx= content (A, x);
825  CanonicalForm contentAy= content (A);
826
827  A= A/(contentAx*contentAy);
828  CFList contentAxFactors, contentAyFactors;
829
830  if (!extension) 
831  { 
832    contentAxFactors= uniFactorizer (contentAx, alpha, GF);
833    contentAyFactors= uniFactorizer (contentAy, alpha, GF);
834  }
835
836  //trivial case
837  CFList factors;
838  if (A.inCoeffDomain()) 
839  {
840    append (factors, contentAxFactors);
841    append (factors, contentAyFactors);
842    decompress (factors, N); 
843    return factors;
844  }
845  else if (A.isUnivariate()) 
846  {
847    factors= uniFactorizer (A, alpha, GF);
848    append (factors, contentAxFactors);
849    append (factors, contentAyFactors);
850    decompress (factors, N); 
851    return factors;
852  }
853
854  // check derivatives
855  bool derivXZero= false;
856  CanonicalForm derivX= deriv (A, x);
857  CanonicalForm gcdDerivX;
858  if (derivX.isZero())
859    derivXZero= true;
860  else
861  {
862    if (GF)
863      gcdDerivX= GCD_GF (A, derivX);
864    else if (alpha == x)
865      gcdDerivX= GCD_small_p (A, derivX);
866    else 
867      gcdDerivX= GCD_Fp_extension (A, derivX, alpha);
868    if (degree (gcdDerivX) > 0)
869    {
870      CanonicalForm g= A/gcdDerivX;
871      CFList factorsG= 
872      Union (biFactorize (g, info), biFactorize (gcdDerivX, info));
873      append (factorsG, contentAxFactors);
874      append (factorsG, contentAyFactors);
875      decompress (factorsG, N);
876      normalize (factors);
877      return factorsG;   
878    }
879  }
880  bool derivYZero= false;
881  CanonicalForm derivY= deriv (A, y);
882  CanonicalForm gcdDerivY;
883  if (derivY.isZero())
884    derivYZero= true;
885  else
886  {
887    if (GF)
888      gcdDerivY= GCD_GF (A, derivY);
889    else if (alpha == x)
890      gcdDerivY= GCD_small_p (A, derivY);
891    else 
892      gcdDerivY= GCD_Fp_extension (A, derivY, alpha);
893    if (degree (gcdDerivY) > 0)
894    {
895      CanonicalForm g= A/gcdDerivY;
896      CFList factorsG= 
897      Union (biFactorize (g, info), biFactorize (gcdDerivY, info));
898      append (factorsG, contentAxFactors);
899      append (factorsG, contentAyFactors);
900      decompress (factorsG, N);
901      normalize (factors);
902      return factorsG;   
903    }
904  }
905  //main variable is chosen s.t. the degree in x is minimal
906  bool swap= false;
907  if ((degree (A) < degree (A, x)) || derivXZero) 
908  {
909    if (!derivYZero) 
910    { 
911      A= swapvar (A, y, x);
912      swap= derivXZero; 
913      derivXZero= derivYZero;
914      derivYZero= swap;
915      swap= true;
916    }
917  }
918
919  bool fail;
920  CanonicalForm Aeval, evaluation, bufAeval, bufEvaluation, buf;
921  CFList uniFactors, list, bufUniFactors;
922  DegreePattern degs;
923  DegreePattern bufDegs;
924
925  bool fail2= false;
926  CanonicalForm Aeval2, evaluation2, bufAeval2, bufEvaluation2;
927  CFList bufUniFactors2, list2, uniFactors2;
928  DegreePattern degs2;
929  DegreePattern bufDegs2;
930  bool swap2= false;
931
932  // several univariate factorizations to obtain more information about the
933  // degree pattern therefore usually less combinations have to be tried during
934  // the recombination process
935  int factorNums= 5;
936  double logarithm= (double) ilog2 (totaldegree (A));
937  logarithm /= log2exp;
938  logarithm= ceil (logarithm);
939  if (factorNums < (int) logarithm)
940    factorNums= (int) logarithm;
941  for (int i= 0; i < factorNums; i++) 
942  {
943    bufAeval= A; 
944    bufEvaluation= evalPoint (A, bufAeval, alpha, list, GF, fail);
945    if (!derivXZero && !fail2)
946    {
947      buf= swapvar (A, x, y);
948      bufAeval2= buf;
949      bufEvaluation2= evalPoint (buf, bufAeval2, alpha, list2, GF, fail2);
950    }
951    // first try to change main variable if there is no valid evaluation point
952    if (fail && (i == 0)) 
953    {
954      if (!derivXZero)
955        bufEvaluation= evalPoint (buf, bufAeval, alpha, list, GF, fail);
956      else
957        fail= true;
958
959      if (!fail) 
960      { 
961        A= buf;
962        swap2= true;
963      }
964    }
965
966    // if there is no valid evaluation point pass to a field extension
967    if (fail && (i == 0)) 
968    {
969      factors= extBiFactorize (A, info);
970      appendSwapDecompress (factors, contentAxFactors, contentAyFactors,
971                            swap, swap2, N);
972      normalize (factors);
973      return factors;
974    }
975
976    // there is at least one valid evaluation point
977    // but we do not compute more univariate factorization over an extension
978    if (fail && (i != 0)) 
979      break;
980
981    // univariate factorization
982    TIMING_START (fac_uni_factorizer);
983    bufUniFactors= uniFactorizer (bufAeval, alpha, GF);
984    TIMING_END_AND_PRINT (fac_uni_factorizer, 
985                          "time for univariate factorization: ");
986    DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " << 
987              prod (bufUniFactors)*Lc (bufAeval) == bufAeval);
988
989    if (!derivXZero && !fail2)
990    {
991      TIMING_START (fac_uni_factorizer);
992      bufUniFactors2= uniFactorizer (bufAeval2, alpha, GF);
993      TIMING_END_AND_PRINT (fac_uni_factorizer, 
994                            "time for univariate factorization in y: ");
995      DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " << 
996                prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2);
997    }
998
999    if (bufUniFactors.length() == 1 || 
1000        (!fail2 && !derivXZero && (bufUniFactors2.length() == 1))) 
1001    {
1002      if (extension)
1003      {
1004        CFList source, dest;
1005        ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k, 
1006                             info.getGFName(), info.isInExtension());
1007        appendMapDown (factors, A, info2, source, dest);
1008      }
1009      else 
1010        factors.append (A);
1011
1012      appendSwapDecompress (factors, contentAxFactors, contentAyFactors, 
1013                            swap, swap2, N);
1014
1015      normalize (factors);
1016      return factors;
1017    }
1018
1019    // degree analysis   
1020    bufDegs = DegreePattern (bufUniFactors); 
1021    if (!derivXZero && !fail2)
1022      bufDegs2= DegreePattern (bufUniFactors2);
1023
1024    if (i == 0) 
1025    {
1026      Aeval= bufAeval;
1027      evaluation= bufEvaluation;
1028      uniFactors= bufUniFactors;
1029      degs= bufDegs;
1030      if (!derivXZero && !fail2)
1031      {
1032        Aeval2= bufAeval2;
1033        evaluation2= bufEvaluation2;
1034        uniFactors2= bufUniFactors2;
1035        degs2= bufDegs2;
1036      }
1037    }
1038    else 
1039    {
1040      degs.intersect (bufDegs);
1041      if (!derivXZero && !fail2) 
1042      {
1043        degs2.intersect (bufDegs2);
1044        if (bufUniFactors2.length() < uniFactors2.length()) 
1045        { 
1046          uniFactors2= bufUniFactors2;
1047          Aeval2= bufAeval2;
1048          evaluation2= bufEvaluation2;
1049        }
1050      }   
1051      if (bufUniFactors.length() < uniFactors.length()) 
1052      { 
1053        uniFactors= bufUniFactors;
1054        Aeval= bufAeval;
1055        evaluation= bufEvaluation;
1056      }
1057    }
1058    list.append (bufEvaluation);
1059    if (!derivXZero && !fail2)
1060      list2.append (bufEvaluation2);
1061  }
1062 
1063  if (!derivXZero && !fail2)
1064  {
1065    if (uniFactors.length() > uniFactors2.length() || 
1066        (uniFactors.length() == uniFactors2.length() 
1067         && degs.getLength() > degs2.getLength())) 
1068    {
1069      degs= degs2;
1070      uniFactors= uniFactors2;
1071      evaluation= evaluation2;
1072      Aeval= Aeval2;
1073      A= buf;
1074      swap2= true;
1075    }
1076  }
1077 
1078  if (degs.getLength() == 1) // A is irreducible
1079  {
1080    if (extension) 
1081    {
1082      CFList source, dest;
1083      ExtensionInfo info2= ExtensionInfo (beta, alpha, delta, gamma, k, 
1084                           info.getGFName(), info.isInExtension());
1085      appendMapDown (factors, A, info2, source, dest);
1086    }
1087    else 
1088      factors.append (A);
1089    appendSwapDecompress (factors, contentAxFactors, contentAyFactors, 
1090                            swap, swap2, N);
1091    normalize (factors);
1092    return factors;
1093  }
1094
1095  A= A (y + evaluation, y);
1096
1097  int liftBound;
1098  if (degree (A, y) == degree (LC (A, x)))
1099    liftBound= degree (LC (A, x)) + 1;
1100  else
1101    liftBound= degree (A, y) + 1 + degree (LC(A, x)); 
1102
1103  DEBOUTLN (cerr, "uniFactors= " << uniFactors);
1104  bool earlySuccess= false;
1105  CFList earlyFactors; 
1106  TIMING_START (fac_hensel_lift);
1107  uniFactors= henselLiftAndEarly
1108               (A, earlySuccess, earlyFactors, degs, liftBound,
1109                uniFactors, info, evaluation); 
1110  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1111  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
1112
1113  CanonicalForm MODl= power (y, liftBound);
1114  TIMING_START (fac_factor_recombination);
1115  if (!extension)
1116    factors= factorRecombination (uniFactors, A, MODl, degs);
1117  else 
1118    factors= extFactorRecombination (uniFactors, A, MODl, info, degs, 
1119                                     evaluation);
1120  TIMING_END_AND_PRINT (fac_factor_recombination,
1121                        "time for factor recombination: ");
1122  if (earlySuccess)
1123    factors= Union (earlyFactors, factors);
1124  else if (!earlySuccess && degs.getLength() == 1)
1125    factors= earlyFactors;
1126 
1127  if (!extension)
1128  {
1129    for (CFListIterator i= factors; i.hasItem(); i++)
1130      i.getItem()= i.getItem() (y - evaluation, y);
1131  }
1132
1133  appendSwapDecompress (factors, contentAxFactors, contentAyFactors, 
1134                        swap, swap2, N);
1135  normalize (factors); 
1136
1137  return factors;
1138}
1139
1140CFList
1141extBiFactorize (const CanonicalForm& F, const ExtensionInfo& info) 
1142{
1143 
1144  CanonicalForm A= F;
1145  Variable alpha= info.getAlpha();
1146  Variable beta= info.getBeta();
1147  int k= info.getGFDegree();
1148  char cGFName= info.getGFName();
1149  CanonicalForm delta= info.getDelta();
1150
1151  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1152  Variable x= Variable (1); 
1153  CFList factors;
1154  if (!GF && alpha == x)  // we are in F_p
1155  {
1156    bool extension= true;         
1157    int p= getCharacteristic();
1158    if (p*p < (1<<16)) // pass to GF if possible
1159    { 
1160      setCharacteristic (getCharacteristic(), 2, 'Z');
1161      A= A.mapinto();
1162      ExtensionInfo info= ExtensionInfo (extension);
1163      factors= biFactorize (A, info);
1164
1165      Variable vBuf= rootOf (gf_mipo);
1166      setCharacteristic (getCharacteristic());
1167      for (CFListIterator j= factors; j.hasItem(); j++) 
1168        j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
1169    }
1170    else // not able to pass to GF, pass to F_p(\alpha)
1171    { 
1172      CanonicalForm mipo= randomIrredpoly (3, Variable (1));
1173      Variable v= rootOf (mipo);
1174      ExtensionInfo info= ExtensionInfo (v);
1175      factors= biFactorize (A, info);
1176    } 
1177    return factors;
1178  }
1179  else if (!GF && (alpha != x)) // we are in F_p(\alpha)
1180  { 
1181    if (k == 1) // need factorization over F_p
1182    {
1183      int extDeg= degree (getMipo (alpha));
1184      extDeg++; 
1185      CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
1186      Variable v= rootOf (mipo);
1187      ExtensionInfo info= ExtensionInfo (v);
1188      factors= biFactorize (A, info);
1189    }
1190    else 
1191    {
1192      if (beta == Variable (1))
1193      {
1194        Variable v= chooseExtension (A, alpha);
1195        CanonicalForm primElem, imPrimElem;
1196        bool primFail= false;
1197        Variable vBuf;
1198        primElem= primitiveElement (alpha, vBuf, primFail);
1199        ASSERT (!primFail, "failure in integer factorizer");
1200        if (primFail)
1201          ; //ERROR
1202        else
1203          imPrimElem= mapPrimElem (primElem, vBuf, v);
1204
1205        CFList source, dest;
1206        CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem, 
1207                                   source, dest);
1208        ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
1209        factors= biFactorize (bufA, info);
1210      }
1211      else
1212      {
1213        Variable v= chooseExtension (A, alpha);
1214        CanonicalForm primElem, imPrimElem;
1215        bool primFail= false;
1216        Variable vBuf;
1217        ASSERT (!primFail, "failure in integer factorizer");
1218        if (primFail)
1219          ; //ERROR
1220        else
1221          imPrimElem= mapPrimElem (delta, beta, v); //oder mapPrimElem (primElem, vBuf, v);
1222
1223        CFList source, dest;
1224        CanonicalForm bufA= mapDown (A, info, source, dest);
1225        source= CFList();
1226        dest= CFList();
1227        bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
1228        ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
1229        factors= biFactorize (bufA, info);
1230      }
1231    }
1232    return factors;
1233  }
1234  else // we are in GF (p^k)
1235  { 
1236    int p= getCharacteristic();
1237    int extensionDeg= getGFDegree();
1238    bool extension= true;
1239    if (k == 1) // need factorization over F_p
1240    {
1241      extensionDeg++;
1242      if (ipower (p, extensionDeg) < (1<<16)) 
1243      // pass to GF(p^k+1)
1244      {
1245        setCharacteristic (p);
1246        Variable vBuf= rootOf (gf_mipo);
1247        A= GF2FalphaRep (A, vBuf);
1248        setCharacteristic (p, extensionDeg, 'Z');
1249        ExtensionInfo info= ExtensionInfo (extension);
1250        factors= biFactorize (A.mapinto(), info);
1251      }
1252      else // not able to pass to another GF, pass to F_p(\alpha)
1253      { 
1254        setCharacteristic (p);
1255        Variable vBuf= rootOf (gf_mipo);
1256        A= GF2FalphaRep (A, vBuf);
1257        Variable v= chooseExtension (A, beta); 
1258        ExtensionInfo info= ExtensionInfo (v, extension);
1259        factors= biFactorize (A, info);
1260      }
1261    }
1262    else // need factorization over GF (p^k)
1263    {
1264      if (ipower (p, 2*extensionDeg) < (1<<16)) 
1265      // pass to GF (p^2k)
1266      { 
1267        setCharacteristic (p, 2*extensionDeg, 'Z');
1268        ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
1269        factors= biFactorize (GFMapUp (A, extensionDeg), info); 
1270        setCharacteristic (p, extensionDeg, cGFName);
1271      }
1272      else // not able to pass to GF (p^2k), pass to F_p (\alpha)
1273      { 
1274        setCharacteristic (p);
1275        Variable v1= rootOf (gf_mipo); 
1276        A= GF2FalphaRep (A, v1);
1277        Variable v2= chooseExtension (A, v1); 
1278        CanonicalForm primElem, imPrimElem;
1279        bool primFail= false;
1280        Variable vBuf;
1281        primElem= primitiveElement (v1, vBuf, primFail);
1282        ASSERT (!primFail, "failure in integer factorizer");
1283        if (primFail)
1284          ; //ERROR
1285        else
1286          imPrimElem= mapPrimElem (primElem, vBuf, v2);
1287
1288        CFList source, dest;
1289        CanonicalForm bufA= mapUp (A, v1, v2, primElem, imPrimElem, 
1290                                     source, dest);
1291        ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
1292        factors= biFactorize (bufA, info); 
1293        setCharacteristic (p, k, cGFName);
1294        for (CFListIterator i= factors; i.hasItem(); i++) 
1295          i.getItem()= Falpha2GFRep (i.getItem());
1296      }
1297    }
1298    return factors;
1299  }
1300}
1301
1302#endif
1303/* HAVE_NTL */
1304
1305
Note: See TracBrowser for help on using the repository browser.