source: git/factory/cf_gcd.cc @ 845303c

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