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