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_sub(union ieee754dp x, union ieee754dp y) 13 { 14 int s; 15 16 COMPXDP; 17 COMPYDP; 18 19 EXPLODEXDP; 20 EXPLODEYDP; 21 22 ieee754_clearcx(); 23 24 FLUSHXDP; 25 FLUSHYDP; 26 27 switch (CLPAIR(xc, yc)) { 28 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): 29 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN): 30 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN): 31 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN): 32 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN): 33 return ieee754dp_nanxcpt(y); 34 35 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN): 36 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): 37 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO): 38 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM): 39 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM): 40 case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF): 41 return ieee754dp_nanxcpt(x); 42 43 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN): 44 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN): 45 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN): 46 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN): 47 return y; 48 49 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN): 50 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO): 51 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM): 52 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM): 53 case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF): 54 return x; 55 56 57 /* 58 * Infinity handling 59 */ 60 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF): 61 if (xs != ys) 62 return x; 63 ieee754_setcx(IEEE754_INVALID_OPERATION); 64 return ieee754dp_indef(); 65 66 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF): 67 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF): 68 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF): 69 return ieee754dp_inf(ys ^ 1); 70 71 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO): 72 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM): 73 case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM): 74 return x; 75 76 /* 77 * Zero handling 78 */ 79 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): 80 if (xs != ys) 81 return x; 82 else 83 return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD); 84 85 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO): 86 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO): 87 return x; 88 89 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM): 90 case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM): 91 /* quick fix up */ 92 DPSIGN(y) ^= 1; 93 return y; 94 95 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): 96 DPDNORMX; 97 /* fall through */ 98 99 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): 100 /* normalize ym,ye */ 101 DPDNORMY; 102 break; 103 104 case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): 105 /* normalize xm,xe */ 106 DPDNORMX; 107 break; 108 109 case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): 110 break; 111 } 112 /* flip sign of y and handle as add */ 113 ys ^= 1; 114 115 assert(xm & DP_HIDDEN_BIT); 116 assert(ym & DP_HIDDEN_BIT); 117 118 119 /* provide guard,round and stick bit dpace */ 120 xm <<= 3; 121 ym <<= 3; 122 123 if (xe > ye) { 124 /* 125 * Have to shift y fraction right to align 126 */ 127 s = xe - ye; 128 ym = XDPSRS(ym, s); 129 ye += s; 130 } else if (ye > xe) { 131 /* 132 * Have to shift x fraction right to align 133 */ 134 s = ye - xe; 135 xm = XDPSRS(xm, s); 136 xe += s; 137 } 138 assert(xe == ye); 139 assert(xe <= DP_EMAX); 140 141 if (xs == ys) { 142 /* generate 28 bit result of adding two 27 bit numbers 143 */ 144 xm = xm + ym; 145 146 if (xm >> (DP_FBITS + 1 + 3)) { /* carry out */ 147 xm = XDPSRS1(xm); /* shift preserving sticky */ 148 xe++; 149 } 150 } else { 151 if (xm >= ym) { 152 xm = xm - ym; 153 } else { 154 xm = ym - xm; 155 xs = ys; 156 } 157 if (xm == 0) { 158 if (ieee754_csr.rm == FPU_CSR_RD) 159 return ieee754dp_zero(1); /* round negative inf. => sign = -1 */ 160 else 161 return ieee754dp_zero(0); /* other round modes => sign = 1 */ 162 } 163 164 /* normalize to rounding precision 165 */ 166 while ((xm >> (DP_FBITS + 3)) == 0) { 167 xm <<= 1; 168 xe--; 169 } 170 } 171 172 return ieee754dp_format(xs, xe, xm); 173 } 174