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 #if defined(ELFOBJ) 31*25c28e83SPiotr Jasiukajtis #pragma weak fmal = __fmal 32*25c28e83SPiotr Jasiukajtis #endif 33*25c28e83SPiotr Jasiukajtis 34*25c28e83SPiotr Jasiukajtis #include "libm.h" 35*25c28e83SPiotr Jasiukajtis #include "fma.h" 36*25c28e83SPiotr Jasiukajtis #include "fenv_inlines.h" 37*25c28e83SPiotr Jasiukajtis 38*25c28e83SPiotr Jasiukajtis #if defined(__sparc) 39*25c28e83SPiotr Jasiukajtis 40*25c28e83SPiotr Jasiukajtis static const union { 41*25c28e83SPiotr Jasiukajtis unsigned i[2]; 42*25c28e83SPiotr Jasiukajtis double d; 43*25c28e83SPiotr Jasiukajtis } C[] = { 44*25c28e83SPiotr Jasiukajtis { 0x3fe00000u, 0 }, 45*25c28e83SPiotr Jasiukajtis { 0x40000000u, 0 }, 46*25c28e83SPiotr Jasiukajtis { 0x3ef00000u, 0 }, 47*25c28e83SPiotr Jasiukajtis { 0x3e700000u, 0 }, 48*25c28e83SPiotr Jasiukajtis { 0x41300000u, 0 }, 49*25c28e83SPiotr Jasiukajtis { 0x3e300000u, 0 }, 50*25c28e83SPiotr Jasiukajtis { 0x3b300000u, 0 }, 51*25c28e83SPiotr Jasiukajtis { 0x38300000u, 0 }, 52*25c28e83SPiotr Jasiukajtis { 0x42300000u, 0 }, 53*25c28e83SPiotr Jasiukajtis { 0x3df00000u, 0 }, 54*25c28e83SPiotr Jasiukajtis { 0x7fe00000u, 0 }, 55*25c28e83SPiotr Jasiukajtis { 0x00100000u, 0 }, 56*25c28e83SPiotr Jasiukajtis { 0x00100001u, 0 }, 57*25c28e83SPiotr Jasiukajtis { 0, 0 }, 58*25c28e83SPiotr Jasiukajtis { 0x7ff00000u, 0 }, 59*25c28e83SPiotr Jasiukajtis { 0x7ff00001u, 0 } 60*25c28e83SPiotr Jasiukajtis }; 61*25c28e83SPiotr Jasiukajtis 62*25c28e83SPiotr Jasiukajtis #define half C[0].d 63*25c28e83SPiotr Jasiukajtis #define two C[1].d 64*25c28e83SPiotr Jasiukajtis #define twom16 C[2].d 65*25c28e83SPiotr Jasiukajtis #define twom24 C[3].d 66*25c28e83SPiotr Jasiukajtis #define two20 C[4].d 67*25c28e83SPiotr Jasiukajtis #define twom28 C[5].d 68*25c28e83SPiotr Jasiukajtis #define twom76 C[6].d 69*25c28e83SPiotr Jasiukajtis #define twom124 C[7].d 70*25c28e83SPiotr Jasiukajtis #define two36 C[8].d 71*25c28e83SPiotr Jasiukajtis #define twom32 C[9].d 72*25c28e83SPiotr Jasiukajtis #define huge C[10].d 73*25c28e83SPiotr Jasiukajtis #define tiny C[11].d 74*25c28e83SPiotr Jasiukajtis #define tiny2 C[12].d 75*25c28e83SPiotr Jasiukajtis #define zero C[13].d 76*25c28e83SPiotr Jasiukajtis #define inf C[14].d 77*25c28e83SPiotr Jasiukajtis #define snan C[15].d 78*25c28e83SPiotr Jasiukajtis 79*25c28e83SPiotr Jasiukajtis static const unsigned int fsr_rm = 0xc0000000u; 80*25c28e83SPiotr Jasiukajtis 81*25c28e83SPiotr Jasiukajtis /* 82*25c28e83SPiotr Jasiukajtis * fmal for SPARC: 128-bit quad precision, big-endian 83*25c28e83SPiotr Jasiukajtis */ 84*25c28e83SPiotr Jasiukajtis long double 85*25c28e83SPiotr Jasiukajtis __fmal(long double x, long double y, long double z) { 86*25c28e83SPiotr Jasiukajtis union { 87*25c28e83SPiotr Jasiukajtis unsigned int i[4]; 88*25c28e83SPiotr Jasiukajtis long double q; 89*25c28e83SPiotr Jasiukajtis } xx, yy, zz; 90*25c28e83SPiotr Jasiukajtis union { 91*25c28e83SPiotr Jasiukajtis unsigned int i[2]; 92*25c28e83SPiotr Jasiukajtis double d; 93*25c28e83SPiotr Jasiukajtis } u; 94*25c28e83SPiotr Jasiukajtis double dx[5], dy[5], dxy[9], c, s; 95*25c28e83SPiotr Jasiukajtis unsigned int xy0, xy1, xy2, xy3, xy4, xy5, xy6, xy7; 96*25c28e83SPiotr Jasiukajtis unsigned int z0, z1, z2, z3, z4, z5, z6, z7; 97*25c28e83SPiotr Jasiukajtis unsigned int rm, sticky; 98*25c28e83SPiotr Jasiukajtis unsigned int fsr; 99*25c28e83SPiotr Jasiukajtis int hx, hy, hz, ex, ey, ez, exy, sxy, sz, e, ibit; 100*25c28e83SPiotr Jasiukajtis int cx, cy, cz; 101*25c28e83SPiotr Jasiukajtis volatile double dummy; 102*25c28e83SPiotr Jasiukajtis 103*25c28e83SPiotr Jasiukajtis /* extract the high order words of the arguments */ 104*25c28e83SPiotr Jasiukajtis xx.q = x; 105*25c28e83SPiotr Jasiukajtis yy.q = y; 106*25c28e83SPiotr Jasiukajtis zz.q = z; 107*25c28e83SPiotr Jasiukajtis hx = xx.i[0] & ~0x80000000; 108*25c28e83SPiotr Jasiukajtis hy = yy.i[0] & ~0x80000000; 109*25c28e83SPiotr Jasiukajtis hz = zz.i[0] & ~0x80000000; 110*25c28e83SPiotr Jasiukajtis 111*25c28e83SPiotr Jasiukajtis /* 112*25c28e83SPiotr Jasiukajtis * distinguish zero, finite nonzero, infinite, and quiet nan 113*25c28e83SPiotr Jasiukajtis * arguments; raise invalid and return for signaling nans 114*25c28e83SPiotr Jasiukajtis */ 115*25c28e83SPiotr Jasiukajtis if (hx >= 0x7fff0000) { 116*25c28e83SPiotr Jasiukajtis if ((hx & 0xffff) | xx.i[1] | xx.i[2] | xx.i[3]) { 117*25c28e83SPiotr Jasiukajtis if (!(hx & 0x8000)) { 118*25c28e83SPiotr Jasiukajtis /* signaling nan, raise invalid */ 119*25c28e83SPiotr Jasiukajtis dummy = snan; 120*25c28e83SPiotr Jasiukajtis dummy += snan; 121*25c28e83SPiotr Jasiukajtis xx.i[0] |= 0x8000; 122*25c28e83SPiotr Jasiukajtis return (xx.q); 123*25c28e83SPiotr Jasiukajtis } 124*25c28e83SPiotr Jasiukajtis cx = 3; /* quiet nan */ 125*25c28e83SPiotr Jasiukajtis } else 126*25c28e83SPiotr Jasiukajtis cx = 2; /* inf */ 127*25c28e83SPiotr Jasiukajtis } else if (hx == 0) { 128*25c28e83SPiotr Jasiukajtis cx = (xx.i[1] | xx.i[2] | xx.i[3]) ? 1 : 0; 129*25c28e83SPiotr Jasiukajtis /* subnormal or zero */ 130*25c28e83SPiotr Jasiukajtis } else 131*25c28e83SPiotr Jasiukajtis cx = 1; /* finite nonzero */ 132*25c28e83SPiotr Jasiukajtis 133*25c28e83SPiotr Jasiukajtis if (hy >= 0x7fff0000) { 134*25c28e83SPiotr Jasiukajtis if ((hy & 0xffff) | yy.i[1] | yy.i[2] | yy.i[3]) { 135*25c28e83SPiotr Jasiukajtis if (!(hy & 0x8000)) { 136*25c28e83SPiotr Jasiukajtis dummy = snan; 137*25c28e83SPiotr Jasiukajtis dummy += snan; 138*25c28e83SPiotr Jasiukajtis yy.i[0] |= 0x8000; 139*25c28e83SPiotr Jasiukajtis return (yy.q); 140*25c28e83SPiotr Jasiukajtis } 141*25c28e83SPiotr Jasiukajtis cy = 3; 142*25c28e83SPiotr Jasiukajtis } else 143*25c28e83SPiotr Jasiukajtis cy = 2; 144*25c28e83SPiotr Jasiukajtis } else if (hy == 0) { 145*25c28e83SPiotr Jasiukajtis cy = (yy.i[1] | yy.i[2] | yy.i[3]) ? 1 : 0; 146*25c28e83SPiotr Jasiukajtis } else 147*25c28e83SPiotr Jasiukajtis cy = 1; 148*25c28e83SPiotr Jasiukajtis 149*25c28e83SPiotr Jasiukajtis if (hz >= 0x7fff0000) { 150*25c28e83SPiotr Jasiukajtis if ((hz & 0xffff) | zz.i[1] | zz.i[2] | zz.i[3]) { 151*25c28e83SPiotr Jasiukajtis if (!(hz & 0x8000)) { 152*25c28e83SPiotr Jasiukajtis dummy = snan; 153*25c28e83SPiotr Jasiukajtis dummy += snan; 154*25c28e83SPiotr Jasiukajtis zz.i[0] |= 0x8000; 155*25c28e83SPiotr Jasiukajtis return (zz.q); 156*25c28e83SPiotr Jasiukajtis } 157*25c28e83SPiotr Jasiukajtis cz = 3; 158*25c28e83SPiotr Jasiukajtis } else 159*25c28e83SPiotr Jasiukajtis cz = 2; 160*25c28e83SPiotr Jasiukajtis } else if (hz == 0) { 161*25c28e83SPiotr Jasiukajtis cz = (zz.i[1] | zz.i[2] | zz.i[3]) ? 1 : 0; 162*25c28e83SPiotr Jasiukajtis } else 163*25c28e83SPiotr Jasiukajtis cz = 1; 164*25c28e83SPiotr Jasiukajtis 165*25c28e83SPiotr Jasiukajtis /* get the fsr and clear current exceptions */ 166*25c28e83SPiotr Jasiukajtis __fenv_getfsr32(&fsr); 167*25c28e83SPiotr Jasiukajtis fsr &= ~FSR_CEXC; 168*25c28e83SPiotr Jasiukajtis 169*25c28e83SPiotr Jasiukajtis /* handle all other zero, inf, and nan cases */ 170*25c28e83SPiotr Jasiukajtis if (cx != 1 || cy != 1 || cz != 1) { 171*25c28e83SPiotr Jasiukajtis /* if x or y is a quiet nan, return it */ 172*25c28e83SPiotr Jasiukajtis if (cx == 3) { 173*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 174*25c28e83SPiotr Jasiukajtis return (x); 175*25c28e83SPiotr Jasiukajtis } 176*25c28e83SPiotr Jasiukajtis if (cy == 3) { 177*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 178*25c28e83SPiotr Jasiukajtis return (y); 179*25c28e83SPiotr Jasiukajtis } 180*25c28e83SPiotr Jasiukajtis 181*25c28e83SPiotr Jasiukajtis /* if x*y is 0*inf, raise invalid and return the default nan */ 182*25c28e83SPiotr Jasiukajtis if ((cx == 0 && cy == 2) || (cx == 2 && cy == 0)) { 183*25c28e83SPiotr Jasiukajtis dummy = zero; 184*25c28e83SPiotr Jasiukajtis dummy *= inf; 185*25c28e83SPiotr Jasiukajtis zz.i[0] = 0x7fffffff; 186*25c28e83SPiotr Jasiukajtis zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff; 187*25c28e83SPiotr Jasiukajtis return (zz.q); 188*25c28e83SPiotr Jasiukajtis } 189*25c28e83SPiotr Jasiukajtis 190*25c28e83SPiotr Jasiukajtis /* if z is a quiet nan, return it */ 191*25c28e83SPiotr Jasiukajtis if (cz == 3) { 192*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 193*25c28e83SPiotr Jasiukajtis return (z); 194*25c28e83SPiotr Jasiukajtis } 195*25c28e83SPiotr Jasiukajtis 196*25c28e83SPiotr Jasiukajtis /* 197*25c28e83SPiotr Jasiukajtis * now none of x, y, or z is nan; handle cases where x or y 198*25c28e83SPiotr Jasiukajtis * is inf 199*25c28e83SPiotr Jasiukajtis */ 200*25c28e83SPiotr Jasiukajtis if (cx == 2 || cy == 2) { 201*25c28e83SPiotr Jasiukajtis /* 202*25c28e83SPiotr Jasiukajtis * if z is also inf, either we have inf-inf or 203*25c28e83SPiotr Jasiukajtis * the result is the same as z depending on signs 204*25c28e83SPiotr Jasiukajtis */ 205*25c28e83SPiotr Jasiukajtis if (cz == 2) { 206*25c28e83SPiotr Jasiukajtis if ((int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) { 207*25c28e83SPiotr Jasiukajtis dummy = inf; 208*25c28e83SPiotr Jasiukajtis dummy -= inf; 209*25c28e83SPiotr Jasiukajtis zz.i[0] = 0x7fffffff; 210*25c28e83SPiotr Jasiukajtis zz.i[1] = zz.i[2] = zz.i[3] = 211*25c28e83SPiotr Jasiukajtis 0xffffffff; 212*25c28e83SPiotr Jasiukajtis return (zz.q); 213*25c28e83SPiotr Jasiukajtis } 214*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 215*25c28e83SPiotr Jasiukajtis return (z); 216*25c28e83SPiotr Jasiukajtis } 217*25c28e83SPiotr Jasiukajtis 218*25c28e83SPiotr Jasiukajtis /* otherwise the result is inf with appropriate sign */ 219*25c28e83SPiotr Jasiukajtis zz.i[0] = ((xx.i[0] ^ yy.i[0]) & 0x80000000) | 220*25c28e83SPiotr Jasiukajtis 0x7fff0000; 221*25c28e83SPiotr Jasiukajtis zz.i[1] = zz.i[2] = zz.i[3] = 0; 222*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 223*25c28e83SPiotr Jasiukajtis return (zz.q); 224*25c28e83SPiotr Jasiukajtis } 225*25c28e83SPiotr Jasiukajtis 226*25c28e83SPiotr Jasiukajtis /* if z is inf, return it */ 227*25c28e83SPiotr Jasiukajtis if (cz == 2) { 228*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 229*25c28e83SPiotr Jasiukajtis return (z); 230*25c28e83SPiotr Jasiukajtis } 231*25c28e83SPiotr Jasiukajtis 232*25c28e83SPiotr Jasiukajtis /* 233*25c28e83SPiotr Jasiukajtis * now x, y, and z are all finite; handle cases where x or y 234*25c28e83SPiotr Jasiukajtis * is zero 235*25c28e83SPiotr Jasiukajtis */ 236*25c28e83SPiotr Jasiukajtis if (cx == 0 || cy == 0) { 237*25c28e83SPiotr Jasiukajtis /* either we have 0-0 or the result is the same as z */ 238*25c28e83SPiotr Jasiukajtis if (cz == 0 && (int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 239*25c28e83SPiotr Jasiukajtis 0) { 240*25c28e83SPiotr Jasiukajtis zz.i[0] = (fsr >> 30) == FSR_RM ? 0x80000000 : 241*25c28e83SPiotr Jasiukajtis 0; 242*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 243*25c28e83SPiotr Jasiukajtis return (zz.q); 244*25c28e83SPiotr Jasiukajtis } 245*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 246*25c28e83SPiotr Jasiukajtis return (z); 247*25c28e83SPiotr Jasiukajtis } 248*25c28e83SPiotr Jasiukajtis 249*25c28e83SPiotr Jasiukajtis /* if we get here, x and y are nonzero finite, z must be zero */ 250*25c28e83SPiotr Jasiukajtis return (x * y); 251*25c28e83SPiotr Jasiukajtis } 252*25c28e83SPiotr Jasiukajtis 253*25c28e83SPiotr Jasiukajtis /* 254*25c28e83SPiotr Jasiukajtis * now x, y, and z are all finite and nonzero; set round-to- 255*25c28e83SPiotr Jasiukajtis * negative-infinity mode 256*25c28e83SPiotr Jasiukajtis */ 257*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr_rm); 258*25c28e83SPiotr Jasiukajtis 259*25c28e83SPiotr Jasiukajtis /* 260*25c28e83SPiotr Jasiukajtis * get the signs and exponents and normalize the significands 261*25c28e83SPiotr Jasiukajtis * of x and y 262*25c28e83SPiotr Jasiukajtis */ 263*25c28e83SPiotr Jasiukajtis sxy = (xx.i[0] ^ yy.i[0]) & 0x80000000; 264*25c28e83SPiotr Jasiukajtis ex = hx >> 16; 265*25c28e83SPiotr Jasiukajtis hx &= 0xffff; 266*25c28e83SPiotr Jasiukajtis if (!ex) { 267*25c28e83SPiotr Jasiukajtis if (hx | (xx.i[1] & 0xfffe0000)) { 268*25c28e83SPiotr Jasiukajtis ex = 1; 269*25c28e83SPiotr Jasiukajtis } else if (xx.i[1] | (xx.i[2] & 0xfffe0000)) { 270*25c28e83SPiotr Jasiukajtis hx = xx.i[1]; 271*25c28e83SPiotr Jasiukajtis xx.i[1] = xx.i[2]; 272*25c28e83SPiotr Jasiukajtis xx.i[2] = xx.i[3]; 273*25c28e83SPiotr Jasiukajtis xx.i[3] = 0; 274*25c28e83SPiotr Jasiukajtis ex = -31; 275*25c28e83SPiotr Jasiukajtis } else if (xx.i[2] | (xx.i[3] & 0xfffe0000)) { 276*25c28e83SPiotr Jasiukajtis hx = xx.i[2]; 277*25c28e83SPiotr Jasiukajtis xx.i[1] = xx.i[3]; 278*25c28e83SPiotr Jasiukajtis xx.i[2] = xx.i[3] = 0; 279*25c28e83SPiotr Jasiukajtis ex = -63; 280*25c28e83SPiotr Jasiukajtis } else { 281*25c28e83SPiotr Jasiukajtis hx = xx.i[3]; 282*25c28e83SPiotr Jasiukajtis xx.i[1] = xx.i[2] = xx.i[3] = 0; 283*25c28e83SPiotr Jasiukajtis ex = -95; 284*25c28e83SPiotr Jasiukajtis } 285*25c28e83SPiotr Jasiukajtis while ((hx & 0x10000) == 0) { 286*25c28e83SPiotr Jasiukajtis hx = (hx << 1) | (xx.i[1] >> 31); 287*25c28e83SPiotr Jasiukajtis xx.i[1] = (xx.i[1] << 1) | (xx.i[2] >> 31); 288*25c28e83SPiotr Jasiukajtis xx.i[2] = (xx.i[2] << 1) | (xx.i[3] >> 31); 289*25c28e83SPiotr Jasiukajtis xx.i[3] <<= 1; 290*25c28e83SPiotr Jasiukajtis ex--; 291*25c28e83SPiotr Jasiukajtis } 292*25c28e83SPiotr Jasiukajtis } else 293*25c28e83SPiotr Jasiukajtis hx |= 0x10000; 294*25c28e83SPiotr Jasiukajtis ey = hy >> 16; 295*25c28e83SPiotr Jasiukajtis hy &= 0xffff; 296*25c28e83SPiotr Jasiukajtis if (!ey) { 297*25c28e83SPiotr Jasiukajtis if (hy | (yy.i[1] & 0xfffe0000)) { 298*25c28e83SPiotr Jasiukajtis ey = 1; 299*25c28e83SPiotr Jasiukajtis } else if (yy.i[1] | (yy.i[2] & 0xfffe0000)) { 300*25c28e83SPiotr Jasiukajtis hy = yy.i[1]; 301*25c28e83SPiotr Jasiukajtis yy.i[1] = yy.i[2]; 302*25c28e83SPiotr Jasiukajtis yy.i[2] = yy.i[3]; 303*25c28e83SPiotr Jasiukajtis yy.i[3] = 0; 304*25c28e83SPiotr Jasiukajtis ey = -31; 305*25c28e83SPiotr Jasiukajtis } else if (yy.i[2] | (yy.i[3] & 0xfffe0000)) { 306*25c28e83SPiotr Jasiukajtis hy = yy.i[2]; 307*25c28e83SPiotr Jasiukajtis yy.i[1] = yy.i[3]; 308*25c28e83SPiotr Jasiukajtis yy.i[2] = yy.i[3] = 0; 309*25c28e83SPiotr Jasiukajtis ey = -63; 310*25c28e83SPiotr Jasiukajtis } else { 311*25c28e83SPiotr Jasiukajtis hy = yy.i[3]; 312*25c28e83SPiotr Jasiukajtis yy.i[1] = yy.i[2] = yy.i[3] = 0; 313*25c28e83SPiotr Jasiukajtis ey = -95; 314*25c28e83SPiotr Jasiukajtis } 315*25c28e83SPiotr Jasiukajtis while ((hy & 0x10000) == 0) { 316*25c28e83SPiotr Jasiukajtis hy = (hy << 1) | (yy.i[1] >> 31); 317*25c28e83SPiotr Jasiukajtis yy.i[1] = (yy.i[1] << 1) | (yy.i[2] >> 31); 318*25c28e83SPiotr Jasiukajtis yy.i[2] = (yy.i[2] << 1) | (yy.i[3] >> 31); 319*25c28e83SPiotr Jasiukajtis yy.i[3] <<= 1; 320*25c28e83SPiotr Jasiukajtis ey--; 321*25c28e83SPiotr Jasiukajtis } 322*25c28e83SPiotr Jasiukajtis } else 323*25c28e83SPiotr Jasiukajtis hy |= 0x10000; 324*25c28e83SPiotr Jasiukajtis exy = ex + ey - 0x3fff; 325*25c28e83SPiotr Jasiukajtis 326*25c28e83SPiotr Jasiukajtis /* convert the significands of x and y to doubles */ 327*25c28e83SPiotr Jasiukajtis c = twom16; 328*25c28e83SPiotr Jasiukajtis dx[0] = (double) ((int) hx) * c; 329*25c28e83SPiotr Jasiukajtis dy[0] = (double) ((int) hy) * c; 330*25c28e83SPiotr Jasiukajtis 331*25c28e83SPiotr Jasiukajtis c *= twom24; 332*25c28e83SPiotr Jasiukajtis dx[1] = (double) ((int) (xx.i[1] >> 8)) * c; 333*25c28e83SPiotr Jasiukajtis dy[1] = (double) ((int) (yy.i[1] >> 8)) * c; 334*25c28e83SPiotr Jasiukajtis 335*25c28e83SPiotr Jasiukajtis c *= twom24; 336*25c28e83SPiotr Jasiukajtis dx[2] = (double) ((int) (((xx.i[1] << 16) | (xx.i[2] >> 16)) & 337*25c28e83SPiotr Jasiukajtis 0xffffff)) * c; 338*25c28e83SPiotr Jasiukajtis dy[2] = (double) ((int) (((yy.i[1] << 16) | (yy.i[2] >> 16)) & 339*25c28e83SPiotr Jasiukajtis 0xffffff)) * c; 340*25c28e83SPiotr Jasiukajtis 341*25c28e83SPiotr Jasiukajtis c *= twom24; 342*25c28e83SPiotr Jasiukajtis dx[3] = (double) ((int) (((xx.i[2] << 8) | (xx.i[3] >> 24)) & 343*25c28e83SPiotr Jasiukajtis 0xffffff)) * c; 344*25c28e83SPiotr Jasiukajtis dy[3] = (double) ((int) (((yy.i[2] << 8) | (yy.i[3] >> 24)) & 345*25c28e83SPiotr Jasiukajtis 0xffffff)) * c; 346*25c28e83SPiotr Jasiukajtis 347*25c28e83SPiotr Jasiukajtis c *= twom24; 348*25c28e83SPiotr Jasiukajtis dx[4] = (double) ((int) (xx.i[3] & 0xffffff)) * c; 349*25c28e83SPiotr Jasiukajtis dy[4] = (double) ((int) (yy.i[3] & 0xffffff)) * c; 350*25c28e83SPiotr Jasiukajtis 351*25c28e83SPiotr Jasiukajtis /* form the "digits" of the product */ 352*25c28e83SPiotr Jasiukajtis dxy[0] = dx[0] * dy[0]; 353*25c28e83SPiotr Jasiukajtis dxy[1] = dx[0] * dy[1] + dx[1] * dy[0]; 354*25c28e83SPiotr Jasiukajtis dxy[2] = dx[0] * dy[2] + dx[1] * dy[1] + dx[2] * dy[0]; 355*25c28e83SPiotr Jasiukajtis dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] + 356*25c28e83SPiotr Jasiukajtis dx[3] * dy[0]; 357*25c28e83SPiotr Jasiukajtis dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] + 358*25c28e83SPiotr Jasiukajtis dx[3] * dy[1] + dx[4] * dy[0]; 359*25c28e83SPiotr Jasiukajtis dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] + 360*25c28e83SPiotr Jasiukajtis dx[4] * dy[1]; 361*25c28e83SPiotr Jasiukajtis dxy[6] = dx[2] * dy[4] + dx[3] * dy[3] + dx[4] * dy[2]; 362*25c28e83SPiotr Jasiukajtis dxy[7] = dx[3] * dy[4] + dx[4] * dy[3]; 363*25c28e83SPiotr Jasiukajtis dxy[8] = dx[4] * dy[4]; 364*25c28e83SPiotr Jasiukajtis 365*25c28e83SPiotr Jasiukajtis /* split odd-numbered terms and combine into even-numbered terms */ 366*25c28e83SPiotr Jasiukajtis c = (dxy[1] + two20) - two20; 367*25c28e83SPiotr Jasiukajtis dxy[0] += c; 368*25c28e83SPiotr Jasiukajtis dxy[1] -= c; 369*25c28e83SPiotr Jasiukajtis c = (dxy[3] + twom28) - twom28; 370*25c28e83SPiotr Jasiukajtis dxy[2] += c + dxy[1]; 371*25c28e83SPiotr Jasiukajtis dxy[3] -= c; 372*25c28e83SPiotr Jasiukajtis c = (dxy[5] + twom76) - twom76; 373*25c28e83SPiotr Jasiukajtis dxy[4] += c + dxy[3]; 374*25c28e83SPiotr Jasiukajtis dxy[5] -= c; 375*25c28e83SPiotr Jasiukajtis c = (dxy[7] + twom124) - twom124; 376*25c28e83SPiotr Jasiukajtis dxy[6] += c + dxy[5]; 377*25c28e83SPiotr Jasiukajtis dxy[8] += (dxy[7] - c); 378*25c28e83SPiotr Jasiukajtis 379*25c28e83SPiotr Jasiukajtis /* propagate carries, adjusting the exponent if need be */ 380*25c28e83SPiotr Jasiukajtis dxy[7] = dxy[6] + dxy[8]; 381*25c28e83SPiotr Jasiukajtis dxy[5] = dxy[4] + dxy[7]; 382*25c28e83SPiotr Jasiukajtis dxy[3] = dxy[2] + dxy[5]; 383*25c28e83SPiotr Jasiukajtis dxy[1] = dxy[0] + dxy[3]; 384*25c28e83SPiotr Jasiukajtis if (dxy[1] >= two) { 385*25c28e83SPiotr Jasiukajtis dxy[0] *= half; 386*25c28e83SPiotr Jasiukajtis dxy[1] *= half; 387*25c28e83SPiotr Jasiukajtis dxy[2] *= half; 388*25c28e83SPiotr Jasiukajtis dxy[3] *= half; 389*25c28e83SPiotr Jasiukajtis dxy[4] *= half; 390*25c28e83SPiotr Jasiukajtis dxy[5] *= half; 391*25c28e83SPiotr Jasiukajtis dxy[6] *= half; 392*25c28e83SPiotr Jasiukajtis dxy[7] *= half; 393*25c28e83SPiotr Jasiukajtis dxy[8] *= half; 394*25c28e83SPiotr Jasiukajtis exy++; 395*25c28e83SPiotr Jasiukajtis } 396*25c28e83SPiotr Jasiukajtis 397*25c28e83SPiotr Jasiukajtis /* extract the significand of x*y */ 398*25c28e83SPiotr Jasiukajtis s = two36; 399*25c28e83SPiotr Jasiukajtis u.d = c = dxy[1] + s; 400*25c28e83SPiotr Jasiukajtis xy0 = u.i[1]; 401*25c28e83SPiotr Jasiukajtis c -= s; 402*25c28e83SPiotr Jasiukajtis dxy[1] -= c; 403*25c28e83SPiotr Jasiukajtis dxy[0] -= c; 404*25c28e83SPiotr Jasiukajtis 405*25c28e83SPiotr Jasiukajtis s *= twom32; 406*25c28e83SPiotr Jasiukajtis u.d = c = dxy[1] + s; 407*25c28e83SPiotr Jasiukajtis xy1 = u.i[1]; 408*25c28e83SPiotr Jasiukajtis c -= s; 409*25c28e83SPiotr Jasiukajtis dxy[2] += (dxy[0] - c); 410*25c28e83SPiotr Jasiukajtis dxy[3] = dxy[2] + dxy[5]; 411*25c28e83SPiotr Jasiukajtis 412*25c28e83SPiotr Jasiukajtis s *= twom32; 413*25c28e83SPiotr Jasiukajtis u.d = c = dxy[3] + s; 414*25c28e83SPiotr Jasiukajtis xy2 = u.i[1]; 415*25c28e83SPiotr Jasiukajtis c -= s; 416*25c28e83SPiotr Jasiukajtis dxy[4] += (dxy[2] - c); 417*25c28e83SPiotr Jasiukajtis dxy[5] = dxy[4] + dxy[7]; 418*25c28e83SPiotr Jasiukajtis 419*25c28e83SPiotr Jasiukajtis s *= twom32; 420*25c28e83SPiotr Jasiukajtis u.d = c = dxy[5] + s; 421*25c28e83SPiotr Jasiukajtis xy3 = u.i[1]; 422*25c28e83SPiotr Jasiukajtis c -= s; 423*25c28e83SPiotr Jasiukajtis dxy[4] -= c; 424*25c28e83SPiotr Jasiukajtis dxy[5] = dxy[4] + dxy[7]; 425*25c28e83SPiotr Jasiukajtis 426*25c28e83SPiotr Jasiukajtis s *= twom32; 427*25c28e83SPiotr Jasiukajtis u.d = c = dxy[5] + s; 428*25c28e83SPiotr Jasiukajtis xy4 = u.i[1]; 429*25c28e83SPiotr Jasiukajtis c -= s; 430*25c28e83SPiotr Jasiukajtis dxy[6] += (dxy[4] - c); 431*25c28e83SPiotr Jasiukajtis dxy[7] = dxy[6] + dxy[8]; 432*25c28e83SPiotr Jasiukajtis 433*25c28e83SPiotr Jasiukajtis s *= twom32; 434*25c28e83SPiotr Jasiukajtis u.d = c = dxy[7] + s; 435*25c28e83SPiotr Jasiukajtis xy5 = u.i[1]; 436*25c28e83SPiotr Jasiukajtis c -= s; 437*25c28e83SPiotr Jasiukajtis dxy[8] += (dxy[6] - c); 438*25c28e83SPiotr Jasiukajtis 439*25c28e83SPiotr Jasiukajtis s *= twom32; 440*25c28e83SPiotr Jasiukajtis u.d = c = dxy[8] + s; 441*25c28e83SPiotr Jasiukajtis xy6 = u.i[1]; 442*25c28e83SPiotr Jasiukajtis c -= s; 443*25c28e83SPiotr Jasiukajtis dxy[8] -= c; 444*25c28e83SPiotr Jasiukajtis 445*25c28e83SPiotr Jasiukajtis s *= twom32; 446*25c28e83SPiotr Jasiukajtis u.d = c = dxy[8] + s; 447*25c28e83SPiotr Jasiukajtis xy7 = u.i[1]; 448*25c28e83SPiotr Jasiukajtis 449*25c28e83SPiotr Jasiukajtis /* extract the sign, exponent, and significand of z */ 450*25c28e83SPiotr Jasiukajtis sz = zz.i[0] & 0x80000000; 451*25c28e83SPiotr Jasiukajtis ez = hz >> 16; 452*25c28e83SPiotr Jasiukajtis z0 = hz & 0xffff; 453*25c28e83SPiotr Jasiukajtis if (!ez) { 454*25c28e83SPiotr Jasiukajtis if (z0 | (zz.i[1] & 0xfffe0000)) { 455*25c28e83SPiotr Jasiukajtis z1 = zz.i[1]; 456*25c28e83SPiotr Jasiukajtis z2 = zz.i[2]; 457*25c28e83SPiotr Jasiukajtis z3 = zz.i[3]; 458*25c28e83SPiotr Jasiukajtis ez = 1; 459*25c28e83SPiotr Jasiukajtis } else if (zz.i[1] | (zz.i[2] & 0xfffe0000)) { 460*25c28e83SPiotr Jasiukajtis z0 = zz.i[1]; 461*25c28e83SPiotr Jasiukajtis z1 = zz.i[2]; 462*25c28e83SPiotr Jasiukajtis z2 = zz.i[3]; 463*25c28e83SPiotr Jasiukajtis z3 = 0; 464*25c28e83SPiotr Jasiukajtis ez = -31; 465*25c28e83SPiotr Jasiukajtis } else if (zz.i[2] | (zz.i[3] & 0xfffe0000)) { 466*25c28e83SPiotr Jasiukajtis z0 = zz.i[2]; 467*25c28e83SPiotr Jasiukajtis z1 = zz.i[3]; 468*25c28e83SPiotr Jasiukajtis z2 = z3 = 0; 469*25c28e83SPiotr Jasiukajtis ez = -63; 470*25c28e83SPiotr Jasiukajtis } else { 471*25c28e83SPiotr Jasiukajtis z0 = zz.i[3]; 472*25c28e83SPiotr Jasiukajtis z1 = z2 = z3 = 0; 473*25c28e83SPiotr Jasiukajtis ez = -95; 474*25c28e83SPiotr Jasiukajtis } 475*25c28e83SPiotr Jasiukajtis while ((z0 & 0x10000) == 0) { 476*25c28e83SPiotr Jasiukajtis z0 = (z0 << 1) | (z1 >> 31); 477*25c28e83SPiotr Jasiukajtis z1 = (z1 << 1) | (z2 >> 31); 478*25c28e83SPiotr Jasiukajtis z2 = (z2 << 1) | (z3 >> 31); 479*25c28e83SPiotr Jasiukajtis z3 <<= 1; 480*25c28e83SPiotr Jasiukajtis ez--; 481*25c28e83SPiotr Jasiukajtis } 482*25c28e83SPiotr Jasiukajtis } else { 483*25c28e83SPiotr Jasiukajtis z0 |= 0x10000; 484*25c28e83SPiotr Jasiukajtis z1 = zz.i[1]; 485*25c28e83SPiotr Jasiukajtis z2 = zz.i[2]; 486*25c28e83SPiotr Jasiukajtis z3 = zz.i[3]; 487*25c28e83SPiotr Jasiukajtis } 488*25c28e83SPiotr Jasiukajtis z4 = z5 = z6 = z7 = 0; 489*25c28e83SPiotr Jasiukajtis 490*25c28e83SPiotr Jasiukajtis /* 491*25c28e83SPiotr Jasiukajtis * now x*y is represented by sxy, exy, and xy[0-7], and z is 492*25c28e83SPiotr Jasiukajtis * represented likewise; swap if need be so |xy| <= |z| 493*25c28e83SPiotr Jasiukajtis */ 494*25c28e83SPiotr Jasiukajtis if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 || 495*25c28e83SPiotr Jasiukajtis (xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 || 496*25c28e83SPiotr Jasiukajtis (xy3 == z3 && (xy4 | xy5 | xy6 | xy7) != 0)))))))))) { 497*25c28e83SPiotr Jasiukajtis e = sxy; sxy = sz; sz = e; 498*25c28e83SPiotr Jasiukajtis e = exy; exy = ez; ez = e; 499*25c28e83SPiotr Jasiukajtis e = xy0; xy0 = z0; z0 = e; 500*25c28e83SPiotr Jasiukajtis e = xy1; xy1 = z1; z1 = e; 501*25c28e83SPiotr Jasiukajtis e = xy2; xy2 = z2; z2 = e; 502*25c28e83SPiotr Jasiukajtis e = xy3; xy3 = z3; z3 = e; 503*25c28e83SPiotr Jasiukajtis z4 = xy4; xy4 = 0; 504*25c28e83SPiotr Jasiukajtis z5 = xy5; xy5 = 0; 505*25c28e83SPiotr Jasiukajtis z6 = xy6; xy6 = 0; 506*25c28e83SPiotr Jasiukajtis z7 = xy7; xy7 = 0; 507*25c28e83SPiotr Jasiukajtis } 508*25c28e83SPiotr Jasiukajtis 509*25c28e83SPiotr Jasiukajtis /* shift the significand of xy keeping a sticky bit */ 510*25c28e83SPiotr Jasiukajtis e = ez - exy; 511*25c28e83SPiotr Jasiukajtis if (e > 236) { 512*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0; 513*25c28e83SPiotr Jasiukajtis xy7 = 1; 514*25c28e83SPiotr Jasiukajtis } else if (e >= 224) { 515*25c28e83SPiotr Jasiukajtis sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 | 516*25c28e83SPiotr Jasiukajtis ((xy0 << 1) << (255 - e)); 517*25c28e83SPiotr Jasiukajtis xy7 = xy0 >> (e - 224); 518*25c28e83SPiotr Jasiukajtis if (sticky) 519*25c28e83SPiotr Jasiukajtis xy7 |= 1; 520*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0; 521*25c28e83SPiotr Jasiukajtis } else if (e >= 192) { 522*25c28e83SPiotr Jasiukajtis sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | 523*25c28e83SPiotr Jasiukajtis ((xy1 << 1) << (223 - e)); 524*25c28e83SPiotr Jasiukajtis xy7 = (xy1 >> (e - 192)) | ((xy0 << 1) << (223 - e)); 525*25c28e83SPiotr Jasiukajtis if (sticky) 526*25c28e83SPiotr Jasiukajtis xy7 |= 1; 527*25c28e83SPiotr Jasiukajtis xy6 = xy0 >> (e - 192); 528*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = 0; 529*25c28e83SPiotr Jasiukajtis } else if (e >= 160) { 530*25c28e83SPiotr Jasiukajtis sticky = xy7 | xy6 | xy5 | xy4 | xy3 | 531*25c28e83SPiotr Jasiukajtis ((xy2 << 1) << (191 - e)); 532*25c28e83SPiotr Jasiukajtis xy7 = (xy2 >> (e - 160)) | ((xy1 << 1) << (191 - e)); 533*25c28e83SPiotr Jasiukajtis if (sticky) 534*25c28e83SPiotr Jasiukajtis xy7 |= 1; 535*25c28e83SPiotr Jasiukajtis xy6 = (xy1 >> (e - 160)) | ((xy0 << 1) << (191 - e)); 536*25c28e83SPiotr Jasiukajtis xy5 = xy0 >> (e - 160); 537*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = xy4 = 0; 538*25c28e83SPiotr Jasiukajtis } else if (e >= 128) { 539*25c28e83SPiotr Jasiukajtis sticky = xy7 | xy6 | xy5 | xy4 | ((xy3 << 1) << (159 - e)); 540*25c28e83SPiotr Jasiukajtis xy7 = (xy3 >> (e - 128)) | ((xy2 << 1) << (159 - e)); 541*25c28e83SPiotr Jasiukajtis if (sticky) 542*25c28e83SPiotr Jasiukajtis xy7 |= 1; 543*25c28e83SPiotr Jasiukajtis xy6 = (xy2 >> (e - 128)) | ((xy1 << 1) << (159 - e)); 544*25c28e83SPiotr Jasiukajtis xy5 = (xy1 >> (e - 128)) | ((xy0 << 1) << (159 - e)); 545*25c28e83SPiotr Jasiukajtis xy4 = xy0 >> (e - 128); 546*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = 0; 547*25c28e83SPiotr Jasiukajtis } else if (e >= 96) { 548*25c28e83SPiotr Jasiukajtis sticky = xy7 | xy6 | xy5 | ((xy4 << 1) << (127 - e)); 549*25c28e83SPiotr Jasiukajtis xy7 = (xy4 >> (e - 96)) | ((xy3 << 1) << (127 - e)); 550*25c28e83SPiotr Jasiukajtis if (sticky) 551*25c28e83SPiotr Jasiukajtis xy7 |= 1; 552*25c28e83SPiotr Jasiukajtis xy6 = (xy3 >> (e - 96)) | ((xy2 << 1) << (127 - e)); 553*25c28e83SPiotr Jasiukajtis xy5 = (xy2 >> (e - 96)) | ((xy1 << 1) << (127 - e)); 554*25c28e83SPiotr Jasiukajtis xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e)); 555*25c28e83SPiotr Jasiukajtis xy3 = xy0 >> (e - 96); 556*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = 0; 557*25c28e83SPiotr Jasiukajtis } else if (e >= 64) { 558*25c28e83SPiotr Jasiukajtis sticky = xy7 | xy6 | ((xy5 << 1) << (95 - e)); 559*25c28e83SPiotr Jasiukajtis xy7 = (xy5 >> (e - 64)) | ((xy4 << 1) << (95 - e)); 560*25c28e83SPiotr Jasiukajtis if (sticky) 561*25c28e83SPiotr Jasiukajtis xy7 |= 1; 562*25c28e83SPiotr Jasiukajtis xy6 = (xy4 >> (e - 64)) | ((xy3 << 1) << (95 - e)); 563*25c28e83SPiotr Jasiukajtis xy5 = (xy3 >> (e - 64)) | ((xy2 << 1) << (95 - e)); 564*25c28e83SPiotr Jasiukajtis xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e)); 565*25c28e83SPiotr Jasiukajtis xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e)); 566*25c28e83SPiotr Jasiukajtis xy2 = xy0 >> (e - 64); 567*25c28e83SPiotr Jasiukajtis xy0 = xy1 = 0; 568*25c28e83SPiotr Jasiukajtis } else if (e >= 32) { 569*25c28e83SPiotr Jasiukajtis sticky = xy7 | ((xy6 << 1) << (63 - e)); 570*25c28e83SPiotr Jasiukajtis xy7 = (xy6 >> (e - 32)) | ((xy5 << 1) << (63 - e)); 571*25c28e83SPiotr Jasiukajtis if (sticky) 572*25c28e83SPiotr Jasiukajtis xy7 |= 1; 573*25c28e83SPiotr Jasiukajtis xy6 = (xy5 >> (e - 32)) | ((xy4 << 1) << (63 - e)); 574*25c28e83SPiotr Jasiukajtis xy5 = (xy4 >> (e - 32)) | ((xy3 << 1) << (63 - e)); 575*25c28e83SPiotr Jasiukajtis xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e)); 576*25c28e83SPiotr Jasiukajtis xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e)); 577*25c28e83SPiotr Jasiukajtis xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e)); 578*25c28e83SPiotr Jasiukajtis xy1 = xy0 >> (e - 32); 579*25c28e83SPiotr Jasiukajtis xy0 = 0; 580*25c28e83SPiotr Jasiukajtis } else if (e) { 581*25c28e83SPiotr Jasiukajtis sticky = (xy7 << 1) << (31 - e); 582*25c28e83SPiotr Jasiukajtis xy7 = (xy7 >> e) | ((xy6 << 1) << (31 - e)); 583*25c28e83SPiotr Jasiukajtis if (sticky) 584*25c28e83SPiotr Jasiukajtis xy7 |= 1; 585*25c28e83SPiotr Jasiukajtis xy6 = (xy6 >> e) | ((xy5 << 1) << (31 - e)); 586*25c28e83SPiotr Jasiukajtis xy5 = (xy5 >> e) | ((xy4 << 1) << (31 - e)); 587*25c28e83SPiotr Jasiukajtis xy4 = (xy4 >> e) | ((xy3 << 1) << (31 - e)); 588*25c28e83SPiotr Jasiukajtis xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e)); 589*25c28e83SPiotr Jasiukajtis xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e)); 590*25c28e83SPiotr Jasiukajtis xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e)); 591*25c28e83SPiotr Jasiukajtis xy0 >>= e; 592*25c28e83SPiotr Jasiukajtis } 593*25c28e83SPiotr Jasiukajtis 594*25c28e83SPiotr Jasiukajtis /* if this is a magnitude subtract, negate the significand of xy */ 595*25c28e83SPiotr Jasiukajtis if (sxy ^ sz) { 596*25c28e83SPiotr Jasiukajtis xy0 = ~xy0; 597*25c28e83SPiotr Jasiukajtis xy1 = ~xy1; 598*25c28e83SPiotr Jasiukajtis xy2 = ~xy2; 599*25c28e83SPiotr Jasiukajtis xy3 = ~xy3; 600*25c28e83SPiotr Jasiukajtis xy4 = ~xy4; 601*25c28e83SPiotr Jasiukajtis xy5 = ~xy5; 602*25c28e83SPiotr Jasiukajtis xy6 = ~xy6; 603*25c28e83SPiotr Jasiukajtis xy7 = -xy7; 604*25c28e83SPiotr Jasiukajtis if (xy7 == 0) 605*25c28e83SPiotr Jasiukajtis if (++xy6 == 0) 606*25c28e83SPiotr Jasiukajtis if (++xy5 == 0) 607*25c28e83SPiotr Jasiukajtis if (++xy4 == 0) 608*25c28e83SPiotr Jasiukajtis if (++xy3 == 0) 609*25c28e83SPiotr Jasiukajtis if (++xy2 == 0) 610*25c28e83SPiotr Jasiukajtis if (++xy1 == 0) 611*25c28e83SPiotr Jasiukajtis xy0++; 612*25c28e83SPiotr Jasiukajtis } 613*25c28e83SPiotr Jasiukajtis 614*25c28e83SPiotr Jasiukajtis /* add, propagating carries */ 615*25c28e83SPiotr Jasiukajtis z7 += xy7; 616*25c28e83SPiotr Jasiukajtis e = (z7 < xy7); 617*25c28e83SPiotr Jasiukajtis z6 += xy6; 618*25c28e83SPiotr Jasiukajtis if (e) { 619*25c28e83SPiotr Jasiukajtis z6++; 620*25c28e83SPiotr Jasiukajtis e = (z6 <= xy6); 621*25c28e83SPiotr Jasiukajtis } else 622*25c28e83SPiotr Jasiukajtis e = (z6 < xy6); 623*25c28e83SPiotr Jasiukajtis z5 += xy5; 624*25c28e83SPiotr Jasiukajtis if (e) { 625*25c28e83SPiotr Jasiukajtis z5++; 626*25c28e83SPiotr Jasiukajtis e = (z5 <= xy5); 627*25c28e83SPiotr Jasiukajtis } else 628*25c28e83SPiotr Jasiukajtis e = (z5 < xy5); 629*25c28e83SPiotr Jasiukajtis z4 += xy4; 630*25c28e83SPiotr Jasiukajtis if (e) { 631*25c28e83SPiotr Jasiukajtis z4++; 632*25c28e83SPiotr Jasiukajtis e = (z4 <= xy4); 633*25c28e83SPiotr Jasiukajtis } else 634*25c28e83SPiotr Jasiukajtis e = (z4 < xy4); 635*25c28e83SPiotr Jasiukajtis z3 += xy3; 636*25c28e83SPiotr Jasiukajtis if (e) { 637*25c28e83SPiotr Jasiukajtis z3++; 638*25c28e83SPiotr Jasiukajtis e = (z3 <= xy3); 639*25c28e83SPiotr Jasiukajtis } else 640*25c28e83SPiotr Jasiukajtis e = (z3 < xy3); 641*25c28e83SPiotr Jasiukajtis z2 += xy2; 642*25c28e83SPiotr Jasiukajtis if (e) { 643*25c28e83SPiotr Jasiukajtis z2++; 644*25c28e83SPiotr Jasiukajtis e = (z2 <= xy2); 645*25c28e83SPiotr Jasiukajtis } else 646*25c28e83SPiotr Jasiukajtis e = (z2 < xy2); 647*25c28e83SPiotr Jasiukajtis z1 += xy1; 648*25c28e83SPiotr Jasiukajtis if (e) { 649*25c28e83SPiotr Jasiukajtis z1++; 650*25c28e83SPiotr Jasiukajtis e = (z1 <= xy1); 651*25c28e83SPiotr Jasiukajtis } else 652*25c28e83SPiotr Jasiukajtis e = (z1 < xy1); 653*25c28e83SPiotr Jasiukajtis z0 += xy0; 654*25c28e83SPiotr Jasiukajtis if (e) 655*25c28e83SPiotr Jasiukajtis z0++; 656*25c28e83SPiotr Jasiukajtis 657*25c28e83SPiotr Jasiukajtis /* postnormalize and collect rounding information into z4 */ 658*25c28e83SPiotr Jasiukajtis if (ez < 1) { 659*25c28e83SPiotr Jasiukajtis /* result is tiny; shift right until exponent is within range */ 660*25c28e83SPiotr Jasiukajtis e = 1 - ez; 661*25c28e83SPiotr Jasiukajtis if (e > 116) { 662*25c28e83SPiotr Jasiukajtis z4 = 1; /* result can't be exactly zero */ 663*25c28e83SPiotr Jasiukajtis z0 = z1 = z2 = z3 = 0; 664*25c28e83SPiotr Jasiukajtis } else if (e >= 96) { 665*25c28e83SPiotr Jasiukajtis sticky = z7 | z6 | z5 | z4 | z3 | z2 | 666*25c28e83SPiotr Jasiukajtis ((z1 << 1) << (127 - e)); 667*25c28e83SPiotr Jasiukajtis z4 = (z1 >> (e - 96)) | ((z0 << 1) << (127 - e)); 668*25c28e83SPiotr Jasiukajtis if (sticky) 669*25c28e83SPiotr Jasiukajtis z4 |= 1; 670*25c28e83SPiotr Jasiukajtis z3 = z0 >> (e - 96); 671*25c28e83SPiotr Jasiukajtis z0 = z1 = z2 = 0; 672*25c28e83SPiotr Jasiukajtis } else if (e >= 64) { 673*25c28e83SPiotr Jasiukajtis sticky = z7 | z6 | z5 | z4 | z3 | 674*25c28e83SPiotr Jasiukajtis ((z2 << 1) << (95 - e)); 675*25c28e83SPiotr Jasiukajtis z4 = (z2 >> (e - 64)) | ((z1 << 1) << (95 - e)); 676*25c28e83SPiotr Jasiukajtis if (sticky) 677*25c28e83SPiotr Jasiukajtis z4 |= 1; 678*25c28e83SPiotr Jasiukajtis z3 = (z1 >> (e - 64)) | ((z0 << 1) << (95 - e)); 679*25c28e83SPiotr Jasiukajtis z2 = z0 >> (e - 64); 680*25c28e83SPiotr Jasiukajtis z0 = z1 = 0; 681*25c28e83SPiotr Jasiukajtis } else if (e >= 32) { 682*25c28e83SPiotr Jasiukajtis sticky = z7 | z6 | z5 | z4 | ((z3 << 1) << (63 - e)); 683*25c28e83SPiotr Jasiukajtis z4 = (z3 >> (e - 32)) | ((z2 << 1) << (63 - e)); 684*25c28e83SPiotr Jasiukajtis if (sticky) 685*25c28e83SPiotr Jasiukajtis z4 |= 1; 686*25c28e83SPiotr Jasiukajtis z3 = (z2 >> (e - 32)) | ((z1 << 1) << (63 - e)); 687*25c28e83SPiotr Jasiukajtis z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e)); 688*25c28e83SPiotr Jasiukajtis z1 = z0 >> (e - 32); 689*25c28e83SPiotr Jasiukajtis z0 = 0; 690*25c28e83SPiotr Jasiukajtis } else { 691*25c28e83SPiotr Jasiukajtis sticky = z7 | z6 | z5 | (z4 << 1) << (31 - e); 692*25c28e83SPiotr Jasiukajtis z4 = (z4 >> e) | ((z3 << 1) << (31 - e)); 693*25c28e83SPiotr Jasiukajtis if (sticky) 694*25c28e83SPiotr Jasiukajtis z4 |= 1; 695*25c28e83SPiotr Jasiukajtis z3 = (z3 >> e) | ((z2 << 1) << (31 - e)); 696*25c28e83SPiotr Jasiukajtis z2 = (z2 >> e) | ((z1 << 1) << (31 - e)); 697*25c28e83SPiotr Jasiukajtis z1 = (z1 >> e) | ((z0 << 1) << (31 - e)); 698*25c28e83SPiotr Jasiukajtis z0 >>= e; 699*25c28e83SPiotr Jasiukajtis } 700*25c28e83SPiotr Jasiukajtis ez = 1; 701*25c28e83SPiotr Jasiukajtis } else if (z0 >= 0x20000) { 702*25c28e83SPiotr Jasiukajtis /* carry out; shift right by one */ 703*25c28e83SPiotr Jasiukajtis sticky = (z4 & 1) | z5 | z6 | z7; 704*25c28e83SPiotr Jasiukajtis z4 = (z4 >> 1) | (z3 << 31); 705*25c28e83SPiotr Jasiukajtis if (sticky) 706*25c28e83SPiotr Jasiukajtis z4 |= 1; 707*25c28e83SPiotr Jasiukajtis z3 = (z3 >> 1) | (z2 << 31); 708*25c28e83SPiotr Jasiukajtis z2 = (z2 >> 1) | (z1 << 31); 709*25c28e83SPiotr Jasiukajtis z1 = (z1 >> 1) | (z0 << 31); 710*25c28e83SPiotr Jasiukajtis z0 >>= 1; 711*25c28e83SPiotr Jasiukajtis ez++; 712*25c28e83SPiotr Jasiukajtis } else { 713*25c28e83SPiotr Jasiukajtis if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7) 714*25c28e83SPiotr Jasiukajtis != 0) { 715*25c28e83SPiotr Jasiukajtis /* 716*25c28e83SPiotr Jasiukajtis * borrow/cancellation; shift left as much as 717*25c28e83SPiotr Jasiukajtis * exponent allows 718*25c28e83SPiotr Jasiukajtis */ 719*25c28e83SPiotr Jasiukajtis while (!(z0 | (z1 & 0xfffe0000)) && ez >= 33) { 720*25c28e83SPiotr Jasiukajtis z0 = z1; 721*25c28e83SPiotr Jasiukajtis z1 = z2; 722*25c28e83SPiotr Jasiukajtis z2 = z3; 723*25c28e83SPiotr Jasiukajtis z3 = z4; 724*25c28e83SPiotr Jasiukajtis z4 = z5; 725*25c28e83SPiotr Jasiukajtis z5 = z6; 726*25c28e83SPiotr Jasiukajtis z6 = z7; 727*25c28e83SPiotr Jasiukajtis z7 = 0; 728*25c28e83SPiotr Jasiukajtis ez -= 32; 729*25c28e83SPiotr Jasiukajtis } 730*25c28e83SPiotr Jasiukajtis while (z0 < 0x10000 && ez > 1) { 731*25c28e83SPiotr Jasiukajtis z0 = (z0 << 1) | (z1 >> 31); 732*25c28e83SPiotr Jasiukajtis z1 = (z1 << 1) | (z2 >> 31); 733*25c28e83SPiotr Jasiukajtis z2 = (z2 << 1) | (z3 >> 31); 734*25c28e83SPiotr Jasiukajtis z3 = (z3 << 1) | (z4 >> 31); 735*25c28e83SPiotr Jasiukajtis z4 = (z4 << 1) | (z5 >> 31); 736*25c28e83SPiotr Jasiukajtis z5 = (z5 << 1) | (z6 >> 31); 737*25c28e83SPiotr Jasiukajtis z6 = (z6 << 1) | (z7 >> 31); 738*25c28e83SPiotr Jasiukajtis z7 <<= 1; 739*25c28e83SPiotr Jasiukajtis ez--; 740*25c28e83SPiotr Jasiukajtis } 741*25c28e83SPiotr Jasiukajtis } 742*25c28e83SPiotr Jasiukajtis if (z5 | z6 | z7) 743*25c28e83SPiotr Jasiukajtis z4 |= 1; 744*25c28e83SPiotr Jasiukajtis } 745*25c28e83SPiotr Jasiukajtis 746*25c28e83SPiotr Jasiukajtis /* get the rounding mode */ 747*25c28e83SPiotr Jasiukajtis rm = fsr >> 30; 748*25c28e83SPiotr Jasiukajtis 749*25c28e83SPiotr Jasiukajtis /* strip off the integer bit, if there is one */ 750*25c28e83SPiotr Jasiukajtis ibit = z0 & 0x10000; 751*25c28e83SPiotr Jasiukajtis if (ibit) 752*25c28e83SPiotr Jasiukajtis z0 -= 0x10000; 753*25c28e83SPiotr Jasiukajtis else { 754*25c28e83SPiotr Jasiukajtis ez = 0; 755*25c28e83SPiotr Jasiukajtis if (!(z0 | z1 | z2 | z3 | z4)) { /* exact zero */ 756*25c28e83SPiotr Jasiukajtis zz.i[0] = rm == FSR_RM ? 0x80000000 : 0; 757*25c28e83SPiotr Jasiukajtis zz.i[1] = zz.i[2] = zz.i[3] = 0; 758*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 759*25c28e83SPiotr Jasiukajtis return (zz.q); 760*25c28e83SPiotr Jasiukajtis } 761*25c28e83SPiotr Jasiukajtis } 762*25c28e83SPiotr Jasiukajtis 763*25c28e83SPiotr Jasiukajtis /* 764*25c28e83SPiotr Jasiukajtis * flip the sense of directed roundings if the result is negative; 765*25c28e83SPiotr Jasiukajtis * the logic below applies to a positive result 766*25c28e83SPiotr Jasiukajtis */ 767*25c28e83SPiotr Jasiukajtis if (sz) 768*25c28e83SPiotr Jasiukajtis rm ^= rm >> 1; 769*25c28e83SPiotr Jasiukajtis 770*25c28e83SPiotr Jasiukajtis /* round and raise exceptions */ 771*25c28e83SPiotr Jasiukajtis if (z4) { 772*25c28e83SPiotr Jasiukajtis fsr |= FSR_NXC; 773*25c28e83SPiotr Jasiukajtis 774*25c28e83SPiotr Jasiukajtis /* decide whether to round the fraction up */ 775*25c28e83SPiotr Jasiukajtis if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u || 776*25c28e83SPiotr Jasiukajtis (z4 == 0x80000000u && (z3 & 1))))) { 777*25c28e83SPiotr Jasiukajtis /* round up and renormalize if necessary */ 778*25c28e83SPiotr Jasiukajtis if (++z3 == 0) 779*25c28e83SPiotr Jasiukajtis if (++z2 == 0) 780*25c28e83SPiotr Jasiukajtis if (++z1 == 0) 781*25c28e83SPiotr Jasiukajtis if (++z0 == 0x10000) { 782*25c28e83SPiotr Jasiukajtis z0 = 0; 783*25c28e83SPiotr Jasiukajtis ez++; 784*25c28e83SPiotr Jasiukajtis } 785*25c28e83SPiotr Jasiukajtis } 786*25c28e83SPiotr Jasiukajtis } 787*25c28e83SPiotr Jasiukajtis 788*25c28e83SPiotr Jasiukajtis /* check for under/overflow */ 789*25c28e83SPiotr Jasiukajtis if (ez >= 0x7fff) { 790*25c28e83SPiotr Jasiukajtis if (rm == FSR_RN || rm == FSR_RP) { 791*25c28e83SPiotr Jasiukajtis zz.i[0] = sz | 0x7fff0000; 792*25c28e83SPiotr Jasiukajtis zz.i[1] = zz.i[2] = zz.i[3] = 0; 793*25c28e83SPiotr Jasiukajtis } else { 794*25c28e83SPiotr Jasiukajtis zz.i[0] = sz | 0x7ffeffff; 795*25c28e83SPiotr Jasiukajtis zz.i[1] = zz.i[2] = zz.i[3] = 0xffffffff; 796*25c28e83SPiotr Jasiukajtis } 797*25c28e83SPiotr Jasiukajtis fsr |= FSR_OFC | FSR_NXC; 798*25c28e83SPiotr Jasiukajtis } else { 799*25c28e83SPiotr Jasiukajtis zz.i[0] = sz | (ez << 16) | z0; 800*25c28e83SPiotr Jasiukajtis zz.i[1] = z1; 801*25c28e83SPiotr Jasiukajtis zz.i[2] = z2; 802*25c28e83SPiotr Jasiukajtis zz.i[3] = z3; 803*25c28e83SPiotr Jasiukajtis 804*25c28e83SPiotr Jasiukajtis /* 805*25c28e83SPiotr Jasiukajtis * !ibit => exact result was tiny before rounding, 806*25c28e83SPiotr Jasiukajtis * z4 nonzero => result delivered is inexact 807*25c28e83SPiotr Jasiukajtis */ 808*25c28e83SPiotr Jasiukajtis if (!ibit) { 809*25c28e83SPiotr Jasiukajtis if (z4) 810*25c28e83SPiotr Jasiukajtis fsr |= FSR_UFC | FSR_NXC; 811*25c28e83SPiotr Jasiukajtis else if (fsr & FSR_UFM) 812*25c28e83SPiotr Jasiukajtis fsr |= FSR_UFC; 813*25c28e83SPiotr Jasiukajtis } 814*25c28e83SPiotr Jasiukajtis } 815*25c28e83SPiotr Jasiukajtis 816*25c28e83SPiotr Jasiukajtis /* restore the fsr and emulate exceptions as needed */ 817*25c28e83SPiotr Jasiukajtis if ((fsr & FSR_CEXC) & (fsr >> 23)) { 818*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 819*25c28e83SPiotr Jasiukajtis if (fsr & FSR_OFC) { 820*25c28e83SPiotr Jasiukajtis dummy = huge; 821*25c28e83SPiotr Jasiukajtis dummy *= huge; 822*25c28e83SPiotr Jasiukajtis } else if (fsr & FSR_UFC) { 823*25c28e83SPiotr Jasiukajtis dummy = tiny; 824*25c28e83SPiotr Jasiukajtis if (fsr & FSR_NXC) 825*25c28e83SPiotr Jasiukajtis dummy *= tiny; 826*25c28e83SPiotr Jasiukajtis else 827*25c28e83SPiotr Jasiukajtis dummy -= tiny2; 828*25c28e83SPiotr Jasiukajtis } else { 829*25c28e83SPiotr Jasiukajtis dummy = huge; 830*25c28e83SPiotr Jasiukajtis dummy += tiny; 831*25c28e83SPiotr Jasiukajtis } 832*25c28e83SPiotr Jasiukajtis } else { 833*25c28e83SPiotr Jasiukajtis fsr |= (fsr & 0x1f) << 5; 834*25c28e83SPiotr Jasiukajtis __fenv_setfsr32(&fsr); 835*25c28e83SPiotr Jasiukajtis } 836*25c28e83SPiotr Jasiukajtis return (zz.q); 837*25c28e83SPiotr Jasiukajtis } 838*25c28e83SPiotr Jasiukajtis 839*25c28e83SPiotr Jasiukajtis #elif defined(__x86) 840*25c28e83SPiotr Jasiukajtis 841*25c28e83SPiotr Jasiukajtis static const union { 842*25c28e83SPiotr Jasiukajtis unsigned i[2]; 843*25c28e83SPiotr Jasiukajtis double d; 844*25c28e83SPiotr Jasiukajtis } C[] = { 845*25c28e83SPiotr Jasiukajtis { 0, 0x3fe00000u }, 846*25c28e83SPiotr Jasiukajtis { 0, 0x40000000u }, 847*25c28e83SPiotr Jasiukajtis { 0, 0x3df00000u }, 848*25c28e83SPiotr Jasiukajtis { 0, 0x3bf00000u }, 849*25c28e83SPiotr Jasiukajtis { 0, 0x41f00000u }, 850*25c28e83SPiotr Jasiukajtis { 0, 0x43e00000u }, 851*25c28e83SPiotr Jasiukajtis { 0, 0x7fe00000u }, 852*25c28e83SPiotr Jasiukajtis { 0, 0x00100000u }, 853*25c28e83SPiotr Jasiukajtis { 0, 0x00100001u } 854*25c28e83SPiotr Jasiukajtis }; 855*25c28e83SPiotr Jasiukajtis 856*25c28e83SPiotr Jasiukajtis #define half C[0].d 857*25c28e83SPiotr Jasiukajtis #define two C[1].d 858*25c28e83SPiotr Jasiukajtis #define twom32 C[2].d 859*25c28e83SPiotr Jasiukajtis #define twom64 C[3].d 860*25c28e83SPiotr Jasiukajtis #define two32 C[4].d 861*25c28e83SPiotr Jasiukajtis #define two63 C[5].d 862*25c28e83SPiotr Jasiukajtis #define huge C[6].d 863*25c28e83SPiotr Jasiukajtis #define tiny C[7].d 864*25c28e83SPiotr Jasiukajtis #define tiny2 C[8].d 865*25c28e83SPiotr Jasiukajtis 866*25c28e83SPiotr Jasiukajtis #if defined(__amd64) 867*25c28e83SPiotr Jasiukajtis #define NI 4 868*25c28e83SPiotr Jasiukajtis #else 869*25c28e83SPiotr Jasiukajtis #define NI 3 870*25c28e83SPiotr Jasiukajtis #endif 871*25c28e83SPiotr Jasiukajtis 872*25c28e83SPiotr Jasiukajtis /* 873*25c28e83SPiotr Jasiukajtis * fmal for x86: 80-bit extended double precision, little-endian 874*25c28e83SPiotr Jasiukajtis */ 875*25c28e83SPiotr Jasiukajtis long double 876*25c28e83SPiotr Jasiukajtis __fmal(long double x, long double y, long double z) { 877*25c28e83SPiotr Jasiukajtis union { 878*25c28e83SPiotr Jasiukajtis unsigned i[NI]; 879*25c28e83SPiotr Jasiukajtis long double e; 880*25c28e83SPiotr Jasiukajtis } xx, yy, zz; 881*25c28e83SPiotr Jasiukajtis long double xhi, yhi, xlo, ylo, t; 882*25c28e83SPiotr Jasiukajtis unsigned xy0, xy1, xy2, xy3, xy4, z0, z1, z2, z3, z4; 883*25c28e83SPiotr Jasiukajtis unsigned oldcwsw, cwsw, rm, sticky, carry; 884*25c28e83SPiotr Jasiukajtis int ex, ey, ez, exy, sxy, sz, e, tinyafter; 885*25c28e83SPiotr Jasiukajtis volatile double dummy; 886*25c28e83SPiotr Jasiukajtis 887*25c28e83SPiotr Jasiukajtis /* extract the exponents of the arguments */ 888*25c28e83SPiotr Jasiukajtis xx.e = x; 889*25c28e83SPiotr Jasiukajtis yy.e = y; 890*25c28e83SPiotr Jasiukajtis zz.e = z; 891*25c28e83SPiotr Jasiukajtis ex = xx.i[2] & 0x7fff; 892*25c28e83SPiotr Jasiukajtis ey = yy.i[2] & 0x7fff; 893*25c28e83SPiotr Jasiukajtis ez = zz.i[2] & 0x7fff; 894*25c28e83SPiotr Jasiukajtis 895*25c28e83SPiotr Jasiukajtis /* dispense with inf, nan, and zero cases */ 896*25c28e83SPiotr Jasiukajtis if (ex == 0x7fff || ey == 0x7fff || (ex | xx.i[1] | xx.i[0]) == 0 || 897*25c28e83SPiotr Jasiukajtis (ey | yy.i[1] | yy.i[0]) == 0) /* x or y is inf, nan, or 0 */ 898*25c28e83SPiotr Jasiukajtis return (x * y + z); 899*25c28e83SPiotr Jasiukajtis 900*25c28e83SPiotr Jasiukajtis if (ez == 0x7fff) /* z is inf or nan */ 901*25c28e83SPiotr Jasiukajtis return (x + z); /* avoid spurious under/overflow in x * y */ 902*25c28e83SPiotr Jasiukajtis 903*25c28e83SPiotr Jasiukajtis if ((ez | zz.i[1] | zz.i[0]) == 0) /* z is zero */ 904*25c28e83SPiotr Jasiukajtis /* 905*25c28e83SPiotr Jasiukajtis * x * y isn't zero but could underflow to zero, 906*25c28e83SPiotr Jasiukajtis * so don't add z, lest we perturb the sign 907*25c28e83SPiotr Jasiukajtis */ 908*25c28e83SPiotr Jasiukajtis return (x * y); 909*25c28e83SPiotr Jasiukajtis 910*25c28e83SPiotr Jasiukajtis /* 911*25c28e83SPiotr Jasiukajtis * now x, y, and z are all finite and nonzero; extract signs and 912*25c28e83SPiotr Jasiukajtis * normalize the significands (this will raise the denormal operand 913*25c28e83SPiotr Jasiukajtis * exception if need be) 914*25c28e83SPiotr Jasiukajtis */ 915*25c28e83SPiotr Jasiukajtis sxy = (xx.i[2] ^ yy.i[2]) & 0x8000; 916*25c28e83SPiotr Jasiukajtis sz = zz.i[2] & 0x8000; 917*25c28e83SPiotr Jasiukajtis if (!ex) { 918*25c28e83SPiotr Jasiukajtis xx.e = x * two63; 919*25c28e83SPiotr Jasiukajtis ex = (xx.i[2] & 0x7fff) - 63; 920*25c28e83SPiotr Jasiukajtis } 921*25c28e83SPiotr Jasiukajtis if (!ey) { 922*25c28e83SPiotr Jasiukajtis yy.e = y * two63; 923*25c28e83SPiotr Jasiukajtis ey = (yy.i[2] & 0x7fff) - 63; 924*25c28e83SPiotr Jasiukajtis } 925*25c28e83SPiotr Jasiukajtis if (!ez) { 926*25c28e83SPiotr Jasiukajtis zz.e = z * two63; 927*25c28e83SPiotr Jasiukajtis ez = (zz.i[2] & 0x7fff) - 63; 928*25c28e83SPiotr Jasiukajtis } 929*25c28e83SPiotr Jasiukajtis 930*25c28e83SPiotr Jasiukajtis /* 931*25c28e83SPiotr Jasiukajtis * save the control and status words, mask all exceptions, and 932*25c28e83SPiotr Jasiukajtis * set rounding to 64-bit precision and toward-zero 933*25c28e83SPiotr Jasiukajtis */ 934*25c28e83SPiotr Jasiukajtis __fenv_getcwsw(&oldcwsw); 935*25c28e83SPiotr Jasiukajtis cwsw = (oldcwsw & 0xf0c0ffff) | 0x0f3f0000; 936*25c28e83SPiotr Jasiukajtis __fenv_setcwsw(&cwsw); 937*25c28e83SPiotr Jasiukajtis 938*25c28e83SPiotr Jasiukajtis /* multiply x*y to 128 bits */ 939*25c28e83SPiotr Jasiukajtis exy = ex + ey - 0x3fff; 940*25c28e83SPiotr Jasiukajtis xx.i[2] = 0x3fff; 941*25c28e83SPiotr Jasiukajtis yy.i[2] = 0x3fff; 942*25c28e83SPiotr Jasiukajtis x = xx.e; 943*25c28e83SPiotr Jasiukajtis y = yy.e; 944*25c28e83SPiotr Jasiukajtis xhi = ((x + twom32) + two32) - two32; 945*25c28e83SPiotr Jasiukajtis yhi = ((y + twom32) + two32) - two32; 946*25c28e83SPiotr Jasiukajtis xlo = x - xhi; 947*25c28e83SPiotr Jasiukajtis ylo = y - yhi; 948*25c28e83SPiotr Jasiukajtis x *= y; 949*25c28e83SPiotr Jasiukajtis y = ((xhi * yhi - x) + xhi * ylo + xlo * yhi) + xlo * ylo; 950*25c28e83SPiotr Jasiukajtis if (x >= two) { 951*25c28e83SPiotr Jasiukajtis x *= half; 952*25c28e83SPiotr Jasiukajtis y *= half; 953*25c28e83SPiotr Jasiukajtis exy++; 954*25c28e83SPiotr Jasiukajtis } 955*25c28e83SPiotr Jasiukajtis 956*25c28e83SPiotr Jasiukajtis /* extract the significands */ 957*25c28e83SPiotr Jasiukajtis xx.e = x; 958*25c28e83SPiotr Jasiukajtis xy0 = xx.i[1]; 959*25c28e83SPiotr Jasiukajtis xy1 = xx.i[0]; 960*25c28e83SPiotr Jasiukajtis yy.e = t = y + twom32; 961*25c28e83SPiotr Jasiukajtis xy2 = yy.i[0]; 962*25c28e83SPiotr Jasiukajtis yy.e = (y - (t - twom32)) + twom64; 963*25c28e83SPiotr Jasiukajtis xy3 = yy.i[0]; 964*25c28e83SPiotr Jasiukajtis xy4 = 0; 965*25c28e83SPiotr Jasiukajtis z0 = zz.i[1]; 966*25c28e83SPiotr Jasiukajtis z1 = zz.i[0]; 967*25c28e83SPiotr Jasiukajtis z2 = z3 = z4 = 0; 968*25c28e83SPiotr Jasiukajtis 969*25c28e83SPiotr Jasiukajtis /* 970*25c28e83SPiotr Jasiukajtis * now x*y is represented by sxy, exy, and xy[0-4], and z is 971*25c28e83SPiotr Jasiukajtis * represented likewise; swap if need be so |xy| <= |z| 972*25c28e83SPiotr Jasiukajtis */ 973*25c28e83SPiotr Jasiukajtis if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && 974*25c28e83SPiotr Jasiukajtis (xy1 > z1 || (xy1 == z1 && (xy2 | xy3) != 0)))))) { 975*25c28e83SPiotr Jasiukajtis e = sxy; sxy = sz; sz = e; 976*25c28e83SPiotr Jasiukajtis e = exy; exy = ez; ez = e; 977*25c28e83SPiotr Jasiukajtis e = xy0; xy0 = z0; z0 = e; 978*25c28e83SPiotr Jasiukajtis e = xy1; xy1 = z1; z1 = e; 979*25c28e83SPiotr Jasiukajtis z2 = xy2; xy2 = 0; 980*25c28e83SPiotr Jasiukajtis z3 = xy3; xy3 = 0; 981*25c28e83SPiotr Jasiukajtis } 982*25c28e83SPiotr Jasiukajtis 983*25c28e83SPiotr Jasiukajtis /* shift the significand of xy keeping a sticky bit */ 984*25c28e83SPiotr Jasiukajtis e = ez - exy; 985*25c28e83SPiotr Jasiukajtis if (e > 130) { 986*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = 0; 987*25c28e83SPiotr Jasiukajtis xy4 = 1; 988*25c28e83SPiotr Jasiukajtis } else if (e >= 128) { 989*25c28e83SPiotr Jasiukajtis sticky = xy3 | xy2 | xy1 | ((xy0 << 1) << (159 - e)); 990*25c28e83SPiotr Jasiukajtis xy4 = xy0 >> (e - 128); 991*25c28e83SPiotr Jasiukajtis if (sticky) 992*25c28e83SPiotr Jasiukajtis xy4 |= 1; 993*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = 0; 994*25c28e83SPiotr Jasiukajtis } else if (e >= 96) { 995*25c28e83SPiotr Jasiukajtis sticky = xy3 | xy2 | ((xy1 << 1) << (127 - e)); 996*25c28e83SPiotr Jasiukajtis xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e)); 997*25c28e83SPiotr Jasiukajtis if (sticky) 998*25c28e83SPiotr Jasiukajtis xy4 |= 1; 999*25c28e83SPiotr Jasiukajtis xy3 = xy0 >> (e - 96); 1000*25c28e83SPiotr Jasiukajtis xy0 = xy1 = xy2 = 0; 1001*25c28e83SPiotr Jasiukajtis } else if (e >= 64) { 1002*25c28e83SPiotr Jasiukajtis sticky = xy3 | ((xy2 << 1) << (95 - e)); 1003*25c28e83SPiotr Jasiukajtis xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e)); 1004*25c28e83SPiotr Jasiukajtis if (sticky) 1005*25c28e83SPiotr Jasiukajtis xy4 |= 1; 1006*25c28e83SPiotr Jasiukajtis xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e)); 1007*25c28e83SPiotr Jasiukajtis xy2 = xy0 >> (e - 64); 1008*25c28e83SPiotr Jasiukajtis xy0 = xy1 = 0; 1009*25c28e83SPiotr Jasiukajtis } else if (e >= 32) { 1010*25c28e83SPiotr Jasiukajtis sticky = (xy3 << 1) << (63 - e); 1011*25c28e83SPiotr Jasiukajtis xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e)); 1012*25c28e83SPiotr Jasiukajtis if (sticky) 1013*25c28e83SPiotr Jasiukajtis xy4 |= 1; 1014*25c28e83SPiotr Jasiukajtis xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e)); 1015*25c28e83SPiotr Jasiukajtis xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e)); 1016*25c28e83SPiotr Jasiukajtis xy1 = xy0 >> (e - 32); 1017*25c28e83SPiotr Jasiukajtis xy0 = 0; 1018*25c28e83SPiotr Jasiukajtis } else if (e) { 1019*25c28e83SPiotr Jasiukajtis xy4 = (xy3 << 1) << (31 - e); 1020*25c28e83SPiotr Jasiukajtis xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e)); 1021*25c28e83SPiotr Jasiukajtis xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e)); 1022*25c28e83SPiotr Jasiukajtis xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e)); 1023*25c28e83SPiotr Jasiukajtis xy0 >>= e; 1024*25c28e83SPiotr Jasiukajtis } 1025*25c28e83SPiotr Jasiukajtis 1026*25c28e83SPiotr Jasiukajtis /* if this is a magnitude subtract, negate the significand of xy */ 1027*25c28e83SPiotr Jasiukajtis if (sxy ^ sz) { 1028*25c28e83SPiotr Jasiukajtis xy0 = ~xy0; 1029*25c28e83SPiotr Jasiukajtis xy1 = ~xy1; 1030*25c28e83SPiotr Jasiukajtis xy2 = ~xy2; 1031*25c28e83SPiotr Jasiukajtis xy3 = ~xy3; 1032*25c28e83SPiotr Jasiukajtis xy4 = -xy4; 1033*25c28e83SPiotr Jasiukajtis if (xy4 == 0) 1034*25c28e83SPiotr Jasiukajtis if (++xy3 == 0) 1035*25c28e83SPiotr Jasiukajtis if (++xy2 == 0) 1036*25c28e83SPiotr Jasiukajtis if (++xy1 == 0) 1037*25c28e83SPiotr Jasiukajtis xy0++; 1038*25c28e83SPiotr Jasiukajtis } 1039*25c28e83SPiotr Jasiukajtis 1040*25c28e83SPiotr Jasiukajtis /* add, propagating carries */ 1041*25c28e83SPiotr Jasiukajtis z4 += xy4; 1042*25c28e83SPiotr Jasiukajtis carry = (z4 < xy4); 1043*25c28e83SPiotr Jasiukajtis z3 += xy3; 1044*25c28e83SPiotr Jasiukajtis if (carry) { 1045*25c28e83SPiotr Jasiukajtis z3++; 1046*25c28e83SPiotr Jasiukajtis carry = (z3 <= xy3); 1047*25c28e83SPiotr Jasiukajtis } else 1048*25c28e83SPiotr Jasiukajtis carry = (z3 < xy3); 1049*25c28e83SPiotr Jasiukajtis z2 += xy2; 1050*25c28e83SPiotr Jasiukajtis if (carry) { 1051*25c28e83SPiotr Jasiukajtis z2++; 1052*25c28e83SPiotr Jasiukajtis carry = (z2 <= xy2); 1053*25c28e83SPiotr Jasiukajtis } else 1054*25c28e83SPiotr Jasiukajtis carry = (z2 < xy2); 1055*25c28e83SPiotr Jasiukajtis z1 += xy1; 1056*25c28e83SPiotr Jasiukajtis if (carry) { 1057*25c28e83SPiotr Jasiukajtis z1++; 1058*25c28e83SPiotr Jasiukajtis carry = (z1 <= xy1); 1059*25c28e83SPiotr Jasiukajtis } else 1060*25c28e83SPiotr Jasiukajtis carry = (z1 < xy1); 1061*25c28e83SPiotr Jasiukajtis z0 += xy0; 1062*25c28e83SPiotr Jasiukajtis if (carry) { 1063*25c28e83SPiotr Jasiukajtis z0++; 1064*25c28e83SPiotr Jasiukajtis carry = (z0 <= xy0); 1065*25c28e83SPiotr Jasiukajtis } else 1066*25c28e83SPiotr Jasiukajtis carry = (z0 < xy0); 1067*25c28e83SPiotr Jasiukajtis 1068*25c28e83SPiotr Jasiukajtis /* for a magnitude subtract, ignore the last carry out */ 1069*25c28e83SPiotr Jasiukajtis if (sxy ^ sz) 1070*25c28e83SPiotr Jasiukajtis carry = 0; 1071*25c28e83SPiotr Jasiukajtis 1072*25c28e83SPiotr Jasiukajtis /* postnormalize and collect rounding information into z2 */ 1073*25c28e83SPiotr Jasiukajtis if (ez < 1) { 1074*25c28e83SPiotr Jasiukajtis /* result is tiny; shift right until exponent is within range */ 1075*25c28e83SPiotr Jasiukajtis e = 1 - ez; 1076*25c28e83SPiotr Jasiukajtis if (e > 67) { 1077*25c28e83SPiotr Jasiukajtis z2 = 1; /* result can't be exactly zero */ 1078*25c28e83SPiotr Jasiukajtis z0 = z1 = 0; 1079*25c28e83SPiotr Jasiukajtis } else if (e >= 64) { 1080*25c28e83SPiotr Jasiukajtis sticky = z4 | z3 | z2 | z1 | ((z0 << 1) << (95 - e)); 1081*25c28e83SPiotr Jasiukajtis z2 = (z0 >> (e - 64)) | ((carry << 1) << (95 - e)); 1082*25c28e83SPiotr Jasiukajtis if (sticky) 1083*25c28e83SPiotr Jasiukajtis z2 |= 1; 1084*25c28e83SPiotr Jasiukajtis z1 = carry >> (e - 64); 1085*25c28e83SPiotr Jasiukajtis z0 = 0; 1086*25c28e83SPiotr Jasiukajtis } else if (e >= 32) { 1087*25c28e83SPiotr Jasiukajtis sticky = z4 | z3 | z2 | ((z1 << 1) << (63 - e)); 1088*25c28e83SPiotr Jasiukajtis z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e)); 1089*25c28e83SPiotr Jasiukajtis if (sticky) 1090*25c28e83SPiotr Jasiukajtis z2 |= 1; 1091*25c28e83SPiotr Jasiukajtis z1 = (z0 >> (e - 32)) | ((carry << 1) << (63 - e)); 1092*25c28e83SPiotr Jasiukajtis z0 = carry >> (e - 32); 1093*25c28e83SPiotr Jasiukajtis } else { 1094*25c28e83SPiotr Jasiukajtis sticky = z4 | z3 | (z2 << 1) << (31 - e); 1095*25c28e83SPiotr Jasiukajtis z2 = (z2 >> e) | ((z1 << 1) << (31 - e)); 1096*25c28e83SPiotr Jasiukajtis if (sticky) 1097*25c28e83SPiotr Jasiukajtis z2 |= 1; 1098*25c28e83SPiotr Jasiukajtis z1 = (z1 >> e) | ((z0 << 1) << (31 - e)); 1099*25c28e83SPiotr Jasiukajtis z0 = (z0 >> e) | ((carry << 1) << (31 - e)); 1100*25c28e83SPiotr Jasiukajtis } 1101*25c28e83SPiotr Jasiukajtis ez = 1; 1102*25c28e83SPiotr Jasiukajtis } else if (carry) { 1103*25c28e83SPiotr Jasiukajtis /* carry out; shift right by one */ 1104*25c28e83SPiotr Jasiukajtis sticky = (z2 & 1) | z3 | z4; 1105*25c28e83SPiotr Jasiukajtis z2 = (z2 >> 1) | (z1 << 31); 1106*25c28e83SPiotr Jasiukajtis if (sticky) 1107*25c28e83SPiotr Jasiukajtis z2 |= 1; 1108*25c28e83SPiotr Jasiukajtis z1 = (z1 >> 1) | (z0 << 31); 1109*25c28e83SPiotr Jasiukajtis z0 = (z0 >> 1) | 0x80000000; 1110*25c28e83SPiotr Jasiukajtis ez++; 1111*25c28e83SPiotr Jasiukajtis } else { 1112*25c28e83SPiotr Jasiukajtis if (z0 < 0x80000000u && (z0 | z1 | z2 | z3 | z4) != 0) { 1113*25c28e83SPiotr Jasiukajtis /* 1114*25c28e83SPiotr Jasiukajtis * borrow/cancellation; shift left as much as 1115*25c28e83SPiotr Jasiukajtis * exponent allows 1116*25c28e83SPiotr Jasiukajtis */ 1117*25c28e83SPiotr Jasiukajtis while (!z0 && ez >= 33) { 1118*25c28e83SPiotr Jasiukajtis z0 = z1; 1119*25c28e83SPiotr Jasiukajtis z1 = z2; 1120*25c28e83SPiotr Jasiukajtis z2 = z3; 1121*25c28e83SPiotr Jasiukajtis z3 = z4; 1122*25c28e83SPiotr Jasiukajtis z4 = 0; 1123*25c28e83SPiotr Jasiukajtis ez -= 32; 1124*25c28e83SPiotr Jasiukajtis } 1125*25c28e83SPiotr Jasiukajtis while (z0 < 0x80000000u && ez > 1) { 1126*25c28e83SPiotr Jasiukajtis z0 = (z0 << 1) | (z1 >> 31); 1127*25c28e83SPiotr Jasiukajtis z1 = (z1 << 1) | (z2 >> 31); 1128*25c28e83SPiotr Jasiukajtis z2 = (z2 << 1) | (z3 >> 31); 1129*25c28e83SPiotr Jasiukajtis z3 = (z3 << 1) | (z4 >> 31); 1130*25c28e83SPiotr Jasiukajtis z4 <<= 1; 1131*25c28e83SPiotr Jasiukajtis ez--; 1132*25c28e83SPiotr Jasiukajtis } 1133*25c28e83SPiotr Jasiukajtis } 1134*25c28e83SPiotr Jasiukajtis if (z3 | z4) 1135*25c28e83SPiotr Jasiukajtis z2 |= 1; 1136*25c28e83SPiotr Jasiukajtis } 1137*25c28e83SPiotr Jasiukajtis 1138*25c28e83SPiotr Jasiukajtis /* get the rounding mode */ 1139*25c28e83SPiotr Jasiukajtis rm = oldcwsw & 0x0c000000; 1140*25c28e83SPiotr Jasiukajtis 1141*25c28e83SPiotr Jasiukajtis /* adjust exponent if result is subnormal */ 1142*25c28e83SPiotr Jasiukajtis tinyafter = 0; 1143*25c28e83SPiotr Jasiukajtis if (!(z0 & 0x80000000)) { 1144*25c28e83SPiotr Jasiukajtis ez = 0; 1145*25c28e83SPiotr Jasiukajtis tinyafter = 1; 1146*25c28e83SPiotr Jasiukajtis if (!(z0 | z1 | z2)) { /* exact zero */ 1147*25c28e83SPiotr Jasiukajtis zz.i[2] = rm == FCW_RM ? 0x8000 : 0; 1148*25c28e83SPiotr Jasiukajtis zz.i[1] = zz.i[0] = 0; 1149*25c28e83SPiotr Jasiukajtis __fenv_setcwsw(&oldcwsw); 1150*25c28e83SPiotr Jasiukajtis return (zz.e); 1151*25c28e83SPiotr Jasiukajtis } 1152*25c28e83SPiotr Jasiukajtis } 1153*25c28e83SPiotr Jasiukajtis 1154*25c28e83SPiotr Jasiukajtis /* 1155*25c28e83SPiotr Jasiukajtis * flip the sense of directed roundings if the result is negative; 1156*25c28e83SPiotr Jasiukajtis * the logic below applies to a positive result 1157*25c28e83SPiotr Jasiukajtis */ 1158*25c28e83SPiotr Jasiukajtis if (sz && (rm == FCW_RM || rm == FCW_RP)) 1159*25c28e83SPiotr Jasiukajtis rm = (FCW_RM + FCW_RP) - rm; 1160*25c28e83SPiotr Jasiukajtis 1161*25c28e83SPiotr Jasiukajtis /* round */ 1162*25c28e83SPiotr Jasiukajtis if (z2) { 1163*25c28e83SPiotr Jasiukajtis if (rm == FCW_RP || (rm == FCW_RN && (z2 > 0x80000000u || 1164*25c28e83SPiotr Jasiukajtis (z2 == 0x80000000u && (z1 & 1))))) { 1165*25c28e83SPiotr Jasiukajtis /* round up and renormalize if necessary */ 1166*25c28e83SPiotr Jasiukajtis if (++z1 == 0) { 1167*25c28e83SPiotr Jasiukajtis if (++z0 == 0) { 1168*25c28e83SPiotr Jasiukajtis z0 = 0x80000000; 1169*25c28e83SPiotr Jasiukajtis ez++; 1170*25c28e83SPiotr Jasiukajtis } else if (z0 == 0x80000000) { 1171*25c28e83SPiotr Jasiukajtis /* rounded up to smallest normal */ 1172*25c28e83SPiotr Jasiukajtis ez = 1; 1173*25c28e83SPiotr Jasiukajtis if ((rm == FCW_RP && z2 > 1174*25c28e83SPiotr Jasiukajtis 0x80000000u) || (rm == FCW_RN && 1175*25c28e83SPiotr Jasiukajtis z2 >= 0xc0000000u)) 1176*25c28e83SPiotr Jasiukajtis /* 1177*25c28e83SPiotr Jasiukajtis * would have rounded up to 1178*25c28e83SPiotr Jasiukajtis * smallest normal even with 1179*25c28e83SPiotr Jasiukajtis * unbounded range 1180*25c28e83SPiotr Jasiukajtis */ 1181*25c28e83SPiotr Jasiukajtis tinyafter = 0; 1182*25c28e83SPiotr Jasiukajtis } 1183*25c28e83SPiotr Jasiukajtis } 1184*25c28e83SPiotr Jasiukajtis } 1185*25c28e83SPiotr Jasiukajtis } 1186*25c28e83SPiotr Jasiukajtis 1187*25c28e83SPiotr Jasiukajtis /* restore the control and status words, check for over/underflow */ 1188*25c28e83SPiotr Jasiukajtis __fenv_setcwsw(&oldcwsw); 1189*25c28e83SPiotr Jasiukajtis if (ez >= 0x7fff) { 1190*25c28e83SPiotr Jasiukajtis if (rm == FCW_RN || rm == FCW_RP) { 1191*25c28e83SPiotr Jasiukajtis zz.i[2] = sz | 0x7fff; 1192*25c28e83SPiotr Jasiukajtis zz.i[1] = 0x80000000; 1193*25c28e83SPiotr Jasiukajtis zz.i[0] = 0; 1194*25c28e83SPiotr Jasiukajtis } else { 1195*25c28e83SPiotr Jasiukajtis zz.i[2] = sz | 0x7ffe; 1196*25c28e83SPiotr Jasiukajtis zz.i[1] = 0xffffffff; 1197*25c28e83SPiotr Jasiukajtis zz.i[0] = 0xffffffff; 1198*25c28e83SPiotr Jasiukajtis } 1199*25c28e83SPiotr Jasiukajtis dummy = huge; 1200*25c28e83SPiotr Jasiukajtis dummy *= huge; 1201*25c28e83SPiotr Jasiukajtis } else { 1202*25c28e83SPiotr Jasiukajtis zz.i[2] = sz | ez; 1203*25c28e83SPiotr Jasiukajtis zz.i[1] = z0; 1204*25c28e83SPiotr Jasiukajtis zz.i[0] = z1; 1205*25c28e83SPiotr Jasiukajtis 1206*25c28e83SPiotr Jasiukajtis /* 1207*25c28e83SPiotr Jasiukajtis * tinyafter => result rounded w/ unbounded range would be tiny, 1208*25c28e83SPiotr Jasiukajtis * z2 nonzero => result delivered is inexact 1209*25c28e83SPiotr Jasiukajtis */ 1210*25c28e83SPiotr Jasiukajtis if (tinyafter) { 1211*25c28e83SPiotr Jasiukajtis dummy = tiny; 1212*25c28e83SPiotr Jasiukajtis if (z2) 1213*25c28e83SPiotr Jasiukajtis dummy *= tiny; 1214*25c28e83SPiotr Jasiukajtis else 1215*25c28e83SPiotr Jasiukajtis dummy -= tiny2; 1216*25c28e83SPiotr Jasiukajtis } else if (z2) { 1217*25c28e83SPiotr Jasiukajtis dummy = huge; 1218*25c28e83SPiotr Jasiukajtis dummy += tiny; 1219*25c28e83SPiotr Jasiukajtis } 1220*25c28e83SPiotr Jasiukajtis } 1221*25c28e83SPiotr Jasiukajtis 1222*25c28e83SPiotr Jasiukajtis return (zz.e); 1223*25c28e83SPiotr Jasiukajtis } 1224*25c28e83SPiotr Jasiukajtis 1225*25c28e83SPiotr Jasiukajtis #else 1226*25c28e83SPiotr Jasiukajtis #error Unknown architecture 1227*25c28e83SPiotr Jasiukajtis #endif 1228