/*
 * CDDL HEADER START
 *
 * The contents of this file are subject to the terms of the
 * Common Development and Distribution License, Version 1.0 only
 * (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 (c) 1994-1997, by Sun Microsystems, Inc.
 * All rights reserved.
 */

/*
 * This file contains __quad_mag_add and __quad_mag_sub, the core
 * of the quad precision add and subtract operations.
 */

#include "quad.h"

/*
 * __quad_mag_add(x, y, z, fsr)
 *
 * Sets *z = *x + *y, rounded according to the rounding mode in *fsr,
 * and updates the current exceptions in *fsr.  This routine assumes
 * *x and *y are finite, with the same sign (i.e., an addition of
 * magnitudes), |*x| >= |*y|, and *z already has its sign bit set.
 */
void
__quad_mag_add(const union longdouble *x, const union longdouble *y,
        union longdouble *z, unsigned int *fsr)
{
        unsigned int    lx, ly, ex, ey, frac2, frac3, frac4;
        unsigned int    round, sticky, carry, rm;
        int             e, uflo;

        /* get the leading significand words and exponents */
        ex = (x->l.msw & 0x7fffffff) >> 16;
        lx = x->l.msw & 0xffff;
        if (ex == 0)
                ex = 1;
        else
                lx |= 0x10000;

        ey = (y->l.msw & 0x7fffffff) >> 16;
        ly = y->l.msw & 0xffff;
        if (ey == 0)
                ey = 1;
        else
                ly |= 0x10000;

        /* prenormalize y */
        e = (int) ex - (int) ey;
        round = sticky = 0;
        if (e >= 114) {
                frac2 = x->l.frac2;
                frac3 = x->l.frac3;
                frac4 = x->l.frac4;
                sticky = ly | y->l.frac2 | y->l.frac3 | y->l.frac4;
        } else {
                frac2 = y->l.frac2;
                frac3 = y->l.frac3;
                frac4 = y->l.frac4;
                if (e >= 96) {
                        sticky = frac4 | frac3 | (frac2 & 0x7fffffff);
                        round = frac2 & 0x80000000;
                        frac4 = ly;
                        frac3 = frac2 = ly = 0;
                        e -= 96;
                } else if (e >= 64) {
                        sticky = frac4 | (frac3 & 0x7fffffff);
                        round = frac3 & 0x80000000;
                        frac4 = frac2;
                        frac3 = ly;
                        frac2 = ly = 0;
                        e -= 64;
                } else if (e >= 32) {
                        sticky = frac4 & 0x7fffffff;
                        round = frac4 & 0x80000000;
                        frac4 = frac3;
                        frac3 = frac2;
                        frac2 = ly;
                        ly = 0;
                        e -= 32;
                }
                if (e) {
                        sticky |= round | (frac4 & ((1 << (e - 1)) - 1));
                        round = frac4 & (1 << (e - 1));
                        frac4 = (frac4 >> e) | (frac3 << (32 - e));
                        frac3 = (frac3 >> e) | (frac2 << (32 - e));
                        frac2 = (frac2 >> e) | (ly << (32 - e));
                        ly >>= e;
                }

                /* add, propagating carries */
                frac4 += x->l.frac4;
                carry = (frac4 < x->l.frac4);
                frac3 += x->l.frac3;
                if (carry) {
                        frac3++;
                        carry = (frac3 <= x->l.frac3);
                } else {
                        carry = (frac3 < x->l.frac3);
                }
                frac2 += x->l.frac2;
                if (carry) {
                        frac2++;
                        carry = (frac2 <= x->l.frac2);
                } else {
                        carry = (frac2 < x->l.frac2);
                }
                lx += ly;
                if (carry)
                        lx++;

                /* postnormalize */
                if (lx >= 0x20000) {
                        sticky |= round;
                        round = frac4 & 1;
                        frac4 = (frac4 >> 1) | (frac3 << 31);
                        frac3 = (frac3 >> 1) | (frac2 << 31);
                        frac2 = (frac2 >> 1) | (lx << 31);
                        lx >>= 1;
                        ex++;
                }
        }

        /* keep track of whether the result before rounding is tiny */
        uflo = (lx < 0x10000);

        /* get the rounding mode, fudging directed rounding modes */
        /* as though the result were positive */
        rm = *fsr >> 30;
        if (z->l.msw)
                rm ^= (rm >> 1);

