source: git/libfac/factor/MVMultiHensel.cc @ 456842

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