c_cpp/etc/calc/qtrans.c

   1 /*
   2  * qtrans - transcendental functions for real numbers
   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: qtrans.c,v 30.2 2008/10/24 07:09:41 chongo Exp $
  24  * @(#) $Source: /usr/local/src/bin/calc/RCS/qtrans.c,v $
  25  *
  26  * Under source code control:   1990/02/15 01:48:22
  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  * Transcendental functions for real numbers.
  34  * These are sin, cos, exp, ln, power, cosh, sinh.
  35  */
  36
  37
  38 #include "qmath.h"
  39
  40 HALF _qlgenum_[] = { 36744 };
  41 HALF _qlgeden_[] = { 25469 };
  42 NUMBER _qlge_ = { { _qlgenum_, 1, 0 }, { _qlgeden_, 1, 0 }, 1, NULL };
  43
  44 /*
  45  * cache the natural logarithm of 10
  46  */
  47 STATIC NUMBER *ln_10 = NULL;
  48 STATIC NUMBER *ln_10_epsilon = NULL;
  49
  50 STATIC NUMBER *pivalue[2];
  51 STATIC NUMBER *qexprel(NUMBER *q, long bitnum);
  52
  53 /*
  54  * Evaluate and store in specified locations the sin and cos of a given
  55  * number to accuracy corresponding to a specified number of binary digits.
  56  */
  57 void
  58 qsincos(NUMBER *q, long bitnum, NUMBER **vs, NUMBER **vc)
  59 {
  60         long n, m, k, h, s, t, d;
  61         NUMBER *qtmp1, *qtmp2;
  62         ZVALUE X, cossum, sinsum, mul, ztmp1, ztmp2, ztmp3;
  63
  64         qtmp1 = qqabs(q);
  65         h = qilog2(qtmp1);
  66         qfree(qtmp1);
  67         k = bitnum + h + 1;
  68         if (k < 0) {
  69                 *vs = qlink(&_qzero_);
  70                 *vc = qlink(&_qone_);
  71                 return;
  72         }
  73         s = k;
  74         if (k) {
  75                 do {
  76                         t = s;
  77                         s = (s + k/s)/2;
  78                 }
  79                 while (t > s);
  80         }               /* s is int(sqrt(k)) */
  81         s++;
  82         if (s < -h)
  83                 s = -h;
  84         n = h + s;      /* n is number of squarings that will be required */
  85         m = bitnum + n;
  86         while (s > 0) {         /* increasing m by ilog2(s) */
  87                 s >>= 1;
  88                 m++;
  89         }                       /* m is working number of bits */
  90         qtmp1 = qscale(q, m - n);
  91         zquo(qtmp1->num, qtmp1->den, &X, 24);
  92         qfree(qtmp1);
  93         if (ziszero(X)) {
  94                 zfree(X);
  95                 *vs = qlink(&_qzero_);
  96                 *vc = qlink(&_qone_);
  97                 return;
  98         }
  99         zbitvalue(m, &cossum);
 100         zcopy(X, &sinsum);
 101         zcopy(X, &mul);
 102         d = 1;
 103         for (;;) {
 104                 X.sign = !X.sign;
 105                 zmul(X, mul, &ztmp1);
 106                 zfree(X);
 107                 zshift(ztmp1, -m, &ztmp2);
 108                 zfree(ztmp1);
 109                 zdivi(ztmp2, ++d, &X);
 110                 zfree(ztmp2);
 111                 if (ziszero(X))
 112                         break;
 113                 zadd(cossum, X, &ztmp1);
 114                 zfree(cossum);
 115                 cossum = ztmp1;
 116                 zmul(X, mul, &ztmp1);
 117                 zfree(X);
 118                 zshift(ztmp1, -m, &ztmp2);
 119                 zfree(ztmp1);
 120                 zdivi(ztmp2, ++d, &X);
 121                 zfree(ztmp2);
 122                 if (ziszero(X))
 123                         break;
 124                 zadd(sinsum, X, &ztmp1);
 125                 zfree(sinsum);
 126                 sinsum = ztmp1;
 127         }
 128         zfree(X);
 129         zfree(mul);
 130         while (n-- > 0) {
 131                 zsquare(cossum, &ztmp1);
 132                 zsquare(sinsum, &ztmp2);
 133                 zsub(ztmp1, ztmp2, &ztmp3);
 134                 zfree(ztmp1);
 135                 zfree(ztmp2);
 136                 zmul(cossum, sinsum, &ztmp1);
 137                 zfree(cossum);
 138                 zfree(sinsum);
 139                 zshift(ztmp3, -m, &cossum);
 140                 zfree(ztmp3);
 141                 zshift(ztmp1, 1 - m, &sinsum);
 142                 zfree(ztmp1);
 143         }
 144         h = zlowbit(cossum);
 145         qtmp1 = qalloc();
 146         if (m > h) {
 147                 zshift(cossum, -h, &qtmp1->num);
 148                 zbitvalue(m - h, &qtmp1->den);
 149         } else {
 150                 zshift(cossum, - m, &qtmp1->num);
 151         }
 152         zfree(cossum);
 153         *vc = qtmp1;
 154         h = zlowbit(sinsum);
 155         qtmp2 = qalloc();
 156         if (m > h) {
 157                 zshift(sinsum, -h, &qtmp2->num);
 158                 zbitvalue(m - h, &qtmp2->den);
 159         } else {
 160                 zshift(sinsum, -m, &qtmp2->num);
 161         }
 162         zfree(sinsum);
 163         *vs = qtmp2;
 164         return;
 165 }
 166
 167 /*
 168  * Calculate the cosine of a number to a near multiple of epsilon.
 169  * This calls qsincos() and discards the value of sin.
 170  */
 171 NUMBER *
 172 qcos(NUMBER *q, NUMBER *epsilon)
 173 {
 174         NUMBER *sin, *cos, *res;
 175         long n;
 176
 177         if (qiszero(epsilon)) {
 178                 math_error("Zero epsilon value for cosine");
 179                 /*NOTREACHED*/
 180         }
 181         if (qiszero(q))
 182                 return qlink(&_qone_);
 183         n = -qilog2(epsilon);
 184         if (n < 0)
 185                 return qlink(&_qzero_);
 186         qsincos(q, n + 2, &sin, &cos);
 187         qfree(sin);
 188         res = qmappr(cos, epsilon, 24);
 189         qfree(cos);
 190         return res;
 191 }
 192
 193
 194
 195 /*
 196  * This calls qsincos() and discards the value of cos.
 197  */
 198 NUMBER *
 199 qsin(NUMBER *q, NUMBER *epsilon)
 200 {
 201         NUMBER *sin, *cos, *res;
 202         long n;
 203
 204         if (qiszero(epsilon)) {
 205                 math_error("Zero epsilon value for sine");
 206                 /*NOTREACHED*/
 207         }
 208         n = -qilog2(epsilon);
 209         if (qiszero(q) || n < 0)
 210                 return qlink(&_qzero_);
 211         qsincos(q, n + 2, &sin, &cos);
 212         qfree(cos);
 213         res = qmappr(sin, epsilon, 24);
 214         qfree(sin);
 215         return res;
 216 }
 217
 218
 219 /*
 220  * Calculate the tangent function.
 221  */
 222 NUMBER *
 223 qtan(NUMBER *q, NUMBER *epsilon)
 224 {
 225         NUMBER *sin, *cos, *tan, *res;
 226         long n, k, m;
 227
 228         if (qiszero(epsilon)) {
 229                 math_error("Zero epsilon value for tangent");
 230                 /*NOTREACHED*/
 231         }
 232         if (qiszero(q))
 233                 return qlink(q);
 234         n = qilog2(epsilon);
 235         k = (n > 0) ? 4 + n/2 : 4;
 236         for (;;) {
 237                 qsincos(q, 2 * k - n, &sin, &cos);
 238                 if (qiszero(cos)) {
 239                         qfree(sin);
 240                         qfree(cos);
 241                         k = 2 * k - n + 4;
 242                         continue;
 243                 }
 244                 m = -qilog2(cos);
 245                 if (m < k)
 246                          break;
 247                 qfree(sin);
 248                 qfree(cos);
 249                 k = m + 1;
 250         }
 251         tan = qqdiv(sin, cos);
 252         qfree(sin);
 253         qfree(cos);
 254         res = qmappr(tan, epsilon, 24);
 255         qfree(tan);
 256         return res;
 257 }
 258
 259
 260 /*
 261  * Calculate the cotangent function.
 262  */
 263 NUMBER *
 264 qcot(NUMBER *q, NUMBER *epsilon)
 265 {
 266         NUMBER *sin, *cos, *cot, *res;
 267         long n, k, m;
 268
 269         if (qiszero(epsilon)) {
 270                 math_error("Zero epsilon value for cotangent");
 271                 /*NOTREACHED*/
 272         }
 273         if (qiszero(q)) {
 274                 math_error("Zero argument for cotangent");
 275                 /*NOTREACHED*/
 276         }
 277         k = -qilog2(q);
 278         n = qilog2(epsilon);
 279         if (k < 0)
 280                 k = (n > 0) ? n/2 : 0;
 281         k += 4;
 282         for (;;) {
 283                 qsincos(q, 2 * k - n, &sin, &cos);
 284                 if (qiszero(sin)) {
 285                         qfree(sin);
 286                         qfree(cos);
 287                         k = 2 * k - n + 4;
 288                         continue;
 289                 }
 290                 m = -qilog2(sin);
 291                 if (m < k)
 292                          break;
 293                 qfree(sin);
 294                 qfree(cos);
 295                 k = m + 1;
 296         }
 297         cot = qqdiv(cos, sin);
 298         qfree(sin);
 299         qfree(cos);
 300         res = qmappr(cot, epsilon, 24);
 301         qfree(cot);
 302         return res;
 303 }
 304
 305
 306 /*
 307  * Calculate the secant function.
 308  */
 309 NUMBER *
 310 qsec(NUMBER *q, NUMBER *epsilon)
 311 {
 312         NUMBER *sin, *cos, *sec, *res;
 313         long n, k, m;
 314
 315         if (qiszero(epsilon)) {
 316                 math_error("Zero epsilon value for secant");
 317                 /*NOTREACHED*/
 318         }
 319         if (qiszero(q))
 320                 return qlink(&_qone_);
 321         n = qilog2(epsilon);
 322         k = (n > 0) ? 4 + n/2 : 4;
 323         for (;;) {
 324                 qsincos(q, 2 * k - n, &sin, &cos);
 325                 qfree(sin);
 326                 if (qiszero(cos)) {
 327                         qfree(cos);
 328                         k = 2 * k - n + 4;
 329                         continue;
 330                 }
 331                 m = -qilog2(cos);
 332                 if (m < k)
 333                          break;
 334                 qfree(cos);
 335                 k = m + 1;
 336         }
 337         sec = qinv(cos);
 338         qfree(cos);
 339         res = qmappr(sec, epsilon, 24);
 340         qfree(sec);
 341         return res;
 342 }
 343
 344
 345 /*
 346  * Calculate the cosecant function.
 347  */
 348 NUMBER *
 349 qcsc(NUMBER *q, NUMBER *epsilon)
 350 {
 351         NUMBER *sin, *cos, *csc, *res;
 352         long n, k, m;
 353
 354         if (qiszero(epsilon)) {
 355                 math_error("Zero epsilon value for cosecant");
 356                 /*NOTREACHED*/
 357         }
 358         if (qiszero(q)) {
 359                 math_error("Zero argument for cosecant");
 360                 /*NOTREACHED*/
 361         }
 362         k = -qilog2(q);
 363         n = qilog2(epsilon);
 364         if (k < 0)
 365                 k = (n > 0) ? n/2 : 0;
 366         k += 4;
 367         for (;;) {
 368                 qsincos(q, 2 * k - n, &sin, &cos);
 369                 qfree(cos);
 370                 if (qiszero(sin)) {
 371                         qfree(sin);
 372                         k = 2 * k - n + 4;
 373                         continue;
 374                 }
 375                 m = -qilog2(sin);
 376                 if (m < k)
 377                          break;
 378                 qfree(sin);
 379                 k = m + 1;
 380         }
 381         csc = qinv(sin);
 382         qfree(sin);
 383         res = qmappr(csc, epsilon, 24);
 384         qfree(csc);
 385         return res;
 386 }
 387 /*
 388  * Calculate the arcsine function.
 389  * The result is in the range -pi/2 to pi/2.
 390  */
 391 NUMBER *
 392 qasin(NUMBER *q, NUMBER *epsilon)
 393 {
 394         NUMBER *qtmp1, *qtmp2, *epsilon1;
 395         ZVALUE ztmp;
 396         BOOL neg;
 397         FLAG r;
 398
 399         if (qiszero(epsilon)) {
 400                 math_error("Zero epsilon value for asin");
 401                 /*NOTREACHED*/
 402         }
 403         if (qiszero(q))
 404                 return qlink(&_qzero_);
 405         ztmp = q->num;
 406         neg = ztmp.sign;
 407         ztmp.sign = 0;
 408         r = zrel(ztmp, q->den);
 409         if (r > 0)
 410                 return NULL;
 411         if (r == 0) {
 412                 epsilon1 = qscale(epsilon, 1L);
 413                 qtmp2 = qpi(epsilon1);
 414                 qtmp1 = qscale(qtmp2, -1L);
 415         } else {
 416                 epsilon1 = qscale(epsilon, -2L);
 417                 qtmp1 = qalloc();
 418                 zsquare(q->num, &qtmp1->num);
 419                 zsquare(q->den, &ztmp);
 420                 zsub(ztmp, qtmp1->num, &qtmp1->den);
 421                 zfree(ztmp);
 422                 qtmp2 = qsqrt(qtmp1, epsilon1, 24);
 423                 qfree(qtmp1);
 424                 qtmp1 = qatan(qtmp2, epsilon);
 425         }
 426         qfree(qtmp2);
 427         qfree(epsilon1);
 428         if (neg) {
 429                 qtmp2 = qneg(qtmp1);
 430                 qfree(qtmp1);
 431                 return(qtmp2);
 432         }
 433         return qtmp1;
 434 }
 435
 436
 437
 438 /*
 439  * Calculate the acos function.
 440  * The result is in the range 0 to pi.
 441  */
 442 NUMBER *
 443 qacos(NUMBER *q, NUMBER *epsilon)
 444 {
 445         NUMBER *q1, *q2, *epsilon1;
 446         ZVALUE z;
 447
 448         if (qiszero(epsilon)) {
 449                 math_error("Zero epsilon value for acos");
 450                 /*NOTREACHED*/
 451         }
 452         if (qisone(q))
 453                 return qlink(&_qzero_);
 454         if (qisnegone(q))
 455                 return qpi(epsilon);
 456
 457         z = q->num;
 458         z.sign = 0;
 459         if (zrel(z, q->den) > 0)
 460                 return NULL;
 461         epsilon1 = qscale(epsilon, -3L);                /* ??? */
 462         q1 = qalloc();
 463         zsub(q->den, q->num, &q1->num);
 464         zadd(q->den, q->num, &q1->den);
 465         q2 = qsqrt(q1, epsilon1, 24L);
 466         qfree(q1);
 467         qfree(epsilon1);
 468         epsilon1 = qscale(epsilon, -1L);
 469         q1 = qatan(q2, epsilon1);
 470         qfree(epsilon1);
 471         qfree(q2);
 472         q2 = qscale(q1, 1L);
 473         qfree(q1)
 474         return q2;
 475 }
 476
 477
 478 /*
 479  * Calculate the arctangent function to the nearest or next to nearest
 480  * multiple of epsilon.  Algorithm uses
 481  *      atan(x) = 2 * atan(x/(1 + sqrt(1+x^2)))
 482  * to reduce x to a small value and then
 483  *      atan(x) = x - x^3/3 + ...
 484  */
 485 NUMBER *
 486 qatan(NUMBER *q, NUMBER *epsilon)
 487 {
 488         long m, k, i;
 489         FULL d;
 490         ZVALUE X, D, DD, sum, mul, term, ztmp1, ztmp2;
 491         NUMBER *qtmp, *res;
 492         BOOL sign;
 493
 494         if (qiszero(epsilon)) {
 495                 math_error("Zero epsilon value for arctangent");
 496                 /*NOTREACHED*/
 497         }
 498         if (qiszero(q))
 499                 return qlink(&_qzero_);
 500         m = 12 - qilog2(epsilon);
 501                 /* 4 bits for 4 doublings; 8 for rounding */
 502         if (m < 8)
 503                 m = 8;          /* m is number of working binary digits */
 504         qtmp = qscale(q, m);
 505         zquo(qtmp->num, qtmp->den, &X, 24);
 506         qfree(qtmp);
 507         zbitvalue(m, &D);       /* q has become X/D */
 508         zsquare(D, &DD);
 509         i = 4;                  /* maybe this should be larger */
 510         while (i-- > 0 && !ziszero(X)) {
 511                 zsquare(X, &ztmp1);
 512                 zadd(ztmp1, DD, &ztmp2);
 513                 zfree(ztmp1);
 514                 zsqrt(ztmp2, &ztmp1, 24L);
 515                 zfree(ztmp2);
 516                 zadd(ztmp1, D, &ztmp2);
 517                 zfree(ztmp1);
 518                 zshift(X, m, &ztmp1);
 519                 zfree(X);
 520                 zquo(ztmp1, ztmp2, &X, 24L);
 521                 zfree(ztmp1);
 522                 zfree(ztmp2);
 523         }
 524         zfree(DD);
 525         zfree(D);
 526         if (ziszero(X)) {
 527                 zfree(X);
 528                 return qlink(&_qzero_);
 529         }
 530         zcopy(X, &sum);
 531         zsquare(X, &ztmp1);
 532         zshift(ztmp1, -m, &mul);
 533         zfree(ztmp1);
 534         d = 3;
 535         sign = !X.sign;
 536         for (;;) {
 537                 if (d > BASE) {
 538                         math_error("Too many terms required for atan");
 539                         /*NOTREACHED*/
 540                 }
 541                 zmul(X, mul, &ztmp1);
 542                 zfree(X);
 543                 zshift(ztmp1, -m, &X);  /* X now (original X)^d */
 544                 zfree(ztmp1);
 545                 zdivi(X, d, &term);
 546                 if (ziszero(term)) {
 547                         zfree(term);
 548                         break;
 549                 }
 550                 term.sign = sign;
 551                 zadd(sum, term, &ztmp1);
 552                 zfree(sum);
 553                 zfree(term);
 554                 sum = ztmp1;
 555                 sign = !sign;
 556                 d += 2;
 557         }
 558         zfree(mul);
 559         zfree(X);
 560         qtmp = qalloc();
 561         k = zlowbit(sum);
 562         if (k) {
 563                 zshift(sum, -k, &qtmp->num);
 564                 zfree(sum);
 565         } else {
 566                 qtmp->num = sum;
 567         }
 568         zbitvalue(m - 4 - k, &qtmp->den);
 569         res = qmappr(qtmp, epsilon, 24L);
 570         qfree(qtmp);
 571         return res;
 572 }
 573
 574 /*
 575  * Inverse secant function
 576  */
 577 NUMBER *
 578 qasec(NUMBER *q, NUMBER *epsilon)
 579 {
 580         NUMBER *tmp, *res;
 581
 582         tmp = qinv(q);
 583         res = qacos(tmp, epsilon);
 584         qfree(tmp);
 585         return res;
 586 }
 587
 588
 589 /*
 590  * Inverse cosecant function
 591  */
 592 NUMBER *
 593 qacsc(NUMBER *q, NUMBER *epsilon)
 594 {
 595         NUMBER *tmp, *res;
 596
 597         tmp = qinv(q);
 598         res = qasin(tmp, epsilon);
 599         qfree(tmp);
 600         return res;
 601 }
 602
 603
 604 /*
 605  * Inverse cotangent function
 606  */
 607 NUMBER *
 608 qacot(NUMBER *q, NUMBER *epsilon)
 609 {
 610         NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
 611
 612         if (qiszero(epsilon)) {
 613                 math_error("Zero epsilon for acot");
 614                 /*NOTREACHED*/
 615         }
 616         if (qiszero(q)) {
 617                 epsilon1 = qscale(epsilon, 1L);
 618                 tmp1 = qpi(epsilon1);
 619                 qfree(epsilon1);
 620                 tmp2 = qscale(tmp1, -1L);
 621                 qfree(tmp1);
 622                 return tmp2;
 623         }
 624         tmp1 = qinv(q);
 625         if (!qisneg(q)) {
 626                 tmp2 = qatan(tmp1, epsilon);
 627                 qfree(tmp1);
 628                 return tmp2;
 629         }
 630         epsilon1 = qscale(epsilon, -2L);
 631         tmp2 = qatan(tmp1, epsilon1);
 632         qfree(tmp1);
 633         tmp1 = qpi(epsilon1);
 634         qfree(epsilon1);
 635         tmp3 = qqadd(tmp1, tmp2);
 636         qfree(tmp1);
 637         qfree(tmp2);
 638         tmp1 = qmappr(tmp3, epsilon, 24L);
 639         qfree(tmp3);
 640         return tmp1;
 641 }
 642
 643
 644 /*
 645  * Calculate the angle which is determined by the point (x,y).
 646  * This is the same as atan(y/x) for positive x, but is continuous
 647  * except for y = 0, x <= 0.  By convention, y is the first argument.
 648  * For all x, y, -pi < atan2 <= pi.  For example, qatan2(1, -1) = 3/4 * pi.
 649  */
 650 NUMBER *
 651 qatan2(NUMBER *qy, NUMBER *qx, NUMBER *epsilon)
 652 {
 653         NUMBER *tmp1, *tmp2, *tmp3, *epsilon2;
 654
 655         if (qiszero(epsilon)) {
 656                 math_error("Zero epsilon value for atan2");
 657                 /*NOTREACHED*/
 658         }
 659         if (qiszero(qy) && qiszero(qx)) {
 660                 /* conform to 4.3BSD ANSI/IEEE 754-1985 math lib */
 661                 return qlink(&_qzero_);
 662         }
 663         /*
 664          * If the point is on the negative real axis, then the answer is pi.
 665          */
 666         if (qiszero(qy) && qisneg(qx))
 667                 return qpi(epsilon);
 668         /*
 669          * If the point is in the right half plane, then use the normal atan.
 670          */
 671         if (!qisneg(qx) && !qiszero(qx)) {
 672                 if (qiszero(qy))
 673                         return qlink(&_qzero_);
 674                 tmp1 = qqdiv(qy, qx);
 675                 tmp2 = qatan(tmp1, epsilon);
 676                 qfree(tmp1);
 677                 return tmp2;
 678         }
 679         /*
 680          * The point is in the left half plane (x <= 0) with nonzero y.
 681          * Calculate the angle by using the formula:
 682          *      atan2(y,x) = 2 * atan(sgn(y) * sqrt((x/y)^2 + 1) - x/y).
 683          */
 684         epsilon2 = qscale(epsilon, -4L);
 685         tmp1 = qqdiv(qx, qy);
 686         tmp2 = qsquare(tmp1);
 687         tmp3 = qqadd(tmp2, &_qone_);
 688         qfree(tmp2);
 689         tmp2 = qsqrt(tmp3, epsilon2, 24L | (qy->num.sign * 64));
 690         qfree(tmp3);
 691         tmp3 = qsub(tmp2, tmp1);
 692         qfree(tmp2);
 693         qfree(tmp1);
 694         qfree(epsilon2);
 695         epsilon2 = qscale(epsilon, -1L);
 696         tmp1 = qatan(tmp3, epsilon2);
 697         qfree(epsilon2);
 698         qfree(tmp3);
 699         tmp2 = qscale(tmp1, 1L);
 700         qfree(tmp1);
 701         return tmp2;
 702 }
 703
 704
 705 /*
 706  * Calculate the value of pi to within the required epsilon.
 707  * This uses the following formula which only needs integer calculations
 708  * except for the final operation:
 709  *      pi = 1 / SUMOF(comb(2 * N, N) ^ 3 * (42 * N + 5) / 2 ^ (12 * N + 4)),
 710  * where the summation runs from N=0.  This formula gives about 6 bits of
 711  * accuracy per term.  Since the denominator for each term is a power of two,
 712  * we can simply use shifts to sum the terms.  The combinatorial numbers
 713  * in the formula are calculated recursively using the formula:
 714  *      comb(2*(N+1), N+1) = 2 * comb(2 * N, N) * (2 * N + 1) / N.
 715  */
 716 NUMBER *
 717 qpi(NUMBER *epsilon)
 718 {
 719         ZVALUE comb;                    /* current combinatorial value */
 720         ZVALUE sum;                     /* current sum */
 721         ZVALUE tmp1, tmp2;
 722         NUMBER *r, *t1, qtmp;
 723         long shift;                     /* current shift of result */
 724         long N;                         /* current term number */
 725         long bits;                      /* needed number of bits of precision */
 726         long t;
 727
 728         if (qiszero(epsilon)) {
 729                 math_error("zero epsilon value for pi");
 730                 /*NOTREACHED*/
 731         }
 732         if (epsilon == pivalue[0])
 733                 return qlink(pivalue[1]);
 734         if (pivalue[0]) {
 735                 qfree(pivalue[0]);
 736                 qfree(pivalue[1]);
 737         }
 738         bits = -qilog2(epsilon) + 4;
 739         if (bits < 4)
 740                 bits = 4;
 741         comb = _one_;
 742         itoz(5L, &sum);
 743         N = 0;
 744         shift = 4;
 745         do {
 746                 t = 1 + (++N & 0x1);
 747                 (void) zdivi(comb, N / (3 - t), &tmp1);
 748                 zfree(comb);
 749                 zmuli(tmp1, t * (2 * N - 1), &comb);
 750                 zfree(tmp1);
 751                 zsquare(comb, &tmp1);
 752                 zmul(comb, tmp1, &tmp2);
 753                 zfree(tmp1);
 754                 zmuli(tmp2, 42 * N + 5, &tmp1);
 755                 zfree(tmp2);
 756                 zshift(sum, 12L, &tmp2);
 757                 zfree(sum);
 758                 zadd(tmp1, tmp2, &sum);
 759                 t = zhighbit(tmp1);
 760                 zfree(tmp1);
 761                 zfree(tmp2);
 762                 shift += 12;
 763         } while ((shift - t) < bits);
 764         zfree(comb);
 765         qtmp.num = _one_;
 766         qtmp.den = sum;
 767         t1 = qscale(&qtmp, shift);
 768         zfree(sum);
 769         r = qmappr(t1, epsilon, 24L);
 770         qfree(t1);
 771         pivalue[0] = qlink(epsilon);
 772         pivalue[1] = qlink(r);
 773         return r;
 774 }
 775
 776 /*
 777  * Calculate the exponential function to the nearest or next to nearest
 778  * multiple of the positive number epsilon.
 779  */
 780 NUMBER *
 781 qexp(NUMBER *q, NUMBER *epsilon)
 782 {
 783         long m, n;
 784         NUMBER *tmp1, *tmp2;
 785
 786         if (qiszero(epsilon)) {
 787                 math_error("Zero epsilon value for exp");
 788                 /*NOTREACHED*/
 789         }
 790         if (qiszero(q))
 791                 return qlink(&_qone_);
 792         tmp1 = qmul(q, &_qlge_);
 793         m = qtoi(tmp1);         /* exp(q) < 2^(m+1) or m == MAXLONG */
 794         qfree(tmp1);
 795
 796         if (m > (1 << 30))
 797                 return NULL;
 798
 799         n = qilog2(epsilon);    /* 2^n <= epsilon < 2^(n+1) */
 800         if (m < n)
 801                 return qlink(&_qzero_);
 802         tmp1 = qqabs(q);
 803         tmp2 = qexprel(tmp1, m - n + 1);
 804         qfree(tmp1);
 805         if (tmp2 == NULL)
 806                 return NULL;
 807         if (qisneg(q)) {
 808                 tmp1 = qinv(tmp2);
 809                 qfree(tmp2);
 810                 tmp2 = tmp1;
 811         }
 812         tmp1 = qmappr(tmp2, epsilon, 24L);
 813         qfree(tmp2);
 814         return tmp1;
 815 }
 816
 817
 818 /*
 819  * Calculate the exponential function with relative error corresponding
 820  * to a specified number of significant bits
 821  * Requires *q >= 0, bitnum >= 0.
 822  * This returns NULL if more than 2^30 working bits would be required.
 823  */
 824 S_FUNC NUMBER *
 825 qexprel(NUMBER *q, long bitnum)
 826 {
 827         long n, m, k, h, s, t, d;
 828         NUMBER *qtmp1;
 829         ZVALUE X, B, sum, term, ztmp1, ztmp2;
 830
 831         h = qilog2(q);
 832         k = bitnum + h + 1;
 833         if (k < 0)
 834                 return qlink(&_qone_);
 835         s = k;
 836         if (k) {
 837                 do {
 838                         t = s;
 839                         s = (s + k/s)/2;
 840                 }
 841                 while (t > s);
 842         }               /* s is int(sqrt(k)) */
 843         s++;
 844         if (s < -h)
 845                 s = -h;
 846         n = h + s;      /* n is number of squarings that will be required */
 847         m = bitnum + n;
 848         if (m > (1 << 30))
 849                 return NULL;
 850         while (s > 0) {         /* increasing m by ilog2(s) */
 851                 s >>= 1;
 852                 m++;
 853         }                       /* m is working number of bits */
 854         qtmp1 = qscale(q, m - n);
 855         zquo(qtmp1->num, qtmp1->den, &X, 24);
 856         qfree(qtmp1);
 857         if (ziszero(X)) {
 858                 zfree(X);
 859                 return qlink(&_qone_);
 860         }
 861         zbitvalue(m, &sum);
 862         zcopy(X, &term);
 863         d = 1;
 864         do {
 865                 zadd(sum, term, &ztmp1);
 866                 zfree(sum);
 867                 sum = ztmp1;
 868                 zmul(term, X, &ztmp1);
 869                 zfree(term);
 870                 zshift(ztmp1, -m, &ztmp2);
 871                 zfree(ztmp1);
 872                 zdivi(ztmp2, ++d, &term);
 873                 zfree(ztmp2);
 874         }
 875         while (!ziszero(term));
 876         zfree(term);
 877         zfree(X);
 878         k = 0;
 879         zbitvalue(2 * m + 1, &B);
 880         while (n-- > 0) {
 881                 k *= 2;
 882                 zsquare(sum, &ztmp1);
 883                 zfree(sum);
 884                 if (zrel(ztmp1, B) >= 0) {
 885                         zshift(ztmp1, -m - 1, &sum);
 886                         k++;
 887                 } else {
 888                         zshift(ztmp1, -m, &sum);
 889                 }
 890                 zfree(ztmp1);
 891         }
 892         zfree(B);
 893         h = zlowbit(sum);
 894         qtmp1 = qalloc();
 895         if (m > h + k) {
 896                 zshift(sum, -h, &qtmp1->num);
 897                 zbitvalue(m - h - k, &qtmp1->den);
 898         }
 899         else
 900                 zshift(sum, k - m, &qtmp1->num);
 901         zfree(sum);
 902         return qtmp1;
 903 }
 904
 905
 906 /*
 907  * Calculate the natural logarithm of a number accurate to the specified
 908  * positive epsilon.
 909  */
 910 NUMBER *
 911 qln(NUMBER *q, NUMBER *epsilon)
 912 {
 913         long m, n, k, h, d;
 914         ZVALUE term, sum, mul, pow, X, D, B, ztmp;
 915         NUMBER *qtmp, *res;
 916         BOOL neg;
 917
 918         if (qiszero(q) || qiszero(epsilon)) {
 919                 math_error("logarithm of 0");
 920                 /*NOTREACHED*/
 921         }
 922         if (qisunit(q))
 923                 return qlink(&_qzero_);
 924         q = qqabs(q);                   /* Ignore sign of q */
 925         neg = (zrel(q->num, q->den) < 0);
 926         if (neg) {
 927                 qtmp = qinv(q);
 928                 qfree(q);
 929                 q = qtmp;
 930         }
 931         k = qilog2(q);
 932         m = -qilog2(epsilon);           /* m will be number of working bits */
 933         if (m < 0)
 934                 m = 0;
 935         h = k;
 936         while (h > 0) {
 937                 h /= 2;
 938                 m++;                    /* Add 1 for each sqrt until X < 2 */
 939         }
 940         m += 18;        /* 8 more sqrts, 8 for rounding, 2 for epsilon/4 */
 941         qtmp = qscale(q, m - k);
 942         zquo(qtmp->num, qtmp->den, &X, 24L);
 943         qfree(q);
 944         qfree(qtmp);
 945
 946         zbitvalue(m, &D);               /* Now "q" = X/D */
 947         zbitvalue(m - 8, &ztmp);
 948         zadd(D, ztmp, &B);              /* Will take sqrts until X <= B */
 949         zfree(ztmp);
 950
 951         n = 1;                  /* n is to count 1 + number of sqrts */
 952
 953         while (k > 0 || zrel(X, B) > 0) {
 954                 n++;
 955                 zshift(X, m + (k & 1), &ztmp);
 956                 zfree(X);
 957                 zsqrt(ztmp, &X, 24);
 958                 zfree(ztmp)
 959                 k /= 2;
 960         }
 961         zfree(B);
 962         zsub(X, D, &pow);       /* pow, mul used as tmps */
 963         zadd(X, D, &mul);
 964         zfree(X);
 965         zfree(D);
 966         zshift(pow, m, &ztmp);
 967         zfree(pow);
 968         zquo(ztmp, mul, &pow, 24);      /* pow now (X - D)/(X + D) */
 969         zfree(ztmp);
 970         zfree(mul);
 971
 972         zcopy(pow, &sum);       /* pow is first term of sum */
 973         zsquare(pow, &ztmp);
 974         zshift(ztmp, -m, &mul); /* mul is now multiplier for powers */
 975         zfree(ztmp);
 976
 977         d = 1;
 978         for (;;) {
 979                 zmul(pow, mul, &ztmp);
 980                 zfree(pow);
 981                 zshift(ztmp, -m, &pow);
 982                 zfree(ztmp);
 983                 d += 2;
 984                 zdivi(pow, d, &term);   /* Round down div should be round off */
 985                 if (ziszero(term)) {
 986                         zfree(term);
 987                         break;
 988                 }
 989                 zadd(sum, term, &ztmp);
 990                 zfree(term);
 991                 zfree(sum);
 992                 sum = ztmp;
 993         }
 994         zfree(pow);
 995         zfree(mul);
 996         qtmp = qalloc();        /* qtmp is to be 2^n * sum / 2^m */
 997         k = zlowbit(sum);
 998         sum.sign = neg;
 999         if (k + n >= m) {
1000                 zshift(sum, n - m, &qtmp->num);
1001         } else {
1002                 if (k) {
1003                         zshift(sum, -k, &qtmp->num);
1004                         zfree(sum);
1005                 } else {
1006                         qtmp->num = sum;
1007                 }
1008                 zbitvalue(m - k - n, &qtmp->den);
1009         }
1010         res = qmappr(qtmp, epsilon, 24L);
1011         qfree(qtmp);
1012         return res;
1013 }
1014
1015
1016 /*
1017  * Calculate the base 10 logarithm
1018  *
1019  *      log(q) = ln(q) / ln(10)
1020  */
1021 NUMBER *
1022 qlog(NUMBER *q, NUMBER *epsilon)
1023 {
1024         int need_new_ln_10 = TRUE;      /* FALSE => use cached ln_10 value */
1025         NUMBER *ln_q;                   /* ln(x) */
1026         NUMBER *ret;                    /* base 10 logarithm of x */
1027
1028         /* firewall */
1029         if (qiszero(q) || qiszero(epsilon)) {
1030                 math_error("logarithm of 0");
1031                 /*NOTREACHED*/
1032         }
1033
1034         /*
1035          * shortcut for small integer powers of 10
1036          */
1037         if (qisint(q) && qispos(q) && !zge8192b(q->num) && ziseven(q->num)) {
1038                 BOOL is_10_power;       /* TRUE ==> q is a power of 10 */
1039                 long ilog_10;           /* integer log base 10 */
1040
1041                 /* try for a quick small power of 10 log */
1042                 ilog_10 = zlog10(q->num, &is_10_power );
1043                 if (is_10_power == TRUE) {
1044                         /* is small power of 10, return log */
1045                         return itoq(ilog_10);
1046                 }
1047                 /* q is an even integer that is not a power of 10 */
1048         }
1049
1050         /*
1051          * compute ln(c) first
1052          */
1053         ln_q = qln(q, epsilon);
1054         /* log(1) == 0 */
1055         if (qiszero(ln_q)) {
1056                 return ln_q;
1057         }
1058
1059         /*
1060          * save epsilon for ln(10) if needed
1061          */
1062         if (ln_10_epsilon == NULL) {
1063                 /* first time call */
1064                 ln_10_epsilon = qcopy(epsilon);
1065         } else if (qcmp(ln_10_epsilon, epsilon) == TRUE) {
1066                 /* replaced cacheed value with epsilon arg */
1067                 qfree(ln_10_epsilon);
1068                 ln_10_epsilon = qcopy(epsilon);
1069         } else if (ln_10 != NULL) {
1070                 /* the previously computed ln(2) is OK to use */
1071                 need_new_ln_10 = FALSE;
1072         }
1073
1074         /*
1075          * compute ln(10) if needed
1076          */
1077         if (need_new_ln_10 == TRUE) {
1078                 if (ln_10 != NULL) {
1079                         qfree(ln_10);
1080                 }
1081                 ln_10 = qln(&_qten_, ln_10_epsilon);
1082         }
1083
1084         /*
1085          * return ln(q) / ln(10)
1086          */
1087         ret = qqdiv(ln_q, ln_10);
1088         qfree(ln_q);
1089         return ret;
1090 }
1091
1092
1093 /*
1094  * Calculate the result of raising one number to the power of another.
1095  * The result is calculated to the nearest or next to nearest multiple of
1096  * epsilon.
1097  */
1098 NUMBER *
1099 qpower(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
1100 {
1101         NUMBER *tmp1, *tmp2, *epsilon2;
1102         NUMBER *q1tmp, *q2tmp;
1103         long m, n;
1104
1105         if (qiszero(epsilon)) {
1106                 math_error("Zero epsilon for power");
1107                 /*NOTREACHED*/
1108         }
1109         if (qiszero(q1) && qisneg(q2)) {
1110                 math_error("Negative power of zero");
1111                 /*NOTREACHED*/
1112         }
1113         if (qiszero(q2) || qisone(q1))
1114                 return qlink(&_qone_);
1115         if (qiszero(q1))
1116                 return qlink(&_qzero_);
1117         if (qisneg(q1)) {
1118                 math_error("Negative base for qpower");
1119                 /*NOTREACHED*/
1120         }
1121         if (qisone(q2))
1122                 return qmappr(q1, epsilon, 24);
1123         if (zrel(q1->num, q1->den) < 0) {
1124                 q1tmp = qinv(q1);
1125                 q2tmp = qneg(q2);
1126         } else {
1127                 q1tmp = qlink(q1);
1128                 q2tmp = qlink(q2);
1129         }
1130         if (qisone(q2tmp)) {
1131                 qfree(q2tmp);
1132                 q2tmp = qmappr(q1tmp, epsilon, 24);
1133                 qfree(q1tmp);
1134                 return q2tmp;
1135         }
1136         m = qilog2(q1tmp);
1137         n = qilog2(epsilon);
1138         if (qisneg(q2tmp)) {
1139                 if (m > 0) {
1140                         tmp1 = itoq(m);
1141                         tmp2 = qmul(tmp1, q2tmp);
1142                         m = qtoi(tmp2);
1143                 } else {
1144                         tmp1 = qdec(q1tmp);
1145                         tmp2 = qqdiv(tmp1, q1tmp);
1146                         qfree(tmp1);
1147                         tmp1 = qmul(tmp2, q2tmp);
1148                         qfree(tmp2);
1149                         tmp2 = qmul(tmp1, &_qlge_);
1150                         m = qtoi(tmp2);
1151                 }
1152         } else {
1153                 if (m > 0) {
1154                         tmp1 = itoq(m + 1);
1155                         tmp2 = qmul(tmp1, q2tmp);
1156                         m = qtoi(tmp2);
1157                 } else {
1158                         tmp1 = qdec(q1tmp);
1159                         tmp2 = qmul(tmp1, q2tmp);
1160                         qfree(tmp1);
1161                         tmp1 = qmul(tmp2, &_qlge_);
1162                         m = qtoi(tmp1);
1163                 }
1164         }
1165         qfree(tmp1);
1166         qfree(tmp2);
1167         if (m > (1 << 30)) {
1168                 qfree(q1tmp);
1169                 qfree(q2tmp);
1170                 return NULL;
1171         }
1172         m += 1;
1173         if (m < n) {
1174                 qfree(q1tmp);
1175                 qfree(q2tmp);
1176                 return qlink(&_qzero_);
1177         }
1178         tmp1 = qqdiv(epsilon, q2tmp);
1179         tmp2 = qscale(tmp1, -m - 4);
1180         epsilon2 = qqabs(tmp2);
1181         qfree(tmp1);
1182         qfree(tmp2);
1183         tmp1 = qln(q1tmp, epsilon2);
1184         qfree(epsilon2);
1185         tmp2 = qmul(tmp1, q2tmp);
1186         qfree(tmp1);
1187         qfree(q1tmp);
1188         qfree(q2tmp);
1189         if (qisneg(tmp2)) {
1190                 tmp1 = qneg(tmp2);
1191                 qfree(tmp2);
1192                 tmp2 = qexprel(tmp1, m - n + 3);
1193                 if (tmp2 == NULL) {
1194                         qfree(tmp1);
1195                         return NULL;
1196                 }
1197                 qfree(tmp1);
1198                 tmp1 = qinv(tmp2);
1199         } else {
1200                 tmp1 = qexprel(tmp2, m - n + 3) ;
1201                 if (tmp1 == NULL) {
1202                         qfree(tmp2);
1203                         return NULL;
1204                 }
1205         }
1206         qfree(tmp2);
1207         tmp2 = qmappr(tmp1, epsilon, 24L);
1208         qfree(tmp1);
1209         return tmp2;
1210 }
1211
1212
1213 /*
1214  * Calculate the Kth root of a number to within the specified accuracy.
1215  */
1216 NUMBER *
1217 qroot(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
1218 {
1219         NUMBER *tmp1, *tmp2;
1220         int neg;
1221
1222         if (qiszero(epsilon)) {
1223                 math_error("Zero epsilon for root");
1224                 /*NOTREACHED*/
1225         }
1226         if (qisneg(q2) || qiszero(q2) || qisfrac(q2)) {
1227                 math_error("Taking bad root of number");
1228                 /*NOTREACHED*/
1229         }
1230         if (qiszero(q1) || qisone(q1) || qisone(q2))
1231                 return qlink(q1);
1232         if (qistwo(q2))
1233                 return qsqrt(q1, epsilon, 24L);
1234         neg = qisneg(q1);
1235         if (neg) {
1236                 if (ziseven(q2->num)) {
1237                         math_error("Taking even root of negative number");
1238                         /*NOTREACHED*/
1239                 }
1240                 q1 = qqabs(q1);
1241         }
1242         tmp2 = qinv(q2);
1243         tmp1 = qpower(q1, tmp2, epsilon);
1244         qfree(tmp2);
1245         if (tmp1 == NULL)
1246                 return NULL;
1247         if (neg) {
1248                 tmp2 = qneg(tmp1);
1249                 qfree(tmp1);
1250                 tmp1 = tmp2;
1251         }
1252         return tmp1;
1253 }
1254
1255
1256 /* Calculate the hyperbolic cosine function to the nearest or next to
1257  * nearest multiple of epsilon.
1258  * This is calculated using cosh(x) = (exp(x) + 1/exp(x))/2;
1259  */
1260 NUMBER *
1261 qcosh(NUMBER *q, NUMBER *epsilon)
1262 {
1263         NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
1264
1265         epsilon1 = qscale(epsilon, -2);
1266         tmp1 = qqabs(q);
1267         tmp2 = qexp(tmp1, epsilon1);
1268         qfree(tmp1);
1269         qfree(epsilon1);
1270         if (tmp2 == NULL)
1271                 return NULL;
1272         tmp1 = qinv(tmp2);
1273         tmp3 = qqadd(tmp1, tmp2);
1274         qfree(tmp1)
1275         qfree(tmp2)
1276         tmp1 = qscale(tmp3, -1);
1277         qfree(tmp3);
1278         tmp2 = qmappr(tmp1, epsilon, 24L);
1279         qfree(tmp1);
1280         return tmp2;
1281 }
1282
1283
1284 /*
1285  * Calculate the hyperbolic sine to the nearest or next to nearest
1286  * multiple of epsilon.
1287  * This is calculated using sinh(x) = (exp(x) - 1/exp(x))/2.
1288  */
1289 NUMBER *
1290 qsinh(NUMBER *q, NUMBER *epsilon)
1291 {
1292         NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
1293
1294         if (qiszero(q))
1295                 return qlink(&_qzero_);
1296         epsilon1 = qscale(epsilon, -3);
1297         tmp1 = qqabs(q);
1298         tmp2 = qexp(tmp1, epsilon1);
1299         qfree(tmp1);
1300         qfree(epsilon1);
1301         if (tmp2 == NULL)
1302                 return NULL;
1303         tmp1 = qinv(tmp2);
1304         tmp3 = qispos(q) ? qsub(tmp2, tmp1) : qsub(tmp1, tmp2);
1305         qfree(tmp1)
1306         qfree(tmp2)
1307         tmp1 = qscale(tmp3, -1);
1308         qfree(tmp3);
1309         tmp2 = qmappr(tmp1, epsilon, 24L);
1310         qfree(tmp1);
1311         return tmp2;
1312 }
1313
1314
1315 /*
1316  * Calculate the hyperbolic tangent to the nearest or next to nearest
1317  * multiple of epsilon.
1318  * This is calculated using the formulae:
1319  *      tanh(x) = 1 or -1 for very large abs(x)
1320  *      tanh(x) = (+ or -)(1 - 2 * exp(2 * x))  for large abx(x)
1321  *      tanh(x) = (exp(2*x) - 1)/(exp(2*x) + 1) otherwise
1322  */
1323 NUMBER *
1324 qtanh(NUMBER *q, NUMBER *epsilon)
1325 {
1326         NUMBER *tmp1, *tmp2, *tmp3;
1327         long n, m;
1328
1329         n = qilog2(epsilon);
1330         if (n > 0 || qiszero(q))
1331                 return qlink(&_qzero_);
1332         n = -n;
1333         tmp1 = qqabs(q);
1334         tmp2 = qmul(tmp1, &_qlge_);
1335         m = qtoi(tmp2);         /* exp(|q|) < 2^(m+1) or m == MAXLONG */
1336         qfree(tmp2);
1337         if (m > 1 + n/2) {
1338                 qfree(tmp1);
1339                 return q->num.sign ? qlink(&_qnegone_) : qlink(&_qone_);
1340         }
1341         tmp2 = qscale(tmp1, 1);
1342         qfree(tmp1);
1343         tmp1 = qexprel(tmp2, 2 + n);
1344         qfree(tmp2);
1345         if (m > 1 + n/4) {
1346                 tmp2 = qqdiv(&_qtwo_, tmp1);
1347                 qfree(tmp1);
1348                 tmp1 = qsub(&_qone_, tmp2);
1349                 qfree(tmp2);
1350         } else {
1351                 tmp2 = qdec(tmp1);
1352                 tmp3 = qinc(tmp1);
1353                 qfree(tmp1);
1354                 tmp1 = qqdiv(tmp2, tmp3);
1355                 qfree(tmp2);
1356                 qfree(tmp3);
1357         }
1358         tmp2 = qmappr(tmp1, epsilon, 24L);
1359         qfree(tmp1);
1360         if (qisneg(q)) {
1361                 tmp1 = qneg(tmp2);
1362                 qfree(tmp2);
1363                 return tmp1;
1364         }
1365         return tmp2;
1366 }
1367
1368
1369 /*
1370  * Hyperbolic cotangent.
1371  * Calculated using coth(x) = 1 + 2/(exp(2*x) - 1)
1372  */
1373 NUMBER *
1374 qcoth(NUMBER *q, NUMBER *epsilon)
1375 {
1376         NUMBER *tmp1, *tmp2, *res;
1377         long n, k;
1378
1379         if (qiszero(epsilon)) {
1380                 math_error("Zero epsilon value for coth");
1381                 /*NOTREACHED*/
1382         }
1383         if (qiszero(q)) {
1384                 math_error("Zero argument for coth");
1385                 /*NOTREACHED*/
1386         }
1387         tmp1 = qscale(q, 1);
1388         tmp2 = qqabs(tmp1);
1389         qfree(tmp1);
1390         k = qilog2(tmp2);
1391         n = qilog2(epsilon);
1392         if (k > 0) {
1393                 tmp1 = qmul(&_qlge_, tmp2);
1394                 k = qtoi(tmp1);
1395                 qfree(tmp1);
1396         } else {
1397                 k = 2 * k;
1398         }
1399         k = 4 - k - n;
1400         if (k < 4)
1401                 k = 4;
1402         tmp1 = qexprel(tmp2,  k);
1403         qfree(tmp2);
1404         tmp2 = qdec(tmp1);
1405         qfree(tmp1);
1406         if (qiszero(tmp2)) {
1407                 math_error("This should not happen ????");
1408                 /*NOTREACHED*/
1409         }
1410         tmp1 = qinv(tmp2);
1411         qfree(tmp2);
1412         tmp2 = qscale(tmp1, 1);
1413         qfree(tmp1);
1414         tmp1 = qinc(tmp2);
1415         qfree(tmp2);
1416         if (qisneg(q)) {
1417                 tmp2 = qneg(tmp1);
1418                 qfree(tmp1);
1419                 tmp1 = tmp2;
1420         }
1421         res = qmappr(tmp1, epsilon, 24L);
1422         qfree(tmp1);
1423         return res;
1424 }
1425
1426
1427 NUMBER *
1428 qsech(NUMBER *q, NUMBER *epsilon)
1429 {
1430         NUMBER *tmp1, *tmp2, *tmp3, *res;
1431         long n, k;
1432
1433         if (qiszero(epsilon)) {
1434                 math_error("Zero epsilon value for sech");
1435                 /*NOTREACHED*/
1436         }
1437         if (qiszero(q))
1438                 return qlink(&_qone_);
1439
1440         tmp1 = qqabs(q);
1441         k = 0;
1442         if (zrel(tmp1->num, tmp1->den) >= 0) {
1443                 tmp2 = qmul(&_qlge_, tmp1);
1444                 k = qtoi(tmp2);
1445                 qfree(tmp2);
1446         }
1447         n = qilog2(epsilon);
1448         if (k + n > 1) {
1449                 qfree(tmp1);
1450                 return qlink(&_qzero_);
1451         }
1452         tmp2 = qexprel(tmp1, 4 - k - n);
1453         qfree(tmp1);
1454         tmp1 = qinv(tmp2);
1455         tmp3 = qqadd(tmp1, tmp2);
1456         qfree(tmp1);
1457         qfree(tmp2);
1458         tmp1 = qinv(tmp3);
1459         qfree(tmp3);
1460         tmp2 = qscale(tmp1, 1);
1461         qfree(tmp1);
1462         res = qmappr(tmp2, epsilon, 24L);
1463         qfree(tmp2);
1464         return res;
1465 }
1466
1467
1468 NUMBER *
1469 qcsch(NUMBER *q, NUMBER *epsilon)
1470 {
1471         NUMBER *tmp1, *tmp2, *tmp3, *res;
1472         long n, k;
1473
1474         if (qiszero(epsilon)) {
1475                 math_error("Zero epsilon value for csch");
1476                 /*NOTREACHED*/
1477         }
1478         if (qiszero(q)) {
1479                 math_error("Zero argument for csch");
1480                 /*NOTREACHED*/
1481         }
1482
1483         n = qilog2(epsilon);
1484         tmp1 = qqabs(q);
1485         if (zrel(tmp1->num, tmp1->den) >= 0) {
1486                 tmp2 = qmul(&_qlge_, tmp1);
1487                 k = qtoi(tmp2);
1488                 qfree(tmp2);
1489         } else {
1490                 k = 2 * qilog2(tmp1);
1491         }
1492         if (k + n >= 1) {
1493                 qfree(tmp1);
1494                 return qlink(&_qzero_);
1495         }
1496         tmp2 = qexprel(tmp1, 4 - k - n);
1497         qfree(tmp1);
1498         tmp1 = qinv(tmp2);
1499         if (qisneg(q))
1500                 tmp3 = qsub(tmp1, tmp2);
1501         else
1502                 tmp3 = qsub(tmp2, tmp1);
1503         qfree(tmp1);
1504         qfree(tmp2);
1505         tmp1 = qinv(tmp3);
1506         qfree(tmp3)
1507         tmp2 = qscale(tmp1, 1);
1508         qfree(tmp1);
1509         res = qmappr(tmp2, epsilon, 24L);
1510         qfree(tmp2);
1511         return res;
1512 }
1513
1514
1515 /*
1516  * Compute the hyperbolic arccosine within the specified accuracy.
1517  * This is calculated using the formula:
1518  *      acosh(x) = ln(x + sqrt(x^2 - 1)).
1519  */
1520 NUMBER *
1521 qacosh(NUMBER *q, NUMBER *epsilon)
1522 {
1523         NUMBER *tmp1, *tmp2, *epsilon1;
1524         long n;
1525
1526         if (qiszero(epsilon)) {
1527                 math_error("Zero epsilon value for acosh");
1528                 /*NOTREACHED*/
1529         }
1530         if (qisone(q))
1531                 return qlink(&_qzero_);
1532         if (zrel(q->num, q->den) < 0)
1533                 return NULL;
1534         n = qilog2(epsilon);
1535         epsilon1 = qbitvalue(n - 3);
1536         tmp1 = qsquare(q);
1537         tmp2 = qdec(tmp1);
1538         qfree(tmp1);
1539         tmp1 = qsqrt(tmp2, epsilon1, 24L);
1540         qfree(tmp2);
1541         tmp2 = qqadd(tmp1, q);
1542         qfree(tmp1);
1543         tmp1 = qln(tmp2, epsilon1);
1544         qfree(tmp2);
1545         qfree(epsilon1);
1546         tmp2 = qmappr(tmp1, epsilon, 24L);
1547         qfree(tmp1);
1548         return tmp2;
1549 }
1550
1551
1552 /*
1553  * Compute the hyperbolic arcsine within the specified accuracy.
1554  * This is calculated using the formula:
1555  *      asinh(x) = ln(x + sqrt(x^2 + 1)).
1556  */
1557 NUMBER *
1558 qasinh(NUMBER *q, NUMBER *epsilon)
1559 {
1560         NUMBER *tmp1, *tmp2, *epsilon1;
1561         long n;
1562         BOOL neg;
1563
1564         if (qiszero(epsilon)) {
1565                 math_error("Zero epsilon value for asinh");
1566                 /*NOTREACHED*/
1567         }
1568         if (qiszero(q))
1569                 return qlink(&_qzero_);
1570         neg = qisneg(q);
1571         q = qqabs(q);
1572         n = qilog2(epsilon);
1573         epsilon1 = qbitvalue(n - 3);
1574         tmp1 = qsquare(q);
1575         tmp2 = qinc(tmp1);
1576         qfree(tmp1);
1577         tmp1 = qsqrt(tmp2, epsilon1, 24L);
1578         qfree(tmp2);
1579         tmp2 = qqadd(tmp1, q);
1580         qfree(tmp1);
1581         tmp1 = qln(tmp2, epsilon1);
1582         qfree(tmp2);
1583         qfree(q);
1584         qfree(epsilon1);
1585         tmp2 = qmappr(tmp1, epsilon, 24L);
1586         if (neg) {
1587                 tmp1 = qneg(tmp2);
1588                 qfree(tmp2);
1589                 return tmp1;
1590         }
1591         return tmp2;
1592 }
1593
1594
1595 /*
1596  * Compute the hyperbolic arctangent within the specified accuracy.
1597  * This is calculated using the formula:
1598  *      atanh(x) = ln((1 + x) / (1 - x)) / 2.
1599  */
1600 NUMBER *
1601 qatanh(NUMBER *q, NUMBER *epsilon)
1602 {
1603         NUMBER *tmp1, *tmp2, *tmp3, *epsilon1;
1604         ZVALUE z;
1605
1606         if (qiszero(epsilon)) {
1607                 math_error("Zero epsilon value for atanh");
1608                 /*NOTREACHED*/
1609         }
1610         if (qiszero(q))
1611                 return qlink(&_qzero_);
1612         z = q->num;
1613         z.sign = 0;
1614         if (zrel(z, q->den) >= 0)
1615                 return NULL;
1616         tmp1 = qinc(q);
1617         tmp2 = qsub(&_qone_, q);
1618         tmp3 = qqdiv(tmp1, tmp2);
1619         qfree(tmp1);
1620         qfree(tmp2);
1621         epsilon1 = qscale(epsilon, 1L);
1622         tmp1 = qln(tmp3, epsilon1);
1623         qfree(tmp3);
1624         tmp2 = qscale(tmp1, -1L);
1625         qfree(tmp1);
1626         qfree(epsilon1);
1627         return tmp2;
1628 }
1629
1630
1631 /*
1632  * Inverse hyperbolic secant function
1633  */
1634 NUMBER *
1635 qasech(NUMBER *q, NUMBER *epsilon)
1636 {
1637         NUMBER *tmp, *res;
1638
1639         tmp = qinv(q);
1640         res = qacosh(tmp, epsilon);
1641         qfree(tmp);
1642         return res;
1643 }
1644
1645
1646 /*
1647  * Inverse hyperbolic cosecant function
1648  */
1649 NUMBER *
1650 qacsch(NUMBER *q, NUMBER *epsilon)
1651 {
1652         NUMBER *tmp, *res;
1653
1654         tmp = qinv(q);
1655         res = qasinh(tmp, epsilon);
1656         qfree(tmp);
1657         return res;
1658 }
1659
1660
1661 /*
1662  * Inverse hyperbolic cotangent function
1663  */
1664 NUMBER *
1665 qacoth(NUMBER *q, NUMBER *epsilon)
1666 {
1667         NUMBER *tmp, *res;
1668
1669         tmp = qinv(q);
1670         res = qatanh(tmp, epsilon);
1671         qfree(tmp);
1672         return res;
1673 }