1 /* mpc_norm -- Square of the norm of 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 /* a <- norm(b) = b * conj(b) 25 (the rounding mode is mpfr_rnd_t here since we return an mpfr number) */ 26 int 27 mpc_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), MPFR_RNDD); /* err <= 1 ulp in prec */ 81 inexact |= mpfr_sqr (v, mpc_imagref(b), MPFR_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, MPFR_RNDD); /* err <= 3 ulp in prec */ 90 } 91 92 } while (loops < max_loops && inexact != 0 93 && !mpfr_can_round (res, prec - 2, MPFR_RNDD, MPFR_RNDU, 94 mpfr_get_prec (a) + (rnd == MPFR_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, MPFR_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, MPFR_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, MPFR_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, MPFR_RNDN); 149 mpfr_sqr (u, u, MPFR_RNDN); 150 } 151 else /* just scale */ 152 mpfr_mul_2ui (u, u, 2*scale, MPFR_RNDN); 153 if (mpfr_zero_p (v)) { 154 mpfr_mul_2ui (v, mpc_imagref (b), scale, MPFR_RNDN); 155 mpfr_sqr (v, v, MPFR_RNDN); 156 } 157 else 158 mpfr_mul_2ui (v, v, 2*scale, MPFR_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
Indexes created Sat Sep 12 00:25:55 UTC 2026