12aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm/* @(#)e_fmod.c 5.1 93/09/24 */ 22aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm/* 32aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * ==================================================== 42aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 52aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * 62aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * Developed at SunPro, a Sun Microsystems, Inc. business. 72aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * Permission to use, copy, modify, and distribute this 82aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * software is freely granted, provided that this notice 92aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * is preserved. 102aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * ==================================================== 112aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm */ 122aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm#include <LibConfig.h> 132aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm#include <sys/EfiCdefs.h> 142aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm#if defined(LIBM_SCCS) && !defined(lint) 152aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm__RCSID("$NetBSD: e_fmod.c,v 1.11 2002/05/26 22:01:49 wiz Exp $"); 162aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm#endif 172aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 182aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm/* 192aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * __ieee754_fmod(x,y) 202aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * Return x mod y in exact arithmetic 212aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm * Method: shift and subtract 222aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm */ 232aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 242aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm#include "math.h" 252aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm#include "math_private.h" 262aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 272aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */ 282aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm // unary minus operator applied to unsigned type, result still unsigned 292aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm #pragma warning ( disable : 4146 ) 302aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm#endif 312aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 322aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylmstatic const double one = 1.0, Zero[] = {0.0, -0.0,}; 332aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 342aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylmdouble 352aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm__ieee754_fmod(double x, double y) 362aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm{ 372aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm int32_t n,hx,hy,hz,ix,iy,sx,i; 382aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm u_int32_t lx,ly,lz; 392aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 402aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm EXTRACT_WORDS(hx,lx,x); 412aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm EXTRACT_WORDS(hy,ly,y); 422aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm sx = hx&0x80000000; /* sign of x */ 432aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hx ^=sx; /* |x| */ 442aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hy &= 0x7fffffff; /* |y| */ 452aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 462aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm /* purge off exception values */ 472aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ 482aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ 492aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm return (x*y)/(x*y); 502aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(hx<=hy) { 512aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */ 522aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(lx==ly) 532aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ 542aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 552aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 562aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm /* determine ix = ilogb(x) */ 572aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(hx<0x00100000) { /* subnormal x */ 582aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(hx==0) { 592aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; 602aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else { 612aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; 622aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 632aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else ix = (hx>>20)-1023; 642aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 652aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm /* determine iy = ilogb(y) */ 662aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(hy<0x00100000) { /* subnormal y */ 672aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(hy==0) { 682aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; 692aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else { 702aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; 712aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 722aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else iy = (hy>>20)-1023; 732aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 742aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm /* set up {hx,lx}, {hy,ly} and align y to x */ 752aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(ix >= -1022) 762aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hx = 0x00100000|(0x000fffff&hx); 772aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm else { /* subnormal x, shift x to normal */ 782aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm n = -1022-ix; 792aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(n<=31) { 802aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hx = (hx<<n)|(lx>>(32-n)); 812aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm lx <<= n; 822aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else { 832aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hx = lx<<(n-32); 842aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm lx = 0; 852aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 862aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 872aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(iy >= -1022) 882aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hy = 0x00100000|(0x000fffff&hy); 892aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm else { /* subnormal y, shift y to normal */ 902aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm n = -1022-iy; 912aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(n<=31) { 922aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hy = (hy<<n)|(ly>>(32-n)); 932aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm ly <<= n; 942aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else { 952aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hy = ly<<(n-32); 962aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm ly = 0; 972aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 982aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 992aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 1002aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm /* fix point fmod */ 1012aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm n = ix - iy; 1022aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm while(n--) { 1032aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 1042aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} 1052aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm else { 1062aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if((hz|lz)==0) /* return sign(x)*0 */ 1072aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm return Zero[(u_int32_t)sx>>31]; 1082aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hx = hz+hz+(lz>>31); lx = lz+lz; 1092aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 1102aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 1112aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 1122aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(hz>=0) {hx=hz;lx=lz;} 1132aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm 1142aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm /* convert back to floating value and restore the sign */ 1152aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if((hx|lx)==0) /* return sign(x)*0 */ 1162aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm return Zero[(u_int32_t)sx>>31]; 1172aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm while(hx<0x00100000) { /* normalize x */ 1182aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hx = hx+hx+(lx>>31); lx = lx+lx; 1192aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm iy -= 1; 1202aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 1212aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(iy>= -1022) { /* normalize output */ 1222aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hx = ((hx-0x00100000)|((iy+1023)<<20)); 1232aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm INSERT_WORDS(x,hx|sx,lx); 1242aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else { /* subnormal output */ 1252aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm n = -1022 - iy; 1262aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm if(n<=20) { 1272aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm lx = (lx>>n)|((u_int32_t)hx<<(32-n)); 1282aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm hx >>= n; 1292aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else if (n<=31) { 1302aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm lx = (hx<<(32-n))|(lx>>n); hx = sx; 1312aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } else { 1322aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm lx = hx>>(n-32); hx = sx; 1332aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 1342aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm INSERT_WORDS(x,hx|sx,lx); 1352aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm x *= one; /* create necessary signal */ 1362aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm } 1372aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm return x; /* exact output */ 1382aa62f2bc9a9654687b377d9ca8a8c2c860a3852darylm} 139