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

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