/*
 * CDDL HEADER START
 *
 * The contents of this file are subject to the terms of the
 * Common Development and Distribution License (the "License").
 * You may not use this file except in compliance with the License.
 *
 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
 * or http://www.opensolaris.org/os/licensing.
 * See the License for the specific language governing permissions
 * and limitations under the License.
 *
 * When distributing Covered Code, include this CDDL HEADER in each
 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
 * If applicable, add the following below this CDDL HEADER, with the
 * fields enclosed by brackets "[]" replaced with your own identifying
 * information: Portions Copyright [yyyy] [name of copyright owner]
 *
 * CDDL HEADER END
 */
/*
 * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
 */
/*
 * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
 * Use is subject to license terms.
 */

#pragma weak __fmod = fmod

#include "libm.h"

static const double zero = 0.0;

/*
 * The following implementation assumes fast 64-bit integer arith-
 * metic.  This is fine for sparc because we build libm in v8plus
 * mode.  It's also fine for sparcv9 and amd64, although we have
 * assembly code on amd64.  For x86, it would be better to use
 * 32-bit code, but we have assembly for x86, too.
 */
double
fmod(double x, double y) {
        double          w;
        long long       hx, ix, iy, iz;
        int             nd, k, ny;

        hx = *(long long *)&x;
        ix = hx & ~0x8000000000000000ull;
        iy = *(long long *)&y & ~0x8000000000000000ull;

        /* handle special cases */
        if (iy == 0ll)
                return (_SVID_libm_err(x, y, 27));

        if (ix >= 0x7ff0000000000000ll || iy > 0x7ff0000000000000ll)
                return ((x * y) * zero);

        if (ix <= iy)
                return ((ix < iy)? x : x * zero);

        /*
         * Set:
         *      ny = true exponent of y
         *      nd = true exponent of x minus true exponent of y
         *      ix = normalized significand of x
         *      iy = normalized significand of y
         */
        ny = iy >> 52;
        k = ix >> 52;
        if (ny == 0) {
                /* y is subnormal, x could be normal or subnormal */
                ny = 1;
                while (iy < 0x0010000000000000ll) {
                        ny -= 1;
                        iy += iy;
                }
                nd = k - ny;
                if (k == 0) {
                        nd += 1;
                        while (ix < 0x0010000000000000ll) {
                                nd -= 1;
                                ix += ix;
                        }
                } else {
                        ix = 0x0010000000000000ll | (ix & 0x000fffffffffffffll);
                }
        } else {
                /* both x and y are normal */
                nd = k - ny;
                ix = 0x0010000000000000ll | (ix & 0x000fffffffffffffll);
                iy = 0x0010000000000000ll | (iy & 0x000fffffffffffffll);
        }

        /* perform fixed point mod */
        while (nd--) {
                iz = ix - iy;
                if (iz >= 0)
                        ix = iz;
                ix += ix;
        }
        iz = ix - iy;
        if (iz >= 0)
                ix = iz;

        /* convert back to floating point and restore the sign */
        if (ix == 0ll)
                return (x * zero);
        while (ix < 0x0010000000000000ll) {
                ix += ix;
                ny -= 1;
        }
        while (ix > 0x0020000000000000ll) {     /* XXX can this ever happen? */
                ny += 1;
                ix >>= 1;
        }
        if (ny <= 0) {
                /* result is subnormal */
                k = -ny + 1;
                ix >>= k;
                *(long long *)&w = (hx & 0x8000000000000000ull) | ix;
                return (w);
        }
        *(long long *)&w = (hx & 0x8000000000000000ull) |
            ((long long)ny << 52) | (ix & 0x000fffffffffffffll);
        return (w);
}