#ifndef lintstatic	char sccsid[] = "@(#)remainder.c 1.1 92/07/30 SMI";#endif/* * Copyright (c) 1987 by Sun Microsystems, Inc. *//* remainder(x,p) * Code originated from 4.3bsd. * Modified by K.C. Ng for SUN 4.0 libm. * Return :                   * 	returns  x REM p  =  x - [x/p]*p as if in infinite precise arithmetic, *	where [x/p] is the (inifinite bit) integer nearest x/p (in half way  *	case choose the even one). * Method :  *	Based on fmod() return x-[x/p]chopped*p exactly. */#include <math.h>static double zero = 0.0;double remainder(x,p)double x,p;{double	one	= 1.0;static	long	hfmaxx[2], twom1021x[2];	double hp;	long sx;	if ((* (int *) &one) != 0) {	  hfmaxx[0] = 0x7fdfffff; hfmaxx[1] = 0xffffffff;	  twom1021x[0] = 0x00200000; twom1021x[1] = 0x0;	}						/* not a i386	*/	else {	  hfmaxx[0] = 0xffffffff; hfmaxx[1] = 0x7fdfffff;	  twom1021x[0] = 0x0; twom1021x[1] = 0x00200000;	}						/* is a i386	*/#define	hfmax (*(double *)hfmaxx)#define	twom1021 (*(double *)twom1021x)	if(p!=p) return x+p;	if(!finite(x)) return x-x;	p  = fabs(p);	if(p<=hfmax) x = fmod(x,p+p);	sx = signbit(x);	x  = fabs(x);	if (p<twom1021) {	    if(x+x>p) {		if(x==p) x=zero; else x-=p;	/* avoid x-x=-0 in RM mode */		if(x+x>=p) x -= p;	    }	} else {	    hp = 0.5*p;	    if(x>hp) {		if(x==p) x=zero; else x-=p;	/* avoid x-x=-0 in RM mode */		if(x>=hp) x -= p;	    }	}	return (sx==0)? x: -x;}double drem(x,p)double x,p;{	return(remainder(x,p));}