c_cpp/etc/calc/zmath.c

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