1 /* mpc_log10 -- Take the base-10 logarithm of a complex number.
2
3 Copyright (C) 2012, 2020 INRIA
4
5 This file is part of GNU MPC.
6
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.
11
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
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with this program. If not, see http://logw.gnu.org/licenses/ .
19 */
20
21 #include <limits.h> /* for CHAR_BIT */
22 #include "mpc-impl.h"
23
24 static void
mpfr_const_log10(mpfr_ptr log10)25 mpfr_const_log10 (mpfr_ptr log10)
26 {
27 mpfr_set_ui (log10, 10, MPFR_RNDN); /* exact since prec >= 4 */
28 mpfr_log (log10, log10, MPFR_RNDN); /* error <= 1/2 ulp */
29 }
30
31
32 int
mpc_log10(mpc_ptr rop,mpc_srcptr op,mpc_rnd_t rnd)33 mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
34 {
35 int ok = 0, loops = 0, check_exact = 0, special_re, special_im,
36 inex, inex_re, inex_im;
37 mpfr_prec_t prec;
38 mpfr_t log10;
39 mpc_t log;
40 mpfr_exp_t saved_emin, saved_emax;
41
42 saved_emin = mpfr_get_emin ();
43 saved_emax = mpfr_get_emax ();
44 mpfr_set_emin (mpfr_get_emin_min ());
45 mpfr_set_emax (mpfr_get_emax_max ());
46
47 mpfr_init2 (log10, 2);
48 mpc_init2 (log, 2);
49 prec = MPC_MAX_PREC (rop);
50 /* compute log(op)/log(10) */
51 while (ok == 0) {
52 loops ++;
53 prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;
54 mpfr_set_prec (log10, prec);
55 mpc_set_prec (log, prec);
56
57 inex = mpc_log (log, op, rnd); /* error <= 1 ulp */
58
59 if (!mpfr_number_p (mpc_imagref (log))
60 || mpfr_zero_p (mpc_imagref (log))) {
61 /* no need to divide by log(10) */
62 special_im = 1;
63 ok = 1;
64 }
65 else {
66 special_im = 0;
67 mpfr_const_log10 (log10);
68 mpfr_div (mpc_imagref (log), mpc_imagref (log), log10, MPFR_RNDN);
69
70 ok = mpfr_can_round (mpc_imagref (log), prec - 2,
71 MPFR_RNDN, MPFR_RNDZ,
72 MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == MPFR_RNDN));
73 }
74
75 if (ok) {
76 if (!mpfr_number_p (mpc_realref (log))
77 || mpfr_zero_p (mpc_realref (log)))
78 special_re = 1;
79 else {
80 special_re = 0;
81 if (special_im)
82 /* log10 not yet computed */
83 mpfr_const_log10 (log10);
84 mpfr_div (mpc_realref (log), mpc_realref (log), log10, MPFR_RNDN);
85 /* error <= 24/7 ulp < 4 ulp for prec >= 4, see algorithms.tex */
86
87 ok = mpfr_can_round (mpc_realref (log), prec - 2,
88 MPFR_RNDN, MPFR_RNDZ,
89 MPC_PREC_RE(rop) + (MPC_RND_RE (rnd) == MPFR_RNDN));
90 }
91
92 /* Special code to deal with cases where the real part of log10(x+i*y)
93 is exact, like x=3 and y=1. Since Re(log10(x+i*y)) = log10(x^2+y^2)/2
94 this happens whenever x^2+y^2 is a nonnegative power of 10.
95 Indeed x^2+y^2 cannot equal 10^(a/2^b) for a, b integers, a odd, b>0,
96 since x^2+y^2 is rational, and 10^(a/2^b) is irrational.
97 Similarly, for b=0, x^2+y^2 cannot equal 10^a for a < 0 since x^2+y^2
98 is a rational with denominator a power of 2.
99 Now let x^2+y^2 = 10^s. Without loss of generality we can assume
100 x = u/2^e and y = v/2^e with u, v, e integers: u^2+v^2 = 10^s*2^(2e)
101 thus u^2+v^2 = 0 mod 2^(2e). By recurrence on e, necessarily
102 u = v = 0 mod 2^e, thus x and y are necessarily integers.
103 */
104 if (!ok && !check_exact && mpfr_integer_p (mpc_realref (op)) &&
105 mpfr_integer_p (mpc_imagref (op))) {
106 mpz_t x, y;
107 unsigned long s, v;
108
109 check_exact = 1;
110 mpz_init (x);
111 mpz_init (y);
112 mpfr_get_z (x, mpc_realref (op), MPFR_RNDN); /* exact */
113 mpfr_get_z (y, mpc_imagref (op), MPFR_RNDN); /* exact */
114 mpz_mul (x, x, x);
115 mpz_mul (y, y, y);
116 mpz_add (x, x, y); /* x^2+y^2 */
117 v = mpz_scan1 (x, 0);
118 /* if x = 10^s then necessarily s = v */
119 s = mpz_sizeinbase (x, 10);
120 /* since s is either the number of digits of x or one more,
121 then x = 10^(s-1) or 10^(s-2) */
122 if (s == v + 1 || s == v + 2) {
123 mpz_div_2exp (x, x, v);
124 mpz_ui_pow_ui (y, 5, v);
125 if (mpz_cmp (y, x) == 0) {
126 /* Re(log10(x+i*y)) is exactly v/2
127 we reset the precision of Re(log) so that v can be
128 represented exactly */
129 mpfr_set_prec (mpc_realref (log),
130 sizeof(unsigned long)*CHAR_BIT);
131 mpfr_set_ui_2exp (mpc_realref (log), v, -1, MPFR_RNDN);
132 /* exact */
133 ok = 1;
134 }
135 }
136 mpz_clear (x);
137 mpz_clear (y);
138 }
139 }
140 }
141
142 inex_re = mpfr_set (mpc_realref(rop), mpc_realref (log), MPC_RND_RE (rnd));
143 if (special_re)
144 inex_re = MPC_INEX_RE (inex);
145 /* recover flag from call to mpc_log above */
146 inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (log), MPC_RND_IM (rnd));
147 if (special_im)
148 inex_im = MPC_INEX_IM (inex);
149 mpfr_clear (log10);
150 mpc_clear (log);
151
152 /* restore the exponent range, and check the range of results */
153 mpfr_set_emin (saved_emin);
154 mpfr_set_emax (saved_emax);
155 inex_re = mpfr_check_range (mpc_realref (rop), inex_re, MPC_RND_RE (rnd));
156 inex_im = mpfr_check_range (mpc_imagref (rop), inex_im, MPC_RND_IM (rnd));
157
158 return MPC_INEX(inex_re, inex_im);
159 }
160