mips/math-emu/dp_mul.c

0001 // SPDX-License-Identifier: GPL-2.0-only
0002 /* IEEE754 floating point arithmetic
0003  * double precision: common utilities
0004  */
0005 /*
0006  * MIPS floating point support
0007  * Copyright (C) 1994-2000 Algorithmics Ltd.
0008  */
0009
0010 #include "ieee754dp.h"
0011
0012 union ieee754dp ieee754dp_mul(union ieee754dp x, union ieee754dp y)
0013 {
0014     int re;
0015     int rs;
0016     u64 rm;
0017     unsigned int lxm;
0018     unsigned int hxm;
0019     unsigned int lym;
0020     unsigned int hym;
0021     u64 lrm;
0022     u64 hrm;
0023     u64 t;
0024     u64 at;
0025
0026     COMPXDP;
0027     COMPYDP;
0028
0029     EXPLODEXDP;
0030     EXPLODEYDP;
0031
0032     ieee754_clearcx();
0033
0034     FLUSHXDP;
0035     FLUSHYDP;
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 ieee754dp_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 ieee754dp_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 ieee754dp_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 ieee754dp_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 ieee754dp_zero(xs ^ ys);
0088
0089
0090     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
0091         DPDNORMX;
0092         fallthrough;
0093     case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
0094         DPDNORMY;
0095         break;
0096
0097     case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
0098         DPDNORMX;
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 & DP_HIDDEN_BIT);
0106     assert(ym & DP_HIDDEN_BIT);
0107
0108     re = xe + ye;
0109     rs = xs ^ ys;
0110
0111     /* shunt to top of word */
0112     xm <<= 64 - (DP_FBITS + 1);
0113     ym <<= 64 - (DP_FBITS + 1);
0114
0115     /*
0116      * Multiply 64 bits xm, ym to give high 64 bits rm with stickness.
0117      */
0118
0119     lxm = xm;
0120     hxm = xm >> 32;
0121     lym = ym;
0122     hym = ym >> 32;
0123
0124     lrm = DPXMULT(lxm, lym);
0125     hrm = DPXMULT(hxm, hym);
0126
0127     t = DPXMULT(lxm, hym);
0128
0129     at = lrm + (t << 32);
0130     hrm += at < lrm;
0131     lrm = at;
0132
0133     hrm = hrm + (t >> 32);
0134
0135     t = DPXMULT(hxm, lym);
0136
0137     at = lrm + (t << 32);
0138     hrm += at < lrm;
0139     lrm = at;
0140
0141     hrm = hrm + (t >> 32);
0142
0143     rm = hrm | (lrm != 0);
0144
0145     /*
0146      * Sticky shift down to normal rounding precision.
0147      */
0148     if ((s64) rm < 0) {
0149         rm = (rm >> (64 - (DP_FBITS + 1 + 3))) |
0150              ((rm << (DP_FBITS + 1 + 3)) != 0);
0151         re++;
0152     } else {
0153         rm = (rm >> (64 - (DP_FBITS + 1 + 3 + 1))) |
0154              ((rm << (DP_FBITS + 1 + 3 + 1)) != 0);
0155     }
0156     assert(rm & (DP_HIDDEN_BIT << 3));
0157
0158     return ieee754dp_format(rs, re, rm);
0159 }