src/rational/elem/cl_RA_mul.cc

   1 // binary operator *
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/rational.h"
   8
   9
  10 // Implementation.
  11
  12 #include "rational/cl_RA.h"
  13 #include "cln/integer.h"
  14 #include "integer/cl_I.h"
  15
  16 namespace cln {
  17
  18 const cl_RA operator* (const cl_RA& r, const cl_RA& s)
  19 {
  20 // Methode (vgl. [Buchberger, Collins, Loos: Computer Algebra, S.201])
  21 // r,s beide Integers -> klar.
  22 // r=a/b, s=c ->
  23 //   Bei c=0 Ergebnis 0.
  24 //   g:=ggT(b,c).
  25 //   Falls g=1: Ergebnis (a*c)/b (mit b>1, ggT(a*c,b)=1).
  26 //   Sonst: b':=b/g, c':=c/g, Ergebnis (a*c')/b' (mit ggT(a*c',b')=1).
  27 // r=a, s=c/d analog.
  28 // r=a/b, s=c/d ->
  29 //   g:=ggT(a,d), h:=ggT(b,c).
  30 //   a':=a/g, d':=d/g (nur bei g>1 bedeutet das Rechnung).
  31 //   b':=b/h, c':=c/h (nur bei h>1 bedeutet das Rechnung).
  32 //   Ergebnis ist = (a'*c')/(b'*d').
  33         if (integerp(s)) {
  34                 // s Integer
  35                 DeclareType(cl_I,s);
  36                 if (integerp(r)) {
  37                         // beides Integer
  38                         DeclareType(cl_I,r);
  39                         return r*s;
  40                 } else {
  41                         DeclareType(cl_RT,r);
  42                         var const cl_I& a = numerator(r);
  43                         var const cl_I& b = denominator(r);
  44                         var const cl_I& c = s;
  45                         // r=a/b, s=c, bilde a/b * c.
  46                         if (zerop(c))
  47                                 { return 0; } // c=0 -> Ergebnis 0
  48                         var cl_I g = gcd(b,c);
  49                         if (eq(g,1))
  50                                 // g=1
  51                                 return I_I_to_RT(a*c,b); // (a*c)/b
  52                         else
  53                                 // g>1
  54                                 return I_I_to_RA(a*exquo(c,g),exquopos(b,g)); // (a*(c/g))/(b/g)
  55                 }
  56         } else {
  57                 // s ist Ratio
  58                 DeclareType(cl_RT,s);
  59                 if (integerp(r)) {
  60                         // r Integer
  61                         DeclareType(cl_I,r);
  62                         var const cl_I& a = r;
  63                         var const cl_I& b = numerator(s);
  64                         var const cl_I& c = denominator(s);
  65                         // r=a, s=b/c, bilde a * b/c.
  66                         if (zerop(a))
  67                                 { return 0; } // a=0 -> Ergebnis 0
  68                         var cl_I g = gcd(a,c);
  69                         if (eq(g,1))
  70                                 // g=1
  71                                 return I_I_to_RT(a*b,c); // (a*b)/c
  72                         else
  73                                 // g>1
  74                                 return I_I_to_RA(exquo(a,g)*b,exquopos(c,g)); // ((a/g)*b)/(c/g)
  75                 } else {
  76                         // r,s beide Ratios
  77                         DeclareType(cl_RT,r);
  78                         var const cl_I& a = numerator(r);
  79                         var const cl_I& b = denominator(r);
  80                         var const cl_I& c = numerator(s);
  81                         var const cl_I& d = denominator(s);
  82                         var cl_I ap, dp;
  83                         {
  84                                 var cl_I g = gcd(a,d);
  85                                 if (eq(g,1))
  86                                         { ap = a; dp = d; }
  87                                 else
  88                                         { ap = exquo(a,g); dp = exquopos(d,g); }
  89                         }
  90                         var cl_I cp, bp;
  91                         {
  92                                 var cl_I h = gcd(b,c);
  93                                 if (eq(h,1))
  94                                         { cp = c; bp = b; }
  95                                 else
  96                                         { cp = exquo(c,h); bp = exquopos(b,h); }
  97                         }
  98                         return I_I_to_RA(ap*cp,bp*dp); // (a'*c')/(b'*d')
  99                 }
 100         }
 101 }
 102
 103 }  // namespace cln