mips/math-emu/sp_mul.c

0001 // SPDX-License-Identifier: GPL-2.0-only
0002 /* IEEE754 floating point arithmetic
0003  * single precision
0004  */
0005 /*
0006  * MIPS floating point support
0007  * Copyright (C) 1994-2000 Algorithmics Ltd.
0008  */
0009
0010 #include "ieee754sp.h"
0011
0012 union ieee754sp ieee754sp_mul(union ieee754sp x, union ieee754sp y)
0013 {
0014     int re;
0015     int rs;
0016     unsigned int rm;
0017     unsigned short lxm;
0018     unsigned short hxm;
0019     unsigned short lym;
0020     unsigned short hym;
0021     unsigned int lrm;
0022     unsigned int hrm;
0023     unsigned int t;
0024     unsigned int at;
0025
0026     COMPXSP;
0027     COMPYSP;
0028
0029     EXPLODEXSP;
0030     EXPLODEYSP;
0031
0032     ieee754_clearcx();
0033
0034     FLUSHXSP;
0035     FLUSHYSP;
0036
0037     switch (CLPAIR(xc, yc)) {
0038     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
0039     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
0040     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
0041     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
0042     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
0043         return ieee754sp_nanxcpt(y);
0044
0045     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
0046     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
0047     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
0048     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
0049     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
0050     case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
0051         return ieee754sp_nanxcpt(x);
0052
0053     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
0054     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
0055     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
0056     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
0057         return y;
0058
0059     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
0060     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
0061     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
0062     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
0063     case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
0064         return x;
0065
0066
0067     /*
0068      * Infinity handling
0069      */
0070     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
0071     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
0072         ieee754_setcx(IEEE754_INVALID_OPERATION);
0073         return ieee754sp_indef();
0074
0075     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
0076     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
0077     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
0078     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
0079     case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
0080         return ieee754sp_inf(xs ^ ys);
0081
0082     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
0083     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
0084     case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
0085     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
0086     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
0087         return ieee754sp_zero(xs ^ ys);
0088
0089
0090     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
0091         SPDNORMX;
0092         fallthrough;
0093     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
0094         SPDNORMY;
0095         break;
0096
0097     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
0098         SPDNORMX;
0099         break;
0100
0101     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
0102         break;
0103     }
0104     /* rm = xm * ym, re = xe+ye basically */
0105     assert(xm & SP_HIDDEN_BIT);
0106     assert(ym & SP_HIDDEN_BIT);
0107
0108     re = xe + ye;
0109     rs = xs ^ ys;
0110
0111     /* shunt to top of word */
0112     xm <<= 32 - (SP_FBITS + 1);
0113     ym <<= 32 - (SP_FBITS + 1);
0114
0115     /*
0116      * Multiply 32 bits xm, ym to give high 32 bits rm with stickness.
0117      */
0118     lxm = xm & 0xffff;
0119     hxm = xm >> 16;
0120     lym = ym & 0xffff;
0121     hym = ym >> 16;
0122
0123     lrm = lxm * lym;    /* 16 * 16 => 32 */
0124     hrm = hxm * hym;    /* 16 * 16 => 32 */
0125
0126     t = lxm * hym; /* 16 * 16 => 32 */
0127     at = lrm + (t << 16);
0128     hrm += at < lrm;
0129     lrm = at;
0130     hrm = hrm + (t >> 16);
0131
0132     t = hxm * lym; /* 16 * 16 => 32 */
0133     at = lrm + (t << 16);
0134     hrm += at < lrm;
0135     lrm = at;
0136     hrm = hrm + (t >> 16);
0137
0138     rm = hrm | (lrm != 0);
0139
0140     /*
0141      * Sticky shift down to normal rounding precision.
0142      */
0143     if ((int) rm < 0) {
0144         rm = (rm >> (32 - (SP_FBITS + 1 + 3))) |
0145             ((rm << (SP_FBITS + 1 + 3)) != 0);
0146         re++;
0147     } else {
0148         rm = (rm >> (32 - (SP_FBITS + 1 + 3 + 1))) |
0149              ((rm << (SP_FBITS + 1 + 3 + 1)) != 0);
0150     }
0151     assert(rm & (SP_HIDDEN_BIT << 3));
0152
0153     return ieee754sp_format(rs, re, rm);
0154 }