1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #pragma weak __remquof = remquof 31 32 /* INDENT OFF */ 33 /* 34 * float remquof(float x, float y, int *quo) return remainderf(x,y) and an 35 * integer pointer quo such that *quo = N mod (2**31), where N is the 36 * exact integeral part of x/y rounded to nearest even. 37 * 38 * remquof call internal fmodquof 39 */ 40 41 #include "libm.h" 42 #include "libm_protos.h" 43 #include <math.h> 44 extern float fabsf(float); 45 46 static const int 47 is = (int) 0x80000000, 48 im = 0x007fffff, 49 ii = 0x7f800000, 50 iu = 0x00800000; 51 52 static const float zero = 0.0F, half = 0.5F; 53 /* INDENT ON */ 54 55 static float 56 fmodquof(float x, float y, int *quo) { 57 float w; 58 int hx, ix, iy, iz, k, ny, nd, m, sq; 59 60 hx = *(int *) &x; 61 ix = hx & 0x7fffffff; 62 iy = *(int *) &y; 63 sq = (iy ^ hx) & is; /* sign of x/y */ 64 iy &= 0x7fffffff; 65 66 /* purge off exception values */ 67 *quo = 0; 68 if (ix >= ii || iy > ii || iy == 0) { 69 w = x * y; 70 w = w / w; 71 } else if (ix <= iy) { 72 if (ix < iy) 73 w = x; /* return x if |x|<|y| */ 74 else { 75 *quo = 1 + (sq >> 30); 76 w = zero * x; /* return sign(x)*0.0 */ 77 } 78 } else { 79 /* INDENT OFF */ 80 /* 81 * scale x,y to "normal" with 82 * ny = exponent of y 83 * nd = exponent of x minus exponent of y 84 */ 85 /* INDENT ON */ 86 ny = iy >> 23; 87 k = ix >> 23; 88 89 /* special case for subnormal y or x */ 90 if (ny == 0) { 91 ny = 1; 92 while (iy < iu) { 93 ny -= 1; 94 iy += iy; 95 } 96 nd = k - ny; 97 if (k == 0) { 98 nd += 1; 99 while (ix < iu) { 100 nd -= 1; 101 ix += ix; 102 } 103 } else 104 ix = iu | (ix & im); 105 } else { 106 nd = k - ny; 107 ix = iu | (ix & im); 108 iy = iu | (iy & im); 109 } 110 /* INDENT OFF */ 111 /* fix point fmod for normalized ix and iy */ 112 /* 113 * while (nd--) { 114 * iz = ix - iy; 115 * if (iz < 0) 116 * ix = ix + ix; 117 * else if (iz == 0) { 118 * *(int *) &w = is & hx; 119 * return w; 120 * } else 121 * ix = iz + iz; 122 * } 123 */ 124 /* INDENT ON */ 125 /* unroll the above loop 4 times to gain performance */ 126 m = 0; 127 k = nd >> 2; 128 nd -= (k << 2); 129 while (k--) { 130 iz = ix - iy; 131 if (iz >= 0) { 132 m += 1; 133 ix = iz + iz; 134 } else 135 ix += ix; 136 m += m; 137 iz = ix - iy; 138 if (iz >= 0) { 139 m += 1; 140 ix = iz + iz; 141 } else 142 ix += ix; 143 m += m; 144 iz = ix - iy; 145 if (iz >= 0) { 146 m += 1; 147 ix = iz + iz; 148 } else 149 ix += ix; 150 m += m; 151 iz = ix - iy; 152 if (iz >= 0) { 153 m += 1; 154 ix = iz + iz; 155 } else 156 ix += ix; 157 m += m; 158 if (iz == 0) { 159 iz = (k << 2) + nd; 160 if (iz < 32) 161 m <<= iz; 162 else 163 m = 0; 164 m &= 0x7fffffff; 165 *quo = sq >= 0 ? m : -m; 166 *(int *) &w = is & hx; 167 return (w); 168 } 169 } 170 while (nd--) { 171 iz = ix - iy; 172 if (iz >= 0) { 173 m += 1; 174 ix = iz + iz; 175 } else 176 ix += ix; 177 m += m; 178 } 179 /* end of unrolling */ 180 181 iz = ix - iy; 182 if (iz >= 0) { 183 m += 1; 184 ix = iz; 185 } 186 m &= 0x7fffffff; 187 *quo = sq >= 0 ? m : -m; 188 189 /* convert back to floating value and restore the sign */ 190 if (ix == 0) { 191 *(int *) &w = is & hx; 192 return (w); 193 } 194 while (ix < iu) { 195 ix += ix; 196 ny -= 1; 197 } 198 while (ix > (iu + iu)) { 199 ny += 1; 200 ix >>= 1; 201 } 202 if (ny > 0) 203 *(int *) &w = (is & hx) | (ix & im) | (ny << 23); 204 else { /* subnormal output */ 205 k = -ny + 1; 206 ix >>= k; 207 *(int *) &w = (is & hx) | ix; 208 } 209 } 210 return (w); 211 } 212 213 float 214 remquof(float x, float y, int *quo) { 215 int hx, hy, sx, sq; 216 float v; 217 218 hx = *(int *) &x; /* high word of x */ 219 hy = *(int *) &y; /* high word of y */ 220 sx = hx & is; /* sign of x */ 221 sq = (hx ^ hy) & is; /* sign of x/y */ 222 hx ^= sx; /* |x| */ 223 hy &= 0x7fffffff; /* |y| */ 224 225 /* purge off exception values: y is 0 or NaN, x is Inf or NaN */ 226 *quo = 0; 227 if (hx >= ii || hy > ii || hy == 0) { 228 v = x * y; 229 return (v / v); 230 } 231 232 y = fabsf(y); 233 x = fabsf(x); 234 if (hy <= 0x7f7fffff) { 235 x = fmodquof(x, y + y, quo); 236 *quo = ((*quo) & 0x3fffffff) << 1; 237 } 238 if (hy < 0x01000000) { 239 if (x + x > y) { 240 *quo += 1; 241 if (x == y) 242 x = zero; 243 else 244 x -= y; 245 if (x + x >= y) { 246 x -= y; 247 *quo += 1; 248 } 249 } 250 } else { 251 v = half * y; 252 if (x > v) { 253 *quo += 1; 254 if (x == y) 255 x = zero; 256 else 257 x -= y; 258 if (x >= v) { 259 x -= y; 260 *quo += 1; 261 } 262 } 263 } 264 if (sq != 0) 265 *quo = -(*quo); 266 return (sx == 0 ? x : -x); 267 } 268