25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
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 *
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 *
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 *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
ddc0e0b53c661f6e439e3b7072b3ef353eadb4afRichard Lowe#pragma weak __remquol = remquol
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm.h"
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__SUNPRO_C)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include <sunmath.h> /* fabsl */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const int
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis is = -0x7fffffff - 1,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis im = 0x0000ffff,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis iu = 0x00010000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const long double zero = 0.0L, one = 1.0L;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(__sparc)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define __H0(x) ((int *) &x)[0]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define __H1(x) ((int *) &x)[1]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define __H2(x) ((int *) &x)[2]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define __H3(x) ((int *) &x)[3]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#error Unsupported architecture
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * On entrance: *quo is initialized to 0, x finite and y non-zero & ordered
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic long double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisfmodquol(long double x, long double y, int *quo) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis long double a, b;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int n, ix, iy, k, sx, sq, m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int hx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int x0, y0, z0, carry;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis hx = __H0(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = __H1(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = __H2(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = __H3(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y0 = __H0(y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y1 = __H1(y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y2 = __H2(y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y3 = __H3(y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sx = hx & is;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sq = (hx ^ y0) & is;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = hx ^ sx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y0 &= ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis a = fabsl(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis b = fabsl(y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (a <= b) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (a < b)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo = 1 + (sq >> 30);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (zero * x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* determine ix = ilogbl(x) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x0 < iu) { /* subnormal x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ix = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ix = -16382;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (x0 == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ix -= 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = x1 >> 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = (x1 << 16) | (x2 >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = (x2 << 16) | (x3 >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = (x3 << 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (x0 < iu) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ix -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = (x0 << 1) | (x1 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = (x1 << 1) | (x2 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = (x2 << 1) | (x3 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 <<= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ix = (x0 >> 16) - 16383;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = iu | (x0 & im);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* determine iy = ilogbl(y) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (y0 < iu) { /* subnormal y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis iy = -16382;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (y0 == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis iy -= 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y0 = y1 >> 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y1 = (y1 << 16) | (y2 >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y2 = (y2 << 16) | (y3 >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y3 = (y3 << 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (y0 < iu) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis iy -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y0 = (y0 << 1) | (y1 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y1 = (y1 << 1) | (y2 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y2 = (y2 << 1) | (y3 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y3 <<= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis iy = (y0 >> 16) - 16383;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y0 = iu | (y0 & im);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* fix point fmod */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n = ix - iy;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (n--) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (x0 == 0 && n >= 16) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m <<= 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n -= 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = x1 >> 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = (x1 << 16) | (x2 >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = (x2 << 16) | (x3 >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = (x3 << 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (x0 < iu && n >= 1) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m += m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = (x0 << 1) | (x1 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = (x1 << 1) | (x2 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = (x2 << 1) | (x3 >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = (x3 << 1);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z3 = x3 - y3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z3 > x3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (carry == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = x2 - y2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z2 > x2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = x2 - y2 - 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z2 >= x2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (carry == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z1 = x1 - y1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z1 > x1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z1 = x1 - y1 - 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z1 >= x1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z0 = x0 - y0 - carry;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (z0 < 0) { /* double x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = x0 + x0 + ((x1 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = x1 + x1 + ((x2 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = x2 + x2 + ((x3 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = x3 + x3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m += m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (z0 == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((z1 | z2 | z3) == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* 0: we are done */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (n < 31)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m <<= (1 + n);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m &= ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo = sq >= 0 ? m : -m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H0(a) = hx & is;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H1(a) = __H2(a) = __H3(a) = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (a);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* x = z << 1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z0 = z0 + z0 + ((z1 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z1 = z1 + z1 + ((z2 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = z2 + z2 + ((z3 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z3 = z3 + z3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = z0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = z1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = z2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = z3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m += m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z3 = x3 - y3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z3 > x3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (carry == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = x2 - y2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z2 > x2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z2 = x2 - y2 - 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z2 >= x2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (carry == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z1 = x1 - y1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z1 > x1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z1 = x1 - y1 - 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis carry = z1 >= x1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis z0 = x0 - y0 - carry;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (z0 >= 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = z0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = z1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = z2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = z3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis m &= ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo = sq >= 0 ? m : -m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* convert back to floating value and restore the sign */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if ((x0 | x1 | x2 | x3) == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H0(a) = hx & is;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H1(a) = __H2(a) = __H3(a) = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (a);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (x0 < iu) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x0 == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis iy -= 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = x1 >> 16;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = (x1 << 16) | (x2 >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = (x2 << 16) | (x3 >> 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = (x3 << 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = x0 + x0 + ((x1 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = x1 + x1 + ((x2 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = x2 + x2 + ((x3 & is) != 0);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = x3 + x3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis iy -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* normalize output */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (iy >= -16382) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H1(a) = x1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H2(a) = x2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H3(a) = x3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else { /* subnormal output */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n = -16382 - iy;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis k = n & 31;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (k <= 16) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = (x2 << (32 - k)) | (x3 >> k);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = (x1 << (32 - k)) | (x2 >> k);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = (x0 << (32 - k)) | (x1 >> k);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 >>= k;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = (x2 << (32 - k)) | (x3 >> k);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = (x1 << (32 - k)) | (x2 >> k);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = (x0 << (32 - k)) | (x1 >> k);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis while (n >= 32) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis n -= 32;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x3 = x2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = x1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x1 = x0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x0 = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H0(a) = x0 | sx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H1(a) = x1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H2(a) = x2;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis __H3(a) = x3;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis a *= one;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (a);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtislong double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisremquol(long double x, long double y, int *quo) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int hx, hy, sx, sq;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis long double v;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis hx = __H0(x); /* high word of x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis hy = __H0(y); /* high word of y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sx = hx & is; /* sign of x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sq = (hx ^ hy) & is; /* sign of x/y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis hx ^= sx; /* |x| */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis hy &= ~0x80000000;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* purge off exception values */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* y=0, y is NaN, x is NaN or inf */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (y == 0.0L || y != y || hx >= 0x7fff0000)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return ((x * y) / (x * y));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y = fabsl(y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x = fabsl(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (hy <= 0x7ffdffff) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x = fmodquol(x, y + y, quo);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo = ((*quo) & 0x3fffffff) << 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (hy < 0x00020000) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x + x > y) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x == y)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x = zero;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x -= y;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x + x >= y) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x -= y;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis v = 0.5L * y;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x > v) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x == y)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x = zero;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x -= y;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (x >= v) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x -= y;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (sq != 0)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *quo = -(*quo);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis return (sx == 0 ? x : -x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}