src/rational/elem/cl_RA_minus.cc

   1 // binary operator -
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/rational.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_RA.h"
  13 #include "cln/integer.h"
  14 #include "cl_I.h"
  15
  16 namespace cln {
  17
  18 const cl_RA operator- (const cl_RA& r, const cl_RA& s)
  19 {
  20 #if 0
  21 // Methode:
  22 // (+ r (- s))
  23         return r + (- s);
  24 #else
  25 // Methode (vgl. [Buchberger, Collins, Loos: Computer Algebra, S.200-201])
  26 // r,s beide Integers -> klar.
  27 // r=a/b, s=c -> Ergebnis (a-b*c)/b
  28 //   (mit b>1 und ggT(a-b*c,b) = ggT(a,b) = 1)
  29 //   Bei c=0 direkt r als Ergebnis.
  30 // r=a, s=c/d -> Ergebnis (a*d-c)/d
  31 //   (mit d>1 und ggT(a*d-c,d) = ggT(-c,d) = ggT(c,d) = 1)
  32 //   Bei a=0 direkt -s = (-c)/d als Ergebnis.
  33 // r=a/b, s=c/d:
  34 //   g:=ggT(b,d)>0.
  35 //   Falls g=1:
  36 //     Ergebnis (a*d-b*c)/(b*d),
  37 //     (mit b*d>1 wegen b>1, d>1, und
  38 //      ggT(a*d-b*c,b*d) = 1
  39 //      wegen ggT(a*d-b*c,b) = ggT(a*d,b) = 1 (wegen ggT(a,b)=1 und ggT(d,b)=1)
  40 //      und   ggT(a*d-b*c,d) = ggT(b*c,d) = 1 (wegen ggT(b,d)=1 und ggT(c,d)=1)
  41 //     )
  42 //   Sonst b' := b/g, d' := d/g. e := a*d'-b'*c, f:= b'*d = b*d'.
  43 //   Es ist g = ggT(g*b',g*d') = g*ggT(b',d'), also ggT(b',d')=1.
  44 //   Es ist r-s = (a*d-b*c)/(b*d) = (nach Kürzen mit g) e/f.
  45 //   Außerdem:
  46 //     ggT(a,b') teilt ggT(a,b)=1, also ggT(a,b')=1. Mit ggT(d',b')=1 folgt
  47 //     1 = ggT(a*d',b') = ggT(a*d'-b'*c,b') = ggT(e,b').
  48 //     ggT(c,d') teilt ggT(c,d)=1, also ggT(c,d')=1. Mit ggT(b',d')=1 folgt
  49 //     1 = ggT(b'*c,d') = ggT(a*d'-b'*c,d') = ggT(e,d').
  50 //     Daher ist ggT(e,f) = ggT(e,b'*d'*g) = ggT(e,g).
  51 //   Errechne daher h=ggT(e,g).
  52 //   Bei h=1 ist e/f das Ergebnis (mit f>1, da d>1, und ggT(e,f)=1),
  53 //   sonst ist (e/h)/(f/h) das Ergebnis.
  54         if (integerp(s)) {
  55                 // s ist Integer
  56                 DeclareType(cl_I,s);
  57                 if (eq(s,0)) { return r; } // s=0 -> r als Ergebnis
  58                 if (integerp(r)) {
  59                         // beides Integers
  60                         DeclareType(cl_I,r);
  61                         return r-s;
  62                 } else {
  63                         DeclareType(cl_RT,r);
  64                         var const cl_I& a = numerator(r);
  65                         var const cl_I& b = denominator(r);
  66                         var const cl_I& c = s;
  67                         // r = a/b, s = c.
  68                         return I_I_to_RT(a-b*c,b);
  69                 }
  70         } else {
  71                 // s ist Ratio
  72                 DeclareType(cl_RT,s);
  73                 if (integerp(r)) {
  74                         // r ist Integer
  75                         DeclareType(cl_I,r);
  76                         if (eq(r,0)) {
  77                                 // r=0 -> -s als Ergebnis
  78                                 var const cl_I& c = numerator(s);
  79                                 var const cl_I& d = denominator(s);
  80                                 return I_I_to_RT(-c,d);
  81                         }
  82                         var const cl_I& a = r;
  83                         var const cl_I& c = numerator(s);
  84                         var const cl_I& d = denominator(s);
  85                         // r = a, s = c/d.
  86                         return I_I_to_RT(a*d-c,d);
  87                 } else {
  88                         // r,s beide Ratios
  89                         DeclareType(cl_RT,r);
  90                         var const cl_I& a = numerator(r);
  91                         var const cl_I& b = denominator(r);
  92                         var const cl_I& c = numerator(s);
  93                         var const cl_I& d = denominator(s);
  94                         var cl_I g = gcd(b,d); // g = ggT(b,d) >0 bilden
  95                         if (eq(g,1))
  96                                 // g=1 -> Ergebnis (a*d-b*c)/(b*d)
  97                                 return I_I_to_RT(a*d-b*c,b*d);
  98                         // g>1
  99                         var cl_I bp = exquopos(b,g); // b' := b/g (b,g>0)
 100                         var cl_I dp = exquopos(d,g); // d' := d/g (d,g>0)
 101                         var cl_I e = a*dp-bp*c; // e := a*d'-b'*c
 102                         var cl_I f = bp*d; // f := b'*d
 103                         var cl_I h = gcd(e,g); // h := ggT(e,g)
 104                         if (eq(h,1))
 105                                 // h=1
 106                                 return I_I_to_RT(e,f);
 107                         // h>1
 108                         return I_I_to_RA(exquo(e,h),exquopos(f,h)); // (e/h)/(f/h) als Ergebnis
 109                 }
 110         }
 111 #endif
 112 }
 113
 114 }  // namespace cln