11dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* @(#)e_fmod.c 1.3 95/01/18 */ 21dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/*- 31dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ==================================================== 41dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 51dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 61dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Developed at SunSoft, a Sun Microsystems, Inc. business. 71dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Permission to use, copy, modify, and distribute this 81dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * software is freely granted, provided that this notice 91dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * is preserved. 101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ==================================================== 111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include <sys/cdefs.h> 14a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes__FBSDID("$FreeBSD$"); 15a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes 16a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes#include <float.h> 171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include "math.h" 191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include "math_private.h" 201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic const double Zero[] = {0.0, -0.0,}; 221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Return the IEEE remainder and set *quo to the last n bits of the 251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * quotient, rounded to the nearest integer. We choose n=31 because 261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * we wind up computing all the integer bits of the quotient anyway as 271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * a side-effect of computing the remainder by the shift and subtract 281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * method. In practice, this is far more bits than are needed to use 291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * remquo in reduction algorithms. 301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectdouble 321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectremquo(double x, double y, int *quo) 331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project int32_t n,hx,hy,hz,ix,iy,sx,i; 351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u_int32_t lx,ly,lz,q,sxy; 361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project EXTRACT_WORDS(hx,lx,x); 381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project EXTRACT_WORDS(hy,ly,y); 391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project sxy = (hx ^ hy) & 0x80000000; 401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project sx = hx&0x80000000; /* sign of x */ 411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hx ^=sx; /* |x| */ 421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hy &= 0x7fffffff; /* |y| */ 431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* purge off exception values */ 451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ 461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ 471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (x*y)/(x*y); 481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hx<=hy) { 491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if((hx<hy)||(lx<ly)) { 501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q = 0; 511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project goto fixup; /* |x|<|y| return x or x-y */ 521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(lx==ly) { 54a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes *quo = (sxy ? -1 : 1); 551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ 561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* determine ix = ilogb(x) */ 601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hx<0x00100000) { /* subnormal x */ 611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hx==0) { 621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; 631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; 651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else ix = (hx>>20)-1023; 671dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 681dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* determine iy = ilogb(y) */ 691dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hy<0x00100000) { /* subnormal y */ 701dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hy==0) { 711dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; 721dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 731dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; 741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 751dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else iy = (hy>>20)-1023; 761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* set up {hx,lx}, {hy,ly} and align y to x */ 781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(ix >= -1022) 791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hx = 0x00100000|(0x000fffff&hx); 801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else { /* subnormal x, shift x to normal */ 811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project n = -1022-ix; 821dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(n<=31) { 831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hx = (hx<<n)|(lx>>(32-n)); 841dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project lx <<= n; 851dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 861dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hx = lx<<(n-32); 871dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project lx = 0; 881dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(iy >= -1022) 911dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hy = 0x00100000|(0x000fffff&hy); 921dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else { /* subnormal y, shift y to normal */ 931dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project n = -1022-iy; 941dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(n<=31) { 951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hy = (hy<<n)|(ly>>(32-n)); 961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ly <<= n; 971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hy = ly<<(n-32); 991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ly = 0; 1001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* fix point fmod */ 1041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project n = ix - iy; 1051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q = 0; 1061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project while(n--) { 1071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 1081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} 1091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} 1101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q <<= 1; 1111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 1131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(hz>=0) {hx=hz;lx=lz;q++;} 1141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* convert back to floating value and restore the sign */ 1161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if((hx|lx)==0) { /* return sign(x)*0 */ 117a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes q &= 0x7fffffff; 118a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes *quo = (sxy ? -q : q); 1191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return Zero[(u_int32_t)sx>>31]; 1201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project while(hx<0x00100000) { /* normalize x */ 1221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hx = hx+hx+(lx>>31); lx = lx+lx; 1231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project iy -= 1; 1241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(iy>= -1022) { /* normalize output */ 1261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hx = ((hx-0x00100000)|((iy+1023)<<20)); 1271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { /* subnormal output */ 1281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project n = -1022 - iy; 1291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(n<=20) { 1301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project lx = (lx>>n)|((u_int32_t)hx<<(32-n)); 1311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project hx >>= n; 1321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else if (n<=31) { 133a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes lx = (hx<<(32-n))|(lx>>n); hx = 0; 1341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 135a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes lx = hx>>(n-32); hx = 0; 1361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectfixup: 1391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project INSERT_WORDS(x,hx,lx); 1401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = fabs(y); 1411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (y < 0x1p-1021) { 1421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (x+x>y || (x+x==y && (q & 1))) { 1431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q++; 1441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project x-=y; 1451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else if (x>0.5*y || (x==0.5*y && (q & 1))) { 1471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q++; 1481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project x-=y; 1491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project GET_HIGH_WORD(hx,x); 1511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project SET_HIGH_WORD(x,hx^sx); 1521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q &= 0x7fffffff; 1531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project *quo = (sxy ? -q : q); 1541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return x; 1551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 156a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes 157a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes#if LDBL_MANT_DIG == 53 158a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes__weak_reference(remquo, remquol); 159a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes#endif 160