1 /* mpc_log -- Take the logarithm of a complex number.
3 Copyright (C) 2008, 2009, 2010, 2011, 2012 INRIA
5 This file is part of GNU MPC.
7 GNU MPC is free software; you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
17 You should have received a copy of the GNU Lesser General Public License
18 along with this program. If not, see http://www.gnu.org/licenses/ .
21 #include <stdio.h> /* for MPC_ASSERT */
25 mpc_log (mpc_ptr rop
, mpc_srcptr op
, mpc_rnd_t rnd
){
26 int ok
, underflow
= 0;
37 /* special values: NaN and infinities */
38 if (!mpc_fin_p (op
)) {
39 if (mpfr_nan_p (mpc_realref (op
))) {
40 if (mpfr_inf_p (mpc_imagref (op
)))
41 mpfr_set_inf (mpc_realref (rop
), +1);
43 mpfr_set_nan (mpc_realref (rop
));
44 mpfr_set_nan (mpc_imagref (rop
));
45 inex_im
= 0; /* Inf/NaN is exact */
47 else if (mpfr_nan_p (mpc_imagref (op
))) {
48 if (mpfr_inf_p (mpc_realref (op
)))
49 mpfr_set_inf (mpc_realref (rop
), +1);
51 mpfr_set_nan (mpc_realref (rop
));
52 mpfr_set_nan (mpc_imagref (rop
));
53 inex_im
= 0; /* Inf/NaN is exact */
55 else /* We have an infinity in at least one part. */ {
56 inex_im
= mpfr_atan2 (mpc_imagref (rop
), mpc_imagref (op
), mpc_realref (op
),
58 mpfr_set_inf (mpc_realref (rop
), +1);
60 return MPC_INEX(0, inex_im
);
63 /* special cases: real and purely imaginary numbers */
64 re_cmp
= mpfr_cmp_ui (mpc_realref (op
), 0);
65 im_cmp
= mpfr_cmp_ui (mpc_imagref (op
), 0);
68 inex_im
= mpfr_atan2 (mpc_imagref (rop
), mpc_imagref (op
), mpc_realref (op
),
70 mpfr_set_inf (mpc_realref (rop
), -1);
71 inex_re
= 0; /* -Inf is exact */
73 else if (re_cmp
> 0) {
74 inex_re
= mpfr_log (mpc_realref (rop
), mpc_realref (op
), MPC_RND_RE (rnd
));
75 inex_im
= mpfr_set (mpc_imagref (rop
), mpc_imagref (op
), MPC_RND_IM (rnd
));
78 /* op = x + 0*y; let w = -x = |x| */
82 negative_zero
= mpfr_signbit (mpc_imagref (op
));
84 rnd_im
= INV_RND (MPC_RND_IM (rnd
));
86 rnd_im
= MPC_RND_IM (rnd
);
87 w
[0] = *mpc_realref (op
);
89 inex_re
= mpfr_log (mpc_realref (rop
), w
, MPC_RND_RE (rnd
));
90 inex_im
= mpfr_const_pi (mpc_imagref (rop
), rnd_im
);
92 mpc_conj (rop
, rop
, MPC_RNDNN
);
96 return MPC_INEX(inex_re
, inex_im
);
98 else if (re_cmp
== 0) {
100 inex_re
= mpfr_log (mpc_realref (rop
), mpc_imagref (op
), MPC_RND_RE (rnd
));
101 inex_im
= mpfr_const_pi (mpc_imagref (rop
), MPC_RND_IM (rnd
));
102 /* division by 2 does not change the ternary flag */
103 mpfr_div_2ui (mpc_imagref (rop
), mpc_imagref (rop
), 1, GMP_RNDN
);
106 w
[0] = *mpc_imagref (op
);
107 MPFR_CHANGE_SIGN (w
);
108 inex_re
= mpfr_log (mpc_realref (rop
), w
, MPC_RND_RE (rnd
));
109 inex_im
= mpfr_const_pi (mpc_imagref (rop
), INV_RND (MPC_RND_IM (rnd
)));
110 /* division by 2 does not change the ternary flag */
111 mpfr_div_2ui (mpc_imagref (rop
), mpc_imagref (rop
), 1, GMP_RNDN
);
112 mpfr_neg (mpc_imagref (rop
), mpc_imagref (rop
), GMP_RNDN
);
113 inex_im
= -inex_im
; /* negate the ternary flag */
115 return MPC_INEX(inex_re
, inex_im
);
118 prec
= MPC_PREC_RE(rop
);
120 /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x) */
121 /* loop for the real part: 1/2 log (x^2 + y^2), fast, but unsafe */
124 for (loops
= 1; !ok
&& loops
<= 2; loops
++) {
125 prec
+= mpc_ceil_log2 (prec
) + 4;
126 mpfr_set_prec (w
, prec
);
128 mpc_abs (w
, op
, GMP_RNDN
);
131 /* intermediate overflow; the logarithm may be representable.
132 Intermediate underflow is impossible. */
135 mpfr_log (w
, w
, GMP_RNDN
);
136 /* generic error of log: (2^(- exp(w)) + 0.5) ulp */
139 /* impossible to round, switch to second algorithm */
142 err
= MPC_MAX (-mpfr_get_exp (w
), 0) + 1;
143 /* number of lost digits */
144 ok
= mpfr_can_round (w
, prec
- err
, GMP_RNDN
, GMP_RNDZ
,
145 mpfr_get_prec (mpc_realref (rop
)) + (MPC_RND_RE (rnd
) == GMP_RNDN
));
149 prec
= MPC_PREC_RE(rop
);
151 /* compute 1/2 log (x^2 + y^2) = log |x| + 1/2 * log (1 + (y/x)^2)
152 if |x| >= |y|; otherwise, exchange x and y */
153 if (mpfr_cmpabs (mpc_realref (op
), mpc_imagref (op
)) >= 0) {
154 x
= mpc_realref (op
);
155 y
= mpc_imagref (op
);
158 x
= mpc_imagref (op
);
159 y
= mpc_realref (op
);
163 prec
+= mpc_ceil_log2 (prec
) + 4;
164 mpfr_set_prec (v
, prec
);
165 mpfr_set_prec (w
, prec
);
167 mpfr_div (v
, y
, x
, GMP_RNDD
); /* error 1 ulp */
168 mpfr_sqr (v
, v
, GMP_RNDD
);
169 /* generic error of multiplication:
170 1 + 2*1*(2+1*2^(1-prec)) <= 5.0625 since prec >= 6 */
171 mpfr_log1p (v
, v
, GMP_RNDD
);
172 /* error 1 + 4*5.0625 = 21.25 , see algorithms.tex */
173 mpfr_div_2ui (v
, v
, 1, GMP_RNDD
);
174 /* If the result is 0, then there has been an underflow somewhere. */
176 mpfr_abs (w
, x
, GMP_RNDN
); /* exact */
177 mpfr_log (w
, w
, GMP_RNDN
); /* error 0.5 ulp */
178 expw
= mpfr_get_exp (w
);
179 sgnw
= mpfr_signbit (w
);
181 mpfr_add (w
, w
, v
, GMP_RNDN
);
182 if (!sgnw
) /* v is positive, so no cancellation;
183 error 22.25 ulp; error counts lost bits */
186 err
= MPC_MAX (5 + mpfr_get_exp (v
),
187 /* 21.25 ulp (v) rewritten in ulp (result, now in w) */
188 -1 + expw
- mpfr_get_exp (w
)
189 /* 0.5 ulp (previous w), rewritten in ulp (result) */
192 /* handle one special case: |x|=1, and (y/x)^2 underflows;
193 then 1/2*log(x^2+y^2) \approx 1/2*y^2 also underflows. */
194 if ( (mpfr_cmp_si (x
, -1) == 0 || mpfr_cmp_ui (x
, 1) == 0)
198 } while (!underflow
&&
199 !mpfr_can_round (w
, prec
- err
, GMP_RNDN
, GMP_RNDZ
,
200 mpfr_get_prec (mpc_realref (rop
)) + (MPC_RND_RE (rnd
) == GMP_RNDN
)));
205 inex_im
= mpfr_atan2 (mpc_imagref (rop
), mpc_imagref (op
), mpc_realref (op
),
208 /* set the real part; cannot be done before if rop==op */
210 /* create underflow in result */
211 inex_re
= mpfr_set_ui_2exp (mpc_realref (rop
), 1,
212 mpfr_get_emin_min () - 2, MPC_RND_RE (rnd
));
214 inex_re
= mpfr_set (mpc_realref (rop
), w
, MPC_RND_RE (rnd
));
216 return MPC_INEX(inex_re
, inex_im
);