1 // SPDX-License-Identifier: GPL-2.0-only 2 /* IEEE754 floating point arithmetic 3 * double precision: common utilities 4 */ 5 /* 6 * MIPS floating point support 7 * Copyright (C) 1994-2000 Algorithmics Ltd. 8 */ 9 10 #include "ieee754dp.h" 11 12 union ieee754dp ieee754dp_mul(union ieee754dp x, union ieee754dp y) 13 { 14 int re; 15 int rs; 16 u64 rm; 17 unsigned int lxm; 18 unsigned int hxm; 19 unsigned int lym; 20 unsigned int hym; 21 u64 lrm; 22 u64 hrm; 23 u64 t; 24 u64 at; 25 26 COMPXDP; 27 COMPYDP; 28 29 EXPLODEXDP; 30 EXPLODEYDP; 31 32 ieee754_clearcx(); 33 34 FLUSHXDP; 35 FLUSHYDP; 36 37 switch (CLPAIR(xc, yc)) { 38 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): 39 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN): 40 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN): 41 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN): 42 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN): 43 return ieee754dp_nanxcpt(y); 44 45 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN): 46 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): 47 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO): 48 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM): 49 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM): 50 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF): 51 return ieee754dp_nanxcpt(x); 52 53 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN): 54 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN): 55 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN): 56 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN): 57 return y; 58 59 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN): 60 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO): 61 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM): 62 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM): 63 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF): 64 return x; 65 66 67 /* 68 * Infinity handling 69 */ 70 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO): 71 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF): 72 ieee754_setcx(IEEE754_INVALID_OPERATION); 73 return ieee754dp_indef(); 74 75 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF): 76 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF): 77 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM): 78 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM): 79 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF): 80 return ieee754dp_inf(xs ^ ys); 81 82 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): 83 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM): 84 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM): 85 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO): 86 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO): 87 return ieee754dp_zero(xs ^ ys); 88 89 90 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): 91 DPDNORMX; 92 fallthrough; 93 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): 94 DPDNORMY; 95 break; 96 97 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): 98 DPDNORMX; 99 break; 100 101 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): 102 break; 103 } 104 /* rm = xm * ym, re = xe+ye basically */ 105 assert(xm & DP_HIDDEN_BIT); 106 assert(ym & DP_HIDDEN_BIT); 107 108 re = xe + ye; 109 rs = xs ^ ys; 110 111 /* shunt to top of word */ 112 xm <<= 64 - (DP_FBITS + 1); 113 ym <<= 64 - (DP_FBITS + 1); 114 115 /* 116 * Multiply 64 bits xm, ym to give high 64 bits rm with stickness. 117 */ 118 119 lxm = xm; 120 hxm = xm >> 32; 121 lym = ym; 122 hym = ym >> 32; 123 124 lrm = DPXMULT(lxm, lym); 125 hrm = DPXMULT(hxm, hym); 126 127 t = DPXMULT(lxm, hym); 128 129 at = lrm + (t << 32); 130 hrm += at < lrm; 131 lrm = at; 132 133 hrm = hrm + (t >> 32); 134 135 t = DPXMULT(hxm, lym); 136 137 at = lrm + (t << 32); 138 hrm += at < lrm; 139 lrm = at; 140 141 hrm = hrm + (t >> 32); 142 143 rm = hrm | (lrm != 0); 144 145 /* 146 * Sticky shift down to normal rounding precision. 147 */ 148 if ((s64) rm < 0) { 149 rm = (rm >> (64 - (DP_FBITS + 1 + 3))) | 150 ((rm << (DP_FBITS + 1 + 3)) != 0); 151 re++; 152 } else { 153 rm = (rm >> (64 - (DP_FBITS + 1 + 3 + 1))) | 154 ((rm << (DP_FBITS + 1 + 3 + 1)) != 0); 155 } 156 assert(rm & (DP_HIDDEN_BIT << 3)); 157 158 return ieee754dp_format(rs, re, rm); 159 } 160