source: git/libfac/factor/MVMultiHensel.cc @ 91b36d

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