remainderf.c File Reference

#include "libm.h"
#include "libm_util.h"
#include "libm_inlines.h"
#include "libm_errno.h"

Include dependency graph for remainderf.c:

Go to the source code of this file.

Macros
#define	USE_NANF_WITH_FLAGS
#define	USE_SCALEDOUBLE_1
#define	USE_GET_FPSW_INLINE
#define	USE_SET_FPSW_INLINE
#define	USE_HANDLE_ERRORF
#define	_CRTBLD_C9X
#define	_FUNCNAME   "remainderf"
#define	_OPERATION   OP_REM

Functions
float	remainderf (float x, float y)

Macro Definition Documentation

_CRTBLD_C9X

#define _CRTBLD_C9X

Definition at line 43 of file remainderf.c.

_FUNCNAME

#define _FUNCNAME   "remainderf"

_OPERATION

#define _OPERATION   OP_REM

USE_GET_FPSW_INLINE

#define USE_GET_FPSW_INLINE

Definition at line 32 of file remainderf.c.

USE_HANDLE_ERRORF

#define USE_HANDLE_ERRORF

Definition at line 34 of file remainderf.c.

USE_NANF_WITH_FLAGS

#define USE_NANF_WITH_FLAGS

Definition at line 30 of file remainderf.c.

USE_SCALEDOUBLE_1

#define USE_SCALEDOUBLE_1

Definition at line 31 of file remainderf.c.

USE_SET_FPSW_INLINE

#define USE_SET_FPSW_INLINE

Definition at line 33 of file remainderf.c.

Function Documentation

remainderf()

float remainderf	(	float	x,
		float	y
	)

Definition at line 60 of file remainderf.c.

