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 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #pragma weak __fmodf = fmodf 30 31 #include "libm.h" 32 33 /* INDENT OFF */ 34 static const int 35 is = (int)0x80000000, 36 im = 0x007fffff, 37 ii = 0x7f800000, 38 iu = 0x00800000; 39 /* INDENT ON */ 40 41 static const float zero = 0.0; 42 43 float 44 fmodf(float x, float y) { 45 float w; 46 int hx, ix, iy, iz, k, ny, nd; 47 48 hx = *(int *)&x; 49 ix = hx & 0x7fffffff; 50 iy = *(int *)&y & 0x7fffffff; 51 52 /* purge off exception values */ 53 if (ix >= ii || iy > ii || iy == 0) { 54 w = x * y; 55 w = w / w; 56 } else if (ix <= iy) { 57 if (ix < iy) 58 w = x; /* return x if |x|<|y| */ 59 else 60 w = zero * x; /* return sign(x)*0.0 */ 61 } else { 62 /* INDENT OFF */ 63 /* 64 * scale x,y to "normal" with 65 * ny = exponent of y 66 * nd = exponent of x minus exponent of y 67 */ 68 /* INDENT ON */ 69 ny = iy >> 23; 70 k = ix >> 23; 71 72 /* special case for subnormal y or x */ 73 if (ny == 0) { 74 ny = 1; 75 while (iy < iu) { 76 ny -= 1; 77 iy += iy; 78 } 79 nd = k - ny; 80 if (k == 0) { 81 nd += 1; 82 while (ix < iu) { 83 nd -= 1; 84 ix += ix; 85 } 86 } else { 87 ix = iu | (ix & im); 88 } 89 } else { 90 nd = k - ny; 91 ix = iu | (ix & im); 92 iy = iu | (iy & im); 93 } 94 95 /* fix point fmod for normalized ix and iy */ 96 /* INDENT OFF */ 97 /* 98 * while (nd--) { 99 * iz = ix - iy; 100 * if (iz < 0) 101 * ix = ix + ix; 102 * else if (iz == 0) { 103 * *(int *) &w = is & hx; 104 * return w; 105 * } 106 * else 107 * ix = iz + iz; 108 * } 109 */ 110 /* INDENT ON */ 111 /* unroll the above loop 4 times to gain performance */ 112 k = nd >> 2; 113 nd -= k << 2; 114 while (k--) { 115 iz = ix - iy; 116 if (iz >= 0) 117 ix = iz + iz; 118 else 119 ix += ix; 120 iz = ix - iy; 121 if (iz >= 0) 122 ix = iz + iz; 123 else 124 ix += ix; 125 iz = ix - iy; 126 if (iz >= 0) 127 ix = iz + iz; 128 else 129 ix += ix; 130 iz = ix - iy; 131 if (iz >= 0) 132 ix = iz + iz; 133 else 134 ix += ix; 135 if (iz == 0) { 136 *(int *)&w = is & hx; 137 return (w); 138 } 139 } 140 while (nd--) { 141 iz = ix - iy; 142 if (iz >= 0) 143 ix = iz + iz; 144 else 145 ix += ix; 146 } 147 /* end of unrolling */ 148 149 iz = ix - iy; 150 if (iz >= 0) 151 ix = iz; 152 153 /* convert back to floating value and restore the sign */ 154 if (ix == 0) { 155 *(int *)&w = is & hx; 156 return (w); 157 } 158 while (ix < iu) { 159 ix += ix; 160 ny -= 1; 161 } 162 while (ix > (iu + iu)) { 163 ny += 1; 164 ix >>= 1; 165 } 166 if (ny > 0) { 167 *(int *)&w = (is & hx) | (ix & im) | (ny << 23); 168 } else { 169 /* subnormal output */ 170 k = -ny + 1; 171 ix >>= k; 172 *(int *)&w = (is & hx) | ix; 173 } 174 } 175 return (w); 176 } 177