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