source: git/factory/sm_sparsemod.cc @ 341696

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