^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_div(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) u64 rm;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 15) int re;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 16) u64 bm;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 17)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 18) COMPXDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 19) COMPYDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 20)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 21) EXPLODEXDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 22) EXPLODEYDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 23)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 24) ieee754_clearcx();
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 25)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 26) FLUSHXDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 27) FLUSHYDP;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 28)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 29) switch (CLPAIR(xc, yc)) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 30) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 31) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 32) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 33) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 34) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 35) return ieee754dp_nanxcpt(y);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 36)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 37) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 38) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 39) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 40) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 41) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 42) case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 43) return ieee754dp_nanxcpt(x);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 44)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 45) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 46) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 47) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 48) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 49) return y;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 50)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 51) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 52) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 53) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 54) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 55) case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 56) return x;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 57)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 58)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 59) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 60) * Infinity handling
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 61) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 62) case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
^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 ieee754dp_zero(xs ^ ys);
^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 ieee754dp_inf(xs ^ ys);
^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) ieee754_setcx(IEEE754_INVALID_OPERATION);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 81) return ieee754dp_indef();
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 82)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 83) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 84) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 85) ieee754_setcx(IEEE754_ZERO_DIVIDE);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 86) return ieee754dp_inf(xs ^ ys);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 87)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 88) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 89) case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 90) return ieee754dp_zero(xs == ys ? 0 : 1);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 91)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 92) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 93) DPDNORMX;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 94) fallthrough;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 95) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 96) DPDNORMY;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 97) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 98)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 99) case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 100) DPDNORMX;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 101) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 102)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 103) case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 104) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 105) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 106) assert(xm & DP_HIDDEN_BIT);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 107) assert(ym & DP_HIDDEN_BIT);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 108)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 109) /* provide rounding space */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 110) xm <<= 3;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 111) ym <<= 3;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 112)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 113) /* now the dirty work */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 114)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 115) rm = 0;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 116) re = xe - ye;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 117)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 118) for (bm = DP_MBIT(DP_FBITS + 2); bm; bm >>= 1) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 119) if (xm >= ym) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 120) xm -= ym;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 121) rm |= bm;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 122) if (xm == 0)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 123) break;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 124) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 125) xm <<= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 126) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 127)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 128) rm <<= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 129) if (xm)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 130) rm |= 1; /* have remainder, set sticky */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 131)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 132) assert(rm);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 133)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 134) /*
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 135) * Normalise rm to rounding precision ?
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 136) */
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 137) while ((rm >> (DP_FBITS + 3)) == 0) {
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 138) rm <<= 1;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 139) re--;
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 140) }
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 141)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 142) return ieee754dp_format(xs == ys ? 0 : 1, re, rm);
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 143) }