s_remquo.c revision 144091
1144091Sdas/* @(#)e_fmod.c 1.3 95/01/18 */ 2144091Sdas/*- 3144091Sdas * ==================================================== 4144091Sdas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5144091Sdas * 6144091Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 7144091Sdas * Permission to use, copy, modify, and distribute this 8144091Sdas * software is freely granted, provided that this notice 9144091Sdas * is preserved. 10144091Sdas * ==================================================== 11144091Sdas */ 12144091Sdas 13144091Sdas#include <sys/cdefs.h> 14144091Sdas__FBSDID("$FreeBSD: head/lib/msun/src/s_remquo.c 144091 2005-03-25 04:40:44Z das $"); 15144091Sdas 16144091Sdas#include "math.h" 17144091Sdas#include "math_private.h" 18144091Sdas 19144091Sdasstatic const double Zero[] = {0.0, -0.0,}; 20144091Sdas 21144091Sdas/* 22144091Sdas * Return the IEEE remainder and set *quo to the last n bits of the 23144091Sdas * quotient, rounded to the nearest integer. We choose n=31 because 24144091Sdas * we wind up computing all the integer bits of the quotient anyway as 25144091Sdas * a side-effect of computing the remainder by the shift and subtract 26144091Sdas * method. In practice, this is far more bits than are needed to use 27144091Sdas * remquo in reduction algorithms. 28144091Sdas */ 29144091Sdasdouble 30144091Sdasremquo(double x, double y, int *quo) 31144091Sdas{ 32144091Sdas int32_t n,hx,hy,hz,ix,iy,sx,i; 33144091Sdas u_int32_t lx,ly,lz,q,sxy; 34144091Sdas 35144091Sdas EXTRACT_WORDS(hx,lx,x); 36144091Sdas EXTRACT_WORDS(hy,ly,y); 37144091Sdas sxy = (hx ^ hy) & 0x80000000; 38144091Sdas sx = hx&0x80000000; /* sign of x */ 39144091Sdas hx ^=sx; /* |x| */ 40144091Sdas hy &= 0x7fffffff; /* |y| */ 41144091Sdas 42144091Sdas /* purge off exception values */ 43144091Sdas if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ 44144091Sdas ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ 45144091Sdas return (x*y)/(x*y); 46144091Sdas if(hx<=hy) { 47144091Sdas if((hx<hy)||(lx<ly)) { 48144091Sdas q = 0; 49144091Sdas goto fixup; /* |x|<|y| return x or x-y */ 50144091Sdas } 51144091Sdas if(lx==ly) { 52144091Sdas *quo = 1; 53144091Sdas return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ 54144091Sdas } 55144091Sdas } 56144091Sdas 57144091Sdas /* determine ix = ilogb(x) */ 58144091Sdas if(hx<0x00100000) { /* subnormal x */ 59144091Sdas if(hx==0) { 60144091Sdas for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; 61144091Sdas } else { 62144091Sdas for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; 63144091Sdas } 64144091Sdas } else ix = (hx>>20)-1023; 65144091Sdas 66144091Sdas /* determine iy = ilogb(y) */ 67144091Sdas if(hy<0x00100000) { /* subnormal y */ 68144091Sdas if(hy==0) { 69144091Sdas for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; 70144091Sdas } else { 71144091Sdas for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; 72144091Sdas } 73144091Sdas } else iy = (hy>>20)-1023; 74144091Sdas 75144091Sdas /* set up {hx,lx}, {hy,ly} and align y to x */ 76144091Sdas if(ix >= -1022) 77144091Sdas hx = 0x00100000|(0x000fffff&hx); 78144091Sdas else { /* subnormal x, shift x to normal */ 79144091Sdas n = -1022-ix; 80144091Sdas if(n<=31) { 81144091Sdas hx = (hx<<n)|(lx>>(32-n)); 82144091Sdas lx <<= n; 83144091Sdas } else { 84144091Sdas hx = lx<<(n-32); 85144091Sdas lx = 0; 86144091Sdas } 87144091Sdas } 88144091Sdas if(iy >= -1022) 89144091Sdas hy = 0x00100000|(0x000fffff&hy); 90144091Sdas else { /* subnormal y, shift y to normal */ 91144091Sdas n = -1022-iy; 92144091Sdas if(n<=31) { 93144091Sdas hy = (hy<<n)|(ly>>(32-n)); 94144091Sdas ly <<= n; 95144091Sdas } else { 96144091Sdas hy = ly<<(n-32); 97144091Sdas ly = 0; 98144091Sdas } 99144091Sdas } 100144091Sdas 101144091Sdas /* fix point fmod */ 102144091Sdas n = ix - iy; 103144091Sdas q = 0; 104144091Sdas while(n--) { 105144091Sdas hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 106144091Sdas if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} 107144091Sdas else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} 108144091Sdas q <<= 1; 109144091Sdas } 110144091Sdas hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 111144091Sdas if(hz>=0) {hx=hz;lx=lz;q++;} 112144091Sdas 113144091Sdas /* convert back to floating value and restore the sign */ 114144091Sdas if((hx|lx)==0) { /* return sign(x)*0 */ 115144091Sdas *quo = (sxy ? -q : q); 116144091Sdas return Zero[(u_int32_t)sx>>31]; 117144091Sdas } 118144091Sdas while(hx<0x00100000) { /* normalize x */ 119144091Sdas hx = hx+hx+(lx>>31); lx = lx+lx; 120144091Sdas iy -= 1; 121144091Sdas } 122144091Sdas if(iy>= -1022) { /* normalize output */ 123144091Sdas hx = ((hx-0x00100000)|((iy+1023)<<20)); 124144091Sdas } else { /* subnormal output */ 125144091Sdas n = -1022 - iy; 126144091Sdas if(n<=20) { 127144091Sdas lx = (lx>>n)|((u_int32_t)hx<<(32-n)); 128144091Sdas hx >>= n; 129144091Sdas } else if (n<=31) { 130144091Sdas lx = (hx<<(32-n))|(lx>>n); hx = sx; 131144091Sdas } else { 132144091Sdas lx = hx>>(n-32); hx = sx; 133144091Sdas } 134144091Sdas } 135144091Sdasfixup: 136144091Sdas INSERT_WORDS(x,hx,lx); 137144091Sdas y = fabs(y); 138144091Sdas if (y < 0x1p-1021) { 139144091Sdas if (x+x>y || (x+x==y && (q & 1))) { 140144091Sdas q++; 141144091Sdas x-=y; 142144091Sdas } 143144091Sdas } else if (x>0.5*y || (x==0.5*y && (q & 1))) { 144144091Sdas q++; 145144091Sdas x-=y; 146144091Sdas } 147144091Sdas GET_HIGH_WORD(hx,x); 148144091Sdas SET_HIGH_WORD(x,hx^sx); 149144091Sdas q &= 0x7fffffff; 150144091Sdas *quo = (sxy ? -q : q); 151144091Sdas return x; 152144091Sdas} 153