1 /*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunSoft, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12 #include "bsd_private.h"
13
bsd__ieee754_remainder(double x,double p)14 double bsd__ieee754_remainder(double x, double p)
15 {
16 int32_t hx,hp;
17 u_int32_t sx,lx,lp;
18 double p_half;
19
20 EXTRACT_WORDS(hx,lx,x);
21 EXTRACT_WORDS(hp,lp,p);
22 sx = hx&0x80000000;
23 hp &= 0x7fffffff;
24 hx &= 0x7fffffff;
25
26 /* purge off exception values */
27 if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */
28 if((hx>=0x7ff00000)|| /* x not finite */
29 ((hp>=0x7ff00000)&& /* p is NaN */
30 (((hp-0x7ff00000)|lp)!=0)))
31 return ((long double)x*p)/((long double)x*p);
32
33
34 if (hp<=0x7fdfffff) x = bsd__ieee754_fmod(x,p+p); /* now x < 2p */
35 if (((hx-hp)|(lx-lp))==0) return zero*x;
36 x = fabs(x);
37 p = fabs(p);
38 if (hp<0x00200000) {
39 if(x+x>p) {
40 x-=p;
41 if(x+x>=p) x -= p;
42 }
43 } else {
44 p_half = 0.5*p;
45 if(x>p_half) {
46 x-=p;
47 if(x>=p_half) x -= p;
48 }
49 }
50 GET_HIGH_WORD(hx,x);
51 if ((hx&0x7fffffff)==0) hx = 0;
52 SET_HIGH_WORD(x,hx^sx);
53 return x;
54 }
55