source: git/factory/cf_gcd.cc @ db1996e

fieker-DuValspielwiese
Last change on this file since db1996e was db1996e, checked in by Wilfred Pohl <pohl@…>, 18 years ago
bugfix git-svn-id: file:///usr/local/Singular/svn/trunk@8923 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 17.9 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: cf_gcd.cc,v 1.38 2006-02-01 09:37:07 pohl Exp $ */
3
4#include <config.h>
5
6#include "assert.h"
7#include "debug.h"
8
9#include "cf_defs.h"
10#include "canonicalform.h"
11#include "cf_iter.h"
12#include "cf_reval.h"
13#include "cf_primes.h"
14#include "cf_algorithm.h"
15#include "cf_map.h"
16#include "sm_sparsemod.h"
17#include "fac_util.h"
18#include "ftmpl_functions.h"
19
20#ifdef HAVE_NTL
21#include <NTL/ZZX.h>
22#include "NTLconvert.h"
23bool isPurePoly(const CanonicalForm & f);
24#endif
25
26static CanonicalForm gcd_poly( const CanonicalForm & f, const CanonicalForm& g, bool modularflag );
27static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g );
28
29bool
30gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap )
31{
32    int count = 0;
33    // assume polys have same level;
34    CFRandom * sample = CFRandomFactory::generate();
35    REvaluation e( 2, tmax( f.level(), g.level() ), *sample );
36    delete sample;
37    CanonicalForm lcf, lcg;
38    if ( swap ) {
39        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
40        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
41    }
42    else {
43        lcf = LC( f, Variable(1) );
44        lcg = LC( g, Variable(1) );
45    }
46    while ( ( e( lcf ).isZero() || e( lcg ).isZero() ) && count < 100 ) {
47        e.nextpoint();
48        count++;
49    }
50    if ( count == 100 )
51        return false;
52    CanonicalForm F, G;
53    if ( swap ) {
54        F=swapvar( f, Variable(1), f.mvar() );
55        G=swapvar( g, Variable(1), g.mvar() );
56    }
57    else {
58        F = f;
59        G = g;
60    }
61    if ( e(F).taildegree() > 0 && e(G).taildegree() > 0 )
62        return false;
63    return gcd( e( F ), e( G ) ).degree() < 1;
64}
65
66//{{{ static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q )
67//{{{ docu
68//
69// balance() - map f from positive to symmetric representation
70//   mod q.
71//
72// This makes sense for univariate polynomials over Z only.
73// q should be an integer.
74//
75// Used by gcd_poly_univar0().
76//
77//}}}
78static CanonicalForm
79balance ( const CanonicalForm & f, const CanonicalForm & q )
80{
81    Variable x = f.mvar();
82    CanonicalForm result = 0, qh = q / 2;
83    CanonicalForm c;
84    CFIterator i;
85    for ( i = f; i.hasTerms(); i++ ) {
86        c = mod( i.coeff(), q );
87        if ( c > qh )
88            result += power( x, i.exp() ) * (c - q);
89        else
90            result += power( x, i.exp() ) * c;
91    }
92    return result;
93}
94//}}}
95
96//{{{ static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
97//{{{ docu
98//
99// icontent() - return gcd of c and all coefficients of f which
100//   are in a coefficient domain.
101//
102// Used by icontent().
103//
104//}}}
105static CanonicalForm
106icontent ( const CanonicalForm & f, const CanonicalForm & c )
107{
108    if ( f.inCoeffDomain() )
109        return gcd( f, c );
110    else {
111        CanonicalForm g = c;
112        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
113            g = icontent( i.coeff(), g );
114        return g;
115    }
116}
117//}}}
118
119//{{{ static CanonicalForm bcontent ( const CanonicalForm & f, const CanonicalForm & b )
120//{{{ docu
121//
122// bcontent() - return gcd of b and all coefficients of f which
123//   are in a basic domain.
124//
125// Used by gcd().
126//
127//}}}
128static CanonicalForm
129bcontent ( const CanonicalForm & f, const CanonicalForm & c )
130{
131    if ( f.inBaseDomain() )
132        return bgcd( f, c );
133    else if ( f.inCoeffDomain() )
134        return f.genOne();
135    else {
136        CanonicalForm g = c;
137        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
138            g = bcontent( i.coeff(), g );
139        if( g.lc().sign() < 0 )
140            return -g;
141        else
142            return g;
143    }
144}
145//}}}
146
147
148//{{{ CanonicalForm icontent ( const CanonicalForm & f )
149//{{{ docu
150//
151// icontent() - return gcd over all coefficients of f which are
152//   in a coefficient domain.
153//
154//}}}
155CanonicalForm
156icontent ( const CanonicalForm & f )
157{
158    return icontent( f, 0 );
159}
160//}}}
161
162//{{{ CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
163//{{{ docu
164//
165// extgcd() - returns polynomial extended gcd of f and g.
166//
167// Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g).
168// The gcd is calculated using an extended euclidean polynomial
169// remainder sequence, so f and g should be polynomials over an
170// euclidean domain.  Normalizes result.
171//
172// Note: be sure that f and g have the same level!
173//
174//}}}
175CanonicalForm
176extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
177{
178#ifdef HAVE_NTL
179  if (isOn(SW_USE_NTL_GCD_P) && ( getCharacteristic() > 0 )
180  && isPurePoly(f) && isPurePoly(g))
181  {
182    zz_pContext ccc(getCharacteristic());
183    ccc.restore();
184    zz_p::init(getCharacteristic());
185    zz_pX F1=convertFacCF2NTLzzpX(f);
186    zz_pX G1=convertFacCF2NTLzzpX(g);
187    zz_pX R;
188    zz_pX A,B;
189    XGCD(R,A,B,F1,G1);
190    a=convertNTLzzpX2CF(A,f.mvar());
191    b=convertNTLzzpX2CF(B,f.mvar());
192    return convertNTLzzpX2CF(R,f.mvar());
193  }
194#endif
195  CanonicalForm contf = content( f );
196  CanonicalForm contg = content( g );
197
198  CanonicalForm p0 = f / contf, p1 = g / contg;
199  CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
200
201  while ( ! p1.isZero() ) {
202      divrem( p0, p1, q, r );
203      p0 = p1; p1 = r;
204      r = g0 - g1 * q;
205      g0 = g1; g1 = r;
206      r = f0 - f1 * q;
207      f0 = f1; f1 = r;
208  }
209  CanonicalForm contp0 = content( p0 );
210  a = f0 / ( contf * contp0 );
211  b = g0 / ( contg * contp0 );
212  p0 /= contp0;
213  if ( p0.sign() < 0 ) {
214      p0 = -p0;
215      a = -a;
216      b = -b;
217  }
218  return p0;
219}
220//}}}
221
222static CanonicalForm
223gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
224{
225#ifdef HAVE_NTL
226  if (isOn(SW_USE_NTL_GCD_P) && isPurePoly(F) && isPurePoly(G))
227  {
228    if ( getCharacteristic() > 0 )
229    {
230      //CanonicalForm cf=F.lc();
231      //CanonicalForm f=F / cf;
232      //CanonicalForm cg=G.lc();
233      //CanonicalForm g= G / cg;
234      zz_pContext ccc(getCharacteristic());
235      ccc.restore();
236      zz_p::init(getCharacteristic());
237      zz_pX F1=convertFacCF2NTLzzpX(F);
238      zz_pX G1=convertFacCF2NTLzzpX(G);
239      zz_pX R=GCD(F1,G1);
240      return  convertNTLzzpX2CF(R,F.mvar());
241    }
242    else
243    {
244      CanonicalForm f=F ;
245      CanonicalForm g=G ;
246      bool rat=isOn( SW_RATIONAL );
247      if ( isOn( SW_RATIONAL ) )
248      {
249         DEBOUTLN( cerr, "NTL_gcd: ..." );
250         CanonicalForm cdF = bCommonDen( F );
251         CanonicalForm cdG = bCommonDen( G );
252         Off( SW_RATIONAL );
253         CanonicalForm l = lcm( cdF, cdG );
254         On( SW_RATIONAL );
255         f *= l, g *= l;
256      }
257      DEBOUTLN( cerr, "NTL_gcd: f=" << f );
258      DEBOUTLN( cerr, "NTL_gcd: g=" << g );
259      ZZX F1=convertFacCF2NTLZZX(f);
260      ZZX G1=convertFacCF2NTLZZX(g);
261      ZZX R=GCD(F1,G1);
262      CanonicalForm r=convertNTLZZX2CF(R,F.mvar());
263      DEBOUTLN( cerr, "NTL_gcd: -> " << r );
264      if (rat)
265      {
266        r /= r.lc();
267        DEBOUTLN( cerr, "NTL_gcd2: -> " << r );
268      }
269      return r;
270    }
271  }
272#endif
273  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
274  int p, i, n;
275
276  if ( primitive )
277  {
278    f = F;
279    g = G;
280    c = 1;
281  }
282  else
283  {
284    CanonicalForm cF = content( F ), cG = content( G );
285    f = F / cF;
286    g = G / cG;
287    c = bgcd( cF, cG );
288  }
289  cg = gcd( f.lc(), g.lc() );
290  cl = ( f.lc() / cg ) * g.lc();
291//     B = 2 * cg * tmin(
292//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
293//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
294//         )+1;
295  M = tmin( maxNorm(f), maxNorm(g) );
296  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
297  q = 0;
298  i = cf_getNumSmallPrimes() - 1;
299  while ( true )
300  {
301    B = BB;
302    while ( i >= 0 && q < B )
303    {
304      p = cf_getSmallPrime( i );
305      i--;
306      while ( i >= 0 && mod( cl, p ) == 0 )
307      {
308        p = cf_getSmallPrime( i );
309        i--;
310      }
311      setCharacteristic( p );
312      Dp = gcd( mapinto( f ), mapinto( g ) );
313      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
314      setCharacteristic( 0 );
315      if ( Dp.degree() == 0 )
316        return c;
317      if ( q.isZero() )
318      {
319        D = mapinto( Dp );
320        q = p;
321        B = power(CanonicalForm(2),D.degree())*M+1;
322      }
323      else
324      {
325        if ( Dp.degree() == D.degree() )
326        {
327          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
328          q = newq;
329          D = newD;
330        }
331        else if ( Dp.degree() < D.degree() )
332        {
333          // all previous p's are bad primes
334          q = p;
335          D = mapinto( Dp );
336          B = power(CanonicalForm(2),D.degree())*M+1;
337        }
338        // else p is a bad prime
339      }
340    }
341    if ( i >= 0 )
342    {
343      // now balance D mod q
344      D = pp( balance( D, q ) );
345      if ( divides( D, f ) && divides( D, g ) )
346        return D * c;
347      else
348        q = 0;
349    }
350    else
351      return gcd_poly( F, G, false );
352    DEBOUTLN( cerr, "another try ..." );
353  }
354}
355
356CanonicalForm
357gcd_poly1( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
358{
359    CanonicalForm C, Ci, Ci1, Hi, bi, pi, pi1, pi2;
360    int delta;
361    Variable v = f.mvar();
362
363    if ( f.degree( v ) >= g.degree( v ) ) {
364        pi = f; pi1 = g;
365    }
366    else {
367        pi = g; pi1 = f;
368    }
369    Ci = content( pi ); Ci1 = content( pi1 );
370    C = gcd( Ci, Ci1 );
371    pi1 = pi1 / Ci1; pi = pi / Ci;
372    if ( pi.isUnivariate() && pi1.isUnivariate() )
373    {
374#ifdef HAVE_NTL
375      if ((isOn(SW_USE_NTL_GCD_P)||isOn(SW_USE_NTL_GCD_0))
376       && isPurePoly(pi) && isPurePoly(pi1))
377         return gcd_poly_univar0(f, g, true);
378#endif
379      if ( modularflag)
380        return gcd_poly_univar0( pi, pi1, true ) * C;
381    }
382    else if ( gcd_test_one( pi1, pi, true ) )
383      return C;
384    delta = degree( pi, v ) - degree( pi1, v );
385    Hi = power( LC( pi1, v ), delta );
386    if ( (delta+1) % 2 )
387        bi = 1;
388    else
389        bi = -1;
390    while ( degree( pi1, v ) > 0 ) {
391        pi2 = psr( pi, pi1, v );
392        pi2 = pi2 / bi;
393        pi = pi1; pi1 = pi2;
394        if ( degree( pi1, v ) > 0 ) {
395            delta = degree( pi, v ) - degree( pi1, v );
396            if ( (delta+1) % 2 )
397                bi = LC( pi, v ) * power( Hi, delta );
398            else
399                bi = -LC( pi, v ) * power( Hi, delta );
400            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
401        }
402    }
403    if ( degree( pi1, v ) == 0 )
404        return C;
405    else {
406        return C * pp( pi );
407    }
408}
409
410//{{{ static CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
411//{{{ docu
412//
413// gcd_poly() - calculate polynomial gcd.
414//
415// This is the dispatcher for polynomial gcd calculation.  We call either
416// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
417// characteristic and settings of SW_USE_EZGCD and SW_USE_SPARSEMOD, resp.
418//
419// modularflag is reached down to gcd_poly1() without change in case of
420// zero characteristic.
421//
422// Used by gcd() and gcd_poly_univar0().
423//
424//}}}
425int si_factor_reminder=1;
426static CanonicalForm
427gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
428{
429    if ( getCharacteristic() != 0 ) {
430        if (! ( f.isUnivariate() && g.isUnivariate() ) ) {
431            CFMap M, N;
432            compress( f, g, M, N );
433            CanonicalForm fM = M(f);
434            CanonicalForm gM = M(g);
435            if ( fM.mvar() != gM.mvar() ) {
436                if ( fM.mvar() > gM.mvar() )
437                    return N( cf_content( fM, gM ) );
438                else
439                    return N( cf_content( gM, fM ) );
440            }
441            else
442                return N( gcd_poly1( fM, gM, false ) );
443        }
444        else
445            return gcd_poly1( f, g, false );
446    }
447    else if ( isOn( SW_USE_EZGCD ) && ! ( f.isUnivariate() && g.isUnivariate() ) ) {
448        CFMap M, N;
449        compress( f, g, M, N );
450#if 0
451        CanonicalForm r=N( ezgcd( M(f), M(g) ) );
452        if ((f%r!=0) || (g % r !=0))
453        {
454           if (si_factor_reminder)
455           printf("ezgcd failed, trying gcd_poly1\n");
456           return gcd_poly1( f, g, modularflag);
457        }
458        else
459           return r;
460#else
461         return N( ezgcd( M(f), M(g) ) );
462#endif
463    }
464    else if ( isOn( SW_USE_SPARSEMOD )
465    && ! ( f.isUnivariate() && g.isUnivariate() ) )
466    {
467#if 0
468        CanonicalForm r=sparsemod( f, g );
469        if ((f%r!=0) || (g % r !=0))
470        {
471           if (si_factor_reminder)
472           printf("sparsemod failed, trying gcd_poly1\n");
473           return r;
474           //return gcd_poly1( f, g, modularflag);
475        }
476        else
477          return r;
478#else
479        return sparsemod( f, g );
480#endif
481    }
482    else
483    {
484        return gcd_poly1( f, g, modularflag );
485    }
486}
487//}}}
488
489//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
490//{{{ docu
491//
492// cf_content() - return gcd(g, content(f)).
493//
494// content(f) is calculated with respect to f's main variable.
495//
496// Used by gcd(), content(), content( CF, Variable ).
497//
498//}}}
499static CanonicalForm
500cf_content ( const CanonicalForm & f, const CanonicalForm & g )
501{
502    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
503        CFIterator i = f;
504        CanonicalForm result = g;
505        while ( i.hasTerms() && ! result.isOne() ) {
506            result = gcd( result, i.coeff() );
507            i++;
508        }
509        return result;
510    }
511    else
512    {
513        ASSERT(g==0,"invalid call of cf_gcd");
514        if ( f.sign() < 0 )
515            return -f;
516        else
517            return f;
518    }
519}
520//}}}
521
522//{{{ CanonicalForm content ( const CanonicalForm & f )
523//{{{ docu
524//
525// content() - return content(f) with respect to main variable.
526//
527// Normalizes result.
528//
529//}}}
530CanonicalForm
531content ( const CanonicalForm & f )
532{
533    return cf_content( f, 0 );
534}
535//}}}
536
537//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
538//{{{ docu
539//
540// content() - return content(f) with respect to x.
541//
542// x should be a polynomial variable.
543//
544//}}}
545CanonicalForm
546content ( const CanonicalForm & f, const Variable & x )
547{
548    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
549    Variable y = f.mvar();
550
551    if ( y == x )
552        return cf_content( f, 0 );
553    else  if ( y < x )
554        return f;
555    else
556        return swapvar( content( swapvar( f, y, x ), y ), y, x );
557}
558//}}}
559
560//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
561//{{{ docu
562//
563// vcontent() - return content of f with repect to variables >= x.
564//
565// The content is recursively calculated over all coefficients in
566// f having level less than x.  x should be a polynomial
567// variable.
568//
569//}}}
570CanonicalForm
571vcontent ( const CanonicalForm & f, const Variable & x )
572{
573    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
574
575    if ( f.mvar() <= x )
576        return content( f, x );
577    else {
578        CFIterator i;
579        CanonicalForm d = 0;
580        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
581            d = gcd( d, vcontent( i.coeff(), x ) );
582        return d;
583    }
584}
585//}}}
586
587//{{{ CanonicalForm pp ( const CanonicalForm & f )
588//{{{ docu
589//
590// pp() - return primitive part of f.
591//
592// Returns zero if f equals zero, otherwise f / content(f).
593//
594//}}}
595CanonicalForm
596pp ( const CanonicalForm & f )
597{
598    if ( f.isZero() )
599        return f;
600    else
601        return f / content( f );
602}
603//}}}
604
605CanonicalForm
606gcd ( const CanonicalForm & f, const CanonicalForm & g )
607{
608    if ( f.isZero() )
609    {
610      if ( g.lc().sign() < 0 )
611        return -g;
612      else
613       return g;
614    }
615    else  if ( g.isZero() )
616    {
617      if ( f.lc().sign() < 0 )
618        return -f;
619      else
620        return f;
621    }
622    else  if ( f.inBaseDomain() )
623        return bcontent( g, f );
624    else  if ( g.inBaseDomain() )
625        return bcontent( f, g );
626    else  if ( f.mvar() == g.mvar() )
627        if ( f.inExtension() && getReduce( f.mvar() ) )
628            return 1;
629        else {
630            if ( divides( f, g ) )
631                if ( f.lc().sign() < 0 )
632                    return -f;
633                else
634                    return f;
635            else  if ( divides( g, f ) )
636                if ( g.lc().sign() < 0 )
637                    return -g;
638                else
639                    return g;
640            if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) {
641                CanonicalForm cdF = bCommonDen( f );
642                CanonicalForm cdG = bCommonDen( g );
643                Off( SW_RATIONAL );
644                CanonicalForm l = lcm( cdF, cdG );
645                On( SW_RATIONAL );
646                CanonicalForm F = f * l, G = g * l;
647                Off( SW_RATIONAL );
648                do { l = gcd_poly( F, G, true ); 
649                //if ((!divides(l,F)) || (!divides(l,G)))
650                //{printf("fehler\n");
651                //cout<<"F="<<F<<"\nG="<<G<<"wrong gcd="<<l<<"\n";
652                //}
653                }
654                while (0); //((!divides(l,F)) || (!divides(l,G)));
655                On( SW_RATIONAL );
656                if ( l.lc().sign() < 0 )
657                    return -l;
658                else
659                    return l;
660            }
661            else {
662                CanonicalForm d;
663                do{ d = gcd_poly( f, g, getCharacteristic()==0 ); 
664                //if ((!divides(d,f)) || (!divides(d,g)))
665                //{printf("fehler2\n"); cout<<"F="<<f<<"\nG="<<g<<"wrong gcd="<<d<<"\n";
666                //}
667                }
668                while (0);//((!divides(d,f)) || (!divides(d,g)));
669                if ( d.lc().sign() < 0 )
670                    return -d;
671                else
672                    return d;
673            }
674        }
675    else  if ( f.mvar() > g.mvar() )
676        return cf_content( f, g );
677    else
678        return cf_content( g, f );
679}
680
681//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
682//{{{ docu
683//
684// lcm() - return least common multiple of f and g.
685//
686// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
687//
688// Returns zero if one of f or g equals zero.
689//
690//}}}
691CanonicalForm
692lcm ( const CanonicalForm & f, const CanonicalForm & g )
693{
694    if ( f.isZero() || g.isZero() )
695        return f;
696    else
697        return ( f / gcd( f, g ) ) * g;
698}
699//}}}
Note: See TracBrowser for help on using the repository browser.