1 /*
2 * qmath - extended precision rational arithmetic primitive routines
3 *
4 * Copyright (C) 1999-2007,2014,2021 David I. Bell and Ernest Bowen
5 *
6 * Primary author: David I. Bell
7 *
8 * Calc is open software; you can redistribute it and/or modify it under
9 * the terms of the version 2.1 of the GNU Lesser General Public License
10 * as published by the Free Software Foundation.
11 *
12 * Calc is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
15 * Public License for more details.
16 *
17 * A copy of version 2.1 of the GNU Lesser General Public License is
18 * distributed with calc under the filename COPYING-LGPL. You should have
19 * received a copy with calc; if not, write to Free Software Foundation, Inc.
20 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
21 *
22 * Under source code control: 1990/02/15 01:48:21
23 * File existed as early as: before 1990
24 *
25 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
26 */
27
28
29 #include <stdio.h>
30 #include "qmath.h"
31 #include "config.h"
32
33
34 #include "banned.h" /* include after system header <> includes */
35
36
37 NUMBER _qzero_ = { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
38 NUMBER _qone_ = { { _oneval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
39 NUMBER _qtwo_ = { { _twoval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
40 NUMBER _qten_ = { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
41 NUMBER _qnegone_ = { { _oneval_, 1, 1 }, { _oneval_, 1, 0 }, 1, NULL };
42 NUMBER _qonehalf_ = { { _oneval_, 1, 0 }, { _twoval_, 1, 0 }, 1, NULL };
43 NUMBER _qneghalf_ = { { _oneval_, 1, 1 }, { _twoval_, 1, 0 }, 1, NULL };
44 NUMBER _qonesqbase_ = { { _oneval_, 1, 0 }, { _sqbaseval_, 2, 0 }, 1, NULL };
45
46 NUMBER * initnumbs[] = {&_qzero_, &_qone_, &_qtwo_,
47 &_qten_, &_qnegone_, &_qonehalf_, &_qneghalf_,
48 &_qonesqbase_,
49 NULL /* must be last */
50 };
51
52
53 /*
54 * Create another copy of a number.
55 * q2 = qcopy(q1);
56 */
57 NUMBER *
qcopy(NUMBER * q)58 qcopy(NUMBER *q)
59 {
60 register NUMBER *r;
61
62 r = qalloc();
63 r->num.sign = q->num.sign;
64 if (!zisunit(q->num)) {
65 r->num.len = q->num.len;
66 r->num.v = alloc(r->num.len);
67 zcopyval(q->num, r->num);
68 }
69 if (!zisunit(q->den)) {
70 r->den.len = q->den.len;
71 r->den.v = alloc(r->den.len);
72 zcopyval(q->den, r->den);
73 }
74 return r;
75 }
76
77
78 /*
79 * Convert a number to a normal integer.
80 * i = qtoi(q);
81 */
82 long
qtoi(NUMBER * q)83 qtoi(NUMBER *q)
84 {
85 long i;
86 ZVALUE res;
87
88 if (qisint(q))
89 return ztoi(q->num);
90 zquo(q->num, q->den, &res, 0);
91 i = ztoi(res);
92 zfree(res);
93 return i;
94 }
95
96
97 /*
98 * Convert a normal integer into a number.
99 * q = itoq(i);
100 */
101 NUMBER *
itoq(long i)102 itoq(long i)
103 {
104 register NUMBER *q;
105
106 if ((i >= -1) && (i <= 10)) {
107 switch ((int) i) {
108 case 0: q = &_qzero_; break;
109 case 1: q = &_qone_; break;
110 case 2: q = &_qtwo_; break;
111 case 10: q = &_qten_; break;
112 case -1: q = &_qnegone_; break;
113 default: q = NULL;
114 }
115 if (q)
116 return qlink(q);
117 }
118 q = qalloc();
119 itoz(i, &q->num);
120 return q;
121 }
122
123
124 /*
125 * Convert a number to a normal unsigned integer.
126 * u = qtou(q);
127 */
128 FULL
qtou(NUMBER * q)129 qtou(NUMBER *q)
130 {
131 FULL i;
132 ZVALUE res;
133
134 if (qisint(q))
135 return ztou(q->num);
136 zquo(q->num, q->den, &res, 0);
137 i = ztou(res);
138 zfree(res);
139 return i;
140 }
141
142
143 /*
144 * Convert a number to a normal signed integer.
145 * s = qtos(q);
146 */
147 SFULL
qtos(NUMBER * q)148 qtos(NUMBER *q)
149 {
150 SFULL i;
151 ZVALUE res;
152
153 if (qisint(q))
154 return ztos(q->num);
155 zquo(q->num, q->den, &res, 0);
156 i = ztos(res);
157 zfree(res);
158 return i;
159 }
160
161
162 /*
163 * Convert a normal unsigned integer into a number.
164 * q = utoq(i);
165 */
166 NUMBER *
utoq(FULL i)167 utoq(FULL i)
168 {
169 register NUMBER *q;
170
171 if (i <= 10) {
172 switch ((int) i) {
173 case 0: q = &_qzero_; break;
174 case 1: q = &_qone_; break;
175 case 2: q = &_qtwo_; break;
176 case 10: q = &_qten_; break;
177 default: q = NULL;
178 }
179 if (q)
180 return qlink(q);
181 }
182 q = qalloc();
183 utoz(i, &q->num);
184 return q;
185 }
186
187
188 /*
189 * Convert a normal signed integer into a number.
190 * q = stoq(s);
191 */
192 NUMBER *
stoq(SFULL i)193 stoq(SFULL i)
194 {
195 register NUMBER *q;
196
197 if (i <= 10) {
198 switch ((int) i) {
199 case 0: q = &_qzero_; break;
200 case 1: q = &_qone_; break;
201 case 2: q = &_qtwo_; break;
202 case 10: q = &_qten_; break;
203 default: q = NULL;
204 }
205 if (q)
206 return qlink(q);
207 }
208 q = qalloc();
209 stoz(i, &q->num);
210 return q;
211 }
212
213
214 /*
215 * Create a number from the given FULL numerator and denominator.
216 * q = uutoq(inum, iden);
217 */
218 NUMBER *
uutoq(FULL inum,FULL iden)219 uutoq(FULL inum, FULL iden)
220 {
221 register NUMBER *q;
222 FULL d;
223 BOOL sign;
224
225 if (iden == 0) {
226 math_error("Division by zero");
227 /*NOTREACHED*/
228 }
229 if (inum == 0)
230 return qlink(&_qzero_);
231 sign = 0;
232 d = uugcd(inum, iden);
233 inum /= d;
234 iden /= d;
235 if (iden == 1)
236 return utoq(inum);
237 q = qalloc();
238 if (inum != 1)
239 utoz(inum, &q->num);
240 utoz(iden, &q->den);
241 q->num.sign = sign;
242 return q;
243 }
244
245
246 /*
247 * Create a number from the given integral numerator and denominator.
248 * q = iitoq(inum, iden);
249 */
250 NUMBER *
iitoq(long inum,long iden)251 iitoq(long inum, long iden)
252 {
253 register NUMBER *q;
254 long d;
255 BOOL sign;
256
257 if (iden == 0) {
258 math_error("Division by zero");
259 /*NOTREACHED*/
260 }
261 if (inum == 0)
262 return qlink(&_qzero_);
263 sign = 0;
264 if (inum < 0) {
265 sign = 1;
266 inum = -inum;
267 }
268 if (iden < 0) {
269 sign = 1 - sign;
270 iden = -iden;
271 }
272 d = iigcd(inum, iden);
273 inum /= d;
274 iden /= d;
275 if (iden == 1)
276 return itoq(sign ? -inum : inum);
277 q = qalloc();
278 if (inum != 1)
279 itoz(inum, &q->num);
280 itoz(iden, &q->den);
281 q->num.sign = sign;
282 return q;
283 }
284
285
286 /*
287 * Add two numbers to each other.
288 * q3 = qqadd(q1, q2);
289 */
290 NUMBER *
qqadd(NUMBER * q1,NUMBER * q2)291 qqadd(NUMBER *q1, NUMBER *q2)
292 {
293 NUMBER *r;
294 ZVALUE t1, t2, temp, d1, d2, vpd1, upd1;
295
296 if (qiszero(q1))
297 return qlink(q2);
298 if (qiszero(q2))
299 return qlink(q1);
300 r = qalloc();
301 /*
302 * If either number is an integer, then the result is easy.
303 */
304 if (qisint(q1) && qisint(q2)) {
305 zadd(q1->num, q2->num, &r->num);
306 return r;
307 }
308 if (qisint(q2)) {
309 zmul(q1->den, q2->num, &temp);
310 zadd(q1->num, temp, &r->num);
311 zfree(temp);
312 zcopy(q1->den, &r->den);
313 return r;
314 }
315 if (qisint(q1)) {
316 zmul(q2->den, q1->num, &temp);
317 zadd(q2->num, temp, &r->num);
318 zfree(temp);
319 zcopy(q2->den, &r->den);
320 return r;
321 }
322 /*
323 * Both arguments are true fractions, so we need more work.
324 * If the denominators are relatively prime, then the answer is the
325 * straightforward cross product result with no need for reduction.
326 */
327 zgcd(q1->den, q2->den, &d1);
328 if (zisunit(d1)) {
329 zfree(d1);
330 zmul(q1->num, q2->den, &t1);
331 zmul(q1->den, q2->num, &t2);
332 zadd(t1, t2, &r->num);
333 zfree(t1);
334 zfree(t2);
335 zmul(q1->den, q2->den, &r->den);
336 return r;
337 }
338 /*
339 * The calculation is now more complicated.
340 * See Knuth Vol 2 for details.
341 */
342 zquo(q2->den, d1, &vpd1, 0);
343 zquo(q1->den, d1, &upd1, 0);
344 zmul(q1->num, vpd1, &t1);
345 zmul(q2->num, upd1, &t2);
346 zadd(t1, t2, &temp);
347 zfree(t1);
348 zfree(t2);
349 zfree(vpd1);
350 zgcd(temp, d1, &d2);
351 zfree(d1);
352 if (zisunit(d2)) {
353 zfree(d2);
354 r->num = temp;
355 zmul(upd1, q2->den, &r->den);
356 zfree(upd1);
357 return r;
358 }
359 zquo(temp, d2, &r->num, 0);
360 zfree(temp);
361 zquo(q2->den, d2, &temp, 0);
362 zfree(d2);
363 zmul(temp, upd1, &r->den);
364 zfree(temp);
365 zfree(upd1);
366 return r;
367 }
368
369
370 /*
371 * Subtract one number from another.
372 * q3 = qsub(q1, q2);
373 */
374 NUMBER *
qsub(NUMBER * q1,NUMBER * q2)375 qsub(NUMBER *q1, NUMBER *q2)
376 {
377 NUMBER *r;
378
379 if (q1 == q2)
380 return qlink(&_qzero_);
381 if (qiszero(q2))
382 return qlink(q1);
383 if (qisint(q1) && qisint(q2)) {
384 r = qalloc();
385 zsub(q1->num, q2->num, &r->num);
386 return r;
387 }
388 q2 = qneg(q2);
389 if (qiszero(q1))
390 return q2;
391 r = qqadd(q1, q2);
392 qfree(q2);
393 return r;
394 }
395
396
397 /*
398 * Increment a number by one.
399 */
400 NUMBER *
qinc(NUMBER * q)401 qinc(NUMBER *q)
402 {
403 NUMBER *r;
404
405 r = qalloc();
406 if (qisint(q)) {
407 zadd(q->num, _one_, &r->num);
408 return r;
409 }
410 zadd(q->num, q->den, &r->num);
411 zcopy(q->den, &r->den);
412 return r;
413 }
414
415
416 /*
417 * Decrement a number by one.
418 */
419 NUMBER *
qdec(NUMBER * q)420 qdec(NUMBER *q)
421 {
422 NUMBER *r;
423
424 r = qalloc();
425 if (qisint(q)) {
426 zsub(q->num, _one_, &r->num);
427 return r;
428 }
429 zsub(q->num, q->den, &r->num);
430 zcopy(q->den, &r->den);
431 return r;
432 }
433
434
435 /*
436 * Add a normal small integer value to an arbitrary number.
437 */
438 NUMBER *
qaddi(NUMBER * q1,long n)439 qaddi(NUMBER *q1, long n)
440 {
441 NUMBER addnum; /* temporary number */
442 HALF addval[2]; /* value of small number */
443 BOOL neg; /* TRUE if number is neg */
444 #if LONG_BITS > BASEB
445 FULL nf;
446 #endif
447
448 if (n == 0)
449 return qlink(q1);
450 if (n == 1)
451 return qinc(q1);
452 if (n == -1)
453 return qdec(q1);
454 if (qiszero(q1))
455 return itoq(n);
456 addnum.num.sign = 0;
457 addnum.num.v = addval;
458 addnum.den = _one_;
459 neg = (n < 0);
460 if (neg)
461 n = -n;
462 addval[0] = (HALF) n;
463 #if LONG_BITS > BASEB
464 nf = (((FULL) n) >> BASEB);
465 if (nf) {
466 addval[1] = (HALF) nf;
467 addnum.num.len = 2;
468 }
469 #else
470 addnum.num.len = 1;
471 #endif
472 if (neg)
473 return qsub(q1, &addnum);
474 else
475 return qqadd(q1, &addnum);
476 }
477
478
479 /*
480 * Multiply two numbers.
481 * q3 = qmul(q1, q2);
482 */
483 NUMBER *
qmul(NUMBER * q1,NUMBER * q2)484 qmul(NUMBER *q1, NUMBER *q2)
485 {
486 NUMBER *r; /* returned value */
487 ZVALUE n1, n2, d1, d2; /* numerators and denominators */
488 ZVALUE tmp;
489
490 if (qiszero(q1) || qiszero(q2))
491 return qlink(&_qzero_);
492 if (qisone(q1))
493 return qlink(q2);
494 if (qisone(q2))
495 return qlink(q1);
496 if (qisint(q1) && qisint(q2)) { /* easy results if integers */
497 r = qalloc();
498 zmul(q1->num, q2->num, &r->num);
499 return r;
500 }
501 n1 = q1->num;
502 n2 = q2->num;
503 d1 = q1->den;
504 d2 = q2->den;
505 if (ziszero(d1) || ziszero(d2)) {
506 math_error("Division by zero");
507 /*NOTREACHED*/
508 }
509 if (ziszero(n1) || ziszero(n2))
510 return qlink(&_qzero_);
511 if (!zisunit(n1) && !zisunit(d2)) { /* possibly reduce */
512 zgcd(n1, d2, &tmp);
513 if (!zisunit(tmp)) {
514 zequo(q1->num, tmp, &n1);
515 zequo(q2->den, tmp, &d2);
516 }
517 zfree(tmp);
518 }
519 if (!zisunit(n2) && !zisunit(d1)) { /* again possibly reduce */
520 zgcd(n2, d1, &tmp);
521 if (!zisunit(tmp)) {
522 zequo(q2->num, tmp, &n2);
523 zequo(q1->den, tmp, &d1);
524 }
525 zfree(tmp);
526 }
527 r = qalloc();
528 zmul(n1, n2, &r->num);
529 zmul(d1, d2, &r->den);
530 if (q1->num.v != n1.v)
531 zfree(n1);
532 if (q1->den.v != d1.v)
533 zfree(d1);
534 if (q2->num.v != n2.v)
535 zfree(n2);
536 if (q2->den.v != d2.v)
537 zfree(d2);
538 return r;
539 }
540
541
542 /*
543 * Multiply a number by a small integer.
544 * q2 = qmuli(q1, n);
545 */
546 NUMBER *
qmuli(NUMBER * q,long n)547 qmuli(NUMBER *q, long n)
548 {
549 NUMBER *r;
550 long d; /* gcd of multiplier and denominator */
551 int sign;
552
553 if ((n == 0) || qiszero(q))
554 return qlink(&_qzero_);
555 if (n == 1)
556 return qlink(q);
557 r = qalloc();
558 if (qisint(q)) {
559 zmuli(q->num, n, &r->num);
560 return r;
561 }
562 sign = 1;
563 if (n < 0) {
564 n = -n;
565 sign = -1;
566 }
567 d = zmodi(q->den, n);
568 d = iigcd(d, n);
569 zmuli(q->num, (n * sign) / d, &r->num);
570 (void) zdivi(q->den, d, &r->den);
571 return r;
572 }
573
574
575 /*
576 * Divide two numbers (as fractions).
577 * q3 = qqdiv(q1, q2);
578 */
579 NUMBER *
qqdiv(NUMBER * q1,NUMBER * q2)580 qqdiv(NUMBER *q1, NUMBER *q2)
581 {
582 NUMBER temp;
583
584 if (qiszero(q2)) {
585 math_error("Division by zero");
586 /*NOTREACHED*/
587 }
588 if ((q1 == q2) || !qcmp(q1, q2))
589 return qlink(&_qone_);
590 if (qisone(q1))
591 return qinv(q2);
592 temp.num = q2->den;
593 temp.den = q2->num;
594 temp.num.sign = temp.den.sign;
595 temp.den.sign = 0;
596 temp.links = 1;
597 return qmul(q1, &temp);
598 }
599
600
601 /*
602 * Divide a number by a small integer.
603 * q2 = qdivi(q1, n);
604 */
605 NUMBER *
qdivi(NUMBER * q,long n)606 qdivi(NUMBER *q, long n)
607 {
608 NUMBER *r;
609 long d; /* gcd of divisor and numerator */
610 int sign;
611
612 if (n == 0) {
613 math_error("Division by zero");
614 /*NOTREACHED*/
615 }
616 if ((n == 1) || qiszero(q))
617 return qlink(q);
618 sign = 1;
619 if (n < 0) {
620 n = -n;
621 sign = -1;
622 }
623 r = qalloc();
624 d = zmodi(q->num, n);
625 d = iigcd(d, n);
626 (void) zdivi(q->num, d * sign, &r->num);
627 zmuli(q->den, n / d, &r->den);
628 return r;
629 }
630
631
632 /*
633 * Return the integer quotient of a pair of numbers
634 * If q1/q2 is an integer qquo(q1, q2) returns this integer
635 * If q2 is zero, zero is returned
636 * In other cases whether rounding is down, up, towards zero, etc.
637 * is determined by rnd.
638 */
639 NUMBER *
qquo(NUMBER * q1,NUMBER * q2,long rnd)640 qquo(NUMBER *q1, NUMBER *q2, long rnd)
641 {
642 ZVALUE tmp, tmp1, tmp2;
643 NUMBER *q;
644
645 if (qiszero(q1) || qiszero(q2))
646 return qlink(&_qzero_);
647 if (qisint(q1) && qisint(q2)) {
648 zquo(q1->num, q2->num, &tmp, rnd);
649 } else {
650 zmul(q1->num, q2->den, &tmp1);
651 zmul(q2->num, q1->den, &tmp2);
652 zquo(tmp1, tmp2, &tmp, rnd);
653 zfree(tmp1);
654 zfree(tmp2);
655 }
656 if (ziszero(tmp)) {
657 zfree(tmp);
658 return qlink(&_qzero_);
659 }
660 q = qalloc();
661 q->num = tmp;
662 return q;
663 }
664
665
666 /*
667 * Return the absolute value of a number.
668 * q2 = qqabs(q1);
669 */
670 NUMBER *
qqabs(NUMBER * q)671 qqabs(NUMBER *q)
672 {
673 register NUMBER *r;
674
675 if (q->num.sign == 0)
676 return qlink(q);
677 r = qalloc();
678 if (!zisunit(q->num))
679 zcopy(q->num, &r->num);
680 if (!zisunit(q->den))
681 zcopy(q->den, &r->den);
682 r->num.sign = 0;
683 return r;
684 }
685
686
687 /*
688 * Negate a number.
689 * q2 = qneg(q1);
690 */
691 NUMBER *
qneg(NUMBER * q)692 qneg(NUMBER *q)
693 {
694 register NUMBER *r;
695
696 if (qiszero(q))
697 return qlink(&_qzero_);
698 r = qalloc();
699 if (!zisunit(q->num))
700 zcopy(q->num, &r->num);
701 if (!zisunit(q->den))
702 zcopy(q->den, &r->den);
703 r->num.sign = !q->num.sign;
704 return r;
705 }
706
707
708 /*
709 * Return the sign of a number (-1, 0 or 1)
710 */
711 NUMBER *
qsign(NUMBER * q)712 qsign(NUMBER *q)
713 {
714 if (qiszero(q))
715 return qlink(&_qzero_);
716 if (qisneg(q))
717 return qlink(&_qnegone_);
718 return qlink(&_qone_);
719 }
720
721
722 /*
723 * Invert a number.
724 * q2 = qinv(q1);
725 */
726 NUMBER *
qinv(NUMBER * q)727 qinv(NUMBER *q)
728 {
729 register NUMBER *r;
730
731 if (qisunit(q)) {
732 r = (qisneg(q) ? &_qnegone_ : &_qone_);
733 return qlink(r);
734 }
735 if (qiszero(q)) {
736 math_error("Division by zero");
737 /*NOTREACHED*/
738 }
739 r = qalloc();
740 if (!zisunit(q->num))
741 zcopy(q->num, &r->den);
742 if (!zisunit(q->den))
743 zcopy(q->den, &r->num);
744 r->num.sign = q->num.sign;
745 r->den.sign = 0;
746 return r;
747 }
748
749
750 /*
751 * Return just the numerator of a number.
752 * q2 = qnum(q1);
753 */
754 NUMBER *
qnum(NUMBER * q)755 qnum(NUMBER *q)
756 {
757 register NUMBER *r;
758
759 if (qisint(q))
760 return qlink(q);
761 if (zisunit(q->num)) {
762 r = (qisneg(q) ? &_qnegone_ : &_qone_);
763 return qlink(r);
764 }
765 r = qalloc();
766 zcopy(q->num, &r->num);
767 return r;
768 }
769
770
771 /*
772 * Return just the denominator of a number.
773 * q2 = qden(q1);
774 */
775 NUMBER *
qden(NUMBER * q)776 qden(NUMBER *q)
777 {
778 register NUMBER *r;
779
780 if (qisint(q))
781 return qlink(&_qone_);
782 r = qalloc();
783 zcopy(q->den, &r->num);
784 return r;
785 }
786
787
788 /*
789 * Return the fractional part of a number.
790 * q2 = qfrac(q1);
791 */
792 NUMBER *
qfrac(NUMBER * q)793 qfrac(NUMBER *q)
794 {
795 register NUMBER *r;
796
797 if (qisint(q))
798 return qlink(&_qzero_);
799 if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
800 (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
801 return qlink(q);
802 r = qalloc();
803 zmod(q->num, q->den, &r->num, 2);
804 zcopy(q->den, &r->den);
805 return r;
806 }
807
808
809 /*
810 * Return the integral part of a number.
811 * q2 = qint(q1);
812 */
813 NUMBER *
qint(NUMBER * q)814 qint(NUMBER *q)
815 {
816 register NUMBER *r;
817
818 if (qisint(q))
819 return qlink(q);
820 if ((q->num.len < q->den.len) || ((q->num.len == q->den.len) &&
821 (q->num.v[q->num.len - 1] < q->den.v[q->num.len - 1])))
822 return qlink(&_qzero_);
823 r = qalloc();
824 zquo(q->num, q->den, &r->num, 2);
825 return r;
826 }
827
828
829 /*
830 * Compute the square of a number.
831 */
832 NUMBER *
qsquare(NUMBER * q)833 qsquare(NUMBER *q)
834 {
835 ZVALUE num, zden;
836
837 if (qiszero(q))
838 return qlink(&_qzero_);
839 if (qisunit(q))
840 return qlink(&_qone_);
841 num = q->num;
842 zden = q->den;
843 q = qalloc();
844 if (!zisunit(num))
845 zsquare(num, &q->num);
846 if (!zisunit(zden))
847 zsquare(zden, &q->den);
848 return q;
849 }
850
851
852 /*
853 * Shift an integer by a given number of bits. This multiplies the number
854 * by the appropriate power of two. Positive numbers shift left, negative
855 * ones shift right. Low bits are truncated when shifting right.
856 */
857 NUMBER *
qshift(NUMBER * q,long n)858 qshift(NUMBER *q, long n)
859 {
860 register NUMBER *r;
861
862 if (qisfrac(q)) {
863 math_error("Shift of non-integer");
864 /*NOTREACHED*/
865 }
866 if (qiszero(q) || (n == 0))
867 return qlink(q);
868 if (n <= -(q->num.len * BASEB))
869 return qlink(&_qzero_);
870 r = qalloc();
871 zshift(q->num, n, &r->num);
872 return r;
873 }
874
875
876 /*
877 * Scale a number by a power of two, as in:
878 * ans = q * 2^n.
879 * This is similar to shifting, except that fractions work.
880 */
881 NUMBER *
qscale(NUMBER * q,long pow)882 qscale(NUMBER *q, long pow)
883 {
884 long numshift, denshift, tmp;
885 NUMBER *r;
886
887 if (qiszero(q) || (pow == 0))
888 return qlink(q);
889 numshift = zisodd(q->num) ? 0 : zlowbit(q->num);
890 denshift = zisodd(q->den) ? 0 : zlowbit(q->den);
891 if (pow > 0) {
892 tmp = pow;
893 if (tmp > denshift)
894 tmp = denshift;
895 denshift = -tmp;
896 numshift = (pow - tmp);
897 } else {
898 pow = -pow;
899 tmp = pow;
900 if (tmp > numshift)
901 tmp = numshift;
902 numshift = -tmp;
903 denshift = (pow - tmp);
904 }
905 r = qalloc();
906 if (numshift)
907 zshift(q->num, numshift, &r->num);
908 else
909 zcopy(q->num, &r->num);
910 if (denshift)
911 zshift(q->den, denshift, &r->den);
912 else
913 zcopy(q->den, &r->den);
914 return r;
915 }
916
917
918 /*
919 * Return the minimum of two numbers.
920 */
921 NUMBER *
qmin(NUMBER * q1,NUMBER * q2)922 qmin(NUMBER *q1, NUMBER *q2)
923 {
924 if (q1 == q2)
925 return qlink(q1);
926 if (qrel(q1, q2) > 0)
927 q1 = q2;
928 return qlink(q1);
929 }
930
931
932 /*
933 * Return the maximum of two numbers.
934 */
935 NUMBER *
qmax(NUMBER * q1,NUMBER * q2)936 qmax(NUMBER *q1, NUMBER *q2)
937 {
938 if (q1 == q2)
939 return qlink(q1);
940 if (qrel(q1, q2) < 0)
941 q1 = q2;
942 return qlink(q1);
943 }
944
945
946 /*
947 * Perform the bitwise OR of two integers.
948 */
949 NUMBER *
qor(NUMBER * q1,NUMBER * q2)950 qor(NUMBER *q1, NUMBER *q2)
951 {
952 register NUMBER *r;
953 NUMBER *q1tmp, *q2tmp, *q;
954
955 if (qisfrac(q1) || qisfrac(q2)) {
956 math_error("Non-integers for bitwise or");
957 /*NOTREACHED*/
958 }
959 if (qcmp(q1,q2) == 0 || qiszero(q2))
960 return qlink(q1);
961 if (qiszero(q1))
962 return qlink(q2);
963 if (qisneg(q1)) {
964 q1tmp = qcomp(q1);
965 if (qisneg(q2)) {
966 q2tmp = qcomp(q2);
967 q = qand(q1tmp,q2tmp);
968 r = qcomp(q);
969 qfree(q1tmp);
970 qfree(q2tmp);
971 qfree(q);
972 return r;
973 }
974 q = qandnot(q1tmp, q2);
975 qfree(q1tmp);
976 r = qcomp(q);
977 qfree(q);
978 return r;
979 }
980 if (qisneg(q2)) {
981 q2tmp = qcomp(q2);
982 q = qandnot(q2tmp, q1);
983 qfree(q2tmp);
984 r = qcomp(q);
985 qfree(q);
986 return r;
987 }
988 r = qalloc();
989 zor(q1->num, q2->num, &r->num);
990 return r;
991 }
992
993
994 /*
995 * Perform the bitwise AND of two integers.
996 */
997 NUMBER *
qand(NUMBER * q1,NUMBER * q2)998 qand(NUMBER *q1, NUMBER *q2)
999 {
1000 register NUMBER *r;
1001 NUMBER *q1tmp, *q2tmp, *q;
1002 ZVALUE res;
1003
1004 if (qisfrac(q1) || qisfrac(q2)) {
1005 math_error("Non-integers for bitwise and");
1006 /*NOTREACHED*/
1007 }
1008 if (qcmp(q1, q2) == 0)
1009 return qlink(q1);
1010 if (qiszero(q1) || qiszero(q2))
1011 return qlink(&_qzero_);
1012 if (qisneg(q1)) {
1013 q1tmp = qcomp(q1);
1014 if (qisneg(q2)) {
1015 q2tmp = qcomp(q2);
1016 q = qor(q1tmp, q2tmp);
1017 qfree(q1tmp);
1018 qfree(q2tmp);
1019 r = qcomp(q);
1020 qfree(q);
1021 return r;
1022 }
1023 r = qandnot(q2, q1tmp);
1024 qfree(q1tmp);
1025 return r;
1026 }
1027 if (qisneg(q2)) {
1028 q2tmp = qcomp(q2);
1029 r = qandnot(q1, q2tmp);
1030 qfree(q2tmp);
1031 return r;
1032 }
1033 zand(q1->num, q2->num, &res);
1034 if (ziszero(res)) {
1035 zfree(res);
1036 return qlink(&_qzero_);
1037 }
1038 r = qalloc();
1039 r->num = res;
1040 return r;
1041 }
1042
1043
1044 /*
1045 * Perform the bitwise XOR of two integers.
1046 */
1047 NUMBER *
qxor(NUMBER * q1,NUMBER * q2)1048 qxor(NUMBER *q1, NUMBER *q2)
1049 {
1050 register NUMBER *r;
1051 NUMBER *q1tmp, *q2tmp, *q;
1052 ZVALUE res;
1053
1054 if (qisfrac(q1) || qisfrac(q2)) {
1055 math_error("Non-integers for bitwise xor");
1056 /*NOTREACHED*/
1057 }
1058 if (qcmp(q1,q2) == 0)
1059 return qlink(&_qzero_);
1060 if (qiszero(q1))
1061 return qlink(q2);
1062 if (qiszero(q2))
1063 return qlink(q1);
1064 if (qisneg(q1)) {
1065 q1tmp = qcomp(q1);
1066 if (qisneg(q2)) {
1067 q2tmp = qcomp(q2);
1068 r = qxor(q1tmp, q2tmp);
1069 qfree(q1tmp);
1070 qfree(q2tmp);
1071 return r;
1072 }
1073 q = qxor(q1tmp, q2);
1074 qfree(q1tmp);
1075 r = qcomp(q);
1076 qfree(q);
1077 return r;
1078 }
1079 if (qisneg(q2)) {
1080 q2tmp = qcomp(q2);
1081 q = qxor(q1, q2tmp);
1082 qfree(q2tmp);
1083 r = qcomp(q);
1084 qfree(q);
1085 return r;
1086 }
1087 zxor(q1->num, q2->num, &res);
1088 if (ziszero(res)) {
1089 zfree(res);
1090 return qlink(&_qzero_);
1091 }
1092 r = qalloc();
1093 r->num = res;
1094 return r;
1095 }
1096
1097
1098 /*
1099 * Perform the bitwise AND-NOT of two integers.
1100 */
1101 NUMBER *
qandnot(NUMBER * q1,NUMBER * q2)1102 qandnot(NUMBER *q1, NUMBER *q2)
1103 {
1104 register NUMBER *r;
1105 NUMBER *q1tmp, *q2tmp, *q;
1106
1107 if (qisfrac(q1) || qisfrac(q2)) {
1108 math_error("Non-integers for bitwise xor");
1109 /*NOTREACHED*/
1110 }
1111 if (qcmp(q1,q2) == 0 || qiszero(q1))
1112 return qlink(&_qzero_);
1113 if (qiszero(q2))
1114 return qlink(q1);
1115 if (qisneg(q1)) {
1116 q1tmp = qcomp(q1);
1117 if (qisneg(q2)) {
1118 q2tmp = qcomp(q2);
1119 r = qandnot(q2tmp, q1tmp);
1120 qfree(q1tmp);
1121 qfree(q2tmp);
1122 return r;
1123 }
1124 q = qor(q1tmp, q2);
1125 qfree(q1tmp);
1126 r = qcomp(q);
1127 qfree(q);
1128 return r;
1129 }
1130 if (qisneg(q2)) {
1131 q2tmp = qcomp(q2);
1132 r = qand(q1, q2tmp);
1133 qfree(q2tmp);
1134 return r;
1135 }
1136 r = qalloc();
1137 zandnot(q1->num, q2->num, &r->num);
1138 return r;
1139 }
1140
1141 /*
1142 * Return the bitwise "complement" of a number. This is - q -1 if q is an
1143 * integer, - q otherwise.
1144 */
1145 NUMBER *
qcomp(NUMBER * q)1146 qcomp(NUMBER *q)
1147 {
1148 NUMBER *qtmp;
1149 NUMBER *res;
1150
1151 if (qiszero(q))
1152 return qlink(&_qnegone_);
1153 if (qisnegone(q))
1154 return qlink(&_qzero_);
1155 qtmp = qneg(q);
1156 if (qisfrac(q))
1157 return qtmp;
1158 res = qdec(qtmp);
1159 qfree(qtmp);
1160 return res;
1161 }
1162
1163
1164 /*
1165 * Return the number whose binary representation only has the specified
1166 * bit set (counting from zero). This thus produces a given power of two.
1167 */
1168 NUMBER *
qbitvalue(long n)1169 qbitvalue(long n)
1170 {
1171 register NUMBER *r;
1172
1173 if (n == 0)
1174 return qlink(&_qone_);
1175 r = qalloc();
1176 if (n > 0)
1177 zbitvalue(n, &r->num);
1178 else
1179 zbitvalue(-n, &r->den);
1180 return r;
1181 }
1182
1183 /*
1184 * Return 10^n
1185 */
1186 NUMBER *
qtenpow(long n)1187 qtenpow(long n)
1188 {
1189 register NUMBER *r;
1190
1191 if (n == 0)
1192 return qlink(&_qone_);
1193 r = qalloc();
1194 if (n > 0)
1195 ztenpow(n, &r->num);
1196 else
1197 ztenpow(-n, &r->den);
1198 return r;
1199 }
1200
1201
1202 /*
1203 * Return the precision of a number (usually for examining an epsilon value).
1204 * The precision of a number e less than 1 is the positive
1205 * integer p for which e = 2^-p * f, where 1 <= f < 2.
1206 * Numbers greater than or equal to one have a precision of zero.
1207 * For example, the precision of e is 6 if 1/64 <= e < 1/32.
1208 */
1209 long
qprecision(NUMBER * q)1210 qprecision(NUMBER *q)
1211 {
1212 long r;
1213
1214 if (qiszero(q) || qisneg(q)) {
1215 math_error("Non-positive number for precision");
1216 /*NOTREACHED*/
1217 }
1218 r = - qilog2(q);
1219 return (r < 0 ? 0 : r);
1220 }
1221
1222
1223 /*
1224 * Determine whether or not one number exactly divides another one.
1225 * Returns TRUE if the first number is an integer multiple of the second one.
1226 */
1227 BOOL
qdivides(NUMBER * q1,NUMBER * q2)1228 qdivides(NUMBER *q1, NUMBER *q2)
1229 {
1230 if (qiszero(q1))
1231 return TRUE;
1232 if (qisint(q1) && qisint(q2)) {
1233 if (qisunit(q2))
1234 return TRUE;
1235 return zdivides(q1->num, q2->num);
1236 }
1237 return zdivides(q1->num, q2->num) && zdivides(q2->den, q1->den);
1238 }
1239
1240
1241 /*
1242 * Compare two numbers and return an integer indicating their relative size.
1243 * i = qrel(q1, q2);
1244 */
1245 FLAG
qrel(NUMBER * q1,NUMBER * q2)1246 qrel(NUMBER *q1, NUMBER *q2)
1247 {
1248 ZVALUE z1, z2;
1249 long wc1, wc2;
1250 int sign;
1251 int z1f = 0, z2f = 0;
1252
1253 if (q1 == q2)
1254 return 0;
1255 sign = q2->num.sign - q1->num.sign;
1256 if (sign)
1257 return sign;
1258 if (qiszero(q2))
1259 return !qiszero(q1);
1260 if (qiszero(q1))
1261 return -1;
1262 /*
1263 * Make a quick comparison by calculating the number of words
1264 * resulting as if we multiplied through by the denominators,
1265 * and then comparing the word counts.
1266 */
1267 sign = 1;
1268 if (qisneg(q1))
1269 sign = -1;
1270 wc1 = q1->num.len + q2->den.len;
1271 wc2 = q2->num.len + q1->den.len;
1272 if (wc1 < wc2 - 1)
1273 return -sign;
1274 if (wc2 < wc1 - 1)
1275 return sign;
1276 /*
1277 * Quick check failed, must actually do the full comparison.
1278 */
1279 if (zisunit(q2->den)) {
1280 z1 = q1->num;
1281 } else if (zisone(q1->num)) {
1282 z1 = q2->den;
1283 } else {
1284 z1f = 1;
1285 zmul(q1->num, q2->den, &z1);
1286 }
1287 if (zisunit(q1->den)) {
1288 z2 = q2->num;
1289 } else if (zisone(q2->num)) {
1290 z2 = q1->den;
1291 } else {
1292 z2f = 1;
1293 zmul(q2->num, q1->den, &z2);
1294 }
1295 sign = zrel(z1, z2);
1296 if (z1f)
1297 zfree(z1);
1298 if (z2f)
1299 zfree(z2);
1300 return sign;
1301 }
1302
1303
1304 /*
1305 * Compare two numbers to see if they are equal.
1306 * This differs from qrel in that the numbers are not ordered.
1307 * Returns TRUE if they differ.
1308 */
1309 BOOL
qcmp(NUMBER * q1,NUMBER * q2)1310 qcmp(NUMBER *q1, NUMBER *q2)
1311 {
1312 if (q1 == q2)
1313 return FALSE;
1314 if ((q1->num.sign != q2->num.sign) || (q1->num.len != q2->num.len) ||
1315 (q1->den.len != q2->den.len) || (*q1->num.v != *q2->num.v) ||
1316 (*q1->den.v != *q2->den.v))
1317 return TRUE;
1318 if (zcmp(q1->num, q2->num))
1319 return TRUE;
1320 if (qisint(q1))
1321 return FALSE;
1322 return zcmp(q1->den, q2->den);
1323 }
1324
1325
1326 /*
1327 * Compare a number against a normal small integer.
1328 * Returns 1, 0, or -1, according to whether the first number is greater,
1329 * equal, or less than the second number.
1330 * res = qreli(q, n);
1331 */
1332 FLAG
qreli(NUMBER * q,long n)1333 qreli(NUMBER *q, long n)
1334 {
1335 ZVALUE z1, z2;
1336 FLAG res;
1337
1338 if (qiszero(q))
1339 return ((n > 0) ? -1 : (n < 0));
1340
1341 if (n == 0)
1342 return (q->num.sign ? -1 : 0);
1343
1344 if (q->num.sign != (n < 0))
1345 return ((n < 0) ? 1 : -1);
1346
1347 itoz(n, &z1);
1348
1349 if (qisfrac(q)) {
1350 zmul(q->den, z1, &z2);
1351 zfree(z1);
1352 z1 = z2;
1353 }
1354
1355 res = zrel(q->num, z1);
1356 zfree(z1);
1357 return res;
1358 }
1359
1360
1361 /*
1362 * Compare a number against a small integer to see if they are equal.
1363 * Returns TRUE if they differ.
1364 */
1365 BOOL
qcmpi(NUMBER * q,long n)1366 qcmpi(NUMBER *q, long n)
1367 {
1368 long len;
1369 #if LONG_BITS > BASEB
1370 FULL nf;
1371 #endif
1372
1373 len = q->num.len;
1374 if (qisfrac(q) || (q->num.sign != (n < 0)))
1375 return TRUE;
1376 if (n < 0)
1377 n = -n;
1378 if (((HALF)(n)) != q->num.v[0])
1379 return TRUE;
1380 #if LONG_BITS > BASEB
1381 nf = ((FULL) n) >> BASEB;
1382 return ((nf != 0 || len > 1) && (len != 2 || nf != q->num.v[1]));
1383 #else
1384 return (len > 1);
1385 #endif
1386 }
1387
1388
1389 /*
1390 * Number node allocation routines
1391 */
1392
1393 #define NNALLOC 1000
1394
1395
1396 STATIC NUMBER *freeNum = NULL;
1397 STATIC NUMBER **firstNums = NULL;
1398 STATIC long blockcount = 0;
1399
1400
1401 NUMBER *
qalloc(void)1402 qalloc(void)
1403 {
1404 NUMBER *temp;
1405 NUMBER ** newfn;
1406
1407 if (freeNum == NULL) {
1408 freeNum = (NUMBER *) malloc(sizeof (NUMBER) * NNALLOC);
1409 if (freeNum == NULL) {
1410 math_error("Not enough memory");
1411 /*NOTREACHED*/
1412 }
1413 freeNum[NNALLOC - 1].next = NULL;
1414 freeNum[NNALLOC - 1].links = 0;
1415
1416 /*
1417 * We prevent the temp pointer from walking behind freeNum
1418 * by stopping one short of the end and running the loop one
1419 * more time.
1420 *
1421 * We would stop the loop with just temp >= freeNum, but
1422 * doing this helps make code checkers such as insure happy.
1423 */
1424 for (temp = freeNum + NNALLOC - 2; temp > freeNum; --temp) {
1425 temp->next = temp + 1;
1426 temp->links = 0;
1427 }
1428 /* run the loop manually one last time */
1429 temp->next = temp + 1;
1430 temp->links = 0;
1431
1432 blockcount++;
1433 if (firstNums == NULL) {
1434 newfn = (NUMBER **) malloc(blockcount * sizeof(NUMBER *));
1435 } else {
1436 newfn = (NUMBER **)
1437 realloc(firstNums, blockcount * sizeof(NUMBER *));
1438 }
1439 if (newfn == NULL) {
1440 math_error("Cannot allocate new number block");
1441 /*NOTREACHED*/
1442 }
1443 firstNums = newfn;
1444 firstNums[blockcount - 1] = freeNum;
1445 }
1446 temp = freeNum;
1447 freeNum = temp->next;
1448 temp->links = 1;
1449 temp->num = _one_;
1450 temp->den = _one_;
1451 return temp;
1452 }
1453
1454
1455 void
qfreenum(NUMBER * q)1456 qfreenum(NUMBER *q)
1457 {
1458 if (q == NULL) {
1459 math_error("Calling qfreenum with null argument!!!");
1460 /*NOTREACHED*/
1461 }
1462 if (q->links != 0) {
1463 math_error("Calling qfreenum with non-zero links!!!");
1464 /*NOTREACHED*/
1465 }
1466 zfree(q->num);
1467 zfree(q->den);
1468 q->next = freeNum;
1469 freeNum = q;
1470 }
1471
1472 void
shownumbers(void)1473 shownumbers(void)
1474 {
1475 NUMBER *vp;
1476 long i, j, k;
1477 long count = 0;
1478
1479 printf("Index Links Digits Value\n");
1480 printf("----- ----- ------ -----\n");
1481
1482 for (i = 0, k = 0; initnumbs[i] != NULL; i++) {
1483 count++;
1484 vp = initnumbs[i];
1485 printf("%6ld %4ld ", k++, vp->links);
1486 fitprint(vp, 40);
1487 printf("\n");
1488 }
1489
1490 for (i = 0; i < blockcount; i++) {
1491 vp = firstNums[i];
1492 for (j = 0; j < NNALLOC; j++, k++, vp++) {
1493 if (vp->links > 0) {
1494 count++;
1495 printf("%6ld %4ld ", k, vp->links);
1496 fitprint(vp, 40);
1497 printf("\n");
1498 }
1499 }
1500 }
1501 printf("\nNumber: %ld\n", count);
1502 }
1503