        /* see if we need to round */
        if (round | sticky) {
                *fsr |= FSR_NXC;

                /* round up if necessary */
                if (rm == FSR_RP || (rm == FSR_RN && round &&
                        (sticky || (frac4 & 1)))) {
                        if (++frac4 == 0)
                                if (++frac3 == 0)
                                        if (++frac2 == 0)
                                                if (++lx >= 0x20000) {
                                                        lx >>= 1;
                                                        ex++;
                                                }
                }
        }

        /* check for overflow */
        if (ex >= 0x7fff) {
                /* store the default overflowed result */
                *fsr |= FSR_OFC | FSR_NXC;
                if (rm == FSR_RN || rm == FSR_RP) {
                        z->l.msw |= 0x7fff0000;
                        z->l.frac2 = z->l.frac3 = z->l.frac4 = 0;
                } else {
                        z->l.msw |= 0x7ffeffff;
                        z->l.frac2 = z->l.frac3 = z->l.frac4 = 0xffffffff;
                }
        } else {
                /* store the result */
                if (lx >= 0x10000)
                        z->l.msw |= (ex << 16);
                z->l.msw |= (lx & 0xffff);
                z->l.frac2 = frac2;
                z->l.frac3 = frac3;
                z->l.frac4 = frac4;

                /* if the pre-rounded result was tiny and underflow trapping */
                /* is enabled, simulate underflow */
                if (uflo && (*fsr & FSR_UFM))
                        *fsr |= FSR_UFC;
        }
}

/*
 * __quad_mag_sub(x, y, z, fsr)
 *
 * Sets *z = *x - *y, rounded according to the rounding mode in *fsr,
 * and updates the current exceptions in *fsr.  This routine assumes
 * *x and *y are finite, with opposite signs (i.e., a subtraction of
 * magnitudes), |*x| >= |*y|, and *z already has its sign bit set.
 */
void
__quad_mag_sub(const union longdouble *x, const union longdouble *y,
        union longdouble *z, unsigned int *fsr)
{
        unsigned int    lx, ly, ex, ey, frac2, frac3, frac4;
        unsigned int    guard, round, sticky, borrow, rm;
        int             e;

        /* get the leading significand words and exponents */
        ex = (x->l.msw & 0x7fffffff) >> 16;
        lx = x->l.msw & 0xffff;
        if (ex == 0)
                ex = 1;
        else
                lx |= 0x10000;

        ey = (y->l.msw & 0x7fffffff) >> 16;
        ly = y->l.msw & 0xffff;
        if (ey == 0)
                ey = 1;
        else
                ly |= 0x10000;

        /* prenormalize y */
        e = (int) ex - (int) ey;
        guard = round = sticky = 0;
        if (e > 114) {
                sticky = ly | y->l.frac2 | y->l.frac3 | y->l.frac4;
                ly = frac2 = frac3 = frac4 = 0;
        } else {
                frac2 = y->l.frac2;
                frac3 = y->l.frac3;
                frac4 = y->l.frac4;
                if (e >= 96) {
                        sticky = frac4 | frac3 | (frac2 & 0x3fffffff);
                        round = frac2 & 0x40000000;
                        guard = frac2 & 0x80000000;
                        frac4 = ly;
                        frac3 = frac2 = ly = 0;
                        e -= 96;
                } else if (e >= 64) {
                        sticky = frac4 | (frac3 & 0x3fffffff);
                        round = frac3 & 0x40000000;
                        guard = frac3 & 0x80000000;
                        frac4 = frac2;
                        frac3 = ly;
                        frac2 = ly = 0;
                        e -= 64;
                } else if (e >= 32) {
                        sticky = frac4 & 0x3fffffff;
                        round = frac4 & 0x40000000;
                        guard = frac4 & 0x80000000;
                        frac4 = frac3;
                        frac3 = frac2;
                        frac2 = ly;
                        ly = 0;
                        e -= 32;
                }
                if (e > 1) {
                        sticky |= guard | round |
                                (frac4 & ((1 << (e - 2)) - 1));
                        round = frac4 & (1 << (e - 2));
                        guard = frac4 & (1 << (e - 1));
                        frac4 = (frac4 >> e) | (frac3 << (32 - e));
                        frac3 = (frac3 >> e) | (frac2 << (32 - e));
                        frac2 = (frac2 >> e) | (ly << (32 - e));
                        ly >>= e;
                } else if (e == 1) {
                        sticky |= round;
                        round = guard;
                        guard = frac4 & 1;
                        frac4 = (frac4 >> 1) | (frac3 << 31);
                        frac3 = (frac3 >> 1) | (frac2 << 31);
                        frac2 = (frac2 >> 1) | (ly << 31);
                        ly >>= 1;
                }
        }

