Subversion Repositories shark

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 pj 1
/* @(#)e_remainder.c 5.1 93/09/24 */
2
/*
3
 * ====================================================
4
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5
 *
6
 * Developed at SunPro, a Sun Microsystems, Inc. business.
7
 * Permission to use, copy, modify, and distribute this
8
 * software is freely granted, provided that this notice
9
 * is preserved.
10
 * ====================================================
11
 */
12
 
13
#ifndef lint
14
static char rcsid[] = "$\Id: e_remainder.c,v 1.2.6.1 1997/02/23 11:03:07 joerg Exp $";
15
#endif
16
 
17
/* __ieee754_remainder(x,p)
18
 * Return :
19
 *      returns  x REM p  =  x - [x/p]*p as if in infinite
20
 *      precise arithmetic, where [x/p] is the (infinite bit)
21
 *      integer nearest x/p (in half way case choose the even one).
22
 * Method :
23
 *      Based on fmod() return x-[x/p]chopped*p exactlp.
24
 */
25
 
26
#include "math.h"
27
#include "math_private.h"
28
 
29
#ifdef __STDC__
30
static const double zero = 0.0;
31
#else
32
static double zero = 0.0;
33
#endif
34
 
35
 
36
#ifdef __STDC__
37
        double __generic___ieee754_remainder(double x, double p)
38
#else
39
        double __generic___ieee754_remainder(x,p)
40
        double x,p;
41
#endif
42
{
43
        int32_t hx,hp;
44
        u_int32_t sx,lx,lp;
45
        double p_half;
46
 
47
        EXTRACT_WORDS(hx,lx,x);
48
        EXTRACT_WORDS(hp,lp,p);
49
        sx = hx&0x80000000;
50
        hp &= 0x7fffffff;
51
        hx &= 0x7fffffff;
52
 
53
    /* purge off exception values */
54
        if((hp|lp)==0) return (x*p)/(x*p);      /* p = 0 */
55
        if((hx>=0x7ff00000)||                   /* x not finite */
56
          ((hp>=0x7ff00000)&&                   /* p is NaN */
57
          (((hp-0x7ff00000)|lp)!=0)))
58
            return (x*p)/(x*p);
59
 
60
 
61
        if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);  /* now x < 2p */
62
        if (((hx-hp)|(lx-lp))==0) return zero*x;
63
        x  = fabs(x);
64
        p  = fabs(p);
65
        if (hp<0x00200000) {
66
            if(x+x>p) {
67
                x-=p;
68
                if(x+x>=p) x -= p;
69
            }
70
        } else {
71
            p_half = 0.5*p;
72
            if(x>p_half) {
73
                x-=p;
74
                if(x>=p_half) x -= p;
75
            }
76
        }
77
        GET_HIGH_WORD(hx,x);
78
        SET_HIGH_WORD(x,hx^sx);
79
        return x;
80
}