19d5a6349SThomas Gleixner // SPDX-License-Identifier: GPL-2.0-only
21da177e4SLinus Torvalds /* IEEE754 floating point arithmetic
31da177e4SLinus Torvalds * double precision: common utilities
41da177e4SLinus Torvalds */
51da177e4SLinus Torvalds /*
61da177e4SLinus Torvalds * MIPS floating point support
71da177e4SLinus Torvalds * Copyright (C) 1994-2000 Algorithmics Ltd.
81da177e4SLinus Torvalds */
91da177e4SLinus Torvalds
101da177e4SLinus Torvalds #include "ieee754dp.h"
111da177e4SLinus Torvalds
ieee754dp_mul(union ieee754dp x,union ieee754dp y)122209bcb1SRalf Baechle union ieee754dp ieee754dp_mul(union ieee754dp x, union ieee754dp y)
131da177e4SLinus Torvalds {
143f7cac41SRalf Baechle int re;
153f7cac41SRalf Baechle int rs;
163f7cac41SRalf Baechle u64 rm;
17a58f85b5SAleksandar Markovic unsigned int lxm;
18a58f85b5SAleksandar Markovic unsigned int hxm;
19a58f85b5SAleksandar Markovic unsigned int lym;
20a58f85b5SAleksandar Markovic unsigned int hym;
213f7cac41SRalf Baechle u64 lrm;
223f7cac41SRalf Baechle u64 hrm;
233f7cac41SRalf Baechle u64 t;
243f7cac41SRalf Baechle u64 at;
253f7cac41SRalf Baechle
261da177e4SLinus Torvalds COMPXDP;
271da177e4SLinus Torvalds COMPYDP;
281da177e4SLinus Torvalds
291da177e4SLinus Torvalds EXPLODEXDP;
301da177e4SLinus Torvalds EXPLODEYDP;
311da177e4SLinus Torvalds
329e8bad1fSRalf Baechle ieee754_clearcx();
331da177e4SLinus Torvalds
341da177e4SLinus Torvalds FLUSHXDP;
351da177e4SLinus Torvalds FLUSHYDP;
361da177e4SLinus Torvalds
371da177e4SLinus Torvalds switch (CLPAIR(xc, yc)) {
381da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
391da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
401da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
411da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
421da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
43d5afa7e9SMaciej W. Rozycki return ieee754dp_nanxcpt(y);
44d5afa7e9SMaciej W. Rozycki
45d5afa7e9SMaciej W. Rozycki case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
46d5afa7e9SMaciej W. Rozycki case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
471da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
481da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
491da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
501da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
51d5afa7e9SMaciej W. Rozycki return ieee754dp_nanxcpt(x);
521da177e4SLinus Torvalds
531da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
541da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
551da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
561da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
571da177e4SLinus Torvalds return y;
581da177e4SLinus Torvalds
591da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
601da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
611da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
621da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
631da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
641da177e4SLinus Torvalds return x;
651da177e4SLinus Torvalds
661da177e4SLinus Torvalds
673f7cac41SRalf Baechle /*
683f7cac41SRalf Baechle * Infinity handling
693f7cac41SRalf Baechle */
701da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
711da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
729e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INVALID_OPERATION);
7390efba36SRalf Baechle return ieee754dp_indef();
741da177e4SLinus Torvalds
751da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
761da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
771da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
781da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
791da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
801da177e4SLinus Torvalds return ieee754dp_inf(xs ^ ys);
811da177e4SLinus Torvalds
821da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
831da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
841da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
851da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
861da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
871da177e4SLinus Torvalds return ieee754dp_zero(xs ^ ys);
881da177e4SLinus Torvalds
891da177e4SLinus Torvalds
901da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
911da177e4SLinus Torvalds DPDNORMX;
92c9b02990SLiangliang Huang fallthrough;
931da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
941da177e4SLinus Torvalds DPDNORMY;
951da177e4SLinus Torvalds break;
961da177e4SLinus Torvalds
971da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
981da177e4SLinus Torvalds DPDNORMX;
991da177e4SLinus Torvalds break;
1001da177e4SLinus Torvalds
1011da177e4SLinus Torvalds case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
1021da177e4SLinus Torvalds break;
1031da177e4SLinus Torvalds }
10425985edcSLucas De Marchi /* rm = xm * ym, re = xe+ye basically */
1051da177e4SLinus Torvalds assert(xm & DP_HIDDEN_BIT);
1061da177e4SLinus Torvalds assert(ym & DP_HIDDEN_BIT);
1073f7cac41SRalf Baechle
1083f7cac41SRalf Baechle re = xe + ye;
1093f7cac41SRalf Baechle rs = xs ^ ys;
1101da177e4SLinus Torvalds
1111da177e4SLinus Torvalds /* shunt to top of word */
112ad8fb553SRalf Baechle xm <<= 64 - (DP_FBITS + 1);
113ad8fb553SRalf Baechle ym <<= 64 - (DP_FBITS + 1);
1141da177e4SLinus Torvalds
1153f7cac41SRalf Baechle /*
11695bff241SPaul Burton * Multiply 64 bits xm, ym to give high 64 bits rm with stickness.
1171da177e4SLinus Torvalds */
1181da177e4SLinus Torvalds
1193f7cac41SRalf Baechle lxm = xm;
1203f7cac41SRalf Baechle hxm = xm >> 32;
1213f7cac41SRalf Baechle lym = ym;
1223f7cac41SRalf Baechle hym = ym >> 32;
1231da177e4SLinus Torvalds
1241da177e4SLinus Torvalds lrm = DPXMULT(lxm, lym);
1251da177e4SLinus Torvalds hrm = DPXMULT(hxm, hym);
1261da177e4SLinus Torvalds
1273f7cac41SRalf Baechle t = DPXMULT(lxm, hym);
1281da177e4SLinus Torvalds
1293f7cac41SRalf Baechle at = lrm + (t << 32);
1301da177e4SLinus Torvalds hrm += at < lrm;
1311da177e4SLinus Torvalds lrm = at;
1323f7cac41SRalf Baechle
1331da177e4SLinus Torvalds hrm = hrm + (t >> 32);
1343f7cac41SRalf Baechle
1353f7cac41SRalf Baechle t = DPXMULT(hxm, lym);
1363f7cac41SRalf Baechle
1373f7cac41SRalf Baechle at = lrm + (t << 32);
1383f7cac41SRalf Baechle hrm += at < lrm;
1393f7cac41SRalf Baechle lrm = at;
1403f7cac41SRalf Baechle
1413f7cac41SRalf Baechle hrm = hrm + (t >> 32);
1423f7cac41SRalf Baechle
1431da177e4SLinus Torvalds rm = hrm | (lrm != 0);
1441da177e4SLinus Torvalds
1451da177e4SLinus Torvalds /*
1463f7cac41SRalf Baechle * Sticky shift down to normal rounding precision.
1471da177e4SLinus Torvalds */
1481da177e4SLinus Torvalds if ((s64) rm < 0) {
1493f7cac41SRalf Baechle rm = (rm >> (64 - (DP_FBITS + 1 + 3))) |
150ad8fb553SRalf Baechle ((rm << (DP_FBITS + 1 + 3)) != 0);
1511da177e4SLinus Torvalds re++;
1521da177e4SLinus Torvalds } else {
1533f7cac41SRalf Baechle rm = (rm >> (64 - (DP_FBITS + 1 + 3 + 1))) |
154ad8fb553SRalf Baechle ((rm << (DP_FBITS + 1 + 3 + 1)) != 0);
1551da177e4SLinus Torvalds }
1561da177e4SLinus Torvalds assert(rm & (DP_HIDDEN_BIT << 3));
15790efba36SRalf Baechle
15890efba36SRalf Baechle return ieee754dp_format(rs, re, rm);
1591da177e4SLinus Torvalds }
160