1 /*
2 * zmath - extended precision integral arithmetic primitives
3 *
4 * Copyright (C) 1999-2007,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:28
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 "alloc.h"
31 #include "zmath.h"
32
33
34 #include "banned.h" /* include after system header <> includes */
35
36
37 HALF _zeroval_[] = { 0 };
38 ZVALUE _zero_ = { _zeroval_, 1, 0};
39
40 HALF _oneval_[] = { 1 };
41 ZVALUE _one_ = { _oneval_, 1, 0 };
42 ZVALUE _neg_one_ = { _oneval_, 1, 1 };
43
44 HALF _twoval_[] = { 2 };
45 ZVALUE _two_ = { _twoval_, 1, 0 };
46
47 HALF _tenval_[] = { 10 };
48 ZVALUE _ten_ = { _tenval_, 1, 0 };
49
50 HALF _sqbaseval_[] = { 0, 1 };
51 ZVALUE _sqbase_ = { _sqbaseval_, 2, 0 };
52
53 HALF _pow4baseval_[] = { 0, 0, 1 };
54 ZVALUE _pow4base_ = { _pow4baseval_, 4, 0 };
55
56 HALF _pow8baseval_[] = { 0, 0, 0, 0, 1 };
57 ZVALUE _pow8base_ = { _pow8baseval_, 4, 0 };
58
59 /*
60 * 2^64 as a ZVALUE
61 */
62 #if BASEB == 32
63 ZVALUE _b32_ = { _sqbaseval_, 2, 0 };
64 ZVALUE _b64_ = { _pow4baseval_, 3, 0 };
65 #elif BASEB == 16
66 ZVALUE _b32_ = { _pow4baseval_, 3, 0 };
67 ZVALUE _b64_ = { _pow8baseval_, 5, 0 };
68 #else
69 -=@=- BASEB not 16 or 32 -=@=-
70 #endif
71
72 /*
73 * ZVALUE - values that should not be freed
74 */
75 HALF *half_tbl[] = {
76 _zeroval_,
77 _oneval_,
78 _twoval_,
79 _tenval_,
80 _sqbaseval_,
81 _pow4baseval_,
82 _pow8baseval_,
83 NULL /* must be last */
84 };
85
86
87 /*
88 * highhalf[i] - masks off the upper i bits of a HALF
89 * rhighhalf[i] - masks off the upper BASEB-i bits of a HALF
90 * lowhalf[i] - masks off the upper i bits of a HALF
91 * rlowhalf[i] - masks off the upper BASEB-i bits of a HALF
92 * bitmask[i] - (1 << i) for 0 <= i <= BASEB*2
93 */
94 HALF highhalf[BASEB+1] = {
95 #if BASEB == 32
96 0x00000000,
97 0x80000000, 0xC0000000, 0xE0000000, 0xF0000000,
98 0xF8000000, 0xFC000000, 0xFE000000, 0xFF000000,
99 0xFF800000, 0xFFC00000, 0xFFE00000, 0xFFF00000,
100 0xFFF80000, 0xFFFC0000, 0xFFFE0000, 0xFFFF0000,
101 0xFFFF8000, 0xFFFFC000, 0xFFFFE000, 0xFFFFF000,
102 0xFFFFF800, 0xFFFFFC00, 0xFFFFFE00, 0xFFFFFF00,
103 0xFFFFFF80, 0xFFFFFFC0, 0xFFFFFFE0, 0xFFFFFFF0,
104 0xFFFFFFF8, 0xFFFFFFFC, 0xFFFFFFFE, 0xFFFFFFFF
105 #elif BASEB == 16
106 0x0000,
107 0x8000, 0xC000, 0xE000, 0xF000,
108 0xF800, 0xFC00, 0xFE00, 0xFF00,
109 0xFF80, 0xFFC0, 0xFFE0, 0xFFF0,
110 0xFFF8, 0xFFFC, 0xFFFE, 0xFFFF
111 #else
112 -=@=- BASEB not 16 or 32 -=@=-
113 #endif
114 };
115 HALF rhighhalf[BASEB+1] = {
116 #if BASEB == 32
117 0xFFFFFFFF, 0xFFFFFFFE, 0xFFFFFFFC, 0xFFFFFFF8,
118 0xFFFFFFF0, 0xFFFFFFE0, 0xFFFFFFC0, 0xFFFFFF80,
119 0xFFFFFF00, 0xFFFFFE00, 0xFFFFFC00, 0xFFFFF800,
120 0xFFFFF000, 0xFFFFE000, 0xFFFFC000, 0xFFFF8000,
121 0xFFFF0000, 0xFFFE0000, 0xFFFC0000, 0xFFF80000,
122 0xFFF00000, 0xFFE00000, 0xFFC00000, 0xFF800000,
123 0xFF000000, 0xFE000000, 0xFC000000, 0xF8000000,
124 0xF0000000, 0xE0000000, 0xC0000000, 0x80000000,
125 0x00000000
126 #elif BASEB == 16
127 0xFFFF, 0xFFFE, 0xFFFC, 0xFFF8,
128 0xFFF0, 0xFFE0, 0xFFC0, 0xFF80,
129 0xFF00, 0xFE00, 0xFC00, 0xF800,
130 0xF000, 0xE000, 0xC000, 0x8000,
131 0x0000
132 #else
133 -=@=- BASEB not 16 or 32 -=@=-
134 #endif
135 };
136 HALF lowhalf[BASEB+1] = {
137 0x0,
138 0x1, 0x3, 0x7, 0xF,
139 0x1F, 0x3F, 0x7F, 0xFF,
140 0x1FF, 0x3FF, 0x7FF, 0xFFF,
141 0x1FFF, 0x3FFF, 0x7FFF, 0xFFFF
142 #if BASEB == 32
143 ,0x1FFFF, 0x3FFFF, 0x7FFFF, 0xFFFFF,
144 0x1FFFFF, 0x3FFFFF, 0x7FFFFF, 0xFFFFFF,
145 0x1FFFFFF, 0x3FFFFFF, 0x7FFFFFF, 0xFFFFFFF,
146 0x1FFFFFFF, 0x3FFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF
147 #endif
148 };
149 HALF rlowhalf[BASEB+1] = {
150 #if BASEB == 32
151 0xFFFFFFFF, 0x7FFFFFFF, 0x3FFFFFFF, 0x1FFFFFFF,
152 0xFFFFFFF, 0x7FFFFFF, 0x3FFFFFF, 0x1FFFFFF,
153 0xFFFFFF, 0x7FFFFF, 0x3FFFFF, 0x1FFFFF,
154 0xFFFFF, 0x7FFFF, 0x3FFFF, 0x1FFFF,
155 #endif
156 0xFFFF, 0x7FFF, 0x3FFF, 0x1FFF,
157 0xFFF, 0x7FF, 0x3FF, 0x1FF,
158 0xFF, 0x7F, 0x3F, 0x1F,
159 0xF, 0x7, 0x3, 0x1,
160 0x0
161 };
162 HALF bitmask[(2*BASEB)+1] = {
163 #if BASEB == 32
164 0x00000001, 0x00000002, 0x00000004, 0x00000008,
165 0x00000010, 0x00000020, 0x00000040, 0x00000080,
166 0x00000100, 0x00000200, 0x00000400, 0x00000800,
167 0x00001000, 0x00002000, 0x00004000, 0x00008000,
168 0x00010000, 0x00020000, 0x00040000, 0x00080000,
169 0x00100000, 0x00200000, 0x00400000, 0x00800000,
170 0x01000000, 0x02000000, 0x04000000, 0x08000000,
171 0x10000000, 0x20000000, 0x40000000, 0x80000000,
172 0x00000001, 0x00000002, 0x00000004, 0x00000008,
173 0x00000010, 0x00000020, 0x00000040, 0x00000080,
174 0x00000100, 0x00000200, 0x00000400, 0x00000800,
175 0x00001000, 0x00002000, 0x00004000, 0x00008000,
176 0x00010000, 0x00020000, 0x00040000, 0x00080000,
177 0x00100000, 0x00200000, 0x00400000, 0x00800000,
178 0x01000000, 0x02000000, 0x04000000, 0x08000000,
179 0x10000000, 0x20000000, 0x40000000, 0x80000000,
180 0x00000001
181 #elif BASEB == 16
182 0x0001, 0x0002, 0x0004, 0x0008,
183 0x0010, 0x0020, 0x0040, 0x0080,
184 0x0100, 0x0200, 0x0400, 0x0800,
185 0x1000, 0x2000, 0x4000, 0x8000,
186 0x0001, 0x0002, 0x0004, 0x0008,
187 0x0010, 0x0020, 0x0040, 0x0080,
188 0x0100, 0x0200, 0x0400, 0x0800,
189 0x1000, 0x2000, 0x4000, 0x8000,
190 0x0001
191 #else
192 -=@=- BASEB not 16 or 32 -=@=-
193 #endif
194 }; /* actual rotation thru 8 cycles */
195
196 BOOL _math_abort_; /* nonzero to abort calculations */
197
198
199 /*
200 * popcnt - popcnt[x] number of 1 bits in 0 <= x < 256
201 */
202 char popcnt[256] = {
203 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
204 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
205 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
206 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
207 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
208 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
209 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
210 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
211 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
212 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
213 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
214 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
215 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
216 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
217 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
218 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
219 };
220
221
222
223 #ifdef ALLOCTEST
224 STATIC long nalloc = 0;
225 STATIC long nfree = 0;
226 #endif
227
228
229 HALF *
alloc(LEN len)230 alloc(LEN len)
231 {
232 HALF *hp;
233
234 if (_math_abort_) {
235 math_error("Calculation aborted");
236 /*NOTREACHED*/
237 }
238 hp = (HALF *) malloc((len+1) * sizeof(HALF));
239 if (hp == 0) {
240 math_error("Not enough memory");
241 /*NOTREACHED*/
242 }
243 #ifdef ALLOCTEST
244 ++nalloc;
245 #endif
246 return hp;
247 }
248
249
250 /*
251 * is_const - determine if a HALF array is an pre-allocated array
252 *
253 * given:
254 * h pointer to the beginning of the HALF array
255 *
256 * returns:
257 * TRUE - h is found in the half_tbl array
258 * FALSE - is is not found in the half_tbl array
259 */
260 int
is_const(HALF * h)261 is_const(HALF* h)
262 {
263 HALF **h_p; /* half_tbl array pointer */
264
265 /* search the half_tbl for h */
266 for (h_p = &half_tbl[0]; *h_p != NULL; ++h_p) {
267 if (h == *h_p) {
268 return TRUE; /* found in the half_tbl array */
269 }
270 }
271
272 /* not found in the half_tbl array */
273 return FALSE;
274 }
275
276
277 #ifdef ALLOCTEST
278 void
freeh(HALF * h)279 freeh(HALF *h)
280 {
281 if (!is_const(h)) {
282 free(h);
283 ++nfree;
284 }
285 }
286
287
288 void
allocStat(void)289 allocStat(void)
290 {
291 fprintf(stderr, "nalloc: %ld nfree: %ld kept: %ld\n",
292 nalloc, nfree, nalloc - nfree);
293 }
294 #endif
295
296
297 /*
298 * Convert a normal integer to a number.
299 */
300 void
itoz(long i,ZVALUE * res)301 itoz(long i, ZVALUE *res)
302 {
303 long diddle, len;
304
305 res->len = 1;
306 res->sign = 0;
307 diddle = 0;
308 if (i == 0) {
309 res->v = _zeroval_;
310 return;
311 }
312 if (i < 0) {
313 res->sign = 1;
314 i = -i;
315 if (i < 0) { /* fix most negative number */
316 diddle = 1;
317 i--;
318 }
319 }
320 if (i == 1) {
321 res->v = _oneval_;
322 return;
323 }
324 len = 1 + (((FULL) i) >= BASE);
325 res->len = (LEN)len;
326 res->v = alloc((LEN)len);
327 res->v[0] = (HALF) (i + diddle);
328 if (len == 2)
329 res->v[1] = (HALF) (i / BASE);
330 }
331
332
333 /*
334 * Convert a number to a normal integer, as far as possible.
335 * If the number is out of range, the largest number is returned.
336 */
337 long
ztoi(ZVALUE z)338 ztoi(ZVALUE z)
339 {
340 long i;
341
342 if (zgtmaxlong(z)) {
343 i = MAXLONG;
344 return (z.sign ? -i : i);
345 }
346 i = ztolong(z);
347 return (z.sign ? -i : i);
348 }
349
350
351 /*
352 * Convert a normal unsigned integer to a number.
353 */
354 void
utoz(FULL i,ZVALUE * res)355 utoz(FULL i, ZVALUE *res)
356 {
357 long len;
358
359 res->len = 1;
360 res->sign = 0;
361 if (i == 0) {
362 res->v = _zeroval_;
363 return;
364 }
365 if (i == 1) {
366 res->v = _oneval_;
367 return;
368 }
369 len = 1 + (((FULL) i) >= BASE);
370 res->len = (LEN)len;
371 res->v = alloc((LEN)len);
372 res->v[0] = (HALF)i;
373 if (len == 2)
374 res->v[1] = (HALF) (i / BASE);
375 }
376
377
378 /*
379 * Convert a normal signed integer to a number.
380 */
381 void
stoz(SFULL i,ZVALUE * res)382 stoz(SFULL i, ZVALUE *res)
383 {
384 long diddle, len;
385
386 res->len = 1;
387 res->sign = 0;
388 diddle = 0;
389 if (i == 0) {
390 res->v = _zeroval_;
391 return;
392 }
393 if (i < 0) {
394 res->sign = 1;
395 i = -i;
396 if (i < 0) { /* fix most negative number */
397 diddle = 1;
398 i--;
399 }
400 }
401 if (i == 1) {
402 res->v = _oneval_;
403 return;
404 }
405 len = 1 + (((FULL) i) >= BASE);
406 res->len = (LEN)len;
407 res->v = alloc((LEN)len);
408 res->v[0] = (HALF) (i + diddle);
409 if (len == 2)
410 res->v[1] = (HALF) (i / BASE);
411 }
412
413
414 /*
415 * Convert a number to a unsigned normal integer, as far as possible.
416 * If the number is out of range, the largest number is returned.
417 * The absolute value of z is converted.
418 */
419 FULL
ztou(ZVALUE z)420 ztou(ZVALUE z)
421 {
422 if (z.len > 2) {
423 return MAXUFULL;
424 }
425 return ztofull(z);
426 }
427
428
429 /*
430 * Convert a number to a signed normal integer, as far as possible.
431 *
432 * If the number is too large to fit into an integer, than the largest
433 * positive or largest negative integer is returned depending on the sign.
434 */
435 SFULL
ztos(ZVALUE z)436 ztos(ZVALUE z)
437 {
438 FULL val; /* absolute value of the return value */
439
440 /* negative value processing */
441 if (z.sign) {
442 if (z.len > 2) {
443 return MINSFULL;
444 }
445 val = ztofull(z);
446 if (val > TOPFULL) {
447 return MINSFULL;
448 }
449 return (SFULL)(-val);
450 }
451
452 /* positive value processing */
453 if (z.len > 2) {
454 return (SFULL)MAXFULL;
455 }
456 val = ztofull(z);
457 if (val > MAXFULL) {
458 return (SFULL)MAXFULL;
459 }
460 return (SFULL)val;
461 }
462
463
464 /*
465 * Make a copy of an integer value
466 */
467 void
zcopy(ZVALUE z,ZVALUE * res)468 zcopy(ZVALUE z, ZVALUE *res)
469 {
470 res->sign = z.sign;
471 res->len = z.len;
472 if (zisabsleone(z)) { /* zero or plus or minus one are easy */
473 res->v = (z.v[0] ? _oneval_ : _zeroval_);
474 return;
475 }
476 res->v = alloc(z.len);
477 zcopyval(z, *res);
478 }
479
480
481 /*
482 * Add together two integers.
483 */
484 void
zadd(ZVALUE z1,ZVALUE z2,ZVALUE * res)485 zadd(ZVALUE z1, ZVALUE z2, ZVALUE *res)
486 {
487 ZVALUE dest;
488 HALF *p1, *p2, *pd;
489 long len;
490 FULL carry;
491 SIUNION sival;
492
493 if (z1.sign && !z2.sign) {
494 z1.sign = 0;
495 zsub(z2, z1, res);
496 return;
497 }
498 if (z2.sign && !z1.sign) {
499 z2.sign = 0;
500 zsub(z1, z2, res);
501 return;
502 }
503 if (z2.len > z1.len) {
504 pd = z1.v; z1.v = z2.v; z2.v = pd;
505 len = z1.len; z1.len = z2.len; z2.len = (LEN)len;
506 }
507 dest.len = z1.len + 1;
508 dest.v = alloc(dest.len);
509 dest.sign = z1.sign;
510 carry = 0;
511 pd = dest.v;
512 p1 = z1.v;
513 p2 = z2.v;
514 len = z2.len;
515 while (len--) {
516 sival.ivalue = ((FULL) *p1++) + ((FULL) *p2++) + carry;
517 /* ignore Saber-C warning #112 - get ushort from uint */
518 /* OK to ignore on name zadd`sival */
519 *pd++ = sival.silow;
520 carry = sival.sihigh;
521 }
522 len = z1.len - z2.len;
523 while (len--) {
524 sival.ivalue = ((FULL) *p1++) + carry;
525 *pd++ = sival.silow;
526 carry = sival.sihigh;
527 }
528 *pd = (HALF)carry;
529 zquicktrim(dest);
530 *res = dest;
531 }
532
533
534 /*
535 * Subtract two integers.
536 */
537 void
zsub(ZVALUE z1,ZVALUE z2,ZVALUE * res)538 zsub(ZVALUE z1, ZVALUE z2, ZVALUE *res)
539 {
540 register HALF *h1, *h2, *hd;
541 long len1, len2;
542 FULL carry;
543 SIUNION sival;
544 ZVALUE dest;
545
546 if (z1.sign != z2.sign) {
547 z2.sign = z1.sign;
548 zadd(z1, z2, res);
549 return;
550 }
551 len1 = z1.len;
552 len2 = z2.len;
553 if (len1 == len2) {
554 h1 = z1.v + len1;
555 h2 = z2.v + len2;
556 while ((len1 > 0) && ((FULL)*--h1 == (FULL)*--h2)) {
557 len1--;
558 }
559 if (len1 == 0) {
560 *res = _zero_;
561 return;
562 }
563 len2 = len1;
564 carry = ((FULL)*h1 < (FULL)*h2);
565 } else {
566 carry = (len1 < len2);
567 }
568 dest.sign = z1.sign;
569 h1 = z1.v;
570 h2 = z2.v;
571 if (carry) {
572 carry = len1;
573 len1 = len2;
574 len2 = (long)carry;
575 h1 = z2.v;
576 h2 = z1.v;
577 dest.sign = !dest.sign;
578 }
579 hd = alloc((LEN)len1);
580 dest.v = hd;
581 dest.len = (LEN)len1;
582 len1 -= len2;
583 carry = 0;
584 while (--len2 >= 0) {
585 /* ignore Saber-C warning #112 - get ushort from uint */
586 /* OK to ignore on name zsub`sival */
587 sival.ivalue = (BASE1 - ((FULL) *h1++)) + *h2++ + carry;
588 *hd++ = (HALF)(BASE1 - sival.silow);
589 carry = sival.sihigh;
590 }
591 while (--len1 >= 0) {
592 sival.ivalue = (BASE1 - ((FULL) *h1++)) + carry;
593 *hd++ = (HALF)(BASE1 - sival.silow);
594 carry = sival.sihigh;
595 }
596 if (hd[-1] == 0)
597 ztrim(&dest);
598 *res = dest;
599 }
600
601
602 /*
603 * Multiply an integer by a small number.
604 */
605 void
zmuli(ZVALUE z,long n,ZVALUE * res)606 zmuli(ZVALUE z, long n, ZVALUE *res)
607 {
608 register HALF *h1, *sd;
609 FULL low;
610 FULL high;
611 FULL carry;
612 long len;
613 SIUNION sival;
614 ZVALUE dest;
615
616 if ((n == 0) || ziszero(z)) {
617 *res = _zero_;
618 return;
619 }
620 if (n < 0) {
621 n = -n;
622 z.sign = !z.sign;
623 }
624 if (n == 1) {
625 zcopy(z, res);
626 return;
627 }
628 #if LONG_BITS > BASEB
629 low = ((FULL) n) & BASE1;
630 high = ((FULL) n) >> BASEB;
631 #else
632 low = (FULL)n;
633 high = 0;
634 #endif
635 dest.len = z.len + 2;
636 dest.v = alloc(dest.len);
637 dest.sign = z.sign;
638 /*
639 * Multiply by the low digit.
640 */
641 h1 = z.v;
642 sd = dest.v;
643 len = z.len;
644 carry = 0;
645 while (len--) {
646 /* ignore Saber-C warning #112 - get ushort from uint */
647 /* OK to ignore on name zmuli`sival */
648 sival.ivalue = ((FULL) *h1++) * low + carry;
649 *sd++ = sival.silow;
650 carry = sival.sihigh;
651 }
652 *sd = (HALF)carry;
653 /*
654 * If there was only one digit, then we are all done except
655 * for trimming the number if there was no last carry.
656 */
657 if (high == 0) {
658 dest.len--;
659 if (carry == 0)
660 dest.len--;
661 *res = dest;
662 return;
663 }
664 /*
665 * Need to multiply by the high digit and add it into the
666 * previous value. Clear the final word of rubbish first.
667 */
668 *(++sd) = 0;
669 h1 = z.v;
670 sd = dest.v + 1;
671 len = z.len;
672 carry = 0;
673 while (len--) {
674 sival.ivalue = ((FULL) *h1++) * high + ((FULL) *sd) + carry;
675 *sd++ = sival.silow;
676 carry = sival.sihigh;
677 }
678 *sd = (HALF)carry;
679 zquicktrim(dest);
680 *res = dest;
681 }
682
683
684 /*
685 * Divide two numbers by their greatest common divisor.
686 * This is useful for reducing the numerator and denominator of
687 * a fraction to its lowest terms.
688 */
689 void
zreduce(ZVALUE z1,ZVALUE z2,ZVALUE * z1res,ZVALUE * z2res)690 zreduce(ZVALUE z1, ZVALUE z2, ZVALUE *z1res, ZVALUE *z2res)
691 {
692 ZVALUE tmp;
693
694 if (zisabsleone(z1) || zisabsleone(z2))
695 tmp = _one_;
696 else
697 zgcd(z1, z2, &tmp);
698 if (zisunit(tmp)) {
699 zcopy(z1, z1res);
700 zcopy(z2, z2res);
701 } else {
702 zequo(z1, tmp, z1res);
703 zequo(z2, tmp, z2res);
704 }
705 zfree(tmp);
706 }
707
708
709
710 /*
711 * Compute the quotient and remainder for division of an integer by an
712 * integer, rounding criteria determined by rnd. Returns the sign of
713 * the remainder.
714 */
715 long
zdiv(ZVALUE z1,ZVALUE z2,ZVALUE * quo,ZVALUE * rem,long rnd)716 zdiv(ZVALUE z1, ZVALUE z2, ZVALUE *quo, ZVALUE *rem, long rnd)
717 {
718 register HALF *a, *b;
719 HALF s, u;
720 HALF *A, *B, *a1, *b0;
721 FULL f, g, h, x;
722 BOOL adjust, onebit;
723 LEN m, n, len, i, p, j1, j2, k;
724 long t, val;
725
726 if (ziszero(z2)) {
727 math_error("Division by zero in zdiv");
728 /*NOTREACHED*/
729 }
730 m = z1.len;
731 n = z2.len;
732 B = z2.v;
733 s = 0;
734 if (m < n) {
735 A = alloc(n + 1);
736 memcpy(A, z1.v, m * sizeof(HALF));
737 for (i = m; i <= n; i++)
738 A[i] = 0;
739 a1 = A + n;
740 len = 1;
741 goto done;
742 }
743 A = alloc(m + 2);
744 memcpy(A, z1.v, m * sizeof(HALF));
745 A[m] = 0;
746 A[m + 1] = 0;
747 len = m - n + 1; /* quotient length will be len or len +/- 1 */
748 a1 = A + n; /* start of digits for quotient */
749 b0 = B;
750 p = n;
751 while (!*b0) { /* b0: working start for divisor */
752 ++b0;
753 --p;
754 }
755 if (p == 1) {
756 u = *b0;
757 if (u == 1) {
758 for (; m >= n; m--)
759 A[m] = A[m - 1];
760 A[m] = 0;
761 goto done;
762 }
763 f = 0;
764 a = A + m;
765 i = len;
766 while (i--) {
767 f = f << BASEB | *--a;
768 a[1] = (HALF)(f / u);
769 f = f % u;
770 }
771 *a = (HALF)f;
772 m = n;
773 goto done;
774 }
775 f = B[n - 1];
776 k = 1;
777 while (f >>= 1)
778 k++; /* k: number of bits in top divisor digit */
779 j1 = BASEB - k;
780 j2 = BASEB + j1;
781 h = j1 ? ((FULL) B[n - 1] << j1 | B[n - 2] >> k) : B[n-1];
782 onebit = (BOOL)((B[n - 2] >> (k - 1)) & 1);
783 m++;
784 while (m > n) {
785 m--;
786 f = (FULL) A[m] << j2 | (FULL) A[m - 1] << j1;
787 if (j1) f |= A[m - 2] >> k;
788 if (s) f = ~f;
789 x = f / h;
790 if (x) {
791 if (onebit && x > TOPHALF + f % h)
792 x--;
793 a = A + m - p;
794 b = b0;
795 u = 0;
796 i = p;
797 if (s) {
798 while (i--) {
799 f = (FULL) *a + u + x * *b++;
800 *a++ = (HALF) f;
801 u = (HALF) (f >> BASEB);
802 }
803 s = *a + u;
804 A[m] = (HALF) (~x + !s);
805 } else {
806 while (i--) {
807 f = (FULL) *a - u - x * *b++;
808 *a++ = (HALF) f;
809 u = -(HALF)(f >> BASEB);
810 }
811 s = *a - u;
812 A[m] = (HALF)(x + s);
813 }
814 }
815 else
816 A[m] = s;
817 }
818 done: while (m > 0 && A[m - 1] == 0)
819 m--;
820 if (m == 0 && s == 0) {
821 *rem = _zero_;
822 val = 0;
823 if (a1[len - 1] == 0)
824 len--;
825 if (len == 0) {
826 *quo = _zero_;
827 } else {
828 quo->len = len;
829 quo->v = alloc(len);
830 memcpy(quo->v, a1, len * sizeof(HALF));
831 quo->sign = z1.sign ^ z2.sign;
832 }
833 freeh(A);
834 return val;
835 }
836 if (rnd & 8)
837 adjust = (((*a1 ^ rnd) & 1) ? TRUE : FALSE);
838 else
839 adjust = (((rnd & 1) ^ z1.sign ^ z2.sign) ? TRUE : FALSE);
840 if (rnd & 2)
841 adjust ^= z1.sign ^ z2.sign;
842 if (rnd & 4)
843 adjust ^= z2.sign;
844 if (rnd & 16) { /* round-off case */
845 a = A + n;
846 b = B + n;
847 i = n + 1;
848 f = g = 0;
849 t = -1;
850 if (s) {
851 while (--i > 0) {
852 g = (FULL) *--a + (*--b >> 1 | f);
853 f = *b & 1 ? TOPHALF : 0;
854 if (g != BASE1)
855 break;
856 }
857 if (g == BASE && f == 0) {
858 while ((--i > 0) && ((*--a | *--b) == 0));
859 t = (i > 0);
860 } else if (g >= BASE) {
861 t = 1;
862 }
863 } else {
864 while (--i > 0) {
865 g = (FULL) *--a - (*--b >> 1 | f);
866 f = *b & 1 ? TOPHALF : 0;
867 if (g != 0)
868 break;
869 }
870 if (g > 0 && g < BASE)
871 t = 1;
872 else if (g == 0 && f == 0)
873 t = 0;
874 }
875 if (t)
876 adjust = (t > 0);
877 }
878 if (adjust) {
879 i = len;
880 a = a1;
881 while (i > 0 && *a == BASE1) {
882 i--;
883 *a++ = 0;
884 }
885 (*a)++;
886 if (i == 0)
887 len++;
888 }
889 if (s && adjust) {
890 i = 0;
891 while (A[i] == 0)
892 i++;
893 A[i] = -A[i];
894 while (++i < n)
895 A[i] = ~A[i];
896 m = n;
897 while (A[m - 1] == 0)
898 m--;
899 }
900 if (!s && adjust) {
901 a = A;
902 b = B;
903 i = n;
904 u = 0;
905 while (i--) {
906 f = (FULL) *b++ - *a - (FULL) u;
907 *a++ = (HALF) f;
908 u = -(HALF)(f >> BASEB);
909 }
910 m = n;
911 while (A[m - 1] == 0)
912 m--;
913 }
914 if (s && !adjust) {
915 a = A;
916 b = B;
917 i = n;
918 f = 0;
919 while (i--) {
920 f = (FULL) *b++ + *a + (f >> BASEB);
921 *a++ = (HALF) f;
922 }
923 m = n;
924 while (A[m-1] == 0)
925 m--;
926 }
927 rem->len = m;
928 rem->v = alloc(m);
929 memcpy(rem->v, A, m * sizeof(HALF));
930 rem->sign = z1.sign ^ adjust;
931 val = rem->sign ? -1 : 1;
932 if (a1[len - 1] == 0)
933 len--;
934 if (len == 0) {
935 *quo = _zero_;
936 } else {
937 quo->len = len;
938 quo->v = alloc(len);
939 memcpy(quo->v, a1, len * sizeof(HALF));
940 quo->sign = z1.sign ^ z2.sign;
941 }
942 freeh(A);
943 return val;
944 }
945
946
947 /*
948 * Compute and store at a specified location the integer quotient
949 * of two integers, the type of rounding being determined by rnd.
950 * Returns the sign of z1/z2 - calculated quotient.
951 */
952 long
zquo(ZVALUE z1,ZVALUE z2,ZVALUE * res,long rnd)953 zquo(ZVALUE z1, ZVALUE z2, ZVALUE *res, long rnd)
954 {
955 ZVALUE tmp;
956 long val;
957
958 val = zdiv(z1, z2, res, &tmp, rnd);
959 if (z2.sign)
960 val = -val;
961 zfree(tmp);
962 return val;
963 }
964
965
966 /*
967 * Compute and store at a specified location the remainder for
968 * division of an integer by an integer, the type of rounding
969 * used being determined by rnd. Returns the sign of the remainder.
970 */
971 long
zmod(ZVALUE z1,ZVALUE z2,ZVALUE * res,long rnd)972 zmod(ZVALUE z1, ZVALUE z2, ZVALUE *res, long rnd)
973 {
974 ZVALUE tmp;
975 long val;
976
977 val = zdiv(z1, z2, &tmp, res, rnd);
978 zfree(tmp);
979 return val;
980 }
981
982
983 /*
984 * Computes the quotient z1/z2 on the assumption that this is exact.
985 * There is no thorough check on the exactness of the division
986 * so this should not be called unless z1/z2 is an integer.
987 */
988 void
zequo(ZVALUE z1,ZVALUE z2,ZVALUE * res)989 zequo(ZVALUE z1, ZVALUE z2, ZVALUE *res)
990 {
991 LEN i, j, k, len, m, n, o, p;
992 HALF *a, *a0, *A, *b, *B, u, v, w, x;
993 FULL f, g;
994
995 if (ziszero(z1)) {
996 *res = _zero_;
997 return;
998 }
999 if (ziszero(z2)) {
1000 math_error("Division by zero");
1001 /*NOTREACHED*/
1002 }
1003 if (zisone(z2)) {
1004 zcopy(z1, res);
1005 return;
1006 }
1007 if (zhighbit(z1) < zhighbit(z2)) {
1008 math_error("Bad call to zequo");
1009 /*NOTREACHED*/
1010 }
1011 B = z2.v;
1012 o = 0;
1013 while (!*B) {
1014 ++B;
1015 ++o;
1016 }
1017 m = z1.len - o;
1018 n = z2.len - o;
1019 len = m - n + 1; /* Maximum length of quotient */
1020 v = *B;
1021 A = alloc(len+1);
1022 memcpy(A, z1.v + o, len * sizeof(HALF));
1023 A[len] = 0;
1024 if (n == 1) {
1025 if (v > 1) {
1026 a = A + len;
1027 i = len;
1028 f = 0;
1029 while (i--) {
1030 f = f << BASEB | *--a;
1031 *a = (HALF)(f / v);
1032 f %= v;
1033 }
1034 }
1035 } else {
1036 k = 0;
1037 while (!(v & 1)) {
1038 k++;
1039 v >>= 1;
1040 }
1041 j = BASEB - k;
1042 if (n > 1 && k > 0)
1043 v |= B[1] << j;
1044 u = v - 1;
1045 w = x = 1;
1046 while (u) { /* To find w = inverse of v modulo BASE */
1047 do {
1048 v <<= 1;
1049 x <<= 1;
1050 }
1051 while (!(u & x));
1052 u += v;
1053 w |= x;
1054 }
1055 a0 = A;
1056 p = len;
1057 while (p > 1) {
1058 if (!*a0) {
1059 while (!*++a0 && p > 1)
1060 p--;
1061 --a0;
1062 }
1063 if (p == 1)
1064 break;
1065 x = k ? w * (*a0 >> k | a0[1] << j) : w * *a0;
1066 g = x;
1067 if (x) {
1068 a = a0;
1069 b = B;
1070 u = 0;
1071 i = n;
1072 if (i > p)
1073 i = p;
1074 while (i--) {
1075 f = (FULL) *a - g * *b++ - (FULL) u;
1076 *a++ = (HALF)f;
1077 u = -(HALF)(f >> BASEB);
1078 }
1079 if (u && p > n) {
1080 i = p - n;
1081 while (u && i--) {
1082 f = (FULL) *a - u;
1083 *a++ = (HALF) f;
1084 u = -(HALF)(f >> BASEB);
1085 }
1086 }
1087 }
1088 *a0++ = x;
1089 p--;
1090 }
1091 if (k == 0) {
1092 *a0 = w * *a0;
1093 } else {
1094 u = (HALF)(w * *a0) >> k;
1095 x = (HALF)(((FULL) z1.v[z1.len - 1] << BASEB
1096 | z1.v[z1.len - 2])
1097 /((FULL) B[n-1] << BASEB | B[n-2]));
1098 if ((x ^ u) & 1) x--;
1099 *a0 = x;
1100 }
1101 }
1102 if (!A[len - 1]) len--;
1103 res->v = A;
1104 res->len = len;
1105 res->sign = z1.sign != z2.sign;
1106 }
1107
1108
1109
1110 /*
1111 * Return the quotient and remainder of an integer divided by a small
1112 * number. A nonzero remainder is only meaningful when both numbers
1113 * are positive.
1114 */
1115 long
zdivi(ZVALUE z,long n,ZVALUE * res)1116 zdivi(ZVALUE z, long n, ZVALUE *res)
1117 {
1118 HALF *h1, *sd;
1119 FULL val;
1120 HALF divval[2];
1121 ZVALUE div;
1122 ZVALUE dest;
1123 LEN len;
1124
1125 if (n == 0) {
1126 math_error("Division by zero");
1127 /*NOTREACHED*/
1128 }
1129 if (ziszero(z)) {
1130 *res = _zero_;
1131 return 0;
1132 }
1133 if (n < 0) {
1134 n = -n;
1135 z.sign = !z.sign;
1136 }
1137 if (n == 1) {
1138 zcopy(z, res);
1139 return 0;
1140 }
1141 /*
1142 * If the division is by a large number, then call the normal
1143 * divide routine.
1144 */
1145 if (n & ~BASE1) {
1146 div.sign = 0;
1147 div.v = divval;
1148 divval[0] = (HALF) n;
1149 #if LONG_BITS > BASEB
1150 divval[1] = (HALF)(((FULL) n) >> BASEB);
1151 div.len = 2;
1152 #else
1153 div.len = 1;
1154 #endif
1155 zdiv(z, div, res, &dest, 0);
1156 n = ztolong(dest);
1157 zfree(dest);
1158 return n;
1159 }
1160 /*
1161 * Division is by a small number, so we can be quick about it.
1162 */
1163 len = z.len;
1164 dest.sign = z.sign;
1165 dest.len = len;
1166 dest.v = alloc(len);
1167 h1 = z.v + len;
1168 sd = dest.v + len;
1169 val = 0;
1170 while (len--) {
1171 val = ((val << BASEB) + ((FULL) *--h1));
1172 *--sd = (HALF)(val / n);
1173 val %= n;
1174 }
1175 zquicktrim(dest);
1176 *res = dest;
1177 return (long)val;
1178 }
1179
1180
1181
1182 /*
1183 * Calculate the mod of an integer by a small number.
1184 * This is only defined for positive moduli.
1185 */
1186 long
zmodi(ZVALUE z,long n)1187 zmodi(ZVALUE z, long n)
1188 {
1189 register HALF *h1;
1190 FULL val;
1191 HALF divval[2];
1192 ZVALUE div;
1193 ZVALUE temp;
1194 long len;
1195
1196 if (n == 0) {
1197 math_error("Division by zero");
1198 /*NOTREACHED*/
1199 }
1200 if (n < 0) {
1201 math_error("Non-positive modulus");
1202 /*NOTREACHED*/
1203 }
1204 if (ziszero(z) || (n == 1))
1205 return 0;
1206 if (zisone(z))
1207 return 1;
1208 /*
1209 * If the modulus is by a large number, then call the normal
1210 * modulo routine.
1211 */
1212 if (n & ~BASE1) {
1213 div.sign = 0;
1214 div.v = divval;
1215 divval[0] = (HALF) n;
1216 #if LONG_BITS > BASEB
1217 divval[1] = (HALF)(((FULL) n) >> BASEB);
1218 div.len = 2;
1219 #else
1220 div.len = 1;
1221 #endif
1222 zmod(z, div, &temp, 0);
1223 n = ztolong(temp);
1224 zfree(temp);
1225 return n;
1226 }
1227 /*
1228 * The modulus is by a small number, so we can do this quickly.
1229 */
1230 len = z.len;
1231 h1 = z.v + len;
1232 val = 0;
1233 while (len-- > 0)
1234 val = ((val << BASEB) + ((FULL) *--h1)) % n;
1235 if (val && z.sign)
1236 val = n - val;
1237 return (long)val;
1238 }
1239
1240
1241 /*
1242 * Return whether or not one number exactly divides another one.
1243 * Returns TRUE if division occurs with no remainder.
1244 * z1 is the number to be divided by z2.
1245 */
1246
1247 BOOL
zdivides(ZVALUE z1,ZVALUE z2)1248 zdivides(ZVALUE z1, ZVALUE z2)
1249 {
1250 LEN i, j, k, m, n;
1251 HALF u, v, w, x;
1252 HALF *a, *a0, *A, *b, *B, *c, *d;
1253 FULL f;
1254 BOOL ans;
1255
1256 if (zisunit(z2) || ziszero(z1)) return TRUE;
1257 if (ziszero(z2)) return FALSE;
1258
1259 m = z1.len;
1260 n = z2.len;
1261 if (m < n) return FALSE;
1262
1263 c = z1.v;
1264 d = z2.v;
1265
1266 while (!*d) {
1267 if (*c++) return FALSE;
1268 d++;
1269 m--;
1270 n--;
1271 }
1272
1273 j = 0;
1274 u = *c;
1275 v = *d;
1276 while (!(v & 1)) { /* Counting and checking zero bits */
1277 if (u & 1) return FALSE;
1278 u >>= 1;
1279 v >>= 1;
1280 j++;
1281 }
1282
1283 if (n == 1 && v == 1) return TRUE;
1284 if (!*c) { /* Removing any further zeros */
1285 while(!*++c)
1286 m--;
1287 c--;
1288 }
1289
1290 if (m < n) return FALSE;
1291
1292 if (j) {
1293 B = alloc((LEN)n); /* Array for shifted z2 */
1294 d += n;
1295 b = B + n;
1296 i = n;
1297 f = 0;
1298 while(i--) {
1299 f = f << BASEB | *--d;
1300 *--b = (HALF)(f >> j);
1301 }
1302 if (!B[n - 1]) n--;
1303 }
1304 else B = d;
1305 u = *B;
1306 v = x = 1;
1307 w = 0;
1308 while (x) { /* Finding minv(*B, BASE) */
1309 if (v & x) {
1310 v -= x * u;
1311 w |= x;
1312 }
1313 x <<= 1;
1314 }
1315
1316 A = alloc((LEN)(m + 2)); /* Main working array */
1317 memcpy(A, c, m * sizeof(HALF));
1318 A[m + 1] = 1;
1319 A[m] = 0;
1320
1321 k = m - n + 1; /* Length of presumed quotient */
1322
1323 a0 = A;
1324
1325 while (k--) {
1326 if (*a0) {
1327 x = w * *a0;
1328 a = a0;
1329 b = B;
1330 i = n;
1331 u = 0;
1332 while (i--) {
1333 f = (FULL) *a - (FULL) x * *b++ - u;
1334 *a++ = (HALF)f;
1335 u = -(HALF)(f >> BASEB);
1336 }
1337 f = (FULL) *a - u;
1338 *a++ = (HALF)f;
1339 u = -(HALF)(f >> BASEB);
1340 if (u) {
1341 while (*a == 0) *a++ = BASE1;
1342 (*a)--;
1343 }
1344 }
1345 a0++;
1346 }
1347 ans = FALSE;
1348 if (A[m + 1]) {
1349 a = A + m;
1350 while (--n && !*--a);
1351 if (!n) ans = TRUE;
1352 }
1353 freeh(A);
1354 if (j) freeh(B);
1355 return ans;
1356 }
1357
1358
1359 /*
1360 * Compute the bitwise OR of two integers
1361 */
1362 void
zor(ZVALUE z1,ZVALUE z2,ZVALUE * res)1363 zor(ZVALUE z1, ZVALUE z2, ZVALUE *res)
1364 {
1365 register HALF *sp, *dp;
1366 long len;
1367 ZVALUE bz, lz, dest;
1368
1369 if (z1.len >= z2.len) {
1370 bz = z1;
1371 lz = z2;
1372 } else {
1373 bz = z2;
1374 lz = z1;
1375 }
1376 dest.len = bz.len;
1377 dest.v = alloc(dest.len);
1378 dest.sign = 0;
1379 zcopyval(bz, dest);
1380 len = lz.len;
1381 sp = lz.v;
1382 dp = dest.v;
1383 while (len--)
1384 *dp++ |= *sp++;
1385 *res = dest;
1386 }
1387
1388
1389 /*
1390 * Compute the bitwise AND of two integers
1391 */
1392 void
zand(ZVALUE z1,ZVALUE z2,ZVALUE * res)1393 zand(ZVALUE z1, ZVALUE z2, ZVALUE *res)
1394 {
1395 HALF *h1, *h2, *hd;
1396 LEN len;
1397 ZVALUE dest;
1398
1399 len = ((z1.len <= z2.len) ? z1.len : z2.len);
1400 h1 = &z1.v[len-1];
1401 h2 = &z2.v[len-1];
1402 while ((len > 1) && ((*h1 & *h2) == 0)) {
1403 h1--;
1404 h2--;
1405 len--;
1406 }
1407 dest.len = len;
1408 dest.v = alloc(len);
1409 dest.sign = 0;
1410 h1 = z1.v;
1411 h2 = z2.v;
1412 hd = dest.v;
1413 while (len--)
1414 *hd++ = (*h1++ & *h2++);
1415 *res = dest;
1416 }
1417
1418
1419 /*
1420 * Compute the bitwise XOR of two integers.
1421 */
1422 void
zxor(ZVALUE z1,ZVALUE z2,ZVALUE * res)1423 zxor(ZVALUE z1, ZVALUE z2, ZVALUE *res)
1424 {
1425 HALF *dp, *h1, *h2;
1426 LEN len, j, k;
1427 ZVALUE dest;
1428
1429 h1 = z1.v;
1430 h2 = z2.v;
1431 len = z1.len;
1432 j = z2.len;
1433 if (z1.len < z2.len) {
1434 len = z2.len;
1435 j = z1.len;
1436 h1 = z2.v;
1437 h2 = z1.v;
1438 } else if (z1.len == z2.len) {
1439 while (len > 1 && z1.v[len-1] == z2.v[len-1])
1440 len--;
1441 j = len;
1442 }
1443 k = len - j;
1444 dest.len = len;
1445 dest.v = alloc(len);
1446 dest.sign = 0;
1447 dp = dest.v;
1448 while (j-- > 0)
1449 *dp++ = *h1++ ^ *h2++;
1450 while (k-- > 0)
1451 *dp++ = *h1++;
1452 *res = dest;
1453 }
1454
1455
1456 /*
1457 * Compute the bitwise AND NOT of two integers.
1458 */
1459 void
zandnot(ZVALUE z1,ZVALUE z2,ZVALUE * res)1460 zandnot(ZVALUE z1, ZVALUE z2, ZVALUE *res)
1461 {
1462 HALF *dp, *h1, *h2;
1463 LEN len, j, k;
1464 ZVALUE dest;
1465
1466 len = z1.len;
1467 if (z2.len >= len) {
1468 while (len > 1 && (z1.v[len-1] & ~z2.v[len-1]) == 0)
1469 len--;
1470 j = len;
1471 k = 0;
1472 } else {
1473 j = z2.len;
1474 k = len - z2.len;
1475 }
1476 dest.len = len;
1477 dest.v = alloc(len);
1478 dest.sign = 0;
1479 dp = dest.v;
1480 h1 = z1.v;
1481 h2 = z2.v;
1482 while (j-- > 0)
1483 *dp++ = *h1++ & ~*h2++;
1484 while (k-- > 0)
1485 *dp++ = *h1++;
1486 *res = dest;
1487 }
1488
1489
1490 /*
1491 * Shift a number left (or right) by the specified number of bits.
1492 * Positive shift means to the left. When shifting right, rightmost
1493 * bits are lost. The sign of the number is preserved.
1494 */
1495 void
zshift(ZVALUE z,long n,ZVALUE * res)1496 zshift(ZVALUE z, long n, ZVALUE *res)
1497 {
1498 ZVALUE ans;
1499 LEN hc; /* number of halfwords shift is by */
1500
1501 if (ziszero(z)) {
1502 *res = _zero_;
1503 return;
1504 }
1505 if (n == 0) {
1506 zcopy(z, res);
1507 return;
1508 }
1509 /*
1510 * If shift value is negative, then shift right.
1511 * Check for large shifts, and handle word-sized shifts quickly.
1512 */
1513 if (n < 0) {
1514 n = -n;
1515 if ((n < 0) || (n >= (z.len * BASEB))) {
1516 *res = _zero_;
1517 return;
1518 }
1519 hc = (LEN)(n / BASEB);
1520 n %= BASEB;
1521 z.v += hc;
1522 z.len -= hc;
1523 ans.len = z.len;
1524 ans.v = alloc(ans.len);
1525 ans.sign = z.sign;
1526 zcopyval(z, ans);
1527 if (n > 0) {
1528 zshiftr(ans, n);
1529 ztrim(&ans);
1530 }
1531 if (ziszero(ans)) {
1532 zfree(ans);
1533 ans = _zero_;
1534 }
1535 *res = ans;
1536 return;
1537 }
1538 /*
1539 * Shift value is positive, so shift leftwards.
1540 * Check specially for a shift of the value 1, since this is common.
1541 * Also handle word-sized shifts quickly.
1542 */
1543 if (zisunit(z)) {
1544 zbitvalue(n, res);
1545 res->sign = z.sign;
1546 return;
1547 }
1548 hc = (LEN)(n / BASEB);
1549 n %= BASEB;
1550 ans.len = z.len + hc + 1;
1551 ans.v = alloc(ans.len);
1552 ans.sign = z.sign;
1553 if (hc > 0)
1554 memset((char *) ans.v, 0, hc * sizeof(HALF));
1555 memcpy((char *) (ans.v + hc),
1556 (char *) z.v, z.len * sizeof(HALF));
1557 ans.v[ans.len - 1] = 0;
1558 if (n > 0) {
1559 ans.v += hc;
1560 ans.len -= hc;
1561 zshiftl(ans, n);
1562 ans.v -= hc;
1563 ans.len += hc;
1564 }
1565 ztrim(&ans);
1566 *res = ans;
1567 }
1568
1569
1570 /*
1571 * Return the position of the lowest bit which is set in the binary
1572 * representation of a number (counting from zero). This is the highest
1573 * power of two which evenly divides the number.
1574 */
1575 long
zlowbit(ZVALUE z)1576 zlowbit(ZVALUE z)
1577 {
1578 register HALF *zp;
1579 long n;
1580 HALF dataval;
1581 HALF *bitval;
1582
1583 n = 0;
1584 for (zp = z.v; *zp == 0; zp++)
1585 if (++n >= z.len)
1586 return 0;
1587 dataval = *zp;
1588 bitval = bitmask;
1589 /* ignore Saber-C warning #530 about empty while statement */
1590 /* OK to ignore in proc zlowbit */
1591 while ((*(bitval++) & dataval) == 0) {
1592 }
1593 return (n*BASEB)+(bitval-bitmask-1);
1594 }
1595
1596
1597 /*
1598 * Return the position of the highest bit which is set in the binary
1599 * representation of a number (counting from zero). This is the highest power
1600 * of two which is less than or equal to the number (which is assumed nonzero).
1601 */
1602 LEN
zhighbit(ZVALUE z)1603 zhighbit(ZVALUE z)
1604 {
1605 HALF dataval;
1606 HALF *bitval;
1607
1608 dataval = z.v[z.len-1];
1609 if (dataval == 0) {
1610 return 0;
1611 }
1612 bitval = bitmask+BASEB;
1613 if (dataval) {
1614 /* ignore Saber-C warning #530 about empty while statement */
1615 /* OK to ignore in proc zhighbit */
1616 while ((*(--bitval) & dataval) == 0) {
1617 }
1618 }
1619 return (z.len*BASEB)+(LEN)(bitval-bitmask-BASEB);
1620 }
1621
1622
1623 /*
1624 * Return whether or not the specified bit number is set in a number.
1625 * Rightmost bit of a number is bit 0.
1626 */
1627 BOOL
zisset(ZVALUE z,long n)1628 zisset(ZVALUE z, long n)
1629 {
1630 if ((n < 0) || ((n / BASEB) >= z.len))
1631 return FALSE;
1632 return ((z.v[n / BASEB] & (((HALF) 1) << (n % BASEB))) != 0);
1633 }
1634
1635
1636 /*
1637 * Check whether or not a number has exactly one bit set, and
1638 * thus is an exact power of two. Returns TRUE if so.
1639 */
1640 BOOL
zisonebit(ZVALUE z)1641 zisonebit(ZVALUE z)
1642 {
1643 register HALF *hp;
1644 register LEN len;
1645
1646 if (ziszero(z) || zisneg(z))
1647 return FALSE;
1648 hp = z.v;
1649 len = z.len;
1650 while (len > 4) {
1651 len -= 4;
1652 if (*hp++ || *hp++ || *hp++ || *hp++)
1653 return FALSE;
1654 }
1655 while (--len > 0) {
1656 if (*hp++)
1657 return FALSE;
1658 }
1659 return ((*hp & -*hp) == *hp); /* NEEDS 2'S COMPLEMENT */
1660 }
1661
1662
1663 /*
1664 * Check whether or not a number has all of its bits set below some
1665 * bit position, and thus is one less than an exact power of two.
1666 * Returns TRUE if so.
1667 */
1668 BOOL
zisallbits(ZVALUE z)1669 zisallbits(ZVALUE z)
1670 {
1671 register HALF *hp;
1672 register LEN len;
1673 HALF digit;
1674
1675 if (ziszero(z) || zisneg(z))
1676 return FALSE;
1677 hp = z.v;
1678 len = z.len;
1679 while (len > 4) {
1680 len -= 4;
1681 if ((*hp++ != BASE1) || (*hp++ != BASE1) ||
1682 (*hp++ != BASE1) || (*hp++ != BASE1))
1683 return FALSE;
1684 }
1685 while (--len > 0) {
1686 if (*hp++ != BASE1)
1687 return FALSE;
1688 }
1689 digit = (HALF)(*hp + 1);
1690 return ((digit & -digit) == digit); /* NEEDS 2'S COMPLEMENT */
1691 }
1692
1693
1694 /*
1695 * Return the number whose binary representation contains only one bit which
1696 * is in the specified position (counting from zero). This is equivalent
1697 * to raising two to the given power.
1698 */
1699 void
zbitvalue(long n,ZVALUE * res)1700 zbitvalue(long n, ZVALUE *res)
1701 {
1702 ZVALUE z;
1703
1704 if (n < 0) n = 0;
1705 z.sign = 0;
1706 z.len = (LEN)((n / BASEB) + 1);
1707 z.v = alloc(z.len);
1708 zclearval(z);
1709 z.v[z.len-1] = (((HALF) 1) << (n % BASEB));
1710 *res = z;
1711 }
1712
1713
1714 /*
1715 * Compare a number against zero.
1716 * Returns the sgn function of the number (-1, 0, or 1).
1717 */
1718 FLAG
ztest(ZVALUE z)1719 ztest(ZVALUE z)
1720 {
1721 register int sign;
1722 register HALF *h;
1723 register long len;
1724
1725 sign = 1;
1726 if (z.sign)
1727 sign = -sign;
1728 h = z.v;
1729 len = z.len;
1730 while (len--)
1731 if (*h++)
1732 return sign;
1733 return 0;
1734 }
1735
1736
1737 /*
1738 * Return the sign of z1 - z2, i.e. 1 if the first integer is greater,
1739 * 0 if they are equal, -1 otherwise.
1740 */
1741 FLAG
zrel(ZVALUE z1,ZVALUE z2)1742 zrel(ZVALUE z1, ZVALUE z2)
1743 {
1744 HALF *h1, *h2;
1745 LEN len;
1746 int sign;
1747
1748 sign = 1;
1749 if (z1.sign < z2.sign)
1750 return 1;
1751 if (z2.sign < z1.sign)
1752 return -1;
1753 if (z2.sign)
1754 sign = -1;
1755 if (z1.len != z2.len)
1756 return (z1.len > z2.len) ? sign : -sign;
1757 len = z1.len;
1758 h1 = z1.v + len;
1759 h2 = z2.v + len;
1760
1761 while (len > 0) {
1762 if (*--h1 != *--h2)
1763 break;
1764 len--;
1765 }
1766 if (len > 0)
1767 return (*h1 > *h2) ? sign : -sign;
1768 return 0;
1769 }
1770
1771
1772 /*
1773 * Return the sign of abs(z1) - abs(z2), i.e. 1 if the first integer
1774 * has greater absolute value, 0 is they have equal absolute value,
1775 * -1 otherwise.
1776 */
1777 FLAG
zabsrel(ZVALUE z1,ZVALUE z2)1778 zabsrel(ZVALUE z1, ZVALUE z2)
1779 {
1780 HALF *h1, *h2;
1781 LEN len;
1782
1783 if (z1.len != z2.len)
1784 return (z1.len > z2.len) ? 1 : -1;
1785
1786 len = z1.len;
1787 h1 = z1.v + len;
1788 h2 = z2.v + len;
1789 while (len > 0) {
1790 if (*--h1 != *--h2)
1791 break;
1792 len--;
1793 }
1794 if (len > 0)
1795 return (*h1 > *h2) ? 1 : -1;
1796 return 0;
1797 }
1798
1799
1800 /*
1801 * Compare two numbers to see if they are equal or not.
1802 * Returns TRUE if they differ.
1803 */
1804 BOOL
zcmp(ZVALUE z1,ZVALUE z2)1805 zcmp(ZVALUE z1, ZVALUE z2)
1806 {
1807 register HALF *h1, *h2;
1808 register long len;
1809
1810 if ((z1.sign != z2.sign) || (z1.len != z2.len) || (*z1.v != *z2.v))
1811 return TRUE;
1812 len = z1.len;
1813 h1 = z1.v;
1814 h2 = z2.v;
1815 while (--len > 0) {
1816 if (*++h1 != *++h2)
1817 return TRUE;
1818 }
1819 return FALSE;
1820 }
1821
1822
1823 /*
1824 * Utility to calculate the gcd of two FULL integers.
1825 */
1826 FULL
uugcd(FULL f1,FULL f2)1827 uugcd(FULL f1, FULL f2)
1828 {
1829 FULL temp;
1830
1831 while (f1) {
1832 temp = f2 % f1;
1833 f2 = f1;
1834 f1 = temp;
1835 }
1836 return (FULL) f2;
1837 }
1838
1839
1840 /*
1841 * Utility to calculate the gcd of two small integers.
1842 */
1843 long
iigcd(long i1,long i2)1844 iigcd(long i1, long i2)
1845 {
1846 FULL f1, f2, temp;
1847
1848 f1 = (FULL) ((i1 >= 0) ? i1 : -i1);
1849 f2 = (FULL) ((i2 >= 0) ? i2 : -i2);
1850 while (f1) {
1851 temp = f2 % f1;
1852 f2 = f1;
1853 f1 = temp;
1854 }
1855 return (long) f2;
1856 }
1857
1858
1859 void
ztrim(ZVALUE * z)1860 ztrim(ZVALUE *z)
1861 {
1862 HALF *h;
1863 LEN len;
1864
1865 h = z->v + z->len - 1;
1866 len = z->len;
1867 while (*h == 0 && len > 1) {
1868 --h;
1869 --len;
1870 }
1871 z->len = len;
1872 }
1873
1874
1875 /*
1876 * Utility routine to shift right.
1877 *
1878 * NOTE: The ZVALUE length is not adjusted instead, the value is
1879 * zero padded from the left. One may need to call ztrim()
1880 * or use zshift() instead.
1881 */
1882 void
zshiftr(ZVALUE z,long n)1883 zshiftr(ZVALUE z, long n)
1884 {
1885 register HALF *h, *lim;
1886 FULL mask, maskt;
1887 long len;
1888
1889 if (n >= BASEB) {
1890 len = n / BASEB;
1891 h = z.v;
1892 lim = z.v + z.len - len;
1893 while (h < lim) {
1894 h[0] = h[len];
1895 ++h;
1896 }
1897 n -= BASEB * len;
1898 lim = z.v + z.len;
1899 while (h < lim)
1900 *h++ = 0;
1901 }
1902 if (n) {
1903 len = z.len;
1904 h = z.v + len;
1905 mask = 0;
1906 while (len--) {
1907 maskt = (((FULL) *--h) << (BASEB - n)) & BASE1;
1908 *h = ((*h >> n) | (HALF)mask);
1909 mask = maskt;
1910 }
1911 }
1912 }
1913
1914
1915 /*
1916 * Utility routine to shift left.
1917 *
1918 * NOTE: The ZVALUE length is not adjusted. The bits in the upper
1919 * HALF are simply tossed. You may want to use zshift() instead.
1920 */
1921 void
zshiftl(ZVALUE z,long n)1922 zshiftl(ZVALUE z, long n)
1923 {
1924 register HALF *h;
1925 FULL mask, i;
1926 long len;
1927
1928 if (n >= BASEB) {
1929 len = n / BASEB;
1930 h = z.v + z.len - 1;
1931 while (!*h)
1932 --h;
1933 while (h >= z.v) {
1934 h[len] = h[0];
1935 --h;
1936 }
1937 n -= BASEB * len;
1938 while (len)
1939 h[len--] = 0;
1940 }
1941 if (n > 0) {
1942 len = z.len;
1943 h = z.v;
1944 mask = 0;
1945 while (len--) {
1946 i = (((FULL) *h) << n) | mask;
1947 if (i > BASE1) {
1948 mask = i >> BASEB;
1949 i &= BASE1;
1950 } else {
1951 mask = 0;
1952 }
1953 *h = (HALF) i;
1954 ++h;
1955 }
1956 }
1957 }
1958
1959
1960 /*
1961 * popcnt - count the number of 0 or 1 bits in an integer
1962 *
1963 * We ignore all 0 bits above the highest bit.
1964 */
1965 long
zpopcnt(ZVALUE z,int bitval)1966 zpopcnt(ZVALUE z, int bitval)
1967 {
1968 long cnt = 0; /* number of times found */
1969 HALF h; /* HALF to count */
1970 int i;
1971
1972 /*
1973 * count 1's
1974 */
1975 if (bitval) {
1976
1977 /*
1978 * count each HALF
1979 */
1980 for (i=0; i < z.len; ++i) {
1981 /* count each octet */
1982 for (h = z.v[i]; h; h >>= 8) {
1983 cnt += (long)popcnt[h & 0xff];
1984 }
1985 }
1986
1987 /*
1988 * count 0's
1989 */
1990 } else {
1991
1992 /*
1993 * count each HALF up until the last
1994 */
1995 for (i=0; i < z.len-1; ++i) {
1996
1997 /* count each octet */
1998 cnt += BASEB;
1999 for (h = z.v[i]; h; h >>= 8) {
2000 cnt -= (long)popcnt[h & 0xff];
2001 }
2002 }
2003
2004 /*
2005 * count the last octet up until the highest 1 bit
2006 */
2007 for (h = z.v[z.len-1]; h; h>>=1) {
2008 /* count each 0 bit */
2009 if ((h & 0x1) == 0) {
2010 ++cnt;
2011 }
2012 }
2013 }
2014
2015 /*
2016 * return count
2017 */
2018 return cnt;
2019 }
2020