source: git/factory/libfac/factor/MVMultiHensel.cc @ d92d71

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