source: git/libfac/factor/MVMultiHensel.cc @ 18500b

jengelh-datetimespielwiese
Last change on this file since 18500b was 18500b, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: moved libfac back
  • Property mode set to 100644
File size: 14.1 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2// emacs edit mode for this file is -*- C++ -*-
3///////////////////////////////////////////////////////////////////////////////
4// FACTORY - Includes
5#include <factory.h>
6#ifndef NOSTREAMIO
7#ifdef HAVE_IOSTREAM
8#include <iostream>
9#define OSTREAM std::ostream
10#define ISTREAM std::istream
11#define CERR std::cerr
12#define CIN std::cin
13#elif defined(HAVE_IOSTREAM_H)
14#include <iostream.h>
15#define OSTREAM ostream
16#define ISTREAM istream
17#define CERR cerr
18#define CIN cin
19#endif
20#endif
21// Factor - Includes
22#include "tmpl_inst.h"
23#include "helpstuff.h"
24// some CC's need this:
25#include "MVMultiHensel.h"
26
27#ifndef NOSTREAMIO
28void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
29#endif
30
31extern int libfac_interruptflag;
32
33#ifdef HENSELDEBUG
34#  define DEBUGOUTPUT
35#else
36#  undef DEBUGOUTPUT
37#endif
38
39#include <libfac/factor/debug.h>
40#include "interrupt.h"
41#include "timing.h"
42
43///////////////////////////////////////////////////////////////
44// some class definition needed in MVMultiHensel
45///////////////////////////////////////////////////////////////
46typedef bool Boolean;
47
48class DiophantForm
49{
50 public:
51   CanonicalForm One;
52   CanonicalForm Two;
53   inline DiophantForm& operator=( const DiophantForm&  value )
54   {
55     if ( this != &value )
56     {
57       One = value.One;
58       Two = value.Two;
59     }
60     return *this;
61   }
62};
63
64// We remember an already calculated value; simple class for RememberArray
65class RememberForm
66{
67public:
68  inline RememberForm operator=( CanonicalForm & value )
69  {
70    this->calculated = true;
71    this->poly = value;
72    return *this;
73  }
74  RememberForm() : poly(0), calculated(false) {}
75  Boolean calculated;
76  CanonicalForm poly;
77};
78
79// Array to remember already calculated values; used for the diophantine
80// equation s*f + t*g = x^i
81class RememberArray
82{
83public:
84// operations performed on arrays
85  RememberArray( int sz )
86  {
87    size = sz;
88    ia = new RememberForm[size];
89//    assert( ia != 0 ); // test if we got the memory
90    init( sz );
91  }
92  ~RememberArray(){ delete [] ia; }
93  inline RememberForm& operator[]( int index )
94  {
95    return ia[index];
96  }
97  bool checksize(int i) {return i<size;}
98protected:
99  void init( int )
100  {
101    for ( int ix=0; ix < size; ++ix)
102    {
103      ia[ix].calculated=false;
104      ia[ix].poly=0;
105    }
106  }
107// internal data representation
108  int size;
109  RememberForm *ia;
110};
111
112///////////////////////////////////////////////////////////////
113// Solve the Diophantine equation: ( levelU == mainvar )     //
114//            s*F1 + t*F2 = (mainvar)^i                      //
115// Returns s and t.                                          //
116///////////////////////////////////////////////////////////////
117static DiophantForm
118diophant( int levelU , const CanonicalForm & F1 , const CanonicalForm & F2 ,
119          int i , RememberArray & A, RememberArray & B,
120          const CanonicalForm &alpha )
121{
122  DiophantForm Retvalue;
123  CanonicalForm s,t,q,r;
124  Variable x(levelU);
125
126  DEBOUT(CERR, "diophant: called with: ", F1);
127  DEBOUT(CERR, "  ", F2); DEBOUTLN(CERR, "  ", i);
128
129  // Did we solve the diophantine equation yet?
130  // If so, return the calculated values
131  if (A.checksize(i) && A[i].calculated && B[i].calculated )
132  {
133    Retvalue.One=A[i].poly;
134    Retvalue.Two=B[i].poly;
135    return Retvalue;
136  }
137
138  // Degrees ok? degree(F1,mainvar) + degree(F2,mainvar) <= i ?
139  if ( (degree(F1,levelU) + degree(F2,levelU) ) <= i )
140  {
141    if (!interrupt_handle()) factoryError("libfac: diophant ERROR: degree too large!"); /*  (%d + %d <= %d)",degree(F1,levelU), degree(F2,levelU), i);*/
142    Retvalue.One=F1;Retvalue.Two=F2;
143    return Retvalue;
144  }
145
146  if ( i == 0 )
147  { // call the extended gcd
148    r=extgcd(F1,F2,s,t);
149    // check if gcd(F1,F2) <> 1 , i.e. F1 and F2 are not relatively prime
150    if ( ! r.isOne() )
151    {
152      if (r.degree()<1) // some constant != 1 ?
153      {
154        Retvalue.One=s/r;Retvalue.Two=t/r;
155        return Retvalue;
156      }
157      else if (alpha!=0)
158      {
159        Variable Alpha=alpha.mvar();
160        if (r.mvar()==Alpha)   // from field extension ?
161        {
162          Variable X=rootOf(alpha);
163          r=replacevar(r,Alpha,X);
164          s=replacevar(s,Alpha,X);
165          t=replacevar(t,Alpha,X);
166          s/=r;
167          t/=r;
168          s=replacevar(s,X,Alpha);
169          t=replacevar(t,X,Alpha);
170          Retvalue.One=s; Retvalue.Two=t;
171          return Retvalue;
172        }
173      } 
174      if (!interrupt_handle())
175        factoryError("libfac: diophant ERROR: F1 and F2 are not relatively prime! ");
176      Retvalue.One=s/r;Retvalue.Two=t/r;
177      return Retvalue;
178    }
179    Retvalue.One = s; Retvalue.Two = t;
180  }
181  else
182  { // recursively call diophant
183    Retvalue=diophant(levelU,F1,F2,i-1,A,B,alpha);
184    Retvalue.One *= x; // createVar(levelU,1);
185    Retvalue.Two *= x; // createVar(levelU,1);
186
187    if (interrupt_handle()) return Retvalue;
188
189    // Check degrees.
190
191    if ( degree(Retvalue.One,levelU) > degree(F2,levelU) )
192    {
193      // Make degree(Retvalue.one,mainvar) < degree(F2,mainvar)
194      divrem(Retvalue.One,F2,q,r);
195      Retvalue.One = r; Retvalue.Two += F1*q;
196    }
197    else
198    {
199      if ( degree(Retvalue.Two,levelU) >= degree(F1,levelU) )
200      {
201        // Make degree(Retvalue.Two,mainvar) <= degree(F1,mainvar)
202        divrem(Retvalue.Two,F1,q,r);
203        Retvalue.One += F2*q; Retvalue.Two = r;
204      }
205    }
206
207  }
208  if (A.checksize(i))
209  {
210    A[i].poly = Retvalue.One ;
211    B[i].poly = Retvalue.Two ;
212    A[i].calculated = true ; B[i].calculated = true ;
213  }
214  DEBOUT(CERR, "diophant: Returnvalue is: ", Retvalue.One);
215  DEBOUTLN(CERR, "  ", Retvalue.Two);
216
217  return  Retvalue;
218}
219
220///////////////////////////////////////////////////////////////
221// A more efficient way to solve s*F1 + t*F2 = W             //
222// as in Wang and Rothschild [Wang&Roth75].                  //
223///////////////////////////////////////////////////////////////
224static CanonicalForm
225make_delta( int levelU, const CanonicalForm & W,
226            const CanonicalForm & F1, const CanonicalForm & F2,
227            RememberArray & A, RememberArray & B,
228            const CanonicalForm &alpha)
229{
230  CanonicalForm Retvalue;
231  DiophantForm intermediate;
232
233  DEBOUT(CERR, "make_delta: W= ", W);
234  DEBOUTLN(CERR, "  degree(W,levelU)= ", degree(W,levelU) );
235
236  if ( levelU == level(W) ) // same level, good
237  {
238    for ( CFIterator i=W; i.hasTerms(); i++)
239    {
240      intermediate=diophant(levelU,F1,F2,i.exp(),A,B,alpha);
241      Retvalue += intermediate.One * i.coeff();
242
243      if (interrupt_handle()) return Retvalue;
244    }
245  }
246  else // level(W) < levelU ; i.e. degree(w,levelU) == 0
247  {
248    intermediate=diophant(levelU,F1,F2,0,A,B,alpha);
249    Retvalue = W * intermediate.One;
250  }
251  DEBOUTLN(CERR, "make_delta: Returnvalue= ", Retvalue);
252  return Retvalue;
253}
254
255static CanonicalForm
256make_square( int levelU, const CanonicalForm & W,
257             const CanonicalForm & F1, const CanonicalForm & F2,
258             RememberArray & A, RememberArray & B,const CanonicalForm &alpha)
259{
260  CanonicalForm Retvalue;
261  DiophantForm intermediate;
262
263  DEBOUT(CERR, "make_square: W= ", W );
264  DEBOUTLN(CERR, "  degree(W,levelU)= ", degree(W,levelU));
265
266  if ( levelU == level(W) ) // same level, good
267  {
268    for ( CFIterator i=W; i.hasTerms(); i++)
269    {
270      intermediate=diophant(levelU,F1,F2,i.exp(),A,B,alpha);
271      Retvalue += i.coeff() * intermediate.Two;
272
273      if (interrupt_handle()) return Retvalue;
274    }
275  }
276  else // level(W) < levelU ; i.e. degree(w,levelU) == 0
277  {
278    intermediate=diophant(levelU,F1,F2,0,A,B,alpha);
279    Retvalue = W * intermediate.Two;
280  }
281  DEBOUTLN(CERR, "make_square: Returnvalue= ", Retvalue);
282
283  return Retvalue;
284}
285
286
287///////////////////////////////////////////////////////////////
288// Multivariat Hensel routine for two factors F and G .      //
289// U is the monic univariat polynomial; we manage two arrays //
290// to remember already calculated values for the diophantine //
291// equation. This is suggested by Joel Moses [Moses71] .     //
292// Return the fully lifted factors.                          //
293///////////////////////////////////////////////////////////////
294static DiophantForm
295mvhensel( const CanonicalForm & U , const CanonicalForm & F ,
296          const CanonicalForm & G , const SFormList & Substitutionlist,
297          const CanonicalForm &alpha)
298{
299  CanonicalForm V,Fk=F,Gk=G,Rk,W,D,S;
300  int  levelU=level(U), degU=subvardegree(U,levelU); // degree(U,{x_1,..,x_(level(U)-1)})
301  DiophantForm Retvalue;
302  RememberArray A(degree(F,levelU)+degree(G,levelU)+1);
303  RememberArray B(degree(F,levelU)+degree(G,levelU)+1);
304
305  DEBOUTLN(CERR, "mvhensel called with: U= ", U);
306  DEBOUTLN(CERR, "                      F= ", F);
307  DEBOUTLN(CERR, "                      G= ", G);
308  DEBOUTLN(CERR, "                   degU= ", degU);
309
310  V=change_poly(U,Substitutionlist,0); // change x_i <- x_i + a_i for all i
311  Rk = F*G-V;
312#ifdef HENSELDEBUG2
313  CERR << "mvhensel: V = " << V << "\n"
314       << "          Fk= " << F << "\n"
315       << "          Gk= " << G << "\n"
316       << "          Rk= " << Rk << "\n";
317#endif
318  for ( int k=2; k<=degU+1; k++)
319  {
320    W = mod_power(Rk,k,levelU);
321#ifdef HENSELDEBUG2
322    CERR << "mvhensel: Iteration: " << k << "\n";
323    CERR << "mvhensel: W= " << W << "\n";
324#endif
325    D = make_delta(levelU,W,F,G,A,B,alpha);
326#ifdef HENSELDEBUG2
327    CERR << "mvhensel: D= " << D << "\n";
328#endif
329    S = make_square(levelU,W,F,G,A,B,alpha);
330#ifdef HENSELDEBUG2
331    CERR << "mvhensel: S= " << S << "\n";
332#endif
333    Rk += S*D - D*Fk - S*Gk;
334#ifdef HENSELDEBUG2
335    CERR << "mvhensel: Rk= " << Rk << "\n";
336#endif
337    Fk -= S;
338#ifdef HENSELDEBUG2
339    CERR << "mvhensel: Fk= " << Fk << "\n";
340#endif
341    Gk -= D;
342#ifdef HENSELDEBUG2
343    CERR << "mvhensel: Gk= " << Gk << "\n";
344#endif
345    if ( Rk.isZero() ) break;
346    if (interrupt_handle()) break;
347  }
348  Retvalue.One = change_poly(Fk,Substitutionlist,1);
349  Retvalue.Two = change_poly(Gk,Substitutionlist,1);
350
351  DEBOUTLN(CERR, "mvhensel: Retvalue: ", Retvalue.One);
352  DEBOUTLN(CERR, "                    ", Retvalue.Two);
353
354  return Retvalue ;
355}
356
357///////////////////////////////////////////////////////////////
358// Recursive Version of MultiHensel.                         //
359///////////////////////////////////////////////////////////////
360CFFList
361multihensel( const CanonicalForm & mF, const CFFList & Factorlist,
362             const SFormList & Substitutionlist,
363             const CanonicalForm &alpha)
364{
365  CFFList Returnlist,factorlist=Factorlist;
366  DiophantForm intermediat;
367  CanonicalForm Pl,Pr;
368  int n = factorlist.length();
369
370  DEBOUT(CERR, "multihensel: called with ", mF);
371  DEBOUTLN(CERR, "  ", factorlist);
372
373  if ( n == 1 )
374  {
375    Returnlist.append(CFFactor(mF,1));
376  }
377  else
378  {
379    if ( n == 2 )
380    {
381      intermediat= mvhensel(mF, factorlist.getFirst().factor(),
382                            Factorlist.getLast().factor(),
383                            Substitutionlist,alpha);
384      Returnlist.append(CFFactor(intermediat.One,1));
385      Returnlist.append(CFFactor(intermediat.Two,1));
386    }
387    else  // more then two factors
388    {
389#ifdef HENSELDEBUG2
390      CERR << "multihensel: more than two factors!" << "\n";
391#endif
392      Pl=factorlist.getFirst().factor(); factorlist.removeFirst();
393      Pr=Pl.genOne();
394      for ( ListIterator<CFFactor> i=factorlist; i.hasItem(); i++ )
395        Pr *=  i.getItem().factor() ;
396#ifdef HENSELDEBUG2
397      CERR << "multihensel: Pl,Pr, factorlist: " << Pl << "  " << Pr
398           << "  " << factorlist << "\n";
399#endif
400      intermediat= mvhensel(mF,Pl,Pr,Substitutionlist,alpha);
401      Returnlist.append(CFFactor(intermediat.One,1));
402      Returnlist=Union( multihensel(intermediat.Two,factorlist,Substitutionlist,alpha), 
403                        Returnlist);
404    }
405  }
406  return Returnlist;
407}
408
409///////////////////////////////////////////////////////////////
410// Generalized Hensel routine.                               //
411// mF is the monic multivariat polynomial, Factorlist is the //
412// list of factors, Substitutionlist represents the ideal    //
413// <x_1+a_1, .. , x_r+a_r>, where r=level(mF)-1 .            //
414// Returns the list of fully lifted factors.                 //
415///////////////////////////////////////////////////////////////
416CFFList
417MultiHensel( const CanonicalForm & mF, const CFFList & Factorlist,
418             const SFormList & Substitutionlist, const CanonicalForm &alpha)
419{
420  CFFList Ll;
421  CFFList Returnlist,Retlistinter,factorlist=Factorlist;
422  CFFListIterator i;
423  DiophantForm intermediat;
424  CanonicalForm Pl,Pr;
425  int n = factorlist.length(),h=n/2, k;
426
427  DEBOUT(CERR, "MultiHensel: called with ", mF);
428  DEBOUTLN(CERR, "  ", factorlist);
429  DEBOUT(CERR,"           : n,h = ", n);
430  DEBOUTLN(CERR,"  ", h);
431
432  if ( n == 1 )
433  {
434    Returnlist.append(CFFactor(mF,1));
435  }
436  else
437  {
438    if ( n == 2 )
439    {
440      intermediat= mvhensel(mF, factorlist.getFirst().factor(),
441                            Factorlist.getLast().factor(),
442                            Substitutionlist,alpha);
443      Returnlist.append(CFFactor(intermediat.One,1));
444      Returnlist.append(CFFactor(intermediat.Two,1));
445    }
446    else  // more then two factors
447    {
448      for ( k=1 ; k<=h; k++)
449      {
450        Ll.append(factorlist.getFirst());
451        factorlist.removeFirst();
452      }
453
454      DEBOUTLN(CERR, "MultiHensel: Ll= ", Ll);
455      DEBOUTLN(CERR, "     factorlist= ", factorlist);
456
457      Pl = 1; Pr = 1;
458      for ( i = Ll; i.hasItem(); i++)
459        Pl *= i.getItem().factor();
460      DEBOUTLN(CERR, "MultiHensel: Pl= ", Pl);
461      for ( i = factorlist; i.hasItem(); i++)
462        Pr *= i.getItem().factor();
463      DEBOUTLN(CERR, "MultiHensel: Pr= ", Pr);
464      intermediat = mvhensel(mF,Pl,Pr,Substitutionlist,alpha);
465      // divison test for intermediat.One and intermediat.Two ?
466      CanonicalForm a,b;
467      // we add a division test now for intermediat.One and intermediat.Two
468      if ( mydivremt (mF,intermediat.One,a,b) && b == mF.genZero() )
469        Retlistinter.append(CFFactor(intermediat.One,1) );
470      if ( mydivremt (mF,intermediat.Two,a,b) && b == mF.genZero() )
471        Retlistinter.append(CFFactor(intermediat.Two,1)  );
472
473      Ll = MultiHensel(intermediat.One, Ll, Substitutionlist,alpha);
474      Returnlist = MultiHensel(intermediat.Two, factorlist, Substitutionlist,alpha);
475      Returnlist = Union(Returnlist,Ll);
476
477      Returnlist = Union(Retlistinter,Returnlist);
478
479    }
480  }
481  return Returnlist;
482}
Note: See TracBrowser for help on using the repository browser.