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 fmodl = __fmodl 31 32 #include "libm.h" 33 34 static const int 35 is = -0x7fffffff - 1, 36 im = 0x0000ffff, 37 iu = 0x00010000; 38 39 static const long double 40 zero = 0.0L, 41 one = 1.0L; 42 43 #ifdef __LITTLE_ENDIAN 44 #define __H0(x) *(3 + (int *) &x) 45 #define __H1(x) *(2 + (int *) &x) 46 #define __H2(x) *(1 + (int *) &x) 47 #define __H3(x) *(0 + (int *) &x) 48 #else 49 #define __H0(x) *(0 + (int *) &x) 50 #define __H1(x) *(1 + (int *) &x) 51 #define __H2(x) *(2 + (int *) &x) 52 #define __H3(x) *(3 + (int *) &x) 53 #endif 54 55 long double 56 fmodl(long double x, long double y) { 57 long double a, b; 58 int n, ix, iy, k, sx; 59 int hx; 60 int x0, y0, z0, carry; 61 unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3; 62 63 hx = __H0(x); 64 x1 = __H1(x); 65 x2 = __H2(x); 66 x3 = __H3(x); 67 y0 = __H0(y); 68 y1 = __H1(y); 69 y2 = __H2(y); 70 y3 = __H3(y); 71 72 sx = hx & 0x80000000; 73 x0 = hx ^ sx; 74 y0 &= 0x7fffffff; 75 76 /* purge off exception values */ 77 if (x0 >= 0x7fff0000 || /* !finitel(x) */ 78 (y0 > 0x7fff0000) || (y0 == 0x7fff0000 && ((y1 | y2 | y3) != 0)) || 79 (y0 | y1 | y2 | y3) == 0) /* isnanl(y) || y = 0 */ 80 return ((x * y) / (x * y)); 81 a = fabsl(x); 82 b = fabsl(y); 83 if (a <= b) { 84 if (a < b) 85 return (x); 86 else 87 return (zero * x); 88 } 89 /* determine ix = ilogbl(x) */ 90 if (x0 < iu) { /* subnormal x */ 91 ix = -16382; 92 while (x0 == 0) { 93 ix -= 16; 94 x0 = x1 >> 16; 95 x1 = (x1 << 16) | (x2 >> 16); 96 x2 = (x2 << 16) | (x3 >> 16); 97 x3 = (x3 << 16); 98 } 99 while (x0 < iu) { 100 ix -= 1; 101 x0 = (x0 << 1) | (x1 >> 31); 102 x1 = (x1 << 1) | (x2 >> 31); 103 x2 = (x2 << 1) | (x3 >> 31); 104 x3 <<= 1; 105 } 106 } else { 107 ix = (x0 >> 16) - 16383; 108 x0 = iu | (x0 & im); 109 } 110 111 /* determine iy = ilogbl(y) */ 112 if (y0 < iu) { /* subnormal y */ 113 iy = -16382; 114 while (y0 == 0) { 115 iy -= 16; 116 y0 = y1 >> 16; 117 y1 = (y1 << 16) | (y2 >> 16); 118 y2 = (y2 << 16) | (y3 >> 16); 119 y3 = (y3 << 16); 120 } 121 while (y0 < iu) { 122 iy -= 1; 123 y0 = (y0 << 1) | (y1 >> 31); 124 y1 = (y1 << 1) | (y2 >> 31); 125 y2 = (y2 << 1) | (y3 >> 31); 126 y3 <<= 1; 127 } 128 } else { 129 iy = (y0 >> 16) - 16383; 130 y0 = iu | (y0 & im); 131 } 132 133 /* fix point fmod */ 134 n = ix - iy; 135 while (n--) { 136 while (x0 == 0 && n >= 16) { 137 n -= 16; 138 x0 = x1 >> 16; 139 x1 = (x1 << 16) | (x2 >> 16); 140 x2 = (x2 << 16) | (x3 >> 16); 141 x3 = (x3 << 16); 142 } 143 while (x0 < iu && n >= 1) { 144 n -= 1; 145 x0 = (x0 << 1) | (x1 >> 31); 146 x1 = (x1 << 1) | (x2 >> 31); 147 x2 = (x2 << 1) | (x3 >> 31); 148 x3 = (x3 << 1); 149 } 150 carry = 0; 151 z3 = x3 - y3; 152 carry = (z3 > x3); 153 if (carry == 0) { 154 z2 = x2 - y2; 155 carry = (z2 > x2); 156 } else { 157 z2 = x2 - y2 - 1; 158 carry = (z2 >= x2); 159 } 160 if (carry == 0) { 161 z1 = x1 - y1; 162 carry = (z1 > x1); 163 } else { 164 z1 = x1 - y1 - 1; 165 carry = (z1 >= x1); 166 } 167 z0 = x0 - y0 - carry; 168 if (z0 < 0) { /* double x */ 169 x0 = x0 + x0 + ((x1 & is) != 0); 170 x1 = x1 + x1 + ((x2 & is) != 0); 171 x2 = x2 + x2 + ((x3 & is) != 0); 172 x3 = x3 + x3; 173 } else { 174 if (z0 == 0) { 175 if ((z1 | z2 | z3) == 0) { /* 0: done */ 176 __H0(a) = hx & is; 177 __H1(a) = __H2(a) = __H3(a) = 0; 178 return (a); 179 } 180 } 181 /* x = z << 1 */ 182 z0 = z0 + z0 + ((z1 & is) != 0); 183 z1 = z1 + z1 + ((z2 & is) != 0); 184 z2 = z2 + z2 + ((z3 & is) != 0); 185 z3 = z3 + z3; 186 x0 = z0; 187 x1 = z1; 188 x2 = z2; 189 x3 = z3; 190 } 191 } 192 193 carry = 0; 194 z3 = x3 - y3; 195 carry = (z3 > x3); 196 if (carry == 0) { 197 z2 = x2 - y2; 198 carry = (z2 > x2); 199 } else { 200 z2 = x2 - y2 - 1; 201 carry = (z2 >= x2); 202 } 203 if (carry == 0) { 204 z1 = x1 - y1; 205 carry = (z1 > x1); 206 } else { 207 z1 = x1 - y1 - 1; 208 carry = (z1 >= x1); 209 } 210 z0 = x0 - y0 - carry; 211 if (z0 >= 0) { 212 x0 = z0; 213 x1 = z1; 214 x2 = z2; 215 x3 = z3; 216 } 217 /* convert back to floating value and restore the sign */ 218 if ((x0 | x1 | x2 | x3) == 0) { 219 __H0(a) = hx & is; 220 __H1(a) = __H2(a) = __H3(a) = 0; 221 return (a); 222 } 223 while (x0 < iu) { 224 if (x0 == 0) { 225 iy -= 16; 226 x0 = x1 >> 16; 227 x1 = (x1 << 16) | (x2 >> 16); 228 x2 = (x2 << 16) | (x3 >> 16); 229 x3 = (x3 << 16); 230 } else { 231 x0 = x0 + x0 + ((x1 & is) != 0); 232 x1 = x1 + x1 + ((x2 & is) != 0); 233 x2 = x2 + x2 + ((x3 & is) != 0); 234 x3 = x3 + x3; 235 iy -= 1; 236 } 237 } 238 239 /* normalize output */ 240 if (iy >= -16382) { 241 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16); 242 __H1(a) = x1; 243 __H2(a) = x2; 244 __H3(a) = x3; 245 } else { /* subnormal output */ 246 n = -16382 - iy; 247 k = n & 31; 248 if (k != 0) { 249 if (k <= 16) { 250 x3 = (x2 << (32 - k)) | (x3 >> k); 251 x2 = (x1 << (32 - k)) | (x2 >> k); 252 x1 = (x0 << (32 - k)) | (x1 >> k); 253 x0 >>= k; 254 } else { 255 x3 = (x2 << (32 - k)) | (x3 >> k); 256 x2 = (x1 << (32 - k)) | (x2 >> k); 257 x1 = (x0 << (32 - k)) | (x1 >> k); 258 x0 = 0; 259 } 260 } 261 while (n >= 32) { 262 n -= 32; 263 x3 = x2; 264 x2 = x1; 265 x1 = x0; 266 x0 = 0; 267 } 268 __H0(a) = x0 | sx; 269 __H1(a) = x1; 270 __H2(a) = x2; 271 __H3(a) = x3; 272 a *= one; 273 } 274 return (a); 275 } 276