![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to b1nuR.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint|
Return to b1nuR.c CVS log | Up to [CSRG BSD Unix] / 43BSD / contrib / B / src / bint |
1.1 root 1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
2:
3: /*
4: $Header: b1nuR.c,v 1.4 85/08/22 16:51:49 timo Exp $
5: */
6:
7: /* Rational arithmetic */
8:
9: #include "b.h"
10: #include "b0con.h"
11: #include "b1obj.h"
12: #include "b1num.h"
13: #include "b3err.h"
14:
15: /* Length calculations used for fraction sizes: */
16: #define Maxlen(u, v) \
17: (Roundsize(u) > Roundsize(v) ? Roundsize(u) : Roundsize(v))
18: #define Sumlen(u, v) (Roundsize(u)+Roundsize(v))
19: #define Difflen(u, v) (Roundsize(u)-Roundsize(v))
20:
21: /* To shut off lint and other warnings: */
22: #undef Copy
23: #define Copy(x) ((integer)copy((value)(x)))
24:
25: /* Globally used constants */
26:
27: rational rat_zero;
28: rational rat_half;
29:
30: /* Make a normalized rational from two integers */
31:
32: Visible rational mk_rat(x, y, len) integer x, y; int len; {
33: rational a;
34: integer u,v;
35:
36: if (y == int_0) {
37: if (interrupted)
38: return rat_zero;
39: syserr(MESS(1200, "mk_rat(x, y) with y=0"));
40: }
41:
42: if (x == int_0 && len <= 0) return (rational) Copy(rat_zero);
43:
44: if (Msd(y) < 0) { /* interchange signs */
45: u = int_neg(x);
46: v = int_neg(y);
47: } else {
48: u = Copy(x);
49: v = Copy(y);
50: }
51:
52: a = (rational) grab_rat();
53: if (len > 0 && len+2 <= Maxintlet) Length(a) = -2 - len;
54:
55: if (u == int_0 || v == int_1) {
56: /* No simplification possible */
57: Numerator(a) = Copy(u);
58: Denominator(a) = int_1;
59: } else {
60: integer g, abs_u;
61:
62: if (Msd(u) < 0) abs_u = int_neg(u);
63: else abs_u = Copy(u);
64: g = int_gcd(abs_u, v);
65: release((value) abs_u);
66:
67: if (g != int_1) {
68: Numerator(a) = int_quot(u, g);
69: Denominator(a) = int_quot(v, g);
70: } else {
71: Numerator(a) = Copy(u);
72: Denominator(a) = Copy(v);
73: }
74: release((value) g);
75: }
76:
77: release((value) u); release((value) v);
78:
79: return a;
80: }
81:
82:
83: /* Arithmetic on rational numbers */
84:
85: /* Shorthands: */
86: #define N(u) Numerator(u)
87: #define D(u) Denominator(u)
88:
89: Visible rational rat_sum(u, v) register rational u, v; {
90: integer t1, t2, t3, t4;
91: rational a;
92:
93: t2= int_prod(N(u), D(v));
94: t3= int_prod(N(v), D(u));
95: t1= int_sum(t2, t3);
96: t4= int_prod(D(u), D(v));
97: a= mk_rat(t1, t4, Maxlen(u, v));
98: release((value) t1); release((value) t2);
99: release((value) t3); release((value) t4);
100:
101: return a;
102: }
103:
104:
105: Visible rational rat_diff(u, v) register rational u, v; {
106: integer t1, t2, t3, t4;
107: rational a;
108:
109: t2= int_prod(N(u), D(v));
110: t3= int_prod(N(v), D(u));
111: t1= int_diff(t2, t3);
112: t4= int_prod(D(u), D(v));
113: a= mk_rat(t1, t4, Maxlen(u, v));
114: release((value) t1); release((value) t2);
115: release((value) t3); release((value) t4);
116:
117: return a;
118: }
119:
120:
121: Visible rational rat_prod(u, v) register rational u, v; {
122: integer t1, t2;
123: rational a;
124:
125: t1= int_prod(N(u), N(v));
126: t2= int_prod(D(u), D(v));
127: a= mk_rat(t1, t2, Sumlen(u, v));
128: release((value) t1); release((value) t2);
129:
130: return a;
131: }
132:
133:
134: Visible rational rat_quot(u, v) register rational u, v; {
135: integer t1, t2;
136: rational a;
137:
138: if (Numerator(v) == int_0) {
139: error(MESS(1201, "in u/v, v is zero"));
140: return (rational) Copy(rat_zero);
141: }
142:
143: t1= int_prod(N(u), D(v));
144: t2= int_prod(D(u), N(v));
145: a= mk_rat(t1, t2, Difflen(u, v));
146: release((value) t1); release((value) t2);
147:
148: return a;
149: }
150:
151:
152: Visible rational rat_neg(u) register rational u; {
153: register rational a;
154:
155: /* Avoid a real subtraction from zero */
156:
157: if (Numerator(u) == int_0) return (rational) Copy(u);
158:
159: a = (rational) grab_rat();
160: N(a) = int_neg(N(u));
161: D(a) = Copy(D(u));
162: Length(a) = Length(u);
163:
164: return a;
165: }
166:
167:
168: /* Rational number to the integral power */
169:
170: Visible rational rat_power(a, n) rational a; integer n; {
171: integer u, v, tu, tv, temp;
172:
173: if (n == int_0) return mk_rat(int_1, int_1, 0);
174:
175: if (Msd(n) < 0) {
176: if (Numerator(a) == int_0) {
177: error(MESS(1202, "in 0**n, n is negative"));
178: return (rational) Copy(a);
179: }
180: if (Msd(Numerator(a)) < 0) {
181: u= int_neg(Denominator(a));
182: v = int_neg(Numerator(a));
183: }
184: else {
185: u = Copy(Denominator(a));
186: v = Copy(Numerator(a));
187: }
188: n = int_neg(n);
189: } else {
190: if (Numerator(a) == int_0) return (rational) Copy(a);
191: /* To avoid necessary simplification later on */
192: u = Copy(Numerator(a));
193: v = Copy(Denominator(a));
194: n = Copy(n);
195: }
196:
197: tu = int_1;
198: tv = int_1;
199:
200: while (n != int_0 && !interrupted) {
201: if (Odd(Lsd(n))) {
202: if (u != int_1) {
203: temp = tu;
204: tu = int_prod(u, tu);
205: release((value) temp);
206: }
207: if (v != int_1) {
208: temp = tv;
209: tv = int_prod(v, tv);
210: release((value) temp);
211: }
212: if (n == int_1)
213: break; /* Avoid useless last squaring */
214: }
215:
216: /* Square u, v */
217:
218: if (u != int_1) {
219: temp = u;
220: u = int_prod(u, u);
221: release((value) temp);
222: }
223: if (v != int_1) {
224: temp = v;
225: v = int_prod(v, v);
226: release((value) temp);
227: }
228:
229: n = int_half(n);
230: } /* while (n!=0) */
231:
232: release((value) n);
233: release((value) u);
234: release((value) v);
235: a = (rational) grab_rat();
236: Numerator(a) = tu;
237: Denominator(a) = tv;
238:
239: return a;
240: }
241:
242:
243: /* Compare two rational numbers */
244:
245: Visible relation rat_comp(u, v) register rational u, v; {
246: int sd, su, sv;
247: integer nu, nv;
248:
249: /* 1. Compare pointers */
250: if (u == v || N(u) == N(v) && D(u) == D(v)) return 0;
251:
252: /* 2. Either zero? */
253: if (N(u) == int_0) return int_comp(int_0, N(v));
254: if (N(v) == int_0) return int_comp(N(u), int_0);
255:
256: /* 3. Compare signs */
257: su = Msd(N(u));
258: sv = Msd(N(v));
259: su = (su>0) - (su<0);
260: sv = (sv>0) - (sv<0);
261: if (su != sv) return su > sv ? 1 : -1;
262:
263: /* 4. Compute numerator of difference and return sign */
264: nu= int_prod(N(u), D(v));
265: nv= int_prod(N(v), D(u));
266: sd= int_comp(nu, nv);
267: release((value) nu); release((value) nv);
268: return sd;
269: }
270:
271: Visible Procedure rat_init() {
272: rat_zero = (rational) grab_rat();
273: Numerator(rat_zero) = int_0;
274: Denominator(rat_zero) = int_1;
275:
276: rat_half = (rational) grab_rat();
277: Numerator(rat_half) = int_1;
278: Denominator(rat_half) = int_2;
279: }
280:
281: Visible Procedure endrat() {
282: release((value) rat_zero);
283: release((value) rat_half);
284: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.