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_synonyms.h" 43 #include "libm_protos.h" 44 #include <math.h> 45 extern float fabsf(float); 46 47 static const int 48 is = (int) 0x80000000, 49 im = 0x007fffff, 50 ii = 0x7f800000, 51 iu = 0x00800000; 52 53 static const float zero = 0.0F, half = 0.5F; 54 /* INDENT ON */ 55 56 static float 57 fmodquof(float x, float y, int *quo) { 58 float w; 59 int hx, ix, iy, iz, k, ny, nd, m, sq; 60 61 hx = *(int *) &x; 62 ix = hx & 0x7fffffff; 63 iy = *(int *) &y; 64 sq = (iy ^ hx) & is; /* sign of x/y */ 65 iy &= 0x7fffffff; 66 67 /* purge off exception values */ 68 *quo = 0; 69 if (ix >= ii || iy > ii || iy == 0) { 70 w = x * y; 71 w = w / w; 72 } else if (ix <= iy) { 73 if (ix < iy) 74 w = x; /* return x if |x|<|y| */ 75 else { 76 *quo = 1 + (sq >> 30); 77 w = zero * x; /* return sign(x)*0.0 */ 78 } 79 } else { 80 /* INDENT OFF */ 81 /* 82 * scale x,y to "normal" with 83 * ny = exponent of y 84 * nd = exponent of x minus exponent of y 85 */ 86 /* INDENT ON */ 87 ny = iy >> 23; 88 k = ix >> 23; 89 90 /* special case for subnormal y or x */ 91 if (ny == 0) { 92 ny = 1; 93 while (iy < iu) { 94 ny -= 1; 95 iy += iy; 96 } 97 nd = k - ny; 98 if (k == 0) { 99 nd += 1; 100 while (ix < iu) { 101 nd -= 1; 102 ix += ix; 103 } 104 } else 105 ix = iu | (ix & im); 106 } else { 107 nd = k - ny; 108 ix = iu | (ix & im); 109 iy = iu | (iy & im); 110 } 111 /* INDENT OFF */ 112 /* fix point fmod for normalized ix and iy */ 113 /* 114 * while (nd--) { 115 * iz = ix - iy; 116 * if (iz < 0) 117 * ix = ix + ix; 118 * else if (iz == 0) { 119 * *(int *) &w = is & hx; 120 * return w; 121 * } else 122 * ix = iz + iz; 123 * } 124 */ 125 /* INDENT ON */ 126 /* unroll the above loop 4 times to gain performance */ 127 m = 0; 128 k = nd >> 2; 129 nd -= (k << 2); 130 while (k--) { 131 iz = ix - iy; 132 if (iz >= 0) { 133 m += 1; 134 ix = iz + iz; 135 } else 136 ix += ix; 137 m += m; 138 iz = ix - iy; 139 if (iz >= 0) { 140 m += 1; 141 ix = iz + iz; 142 } else 143 ix += ix; 144 m += m; 145 iz = ix - iy; 146 if (iz >= 0) { 147 m += 1; 148 ix = iz + iz; 149 } else 150 ix += ix; 151 m += m; 152 iz = ix - iy; 153 if (iz >= 0) { 154 m += 1; 155 ix = iz + iz; 156 } else 157 ix += ix; 158 m += m; 159 if (iz == 0) { 160 iz = (k << 2) + nd; 161 if (iz < 32) 162 m <<= iz; 163 else 164 m = 0; 165 m &= 0x7fffffff; 166 *quo = sq >= 0 ? m : -m; 167 *(int *) &w = is & hx; 168 return (w); 169 } 170 } 171 while (nd--) { 172 iz = ix - iy; 173 if (iz >= 0) { 174 m += 1; 175 ix = iz + iz; 176 } else 177 ix += ix; 178 m += m; 179 } 180 /* end of unrolling */ 181 182 iz = ix - iy; 183 if (iz >= 0) { 184 m += 1; 185 ix = iz; 186 } 187 m &= 0x7fffffff; 188 *quo = sq >= 0 ? m : -m; 189 190 /* convert back to floating value and restore the sign */ 191 if (ix == 0) { 192 *(int *) &w = is & hx; 193 return (w); 194 } 195 while (ix < iu) { 196 ix += ix; 197 ny -= 1; 198 } 199 while (ix > (iu + iu)) { 200 ny += 1; 201 ix >>= 1; 202 } 203 if (ny > 0) 204 *(int *) &w = (is & hx) | (ix & im) | (ny << 23); 205 else { /* subnormal output */ 206 k = -ny + 1; 207 ix >>= k; 208 *(int *) &w = (is & hx) | ix; 209 } 210 } 211 return (w); 212 } 213 214 float 215 remquof(float x, float y, int *quo) { 216 int hx, hy, sx, sq; 217 float v; 218 219 hx = *(int *) &x; /* high word of x */ 220 hy = *(int *) &y; /* high word of y */ 221 sx = hx & is; /* sign of x */ 222 sq = (hx ^ hy) & is; /* sign of x/y */ 223 hx ^= sx; /* |x| */ 224 hy &= 0x7fffffff; /* |y| */ 225 226 /* purge off exception values: y is 0 or NaN, x is Inf or NaN */ 227 *quo = 0; 228 if (hx >= ii || hy > ii || hy == 0) { 229 v = x * y; 230 return (v / v); 231 } 232 233 y = fabsf(y); 234 x = fabsf(x); 235 if (hy <= 0x7f7fffff) { 236 x = fmodquof(x, y + y, quo); 237 *quo = ((*quo) & 0x3fffffff) << 1; 238 } 239 if (hy < 0x01000000) { 240 if (x + x > y) { 241 *quo += 1; 242 if (x == y) 243 x = zero; 244 else 245 x -= y; 246 if (x + x >= y) { 247 x -= y; 248 *quo += 1; 249 } 250 } 251 } else { 252 v = half * y; 253 if (x > v) { 254 *quo += 1; 255 if (x == y) 256 x = zero; 257 else 258 x -= y; 259 if (x >= v) { 260 x -= y; 261 *quo += 1; 262 } 263 } 264 } 265 if (sq != 0) 266 *quo = -(*quo); 267 return (sx == 0 ? x : -x); 268 } 269