1  /* e_remainderf.c -- float version of e_remainder.c.
       2   */
       3  
       4  /*
       5   * ====================================================
       6   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       7   *
       8   * Developed at SunPro, a Sun Microsystems, Inc. business.
       9   * Permission to use, copy, modify, and distribute this
      10   * software is freely granted, provided that this notice
      11   * is preserved.
      12   * ====================================================
      13   */
      14  
      15  #include <math.h>
      16  #include <math_private.h>
      17  #include <libm-alias-finite.h>
      18  
      19  static const float zero = 0.0;
      20  
      21  
      22  float
      23  __ieee754_remainderf(float x, float p)
      24  {
      25  	int32_t hx,hp;
      26  	uint32_t sx;
      27  	float p_half;
      28  
      29  	GET_FLOAT_WORD(hx,x);
      30  	GET_FLOAT_WORD(hp,p);
      31  	sx = hx&0x80000000;
      32  	hp &= 0x7fffffff;
      33  	hx &= 0x7fffffff;
      34  
      35      /* purge off exception values */
      36  	if(hp==0) return (x*p)/(x*p);		/* p = 0 */
      37  	if((hx>=0x7f800000)||			/* x not finite */
      38  	  ((hp>0x7f800000)))			/* p is NaN */
      39  	    return (x*p)/(x*p);
      40  
      41  
      42  	if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p);	/* now x < 2p */
      43  	if ((hx-hp)==0) return zero*x;
      44  	x  = fabsf(x);
      45  	p  = fabsf(p);
      46  	if (hp<0x01000000) {
      47  	    if(x+x>p) {
      48  		x-=p;
      49  		if(x+x>=p) x -= p;
      50  	    }
      51  	} else {
      52  	    p_half = (float)0.5*p;
      53  	    if(x>p_half) {
      54  		x-=p;
      55  		if(x>=p_half) x -= p;
      56  	    }
      57  	}
      58  	GET_FLOAT_WORD(hx,x);
      59  	SET_FLOAT_WORD(x,hx^sx);
      60  	return x;
      61  }
      62  libm_alias_finite (__ieee754_remainderf, __remainderf)