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$"); 15144091Sdas 16177764Sdas#include <float.h> 17177764Sdas 18144091Sdas#include "math.h" 19144091Sdas#include "math_private.h" 20144091Sdas 21144091Sdasstatic const double Zero[] = {0.0, -0.0,}; 22144091Sdas 23144091Sdas/* 24144091Sdas * Return the IEEE remainder and set *quo to the last n bits of the 25144091Sdas * quotient, rounded to the nearest integer. We choose n=31 because 26144091Sdas * we wind up computing all the integer bits of the quotient anyway as 27144091Sdas * a side-effect of computing the remainder by the shift and subtract 28144091Sdas * method. In practice, this is far more bits than are needed to use 29144091Sdas * remquo in reduction algorithms. 30144091Sdas */ 31144091Sdasdouble 32144091Sdasremquo(double x, double y, int *quo) 33144091Sdas{ 34144091Sdas int32_t n,hx,hy,hz,ix,iy,sx,i; 35144091Sdas u_int32_t lx,ly,lz,q,sxy; 36144091Sdas 37144091Sdas EXTRACT_WORDS(hx,lx,x); 38144091Sdas EXTRACT_WORDS(hy,ly,y); 39144091Sdas sxy = (hx ^ hy) & 0x80000000; 40144091Sdas sx = hx&0x80000000; /* sign of x */ 41144091Sdas hx ^=sx; /* |x| */ 42144091Sdas hy &= 0x7fffffff; /* |y| */ 43144091Sdas 44144091Sdas /* purge off exception values */ 45144091Sdas if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ 46144091Sdas ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ 47144091Sdas return (x*y)/(x*y); 48144091Sdas if(hx<=hy) { 49144091Sdas if((hx<hy)||(lx<ly)) { 50144091Sdas q = 0; 51144091Sdas goto fixup; /* |x|<|y| return x or x-y */ 52144091Sdas } 53144091Sdas if(lx==ly) { 54233973Sdas *quo = (sxy ? -1 : 1); 55144091Sdas return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ 56144091Sdas } 57144091Sdas } 58144091Sdas 59144091Sdas /* determine ix = ilogb(x) */ 60144091Sdas if(hx<0x00100000) { /* subnormal x */ 61144091Sdas if(hx==0) { 62144091Sdas for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; 63144091Sdas } else { 64144091Sdas for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; 65144091Sdas } 66144091Sdas } else ix = (hx>>20)-1023; 67144091Sdas 68144091Sdas /* determine iy = ilogb(y) */ 69144091Sdas if(hy<0x00100000) { /* subnormal y */ 70144091Sdas if(hy==0) { 71144091Sdas for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; 72144091Sdas } else { 73144091Sdas for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; 74144091Sdas } 75144091Sdas } else iy = (hy>>20)-1023; 76144091Sdas 77144091Sdas /* set up {hx,lx}, {hy,ly} and align y to x */ 78144091Sdas if(ix >= -1022) 79144091Sdas hx = 0x00100000|(0x000fffff&hx); 80144091Sdas else { /* subnormal x, shift x to normal */ 81144091Sdas n = -1022-ix; 82144091Sdas if(n<=31) { 83144091Sdas hx = (hx<<n)|(lx>>(32-n)); 84144091Sdas lx <<= n; 85144091Sdas } else { 86144091Sdas hx = lx<<(n-32); 87144091Sdas lx = 0; 88144091Sdas } 89144091Sdas } 90144091Sdas if(iy >= -1022) 91144091Sdas hy = 0x00100000|(0x000fffff&hy); 92144091Sdas else { /* subnormal y, shift y to normal */ 93144091Sdas n = -1022-iy; 94144091Sdas if(n<=31) { 95144091Sdas hy = (hy<<n)|(ly>>(32-n)); 96144091Sdas ly <<= n; 97144091Sdas } else { 98144091Sdas hy = ly<<(n-32); 99144091Sdas ly = 0; 100144091Sdas } 101144091Sdas } 102144091Sdas 103144091Sdas /* fix point fmod */ 104144091Sdas n = ix - iy; 105144091Sdas q = 0; 106144091Sdas while(n--) { 107144091Sdas hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 108144091Sdas if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} 109144091Sdas else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} 110144091Sdas q <<= 1; 111144091Sdas } 112144091Sdas hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 113144091Sdas if(hz>=0) {hx=hz;lx=lz;q++;} 114144091Sdas 115144091Sdas /* convert back to floating value and restore the sign */ 116144091Sdas if((hx|lx)==0) { /* return sign(x)*0 */ 117233973Sdas q &= 0x7fffffff; 118144091Sdas *quo = (sxy ? -q : q); 119144091Sdas return Zero[(u_int32_t)sx>>31]; 120144091Sdas } 121144091Sdas while(hx<0x00100000) { /* normalize x */ 122144091Sdas hx = hx+hx+(lx>>31); lx = lx+lx; 123144091Sdas iy -= 1; 124144091Sdas } 125144091Sdas if(iy>= -1022) { /* normalize output */ 126144091Sdas hx = ((hx-0x00100000)|((iy+1023)<<20)); 127144091Sdas } else { /* subnormal output */ 128144091Sdas n = -1022 - iy; 129144091Sdas if(n<=20) { 130144091Sdas lx = (lx>>n)|((u_int32_t)hx<<(32-n)); 131144091Sdas hx >>= n; 132144091Sdas } else if (n<=31) { 133233973Sdas lx = (hx<<(32-n))|(lx>>n); hx = 0; 134144091Sdas } else { 135233973Sdas lx = hx>>(n-32); hx = 0; 136144091Sdas } 137144091Sdas } 138144091Sdasfixup: 139144091Sdas INSERT_WORDS(x,hx,lx); 140144091Sdas y = fabs(y); 141144091Sdas if (y < 0x1p-1021) { 142144091Sdas if (x+x>y || (x+x==y && (q & 1))) { 143144091Sdas q++; 144144091Sdas x-=y; 145144091Sdas } 146144091Sdas } else if (x>0.5*y || (x==0.5*y && (q & 1))) { 147144091Sdas q++; 148144091Sdas x-=y; 149144091Sdas } 150144091Sdas GET_HIGH_WORD(hx,x); 151144091Sdas SET_HIGH_WORD(x,hx^sx); 152144091Sdas q &= 0x7fffffff; 153144091Sdas *quo = (sxy ? -q : q); 154144091Sdas return x; 155144091Sdas} 156177764Sdas 157177764Sdas#if LDBL_MANT_DIG == 53 158177764Sdas__weak_reference(remquo, remquol); 159177764Sdas#endif 160