1 
2 /* @(#)e_remainder.c 5.1 93/09/24 */
3 /*
4  * ====================================================
5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6  *
7  * Developed at SunPro, a Sun Microsystems, Inc. business.
8  * Permission to use, copy, modify, and distribute this
9  * software is freely granted, provided that this notice
10  * is preserved.
11  * ====================================================
12  */
13 
14 /*
15 FUNCTION
16 <<remainder>>, <<remainderf>>---round and  remainder
17 INDEX
18         remainder
19 INDEX
20         remainderf
21 
22 ANSI_SYNOPSIS
23         #include <math.h>
24         double remainder(double <[x]>, double <[y]>);
25         float remainderf(float <[x]>, float <[y]>);
26 
27 TRAD_SYNOPSIS
28         #include <math.h>
29         double remainder(<[x]>,<[y]>)
30         double <[x]>, <[y]>;
31         float remainderf(<[x]>,<[y]>)
32         float <[x]>, <[y]>;
33 
34 DESCRIPTION
35 <<remainder>> and <<remainderf>> find the remainder of
36 <[x]>/<[y]>; this value is in the range -<[y]>/2 .. +<[y]>/2.
37 
38 RETURNS
39 <<remainder>> returns the integer result as a double.
40 
41 PORTABILITY
42 <<remainder>> is a System V release 4.
43 <<remainderf>> is an extension.
44 
45 */
46 
47 /* remainder(x,p)
48  * Return :
49  * 	returns  x REM p  =  x - [x/p]*p as if in infinite
50  * 	precise arithmetic, where [x/p] is the (infinite bit)
51  *	integer nearest x/p (in half way case choose the even one).
52  * Method :
53  *	Based on fmod() return x-[x/p]chopped*p exactlp.
54  */
55 
56 #include "fdlibm.h"
57 
58 #ifndef _DOUBLE_IS_32BITS
59 
60 #ifdef __STDC__
61 static const double zero = 0.0;
62 #else
63 static double zero = 0.0;
64 #endif
65 
66 
67 #ifdef __STDC__
remainder(double x,double p)68 	double remainder(double x, double p)
69 #else
70 	double remainder(x,p)
71 	double x,p;
72 #endif
73 {
74 	__int32_t hx,hp;
75 	__uint32_t sx,lx,lp;
76 	double p_half;
77 
78 	EXTRACT_WORDS(hx,lx,x);
79 	EXTRACT_WORDS(hp,lp,p);
80 	sx = hx&0x80000000;
81 	hp &= 0x7fffffff;
82 	hx &= 0x7fffffff;
83 
84     /* purge off exception values */
85 	if((hp|lp)==0) return (x*p)/(x*p); 	/* p = 0 */
86 	if((hx>=0x7ff00000)||			/* x not finite */
87 	  ((hp>=0x7ff00000)&&			/* p is NaN */
88 	  (((hp-0x7ff00000)|lp)!=0)))
89 	    return (x*p)/(x*p);
90 
91 
92 	if (hp<=0x7fdfffff) x = fmod(x,p+p);	/* now x < 2p */
93 	if (((hx-hp)|(lx-lp))==0) return zero*x;
94 	x  = fabs(x);
95 	p  = fabs(p);
96 	if (hp<0x00200000) {
97 	    if(x+x>p) {
98 		x-=p;
99 		if(x+x>=p) x -= p;
100 	    }
101 	} else {
102 	    p_half = 0.5*p;
103 	    if(x>p_half) {
104 		x-=p;
105 		if(x>=p_half) x -= p;
106 	    }
107 	}
108 	GET_HIGH_WORD(hx,x);
109 	SET_HIGH_WORD(x,hx^sx);
110 	return x;
111 }
112 
113 #endif /* defined(_DOUBLE_IS_32BITS) */
114