1 /*
2 Copyright (C) 1988 Free Software Foundation
3 written by Doug Lea (dl@rocky.oswego.edu)
4
5 This file is part of GNU CC.
6
7 GNU CC is distributed in the hope that it will be useful,
8 but WITHOUT ANY WARRANTY. No author or distributor
9 accepts responsibility to anyone for the consequences of using it
10 or for whether it serves any particular purpose or works at all,
11 unless he says so in writing. Refer to the GNU CC General Public
12 License for full details.
13
14 Everyone is granted permission to copy, modify and redistribute
15 GNU CC, but only under the conditions described in the
16 GNU CC General Public License. A copy of this license is
17 supposed to have been given to you along with GNU CC so you
18 can know your rights and responsibilities. It should be in a
19 file named COPYING. Among other things, the copyright notice
20 and this notice must be preserved on all copies.
21 */
22
23 #ifdef __GNUG__
24 #pragma implementation
25 #endif
26 #include <Rational.h>
27 #include <std.h>
28 #include <math.h>
29 #include <values.h>
30
31
32
error(const char * msg) const33 volatile void Rational::error(const char* msg) const
34 {
35 (*lib_error_handler)("Rational", msg);
36 }
37
38 static const Integer _Int_One(1);
39
normalize()40 void Rational::normalize()
41 {
42 int s = sign(den);
43 if (s == 0)
44 error("Zero denominator.");
45 else if (s < 0)
46 {
47 den.negate();
48 num.negate();
49 }
50
51 Integer g = gcd(num, den);
52 if (ucompare(g, _Int_One) != 0)
53 {
54 num /= g;
55 den /= g;
56 }
57 }
58
add(const Rational & x,const Rational & y,Rational & r)59 void add(const Rational& x, const Rational& y, Rational& r)
60 {
61 if (&r != &x && &r != &y)
62 {
63 mul(x.num, y.den, r.num);
64 mul(x.den, y.num, r.den);
65 add(r.num, r.den, r.num);
66 mul(x.den, y.den, r.den);
67 }
68 else
69 {
70 Integer tmp;
71 mul(x.den, y.num, tmp);
72 mul(x.num, y.den, r.num);
73 add(r.num, tmp, r.num);
74 mul(x.den, y.den, r.den);
75 }
76 r.normalize();
77 }
78
sub(const Rational & x,const Rational & y,Rational & r)79 void sub(const Rational& x, const Rational& y, Rational& r)
80 {
81 if (&r != &x && &r != &y)
82 {
83 mul(x.num, y.den, r.num);
84 mul(x.den, y.num, r.den);
85 sub(r.num, r.den, r.num);
86 mul(x.den, y.den, r.den);
87 }
88 else
89 {
90 Integer tmp;
91 mul(x.den, y.num, tmp);
92 mul(x.num, y.den, r.num);
93 sub(r.num, tmp, r.num);
94 mul(x.den, y.den, r.den);
95 }
96 r.normalize();
97 }
98
mul(const Rational & x,const Rational & y,Rational & r)99 void mul(const Rational& x, const Rational& y, Rational& r)
100 {
101 mul(x.num, y.num, r.num);
102 mul(x.den, y.den, r.den);
103 r.normalize();
104 }
105
div(const Rational & x,const Rational & y,Rational & r)106 void div(const Rational& x, const Rational& y, Rational& r)
107 {
108 if (&r != &x && &r != &y)
109 {
110 mul(x.num, y.den, r.num);
111 mul(x.den, y.num, r.den);
112 }
113 else
114 {
115 Integer tmp;
116 mul(x.num, y.den, tmp);
117 mul(y.num, x.den, r.den);
118 r.num = tmp;
119 }
120 r.normalize();
121 }
122
123
124
125
invert()126 void Rational::invert()
127 {
128 Integer tmp = num;
129 num = den;
130 den = tmp;
131 int s = sign(den);
132 if (s == 0)
133 error("Zero denominator.");
134 else if (s < 0)
135 {
136 den.negate();
137 num.negate();
138 }
139 }
140
compare(const Rational & x,const Rational & y)141 int compare(const Rational& x, const Rational& y)
142 {
143 int xsgn = sign(x.num);
144 int ysgn = sign(y.num);
145 int d = xsgn - ysgn;
146 if (d == 0 && xsgn != 0) d = compare(x.num * y.den, x.den * y.num);
147 return d;
148 }
149
Rational(double x)150 Rational::Rational(double x)
151 {
152 num = 0;
153 den = 1;
154 if (x != 0.0)
155 {
156 int neg = x < 0;
157 if (neg)
158 x = -x;
159
160 const long shift = 15; // a safe shift per step
161 const double width = 32768.0; // = 2^shift
162 const int maxiter = 20; // ought not be necessary, but just in case,
163 // max 300 bits of precision
164 int expt;
165 double mantissa = frexp(x, &expt);
166 long exponent = expt;
167 double intpart;
168 int k = 0;
169 while (mantissa != 0.0 && k++ < maxiter)
170 {
171 mantissa *= width;
172 mantissa = modf(mantissa, &intpart);
173 num <<= shift;
174 num += (long)intpart;
175 exponent -= shift;
176 }
177 if (exponent > 0)
178 num <<= exponent;
179 else if (exponent < 0)
180 den <<= -exponent;
181 if (neg)
182 num.negate();
183 }
184 normalize();
185 }
186
187
trunc(const Rational & x)188 Integer trunc(const Rational& x)
189 {
190 return x.num / x.den ;
191 }
192
193
pow(const Rational & x,const Integer & y)194 Rational pow(const Rational& x, const Integer& y)
195 {
196 long yy = long(y);
197 return pow(x, yy);
198 }
199
200 #if defined(__GNUG__) && !defined(NO_NRV)
201
r(x)202 Rational operator - (const Rational& x) return r(x)
203 {
204 r.negate();
205 }
206
r(x)207 Rational abs(const Rational& x) return r(x)
208 {
209 if (sign(r.num) < 0) r.negate();
210 }
211
212
213 Rational sqr(const Rational& x) return r
214 {
215 mul(x.num, x.num, r.num);
216 mul(x.den, x.den, r.den);
217 r.normalize();
218 }
219
220 Integer floor(const Rational& x) return q
221 {
222 Integer r;
223 divide(x.num, x.den, q, r);
224 if (sign(x.num) < 0 && sign(r) != 0) q--;
225 }
226
227 Integer ceil(const Rational& x) return q
228 {
229 Integer r;
230 divide(x.num, x.den, q, r);
231 if (sign(x.num) >= 0 && sign(r) != 0) q++;
232 }
233
234 Integer round(const Rational& x) return q
235 {
236 Integer r;
237 divide(x.num, x.den, q, r);
238 r <<= 1;
239 if (ucompare(r, x.den) >= 0)
240 {
241 if (sign(x.num) >= 0)
242 q++;
243 else
244 q--;
245 }
246 }
247
248 // power: no need to normalize since num & den already relatively prime
249
250 Rational pow(const Rational& x, long y) return r
251 {
252 if (y >= 0)
253 {
254 pow(x.num, y, r.num);
255 pow(x.den, y, r.den);
256 }
257 else
258 {
259 y = -y;
260 pow(x.num, y, r.den);
261 pow(x.den, y, r.num);
262 if (sign(r.den) < 0)
263 {
264 r.num.negate();
265 r.den.negate();
266 }
267 }
268 }
269
270 #else
271
272 Rational operator - (const Rational& x)
273 {
274 Rational r(x); r.negate(); return r;
275 }
276
277 Rational abs(const Rational& x)
278 {
279 Rational r(x);
280 if (sign(r.num) < 0) r.negate();
281 return r;
282 }
283
284
285 Rational sqr(const Rational& x)
286 {
287 Rational r;
288 mul(x.num, x.num, r.num);
289 mul(x.den, x.den, r.den);
290 r.normalize();
291 return r;
292 }
293
294 Integer floor(const Rational& x)
295 {
296 Integer q;
297 Integer r;
298 divide(x.num, x.den, q, r);
299 if (sign(x.num) < 0 && sign(r) != 0) q--;
300 return q;
301 }
302
303 Integer ceil(const Rational& x)
304 {
305 Integer q;
306 Integer r;
307 divide(x.num, x.den, q, r);
308 if (sign(x.num) >= 0 && sign(r) != 0) q++;
309 return q;
310 }
311
312 Integer round(const Rational& x)
313 {
314 Integer q;
315 Integer r;
316 divide(x.num, x.den, q, r);
317 r <<= 1;
318 if (ucompare(r, x.den) >= 0)
319 {
320 if (sign(x.num) >= 0)
321 q++;
322 else
323 q--;
324 }
325 return q;
326 }
327
328 Rational pow(const Rational& x, long y)
329 {
330 Rational r;
331 if (y >= 0)
332 {
333 pow(x.num, y, r.num);
334 pow(x.den, y, r.den);
335 }
336 else
337 {
338 y = -y;
339 pow(x.num, y, r.den);
340 pow(x.den, y, r.num);
341 if (sign(r.den) < 0)
342 {
343 r.num.negate();
344 r.den.negate();
345 }
346 }
347 return r;
348 }
349
350 #endif
351
352 ostream& operator << (ostream& s, const Rational& y)
353 {
354 if (y.den == 1)
355 s << Itoa(y.num);
356 else
357 {
358 s << Itoa(y.num);
359 s << "/";
360 s << Itoa(y.den);
361 }
362 return s;
363 }
364
operator >>(istream & s,Rational & y)365 istream& operator >> (istream& s, Rational& y)
366 {
367 s >> y.num;
368 if (s)
369 {
370 char ch = 0;
371 s.get(ch);
372 if (ch == '/')
373 {
374 s >> y.den;
375 y.normalize();
376 }
377 else
378 {
379 s.unget(ch);
380 y.den = 1;
381 }
382 }
383 return s;
384 }
385
OK() const386 int Rational::OK() const
387 {
388 int v = num.OK() && den.OK(); // have valid num and denom
389 v &= sign(den) > 0; // denominator positive;
390 v &= ucompare(gcd(num, den), _Int_One) == 0; // relatively prime
391 if (!v) error("invariant failure");
392 return v;
393 }
394