^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 1) // SPDX-License-Identifier: GPL-2.0-only
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 2) /* IEEE754 floating point arithmetic
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 3) * double precision: common utilities
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 4) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 5) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 6) * MIPS floating point support
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 7) * Copyright (C) 1994-2000 Algorithmics Ltd.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 8) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 9)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 10) #include "ieee754dp.h"
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 11)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 12) union ieee754dp ieee754dp_add(union ieee754dp x, union ieee754dp y)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 13) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 14) int s;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 15)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 16) COMPXDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 17) COMPYDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 18)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 19) EXPLODEXDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 20) EXPLODEYDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 21)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 22) ieee754_clearcx();
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 23)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 24) FLUSHXDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 25) FLUSHYDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 26)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 27) switch (CLPAIR(xc, yc)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 28) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 29) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 30) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 31) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 32) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 33) return ieee754dp_nanxcpt(y);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 34)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 35) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 36) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 37) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 38) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 39) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 40) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 41) return ieee754dp_nanxcpt(x);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 42)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 43) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 44) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 45) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 46) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 47) return y;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 48)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 49) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 50) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 51) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 52) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 53) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 54) return x;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 55)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 56)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 57) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 58) * Infinity handling
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 59) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 60) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 61) if (xs == ys)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 62) return x;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 63) ieee754_setcx(IEEE754_INVALID_OPERATION);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 64) return ieee754dp_indef();
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 65)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 66) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 67) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 68) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 69) return y;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 70)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 71) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 72) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 73) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 74) return x;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 75)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 76) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 77) * Zero handling
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 78) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 79) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 80) if (xs == ys)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 81) return x;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 82) else
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 83) return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 84)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 85) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 86) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 87) return x;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 88)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 89) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 90) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 91) return y;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 92)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 93) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 94) DPDNORMX;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 95) fallthrough;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 96) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 97) DPDNORMY;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 98) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 99)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 100) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 101) DPDNORMX;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 102) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 103)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 104) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 105) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 106) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 107) assert(xm & DP_HIDDEN_BIT);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 108) assert(ym & DP_HIDDEN_BIT);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 109)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 110) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 111) * Provide guard,round and stick bit space.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 112) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 113) xm <<= 3;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 114) ym <<= 3;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 115)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 116) if (xe > ye) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 117) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 118) * Have to shift y fraction right to align.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 119) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 120) s = xe - ye;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 121) ym = XDPSRS(ym, s);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 122) ye += s;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 123) } else if (ye > xe) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 124) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 125) * Have to shift x fraction right to align.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 126) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 127) s = ye - xe;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 128) xm = XDPSRS(xm, s);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 129) xe += s;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 130) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 131) assert(xe == ye);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 132) assert(xe <= DP_EMAX);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 133)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 134) if (xs == ys) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 135) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 136) * Generate 28 bit result of adding two 27 bit numbers
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 137) * leaving result in xm, xs and xe.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 138) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 139) xm = xm + ym;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 140)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 141) if (xm >> (DP_FBITS + 1 + 3)) { /* carry out */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 142) xm = XDPSRS1(xm);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 143) xe++;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 144) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 145) } else {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 146) if (xm >= ym) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 147) xm = xm - ym;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 148) } else {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 149) xm = ym - xm;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 150) xs = ys;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 151) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 152) if (xm == 0)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 153) return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 154)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 155) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 156) * Normalize to rounding precision.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 157) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 158) while ((xm >> (DP_FBITS + 3)) == 0) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 159) xm <<= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 160) xe--;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 161) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 162) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 163)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 164) return ieee754dp_format(xs, xe, xm);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 165) }