1 /* mpc_log10 -- Take the base-10 logarithm of a complex number.
3 Copyright (C) 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 <limits.h> /* for CHAR_BIT */
24 /* Auxiliary functions which implement Ziv's strategy for special cases.
25 if flag = 0: compute only real part
26 if flag = 1: compute only imaginary
27 Exact cases should be dealt with separately. */
29 mpc_log10_aux (mpc_ptr rop
, mpc_srcptr op
, mpc_rnd_t rnd
, int flag
, int nb
)
31 mp_prec_t prec
= (MPFR_PREC_MIN
> 4) ? MPFR_PREC_MIN
: 4;
36 prec
= mpfr_get_prec ((flag
== 0) ? mpc_realref (rop
) : mpc_imagref (rop
));
38 mpc_init2 (tmp
, prec
);
39 mpfr_init2 (log10
, prec
);
42 mpfr_set_ui (log10
, 10, GMP_RNDN
); /* exact since prec >= 4 */
43 mpfr_log (log10
, log10
, GMP_RNDN
);
44 /* In each case we have two roundings, thus the final value is
45 x * (1+u)^2 where x is the exact value, and |u| <= 2^(-prec-1).
46 Thus the error is always less than 3 ulps. */
49 case 0: /* imag <- atan2(y/x) */
50 mpfr_atan2 (mpc_imagref (tmp
), mpc_imagref (op
), mpc_realref (op
),
52 mpfr_div (mpc_imagref (tmp
), mpc_imagref (tmp
), log10
, GMP_RNDN
);
53 ok
= mpfr_can_round (mpc_imagref (tmp
), prec
- 2, GMP_RNDN
,
54 GMP_RNDZ
, MPC_PREC_IM(rop
) +
55 (MPC_RND_IM (rnd
) == GMP_RNDN
));
57 ret
= mpfr_set (mpc_imagref (rop
), mpc_imagref (tmp
),
60 case 1: /* real <- log(x) */
61 mpfr_log (mpc_realref (tmp
), mpc_realref (op
), MPC_RND_RE (rnd
));
62 mpfr_div (mpc_realref (tmp
), mpc_realref (tmp
), log10
, GMP_RNDN
);
63 ok
= mpfr_can_round (mpc_realref (tmp
), prec
- 2, GMP_RNDN
,
64 GMP_RNDZ
, MPC_PREC_RE(rop
) +
65 (MPC_RND_RE (rnd
) == GMP_RNDN
));
67 ret
= mpfr_set (mpc_realref (rop
), mpc_realref (tmp
),
70 case 2: /* imag <- pi */
71 mpfr_const_pi (mpc_imagref (tmp
), MPC_RND_IM (rnd
));
72 mpfr_div (mpc_imagref (tmp
), mpc_imagref (tmp
), log10
, GMP_RNDN
);
73 ok
= mpfr_can_round (mpc_imagref (tmp
), prec
- 2, GMP_RNDN
,
74 GMP_RNDZ
, MPC_PREC_IM(rop
) +
75 (MPC_RND_IM (rnd
) == GMP_RNDN
));
77 ret
= mpfr_set (mpc_imagref (rop
), mpc_imagref (tmp
),
82 mpc_set_prec (tmp
, prec
);
83 mpfr_set_prec (log10
, prec
);
91 mpc_log10 (mpc_ptr rop
, mpc_srcptr op
, mpc_rnd_t rnd
)
93 int ok
= 0, loops
= 0, re_cmp
, im_cmp
, inex_re
, inex_im
, negative_zero
;
100 /* special values: NaN and infinities: same as mpc_log */
101 if (!mpc_fin_p (op
)) /* real or imaginary parts are NaN or Inf */
103 if (mpfr_nan_p (mpc_realref (op
)))
105 if (mpfr_inf_p (mpc_imagref (op
)))
106 /* (NaN, Inf) -> (+Inf, NaN) */
107 mpfr_set_inf (mpc_realref (rop
), +1);
109 /* (NaN, xxx) -> (NaN, NaN) */
110 mpfr_set_nan (mpc_realref (rop
));
111 mpfr_set_nan (mpc_imagref (rop
));
112 inex_im
= 0; /* Inf/NaN is exact */
114 else if (mpfr_nan_p (mpc_imagref (op
)))
116 if (mpfr_inf_p (mpc_realref (op
)))
117 /* (Inf, NaN) -> (+Inf, NaN) */
118 mpfr_set_inf (mpc_realref (rop
), +1);
120 /* (xxx, NaN) -> (NaN, NaN) */
121 mpfr_set_nan (mpc_realref (rop
));
122 mpfr_set_nan (mpc_imagref (rop
));
123 inex_im
= 0; /* Inf/NaN is exact */
125 else /* We have an infinity in at least one part. */
127 /* (+Inf, y) -> (+Inf, 0) for finite positive-signed y */
128 if (mpfr_inf_p (mpc_realref (op
)) && mpfr_signbit (mpc_realref (op
))
129 == 0 && mpfr_number_p (mpc_imagref (op
)))
130 inex_im
= mpfr_atan2 (mpc_imagref (rop
), mpc_imagref (op
),
131 mpc_realref (op
), MPC_RND_IM (rnd
));
133 /* (xxx, Inf) -> (+Inf, atan2(Inf/xxx))
134 (Inf, yyy) -> (+Inf, atan2(yyy/Inf)) */
135 inex_im
= mpc_log10_aux (rop
, op
, rnd
, 1, 0);
136 mpfr_set_inf (mpc_realref (rop
), +1);
138 return MPC_INEX(0, inex_im
);
141 /* special cases: real and purely imaginary numbers */
142 re_cmp
= mpfr_cmp_ui (mpc_realref (op
), 0);
143 im_cmp
= mpfr_cmp_ui (mpc_imagref (op
), 0);
144 if (im_cmp
== 0) /* Im(op) = 0 */
146 if (re_cmp
== 0) /* Re(op) = 0 */
148 if (mpfr_signbit (mpc_realref (op
)) == 0)
149 inex_im
= mpfr_atan2 (mpc_imagref (rop
), mpc_imagref (op
),
150 mpc_realref (op
), MPC_RND_IM (rnd
));
152 inex_im
= mpc_log10_aux (rop
, op
, rnd
, 1, 0);
153 mpfr_set_inf (mpc_realref (rop
), -1);
154 inex_re
= 0; /* -Inf is exact */
158 inex_re
= mpfr_log10 (mpc_realref (rop
), mpc_realref (op
),
160 inex_im
= mpfr_set (mpc_imagref (rop
), mpc_imagref (op
),
163 else /* log10(x + 0*i) for negative x */
164 { /* op = x + 0*i; let w = -x = |x| */
165 negative_zero
= mpfr_signbit (mpc_imagref (op
));
167 rnd_im
= INV_RND (MPC_RND_IM (rnd
));
169 rnd_im
= MPC_RND_IM (rnd
);
170 ww
->re
[0] = *mpc_realref (op
);
171 MPFR_CHANGE_SIGN (ww
->re
);
172 ww
->im
[0] = *mpc_imagref (op
);
173 if (mpfr_cmp_ui (ww
->re
, 1) == 0)
174 inex_re
= mpfr_set_ui (mpc_realref (rop
), 0, MPC_RND_RE (rnd
));
176 inex_re
= mpc_log10_aux (rop
, ww
, rnd
, 0, 1);
177 inex_im
= mpc_log10_aux (rop
, op
, MPC_RND (0,rnd_im
), 1, 2);
180 mpc_conj (rop
, rop
, MPC_RNDNN
);
184 return MPC_INEX(inex_re
, inex_im
);
186 else if (re_cmp
== 0)
190 inex_re
= mpfr_log10 (mpc_realref (rop
), mpc_imagref (op
), MPC_RND_RE (rnd
));
191 inex_im
= mpc_log10_aux (rop
, op
, rnd
, 1, 2);
192 /* division by 2 does not change the ternary flag */
193 mpfr_div_2ui (mpc_imagref (rop
), mpc_imagref (rop
), 1, GMP_RNDN
);
197 w
[0] = *mpc_imagref (op
);
198 MPFR_CHANGE_SIGN (w
);
199 inex_re
= mpfr_log10 (mpc_realref (rop
), w
, MPC_RND_RE (rnd
));
200 invrnd
= MPC_RND (0, INV_RND (MPC_RND_IM (rnd
)));
201 inex_im
= mpc_log10_aux (rop
, op
, invrnd
, 1, 2);
202 /* division by 2 does not change the ternary flag */
203 mpfr_div_2ui (mpc_imagref (rop
), mpc_imagref (rop
), 1, GMP_RNDN
);
204 mpfr_neg (mpc_imagref (rop
), mpc_imagref (rop
), GMP_RNDN
);
205 inex_im
= -inex_im
; /* negate the ternary flag */
207 return MPC_INEX(inex_re
, inex_im
);
210 /* generic case: neither Re(op) nor Im(op) is NaN, Inf or zero */
211 prec
= MPC_PREC_RE(rop
);
212 mpfr_init2 (w
, prec
);
213 mpc_init2 (ww
, prec
);
214 /* let op = x + iy; compute log(op)/log(10) */
218 prec
+= (loops
<= 2) ? mpc_ceil_log2 (prec
) + 4 : prec
/ 2;
219 mpfr_set_prec (w
, prec
);
220 mpc_set_prec (ww
, prec
);
222 mpc_log (ww
, op
, MPC_RNDNN
);
223 mpfr_set_ui (w
, 10, GMP_RNDN
); /* exact since prec >= 4 */
224 mpfr_log (w
, w
, GMP_RNDN
);
225 mpc_div_fr (ww
, ww
, w
, MPC_RNDNN
);
227 ok
= mpfr_can_round (mpc_realref (ww
), prec
- 2, GMP_RNDN
, GMP_RNDZ
,
228 MPC_PREC_RE(rop
) + (MPC_RND_RE (rnd
) == GMP_RNDN
));
230 /* Special code to deal with cases where the real part of log10(x+i*y)
231 is exact, like x=3 and y=1. Since Re(log10(x+i*y)) = log10(x^2+y^2)/2
232 this happens whenever x^2+y^2 is a nonnegative power of 10.
233 Indeed x^2+y^2 cannot equal 10^(a/2^b) for a, b integers, a odd, b>0,
234 since x^2+y^2 is rational, and 10^(a/2^b) is irrational.
235 Similarly, for b=0, x^2+y^2 cannot equal 10^a for a < 0 since x^2+y^2
236 is a rational with denominator a power of 2.
237 Now let x^2+y^2 = 10^s. Without loss of generality we can assume
238 x = u/2^e and y = v/2^e with u, v, e integers: u^2+v^2 = 10^s*2^(2e)
239 thus u^2+v^2 = 0 mod 2^(2e). By recurrence on e, necessarily
240 u = v = 0 mod 2^e, thus x and y are necessarily integers.
242 if ((ok
== 0) && (loops
== 1) && mpfr_integer_p (mpc_realref (op
)) &&
243 mpfr_integer_p (mpc_imagref (op
)))
250 mpfr_get_z (x
, mpc_realref (op
), GMP_RNDN
); /* exact */
251 mpfr_get_z (y
, mpc_imagref (op
), GMP_RNDN
); /* exact */
254 mpz_add (x
, x
, y
); /* x^2+y^2 */
255 v
= mpz_scan1 (x
, 0);
256 /* if x = 10^s then necessarily s = v */
257 s
= mpz_sizeinbase (x
, 10);
258 /* since s is either the number of digits of x or one more,
259 then x = 10^(s-1) or 10^(s-2) */
260 if (s
== v
+ 1 || s
== v
+ 2)
262 mpz_div_2exp (x
, x
, v
);
263 mpz_ui_pow_ui (y
, 5, v
);
264 if (mpz_cmp (y
, x
) == 0) /* Re(log10(x+i*y)) is exactly v/2 */
266 /* we reset the precision of Re(ww) so that v can be
267 represented exactly */
268 mpfr_set_prec (mpc_realref (ww
), sizeof(unsigned long)*CHAR_BIT
);
269 mpfr_set_ui_2exp (mpc_realref (ww
), v
, -1, GMP_RNDN
); /* exact */
277 ok
= ok
&& mpfr_can_round (mpc_imagref (ww
), prec
-2, GMP_RNDN
, GMP_RNDZ
,
278 MPC_PREC_IM(rop
) + (MPC_RND_IM (rnd
) == GMP_RNDN
));
281 inex_re
= mpfr_set (mpc_realref(rop
), mpc_realref (ww
), MPC_RND_RE (rnd
));
282 inex_im
= mpfr_set (mpc_imagref(rop
), mpc_imagref (ww
), MPC_RND_IM (rnd
));
285 return MPC_INEX(inex_re
, inex_im
);