source: git/libfac/factor/MVMultiHensel.cc @ 14b1e65

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