1 /* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*- 2 * 3 * ***** BEGIN LICENSE BLOCK ***** 4 * Version: MPL 1.1/GPL 2.0/LGPL 2.1 5 * 6 * The contents of this file are subject to the Mozilla Public License Version 7 * 1.1 (the "License"); you may not use this file except in compliance with 8 * the License. You may obtain a copy of the License at 9 * http://www.mozilla.org/MPL/ 10 * 11 * Software distributed under the License is distributed on an "AS IS" basis, 12 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 13 * for the specific language governing rights and limitations under the 14 * License. 15 * 16 * The Original Code is Mozilla Communicator client code, released 17 * March 31, 1998. 18 * 19 * The Initial Developer of the Original Code is 20 * Sun Microsystems, Inc. 21 * Portions created by the Initial Developer are Copyright (C) 1998 22 * the Initial Developer. All Rights Reserved. 23 * 24 * Contributor(s): 25 * 26 * Alternatively, the contents of this file may be used under the terms of 27 * either of the GNU General Public License Version 2 or later (the "GPL"), 28 * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 29 * in which case the provisions of the GPL or the LGPL are applicable instead 30 * of those above. If you wish to allow use of your version of this file only 31 * under the terms of either the GPL or the LGPL, and not to allow others to 32 * use your version of this file under the terms of the MPL, indicate your 33 * decision by deleting the provisions above and replace them with the notice 34 * and other provisions required by the GPL or the LGPL. If you do not delete 35 * the provisions above, a recipient may use your version of this file under 36 * the terms of any one of the MPL, the GPL or the LGPL. 37 * 38 * ***** END LICENSE BLOCK ***** */ 39 40 /* @(#)e_fmod.c 1.3 95/01/18 */ 41 /* 42 * ==================================================== 43 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 44 * 45 * Developed at SunSoft, a Sun Microsystems, Inc. business. 46 * Permission to use, copy, modify, and distribute this 47 * software is freely granted, provided that this notice 48 * is preserved. 49 * ==================================================== 50 */ 51 52 /* 53 * __ieee754_fmod(x,y) 54 * Return x mod y in exact arithmetic 55 * Method: shift and subtract 56 */ 57 58 #include "fdlibm.h" 59 60 #ifdef __STDC__ 61 static const double one = 1.0, Zero[] = {0.0, -0.0,}; 62 #else 63 static double one = 1.0, Zero[] = {0.0, -0.0,}; 64 #endif 65 66 #ifdef __STDC__ __ieee754_fmod(double x,double y)67 double __ieee754_fmod(double x, double y) 68 #else 69 double __ieee754_fmod(x,y) 70 double x,y ; 71 #endif 72 { 73 fd_twoints ux, uy; 74 int n,hx,hy,hz,ix,iy,sx,i; 75 unsigned lx,ly,lz; 76 77 ux.d = x; uy.d = y; 78 hx = __HI(ux); /* high word of x */ 79 lx = __LO(ux); /* low word of x */ 80 hy = __HI(uy); /* high word of y */ 81 ly = __LO(uy); /* low word of y */ 82 sx = hx&0x80000000; /* sign of x */ 83 hx ^=sx; /* |x| */ 84 hy &= 0x7fffffff; /* |y| */ 85 86 /* purge off exception values */ 87 if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ 88 ((hy|((ly|-(int)ly)>>31))>0x7ff00000)) /* or y is NaN */ 89 return (x*y)/(x*y); 90 if(hx<=hy) { 91 if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */ 92 if(lx==ly) 93 return Zero[(unsigned)sx>>31]; /* |x|=|y| return x*0*/ 94 } 95 96 /* determine ix = ilogb(x) */ 97 if(hx<0x00100000) { /* subnormal x */ 98 if(hx==0) { 99 for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; 100 } else { 101 for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; 102 } 103 } else ix = (hx>>20)-1023; 104 105 /* determine iy = ilogb(y) */ 106 if(hy<0x00100000) { /* subnormal y */ 107 if(hy==0) { 108 for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; 109 } else { 110 for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; 111 } 112 } else iy = (hy>>20)-1023; 113 114 /* set up {hx,lx}, {hy,ly} and align y to x */ 115 if(ix >= -1022) 116 hx = 0x00100000|(0x000fffff&hx); 117 else { /* subnormal x, shift x to normal */ 118 n = -1022-ix; 119 if(n<=31) { 120 hx = (hx<<n)|(lx>>(32-n)); 121 lx <<= n; 122 } else { 123 hx = lx<<(n-32); 124 lx = 0; 125 } 126 } 127 if(iy >= -1022) 128 hy = 0x00100000|(0x000fffff&hy); 129 else { /* subnormal y, shift y to normal */ 130 n = -1022-iy; 131 if(n<=31) { 132 hy = (hy<<n)|(ly>>(32-n)); 133 ly <<= n; 134 } else { 135 hy = ly<<(n-32); 136 ly = 0; 137 } 138 } 139 140 /* fix point fmod */ 141 n = ix - iy; 142 while(n--) { 143 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 144 if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} 145 else { 146 if((hz|lz)==0) /* return sign(x)*0 */ 147 return Zero[(unsigned)sx>>31]; 148 hx = hz+hz+(lz>>31); lx = lz+lz; 149 } 150 } 151 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 152 if(hz>=0) {hx=hz;lx=lz;} 153 154 /* convert back to floating value and restore the sign */ 155 if((hx|lx)==0) /* return sign(x)*0 */ 156 return Zero[(unsigned)sx>>31]; 157 while(hx<0x00100000) { /* normalize x */ 158 hx = hx+hx+(lx>>31); lx = lx+lx; 159 iy -= 1; 160 } 161 if(iy>= -1022) { /* normalize output */ 162 hx = ((hx-0x00100000)|((iy+1023)<<20)); 163 ux.d = x; 164 __HI(ux) = hx|sx; 165 __LO(ux) = lx; 166 x = ux.d; 167 } else { /* subnormal output */ 168 n = -1022 - iy; 169 if(n<=20) { 170 lx = (lx>>n)|((unsigned)hx<<(32-n)); 171 hx >>= n; 172 } else if (n<=31) { 173 lx = (hx<<(32-n))|(lx>>n); hx = sx; 174 } else { 175 lx = hx>>(n-32); hx = sx; 176 } 177 ux.d = x; 178 __HI(ux) = hx|sx; 179 __LO(ux) = lx; 180 x = ux.d; 181 x *= one; /* create necessary signal */ 182 } 183 return x; /* exact output */ 184 } 185