Subversion Repositories shark

Rev

Rev 1618 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 1
/* e_fmodf.c -- float version of e_fmod.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_fmodf.c,v 1.2 1995/05/30 05:48:11 rgrimes Exp $";
18
#endif
19
 
20
/*
21
 * __ieee754_fmodf(x,y)
22
 * Return x mod y in exact arithmetic
23
 * Method: shift and subtract
24
 */
25
 
26
#include "math.h"
27
#include "math_private.h"
28
 
29
#ifdef __STDC__
30
static const float one = 1.0, Zero[] = {0.0, -0.0,};
31
#else
32
static float one = 1.0, Zero[] = {0.0, -0.0,};
33
#endif
34
 
35
#ifdef __STDC__
36
        float __ieee754_fmodf(float x, float y)
37
#else
38
        float __ieee754_fmodf(x,y)
39
        float x,y ;
40
#endif
41
{
42
        int32_t n,hx,hy,hz,ix,iy,sx,i;
43
 
44
        GET_FLOAT_WORD(hx,x);
45
        GET_FLOAT_WORD(hy,y);
46
        sx = hx&0x80000000;             /* sign of x */
47
        hx ^=sx;                /* |x| */
48
        hy &= 0x7fffffff;       /* |y| */
49
 
50
    /* purge off exception values */
51
        if(hy==0||(hx>=0x7f800000)||            /* y=0,or x not finite */
52
           (hy>0x7f800000))                     /* or y is NaN */
53
            return (x*y)/(x*y);
54
        if(hx<hy) return x;                     /* |x|<|y| return x */
55
        if(hx==hy)
56
            return Zero[(u_int32_t)sx>>31];     /* |x|=|y| return x*0*/
57
 
58
    /* determine ix = ilogb(x) */
59
        if(hx<0x00800000) {     /* subnormal x */
60
            for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
61
        } else ix = (hx>>23)-127;
62
 
63
    /* determine iy = ilogb(y) */
64
        if(hy<0x00800000) {     /* subnormal y */
65
            for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
66
        } else iy = (hy>>23)-127;
67
 
68
    /* set up {hx,lx}, {hy,ly} and align y to x */
69
        if(ix >= -126)
70
            hx = 0x00800000|(0x007fffff&hx);
71
        else {          /* subnormal x, shift x to normal */
72
            n = -126-ix;
73
            hx = hx<<n;
74
        }
75
        if(iy >= -126)
76
            hy = 0x00800000|(0x007fffff&hy);
77
        else {          /* subnormal y, shift y to normal */
78
            n = -126-iy;
79
            hy = hy<<n;
80
        }
81
 
82
    /* fix point fmod */
83
        n = ix - iy;
84
        while(n--) {
85
            hz=hx-hy;
86
            if(hz<0){hx = hx+hx;}
87
            else {
88
                if(hz==0)               /* return sign(x)*0 */
89
                    return Zero[(u_int32_t)sx>>31];
90
                hx = hz+hz;
91
            }
92
        }
93
        hz=hx-hy;
94
        if(hz>=0) {hx=hz;}
95
 
96
    /* convert back to floating value and restore the sign */
97
        if(hx==0)                       /* return sign(x)*0 */
98
            return Zero[(u_int32_t)sx>>31];
99
        while(hx<0x00800000) {          /* normalize x */
100
            hx = hx+hx;
101
            iy -= 1;
102
        }
103
        if(iy>= -126) {         /* normalize output */
104
            hx = ((hx-0x00800000)|((iy+127)<<23));
105
            SET_FLOAT_WORD(x,hx|sx);
106
        } else {                /* subnormal output */
107
            n = -126 - iy;
108
            hx >>= n;
109
            SET_FLOAT_WORD(x,hx|sx);
110
            x *= one;           /* create necessary signal */
111
        }
112
        return x;               /* exact output */
113
}