contrib/mpfr/fms.c

   1 /* mpfr_fms -- Floating multiply-subtract
   2
   3 Copyright 2001, 2002, 2004, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute
   9 it and/or modify it under the terms of the GNU Lesser
  10 General Public License as published by the Free Software
  11 Foundation; either version 2.1 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MPFR Library is distributed in the hope that it will
  15 be useful, but WITHOUT ANY WARRANTY; without even the
  16 implied warranty of MERCHANTABILITY or FITNESS FOR A
  17 PARTICULAR PURPOSE.  See the GNU Lesser General Public
  18 License for more details.
  19
  20 You should have received a copy of the GNU Lesser
  21 General Public License along with the GNU MPFR Library; see
  22 the file COPYING.LIB.  If not, write to the Free Software
  23 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  24 MA 02110-1301, USA. */
  25
  26 #include "mpfr-impl.h"
  27
  28 /* The fused-multiply-subtract (fms) of x, y and z is defined by:
  29    fms(x,y,z)= x*y - z
  30    Note: this is neither in IEEE754R, nor in LIA-2, but both the
  31    PowerPC and the Itanium define fms as x*y - z.
  32 */
  33
  34 int
  35 mpfr_fms (mpfr_ptr s, mpfr_srcptr x, mpfr_srcptr y, mpfr_srcptr z,
  36           mp_rnd_t rnd_mode)
  37 {
  38   int inexact;
  39   mpfr_t u;
  40   MPFR_SAVE_EXPO_DECL (expo);
  41
  42   /* particular cases */
  43   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ||
  44                      MPFR_IS_SINGULAR(y) ||
  45                      MPFR_IS_SINGULAR(z) ))
  46     {
  47       if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y) || MPFR_IS_NAN(z))
  48         {
  49           MPFR_SET_NAN(s);
  50           MPFR_RET_NAN;
  51         }
  52       /* now neither x, y or z is NaN */
  53       else if (MPFR_IS_INF(x) || MPFR_IS_INF(y))
  54         {
  55           /* cases Inf*0-z, 0*Inf-z, Inf-Inf */
  56           if ((MPFR_IS_ZERO(y)) ||
  57               (MPFR_IS_ZERO(x)) ||
  58               (MPFR_IS_INF(z) &&
  59                ((MPFR_MULT_SIGN(MPFR_SIGN(x), MPFR_SIGN(y))) == MPFR_SIGN(z))))
  60             {
  61               MPFR_SET_NAN(s);
  62               MPFR_RET_NAN;
  63             }
  64           else if (MPFR_IS_INF(z)) /* case Inf-Inf already checked above */
  65             {
  66               MPFR_SET_INF(s);
  67               MPFR_SET_OPPOSITE_SIGN(s, z);
  68               MPFR_RET(0);
  69             }
  70           else /* z is finite */
  71             {
  72               MPFR_SET_INF(s);
  73               MPFR_SET_SIGN(s, MPFR_MULT_SIGN(MPFR_SIGN(x) , MPFR_SIGN(y)));
  74               MPFR_RET(0);
  75             }
  76         }
  77       /* now x and y are finite */
  78       else if (MPFR_IS_INF(z))
  79         {
  80           MPFR_SET_INF(s);
  81           MPFR_SET_OPPOSITE_SIGN(s, z);
  82           MPFR_RET(0);
  83         }
  84       else if (MPFR_IS_ZERO(x) || MPFR_IS_ZERO(y))
  85         {
  86           if (MPFR_IS_ZERO(z))
  87             {
  88               int sign_p;
  89               sign_p = MPFR_MULT_SIGN( MPFR_SIGN(x) , MPFR_SIGN(y) );
  90               MPFR_SET_SIGN(s,(rnd_mode != GMP_RNDD ?
  91                                ((MPFR_IS_NEG_SIGN(sign_p) && MPFR_IS_POS(z))
  92                                 ? -1 : 1) :
  93                                ((MPFR_IS_POS_SIGN(sign_p) && MPFR_IS_NEG(z))
  94                                 ? 1 : -1)));
  95               MPFR_SET_ZERO(s);
  96               MPFR_RET(0);
  97             }
  98           else
  99             return mpfr_neg (s, z, rnd_mode);
 100         }
 101       else /* necessarily z is zero here */
 102         {
 103           MPFR_ASSERTD(MPFR_IS_ZERO(z));
 104           return mpfr_mul (s, x, y, rnd_mode);
 105         }
 106     }
 107
 108   /* If we take prec(u) >= prec(x) + prec(y), the product u <- x*y
 109      is exact, except in case of overflow or underflow. */
 110   MPFR_SAVE_EXPO_MARK (expo);
 111   mpfr_init2 (u, MPFR_PREC(x) + MPFR_PREC(y));
 112
 113   if (MPFR_UNLIKELY (mpfr_mul (u, x, y, GMP_RNDN)))
 114     {
 115       /* overflow or underflow - this case is regarded as rare, thus
 116          does not need to be very efficient (even if some tests below
 117          could have been done earlier).
 118          It is an overflow iff u is an infinity (since GMP_RNDN was used).
 119          Alternatively, we could test the overflow flag, but in this case,
 120          mpfr_clear_flags would have been necessary. */
 121       if (MPFR_IS_INF (u))  /* overflow */
 122         {
 123           /* Let's eliminate the obvious case where x*y and z have the
 124              same sign. No possible cancellation -> real overflow.
 125              Also, we know that |z| < 2^emax. If E(x) + E(y) >= emax+3,
 126              then |x*y| >= 2^(emax+1), and |x*y - z| >= 2^emax. This case
 127              is also an overflow. */
 128           if (MPFR_SIGN (u) != MPFR_SIGN (z) ||
 129               MPFR_GET_EXP (x) + MPFR_GET_EXP (y) >= __gmpfr_emax + 3)
 130             {
 131               mpfr_clear (u);
 132               MPFR_SAVE_EXPO_FREE (expo);
 133               return mpfr_overflow (s, rnd_mode, - MPFR_SIGN (z));
 134             }
 135
 136           /* E(x) + E(y) <= emax+2, therefore |x*y| < 2^(emax+2), and
 137              (x/4)*y does not overflow (let's recall that the result
 138              is exact with an unbounded exponent range). It does not
 139              underflow either, because x*y overflows and the exponent
 140              range is large enough. */
 141           inexact = mpfr_div_2ui (u, x, 2, GMP_RNDN);
 142           MPFR_ASSERTN (inexact == 0);
 143           inexact = mpfr_mul (u, u, y, GMP_RNDN);
 144           MPFR_ASSERTN (inexact == 0);
 145
 146           /* Now, we need to subtract z/4... But it may underflow! */
 147           {
 148             mpfr_t zo4;
 149             mpfr_srcptr zz;
 150             MPFR_BLOCK_DECL (flags);
 151
 152             if (MPFR_GET_EXP (u) > MPFR_GET_EXP (z) &&
 153                 MPFR_GET_EXP (u) - MPFR_GET_EXP (z) > MPFR_PREC (u))
 154               {
 155                 /* |z| < ulp(u)/2, therefore one can use z instead of z/4. */
 156                 zz = z;
 157               }
 158             else
 159               {
 160                 mpfr_init2 (zo4, MPFR_PREC (z));
 161                 if (mpfr_div_2ui (zo4, z, 2, GMP_RNDZ))
 162                   {
 163                     /* The division by 4 underflowed! */
 164                     MPFR_ASSERTN (0); /* TODO... */
 165                   }
 166                 zz = zo4;
 167               }
 168
 169             /* Let's recall that u = x*y/4 and zz = z/4 (or z if the
 170                following subtraction would give the same result). */
 171             MPFR_BLOCK (flags, inexact = mpfr_sub (s, u, zz, rnd_mode));
 172             /* u and zz have the same sign, so that an overflow
 173                is not possible. But an underflow is theoretically
 174                possible! */
 175             if (MPFR_UNDERFLOW (flags))
 176               {
 177                 MPFR_ASSERTN (zz != z);
 178                 MPFR_ASSERTN (0); /* TODO... */
 179                 mpfr_clears (zo4, u, (mpfr_ptr) 0);
 180               }
 181             else
 182               {
 183                 int inex2;
 184
 185                 if (zz != z)
 186                   mpfr_clear (zo4);
 187                 mpfr_clear (u);
 188                 MPFR_ASSERTN (! MPFR_OVERFLOW (flags));
 189                 inex2 = mpfr_mul_2ui (s, s, 2, rnd_mode);
 190                 if (inex2)  /* overflow */
 191                   {
 192                     inexact = inex2;
 193                     MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 194                   }
 195                 goto end;
 196               }
 197           }
 198         }
 199       else  /* underflow: one has |xy| < 2^(emin-1). */
 200         {
 201           unsigned long scale = 0;
 202           mpfr_t scaled_z;
 203           mpfr_srcptr new_z;
 204           mp_exp_t diffexp;
 205           mp_prec_t pzs;
 206           int xy_underflows;
 207
 208           /* Let's scale z so that ulp(z) > 2^emin and ulp(s) > 2^emin
 209              (the + 1 on MPFR_PREC (s) is necessary because the exponent
 210              of the result can be EXP(z) - 1). */
 211           diffexp = MPFR_GET_EXP (z) - __gmpfr_emin;
 212           pzs = MAX (MPFR_PREC (z), MPFR_PREC (s) + 1);
 213           if (diffexp <= pzs)
 214             {
 215               mp_exp_unsigned_t uscale;
 216               mpfr_t scaled_v;
 217               MPFR_BLOCK_DECL (flags);
 218
 219               uscale = (mp_exp_unsigned_t) pzs - diffexp + 1;
 220               MPFR_ASSERTN (uscale > 0);
 221               MPFR_ASSERTN (uscale <= ULONG_MAX);
 222               scale = uscale;
 223               mpfr_init2 (scaled_z, MPFR_PREC (z));
 224               inexact = mpfr_mul_2ui (scaled_z, z, scale, GMP_RNDN);
 225               MPFR_ASSERTN (inexact == 0);  /* TODO: overflow case */
 226               new_z = scaled_z;
 227               /* Now we need to recompute u = xy * 2^scale. */
 228               MPFR_BLOCK (flags,
 229                           if (MPFR_GET_EXP (x) < MPFR_GET_EXP (y))
 230                             {
 231                               mpfr_init2 (scaled_v, MPFR_PREC (x));
 232                               mpfr_mul_2ui (scaled_v, x, scale, GMP_RNDN);
 233                               mpfr_mul (u, scaled_v, y, GMP_RNDN);
 234                             }
 235                           else
 236                             {
 237                               mpfr_init2 (scaled_v, MPFR_PREC (y));
 238                               mpfr_mul_2ui (scaled_v, y, scale, GMP_RNDN);
 239                               mpfr_mul (u, x, scaled_v, GMP_RNDN);
 240                             });
 241               mpfr_clear (scaled_v);
 242               MPFR_ASSERTN (! MPFR_OVERFLOW (flags));
 243               xy_underflows = MPFR_UNDERFLOW (flags);
 244             }
 245           else
 246             {
 247               new_z = z;
 248               xy_underflows = 1;
 249             }
 250
 251           if (xy_underflows)
 252             {
 253               /* Let's replace xy by sign(xy) * 2^(emin-1). */
 254               mpfr_set_prec (u, MPFR_PREC_MIN);
 255               mpfr_setmin (u, __gmpfr_emin);
 256               MPFR_SET_SIGN (u, MPFR_MULT_SIGN (MPFR_SIGN (x),
 257                                                 MPFR_SIGN (y)));
 258             }
 259
 260           {
 261             MPFR_BLOCK_DECL (flags);
 262
 263             MPFR_BLOCK (flags, inexact = mpfr_sub (s, u, new_z, rnd_mode));
 264             mpfr_clear (u);
 265
 266             if (scale != 0)
 267               {
 268                 int inex2;
 269
 270                 mpfr_clear (scaled_z);
 271                 /* Here an overflow is theoretically possible, in which case
 272                    the result may be wrong, hence the assert. An underflow
 273                    is not possible, but let's check that anyway. */
 274                 MPFR_ASSERTN (! MPFR_OVERFLOW (flags));  /* TODO... */
 275                 MPFR_ASSERTN (! MPFR_UNDERFLOW (flags));  /* not possible */
 276                 inex2 = mpfr_div_2ui (s, s, scale, GMP_RNDN);
 277                 /* FIXME: this seems incorrect. GMP_RNDN -> rnd_mode?
 278                    Also, handle the double rounding case:
 279                    s / 2^scale = 2^(emin - 2) in GMP_RNDN. */
 280                 if (inex2)  /* underflow */
 281                   inexact = inex2;
 282               }
 283           }
 284
 285           /* FIXME/TODO: I'm not sure that the following is correct.
 286              Check for possible spurious exceptions due to intermediate
 287              computations. */
 288           MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 289           goto end;
 290         }
 291     }
 292
 293   inexact = mpfr_sub (s, u, z, rnd_mode);
 294   mpfr_clear (u);
 295   MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags);
 296  end:
 297   MPFR_SAVE_EXPO_FREE (expo);
 298   return mpfr_check_range (s, inexact, rnd_mode);
 299 }