1 /* mpc_norm -- Square of the norm of a complex number.
3 Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011 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 */
24 /* a <- norm(b) = b * conj(b)
25 (the rounding mode is mpfr_rnd_t here since we return an mpfr number) */
27 mpc_norm (mpfr_ptr a
, mpc_srcptr b
, mpfr_rnd_t rnd
)
30 int saved_underflow
, saved_overflow
;
32 /* handling of special values; consistent with abs in that
33 norm = abs^2; so norm (+-inf, xxx) = norm (xxx, +-inf) = +inf */
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 */
40 return mpfr_sqr (a
, mpc_imagref (b
), rnd
);
42 else if (mpfr_zero_p (mpc_imagref (b
)))
43 return mpfr_sqr (a
, mpc_realref (b
), rnd
); /* Re(b) <> 0 */
45 else /* everything finite and non-zero */ {
47 mpfr_prec_t prec
, prec_u
, prec_v
;
49 const int max_loops
= 2;
50 /* switch to exact squarings when loops==max_loops */
52 prec
= mpfr_get_prec (a
);
58 /* save the underflow or overflow flags from MPFR */
59 saved_underflow
= mpfr_underflow_p ();
60 saved_overflow
= mpfr_overflow_p ();
63 mpfr_clear_underflow ();
64 mpfr_clear_overflow ();
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
);
73 prec_u
= MPC_MIN (prec
, 2 * MPC_PREC_RE (b
));
74 prec_v
= MPC_MIN (prec
, 2 * MPC_PREC_IM (b
));
77 mpfr_set_prec (u
, prec_u
);
78 mpfr_set_prec (v
, prec_v
);
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 */
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. */
88 mpfr_set_prec (res
, prec
);
89 mpfr_add (res
, u
, v
, GMP_RNDD
); /* err <= 3 ulp in prec */
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
)));
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
);
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 ();
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. */
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
);
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
);
137 unsigned long int scale
, exp_re
, exp_im
;
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
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
);
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
);
158 mpfr_mul_2ui (v
, v
, 2*scale
, GMP_RNDN
);
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
;
167 else /* no problems, ternary value due to mpfr_can_round trick */
168 inexact
= mpfr_set (a
, res
, rnd
);
170 /* restore underflow and overflow flags from MPFR */
172 mpfr_set_underflow ();
174 mpfr_set_overflow ();