contrib/mpc/src/norm.c

   1 /* mpc_norm -- Square of the norm of a complex number.
   2
   3 Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011 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), 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 }