priv/host_generic_maddf.c

/*---------------------------------------------------------------*/
/*--- begin                              host_generic_maddf.c ---*/
/*---------------------------------------------------------------*/

/* 
   Compute x * y + z as ternary operation.
   Copyright (C) 2010-2013 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, see
   <http://www.gnu.org/licenses/>.
*/

/* Generic helper functions for doing FMA, i.e. compute x * y + z
   as ternary operation.
   These are purely back-end entities and cannot be seen/referenced
   from IR. */

#include "libvex_basictypes.h"
#include "host_generic_maddf.h"
#include "main_util.h"

/* This implementation relies on Double being more than twice as
   precise as Float and uses rounding to odd in order to avoid problems
   with double rounding.
   See a paper by Boldo and Melquiond:
   http://www.lri.fr/~melquion/doc/08-tc.pdf  */

#define FORCE_EVAL(X) __asm __volatile__ ("" : : "m" (X))

#if defined(__x86_64__) && defined(__SSE2_MATH__)
# define ENV_TYPE unsigned int
/* Save current rounding mode into ENV, hold exceptions, set rounding
   mode to rounding toward zero.  */
# define ROUNDTOZERO(env) \
   do {							\
      unsigned int mxcsr;				\
      __asm __volatile__ ("stmxcsr %0" : "=m" (mxcsr));	\
      (env) = mxcsr;					\
      mxcsr = (mxcsr | 0x7f80) & ~0x3f;			\
      __asm __volatile__ ("ldmxcsr %0" : : "m" (mxcsr));\
   } while (0)
/* Restore exceptions from ENV, return if inexact exception has been raised
   since ROUNDTOZERO.  */
# define RESET_TESTINEXACT(env) \
   ({							\
      unsigned int mxcsr, ret;				\
      __asm __volatile__ ("stmxcsr %0" : "=m" (mxcsr));	\
      ret = (mxcsr >> 5) & 1;				\
      mxcsr = (mxcsr & 0x3d) | (env);			\
      __asm __volatile__ ("ldmxcsr %0" : : "m" (mxcsr));\
      ret;						\
   })
/* Return if inexact exception has been raised since ROUNDTOZERO.  */
# define TESTINEXACT() \
   ({							\
      unsigned int mxcsr;				\
      __asm __volatile__ ("stmxcsr %0" : "=m" (mxcsr));	\
      (mxcsr >> 5) & 1;					\
   })
#endif

#define DBL_MANT_DIG 53
#define IEEE754_DOUBLE_BIAS 0x3ff

union vg_ieee754_double {
   Double d;

   /* This is the IEEE 754 double-precision format.  */
   struct {
#ifdef VKI_BIG_ENDIAN
      unsigned int negative:1;
      unsigned int exponent:11;
      unsigned int mantissa0:20;
      unsigned int mantissa1:32;
#else
      unsigned int mantissa1:32;
      unsigned int mantissa0:20;
      unsigned int exponent:11;
      unsigned int negative:1;
#endif
   } ieee;
};

void VEX_REGPARM(3)
     h_generic_calc_MAddF32 ( /*OUT*/Float* res,
                               Float* argX, Float* argY, Float* argZ )
{
#ifndef ENV_TYPE
   /* Lame fallback implementation.  */
   *res = *argX * *argY + *argZ;
#else
   ENV_TYPE env;
   /* Multiplication is always exact.  */
   Double temp = (Double) *argX * (Double) *argY;
   union vg_ieee754_double u;

   ROUNDTOZERO (env);

   /* Perform addition with round to odd.  */
   u.d = temp + (Double) *argZ;
   /* Ensure the addition is not scheduled after fetestexcept call.  */
   FORCE_EVAL (u.d);

   /* Reset rounding mode and test for inexact simultaneously.  */
   int j = RESET_TESTINEXACT (env);

   if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
      u.ieee.mantissa1 |= j;

   /* And finally truncation with round to nearest.  */
   *res = (Float) u.d;
#endif
}


