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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #pragma weak __fmod = fmod 30 31 #include "libm.h" 32 33 static const double zero = 0.0; 34 35 /* 36 * The following implementation assumes fast 64-bit integer arith- 37 * metic. This is fine for sparc because we build libm in v8plus 38 * mode. It's also fine for sparcv9 and amd64, although we have 39 * assembly code on amd64. For x86, it would be better to use 40 * 32-bit code, but we have assembly for x86, too. 41 */ 42 double 43 fmod(double x, double y) { 44 double w; 45 long long hx, ix, iy, iz; 46 int nd, k, ny; 47 48 hx = *(long long *)&x; 49 ix = hx & ~0x8000000000000000ull; 50 iy = *(long long *)&y & ~0x8000000000000000ull; 51 52 /* handle special cases */ 53 if (iy == 0ll) 54 return (_SVID_libm_err(x, y, 27)); 55 56 if (ix >= 0x7ff0000000000000ll || iy > 0x7ff0000000000000ll) 57 return ((x * y) * zero); 58 59 if (ix <= iy) 60 return ((ix < iy)? x : x * zero); 61 62 /* 63 * Set: 64 * ny = true exponent of y 65 * nd = true exponent of x minus true exponent of y 66 * ix = normalized significand of x 67 * iy = normalized significand of y 68 */ 69 ny = iy >> 52; 70 k = ix >> 52; 71 if (ny == 0) { 72 /* y is subnormal, x could be normal or subnormal */ 73 ny = 1; 74 while (iy < 0x0010000000000000ll) { 75 ny -= 1; 76 iy += iy; 77 } 78 nd = k - ny; 79 if (k == 0) { 80 nd += 1; 81 while (ix < 0x0010000000000000ll) { 82 nd -= 1; 83 ix += ix; 84 } 85 } else { 86 ix = 0x0010000000000000ll | (ix & 0x000fffffffffffffll); 87 } 88 } else { 89 /* both x and y are normal */ 90 nd = k - ny; 91 ix = 0x0010000000000000ll | (ix & 0x000fffffffffffffll); 92 iy = 0x0010000000000000ll | (iy & 0x000fffffffffffffll); 93 } 94 95 /* perform fixed point mod */ 96 while (nd--) { 97 iz = ix - iy; 98 if (iz >= 0) 99 ix = iz; 100 ix += ix; 101 } 102 iz = ix - iy; 103 if (iz >= 0) 104 ix = iz; 105 106 /* convert back to floating point and restore the sign */ 107 if (ix == 0ll) 108 return (x * zero); 109 while (ix < 0x0010000000000000ll) { 110 ix += ix; 111 ny -= 1; 112 } 113 while (ix > 0x0020000000000000ll) { /* XXX can this ever happen? */ 114 ny += 1; 115 ix >>= 1; 116 } 117 if (ny <= 0) { 118 /* result is subnormal */ 119 k = -ny + 1; 120 ix >>= k; 121 *(long long *)&w = (hx & 0x8000000000000000ull) | ix; 122 return (w); 123 } 124 *(long long *)&w = (hx & 0x8000000000000000ull) | 125 ((long long)ny << 52) | (ix & 0x000fffffffffffffll); 126 return (w); 127 } 128