c_cpp/etc/calc/comfunc.c

   1 /*
   2  * comfunc - extended precision complex arithmetic non-primitive routines
   3  *
   4  * Copyright (C) 1999-2007  David I. Bell and Ernest Bowen
   5  *
   6  * Primary author:  David I. Bell
   7  *
   8  * Calc is open software; you can redistribute it and/or modify it under
   9  * the terms of the version 2.1 of the GNU Lesser General Public License
  10  * as published by the Free Software Foundation.
  11  *
  12  * Calc is distributed in the hope that it will be useful, but WITHOUT
  13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  15  * Public License for more details.
  16  *
  17  * A copy of version 2.1 of the GNU Lesser General Public License is
  18  * distributed with calc under the filename COPYING-LGPL.  You should have
  19  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  20  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  21  *
  22  * @(#) $Revision: 30.3 $
  23  * @(#) $Id: comfunc.c,v 30.3 2008/10/24 07:09:41 chongo Exp $
  24  * @(#) $Source: /usr/local/src/bin/calc/RCS/comfunc.c,v $
  25  *
  26  * Under source code control:   1990/02/15 01:48:13
  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 "config.h"
  34 #include "cmath.h"
  35
  36 /*
  37  * cache the natural logarithm of 10
  38  */
  39 STATIC COMPLEX *cln_10 = NULL;
  40 STATIC NUMBER *cln_10_epsilon = NULL;
  41 STATIC NUMBER _q10_ = { { _tenval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
  42 STATIC NUMBER _q0_ = { { _zeroval_, 1, 0 }, { _oneval_, 1, 0 }, 1, NULL };
  43 COMPLEX _cten_ = { &_q10_, &_q0_, 1 };
  44
  45
  46 /*
  47  * Compute the result of raising a complex number to an integer power.
  48  *
  49  * given:
  50  *      c               complex number to be raised
  51  *      q               power to raise it to
  52  */
  53 COMPLEX *
  54 c_powi(COMPLEX *c, NUMBER *q)
  55 {
  56         COMPLEX *tmp, *res;     /* temporary values */
  57         long power;             /* power to raise to */
  58         FULL bit;               /* current bit value */
  59         int sign;
  60
  61         if (qisfrac(q)) {
  62                 math_error("Raising number to non-integral power");
  63                 /*NOTREACHED*/
  64         }
  65         if (zge31b(q->num)) {
  66                 math_error("Raising number to very large power");
  67                 /*NOTREACHED*/
  68         }
  69         power = ztolong(q->num);
  70         if (ciszero(c) && (power == 0)) {
  71                 math_error("Raising zero to zeroth power");
  72                 /*NOTREACHED*/
  73         }
  74         sign = 1;
  75         if (qisneg(q))
  76                 sign = -1;
  77         /*
  78          * Handle some low powers specially
  79          */
  80         if (power <= 4) {
  81                 switch ((int) (power * sign)) {
  82                 case 0:
  83                         return clink(&_cone_);
  84                 case 1:
  85                         return clink(c);
  86                 case -1:
  87                         return c_inv(c);
  88                 case 2:
  89                         return c_square(c);
  90                 case -2:
  91                         tmp = c_square(c);
  92                         res = c_inv(tmp);
  93                         comfree(tmp);
  94                         return res;
  95                 case 3:
  96                         tmp = c_square(c);
  97                         res = c_mul(c, tmp);
  98                         comfree(tmp);
  99                         return res;
 100                 case 4:
 101                         tmp = c_square(c);
 102                         res = c_square(tmp);
 103                         comfree(tmp);
 104                         return res;
 105                 }
 106         }
 107         /*
 108          * Compute the power by squaring and multiplying.
 109          * This uses the left to right method of power raising.
 110          */
 111         bit = TOPFULL;
 112         while ((bit & power) == 0)
 113                 bit >>= 1L;
 114         bit >>= 1L;
 115         res = c_square(c);
 116         if (bit & power) {
 117                 tmp = c_mul(res, c);
 118                 comfree(res);
 119                 res = tmp;
 120         }
 121         bit >>= 1L;
 122         while (bit) {
 123                 tmp = c_square(res);
 124                 comfree(res);
 125                 res = tmp;
 126                 if (bit & power) {
 127                         tmp = c_mul(res, c);
 128                         comfree(res);
 129                         res = tmp;
 130                 }
 131                 bit >>= 1L;
 132         }
 133         if (sign < 0) {
 134                 tmp = c_inv(res);
 135                 comfree(res);
 136                 res = tmp;
 137         }
 138         return res;
 139 }
 140
 141
 142 /*
 143  * Calculate the square root of a complex number to specified accuracy.
 144  * Type of rounding of each component specified by R as for qsqrt().
 145  */
 146 COMPLEX *
 147 c_sqrt(COMPLEX *c, NUMBER *epsilon, long R)
 148 {
 149         COMPLEX *r;
 150         NUMBER *es, *aes, *bes, *u, *v, qtemp;
 151         NUMBER *ntmp;
 152         ZVALUE g, a, b, d, aa, cc;
 153         ZVALUE tmp1, tmp2, tmp3, mul1, mul2;
 154         long s1, s2, s3, up1, up2;
 155         int imsign, sign;
 156
 157         if (ciszero(c))
 158                 return clink(c);
 159         if (cisreal(c)) {
 160                 r = comalloc();
 161                 if (!qisneg(c->real)) {
 162                         qfree(r->real);
 163                         r->real = qsqrt(c->real, epsilon, R);
 164                         return r;
 165                 }
 166                 ntmp = qneg(c->real);
 167                 qfree(r->imag);
 168                 r->imag = qsqrt(ntmp, epsilon, R);
 169                 qfree(ntmp);
 170                 return r;
 171         }
 172
 173         up1 = up2 = 0;
 174         sign = (R & 64) != 0;
 175         imsign = c->imag->num.sign;
 176         es = qsquare(epsilon);
 177         aes = qqdiv(c->real, es);
 178         bes = qqdiv(c->imag, es);
 179         qfree(es);
 180         zgcd(aes->den, bes->den, &g);
 181         zequo(bes->den, g, &tmp1);
 182         zmul(aes->num, tmp1, &a);
 183         zmul(aes->den, tmp1, &tmp2);
 184         zshift(tmp2, 1, &d);
 185         zfree(tmp1);
 186         zfree(tmp2);
 187         zequo(aes->den, g, &tmp1);
 188         zmul(bes->num, tmp1, &b);
 189         zfree(tmp1);
 190         zfree(g);
 191         qfree(aes);
 192         qfree(bes);
 193         zsquare(a, &tmp1);
 194         zsquare(b, &tmp2);
 195         zfree(b);
 196         zadd(tmp1, tmp2, &tmp3);
 197         zfree(tmp1);
 198         zfree(tmp2);
 199         if (R & 16) {
 200                 zshift(tmp3, 4, &tmp1);
 201                 zfree(tmp3);
 202                 zshift(a, 2, &aa);
 203                 zfree(a);
 204                 s1 = zsqrt(tmp1, &cc, 16);
 205                 zfree(tmp1);
 206                 zadd(cc, aa, &tmp1);
 207                 if (s1 == 0 && R & 32) {
 208                         zmul(tmp1, d, &tmp2);
 209                         zfree(tmp1);
 210                         s2 = zsqrt(tmp2, &tmp3, 16);
 211                         zfree(tmp2);
 212                         if (s2 == 0) {
 213                                 aes = qalloc();
 214                                 zshift(d, 1, &tmp1);
 215                                 zreduce(tmp3, tmp1, &aes->num, &aes->den);
 216                                 zfree(tmp1);
 217                                 zfree(tmp3);
 218                                 zfree(aa);
 219                                 zfree(cc);
 220                                 zfree(d);
 221                                 r = comalloc();
 222                                 qtemp = *aes;
 223                                 qtemp.num.sign = sign;
 224                                 qfree(r->real);
 225                                 r->real = qmul(&qtemp, epsilon);
 226                                 qfree(aes);
 227                                 bes = qscale(r->real, 1);
 228                                 qtemp = *bes;
 229                                 qtemp.num.sign = sign ^ imsign;
 230                                 qfree(r->imag);
 231                                 r->imag = qqdiv(c->imag, &qtemp);
 232                                 qfree(bes);
 233                                 return r;
 234                         }
 235                         s3 = zquo(tmp3, d, &tmp1, s2 < 0);
 236                 } else {
 237                         s2 = zquo(tmp1, d, &tmp3, s1 ? (s1 < 0) : 16);
 238                         zfree(tmp1);
 239                         s3 = zsqrt(tmp3,&tmp1,(s1||s2) ? (s1<0 || s2<0) : 16);
 240                 }
 241                 zfree(tmp3);
 242                 zshift(tmp1, -1, &mul1);
 243                 if (*tmp1.v & 1)
 244                         up1 = s1 + s2 + s3;
 245                 else
 246                         up1 = -1;
 247                 zfree(tmp1);
 248                 zsub(cc, aa, &tmp1);
 249                 s2 = zquo(tmp1, d, &tmp2, s1 ? (s1 < 0) : 16);
 250                 zfree(tmp1);
 251                 s3 = zsqrt(tmp2, &tmp1, (s1 || s2) ? (s1 < 0 || s2 < 0) : 16);
 252                 zfree(tmp2);
 253                 zshift(tmp1, -1, &mul2);
 254                 if (*tmp1.v & 1)
 255                         up2 = s1 + s2 + s3;
 256                 else
 257                         up2 = -1;
 258                 zfree(tmp1);
 259                 zfree(aa);
 260         } else {
 261                 s1 = zsqrt(tmp3, &cc, 0);
 262                 zfree(tmp3);
 263                 zadd(cc, a, &tmp1);
 264                 if (s1 == 0 && R & 32) {
 265                         zmul(tmp1, d, &tmp2);
 266                         zfree(tmp1);
 267                         s2 = zsqrt(tmp2, &tmp3, 0);
 268                         zfree(tmp2);
 269                         if (s2 == 0) {
 270                                 aes = qalloc();
 271                                 zreduce(tmp3, d, &aes->num, &aes->den);
 272                                 zfree(tmp3);
 273                                 zfree(a);
 274                                 zfree(cc);
 275                                 zfree(d);
 276                                 r = comalloc();
 277                                 qtemp = *aes;
 278                                 qtemp.num.sign = sign;
 279                                 qfree(r->real);
 280                                 r->real = qmul(&qtemp, epsilon);
 281                                 qfree(aes);
 282                                 bes = qscale(r->real, 1);
 283                                 qtemp = *bes;
 284                                 qtemp.num.sign = sign ^ imsign;
 285                                 qfree(r->imag);
 286                                 r->imag = qqdiv(c->imag, &qtemp);
 287                                 qfree(bes);
 288                                 return r;
 289                         }
 290                         s3 = zquo(tmp3, d, &mul1, 0);
 291                 } else {
 292                         s2 = zquo(tmp1, d, &tmp3, 0);
 293                         zfree(tmp1);
 294                         s3 = zsqrt(tmp3, &mul1, 0);
 295                 }
 296                 up1 = (s1 + s2 + s3) ? 0 : -1;
 297                 zfree(tmp3);
 298                 zsub(cc, a, &tmp1);
 299                 s2 = zquo(tmp1, d, &tmp2, 0);
 300                 zfree(tmp1);
 301                 s3 = zsqrt(tmp2, &mul2, 0);
 302                 up2 = (s1 + s2 + s3) ? 0 : -1;
 303                 zfree(tmp2);
 304                 zfree(a);
 305         }
 306         zfree(cc); zfree(d);
 307         if (up1 == 0) {
 308                 if (R & 8)
 309                         up1 = (long)((R ^ *mul1.v) & 1);
 310                 else
 311                         up1 = (R ^ epsilon->num.sign ^ sign) & 1;
 312                 if (R & 2)
 313                         up1 ^= epsilon->num.sign ^ sign;
 314                 if (R & 4)
 315                         up1 ^= epsilon->num.sign;
 316         }
 317         if (up2 == 0) {
 318                 if (R & 8)
 319                         up2 = (long)((R ^ *mul2.v) & 1);
 320                 else
 321                         up2 = (R ^ epsilon->num.sign ^ sign ^ imsign) & 1;
 322                 if (R & 2)
 323                         up2 ^= epsilon->num.sign ^ imsign ^ sign;
 324                 if (R & 4)
 325                         up2 ^= epsilon->num.sign;
 326         }
 327         if (up1 > 0) {
 328                 zadd(mul1, _one_, &tmp1);
 329                 zfree(mul1);
 330                 mul1 = tmp1;
 331         }
 332         if (up2 > 0) {
 333                 zadd(mul2, _one_, &tmp2);
 334                 zfree(mul2);
 335                 mul2 = tmp2;
 336         }
 337         if (ziszero(mul1)) {
 338                 u = qlink(&_qzero_);
 339         } else {
 340                 mul1.sign = sign ^ epsilon->num.sign;
 341                 u = qalloc();
 342                 zreduce(mul1, epsilon->den, &tmp2, &u->den);
 343                 zmul(tmp2, epsilon->num, &u->num);
 344                 zfree(tmp2);
 345         }
 346         zfree(mul1);
 347         if (ziszero(mul2)) {
 348                 v = qlink(&_qzero_);
 349         } else {
 350                 mul2.sign = imsign ^ sign ^ epsilon->num.sign;
 351                 v = qalloc();
 352                 zreduce(mul2, epsilon->den, &tmp2, &v->den);
 353                 zmul(tmp2, epsilon->num, &v->num);
 354                 zfree(tmp2);
 355         }
 356         zfree(mul2);
 357         if (qiszero(u) && qiszero(v)) {
 358                 qfree(u);
 359                 qfree(v);
 360                 return clink(&_czero_);
 361         }
 362         r = comalloc();
 363         qfree(r->real);
 364         qfree(r->imag);
 365         r->real = u;
 366         r->imag = v;
 367         return r;
 368 }
 369
 370
 371 /*
 372  * Take the Nth root of a complex number, where N is a positive integer.
 373  * Each component of the result is within the specified error.
 374  */
 375 COMPLEX *
 376 c_root(COMPLEX *c, NUMBER *q, NUMBER *epsilon)
 377 {
 378         COMPLEX *r;
 379         NUMBER *a2pb2, *root, *tmp1, *tmp2, *epsilon2;
 380         long n, m;
 381
 382         if (qisneg(q) || qiszero(q) || qisfrac(q)) {
 383                 math_error("Taking bad root of complex number");
 384                 /*NOTREACHED*/
 385         }
 386         if (cisone(c) || qisone(q))
 387                 return clink(c);
 388         if (qistwo(q))
 389                 return c_sqrt(c, epsilon, 24L);
 390         if (cisreal(c) && !qisneg(c->real)) {
 391                 tmp1 = qroot(c->real, q, epsilon);
 392                 if (tmp1 == NULL)
 393                         return NULL;
 394                 r = comalloc();
 395                 qfree(r->real);
 396                 r->real = tmp1;
 397                 return r;
 398         }
 399         /*
 400          * Calculate the root using the formula:
 401          *      c_root(a + bi, n) =
 402          *              c_polar(qroot(a^2 + b^2, 2 * n), qatan2(b, a) / n).
 403          */
 404         n = qilog2(epsilon);
 405         epsilon2 = qbitvalue(n - 4);
 406         tmp1 = qsquare(c->real);
 407         tmp2 = qsquare(c->imag);
 408         a2pb2 = qqadd(tmp1, tmp2);
 409         qfree(tmp1);
 410         qfree(tmp2);
 411         tmp1 = qscale(q, 1L);
 412         root = qroot(a2pb2, tmp1, epsilon2);
 413         qfree(a2pb2);
 414         qfree(tmp1);
 415         qfree(epsilon2);
 416         if (root == NULL)
 417                 return NULL;
 418         m = qilog2(root);
 419         if (m < n) {
 420                 qfree(root);
 421                 return clink(&_czero_);
 422         }
 423         epsilon2 = qbitvalue(n - m - 4);
 424         tmp1 = qatan2(c->imag, c->real, epsilon2);
 425         qfree(epsilon2);
 426         tmp2 = qqdiv(tmp1, q);
 427         qfree(tmp1);
 428         r = c_polar(root, tmp2, epsilon);
 429         qfree(root);
 430         qfree(tmp2);
 431         return r;
 432 }
 433
 434
 435 /*
 436  * Calculate the complex exponential function to the desired accuracy.
 437  * We use the formula:
 438  *      exp(a + bi) = exp(a) * (cos(b) + i * sin(b)).
 439  */
 440 COMPLEX *
 441 c_exp(COMPLEX *c, NUMBER *epsilon)
 442 {
 443         COMPLEX *r;
 444         NUMBER *sin, *cos, *tmp1, *tmp2, *epsilon1;
 445         long k, n;
 446
 447         if (qiszero(epsilon)) {
 448                 math_error("Zero epsilon for cexp");
 449                 /*NOTREACHED*/
 450         }
 451         if (cisreal(c)) {
 452                 tmp1 = qexp(c->real, epsilon);
 453                 if (tmp1 == NULL)
 454                         return NULL;
 455                 r = comalloc();
 456                 qfree(r->real);
 457                 r->real = qexp(c->real, epsilon);
 458                 return r;
 459         }
 460         n = qilog2(epsilon);
 461         epsilon1 = qbitvalue(n - 2);
 462         tmp1 = qexp(c->real, epsilon1);
 463         qfree(epsilon1);
 464         if (tmp1 == NULL)
 465                 return NULL;
 466         if (qiszero(tmp1)) {
 467                 qfree(tmp1);
 468                 return clink(&_czero_);
 469         }
 470         k = qilog2(tmp1) + 1;
 471         if (k < n) {
 472                 qfree(tmp1);
 473                 return clink(&_czero_);
 474         }
 475         qsincos(c->imag, k - n + 2, &sin, &cos);
 476         tmp2 = qmul(tmp1, cos);
 477         qfree(cos);
 478         r = comalloc();
 479         qfree(r->real);
 480         r->real = qmappr(tmp2, epsilon, 24L);
 481         qfree(tmp2);
 482         tmp2 = qmul(tmp1, sin);
 483         qfree(tmp1);
 484         qfree(sin);
 485         qfree(r->imag);
 486         r->imag = qmappr(tmp2, epsilon, 24L);
 487         qfree(tmp2);
 488         return r;
 489 }
 490
 491
 492 /*
 493  * Calculate the natural logarithm of a complex number within the specified
 494  * error.  We use the formula:
 495  *      ln(a + bi) = ln(a^2 + b^2) / 2 + i * atan2(b, a).
 496  */
 497 COMPLEX *
 498 c_ln(COMPLEX *c, NUMBER *epsilon)
 499 {
 500         COMPLEX *r;
 501         NUMBER *a2b2, *tmp1, *tmp2, *epsilon1;
 502
 503         if (ciszero(c)) {
 504                 math_error("logarithm of zero");
 505                 /*NOTREACHED*/
 506         }
 507         if (cisone(c))
 508                 return clink(&_czero_);
 509         r = comalloc();
 510         if (cisreal(c) && !qisneg(c->real)) {
 511                 qfree(r->real);
 512                 r->real = qln(c->real, epsilon);
 513                 return r;
 514         }
 515         tmp1 = qsquare(c->real);
 516         tmp2 = qsquare(c->imag);
 517         a2b2 = qqadd(tmp1, tmp2);
 518         qfree(tmp1);
 519         qfree(tmp2);
 520         epsilon1 = qscale(epsilon, 1L);
 521         tmp1 = qln(a2b2, epsilon1);
 522         qfree(a2b2);
 523         qfree(epsilon1);
 524         qfree(r->real);
 525         r->real = qscale(tmp1, -1L);
 526         qfree(tmp1);
 527         qfree(r->imag);
 528         r->imag = qatan2(c->imag, c->real, epsilon);
 529         return r;
 530 }
 531
 532 /*
 533  * Calculate base 10 logarithm by:
 534  *
 535  *      log(c) = ln(c) / ln(10)
 536  */
 537 COMPLEX *
 538 c_log(COMPLEX *c, NUMBER *epsilon)
 539 {
 540         int need_new_cln_10 = TRUE;     /* FALSE => use cached cln_10 value */
 541         COMPLEX *ln_c;                  /* ln(x) */
 542         COMPLEX *ret;                   /* base 10 logarithm of x */
 543
 544         /*
 545          * compute ln(c) first
 546          */
 547         ln_c = c_ln(c, epsilon);
 548         /* log(1) == 0 */
 549         if (ciszero(ln_c)) {
 550                 return ln_c;
 551         }
 552
 553         /*
 554          * save epsilon for ln(10) if needed
 555          */
 556         if (cln_10_epsilon == NULL) {
 557                 /* first time call */
 558                 cln_10_epsilon = qcopy(epsilon);
 559         } else if (qcmp(cln_10_epsilon, epsilon) == TRUE) {
 560                 /* replaced cacheed value with epsilon arg */
 561                 qfree(cln_10_epsilon);
 562                 cln_10_epsilon = qcopy(epsilon);
 563         } else if (cln_10 != NULL) {
 564                 /* the previously computed ln(2) is OK to use */
 565                 need_new_cln_10 = FALSE;
 566         }
 567
 568         /*
 569          * compute ln(10) if needed
 570          */
 571         if (need_new_cln_10 == TRUE) {
 572                 if (cln_10 != NULL) {
 573                         comfree(cln_10);
 574                 }
 575                 cln_10 = c_ln(&_cten_, cln_10_epsilon);
 576         }
 577
 578         /*
 579          * return ln(c) / ln(10)
 580          */
 581         ret = c_div(ln_c, cln_10);
 582         comfree(ln_c);
 583         return ret;
 584 }
 585
 586
 587 /*
 588  * Calculate the complex cosine within the specified accuracy.
 589  * This uses the formula:
 590  *      cos(x) = (exp(1i * x) + exp(-1i * x))/2;
 591  */
 592 COMPLEX *
 593 c_cos(COMPLEX *c, NUMBER *epsilon)
 594 {
 595         COMPLEX *r, *ctmp1, *ctmp2, *ctmp3;
 596         NUMBER *epsilon1;
 597         long n;
 598         BOOL neg;
 599
 600         if (qiszero(epsilon)) {
 601                 math_error("Zero epsilon for ccos");
 602                 /*NOTREACHED*/
 603         }
 604         n = qilog2(epsilon);
 605         ctmp1 = comalloc();
 606         qfree(ctmp1->real);
 607         qfree(ctmp1->imag);
 608         neg = qisneg(c->imag);
 609         ctmp1->real = neg ? qneg(c->imag) : qlink(c->imag);
 610         ctmp1->imag = neg ? qlink(c->real) : qneg(c->real);
 611         epsilon1 = qbitvalue(n - 2);
 612         ctmp2 = c_exp(ctmp1, epsilon1);
 613         comfree(ctmp1);
 614         qfree(epsilon1);
 615         if (ctmp2 == NULL)
 616                 return NULL;
 617         if (ciszero(ctmp2)) {
 618                 comfree(ctmp2);
 619                 return clink(&_czero_);
 620         }
 621         ctmp1 = c_inv(ctmp2);
 622         ctmp3 = c_add(ctmp2, ctmp1);
 623         comfree(ctmp1);
 624         comfree(ctmp2);
 625         ctmp1 = c_scale(ctmp3, -1);
 626         comfree(ctmp3);
 627         r = comalloc();
 628         qfree(r->real);
 629         r->real = qmappr(ctmp1->real, epsilon, 24L);
 630         qfree(r->imag);
 631         r->imag = qmappr(ctmp1->imag, epsilon, 24L);
 632         comfree(ctmp1);
 633         return r;
 634 }
 635
 636
 637 /*
 638  * Calculate the complex sine within the specified accuracy.
 639  * This uses the formula:
 640  *      sin(x) = (exp(1i * x) - exp(-i1*x))/(2i).
 641  */
 642 COMPLEX *
 643 c_sin(COMPLEX *c, NUMBER *epsilon)
 644 {
 645         COMPLEX *r, *ctmp1, *ctmp2, *ctmp3;
 646         NUMBER *qtmp, *epsilon1;
 647         long n;
 648         BOOL neg;
 649
 650         if (qiszero(epsilon)) {
 651                 math_error("Zero epsilon for csin");
 652                 /*NOTREACHED*/
 653         }
 654         if (ciszero(c))
 655                 return clink(&_czero_);
 656         n = qilog2(epsilon);
 657         ctmp1 = comalloc();
 658         neg = qisneg(c->imag);
 659         qfree(ctmp1->real);
 660         qfree(ctmp1->imag);
 661         ctmp1->real = neg ? qneg(c->imag) : qlink(c->imag);
 662         ctmp1->imag = neg ? qlink(c->real) : qneg(c->real);
 663         epsilon1 = qbitvalue(n - 2);
 664         ctmp2 = c_exp(ctmp1, epsilon1);
 665         comfree(ctmp1);
 666         qfree(epsilon1);
 667         if (ctmp2 == NULL)
 668                 return NULL;
 669         if (ciszero(ctmp2)) {
 670                 comfree(ctmp2);
 671                 return clink(&_czero_);
 672         }
 673         ctmp1 = c_inv(ctmp2);
 674         ctmp3 = c_sub(ctmp2, ctmp1);
 675         comfree(ctmp1);
 676         comfree(ctmp2);
 677         ctmp1 = c_scale(ctmp3, -1);
 678         comfree(ctmp3);
 679         r = comalloc();
 680         qtmp = neg ? qlink(ctmp1->imag) : qneg(ctmp1->imag);
 681         qfree(r->real);
 682         r->real = qmappr(qtmp, epsilon, 24L);
 683         qfree(qtmp);
 684         qtmp = neg ? qneg(ctmp1->real) : qlink(ctmp1->real);
 685         qfree(r->imag);
 686         r->imag = qmappr(qtmp, epsilon, 24L);
 687         qfree(qtmp);
 688         comfree(ctmp1);
 689         return r;
 690 }
 691
 692
 693 COMPLEX *
 694 c_cosh(COMPLEX *c, NUMBER *epsilon)
 695 {
 696         COMPLEX *tmp1, *tmp2, *tmp3;
 697
 698         tmp1 = c_exp(c, epsilon);
 699         if (tmp1 == NULL)
 700                 return NULL;
 701         tmp2 = c_neg(c);
 702         tmp3 = c_exp(tmp2, epsilon);
 703         comfree(tmp2);
 704         if (tmp3 == NULL)
 705                 return NULL;
 706         tmp2 = c_add(tmp1, tmp3);
 707         comfree(tmp1);
 708         comfree(tmp3);
 709         tmp1 = c_scale(tmp2, -1);
 710         comfree(tmp2);
 711         return tmp1;
 712 }
 713
 714
 715 COMPLEX *
 716 c_sinh(COMPLEX *c, NUMBER *epsilon)
 717 {
 718         COMPLEX *tmp1, *tmp2, *tmp3;
 719
 720         tmp1 = c_exp(c, epsilon);
 721         if (tmp1 == NULL)
 722                 return NULL;
 723         tmp2 = c_neg(c);
 724         tmp3 = c_exp(tmp2, epsilon);
 725         comfree(tmp2);
 726         if (tmp3 == NULL)
 727                 return NULL;
 728         tmp2 = c_sub(tmp1, tmp3);
 729         comfree(tmp1);
 730         comfree(tmp3);
 731         tmp1 = c_scale(tmp2, -1);
 732         comfree(tmp2);
 733         return tmp1;
 734 }
 735
 736
 737 COMPLEX *
 738 c_asin(COMPLEX *c, NUMBER *epsilon)
 739 {
 740         COMPLEX *tmp1, *tmp2;
 741
 742         tmp1 = c_mul(&_conei_, c);
 743         tmp2 = c_asinh(tmp1, epsilon);
 744         comfree(tmp1);
 745         tmp1 = c_div(tmp2, &_conei_);
 746         comfree(tmp2);
 747         return tmp1;
 748 }
 749
 750
 751 COMPLEX *
 752 c_acos(COMPLEX *c, NUMBER *epsilon)
 753 {
 754         COMPLEX *tmp1, *tmp2;
 755
 756         tmp1 = c_square(c);
 757         tmp2 = c_sub(&_cone_, tmp1);
 758         comfree(tmp1);
 759         tmp1 = c_sqrt(tmp2, epsilon, 24);
 760         comfree(tmp2);
 761         tmp2 = c_mul(&_conei_, tmp1);
 762         comfree(tmp1);
 763         tmp1 = c_add(c, tmp2);
 764         comfree(tmp2);
 765         tmp2 = c_ln(tmp1, epsilon);
 766         comfree(tmp1);
 767         tmp1 = c_div(tmp2, &_conei_);
 768         comfree(tmp2);
 769         return tmp1;
 770 }
 771
 772
 773 COMPLEX *
 774 c_asinh(COMPLEX *c, NUMBER *epsilon)
 775 {
 776         COMPLEX *tmp1, *tmp2, *tmp3;
 777         BOOL neg;
 778
 779         neg = qisneg(c->real);
 780         tmp1 = neg ? c_neg(c) : clink(c);
 781         tmp2 = c_square(tmp1);
 782         tmp3 = c_add(&_cone_, tmp2);
 783         comfree(tmp2);
 784         tmp2 = c_sqrt(tmp3, epsilon, 24);
 785         comfree(tmp3);
 786         tmp3 = c_add(tmp2, tmp1);
 787         comfree(tmp1);
 788         comfree(tmp2);
 789         tmp1 = c_ln(tmp3, epsilon);
 790         comfree(tmp3);
 791         if (neg) {
 792                 tmp2 = c_neg(tmp1);
 793                 comfree(tmp1);
 794                 return tmp2;
 795         }
 796         return tmp1;
 797 }
 798
 799
 800 COMPLEX *
 801 c_acosh(COMPLEX *c, NUMBER *epsilon)
 802 {
 803         COMPLEX *tmp1, *tmp2;
 804
 805         tmp1 = c_square(c);
 806         tmp2 = c_sub(tmp1, &_cone_);
 807         comfree(tmp1);
 808         tmp1 = c_sqrt(tmp2, epsilon, 24);
 809         comfree(tmp2);
 810         tmp2 = c_add(c, tmp1);
 811         comfree(tmp1);
 812         tmp1 = c_ln(tmp2, epsilon);
 813         comfree(tmp2);
 814         return tmp1;
 815 }
 816
 817
 818 COMPLEX *
 819 c_atan(COMPLEX *c, NUMBER *epsilon)
 820 {
 821         COMPLEX *tmp1, *tmp2, *tmp3;
 822
 823         if (qiszero(c->real) && qisunit(c->imag))
 824                 return NULL;
 825         tmp1 = c_sub(&_conei_, c);
 826         tmp2 = c_add(&_conei_, c);
 827         tmp3 = c_div(tmp1, tmp2);
 828         comfree(tmp1);
 829         comfree(tmp2);
 830         tmp1 = c_ln(tmp3, epsilon);
 831         comfree(tmp3);
 832         tmp2 = c_scale(tmp1, -1);
 833         comfree(tmp1);
 834         tmp1 = c_div(tmp2, &_conei_);
 835         comfree(tmp2);
 836         return tmp1;
 837 }
 838
 839
 840 COMPLEX *
 841 c_acot(COMPLEX *c, NUMBER *epsilon)
 842 {
 843         COMPLEX *tmp1, *tmp2, *tmp3;
 844
 845         if (qiszero(c->real) && qisunit(c->imag))
 846                 return NULL;
 847         tmp1 = c_add(c, &_conei_);
 848         tmp2 = c_sub(c, &_conei_);
 849         tmp3 = c_div(tmp1, tmp2);
 850         comfree(tmp1);
 851         comfree(tmp2);
 852         tmp1 = c_ln(tmp3, epsilon);
 853         comfree(tmp3);
 854         tmp2 = c_scale(tmp1, -1);
 855         comfree(tmp1);
 856         tmp1 = c_div(tmp2, &_conei_);
 857         comfree(tmp2);
 858         return tmp1;
 859 }
 860
 861 COMPLEX *
 862 c_asec(COMPLEX *c, NUMBER *epsilon)
 863 {
 864         COMPLEX *tmp1, *tmp2;
 865
 866         tmp1 = c_inv(c);
 867         tmp2 = c_acos(tmp1, epsilon);
 868         comfree(tmp1);
 869         return tmp2;
 870 }
 871
 872 COMPLEX *
 873 c_acsc(COMPLEX *c, NUMBER *epsilon)
 874 {
 875         COMPLEX *tmp1, *tmp2;
 876
 877         tmp1 = c_inv(c);
 878         tmp2 = c_asin(tmp1, epsilon);
 879         comfree(tmp1);
 880         return tmp2;
 881 }
 882
 883
 884 COMPLEX *
 885 c_atanh(COMPLEX *c, NUMBER *epsilon)
 886 {
 887         COMPLEX *tmp1, *tmp2, *tmp3;
 888
 889         if (qiszero(c->imag) && qisunit(c->real))
 890                 return NULL;
 891         tmp1 = c_add(&_cone_, c);
 892         tmp2 = c_sub(&_cone_, c);
 893         tmp3 = c_div(tmp1, tmp2);
 894         comfree(tmp1);
 895         comfree(tmp2);
 896         tmp1 = c_ln(tmp3, epsilon);
 897         comfree(tmp3);
 898         tmp2 = c_scale(tmp1, -1);
 899         comfree(tmp1);
 900         return tmp2;
 901 }
 902
 903
 904 COMPLEX *
 905 c_acoth(COMPLEX *c, NUMBER *epsilon)
 906 {
 907         COMPLEX *tmp1, *tmp2, *tmp3;
 908
 909         if (qiszero(c->imag) && qisunit(c->real))
 910                 return NULL;
 911         tmp1 = c_add(c, &_cone_);
 912         tmp2 = c_sub(c, &_cone_);
 913         tmp3 = c_div(tmp1, tmp2);
 914         comfree(tmp1);
 915         comfree(tmp2);
 916         tmp1 = c_ln(tmp3, epsilon);
 917         comfree(tmp3);
 918         tmp2 = c_scale(tmp1, -1);
 919         comfree(tmp1);
 920         return tmp2;
 921 }
 922
 923 COMPLEX *
 924 c_asech(COMPLEX *c, NUMBER *epsilon)
 925 {
 926         COMPLEX *tmp1, *tmp2;
 927
 928         tmp1 = c_inv(c);
 929         tmp2 = c_acosh(tmp1, epsilon);
 930         comfree(tmp1);
 931         return tmp2;
 932 }
 933
 934 COMPLEX *
 935 c_acsch(COMPLEX *c, NUMBER *epsilon)
 936 {
 937         COMPLEX *tmp1, *tmp2;
 938
 939         tmp1 = c_inv(c);
 940         tmp2 = c_asinh(tmp1, epsilon);
 941         comfree(tmp1);
 942         return tmp2;
 943 }
 944
 945
 946 COMPLEX *
 947 c_gd(COMPLEX *c, NUMBER *epsilon)
 948 {
 949         COMPLEX *tmp1, *tmp2, *tmp3;
 950         NUMBER *q1, *q2;
 951         NUMBER *sin, *cos;
 952         NUMBER *eps;
 953         int n, n1;
 954         BOOL neg;
 955
 956         if (cisreal(c)) {
 957                 q1 = qscale(c->real, -1);
 958                 eps = qscale(epsilon, -1);
 959                 q2 = qtanh(q1, eps);
 960                 qfree(q1);
 961                 q1 = qatan(q2, eps);
 962                 qfree(eps);
 963                 qfree(q2);
 964                 tmp1 = comalloc();
 965                 qfree(tmp1->real);
 966                 tmp1->real = qscale(q1, 1);
 967                 qfree(q1);
 968                 return tmp1;
 969         }
 970         if (qiszero(c->real)) {
 971                 n = - qilog2(epsilon);
 972                 qsincos(c->imag, n + 8, &sin, &cos);
 973                 if (qiszero(cos) || (n1 = -qilog2(cos)) > n) {
 974                         qfree(sin);
 975                         qfree(cos);
 976                         return NULL;
 977                 }
 978                 neg = qisneg(sin);
 979                 q1 = neg ? qsub(&_qone_, sin) : qqadd(&_qone_, sin);
 980                 qfree(sin);
 981                 if (n1 > 8) {
 982                         qfree(q1);
 983                         qfree(cos);
 984                         qsincos(c->imag, n + n1, &sin, &cos);
 985                         q1 = neg ? qsub(&_qone_, sin) : qqadd(&_qone_, sin);
 986                         qfree(sin);
 987                 }
 988                 q2 = qqdiv(q1, cos);
 989                 qfree(q1);
 990                 q1 = qln(q2, epsilon);
 991                 qfree(q2);
 992                 if (neg) {
 993                         q2 = qneg(q1);
 994                         qfree(q1);
 995                         q1 = q2;
 996                 }
 997                 tmp1 = comalloc();
 998                 qfree(tmp1->imag);
 999                 tmp1->imag = q1;
1000                 if (qisneg(cos)) {
1001                         qfree(tmp1->real);
1002                         q1 = qpi(epsilon);
1003                         if (qisneg(c->imag)) {
1004                                 q2 = qneg(q1);
1005                                 qfree(q1);
1006                                 q1 = q2;
1007                         }
1008                         tmp1->real = q1;
1009                 }
1010                 qfree(cos);
1011                 return tmp1;
1012         }
1013         neg = qisneg(c->real);
1014         tmp1 = neg ? c_neg(c) : clink(c);
1015         tmp2 = c_exp(tmp1, epsilon);
1016         comfree(tmp1);
1017         if (tmp2 == NULL)
1018                 return NULL;
1019         tmp1 = c_mul(&_conei_, tmp2);
1020         tmp3 = c_add(&_conei_, tmp2);
1021         comfree(tmp2);
1022         tmp2 = c_add(tmp1, &_cone_);
1023         comfree(tmp1);
1024         if (ciszero(tmp2) || ciszero(tmp3)) {
1025                 comfree(tmp2);
1026                 comfree(tmp3);
1027                 return NULL;
1028         }
1029         tmp1 = c_div(tmp2, tmp3);
1030         comfree(tmp2);
1031         comfree(tmp3);
1032         tmp2 = c_ln(tmp1, epsilon);
1033         comfree(tmp1);
1034         tmp1 = c_div(tmp2, &_conei_);
1035         comfree(tmp2);
1036         if (neg) {
1037                 tmp2 = c_neg(tmp1);
1038                 comfree(tmp1);
1039                 return tmp2;
1040         }
1041         return tmp1;
1042 }
1043
1044
1045 COMPLEX *
1046 c_agd(COMPLEX *c, NUMBER *epsilon)
1047 {
1048         COMPLEX *tmp1, *tmp2;
1049
1050         tmp1 = c_mul(&_conei_, c);
1051         tmp2 = c_gd(tmp1, epsilon);
1052         comfree(tmp1);
1053         if (tmp2 == NULL)
1054                 return NULL;
1055         tmp1 = c_div(tmp2, &_conei_);
1056         comfree(tmp2);
1057         return tmp1;
1058 }
1059
1060
1061 /*
1062  * Convert a number from polar coordinates to normal complex number form
1063  * within the specified accuracy.  This produces the value:
1064  *      q1 * cos(q2) + q1 * sin(q2) * i.
1065  */
1066 COMPLEX *
1067 c_polar(NUMBER *q1, NUMBER *q2, NUMBER *epsilon)
1068 {
1069         COMPLEX *r;
1070         NUMBER *tmp, *cos, *sin;
1071         long m, n;
1072
1073         if (qiszero(epsilon)) {
1074                 math_error("Zero epsilon for cpolar");
1075                 /*NOTREACHED*/
1076         }
1077         if (qiszero(q1))
1078                 return clink(&_czero_);
1079         m = qilog2(q1) + 1;
1080         n = qilog2(epsilon);
1081         if (m < n)
1082                 return clink(&_czero_);
1083         r = comalloc();
1084         if (qiszero(q2)) {
1085                 qfree(r->real);
1086                 r->real = qlink(q1);
1087                 return r;
1088         }
1089         qsincos(q2, m - n + 2, &sin, &cos);
1090         tmp = qmul(q1, cos);
1091         qfree(cos);
1092         qfree(r->real);
1093         r->real = qmappr(tmp, epsilon, 24L);
1094         qfree(tmp);
1095         tmp = qmul(q1, sin);
1096         qfree(sin);
1097         qfree(r->imag);
1098         r->imag = qmappr(tmp, epsilon, 24L);
1099         qfree(tmp);
1100         return r;
1101 }
1102
1103
1104 /*
1105  * Raise one complex number to the power of another one to within the
1106  * specified error.
1107  */
1108 COMPLEX *
1109 c_power(COMPLEX *c1, COMPLEX *c2, NUMBER *epsilon)
1110 {
1111         COMPLEX *ctmp1, *ctmp2;
1112         long k1, k2, k, m1, m2, m, n;
1113         NUMBER *a2b2, *qtmp1, *qtmp2, *epsilon1;
1114
1115         if (qiszero(epsilon)) {
1116                 math_error("Zero epsilon for cpower");
1117                 /*NOTREACHED*/
1118         }
1119         if (ciszero(c1)) {
1120                 if (cisreal(c2) && qisneg(c2->real)) {
1121                         math_error ("Non-positive real exponent of zero");
1122                         /*NOTREACHED*/
1123                 }
1124                 return clink(&_czero_);
1125         }
1126         n = qilog2(epsilon);
1127         m1 = m2 = -1000000;
1128         k1 = k2 = 0;
1129         if (!qiszero(c2->real)) {
1130                 qtmp1 = qsquare(c1->real);
1131                 qtmp2 = qsquare(c1->imag);
1132                 a2b2 = qqadd(qtmp1, qtmp2);
1133                 qfree(qtmp1);
1134                 qfree(qtmp2);
1135                 m1 = qilog2(c2->real);
1136                 epsilon1 = qbitvalue(-m1 - 1);
1137                 qtmp1 = qln(a2b2, epsilon1);
1138                 qfree(epsilon1);
1139                 qfree(a2b2);
1140                 qtmp2 = qmul(qtmp1, c2->real);
1141                 qfree(qtmp1);
1142                 qtmp1 = qmul(qtmp2, &_qlge_);
1143                 qfree(qtmp2);
1144                 k1 = qtoi(qtmp1);
1145                 qfree(qtmp1);
1146         }
1147         if (!qiszero(c2->imag)) {
1148                 m2 = qilog2(c2->imag);
1149                 epsilon1 = qbitvalue(-m2 - 1);
1150                 qtmp1 = qatan2(c1->imag, c1->real, epsilon1);
1151                 qfree(epsilon1);
1152                 qtmp2 = qmul(qtmp1, c2->imag);
1153                 qfree(qtmp1);
1154                 qtmp1 = qscale(qtmp2, -1);
1155                 qfree(qtmp2);
1156                 qtmp2 = qmul(qtmp1, &_qlge_);
1157                 qfree(qtmp1);
1158                 k2 = qtoi(qtmp2);
1159                 qfree(qtmp2);
1160         }
1161         m = (m2 > m1) ? m2 : m1;
1162         k = k1 - k2 + 1;
1163         if (k < n)
1164                 return clink(&_czero_);
1165         epsilon1 = qbitvalue(n - k - m - 2);
1166         ctmp1 = c_ln(c1, epsilon1);
1167         qfree(epsilon1);
1168         ctmp2 = c_mul(ctmp1, c2);
1169         comfree(ctmp1);
1170         ctmp1 = c_exp(ctmp2, epsilon);
1171         comfree(ctmp2);
1172         return ctmp1;
1173 }
1174
1175
1176 /*
1177  * Print a complex number in the current output mode.
1178  */
1179 void
1180 comprint(COMPLEX *c)
1181 {
1182         NUMBER qtmp;
1183
1184         if (conf->outmode == MODE_FRAC) {
1185                 cprintfr(c);
1186                 return;
1187         }
1188         if (!qiszero(c->real) || qiszero(c->imag))
1189                 qprintnum(c->real, MODE_DEFAULT);
1190         qtmp = c->imag[0];
1191         if (qiszero(&qtmp))
1192                 return;
1193         if (!qiszero(c->real) && !qisneg(&qtmp))
1194                 math_chr('+');
1195         if (qisneg(&qtmp)) {
1196                 math_chr('-');
1197                 qtmp.num.sign = 0;
1198         }
1199         qprintnum(&qtmp, MODE_DEFAULT);
1200         math_chr('i');
1201 }
1202
1203
1204 /*
1205  * Print a complex number in rational representation.
1206  * Example:  2/3-4i/5
1207  */
1208 void
1209 cprintfr(COMPLEX *c)
1210 {
1211         NUMBER *r;
1212         NUMBER *i;
1213
1214         r = c->real;
1215         i = c->imag;
1216         if (!qiszero(r) || qiszero(i))
1217                 qprintfr(r, 0L, FALSE);
1218         if (qiszero(i))
1219                 return;
1220         if (!qiszero(r) && !qisneg(i))
1221                 math_chr('+');
1222         zprintval(i->num, 0L, 0L);
1223         math_chr('i');
1224         if (qisfrac(i)) {
1225                 math_chr('/');
1226                 zprintval(i->den, 0L, 0L);
1227         }
1228 }
1229
1230
1231 NUMBER *
1232 c_ilog(COMPLEX *c, ZVALUE base)
1233 {
1234         NUMBER *qr, *qi;
1235
1236         qr = qilog(c->real, base);
1237         qi = qilog(c->imag, base);
1238
1239         if (qr == NULL) {
1240                 if (qi == NULL)
1241                         return NULL;
1242                 return qi;
1243         }
1244         if (qi == NULL)
1245                 return qr;
1246         if (qrel(qr, qi) >= 0) {
1247                 qfree(qi);
1248                 return qr;
1249         }
1250         qfree(qr);
1251         return qi;
1252 }