1 
/* * Mode: C; tabwidth: 8; indenttabsmode: nil; cbasicoffset: 4 * 
2 
* 
3 
* ***** BEGIN LICENSE BLOCK ***** 
4 
* Version: MPL 1.1/GPL 2.0/LGPL 2.1 
5 
* 
6 
* The contents of this file are subject to the Mozilla Public License Version 
7 
* 1.1 (the "License"); you may not use this file except in compliance with 
8 
* the License. You may obtain a copy of the License at 
9 
* http://www.mozilla.org/MPL/ 
10 
* 
11 
* Software distributed under the License is distributed on an "AS IS" basis, 
12 
* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 
13 
* for the specific language governing rights and limitations under the 
14 
* License. 
15 
* 
16 
* The Original Code is Mozilla Communicator client code, released 
17 
* March 31, 1998. 
18 
* 
19 
* The Initial Developer of the Original Code is 
20 
* Sun Microsystems, Inc. 
21 
* Portions created by the Initial Developer are Copyright (C) 1998 
22 
* the Initial Developer. All Rights Reserved. 
23 
* 
24 
* Contributor(s): 
25 
* 
26 
* Alternatively, the contents of this file may be used under the terms of 
27 
* either of the GNU General Public License Version 2 or later (the "GPL"), 
28 
* or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 
29 
* in which case the provisions of the GPL or the LGPL are applicable instead 
30 
* of those above. If you wish to allow use of your version of this file only 
31 
* under the terms of either the GPL or the LGPL, and not to allow others to 
32 
* use your version of this file under the terms of the MPL, indicate your 
33 
* decision by deleting the provisions above and replace them with the notice 
34 
* and other provisions required by the GPL or the LGPL. If you do not delete 
35 
* the provisions above, a recipient may use your version of this file under 
36 
* the terms of any one of the MPL, the GPL or the LGPL. 
37 
* 
38 
* ***** END LICENSE BLOCK ***** */ 
39 

40 
/* @(#)e_remainder.c 1.3 95/01/18 */ 
41 
/* 
42 
* ==================================================== 
43 
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 
44 
* 
45 
* Developed at SunSoft, a Sun Microsystems, Inc. business. 
46 
* Permission to use, copy, modify, and distribute this 
47 
* software is freely granted, provided that this notice 
48 
* is preserved. 
49 
* ==================================================== 
50 
*/ 
51 

52 
/* __ieee754_remainder(x,p) 
53 
* Return : 
54 
* returns x REM p = x  [x/p]*p as if in infinite 
55 
* precise arithmetic, where [x/p] is the (infinite bit) 
56 
* integer nearest x/p (in half way case choose the even one). 
57 
* Method : 
58 
* Based on fmod() return x[x/p]chopped*p exactlp. 
59 
*/ 
60 

61 
#include "fdlibm.h" 
62 

63 
#ifdef __STDC__ 
64 
static const double zero = 0.0; 
65 
#else 
66 
static double zero = 0.0; 
67 
#endif 
68 

69 

70 
#ifdef __STDC__ 
71 
double __ieee754_remainder(double x, double p) 
72 
#else 
73 
double __ieee754_remainder(x,p) 
74 
double x,p; 
75 
#endif 
76 
{ 
77 
fd_twoints u; 
78 
int hx,hp; 
79 
unsigned sx,lx,lp; 
80 
double p_half; 
81 

82 
u.d = x; 
83 
hx = __HI(u); /* high word of x */ 
84 
lx = __LO(u); /* low word of x */ 
85 
u.d = p; 
86 
hp = __HI(u); /* high word of p */ 
87 
lp = __LO(u); /* low word of p */ 
88 
sx = hx&0x80000000; 
89 
hp &= 0x7fffffff; 
90 
hx &= 0x7fffffff; 
91 

92 
/* purge off exception values */ 
93 
if((hplp)==0) return (x*p)/(x*p); /* p = 0 */ 
94 
if((hx>=0x7ff00000) /* x not finite */ 
95 
((hp>=0x7ff00000)&& /* p is NaN */ 
96 
(((hp0x7ff00000)lp)!=0))) 
97 
return (x*p)/(x*p); 
98 

99 

100 
if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */ 
101 
if (((hxhp)(lxlp))==0) return zero*x; 
102 
x = fd_fabs(x); 
103 
p = fd_fabs(p); 
104 
if (hp<0x00200000) { 
105 
if(x+x>p) { 
106 
x=p; 
107 
if(x+x>=p) x = p; 
108 
} 
109 
} else { 
110 
p_half = 0.5*p; 
111 
if(x>p_half) { 
112 
x=p; 
113 
if(x>=p_half) x = p; 
114 
} 
115 
} 
116 
u.d = x; 
117 
__HI(u) ^= sx; 
118 
x = u.d; 
119 
return x; 
120 
} 