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