c_cpp/etc/calc/qmath.c

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