Back to index

glibc  2.9
e_remainderf.c
Go to the documentation of this file.
00001 /* e_remainderf.c -- float version of e_remainder.c.
00002  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice 
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: e_remainderf.c,v 1.4 1995/05/10 20:46:08 jtc Exp $";
00018 #endif
00019 
00020 #include "math.h"
00021 #include "math_private.h"
00022 
00023 #ifdef __STDC__
00024 static const float zero = 0.0;
00025 #else
00026 static float zero = 0.0;
00027 #endif
00028 
00029 
00030 #ifdef __STDC__
00031        float __ieee754_remainderf(float x, float p)
00032 #else
00033        float __ieee754_remainderf(x,p)
00034        float x,p;
00035 #endif
00036 {
00037        int32_t hx,hp;
00038        u_int32_t sx;
00039        float p_half;
00040 
00041        GET_FLOAT_WORD(hx,x);
00042        GET_FLOAT_WORD(hp,p);
00043        sx = hx&0x80000000;
00044        hp &= 0x7fffffff;
00045        hx &= 0x7fffffff;
00046 
00047     /* purge off exception values */
00048        if(hp==0) return (x*p)/(x*p);             /* p = 0 */
00049        if((hx>=0x7f800000)||                     /* x not finite */
00050          ((hp>0x7f800000)))               /* p is NaN */
00051            return (x*p)/(x*p);
00052 
00053 
00054        if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p);  /* now x < 2p */
00055        if ((hx-hp)==0) return zero*x;
00056        x  = fabsf(x);
00057        p  = fabsf(p);
00058        if (hp<0x01000000) {
00059            if(x+x>p) {
00060               x-=p;
00061               if(x+x>=p) x -= p;
00062            }
00063        } else {
00064            p_half = (float)0.5*p;
00065            if(x>p_half) {
00066               x-=p;
00067               if(x>=p_half) x -= p;
00068            }
00069        }
00070        GET_FLOAT_WORD(hx,x);
00071        SET_FLOAT_WORD(x,hx^sx);
00072        return x;
00073 }