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