Subversion Repositories shark

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 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[] = "$\Id: e_remainderf.c,v 1.2 1995/05/30 05:48:40 rgrimes Exp $";
18
#endif
19
 
20
#include "math.h"
21
#include "math_private.h"
22
 
23
#ifdef __STDC__
24
static const float zero = 0.0;
25
#else
26
static float zero = 0.0;
27
#endif
28
 
29
 
30
#ifdef __STDC__
31
        float __ieee754_remainderf(float x, float p)
32
#else
33
        float __ieee754_remainderf(x,p)
34
        float x,p;
35
#endif
36
{
37
        int32_t hx,hp;
38
        u_int32_t sx;
39
        float p_half;
40
 
41
        GET_FLOAT_WORD(hx,x);
42
        GET_FLOAT_WORD(hp,p);
43
        sx = hx&0x80000000;
44
        hp &= 0x7fffffff;
45
        hx &= 0x7fffffff;
46
 
47
    /* purge off exception values */
48
        if(hp==0) return (x*p)/(x*p);           /* p = 0 */
49
        if((hx>=0x7f800000)||                   /* x not finite */
50
          ((hp>0x7f800000)))                    /* p is NaN */
51
            return (x*p)/(x*p);
52
 
53
 
54
        if (hp<=0x7effffff) x = __ieee754_fmodf(x,p+p); /* now x < 2p */
55
        if ((hx-hp)==0) return zero*x;
56
        x  = fabsf(x);
57
        p  = fabsf(p);
58
        if (hp<0x01000000) {
59
            if(x+x>p) {
60
                x-=p;
61
                if(x+x>=p) x -= p;
62
            }
63
        } else {
64
            p_half = (float)0.5*p;
65
            if(x>p_half) {
66
                x-=p;
67
                if(x>=p_half) x -= p;
68
            }
69
        }
70
        GET_FLOAT_WORD(hx,x);
71
        SET_FLOAT_WORD(x,hx^sx);
72
        return x;
73
}