source: git/factory/cf_gcd.cc @ 767da0

spielwiese
Last change on this file since 767da0 was 4b24c19, checked in by Hans Schoenemann <hannes@…>, 4 years ago
Revert "use: chineseRemainderCached" This reverts commit 0aeeeeaed39a3ac804f59d6a4f58df7e1572b687.
  • Property mode set to 100644
File size: 18.8 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3/**
4 * @file cf_gcd.cc
5 *
6 * gcd/content/lcm of polynomials
7 *
8 * To compute the GCD different variants are chosen automatically
9**/
10
11#include "config.h"
12
13
14#include "timing.h"
15#include "cf_assert.h"
16#include "debug.h"
17
18#include "cf_defs.h"
19#include "canonicalform.h"
20#include "cf_iter.h"
21#include "cf_reval.h"
22#include "cf_primes.h"
23#include "cf_algorithm.h"
24#include "cfEzgcd.h"
25#include "cfGcdAlgExt.h"
26#include "cfSubResGcd.h"
27#include "cfModGcd.h"
28#include "FLINTconvert.h"
29#include "facAlgFuncUtil.h"
30#include "templates/ftmpl_functions.h"
31
32#ifdef HAVE_NTL
33#include <NTL/ZZX.h>
34#include "NTLconvert.h"
35bool isPurePoly(const CanonicalForm & );
36#endif
37
38void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
39
40/** static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
41 *
42 * icontent() - return gcd of c and all coefficients of f which
43 *   are in a coefficient domain.
44 *
45 * @sa icontent().
46 *
47**/
48static CanonicalForm
49icontent ( const CanonicalForm & f, const CanonicalForm & c )
50{
51    if ( f.inBaseDomain() )
52    {
53      if (c.isZero()) return abs(f);
54      return bgcd( f, c );
55    }
56    //else if ( f.inCoeffDomain() )
57    //   return gcd(f,c);
58    else
59    {
60        CanonicalForm g = c;
61        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
62            g = icontent( i.coeff(), g );
63        return g;
64    }
65}
66
67/** CanonicalForm icontent ( const CanonicalForm & f )
68 *
69 * icontent() - return gcd over all coefficients of f which are
70 *   in a coefficient domain.
71 *
72**/
73CanonicalForm
74icontent ( const CanonicalForm & f )
75{
76    return icontent( f, 0 );
77}
78
79#ifdef HAVE_FLINT
80static CanonicalForm
81gcd_univar_flintp (const CanonicalForm& F, const CanonicalForm& G)
82{
83  nmod_poly_t F1, G1;
84  convertFacCF2nmod_poly_t (F1, F);
85  convertFacCF2nmod_poly_t (G1, G);
86  nmod_poly_gcd (F1, F1, G1);
87  CanonicalForm result= convertnmod_poly_t2FacCF (F1, F.mvar());
88  nmod_poly_clear (F1);
89  nmod_poly_clear (G1);
90  return result;
91}
92
93static CanonicalForm
94gcd_univar_flint0( const CanonicalForm & F, const CanonicalForm & G )
95{
96  fmpz_poly_t F1, G1;
97  convertFacCF2Fmpz_poly_t(F1, F);
98  convertFacCF2Fmpz_poly_t(G1, G);
99  fmpz_poly_gcd (F1, F1, G1);
100  CanonicalForm result= convertFmpz_poly_t2FacCF (F1, F.mvar());
101  fmpz_poly_clear (F1);
102  fmpz_poly_clear (G1);
103  return result;
104}
105#endif
106
107#ifdef HAVE_NTL
108#ifndef HAVE_FLINT
109static CanonicalForm
110gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
111{
112    ZZX F1=convertFacCF2NTLZZX(F);
113    ZZX G1=convertFacCF2NTLZZX(G);
114    ZZX R=GCD(F1,G1);
115    return convertNTLZZX2CF(R,F.mvar());
116}
117
118static CanonicalForm
119gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
120{
121  int ch=getCharacteristic();
122  if (fac_NTL_char!=ch)
123  {
124    fac_NTL_char=ch;
125    zz_p::init(ch);
126  }
127  zz_pX F1=convertFacCF2NTLzzpX(F);
128  zz_pX G1=convertFacCF2NTLzzpX(G);
129  zz_pX R=GCD(F1,G1);
130  return  convertNTLzzpX2CF(R,F.mvar());
131}
132#endif
133#endif
134
135//{{{ static CanonicalForm balance_p ( const CanonicalForm & f, const CanonicalForm & q )
136//{{{ docu
137//
138// balance_p() - map f from positive to symmetric representation
139//   mod q.
140//
141// This makes sense for univariate polynomials over Z only.
142// q should be an integer.
143//
144// Used by gcd_poly_univar0().
145//
146//}}}
147
148static CanonicalForm
149balance_p ( const CanonicalForm & f, const CanonicalForm & q, const CanonicalForm & qh )
150{
151    Variable x = f.mvar();
152    CanonicalForm result = 0;
153    CanonicalForm c;
154    CFIterator i;
155    for ( i = f; i.hasTerms(); i++ )
156    {
157        c = i.coeff();
158        if ( c.inCoeffDomain())
159        {
160          if ( c > qh )
161            result += power( x, i.exp() ) * (c - q);
162          else
163            result += power( x, i.exp() ) * c;
164        }
165        else
166          result += power( x, i.exp() ) * balance_p(c,q,qh);
167    }
168    return result;
169}
170
171static CanonicalForm
172balance_p ( const CanonicalForm & f, const CanonicalForm & q )
173{
174    CanonicalForm qh = q / 2;
175    return balance_p (f, q, qh);
176}
177
178static CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
179{
180  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
181  int p, i;
182
183  if ( primitive )
184  {
185    f = F;
186    g = G;
187    c = 1;
188  }
189  else
190  {
191    CanonicalForm cF = content( F ), cG = content( G );
192    f = F / cF;
193    g = G / cG;
194    c = bgcd( cF, cG );
195  }
196  cg = gcd( f.lc(), g.lc() );
197  cl = ( f.lc() / cg ) * g.lc();
198//     B = 2 * cg * tmin(
199//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
200//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
201//         )+1;
202  M = tmin( maxNorm(f), maxNorm(g) );
203  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
204  q = 0;
205  i = cf_getNumSmallPrimes() - 1;
206  while ( true )
207  {
208    B = BB;
209    while ( i >= 0 && q < B )
210    {
211      p = cf_getSmallPrime( i );
212      i--;
213      while ( i >= 0 && mod( cl, p ) == 0 )
214      {
215        p = cf_getSmallPrime( i );
216        i--;
217      }
218      setCharacteristic( p );
219      Dp = gcd( mapinto( f ), mapinto( g ) );
220      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
221      setCharacteristic( 0 );
222      if ( Dp.degree() == 0 )
223        return c;
224      if ( q.isZero() )
225      {
226        D = mapinto( Dp );
227        q = p;
228        B = power(CanonicalForm(2),D.degree())*M+1;
229      }
230      else
231      {
232        if ( Dp.degree() == D.degree() )
233        {
234          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
235          q = newq;
236          D = newD;
237        }
238        else if ( Dp.degree() < D.degree() )
239        {
240          // all previous p's are bad primes
241          q = p;
242          D = mapinto( Dp );
243          B = power(CanonicalForm(2),D.degree())*M+1;
244        }
245        // else p is a bad prime
246      }
247    }
248    if ( i >= 0 )
249    {
250      // now balance D mod q
251      D = pp( balance_p( D, q ) );
252      if ( fdivides( D, f ) && fdivides( D, g ) )
253        return D * c;
254      else
255        q = 0;
256    }
257    else
258      return gcd_poly( F, G );
259    DEBOUTLN( cerr, "another try ..." );
260  }
261}
262
263static CanonicalForm
264gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
265{
266    if (f.inCoeffDomain() || g.inCoeffDomain()) //zero case should be caught by gcd
267      return 1;
268    CanonicalForm pi, pi1;
269    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
270    bool bpure, ezgcdon= isOn (SW_USE_EZGCD_P);
271    int delta = degree( f ) - degree( g );
272
273    if ( delta >= 0 )
274    {
275        pi = f; pi1 = g;
276    }
277    else
278    {
279        pi = g; pi1 = f; delta = -delta;
280    }
281    if (pi.isUnivariate())
282      Ci= 1;
283    else
284    {
285      if (!ezgcdon)
286        On (SW_USE_EZGCD_P);
287      Ci = content( pi );
288      if (!ezgcdon)
289        Off (SW_USE_EZGCD_P);
290      pi = pi / Ci;
291    }
292    if (pi1.isUnivariate())
293      Ci1= 1;
294    else
295    {
296      if (!ezgcdon)
297        On (SW_USE_EZGCD_P);
298      Ci1 = content( pi1 );
299      if (!ezgcdon)
300        Off (SW_USE_EZGCD_P);
301      pi1 = pi1 / Ci1;
302    }
303    C = gcd( Ci, Ci1 );
304    int d= 0;
305    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
306    {
307        if ( gcd_test_one( pi1, pi, true, d ) )
308        {
309          C=abs(C);
310          //out_cf("GCD:",C,"\n");
311          return C;
312        }
313        bpure = false;
314    }
315    else
316    {
317        bpure = isPurePoly(pi) && isPurePoly(pi1);
318#ifdef HAVE_FLINT
319        if (bpure && (CFFactory::gettype() != GaloisFieldDomain))
320          return gcd_univar_flintp(pi,pi1)*C;
321#else
322#ifdef HAVE_NTL
323        if ( bpure && (CFFactory::gettype() != GaloisFieldDomain))
324            return gcd_univar_ntlp(pi, pi1 ) * C;
325#endif
326#endif
327    }
328    Variable v = f.mvar();
329    Hi = power( LC( pi1, v ), delta );
330    int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
331
332    if (!(pi.isUnivariate() && pi1.isUnivariate()))
333    {
334      if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning
335      {
336        On (SW_USE_FF_MOD_GCD);
337        C *= gcd (pi, pi1);
338        Off (SW_USE_FF_MOD_GCD);
339        return C;
340      }
341    }
342
343    if ( (delta+1) % 2 )
344        bi = 1;
345    else
346        bi = -1;
347    CanonicalForm oldPi= pi, oldPi1= pi1, powHi;
348    while ( degree( pi1, v ) > 0 )
349    {
350        if (!(pi.isUnivariate() && pi1.isUnivariate()))
351        {
352          if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500)
353          {
354            On (SW_USE_FF_MOD_GCD);
355            C *= gcd (oldPi, oldPi1);
356            Off (SW_USE_FF_MOD_GCD);
357            return C;
358          }
359        }
360        pi2 = psr( pi, pi1, v );
361        pi2 = pi2 / bi;
362        pi = pi1; pi1 = pi2;
363        maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
364        if (!pi1.isUnivariate() && (size (pi1)/maxNumVars > 500))
365        {
366            On (SW_USE_FF_MOD_GCD);
367            C *= gcd (oldPi, oldPi1);
368            Off (SW_USE_FF_MOD_GCD);
369            return C;
370        }
371        if ( degree( pi1, v ) > 0 )
372        {
373            delta = degree( pi, v ) - degree( pi1, v );
374            powHi= power (Hi, delta-1);
375            if ( (delta+1) % 2 )
376                bi = LC( pi, v ) * powHi*Hi;
377            else
378                bi = -LC( pi, v ) * powHi*Hi;
379            Hi = power( LC( pi1, v ), delta ) / powHi;
380            if (!(pi.isUnivariate() && pi1.isUnivariate()))
381            {
382              if (size (Hi)*size (pi)/(maxNumVars*3) > 1500) //maybe this needs more tuning
383              {
384                On (SW_USE_FF_MOD_GCD);
385                C *= gcd (oldPi, oldPi1);
386                Off (SW_USE_FF_MOD_GCD);
387                return C;
388              }
389            }
390        }
391    }
392    if ( degree( pi1, v ) == 0 )
393    {
394      C=abs(C);
395      //out_cf("GCD:",C,"\n");
396      return C;
397    }
398    if (!pi.isUnivariate())
399    {
400      if (!ezgcdon)
401        On (SW_USE_EZGCD_P);
402      Ci= gcd (LC (oldPi,v), LC (oldPi1,v));
403      pi /= LC (pi,v)/Ci;
404      Ci= content (pi);
405      pi /= Ci;
406      if (!ezgcdon)
407        Off (SW_USE_EZGCD_P);
408    }
409    if ( bpure )
410        pi /= pi.lc();
411    C=abs(C*pi);
412    //out_cf("GCD:",C,"\n");
413    return C;
414}
415
416static CanonicalForm
417gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
418{
419    CanonicalForm pi, pi1;
420    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
421    int delta = degree( f ) - degree( g );
422
423    if ( delta >= 0 )
424    {
425        pi = f; pi1 = g;
426    }
427    else
428    {
429        pi = g; pi1 = f; delta = -delta;
430    }
431    Ci = content( pi ); Ci1 = content( pi1 );
432    pi1 = pi1 / Ci1; pi = pi / Ci;
433    C = gcd( Ci, Ci1 );
434    int d= 0;
435    if ( pi.isUnivariate() && pi1.isUnivariate() )
436    {
437#ifdef HAVE_FLINT
438        if (isPurePoly(pi) && isPurePoly(pi1) )
439            return gcd_univar_flint0(pi, pi1 ) * C;
440#else
441#ifdef HAVE_NTL
442        if ( isPurePoly(pi) && isPurePoly(pi1) )
443            return gcd_univar_ntl0(pi, pi1 ) * C;
444#endif
445#endif
446        return gcd_poly_univar0( pi, pi1, true ) * C;
447    }
448    else if ( gcd_test_one( pi1, pi, true, d ) )
449      return C;
450    Variable v = f.mvar();
451    Hi = power( LC( pi1, v ), delta );
452    if ( (delta+1) % 2 )
453        bi = 1;
454    else
455        bi = -1;
456    while ( degree( pi1, v ) > 0 )
457    {
458        pi2 = psr( pi, pi1, v );
459        pi2 = pi2 / bi;
460        pi = pi1; pi1 = pi2;
461        if ( degree( pi1, v ) > 0 )
462        {
463            delta = degree( pi, v ) - degree( pi1, v );
464            if ( (delta+1) % 2 )
465                bi = LC( pi, v ) * power( Hi, delta );
466            else
467                bi = -LC( pi, v ) * power( Hi, delta );
468            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
469        }
470    }
471    if ( degree( pi1, v ) == 0 )
472        return C;
473    else
474        return C * pp( pi );
475}
476
477/** CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
478 *
479 * gcd_poly() - calculate polynomial gcd.
480 *
481 * This is the dispatcher for polynomial gcd calculation.
482 * Different gcd variants get called depending the input, characteristic, and
483 * on switches (cf_defs.h)
484 *
485 * With the current settings from Singular (i.e. SW_USE_EZGCD= on,
486 * SW_USE_EZGCD_P= on, SW_USE_CHINREM_GCD= on, the EZ GCD variants are the
487 * default algorithms for multivariate polynomial GCD computations)
488 *
489 * @sa gcd(), cf_defs.h
490 *
491**/
492CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
493{
494  CanonicalForm fc, gc;
495  bool fc_isUnivariate=f.isUnivariate();
496  bool gc_isUnivariate=g.isUnivariate();
497  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
498  fc = f;
499  gc = g;
500  int ch=getCharacteristic();
501  if ( ch != 0 )
502  {
503    if (0) {} // dummy, to be able to build without NTL and FLINT
504    #if defined(HAVE_FLINT) && ( __FLINT_RELEASE >= 20503)
505    if ( isOn( SW_USE_FL_GCD_P)
506    && (CFFactory::gettype() != GaloisFieldDomain)
507    #ifdef HAVE_NTL
508    && (ch>10) // if we have NTL: it is better for char <11
509    #endif
510    &&(!hasAlgVar(fc)) && (!hasAlgVar(gc)))
511    {
512      return gcdFlintMP_Zp(fc,gc);
513    }
514    #endif
515    #ifdef HAVE_NTL
516    if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P )))
517    {
518      fc= EZGCD_P (fc, gc);
519    }
520    #endif
521    #if defined(HAVE_NTL) || defined(HAVE_FLINT)
522    else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate)
523    {
524      Variable a;
525      if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a))
526        fc=modGCDFq (fc, gc, a);
527      else if (CFFactory::gettype() == GaloisFieldDomain)
528        fc=modGCDGF (fc, gc);
529      else
530        fc=modGCDFp (fc, gc);
531    }
532    #endif
533    else
534    fc = gcd_poly_p( fc, gc );
535  }
536  else if (!fc_and_gc_Univariate) /* && char==0*/
537  {
538    #if defined(HAVE_FLINT) && ( __FLINT_RELEASE >= 20503)
539    if (( isOn( SW_USE_FL_GCD_0) )
540    &&(!hasAlgVar(fc)) && (!hasAlgVar(gc)))
541    {
542      return gcdFlintMP_QQ(fc,gc);
543    }
544    else
545    #endif
546    #ifdef HAVE_NTL
547    if ( isOn( SW_USE_EZGCD ) )
548      fc= ezgcd (fc, gc);
549    else
550    #endif
551    #if defined(HAVE_NTL) || defined(HAVE_FLINT)
552    if (isOn(SW_USE_CHINREM_GCD))
553      fc = modGCDZ( fc, gc);
554    else
555    #endif
556    {
557       fc = gcd_poly_0( fc, gc );
558    }
559  }
560  else
561  {
562    fc = gcd_poly_0( fc, gc );
563  }
564  if ((ch>0)&&(!hasAlgVar(fc))) fc/=fc.lc();
565  return fc;
566}
567
568/** static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
569 *
570 * cf_content() - return gcd(g, content(f)).
571 *
572 * content(f) is calculated with respect to f's main variable.
573 *
574 * @sa gcd(), content(), content( CF, Variable ).
575 *
576**/
577static CanonicalForm
578cf_content ( const CanonicalForm & f, const CanonicalForm & g )
579{
580    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
581    {
582        CFIterator i = f;
583        CanonicalForm result = g;
584        while ( i.hasTerms() && ! result.isOne() )
585        {
586            result = gcd( i.coeff(), result );
587            i++;
588        }
589        return result;
590    }
591    else
592        return abs( f );
593}
594
595/** CanonicalForm content ( const CanonicalForm & f )
596 *
597 * content() - return content(f) with respect to main variable.
598 *
599 * Normalizes result.
600 *
601**/
602CanonicalForm
603content ( const CanonicalForm & f )
604{
605    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
606    {
607        CFIterator i = f;
608        CanonicalForm result = abs( i.coeff() );
609        i++;
610        while ( i.hasTerms() && ! result.isOne() )
611        {
612            result = gcd( i.coeff(), result );
613            i++;
614        }
615        return result;
616    }
617    else
618        return abs( f );
619}
620
621/** CanonicalForm content ( const CanonicalForm & f, const Variable & x )
622 *
623 * content() - return content(f) with respect to x.
624 *
625 * x should be a polynomial variable.
626 *
627**/
628CanonicalForm
629content ( const CanonicalForm & f, const Variable & x )
630{
631    if (f.inBaseDomain()) return f;
632    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
633    Variable y = f.mvar();
634
635    if ( y == x )
636        return cf_content( f, 0 );
637    else  if ( y < x )
638        return f;
639    else
640        return swapvar( content( swapvar( f, y, x ), y ), y, x );
641}
642
643/** CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
644 *
645 * vcontent() - return content of f with repect to variables >= x.
646 *
647 * The content is recursively calculated over all coefficients in
648 * f having level less than x.  x should be a polynomial
649 * variable.
650 *
651**/
652CanonicalForm
653vcontent ( const CanonicalForm & f, const Variable & x )
654{
655    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
656
657    if ( f.mvar() <= x )
658        return content( f, x );
659    else {
660        CFIterator i;
661        CanonicalForm d = 0;
662        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
663            d = gcd( d, vcontent( i.coeff(), x ) );
664        return d;
665    }
666}
667
668/** CanonicalForm pp ( const CanonicalForm & f )
669 *
670 * pp() - return primitive part of f.
671 *
672 * Returns zero if f equals zero, otherwise f / content(f).
673 *
674**/
675CanonicalForm
676pp ( const CanonicalForm & f )
677{
678    if ( f.isZero() )
679        return f;
680    else
681        return f / content( f );
682}
683
684CanonicalForm
685gcd ( const CanonicalForm & f, const CanonicalForm & g )
686{
687    bool b = f.isZero();
688    if ( b || g.isZero() )
689    {
690        if ( b )
691            return abs( g );
692        else
693            return abs( f );
694    }
695    if ( f.inPolyDomain() || g.inPolyDomain() )
696    {
697        if ( f.mvar() != g.mvar() )
698        {
699            if ( f.mvar() > g.mvar() )
700                return cf_content( f, g );
701            else
702                return cf_content( g, f );
703        }
704        if (isOn(SW_USE_QGCD))
705        {
706          Variable m;
707          if (
708          (getCharacteristic() == 0) &&
709          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
710          )
711          {
712            bool on_rational = isOn(SW_RATIONAL);
713            CanonicalForm r=QGCD(f,g);
714            On(SW_RATIONAL);
715            CanonicalForm cdF = bCommonDen( r );
716            if (!on_rational) Off(SW_RATIONAL);
717            return cdF*r;
718          }
719        }
720
721        if ( f.inExtension() && getReduce( f.mvar() ) )
722            return CanonicalForm(1);
723        else
724        {
725            if ( fdivides( f, g ) )
726                return abs( f );
727            else  if ( fdivides( g, f ) )
728                return abs( g );
729            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
730            {
731                CanonicalForm d;
732                d = gcd_poly( f, g );
733                return abs( d );
734            }
735            else
736            {
737                CanonicalForm cdF = bCommonDen( f );
738                CanonicalForm cdG = bCommonDen( g );
739                CanonicalForm F = f * cdF, G = g * cdG;
740                Off( SW_RATIONAL );
741                CanonicalForm l = gcd_poly( F, G );
742                On( SW_RATIONAL );
743                return abs( l );
744            }
745        }
746    }
747    if ( f.inBaseDomain() && g.inBaseDomain() )
748        return bgcd( f, g );
749    else
750        return 1;
751}
752
753/** CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
754 *
755 * lcm() - return least common multiple of f and g.
756 *
757 * The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
758 *
759 * Returns zero if one of f or g equals zero.
760 *
761**/
762CanonicalForm
763lcm ( const CanonicalForm & f, const CanonicalForm & g )
764{
765    if ( f.isZero() || g.isZero() )
766        return 0;
767    else
768        return ( f / gcd( f, g ) ) * g;
769}
770
Note: See TracBrowser for help on using the repository browser.