1 // SPDX-License-Identifier: GPL-2.0-only 2 /* 3 * IEEE754 floating point arithmetic 4 * double precision: MADDF.f (Fused Multiply Add) 5 * MADDF.fmt: FPR[fd] = FPR[fd] + (FPR[fs] x FPR[ft]) 6 * 7 * MIPS floating point support 8 * Copyright (C) 2015 Imagination Technologies, Ltd. 9 * Author: Markos Chandras <markos.chandras@imgtec.com> 10 */ 11 12 #include "ieee754dp.h" 13 14 15 /* 128 bits shift right logical with rounding. */ 16 static void srl128(u64 *hptr, u64 *lptr, int count) 17 { 18 u64 low; 19 20 if (count >= 128) { 21 *lptr = *hptr != 0 || *lptr != 0; 22 *hptr = 0; 23 } else if (count >= 64) { 24 if (count == 64) { 25 *lptr = *hptr | (*lptr != 0); 26 } else { 27 low = *lptr; 28 *lptr = *hptr >> (count - 64); 29 *lptr |= (*hptr << (128 - count)) != 0 || low != 0; 30 } 31 *hptr = 0; 32 } else { 33 low = *lptr; 34 *lptr = low >> count | *hptr << (64 - count); 35 *lptr |= (low << (64 - count)) != 0; 36 *hptr = *hptr >> count; 37 } 38 } 39 40 static union ieee754dp _dp_maddf(union ieee754dp z, union ieee754dp x, 41 union ieee754dp y, enum maddf_flags flags) 42 { 43 int re; 44 int rs; 45 unsigned int lxm; 46 unsigned int hxm; 47 unsigned int lym; 48 unsigned int hym; 49 u64 lrm; 50 u64 hrm; 51 u64 lzm; 52 u64 hzm; 53 u64 t; 54 u64 at; 55 int s; 56 57 COMPXDP; 58 COMPYDP; 59 COMPZDP; 60 61 EXPLODEXDP; 62 EXPLODEYDP; 63 EXPLODEZDP; 64 65 FLUSHXDP; 66 FLUSHYDP; 67 FLUSHZDP; 68 69 ieee754_clearcx(); 70 71 rs = xs ^ ys; 72 if (flags & MADDF_NEGATE_PRODUCT) 73 rs ^= 1; 74 if (flags & MADDF_NEGATE_ADDITION) 75 zs ^= 1; 76 77 /* 78 * Handle the cases when at least one of x, y or z is a NaN. 79 * Order of precedence is sNaN, qNaN and z, x, y. 80 */ 81 if (zc == IEEE754_CLASS_SNAN) 82 return ieee754dp_nanxcpt(z); 83 if (xc == IEEE754_CLASS_SNAN) 84 return ieee754dp_nanxcpt(x); 85 if (yc == IEEE754_CLASS_SNAN) 86 return ieee754dp_nanxcpt(y); 87 if (zc == IEEE754_CLASS_QNAN) 88 return z; 89 if (xc == IEEE754_CLASS_QNAN) 90 return x; 91 if (yc == IEEE754_CLASS_QNAN) 92 return y; 93 94 if (zc == IEEE754_CLASS_DNORM) 95 DPDNORMZ; 96 /* ZERO z cases are handled separately below */ 97 98 switch (CLPAIR(xc, yc)) { 99 100 /* 101 * Infinity handling 102 */ 103 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO): 104 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF): 105 ieee754_setcx(IEEE754_INVALID_OPERATION); 106 return ieee754dp_indef(); 107 108 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF): 109 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF): 110 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM): 111 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM): 112 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF): 113 if ((zc == IEEE754_CLASS_INF) && (zs != rs)) { 114 /* 115 * Cases of addition of infinities with opposite signs 116 * or subtraction of infinities with same signs. 117 */ 118 ieee754_setcx(IEEE754_INVALID_OPERATION); 119 return ieee754dp_indef(); 120 } 121 /* 122 * z is here either not an infinity, or an infinity having the 123 * same sign as product (x*y). The result must be an infinity, 124 * and its sign is determined only by the sign of product (x*y). 125 */ 126 return ieee754dp_inf(rs); 127 128 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): 129 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM): 130 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM): 131 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO): 132 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO): 133 if (zc == IEEE754_CLASS_INF) 134 return ieee754dp_inf(zs); 135 if (zc == IEEE754_CLASS_ZERO) { 136 /* Handle cases +0 + (-0) and similar ones. */ 137 if (zs == rs) 138 /* 139 * Cases of addition of zeros of equal signs 140 * or subtraction of zeroes of opposite signs. 141 * The sign of the resulting zero is in any 142 * such case determined only by the sign of z. 143 */ 144 return z; 145 146 return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD); 147 } 148 /* x*y is here 0, and z is not 0, so just return z */ 149 return z; 150 151 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): 152 DPDNORMX; 153 fallthrough; 154 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): 155 if (zc == IEEE754_CLASS_INF) 156 return ieee754dp_inf(zs); 157 DPDNORMY; 158 break; 159 160 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): 161 if (zc == IEEE754_CLASS_INF) 162 return ieee754dp_inf(zs); 163 DPDNORMX; 164 break; 165 166 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): 167 if (zc == IEEE754_CLASS_INF) 168 return ieee754dp_inf(zs); 169 /* continue to real computations */ 170 } 171 172 /* Finally get to do some computation */ 173 174 /* 175 * Do the multiplication bit first 176 * 177 * rm = xm * ym, re = xe + ye basically 178 * 179 * At this point xm and ym should have been normalized. 180 */ 181 assert(xm & DP_HIDDEN_BIT); 182 assert(ym & DP_HIDDEN_BIT); 183 184 re = xe + ye; 185 186 /* shunt to top of word */ 187 xm <<= 64 - (DP_FBITS + 1); 188 ym <<= 64 - (DP_FBITS + 1); 189 190 /* 191 * Multiply 64 bits xm and ym to give 128 bits result in hrm:lrm. 192 */ 193 194 lxm = xm; 195 hxm = xm >> 32; 196 lym = ym; 197 hym = ym >> 32; 198 199 lrm = DPXMULT(lxm, lym); 200 hrm = DPXMULT(hxm, hym); 201 202 t = DPXMULT(lxm, hym); 203 204 at = lrm + (t << 32); 205 hrm += at < lrm; 206 lrm = at; 207 208 hrm = hrm + (t >> 32); 209 210 t = DPXMULT(hxm, lym); 211 212 at = lrm + (t << 32); 213 hrm += at < lrm; 214 lrm = at; 215 216 hrm = hrm + (t >> 32); 217 218 /* Put explicit bit at bit 126 if necessary */ 219 if ((int64_t)hrm < 0) { 220 lrm = (hrm << 63) | (lrm >> 1); 221 hrm = hrm >> 1; 222 re++; 223 } 224 225 assert(hrm & (1 << 62)); 226 227 if (zc == IEEE754_CLASS_ZERO) { 228 /* 229 * Move explicit bit from bit 126 to bit 55 since the 230 * ieee754dp_format code expects the mantissa to be 231 * 56 bits wide (53 + 3 rounding bits). 232 */ 233 srl128(&hrm, &lrm, (126 - 55)); 234 return ieee754dp_format(rs, re, lrm); 235 } 236 237 /* Move explicit bit from bit 52 to bit 126 */ 238 lzm = 0; 239 hzm = zm << 10; 240 assert(hzm & (1 << 62)); 241 242 /* Make the exponents the same */ 243 if (ze > re) { 244 /* 245 * Have to shift y fraction right to align. 246 */ 247 s = ze - re; 248 srl128(&hrm, &lrm, s); 249 re += s; 250 } else if (re > ze) { 251 /* 252 * Have to shift x fraction right to align. 253 */ 254 s = re - ze; 255 srl128(&hzm, &lzm, s); 256 ze += s; 257 } 258 assert(ze == re); 259 assert(ze <= DP_EMAX); 260 261 /* Do the addition */ 262 if (zs == rs) { 263 /* 264 * Generate 128 bit result by adding two 127 bit numbers 265 * leaving result in hzm:lzm, zs and ze. 266 */ 267 hzm = hzm + hrm + (lzm > (lzm + lrm)); 268 lzm = lzm + lrm; 269 if ((int64_t)hzm < 0) { /* carry out */ 270 srl128(&hzm, &lzm, 1); 271 ze++; 272 } 273 } else { 274 if (hzm > hrm || (hzm == hrm && lzm >= lrm)) { 275 hzm = hzm - hrm - (lzm < lrm); 276 lzm = lzm - lrm; 277 } else { 278 hzm = hrm - hzm - (lrm < lzm); 279 lzm = lrm - lzm; 280 zs = rs; 281 } 282 if (lzm == 0 && hzm == 0) 283 return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD); 284 285 /* 286 * Put explicit bit at bit 126 if necessary. 287 */ 288 if (hzm == 0) { 289 /* left shift by 63 or 64 bits */ 290 if ((int64_t)lzm < 0) { 291 /* MSB of lzm is the explicit bit */ 292 hzm = lzm >> 1; 293 lzm = lzm << 63; 294 ze -= 63; 295 } else { 296 hzm = lzm; 297 lzm = 0; 298 ze -= 64; 299 } 300 } 301 302 t = 0; 303 while ((hzm >> (62 - t)) == 0) 304 t++; 305 306 assert(t <= 62); 307 if (t) { 308 hzm = hzm << t | lzm >> (64 - t); 309 lzm = lzm << t; 310 ze -= t; 311 } 312 } 313 314 /* 315 * Move explicit bit from bit 126 to bit 55 since the 316 * ieee754dp_format code expects the mantissa to be 317 * 56 bits wide (53 + 3 rounding bits). 318 */ 319 srl128(&hzm, &lzm, (126 - 55)); 320 321 return ieee754dp_format(zs, ze, lzm); 322 } 323 324 union ieee754dp ieee754dp_maddf(union ieee754dp z, union ieee754dp x, 325 union ieee754dp y) 326 { 327 return _dp_maddf(z, x, y, 0); 328 } 329 330 union ieee754dp ieee754dp_msubf(union ieee754dp z, union ieee754dp x, 331 union ieee754dp y) 332 { 333 return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT); 334 } 335 336 union ieee754dp ieee754dp_madd(union ieee754dp z, union ieee754dp x, 337 union ieee754dp y) 338 { 339 return _dp_maddf(z, x, y, 0); 340 } 341 342 union ieee754dp ieee754dp_msub(union ieee754dp z, union ieee754dp x, 343 union ieee754dp y) 344 { 345 return _dp_maddf(z, x, y, MADDF_NEGATE_ADDITION); 346 } 347 348 union ieee754dp ieee754dp_nmadd(union ieee754dp z, union ieee754dp x, 349 union ieee754dp y) 350 { 351 return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT|MADDF_NEGATE_ADDITION); 352 } 353 354 union ieee754dp ieee754dp_nmsub(union ieee754dp z, union ieee754dp x, 355 union ieee754dp y) 356 { 357 return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT); 358 } 359