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 /* fall through */ 93 94 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): 95 DPDNORMY; 96 break; 97 98 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): 99 DPDNORMX; 100 break; 101 102 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): 103 break; 104 } 105 /* rm = xm * ym, re = xe+ye basically */ 106 assert(xm & DP_HIDDEN_BIT); 107 assert(ym & DP_HIDDEN_BIT); 108 109 re = xe + ye; 110 rs = xs ^ ys; 111 112 /* shunt to top of word */ 113 xm <<= 64 - (DP_FBITS + 1); 114 ym <<= 64 - (DP_FBITS + 1); 115 116 /* 117 * Multiply 64 bits xm, ym to give high 64 bits rm with stickness. 118 */ 119 120 lxm = xm; 121 hxm = xm >> 32; 122 lym = ym; 123 hym = ym >> 32; 124 125 lrm = DPXMULT(lxm, lym); 126 hrm = DPXMULT(hxm, hym); 127 128 t = DPXMULT(lxm, hym); 129 130 at = lrm + (t << 32); 131 hrm += at < lrm; 132 lrm = at; 133 134 hrm = hrm + (t >> 32); 135 136 t = DPXMULT(hxm, lym); 137 138 at = lrm + (t << 32); 139 hrm += at < lrm; 140 lrm = at; 141 142 hrm = hrm + (t >> 32); 143 144 rm = hrm | (lrm != 0); 145 146 /* 147 * Sticky shift down to normal rounding precision. 148 */ 149 if ((s64) rm < 0) { 150 rm = (rm >> (64 - (DP_FBITS + 1 + 3))) | 151 ((rm << (DP_FBITS + 1 + 3)) != 0); 152 re++; 153 } else { 154 rm = (rm >> (64 - (DP_FBITS + 1 + 3 + 1))) | 155 ((rm << (DP_FBITS + 1 + 3 + 1)) != 0); 156 } 157 assert(rm & (DP_HIDDEN_BIT << 3)); 158 159 return ieee754dp_format(rs, re, rm); 160 } 161