1 /* mpc_sqr -- Square a complex number.
2
3 Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011, 2012 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://www.gnu.org/licenses/ .
19 */
20
21 #include <stdio.h> /* for MPC_ASSERT */
22 #include "mpc-impl.h"
23
24
25 static int
mpfr_fsss(mpfr_ptr z,mpfr_srcptr a,mpfr_srcptr c,mpfr_rnd_t rnd)26 mpfr_fsss (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr c, mpfr_rnd_t rnd)
27 {
28 /* Computes z = a^2 - c^2.
29 Assumes that a and c are finite and non-zero; so a squaring yielding
30 an infinity is an overflow, and a squaring yielding 0 is an underflow.
31 Assumes further that z is distinct from a and c. */
32
33 int inex;
34 mpfr_t u, v;
35
36 /* u=a^2, v=c^2 exactly */
37 mpfr_init2 (u, 2*mpfr_get_prec (a));
38 mpfr_init2 (v, 2*mpfr_get_prec (c));
39 mpfr_sqr (u, a, MPFR_RNDN);
40 mpfr_sqr (v, c, MPFR_RNDN);
41
42 /* tentatively compute z as u-v; here we need z to be distinct
43 from a and c to not lose the latter */
44 inex = mpfr_sub (z, u, v, rnd);
45
46 if (mpfr_inf_p (z)) {
47 /* replace by "correctly rounded overflow" */
48 mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), MPFR_RNDN);
49 inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd);
50 }
51 else if (mpfr_zero_p (u) && !mpfr_zero_p (v)) {
52 /* exactly u underflowed, determine inexact flag */
53 inex = (mpfr_signbit (u) ? 1 : -1);
54 }
55 else if (mpfr_zero_p (v) && !mpfr_zero_p (u)) {
56 /* exactly v underflowed, determine inexact flag */
57 inex = (mpfr_signbit (v) ? -1 : 1);
58 }
59 else if (mpfr_nan_p (z) || (mpfr_zero_p (u) && mpfr_zero_p (v))) {
60 /* In the first case, u and v are +inf.
61 In the second case, u and v are zeroes; their difference may be 0
62 or the least representable number, with a sign to be determined.
63 Redo the computations with mpz_t exponents */
64 mpfr_exp_t ea, ec;
65 mpz_t eu, ev;
66 /* cheat to work around the const qualifiers */
67
68 /* Normalise the input by shifting and keep track of the shifts in
69 the exponents of u and v */
70 ea = mpfr_get_exp (a);
71 ec = mpfr_get_exp (c);
72
73 mpfr_set_exp ((mpfr_ptr) a, (mpfr_prec_t) 0);
74 mpfr_set_exp ((mpfr_ptr) c, (mpfr_prec_t) 0);
75
76 mpz_init (eu);
77 mpz_init (ev);
78 mpz_set_si (eu, (long int) ea);
79 mpz_mul_2exp (eu, eu, 1);
80 mpz_set_si (ev, (long int) ec);
81 mpz_mul_2exp (ev, ev, 1);
82
83 /* recompute u and v and move exponents to eu and ev */
84 mpfr_sqr (u, a, MPFR_RNDN);
85 /* exponent of u is non-positive */
86 mpz_sub_ui (eu, eu, (unsigned long int) (-mpfr_get_exp (u)));
87 mpfr_set_exp (u, (mpfr_prec_t) 0);
88 mpfr_sqr (v, c, MPFR_RNDN);
89 mpz_sub_ui (ev, ev, (unsigned long int) (-mpfr_get_exp (v)));
90 mpfr_set_exp (v, (mpfr_prec_t) 0);
91 if (mpfr_nan_p (z)) {
92 mpfr_exp_t emax = mpfr_get_emax ();
93 int overflow;
94 /* We have a = ma * 2^ea with 1/2 <= |ma| < 1 and ea <= emax.
95 So eu <= 2*emax, and eu > emax since we have
96 an overflow. The same holds for ev. Shift u and v by as much as
97 possible so that one of them has exponent emax and the
98 remaining exponents in eu and ev are the same. Then carry out
99 the addition. Shifting u and v prevents an underflow. */
100 if (mpz_cmp (eu, ev) >= 0) {
101 mpfr_set_exp (u, emax);
102 mpz_sub_ui (eu, eu, (long int) emax);
103 mpz_sub (ev, ev, eu);
104 mpfr_set_exp (v, (mpfr_exp_t) mpz_get_ui (ev));
105 /* remaining common exponent is now in eu */
106 }
107 else {
108 mpfr_set_exp (v, emax);
109 mpz_sub_ui (ev, ev, (long int) emax);
110 mpz_sub (eu, eu, ev);
111 mpfr_set_exp (u, (mpfr_exp_t) mpz_get_ui (eu));
112 mpz_set (eu, ev);
113 /* remaining common exponent is now also in eu */
114 }
115 inex = mpfr_sub (z, u, v, rnd);
116 /* Result is finite since u and v have the same sign. */
117 overflow = mpfr_mul_2ui (z, z, mpz_get_ui (eu), rnd);
118 if (overflow)
119 inex = overflow;
120 }
121 else {
122 int underflow;
123 /* Subtraction of two zeroes. We have a = ma * 2^ea
124 with 1/2 <= |ma| < 1 and ea >= emin and similarly for b.
125 So 2*emin < 2*emin+1 <= eu < emin < 0, and analogously for v. */
126 mpfr_exp_t emin = mpfr_get_emin ();
127 if (mpz_cmp (eu, ev) <= 0) {
128 mpfr_set_exp (u, emin);
129 mpz_add_ui (eu, eu, (unsigned long int) (-emin));
130 mpz_sub (ev, ev, eu);
131 mpfr_set_exp (v, (mpfr_exp_t) mpz_get_si (ev));
132 }
133 else {
134 mpfr_set_exp (v, emin);
135 mpz_add_ui (ev, ev, (unsigned long int) (-emin));
136 mpz_sub (eu, eu, ev);
137 mpfr_set_exp (u, (mpfr_exp_t) mpz_get_si (eu));
138 mpz_set (eu, ev);
139 }
140 inex = mpfr_sub (z, u, v, rnd);
141 mpz_neg (eu, eu);
142 underflow = mpfr_div_2ui (z, z, mpz_get_ui (eu), rnd);
143 if (underflow)
144 inex = underflow;
145 }
146
147 mpz_clear (eu);
148 mpz_clear (ev);
149
150 mpfr_set_exp ((mpfr_ptr) a, ea);
151 mpfr_set_exp ((mpfr_ptr) c, ec);
152 /* works also when a == c */
153 }
154
155 mpfr_clear (u);
156 mpfr_clear (v);
157
158 return inex;
159 }
160
161
162 int
mpc_sqr(mpc_ptr rop,mpc_srcptr op,mpc_rnd_t rnd)163 mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
164 {
165 int ok;
166 mpfr_t u, v;
167 mpfr_t x;
168 /* temporary variable to hold the real part of op,
169 needed in the case rop==op */
170 mpfr_prec_t prec;
171 int inex_re, inex_im, inexact;
172 mpfr_exp_t emin;
173 int saved_underflow;
174
175 /* special values: NaN and infinities */
176 if (!mpc_fin_p (op)) {
177 if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) {
178 mpfr_set_nan (mpc_realref (rop));
179 mpfr_set_nan (mpc_imagref (rop));
180 }
181 else if (mpfr_inf_p (mpc_realref (op))) {
182 if (mpfr_inf_p (mpc_imagref (op))) {
183 mpfr_set_inf (mpc_imagref (rop),
184 MPFR_SIGN (mpc_realref (op)) * MPFR_SIGN (mpc_imagref (op)));
185 mpfr_set_nan (mpc_realref (rop));
186 }
187 else {
188 if (mpfr_zero_p (mpc_imagref (op)))
189 mpfr_set_nan (mpc_imagref (rop));
190 else
191 mpfr_set_inf (mpc_imagref (rop),
192 MPFR_SIGN (mpc_realref (op)) * MPFR_SIGN (mpc_imagref (op)));
193 mpfr_set_inf (mpc_realref (rop), +1);
194 }
195 }
196 else /* IM(op) is infinity, RE(op) is not */ {
197 if (mpfr_zero_p (mpc_realref (op)))
198 mpfr_set_nan (mpc_imagref (rop));
199 else
200 mpfr_set_inf (mpc_imagref (rop),
201 MPFR_SIGN (mpc_realref (op)) * MPFR_SIGN (mpc_imagref (op)));
202 mpfr_set_inf (mpc_realref (rop), -1);
203 }
204 return MPC_INEX (0, 0); /* exact */
205 }
206
207 prec = MPC_MAX_PREC(rop);
208
209 /* Check for real resp. purely imaginary number */
210 if (mpfr_zero_p (mpc_imagref(op))) {
211 int same_sign = mpfr_signbit (mpc_realref (op)) == mpfr_signbit (mpc_imagref (op));
212 inex_re = mpfr_sqr (mpc_realref(rop), mpc_realref(op), MPC_RND_RE(rnd));
213 inex_im = mpfr_set_ui (mpc_imagref(rop), 0ul, MPFR_RNDN);
214 if (!same_sign)
215 mpc_conj (rop, rop, MPC_RNDNN);
216 return MPC_INEX(inex_re, inex_im);
217 }
218 if (mpfr_zero_p (mpc_realref(op))) {
219 int same_sign = mpfr_signbit (mpc_realref (op)) == mpfr_signbit (mpc_imagref (op));
220 inex_re = -mpfr_sqr (mpc_realref(rop), mpc_imagref(op), INV_RND (MPC_RND_RE(rnd)));
221 mpfr_neg (mpc_realref(rop), mpc_realref(rop), MPFR_RNDN);
222 inex_im = mpfr_set_ui (mpc_imagref(rop), 0ul, MPFR_RNDN);
223 if (!same_sign)
224 mpc_conj (rop, rop, MPC_RNDNN);
225 return MPC_INEX(inex_re, inex_im);
226 }
227
228 if (rop == op)
229 {
230 mpfr_init2 (x, MPC_PREC_RE (op));
231 mpfr_set (x, op->re, MPFR_RNDN);
232 }
233 else
234 x [0] = op->re [0];
235 /* From here on, use x instead of op->re and safely overwrite rop->re. */
236
237 /* Compute real part of result. */
238 if (SAFE_ABS (mpfr_exp_t,
239 mpfr_get_exp (mpc_realref (op)) - mpfr_get_exp (mpc_imagref (op)))
240 > (mpfr_exp_t) MPC_MAX_PREC (op) / 2) {
241 /* If the real and imaginary parts of the argument have very different
242 exponents, it is not reasonable to use Karatsuba squaring; compute
243 exactly with the standard formulae instead, even if this means an
244 additional multiplication. Using the approach copied from mul, over-
245 and underflows are also handled correctly. */
246
247 inex_re = mpfr_fsss (rop->re, x, op->im, MPC_RND_RE (rnd));
248 }
249 else {
250 /* Karatsuba squaring: we compute the real part as (x+y)*(x-y) and the
251 imaginary part as 2*x*y, with a total of 2M instead of 2S+1M for the
252 naive algorithm, which computes x^2-y^2 and 2*y*y */
253 mpfr_init (u);
254 mpfr_init (v);
255
256 emin = mpfr_get_emin ();
257
258 do
259 {
260 prec += mpc_ceil_log2 (prec) + 5;
261
262 mpfr_set_prec (u, prec);
263 mpfr_set_prec (v, prec);
264
265 /* Let op = x + iy. We need u = x+y and v = x-y, rounded away. */
266 /* The error is bounded above by 1 ulp. */
267 /* We first let inexact be 1 if the real part is not computed */
268 /* exactly and determine the sign later. */
269 inexact = mpfr_add (u, x, mpc_imagref (op), MPFR_RNDA)
270 | mpfr_sub (v, x, mpc_imagref (op), MPFR_RNDA);
271
272 /* compute the real part as u*v, rounded away */
273 /* determine also the sign of inex_re */
274
275 if (mpfr_sgn (u) == 0 || mpfr_sgn (v) == 0) {
276 /* as we have rounded away, the result is exact */
277 mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN);
278 inex_re = 0;
279 ok = 1;
280 }
281 else {
282 inexact |= mpfr_mul (u, u, v, MPFR_RNDA); /* error 5 */
283 if (mpfr_get_exp (u) == emin || mpfr_inf_p (u)) {
284 /* under- or overflow */
285 inex_re = mpfr_fsss (rop->re, x, op->im, MPC_RND_RE (rnd));
286 ok = 1;
287 }
288 else {
289 ok = (!inexact) | mpfr_can_round (u, prec - 3,
290 MPFR_RNDA, MPFR_RNDZ,
291 MPC_PREC_RE (rop) + (MPC_RND_RE (rnd) == MPFR_RNDN));
292 if (ok) {
293 inex_re = mpfr_set (mpc_realref (rop), u, MPC_RND_RE (rnd));
294 if (inex_re == 0)
295 /* remember that u was already rounded */
296 inex_re = inexact;
297 }
298 }
299 }
300 }
301 while (!ok);
302
303 mpfr_clear (u);
304 mpfr_clear (v);
305 }
306
307 saved_underflow = mpfr_underflow_p ();
308 mpfr_clear_underflow ();
309 inex_im = mpfr_mul (rop->im, x, op->im, MPC_RND_IM (rnd));
310 if (!mpfr_underflow_p ())
311 inex_im |= mpfr_mul_2ui (rop->im, rop->im, 1, MPC_RND_IM (rnd));
312 /* We must not multiply by 2 if rop->im has been set to the smallest
313 representable number. */
314 if (saved_underflow)
315 mpfr_set_underflow ();
316
317 if (rop == op)
318 mpfr_clear (x);
319
320 return MPC_INEX (inex_re, inex_im);
321 }
322