void VEX_REGPARM(3)
     h_generic_calc_MAddF64 ( /*OUT*/Double* res,
                               Double* argX, Double* argY, Double* argZ )
{
#ifndef ENV_TYPE
   /* Lame fallback implementation.  */
   *res = *argX * *argY + *argZ;
#else
   Double x = *argX, y = *argY, z = *argZ;
   union vg_ieee754_double u, v, w;
   int adjust = 0;
   u.d = x;
   v.d = y;
   w.d = z;
   if (UNLIKELY (u.ieee.exponent + v.ieee.exponent
                 >= 0x7ff + IEEE754_DOUBLE_BIAS - DBL_MANT_DIG)
       || UNLIKELY (u.ieee.exponent >= 0x7ff - DBL_MANT_DIG)
       || UNLIKELY (v.ieee.exponent >= 0x7ff - DBL_MANT_DIG)
       || UNLIKELY (w.ieee.exponent >= 0x7ff - DBL_MANT_DIG)
       || UNLIKELY (u.ieee.exponent + v.ieee.exponent
                    <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG)) {
      /* If z is Inf, but x and y are finite, the result should be
         z rather than NaN.  */
      if (w.ieee.exponent == 0x7ff
          && u.ieee.exponent != 0x7ff
          && v.ieee.exponent != 0x7ff) {
         *res = (z + x) + y;
         return;
      }
      /* If x or y or z is Inf/NaN, or if fma will certainly overflow,
         or if x * y is less than half of DBL_DENORM_MIN,
         compute as x * y + z.  */
      if (u.ieee.exponent == 0x7ff
          || v.ieee.exponent == 0x7ff
          || w.ieee.exponent == 0x7ff
          || u.ieee.exponent + v.ieee.exponent > 0x7ff + IEEE754_DOUBLE_BIAS
          || u.ieee.exponent + v.ieee.exponent
             < IEEE754_DOUBLE_BIAS - DBL_MANT_DIG - 2) {
         *res = x * y + z;
         return;
      }
      if (u.ieee.exponent + v.ieee.exponent
          >= 0x7ff + IEEE754_DOUBLE_BIAS - DBL_MANT_DIG) {
         /* Compute 1p-53 times smaller result and multiply
            at the end.  */
         if (u.ieee.exponent > v.ieee.exponent)
            u.ieee.exponent -= DBL_MANT_DIG;
         else
            v.ieee.exponent -= DBL_MANT_DIG;
         /* If x + y exponent is very large and z exponent is very small,
            it doesn't matter if we don't adjust it.  */
         if (w.ieee.exponent > DBL_MANT_DIG)
            w.ieee.exponent -= DBL_MANT_DIG;
         adjust = 1;
      } else if (w.ieee.exponent >= 0x7ff - DBL_MANT_DIG) {
         /* Similarly.
            If z exponent is very large and x and y exponents are
            very small, it doesn't matter if we don't adjust it.  */
         if (u.ieee.exponent > v.ieee.exponent) {
            if (u.ieee.exponent > DBL_MANT_DIG)
               u.ieee.exponent -= DBL_MANT_DIG;
         } else if (v.ieee.exponent > DBL_MANT_DIG)
            v.ieee.exponent -= DBL_MANT_DIG;
         w.ieee.exponent -= DBL_MANT_DIG;
         adjust = 1;
      } else if (u.ieee.exponent >= 0x7ff - DBL_MANT_DIG) {
         u.ieee.exponent -= DBL_MANT_DIG;
         if (v.ieee.exponent)
            v.ieee.exponent += DBL_MANT_DIG;
         else
            v.d *= 0x1p53;
      } else if (v.ieee.exponent >= 0x7ff - DBL_MANT_DIG) {
         v.ieee.exponent -= DBL_MANT_DIG;
         if (u.ieee.exponent)
            u.ieee.exponent += DBL_MANT_DIG;
         else
            u.d *= 0x1p53;
      } else /* if (u.ieee.exponent + v.ieee.exponent
                    <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG) */ {
         if (u.ieee.exponent > v.ieee.exponent)
            u.ieee.exponent += 2 * DBL_MANT_DIG;
         else
            v.ieee.exponent += 2 * DBL_MANT_DIG;
         if (w.ieee.exponent <= 4 * DBL_MANT_DIG + 4) {
            if (w.ieee.exponent)
               w.ieee.exponent += 2 * DBL_MANT_DIG;
            else
               w.d *= 0x1p106;
            adjust = -1;
         }
         /* Otherwise x * y should just affect inexact
            and nothing else.  */
      }
      x = u.d;
      y = v.d;
      z = w.d;
   }
   /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
#  define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
   Double x1 = x * C;
   Double y1 = y * C;
   Double m1 = x * y;
   x1 = (x - x1) + x1;
   y1 = (y - y1) + y1;
   Double x2 = x - x1;
   Double y2 = y - y1;
   Double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
#  undef C

   /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
   Double a1 = z + m1;
   Double t1 = a1 - z;
   Double t2 = a1 - t1;
   t1 = m1 - t1;
   t2 = z - t2;
   Double a2 = t1 + t2;

   ENV_TYPE env;
   ROUNDTOZERO (env);

   /* Perform m2 + a2 addition with round to odd.  */
   u.d = a2 + m2;

   if (UNLIKELY (adjust < 0)) {
      if ((u.ieee.mantissa1 & 1) == 0)
         u.ieee.mantissa1 |= TESTINEXACT ();
      v.d = a1 + u.d;
      /* Ensure the addition is not scheduled after fetestexcept call.  */
      FORCE_EVAL (v.d);
   }

   /* Reset rounding mode and test for inexact simultaneously.  */
   int j = RESET_TESTINEXACT (env) != 0;

   if (LIKELY (adjust == 0)) {
      if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
         u.ieee.mantissa1 |= j;
      /* Result is a1 + u.d.  */
      *res = a1 + u.d;
   } else if (LIKELY (adjust > 0)) {
      if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7ff)
         u.ieee.mantissa1 |= j;
      /* Result is a1 + u.d, scaled up.  */
      *res = (a1 + u.d) * 0x1p53;
   } else {
      /* If a1 + u.d is exact, the only rounding happens during
         scaling down.  */
      if (j == 0) {
         *res = v.d * 0x1p-106;
         return;
      }
      /* If result rounded to zero is not subnormal, no double
         rounding will occur.  */
      if (v.ieee.exponent > 106) {
         *res = (a1 + u.d) * 0x1p-106;
         return;
      }
      /* If v.d * 0x1p-106 with round to zero is a subnormal above
         or equal to DBL_MIN / 2, then v.d * 0x1p-106 shifts mantissa
         down just by 1 bit, which means v.ieee.mantissa1 |= j would
         change the round bit, not sticky or guard bit.
         v.d * 0x1p-106 never normalizes by shifting up,
         so round bit plus sticky bit should be already enough
         for proper rounding.  */
      if (v.ieee.exponent == 106) {
         /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
            v.ieee.mantissa1 & 1 is the round bit and j is our sticky
            bit.  In round-to-nearest 001 rounds down like 00,
            011 rounds up, even though 01 rounds down (thus we need
            to adjust), 101 rounds down like 10 and 111 rounds up
            like 11.  */
         if ((v.ieee.mantissa1 & 3) == 1) {
            v.d *= 0x1p-106;
            if (v.ieee.negative)
               *res = v.d - 0x1p-1074;
            else
               *res = v.d + 0x1p-1074;
         } else
            *res = v.d * 0x1p-106;
         return;
      }
      v.ieee.mantissa1 |= j;
      *res = v.d * 0x1p-106;
      return;
    }
#endif
}

/*---------------------------------------------------------------*/
/*--- end                                 host_generic_maddf.c --*/
/*---------------------------------------------------------------*/