scout/source/i64.c

   1 /* i64.c */
   2 /*--------------------------------------------------------------------------*\
   3   Copyright (C) 1999 Douglas W. Sauder
   4
   5   This software is provided "as is," without any express or implied
   6   warranty.  In no event will the author be held liable for any damages
   7   arising from the use of this software.
   8
   9   Permission is granted to anyone to use this software for any purpose,
  10   including commercial applications, and to alter it and redistribute it
  11   freely, subject to the following restrictions:
  12
  13   1. The origin of this software must not be misrepresented; you must not
  14      claim that you wrote the original software. If you use this software
  15      in a product, an acknowledgment in the product documentation would be
  16      appreciated but is not required.
  17   2. Altered source versions must be plainly marked as such, and must not be
  18      misrepresented as being the original software.
  19   3. This notice may not be removed or altered from any source distribution.
  20
  21   The original distribution can be obtained from www.hunnysoft.com.
  22   You can email the author, Doug Sauder, at dwsauder@erols.com.
  23
  24   $RCSfile: i64.c,v $
  25   $Revision$
  26   $Date$
  27 \*--------------------------------------------------------------------------*/
  28
  29 #include <ctype.h>
  30 #include <string.h>
  31 #include "i64.h"
  32
  33 static void i64_udiv_impl(bigint dividend, bigint divisor, bigint *quotient,
  34     bigint *remainder);
  35
  36
  37 /*
  38  * Compare two values
  39  *
  40  *    n1 < n2  -->  -1
  41  *    n1 = n2  -->   0
  42  *    n1 > n2  -->   1
  43  */
  44 int i64_cmp(bigint n1, bigint n2)
  45 {
  46     int cmp;
  47
  48     cmp = 1;
  49     if ((int) n1.hi < (int) n2.hi) {
  50         cmp = -1;
  51     }
  52     else if (n1.hi == n2.hi) {
  53         if (n1.lo < n2.lo) {
  54             cmp = -1;
  55         }
  56         else if (n1.lo == n2.lo) {
  57            cmp = 0;
  58         }
  59     }
  60     return cmp;
  61 }
  62
  63
  64 /*
  65  * "Sign" operation
  66  *
  67  *    n1 < 0  -->  -1
  68  *    n1 = 0  -->   0
  69  *    n1 > 0  -->   1
  70  */
  71 int i64_sgn(bigint n1)
  72 {
  73     int sgn;
  74
  75     sgn = 1;
  76     if (n1.hi & 0x80000000) {
  77        sgn = -1;
  78     }
  79     else if ((n1.hi | n1.lo) == 0) {
  80        sgn = 0;
  81     }
  82     return sgn;
  83 }
  84
  85
  86 /*
  87  * Left shift n by b bits.  Undefined if b is negative.
  88  */
  89 bigint i64_lshift(bigint n, int b)
  90 {
  91     if (b > 0) {
  92         if (b < 32) {
  93             n.hi <<= b;
  94             n.hi += n.lo >> (32 - b);
  95             n.lo <<= b;
  96         }
  97         else if (b < 64) {
  98             n.hi = n.lo << (b - 32);
  99             n.lo = 0;
 100         }
 101         else {
 102             n.lo = 0;
 103             n.hi = 0;
 104         }
 105     }
 106     return n;
 107 }
 108
 109
 110 /*
 111  * Unsigned right shift n by b bits.  Zeros are always shifted in.
 112  * Undefined if b is negative.
 113  */
 114 bigint i64_urshift(bigint n, int b)
 115 {
 116     if (b > 0) {
 117         if (b < 32) {
 118             n.lo >>= b;
 119             n.lo += n.hi << (32 - b);
 120             n.hi >>= b;
 121         }
 122         else if (b < 64) {
 123             n.lo = n.hi >> (b - 32);
 124             n.hi = 0;
 125         }
 126         else {
 127             n.lo = 0;
 128             n.hi = 0;
 129         }
 130     }
 131     return n;
 132 }
 133
 134
 135 /*
 136  * Signed right shift n by b bits.  Bit shifted in matches most significant
 137  * bit in n.  Undefined if b is negative.
 138  */
 139 bigint i64_srshift(bigint n, int b)
 140 {
 141     unsigned m;
 142
 143     if (b > 0) {
 144         m = 0;
 145         if (n.hi & 0x80000000) {
 146             m = 0xffffffff;
 147         }
 148         if (b < 32) {
 149             n.lo >>= b;
 150             n.lo += n.hi << (32 - b);
 151             n.hi >>= b;
 152             n.hi += m << (32 - b);
 153         }
 154         else if (b < 64) {
 155             b -= 32;
 156             n.lo = n.hi >> b;
 157             n.lo += m << (32 - b);
 158             n.hi = m;
 159         }
 160         else {
 161             n.lo = m;
 162             n.hi = m;
 163         }
 164     }
 165     return n;
 166 }
 167
 168
 169 /*
 170  * Change sign operation (additive inverse)
 171  *
 172  *   n  -->  -n
 173  */
 174 bigint i64_inv(bigint n)
 175 {
 176     bigint r;
 177
 178     r.lo = (unsigned) - (int) n.lo;
 179     if (r.lo == 0) {
 180         r.hi = (unsigned) - (int) n.hi;
 181     }
 182     else /* if (r.lo != 0) */ {
 183         r.hi = ~n.hi;
 184     }
 185     return r;
 186 }
 187
 188
 189 /*
 190  * Add n1 and n2
 191  */
 192 bigint i64_add(bigint n1, bigint n2)
 193 {
 194     bigint sum;
 195     unsigned c;
 196
 197     sum.lo = n1.lo + n2.lo;
 198     c = (n1.lo & n2.lo) | ((~sum.lo) & (n1.lo | n2.lo));
 199     c >>= 31;
 200     sum.hi = n1.hi + n2.hi + c;
 201     return sum;
 202 }
 203
 204
 205 /*
 206  * Subtract n2 from n1
 207  */
 208 bigint i64_sub(bigint n1, bigint n2)
 209 {
 210     bigint diff;
 211     unsigned c;
 212
 213     diff.lo = n1.lo - n2.lo;
 214     c = (diff.lo & n2.lo) | ((~n1.lo) & (diff.lo | n2.lo));
 215     c >>= 31;
 216     diff.hi = n1.hi - n2.hi - c;
 217     return diff;
 218 }
 219
 220
 221 /*
 222  * Multiply n1 and n2
 223  */
 224 bigint i64_mul(bigint n1, bigint n2)
 225 {
 226     int prodIsMinus;
 227     bigint prod;
 228     unsigned m1_0, m1_1, m1_2, m1_3;
 229     unsigned m2_0, m2_1, m2_2, m2_3;
 230     unsigned ms;
 231
 232     prodIsMinus = 0;
 233     if (n1.hi & 0x80000000) {
 234         n1 = i64_inv(n1);
 235         prodIsMinus = ~prodIsMinus;
 236     }
 237     if (n2.hi & 0x80000000) {
 238         n2 = i64_inv(n2);
 239         prodIsMinus = ~prodIsMinus;
 240     }
 241
 242     m1_3 = n1.hi >> 16;
 243     m1_2 = n1.hi & 0xffff;
 244     m1_1 = n1.lo >> 16;
 245     m1_0 = n1.lo & 0xffff;
 246
 247     m2_3 = n2.hi >> 16;
 248     m2_2 = n2.hi & 0xffff;
 249     m2_1 = n2.lo >> 16;
 250     m2_0 = n2.lo & 0xffff;
 251
 252     ms = m1_0 * m2_0;
 253     prod.lo = ms & 0xffff;
 254     ms >>= 16;
 255
 256     ms += m1_0 * m2_1;
 257     prod.hi = ms >> 16;
 258     ms &= 0xffff;
 259
 260     ms += m1_1 * m2_0;
 261     prod.hi += ms >> 16;
 262     prod.lo += ms << 16;
 263
 264     prod.hi += (m1_0 * m2_2) + (m1_1 * m2_1) + (m1_2 * m2_0);
 265
 266     ms = (m1_0 * m2_3) + (m1_1 * m2_2) + (m1_2 * m2_1) + (m1_3 * m2_0);
 267     prod.hi += ms << 16;
 268
 269     if (prodIsMinus) {
 270         prod = i64_inv(prod);
 271     }
 272
 273     return prod;
 274 }
 275
 276
 277 static char leftZerosInByte[] = {
 278 /* 00  01  02  03  04  05  06  07  08  09  0a  0b  0c  0d  0e  0f */
 279     8,  7,  6,  6,  5,  5,  5,  5,  4,  4,  4,  4,  4,  4,  4,  4,
 280 /* 10  11  12  13  14  15  16  17  18  19  1a  1b  1c  1d  1e  1f */
 281     3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,  3,
 282 /* 20  21  22  23  24  25  26  27  28  29  2a  2b  2c  2d  2e  2f */
 283     2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,
 284 /* 30  31  32  33  34  35  36  37  38  39  3a  3b  3c  3d  3e  3f */
 285     2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,  2,
 286 /* 40  41  42  43  44  45  46  47  48  49  4a  4b  4c  4d  4e  4f */
 287     1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
 288 /* 50  51  52  53  54  55  56  57  58  59  5a  5b  5c  5d  5e  5f */
 289     1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
 290 /* 60  61  62  63  64  65  66  67  68  69  6a  6b  6c  6d  6e  6f */
 291     1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
 292 /* 70  71  72  73  74  75  76  77  78  79  7a  7b  7c  7d  7e  7f */
 293     1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
 294 /* 80  81  82  83  84  78  86  87  88  89  8a  8b  8c  8d  8e  8f */
 295     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
 296 /* 90  91  92  93  94  98  96  97  98  99  9a  9b  9c  9d  9e  9f */
 297     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
 298 /* a0  a1  a2  a3  a4  a8  a6  a7  a8  a9  aa  ab  ac  ad  ae  af */
 299     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
 300 /* b0  b1  b2  b3  b4  b8  b6  b7  b8  b9  ba  bb  bc  bd  be  bf */
 301     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
 302 /* c0  c1  c2  c3  c4  c8  c6  c7  c8  c9  ca  cb  cc  cd  ce  cf */
 303     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
 304 /* d0  d1  d2  d3  d4  d8  d6  d7  d8  d9  da  db  dc  dd  de  df */
 305     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
 306 /* e0  e1  e2  e3  e4  e8  e6  e7  e8  e9  ea  eb  ec  ed  ee  ef */
 307     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
 308 /* f0  f1  f2  f3  f4  f8  f6  f7  f8  f9  fa  fb  fc  fd  fe  ff */
 309     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0
 310 };
 311
 312
 313 /*
 314  * Signed division
 315  */
 316 void i64_div(bigint dividend, bigint divisor, bigint *quotient, bigint *remainder)
 317 {
 318     int quotientIsMinus, remainderIsMinus;
 319
 320     if ((divisor.hi | divisor.lo) == 0) {
 321         return;
 322     }
 323     if ((dividend.hi | dividend.lo) == 0) {
 324         quotient->hi = 0;
 325         quotient->lo = 0;
 326         remainder->hi = 0;
 327         remainder->lo = 0;
 328         return;
 329     }
 330     quotientIsMinus = 0;
 331     remainderIsMinus = dividend.hi & 0x80000000;
 332     if (remainderIsMinus) {
 333         dividend = i64_inv(dividend);
 334         quotientIsMinus = ~quotientIsMinus;
 335     }
 336     if (divisor.hi & 0x80000000) {
 337         divisor = i64_inv(divisor);
 338         quotientIsMinus = ~quotientIsMinus;
 339     }
 340     i64_udiv_impl(dividend, divisor, quotient, remainder);
 341     if (quotientIsMinus) {
 342         *quotient = i64_inv(*quotient);
 343     }
 344     if (remainderIsMinus) {
 345         *remainder = i64_inv(*remainder);
 346     }
 347 }
 348
 349
 350 /*
 351  * Unsigned division
 352  */
 353 void i64_udiv(bigint dividend, bigint divisor, bigint *quotient, bigint *remainder)
 354 {
 355     if ((divisor.hi | divisor.lo) == 0) {
 356         return;
 357     }
 358     if ((dividend.hi | dividend.lo) == 0) {
 359         quotient->hi = 0;
 360         quotient->lo = 0;
 361         remainder->hi = 0;
 362         remainder->lo = 0;
 363         return;
 364     }
 365     i64_udiv_impl(dividend, divisor, quotient, remainder);
 366 }
 367
 368
 369 /*
 370  * Unsigned division implementation.
 371  *
 372  * Note: does not check argument values
 373  */
 374 static void i64_udiv_impl(bigint dividend, bigint divisor, bigint *quotient,
 375     bigint *remainder)
 376 {
 377     int i, j, m, n, d, divisorLeftZeros, dividendLeftZeros;
 378     unsigned v[4], u[5], q[4], qh, rh, a, b, c;
 379     bigint mm;
 380
 381     divisorLeftZeros = 0;
 382     mm = divisor;
 383     n = mm.hi >> 24;
 384     while (n == 0) {
 385         divisorLeftZeros += leftZerosInByte[n];
 386         mm.hi <<= 8;
 387         mm.hi += mm.lo >> 24;
 388         mm.lo <<= 8;
 389         n = mm.hi >> 24;
 390     };
 391     divisorLeftZeros += leftZerosInByte[n];
 392
 393     dividendLeftZeros = 0;
 394     mm = dividend;
 395     n = mm.hi >> 24;
 396     while (n == 0) {
 397         dividendLeftZeros += leftZerosInByte[n];
 398         mm.hi <<= 8;
 399         mm.hi += mm.lo >> 24;
 400         mm.lo <<= 8;
 401         n = mm.hi >> 24;
 402     };
 403     dividendLeftZeros += leftZerosInByte[n];
 404
 405     n = (79 - divisorLeftZeros) / 16;
 406     m = (79 - dividendLeftZeros) / 16 - n;
 407
 408     if (n == 1) {
 409         quotient->hi = dividend.hi / divisor.lo;
 410         rh = dividend.hi % divisor.lo;
 411         a = (rh << 16) + (dividend.lo >> 16);
 412         quotient->lo = (a / divisor.lo) << 16;
 413         rh = a % divisor.lo;
 414         a = (rh << 16) + (dividend.lo & 0xffff);
 415         quotient->lo += a / divisor.lo;
 416         remainder->hi = 0;
 417         remainder->lo = a % divisor.lo;
 418     }
 419     else /* if (n > 1) */ {
 420         d = divisorLeftZeros & 0x0f;
 421
 422         i = 0;
 423         u[i++] = (dividend.lo << d) & 0xffff;
 424         u[i++] = dividend.lo >> (16 - d) & 0xffff;
 425         u[i++] = ((dividend.hi << d) +
 426             ((dividend.lo >> 16) >> (16 - d))) & 0xffff;
 427         u[i++] = (dividend.hi >> (16 - d)) & 0xffff;
 428         u[i]   = (dividend.hi >> 16) >> (16 - d);
 429
 430         i = 0;
 431         v[i++] = (divisor.lo << d) & 0xffff;
 432         v[i++] = (divisor.lo >> (16 - d)) & 0xffff;
 433         v[i++] = ((divisor.hi << d) +
 434             ((divisor.lo >> 16) >> (16 - d))) & 0xffff;
 435         v[i]   = divisor.hi >> (16 - d);
 436
 437         memset(q, 0, sizeof(q));
 438
 439         for (j=m; j >= 0; --j) {
 440             a = (u[j+n] << 16) + u[j+n-1];
 441             qh = a / v[n-1];
 442             rh = a % v[n-1];
 443             if ((qh & 0x30000) || qh*v[n-2] > (rh << 16) + u[j+n-2]) {
 444                 --qh;
 445                 rh += v[n-1];
 446                 if (! (rh & 0x10000)) {
 447                     if ((qh & 0x30000) || qh*v[n-2] > (rh << 16) + u[j+n-2]) {
 448                         --qh;
 449                     }
 450                 }
 451             }
 452             b = 0;
 453             c = 1;
 454             for (i=0; i < n; ++i) {
 455                 b = qh * v[i] + b;
 456                 c = u[j+i] + c + ((~b) & 0xffff);
 457                 u[j+i] = c & 0xffff;
 458                 b = b >> 16;
 459                 c = c >> 16;
 460             }
 461             c = u[j+i] + c + ((~b) & 0xffff);
 462             u[j+i] = c & 0xffff;
 463             if (u[j+n]) {
 464                 c = 0;
 465                 for (i=0; i < n; ++i) {
 466                     c = u[j+i] + v[i] + c;
 467                     u[j+i] = c & 0xffff;
 468                     c = c >> 16;
 469                 }
 470                 c = u[j+i] + c;
 471                 u[j+i] = c & 0xffff;
 472                 --qh;
 473             }
 474             q[j] = qh;
 475         }
 476         quotient->hi = (q[3] << 16) + q[2];
 477         quotient->lo = (q[1] << 16) + q[0];
 478         remainder->hi = (u[3] << (16 - d)) + (u[2] >> d);
 479         remainder->lo = ((u[2] << 16) << (16 - d)) + (u[1] << (16 - d))
 480             + (u[0] >> d);
 481     }
 482 }
 483
 484
 485 /*
 486  * Convert string to integer
 487  */
 488 bigint i64_atoi(const char *str)
 489 {
 490     int i, isNegative;
 491     bigint ten, digit, n;
 492
 493     ten.hi = 0;
 494     ten.lo = 10;
 495     digit.hi = 0;
 496     digit.lo = 0;
 497     n.hi = 0;
 498     n.lo = 0;
 499     isNegative = 0;
 500     i = 0;
 501     while (isspace(str[i])) {
 502         ++i;
 503     }
 504     if (str[i] == '-') {
 505         isNegative = 1;
 506         ++i;
 507     }
 508     else if (str[i] == '+') {
 509         ++i;
 510     }
 511     while (isdigit(str[i])) {
 512         digit.lo = str[i] - '0';
 513         n = i64_mul(ten, n);
 514         n = i64_add(digit, n);
 515         ++i;
 516     }
 517     if (isNegative) {
 518         n = i64_inv(n);
 519     }
 520     return n;
 521 }
 522
 523
 524 /*
 525  * Convert integer to string
 526  */
 527 void i64_itoa(bigint n, char *str, int strSize)
 528 {
 529     char s[24];
 530     int pos;
 531     bigint ten, digit;
 532
 533     ten.hi = 0;
 534     ten.lo = 10;
 535     pos = 20;
 536     s[pos] = 0;
 537     /* positive number */
 538     if (i64_sgn(n) > 0) {
 539         while (n.hi != 0 || n.lo != 0) {
 540             i64_div(n, ten, &n, &digit);
 541             --pos;
 542             s[pos] = digit.lo + '0';
 543         }
 544     }
 545     /* negative number */
 546     else if (i64_sgn(n) < 0) {
 547         n = i64_inv(n);
 548         while (n.hi != 0 || n.lo != 0) {
 549             i64_udiv(n, ten, &n, &digit);
 550             --pos;
 551             s[pos] = digit.lo + '0';
 552         }
 553         --pos;
 554         s[pos] = '-';
 555     }
 556     /* zero */
 557     else /* if (i64_sgn(n) == 0) */ {
 558         --pos;
 559         s[pos] = '0';
 560     }
 561     strncpy(str, &s[pos], strSize);
 562 }
 563
 564 /* Added by Paul Huxham 30/12/2000 */
 565
 566 /* Set the value of a bigint from a long */
 567 bigint i64_set( long n1 )
 568 {
 569     bigint ret;
 570
 571     ret.hi = 0;
 572     ret.lo = n1;
 573
 574     if ( n1 < 1 ) i64_inv( ret );
 575
 576     return ret;
 577 }
 578
 579 /* Set the value of a bigint from a ulong */
 580 bigint i64_uset( unsigned long n1 )
 581 {
 582     bigint ret;
 583
 584     ret.hi = 0;
 585     ret.lo = n1;
 586
 587     return ret;
 588 }
 589
 590 /* Multiply two ulongs together and return a bigint */
 591 bigint i64_uumul( unsigned long n1, unsigned long n2 )
 592 {
 593     bigint a1, a2;
 594
 595     a1 = i64_uset( n1 );
 596     a2 = i64_uset( n2 );
 597
 598     return i64_mul( a1, a2 );
 599 }