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 remquo = __remquo 31*25c28e83SPiotr Jasiukajtis 32*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 33*25c28e83SPiotr Jasiukajtis /* 34*25c28e83SPiotr Jasiukajtis * double remquo(double x, double y, int *quo) return remainder(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 integral part of x/y rounded to nearest even. 37*25c28e83SPiotr Jasiukajtis * 38*25c28e83SPiotr Jasiukajtis * remquo call internal fmodquo 39*25c28e83SPiotr Jasiukajtis */ 40*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 41*25c28e83SPiotr Jasiukajtis 42*25c28e83SPiotr Jasiukajtis #include "libm.h" 43*25c28e83SPiotr Jasiukajtis #include "libm_synonyms.h" 44*25c28e83SPiotr Jasiukajtis #include "libm_protos.h" 45*25c28e83SPiotr Jasiukajtis #include <math.h> /* fabs() */ 46*25c28e83SPiotr Jasiukajtis #include <sys/isa_defs.h> 47*25c28e83SPiotr Jasiukajtis 48*25c28e83SPiotr Jasiukajtis #if defined(_BIG_ENDIAN) 49*25c28e83SPiotr Jasiukajtis #define HIWORD 0 50*25c28e83SPiotr Jasiukajtis #define LOWORD 1 51*25c28e83SPiotr Jasiukajtis #else 52*25c28e83SPiotr Jasiukajtis #define HIWORD 1 53*25c28e83SPiotr Jasiukajtis #define LOWORD 0 54*25c28e83SPiotr Jasiukajtis #endif 55*25c28e83SPiotr Jasiukajtis #define __HI(x) ((int *) &x)[HIWORD] 56*25c28e83SPiotr Jasiukajtis #define __LO(x) ((int *) &x)[LOWORD] 57*25c28e83SPiotr Jasiukajtis 58*25c28e83SPiotr Jasiukajtis static const double one = 1.0, Zero[] = {0.0, -0.0}; 59*25c28e83SPiotr Jasiukajtis 60*25c28e83SPiotr Jasiukajtis static double 61*25c28e83SPiotr Jasiukajtis fmodquo(double x, double y, int *quo) { 62*25c28e83SPiotr Jasiukajtis int n, hx, hy, hz, ix, iy, sx, sq, i, m; 63*25c28e83SPiotr Jasiukajtis unsigned lx, ly, lz; 64*25c28e83SPiotr Jasiukajtis 65*25c28e83SPiotr Jasiukajtis hx = __HI(x); /* high word of x */ 66*25c28e83SPiotr Jasiukajtis lx = __LO(x); /* low word of x */ 67*25c28e83SPiotr Jasiukajtis hy = __HI(y); /* high word of y */ 68*25c28e83SPiotr Jasiukajtis ly = __LO(y); /* low word of y */ 69*25c28e83SPiotr Jasiukajtis sx = hx & 0x80000000; /* sign of x */ 70*25c28e83SPiotr Jasiukajtis sq = (hx ^ hy) & 0x80000000; /* sign of x/y */ 71*25c28e83SPiotr Jasiukajtis hx ^= sx; /* |x| */ 72*25c28e83SPiotr Jasiukajtis hy &= 0x7fffffff; /* |y| */ 73*25c28e83SPiotr Jasiukajtis 74*25c28e83SPiotr Jasiukajtis /* purge off exception values */ 75*25c28e83SPiotr Jasiukajtis *quo = 0; 76*25c28e83SPiotr Jasiukajtis if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */ 77*25c28e83SPiotr Jasiukajtis (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */ 78*25c28e83SPiotr Jasiukajtis return ((x * y) / (x * y)); 79*25c28e83SPiotr Jasiukajtis if (hx <= hy) { 80*25c28e83SPiotr Jasiukajtis if (hx < hy || lx < ly) 81*25c28e83SPiotr Jasiukajtis return (x); /* |x|<|y| return x */ 82*25c28e83SPiotr Jasiukajtis if (lx == ly) { 83*25c28e83SPiotr Jasiukajtis *quo = 1 + (sq >> 30); 84*25c28e83SPiotr Jasiukajtis /* |x|=|y| return x*0 */ 85*25c28e83SPiotr Jasiukajtis return (Zero[(unsigned) sx >> 31]); 86*25c28e83SPiotr Jasiukajtis } 87*25c28e83SPiotr Jasiukajtis } 88*25c28e83SPiotr Jasiukajtis 89*25c28e83SPiotr Jasiukajtis /* determine ix = ilogb(x) */ 90*25c28e83SPiotr Jasiukajtis if (hx < 0x00100000) { /* subnormal x */ 91*25c28e83SPiotr Jasiukajtis if (hx == 0) { 92*25c28e83SPiotr Jasiukajtis for (ix = -1043, i = lx; i > 0; i <<= 1) 93*25c28e83SPiotr Jasiukajtis ix -= 1; 94*25c28e83SPiotr Jasiukajtis } else { 95*25c28e83SPiotr Jasiukajtis for (ix = -1022, i = (hx << 11); i > 0; i <<= 1) 96*25c28e83SPiotr Jasiukajtis ix -= 1; 97*25c28e83SPiotr Jasiukajtis } 98*25c28e83SPiotr Jasiukajtis } else 99*25c28e83SPiotr Jasiukajtis ix = (hx >> 20) - 1023; 100*25c28e83SPiotr Jasiukajtis 101*25c28e83SPiotr Jasiukajtis /* determine iy = ilogb(y) */ 102*25c28e83SPiotr Jasiukajtis if (hy < 0x00100000) { /* subnormal y */ 103*25c28e83SPiotr Jasiukajtis if (hy == 0) { 104*25c28e83SPiotr Jasiukajtis for (iy = -1043, i = ly; i > 0; i <<= 1) 105*25c28e83SPiotr Jasiukajtis iy -= 1; 106*25c28e83SPiotr Jasiukajtis } else { 107*25c28e83SPiotr Jasiukajtis for (iy = -1022, i = (hy << 11); i > 0; i <<= 1) 108*25c28e83SPiotr Jasiukajtis iy -= 1; 109*25c28e83SPiotr Jasiukajtis } 110*25c28e83SPiotr Jasiukajtis } else 111*25c28e83SPiotr Jasiukajtis iy = (hy >> 20) - 1023; 112*25c28e83SPiotr Jasiukajtis 113*25c28e83SPiotr Jasiukajtis /* set up {hx,lx}, {hy,ly} and align y to x */ 114*25c28e83SPiotr Jasiukajtis if (ix >= -1022) 115*25c28e83SPiotr Jasiukajtis hx = 0x00100000 | (0x000fffff & hx); 116*25c28e83SPiotr Jasiukajtis else { /* subnormal x, shift x to normal */ 117*25c28e83SPiotr Jasiukajtis n = -1022 - ix; 118*25c28e83SPiotr Jasiukajtis if (n <= 31) { 119*25c28e83SPiotr Jasiukajtis hx = (hx << n) | (lx >> (32 - n)); 120*25c28e83SPiotr Jasiukajtis lx <<= n; 121*25c28e83SPiotr Jasiukajtis } else { 122*25c28e83SPiotr Jasiukajtis hx = lx << (n - 32); 123*25c28e83SPiotr Jasiukajtis lx = 0; 124*25c28e83SPiotr Jasiukajtis } 125*25c28e83SPiotr Jasiukajtis } 126*25c28e83SPiotr Jasiukajtis if (iy >= -1022) 127*25c28e83SPiotr Jasiukajtis hy = 0x00100000 | (0x000fffff & hy); 128*25c28e83SPiotr Jasiukajtis else { /* subnormal y, shift y to normal */ 129*25c28e83SPiotr Jasiukajtis n = -1022 - iy; 130*25c28e83SPiotr Jasiukajtis if (n <= 31) { 131*25c28e83SPiotr Jasiukajtis hy = (hy << n) | (ly >> (32 - n)); 132*25c28e83SPiotr Jasiukajtis ly <<= n; 133*25c28e83SPiotr Jasiukajtis } else { 134*25c28e83SPiotr Jasiukajtis hy = ly << (n - 32); 135*25c28e83SPiotr Jasiukajtis ly = 0; 136*25c28e83SPiotr Jasiukajtis } 137*25c28e83SPiotr Jasiukajtis } 138*25c28e83SPiotr Jasiukajtis 139*25c28e83SPiotr Jasiukajtis /* fix point fmod */ 140*25c28e83SPiotr Jasiukajtis n = ix - iy; 141*25c28e83SPiotr Jasiukajtis m = 0; 142*25c28e83SPiotr Jasiukajtis while (n--) { 143*25c28e83SPiotr Jasiukajtis hz = hx - hy; 144*25c28e83SPiotr Jasiukajtis lz = lx - ly; 145*25c28e83SPiotr Jasiukajtis if (lx < ly) 146*25c28e83SPiotr Jasiukajtis hz -= 1; 147*25c28e83SPiotr Jasiukajtis if (hz < 0) { 148*25c28e83SPiotr Jasiukajtis hx = hx + hx + (lx >> 31); 149*25c28e83SPiotr Jasiukajtis lx = lx + lx; 150*25c28e83SPiotr Jasiukajtis } else { 151*25c28e83SPiotr Jasiukajtis m += 1; 152*25c28e83SPiotr Jasiukajtis if ((hz | lz) == 0) { /* return sign(x)*0 */ 153*25c28e83SPiotr Jasiukajtis if (n < 31) 154*25c28e83SPiotr Jasiukajtis m <<= 1 + n; 155*25c28e83SPiotr Jasiukajtis else 156*25c28e83SPiotr Jasiukajtis m = 0; 157*25c28e83SPiotr Jasiukajtis m &= 0x7fffffff; 158*25c28e83SPiotr Jasiukajtis *quo = sq >= 0 ? m : -m; 159*25c28e83SPiotr Jasiukajtis return (Zero[(unsigned) sx >> 31]); 160*25c28e83SPiotr Jasiukajtis } 161*25c28e83SPiotr Jasiukajtis hx = hz + hz + (lz >> 31); 162*25c28e83SPiotr Jasiukajtis lx = lz + lz; 163*25c28e83SPiotr Jasiukajtis } 164*25c28e83SPiotr Jasiukajtis m += m; 165*25c28e83SPiotr Jasiukajtis } 166*25c28e83SPiotr Jasiukajtis hz = hx - hy; 167*25c28e83SPiotr Jasiukajtis lz = lx - ly; 168*25c28e83SPiotr Jasiukajtis if (lx < ly) 169*25c28e83SPiotr Jasiukajtis hz -= 1; 170*25c28e83SPiotr Jasiukajtis if (hz >= 0) { 171*25c28e83SPiotr Jasiukajtis hx = hz; 172*25c28e83SPiotr Jasiukajtis lx = lz; 173*25c28e83SPiotr Jasiukajtis m += 1; 174*25c28e83SPiotr Jasiukajtis } 175*25c28e83SPiotr Jasiukajtis m &= 0x7fffffff; 176*25c28e83SPiotr Jasiukajtis *quo = sq >= 0 ? m : -m; 177*25c28e83SPiotr Jasiukajtis 178*25c28e83SPiotr Jasiukajtis /* convert back to floating value and restore the sign */ 179*25c28e83SPiotr Jasiukajtis if ((hx | lx) == 0) { /* return sign(x)*0 */ 180*25c28e83SPiotr Jasiukajtis return (Zero[(unsigned) sx >> 31]); 181*25c28e83SPiotr Jasiukajtis } 182*25c28e83SPiotr Jasiukajtis while (hx < 0x00100000) { /* normalize x */ 183*25c28e83SPiotr Jasiukajtis hx = hx + hx + (lx >> 31); 184*25c28e83SPiotr Jasiukajtis lx = lx + lx; 185*25c28e83SPiotr Jasiukajtis iy -= 1; 186*25c28e83SPiotr Jasiukajtis } 187*25c28e83SPiotr Jasiukajtis if (iy >= -1022) { /* normalize output */ 188*25c28e83SPiotr Jasiukajtis hx = (hx - 0x00100000) | ((iy + 1023) << 20); 189*25c28e83SPiotr Jasiukajtis __HI(x) = hx | sx; 190*25c28e83SPiotr Jasiukajtis __LO(x) = lx; 191*25c28e83SPiotr Jasiukajtis } else { /* subnormal output */ 192*25c28e83SPiotr Jasiukajtis n = -1022 - iy; 193*25c28e83SPiotr Jasiukajtis if (n <= 20) { 194*25c28e83SPiotr Jasiukajtis lx = (lx >> n) | ((unsigned) hx << (32 - n)); 195*25c28e83SPiotr Jasiukajtis hx >>= n; 196*25c28e83SPiotr Jasiukajtis } else if (n <= 31) { 197*25c28e83SPiotr Jasiukajtis lx = (hx << (32 - n)) | (lx >> n); 198*25c28e83SPiotr Jasiukajtis hx = sx; 199*25c28e83SPiotr Jasiukajtis } else { 200*25c28e83SPiotr Jasiukajtis lx = hx >> (n - 32); 201*25c28e83SPiotr Jasiukajtis hx = sx; 202*25c28e83SPiotr Jasiukajtis } 203*25c28e83SPiotr Jasiukajtis __HI(x) = hx | sx; 204*25c28e83SPiotr Jasiukajtis __LO(x) = lx; 205*25c28e83SPiotr Jasiukajtis x *= one; /* create necessary signal */ 206*25c28e83SPiotr Jasiukajtis } 207*25c28e83SPiotr Jasiukajtis return (x); /* exact output */ 208*25c28e83SPiotr Jasiukajtis } 209*25c28e83SPiotr Jasiukajtis 210*25c28e83SPiotr Jasiukajtis #define zero Zero[0] 211*25c28e83SPiotr Jasiukajtis 212*25c28e83SPiotr Jasiukajtis double 213*25c28e83SPiotr Jasiukajtis remquo(double x, double y, int *quo) { 214*25c28e83SPiotr Jasiukajtis int hx, hy, sx, sq; 215*25c28e83SPiotr Jasiukajtis double v; 216*25c28e83SPiotr Jasiukajtis unsigned ly; 217*25c28e83SPiotr Jasiukajtis 218*25c28e83SPiotr Jasiukajtis hx = __HI(x); /* high word of x */ 219*25c28e83SPiotr Jasiukajtis hy = __HI(y); /* high word of y */ 220*25c28e83SPiotr Jasiukajtis ly = __LO(y); /* low word of y */ 221*25c28e83SPiotr Jasiukajtis sx = hx & 0x80000000; /* sign of x */ 222*25c28e83SPiotr Jasiukajtis sq = (hx ^ hy) & 0x80000000; /* 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 */ 227*25c28e83SPiotr Jasiukajtis *quo = 0; 228*25c28e83SPiotr Jasiukajtis if ((hy | ly) == 0 || hx >= 0x7ff00000 || /* y=0, or x !finite */ 229*25c28e83SPiotr Jasiukajtis (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */ 230*25c28e83SPiotr Jasiukajtis return ((x * y) / (x * y)); 231*25c28e83SPiotr Jasiukajtis 232*25c28e83SPiotr Jasiukajtis y = fabs(y); 233*25c28e83SPiotr Jasiukajtis x = fabs(x); 234*25c28e83SPiotr Jasiukajtis if (hy <= 0x7fdfffff) { 235*25c28e83SPiotr Jasiukajtis x = fmodquo(x, y + y, quo); 236*25c28e83SPiotr Jasiukajtis *quo = ((*quo) & 0x3fffffff) << 1; 237*25c28e83SPiotr Jasiukajtis } 238*25c28e83SPiotr Jasiukajtis if (hy < 0x00200000) { 239*25c28e83SPiotr Jasiukajtis if (x + x > y) { 240*25c28e83SPiotr Jasiukajtis *quo += 1; 241*25c28e83SPiotr Jasiukajtis if (x == y) 242*25c28e83SPiotr Jasiukajtis x = zero; 243*25c28e83SPiotr Jasiukajtis else 244*25c28e83SPiotr Jasiukajtis x -= y; 245*25c28e83SPiotr Jasiukajtis if (x + x >= y) { 246*25c28e83SPiotr Jasiukajtis x -= y; 247*25c28e83SPiotr Jasiukajtis *quo += 1; 248*25c28e83SPiotr Jasiukajtis } 249*25c28e83SPiotr Jasiukajtis } 250*25c28e83SPiotr Jasiukajtis } else { 251*25c28e83SPiotr Jasiukajtis v = 0.5 * y; 252*25c28e83SPiotr Jasiukajtis if (x > v) { 253*25c28e83SPiotr Jasiukajtis *quo += 1; 254*25c28e83SPiotr Jasiukajtis if (x == y) 255*25c28e83SPiotr Jasiukajtis x = zero; 256*25c28e83SPiotr Jasiukajtis else 257*25c28e83SPiotr Jasiukajtis x -= y; 258*25c28e83SPiotr Jasiukajtis if (x >= v) { 259*25c28e83SPiotr Jasiukajtis x -= y; 260*25c28e83SPiotr Jasiukajtis *quo += 1; 261*25c28e83SPiotr Jasiukajtis } 262*25c28e83SPiotr Jasiukajtis } 263*25c28e83SPiotr Jasiukajtis } 264*25c28e83SPiotr Jasiukajtis if (sq != 0) 265*25c28e83SPiotr Jasiukajtis *quo = -(*quo); 266*25c28e83SPiotr Jasiukajtis return (sx == 0 ? x : -x); 267*25c28e83SPiotr Jasiukajtis } 268