c_cpp/etc/calc/qfunc.c

   1 /*
   2  * qfunc - extended precision rational arithmetic non-primitive functions
   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: qfunc.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
  24  * @(#) $Source: /usr/local/src/bin/calc/RCS/qfunc.c,v $
  25  *
  26  * Under source code control:   1990/02/15 01:48:20
  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 "qmath.h"
  34 #include "config.h"
  35 #include "prime.h"
  36
  37 STATIC NUMBER **B_table;
  38 STATIC long B_num;
  39 STATIC long B_allocnum;
  40 STATIC NUMBER **E_table;
  41 STATIC long E_num;
  42
  43 #define QALLOCNUM 64
  44
  45 /*
  46  * Set the default epsilon for approximate calculations.
  47  * This must be greater than zero.
  48  *
  49  * given:
  50  *      q               number to be set as the new epsilon
  51  */
  52 void
  53 setepsilon(NUMBER *q)
  54 {
  55         NUMBER *old;
  56
  57         if (qisneg(q) || qiszero(q)) {
  58                 math_error("Epsilon value must be greater than zero");
  59                 /*NOTREACHED*/
  60         }
  61         old = conf->epsilon;
  62         conf->epsilonprec = qprecision(q);
  63         conf->epsilon = qlink(q);
  64         if (old)
  65                 qfree(old);
  66 }
  67
  68
  69 /*
  70  * Return the inverse of one number modulo another.
  71  * That is, find x such that:
  72  *      Ax = 1 (mod B)
  73  * Returns zero if the numbers are not relatively prime (temporary hack).
  74  */
  75 NUMBER *
  76 qminv(NUMBER *q1, NUMBER *q2)
  77 {
  78         NUMBER *r;
  79         ZVALUE z1, z2, tmp;
  80         int s, t;
  81         long rnd;
  82
  83         if (qisfrac(q1) || qisfrac(q2)) {
  84                 math_error("Non-integers for minv");
  85                 /*NOTREACHED*/
  86         }
  87         if (qiszero(q2)) {
  88                 if (qisunit(q1))
  89                         return qlink(q1);
  90                 return qlink(&_qzero_);
  91         }
  92         if (qisunit(q2))
  93                 return qlink(&_qzero_);
  94         rnd = conf->mod;
  95         s = (rnd & 4) ? 0 : q2->num.sign;
  96         if (rnd & 1)
  97                 s^= 1;
  98
  99         tmp = q2->num;
 100         tmp.sign = 0;
 101         if (zmodinv(q1->num, tmp, &z1))
 102                 return qlink(&_qzero_);
 103         zsub(tmp, z1, &z2);
 104         if (rnd & 16) {
 105                 t = zrel(z1, z2);
 106                 if (t < 0)
 107                         s = 0;
 108                 else if (t > 0)
 109                         s = 1;
 110         }
 111         r = qalloc();
 112         if (s) {
 113                 zfree(z1);
 114                 z2.sign = TRUE;
 115                 r->num = z2;
 116                 return r;
 117         }
 118         zfree(z2);
 119         r->num = z1;
 120         return r;
 121 }
 122
 123
 124 /*
 125  * Return the residue modulo an integer (q3) of an integer (q1)
 126  * raised to a positive integer (q2) power.
 127  */
 128 NUMBER *
 129 qpowermod(NUMBER *q1, NUMBER *q2, NUMBER *q3)
 130 {
 131         NUMBER *r;
 132         ZVALUE z1, z2, tmp;
 133         int s, t;
 134         long rnd;
 135
 136         if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) {
 137                 math_error("Non-integers for pmod");
 138                 /*NOTREACHED*/
 139         }
 140         if (qisneg(q2)) {
 141                 math_error("Negative power for pmod");
 142                 /*NOTREACHED*/
 143         }
 144         if (qiszero(q3))
 145                 return qpowi(q1, q2);
 146         if (qisunit(q3))
 147                 return qlink(&_qzero_);
 148         rnd = conf->mod;
 149         s = (rnd & 4) ? 0 : q3->num.sign;
 150         if (rnd & 1)
 151                 s^= 1;
 152         tmp = q3->num;
 153         tmp.sign = 0;
 154         zpowermod(q1->num, q2->num, tmp, &z1);
 155         if (ziszero(z1)) {
 156                 zfree(z1);
 157                 return qlink(&_qzero_);
 158         }
 159         zsub(tmp, z1, &z2);
 160         if (rnd & 16) {
 161                 t = zrel(z1, z2);
 162                 if (t < 0)
 163                         s = 0;
 164                 else if (t > 0)
 165                         s = 1;
 166         }
 167         r = qalloc();
 168         if (s) {
 169                 zfree(z1);
 170                 z2.sign = TRUE;
 171                 r->num = z2;
 172                 return r;
 173         }
 174         zfree(z2);
 175         r->num = z1;
 176         return r;
 177 }
 178
 179
 180 /*
 181  * Return the power function of one number with another.
 182  * The power must be integral.
 183  *      q3 = qpowi(q1, q2);
 184  */
 185 NUMBER *
 186 qpowi(NUMBER *q1, NUMBER *q2)
 187 {
 188         register NUMBER *r;
 189         BOOL invert, sign;
 190         ZVALUE num, zden, z2;
 191
 192         if (qisfrac(q2)) {
 193                 math_error("Raising number to fractional power");
 194                 /*NOTREACHED*/
 195         }
 196         num = q1->num;
 197         zden = q1->den;
 198         z2 = q2->num;
 199         sign = (num.sign && zisodd(z2));
 200         invert = z2.sign;
 201         num.sign = 0;
 202         z2.sign = 0;
 203         /*
 204         * Check for trivial cases first.
 205         */
 206         if (ziszero(num) && !ziszero(z2)) {     /* zero raised to a power */
 207                 if (invert) {
 208                         math_error("Zero raised to negative power");
 209                         /*NOTREACHED*/
 210                 }
 211                 return qlink(&_qzero_);
 212         }
 213         if (zisunit(num) && zisunit(zden)) {    /* 1 or -1 raised to a power */
 214                 r = (sign ? q1 : &_qone_);
 215                 r->links++;
 216                 return r;
 217         }
 218         if (ziszero(z2))        /* raising to zeroth power */
 219                 return qlink(&_qone_);
 220         if (zisunit(z2)) {      /* raising to power 1 or -1 */
 221                 if (invert)
 222                         return qinv(q1);
 223                 return qlink(q1);
 224         }
 225         /*
 226          * Not a trivial case.  Do the real work.
 227          */
 228         r = qalloc();
 229         if (!zisunit(num))
 230                 zpowi(num, z2, &r->num);
 231         if (!zisunit(zden))
 232                 zpowi(zden, z2, &r->den);
 233         if (invert) {
 234                 z2 = r->num;
 235                 r->num = r->den;
 236                 r->den = z2;
 237         }
 238         r->num.sign = sign;
 239         return r;
 240 }
 241
 242
 243 /*
 244  * Given the legs of a right triangle, compute its hypotenuse within
 245  * the specified error.  This is sqrt(a^2 + b^2).
 246  */
 247 NUMBER *
 248 qhypot(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
 249 {
 250         NUMBER *tmp1, *tmp2, *tmp3;
 251
 252         if (qiszero(epsilon)) {
 253                 math_error("Zero epsilon value for hypot");
 254                 /*NOTREACHED*/
 255         }
 256         if (qiszero(q1))
 257                 return qqabs(q2);
 258         if (qiszero(q2))
 259                 return qqabs(q1);
 260         tmp1 = qsquare(q1);
 261         tmp2 = qsquare(q2);
 262         tmp3 = qqadd(tmp1, tmp2);
 263         qfree(tmp1);
 264         qfree(tmp2);
 265         tmp1 = qsqrt(tmp3, epsilon, 24L);
 266         qfree(tmp3);
 267         return tmp1;
 268 }
 269
 270
 271 /*
 272  * Given one leg of a right triangle with unit hypotenuse, calculate
 273  * the other leg within the specified error.  This is sqrt(1 - a^2).
 274  * If the wantneg flag is nonzero, then negative square root is returned.
 275  */
 276 NUMBER *
 277 qlegtoleg(NUMBER *q, NUMBER *epsilon, BOOL wantneg)
 278 {
 279         NUMBER *res, *qtmp1, *qtmp2;
 280         ZVALUE num;
 281
 282         if (qiszero(epsilon)) {
 283                 math_error("Zero epsilon value for legtoleg");
 284                 /*NOTREACHED*/
 285         }
 286         if (qisunit(q))
 287                 return qlink(&_qzero_);
 288         if (qiszero(q)) {
 289                 if (wantneg)
 290                         return qlink(&_qnegone_);
 291                 return qlink(&_qone_);
 292         }
 293         num = q->num;
 294         num.sign = 0;
 295         if (zrel(num, q->den) >= 0) {
 296                 math_error("Leg too large in legtoleg");
 297                 /*NOTREACHED*/
 298         }
 299         qtmp1 = qsquare(q);
 300         qtmp2 = qsub(&_qone_, qtmp1);
 301         qfree(qtmp1);
 302         res = qsqrt(qtmp2, epsilon, 24L);
 303         qfree(qtmp2);
 304         if (wantneg) {
 305                 qtmp1 = qneg(res);
 306                 qfree(res);
 307                 res = qtmp1;
 308         }
 309         return res;
 310 }
 311
 312
 313 /*
 314  * Compute the square root of a real number.
 315  * Type of rounding if any depends on rnd.
 316  * If rnd & 32 is nonzero, result is exact for square numbers;
 317  * If rnd & 64 is nonzero, the negative square root is returned;
 318  * If rnd < 32, result is rounded to a multiple of epsilon
 319  *      up, down, etc. depending on bits 0, 2, 4 of v.
 320  */
 321
 322 NUMBER *
 323 qsqrt(NUMBER *q1, NUMBER *epsilon, long rnd)
 324 {
 325         NUMBER *r, etemp;
 326         ZVALUE tmp1, tmp2, quo, mul, divisor;
 327         long  s1, s2, up, RR, RS;
 328         int sign;
 329
 330         if (qisneg(q1)) {
 331                 math_error("Square root of negative number");
 332                 /*NOTREACHED*/
 333         }
 334         if (qiszero(q1))
 335                 return qlink(&_qzero_);
 336         sign = (rnd & 64) != 0;
 337         if (qiszero(epsilon)) {
 338                 math_error("Zero epsilon for qsqrt");
 339                 /*NOTREACHED*/
 340         }
 341
 342         etemp = *epsilon;
 343         etemp.num.sign = 0;
 344         RS = rnd & 25;
 345         if (!(RS & 8))
 346                 RS ^= epsilon->num.sign;
 347         if (rnd & 2)
 348                 RS ^= sign ^ epsilon->num.sign;
 349         if (rnd & 4)
 350                 RS ^= epsilon->num.sign;
 351         RR = zisunit(q1->den) && qisunit(epsilon);
 352         if (rnd & 32 || RR) {
 353                 s1 = zsqrt(q1->num, &tmp1, RS);
 354                 if (RR) {
 355                         if (ziszero(tmp1)) {
 356                                 zfree(tmp1);
 357                                 return qlink(&_qzero_);
 358                         }
 359                         r = qalloc();
 360                         tmp1.sign = sign;
 361                         r->num = tmp1;
 362                         return r;
 363                 }
 364                 if (!s1) {
 365                         s2 = zsqrt(q1->den, &tmp2, 0);
 366                         if (!s2) {
 367                                 r = qalloc();
 368                                 tmp1.sign = sign;
 369                                 r->num = tmp1;
 370                                 r->den = tmp2;
 371                                 return r;
 372                         }
 373                         zfree(tmp2);
 374                 }
 375                 zfree(tmp1);
 376         }
 377         up = 0;
 378         zsquare(epsilon->den, &tmp1);
 379         zmul(tmp1, q1->num, &tmp2);
 380         zfree(tmp1);
 381         zsquare(epsilon->num, &tmp1);
 382         zmul(tmp1, q1->den, &divisor);
 383         zfree(tmp1);
 384         if (rnd & 16) {
 385                 zshift(tmp2, 2, &tmp1);
 386                 zfree(tmp2);
 387                 s1 = zquo(tmp1, divisor, &quo, 16);
 388                 zfree(tmp1);
 389                 s2 = zsqrt(quo, &tmp1, s1 ? s1 < 0 : 16);
 390                 zshift(tmp1, -1, &mul);
 391                 up = (*tmp1.v & 1) ? s1 + s2 : -1;
 392                 zfree(tmp1);
 393         } else {
 394                 s1 = zquo(tmp2, divisor, &quo, 0);
 395                 zfree(tmp2);
 396                 s2 = zsqrt(quo, &mul, 0);
 397                 up = (s1 + s2) ? 0 : -1;
 398         }
 399         if (up == 0) {
 400                 if (rnd & 8)
 401                         up = (long)((RS ^ *mul.v) & 1);
 402                 else
 403                         up = RS ^ sign;
 404         }
 405         if (up > 0) {
 406                 zadd(mul, _one_, &tmp2);
 407                 zfree(mul);
 408                 mul = tmp2;
 409         }
 410         zfree(divisor);
 411         zfree(quo);
 412         if (ziszero(mul)) {
 413                 zfree(mul);
 414                 return qlink(&_qzero_);
 415         }
 416         r = qalloc();
 417         zreduce(mul, etemp.den, &tmp1, &r->den);
 418         zfree(mul);
 419         tmp1.sign = sign;
 420         zmul(tmp1, etemp.num, &r->num);
 421         zfree(tmp1);
 422         return r;
 423 }
 424
 425
 426 /*
 427  * Calculate the integral part of the square root of a number.
 428  * Example:  qisqrt(13) = 3.
 429  */
 430 NUMBER *
 431 qisqrt(NUMBER *q)
 432 {
 433         NUMBER *r;
 434         ZVALUE tmp;
 435
 436         if (qisneg(q)) {
 437                 math_error("Square root of negative number");
 438                 /*NOTREACHED*/
 439         }
 440         if (qiszero(q))
 441                 return qlink(&_qzero_);
 442         r = qalloc();
 443         if (qisint(q)) {
 444                 (void) zsqrt(q->num, &r->num,0);
 445                 return r;
 446         }
 447         zquo(q->num, q->den, &tmp, 0);
 448         (void) zsqrt(tmp, &r->num,0);
 449         freeh(tmp.v);
 450         return r;
 451 }
 452
 453 /*
 454  * Return whether or not a number is an exact square.
 455  */
 456 BOOL
 457 qissquare(NUMBER *q)
 458 {
 459         BOOL flag;
 460
 461         flag = zissquare(q->num);
 462         if (qisint(q) || !flag)
 463                 return flag;
 464         return zissquare(q->den);
 465 }
 466
 467
 468 /*
 469  * Compute the greatest integer of the Kth root of a number.
 470  * Example:  qiroot(85, 3) = 4.
 471  */
 472 NUMBER *
 473 qiroot(NUMBER *q1, NUMBER *q2)
 474 {
 475         NUMBER *r;
 476         ZVALUE tmp;
 477
 478         if (qisneg(q2) || qiszero(q2) || qisfrac(q2)) {
 479                 math_error("Taking number to bad root value");
 480                 /*NOTREACHED*/
 481         }
 482         if (qiszero(q1))
 483                 return qlink(&_qzero_);
 484         if (qisone(q1) || qisone(q2))
 485                 return qlink(q1);
 486         if (qistwo(q2))
 487                 return qisqrt(q1);
 488         r = qalloc();
 489         if (qisint(q1)) {
 490                 zroot(q1->num, q2->num, &r->num);
 491                 return r;
 492         }
 493         zquo(q1->num, q1->den, &tmp, 0);
 494         zroot(tmp, q2->num, &r->num);
 495         zfree(tmp);
 496         return r;
 497 }
 498
 499
 500 /*
 501  * Return the greatest integer of the base 2 log of a number.
 502  * This is the number such that  1 <= x / log2(x) < 2.
 503  * Examples:  qilog2(8) = 3, qilog2(1.3) = 1, qilog2(1/7) = -3.
 504  *
 505  * given:
 506  *      q               number to take log of
 507  */
 508 long
 509 qilog2(NUMBER *q)
 510 {
 511         long n;                 /* power of two */
 512         int c;                  /* result of comparison */
 513         ZVALUE tmp1, tmp2;      /* temporary values */
 514
 515         if (qiszero(q)) {
 516                 math_error("Zero argument for ilog2");
 517                 /*NOTREACHED*/
 518         }
 519         if (qisint(q))
 520                 return zhighbit(q->num);
 521         tmp1 = q->num;
 522         tmp1.sign = 0;
 523         n = zhighbit(tmp1) - zhighbit(q->den);
 524         if (n == 0)
 525                 c = zrel(tmp1, q->den);
 526         else if (n > 0) {
 527                 zshift(q->den, n, &tmp2);
 528                 c = zrel(tmp1, tmp2);
 529                 zfree(tmp2);
 530         } else {
 531                 zshift(tmp1, -n, &tmp2);
 532                 c = zrel(tmp2, q->den);
 533                 zfree(tmp2);
 534         }
 535         if (c < 0)
 536                 n--;
 537         return n;
 538 }
 539
 540
 541 /*
 542  * Return the greatest integer of the base 10 log of a number.
 543  * This is the number such that  1 <= x / log10(x) < 10.
 544  * Examples:  qilog10(100) = 2, qilog10(12.3) = 1, qilog10(.023) = -2.
 545  *
 546  * given:
 547  *      q               number to take log of
 548  */
 549 long
 550 qilog10(NUMBER *q)
 551 {
 552         long n;                 /* log value */
 553         ZVALUE tmp1, tmp2;      /* temporary values */
 554
 555         if (qiszero(q)) {
 556                 math_error("Zero argument for ilog10");
 557                 /*NOTREACHED*/
 558         }
 559         tmp1 = q->num;
 560         tmp1.sign = 0;
 561         if (qisint(q))
 562                 return zlog10(tmp1, NULL);
 563         /*
 564          * The number is not an integer.
 565          * Compute the result if the number is greater than one.
 566          */
 567         if (zrel(tmp1, q->den) > 0) {
 568                 zquo(tmp1, q->den, &tmp2, 0);
 569                 n = zlog10(tmp2, NULL);
 570                 zfree(tmp2);
 571                 return n;
 572         }
 573         /*
 574          * Here if the number is less than one.
 575          * If the number is the inverse of a power of ten, then the
 576          * obvious answer will be off by one.  Subtracting one if the
 577          * number is the inverse of an integer will fix it.
 578          */
 579         if (zisunit(tmp1))
 580                 zsub(q->den, _one_, &tmp2);
 581         else
 582                 zquo(q->den, tmp1, &tmp2, 0);
 583         n = -zlog10(tmp2, NULL) - 1;
 584         zfree(tmp2);
 585         return n;
 586 }
 587
 588 /*
 589  * Return the integer floor of the logarithm of a number relative to
 590  * a specified integral base.
 591  */
 592 NUMBER *
 593 qilog(NUMBER *q, ZVALUE base)
 594 {
 595         long n;
 596         ZVALUE tmp1, tmp2;
 597
 598         if (qiszero(q))
 599                 return NULL;
 600
 601         if (qisunit(q))
 602                 return qlink(&_qzero_);
 603         if (qisint(q))
 604                 return itoq(zlog(q->num, base));
 605         tmp1 = q->num;
 606         tmp1.sign = 0;
 607         if (zrel(tmp1, q->den) > 0) {
 608                 zquo(tmp1, q->den, &tmp2, 0);
 609                 n = zlog(tmp2, base);
 610                 zfree(tmp2);
 611                 return itoq(n);
 612         }
 613         if (zisunit(tmp1))
 614                 zsub(q->den, _one_, &tmp2);
 615         else
 616                 zquo(q->den, tmp1, &tmp2, 0);
 617         n = -zlog(tmp2, base) - 1;
 618         zfree(tmp2);
 619         return itoq(n);
 620 }
 621
 622
 623 /*
 624  * Return the number of digits in the representation to a specified
 625  * base of the integral part of a number.
 626  *
 627  * Examples: qdigits(3456,10) = 4, qdigits(-23.45, 10) = 2.
 628  *
 629  * One should remember these special cases:
 630  *
 631  *      digits(12.3456) == 2    computes with the integer part only
 632  *      digits(-1234) == 4      computes with the absolute value only
 633  *      digits(0) == 1          specical case
 634  *      digits(-0.123) == 1     combination of all of the above
 635  *
 636  * given:
 637  *      q               number to count digits of
 638  */
 639 long
 640 qdigits(NUMBER *q, ZVALUE base)
 641 {
 642         long n;                 /* number of digits */
 643         ZVALUE temp;            /* temporary value */
 644
 645         if (zabsrel(q->num, q->den) < 0)
 646                 return 1;
 647         if (qisint(q))
 648                 return 1 + zlog(q->num, base);
 649         zquo(q->num, q->den, &temp, 2);
 650         n = 1 + zlog(temp, base);
 651         zfree(temp);
 652         return n;
 653 }
 654
 655
 656 /*
 657  * Return the digit at the specified place in the expansion to specified
 658  * base of a rational number.  The places specified by dpos are numbered from
 659  * the "units" place just before the "decimal" point, so that negative
 660  * dpos indicates the (-dpos)th place to the right of the point.
 661  * Examples: qdigit(1234.5678, 1, 10) = 3, qdigit(1234.5678, -3, 10) = 7.
 662  * The signs of the number and the base are ignored.
 663  */
 664 NUMBER *
 665 qdigit(NUMBER *q, ZVALUE dpos, ZVALUE base)
 666 {
 667         ZVALUE N, D;
 668         ZVALUE K;
 669         long k;
 670         ZVALUE A, B, C;         /* temporary integers */
 671         NUMBER *res;
 672
 673         /*
 674          * In the first stage, q is expressed as base^k * N/D where
 675          *      gcd(D, base) = 1
 676          * K is k as a ZVALUE
 677          */
 678         base.sign = 0;
 679         if (ziszero(base) || zisunit(base))
 680                         return NULL;
 681         if (qiszero(q) || (qisint(q) && zisneg(dpos)) ||
 682                         (zge31b(dpos) && !zisneg(dpos)))
 683                 return qlink(&_qzero_);
 684         k = zfacrem(q->num, base, &N);
 685         if (k == 0) {
 686                 k = zgcdrem(q->den, base, &D);
 687                 if (k > 0) {
 688                         zequo(q->den, D, &A);
 689                         itoz(k, &K);
 690                         zpowi(base, K, &B);
 691                         zfree(K);
 692                         zequo(B, A, &C);
 693                         zfree(A);
 694                         zfree(B);
 695                         zmul(C, q->num, &N);
 696                         zfree(C);
 697                         k = -k;
 698                 }
 699                 else
 700                         N = q->num;
 701         }
 702         if (k >= 0)
 703                 D = q->den;
 704
 705         itoz(k, &K);
 706         if (zrel(dpos, K) >= 0) {
 707                 zsub(dpos, K, &A);
 708                 zfree(K);
 709                 zpowi(base, A, &B);
 710                 zfree(A);
 711                 zmul(D, B, &A);
 712                 zfree(B);
 713                 zquo(N, A, &B, 0);
 714         } else {
 715                 if (zisunit(D)) {
 716                         if (k != 0)
 717                                 zfree(N);
 718                         zfree(K);
 719                         if (k < 0)
 720                                 zfree(D);
 721                         return qlink(&_qzero_);
 722                 }
 723                 zsub(K, dpos, &A);
 724                 zfree(K);
 725                 zpowermod(base, A, D, &C);
 726                 zfree(A);
 727                 zmod(N, D, &A, 0);
 728                 zmul(C, A, &B);
 729                 zfree(A);
 730                 zfree(C);
 731                 zmod(B, D, &A, 0);
 732                 zfree(B);
 733                 zmodinv(D, base, &B);
 734                 zsub(base, B, &C);
 735                 zfree(B);
 736                 zmul(C, A, &B);
 737                 zfree(C);
 738         }
 739         zfree(A);
 740         if (k != 0)
 741                 zfree(N);
 742         if (k < 0)
 743                 zfree(D);
 744         zmod(B, base, &A, 0);
 745         zfree(B);
 746         if (ziszero(A)) {
 747                 zfree(A);
 748                 return qlink(&_qzero_);
 749         }
 750         if (zisone(A)) {
 751                 zfree(A);
 752                 return qlink(&_qone_);
 753         }
 754         res = qalloc();
 755         res->num = A;
 756         return res;
 757 }
 758
 759
 760 /*
 761  * Return whether or not a bit is set at the specified bit position in a
 762  * number.  The lowest bit of the integral part of a number is the zeroth
 763  * bit position.  Negative bit positions indicate bits to the right of the
 764  * binary decimal point.  Examples: qdigit(17.1, 0) = 1, qdigit(17.1, -1) = 0.
 765  */
 766 BOOL
 767 qisset(NUMBER *q, long n)
 768 {
 769         NUMBER *qtmp1, *qtmp2;
 770         ZVALUE ztmp;
 771         BOOL res;
 772
 773         /*
 774          * Zero number or negative bit position place of integer is trivial.
 775          */
 776         if (qiszero(q) || (qisint(q) && (n < 0)))
 777                 return FALSE;
 778         /*
 779          * For non-negative bit positions, answer is easy.
 780          */
 781         if (n >= 0) {
 782                 if (qisint(q))
 783                         return zisset(q->num, n);
 784                 zquo(q->num, q->den, &ztmp, 2);
 785                 res = zisset(ztmp, n);
 786                 zfree(ztmp);
 787                 return res;
 788         }
 789         /*
 790          * Fractional value and want negative bit position, must work harder.
 791          */
 792         qtmp1 = qscale(q, -n);
 793         qtmp2 = qint(qtmp1);
 794         qfree(qtmp1);
 795         res = ((qtmp2->num.v[0] & 0x01) != 0);
 796         qfree(qtmp2);
 797         return res;
 798 }
 799
 800
 801 /*
 802  * Compute the factorial of an integer.
 803  *      q2 = qfact(q1);
 804  */
 805 NUMBER *
 806 qfact(NUMBER *q)
 807 {
 808         register NUMBER *r;
 809
 810         if (qisfrac(q)) {
 811                 math_error("Non-integral factorial");
 812                 /*NOTREACHED*/
 813         }
 814         if (qiszero(q) || zisone(q->num))
 815                 return qlink(&_qone_);
 816         r = qalloc();
 817         zfact(q->num, &r->num);
 818         return r;
 819 }
 820
 821
 822 /*
 823  * Compute the product of the primes less than or equal to a number.
 824  *      q2 = qpfact(q1);
 825  */
 826 NUMBER *
 827 qpfact(NUMBER *q)
 828 {
 829         NUMBER *r;
 830
 831         if (qisfrac(q)) {
 832                 math_error("Non-integral factorial");
 833                 /*NOTREACHED*/
 834         }
 835         r = qalloc();
 836         zpfact(q->num, &r->num);
 837         return r;
 838 }
 839
 840
 841 /*
 842  * Compute the lcm of all the numbers less than or equal to a number.
 843  *      q2 = qlcmfact(q1);
 844  */
 845 NUMBER *
 846 qlcmfact(NUMBER *q)
 847 {
 848         NUMBER *r;
 849
 850         if (qisfrac(q)) {
 851                 math_error("Non-integral lcmfact");
 852                 /*NOTREACHED*/
 853         }
 854         r = qalloc();
 855         zlcmfact(q->num, &r->num);
 856         return r;
 857 }
 858
 859
 860 /*
 861  * Compute the permutation function  q1 * (q1-1) * ... * (q1-q2+1).
 862  */
 863 NUMBER *
 864 qperm(NUMBER *q1, NUMBER *q2)
 865 {
 866         NUMBER *r;
 867         NUMBER *qtmp1, *qtmp2;
 868         long i;
 869
 870         if (qisfrac(q2)) {
 871                 math_error("Non-integral second arg for permutation");
 872                 /*NOTREACHED*/
 873         }
 874         if (qiszero(q2))
 875                 return qlink(&_qone_);
 876         if (qisone(q2))
 877                 return qlink(q1);
 878         if (qisint(q1) && !qisneg(q1)) {
 879                 if (qrel(q2, q1) > 0)
 880                         return qlink(&_qzero_);
 881                 if (qispos(q2)) {
 882                         r = qalloc();
 883                         zperm(q1->num, q2->num, &r->num);
 884                         return r;
 885                 }
 886         }
 887         if (zge31b(q2->num)) {
 888                 math_error("Too large arg2 for permutation");
 889                 /*NOTREACHED*/
 890         }
 891         i = qtoi(q2);
 892         if (i > 0) {
 893                 q1 = qlink(q1);
 894                 r = qlink(q1);
 895                 while (--i > 0) {
 896                         qtmp1 = qdec(q1);
 897                         qtmp2 = qmul(r, qtmp1);
 898                         qfree(q1);
 899                         q1 = qtmp1;
 900                         qfree(r);
 901                         r = qtmp2;
 902                 }
 903                 qfree(q1);
 904                 return r;
 905         }
 906         i = -i;
 907         qtmp1 = qinc(q1);
 908         r = qinv(qtmp1);
 909         while (--i > 0) {
 910                 qtmp2 = qinc(qtmp1);
 911                 qfree(qtmp1);
 912                 qtmp1 = qqdiv(r, qtmp2);
 913                 qfree(r);
 914                 r = qtmp1;
 915                 qtmp1 = qtmp2;
 916         }
 917         qfree(qtmp1);
 918         return r;
 919 }
 920
 921
 922 /*
 923  * Compute the combinatorial function  q(q - 1) ...(q - n + 1)/n!
 924  * n is to be a nonnegative integer
 925  */
 926 NUMBER *
 927 qcomb(NUMBER *q, NUMBER *n)
 928 {
 929         NUMBER *r;
 930         NUMBER *q1, *q2;
 931         long i, j;
 932         ZVALUE z;
 933
 934         if (!qisint(n) || qisneg(n)) {
 935                 math_error("Bad second arg in call to qcomb!");
 936                 /*NOTREACHED*/
 937         }
 938         if (qisint(q)) {
 939                 switch (zcomb(q->num, n->num, &z)) {
 940                 case 0:
 941                         return qlink(&_qzero_);
 942                 case 1:
 943                         return qlink(&_qone_);
 944                 case -1:
 945                         return qlink(&_qnegone_);
 946                 case 2:
 947                         return qlink(q);
 948                 case -2:
 949                         return NULL;
 950                 default:
 951                         r = qalloc();
 952                         r->num = z;
 953                         return r;
 954                 }
 955         }
 956         if (zge31b(n->num))
 957                 return NULL;
 958         i = ztoi(n->num);
 959         q = qlink(q);
 960         r = qlink(q);
 961         j = 1;
 962         while (--i > 0) {
 963                 q1 = qdec(q);
 964                 qfree(q);
 965                 q = q1;
 966                 q2 = qmul(r, q);
 967                 qfree(r);
 968                 r = qdivi(q2, ++j);
 969                 qfree(q2);
 970         }
 971         qfree(q);
 972         return r;
 973 }
 974
 975
 976 /*
 977  * Compute the Bernoulli number with index n
 978  * For even positive n, newly calculated values for all even indices up
 979  * to n are stored in table at B_table.
 980  */
 981 NUMBER *
 982 qbern(ZVALUE z)
 983 {
 984         long n, i, k, m, nn, dd;
 985         NUMBER **p;
 986         NUMBER *s, *s1, *c, *c1, *t;
 987         size_t sz;
 988
 989         if (zisone(z))
 990                 return qlink(&_qneghalf_);
 991
 992         if (zisodd(z) || z.sign)
 993                 return qlink(&_qzero_);
 994
 995         if (zge31b(z))
 996                 return NULL;
 997
 998         n = ztoi(z);
 999
1000         if (n == 0)
1001
1002                 return qlink(&_qone_);
1003
1004         m = (n >> 1) - 1;
1005
1006         if (m < B_num)
1007                 return qlink(B_table[m]);
1008
1009         if (m >= B_allocnum) {
1010                 k = (m/QALLOCNUM + 1) * QALLOCNUM;
1011                 sz = k * sizeof(NUMBER *);
1012                 if (sz < (size_t) k)
1013                         return NULL;
1014                 if (B_allocnum == 0)
1015                         p = (NUMBER **) malloc(sz);
1016                 else
1017                         p = (NUMBER **) realloc(B_table, sz);
1018                 if (p == NULL)
1019                         return NULL;
1020                 B_allocnum = k;
1021                 B_table = p;
1022         }
1023         for (k = B_num; k <= m; k++) {
1024                 nn = 2 * k + 3;
1025                 dd = 1;
1026                 c1 = itoq(nn);
1027                 c = qinv(c1);
1028                 qfree(c1);
1029                 s = qsub(&_qonehalf_, c);
1030                 i = k;
1031                 for (i = 0; i < k; i++) {
1032                         c1 = qmuli(c, nn--);
1033                         qfree(c);
1034                         c = qdivi(c1, dd++);
1035                         qfree(c1);
1036                         c1 = qmuli(c, nn--);
1037                         qfree(c);
1038                         c = qdivi(c1, dd++);
1039                         qfree(c1);
1040                         t = qmul(c, B_table[i]);
1041                         s1 = qsub(s, t);
1042                         qfree(t);
1043                         qfree(s);
1044                         s = s1;
1045                 }
1046                 B_table[k] = s;
1047                 qfree(c);
1048         }
1049         B_num = k;
1050         return qlink(B_table[m]);
1051 }
1052
1053
1054 void
1055 qfreebern(void)
1056 {
1057         long k;
1058
1059         if (B_num > 0) {
1060                 for (k = 0; k < B_num; k++)
1061                         qfree(B_table[k]);
1062                 free(B_table);
1063         }
1064         B_table = NULL;
1065         B_allocnum = 0;
1066         B_num = 0;
1067 }
1068
1069
1070 /*
1071  * Compute the Euler number with index n = z
1072  * For even positive n, newly calculated values with all even indices up
1073  * to n are stored in E_table for later quick access.
1074  */
1075 NUMBER *
1076 qeuler(ZVALUE z)
1077 {
1078         long i, k, m, n, nn, dd;
1079         NUMBER **p;
1080         NUMBER *s, *s1, *c, *c1, *t;
1081         size_t sz;
1082
1083
1084         if (ziszero(z))
1085                 return qlink(&_qone_);
1086         if (zisodd(z) || zisneg(z))
1087                 return qlink(&_qzero_);
1088         if (zge31b(z))
1089                 return NULL;
1090         n = ztoi(z);
1091         m = (n >> 1) - 1;
1092         if (m < E_num)
1093                 return qlink(E_table[m]);
1094         sz = (m + 1) * sizeof(NUMBER *);
1095         if (sz < (size_t) m + 1)
1096                 return NULL;
1097         if (E_num)
1098                 p = (NUMBER **) realloc(E_table, sz);
1099         else
1100                 p = (NUMBER **) malloc(sz);
1101         if (p == NULL)
1102                 return NULL;
1103         E_table = p;
1104         for (k = E_num; k <= m; k++) {
1105                 nn = 2 * k + 2;
1106                 dd = 1;
1107                 c = qlink(&_qone_);
1108                 s = qlink(&_qnegone_);
1109                 i = k;
1110                 for (i = 0; i < k; i++) {
1111                         c1 = qmuli(c, nn--);
1112                         qfree(c);
1113                         c = qdivi(c1, dd++);
1114                         qfree(c1);
1115                         c1 = qmuli(c, nn--);
1116                         qfree(c);
1117                         c = qdivi(c1, dd++);
1118                         qfree(c1);
1119                         t = qmul(c, E_table[i]);
1120                         s1 = qsub(s, t);
1121                         qfree(t);
1122                         qfree(s);
1123                         s = s1;
1124                 }
1125                 E_table[k] = s;
1126                 qfree(c);
1127         }
1128         E_num = k;
1129         return qlink(E_table[m]);
1130 }
1131
1132
1133 void
1134 qfreeeuler(void)
1135 {
1136         long k;
1137
1138         if (E_num > 0) {
1139                 for (k = 0; k < E_num; k++)
1140                         qfree(E_table[k]);
1141                 free(E_table);
1142         }
1143         E_table = NULL;
1144         E_num = 0;
1145 }
1146
1147
1148 /*
1149  * Catalan numbers: catalan(n) = comb(2*n, n)/(n+1)
1150  * To be called only with integer q
1151  */
1152 NUMBER *
1153 qcatalan(NUMBER *q)
1154 {
1155         NUMBER *A, *B;
1156         NUMBER *res;
1157
1158         if (qisneg(q))
1159                 return qlink(&_qzero_);
1160         A = qscale(q, 1);
1161         B = qcomb(A, q);
1162         if (B == NULL)
1163                 return NULL;
1164         qfree(A);
1165         A = qinc(q);
1166         res = qqdiv(B, A);
1167         qfree(A);
1168         qfree(B);
1169         return res;
1170 }
1171
1172 /*
1173  * Compute the Jacobi function (a / b).
1174  * -1 => a is not quadratic residue mod b
1175  * 1 => b is composite, or a is quad residue of b
1176  * 0 => b is even or b < 0
1177  */
1178 NUMBER *
1179 qjacobi(NUMBER *q1, NUMBER *q2)
1180 {
1181         if (qisfrac(q1) || qisfrac(q2)) {
1182                 math_error("Non-integral arguments for jacobi");
1183                 /*NOTREACHED*/
1184         }
1185         return itoq((long) zjacobi(q1->num, q2->num));
1186 }
1187
1188
1189 /*
1190  * Compute the Fibonacci number F(n).
1191  */
1192 NUMBER *
1193 qfib(NUMBER *q)
1194 {
1195         register NUMBER *r;
1196
1197         if (qisfrac(q)) {
1198                 math_error("Non-integral Fibonacci number");
1199                 /*NOTREACHED*/
1200         }
1201         r = qalloc();
1202         zfib(q->num, &r->num);
1203         return r;
1204 }
1205
1206
1207 /*
1208  * Truncate a number to the specified number of decimal places.
1209  */
1210 NUMBER *
1211 qtrunc(NUMBER *q1, NUMBER *q2)
1212 {
1213         long places;
1214         NUMBER *r, *e;
1215
1216         if (qisfrac(q2) || !zistiny(q2->num)) {
1217                 math_error("Bad number of places for qtrunc");
1218                 /*NOTREACHED*/
1219         }
1220         places = qtoi(q2);
1221         e = qtenpow(-places);
1222         r = qmappr(q1, e, 2);
1223         qfree(e);
1224         return r;
1225 }
1226
1227
1228
1229
1230 /*
1231  * Truncate a number to the specified number of binary places.
1232  * Specifying zero places makes the result identical to qint.
1233  */
1234 NUMBER *
1235 qbtrunc(NUMBER *q1, NUMBER *q2)
1236 {
1237         long places;
1238         NUMBER *r, *e;
1239
1240         if (qisfrac(q2) || !zistiny(q2->num)) {
1241                 math_error("Bad number of places for qtrunc");
1242                 /*NOTREACHED*/
1243         }
1244         places = qtoi(q2);
1245         e = qbitvalue(-places);
1246         r = qmappr(q1, e, 2);
1247         qfree(e);
1248         return r;
1249 }
1250
1251
1252 /*
1253  * Round a number to a specified number of binary places.
1254  */
1255 NUMBER *
1256 qbround(NUMBER *q, long places, long rnd)
1257 {
1258         NUMBER *e, *r;
1259
1260         if (qiszero(q))
1261                 return qlink(&_qzero_);
1262         if (rnd & 32)
1263                 places -= qilog2(q) + 1;
1264         e = qbitvalue(-places);
1265         r = qmappr(q, e, rnd & 31);
1266         qfree(e);
1267         return r;
1268 }
1269
1270 /*
1271  * Round a number to a specified number of decimal digits.
1272  */
1273 NUMBER *
1274 qround(NUMBER *q, long places, long rnd)
1275 {
1276         NUMBER *e, *r;
1277
1278         if (qiszero(q))
1279                 return qlink(&_qzero_);
1280         if (rnd & 32)
1281                 places -= qilog10(q) + 1;
1282         e = qtenpow(-places);
1283         r = qmappr(q, e, rnd & 31);
1284         qfree(e);
1285         return r;
1286 }
1287
1288 /*
1289  * Approximate a number to nearest multiple of a given number.  Whether
1290  * rounding is down, up, etc. is determined by rnd.
1291  */
1292 NUMBER *
1293 qmappr(NUMBER *q, NUMBER *e, long rnd)
1294 {
1295         NUMBER *r;
1296         ZVALUE tmp1, tmp2, mul;
1297
1298         if (qiszero(e))
1299                 return qlink(q);
1300         if (qiszero(q))
1301                 return qlink(&_qzero_);
1302         zmul(q->num, e->den, &tmp1);
1303         zmul(q->den, e->num, &tmp2);
1304         zquo(tmp1, tmp2, &mul, rnd);
1305         zfree(tmp1);
1306         zfree(tmp2);
1307         if (ziszero(mul)) {
1308                 zfree(mul);
1309                 return qlink(&_qzero_);
1310         }
1311         r = qalloc();
1312         zreduce(mul, e->den, &tmp1, &r->den);
1313         zmul(tmp1, e->num, &r->num);
1314         zfree(tmp1);
1315         zfree(mul);
1316         return r;
1317 }
1318
1319
1320 /*
1321  * Determine the smallest-denominator rational number in the interval of
1322  * length abs(epsilon) (< 1) with centre or one end at q, or to determine
1323  * the number nearest above or nearest below q with denominator
1324  * not exceeding abs(epsilon).
1325  * Whether the approximation is nearest above or nearest below is
1326  * determined by rnd and the signs of epsilon and q.
1327  */
1328
1329 NUMBER *
1330 qcfappr(NUMBER *q, NUMBER *epsilon, long rnd)
1331 {
1332         NUMBER *res, etemp, *epsilon1;
1333         ZVALUE num, zden, oldnum, oldden;
1334         ZVALUE rem, oldrem, quot;
1335         ZVALUE tmp1, tmp2, tmp3, tmp4;
1336         ZVALUE denbnd;
1337         ZVALUE f, g, k;
1338         BOOL esign;
1339         int s;
1340         BOOL bnddencase;
1341         BOOL useold = FALSE;
1342
1343         if (qiszero(epsilon) || qisint(q))
1344                 return qlink(q);
1345
1346         esign = epsilon->num.sign;
1347         etemp = *epsilon;
1348         etemp.num.sign = 0;
1349         bnddencase = (zrel(etemp.num, etemp.den) >= 0);
1350         if (bnddencase) {
1351                 zquo(etemp.num, etemp.den, &denbnd, 0);
1352                 if (zrel(q->den, denbnd) <= 0) {
1353                         zfree(denbnd);
1354                         return qlink(q);
1355                 }
1356         } else {
1357                 if (rnd & 16)
1358                         epsilon1 = qscale(epsilon, -1);
1359                 else
1360                         epsilon1 = qlink(epsilon);
1361                 zreduce(q->den, epsilon1->den, &tmp1, &g);
1362                 zmul(epsilon1->num, tmp1, &f);
1363                 f.sign = 0;
1364                 zfree(tmp1);
1365                 qfree(epsilon1);
1366         }
1367         if (rnd & 16 && !zistwo(q->den)) {
1368                 s = 0;
1369         } else {
1370                 s = esign ? -1 : 1;
1371                 if (rnd & 1)
1372                         s = -s;
1373                 if (rnd & 2 && q->num.sign ^ esign)
1374                         s = -s;
1375                 if (rnd & 4 && esign)
1376                         s = -s;
1377         }
1378         oldnum = _one_;
1379         oldden = _zero_;
1380         zcopy(q->den, &oldrem);
1381         zdiv(q->num, q->den, &num, &rem, 0);
1382         zden = _one_;
1383         for (;;) {
1384                 if (!bnddencase) {
1385                         zmul(f, zden, &tmp1);
1386                         zmul(g, rem, &tmp2);
1387                         if (ziszero(rem) || (s >= 0 && zrel(tmp1,tmp2) >= 0))
1388                                 break;
1389                         zfree(tmp1);
1390                         zfree(tmp2);
1391                 }
1392                 zdiv(oldrem, rem, &quot, &tmp1, 0);
1393                 zfree(oldrem);
1394                 oldrem = rem;
1395                 rem = tmp1;
1396                 zmul(quot, zden, &tmp1);
1397                 zadd(tmp1, oldden, &tmp2);
1398                 zfree(tmp1);
1399                 zfree(oldden);
1400                 oldden = zden;
1401                 zden = tmp2;
1402                 zmul(quot, num, &tmp1);
1403                 zadd(tmp1, oldnum, &tmp2);
1404                 zfree(tmp1);
1405                 zfree(oldnum);
1406                 oldnum = num;
1407                 num = tmp2;
1408                 zfree(quot);
1409                 if (bnddencase) {
1410                         if (zrel(zden, denbnd) >= 0)
1411                                 break;
1412                 }
1413                 s = -s;
1414         }
1415         if (bnddencase) {
1416                 if (s > 0) {
1417                         useold = TRUE;
1418                 } else {
1419                         zsub(zden, denbnd, &tmp1);
1420                         zquo(tmp1, oldden, &k, 1);
1421                         zfree(tmp1);
1422                 }
1423                 zfree(denbnd);
1424         } else {
1425                 if (s < 0) {
1426                         zfree(tmp1);
1427                         zfree(tmp2);
1428                         zfree(f);
1429                         zfree(g);
1430                         zfree(oldnum);
1431                         zfree(oldden);
1432                         zfree(num);
1433                         zfree(zden);
1434                         zfree(oldrem);
1435                         zfree(rem);
1436                         return qlink(q);
1437                 }
1438                 zsub(tmp1, tmp2, &tmp3);
1439                 zfree(tmp1);
1440                 zfree(tmp2);
1441                 zmul(f, oldden, &tmp1);
1442                 zmul(g, oldrem, &tmp2);
1443                 zfree(f);
1444                 zfree(g);
1445                 zadd(tmp1, tmp2, &tmp4);
1446                 zfree(tmp1);
1447                 zfree(tmp2);
1448                 zquo(tmp3, tmp4, &k, 0);
1449                 zfree(tmp3);
1450                 zfree(tmp4);
1451         }
1452         if (!useold && !ziszero(k)) {
1453                 zmul(k, oldnum, &tmp1);
1454                 zsub(num, tmp1, &tmp2);
1455                 zfree(tmp1);
1456                 zfree(num);
1457                 num = tmp2;
1458                 zmul(k, oldden, &tmp1);
1459                 zsub(zden, tmp1, &tmp2);
1460                 zfree(tmp1);
1461                 zfree(zden);
1462                 zden = tmp2;
1463         }
1464         if (bnddencase && s == 0) {
1465                 zmul(k, oldrem, &tmp1);
1466                 zadd(rem, tmp1, &tmp2);
1467                 zfree(tmp1);
1468                 zfree(rem);
1469                 rem = tmp2;
1470                 zmul(rem, oldden, &tmp1);
1471                 zmul(zden, oldrem, &tmp2);
1472                 useold = (zrel(tmp1, tmp2) >= 0);
1473                 zfree(tmp1);
1474                 zfree(tmp2);
1475         }
1476         if (!bnddencase || s <= 0)
1477                 zfree(k);
1478         zfree(rem);
1479         zfree(oldrem);
1480         res = qalloc();
1481         if (useold) {
1482                 zfree(num);
1483                 zfree(zden);
1484                 res->num = oldnum;
1485                 res->den = oldden;
1486                 return res;
1487         }
1488         zfree(oldnum);
1489         zfree(oldden);
1490         res->num = num;
1491         res->den = zden;
1492         return res;
1493 }
1494
1495
1496 /*
1497  * Calculate the nearest-above, or nearest-below, or nearest, number
1498  * with denominator less than the given number, the choice between
1499  * possibilities being determined by the parameter rnd.
1500  */
1501 NUMBER *
1502 qcfsim(NUMBER *q, long rnd)
1503 {
1504         NUMBER *res;
1505         ZVALUE tmp1, tmp2, den1, den2;
1506         int s;
1507
1508         if (qiszero(q) && rnd & 26)
1509                 return qlink(&_qzero_);
1510         if (rnd & 24) {
1511                 s = q->num.sign;
1512         } else {
1513                 s = rnd & 1;
1514                 if (rnd & 2)
1515                         s ^= q->num.sign;
1516         }
1517         if (qisint(q)) {
1518                 if ((rnd & 8) && !(rnd & 16))
1519                         return qlink(&_qzero_);
1520                 if (s)
1521                         return qinc(q);
1522                 return qdec(q);
1523         }
1524         if (zistwo(q->den)) {
1525                 if (rnd & 16)
1526                         s ^= 1;
1527                 if (s)
1528                         zadd(q->num, _one_, &tmp1);
1529                 else
1530                         zsub(q->num, _one_, &tmp1);
1531                 res = qalloc();
1532                 zshift(tmp1, -1, &res->num);
1533                 zfree(tmp1);
1534                 return res;
1535         }
1536         s = s ? 1 : -1;
1537         if (rnd & 24)
1538                 s = 0;
1539         res = qalloc();
1540         zmodinv(q->num, q->den, &den1);
1541         if (s >= 0) {
1542                 zsub(q->den, den1, &den2);
1543                 if (s > 0 || ((zrel(den1, den2) < 0) ^ (!(rnd & 16)))) {
1544                         zfree(den1);
1545                         res->den = den2;
1546                         zmul(den2, q->num, &tmp1);
1547                         zadd(tmp1, _one_, &tmp2);
1548                         zfree(tmp1);
1549                         zequo(tmp2, q->den, &res->num);
1550                         zfree(tmp2);
1551                         return res;
1552                 }
1553                 zfree(den2);
1554         }
1555         res->den = den1;
1556         zmul(den1, q->num, &tmp1);
1557         zsub(tmp1, _one_, &tmp2);
1558         zfree(tmp1);
1559         zequo(tmp2, q->den, &res->num);
1560         zfree(tmp2);
1561         return res;
1562 }
1563
1564
1565
1566 /*
1567  * Return an indication on whether or not two fractions are approximately
1568  * equal within the specified epsilon. Returns negative if the absolute value
1569  * of the difference between the two numbers is less than epsilon, zero if
1570  * the difference is equal to epsilon, and positive if the difference is
1571  * greater than epsilon.
1572  */
1573 FLAG
1574 qnear(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
1575 {
1576         int res;
1577         NUMBER qtemp, etemp, *qq;
1578
1579         etemp = *epsilon;
1580         etemp.num.sign = 0;
1581         if (q1 == q2) {
1582                 if (qiszero(epsilon))
1583                         return 0;
1584                 return -1;
1585         }
1586         if (qiszero(epsilon))
1587                 return qcmp(q1, q2);
1588         if (qiszero(q2)) {
1589                 qtemp = *q1;
1590                 qtemp.num.sign = 0;
1591                 return qrel(&qtemp, &etemp);
1592         }
1593         if (qiszero(q1)) {
1594                 qtemp = *q2;
1595                 qtemp.num.sign = 0;
1596                 return qrel(&qtemp, &etemp);
1597         }
1598         qq = qsub(q1, q2);
1599         qtemp = *qq;
1600         qtemp.num.sign = 0;
1601         res = qrel(&qtemp, &etemp);
1602         qfree(qq);
1603         return res;
1604 }
1605
1606
1607 /*
1608  * Compute the gcd (greatest common divisor) of two numbers.
1609  *      q3 = qgcd(q1, q2);
1610  */
1611 NUMBER *
1612 qgcd(NUMBER *q1, NUMBER *q2)
1613 {
1614         ZVALUE z;
1615         NUMBER *q;
1616
1617         if (q1 == q2)
1618                 return qqabs(q1);
1619         if (qisfrac(q1) || qisfrac(q2)) {
1620                 q = qalloc();
1621                 zgcd(q1->num, q2->num, &q->num);
1622                 zlcm(q1->den, q2->den, &q->den);
1623                 return q;
1624         }
1625         if (qiszero(q1))
1626                 return qqabs(q2);
1627         if (qiszero(q2))
1628                 return qqabs(q1);
1629         if (qisunit(q1) || qisunit(q2))
1630                 return qlink(&_qone_);
1631         zgcd(q1->num, q2->num, &z);
1632         if (zisunit(z)) {
1633                 zfree(z);
1634                 return qlink(&_qone_);
1635         }
1636         q = qalloc();
1637         q->num = z;
1638         return q;
1639 }
1640
1641
1642 /*
1643  * Compute the lcm (least common multiple) of two numbers.
1644  *      q3 = qlcm(q1, q2);
1645  */
1646 NUMBER *
1647 qlcm(NUMBER *q1, NUMBER *q2)
1648 {
1649         NUMBER *q;
1650
1651         if (qiszero(q1) || qiszero(q2))
1652                 return qlink(&_qzero_);
1653         if (q1 == q2)
1654                 return qqabs(q1);
1655         if (qisunit(q1))
1656                 return qqabs(q2);
1657         if (qisunit(q2))
1658                 return qqabs(q1);
1659         q = qalloc();
1660         zlcm(q1->num, q2->num, &q->num);
1661         if (qisfrac(q1) || qisfrac(q2))
1662                 zgcd(q1->den, q2->den, &q->den);
1663         return q;
1664 }
1665
1666
1667 /*
1668  * Remove all occurrences of the specified factor from a number.
1669  * Returned number is always positive or zero.
1670  */
1671 NUMBER *
1672 qfacrem(NUMBER *q1, NUMBER *q2)
1673 {
1674         long count;
1675         ZVALUE tmp;
1676         NUMBER *r;
1677
1678         if (qisfrac(q1) || qisfrac(q2)) {
1679                 math_error("Non-integers for factor removal");
1680                 /*NOTREACHED*/
1681         }
1682         if (qiszero(q2))
1683                 return qqabs(q1);
1684         if (qiszero(q1))
1685                 return qlink(&_qzero_);
1686         count = zfacrem(q1->num, q2->num, &tmp);
1687         if (zisunit(tmp)) {
1688                 zfree(tmp);
1689                 return qlink(&_qone_);
1690         }
1691         if (count == 0 && !qisneg(q1)) {
1692                 zfree(tmp);
1693                 return qlink(q1);
1694         }
1695         r = qalloc();
1696         r->num = tmp;
1697         return r;
1698 }
1699
1700
1701 /*
1702  * Divide one number by the gcd of it with another number repeatedly until
1703  * the number is relatively prime.
1704  */
1705 NUMBER *
1706 qgcdrem(NUMBER *q1, NUMBER *q2)
1707 {
1708         ZVALUE tmp;
1709         NUMBER *r;
1710
1711         if (qisfrac(q1) || qisfrac(q2)) {
1712                 math_error("Non-integers for gcdrem");
1713                 /*NOTREACHED*/
1714         }
1715         if (qiszero(q2))
1716                 return qlink(&_qone_);
1717         if (qiszero(q1))
1718                 return qlink(&_qzero_);
1719         if (zgcdrem(q1->num, q2->num, &tmp) == 0)
1720                 return qqabs(q1);
1721         if (zisunit(tmp)) {
1722                 zfree(tmp);
1723                 return qlink(&_qone_);
1724         }
1725         if (zcmp(q1->num, tmp) == 0) {
1726                 zfree(tmp);
1727                 return qlink(q1);
1728         }
1729         r = qalloc();
1730         r->num = tmp;
1731         return r;
1732 }
1733
1734
1735 /*
1736  * Return the lowest prime factor of a number.
1737  * Search is conducted for the specified number of primes.
1738  * Returns one if no factor was found.
1739  */
1740 NUMBER *
1741 qlowfactor(NUMBER *q1, NUMBER *q2)
1742 {
1743         unsigned long count;
1744
1745         if (qisfrac(q1) || qisfrac(q2)) {
1746                 math_error("Non-integers for lowfactor");
1747                 /*NOTREACHED*/
1748         }
1749         count = ztoi(q2->num);
1750         if (count > PIX_32B) {
1751                 math_error("lowfactor count is too large");
1752                 /*NOTREACHED*/
1753         }
1754         return utoq(zlowfactor(q1->num, count));
1755 }
1756
1757 /*
1758  * Return the number of places after the decimal point needed to exactly
1759  * represent the specified number as a real number.  Integers return zero,
1760  * and non-terminating decimals return minus one.  Examples:
1761  * qdecplaces(1.23) = 2, qdecplaces(3) = 0, qdecplaces(1/7) = -1.
1762  */
1763 long
1764 qdecplaces(NUMBER *q)
1765 {
1766         long twopow, fivepow;
1767         HALF fiveval[2];
1768         ZVALUE five;
1769         ZVALUE tmp;
1770
1771         if (qisint(q))  /* no decimal places if number is integer */
1772                 return 0;
1773         /*
1774          * The number of decimal places of a fraction in lowest terms is finite
1775          * if an only if the denominator is of the form 2^A * 5^B, and then the
1776          * number of decimal places is equal to MAX(A, B).
1777          */
1778         five.sign = 0;
1779         five.len = 1;
1780         five.v = fiveval;
1781         fiveval[0] = 5;
1782         fivepow = zfacrem(q->den, five, &tmp);
1783         if (!zisonebit(tmp)) {
1784                 zfree(tmp);
1785                 return -1;
1786         }
1787         twopow = zlowbit(tmp);
1788         zfree(tmp);
1789         if (twopow < fivepow)
1790                 twopow = fivepow;
1791         return twopow;
1792 }
1793
1794
1795 /*
1796  * Return, if possible, the minimum number of places after the decimal
1797  * point needed to exactly represent q with the specified base.
1798  * Integers return 0 and numbers with non-terminating expansions -1.
1799  * Returns -2 if base inadmissible
1800  */
1801 long
1802 qplaces(NUMBER *q, ZVALUE base)
1803 {
1804         long count;
1805         ZVALUE tmp;
1806
1807         if (base.len == 1 && base.v[0] == 10)
1808                 return qdecplaces(q);
1809         if (ziszero(base) || zisunit(base))
1810                         return -2;
1811         if (qisint(q))
1812                 return 0;
1813         if (zisonebit(base)) {
1814                 if (!zisonebit(q->den))
1815                         return -1;
1816                 return 1 + (zlowbit(q->den) - 1)/zlowbit(base);
1817         }
1818         count = zgcdrem(q->den, base, &tmp);
1819         if (count == 0)
1820                 return -1;
1821         if (!zisunit(tmp))
1822                 count = -1;
1823         zfree(tmp);
1824         return count;
1825 }
1826
1827
1828 /*
1829  * Perform a probabilistic primality test (algorithm P in Knuth).
1830  * Returns FALSE if definitely not prime, or TRUE if probably prime.
1831  *
1832  * The absolute value of the 2nd arg determines how many times
1833  * to check for primality.  If 2nd arg < 0, then the trivial
1834  * check is omitted.  The 3rd arg determines how tests to
1835  * initially skip.
1836  */
1837 BOOL
1838 qprimetest(NUMBER *q1, NUMBER *q2, NUMBER *q3)
1839 {
1840         if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) {
1841                 math_error("Bad arguments for ptest");
1842                 /*NOTREACHED*/
1843         }
1844         if (zge24b(q2->num)) {
1845                 math_error("ptest count >= 2^24");
1846                 /*NOTREACHED*/
1847         }
1848         return zprimetest(q1->num, ztoi(q2->num), q3->num);
1849 }