        /* complement guard, round, and sticky as need be */
        if (sticky) {
                round = !round;
                guard = !guard;
        } else if (round) {
                guard = !guard;
        }
        borrow = (guard | round | sticky);

        /* subtract, propagating borrows */
        frac4 = x->l.frac4 - frac4;
        if (borrow) {
                frac4--;
                borrow = (frac4 >= x->l.frac4);
        } else {
                borrow = (frac4 > x->l.frac4);
        }
        frac3 = x->l.frac3 - frac3;
        if (borrow) {
                frac3--;
                borrow = (frac3 >= x->l.frac3);
        } else {
                borrow = (frac3 > x->l.frac3);
        }
        frac2 = x->l.frac2 - frac2;
        if (borrow) {
                frac2--;
                borrow = (frac2 >= x->l.frac2);
        } else {
                borrow = (frac2 > x->l.frac2);
        }
        lx -= ly;
        if (borrow)
                lx--;

        /* get the rounding mode */
        rm = *fsr >> 30;

        /* handle zero result */
        if (!(lx | frac2 | frac3 | frac4 | guard)) {
                z->l.msw = ((rm == FSR_RM)? 0x80000000 : 0);
                z->l.frac2 = z->l.frac3 = z->l.frac4 = 0;
                return;
        }

        /* postnormalize */
        if (lx < 0x10000) {
                /* if cancellation occurred or the exponent is 1, */
                /* the result is exact */
                if (lx < 0x8000 || ex == 1) {
                        while ((lx | (frac2 & 0xfffe0000)) == 0 && ex > 32) {
                                lx = frac2;
                                frac2 = frac3;
                                frac3 = frac4;
                                frac4 = ((guard)? 0x80000000 : 0);
                                guard = 0;
                                ex -= 32;
                        }
                        while (lx < 0x10000 && ex > 1) {
                                lx = (lx << 1) | (frac2 >> 31);
                                frac2 = (frac2 << 1) | (frac3 >> 31);
                                frac3 = (frac3 << 1) | (frac4 >> 31);
                                frac4 <<= 1;
                                if (guard) {
                                        frac4 |= 1;
                                        guard = 0;
                                }
                                ex--;
                        }
                        if (lx >= 0x10000)
                                z->l.msw |= (ex << 16);
                        z->l.msw |= (lx & 0xffff);
                        z->l.frac2 = frac2;
                        z->l.frac3 = frac3;
                        z->l.frac4 = frac4;

                        /* if the result is tiny and underflow trapping is */
                        /* enabled, simulate underflow */
                        if (lx < 0x10000 && (*fsr & FSR_UFM))
                                *fsr |= FSR_UFC;
                        return;
                }

                /* otherwise we only borrowed one place */
                lx = (lx << 1) | (frac2 >> 31);
                frac2 = (frac2 << 1) | (frac3 >> 31);
                frac3 = (frac3 << 1) | (frac4 >> 31);
                frac4 <<= 1;
                if (guard)
                        frac4 |= 1;
                ex--;
        } else {
                sticky |= round;
                round = guard;
        }

        /* fudge directed rounding modes as though the result were positive */
        if (z->l.msw)
                rm ^= (rm >> 1);

        /* see if we need to round */
        if (round | sticky) {
                *fsr |= FSR_NXC;

                /* round up if necessary */
                if (rm == FSR_RP || (rm == FSR_RN && round &&
                        (sticky || (frac4 & 1)))) {
                        if (++frac4 == 0)
                                if (++frac3 == 0)
                                        if (++frac2 == 0)
                                                if (++lx >= 0x20000) {
                                                        lx >>= 1;
                                                        ex++;
                                                }
                }
        }

        /* store the result */
        z->l.msw |= (ex << 16) | (lx & 0xffff);
        z->l.frac2 = frac2;
        z->l.frac3 = frac3;
        z->l.frac4 = frac4;
}