1/* mpc_norm -- Square of the norm of a complex number.
2
3Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011 INRIA
4
5This file is part of GNU MPC.
6
7GNU MPC is free software; you can redistribute it and/or modify it under
8the terms of the GNU Lesser General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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/* a <- norm(b) = b * conj(b)
25   (the rounding mode is mpfr_rnd_t here since we return an mpfr number) */
26int
27mpc_norm (mpfr_ptr a, mpc_srcptr b, mpfr_rnd_t rnd)
28{
29   int inexact;
30   int saved_underflow, saved_overflow;
31
32   /* handling of special values; consistent with abs in that
33      norm = abs^2; so norm (+-inf, xxx) = norm (xxx, +-inf) = +inf */
34   if (!mpc_fin_p (b))
35         return mpc_abs (a, b, rnd);
36   else if (mpfr_zero_p (mpc_realref (b))) {
37      if (mpfr_zero_p (mpc_imagref (b)))
38         return mpfr_set_ui (a, 0, rnd); /* +0 */
39      else
40         return mpfr_sqr (a, mpc_imagref (b), rnd);
41   }
42   else if (mpfr_zero_p (mpc_imagref (b)))
43     return mpfr_sqr (a, mpc_realref (b), rnd); /* Re(b) <> 0 */
44
45   else /* everything finite and non-zero */ {
46      mpfr_t u, v, res;
47      mpfr_prec_t prec, prec_u, prec_v;
48      int loops;
49      const int max_loops = 2;
50         /* switch to exact squarings when loops==max_loops */
51
52      prec = mpfr_get_prec (a);
53
54      mpfr_init (u);
55      mpfr_init (v);
56      mpfr_init (res);
57
58      /* save the underflow or overflow flags from MPFR */
59      saved_underflow = mpfr_underflow_p ();
60      saved_overflow = mpfr_overflow_p ();
61
62      loops = 0;
63      mpfr_clear_underflow ();
64      mpfr_clear_overflow ();
65      do {
66         loops++;
67         prec += mpc_ceil_log2 (prec) + 3;
68         if (loops >= max_loops) {
69            prec_u = 2 * MPC_PREC_RE (b);
70            prec_v = 2 * MPC_PREC_IM (b);
71         }
72         else {
73            prec_u = MPC_MIN (prec, 2 * MPC_PREC_RE (b));
74            prec_v = MPC_MIN (prec, 2 * MPC_PREC_IM (b));
75         }
76
77         mpfr_set_prec (u, prec_u);
78         mpfr_set_prec (v, prec_v);
79
80         inexact  = mpfr_sqr (u, mpc_realref(b), GMP_RNDD); /* err <= 1 ulp in prec */
81         inexact |= mpfr_sqr (v, mpc_imagref(b), GMP_RNDD); /* err <= 1 ulp in prec */
82
83         /* If loops = max_loops, inexact should be 0 here, except in case
84               of underflow or overflow.
85            If loops < max_loops and inexact is zero, we can exit the
86            while-loop since it only remains to add u and v into a. */
87         if (inexact) {
88             mpfr_set_prec (res, prec);
89             mpfr_add (res, u, v, GMP_RNDD); /* err <= 3 ulp in prec */
90         }
91
92      } while (loops < max_loops && inexact != 0
93               && !mpfr_can_round (res, prec - 2, GMP_RNDD, GMP_RNDU,
94                                   mpfr_get_prec (a) + (rnd == GMP_RNDN)));
95
96      if (!inexact)
97         /* squarings were exact, neither underflow nor overflow */
98         inexact = mpfr_add (a, u, v, rnd);
99      /* if there was an overflow in Re(b)^2 or Im(b)^2 or their sum,
100         since the norm is larger, there is an overflow for the norm */
101      else if (mpfr_overflow_p ()) {
102         /* replace by "correctly rounded overflow" */
103         mpfr_set_ui (a, 1ul, GMP_RNDN);
104         inexact = mpfr_mul_2ui (a, a, mpfr_get_emax (), rnd);
105      }
106      else if (mpfr_underflow_p ()) {
107         /* necessarily one of the squarings did underflow (otherwise their
108            sum could not underflow), thus one of u, v is zero. */
109         mpfr_exp_t emin = mpfr_get_emin ();
110
111         /* Now either both u and v are zero, or u is zero and v exact,
112            or v is zero and u exact.
113            In the latter case, Im(b)^2 < 2^(emin-1).
114            If ulp(u) >= 2^(emin+1) and norm(b) is not exactly
115            representable at the target precision, then rounding u+Im(b)^2
116            is equivalent to rounding u+2^(emin-1).
117            For instance, if exp(u)>0 and the target precision is smaller
118            than about |emin|, the norm is not representable. To make the
119            scaling in the "else" case work without underflow, we test
120            whether exp(u) is larger than a small negative number instead.
121            The second case is handled analogously.                        */
122         if (!mpfr_zero_p (u)
123             && mpfr_get_exp (u) - 2 * (mpfr_exp_t) prec_u > emin
124             && mpfr_get_exp (u) > -10) {
125               mpfr_set_prec (v, MPFR_PREC_MIN);
126               mpfr_set_ui_2exp (v, 1, emin - 1, GMP_RNDZ);
127               inexact = mpfr_add (a, u, v, rnd);
128         }
129         else if (!mpfr_zero_p (v)
130             && mpfr_get_exp (v) - 2 * (mpfr_exp_t) prec_v > emin
131             && mpfr_get_exp (v) > -10) {
132               mpfr_set_prec (u, MPFR_PREC_MIN);
133               mpfr_set_ui_2exp (u, 1, emin - 1, GMP_RNDZ);
134               inexact = mpfr_add (a, u, v, rnd);
135         }
136         else {
137            unsigned long int scale, exp_re, exp_im;
138            int inex_underflow;
139
140            /* scale the input to an average exponent close to 0 */
141            exp_re = (unsigned long int) (-mpfr_get_exp (mpc_realref (b)));
142            exp_im = (unsigned long int) (-mpfr_get_exp (mpc_imagref (b)));
143            scale = exp_re / 2 + exp_im / 2 + (exp_re % 2 + exp_im % 2) / 2;
144               /* (exp_re + exp_im) / 2, computed in a way avoiding
145                  integer overflow                                  */
146            if (mpfr_zero_p (u)) {
147               /* recompute the scaled value exactly */
148               mpfr_mul_2ui (u, mpc_realref (b), scale, GMP_RNDN);
149               mpfr_sqr (u, u, GMP_RNDN);
150            }
151            else /* just scale */
152               mpfr_mul_2ui (u, u, 2*scale, GMP_RNDN);
153            if (mpfr_zero_p (v)) {
154               mpfr_mul_2ui (v, mpc_imagref (b), scale, GMP_RNDN);
155               mpfr_sqr (v, v, GMP_RNDN);
156            }
157            else
158               mpfr_mul_2ui (v, v, 2*scale, GMP_RNDN);
159
160            inexact = mpfr_add (a, u, v, rnd);
161            mpfr_clear_underflow ();
162            inex_underflow = mpfr_div_2ui (a, a, 2*scale, rnd);
163            if (mpfr_underflow_p ())
164               inexact = inex_underflow;
165         }
166      }
167      else /* no problems, ternary value due to mpfr_can_round trick */
168         inexact = mpfr_set (a, res, rnd);
169
170      /* restore underflow and overflow flags from MPFR */
171      if (saved_underflow)
172        mpfr_set_underflow ();
173      if (saved_overflow)
174        mpfr_set_overflow ();
175
176      mpfr_clear (u);
177      mpfr_clear (v);
178      mpfr_clear (res);
179   }
180
181   return inexact;
182}
183