source: git/factory/sm_sparsemod.cc @ e4fe2b

fieker-DuValspielwiese
Last change on this file since e4fe2b 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: 15.9 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4//{{{ docu
5//
6// sm_sparsemod.cc - zippel's sparse modular gcd.
7//
8// Contributed by Marion Bruder <bruder@math.uni-sb.de>.
9//
10// Used by: cf_gcf.cc
11//
12//}}}
13
14#include "config.h"
15
16#include "cf_assert.h"
17#include "debug.h"
18
19#include "canonicalform.h"
20#include "cf_map.h"
21#include "cf_reval.h"
22#include "cf_irred.h"
23#include "cf_iter.h"
24#include "cf_primes.h"
25#include "cf_random.h"
26#include "fac_util.h"
27#include "cf_algorithm.h"
28#include "sm_util.h"
29#include "sm_sparsemod.h"
30#include "templates/ftmpl_array.h"
31#include "templates/ftmpl_functions.h"
32
33static CanonicalForm
34smodgcd( const CanonicalForm & u, const CanonicalForm & v, const CanonicalForm & lcggt, const REvaluation & alpha, CFRandom & gen, int CHAR, const Variable & extension )
35{
36  DEBOUTLN( cerr, " ALPHA = " << alpha );
37
38  DEBOUTLN( cerr, " u = " << u );
39  DEBOUTLN( cerr, " v = " << v );
40
41  //    diverse Fallunterscheidungen  //
42
43  if ( u.isZero() )
44    return v;
45  else if ( v.isZero() )
46    return u;
47  else if ( u.inBaseDomain() )
48    if ( v.inBaseDomain() )
49      return gcd( u, v );
50    else
51      return gcd( u, icontent( u ) );
52  else if ( v.inBaseDomain() )
53    return gcd( icontent( u ), v);
54  else if ( u.isUnivariate() )
55  {
56    if ( v.isUnivariate() )
57      return gcd( u, v );
58    else
59      return gcd( v, u );
60  }
61
62  //   u und v Polynome in levU - Variablen x1, ..., xlevU
63  //   mit den gleichen Variablen, welches in sparsemod gesichert wurde
64
65  int levU = level( u );
66  CFArray G( 1, levU );
67  Variable x( 1 );
68  CanonicalForm alphau = alpha( u, 2, levU );
69  CanonicalForm alphav = alpha( v, 2, levU );
70
71  G[1] = gcd( alphau, alphav );       // univariater gcd berechnen in x1
72
73  DEBOUTLN( cerr, " G[1] ohne lc = " << G[1] );
74  //  Leitkoeffizienten auf univariaten gcd aufsetzen   //
75
76  if( alpha(lcggt) != 0 && degree(G[1]) != 0 )
77  {
78    //G[1] =( alpha( lcggt ))*( ( 1/lc( G[1] ) ) * G[1] );
79    G[1] =( alpha( lcggt ))*( ( 1/Leitkoeffizient( G[1] ) ) * G[1] );
80  }
81
82  if ( degree( G[1]) < 0 )
83  {
84   return 1 ;
85  }
86
87
88  if ( degree( G[1] ) == 0 )
89  {
90    // zwei Auswertungen um sicherzugehen, dass gcd Eins ist!
91    return 1;
92  }
93
94  CFIterator J = G[ 1 ];
95  CFArray g( 1, degree( G[ 1 ]) + 1 );
96  int * emma = new int[ degree( G[1] ) + 2  ];
97  int s, m, i = 1;
98
99  //  g[i] enthaelt die momentan berechneten Monome bzgl x^emma[i] //
100
101  while ( J.hasTerms() )
102  {
103    g[ i ] = J.coeff() ;
104    emma[ i ] = J.exp();
105    J++;
106    i++;
107  }
108
109  m = i-1;           /// Anzahl der Terme bzgl x1
110  int * n = new int[ m+1  ];
111  CFArray U( 2, levU ), V( 2, levU );
112
113  int  N, d   ;
114  int k, j, zip ;
115
116  // fuer jede Variabel interpolieren, um g[s] zu erhalten  //
117
118  DEBOUTLN( cerr, " unigcd mit alpha mit lc = " << G[1] );
119
120  for( s = 2 ;s <= levU ; s++ )
121  {
122    U[ s ] = alpha( u, s+1, levU );  // s Variablen stehen lassen, fuer
123    V[ s ] = alpha( v, s+1, levU );  // den Rest alpha einsetzen
124    d = tmin( degree( U[ s ], s ), degree( V[ s ], s ));
125
126    DEBOUTLN( cerr, " U["<< s << "] = " << U[s] );
127    DEBOUTLN( cerr, " V["<< s << "] = " << V[s] );
128    DEBOUTLN( cerr, " d = " << d );
129
130    for ( i = 1; i <= m ; i++ )
131    {
132      // Anzahl der Monome berechnen pro gi liefert das Skeletonpolynom //
133      n[ i ] = countmonome( g[ i ] );
134      //cout << "n["<<i<<"] = "<< n[i] << endl;
135    }
136
137    for ( i = 1; i <= m    ; i++ )
138    {
139      if ( i ==1 )
140        N = n[i];
141      else
142      {
143        if ( n[i]> N )
144          N = n[i];
145      }
146    }
147
148    //      int tau[d+1][m+1][N+1];               /// keine Integers !!!
149    typedef CanonicalForm** CF_ptr_ptr;
150    typedef CanonicalForm* CF_ptr;
151
152    CanonicalForm ***tau = new CF_ptr_ptr[m+1];
153    for ( i = 1; i <= m; i++ )
154    {
155      tau[i] = new CF_ptr[d+1];
156      for ( j = 1; j <= d; j++ )
157        tau[i][j] = new CanonicalForm[N+1];
158    }
159
160    CFArray beta( 1, d );
161
162    for ( i = 1; i <= d ; i++ )
163    {
164      beta[ i ] =  gen.generate();   // verschiedene Elemente aus Zp
165      //  cout << "  beta["<<i << "] = " << beta[i];
166    }
167
168    Array<REvaluation> xi( 1, N );
169    REvaluation a( 2, s-1, gen  );
170
171    for ( i = 1 ; i <= N ; i++ )
172    {
173      a.nextpoint();
174      xi[ i ] = a;//REvaluation( 2, s-1, gen );
175      // cout << "  xi["<<i<<"] = "<< xi[i];
176    }
177
178    //CFArray T(1, d*N );
179
180    if ( d == 0 )
181    {
182      ASSERT( 0, "alpha bad point! try some other gcd algorithm" );
183      return gcd(u, v);
184    }
185
186    CFMatrix T( d, N ) ;  // diese Moeglichkeit laeuft!!!
187                    //help = Koeff( T( j, k ), emma[ i ] );
188                    //tau[i][j][k] = help; //.intval();
189
190    for ( j = 1 ; j <= d ; j++ ) // jedesmal andere REvaluation??
191    {
192      for ( k = 1 ; k <= N ; k++ )
193      {
194        CanonicalForm zwischenu, zwischenv, help, nfa ;
195
196        zwischenu = U[ s ]( beta[ j ], s );
197        zwischenv = V[ s ]( beta[ j ], s );
198
199        T( j, k) = gcd ( xi[k]( zwischenu, 2, s-1 ), xi[k]( zwischenv, 2, s-1 ));
200
201        nfa = lcggt( beta[j], s );
202        nfa =  alpha( nfa, s+1, levU );
203
204        //T(j, k ) = (xi[k]( nfa, 2, s-1 ))*((1/lc(T( j, k ))) * T( j, k ));
205        T(j, k ) = (xi[k]( nfa, 2, s-1 ))*((1/Leitkoeffizient(T( j, k ))) * T( j, k ));
206
207        //cout <<"T("<<j<<", "<< k <<") = " << T(j, k) << endl;
208
209        for ( i = 1 ; i <= m ; i++ )
210        {
211          // diese Moeglichkeit laeuft!!!
212          //help = Koeff( T( j, k ), emma[ i ] );
213          //tau[i][j][k] = help; //.intval();
214          if ( T( j, k).inBaseDomain() )
215          {
216            if ( emma[i] == 0 )
217              tau[i][j][k] = T( j, k );
218            else
219              tau[i][j][k] =  0 ;
220          }
221          else
222          {
223            tau[i][j][k] = T(j, k)[emma[i]];
224          }
225        }
226      }
227
228      CFMatrix h( m, d );
229
230      for ( i = 1; i <= m ; i++ )
231      {
232        for ( j = 1 ; j <= d ; j++ )
233        {
234          zip = n[i] +1;
235          CanonicalForm * zwischen = new CanonicalForm[ zip ];//n[i]+1 ];
236
237          for ( k = 1 ; k <= n[i] ; k++ )
238          {
239            zwischen[ k ] = tau[i][j][k];
240          }
241
242          //cout <<" werte fuer sinterpol : " << endl;
243          //cout <<" g["<<i<<"] = " << g[i] << " xi = " << xi << endl;
244          //cout << " zwischen = " << zwischen << " ni = " << n[i]<<endl;
245          h[ i ][ j ] = sinterpol( g[i], xi, zwischen, n[i] );
246          DEBOUTLN( cerr, " Ergebnis von sinterpol h["<<i<<"]["<<j<< "] = " << h[i][j] );
247          delete [] zwischen;
248        }
249      }
250      for ( i = 1 ; i <= m ; i++ )
251      {
252        CFArray zwitscher( 1, d );
253        for ( j = 1 ; j <= d ; j++ )
254        {
255          zwitscher[ j ] = h[ i ][ j ];
256        }
257
258        DEBOUTLN( cerr, "Werte fuer dinterpol : " );
259        DEBOUTLN( cerr, " d = " << d << " g["<<i<<"] = "<< g[i] );
260        DEBOUTLN( cerr, " zwitscher = " << zwitscher );
261        DEBOUTLN( cerr, " alpha["<<s<< "] = "<< alpha[s] << " beta = "<<beta << " n["<<i<<"] = " << n[i] );
262
263        g[ i ] = dinterpol( d, g[i], zwitscher, alpha, s, beta, n[ i ], CHAR );
264        DEBOUTLN( cerr, " Ergebnis von dinterpol g["<<i<<"] = " << g[i] );
265
266      }
267
268      for ( i = 1 ; i <= m ; i++ )
269      {
270        G[ s ] += g[ i ] * ( power( x, emma[i] ) );
271      }
272      DEBOUTLN( cerr, "G["<<s<<"] = " << G[s] );
273      for ( i = 1; i <= m; i++ )
274      {
275        for ( j = 1; j <= d; j++ )
276          delete [] tau[i][j];
277        delete [] tau[i];
278      }
279      delete [] tau;
280    }
281  }
282  delete [] emma;
283
284  return G[ levU ];
285}
286
287// In sparsemod testen, ob G[levU] | u und G[levU] | v ///
288static CanonicalForm
289internalSparsemod( const CanonicalForm & F, const CanonicalForm & G )
290{
291  // On(SW_USE_EZGCD );
292  On( SW_SYMMETRIC_FF );
293  //cout << " in sparse " << endl;
294  if( F.inCoeffDomain() &&  G.inCoeffDomain() )
295  {
296    return gcd( F, G );
297  }
298
299  CanonicalForm f, g, ff, gg, ggt, res, fmodp, gmodp ;
300  int i, count = 10;
301
302  //   COMPRESS    ///////////////////
303
304  CFArray A(1, 2);
305  A[1] = F;
306  A[2] = G;
307  CFMap M, N;
308  compress(A, M, N);
309  f = M(A[1]);
310  g = M(A[2]);
311
312  // POLYNOME PRIMITIV BZGL DER ERSTEN VARIABLE  /////
313
314  CanonicalForm primif, primig, lcf, lcg, lcggt, lcggtmodp, result ;
315  ff = content( f, Variable(1) );//contentsparse( f, 1  );
316  gg = content( g, Variable(1) );//contentsparse( g, 1  );
317
318  primif = f/ff;
319  primig = g/gg;
320  ggt = gcd( ff, gg );
321
322  if( primif.inCoeffDomain() &&  primig.inCoeffDomain() )
323  {
324    return N( gcd( primif, primig ) ) * N( ggt );
325  }
326
327  // Variablen, die in beiden Polynomen auftreten /////
328
329  int levis, m, levelprimif, levelprimig;
330  int  ABFRAGE = 1 ;
331
332  levelprimif = level( primif );
333  levelprimig = level( primig );
334
335  if ( levelprimif > levelprimig )
336    levis = levelprimif;
337  else
338    levis = levelprimig;
339
340  if ( levis < 0 )
341    return N( gcd( primif, primig ) );
342
343  int * varf = new int[levis];
344  int * varg = new int[levis];
345  int * schnitt = new int[levis];
346
347  while ( ABFRAGE == 1 )
348  {
349    levelprimif = level(primif);
350    levelprimig = level(primig);
351
352    if ( levelprimif > levelprimig )
353      levis = levelprimif;
354    else
355      levis = levelprimig;
356
357    if ( levis < 0 )
358      return N( gcd(primif, primig ));
359
360    for( i = 1; i <= levis ; i++ )
361    {
362      if ( degree( primif, i ) != 0 )
363        varf[i] = 1;
364      else
365        varf[i] = 0;
366      if ( degree( primig, i ) != 0 )
367        varg[i] = 1;
368      else
369        varg[i] = 0;
370      if ( varg[i] == 1 && varf[i] == 1 )
371        schnitt[i] = 1;
372      else
373        schnitt[i] = 0;
374    }
375
376    levelprimif = level(primif);
377    levelprimig = level(primig);
378
379    for ( m  = 1; m <= levis ; m++)
380    {
381      if ( schnitt[m] == 0 )
382      {
383        if ( varf[m] == 1)
384        {
385          primif = content( primif, m ); //contentsparse( primif, m  );
386        }
387        else
388        {
389          primig = content( primig, m ); //contentsparse( primig, m  );
390        }
391      }
392    }
393
394    if ( level( primif ) == level( primig ) )
395      ABFRAGE = 0 ;
396
397  }
398
399  delete [] varf; delete [] varg; delete [] schnitt;
400
401  //  Nochmal compress fuer den Fall, dass eine Variable rausfliegt //
402
403  CFArray C(1, 2);
404  C[1] = primif;
405  C[2] = primig;
406  CFMap MM, NN ;
407  compress(C, MM, NN);
408  primif = MM(C[1]);
409  primig = MM(C[2]);
410
411  if ( level( primif) != level( primig ) )
412    ASSERT( 0, "this should n e v e r happen !!!!" );
413
414  //cout << " primif = " << primif << endl;
415  //cout << " primig = " << primig << endl;
416
417  if( primif.inCoeffDomain() &&  primig.inCoeffDomain() )
418  {
419    return N( NN( gcd( primif, primig ) ) ) * N( ggt );
420  }
421
422  // erst hier Leitkoeffizienten updaten  /////////
423
424  //if( primif.inCoeffDomain() || primif.isUnivariate
425
426  lcf = LC(primif, 1 );
427  lcg = LC(primig, 1 );
428  lcggt = gcd ( lcf, lcg );
429
430  //   BOUND BESTIMMEN fuer Charakteristik Null   /////////
431
432  int CHAR = getCharacteristic(), deg  ;
433  if ( CHAR < 10 )
434    deg = 3; // Default Erweiterung bei kleiner Primzahl;
435  else
436    deg = 1;
437  CanonicalForm B, L, Bound, lcF = Leitkoeffizient( primif); //lc( primif ) ;
438  B =  2 * maxNorm( primif ) * maxNorm( g ) ;
439  L = lcF ;
440  Bound = abs( 2 * B * L + 1 );
441  int length = 1;
442
443  CFArray p( 1, 10 ) ;
444
445  if( CHAR == 0 )
446  {
447    p[ 1 ]  = cf_getBigPrime( count );
448    CanonicalForm q = p[ 1 ];
449
450    while ( q < Bound )
451    {
452      count++;
453      length++;
454      p[ length ] = cf_getBigPrime( count );
455      q     *= p[ length ];  //G[levU] = ( 1/lc(G[levU] ))* G[levU];// sinnig?
456      //cout << " lcggt = " << lcggt << endl;
457      //bool ja;
458      //ja = fdivides( lcggt, lc( G[levU] ) );
459      //cout << " ja = " << ja << endl;
460      //cout << " vor Verlassen " << endl;
461    }
462  }
463  else
464  {
465    //int q  = CHAR;
466    //setCharacteristic( 0 );
467    //CanonicalForm Bound2 = mapinto( Bound ) ;
468    //cout << " Bound2 " << Bound2 << endl;
469    //cout << " CHAR =  " << q << endl;
470
471    // erstmal ohne Erweiterung !!!
472    //deg = 3; // Default Erweiterung
473    //q *= ( q * q ) ;cout << "q = " << q << endl;
474
475    //Bestimme Grad der Koerpererweiterung voellig unnuetz!!!
476    //while ( q < abs( Bound2 ) )
477    //        {
478    //  q *= CHAR;cout << " q = " << q << endl;
479    //deg++;cout << " degchar = " << deg << endl;
480    //cout << " in Graderhoehung? " << endl;
481
482    //}
483    //setCharacteristic( CHAR );
484    //cerr << " JUHU " << endl;
485  }
486
487  //        ENDE BOUNDBESTIMMUNG       /////////////////
488
489  FFRandom gen ;
490  levelprimif = level( primif );
491  REvaluation am( 2, levelprimif, FFRandom ( ));
492  int k, help ;
493  CFArray resultat( 1, length );//cout << " nach Resultat " << endl;
494  CanonicalForm zwischen;
495  Variable alpha1( levelprimif + 1 );
496  ABFRAGE = 0;
497
498
499  for ( i = 1 ; i <= 10; i++ )               // 10 Versuche mit sparsemod
500  {
501    if ( CHAR == 0 )
502    {
503      for( k = 1; k <= length ; k++)
504      {
505        help = p[ k ].intval();
506        setCharacteristic( help );
507        FFRandom mache;
508        am = REvaluation( 2, levelprimif, mache  );
509        am.nextpoint();
510        fmodp  = mapinto( primif );
511        gmodp = mapinto( primig );
512        lcggtmodp = mapinto ( lcggt );
513
514        zwischen = smodgcd( fmodp, gmodp, lcggtmodp, am, mache, CHAR, Variable( 1 ));
515        // Variable ( 1 ) Interface fuer endliche Koerper
516
517        // Char auf 0 wegen Chinese Remainder ////////
518
519        setCharacteristic( 0 );
520
521        resultat[ k ] = mapinto( zwischen );
522      }
523
524      ////////////////////////////////////////////////////////////
525      // result[k] mod p[k] via Chinese Remainder zu Resultat //////
526      // ueber Z hochziehen                                //////
527      ////////////////////////////////////////////////////////////
528
529      if( length != 1 )
530        ChinesePoly( length, resultat, p, result );
531      else
532        result = resultat[1];
533
534      CanonicalForm contentresult = content( result, 1 );
535
536      if ( contentresult != 0 )
537        res = result/contentresult;
538      else
539        res = result;
540
541      if ( fdivides( res, primif ) && fdivides( res, primig ) )
542        return  N(NN(res))*N(ggt) ;     /// compress rueckgaengig machen!
543      else
544      {
545
546        // Start all over again ///
547
548        count++;
549        for ( k = 1; k <= length ; k++)
550        {
551          p[k] = cf_getBigPrime( count );
552          count++;
553        }
554      }
555    }
556    else
557    {
558      // Fall Char != 0 ///
559      // in algebraische Erweiterung vom Grad deg gehen //
560      //cout << " IN CHAR != 0 " << endl;
561      //cout << " degree = " << deg << endl;
562      CanonicalForm minimalpoly;
563      //cout << " vor mini " << endl;
564      minimalpoly = find_irreducible( deg, gen, alpha1 );
565      //cout << " nach mini " << endl;
566      Variable alpha2 = rootOf( minimalpoly, 'a' ) ;
567      AlgExtRandomF hallo( alpha2 );
568      //cout << " vor am " << endl;
569      REvaluation am (  2, levelprimif, hallo );
570      //cout << " nach ma " << endl;
571      am.nextpoint();
572      //cout << "vor smodgcd " << endl;
573      result = smodgcd( primif, primig, lcggt, am, hallo, CHAR, alpha2  );
574      if ( result == 1 && ABFRAGE == 0)
575      {
576        // zwei Auswertungen fuer gcd Eins
577        am.nextpoint();
578        ABFRAGE = 1;
579      }
580      //CanonicalForm contentresult = contentsparse(  result, 1 );
581      //zuerst mal auf Nummer sicher gehen ...
582      else
583      {
584        CanonicalForm contentresult = content(  result, 1 );
585        //cout << "RESULT = " << result << endl;
586        if ( contentresult != 0 )
587          res = result/contentresult;
588        else
589          res = result;
590
591        if ( ( fdivides( res, primif )) && ( fdivides ( res, primig ) ))
592        {
593          // make monic ////
594          res = (1/lc(res)) * res;
595          // eventuell auch hier Leitkoeffizient anstatt lc ?
596
597          return  N( NN( res ) ) * N( ggt ) ;
598        }
599        else
600        {
601          // Grad der Erweiterung sukzessive um eins erhoehen
602          deg++;
603          //cout << " deg = " << deg << endl;
604          am.nextpoint();
605          // nextpoint() unnoetig?
606        }
607      }
608    }
609  }
610
611
612  //Fuer den Fall der unwahrscheinlichen Faelle, dass die Versuche
613  //nicht ausreichen :
614
615  ASSERT( 0, "sparsemod failed! try some other gcd algorithm" );
616  if( CHAR == 0 )
617    setCharacteristic( 0 );
618
619  return 0;
620}
621
622CanonicalForm
623sparsemod( const CanonicalForm & F, const CanonicalForm & G )
624{
625    Off( SW_USE_SPARSEMOD );
626    CanonicalForm result = internalSparsemod( F, G );
627    On( SW_USE_SPARSEMOD );
628    return result;
629}
Note: See TracBrowser for help on using the repository browser.