source: git/libfac/factor/MVMultiHensel.cc @ 7d889c

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