source: git/factory/bifac.cc @ 69fdf90

spielwiese
Last change on this file since 69fdf90 was e4fe2b, checked in by Oleksandr Motsak <motsak@…>, 13 years ago
FIX: Fixed huge BUG in cf_gmp.h CHG: starting to cleanup factory
  • Property mode set to 100644
File size: 32.6 KB
Line 
1#include "config.h"
2
3#include "canonicalform.h"
4
5#ifdef HAVE_BIFAC
6
7# include "lgs.h"
8#include "bifacConfig.h"
9
10#define BIFAC_BASIS_OF_G_CHECK        1
11
12void Reduce( bool );
13CanonicalForm Bigcd( const CanonicalForm& f, const CanonicalForm& g);
14
15CanonicalForm MyContent( const CanonicalForm& F, const Variable& x);
16CFFList Mysqrfree( const CanonicalForm& F );
17
18
19CanonicalForm MyGCDmod( const CanonicalForm & a,const CanonicalForm & b);
20CFFList RelFactorize(const CanonicalForm & h);
21
22//====== global definitions ===================
23Variable x( 'x' );
24Variable y( 'y' );
25Variable z( 'z' );
26Variable e( 'e' );
27
28
29///////////////////////////////////////////////////////
30// Class for storing polynomials as vectors.
31// Enables fast access to a certain degree.
32///////////////////////////////////////////////////////
33
34//==================================================
35class PolyVector
36//==================================================
37{
38public:
39  PolyVector ( void ){
40    m = -1;
41  }
42  virtual ~PolyVector( void ){
43    if( m!= -1) delete[] value;
44  }
45  void init (CanonicalForm f){
46    if( f.level()<0 )
47    {
48      m = 0;
49      n = 0;
50      value = new CanonicalForm[1];
51      value[0] = f;
52    }
53    else
54    {
55      m = degree(f,x);
56      n = degree(f,y);
57      ASSERT( m>0 || n>0, "Input is not a polynomial");
58      int correction = 1;  // univariate polynomials
59      if( n==0) correction = n+1;
60
61      value = new CanonicalForm[m*(n+1)+n+1];
62      for(int i=0; i<=m*(n+1)+n; i++) value[i]=0;
63
64
65      for ( CFIterator i = f; i.hasTerms(); i++ ) {
66        for ( CFIterator j = i.coeff(); j.hasTerms(); j++ ){
67                if( i.coeff().mvar().level()< 0 ){
68                        value[ 0*(n+1) + i.exp()*correction ] = j.coeff();}
69            else{
70                        value[ j.exp()*(n+1) + i.exp()*correction ] = j.coeff();}}}
71    }
72  }
73
74  void  push(int mm, int nn, CanonicalForm v){
75    ASSERT( 0<=mm<=m && 0<=nn<=n, "Wrong Index in PolyVector");
76    value[mm*(n+1)+nn] = v;
77  }
78  CanonicalForm get(int mm, int nn){
79
80    if( 0<=mm && mm<=m && 0<=nn && nn<=n )
81      return value[mm*(n+1)+nn];
82    else
83      return 0;
84  }
85#ifndef NOSTREAMIO
86  friend OSTREAM & operator<< ( OSTREAM & s, const PolyVector& V ){
87    for (int i=0;i<=V.m;i++)
88    {
89      s << "[";
90      for (int j=0;j<=V.n;j++)
91        s << V.value[i*(V.n+1)+j] << ", ";
92      s << "]\n";
93    }
94    return s;
95  }
96#endif /* NOSTREAMIO */
97
98
99private:
100  int            m;       // Degree in x
101  int            n;       // Degree in y
102  CanonicalForm* value;   // Value: index = m*(n+1)+n
103};
104////////// END of PolyVector ///////////////////////////
105
106
107
108
109/////////////////////////////////////////////////////////
110//
111//  Default class declarations
112//
113/////////////////////////////////////////////////////////
114
115
116
117//--<>---------------------------------
118BIFAC::BIFAC( void )// KONSTRUKTOR
119//--<>---------------------------------
120{
121}
122
123//--<>---------------------------------
124BIFAC::~BIFAC( void )// DESTRUKTOR
125//--<>---------------------------------
126{
127}
128
129
130/////////////////////////////////////////////////////////
131//
132//  Auxiliary functions
133//
134/////////////////////////////////////////////////////////
135
136//  //--<>---------------------------------
137//  void BIFAC::matrix_drucken( CFMatrix M )
138//  //--<>---------------------------------
139//  {
140//    int i,j;
141//    char* name="matrix.ppm";
142
143//    // === Datei löschen ===
144
145//    ofstream* aus = new ofstream(name, ios::out);
146//    delete aus;
147
148
149//    // === Jetzt immer nur anhängen ===
150
151//    aus  = new ofstream(name, ios::app);
152//    *aus << "// Zeilen Spalten\n"
153//         << "// x-Koord. y-Koord.  Wert\n";
154
155//    *aus  << M.rows()  << " " << M.columns() << endl;
156
157
158//    // === Noch nicht bearbeitet Teile ===
159
160//    for( i=0; i<M.rows(); i++)
161//      for( j=0; j<M.columns(); j++)
162//        *aus << i << " " << j << " " << M(i+1,j+1) << endl;;
163//    delete aus;
164//  }
165
166//=======================================================
167void  BIFAC::passedTime()
168//=======================================================
169{
170        ;
171}
172
173
174//=======================================================
175long int  BIFAC::anz_terme(  CanonicalForm & f )
176//=======================================================
177{
178  long int z=0;
179
180  for ( CFIterator i = f; i.hasTerms(); i++ )
181    for ( CFIterator j = i.coeff(); j.hasTerms(); j++ )
182      z++;
183  return( z );
184}
185
186//=======================================================
187void BIFAC::biGanzMachen(  CanonicalForm & f )
188//=======================================================
189{
190  CanonicalForm ggT;
191  bool init = false;
192  Off( SW_RATIONAL );
193
194  for ( CFIterator i = f; i.hasTerms(); i++ )
195    for ( CFIterator j = i.coeff(); j.hasTerms(); j++ )
196    {
197      if( !init )
198      {
199        ggT = j.coeff();
200        init = true;
201      }
202      else
203        ggT = gcd(j.coeff(), ggT);
204    }
205  f /= ggT;
206  On( SW_RATIONAL );
207}
208
209//=======================================================
210void  BIFAC::biNormieren( CanonicalForm & f )
211//=======================================================
212{
213  if ( getCharacteristic() == 0 )
214  {
215    for ( CFIterator i = f; i.hasTerms(); i++ )
216      for ( CFIterator j = i.coeff(); j.hasTerms(); j++ )
217        if( j.coeff().den() != 1 )
218        {
219          f  *= j.coeff().den();
220          biNormieren( f );
221        }
222    biGanzMachen( f );
223  }
224  else
225  {
226    f /= LC(f);
227  }
228}
229
230
231//=======================================================
232//   * Convert the basis vectors of G into polynomials
233//   * Validate the solutions
234//=======================================================
235CFList BIFAC::matrix2basis(CFMatrix A, int dim, int m, int n, CanonicalForm f)
236//=======================================================
237{
238  Variable x('x'), y('y');
239  int i,j,k;
240  CanonicalForm g,h,ff;
241  CFList Lg, Lh;
242
243  // === Construction of the 'g's ====
244  for(k=1; k<=dim; k++)
245  {
246    g=0;
247    for(i=0; i<=m-1; i++)
248      for(j=0; j<=n; j++)
249        g += A(k, i*(n+1)+j+1)* power(x,i) * power(y,j);
250    Lg.append(g);
251  }
252
253  ///////////  START VALIDATION ////////////////////////////////////
254  if (BIFAC_BASIS_OF_G_CHECK)
255  {
256
257    // === Construction of the 'h's ====
258    for(k=1; k<=dim; k++)
259    {
260      h=0;
261      for(i=0; i<=m; i++)
262        for(j=0; j<n; j++)
263          h += A(k, i*n+j+1 +m*(n+1))* power(x,i) * power(y,j);
264      Lh.append(h);
265    }
266
267    // === Is the solution correct? ===
268    CFListIterator itg=Lg;
269    CFListIterator ith=Lh;
270    for( ; itg.hasItem(); itg++, ith++)
271    {
272      g = itg.getItem();
273      h = ith.getItem();
274      ff = f*(deriv(g,y)-deriv(h,x)) +h*deriv(f,x) -g*deriv(f,y);
275      if( !ff.isZero()) {
276      #ifndef NOSTREAMIO
277        AUSGABE_ERR("* Falsche Polynome!");
278        exit (1);
279      #else
280        printf("wrong polys\n");
281        break;
282      #endif
283      }
284    }
285  }
286  ///////////  END VALIDATION ////////////////////////////////////
287
288  return (Lg);
289}
290
291/////////////////////////////////////////////////////////
292//
293//  Main functions
294//
295/////////////////////////////////////////////////////////
296
297//=======================================================
298//   * Create the matrix belonging to G
299//   * Compute a basis of the kernel
300//=======================================================
301CFList BIFAC::basisOfG(CanonicalForm f)
302//=======================================================
303{
304
305
306  int m = degree(f,x);
307  int n = degree(f,y);
308  int r,s, ii,jj;
309
310
311  // ======= Creation of the system of linear equations for G =============
312  int rows    = 4*m*n;
313  int columns = m*(n+1) + (m+1)*n;
314
315  CFMatrix M(rows, columns); // Remember: The first index is (1,1) -- not (0,0)!
316
317  for ( CFIterator i = f; i.hasTerms(); i++ )  // All coeffizients of y
318  {
319    for ( CFIterator j = i.coeff(); j.hasTerms(); j++ )  // All coeffizients of x
320    {
321      r = j.exp();  // x^r
322      s = i.exp();  // y^s
323
324      // Now we regard g_{ii,jj)
325      for( ii=0; ii<m; ii++)
326        for( jj=0; jj<=n; jj++)
327        {
328          if(  s>= 1) M( (r+ii)*2*n +(jj+s-1)+1, ii*(n+1)+jj +1) += -j.coeff() * s;
329          if( jj>= 1) M( (r+ii)*2*n +(jj+s-1)+1, ii*(n+1)+jj +1) +=  j.coeff() * jj;
330        }
331
332      // Now we regard h_{ii,jj}
333      for( ii=0; ii<=m; ii++)
334        for( jj=0; jj<n; jj++)
335        {
336          if(  r>= 1) M( (r+ii-1)*2*n +(jj+s)+1, (ii*n)+jj +m*(n+1) +1) += j.coeff() * r;
337          if( ii>= 1) M( (r+ii-1)*2*n +(jj+s)+1, (ii*n) +jj +m*(n+1) +1) +=  -j.coeff() * ii;
338        }
339    }
340  }
341  // ========= Solving the  system of linear equations for G =============
342
343//  matrix_drucken(M); // **********************************
344
345  LGS L(rows,columns);
346
347  CFMatrix Z(1,columns);
348  for( ii=1; ii<=rows; ii++)
349  {
350    for( jj=1; jj<=columns; jj++)
351      Z(1,jj) = M(ii,jj);  // Copy the ii-th row
352    L.new_row(Z);
353  }
354
355  if( L.corank() == 1 )
356  {
357    CFList Lg;
358    Lg.append(f);
359    return(Lg);
360  }
361//  L.print();
362  CFMatrix basis = L.GetKernelBasis();
363
364  // ============= TEST AUF KORREKTHEIT /start) ====
365  CanonicalForm tmp;
366  for(int k=1; k<= L.corank(); k++)
367    for(int i=1; i<=rows; i++)
368    {
369      tmp =0;
370      for(int j=1; j<=columns; j++)
371        tmp += M(i,j) * basis(k,j);
372      if( tmp!= 0) {
373        exit(17);
374      }
375    }
376  // ============= TEST AUF KORREKTHEIT (ende) ====
377  return (  matrix2basis( basis, L.corank(), m,n,f ) );
378}
379
380//=======================================================
381//   Compute a   r x r - matrix A=(a_ij) for
382//     gg_i = SUM a_ij * g_j * f_x (mod f)
383//  Return a list consisting of
384//    r x (r+1) Matrix A
385//    the last columns contains only the indices of the
386//    first r linear independent lines
387// REMARK: this is used by BIFAC::createEg but NOT by createEgUni!!
388//=======================================================
389CFMatrix BIFAC::createA (CFList G, CanonicalForm f)
390//=======================================================
391{
392  // === Declarations ===
393  int m,n;
394  int i,j,e;
395  int r = G.length();  // number of factors
396
397  LGS       L(r,r,true);
398//  LGS       L(r,r);
399  CFMatrix  Z(1,r);
400  CFMatrix  A(r,r+2);  // the last two column contain the bi-degree
401
402  CanonicalForm fx   = deriv(f,x);
403  PolyVector*   gifx = new  PolyVector[r];
404
405  // === Convert polynomials into vectors ===
406  i=0;
407  CanonicalForm q;
408
409  for( CFListIterator it=G; it.hasItem(); it++, i++){
410
411    gifx[i].init( (it.getItem()*fx)%f );
412  }
413
414  // === Search linear independent lines ===
415
416  e=1; // row number of A
417  n=0; //
418  m=0; //
419  while (L.rank() != r )
420  {
421    for(j=0;j<r;j++)
422      Z(1,j+1) = gifx[j].get(m,n);
423    if( L.new_row(Z,0) )  // linear independent row?
424    {
425      ASSERT( e<=r, "Wrong index in matrix");
426      A(e,r+1) = m;      // Degree in x
427      A(e,r+2) = n;      // Degree in y
428      e++;
429    }
430    if (m>n) n++;
431    else     { m++; n=0; }
432  }
433  L.print();
434
435  L.inverse(A);
436
437  // === Clean up ==
438  delete[] gifx;
439
440  return A;
441}
442
443//=======================================================
444CanonicalForm BIFAC::create_g (CFList G)
445//=======================================================
446{
447  CanonicalForm g = 0;
448  int i = 0;
449
450  int r    = G.length();  // number of factors
451  float SS =  10*( r*(r-1) / ( 2*( (100- (float) EgSeparable)/100)) );
452  int S    = (int) SS +1;
453
454  IntRandom RANDOM(S);
455
456  int*        rand_coeff1 = new int[r];
457
458
459  //  static for debugging
460  // rand_coeff1[0] = 12; rand_coeff1[1] = 91; rand_coeff1[2] = 42;
461
462  for( CFListIterator it=G; it.hasItem(); it++, i++)
463    {
464      rand_coeff1[i] =  RANDOM.generate().intval();
465
466      g += rand_coeff1[i]  * it.getItem();
467    }
468
469  delete[] rand_coeff1;
470
471  return g;
472}
473
474/////////////////////////////////////////////////////////////
475// This subroutine creates the polynomials Eg(x) and g
476// by using the 'bivariate' methode'.
477// REMARK: There is a 'univariate methode' as well
478//         which ought to be faster!
479////////////////////////////////////////////////////////////
480//=======================================================
481CFList BIFAC::createEg (CFList G, CanonicalForm f)
482//=======================================================
483{
484
485
486  CFMatrix NEU = createA(G,f);
487//  passedTime();
488
489  bool suitable1 = false; // Is Eg by chance unsuitable?
490  bool suitable2 = false;  // Is on of g*g_i or g_i*f_x zero?
491
492  // === (0) Preparation ===
493  CanonicalForm g;
494  CanonicalForm Eg;
495  CanonicalForm fx = deriv(f,x);
496
497  int i,j,e;
498  int r = G.length();  // number of factors
499//    float SS =  ( r*(r-1) / ( 2*( (100- (float) EgSeparable)/100)) );
500//    int S = (int) SS +1;
501
502//    IntRandom RANDOM(S);
503//    int*        rand_coeff = new int[r];
504  CFMatrix  A(r,r);
505  CanonicalForm*     gi  = new CanonicalForm[r];
506  CanonicalForm*    ggi  = new CanonicalForm[r];
507  PolyVector*     v_ggi  = new PolyVector   [r];
508
509
510
511  i=0;
512  for( CFListIterator it=G; it.hasItem(); it++, i++)
513    gi[i] =  it.getItem();
514
515  while ( !suitable1 )
516  {
517
518     suitable2 = false;
519    // === (1) Creating g ===
520    while ( !suitable2 )
521    {
522//        i=0;
523//        g=0;
524//        for( CFListIterator it=G; it.hasItem(); it++, i++)
525//        {
526//          gi[i] =  it.getItem();
527//          rand_coeff[i] =  RANDOM.generate().intval();
528//          g += rand_coeff[i] * it.getItem();
529//        }
530      g = create_g( G );
531
532      // === (2) Computing g_i * g ===
533          //
534      for(i=0; i<r; i++){
535
536          ggi[i]  = (g*gi[i])%f;   // seite 10
537      }
538
539      // ===  Check if all polynomials are <> 0  ===
540      suitable2 = true;    // It should be fine, but ...
541      if( g.isZero() )
542        suitable2 = false;
543//        else
544//          for(i=0; i<r; i++)
545//            if(  ggi[i].isZero() )
546//              suitable2 = false;
547
548    } // end of  Žwhile ( !suitable2 )Ž
549
550    // === (3) Computing Eg(x) ===
551
552    for(i=0;i<r;i++)  // Get Polynomials as vectors
553      v_ggi[i].init(ggi[i]);
554
555    // Matrix A
556    for(i=1; i<=r; i++)
557      for( j=1; j<=r; j++)
558      {
559        A(i,j) = 0;
560        for( e=1; e<=r; e++)
561        {
562
563
564   A(i,j) += ( NEU(j,e ) * v_ggi[i-1].get(NEU(e,r+1).intval(),(NEU(e,r+2).intval() )));
565
566
567//
568
569        }
570      }
571
572    for(j=1; j<=r; j++)
573      A(j,j) -= x;
574    Eg = determinant(A,r);
575//    exit(1);
576    // === (4) Is Eg(x) suitable? ===
577    if( MyGCDmod(Eg, deriv(Eg,x)) == 1 )
578      suitable1 = true;
579    else
580    {
581    }
582  } // end of  Žwhile ( !suitable1 )Ž
583
584  // Delete trash
585
586
587
588  delete[] v_ggi;
589  delete[] gi;
590  delete[] ggi;
591  //  delete[] rand_coeff;
592
593  CFList LL;
594  LL.append(Eg);
595  LL.append(g);
596  return (LL);
597}
598//  /////////////////////////////////////////////////////////////
599//  // It is possible to take univariate polynomials
600//  // with y:=c for a suitable c.
601//  // c is suitable iff gcd( f(x,c), f_x(x,c)) = 1.
602//  ////////////////////////////////////////////////////////////
603//
604//=======================================================
605CFList BIFAC::createEgUni (CFList G, CanonicalForm f)
606//=======================================================
607{
608
609  int i,ii,k;
610  CanonicalForm ff, ffx,g, gg, Eg;
611
612
613  bool suitable1 = false;  // Is Eg unsuitable?
614  bool suitable2 = false;  // Is on of g*g_i or g_i*f_x zero?
615  bool suitable3 = false;  // Is 'konst' unsuitable?
616
617  // ========================
618  // =   (0) Preparation    =
619  // ========================
620  int konst = 0;
621  CanonicalForm fx = deriv(f,x);
622  int m = degree(f,x);
623  int r = G.length();            // number of factors
624  int S =  (int) ((float) ( r*(r-1) / ( 2*( (100- (float) EgSeparable)/100)) )+1);
625
626
627  int*  rand_coeff      = new int[r];
628  CanonicalForm*     gi  = new CanonicalForm[r];
629  CanonicalForm*    ggi  = new CanonicalForm[r];
630
631  CFMatrix  A(r,r);     // We have to find the matrix A,
632  CFMatrix  Z(1,r);     // `Vector` for data transportation
633  CFMatrix  AA(m,r);    // but first we generate AA.
634  CFMatrix  AI(r,r+1);  //
635  LGS       L(r,r,true);
636  IntRandom RANDOM(S);
637
638
639  // ==========================================================
640  // =  (1) Find a suitable constant to make bivariate        =
641  // =      polynomials univariate. Try the following numbers =
642  // =      0, 1, -1, 2, -2, 3,...                            =
643  // ==========================================================
644
645  while ( !suitable3 )
646  {
647    ff  = f(konst,'y');
648    ffx = fx(konst,'y');
649
650    if( gcd(ff, ffx) == 1)
651      suitable3 = true;
652    else
653    {
654      konst *= -1;
655      if( konst >= 0 )
656        konst++;
657    }
658  }
659
660
661  // ===============================================
662  // =    (2) Make g_i univariate                  =
663  // ===============================================
664  i=0;
665  for( CFListIterator it=G; it.hasItem(); it++, i++)
666  {
667    gi[i] =  it.getItem()(konst,'y');
668  }
669
670  // ===============================================
671  // =   (3) Compute the matrices 'AA' and 'AI'    =
672  // ===============================================
673
674
675  for( i=0; i<r; i++) // First store all coeffizients in AA.
676  {
677    ggi[i] = (gi[i]*ffx)%ff;   // now we have degree < m.
678    //biNormieren(ggi[i]);
679    for ( CFIterator j = ggi[i]; j.hasTerms(); j++ )
680      AA( j.exp()+1, i+1) = j.coeff();
681  }
682
683
684  // Now find the lin. indep. rows.
685  i  = 1;
686  ii = 1; // row number of A
687  while (L.rank() != r )
688  {
689    ASSERT( i<=m, "Too few linear independent rows!");
690
691    for (k=1; k<=r; k++)
692      Z(1,k) =  AA(i,k);
693    if( L.new_row(Z,0) )  // linear independent row?
694    {
695      ASSERT( ii<=r, "Wrong index in matrix");
696      AI(ii,r+1) = i;      // Degree in x
697      ii++;
698    }
699    i++;
700    L.print();
701  }
702  L.inverse(AI);
703
704
705  // ==============================================
706  // =   (4) Big loop to find a suitable 'Eg(x)   =
707  // ==============================================
708
709  while ( !suitable1 )    // Is Eg(x) suitable? -> Check at the end of this procedure!
710  {
711    suitable2 = false;   // In case we need a second loop
712
713    // ================================================
714    // =    (4a) Find a suitable 'g'                  =
715    // ================================================
716//    rand_coeff[0] = 0;
717//    rand_coeff[1] = 4;
718
719
720    while ( !suitable2 )
721    {
722      // === (i) Creating g ===
723      i=0;
724      g=0;
725      for( CFListIterator it=G; it.hasItem(); it++, i++)
726      {
727         rand_coeff[i] =  RANDOM.generate().intval();
728        g += rand_coeff[i] * it.getItem();
729      }
730      gg = g(konst,'y');   // univariate!
731      for(i=0; i<r; i++)  ggi[i] = (gi[i]*gg)%ff; // !! Redefinition of ggi !!
732
733      // ===  (ii) Check if all polynomials are <> 0  ===
734      suitable2 = true;    // It should be fine, but ...
735      if( gg.isZero() )
736        suitable2 = false;
737//        else
738//          for(i=0; i<r; i++)
739//            if(  ggi[i].isZero() )
740//              suitable2 = false;
741    } // end of  Žwhile ( !suitable2 )Ž
742
743//    createRg(g,f);
744
745    // ===============================================
746    // =    (b) Compute matrix 'A'                   =
747    // ===============================================
748    for(i=1; i<=r; i++)
749    {
750      for( ii=1; ii<=m; ii++)
751        AA (ii,1) = 0;  // !! Redefinition of AA !!
752      for ( CFIterator j = ggi[i-1]; j.hasTerms(); j++ )
753        AA( j.exp()+1, 1) = j.coeff();
754
755      for( ii=1; ii<=r; ii++)
756      {
757        A(i,ii) = 0;
758        for( k=1; k<=r; k++)
759          A(i,ii) += ( AI(ii,k ) *  AA( AI(k, r+1 ).intval(),1) );
760      }
761    }
762    for(i=1; i<=r; i++)
763      A(i,i) -= x;
764
765    // ===============================================
766    // =    (c) Compute Eg(x) and check it           =
767    // ===============================================
768
769    Eg = determinant(A,r);
770    if( gcd(Eg, deriv(Eg,x)) == 1 )
771    {
772      suitable1 = true;
773    }
774  } // end of  Žwhile ( !suitable1 )Ž
775
776
777  // ==============================================
778  // =   (5) Prepare for leaving                  =
779  // ==============================================
780
781  delete[] gi;
782  delete[] ggi;
783  delete[] rand_coeff;
784
785  CFList LL;
786  LL.append(Eg);
787  LL.append(g);
788
789  return (LL);
790}
791/////////////////////////////////////////////////////////////
792// This subroutine creates the polynomials Rg(x)
793// which can be used instead of Eg(x).
794// No basis of G is neccessary, only one element
795////////////////////////////////////////////////////////////
796//=======================================================
797CFList BIFAC::createRg (CFList G, CanonicalForm f)
798//=======================================================
799{
800
801//  cerr << "* Was ist wenn g versagt???? -> Ausbauen\n";
802
803  CanonicalForm fx = deriv(f,x);
804  CanonicalForm Rg;
805  CanonicalForm g = create_g(G);
806
807
808  // ===============================================
809  // =   (1) Find a suitable constant              =
810  // ===============================================
811
812  CanonicalForm alpha=1;
813
814  while(  resultant( f, fx, x)(alpha) == 0 )
815  {
816        //while( resultant( f, fx, x)(alpha).inCoeffDomain() != true )
817    //alpha +=1;
818  }
819
820
821  // ===============================================
822  // =   (2) Find a suitable constant              =
823  // ===============================================
824
825  Rg = resultant( f(alpha,y), g(alpha,y)-z*fx(alpha,y), x);
826
827
828  CFList LL;
829  LL.append(Rg(x,z));
830  LL.append(g);
831  return (LL);
832}
833/////////////////////////////////////////////////////////
834// Compute the absolute and rational factorization of
835// the univariate polynomial 'ff^grad'.
836//=======================================================
837void BIFAC::unifac (CanonicalForm ff, int grad)
838//=======================================================
839{
840
841  CFFList factorsUni;
842  CFFList factorsAbs;
843  CanonicalForm tmp;
844
845  factorsUni = AbsFactorize(ff);
846
847  for( CFFListIterator l=factorsUni; l.hasItem(); l++)
848    if( ! l.getItem().factor().inBaseDomain() )
849    {
850      gl_RL.append( CFFactor( l.getItem().factor(),l.getItem().exp()*grad) );
851    }
852
853
854}
855
856
857///////////////////////////////////////////////////////
858// Compute the rational factor of f belonging to phi
859//=======================================================
860CanonicalForm BIFAC::RationalFactor (CanonicalForm phi, CanonicalForm ff, \
861                                     CanonicalForm fx, CanonicalForm g)
862//=======================================================
863{
864
865  CanonicalForm h,hh;
866//  CanonicalForm fx = deriv(f,x);
867
868  for ( CFIterator it = phi; it.hasTerms(); it++ )
869    h += it.coeff() * power(fx,phi.degree()-it.exp())*power(g,it.exp());
870
871
872  hh = Bigcd(ff,  h);
873
874  return(hh);
875}
876//=======================================================
877void BIFAC::RationalFactorizationOnly (CFFList Phis, CanonicalForm f0, CanonicalForm g)
878//=======================================================
879{
880  CanonicalForm h,ff;
881  CanonicalForm fx = deriv(f0,x);
882
883  for( CFFListIterator i=Phis; i.hasItem(); i++)
884  {
885    ASSERT( i.getItem().exp() == 1 , "Wrong factor of Eg"); // degree must be 1
886    CanonicalForm phi = i.getItem().factor();
887
888    if( ! phi.inBaseDomain())
889    {
890      h = RationalFactor(phi,f0,fx,g);
891      gl_RL.append( CFFactor(h,exponent ));
892      ff = f0;
893      f0 /= h;
894      ASSERT( f0*h==ff, "Wrong factor found");
895    }
896  }
897}
898
899//=======================================================
900CFList BIFAC::getAbsoluteFactors (CanonicalForm f1, CanonicalForm phi)
901//=======================================================
902{
903  CanonicalForm fac;
904  CanonicalForm root;
905  CFList AbsFac;
906
907  CFFList Fac = factorize(phi,e);
908  for( CFFListIterator i=Fac; i.hasItem(); i++)
909  {
910    fac = i.getItem().factor();
911    if( taildegree(fac) > 0 )  // case: phi = a * x
912      root = 0;
913    else
914      root = -tailcoeff(fac)/lc(fac);
915
916
917    AbsFac.append( f1(root,e) );
918    AbsFac.append( i.getItem().exp() * exponent);
919    AbsFac.append( phi ); // Polynomial of the field extension
920  }
921  return AbsFac;
922}
923//=======================================================
924void BIFAC::AbsoluteFactorization (CFFList Phis, CanonicalForm ff, CanonicalForm g)
925//=======================================================
926{
927
928  int ii;
929  if( getCharacteristic() == 0 )
930  {
931    //cerr << "* Charcteristic 0 is not yet implemented! => Aborting!\n";
932    exit(1);
933  }
934
935
936  CFList  AbsFac;
937  CanonicalForm phi;
938  CanonicalForm h, h_abs, h_res, h_rat;
939  CanonicalForm fx = deriv(ff,x);
940
941
942  for( CFFListIterator i=Phis; i.hasItem(); i++)
943  {
944    ASSERT( i.getItem().exp() == 1 , "Wrong factor of Eg"); // degree must be 1
945    phi = i.getItem().factor();
946
947    if( ! phi.inBaseDomain())
948    {
949
950      // === Case 1:  phi has degree 1 ===
951      if( phi.degree() == 1 )
952      {
953        if( taildegree(phi) > 0 )  // case: phi = a * x
954          h = gcd( ff,g );
955        else                       // case: phi = a * x + c
956        {
957          h =  gcd( ff, g+tailcoeff(phi)/lc(phi)*fx);
958        }
959
960        //biNormieren( h );
961        gl_AL.append(h); // Factor of degree 1
962         gl_AL.append(exponent); // Multiplicity (exponent)
963        gl_AL.append(0); // No field extension
964      } else
965      {
966        // === Case 2:  phi has degree > 1 ===
967        e=rootOf(phi, 'e');
968        h =  gcd( ff, g-e*fx);
969        //biNormieren( h );
970
971        AbsFac = getAbsoluteFactors(h, phi);
972        for( CFListIterator l=AbsFac; l.hasItem(); l++)
973          gl_AL.append( l.getItem() );
974
975
976        // ===  (1) Get the rational factor by multi-  ===
977        // ===      plication of the absolute factor.  ===
978        h_abs=1;
979        ii = 0;
980
981        for( CFListIterator l=AbsFac; l.hasItem(); l++)
982        {
983          ii++;
984          if (ii%3 == 1 )
985            h_abs *= l.getItem();
986        }
987        //biNormieren( h_abs );
988
989
990        // === (2) Compute the rational factor  ===
991        // ===     by using the resultant.      ===
992        h_res =  resultant(phi(z,x), h(z,e), z);
993        //biNormieren( h_res );
994
995
996        // === (3) Compute the rational factor by ignoring  ===
997        // ===     all knowledge of absolute factors.       ===
998        h_rat = RationalFactor(phi, ff,fx, g);
999        //biNormieren( h_rat );
1000
1001        ASSERT(  (h_abs == h_res) && (h_res == h_rat), "Wrong rational factor ?!?");
1002        h = h_abs;
1003      }
1004      // End of absolute factorization.
1005      gl_RL.append(CFFactor( h,exponent )); // Save the rational factor
1006      ff/=h;
1007    }
1008  }
1009}
1010
1011
1012//======================================================
1013//  Factorization of a squarefree bivariate polynomial
1014//  in which every factor appears only once.
1015//  Do we need a complete factorization ('absolute' is true)
1016//  or only a rational factorization ('absolute' false)?
1017//======================================================
1018void BIFAC::bifacSqrFree(CanonicalForm ff)
1019//=======================================================
1020{
1021
1022  int anz=0;                  // number of factors without field elements
1023
1024  CFList G   = basisOfG(ff);
1025
1026  CFList LL;
1027  CanonicalForm Eg,g;
1028
1029
1030
1031  // Case 1: There is only one rational & absolute factor ===
1032  if( G.length() == 1 ){                        // There is only one
1033    gl_RL.append( CFFactor(ff, exponent));      // rational factor
1034    gl_AL.append( ff );
1035    gl_AL.append( exponent );
1036    gl_AL.append( 0 );
1037  }
1038  else // Case 2: There is more than  one absolute factor ===
1039  {
1040//    LL  = createEg(G,ff);
1041//   LL = createEgUni(G,ff); // Hier ist noch ein FEHLER !!!!
1042
1043   LL = createRg( G, ff);  // viel langsamer als EgUni
1044
1045
1046    Eg  =  LL.getFirst();
1047        Eg  =  Eg/LC(Eg);
1048
1049   g   =  LL.getLast();
1050
1051//      g = G.getFirst();
1052
1053
1054    CFFList PHI = AbsFactorize( Eg );
1055
1056        CFFListIterator J=PHI;
1057        CanonicalForm Eg2=1;
1058         for ( ; J.hasItem(); J++)
1059        { Eg2 = Eg2 * J.getItem().factor(); }
1060
1061    // === Is Eg(x) irreducible ? ===
1062    anz=0;
1063
1064        // PHI =  AbsFactorize( Eg) ;
1065        //
1066
1067    for( CFFListIterator i=PHI; i.hasItem(); i++) {
1068      if( !i.getItem().factor().inBaseDomain())
1069        anz++;
1070         }
1071
1072    /* if( absolute ) // Only for a absolute factorization
1073      AbsoluteFactorization( PHI,ff, g);
1074    else         // only for a rational factorization
1075    { */
1076      if( anz==1 ){ ;
1077        gl_RL.append( CFFactor(ff,exponent));}
1078      else
1079        RationalFactorizationOnly( PHI,ff, g);
1080   /* } */
1081  }
1082}
1083
1084/////////////////////////////////////////////
1085// Main procedure for the factorization
1086// of the bivariate polynomial 'f'.
1087// REMARK: 'f' might be univariate, too.
1088//--<>---------------------------------
1089void BIFAC::bifacMain(CanonicalForm  f)
1090//--<>---------------------------------
1091{
1092
1093
1094  CanonicalForm ff, ggT;
1095
1096  // ===============================================
1097  // =    (1) Trivial case: Input is a constant    =
1098  // ===============================================
1099
1100  if( f.inBaseDomain() )
1101  {
1102    gl_AL.append(f); // store polynomial
1103    gl_AL.append(1); // store exponent
1104    gl_AL.append(0); // store ŽpolynomialŽ for field extension
1105
1106    gl_RL.append( CFFactor(f,1) ); // store polynomial
1107    return;
1108  }
1109
1110  // ===============================================
1111  // =       STEP: Squarefree decomposition        =
1112  // ===============================================
1113
1114
1115        CFFList Q =Mysqrfree(f);
1116//
1117//        cout << Q << endl;
1118//
1119
1120
1121
1122  // =========================================================
1123  // =  STEP: Factorization of the squarefree decomposition  =
1124  // =========================================================
1125
1126
1127  for( CFFListIterator i=Q; i.hasItem(); i++)
1128  {
1129
1130        if( i.getItem().factor().level() < 0 ) ;
1131        else
1132        {
1133    if( ( degree(i.getItem().factor(),x) == 0 || degree( i.getItem().factor(),y) == 0) ) {
1134      // case: univariate
1135      unifac( i.getItem().factor(), i.getItem().exp()  ); }
1136    else // case: bivariate
1137    {
1138      exponent =  i.getItem().exp();       // global variable
1139          CanonicalForm dumm = i.getItem().factor();
1140          dumm = dumm.LC();
1141          if( dumm.level() > 0 ){ dumm =  1;  }
1142      bifacSqrFree(i.getItem().factor()/dumm );
1143    }
1144  }}
1145
1146
1147}
1148
1149
1150///////////////////////////////////////////////////////
1151// Find the least prime so that the factorization
1152// works.
1153///////////////////////////////////////////////////////
1154
1155//=======================================================
1156int BIFAC::findCharacteristic(CanonicalForm f)
1157//=======================================================
1158{
1159  int min = (2*degree(f,'x')-1)*degree(f,'y');
1160  int nr=0;
1161
1162  if( min >= 32003 ) return ( 32003 ); // this is the maximum
1163
1164  // Find the smallest poosible prime
1165  while ( cf_getPrime(nr) < min)  { nr++;  }
1166  return ( cf_getPrime(nr) );
1167}
1168
1169/////////////////////////////////////////////////////////
1170//
1171//  PUBLIC functions
1172//
1173/////////////////////////////////////////////////////////
1174
1175// convert the result of the factorization from
1176// the intern storage type into the public one.
1177// Also, check the correctness of the solution
1178// and, if neccessary, change the characteristic.
1179//--<>---------------------------------
1180void BIFAC::convertResult(CanonicalForm & f, int ch, int sw)
1181//--<>---------------------------------
1182{
1183
1184  CanonicalForm ff = 1;
1185  CanonicalForm c;
1186
1187  CFFList aL;
1188
1189  //cout << gl_RL<<endl;
1190
1191        if( sw )
1192        {
1193                Variable W('W');
1194                for( CFFListIterator i=gl_RL; i.hasItem(); i++)
1195            {
1196                        c = i.getItem().factor();
1197                        c = c(W,y);
1198                        c = c(y,x);
1199                        c = c(x,W);
1200                        aL.append( CFFactor( c, i.getItem().exp() ));
1201                }
1202
1203                f = f(W,y); f=f(y,x); f=f(x,W);
1204        }
1205        else aL = gl_RL;
1206
1207        gl_RL = aL;
1208
1209        //cout << aL;
1210
1211
1212
1213  // ==========  OUTPUT  =====================
1214/*  for( CFFListIterator i=aL; i.hasItem(); i++)
1215  {
1216    cout << "(" << i.getItem().factor() << ")";
1217    if( i.getItem().exp() != 1 )
1218      cout << "^" << i.getItem().exp();
1219    cout << " * ";
1220  } */
1221
1222
1223//  cout << "\n* Test auf Korrektheit ...";
1224
1225
1226  for( CFFListIterator i=aL; i.hasItem(); i++)
1227    {
1228      ff *= power(i.getItem().factor(),  i.getItem().exp() );
1229      //      cout << " ff = " << ff
1230      //           << "\n a^b = " << i.getItem().factor() << "  ^ " <<   i.getItem().exp() << endl;
1231    }
1232  c = f.LC()/ff.LC();
1233
1234  ff *= c;
1235
1236
1237//   cout << "\n\nOriginal f = " << f << "\n\nff = " << ff
1238//           << "\n\nDiff = " << f-ff << endl << "Quot "<< f/ff <<endl;
1239//  cout << "degree 0: " << c << endl;
1240
1241
1242#ifndef NOSTREAMIO
1243  if( f != ff ) cout << "\n\nOriginal f = " << f << "\n\nff = " << ff
1244                     << "\n\nDiff = " << f-ff << endl << "Quot "<< f/ff <<endl;
1245#endif
1246  ASSERT( f==ff, "Wrong rational factorization. Abborting!");
1247//  cout << "  [OK]\n";
1248
1249}
1250//--<>---------------------------------
1251void BIFAC::bifac(CanonicalForm f, bool abs)
1252//--<>---------------------------------
1253{
1254  absolute = 1;      // global variables
1255  CFList  factors;
1256  int ch =  getCharacteristic();
1257  int ch2;
1258
1259
1260  ASSERT( ch==0 && !isOn(SW_RATIONAL), "Integer numbers not allowed" );
1261
1262
1263  // === Check the characteristic ===
1264  if( ch != 0 )
1265  {
1266    ch2 = findCharacteristic(f);
1267    if( ch <  ch2 )
1268    {
1269//      setCharacteristic( ch2 );
1270      f = mapinto(f);
1271
1272      // PROVISORISCH
1273      //cerr << "\n Characteristic is too small!"
1274         //  << "\n The result might be wrong!\n\n";
1275      exit(1);
1276
1277    } else ;
1278  }
1279
1280         Variable W('W');
1281          CanonicalForm l;
1282        int sw = 0;
1283
1284        if( degree(f,x) < degree(f,y) ) {
1285                f = f(W,x);   f = f(x,y); f=f(y,W);
1286                sw = 1;
1287        }
1288                l = f.LC();
1289
1290                if( l.level()<0 ) { f = f/f.LC(); gl_RL.append( CFFactor(l,1) ); }
1291
1292
1293  bifacMain(f);                   // start the computation
1294
1295  convertResult(f,ch, sw) ; // and convert the result
1296}
1297
1298// ==============  end of 'bifac.cc'  ==================
1299#endif /* #ifdef HAVE_BIFAC */
Note: See TracBrowser for help on using the repository browser.