1141296Sdas 2141296Sdas/* @(#)e_remainder.c 1.3 95/01/18 */ 32116Sjkh/* 42116Sjkh * ==================================================== 52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 62116Sjkh * 7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 82116Sjkh * Permission to use, copy, modify, and distribute this 9141296Sdas * software is freely granted, provided that this notice 102116Sjkh * is preserved. 112116Sjkh * ==================================================== 122116Sjkh */ 132116Sjkh 14176207Sbde#include <sys/cdefs.h> 15176207Sbde__FBSDID("$FreeBSD$"); 162116Sjkh 172116Sjkh/* __ieee754_remainder(x,p) 18141296Sdas * Return : 19141296Sdas * returns x REM p = x - [x/p]*p as if in infinite 20141296Sdas * precise arithmetic, where [x/p] is the (infinite bit) 212116Sjkh * integer nearest x/p (in half way case choose the even one). 22141296Sdas * Method : 232116Sjkh * Based on fmod() return x-[x/p]chopped*p exactlp. 242116Sjkh */ 252116Sjkh 26177765Sdas#include <float.h> 27177765Sdas 282116Sjkh#include "math.h" 292116Sjkh#include "math_private.h" 302116Sjkh 312116Sjkhstatic const double zero = 0.0; 322116Sjkh 332116Sjkh 3497413Salfreddouble 35117912Speter__ieee754_remainder(double x, double p) 362116Sjkh{ 372116Sjkh int32_t hx,hp; 382116Sjkh u_int32_t sx,lx,lp; 392116Sjkh double p_half; 402116Sjkh 412116Sjkh EXTRACT_WORDS(hx,lx,x); 422116Sjkh EXTRACT_WORDS(hp,lp,p); 432116Sjkh sx = hx&0x80000000; 442116Sjkh hp &= 0x7fffffff; 452116Sjkh hx &= 0x7fffffff; 462116Sjkh 472116Sjkh /* purge off exception values */ 482116Sjkh if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */ 492116Sjkh if((hx>=0x7ff00000)|| /* x not finite */ 502116Sjkh ((hp>=0x7ff00000)&& /* p is NaN */ 512116Sjkh (((hp-0x7ff00000)|lp)!=0))) 52176207Sbde return ((long double)x*p)/((long double)x*p); 532116Sjkh 542116Sjkh 552116Sjkh if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */ 562116Sjkh if (((hx-hp)|(lx-lp))==0) return zero*x; 572116Sjkh x = fabs(x); 582116Sjkh p = fabs(p); 592116Sjkh if (hp<0x00200000) { 602116Sjkh if(x+x>p) { 612116Sjkh x-=p; 622116Sjkh if(x+x>=p) x -= p; 632116Sjkh } 642116Sjkh } else { 652116Sjkh p_half = 0.5*p; 662116Sjkh if(x>p_half) { 672116Sjkh x-=p; 682116Sjkh if(x>=p_half) x -= p; 692116Sjkh } 702116Sjkh } 712116Sjkh GET_HIGH_WORD(hx,x); 72176207Sbde if ((hx&0x7fffffff)==0) hx = 0; 732116Sjkh SET_HIGH_WORD(x,hx^sx); 742116Sjkh return x; 752116Sjkh} 76177765Sdas 77177765Sdas#if LDBL_MANT_DIG == 53 78177765Sdas__weak_reference(remainder, remainderl); 79177765Sdas#endif 80