25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const union {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned i[2];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x3fe00000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x40000000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x3ef00000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x3e700000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x41300000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x3e300000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x3b300000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x38300000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x42300000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x3df00000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x7fe00000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x00100000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x00100001u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x7ff00000u, 0 },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0x7ff00001u, 0 }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const unsigned int fsr_rm = 0xc0000000u;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fmal for SPARC: 128-bit quad precision, big-endian
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__fmal(long double x, long double y, long double z) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned int i[4];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis long double q;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned int i[2];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned int xy0, xy1, xy2, xy3, xy4, xy5, xy6, xy7;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned int z0, z1, z2, z3, z4, z5, z6, z7;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int hx, hy, hz, ex, ey, ez, exy, sxy, sz, e, ibit;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* extract the high order words of the arguments */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * distinguish zero, finite nonzero, infinite, and quiet nan
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * arguments; raise invalid and return for signaling nans
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((hx & 0xffff) | xx.i[1] | xx.i[2] | xx.i[3]) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* signaling nan, raise invalid */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (hx == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis cx = (xx.i[1] | xx.i[2] | xx.i[3]) ? 1 : 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* subnormal or zero */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((hy & 0xffff) | yy.i[1] | yy.i[2] | yy.i[3]) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (hy == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis cy = (yy.i[1] | yy.i[2] | yy.i[3]) ? 1 : 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((hz & 0xffff) | zz.i[1] | zz.i[2] | zz.i[3]) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (hz == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis cz = (zz.i[1] | zz.i[2] | zz.i[3]) ? 1 : 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* get the fsr and clear current exceptions */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* handle all other zero, inf, and nan cases */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* if x or y is a quiet nan, return it */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* if x*y is 0*inf, raise invalid and return the default nan */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((cx == 0 && cy == 2) || (cx == 2 && cy == 0)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* if z is a quiet nan, return it */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * now none of x, y, or z is nan; handle cases where x or y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * if z is also inf, either we have inf-inf or
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * the result is the same as z depending on signs
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) < 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* otherwise the result is inf with appropriate sign */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis zz.i[0] = ((xx.i[0] ^ yy.i[0]) & 0x80000000) |
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* if z is inf, return it */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * now x, y, and z are all finite; handle cases where x or y
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* either we have 0-0 or the result is the same as z */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (cz == 0 && (int) ((xx.i[0] ^ yy.i[0]) ^ zz.i[0]) <
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis zz.i[0] = (fsr >> 30) == FSR_RM ? 0x80000000 :
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* if we get here, x and y are nonzero finite, z must be zero */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (x * y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * now x, y, and z are all finite and nonzero; set round-to-
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * negative-infinity mode
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * get the signs and exponents and normalize the significands
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (xx.i[1] | (xx.i[2] & 0xfffe0000)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (xx.i[2] | (xx.i[3] & 0xfffe0000)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xx.i[1] = (xx.i[1] << 1) | (xx.i[2] >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xx.i[2] = (xx.i[2] << 1) | (xx.i[3] >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (yy.i[1] | (yy.i[2] & 0xfffe0000)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (yy.i[2] | (yy.i[3] & 0xfffe0000)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis yy.i[1] = (yy.i[1] << 1) | (yy.i[2] >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis yy.i[2] = (yy.i[2] << 1) | (yy.i[3] >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* convert the significands of x and y to doubles */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dx[1] = (double) ((int) (xx.i[1] >> 8)) * c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dy[1] = (double) ((int) (yy.i[1] >> 8)) * c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dx[2] = (double) ((int) (((xx.i[1] << 16) | (xx.i[2] >> 16)) &
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 0xffffff)) * c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dy[2] = (double) ((int) (((yy.i[1] << 16) | (yy.i[2] >> 16)) &
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 0xffffff)) * c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dx[3] = (double) ((int) (((xx.i[2] << 8) | (xx.i[3] >> 24)) &
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 0xffffff)) * c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dy[3] = (double) ((int) (((yy.i[2] << 8) | (yy.i[3] >> 24)) &
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis 0xffffff)) * c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dx[4] = (double) ((int) (xx.i[3] & 0xffffff)) * c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dy[4] = (double) ((int) (yy.i[3] & 0xffffff)) * c;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* form the "digits" of the product */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dxy[2] = dx[0] * dy[2] + dx[1] * dy[1] + dx[2] * dy[0];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dxy[3] = dx[0] * dy[3] + dx[1] * dy[2] + dx[2] * dy[1] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dxy[4] = dx[0] * dy[4] + dx[1] * dy[3] + dx[2] * dy[2] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dxy[5] = dx[1] * dy[4] + dx[2] * dy[3] + dx[3] * dy[2] +
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dxy[6] = dx[2] * dy[4] + dx[3] * dy[3] + dx[4] * dy[2];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* split odd-numbered terms and combine into even-numbered terms */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* propagate carries, adjusting the exponent if need be */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* extract the significand of x*y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* extract the sign, exponent, and significand of z */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (zz.i[1] | (zz.i[2] & 0xfffe0000)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (zz.i[2] | (zz.i[3] & 0xfffe0000)) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * now x*y is represented by sxy, exy, and xy[0-7], and z is
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * represented likewise; swap if need be so |xy| <= |z|
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 && (xy1 > z1 ||
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis (xy1 == z1 && (xy2 > z2 || (xy2 == z2 && (xy3 > z3 ||
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis (xy3 == z3 && (xy4 | xy5 | xy6 | xy7) != 0)))))))))) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* shift the significand of xy keeping a sticky bit */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (e > 236) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 224) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 | xy1 |
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy0 = xy1 = xy2 = xy3 = xy4 = xy5 = xy6 = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 192) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = xy7 | xy6 | xy5 | xy4 | xy3 | xy2 |
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy7 = (xy1 >> (e - 192)) | ((xy0 << 1) << (223 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 160) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy7 = (xy2 >> (e - 160)) | ((xy1 << 1) << (191 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy6 = (xy1 >> (e - 160)) | ((xy0 << 1) << (191 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 128) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = xy7 | xy6 | xy5 | xy4 | ((xy3 << 1) << (159 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy7 = (xy3 >> (e - 128)) | ((xy2 << 1) << (159 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy6 = (xy2 >> (e - 128)) | ((xy1 << 1) << (159 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy5 = (xy1 >> (e - 128)) | ((xy0 << 1) << (159 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 96) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = xy7 | xy6 | xy5 | ((xy4 << 1) << (127 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy7 = (xy4 >> (e - 96)) | ((xy3 << 1) << (127 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy6 = (xy3 >> (e - 96)) | ((xy2 << 1) << (127 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy5 = (xy2 >> (e - 96)) | ((xy1 << 1) << (127 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 64) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = xy7 | xy6 | ((xy5 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy7 = (xy5 >> (e - 64)) | ((xy4 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy6 = (xy4 >> (e - 64)) | ((xy3 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy5 = (xy3 >> (e - 64)) | ((xy2 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 32) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy7 = (xy6 >> (e - 32)) | ((xy5 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy6 = (xy5 >> (e - 32)) | ((xy4 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy5 = (xy4 >> (e - 32)) | ((xy3 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy7 = (xy7 >> e) | ((xy6 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy6 = (xy6 >> e) | ((xy5 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy5 = (xy5 >> e) | ((xy4 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy4 = (xy4 >> e) | ((xy3 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* if this is a magnitude subtract, negate the significand of xy */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* add, propagating carries */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* postnormalize and collect rounding information into z4 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* result is tiny; shift right until exponent is within range */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (e > 116) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z4 = 1; /* result can't be exactly zero */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 96) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z4 = (z1 >> (e - 96)) | ((z0 << 1) << (127 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 64) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z4 = (z2 >> (e - 64)) | ((z1 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z3 = (z1 >> (e - 64)) | ((z0 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 32) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = z7 | z6 | z5 | z4 | ((z3 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z4 = (z3 >> (e - 32)) | ((z2 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z3 = (z2 >> (e - 32)) | ((z1 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = z7 | z6 | z5 | (z4 << 1) << (31 - e);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* carry out; shift right by one */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (z0 < 0x10000 && (z0 | z1 | z2 | z3 | z4 | z5 | z6 | z7)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * borrow/cancellation; shift left as much as
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * exponent allows
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (!(z0 | (z1 & 0xfffe0000)) && ez >= 33) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* get the rounding mode */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* strip off the integer bit, if there is one */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (!(z0 | z1 | z2 | z3 | z4)) { /* exact zero */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * flip the sense of directed roundings if the result is negative;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * the logic below applies to a positive result
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* round and raise exceptions */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* decide whether to round the fraction up */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (rm == FSR_RP || (rm == FSR_RN && (z4 > 0x80000000u ||
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* round up and renormalize if necessary */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* check for under/overflow */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * !ibit => exact result was tiny before rounding,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * z4 nonzero => result delivered is inexact
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* restore the fsr and emulate exceptions as needed */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const union {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned i[2];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x3fe00000u },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x40000000u },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x3df00000u },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x3bf00000u },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x41f00000u },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x43e00000u },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x7fe00000u },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x00100000u },
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis { 0, 0x00100001u }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fmal for x86: 80-bit extended double precision, little-endian
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis__fmal(long double x, long double y, long double z) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis long double e;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned xy0, xy1, xy2, xy3, xy4, z0, z1, z2, z3, z4;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned oldcwsw, cwsw, rm, sticky, carry;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int ex, ey, ez, exy, sxy, sz, e, tinyafter;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* extract the exponents of the arguments */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* dispense with inf, nan, and zero cases */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ex == 0x7fff || ey == 0x7fff || (ex | xx.i[1] | xx.i[0]) == 0 ||
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis (ey | yy.i[1] | yy.i[0]) == 0) /* x or y is inf, nan, or 0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (x * y + z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (x + z); /* avoid spurious under/overflow in x * y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((ez | zz.i[1] | zz.i[0]) == 0) /* z is zero */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x * y isn't zero but could underflow to zero,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * so don't add z, lest we perturb the sign
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (x * y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * now x, y, and z are all finite and nonzero; extract signs and
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * normalize the significands (this will raise the denormal operand
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * exception if need be)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * save the control and status words, mask all exceptions, and
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * set rounding to 64-bit precision and toward-zero
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis cwsw = (oldcwsw & 0xf0c0ffff) | 0x0f3f0000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* multiply x*y to 128 bits */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y = ((xhi * yhi - x) + xhi * ylo + xlo * yhi) + xlo * ylo;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* extract the significands */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * now x*y is represented by sxy, exy, and xy[0-4], and z is
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * represented likewise; swap if need be so |xy| <= |z|
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (exy > ez || (exy == ez && (xy0 > z0 || (xy0 == z0 &&
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis (xy1 > z1 || (xy1 == z1 && (xy2 | xy3) != 0)))))) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* shift the significand of xy keeping a sticky bit */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (e > 130) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 128) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = xy3 | xy2 | xy1 | ((xy0 << 1) << (159 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 96) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = xy3 | xy2 | ((xy1 << 1) << (127 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy4 = (xy1 >> (e - 96)) | ((xy0 << 1) << (127 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 64) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy4 = (xy2 >> (e - 64)) | ((xy1 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy3 = (xy1 >> (e - 64)) | ((xy0 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 32) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy4 = (xy3 >> (e - 32)) | ((xy2 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy3 = (xy2 >> (e - 32)) | ((xy1 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy2 = (xy1 >> (e - 32)) | ((xy0 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy3 = (xy3 >> e) | ((xy2 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy2 = (xy2 >> e) | ((xy1 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis xy1 = (xy1 >> e) | ((xy0 << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* if this is a magnitude subtract, negate the significand of xy */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* add, propagating carries */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* for a magnitude subtract, ignore the last carry out */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* postnormalize and collect rounding information into z2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* result is tiny; shift right until exponent is within range */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = 1; /* result can't be exactly zero */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 64) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = z4 | z3 | z2 | z1 | ((z0 << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = (z0 >> (e - 64)) | ((carry << 1) << (95 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (e >= 32) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sticky = z4 | z3 | z2 | ((z1 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = (z1 >> (e - 32)) | ((z0 << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z1 = (z0 >> (e - 32)) | ((carry << 1) << (63 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z0 = (z0 >> e) | ((carry << 1) << (31 - e));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* carry out; shift right by one */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (z0 < 0x80000000u && (z0 | z1 | z2 | z3 | z4) != 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * borrow/cancellation; shift left as much as
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * exponent allows
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* get the rounding mode */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* adjust exponent if result is subnormal */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * flip the sense of directed roundings if the result is negative;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * the logic below applies to a positive result
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (rm == FCW_RP || (rm == FCW_RN && (z2 > 0x80000000u ||
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* round up and renormalize if necessary */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (++z1 == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (++z0 == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* rounded up to smallest normal */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * would have rounded up to
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * smallest normal even with
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * unbounded range
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* restore the control and status words, check for over/underflow */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * tinyafter => result rounded w/ unbounded range would be tiny,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * z2 nonzero => result delivered is inexact
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (z2) {