// SPDX-License-Identifier: GPL-2.0-only
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.
 */

#include "ieee754dp.h"

union ieee754dp ieee754dp_mul(union ieee754dp x, union ieee754dp y)
{
        int re;
        int rs;
        u64 rm;
        unsigned int lxm;
        unsigned int hxm;
        unsigned int lym;
        unsigned int hym;
        u64 lrm;
        u64 hrm;
        u64 t;
        u64 at;

        COMPXDP;
        COMPYDP;

        EXPLODEXDP;
        EXPLODEYDP;

        ieee754_clearcx();

        FLUSHXDP;
        FLUSHYDP;

        switch (CLPAIR(xc, yc)) {
        case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
        case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
        case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
        case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
        case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
                return ieee754dp_nanxcpt(y);

        case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
        case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
        case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
        case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
        case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
        case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
                return ieee754dp_nanxcpt(x);

        case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
        case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
        case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
        case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
                return y;

        case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
        case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
        case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
        case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
        case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
                return x;


        /*
         * Infinity handling
         */
        case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
        case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
                ieee754_setcx(IEEE754_INVALID_OPERATION);
                return ieee754dp_indef();

        case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
        case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
        case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
        case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
        case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
                return ieee754dp_inf(xs ^ ys);

        case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
        case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
        case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
        case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
        case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
                return ieee754dp_zero(xs ^ ys);


        case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
                DPDNORMX;
                fallthrough;
        case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
                DPDNORMY;
                break;

        case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
                DPDNORMX;
                break;

        case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
                break;
        }
        /* rm = xm * ym, re = xe+ye basically */
        assert(xm & DP_HIDDEN_BIT);
        assert(ym & DP_HIDDEN_BIT);

        re = xe + ye;
        rs = xs ^ ys;

        /* shunt to top of word */
        xm <<= 64 - (DP_FBITS + 1);
        ym <<= 64 - (DP_FBITS + 1);

        /*
         * Multiply 64 bits xm, ym to give high 64 bits rm with stickness.
         */

        lxm = xm;
        hxm = xm >> 32;
        lym = ym;
        hym = ym >> 32;

        lrm = DPXMULT(lxm, lym);
        hrm = DPXMULT(hxm, hym);

        t = DPXMULT(lxm, hym);

        at = lrm + (t << 32);
        hrm += at < lrm;
        lrm = at;

        hrm = hrm + (t >> 32);

        t = DPXMULT(hxm, lym);

        at = lrm + (t << 32);
        hrm += at < lrm;
        lrm = at;

        hrm = hrm + (t >> 32);

        rm = hrm | (lrm != 0);

        /*
         * Sticky shift down to normal rounding precision.
         */
        if ((s64) rm < 0) {
                rm = (rm >> (64 - (DP_FBITS + 1 + 3))) |
                     ((rm << (DP_FBITS + 1 + 3)) != 0);
                re++;
        } else {
                rm = (rm >> (64 - (DP_FBITS + 1 + 3 + 1))) |
                     ((rm << (DP_FBITS + 1 + 3 + 1)) != 0);
        }
        assert(rm & (DP_HIDDEN_BIT << 3));

        return ieee754dp_format(rs, re, rm);
}