{
  double dx, dy, scale, w, t;
  int i, ntimes, xexp, yexp;
  unsigned long long ux, uy, ax, ay;
 
  unsigned int sw;
 
  dx = x;
  dy = y;
 
 
  GET_BITS_DP64(dx, ux);
  GET_BITS_DP64(dy, uy);
  ax = ux & ~SIGNBIT_DP64;
  ay = uy & ~SIGNBIT_DP64;
  xexp = (int)((ux & EXPBITS_DP64) >> EXPSHIFTBITS_DP64);
  yexp = (int)((uy & EXPBITS_DP64) >> EXPSHIFTBITS_DP64);
 
  if (xexp < 1 || xexp > BIASEDEMAX_DP64 ||
      yexp < 1 || yexp > BIASEDEMAX_DP64)
    {
      /* x or y is zero, NaN or infinity (neither x nor y can be
         denormalized because we promoted from float to double) */
      if (xexp > BIASEDEMAX_DP64)
        {
          /* x is NaN or infinity */
          if (ux & MANTBITS_DP64)
            {
              /* x is NaN */
              unsigned int ufx;
              GET_BITS_SP32(x, ufx);
              return _handle_errorf(_FUNCNAME, _OPERATION, ufx|0x00400000, _DOMAIN, 0,
                                   EDOM, x, y, 2);
            }
          else
            {
              /* x is infinity; result is NaN */
              return _handle_errorf(_FUNCNAME, _OPERATION, INDEFBITPATT_SP32, _DOMAIN,
                                   AMD_F_INVALID, EDOM, x, y, 2);
            }
        }
      else if (yexp > BIASEDEMAX_DP64)
        {
          /* y is NaN or infinity */
          if (uy & MANTBITS_DP64)
            {
              /* y is NaN */
              unsigned int ufy;
              GET_BITS_SP32(y, ufy);
              return _handle_errorf(_FUNCNAME, _OPERATION, ufy|0x00400000, _DOMAIN, 0,
                                   EDOM, x, y, 2);
            }
          else
            {
#ifdef _CRTBLD_C9X
              /* C99 return for y = +-inf is x */
              return x;
#else
              /* y is infinity; result is indefinite */
              return _handle_errorf(_FUNCNAME, _OPERATION, INDEFBITPATT_SP32, _DOMAIN,
                                   AMD_F_INVALID, EDOM, x, y, 2);
#endif
            }
        }
      else if (xexp < 1)
        {
          /* x must be zero (cannot be denormalized) */
          if (yexp < 1)
            {
              /* y must be zero (cannot be denormalized) */
              return _handle_errorf(_FUNCNAME, _OPERATION, INDEFBITPATT_SP32, _DOMAIN,
                                   AMD_F_INVALID, EDOM, x, y, 2);
            }
          else
              /* C99 return for x = 0 must preserve sign */
              return x;
        }
      else
        {
          /* y must be zero */
          return _handle_errorf(_FUNCNAME, _OPERATION, INDEFBITPATT_SP32, _DOMAIN,
                               AMD_F_INVALID, EDOM, x, y, 2);
        }
    }
  else if (ax == ay)
    {
      /* abs(x) == abs(y); return zero with the sign of x */
      PUT_BITS_DP64(ux & SIGNBIT_DP64, dx);
      return (float)dx;
    }
 
  /* Set dx = abs(x), dy = abs(y) */
  PUT_BITS_DP64(ax, dx);
  PUT_BITS_DP64(ay, dy);
 
  if (ax < ay)
    {
      /* abs(x) < abs(y) */
#if !defined(COMPILING_FMOD)
      if (dx > 0.5*dy)
        dx -= dy;
#endif
      return (float)(x < 0.0? -dx : dx);
    }
 
  /* Save the current floating-point status word. We need
     to do this because the remainder function is always
     exact for finite arguments, but our algorithm causes
     the inexact flag to be raised. We therefore need to
     restore the entry status before exiting. */
  sw = get_fpsw_inline();
 
  /* Set ntimes to the number of times we need to do a
     partial remainder. If the exponent of x is an exact multiple
     of 24 larger than the exponent of y, and the mantissa of x is
     less than the mantissa of y, ntimes will be one too large
     but it doesn't matter - it just means that we'll go round
     the loop below one extra time. */
  if (xexp <= yexp)
    {
      ntimes = 0;
      w = dy;
      scale = 1.0;
    }
  else
    {
      ntimes = (xexp - yexp) / 24;
 
      /* Set w = y * 2^(24*ntimes) */
      PUT_BITS_DP64((unsigned long long)(ntimes * 24 + EXPBIAS_DP64) << EXPSHIFTBITS_DP64,
                    scale);
      w = scale * dy;
      /* Set scale = 2^(-24) */
      PUT_BITS_DP64((unsigned long long)(-24 + EXPBIAS_DP64) << EXPSHIFTBITS_DP64,
                    scale);
    }
 
 
  /* Each time round the loop we compute a partial remainder.
     This is done by subtracting a large multiple of w
     from x each time, where w is a scaled up version of y.
     The subtraction can be performed exactly when performed
     in double precision, and the result at each stage can
     fit exactly in a single precision number. */
  for (i = 0; i < ntimes; i++)
    {
      /* t is the integer multiple of w that we will subtract.
         We use a truncated value for t. */
      t = (double)((int)(dx / w));
      dx -= w * t;
      /* Scale w down by 2^(-24) for the next iteration */
      w *= scale;
    }
 
  /* One more time */
#if defined(COMPILING_FMOD)
  t = (double)((int)(dx / w));
  dx -= w * t;
#else
 {
  unsigned int todd;
  /* Variable todd says whether the integer t is odd or not */
  t = (double)((int)(dx / w));
  todd = ((int)(dx / w)) & 1;
  dx -= w * t;
 
  /* At this point, dx lies in the range [0,dy) */
  /* For the remainder function, we need to adjust dx
     so that it lies in the range (-y/2, y/2] by carefully
     subtracting w (== dy == y) if necessary. */
  if (dx > 0.5 * w || ((dx == 0.5 * w) && todd))
    dx -= w;
 }
#endif
 
  /* **** N.B. for some reason this breaks the 32 bit version
     of remainder when compiling with optimization. */
  /* Restore the entry status flags */
  set_fpsw_inline(sw);
 
  /* Set the result sign according to input argument x */
  return (float)(x < 0.0? -dx : dx);
 
}