e_remainderf.c - vx32 - Local 9vx git repository for patches.
 (HTM) git clone git://r-36.net/vx32
 (DIR) Log
 (DIR) Files
 (DIR) Refs
       ---
       e_remainderf.c (1479B)
       ---
            1 /* e_remainderf.c -- float version of e_remainder.c.
            2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@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 #ifndef lint
           17 static char rcsid[] = "$FreeBSD: src/lib/msun/src/e_remainderf.c,v 1.7 2002/05/28 18:15:04 alfred Exp $";
           18 #endif
           19 
           20 #include "math.h"
           21 #include "math_private.h"
           22 
           23 static const float zero = 0.0;
           24 
           25 
           26 float
           27 __ieee754_remainderf(float x, float p)
           28 {
           29         int32_t hx,hp;
           30         u_int32_t sx;
           31         float p_half;
           32 
           33         GET_FLOAT_WORD(hx,x);
           34         GET_FLOAT_WORD(hp,p);
           35         sx = hx&0x80000000;
           36         hp &= 0x7fffffff;
           37         hx &= 0x7fffffff;
           38 
           39     /* purge off exception values */
           40         if(hp==0) return (x*p)/(x*p);                 /* p = 0 */
           41         if((hx>=0x7f800000)||                        /* x not finite */
           42           ((hp>0x7f800000)))                        /* p is NaN */
           43             return (x*p)/(x*p);
           44 
           45 
           46         if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p);        /* now x < 2p */
           47         if ((hx-hp)==0) return zero*x;
           48         x  = fabsf(x);
           49         p  = fabsf(p);
           50         if (hp<0x01000000) {
           51             if(x+x>p) {
           52                 x-=p;
           53                 if(x+x>=p) x -= p;
           54             }
           55         } else {
           56             p_half = (float)0.5*p;
           57             if(x>p_half) {
           58                 x-=p;
           59                 if(x>=p_half) x -= p;
           60             }
           61         }
           62         GET_FLOAT_WORD(hx,x);
           63         SET_FLOAT_WORD(x,hx^sx);
           64         return x;
           65 }