src/sqr.c

   1 /* mpc_sqr -- Square a complex number.
   2
   3 Copyright (C) 2002, 2005, 2008, 2009 Andreas Enge, Paul Zimmermann, Philippe Th\'eveny
   4
   5 This file is part of the MPC Library.
   6
   7 The MPC Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 2.1 of the License, or (at your
  10 option) any later version.
  11
  12 The MPC Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the MPC Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include <stdio.h>    /* for MPC_ASSERT */
  23 #include "mpc-impl.h"
  24
  25 int
  26 mpc_sqr (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  27 {
  28    int ok;
  29    mpfr_t u, v;
  30    mpfr_t x;
  31       /* rop temporary variable to hold the real part of op,
  32          needed in the case rop==op */
  33    mp_prec_t prec;
  34    int inex_re, inex_im, inexact;
  35    mp_exp_t old_emax, old_emin, emin, emax;
  36
  37    /* special values: NaN and infinities */
  38    if (!mpfr_number_p (MPC_RE (op)) || !mpfr_number_p (MPC_IM (op))) {
  39       if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) {
  40          mpfr_set_nan (MPC_RE (rop));
  41          mpfr_set_nan (MPC_IM (rop));
  42       }
  43       else if (mpfr_inf_p (MPC_RE (op))) {
  44          if (mpfr_inf_p (MPC_IM (op))) {
  45             mpfr_set_nan (MPC_RE (rop));
  46             mpfr_set_inf (MPC_IM (rop),
  47                           MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op)));
  48          }
  49          else {
  50             mpfr_set_inf (MPC_RE (rop), +1);
  51             if (mpfr_zero_p (MPC_IM (op)))
  52                mpfr_set_nan (MPC_IM (rop));
  53             else
  54                mpfr_set_inf (MPC_IM (rop),
  55                              MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op)));
  56          }
  57       }
  58       else /* IM(op) is infinity, RE(op) is not */ {
  59          mpfr_set_inf (MPC_RE (rop), -1);
  60          if (mpfr_zero_p (MPC_RE (op)))
  61             mpfr_set_nan (MPC_IM (rop));
  62          else
  63             mpfr_set_inf (MPC_IM (rop),
  64                           MPFR_SIGN (MPC_RE (op)) * MPFR_SIGN (MPC_IM (op)));
  65       }
  66       return MPC_INEX (0, 0); /* exact */
  67    }
  68
  69    prec = MPC_MAX_PREC(rop);
  70
  71    /* first check for real resp. purely imaginary number */
  72    if (mpfr_zero_p (MPC_IM(op)))
  73    {
  74       int same_sign = mpfr_signbit (MPC_RE (op)) == mpfr_signbit (MPC_IM (op));
  75       inex_re = mpfr_sqr (MPC_RE(rop), MPC_RE(op), MPC_RND_RE(rnd));
  76       inex_im = mpfr_set_ui (MPC_IM(rop), 0ul, GMP_RNDN);
  77       if (!same_sign)
  78         mpc_conj (rop, rop, MPC_RNDNN);
  79       return MPC_INEX(inex_re, inex_im);
  80    }
  81    if (mpfr_zero_p (MPC_RE(op)))
  82    {
  83       int same_sign = mpfr_signbit (MPC_RE (op)) == mpfr_signbit (MPC_IM (op));
  84       inex_re = -mpfr_sqr (MPC_RE(rop), MPC_IM(op), INV_RND (MPC_RND_RE(rnd)));
  85       mpfr_neg (MPC_RE(rop), MPC_RE(rop), GMP_RNDN);
  86       inex_im = mpfr_set_ui (MPC_IM(rop), 0ul, GMP_RNDN);
  87       if (!same_sign)
  88         mpc_conj (rop, rop, MPC_RNDNN);
  89       return MPC_INEX(inex_re, inex_im);
  90    }
  91    /* If the real and imaginary part of the argument have rop very different */
  92    /* exponent, it is not reasonable to use Karatsuba squaring; compute    */
  93    /* exactly with the standard formulae instead, even if this means an    */
  94    /* additional multiplication.                                           */
  95    if (SAFE_ABS (mp_exp_t, MPFR_EXP (MPC_RE (op)) - MPFR_EXP (MPC_IM (op)))
  96        > (mp_exp_t) MPC_MAX_PREC (op) / 2)
  97    {
  98       mpfr_init2 (u, 2*MPFR_PREC (MPC_RE (op)));
  99       mpfr_init2 (v, 2*MPFR_PREC (MPC_IM (op)));
 100
 101       mpfr_sqr (u, MPC_RE (op), GMP_RNDN);
 102       mpfr_sqr (v, MPC_IM (op), GMP_RNDN); /* both are exact */
 103       inex_im = mpfr_mul (rop->im, op->re, op->im, MPC_RND_IM (rnd));
 104       mpfr_mul_2exp (rop->im, rop->im, 1, GMP_RNDN);
 105       inex_re = mpfr_sub (rop->re, u, v, MPC_RND_RE (rnd));
 106
 107       mpfr_clear (u);
 108       mpfr_clear (v);
 109       return MPC_INEX (inex_re, inex_im);
 110    }
 111
 112
 113    mpfr_init (u);
 114    mpfr_init (v);
 115
 116    if (rop == op)
 117    {
 118       mpfr_init2 (x, MPFR_PREC (op->re));
 119       mpfr_set (x, op->re, GMP_RNDN);
 120    }
 121    else
 122       x [0] = op->re [0];
 123
 124    /* store the maximal exponent */
 125    old_emax = mpfr_get_emax ();
 126    old_emin = mpfr_get_emin ();
 127    mpfr_set_emax (emax = mpfr_get_emax_max ());
 128    mpfr_set_emin (emin = mpfr_get_emin_min ());
 129
 130    do
 131    {
 132       prec += mpc_ceil_log2 (prec) + 5;
 133
 134       mpfr_set_prec (u, prec);
 135       mpfr_set_prec (v, prec);
 136
 137       /* Let op = x + iy. We need u = x+y and v = x-y, rounded away.       */
 138       /* As this is not implemented in mpfr, we round towards zero and    */
 139       /* add one ulp if the result is not exact. The error is bounded     */
 140       /* above by 1 ulp.                                                  */
 141       /* We first let inexact be 1 if the real part is not computed       */
 142       /* exactly and determine the sign later.                            */
 143       inexact = 0;
 144       if (mpfr_add (u, x, MPC_IM (op), GMP_RNDZ))
 145          /* we have to use x here instead of MPC_RE (op), as MPC_RE (rop)
 146             may be modified later in the attempt to assign u to it */
 147       {
 148          mpfr_add_one_ulp (u, GMP_RNDN);
 149          inexact = 1;
 150       }
 151       if (mpfr_sub (v, x, MPC_IM (op), GMP_RNDZ))
 152       {
 153          mpfr_add_one_ulp (v, GMP_RNDN);
 154          inexact = 1;
 155       }
 156
 157       /* compute the real part as u*v, rounded away                    */
 158       /* determine also the sign of inex_re                            */
 159       if (mpfr_sgn (u) == 0 || mpfr_sgn (v) == 0)
 160         {
 161           /* as we have rounded away, the result is exact */
 162           mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN);
 163           inex_re = 0;
 164           ok = 1;
 165         }
 166       else if (mpfr_sgn (u) * mpfr_sgn (v) > 0)
 167         {
 168           inexact |= mpfr_mul (u, u, v, GMP_RNDU); /* error 5 */
 169           /* checks that no overflow nor underflow occurs: since u*v > 0
 170              and we round up, an overflow will give +Inf, but an underflow
 171              will give 0.5*2^emin */
 172           MPC_ASSERT (mpfr_get_exp (u) != emin);
 173           if (mpfr_inf_p (u))
 174             {
 175               mp_rnd_t rnd_re = MPC_RND_RE (rnd);
 176               if (rnd_re == GMP_RNDZ || rnd_re == GMP_RNDD)
 177                 {
 178                   mpfr_set_ui_2exp (MPC_RE (rop), 1, emax, rnd_re);
 179                   inex_re = -1;
 180                 }
 181               else /* round up or away from zero */
 182                 inex_re = 1;
 183               break;
 184             }
 185           ok = (!inexact) | mpfr_can_round (u, prec - 3, GMP_RNDU, GMP_RNDZ,
 186                MPFR_PREC (MPC_RE (rop)) + (MPC_RND_RE (rnd) == GMP_RNDN));
 187           if (ok)
 188             {
 189               inex_re = mpfr_set (MPC_RE (rop), u, MPC_RND_RE (rnd));
 190               if (inex_re == 0)
 191                 /* remember that u was already rounded */
 192                 inex_re = inexact;
 193             }
 194         }
 195       else
 196         {
 197           inexact |= mpfr_mul (u, u, v, GMP_RNDD); /* error 5 */
 198           /* checks that no overflow occurs: since u*v < 0 and we round down,
 199              an overflow will give -Inf */
 200           MPC_ASSERT (mpfr_inf_p (u) == 0);
 201           /* if an underflow happens (i.e., u = -0.5*2^emin since we round
 202              away from zero), the result will be an underflow */
 203           if (mpfr_get_exp (u) == emin)
 204             {
 205               mp_rnd_t rnd_re = MPC_RND_RE (rnd);
 206               if (rnd_re == GMP_RNDZ || rnd_re == GMP_RNDN ||
 207                   rnd_re == GMP_RNDU)
 208                 {
 209                   mpfr_set_ui (MPC_RE (rop), 0, rnd_re);
 210                   inex_re = 1;
 211                 }
 212               else /* round down or away from zero */
 213                 inex_re = -1;
 214               break;
 215             }
 216           ok = (!inexact) | mpfr_can_round (u, prec - 3, GMP_RNDD, GMP_RNDZ,
 217                MPFR_PREC (MPC_RE (rop)) + (MPC_RND_RE (rnd) == GMP_RNDN));
 218           if (ok)
 219             {
 220               inex_re = mpfr_set (MPC_RE (rop), u, MPC_RND_RE (rnd));
 221               if (inex_re == 0)
 222                 inex_re = inexact;
 223             }
 224         }
 225    }
 226    while (!ok);
 227
 228    /* compute the imaginary part as 2*x*y, which is always possible */
 229    if (mpfr_get_exp (x) + mpfr_get_exp(MPC_IM (op)) <= emin + 1)
 230      {
 231        mpfr_mul_2ui (x, x, 1, GMP_RNDN);
 232        inex_im = mpfr_mul (MPC_IM (rop), x, MPC_IM (op), MPC_RND_IM (rnd));
 233      }
 234    else
 235      {
 236        inex_im = mpfr_mul (MPC_IM (rop), x, MPC_IM (op), MPC_RND_IM (rnd));
 237        /* checks that no underflow occurs (an overflow can occur here) */
 238        MPC_ASSERT (mpfr_zero_p (MPC_IM (rop)) == 0);
 239        mpfr_mul_2ui (MPC_IM (rop), MPC_IM (rop), 1, GMP_RNDN);
 240      }
 241
 242    mpfr_clear (u);
 243    mpfr_clear (v);
 244
 245    if (rop == op)
 246       mpfr_clear (x);
 247
 248    /* restore the exponent range */
 249    mpfr_set_emax (old_emax);
 250    mpfr_set_emin (old_emin);
 251
 252    return MPC_INEX (inex_re, inex_im);
 253 }