__rem_pio2f.c raw

   1  /* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2f.c */
   2  /*
   3   * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   4   * Debugged and optimized by Bruce D. Evans.
   5   */
   6  /*
   7   * ====================================================
   8   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   9   *
  10   * Developed at SunPro, a Sun Microsystems, Inc. business.
  11   * Permission to use, copy, modify, and distribute this
  12   * software is freely granted, provided that this notice
  13   * is preserved.
  14   * ====================================================
  15   */
  16  /* __rem_pio2f(x,y)
  17   *
  18   * return the remainder of x rem pi/2 in *y
  19   * use double precision for everything except passing x
  20   * use __rem_pio2_large() for large x
  21   */
  22  
  23  #include "libm.h"
  24  
  25  #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1
  26  #define EPS DBL_EPSILON
  27  #elif FLT_EVAL_METHOD==2
  28  #define EPS LDBL_EPSILON
  29  #endif
  30  
  31  /*
  32   * invpio2:  53 bits of 2/pi
  33   * pio2_1:   first 25 bits of pi/2
  34   * pio2_1t:  pi/2 - pio2_1
  35   */
  36  static const double
  37  toint   = 1.5/EPS,
  38  pio4    = 0x1.921fb6p-1,
  39  invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
  40  pio2_1  = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */
  41  pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
  42  
  43  int __rem_pio2f(float x, double *y)
  44  {
  45  	union {float f; uint32_t i;} u = {x};
  46  	double tx[1],ty[1];
  47  	double_t fn;
  48  	uint32_t ix;
  49  	int n, sign, e0;
  50  
  51  	ix = u.i & 0x7fffffff;
  52  	/* 25+53 bit pi is good enough for medium size */
  53  	if (ix < 0x4dc90fdb) {  /* |x| ~< 2^28*(pi/2), medium size */
  54  		/* Use a specialized rint() to get fn. */
  55  		fn = (double_t)x*invpio2 + toint - toint;
  56  		n  = (int32_t)fn;
  57  		*y = x - fn*pio2_1 - fn*pio2_1t;
  58  		/* Matters with directed rounding. */
  59  		if (predict_false(*y < -pio4)) {
  60  			n--;
  61  			fn--;
  62  			*y = x - fn*pio2_1 - fn*pio2_1t;
  63  		} else if (predict_false(*y > pio4)) {
  64  			n++;
  65  			fn++;
  66  			*y = x - fn*pio2_1 - fn*pio2_1t;
  67  		}
  68  		return n;
  69  	}
  70  	if(ix>=0x7f800000) {  /* x is inf or NaN */
  71  		*y = x-x;
  72  		return 0;
  73  	}
  74  	/* scale x into [2^23, 2^24-1] */
  75  	sign = u.i>>31;
  76  	e0 = (ix>>23) - (0x7f+23);  /* e0 = ilogb(|x|)-23, positive */
  77  	u.i = ix - (e0<<23);
  78  	tx[0] = u.f;
  79  	n  =  __rem_pio2_large(tx,ty,e0,1,0);
  80  	if (sign) {
  81  		*y = -ty[0];
  82  		return -n;
  83  	}
  84  	*y = ty[0];
  85  	return n;
  86  }
  87