Home | History | Annotate | Line # | Download | only in src 
      1 /* e_remainderf.c -- float version of e_remainder.c.
      2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian (at) cygnus.com.
      3  */
      4 
      5 /*
      6  * ====================================================
      7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      8  *
      9  * Developed at SunPro, a Sun Microsystems, Inc. business.
     10  * Permission to use, copy, modify, and distribute this
     11  * software is freely granted, provided that this notice
     12  * is preserved.
     13  * ====================================================
     14  */
     15 
     16 #include <sys/cdefs.h>
     17 #if defined(LIBM_SCCS) && !defined(lint)
     18 __RCSID("$NetBSD: e_remainderf.c,v 1.7 2002/05/26 22:01:52 wiz Exp $");
     19 #endif
     20 
     21 #include "math.h"
     22 #include "math_private.h"
     23 
     24 static const float zero = 0.0;
     25 
     26 
     27 float
     28 __ieee754_remainderf(float x, float p)
     29 {
     30 	int32_t hx,hp;
     31 	u_int32_t sx;
     32 	float p_half;
     33 
     34 	GET_FLOAT_WORD(hx,x);
     35 	GET_FLOAT_WORD(hp,p);
     36 	sx = hx&0x80000000;
     37 	hp &= 0x7fffffff;
     38 	hx &= 0x7fffffff;
     39 
     40     /* purge off exception values */
     41 	if(hp==0) return (x*p)/(x*p);	 	/* p = 0 */
     42 	if((hx>=0x7f800000)||			/* x not finite */
     43 	  ((hp>0x7f800000)))			/* p is NaN */
     44 	    return (x*p)/(x*p);
     45 
     46 
     47 	if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p);	/* now x < 2p */
     48 	if ((hx-hp)==0) return zero*x;
     49 	x  = fabsf(x);
     50 	p  = fabsf(p);
     51 	if (hp<0x01000000) {
     52 	    if(x+x>p) {
     53 		x-=p;
     54 		if(x+x>=p) x -= p;
     55 	    }
     56 	} else {
     57 	    p_half = (float)0.5*p;
     58 	    if(x>p_half) {
     59 		x-=p;
     60 		if(x>=p_half) x -= p;
     61 	    }
     62 	}
     63 	GET_FLOAT_WORD(hx,x);
     64 	SET_FLOAT_WORD(x,hx^sx);
     65 	return x;
     66 }
     67