2 * comfunc - extended precision complex arithmetic non-primitive routines
4 * Copyright (C) 1999-2007 David I. Bell and Ernest Bowen
6 * Primary author: David I. Bell
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.
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.
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.
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 $
26 * Under source code control: 1990/02/15 01:48:13
27 * File existed as early as: before 1990
29 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
37 * cache the natural logarithm of 10
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 };
47 * Compute the result of raising a complex number to an integer power.
50 * c complex number to be raised
51 * q power to raise it to
54 c_powi(COMPLEX
*c
, NUMBER
*q
)
56 COMPLEX
*tmp
, *res
; /* temporary values */
57 long power
; /* power to raise to */
58 FULL bit
; /* current bit value */
62 math_error("Raising number to non-integral power");
66 math_error("Raising number to very large power");
69 power
= ztolong(q
->num
);
70 if (ciszero(c
) && (power
== 0)) {
71 math_error("Raising zero to zeroth power");
78 * Handle some low powers specially
81 switch ((int) (power
* sign
)) {
83 return clink(&_cone_
);
108 * Compute the power by squaring and multiplying.
109 * This uses the left to right method of power raising.
112 while ((bit
& power
) == 0)
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().
147 c_sqrt(COMPLEX
*c
, NUMBER
*epsilon
, long R
)
150 NUMBER
*es
, *aes
, *bes
, *u
, *v
, qtemp
;
152 ZVALUE g
, a
, b
, d
, aa
, cc
;
153 ZVALUE tmp1
, tmp2
, tmp3
, mul1
, mul2
;
154 long s1
, s2
, s3
, up1
, up2
;
161 if (!qisneg(c
->real
)) {
163 r
->real
= qsqrt(c
->real
, epsilon
, R
);
166 ntmp
= qneg(c
->real
);
168 r
->imag
= qsqrt(ntmp
, epsilon
, R
);
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
);
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
);
187 zequo(aes
->den
, g
, &tmp1
);
188 zmul(bes
->num
, tmp1
, &b
);
196 zadd(tmp1
, tmp2
, &tmp3
);
200 zshift(tmp3
, 4, &tmp1
);
204 s1
= zsqrt(tmp1
, &cc
, 16);
207 if (s1
== 0 && R
& 32) {
208 zmul(tmp1
, d
, &tmp2
);
210 s2
= zsqrt(tmp2
, &tmp3
, 16);
215 zreduce(tmp3
, tmp1
, &aes
->num
, &aes
->den
);
223 qtemp
.num
.sign
= sign
;
225 r
->real
= qmul(&qtemp
, epsilon
);
227 bes
= qscale(r
->real
, 1);
229 qtemp
.num
.sign
= sign
^ imsign
;
231 r
->imag
= qqdiv(c
->imag
, &qtemp
);
235 s3
= zquo(tmp3
, d
, &tmp1
, s2
< 0);
237 s2
= zquo(tmp1
, d
, &tmp3
, s1
? (s1
< 0) : 16);
239 s3
= zsqrt(tmp3
,&tmp1
,(s1
||s2
) ? (s1
<0 || s2
<0) : 16);
242 zshift(tmp1
, -1, &mul1
);
249 s2
= zquo(tmp1
, d
, &tmp2
, s1
? (s1
< 0) : 16);
251 s3
= zsqrt(tmp2
, &tmp1
, (s1
|| s2
) ? (s1
< 0 || s2
< 0) : 16);
253 zshift(tmp1
, -1, &mul2
);
261 s1
= zsqrt(tmp3
, &cc
, 0);
264 if (s1
== 0 && R
& 32) {
265 zmul(tmp1
, d
, &tmp2
);
267 s2
= zsqrt(tmp2
, &tmp3
, 0);
271 zreduce(tmp3
, d
, &aes
->num
, &aes
->den
);
278 qtemp
.num
.sign
= sign
;
280 r
->real
= qmul(&qtemp
, epsilon
);
282 bes
= qscale(r
->real
, 1);
284 qtemp
.num
.sign
= sign
^ imsign
;
286 r
->imag
= qqdiv(c
->imag
, &qtemp
);
290 s3
= zquo(tmp3
, d
, &mul1
, 0);
292 s2
= zquo(tmp1
, d
, &tmp3
, 0);
294 s3
= zsqrt(tmp3
, &mul1
, 0);
296 up1
= (s1
+ s2
+ s3
) ? 0 : -1;
299 s2
= zquo(tmp1
, d
, &tmp2
, 0);
301 s3
= zsqrt(tmp2
, &mul2
, 0);
302 up2
= (s1
+ s2
+ s3
) ? 0 : -1;
309 up1
= (long)((R
^ *mul1
.v
) & 1);
311 up1
= (R
^ epsilon
->num
.sign
^ sign
) & 1;
313 up1
^= epsilon
->num
.sign
^ sign
;
315 up1
^= epsilon
->num
.sign
;
319 up2
= (long)((R
^ *mul2
.v
) & 1);
321 up2
= (R
^ epsilon
->num
.sign
^ sign
^ imsign
) & 1;
323 up2
^= epsilon
->num
.sign
^ imsign
^ sign
;
325 up2
^= epsilon
->num
.sign
;
328 zadd(mul1
, _one_
, &tmp1
);
333 zadd(mul2
, _one_
, &tmp2
);
340 mul1
.sign
= sign
^ epsilon
->num
.sign
;
342 zreduce(mul1
, epsilon
->den
, &tmp2
, &u
->den
);
343 zmul(tmp2
, epsilon
->num
, &u
->num
);
350 mul2
.sign
= imsign
^ sign
^ epsilon
->num
.sign
;
352 zreduce(mul2
, epsilon
->den
, &tmp2
, &v
->den
);
353 zmul(tmp2
, epsilon
->num
, &v
->num
);
357 if (qiszero(u
) && qiszero(v
)) {
360 return clink(&_czero_
);
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.
376 c_root(COMPLEX
*c
, NUMBER
*q
, NUMBER
*epsilon
)
379 NUMBER
*a2pb2
, *root
, *tmp1
, *tmp2
, *epsilon2
;
382 if (qisneg(q
) || qiszero(q
) || qisfrac(q
)) {
383 math_error("Taking bad root of complex number");
386 if (cisone(c
) || qisone(q
))
389 return c_sqrt(c
, epsilon
, 24L);
390 if (cisreal(c
) && !qisneg(c
->real
)) {
391 tmp1
= qroot(c
->real
, q
, epsilon
);
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).
405 epsilon2
= qbitvalue(n
- 4);
406 tmp1
= qsquare(c
->real
);
407 tmp2
= qsquare(c
->imag
);
408 a2pb2
= qqadd(tmp1
, tmp2
);
411 tmp1
= qscale(q
, 1L);
412 root
= qroot(a2pb2
, tmp1
, epsilon2
);
421 return clink(&_czero_
);
423 epsilon2
= qbitvalue(n
- m
- 4);
424 tmp1
= qatan2(c
->imag
, c
->real
, epsilon2
);
426 tmp2
= qqdiv(tmp1
, q
);
428 r
= c_polar(root
, tmp2
, epsilon
);
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)).
441 c_exp(COMPLEX
*c
, NUMBER
*epsilon
)
444 NUMBER
*sin
, *cos
, *tmp1
, *tmp2
, *epsilon1
;
447 if (qiszero(epsilon
)) {
448 math_error("Zero epsilon for cexp");
452 tmp1
= qexp(c
->real
, epsilon
);
457 r
->real
= qexp(c
->real
, epsilon
);
461 epsilon1
= qbitvalue(n
- 2);
462 tmp1
= qexp(c
->real
, epsilon1
);
468 return clink(&_czero_
);
470 k
= qilog2(tmp1
) + 1;
473 return clink(&_czero_
);
475 qsincos(c
->imag
, k
- n
+ 2, &sin
, &cos
);
476 tmp2
= qmul(tmp1
, cos
);
480 r
->real
= qmappr(tmp2
, epsilon
, 24L);
482 tmp2
= qmul(tmp1
, sin
);
486 r
->imag
= qmappr(tmp2
, epsilon
, 24L);
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).
498 c_ln(COMPLEX
*c
, NUMBER
*epsilon
)
501 NUMBER
*a2b2
, *tmp1
, *tmp2
, *epsilon1
;
504 math_error("logarithm of zero");
508 return clink(&_czero_
);
510 if (cisreal(c
) && !qisneg(c
->real
)) {
512 r
->real
= qln(c
->real
, epsilon
);
515 tmp1
= qsquare(c
->real
);
516 tmp2
= qsquare(c
->imag
);
517 a2b2
= qqadd(tmp1
, tmp2
);
520 epsilon1
= qscale(epsilon
, 1L);
521 tmp1
= qln(a2b2
, epsilon1
);
525 r
->real
= qscale(tmp1
, -1L);
528 r
->imag
= qatan2(c
->imag
, c
->real
, epsilon
);
533 * Calculate base 10 logarithm by:
535 * log(c) = ln(c) / ln(10)
538 c_log(COMPLEX
*c
, NUMBER
*epsilon
)
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 */
545 * compute ln(c) first
547 ln_c
= c_ln(c
, epsilon
);
554 * save epsilon for ln(10) if needed
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
;
569 * compute ln(10) if needed
571 if (need_new_cln_10
== TRUE
) {
572 if (cln_10
!= NULL
) {
575 cln_10
= c_ln(&_cten_
, cln_10_epsilon
);
579 * return ln(c) / ln(10)
581 ret
= c_div(ln_c
, cln_10
);
588 * Calculate the complex cosine within the specified accuracy.
589 * This uses the formula:
590 * cos(x) = (exp(1i * x) + exp(-1i * x))/2;
593 c_cos(COMPLEX
*c
, NUMBER
*epsilon
)
595 COMPLEX
*r
, *ctmp1
, *ctmp2
, *ctmp3
;
600 if (qiszero(epsilon
)) {
601 math_error("Zero epsilon for ccos");
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
);
617 if (ciszero(ctmp2
)) {
619 return clink(&_czero_
);
621 ctmp1
= c_inv(ctmp2
);
622 ctmp3
= c_add(ctmp2
, ctmp1
);
625 ctmp1
= c_scale(ctmp3
, -1);
629 r
->real
= qmappr(ctmp1
->real
, epsilon
, 24L);
631 r
->imag
= qmappr(ctmp1
->imag
, epsilon
, 24L);
638 * Calculate the complex sine within the specified accuracy.
639 * This uses the formula:
640 * sin(x) = (exp(1i * x) - exp(-i1*x))/(2i).
643 c_sin(COMPLEX
*c
, NUMBER
*epsilon
)
645 COMPLEX
*r
, *ctmp1
, *ctmp2
, *ctmp3
;
646 NUMBER
*qtmp
, *epsilon1
;
650 if (qiszero(epsilon
)) {
651 math_error("Zero epsilon for csin");
655 return clink(&_czero_
);
658 neg
= qisneg(c
->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
);
669 if (ciszero(ctmp2
)) {
671 return clink(&_czero_
);
673 ctmp1
= c_inv(ctmp2
);
674 ctmp3
= c_sub(ctmp2
, ctmp1
);
677 ctmp1
= c_scale(ctmp3
, -1);
680 qtmp
= neg
? qlink(ctmp1
->imag
) : qneg(ctmp1
->imag
);
682 r
->real
= qmappr(qtmp
, epsilon
, 24L);
684 qtmp
= neg
? qneg(ctmp1
->real
) : qlink(ctmp1
->real
);
686 r
->imag
= qmappr(qtmp
, epsilon
, 24L);
694 c_cosh(COMPLEX
*c
, NUMBER
*epsilon
)
696 COMPLEX
*tmp1
, *tmp2
, *tmp3
;
698 tmp1
= c_exp(c
, epsilon
);
702 tmp3
= c_exp(tmp2
, epsilon
);
706 tmp2
= c_add(tmp1
, tmp3
);
709 tmp1
= c_scale(tmp2
, -1);
716 c_sinh(COMPLEX
*c
, NUMBER
*epsilon
)
718 COMPLEX
*tmp1
, *tmp2
, *tmp3
;
720 tmp1
= c_exp(c
, epsilon
);
724 tmp3
= c_exp(tmp2
, epsilon
);
728 tmp2
= c_sub(tmp1
, tmp3
);
731 tmp1
= c_scale(tmp2
, -1);
738 c_asin(COMPLEX
*c
, NUMBER
*epsilon
)
740 COMPLEX
*tmp1
, *tmp2
;
742 tmp1
= c_mul(&_conei_
, c
);
743 tmp2
= c_asinh(tmp1
, epsilon
);
745 tmp1
= c_div(tmp2
, &_conei_
);
752 c_acos(COMPLEX
*c
, NUMBER
*epsilon
)
754 COMPLEX
*tmp1
, *tmp2
;
757 tmp2
= c_sub(&_cone_
, tmp1
);
759 tmp1
= c_sqrt(tmp2
, epsilon
, 24);
761 tmp2
= c_mul(&_conei_
, tmp1
);
763 tmp1
= c_add(c
, tmp2
);
765 tmp2
= c_ln(tmp1
, epsilon
);
767 tmp1
= c_div(tmp2
, &_conei_
);
774 c_asinh(COMPLEX
*c
, NUMBER
*epsilon
)
776 COMPLEX
*tmp1
, *tmp2
, *tmp3
;
779 neg
= qisneg(c
->real
);
780 tmp1
= neg
? c_neg(c
) : clink(c
);
781 tmp2
= c_square(tmp1
);
782 tmp3
= c_add(&_cone_
, tmp2
);
784 tmp2
= c_sqrt(tmp3
, epsilon
, 24);
786 tmp3
= c_add(tmp2
, tmp1
);
789 tmp1
= c_ln(tmp3
, epsilon
);
801 c_acosh(COMPLEX
*c
, NUMBER
*epsilon
)
803 COMPLEX
*tmp1
, *tmp2
;
806 tmp2
= c_sub(tmp1
, &_cone_
);
808 tmp1
= c_sqrt(tmp2
, epsilon
, 24);
810 tmp2
= c_add(c
, tmp1
);
812 tmp1
= c_ln(tmp2
, epsilon
);
819 c_atan(COMPLEX
*c
, NUMBER
*epsilon
)
821 COMPLEX
*tmp1
, *tmp2
, *tmp3
;
823 if (qiszero(c
->real
) && qisunit(c
->imag
))
825 tmp1
= c_sub(&_conei_
, c
);
826 tmp2
= c_add(&_conei_
, c
);
827 tmp3
= c_div(tmp1
, tmp2
);
830 tmp1
= c_ln(tmp3
, epsilon
);
832 tmp2
= c_scale(tmp1
, -1);
834 tmp1
= c_div(tmp2
, &_conei_
);
841 c_acot(COMPLEX
*c
, NUMBER
*epsilon
)
843 COMPLEX
*tmp1
, *tmp2
, *tmp3
;
845 if (qiszero(c
->real
) && qisunit(c
->imag
))
847 tmp1
= c_add(c
, &_conei_
);
848 tmp2
= c_sub(c
, &_conei_
);
849 tmp3
= c_div(tmp1
, tmp2
);
852 tmp1
= c_ln(tmp3
, epsilon
);
854 tmp2
= c_scale(tmp1
, -1);
856 tmp1
= c_div(tmp2
, &_conei_
);
862 c_asec(COMPLEX
*c
, NUMBER
*epsilon
)
864 COMPLEX
*tmp1
, *tmp2
;
867 tmp2
= c_acos(tmp1
, epsilon
);
873 c_acsc(COMPLEX
*c
, NUMBER
*epsilon
)
875 COMPLEX
*tmp1
, *tmp2
;
878 tmp2
= c_asin(tmp1
, epsilon
);
885 c_atanh(COMPLEX
*c
, NUMBER
*epsilon
)
887 COMPLEX
*tmp1
, *tmp2
, *tmp3
;
889 if (qiszero(c
->imag
) && qisunit(c
->real
))
891 tmp1
= c_add(&_cone_
, c
);
892 tmp2
= c_sub(&_cone_
, c
);
893 tmp3
= c_div(tmp1
, tmp2
);
896 tmp1
= c_ln(tmp3
, epsilon
);
898 tmp2
= c_scale(tmp1
, -1);
905 c_acoth(COMPLEX
*c
, NUMBER
*epsilon
)
907 COMPLEX
*tmp1
, *tmp2
, *tmp3
;
909 if (qiszero(c
->imag
) && qisunit(c
->real
))
911 tmp1
= c_add(c
, &_cone_
);
912 tmp2
= c_sub(c
, &_cone_
);
913 tmp3
= c_div(tmp1
, tmp2
);
916 tmp1
= c_ln(tmp3
, epsilon
);
918 tmp2
= c_scale(tmp1
, -1);
924 c_asech(COMPLEX
*c
, NUMBER
*epsilon
)
926 COMPLEX
*tmp1
, *tmp2
;
929 tmp2
= c_acosh(tmp1
, epsilon
);
935 c_acsch(COMPLEX
*c
, NUMBER
*epsilon
)
937 COMPLEX
*tmp1
, *tmp2
;
940 tmp2
= c_asinh(tmp1
, epsilon
);
947 c_gd(COMPLEX
*c
, NUMBER
*epsilon
)
949 COMPLEX
*tmp1
, *tmp2
, *tmp3
;
957 q1
= qscale(c
->real
, -1);
958 eps
= qscale(epsilon
, -1);
966 tmp1
->real
= qscale(q1
, 1);
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
) {
979 q1
= neg
? qsub(&_qone_
, sin
) : qqadd(&_qone_
, sin
);
984 qsincos(c
->imag
, n
+ n1
, &sin
, &cos
);
985 q1
= neg
? qsub(&_qone_
, sin
) : qqadd(&_qone_
, sin
);
990 q1
= qln(q2
, epsilon
);
1003 if (qisneg(c
->imag
)) {
1013 neg
= qisneg(c
->real
);
1014 tmp1
= neg
? c_neg(c
) : clink(c
);
1015 tmp2
= c_exp(tmp1
, epsilon
);
1019 tmp1
= c_mul(&_conei_
, tmp2
);
1020 tmp3
= c_add(&_conei_
, tmp2
);
1022 tmp2
= c_add(tmp1
, &_cone_
);
1024 if (ciszero(tmp2
) || ciszero(tmp3
)) {
1029 tmp1
= c_div(tmp2
, tmp3
);
1032 tmp2
= c_ln(tmp1
, epsilon
);
1034 tmp1
= c_div(tmp2
, &_conei_
);
1046 c_agd(COMPLEX
*c
, NUMBER
*epsilon
)
1048 COMPLEX
*tmp1
, *tmp2
;
1050 tmp1
= c_mul(&_conei_
, c
);
1051 tmp2
= c_gd(tmp1
, epsilon
);
1055 tmp1
= c_div(tmp2
, &_conei_
);
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.
1067 c_polar(NUMBER
*q1
, NUMBER
*q2
, NUMBER
*epsilon
)
1070 NUMBER
*tmp
, *cos
, *sin
;
1073 if (qiszero(epsilon
)) {
1074 math_error("Zero epsilon for cpolar");
1078 return clink(&_czero_
);
1080 n
= qilog2(epsilon
);
1082 return clink(&_czero_
);
1086 r
->real
= qlink(q1
);
1089 qsincos(q2
, m
- n
+ 2, &sin
, &cos
);
1090 tmp
= qmul(q1
, cos
);
1093 r
->real
= qmappr(tmp
, epsilon
, 24L);
1095 tmp
= qmul(q1
, sin
);
1098 r
->imag
= qmappr(tmp
, epsilon
, 24L);
1105 * Raise one complex number to the power of another one to within the
1109 c_power(COMPLEX
*c1
, COMPLEX
*c2
, NUMBER
*epsilon
)
1111 COMPLEX
*ctmp1
, *ctmp2
;
1112 long k1
, k2
, k
, m1
, m2
, m
, n
;
1113 NUMBER
*a2b2
, *qtmp1
, *qtmp2
, *epsilon1
;
1115 if (qiszero(epsilon
)) {
1116 math_error("Zero epsilon for cpower");
1120 if (cisreal(c2
) && qisneg(c2
->real
)) {
1121 math_error ("Non-positive real exponent of zero");
1124 return clink(&_czero_
);
1126 n
= qilog2(epsilon
);
1129 if (!qiszero(c2
->real
)) {
1130 qtmp1
= qsquare(c1
->real
);
1131 qtmp2
= qsquare(c1
->imag
);
1132 a2b2
= qqadd(qtmp1
, qtmp2
);
1135 m1
= qilog2(c2
->real
);
1136 epsilon1
= qbitvalue(-m1
- 1);
1137 qtmp1
= qln(a2b2
, epsilon1
);
1140 qtmp2
= qmul(qtmp1
, c2
->real
);
1142 qtmp1
= qmul(qtmp2
, &_qlge_
);
1147 if (!qiszero(c2
->imag
)) {
1148 m2
= qilog2(c2
->imag
);
1149 epsilon1
= qbitvalue(-m2
- 1);
1150 qtmp1
= qatan2(c1
->imag
, c1
->real
, epsilon1
);
1152 qtmp2
= qmul(qtmp1
, c2
->imag
);
1154 qtmp1
= qscale(qtmp2
, -1);
1156 qtmp2
= qmul(qtmp1
, &_qlge_
);
1161 m
= (m2
> m1
) ? m2
: m1
;
1164 return clink(&_czero_
);
1165 epsilon1
= qbitvalue(n
- k
- m
- 2);
1166 ctmp1
= c_ln(c1
, epsilon1
);
1168 ctmp2
= c_mul(ctmp1
, c2
);
1170 ctmp1
= c_exp(ctmp2
, epsilon
);
1177 * Print a complex number in the current output mode.
1180 comprint(COMPLEX
*c
)
1184 if (conf
->outmode
== MODE_FRAC
) {
1188 if (!qiszero(c
->real
) || qiszero(c
->imag
))
1189 qprintnum(c
->real
, MODE_DEFAULT
);
1193 if (!qiszero(c
->real
) && !qisneg(&qtmp
))
1195 if (qisneg(&qtmp
)) {
1199 qprintnum(&qtmp
, MODE_DEFAULT
);
1205 * Print a complex number in rational representation.
1209 cprintfr(COMPLEX
*c
)
1216 if (!qiszero(r
) || qiszero(i
))
1217 qprintfr(r
, 0L, FALSE
);
1220 if (!qiszero(r
) && !qisneg(i
))
1222 zprintval(i
->num
, 0L, 0L);
1226 zprintval(i
->den
, 0L, 0L);
1232 c_ilog(COMPLEX
*c
, ZVALUE base
)
1236 qr
= qilog(c
->real
, base
);
1237 qi
= qilog(c
->imag
, base
);
1246 if (qrel(qr
, qi
) >= 0) {