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