Commit | Line | Data |
---|---|---|
e1a1d2b4 C |
1 | /* |
2 | * Copyright (c) 1994 David I. Bell | |
3 | * Permission is granted to use, distribute, or modify this source, | |
4 | * provided that this copyright notice remains intact. | |
5 | * | |
6 | * Extended precision rational arithmetic primitive routines | |
7 | */ | |
8 | ||
9 | #include "qmath.h" | |
10 | ||
11 | ||
12 | NUMBER _qzero_ = { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1 }; | |
13 | NUMBER _qone_ = { { _oneval_, 1, 0 }, { _oneval_, 1, 0 }, 1 }; | |
14 | static NUMBER _qtwo_ = { { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1 }; | |
15 | static NUMBER _qten_ = { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1 }; | |
16 | NUMBER _qnegone_ = { { _oneval_, 1, 1 }, { _oneval_, 1, 0 }, 1 }; | |
17 | NUMBER _qonehalf_ = { { _oneval_, 1, 0 }, { _twoval_, 1, 0 }, 1 }; | |
18 | ||
19 | ||
20 | /* | |
21 | * Create another copy of a number. | |
22 | * q2 = qcopy(q1); | |
23 | */ | |
24 | NUMBER * | |
25 | qcopy(q) | |
26 | register NUMBER *q; | |
27 | { | |
28 | register NUMBER *r; | |
29 | ||
30 | r = qalloc(); | |
31 | r->num.sign = q->num.sign; | |
32 | if (!zisunit(q->num)) { | |
33 | r->num.len = q->num.len; | |
34 | r->num.v = alloc(r->num.len); | |
35 | zcopyval(q->num, r->num); | |
36 | } | |
37 | if (!zisunit(q->den)) { | |
38 | r->den.len = q->den.len; | |
39 | r->den.v = alloc(r->den.len); | |
40 | zcopyval(q->den, r->den); | |
41 | } | |
42 | return r; | |
43 | } | |
44 | ||
45 | ||
46 | /* | |
47 | * Convert a number to a normal integer. | |
48 | * i = qtoi(q); | |
49 | */ | |
50 | long | |
51 | qtoi(q) | |
52 | register NUMBER *q; | |
53 | { | |
54 | long i; | |
55 | ZVALUE res; | |
56 | ||
57 | if (qisint(q)) | |
58 | return ztoi(q->num); | |
59 | zquo(q->num, q->den, &res); | |
60 | i = ztoi(res); | |
61 | zfree(res); | |
62 | return i; | |
63 | } | |
64 | ||
65 | ||
66 | /* | |
67 | * Convert a normal integer into a number. | |
68 | * q = itoq(i); | |
69 | */ | |
70 | NUMBER * | |
71 | itoq(i) | |
72 | long i; | |
73 | { | |
74 | register NUMBER *q; | |
75 | ||
76 | if ((i >= -1) && (i <= 10)) { | |
77 | switch ((int) i) { | |
78 | case 0: q = &_qzero_; break; | |
79 | case 1: q = &_qone_; break; | |
80 | case 2: q = &_qtwo_; break; | |
81 | case 10: q = &_qten_; break; | |
82 | case -1: q = &_qnegone_; break; | |
83 | default: q = NULL; | |
84 | } | |
85 | if (q) | |
86 | return qlink(q); | |
87 | } | |
88 | q = qalloc(); | |
89 | itoz(i, &q->num); | |
90 | return q; | |
91 | } | |
92 | ||
93 | ||
94 | /* | |
95 | * Create a number from the given integral numerator and denominator. | |
96 | * q = iitoq(inum, iden); | |
97 | */ | |
98 | NUMBER * | |
99 | iitoq(inum, iden) | |
100 | long inum, iden; | |
101 | { | |
102 | register NUMBER *q; | |
103 | long d; | |
104 | BOOL sign; | |
105 | ||
106 | if (iden == 0) | |
107 | math_error("Division by zero"); | |
108 | if (inum == 0) | |
109 | return qlink(&_qzero_); | |
110 | sign = 0; | |
111 | if (inum < 0) { | |
112 | sign = 1; | |
113 | inum = -inum; | |
114 | } | |
115 | if (iden < 0) { | |
116 | sign = 1 - sign; | |
117 | iden = -iden; | |
118 | } | |
119 | d = iigcd(inum, iden); | |
120 | inum /= d; | |
121 | iden /= d; | |
122 | if (iden == 1) | |
123 | return itoq(sign ? -inum : inum); | |
124 | q = qalloc(); | |
125 | if (inum != 1) | |
126 | itoz(inum, &q->num); | |
127 | itoz(iden, &q->den); | |
128 | q->num.sign = sign; | |
129 | return q; | |
130 | } | |
131 | ||
132 | ||
133 | /* | |
134 | * Add two numbers to each other. | |
135 | * q3 = qadd(q1, q2); | |
136 | */ | |
137 | NUMBER * | |
138 | qadd(q1, q2) | |
139 | register NUMBER *q1, *q2; | |
140 | { | |
141 | NUMBER *r; | |
142 | ZVALUE t1, t2, temp, d1, d2, vpd1, upd1; | |
143 | ||
144 | if (qiszero(q1)) | |
145 | return qlink(q2); | |
146 | if (qiszero(q2)) | |
147 | return qlink(q1); | |
148 | r = qalloc(); | |
149 | /* | |
150 | * If either number is an integer, then the result is easy. | |
151 | */ | |
152 | if (qisint(q1) && qisint(q2)) { | |
153 | zadd(q1->num, q2->num, &r->num); | |
154 | return r; | |
155 | } | |
156 | if (qisint(q2)) { | |
157 | zmul(q1->den, q2->num, &temp); | |
158 | zadd(q1->num, temp, &r->num); | |
159 | zfree(temp); | |
160 | zcopy(q1->den, &r->den); | |
161 | return r; | |
162 | } | |
163 | if (qisint(q1)) { | |
164 | zmul(q2->den, q1->num, &temp); | |
165 | zadd(q2->num, temp, &r->num); | |
166 | zfree(temp); | |
167 | zcopy(q2->den, &r->den); | |
168 | return r; | |
169 | } | |
170 | /* | |
171 | * Both arguments are true fractions, so we need more work. | |
172 | * If the denominators are relatively prime, then the answer is the | |
173 | * straightforward cross product result with no need for reduction. | |
174 | */ | |
175 | zgcd(q1->den, q2->den, &d1); | |
176 | if (zisunit(d1)) { | |
177 | zfree(d1); | |
178 | zmul(q1->num, q2->den, &t1); | |
179 | zmul(q1->den, q2->num, &t2); | |
180 | zadd(t1, t2, &r->num); | |
181 | zfree(t1); | |
182 | zfree(t2); | |
183 | zmul(q1->den, q2->den, &r->den); | |
184 | return r; | |
185 | } | |
186 | /* | |
187 | * The calculation is now more complicated. | |
188 | * See Knuth Vol 2 for details. | |
189 | */ | |
190 | zquo(q2->den, d1, &vpd1); | |
191 | zquo(q1->den, d1, &upd1); | |
192 | zmul(q1->num, vpd1, &t1); | |
193 | zmul(q2->num, upd1, &t2); | |
194 | zadd(t1, t2, &temp); | |
195 | zfree(t1); | |
196 | zfree(t2); | |
197 | zfree(vpd1); | |
198 | zgcd(temp, d1, &d2); | |
199 | zfree(d1); | |
200 | if (zisunit(d2)) { | |
201 | zfree(d2); | |
202 | r->num = temp; | |
203 | zmul(upd1, q2->den, &r->den); | |
204 | zfree(upd1); | |
205 | return r; | |
206 | } | |
207 | zquo(temp, d2, &r->num); | |
208 | zfree(temp); | |
209 | zquo(q2->den, d2, &temp); | |
210 | zfree(d2); | |
211 | zmul(temp, upd1, &r->den); | |
212 | zfree(temp); | |
213 | zfree(upd1); | |
214 | return r; | |
215 | } | |
216 | ||
217 | ||
218 | /* | |
219 | * Subtract one number from another. | |
220 | * q3 = qsub(q1, q2); | |
221 | */ | |
222 | NUMBER * | |
223 | qsub(q1, q2) | |
224 | register NUMBER *q1, *q2; | |
225 | { | |
226 | NUMBER *r; | |
227 | ||
228 | if (q1 == q2) | |
229 | return qlink(&_qzero_); | |
230 | if (qiszero(q2)) | |
231 | return qlink(q1); | |
232 | if (qisint(q1) && qisint(q2)) { | |
233 | r = qalloc(); | |
234 | zsub(q1->num, q2->num, &r->num); | |
235 | return r; | |
236 | } | |
237 | q2 = qneg(q2); | |
238 | if (qiszero(q1)) | |
239 | return q2; | |
240 | r = qadd(q1, q2); | |
241 | qfree(q2); | |
242 | return r; | |
243 | } | |
244 | ||
245 | ||
246 | /* | |
247 | * Increment a number by one. | |
248 | */ | |
249 | NUMBER * | |
250 | qinc(q) | |
251 | NUMBER *q; | |
252 | { | |
253 | NUMBER *r; | |
254 | ||
255 | r = qalloc(); | |
256 | if (qisint(q)) { | |
257 | zadd(q->num, _one_, &r->num); | |
258 | return r; | |
259 | } | |
260 | zadd(q->num, q->den, &r->num); | |
261 | zcopy(q->den, &r->den); | |
262 | return r; | |
263 | } | |
264 | ||
265 | ||
266 | /* | |
267 | * Decrement a number by one. | |
268 | */ | |
269 | NUMBER * | |
270 | qdec(q) | |
271 | NUMBER *q; | |
272 | { | |
273 | NUMBER *r; | |
274 | ||
275 | r = qalloc(); | |
276 | if (qisint(q)) { | |
277 | zsub(q->num, _one_, &r->num); | |
278 | return r; | |
279 | } | |
280 | zsub(q->num, q->den, &r->num); | |
281 | zcopy(q->den, &r->den); | |
282 | return r; | |
283 | } | |
284 | ||
285 | ||
286 | /* | |
287 | * Add a normal small integer value to an arbitrary number. | |
288 | */ | |
289 | NUMBER * | |
290 | qaddi(q1, n) | |
291 | NUMBER *q1; | |
292 | long n; | |
293 | { | |
294 | NUMBER addnum; /* temporary number */ | |
295 | HALF addval[2]; /* value of small number */ | |
296 | BOOL neg; /* TRUE if number is neg */ | |
297 | ||
298 | if (n == 0) | |
299 | return qlink(q1); | |
300 | if (n == 1) | |
301 | return qinc(q1); | |
302 | if (n == -1) | |
303 | return qdec(q1); | |
304 | if (qiszero(q1)) | |
305 | return itoq(n); | |
306 | addnum.num.sign = 0; | |
307 | addnum.num.len = 1; | |
308 | addnum.num.v = addval; | |
309 | addnum.den = _one_; | |
310 | neg = (n < 0); | |
311 | if (neg) | |
312 | n = -n; | |
313 | addval[0] = (HALF) n; | |
314 | n = (((FULL) n) >> BASEB); | |
315 | if (n) { | |
316 | addval[1] = (HALF) n; | |
317 | addnum.num.len = 2; | |
318 | } | |
319 | if (neg) | |
320 | return qsub(q1, &addnum); | |
321 | else | |
322 | return qadd(q1, &addnum); | |
323 | } | |
324 | ||
325 | ||
326 | /* | |
327 | * Multiply two numbers. | |
328 | * q3 = qmul(q1, q2); | |
329 | */ | |
330 | NUMBER * | |
331 | qmul(q1, q2) | |
332 | register NUMBER *q1, *q2; | |
333 | { | |
334 | NUMBER *r; /* returned value */ | |
335 | ZVALUE n1, n2, d1, d2; /* numerators and denominators */ | |
336 | ZVALUE tmp; | |
337 | ||
338 | if (qiszero(q1) || qiszero(q2)) | |
339 | return qlink(&_qzero_); | |
340 | if (qisone(q1)) | |
341 | return qlink(q2); | |
342 | if (qisone(q2)) | |
343 | return qlink(q1); | |
344 | if (qisint(q1) && qisint(q2)) { /* easy results if integers */ | |
345 | r = qalloc(); | |
346 | zmul(q1->num, q2->num, &r->num); | |
347 | return r; | |
348 | } | |
349 | n1 = q1->num; | |
350 | n2 = q2->num; | |
351 | d1 = q1->den; | |
352 | d2 = q2->den; | |
353 | if (ziszero(d1) || ziszero(d2)) | |
354 | math_error("Division by zero"); | |
355 | if (ziszero(n1) || ziszero(n2)) | |
356 | return qlink(&_qzero_); | |
357 | if (!zisunit(n1) && !zisunit(d2)) { /* possibly reduce */ | |
358 | zgcd(n1, d2, &tmp); | |
359 | if (!zisunit(tmp)) { | |
360 | zquo(q1->num, tmp, &n1); | |
361 | zquo(q2->den, tmp, &d2); | |
362 | } | |
363 | zfree(tmp); | |
364 | } | |
365 | if (!zisunit(n2) && !zisunit(d1)) { /* again possibly reduce */ | |
366 | zgcd(n2, d1, &tmp); | |
367 | if (!zisunit(tmp)) { | |
368 | zquo(q2->num, tmp, &n2); | |
369 | zquo(q1->den, tmp, &d1); | |
370 | } | |
371 | zfree(tmp); | |
372 | } | |
373 | r = qalloc(); | |
374 | zmul(n1, n2, &r->num); | |
375 | zmul(d1, d2, &r->den); | |
376 | if (q1->num.v != n1.v) | |
377 | zfree(n1); | |
378 | if (q1->den.v != d1.v) | |
379 | zfree(d1); | |
380 | if (q2->num.v != n2.v) | |
381 | zfree(n2); | |
382 | if (q2->den.v != d2.v) | |
383 | zfree(d2); | |
384 | return r; | |
385 | } | |
386 | ||
387 | ||
388 | /* | |
389 | * Multiply a number by a small integer. | |
390 | * q2 = qmuli(q1, n); | |
391 | */ | |
392 | NUMBER * | |
393 | qmuli(q, n) | |
394 | NUMBER *q; | |
395 | long n; | |
396 | { | |
397 | NUMBER *r; | |
398 | long d; /* gcd of multiplier and denominator */ | |
399 | int sign; | |
400 | ||
401 | if ((n == 0) || qiszero(q)) | |
402 | return qlink(&_qzero_); | |
403 | if (n == 1) | |
404 | return qlink(q); | |
405 | r = qalloc(); | |
406 | if (qisint(q)) { | |
407 | zmuli(q->num, n, &r->num); | |
408 | return r; | |
409 | } | |
410 | sign = 1; | |
411 | if (n < 0) { | |
412 | n = -n; | |
413 | sign = -1; | |
414 | } | |
415 | d = zmodi(q->den, n); | |
416 | d = iigcd(d, n); | |
417 | zmuli(q->num, (n * sign) / d, &r->num); | |
418 | (void) zdivi(q->den, d, &r->den); | |
419 | return r; | |
420 | } | |
421 | ||
422 | ||
423 | /* | |
424 | * Divide two numbers (as fractions). | |
425 | * q3 = qdiv(q1, q2); | |
426 | */ | |
427 | NUMBER * | |
428 | qdiv(q1, q2) | |
429 | register NUMBER *q1, *q2; | |
430 | { | |
431 | NUMBER temp; | |
432 | ||
433 | if (qiszero(q2)) | |
434 | math_error("Division by zero"); | |
435 | if ((q1 == q2) || !qcmp(q1, q2)) | |
436 | return qlink(&_qone_); | |
437 | if (qisone(q1)) | |
438 | return qinv(q2); | |
439 | temp.num = q2->den; | |
440 | temp.den = q2->num; | |
441 | temp.num.sign = temp.den.sign; | |
442 | temp.den.sign = 0; | |
443 | temp.links = 1; | |
444 | return qmul(q1, &temp); | |
445 | } | |
446 | ||
447 | ||
448 | /* | |
449 | * Divide a number by a small integer. | |
450 | * q2 = qdivi(q1, n); | |
451 | */ | |
452 | NUMBER * | |
453 | qdivi(q, n) | |
454 | NUMBER *q; | |
455 | long n; | |
456 | { | |
457 | NUMBER *r; | |
458 | long d; /* gcd of divisor and numerator */ | |
459 | int sign; | |
460 | ||
461 | if (n == 0) | |
462 | math_error("Division by zero"); | |
463 | if ((n == 1) || qiszero(q)) | |
464 | return qlink(q); | |
465 | sign = 1; | |
466 | if (n < 0) { | |
467 | n = -n; | |
468 | sign = -1; | |
469 | } | |
470 | r = qalloc(); | |
471 | d = zmodi(q->num, n); | |
472 | d = iigcd(d, n); | |
473 | (void) zdivi(q->num, d * sign, &r->num); | |
474 | zmuli(q->den, n / d, &r->den); | |
475 | return r; | |
476 | } | |
477 | ||
478 | ||
479 | /* | |
480 | * Return the quotient when one number is divided by another. | |
481 | * This works for fractional values also, and in all cases: | |
482 | * qquo(q1, q2) = int(q1 / q2). | |
483 | */ | |
484 | NUMBER * | |
485 | qquo(q1, q2) | |
486 | register NUMBER *q1, *q2; | |
487 | { | |
488 | ZVALUE num, den, res; | |
489 | NUMBER *q; | |
490 | ||
491 | if (zisunit(q1->num)) | |
492 | num = q2->den; | |
493 | else if (zisunit(q2->den)) | |
494 | num = q1->num; | |
495 | else | |
496 | zmul(q1->num, q2->den, &num); | |
497 | if (zisunit(q1->den)) | |
498 | den = q2->num; | |
499 | else if (zisunit(q2->num)) | |
500 | den = q1->den; | |
501 | else | |
502 | zmul(q1->den, q2->num, &den); | |
503 | zquo(num, den, &res); | |
504 | if ((num.v != q2->den.v) && (num.v != q1->num.v)) | |
505 | zfree(num); | |
506 | if ((den.v != q2->num.v) && (den.v != q1->den.v)) | |
507 | zfree(den); | |
508 | if (ziszero(res)) { | |
509 | zfree(res); | |
510 | return qlink(&_qzero_); | |
511 | } | |
512 | res.sign = (q1->num.sign != q2->num.sign); | |
513 | if (zisunit(res)) { | |
514 | q = (res.sign ? &_qnegone_ : &_qone_); | |
515 | zfree(res); | |
516 | return qlink(q); | |
517 | } | |
518 | q = qalloc(); | |
519 | q->num = res; | |
520 | return q; | |
521 | } | |
522 | ||
523 | ||
524 | /* | |
525 | * Return the absolute value of a number. | |
526 | * q2 = qabs(q1); | |
527 | */ | |
528 | NUMBER * | |
529 | qabs(q) | |
530 | register NUMBER *q; | |
531 | { | |
532 | register NUMBER *r; | |
533 | ||
534 | if (q->num.sign == 0) | |
535 | return qlink(q); | |
536 | r = qalloc(); | |
537 | if (!zisunit(q->num)) | |
538 | zcopy(q->num, &r->num); | |
539 | if (!zisunit(q->den)) | |
540 | zcopy(q->den, &r->den); | |
541 | r->num.sign = 0; | |
542 | return r; | |
543 | } | |
544 | ||
545 | ||
546 | /* | |
547 | * Negate a number. | |
548 | * q2 = qneg(q1); | |
549 | */ | |
550 | NUMBER * | |
551 | qneg(q) | |
552 | register NUMBER *q; | |
553 | { | |
554 | register NUMBER *r; | |
555 | ||
556 | if (qiszero(q)) | |
557 | return qlink(&_qzero_); | |
558 | r = qalloc(); | |
559 | if (!zisunit(q->num)) | |
560 | zcopy(q->num, &r->num); | |
561 | if (!zisunit(q->den)) | |
562 | zcopy(q->den, &r->den); | |
563 | r->num.sign = !q->num.sign; | |
564 | return r; | |
565 | } | |
566 | ||
567 | ||
568 | /* | |
569 | * Return the sign of a number (-1, 0 or 1) | |
570 | */ | |
571 | NUMBER * | |
572 | qsign(q) | |
573 | NUMBER *q; | |
574 | { | |
575 | if (qiszero(q)) | |
576 | return qlink(&_qzero_); | |
577 | if (qisneg(q)) | |
578 | return qlink(&_qnegone_); | |
579 | return qlink(&_qone_); | |
580 | } | |
581 | ||
582 | ||
583 | /* | |
584 | * Invert a number. | |
585 | * q2 = qinv(q1); | |
586 | */ | |
587 | NUMBER * | |
588 | qinv(q) | |
589 | register NUMBER *q; | |
590 | { | |
591 | register NUMBER *r; | |
592 | ||
593 | if (qisunit(q)) { | |
594 | r = (qisneg(q) ? &_qnegone_ : &_qone_); | |
595 | return qlink(r); | |
596 | } | |
597 | if (qiszero(q)) | |
598 | math_error("Division by zero"); | |
599 | r = qalloc(); | |
600 | if (!zisunit(q->num)) | |
601 | zcopy(q->num, &r->den); | |
602 | if (!zisunit(q->den)) | |
603 | zcopy(q->den, &r->num); | |
604 | r->num.sign = q->num.sign; | |
605 | r->den.sign = 0; | |
606 | return r; | |
607 | } | |
608 | ||
609 | ||
610 | /* | |
611 | * Return just the numerator of a number. | |
612 | * q2 = qnum(q1); | |
613 | */ | |
614 | NUMBER * | |
615 | qnum(q) | |
616 | register NUMBER *q; | |
617 | { | |
618 | register NUMBER *r; | |
619 | ||
620 | if (qisint(q)) | |
621 | return qlink(q); | |
622 | if (zisunit(q->num)) { | |
623 | r = (qisneg(q) ? &_qnegone_ : &_qone_); | |
624 | return qlink(r); | |
625 | } | |
626 | r = qalloc(); | |
627 | zcopy(q->num, &r->num); | |
628 | return r; | |
629 | } | |
630 | ||
631 | ||
632 | /* | |
633 | * Return just the denominator of a number. | |
634 | * q2 = qden(q1); | |
635 | */ | |
636 | NUMBER * | |
637 | qden(q) | |
638 | register NUMBER *q; | |
639 | { | |
640 | register NUMBER *r; | |
641 | ||
642 | if (qisint(q)) | |
643 | return qlink(&_qone_); | |
644 | r = qalloc(); | |
645 | zcopy(q->den, &r->num); | |
646 | return r; | |
647 | } | |
648 | ||
649 | ||
650 | /* | |
651 | * Return the fractional part of a number. | |
652 | * q2 = qfrac(q1); | |
653 | */ | |
654 | NUMBER * | |
655 | qfrac(q) | |
656 | register NUMBER *q; | |
657 | { | |
658 | register NUMBER *r; | |
659 | ZVALUE z; | |
660 | ||
661 | if (qisint(q)) | |
662 | return qlink(&_qzero_); | |
663 | if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) && | |
664 | (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1]))) | |
665 | return qlink(q); | |
666 | r = qalloc(); | |
667 | if (qisneg(q)) { | |
668 | zmod(q->num, q->den, &z); | |
669 | zsub(q->den, z, &r->num); | |
670 | zfree(z); | |
671 | } else { | |
672 | zmod(q->num, q->den, &r->num); | |
673 | } | |
674 | zcopy(q->den, &r->den); | |
675 | r->num.sign = q->num.sign; | |
676 | return r; | |
677 | } | |
678 | ||
679 | ||
680 | /* | |
681 | * Return the integral part of a number. | |
682 | * q2 = qint(q1); | |
683 | */ | |
684 | NUMBER * | |
685 | qint(q) | |
686 | register NUMBER *q; | |
687 | { | |
688 | register NUMBER *r; | |
689 | ||
690 | if (qisint(q)) | |
691 | return qlink(q); | |
692 | if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) && | |
693 | (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1]))) | |
694 | return qlink(&_qzero_); | |
695 | r = qalloc(); | |
696 | zquo(q->num, q->den, &r->num); | |
697 | return r; | |
698 | } | |
699 | ||
700 | ||
701 | /* | |
702 | * Compute the square of a number. | |
703 | */ | |
704 | NUMBER * | |
705 | qsquare(q) | |
706 | register NUMBER *q; | |
707 | { | |
708 | ZVALUE num, den; | |
709 | ||
710 | if (qiszero(q)) | |
711 | return qlink(&_qzero_); | |
712 | if (qisunit(q)) | |
713 | return qlink(&_qone_); | |
714 | num = q->num; | |
715 | den = q->den; | |
716 | q = qalloc(); | |
717 | if (!zisunit(num)) | |
718 | zsquare(num, &q->num); | |
719 | if (!zisunit(den)) | |
720 | zsquare(den, &q->den); | |
721 | return q; | |
722 | } | |
723 | ||
724 | ||
725 | /* | |
726 | * Shift an integer by a given number of bits. This multiplies the number | |
727 | * by the appropriate power of two. Positive numbers shift left, negative | |
728 | * ones shift right. Low bits are truncated when shifting right. | |
729 | */ | |
730 | NUMBER * | |
731 | qshift(q, n) | |
732 | NUMBER *q; | |
733 | long n; | |
734 | { | |
735 | register NUMBER *r; | |
736 | ||
737 | if (qisfrac(q)) | |
738 | math_error("Shift of non-integer"); | |
739 | if (qiszero(q) || (n == 0)) | |
740 | return qlink(q); | |
741 | if (n <= -(q->num.len * BASEB)) | |
742 | return qlink(&_qzero_); | |
743 | r = qalloc(); | |
744 | zshift(q->num, n, &r->num); | |
745 | return r; | |
746 | } | |
747 | ||
748 | ||
749 | /* | |
750 | * Scale a number by a power of two, as in: | |
751 | * ans = q * 2^n. | |
752 | * This is similar to shifting, except that fractions work. | |
753 | */ | |
754 | NUMBER * | |
755 | qscale(q, pow) | |
756 | NUMBER *q; | |
757 | long pow; | |
758 | { | |
759 | long numshift, denshift, tmp; | |
760 | NUMBER *r; | |
761 | ||
762 | if (qiszero(q) || (pow == 0)) | |
763 | return qlink(q); | |
764 | if ((pow > 1000000L) || (pow < -1000000L)) | |
765 | math_error("Very large scale value"); | |
766 | numshift = zisodd(q->num) ? 0 : zlowbit(q->num); | |
767 | denshift = zisodd(q->den) ? 0 : zlowbit(q->den); | |
768 | if (pow > 0) { | |
769 | tmp = pow; | |
770 | if (tmp > denshift) | |
771 | tmp = denshift; | |
772 | denshift = -tmp; | |
773 | numshift = (pow - tmp); | |
774 | } else { | |
775 | pow = -pow; | |
776 | tmp = pow; | |
777 | if (tmp > numshift) | |
778 | tmp = numshift; | |
779 | numshift = -tmp; | |
780 | denshift = (pow - tmp); | |
781 | } | |
782 | r = qalloc(); | |
783 | if (numshift) | |
784 | zshift(q->num, numshift, &r->num); | |
785 | else | |
786 | zcopy(q->num, &r->num); | |
787 | if (denshift) | |
788 | zshift(q->den, denshift, &r->den); | |
789 | else | |
790 | zcopy(q->den, &r->den); | |
791 | return r; | |
792 | } | |
793 | ||
794 | ||
795 | /* | |
796 | * Return the minimum of two numbers. | |
797 | */ | |
798 | NUMBER * | |
799 | qmin(q1, q2) | |
800 | NUMBER *q1, *q2; | |
801 | { | |
802 | if (q1 == q2) | |
803 | return qlink(q1); | |
804 | if (qrel(q1, q2) > 0) | |
805 | q1 = q2; | |
806 | return qlink(q1); | |
807 | } | |
808 | ||
809 | ||
810 | /* | |
811 | * Return the maximum of two numbers. | |
812 | */ | |
813 | NUMBER * | |
814 | qmax(q1, q2) | |
815 | NUMBER *q1, *q2; | |
816 | { | |
817 | if (q1 == q2) | |
818 | return qlink(q1); | |
819 | if (qrel(q1, q2) < 0) | |
820 | q1 = q2; | |
821 | return qlink(q1); | |
822 | } | |
823 | ||
824 | ||
825 | /* | |
826 | * Perform the logical OR of two integers. | |
827 | */ | |
828 | NUMBER * | |
829 | qor(q1, q2) | |
830 | NUMBER *q1, *q2; | |
831 | { | |
832 | register NUMBER *r; | |
833 | ||
834 | if (qisfrac(q1) || qisfrac(q2)) | |
835 | math_error("Non-integers for logical or"); | |
836 | if ((q1 == q2) || qiszero(q2)) | |
837 | return qlink(q1); | |
838 | if (qiszero(q1)) | |
839 | return qlink(q2); | |
840 | r = qalloc(); | |
841 | zor(q1->num, q2->num, &r->num); | |
842 | return r; | |
843 | } | |
844 | ||
845 | ||
846 | /* | |
847 | * Perform the logical AND of two integers. | |
848 | */ | |
849 | NUMBER * | |
850 | qand(q1, q2) | |
851 | NUMBER *q1, *q2; | |
852 | { | |
853 | register NUMBER *r; | |
854 | ZVALUE res; | |
855 | ||
856 | if (qisfrac(q1) || qisfrac(q2)) | |
857 | math_error("Non-integers for logical and"); | |
858 | if (q1 == q2) | |
859 | return qlink(q1); | |
860 | if (qiszero(q1) || qiszero(q2)) | |
861 | return qlink(&_qzero_); | |
862 | zand(q1->num, q2->num, &res); | |
863 | if (ziszero(res)) { | |
864 | zfree(res); | |
865 | return qlink(&_qzero_); | |
866 | } | |
867 | r = qalloc(); | |
868 | r->num = res; | |
869 | return r; | |
870 | } | |
871 | ||
872 | ||
873 | /* | |
874 | * Perform the logical XOR of two integers. | |
875 | */ | |
876 | NUMBER * | |
877 | qxor(q1, q2) | |
878 | NUMBER *q1, *q2; | |
879 | { | |
880 | register NUMBER *r; | |
881 | ZVALUE res; | |
882 | ||
883 | if (qisfrac(q1) || qisfrac(q2)) | |
884 | math_error("Non-integers for logical xor"); | |
885 | if (q1 == q2) | |
886 | return qlink(&_qzero_); | |
887 | if (qiszero(q1)) | |
888 | return qlink(q2); | |
889 | if (qiszero(q2)) | |
890 | return qlink(q1); | |
891 | zxor(q1->num, q2->num, &res); | |
892 | if (ziszero(res)) { | |
893 | zfree(res); | |
894 | return qlink(&_qzero_); | |
895 | } | |
896 | r = qalloc(); | |
897 | r->num = res; | |
898 | return r; | |
899 | } | |
900 | ||
901 | ||
902 | #if 0 | |
903 | /* | |
904 | * Return the number whose binary representation only has the specified | |
905 | * bit set (counting from zero). This thus produces a given power of two. | |
906 | */ | |
907 | NUMBER * | |
908 | qbitvalue(n) | |
909 | long n; | |
910 | { | |
911 | register NUMBER *r; | |
912 | ||
913 | if (n <= 0) | |
914 | return qlink(&_qone_); | |
915 | r = qalloc(); | |
916 | zbitvalue(n, &r->num); | |
917 | return r; | |
918 | } | |
919 | ||
920 | ||
921 | /* | |
922 | * Test to see if the specified bit of a number is on (counted from zero). | |
923 | * Returns TRUE if the bit is set, or FALSE if it is not. | |
924 | * i = qbittest(q, n); | |
925 | */ | |
926 | BOOL | |
927 | qbittest(q, n) | |
928 | register NUMBER *q; | |
929 | long n; | |
930 | { | |
931 | int x, y; | |
932 | ||
933 | if ((n < 0) || (n >= (q->num.len * BASEB))) | |
934 | return FALSE; | |
935 | x = q->num.v[n / BASEB]; | |
936 | y = (1 << (n % BASEB)); | |
937 | return ((x & y) != 0); | |
938 | } | |
939 | #endif | |
940 | ||
941 | ||
942 | /* | |
943 | * Return the precision of a number (usually for examining an epsilon value). | |
944 | * This is the largest power of two whose reciprocal is not smaller in absolute | |
945 | * value than the specified number. For example, qbitprec(1/100) = 6. | |
946 | * Numbers larger than one have a precision of zero. | |
947 | */ | |
948 | long | |
949 | qprecision(q) | |
950 | NUMBER *q; | |
951 | { | |
952 | long r; | |
953 | ||
954 | if (qisint(q)) | |
955 | return 0; | |
956 | if (zisunit(q->num)) | |
957 | return zhighbit(q->den); | |
958 | r = zhighbit(q->den) - zhighbit(q->num) - 1; | |
959 | if (r < 0) | |
960 | r = 0; | |
961 | return r; | |
962 | } | |
963 | ||
964 | ||
965 | #if 0 | |
966 | /* | |
967 | * Return an integer indicating the sign of a number (-1, 0, or 1). | |
968 | * i = qtst(q); | |
969 | */ | |
970 | FLAG | |
971 | qtest(q) | |
972 | register NUMBER *q; | |
973 | { | |
974 | if (!ztest(q->num)) | |
975 | return 0; | |
976 | if (q->num.sign) | |
977 | return -1; | |
978 | return 1; | |
979 | } | |
980 | #endif | |
981 | ||
982 | ||
983 | /* | |
984 | * Determine whether or not one number exactly divides another one. | |
985 | * Returns TRUE if the first number is an integer multiple of the second one. | |
986 | */ | |
987 | BOOL | |
988 | qdivides(q1, q2) | |
989 | NUMBER *q1, *q2; | |
990 | { | |
991 | if (qiszero(q1)) | |
992 | return TRUE; | |
993 | if (qisint(q1) && qisint(q2)) { | |
994 | if (qisunit(q2)) | |
995 | return TRUE; | |
996 | return zdivides(q1->num, q2->num); | |
997 | } | |
998 | return zdivides(q1->num, q2->num) && zdivides(q2->den, q1->den); | |
999 | } | |
1000 | ||
1001 | ||
1002 | /* | |
1003 | * Compare two numbers and return an integer indicating their relative size. | |
1004 | * i = qrel(q1, q2); | |
1005 | */ | |
1006 | FLAG | |
1007 | qrel(q1, q2) | |
1008 | register NUMBER *q1, *q2; | |
1009 | { | |
1010 | ZVALUE z1, z2; | |
1011 | long wc1, wc2; | |
1012 | int sign; | |
1013 | int z1f = 0, z2f = 0; | |
1014 | ||
1015 | if (q1 == q2) | |
1016 | return 0; | |
1017 | sign = q2->num.sign - q1->num.sign; | |
1018 | if (sign) | |
1019 | return sign; | |
1020 | if (qiszero(q2)) | |
1021 | return !qiszero(q1); | |
1022 | if (qiszero(q1)) | |
1023 | return -1; | |
1024 | /* | |
1025 | * Make a quick comparison by calculating the number of words resulting as | |
1026 | * if we multiplied through by the denominators, and then comparing the | |
1027 | * word counts. | |
1028 | */ | |
1029 | sign = 1; | |
1030 | if (qisneg(q1)) | |
1031 | sign = -1; | |
1032 | wc1 = q1->num.len + q2->den.len; | |
1033 | wc2 = q2->num.len + q1->den.len; | |
1034 | if (wc1 < wc2 - 1) | |
1035 | return -sign; | |
1036 | if (wc2 < wc1 - 1) | |
1037 | return sign; | |
1038 | /* | |
1039 | * Quick check failed, must actually do the full comparison. | |
1040 | */ | |
1041 | if (zisunit(q2->den)) | |
1042 | z1 = q1->num; | |
1043 | else if (zisone(q1->num)) | |
1044 | z1 = q2->den; | |
1045 | else { | |
1046 | z1f = 1; | |
1047 | zmul(q1->num, q2->den, &z1); | |
1048 | } | |
1049 | if (zisunit(q1->den)) | |
1050 | z2 = q2->num; | |
1051 | else if (zisone(q2->num)) | |
1052 | z2 = q1->den; | |
1053 | else { | |
1054 | z2f = 1; | |
1055 | zmul(q2->num, q1->den, &z2); | |
1056 | } | |
1057 | sign = zrel(z1, z2); | |
1058 | if (z1f) | |
1059 | zfree(z1); | |
1060 | if (z2f) | |
1061 | zfree(z2); | |
1062 | return sign; | |
1063 | } | |
1064 | ||
1065 | ||
1066 | /* | |
1067 | * Compare two numbers to see if they are equal. | |
1068 | * This differs from qrel in that the numbers are not ordered. | |
1069 | * Returns TRUE if they differ. | |
1070 | */ | |
1071 | BOOL | |
1072 | qcmp(q1, q2) | |
1073 | register NUMBER *q1, *q2; | |
1074 | { | |
1075 | if (q1 == q2) | |
1076 | return FALSE; | |
1077 | if ((q1->num.sign != q2->num.sign) || (q1->num.len != q2->num.len) || | |
1078 | (q2->den.len != q2->den.len) || (*q1->num.v != *q2->num.v) || | |
1079 | (*q1->den.v != *q2->den.v)) | |
1080 | return TRUE; | |
1081 | if (zcmp(q1->num, q2->num)) | |
1082 | return TRUE; | |
1083 | if (qisint(q1)) | |
1084 | return FALSE; | |
1085 | return zcmp(q1->den, q2->den); | |
1086 | } | |
1087 | ||
1088 | ||
1089 | /* | |
1090 | * Compare a number against a normal small integer. | |
1091 | * Returns 1, 0, or -1, according to whether the first number is greater, | |
1092 | * equal, or less than the second number. | |
1093 | * n = qreli(q, n); | |
1094 | */ | |
1095 | FLAG | |
1096 | qreli(q, n) | |
1097 | NUMBER *q; | |
1098 | long n; | |
1099 | { | |
1100 | int sign; | |
1101 | ZVALUE num; | |
1102 | HALF h2[2]; | |
1103 | NUMBER q2; | |
1104 | ||
1105 | sign = ztest(q->num); /* do trivial sign checks */ | |
1106 | if (sign == 0) { | |
1107 | if (n > 0) | |
1108 | return -1; | |
1109 | return (n < 0); | |
1110 | } | |
1111 | if ((sign < 0) && (n >= 0)) | |
1112 | return -1; | |
1113 | if ((sign > 0) && (n <= 0)) | |
1114 | return 1; | |
1115 | n *= sign; | |
1116 | if (n == 1) { /* quick check against 1 or -1 */ | |
1117 | num = q->num; | |
1118 | num.sign = 0; | |
1119 | return (sign * zrel(num, q->den)); | |
1120 | } | |
1121 | num.sign = (sign < 0); | |
1122 | num.len = 1 + (n >= BASE); | |
1123 | num.v = h2; | |
1124 | h2[0] = (n & BASE1); | |
1125 | h2[1] = (n >> BASEB); | |
1126 | if (zisunit(q->den)) /* integer compare if no denominator */ | |
1127 | return zrel(q->num, num); | |
1128 | q2.num = num; | |
1129 | q2.den = _one_; | |
1130 | q2.links = 1; | |
1131 | return qrel(q, &q2); /* full fractional compare */ | |
1132 | } | |
1133 | ||
1134 | ||
1135 | /* | |
1136 | * Compare a number against a small integer to see if they are equal. | |
1137 | * Returns TRUE if they differ. | |
1138 | */ | |
1139 | BOOL | |
1140 | qcmpi(q, n) | |
1141 | NUMBER *q; | |
1142 | long n; | |
1143 | { | |
1144 | long len; | |
1145 | ||
1146 | len = q->num.len; | |
1147 | if ((len > 2) || qisfrac(q) || (q->num.sign != (n < 0))) | |
1148 | return TRUE; | |
1149 | if (n < 0) | |
1150 | n = -n; | |
1151 | if (((HALF)(n)) != q->num.v[0]) | |
1152 | return TRUE; | |
1153 | n = ((FULL) n) >> BASEB; | |
1154 | return (((n != 0) != (len == 2)) || (n != q->num.v[1])); | |
1155 | } | |
1156 | ||
1157 | ||
1158 | /* | |
1159 | * Number node allocation routines | |
1160 | */ | |
1161 | ||
1162 | #define NNALLOC 1000 | |
1163 | ||
1164 | union allocNode { | |
1165 | NUMBER num; | |
1166 | union allocNode *link; | |
1167 | }; | |
1168 | ||
1169 | static union allocNode *freeNum; | |
1170 | ||
1171 | ||
1172 | NUMBER * | |
1173 | qalloc() | |
1174 | { | |
1175 | register union allocNode *temp; | |
1176 | ||
1177 | if (freeNum == NULL) { | |
1178 | freeNum = (union allocNode *) | |
1179 | malloc(sizeof (NUMBER) * NNALLOC); | |
1180 | if (freeNum == NULL) | |
1181 | math_error("Not enough memory"); | |
1182 | freeNum[NNALLOC-1].link = NULL; | |
1183 | for (temp=freeNum+NNALLOC-2; temp >= freeNum; --temp) { | |
1184 | temp->link = temp+1; | |
1185 | } | |
1186 | } | |
1187 | temp = freeNum; | |
1188 | freeNum = temp->link; | |
1189 | temp->num.links = 1; | |
1190 | temp->num.num = _one_; | |
1191 | temp->num.den = _one_; | |
1192 | return &temp->num; | |
1193 | } | |
1194 | ||
1195 | ||
1196 | void | |
1197 | qfreenum(q) | |
1198 | register NUMBER *q; | |
1199 | { | |
1200 | union allocNode *a; | |
1201 | ||
1202 | if (q == NULL) | |
1203 | return; | |
1204 | zfree(q->num); | |
1205 | zfree(q->den); | |
1206 | a = (union allocNode *) q; | |
1207 | a->link = freeNum; | |
1208 | freeNum = a; | |
1209 | } | |
1210 | ||
1211 | /* END CODE */ |