common/m9x/remquo.c

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 __remquo = remquo
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/*
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * double remquo(double x, double y, int *quo) return remainder(x,y) and an
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * integer pointer quo such that *quo = N mod {2**31}, where N is the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * exact integral part of x/y rounded to nearest even.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis *
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * remquo call internal fmodquo
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm.h"
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm_protos.h"
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include <math.h>       /* fabs() */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include <sys/isa_defs.h>
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#if defined(_BIG_ENDIAN)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define HIWORD  0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define LOWORD  1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define HIWORD  1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define LOWORD  0
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#endif
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define __HI(x) ((int *) &x)[HIWORD]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define __LO(x) ((int *) &x)[LOWORD]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double one = 1.0, Zero[] = {0.0, -0.0};
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisfmodquo(double x, double y, int *quo) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int n, hx, hy, hz, ix, iy, sx, sq, i, m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    unsigned lx, ly, lz;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hx = __HI(x);       /* high word of x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    lx = __LO(x);       /* low  word of x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hy = __HI(y);       /* high word of y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ly = __LO(y);       /* low  word of y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    sx = hx & 0x80000000;   /* sign of x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    sq = (hx ^ hy) & 0x80000000;    /* sign of x/y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hx ^= sx;       /* |x| */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hy &= 0x7fffffff;   /* |y| */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    /* purge off exception values */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    *quo = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if ((hy | ly) == 0 || hx >= 0x7ff00000 ||   /* y=0, or x !finite */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return ((x * y) / (x * y));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (hx <= hy) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (hx < hy || lx < ly)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (x); /* |x|<|y| return x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (lx == ly) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            *quo = 1 + (sq >> 30);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            /* |x|=|y| return x*0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            return (Zero[(unsigned) sx >> 31]);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    /* determine ix = ilogb(x) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (hx < 0x00100000) {  /* subnormal x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (hx == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            for (ix = -1043, i = lx; i > 0; i <<= 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                ix -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                ix -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    } else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        ix = (hx >> 20) - 1023;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    /* determine iy = ilogb(y) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (hy < 0x00100000) {  /* subnormal y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (hy == 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            for (iy = -1043, i = ly; i > 0; i <<= 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                iy -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                iy -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    } else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        iy = (hy >> 20) - 1023;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    /* set up {hx,lx}, {hy,ly} and align y to x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (ix >= -1022)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        hx = 0x00100000 | (0x000fffff & hx);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    else {          /* subnormal x, shift x to normal */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        n = -1022 - ix;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (n <= 31) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hx = (hx << n) | (lx >> (32 - n));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            lx <<= n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hx = lx << (n - 32);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            lx = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (iy >= -1022)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        hy = 0x00100000 | (0x000fffff & hy);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    else {          /* subnormal y, shift y to normal */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        n = -1022 - iy;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (n <= 31) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hy = (hy << n) | (ly >> (32 - n));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ly <<= n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hy = ly << (n - 32);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            ly = 0;
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        hz = hx - hy;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        lz = lx - ly;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (lx < ly)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hz -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (hz < 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hx = hx + hx + (lx >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            lx = lx + lx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            m += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            if ((hz | lz) == 0) {   /* return sign(x)*0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                if (n < 31)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                    m <<= 1 + n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                else
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                    m = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                m &= 0x7fffffff;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                *quo = sq >= 0 ? m : -m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis                return (Zero[(unsigned) sx >> 31]);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hx = hz + hz + (lz >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            lx = lz + lz;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        m += m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hz = hx - hy;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    lz = lx - ly;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (lx < ly)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        hz -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (hz >= 0) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        hx = hz;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        lx = lz;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        m += 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    m &= 0x7fffffff;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    *quo = sq >= 0 ? m : -m;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    /* convert back to floating value and restore the sign */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if ((hx | lx) == 0) {   /* return sign(x)*0 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return (Zero[(unsigned) sx >> 31]);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    while (hx < 0x00100000) {   /* normalize x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        hx = hx + hx + (lx >> 31);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        lx = lx + lx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        iy -= 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (iy >= -1022) {  /* normalize output */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        hx = (hx - 0x00100000) | ((iy + 1023) << 20);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        __HI(x) = hx | sx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        __LO(x) = lx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    } else {            /* subnormal output */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        n = -1022 - iy;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        if (n <= 20) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            lx = (lx >> n) | ((unsigned) hx << (32 - n));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hx >>= n;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else if (n <= 31) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            lx = (hx << (32 - n)) | (lx >> n);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hx = sx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        } else {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            lx = hx >> (n - 32);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis            hx = sx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        __HI(x) = hx | sx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        __LO(x) = lx;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        x *= one;   /* create necessary signal */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    return (x);     /* exact output */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis}
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#define zero    Zero[0]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisdouble
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisremquo(double x, double y, int *quo) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    int hx, hy, sx, sq;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    double v;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    unsigned ly;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hx = __HI(x);       /* high word of x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hy = __HI(y);       /* high word of y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    ly = __LO(y);       /* low  word of y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    sx = hx & 0x80000000;   /* sign of x */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    sq = (hx ^ hy) & 0x80000000;    /* sign of x/y */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hx ^= sx;       /* |x| */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    hy &= 0x7fffffff;   /* |y| */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    /* purge off exception values */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    *quo = 0;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if ((hy | ly) == 0 || hx >= 0x7ff00000 ||   /* y=0, or x !finite */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        (hy | ((ly | -ly) >> 31)) > 0x7ff00000) /* or y is NaN */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        return ((x * y) / (x * y));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    y = fabs(y);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    x = fabs(x);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (hy <= 0x7fdfffff) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        x = fmodquo(x, y + y, quo);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis        *quo = ((*quo) & 0x3fffffff) << 1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    }
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis    if (hy < 0x00200000) {
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.5 * 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}