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

#include "ieee754sp.h"

union ieee754sp ieee754sp_mul(union ieee754sp x, union ieee754sp y)
{
        int re;
        int rs;
        unsigned int rm;
        unsigned short lxm;
        unsigned short hxm;
        unsigned short lym;
        unsigned short hym;
        unsigned int lrm;
        unsigned int hrm;
        unsigned int t;
        unsigned int at;

        COMPXSP;
        COMPYSP;

        EXPLODEXSP;
        EXPLODEYSP;

        ieee754_clearcx();

        FLUSHXSP;
        FLUSHYSP;

        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 ieee754sp_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 ieee754sp_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 ieee754sp_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 ieee754sp_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 ieee754sp_zero(xs ^ ys);


        case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
                SPDNORMX;
                fallthrough;
        case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
                SPDNORMY;
                break;

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

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

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

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

        /*
         * Multiply 32 bits xm, ym to give high 32 bits rm with stickness.
         */
        lxm = xm & 0xffff;
        hxm = xm >> 16;
        lym = ym & 0xffff;
        hym = ym >> 16;

        lrm = lxm * lym;        /* 16 * 16 => 32 */
        hrm = hxm * hym;        /* 16 * 16 => 32 */

        t = lxm * hym; /* 16 * 16 => 32 */
        at = lrm + (t << 16);
        hrm += at < lrm;
        lrm = at;
        hrm = hrm + (t >> 16);

        t = hxm * lym; /* 16 * 16 => 32 */
        at = lrm + (t << 16);
        hrm += at < lrm;
        lrm = at;
        hrm = hrm + (t >> 16);

        rm = hrm | (lrm != 0);

        /*
         * Sticky shift down to normal rounding precision.
         */
        if ((int) rm < 0) {
                rm = (rm >> (32 - (SP_FBITS + 1 + 3))) |
                    ((rm << (SP_FBITS + 1 + 3)) != 0);
                re++;
        } else {
                rm = (rm >> (32 - (SP_FBITS + 1 + 3 + 1))) |
                     ((rm << (SP_FBITS + 1 + 3 + 1)) != 0);
        }
        assert(rm & (SP_HIDDEN_BIT << 3));

        return ieee754sp_format(rs, re, rm);
}