/cln-1.3.2/src/rational/elem/cl_RA_minus.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
 16namespace cln {
 17
 18const 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?&#x;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