Resguardo inicial

linux-7.2.3/arch/mips/math-emu/Makefile

@@ -0,0 +1,19 @@
+# SPDX-License-Identifier: GPL-2.0
+#
+# Makefile for the Linux/MIPS kernel FPU emulation.
+#
+
+obj-y	+= cp1emu.o ieee754dp.o ieee754sp.o ieee754.o \
+	   dp_div.o dp_mul.o dp_sub.o dp_add.o dp_fsp.o dp_cmp.o dp_simple.o \
+	   dp_tint.o dp_fint.o dp_rint.o dp_maddf.o dp_2008class.o dp_fmin.o \
+	   dp_fmax.o							     \
+	   sp_div.o sp_mul.o sp_sub.o sp_add.o sp_fdp.o sp_cmp.o sp_simple.o \
+	   sp_tint.o sp_fint.o sp_rint.o sp_maddf.o sp_2008class.o sp_fmin.o \
+	   sp_fmax.o							     \
+	   dsemul.o
+
+lib-y	+= ieee754d.o \
+	   dp_tlong.o dp_flong.o dp_sqrt.o \
+	   sp_tlong.o sp_flong.o sp_sqrt.o
+
+obj-$(CONFIG_DEBUG_FS) += me-debugfs.o

linux-7.2.3/arch/mips/math-emu/dp_2008class.c

@@ -0,0 +1,52 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * IEEE754 floating point arithmetic
+ * double precision: CLASS.f
+ * FPR[fd] = class(FPR[fs])
+ *
+ * MIPS floating point support
+ * Copyright (C) 2015 Imagination Technologies, Ltd.
+ * Author: Markos Chandras <markos.chandras@imgtec.com>
+ */
+
+#include "ieee754dp.h"
+
+int ieee754dp_2008class(union ieee754dp x)
+{
+	COMPXDP;
+
+	EXPLODEXDP;
+
+	/*
+	 * 10 bit mask as follows:
+	 *
+	 * bit0 = SNAN
+	 * bit1 = QNAN
+	 * bit2 = -INF
+	 * bit3 = -NORM
+	 * bit4 = -DNORM
+	 * bit5 = -ZERO
+	 * bit6 = INF
+	 * bit7 = NORM
+	 * bit8 = DNORM
+	 * bit9 = ZERO
+	 */
+
+	switch(xc) {
+	case IEEE754_CLASS_SNAN:
+		return 0x01;
+	case IEEE754_CLASS_QNAN:
+		return 0x02;
+	case IEEE754_CLASS_INF:
+		return 0x04 << (xs ? 0 : 4);
+	case IEEE754_CLASS_NORM:
+		return 0x08 << (xs ? 0 : 4);
+	case IEEE754_CLASS_DNORM:
+		return 0x10 << (xs ? 0 : 4);
+	case IEEE754_CLASS_ZERO:
+		return 0x20 << (xs ? 0 : 4);
+	default:
+		pr_err("Unknown class: %d\n", xc);
+		return 0;
+	}
+}

linux-7.2.3/arch/mips/math-emu/dp_add.c

@@ -0,0 +1,165 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_add(union ieee754dp x, union ieee754dp y)
+{
+	int s;
+
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+
+	ieee754_clearcx();
+
+	FLUSHXDP;
+	FLUSHYDP;
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754dp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754dp_nanxcpt(x);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return x;
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		if (xs == ys)
+			return x;
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754dp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+		return x;
+
+	/*
+	 * Zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		if (xs == ys)
+			return x;
+		else
+			return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		DPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		break;
+	}
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+	/*
+	 * Provide guard,round and stick bit space.
+	 */
+	xm <<= 3;
+	ym <<= 3;
+
+	if (xe > ye) {
+		/*
+		 * Have to shift y fraction right to align.
+		 */
+		s = xe - ye;
+		ym = XDPSRS(ym, s);
+		ye += s;
+	} else if (ye > xe) {
+		/*
+		 * Have to shift x fraction right to align.
+		 */
+		s = ye - xe;
+		xm = XDPSRS(xm, s);
+		xe += s;
+	}
+	assert(xe == ye);
+	assert(xe <= DP_EMAX);
+
+	if (xs == ys) {
+		/*
+		 * Generate 28 bit result of adding two 27 bit numbers
+		 * leaving result in xm, xs and xe.
+		 */
+		xm = xm + ym;
+
+		if (xm >> (DP_FBITS + 1 + 3)) { /* carry out */
+			xm = XDPSRS1(xm);
+			xe++;
+		}
+	} else {
+		if (xm >= ym) {
+			xm = xm - ym;
+		} else {
+			xm = ym - xm;
+			xs = ys;
+		}
+		if (xm == 0)
+			return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
+
+		/*
+		 * Normalize to rounding precision.
+		 */
+		while ((xm >> (DP_FBITS + 3)) == 0) {
+			xm <<= 1;
+			xe--;
+		}
+	}
+
+	return ieee754dp_format(xs, xe, xm);
+}

linux-7.2.3/arch/mips/math-emu/dp_cmp.c

@@ -0,0 +1,47 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+int ieee754dp_cmp(union ieee754dp x, union ieee754dp y, int cmp, int sig)
+{
+	s64 vx;
+	s64 vy;
+
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+	FLUSHXDP;
+	FLUSHYDP;
+	ieee754_clearcx();	/* Even clear inexact flag here */
+
+	if (ieee754_class_nan(xc) || ieee754_class_nan(yc)) {
+		if (sig ||
+		    xc == IEEE754_CLASS_SNAN || yc == IEEE754_CLASS_SNAN)
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return (cmp & IEEE754_CUN) != 0;
+	} else {
+		vx = x.bits;
+		vy = y.bits;
+
+		if (vx < 0)
+			vx = -vx ^ DP_SIGN_BIT;
+		if (vy < 0)
+			vy = -vy ^ DP_SIGN_BIT;
+
+		if (vx < vy)
+			return (cmp & IEEE754_CLT) != 0;
+		else if (vx == vy)
+			return (cmp & IEEE754_CEQ) != 0;
+		else
+			return (cmp & IEEE754_CGT) != 0;
+	}
+}

linux-7.2.3/arch/mips/math-emu/dp_div.c

@@ -0,0 +1,143 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_div(union ieee754dp x, union ieee754dp y)
+{
+	u64 rm;
+	int re;
+	u64 bm;
+
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+
+	ieee754_clearcx();
+
+	FLUSHXDP;
+	FLUSHYDP;
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754dp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754dp_nanxcpt(x);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return x;
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754dp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+		return ieee754dp_zero(xs ^ ys);
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+		return ieee754dp_inf(xs ^ ys);
+
+	/*
+	 * Zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754dp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		ieee754_setcx(IEEE754_ZERO_DIVIDE);
+		return ieee754dp_inf(xs ^ ys);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return ieee754dp_zero(xs == ys ? 0 : 1);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		DPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		break;
+	}
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+	/* provide rounding space */
+	xm <<= 3;
+	ym <<= 3;
+
+	/* now the dirty work */
+
+	rm = 0;
+	re = xe - ye;
+
+	for (bm = DP_MBIT(DP_FBITS + 2); bm; bm >>= 1) {
+		if (xm >= ym) {
+			xm -= ym;
+			rm |= bm;
+			if (xm == 0)
+				break;
+		}
+		xm <<= 1;
+	}
+
+	rm <<= 1;
+	if (xm)
+		rm |= 1;	/* have remainder, set sticky */
+
+	assert(rm);
+
+	/*
+	 * Normalise rm to rounding precision ?
+	 */
+	while ((rm >> (DP_FBITS + 3)) == 0) {
+		rm <<= 1;
+		re--;
+	}
+
+	return ieee754dp_format(xs == ys ? 0 : 1, re, rm);
+}

linux-7.2.3/arch/mips/math-emu/dp_fint.c

@@ -0,0 +1,44 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_fint(int x)
+{
+	u64 xm;
+	int xe;
+	int xs;
+
+	ieee754_clearcx();
+
+	if (x == 0)
+		return ieee754dp_zero(0);
+	if (x == 1 || x == -1)
+		return ieee754dp_one(x < 0);
+	if (x == 10 || x == -10)
+		return ieee754dp_ten(x < 0);
+
+	xs = (x < 0);
+	if (xs) {
+		if (x == (1 << 31))
+			xm = ((unsigned) 1 << 31);	/* max neg can't be safely negated */
+		else
+			xm = -x;
+	} else {
+		xm = x;
+	}
+
+	/* normalize - result can never be inexact or overflow */
+	xe = DP_FBITS;
+	while ((xm >> DP_FBITS) == 0) {
+		xm <<= 1;
+		xe--;
+	}
+	return builddp(xs, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
+}

linux-7.2.3/arch/mips/math-emu/dp_flong.c

@@ -0,0 +1,53 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_flong(s64 x)
+{
+	u64 xm;
+	int xe;
+	int xs;
+
+	ieee754_clearcx();
+
+	if (x == 0)
+		return ieee754dp_zero(0);
+	if (x == 1 || x == -1)
+		return ieee754dp_one(x < 0);
+	if (x == 10 || x == -10)
+		return ieee754dp_ten(x < 0);
+
+	xs = (x < 0);
+	if (xs) {
+		if (x == (1ULL << 63))
+			xm = (1ULL << 63);	/* max neg can't be safely negated */
+		else
+			xm = -x;
+	} else {
+		xm = x;
+	}
+
+	/* normalize */
+	xe = DP_FBITS + 3;
+	if (xm >> (DP_FBITS + 1 + 3)) {
+		/* shunt out overflow bits */
+		while (xm >> (DP_FBITS + 1 + 3)) {
+			XDPSRSX1();
+		}
+	} else {
+		/* normalize in grs extended double precision */
+		while ((xm >> (DP_FBITS + 3)) == 0) {
+			xm <<= 1;
+			xe--;
+		}
+	}
+
+	return ieee754dp_format(xs, xe, xm);
+}

linux-7.2.3/arch/mips/math-emu/dp_fmax.c

@@ -0,0 +1,252 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * IEEE754 floating point arithmetic
+ * double precision: MIN{,A}.f
+ * MIN : Scalar Floating-Point Minimum
+ * MINA: Scalar Floating-Point argument with Minimum Absolute Value
+ *
+ * MIN.D : FPR[fd] = minNum(FPR[fs],FPR[ft])
+ * MINA.D: FPR[fd] = maxNumMag(FPR[fs],FPR[ft])
+ *
+ * MIPS floating point support
+ * Copyright (C) 2015 Imagination Technologies, Ltd.
+ * Author: Markos Chandras <markos.chandras@imgtec.com>
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_fmax(union ieee754dp x, union ieee754dp y)
+{
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+
+	FLUSHXDP;
+	FLUSHYDP;
+
+	ieee754_clearcx();
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754dp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754dp_nanxcpt(x);
+
+	/*
+	 * Quiet NaN handling
+	 */
+
+	/*
+	 *    The case of both inputs quiet NaNs
+	 */
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+		return x;
+
+	/*
+	 *    The cases of exactly one input quiet NaN (numbers
+	 *    are here preferred as returned values to NaNs)
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return y;
+
+	/*
+	 * Infinity and zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return xs ? y : x;
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return ys ? x : y;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		return ieee754dp_zero(xs & ys);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		DPDNORMX;
+	}
+
+	/* Finally get to do some computation */
+
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+	/* Compare signs */
+	if (xs > ys)
+		return y;
+	else if (xs < ys)
+		return x;
+
+	/* Signs of inputs are equal, let's compare exponents */
+	if (xs == 0) {
+		/* Inputs are both positive */
+		if (xe > ye)
+			return x;
+		else if (xe < ye)
+			return y;
+	} else {
+		/* Inputs are both negative */
+		if (xe > ye)
+			return y;
+		else if (xe < ye)
+			return x;
+	}
+
+	/* Signs and exponents of inputs are equal, let's compare mantissas */
+	if (xs == 0) {
+		/* Inputs are both positive, with equal signs and exponents */
+		if (xm <= ym)
+			return y;
+		return x;
+	}
+	/* Inputs are both negative, with equal signs and exponents */
+	if (xm <= ym)
+		return x;
+	return y;
+}
+
+union ieee754dp ieee754dp_fmaxa(union ieee754dp x, union ieee754dp y)
+{
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+
+	FLUSHXDP;
+	FLUSHYDP;
+
+	ieee754_clearcx();
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754dp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754dp_nanxcpt(x);
+
+	/*
+	 * Quiet NaN handling
+	 */
+
+	/*
+	 *    The case of both inputs quiet NaNs
+	 */
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+		return x;
+
+	/*
+	 *    The cases of exactly one input quiet NaN (numbers
+	 *    are here preferred as returned values to NaNs)
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return y;
+
+	/*
+	 * Infinity and zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		return ieee754dp_inf(xs & ys);
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		return ieee754dp_zero(xs & ys);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		DPDNORMX;
+	}
+
+	/* Finally get to do some computation */
+
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+	/* Compare exponent */
+	if (xe > ye)
+		return x;
+	else if (xe < ye)
+		return y;
+
+	/* Compare mantissa */
+	if (xm < ym)
+		return y;
+	else if (xm > ym)
+		return x;
+	else if (xs == 0)
+		return x;
+	return y;
+}

linux-7.2.3/arch/mips/math-emu/dp_fmin.c

@@ -0,0 +1,252 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * IEEE754 floating point arithmetic
+ * double precision: MIN{,A}.f
+ * MIN : Scalar Floating-Point Minimum
+ * MINA: Scalar Floating-Point argument with Minimum Absolute Value
+ *
+ * MIN.D : FPR[fd] = minNum(FPR[fs],FPR[ft])
+ * MINA.D: FPR[fd] = maxNumMag(FPR[fs],FPR[ft])
+ *
+ * MIPS floating point support
+ * Copyright (C) 2015 Imagination Technologies, Ltd.
+ * Author: Markos Chandras <markos.chandras@imgtec.com>
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_fmin(union ieee754dp x, union ieee754dp y)
+{
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+
+	FLUSHXDP;
+	FLUSHYDP;
+
+	ieee754_clearcx();
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754dp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754dp_nanxcpt(x);
+
+	/*
+	 * Quiet NaN handling
+	 */
+
+	/*
+	 *    The case of both inputs quiet NaNs
+	 */
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+		return x;
+
+	/*
+	 *    The cases of exactly one input quiet NaN (numbers
+	 *    are here preferred as returned values to NaNs)
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return y;
+
+	/*
+	 * Infinity and zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return xs ? x : y;
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return ys ? y : x;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		return ieee754dp_zero(xs | ys);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		DPDNORMX;
+	}
+
+	/* Finally get to do some computation */
+
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+	/* Compare signs */
+	if (xs > ys)
+		return x;
+	else if (xs < ys)
+		return y;
+
+	/* Signs of inputs are the same, let's compare exponents */
+	if (xs == 0) {
+		/* Inputs are both positive */
+		if (xe > ye)
+			return y;
+		else if (xe < ye)
+			return x;
+	} else {
+		/* Inputs are both negative */
+		if (xe > ye)
+			return x;
+		else if (xe < ye)
+			return y;
+	}
+
+	/* Signs and exponents of inputs are equal, let's compare mantissas */
+	if (xs == 0) {
+		/* Inputs are both positive, with equal signs and exponents */
+		if (xm <= ym)
+			return x;
+		return y;
+	}
+	/* Inputs are both negative, with equal signs and exponents */
+	if (xm <= ym)
+		return y;
+	return x;
+}
+
+union ieee754dp ieee754dp_fmina(union ieee754dp x, union ieee754dp y)
+{
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+
+	FLUSHXDP;
+	FLUSHYDP;
+
+	ieee754_clearcx();
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754dp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754dp_nanxcpt(x);
+
+	/*
+	 * Quiet NaN handling
+	 */
+
+	/*
+	 *    The case of both inputs quiet NaNs
+	 */
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+		return x;
+
+	/*
+	 *    The cases of exactly one input quiet NaN (numbers
+	 *    are here preferred as returned values to NaNs)
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return y;
+
+	/*
+	 * Infinity and zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		return ieee754dp_inf(xs | ys);
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		return ieee754dp_zero(xs | ys);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		DPDNORMX;
+	}
+
+	/* Finally get to do some computation */
+
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+	/* Compare exponent */
+	if (xe > ye)
+		return y;
+	else if (xe < ye)
+		return x;
+
+	/* Compare mantissa */
+	if (xm < ym)
+		return x;
+	else if (xm > ym)
+		return y;
+	else if (xs == 1)
+		return x;
+	return y;
+}

linux-7.2.3/arch/mips/math-emu/dp_fsp.c

@@ -0,0 +1,63 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+#include "ieee754dp.h"
+
+static inline union ieee754dp ieee754dp_nan_fsp(int xs, u64 xm)
+{
+	return builddp(xs, DP_EMAX + 1 + DP_EBIAS,
+		       xm << (DP_FBITS - SP_FBITS));
+}
+
+union ieee754dp ieee754dp_fsp(union ieee754sp x)
+{
+	COMPXSP;
+
+	EXPLODEXSP;
+
+	ieee754_clearcx();
+
+	FLUSHXSP;
+
+	switch (xc) {
+	case IEEE754_CLASS_SNAN:
+		return ieee754dp_nanxcpt(ieee754dp_nan_fsp(xs, xm));
+
+	case IEEE754_CLASS_QNAN:
+		return ieee754dp_nan_fsp(xs, xm);
+
+	case IEEE754_CLASS_INF:
+		return ieee754dp_inf(xs);
+
+	case IEEE754_CLASS_ZERO:
+		return ieee754dp_zero(xs);
+
+	case IEEE754_CLASS_DNORM:
+		/* normalize */
+		while ((xm >> SP_FBITS) == 0) {
+			xm <<= 1;
+			xe--;
+		}
+		break;
+
+	case IEEE754_CLASS_NORM:
+		break;
+	}
+
+	/*
+	 * Can't possibly overflow,underflow, or need rounding
+	 */
+
+	/* drop the hidden bit */
+	xm &= ~SP_HIDDEN_BIT;
+
+	return builddp(xs, xe + DP_EBIAS,
+		       (u64) xm << (DP_FBITS - SP_FBITS));
+}

linux-7.2.3/arch/mips/math-emu/dp_maddf.c

@@ -0,0 +1,358 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * IEEE754 floating point arithmetic
+ * double precision: MADDF.f (Fused Multiply Add)
+ * MADDF.fmt: FPR[fd] = FPR[fd] + (FPR[fs] x FPR[ft])
+ *
+ * MIPS floating point support
+ * Copyright (C) 2015 Imagination Technologies, Ltd.
+ * Author: Markos Chandras <markos.chandras@imgtec.com>
+ */
+
+#include "ieee754dp.h"
+
+
+/* 128 bits shift right logical with rounding. */
+static void srl128(u64 *hptr, u64 *lptr, int count)
+{
+	u64 low;
+
+	if (count >= 128) {
+		*lptr = *hptr != 0 || *lptr != 0;
+		*hptr = 0;
+	} else if (count >= 64) {
+		if (count == 64) {
+			*lptr = *hptr | (*lptr != 0);
+		} else {
+			low = *lptr;
+			*lptr = *hptr >> (count - 64);
+			*lptr |= (*hptr << (128 - count)) != 0 || low != 0;
+		}
+		*hptr = 0;
+	} else {
+		low = *lptr;
+		*lptr = low >> count | *hptr << (64 - count);
+		*lptr |= (low << (64 - count)) != 0;
+		*hptr = *hptr >> count;
+	}
+}
+
+static union ieee754dp _dp_maddf(union ieee754dp z, union ieee754dp x,
+				 union ieee754dp y, enum maddf_flags flags)
+{
+	int re;
+	int rs;
+	unsigned int lxm;
+	unsigned int hxm;
+	unsigned int lym;
+	unsigned int hym;
+	u64 lrm;
+	u64 hrm;
+	u64 lzm;
+	u64 hzm;
+	u64 t;
+	u64 at;
+	int s;
+
+	COMPXDP;
+	COMPYDP;
+	COMPZDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+	EXPLODEZDP;
+
+	FLUSHXDP;
+	FLUSHYDP;
+	FLUSHZDP;
+
+	ieee754_clearcx();
+
+	rs = xs ^ ys;
+	if (flags & MADDF_NEGATE_PRODUCT)
+		rs ^= 1;
+	if (flags & MADDF_NEGATE_ADDITION)
+		zs ^= 1;
+
+	/*
+	 * Handle the cases when at least one of x, y or z is a NaN.
+	 * Order of precedence is sNaN, qNaN and z, x, y.
+	 */
+	if (zc == IEEE754_CLASS_SNAN)
+		return ieee754dp_nanxcpt(z);
+	if (xc == IEEE754_CLASS_SNAN)
+		return ieee754dp_nanxcpt(x);
+	if (yc == IEEE754_CLASS_SNAN)
+		return ieee754dp_nanxcpt(y);
+	if (zc == IEEE754_CLASS_QNAN)
+		return z;
+	if (xc == IEEE754_CLASS_QNAN)
+		return x;
+	if (yc == IEEE754_CLASS_QNAN)
+		return y;
+
+	if (zc == IEEE754_CLASS_DNORM)
+		DPDNORMZ;
+	/* ZERO z cases are handled separately below */
+
+	switch (CLPAIR(xc, yc)) {
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754dp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		if ((zc == IEEE754_CLASS_INF) && (zs != rs)) {
+			/*
+			 * Cases of addition of infinities with opposite signs
+			 * or subtraction of infinities with same signs.
+			 */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754dp_indef();
+		}
+		/*
+		 * z is here either not an infinity, or an infinity having the
+		 * same sign as product (x*y). The result must be an infinity,
+		 * and its sign is determined only by the sign of product (x*y).
+		 */
+		return ieee754dp_inf(rs);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		if (zc == IEEE754_CLASS_INF)
+			return ieee754dp_inf(zs);
+		if (zc == IEEE754_CLASS_ZERO) {
+			/* Handle cases +0 + (-0) and similar ones. */
+			if (zs == rs)
+				/*
+				 * Cases of addition of zeros of equal signs
+				 * or subtraction of zeroes of opposite signs.
+				 * The sign of the resulting zero is in any
+				 * such case determined only by the sign of z.
+				 */
+				return z;
+
+			return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
+		}
+		/* x*y is here 0, and z is not 0, so just return z */
+		return z;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		if (zc == IEEE754_CLASS_INF)
+			return ieee754dp_inf(zs);
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		if (zc == IEEE754_CLASS_INF)
+			return ieee754dp_inf(zs);
+		DPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		if (zc == IEEE754_CLASS_INF)
+			return ieee754dp_inf(zs);
+		/* continue to real computations */
+	}
+
+	/* Finally get to do some computation */
+
+	/*
+	 * Do the multiplication bit first
+	 *
+	 * rm = xm * ym, re = xe + ye basically
+	 *
+	 * At this point xm and ym should have been normalized.
+	 */
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+	re = xe + ye;
+
+	/* shunt to top of word */
+	xm <<= 64 - (DP_FBITS + 1);
+	ym <<= 64 - (DP_FBITS + 1);
+
+	/*
+	 * Multiply 64 bits xm and ym to give 128 bits result in hrm:lrm.
+	 */
+
+	lxm = xm;
+	hxm = xm >> 32;
+	lym = ym;
+	hym = ym >> 32;
+
+	lrm = DPXMULT(lxm, lym);
+	hrm = DPXMULT(hxm, hym);
+
+	t = DPXMULT(lxm, hym);
+
+	at = lrm + (t << 32);
+	hrm += at < lrm;
+	lrm = at;
+
+	hrm = hrm + (t >> 32);
+
+	t = DPXMULT(hxm, lym);
+
+	at = lrm + (t << 32);
+	hrm += at < lrm;
+	lrm = at;
+
+	hrm = hrm + (t >> 32);
+
+	/* Put explicit bit at bit 126 if necessary */
+	if ((int64_t)hrm < 0) {
+		lrm = (hrm << 63) | (lrm >> 1);
+		hrm = hrm >> 1;
+		re++;
+	}
+
+	assert(hrm & (1 << 62));
+
+	if (zc == IEEE754_CLASS_ZERO) {
+		/*
+		 * Move explicit bit from bit 126 to bit 55 since the
+		 * ieee754dp_format code expects the mantissa to be
+		 * 56 bits wide (53 + 3 rounding bits).
+		 */
+		srl128(&hrm, &lrm, (126 - 55));
+		return ieee754dp_format(rs, re, lrm);
+	}
+
+	/* Move explicit bit from bit 52 to bit 126 */
+	lzm = 0;
+	hzm = zm << 10;
+	assert(hzm & (1 << 62));
+
+	/* Make the exponents the same */
+	if (ze > re) {
+		/*
+		 * Have to shift y fraction right to align.
+		 */
+		s = ze - re;
+		srl128(&hrm, &lrm, s);
+		re += s;
+	} else if (re > ze) {
+		/*
+		 * Have to shift x fraction right to align.
+		 */
+		s = re - ze;
+		srl128(&hzm, &lzm, s);
+		ze += s;
+	}
+	assert(ze == re);
+	assert(ze <= DP_EMAX);
+
+	/* Do the addition */
+	if (zs == rs) {
+		/*
+		 * Generate 128 bit result by adding two 127 bit numbers
+		 * leaving result in hzm:lzm, zs and ze.
+		 */
+		hzm = hzm + hrm + (lzm > (lzm + lrm));
+		lzm = lzm + lrm;
+		if ((int64_t)hzm < 0) {        /* carry out */
+			srl128(&hzm, &lzm, 1);
+			ze++;
+		}
+	} else {
+		if (hzm > hrm || (hzm == hrm && lzm >= lrm)) {
+			hzm = hzm - hrm - (lzm < lrm);
+			lzm = lzm - lrm;
+		} else {
+			hzm = hrm - hzm - (lrm < lzm);
+			lzm = lrm - lzm;
+			zs = rs;
+		}
+		if (lzm == 0 && hzm == 0)
+			return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
+
+		/*
+		 * Put explicit bit at bit 126 if necessary.
+		 */
+		if (hzm == 0) {
+			/* left shift by 63 or 64 bits */
+			if ((int64_t)lzm < 0) {
+				/* MSB of lzm is the explicit bit */
+				hzm = lzm >> 1;
+				lzm = lzm << 63;
+				ze -= 63;
+			} else {
+				hzm = lzm;
+				lzm = 0;
+				ze -= 64;
+			}
+		}
+
+		t = 0;
+		while ((hzm >> (62 - t)) == 0)
+			t++;
+
+		assert(t <= 62);
+		if (t) {
+			hzm = hzm << t | lzm >> (64 - t);
+			lzm = lzm << t;
+			ze -= t;
+		}
+	}
+
+	/*
+	 * Move explicit bit from bit 126 to bit 55 since the
+	 * ieee754dp_format code expects the mantissa to be
+	 * 56 bits wide (53 + 3 rounding bits).
+	 */
+	srl128(&hzm, &lzm, (126 - 55));
+
+	return ieee754dp_format(zs, ze, lzm);
+}
+
+union ieee754dp ieee754dp_maddf(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y)
+{
+	return _dp_maddf(z, x, y, 0);
+}
+
+union ieee754dp ieee754dp_msubf(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y)
+{
+	return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT);
+}
+
+union ieee754dp ieee754dp_madd(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y)
+{
+	return _dp_maddf(z, x, y, 0);
+}
+
+union ieee754dp ieee754dp_msub(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y)
+{
+	return _dp_maddf(z, x, y, MADDF_NEGATE_ADDITION);
+}
+
+union ieee754dp ieee754dp_nmadd(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y)
+{
+	return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT|MADDF_NEGATE_ADDITION);
+}
+
+union ieee754dp ieee754dp_nmsub(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y)
+{
+	return _dp_maddf(z, x, y, MADDF_NEGATE_PRODUCT);
+}

linux-7.2.3/arch/mips/math-emu/dp_mul.c

@@ -0,0 +1,159 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_mul(union ieee754dp x, union ieee754dp y)
+{
+	int re;
+	int rs;
+	u64 rm;
+	unsigned int lxm;
+	unsigned int hxm;
+	unsigned int lym;
+	unsigned int hym;
+	u64 lrm;
+	u64 hrm;
+	u64 t;
+	u64 at;
+
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+
+	ieee754_clearcx();
+
+	FLUSHXDP;
+	FLUSHYDP;
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754dp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754dp_nanxcpt(x);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return x;
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754dp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		return ieee754dp_inf(xs ^ ys);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return ieee754dp_zero(xs ^ ys);
+
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		DPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		break;
+	}
+	/* rm = xm * ym, re = xe+ye basically */
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+	re = xe + ye;
+	rs = xs ^ ys;
+
+	/* shunt to top of word */
+	xm <<= 64 - (DP_FBITS + 1);
+	ym <<= 64 - (DP_FBITS + 1);
+
+	/*
+	 * Multiply 64 bits xm, ym to give high 64 bits rm with stickness.
+	 */
+
+	lxm = xm;
+	hxm = xm >> 32;
+	lym = ym;
+	hym = ym >> 32;
+
+	lrm = DPXMULT(lxm, lym);
+	hrm = DPXMULT(hxm, hym);
+
+	t = DPXMULT(lxm, hym);
+
+	at = lrm + (t << 32);
+	hrm += at < lrm;
+	lrm = at;
+
+	hrm = hrm + (t >> 32);
+
+	t = DPXMULT(hxm, lym);
+
+	at = lrm + (t << 32);
+	hrm += at < lrm;
+	lrm = at;
+
+	hrm = hrm + (t >> 32);
+
+	rm = hrm | (lrm != 0);
+
+	/*
+	 * Sticky shift down to normal rounding precision.
+	 */
+	if ((s64) rm < 0) {
+		rm = (rm >> (64 - (DP_FBITS + 1 + 3))) |
+		     ((rm << (DP_FBITS + 1 + 3)) != 0);
+		re++;
+	} else {
+		rm = (rm >> (64 - (DP_FBITS + 1 + 3 + 1))) |
+		     ((rm << (DP_FBITS + 1 + 3 + 1)) != 0);
+	}
+	assert(rm & (DP_HIDDEN_BIT << 3));
+
+	return ieee754dp_format(rs, re, rm);
+}

linux-7.2.3/arch/mips/math-emu/dp_rint.c

@@ -0,0 +1,78 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ * Copyright (C) 2017 Imagination Technologies, Ltd.
+ * Author: Aleksandar Markovic <aleksandar.markovic@imgtec.com>
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_rint(union ieee754dp x)
+{
+	union ieee754dp ret;
+	u64 residue;
+	int sticky;
+	int round;
+	int odd;
+
+	COMPXDP;
+
+	ieee754_clearcx();
+
+	EXPLODEXDP;
+	FLUSHXDP;
+
+	if (xc == IEEE754_CLASS_SNAN)
+		return ieee754dp_nanxcpt(x);
+
+	if ((xc == IEEE754_CLASS_QNAN) ||
+	    (xc == IEEE754_CLASS_INF) ||
+	    (xc == IEEE754_CLASS_ZERO))
+		return x;
+
+	if (xe >= DP_FBITS)
+		return x;
+
+	if (xe < -1) {
+		residue = xm;
+		round = 0;
+		sticky = residue != 0;
+		xm = 0;
+	} else {
+		residue = xm << (64 - DP_FBITS + xe);
+		round = (residue >> 63) != 0;
+		sticky = (residue << 1) != 0;
+		xm >>= DP_FBITS - xe;
+	}
+
+	odd = (xm & 0x1) != 0x0;
+
+	switch (ieee754_csr.rm) {
+	case FPU_CSR_RN:	/* toward nearest */
+		if (round && (sticky || odd))
+			xm++;
+		break;
+	case FPU_CSR_RZ:	/* toward zero */
+		break;
+	case FPU_CSR_RU:	/* toward +infinity */
+		if ((round || sticky) && !xs)
+			xm++;
+		break;
+	case FPU_CSR_RD:	/* toward -infinity */
+		if ((round || sticky) && xs)
+			xm++;
+		break;
+	}
+
+	if (round || sticky)
+		ieee754_setcx(IEEE754_INEXACT);
+
+	ret = ieee754dp_flong(xm);
+	DPSIGN(ret) = xs;
+
+	return ret;
+}

linux-7.2.3/arch/mips/math-emu/dp_simple.c

@@ -0,0 +1,49 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_neg(union ieee754dp x)
+{
+	union ieee754dp y;
+
+	if (ieee754_csr.abs2008) {
+		y = x;
+		DPSIGN(y) = !DPSIGN(x);
+	} else {
+		unsigned int oldrm;
+
+		oldrm = ieee754_csr.rm;
+		ieee754_csr.rm = FPU_CSR_RD;
+		y = ieee754dp_sub(ieee754dp_zero(0), x);
+		ieee754_csr.rm = oldrm;
+	}
+	return y;
+}
+
+union ieee754dp ieee754dp_abs(union ieee754dp x)
+{
+	union ieee754dp y;
+
+	if (ieee754_csr.abs2008) {
+		y = x;
+		DPSIGN(y) = 0;
+	} else {
+		unsigned int oldrm;
+
+		oldrm = ieee754_csr.rm;
+		ieee754_csr.rm = FPU_CSR_RD;
+		if (DPSIGN(x))
+			y = ieee754dp_sub(ieee754dp_zero(0), x);
+		else
+			y = ieee754dp_add(ieee754dp_zero(0), x);
+		ieee754_csr.rm = oldrm;
+	}
+	return y;
+}

linux-7.2.3/arch/mips/math-emu/dp_sqrt.c

@@ -0,0 +1,152 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision square root
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+static const unsigned int table[] = {
+	0, 1204, 3062, 5746, 9193, 13348, 18162, 23592,
+	29598, 36145, 43202, 50740, 58733, 67158, 75992,
+	85215, 83599, 71378, 60428, 50647, 41945, 34246,
+	27478, 21581, 16499, 12183, 8588, 5674, 3403,
+	1742, 661, 130
+};
+
+union ieee754dp ieee754dp_sqrt(union ieee754dp x)
+{
+	struct _ieee754_csr oldcsr;
+	union ieee754dp y, z, t;
+	unsigned int scalx, yh;
+	COMPXDP;
+
+	EXPLODEXDP;
+	ieee754_clearcx();
+	FLUSHXDP;
+
+	/* x == INF or NAN? */
+	switch (xc) {
+	case IEEE754_CLASS_SNAN:
+		return ieee754dp_nanxcpt(x);
+
+	case IEEE754_CLASS_QNAN:
+		/* sqrt(Nan) = Nan */
+		return x;
+
+	case IEEE754_CLASS_ZERO:
+		/* sqrt(0) = 0 */
+		return x;
+
+	case IEEE754_CLASS_INF:
+		if (xs) {
+			/* sqrt(-Inf) = Nan */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754dp_indef();
+		}
+		/* sqrt(+Inf) = Inf */
+		return x;
+
+	case IEEE754_CLASS_DNORM:
+		DPDNORMX;
+		fallthrough;
+	case IEEE754_CLASS_NORM:
+		if (xs) {
+			/* sqrt(-x) = Nan */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754dp_indef();
+		}
+		break;
+	}
+
+	/* save old csr; switch off INX enable & flag; set RN rounding */
+	oldcsr = ieee754_csr;
+	ieee754_csr.mx &= ~IEEE754_INEXACT;
+	ieee754_csr.sx &= ~IEEE754_INEXACT;
+	ieee754_csr.rm = FPU_CSR_RN;
+
+	/* adjust exponent to prevent overflow */
+	scalx = 0;
+	if (xe > 512) {		/* x > 2**-512? */
+		xe -= 512;	/* x = x / 2**512 */
+		scalx += 256;
+	} else if (xe < -512) { /* x < 2**-512? */
+		xe += 512;	/* x = x * 2**512 */
+		scalx -= 256;
+	}
+
+	x = builddp(0, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
+	y = x;
+
+	/* magic initial approximation to almost 8 sig. bits */
+	yh = y.bits >> 32;
+	yh = (yh >> 1) + 0x1ff80000;
+	yh = yh - table[(yh >> 15) & 31];
+	y.bits = ((u64) yh << 32) | (y.bits & 0xffffffff);
+
+	/* Heron's rule once with correction to improve to ~18 sig. bits */
+	/* t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; */
+	t = ieee754dp_div(x, y);
+	y = ieee754dp_add(y, t);
+	y.bits -= 0x0010000600000000LL;
+	y.bits &= 0xffffffff00000000LL;
+
+	/* triple to almost 56 sig. bits: y ~= sqrt(x) to within 1 ulp */
+	/* t=y*y; z=t;	pt[n0]+=0x00100000; t+=z; z=(x-z)*y; */
+	t = ieee754dp_mul(y, y);
+	z = t;
+	t.bexp += 0x001;
+	t = ieee754dp_add(t, z);
+	z = ieee754dp_mul(ieee754dp_sub(x, z), y);
+
+	/* t=z/(t+x) ;	pt[n0]+=0x00100000; y+=t; */
+	t = ieee754dp_div(z, ieee754dp_add(t, x));
+	t.bexp += 0x001;
+	y = ieee754dp_add(y, t);
+
+	/* twiddle last bit to force y correctly rounded */
+
+	/* set RZ, clear INEX flag */
+	ieee754_csr.rm = FPU_CSR_RZ;
+	ieee754_csr.sx &= ~IEEE754_INEXACT;
+
+	/* t=x/y; ...chopped quotient, possibly inexact */
+	t = ieee754dp_div(x, y);
+
+	if (ieee754_csr.sx & IEEE754_INEXACT || t.bits != y.bits) {
+
+		if (!(ieee754_csr.sx & IEEE754_INEXACT))
+			/* t = t-ulp */
+			t.bits -= 1;
+
+		/* add inexact to result status */
+		oldcsr.cx |= IEEE754_INEXACT;
+		oldcsr.sx |= IEEE754_INEXACT;
+
+		switch (oldcsr.rm) {
+		case FPU_CSR_RU:
+			y.bits += 1;
+			fallthrough;
+		case FPU_CSR_RN:
+			t.bits += 1;
+			break;
+		}
+
+		/* y=y+t; ...chopped sum */
+		y = ieee754dp_add(y, t);
+
+		/* adjust scalx for correctly rounded sqrt(x) */
+		scalx -= 1;
+	}
+
+	/* py[n0]=py[n0]+scalx; ...scale back y */
+	y.bexp += scalx;
+
+	/* restore rounding mode, possibly set inexact */
+	ieee754_csr = oldcsr;
+
+	return y;
+}

linux-7.2.3/arch/mips/math-emu/dp_sub.c

@@ -0,0 +1,172 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_sub(union ieee754dp x, union ieee754dp y)
+{
+	int s;
+
+	COMPXDP;
+	COMPYDP;
+
+	EXPLODEXDP;
+	EXPLODEYDP;
+
+	ieee754_clearcx();
+
+	FLUSHXDP;
+	FLUSHYDP;
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754dp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754dp_nanxcpt(x);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return x;
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		if (xs != ys)
+			return x;
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754dp_indef();
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+		return ieee754dp_inf(ys ^ 1);
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+		return x;
+
+	/*
+	 * Zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		if (xs != ys)
+			return x;
+		else
+			return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		/* quick fix up */
+		DPSIGN(y) ^= 1;
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		DPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		/* normalize ym,ye */
+		DPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		/* normalize xm,xe */
+		DPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		break;
+	}
+	/* flip sign of y and handle as add */
+	ys ^= 1;
+
+	assert(xm & DP_HIDDEN_BIT);
+	assert(ym & DP_HIDDEN_BIT);
+
+
+	/* provide guard,round and stick bit dpace */
+	xm <<= 3;
+	ym <<= 3;
+
+	if (xe > ye) {
+		/*
+		 * Have to shift y fraction right to align
+		 */
+		s = xe - ye;
+		ym = XDPSRS(ym, s);
+		ye += s;
+	} else if (ye > xe) {
+		/*
+		 * Have to shift x fraction right to align
+		 */
+		s = ye - xe;
+		xm = XDPSRS(xm, s);
+		xe += s;
+	}
+	assert(xe == ye);
+	assert(xe <= DP_EMAX);
+
+	if (xs == ys) {
+		/* generate 28 bit result of adding two 27 bit numbers
+		 */
+		xm = xm + ym;
+
+		if (xm >> (DP_FBITS + 1 + 3)) { /* carry out */
+			xm = XDPSRS1(xm);	/* shift preserving sticky */
+			xe++;
+		}
+	} else {
+		if (xm >= ym) {
+			xm = xm - ym;
+		} else {
+			xm = ym - xm;
+			xs = ys;
+		}
+		if (xm == 0) {
+			if (ieee754_csr.rm == FPU_CSR_RD)
+				return ieee754dp_zero(1);	/* round negative inf. => sign = -1 */
+			else
+				return ieee754dp_zero(0);	/* other round modes   => sign = 1 */
+		}
+
+		/* normalize to rounding precision
+		 */
+		while ((xm >> (DP_FBITS + 3)) == 0) {
+			xm <<= 1;
+			xe--;
+		}
+	}
+
+	return ieee754dp_format(xs, xe, xm);
+}

linux-7.2.3/arch/mips/math-emu/dp_tint.c

@@ -0,0 +1,96 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+int ieee754dp_tint(union ieee754dp x)
+{
+	u64 residue;
+	int round;
+	int sticky;
+	int odd;
+
+	COMPXDP;
+
+	ieee754_clearcx();
+
+	EXPLODEXDP;
+	FLUSHXDP;
+
+	switch (xc) {
+	case IEEE754_CLASS_SNAN:
+	case IEEE754_CLASS_QNAN:
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754si_indef();
+
+	case IEEE754_CLASS_INF:
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754si_overflow(xs);
+
+	case IEEE754_CLASS_ZERO:
+		return 0;
+
+	case IEEE754_CLASS_DNORM:
+	case IEEE754_CLASS_NORM:
+		break;
+	}
+	if (xe > 31) {
+		/* Set invalid. We will only use overflow for floating
+		   point overflow */
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754si_overflow(xs);
+	}
+	/* oh gawd */
+	if (xe > DP_FBITS) {
+		xm <<= xe - DP_FBITS;
+	} else if (xe < DP_FBITS) {
+		if (xe < -1) {
+			residue = xm;
+			round = 0;
+			sticky = residue != 0;
+			xm = 0;
+		} else {
+			residue = xm << (64 - DP_FBITS + xe);
+			round = (residue >> 63) != 0;
+			sticky = (residue << 1) != 0;
+			xm >>= DP_FBITS - xe;
+		}
+		/* Note: At this point upper 32 bits of xm are guaranteed
+		   to be zero */
+		odd = (xm & 0x1) != 0x0;
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RN:
+			if (round && (sticky || odd))
+				xm++;
+			break;
+		case FPU_CSR_RZ:
+			break;
+		case FPU_CSR_RU:	/* toward +Infinity */
+			if ((round || sticky) && !xs)
+				xm++;
+			break;
+		case FPU_CSR_RD:	/* toward -Infinity */
+			if ((round || sticky) && xs)
+				xm++;
+			break;
+		}
+		/* look for valid corner case 0x80000000 */
+		if ((xm >> 31) != 0 && (xs == 0 || xm != 0x80000000)) {
+			/* This can happen after rounding */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754si_overflow(xs);
+		}
+		if (round || sticky)
+			ieee754_setcx(IEEE754_INEXACT);
+	}
+	if (xs)
+		return -xm;
+	else
+		return xm;
+}

linux-7.2.3/arch/mips/math-emu/dp_tlong.c

@@ -0,0 +1,100 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754dp.h"
+
+s64 ieee754dp_tlong(union ieee754dp x)
+{
+	u64 residue;
+	int round;
+	int sticky;
+	int odd;
+
+	COMPXDP;
+
+	ieee754_clearcx();
+
+	EXPLODEXDP;
+	FLUSHXDP;
+
+	switch (xc) {
+	case IEEE754_CLASS_SNAN:
+	case IEEE754_CLASS_QNAN:
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754di_indef();
+
+	case IEEE754_CLASS_INF:
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754di_overflow(xs);
+
+	case IEEE754_CLASS_ZERO:
+		return 0;
+
+	case IEEE754_CLASS_DNORM:
+	case IEEE754_CLASS_NORM:
+		break;
+	}
+	if (xe >= 63) {
+		/* look for valid corner case */
+		if (xe == 63 && xs && xm == DP_HIDDEN_BIT)
+			return -0x8000000000000000LL;
+		/* Set invalid. We will only use overflow for floating
+		   point overflow */
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754di_overflow(xs);
+	}
+	/* oh gawd */
+	if (xe > DP_FBITS) {
+		xm <<= xe - DP_FBITS;
+	} else if (xe < DP_FBITS) {
+		if (xe < -1) {
+			residue = xm;
+			round = 0;
+			sticky = residue != 0;
+			xm = 0;
+		} else {
+			/* Shifting a u64 64 times does not work,
+			* so we do it in two steps. Be aware that xe
+			* may be -1 */
+			residue = xm << (xe + 1);
+			residue <<= 63 - DP_FBITS;
+			round = (residue >> 63) != 0;
+			sticky = (residue << 1) != 0;
+			xm >>= DP_FBITS - xe;
+		}
+		odd = (xm & 0x1) != 0x0;
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RN:
+			if (round && (sticky || odd))
+				xm++;
+			break;
+		case FPU_CSR_RZ:
+			break;
+		case FPU_CSR_RU:	/* toward +Infinity */
+			if ((round || sticky) && !xs)
+				xm++;
+			break;
+		case FPU_CSR_RD:	/* toward -Infinity */
+			if ((round || sticky) && xs)
+				xm++;
+			break;
+		}
+		if ((xm >> 63) != 0) {
+			/* This can happen after rounding */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754di_overflow(xs);
+		}
+		if (round || sticky)
+			ieee754_setcx(IEEE754_INEXACT);
+	}
+	if (xs)
+		return -xm;
+	else
+		return xm;
+}

linux-7.2.3/arch/mips/math-emu/dsemul.c

@@ -0,0 +1,305 @@
+// SPDX-License-Identifier: GPL-2.0
+#include <linux/err.h>
+#include <linux/slab.h>
+#include <linux/mm_types.h>
+#include <linux/sched/task.h>
+
+#include <asm/branch.h>
+#include <asm/cacheflush.h>
+#include <asm/fpu_emulator.h>
+#include <asm/inst.h>
+#include <asm/mipsregs.h>
+#include <linux/uaccess.h>
+
+/**
+ * struct emuframe - The 'emulation' frame structure
+ * @emul:	The instruction to 'emulate'.
+ * @badinst:	A break instruction to cause a return to the kernel.
+ *
+ * This structure defines the frames placed within the delay slot emulation
+ * page in response to a call to mips_dsemul(). Each thread may be allocated
+ * only one frame at any given time. The kernel stores within it the
+ * instruction to be 'emulated' followed by a break instruction, then
+ * executes the frame in user mode. The break causes a trap to the kernel
+ * which leads to do_dsemulret() being called unless the instruction in
+ * @emul causes a trap itself, is a branch, or a signal is delivered to
+ * the thread. In these cases the allocated frame will either be reused by
+ * a subsequent delay slot 'emulation', or be freed during signal delivery or
+ * upon thread exit.
+ *
+ * This approach is used because:
+ *
+ * - Actually emulating all instructions isn't feasible. We would need to
+ *   be able to handle instructions from all revisions of the MIPS ISA,
+ *   all ASEs & all vendor instruction set extensions. This would be a
+ *   whole lot of work & continual maintenance burden as new instructions
+ *   are introduced, and in the case of some vendor extensions may not
+ *   even be possible. Thus we need to take the approach of actually
+ *   executing the instruction.
+ *
+ * - We must execute the instruction within user context. If we were to
+ *   execute the instruction in kernel mode then it would have access to
+ *   kernel resources without very careful checks, leaving us with a
+ *   high potential for security or stability issues to arise.
+ *
+ * - We used to place the frame on the users stack, but this requires
+ *   that the stack be executable. This is bad for security so the
+ *   per-process page is now used instead.
+ *
+ * - The instruction in @emul may be something entirely invalid for a
+ *   delay slot. The user may (intentionally or otherwise) place a branch
+ *   in a delay slot, or a kernel mode instruction, or something else
+ *   which generates an exception. Thus we can't rely upon the break in
+ *   @badinst always being hit. For this reason we track the index of the
+ *   frame allocated to each thread, allowing us to clean it up at later
+ *   points such as signal delivery or thread exit.
+ *
+ * - The user may generate a fake struct emuframe if they wish, invoking
+ *   the BRK_MEMU break instruction themselves. We must therefore not
+ *   trust that BRK_MEMU means there's actually a valid frame allocated
+ *   to the thread, and must not allow the user to do anything they
+ *   couldn't already.
+ */
+struct emuframe {
+	mips_instruction	emul;
+	mips_instruction	badinst;
+};
+
+static const int emupage_frame_count = PAGE_SIZE / sizeof(struct emuframe);
+
+static inline __user struct emuframe *dsemul_page(void)
+{
+	return (__user struct emuframe *)STACK_TOP;
+}
+
+static int alloc_emuframe(void)
+{
+	mm_context_t *mm_ctx = &current->mm->context;
+	int idx;
+
+retry:
+	spin_lock(&mm_ctx->bd_emupage_lock);
+
+	/* Ensure we have an allocation bitmap */
+	if (!mm_ctx->bd_emupage_allocmap) {
+		mm_ctx->bd_emupage_allocmap = bitmap_zalloc(emupage_frame_count,
+							    GFP_ATOMIC);
+		if (!mm_ctx->bd_emupage_allocmap) {
+			idx = BD_EMUFRAME_NONE;
+			goto out_unlock;
+		}
+	}
+
+	/* Attempt to allocate a single bit/frame */
+	idx = bitmap_find_free_region(mm_ctx->bd_emupage_allocmap,
+				      emupage_frame_count, 0);
+	if (idx < 0) {
+		/*
+		 * Failed to allocate a frame. We'll wait until one becomes
+		 * available. We unlock the page so that other threads actually
+		 * get the opportunity to free their frames, which means
+		 * technically the result of bitmap_full may be incorrect.
+		 * However the worst case is that we repeat all this and end up
+		 * back here again.
+		 */
+		spin_unlock(&mm_ctx->bd_emupage_lock);
+		if (!wait_event_killable(mm_ctx->bd_emupage_queue,
+			!bitmap_full(mm_ctx->bd_emupage_allocmap,
+				     emupage_frame_count)))
+			goto retry;
+
+		/* Received a fatal signal - just give in */
+		return BD_EMUFRAME_NONE;
+	}
+
+	/* Success! */
+	pr_debug("allocate emuframe %d to %d\n", idx, current->pid);
+out_unlock:
+	spin_unlock(&mm_ctx->bd_emupage_lock);
+	return idx;
+}
+
+static void free_emuframe(int idx, struct mm_struct *mm)
+{
+	mm_context_t *mm_ctx = &mm->context;
+
+	spin_lock(&mm_ctx->bd_emupage_lock);
+
+	pr_debug("free emuframe %d from %d\n", idx, current->pid);
+	bitmap_clear(mm_ctx->bd_emupage_allocmap, idx, 1);
+
+	/* If some thread is waiting for a frame, now's its chance */
+	wake_up(&mm_ctx->bd_emupage_queue);
+
+	spin_unlock(&mm_ctx->bd_emupage_lock);
+}
+
+static bool within_emuframe(struct pt_regs *regs)
+{
+	unsigned long base = (unsigned long)dsemul_page();
+
+	if (regs->cp0_epc < base)
+		return false;
+	if (regs->cp0_epc >= (base + PAGE_SIZE))
+		return false;
+
+	return true;
+}
+
+bool dsemul_thread_cleanup(struct task_struct *tsk)
+{
+	int fr_idx;
+
+	/* Clear any allocated frame, retrieving its index */
+	fr_idx = atomic_xchg(&tsk->thread.bd_emu_frame, BD_EMUFRAME_NONE);
+
+	/* If no frame was allocated, we're done */
+	if (fr_idx == BD_EMUFRAME_NONE)
+		return false;
+
+	task_lock(tsk);
+
+	/* Free the frame that this thread had allocated */
+	if (tsk->mm)
+		free_emuframe(fr_idx, tsk->mm);
+
+	task_unlock(tsk);
+	return true;
+}
+
+bool dsemul_thread_rollback(struct pt_regs *regs)
+{
+	struct emuframe __user *fr;
+	int fr_idx;
+
+	/* Do nothing if we're not executing from a frame */
+	if (!within_emuframe(regs))
+		return false;
+
+	/* Find the frame being executed */
+	fr_idx = atomic_read(&current->thread.bd_emu_frame);
+	if (fr_idx == BD_EMUFRAME_NONE)
+		return false;
+	fr = &dsemul_page()[fr_idx];
+
+	/*
+	 * If the PC is at the emul instruction, roll back to the branch. If
+	 * PC is at the badinst (break) instruction, we've already emulated the
+	 * instruction so progress to the continue PC. If it's anything else
+	 * then something is amiss & the user has branched into some other area
+	 * of the emupage - we'll free the allocated frame anyway.
+	 */
+	if (msk_isa16_mode(regs->cp0_epc) == (unsigned long)&fr->emul)
+		regs->cp0_epc = current->thread.bd_emu_branch_pc;
+	else if (msk_isa16_mode(regs->cp0_epc) == (unsigned long)&fr->badinst)
+		regs->cp0_epc = current->thread.bd_emu_cont_pc;
+
+	atomic_set(&current->thread.bd_emu_frame, BD_EMUFRAME_NONE);
+	free_emuframe(fr_idx, current->mm);
+	return true;
+}
+
+void dsemul_mm_cleanup(struct mm_struct *mm)
+{
+	mm_context_t *mm_ctx = &mm->context;
+
+	bitmap_free(mm_ctx->bd_emupage_allocmap);
+}
+
+int mips_dsemul(struct pt_regs *regs, mips_instruction ir,
+		unsigned long branch_pc, unsigned long cont_pc)
+{
+	int isa16 = get_isa16_mode(regs->cp0_epc);
+	mips_instruction break_math;
+	unsigned long fr_uaddr;
+	struct emuframe fr;
+	int fr_idx, ret;
+
+	/* NOP is easy */
+	if (ir == 0)
+		return -1;
+
+	/* microMIPS instructions */
+	if (isa16) {
+		union mips_instruction insn = { .word = ir };
+
+		/* NOP16 aka MOVE16 $0, $0 */
+		if ((ir >> 16) == MM_NOP16)
+			return -1;
+
+		/* ADDIUPC */
+		if (insn.mm_a_format.opcode == mm_addiupc_op) {
+			unsigned int rs;
+			s32 v;
+
+			rs = (((insn.mm_a_format.rs + 0xe) & 0xf) + 2);
+			v = regs->cp0_epc & ~3;
+			v += insn.mm_a_format.simmediate << 2;
+			regs->regs[rs] = (long)v;
+			return -1;
+		}
+	}
+
+	pr_debug("dsemul 0x%08lx cont at 0x%08lx\n", regs->cp0_epc, cont_pc);
+
+	/* Allocate a frame if we don't already have one */
+	fr_idx = atomic_read(&current->thread.bd_emu_frame);
+	if (fr_idx == BD_EMUFRAME_NONE)
+		fr_idx = alloc_emuframe();
+	if (fr_idx == BD_EMUFRAME_NONE)
+		return SIGBUS;
+
+	/* Retrieve the appropriately encoded break instruction */
+	break_math = BREAK_MATH(isa16);
+
+	/* Write the instructions to the frame */
+	if (isa16) {
+		union mips_instruction _emul = {
+			.halfword = { ir >> 16, ir }
+		};
+		union mips_instruction _badinst = {
+			.halfword = { break_math >> 16, break_math }
+		};
+
+		fr.emul = _emul.word;
+		fr.badinst = _badinst.word;
+	} else {
+		fr.emul = ir;
+		fr.badinst = break_math;
+	}
+
+	/* Write the frame to user memory */
+	fr_uaddr = (unsigned long)&dsemul_page()[fr_idx];
+	ret = access_process_vm(current, fr_uaddr, &fr, sizeof(fr),
+				FOLL_FORCE | FOLL_WRITE);
+	if (unlikely(ret != sizeof(fr))) {
+		MIPS_FPU_EMU_INC_STATS(errors);
+		free_emuframe(fr_idx, current->mm);
+		return SIGBUS;
+	}
+
+	/* Record the PC of the branch, PC to continue from & frame index */
+	current->thread.bd_emu_branch_pc = branch_pc;
+	current->thread.bd_emu_cont_pc = cont_pc;
+	atomic_set(&current->thread.bd_emu_frame, fr_idx);
+
+	/* Change user register context to execute the frame */
+	regs->cp0_epc = fr_uaddr | isa16;
+
+	return 0;
+}
+
+bool do_dsemulret(struct pt_regs *xcp)
+{
+	/* Cleanup the allocated frame, returning if there wasn't one */
+	if (!dsemul_thread_cleanup(current)) {
+		MIPS_FPU_EMU_INC_STATS(errors);
+		return false;
+	}
+
+	/* Set EPC to return to post-branch instruction */
+	xcp->cp0_epc = current->thread.bd_emu_cont_pc;
+	pr_debug("dsemulret to 0x%08lx\n", xcp->cp0_epc);
+	MIPS_FPU_EMU_INC_STATS(ds_emul);
+	return true;
+}

linux-7.2.3/arch/mips/math-emu/ieee754.c

@@ -0,0 +1,83 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* ieee754 floating point arithmetic
+ * single and double precision
+ *
+ * BUGS
+ * not much dp done
+ * doesn't generate IEEE754_INEXACT
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include <linux/compiler.h>
+
+#include "ieee754.h"
+#include "ieee754sp.h"
+#include "ieee754dp.h"
+
+/*
+ * Special constants
+ */
+
+/*
+ * Older GCC requires the inner braces for initialization of union ieee754dp's
+ * anonymous struct member.  Without an error will result.
+ */
+#define xPCNST(s, b, m, ebias)						\
+{									\
+	{								\
+		.sign	= (s),						\
+		.bexp	= (b) + ebias,					\
+		.mant	= (m)						\
+	}								\
+}
+
+#define DPCNST(s, b, m)							\
+	xPCNST(s, b, m, DP_EBIAS)
+
+const union ieee754dp __ieee754dp_spcvals[] = {
+	DPCNST(0, DP_EMIN - 1, 0x0000000000000ULL),	/* + zero   */
+	DPCNST(1, DP_EMIN - 1, 0x0000000000000ULL),	/* - zero   */
+	DPCNST(0, 0,	       0x0000000000000ULL),	/* + 1.0   */
+	DPCNST(1, 0,	       0x0000000000000ULL),	/* - 1.0   */
+	DPCNST(0, 3,           0x4000000000000ULL),	/* + 10.0   */
+	DPCNST(1, 3,           0x4000000000000ULL),	/* - 10.0   */
+	DPCNST(0, DP_EMAX + 1, 0x0000000000000ULL),	/* + infinity */
+	DPCNST(1, DP_EMAX + 1, 0x0000000000000ULL),	/* - infinity */
+	DPCNST(0, DP_EMAX + 1, 0x7FFFFFFFFFFFFULL),	/* + ind legacy qNaN */
+	DPCNST(0, DP_EMAX + 1, 0x8000000000000ULL),	/* + indef 2008 qNaN */
+	DPCNST(0, DP_EMAX,     0xFFFFFFFFFFFFFULL),	/* + max */
+	DPCNST(1, DP_EMAX,     0xFFFFFFFFFFFFFULL),	/* - max */
+	DPCNST(0, DP_EMIN,     0x0000000000000ULL),	/* + min normal */
+	DPCNST(1, DP_EMIN,     0x0000000000000ULL),	/* - min normal */
+	DPCNST(0, DP_EMIN - 1, 0x0000000000001ULL),	/* + min denormal */
+	DPCNST(1, DP_EMIN - 1, 0x0000000000001ULL),	/* - min denormal */
+	DPCNST(0, 31,          0x0000000000000ULL),	/* + 1.0e31 */
+	DPCNST(0, 63,          0x0000000000000ULL),	/* + 1.0e63 */
+};
+
+#define SPCNST(s, b, m)							\
+	xPCNST(s, b, m, SP_EBIAS)
+
+const union ieee754sp __ieee754sp_spcvals[] = {
+	SPCNST(0, SP_EMIN - 1, 0x000000),	/* + zero   */
+	SPCNST(1, SP_EMIN - 1, 0x000000),	/* - zero   */
+	SPCNST(0, 0,	       0x000000),	/* + 1.0   */
+	SPCNST(1, 0,	       0x000000),	/* - 1.0   */
+	SPCNST(0, 3,	       0x200000),	/* + 10.0   */
+	SPCNST(1, 3,	       0x200000),	/* - 10.0   */
+	SPCNST(0, SP_EMAX + 1, 0x000000),	/* + infinity */
+	SPCNST(1, SP_EMAX + 1, 0x000000),	/* - infinity */
+	SPCNST(0, SP_EMAX + 1, 0x3FFFFF),	/* + indef legacy quiet NaN */
+	SPCNST(0, SP_EMAX + 1, 0x400000),	/* + indef 2008 quiet NaN */
+	SPCNST(0, SP_EMAX,     0x7FFFFF),	/* + max normal */
+	SPCNST(1, SP_EMAX,     0x7FFFFF),	/* - max normal */
+	SPCNST(0, SP_EMIN,     0x000000),	/* + min normal */
+	SPCNST(1, SP_EMIN,     0x000000),	/* - min normal */
+	SPCNST(0, SP_EMIN - 1, 0x000001),	/* + min denormal */
+	SPCNST(1, SP_EMIN - 1, 0x000001),	/* - min denormal */
+	SPCNST(0, 31,	       0x000000),	/* + 1.0e31 */
+	SPCNST(0, 63,	       0x000000),	/* + 1.0e63 */
+};

linux-7.2.3/arch/mips/math-emu/ieee754.h

@@ -0,0 +1,311 @@
+/* SPDX-License-Identifier: GPL-2.0-only */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ *
+ *  Nov 7, 2000
+ *  Modification to allow integration with Linux kernel
+ *
+ *  Kevin D. Kissell, kevink@mips.com and Carsten Langgard, carstenl@mips.com
+ *  Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
+ */
+#ifndef __ARCH_MIPS_MATH_EMU_IEEE754_H
+#define __ARCH_MIPS_MATH_EMU_IEEE754_H
+
+#include <linux/compiler.h>
+#include <asm/byteorder.h>
+#include <linux/kernel.h>
+#include <linux/types.h>
+#include <linux/sched.h>
+#include <asm/bitfield.h>
+
+union ieee754dp {
+	struct {
+		__BITFIELD_FIELD(unsigned int sign:1,
+		__BITFIELD_FIELD(unsigned int bexp:11,
+		__BITFIELD_FIELD(u64 mant:52,
+		;)))
+	};
+	u64 bits;
+};
+
+union ieee754sp {
+	struct {
+		__BITFIELD_FIELD(unsigned sign:1,
+		__BITFIELD_FIELD(unsigned bexp:8,
+		__BITFIELD_FIELD(unsigned mant:23,
+		;)))
+	};
+	u32 bits;
+};
+
+/*
+ * single precision (often aka float)
+*/
+int ieee754sp_class(union ieee754sp x);
+
+union ieee754sp ieee754sp_abs(union ieee754sp x);
+union ieee754sp ieee754sp_neg(union ieee754sp x);
+
+union ieee754sp ieee754sp_add(union ieee754sp x, union ieee754sp y);
+union ieee754sp ieee754sp_sub(union ieee754sp x, union ieee754sp y);
+union ieee754sp ieee754sp_mul(union ieee754sp x, union ieee754sp y);
+union ieee754sp ieee754sp_div(union ieee754sp x, union ieee754sp y);
+
+union ieee754sp ieee754sp_fint(int x);
+union ieee754sp ieee754sp_flong(s64 x);
+union ieee754sp ieee754sp_fdp(union ieee754dp x);
+union ieee754sp ieee754sp_rint(union ieee754sp x);
+
+int ieee754sp_tint(union ieee754sp x);
+s64 ieee754sp_tlong(union ieee754sp x);
+
+int ieee754sp_cmp(union ieee754sp x, union ieee754sp y, int cop, int sig);
+
+union ieee754sp ieee754sp_sqrt(union ieee754sp x);
+
+union ieee754sp ieee754sp_maddf(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y);
+union ieee754sp ieee754sp_msubf(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y);
+union ieee754sp ieee754sp_madd(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y);
+union ieee754sp ieee754sp_msub(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y);
+union ieee754sp ieee754sp_nmadd(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y);
+union ieee754sp ieee754sp_nmsub(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y);
+int ieee754sp_2008class(union ieee754sp x);
+union ieee754sp ieee754sp_fmin(union ieee754sp x, union ieee754sp y);
+union ieee754sp ieee754sp_fmina(union ieee754sp x, union ieee754sp y);
+union ieee754sp ieee754sp_fmax(union ieee754sp x, union ieee754sp y);
+union ieee754sp ieee754sp_fmaxa(union ieee754sp x, union ieee754sp y);
+
+/*
+ * double precision (often aka double)
+*/
+int ieee754dp_class(union ieee754dp x);
+
+union ieee754dp ieee754dp_add(union ieee754dp x, union ieee754dp y);
+union ieee754dp ieee754dp_sub(union ieee754dp x, union ieee754dp y);
+union ieee754dp ieee754dp_mul(union ieee754dp x, union ieee754dp y);
+union ieee754dp ieee754dp_div(union ieee754dp x, union ieee754dp y);
+
+union ieee754dp ieee754dp_abs(union ieee754dp x);
+union ieee754dp ieee754dp_neg(union ieee754dp x);
+
+union ieee754dp ieee754dp_fint(int x);
+union ieee754dp ieee754dp_flong(s64 x);
+union ieee754dp ieee754dp_fsp(union ieee754sp x);
+union ieee754dp ieee754dp_rint(union ieee754dp x);
+
+int ieee754dp_tint(union ieee754dp x);
+s64 ieee754dp_tlong(union ieee754dp x);
+
+int ieee754dp_cmp(union ieee754dp x, union ieee754dp y, int cop, int sig);
+
+union ieee754dp ieee754dp_sqrt(union ieee754dp x);
+
+union ieee754dp ieee754dp_maddf(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y);
+union ieee754dp ieee754dp_msubf(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y);
+union ieee754dp ieee754dp_madd(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y);
+union ieee754dp ieee754dp_msub(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y);
+union ieee754dp ieee754dp_nmadd(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y);
+union ieee754dp ieee754dp_nmsub(union ieee754dp z, union ieee754dp x,
+				union ieee754dp y);
+int ieee754dp_2008class(union ieee754dp x);
+union ieee754dp ieee754dp_fmin(union ieee754dp x, union ieee754dp y);
+union ieee754dp ieee754dp_fmina(union ieee754dp x, union ieee754dp y);
+union ieee754dp ieee754dp_fmax(union ieee754dp x, union ieee754dp y);
+union ieee754dp ieee754dp_fmaxa(union ieee754dp x, union ieee754dp y);
+
+
+/* 5 types of floating point number
+*/
+enum {
+	IEEE754_CLASS_NORM	= 0x00,
+	IEEE754_CLASS_ZERO	= 0x01,
+	IEEE754_CLASS_DNORM	= 0x02,
+	IEEE754_CLASS_INF	= 0x03,
+	IEEE754_CLASS_SNAN	= 0x04,
+	IEEE754_CLASS_QNAN	= 0x05,
+};
+
+/* exception numbers */
+#define IEEE754_INEXACT			0x01
+#define IEEE754_UNDERFLOW		0x02
+#define IEEE754_OVERFLOW		0x04
+#define IEEE754_ZERO_DIVIDE		0x08
+#define IEEE754_INVALID_OPERATION	0x10
+
+/* cmp operators
+*/
+#define IEEE754_CLT	0x01
+#define IEEE754_CEQ	0x02
+#define IEEE754_CGT	0x04
+#define IEEE754_CUN	0x08
+
+/*
+ * The control status register
+ */
+struct _ieee754_csr {
+	__BITFIELD_FIELD(unsigned fcc:7,	/* condition[7:1] */
+	__BITFIELD_FIELD(unsigned nod:1,	/* set 1 for no denormals */
+	__BITFIELD_FIELD(unsigned c:1,		/* condition[0] */
+	__BITFIELD_FIELD(unsigned pad0:3,
+	__BITFIELD_FIELD(unsigned abs2008:1,	/* IEEE 754-2008 ABS/NEG.fmt */
+	__BITFIELD_FIELD(unsigned nan2008:1,	/* IEEE 754-2008 NaN mode */
+	__BITFIELD_FIELD(unsigned cx:6,		/* exceptions this operation */
+	__BITFIELD_FIELD(unsigned mx:5,		/* exception enable  mask */
+	__BITFIELD_FIELD(unsigned sx:5,		/* exceptions total */
+	__BITFIELD_FIELD(unsigned rm:2,		/* current rounding mode */
+	;))))))))))
+};
+#define ieee754_csr (*(struct _ieee754_csr *)(&current->thread.fpu.fcr31))
+
+static inline unsigned int ieee754_getrm(void)
+{
+	return (ieee754_csr.rm);
+}
+
+static inline unsigned int ieee754_setrm(unsigned int rm)
+{
+	return (ieee754_csr.rm = rm);
+}
+
+/*
+ * get current exceptions
+ */
+static inline unsigned int ieee754_getcx(void)
+{
+	return (ieee754_csr.cx);
+}
+
+/* test for current exception condition
+ */
+static inline int ieee754_cxtest(unsigned int n)
+{
+	return (ieee754_csr.cx & n);
+}
+
+/*
+ * get sticky exceptions
+ */
+static inline unsigned int ieee754_getsx(void)
+{
+	return (ieee754_csr.sx);
+}
+
+/* clear sticky conditions
+*/
+static inline unsigned int ieee754_clrsx(void)
+{
+	return (ieee754_csr.sx = 0);
+}
+
+/* test for sticky exception condition
+ */
+static inline int ieee754_sxtest(unsigned int n)
+{
+	return (ieee754_csr.sx & n);
+}
+
+/* debugging */
+union ieee754sp ieee754sp_dump(char *s, union ieee754sp x);
+union ieee754dp ieee754dp_dump(char *s, union ieee754dp x);
+
+#define IEEE754_SPCVAL_PZERO		0	/* +0.0 */
+#define IEEE754_SPCVAL_NZERO		1	/* -0.0 */
+#define IEEE754_SPCVAL_PONE		2	/* +1.0 */
+#define IEEE754_SPCVAL_NONE		3	/* -1.0 */
+#define IEEE754_SPCVAL_PTEN		4	/* +10.0 */
+#define IEEE754_SPCVAL_NTEN		5	/* -10.0 */
+#define IEEE754_SPCVAL_PINFINITY	6	/* +inf */
+#define IEEE754_SPCVAL_NINFINITY	7	/* -inf */
+#define IEEE754_SPCVAL_INDEF_LEG	8	/* legacy quiet NaN */
+#define IEEE754_SPCVAL_INDEF_2008	9	/* IEEE 754-2008 quiet NaN */
+#define IEEE754_SPCVAL_PMAX		10	/* +max norm */
+#define IEEE754_SPCVAL_NMAX		11	/* -max norm */
+#define IEEE754_SPCVAL_PMIN		12	/* +min norm */
+#define IEEE754_SPCVAL_NMIN		13	/* -min norm */
+#define IEEE754_SPCVAL_PMIND		14	/* +min denorm */
+#define IEEE754_SPCVAL_NMIND		15	/* -min denorm */
+#define IEEE754_SPCVAL_P1E31		16	/* + 1.0e31 */
+#define IEEE754_SPCVAL_P1E63		17	/* + 1.0e63 */
+
+extern const union ieee754dp __ieee754dp_spcvals[];
+extern const union ieee754sp __ieee754sp_spcvals[];
+#define ieee754dp_spcvals ((const union ieee754dp *)__ieee754dp_spcvals)
+#define ieee754sp_spcvals ((const union ieee754sp *)__ieee754sp_spcvals)
+
+/*
+ * Return infinity with given sign
+ */
+#define ieee754dp_inf(sn)     (ieee754dp_spcvals[IEEE754_SPCVAL_PINFINITY+(sn)])
+#define ieee754dp_zero(sn)	(ieee754dp_spcvals[IEEE754_SPCVAL_PZERO+(sn)])
+#define ieee754dp_one(sn)	(ieee754dp_spcvals[IEEE754_SPCVAL_PONE+(sn)])
+#define ieee754dp_ten(sn)	(ieee754dp_spcvals[IEEE754_SPCVAL_PTEN+(sn)])
+#define ieee754dp_indef()	(ieee754dp_spcvals[IEEE754_SPCVAL_INDEF_LEG + \
+						   ieee754_csr.nan2008])
+#define ieee754dp_max(sn)	(ieee754dp_spcvals[IEEE754_SPCVAL_PMAX+(sn)])
+#define ieee754dp_min(sn)	(ieee754dp_spcvals[IEEE754_SPCVAL_PMIN+(sn)])
+#define ieee754dp_mind(sn)	(ieee754dp_spcvals[IEEE754_SPCVAL_PMIND+(sn)])
+#define ieee754dp_1e31()	(ieee754dp_spcvals[IEEE754_SPCVAL_P1E31])
+#define ieee754dp_1e63()	(ieee754dp_spcvals[IEEE754_SPCVAL_P1E63])
+
+#define ieee754sp_inf(sn)     (ieee754sp_spcvals[IEEE754_SPCVAL_PINFINITY+(sn)])
+#define ieee754sp_zero(sn)	(ieee754sp_spcvals[IEEE754_SPCVAL_PZERO+(sn)])
+#define ieee754sp_one(sn)	(ieee754sp_spcvals[IEEE754_SPCVAL_PONE+(sn)])
+#define ieee754sp_ten(sn)	(ieee754sp_spcvals[IEEE754_SPCVAL_PTEN+(sn)])
+#define ieee754sp_indef()	(ieee754sp_spcvals[IEEE754_SPCVAL_INDEF_LEG + \
+						   ieee754_csr.nan2008])
+#define ieee754sp_max(sn)	(ieee754sp_spcvals[IEEE754_SPCVAL_PMAX+(sn)])
+#define ieee754sp_min(sn)	(ieee754sp_spcvals[IEEE754_SPCVAL_PMIN+(sn)])
+#define ieee754sp_mind(sn)	(ieee754sp_spcvals[IEEE754_SPCVAL_PMIND+(sn)])
+#define ieee754sp_1e31()	(ieee754sp_spcvals[IEEE754_SPCVAL_P1E31])
+#define ieee754sp_1e63()	(ieee754sp_spcvals[IEEE754_SPCVAL_P1E63])
+
+/*
+ * Indefinite integer value
+ */
+static inline int ieee754si_indef(void)
+{
+	return ieee754_csr.nan2008 ? 0 : INT_MAX;
+}
+
+static inline s64 ieee754di_indef(void)
+{
+	return ieee754_csr.nan2008 ? 0 : S64_MAX;
+}
+
+/*
+ * Overflow integer value
+ */
+static inline int ieee754si_overflow(int xs)
+{
+	return ieee754_csr.nan2008 && xs ? INT_MIN : INT_MAX;
+}
+
+static inline s64 ieee754di_overflow(int xs)
+{
+	return ieee754_csr.nan2008 && xs ? S64_MIN : S64_MAX;
+}
+
+/* result types for xctx.rt */
+#define IEEE754_RT_SP	0
+#define IEEE754_RT_DP	1
+#define IEEE754_RT_XP	2
+#define IEEE754_RT_SI	3
+#define IEEE754_RT_DI	4
+
+/* compat */
+#define ieee754dp_fix(x)	ieee754dp_tint(x)
+#define ieee754sp_fix(x)	ieee754sp_tint(x)
+
+#endif /* __ARCH_MIPS_MATH_EMU_IEEE754_H */

linux-7.2.3/arch/mips/math-emu/ieee754d.c

@@ -0,0 +1,99 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * Some debug functions
+ *
+ * MIPS floating point support
+ *
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ *
+ *  Nov 7, 2000
+ *  Modified to build and operate in Linux kernel environment.
+ *
+ *  Kevin D. Kissell, kevink@mips.com and Carsten Langgaard, carstenl@mips.com
+ *  Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
+ */
+
+#include <linux/types.h>
+#include <linux/printk.h>
+#include "ieee754.h"
+#include "ieee754sp.h"
+#include "ieee754dp.h"
+
+union ieee754dp ieee754dp_dump(char *m, union ieee754dp x)
+{
+	int i;
+
+	printk("%s", m);
+	printk("<%08x,%08x>\n", (unsigned) (x.bits >> 32),
+	       (unsigned) x.bits);
+	printk("\t=");
+	switch (ieee754dp_class(x)) {
+	case IEEE754_CLASS_QNAN:
+	case IEEE754_CLASS_SNAN:
+		printk("Nan %c", DPSIGN(x) ? '-' : '+');
+		for (i = DP_FBITS - 1; i >= 0; i--)
+			printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
+		break;
+	case IEEE754_CLASS_INF:
+		printk("%cInfinity", DPSIGN(x) ? '-' : '+');
+		break;
+	case IEEE754_CLASS_ZERO:
+		printk("%cZero", DPSIGN(x) ? '-' : '+');
+		break;
+	case IEEE754_CLASS_DNORM:
+		printk("%c0.", DPSIGN(x) ? '-' : '+');
+		for (i = DP_FBITS - 1; i >= 0; i--)
+			printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
+		printk("e%d", DPBEXP(x) - DP_EBIAS);
+		break;
+	case IEEE754_CLASS_NORM:
+		printk("%c1.", DPSIGN(x) ? '-' : '+');
+		for (i = DP_FBITS - 1; i >= 0; i--)
+			printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
+		printk("e%d", DPBEXP(x) - DP_EBIAS);
+		break;
+	default:
+		printk("Illegal/Unknown IEEE754 value class");
+	}
+	printk("\n");
+	return x;
+}
+
+union ieee754sp ieee754sp_dump(char *m, union ieee754sp x)
+{
+	int i;
+
+	printk("%s=", m);
+	printk("<%08x>\n", (unsigned) x.bits);
+	printk("\t=");
+	switch (ieee754sp_class(x)) {
+	case IEEE754_CLASS_QNAN:
+	case IEEE754_CLASS_SNAN:
+		printk("Nan %c", SPSIGN(x) ? '-' : '+');
+		for (i = SP_FBITS - 1; i >= 0; i--)
+			printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
+		break;
+	case IEEE754_CLASS_INF:
+		printk("%cInfinity", SPSIGN(x) ? '-' : '+');
+		break;
+	case IEEE754_CLASS_ZERO:
+		printk("%cZero", SPSIGN(x) ? '-' : '+');
+		break;
+	case IEEE754_CLASS_DNORM:
+		printk("%c0.", SPSIGN(x) ? '-' : '+');
+		for (i = SP_FBITS - 1; i >= 0; i--)
+			printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
+		printk("e%d", SPBEXP(x) - SP_EBIAS);
+		break;
+	case IEEE754_CLASS_NORM:
+		printk("%c1.", SPSIGN(x) ? '-' : '+');
+		for (i = SP_FBITS - 1; i >= 0; i--)
+			printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
+		printk("e%d", SPBEXP(x) - SP_EBIAS);
+		break;
+	default:
+		printk("Illegal/Unknown IEEE754 value class");
+	}
+	printk("\n");
+	return x;
+}

linux-7.2.3/arch/mips/math-emu/ieee754dp.c

@@ -0,0 +1,197 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * double precision: common utilities
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include <linux/compiler.h>
+
+#include "ieee754dp.h"
+
+int ieee754dp_class(union ieee754dp x)
+{
+	COMPXDP;
+	EXPLODEXDP;
+	return xc;
+}
+
+static inline int ieee754dp_isnan(union ieee754dp x)
+{
+	return ieee754_class_nan(ieee754dp_class(x));
+}
+
+static inline int ieee754dp_issnan(union ieee754dp x)
+{
+	int qbit;
+
+	assert(ieee754dp_isnan(x));
+	qbit = (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_FBITS - 1);
+	return ieee754_csr.nan2008 ^ qbit;
+}
+
+
+/*
+ * Raise the Invalid Operation IEEE 754 exception
+ * and convert the signaling NaN supplied to a quiet NaN.
+ */
+union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r)
+{
+	assert(ieee754dp_issnan(r));
+
+	ieee754_setcx(IEEE754_INVALID_OPERATION);
+	if (ieee754_csr.nan2008) {
+		DPMANT(r) |= DP_MBIT(DP_FBITS - 1);
+	} else {
+		DPMANT(r) &= ~DP_MBIT(DP_FBITS - 1);
+		if (!ieee754dp_isnan(r))
+			DPMANT(r) |= DP_MBIT(DP_FBITS - 2);
+	}
+
+	return r;
+}
+
+static u64 ieee754dp_get_rounding(int sn, u64 xm)
+{
+	/* inexact must round of 3 bits
+	 */
+	if (xm & (DP_MBIT(3) - 1)) {
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RZ:
+			break;
+		case FPU_CSR_RN:
+			xm += 0x3 + ((xm >> 3) & 1);
+			/* xm += (xm&0x8)?0x4:0x3 */
+			break;
+		case FPU_CSR_RU:	/* toward +Infinity */
+			if (!sn)	/* ?? */
+				xm += 0x8;
+			break;
+		case FPU_CSR_RD:	/* toward -Infinity */
+			if (sn) /* ?? */
+				xm += 0x8;
+			break;
+		}
+	}
+	return xm;
+}
+
+
+/* generate a normal/denormal number with over,under handling
+ * sn is sign
+ * xe is an unbiased exponent
+ * xm is 3bit extended precision value.
+ */
+union ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
+{
+	assert(xm);		/* we don't gen exact zeros (probably should) */
+
+	assert((xm >> (DP_FBITS + 1 + 3)) == 0);	/* no excess */
+	assert(xm & (DP_HIDDEN_BIT << 3));
+
+	if (xe < DP_EMIN) {
+		/* strip lower bits */
+		int es = DP_EMIN - xe;
+
+		if (ieee754_csr.nod) {
+			ieee754_setcx(IEEE754_UNDERFLOW);
+			ieee754_setcx(IEEE754_INEXACT);
+
+			switch(ieee754_csr.rm) {
+			case FPU_CSR_RN:
+			case FPU_CSR_RZ:
+				return ieee754dp_zero(sn);
+			case FPU_CSR_RU:    /* toward +Infinity */
+				if (sn == 0)
+					return ieee754dp_min(0);
+				else
+					return ieee754dp_zero(1);
+			case FPU_CSR_RD:    /* toward -Infinity */
+				if (sn == 0)
+					return ieee754dp_zero(0);
+				else
+					return ieee754dp_min(1);
+			}
+		}
+
+		if (xe == DP_EMIN - 1 &&
+		    ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3))
+		{
+			/* Not tiny after rounding */
+			ieee754_setcx(IEEE754_INEXACT);
+			xm = ieee754dp_get_rounding(sn, xm);
+			xm >>= 1;
+			/* Clear grs bits */
+			xm &= ~(DP_MBIT(3) - 1);
+			xe++;
+		}
+		else {
+			/* sticky right shift es bits
+			 */
+			xm = XDPSRS(xm, es);
+			xe += es;
+			assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
+			assert(xe == DP_EMIN);
+		}
+	}
+	if (xm & (DP_MBIT(3) - 1)) {
+		ieee754_setcx(IEEE754_INEXACT);
+		if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
+			ieee754_setcx(IEEE754_UNDERFLOW);
+		}
+
+		/* inexact must round of 3 bits
+		 */
+		xm = ieee754dp_get_rounding(sn, xm);
+		/* adjust exponent for rounding add overflowing
+		 */
+		if (xm >> (DP_FBITS + 3 + 1)) {
+			/* add causes mantissa overflow */
+			xm >>= 1;
+			xe++;
+		}
+	}
+	/* strip grs bits */
+	xm >>= 3;
+
+	assert((xm >> (DP_FBITS + 1)) == 0);	/* no excess */
+	assert(xe >= DP_EMIN);
+
+	if (xe > DP_EMAX) {
+		ieee754_setcx(IEEE754_OVERFLOW);
+		ieee754_setcx(IEEE754_INEXACT);
+		/* -O can be table indexed by (rm,sn) */
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RN:
+			return ieee754dp_inf(sn);
+		case FPU_CSR_RZ:
+			return ieee754dp_max(sn);
+		case FPU_CSR_RU:	/* toward +Infinity */
+			if (sn == 0)
+				return ieee754dp_inf(0);
+			else
+				return ieee754dp_max(1);
+		case FPU_CSR_RD:	/* toward -Infinity */
+			if (sn == 0)
+				return ieee754dp_max(0);
+			else
+				return ieee754dp_inf(1);
+		}
+	}
+	/* gen norm/denorm/zero */
+
+	if ((xm & DP_HIDDEN_BIT) == 0) {
+		/* we underflow (tiny/zero) */
+		assert(xe == DP_EMIN);
+		if (ieee754_csr.mx & IEEE754_UNDERFLOW)
+			ieee754_setcx(IEEE754_UNDERFLOW);
+		return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
+	} else {
+		assert((xm >> (DP_FBITS + 1)) == 0);	/* no excess */
+		assert(xm & DP_HIDDEN_BIT);
+
+		return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
+	}
+}

linux-7.2.3/arch/mips/math-emu/ieee754dp.h

@@ -0,0 +1,73 @@
+/* SPDX-License-Identifier: GPL-2.0-only */
+/*
+ * IEEE754 floating point
+ * double precision internal header file
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include <linux/compiler.h>
+
+#include "ieee754int.h"
+
+#define assert(expr) ((void)0)
+
+#define DP_EBIAS	1023
+#define DP_EMIN		(-1022)
+#define DP_EMAX		1023
+#define DP_FBITS	52
+#define DP_MBITS	52
+
+#define DP_MBIT(x)	((u64)1 << (x))
+#define DP_HIDDEN_BIT	DP_MBIT(DP_FBITS)
+#define DP_SIGN_BIT	DP_MBIT(63)
+
+#define DPSIGN(dp)	(dp.sign)
+#define DPBEXP(dp)	(dp.bexp)
+#define DPMANT(dp)	(dp.mant)
+
+static inline int ieee754dp_finite(union ieee754dp x)
+{
+	return DPBEXP(x) != DP_EMAX + 1 + DP_EBIAS;
+}
+
+/* 3bit extended double precision sticky right shift */
+#define XDPSRS(v,rs)	\
+	((rs > (DP_FBITS+3))?1:((v) >> (rs)) | ((v) << (64-(rs)) != 0))
+
+#define XDPSRSX1() \
+	(xe++, (xm = (xm >> 1) | (xm & 1)))
+
+#define XDPSRS1(v)	\
+	(((v) >> 1) | ((v) & 1))
+
+/* 32bit * 32bit => 64bit unsigned integer multiplication */
+#define DPXMULT(x, y)	((u64)(x) * (u64)y)
+
+/* convert denormal to normalized with extended exponent */
+#define DPDNORMx(m,e) \
+	while ((m >> DP_FBITS) == 0) { m <<= 1; e--; }
+#define DPDNORMX	DPDNORMx(xm, xe)
+#define DPDNORMY	DPDNORMx(ym, ye)
+#define DPDNORMZ	DPDNORMx(zm, ze)
+
+static inline union ieee754dp builddp(int s, int bx, u64 m)
+{
+	union ieee754dp r;
+
+	assert((s) == 0 || (s) == 1);
+	assert((bx) >= DP_EMIN - 1 + DP_EBIAS
+	       && (bx) <= DP_EMAX + 1 + DP_EBIAS);
+	assert(((m) >> DP_FBITS) == 0);
+
+	r.sign = s;
+	r.bexp = bx;
+	r.mant = m;
+
+	return r;
+}
+
+extern union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp);
+extern union ieee754dp ieee754dp_format(int, int, u64);

linux-7.2.3/arch/mips/math-emu/ieee754int.h

@@ -0,0 +1,149 @@
+/* SPDX-License-Identifier: GPL-2.0-only */
+/*
+ * IEEE754 floating point
+ * common internal header file
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+#ifndef __IEEE754INT_H
+#define __IEEE754INT_H
+
+#include "ieee754.h"
+
+#define CLPAIR(x, y)	((x)*6+(y))
+
+enum maddf_flags {
+	MADDF_NEGATE_PRODUCT	= 1 << 0,
+	MADDF_NEGATE_ADDITION	= 1 << 1,
+};
+
+static inline void ieee754_clearcx(void)
+{
+	ieee754_csr.cx = 0;
+}
+
+static inline void ieee754_setcx(const unsigned int flags)
+{
+	ieee754_csr.cx |= flags;
+	ieee754_csr.sx |= flags;
+}
+
+static inline int ieee754_setandtestcx(const unsigned int x)
+{
+	ieee754_setcx(x);
+
+	return ieee754_csr.mx & x;
+}
+
+static inline int ieee754_class_nan(int xc)
+{
+	return xc >= IEEE754_CLASS_SNAN;
+}
+
+#define COMPXSP \
+	unsigned int xm; int xe; int xs __maybe_unused; int xc
+
+#define COMPYSP \
+	unsigned int ym; int ye; int ys; int yc
+
+#define COMPZSP \
+	unsigned int zm; int ze; int zs; int zc
+
+#define EXPLODESP(v, vc, vs, ve, vm)					\
+{									\
+	vs = SPSIGN(v);							\
+	ve = SPBEXP(v);							\
+	vm = SPMANT(v);							\
+	if (ve == SP_EMAX+1+SP_EBIAS) {					\
+		if (vm == 0)						\
+			vc = IEEE754_CLASS_INF;				\
+		else if (ieee754_csr.nan2008 ^ !(vm & SP_MBIT(SP_FBITS - 1))) \
+			vc = IEEE754_CLASS_QNAN;			\
+		else							\
+			vc = IEEE754_CLASS_SNAN;			\
+	} else if (ve == SP_EMIN-1+SP_EBIAS) {				\
+		if (vm) {						\
+			ve = SP_EMIN;					\
+			vc = IEEE754_CLASS_DNORM;			\
+		} else							\
+			vc = IEEE754_CLASS_ZERO;			\
+	} else {							\
+		ve -= SP_EBIAS;						\
+		vm |= SP_HIDDEN_BIT;					\
+		vc = IEEE754_CLASS_NORM;				\
+	}								\
+}
+#define EXPLODEXSP EXPLODESP(x, xc, xs, xe, xm)
+#define EXPLODEYSP EXPLODESP(y, yc, ys, ye, ym)
+#define EXPLODEZSP EXPLODESP(z, zc, zs, ze, zm)
+
+
+#define COMPXDP \
+	u64 xm; int xe; int xs __maybe_unused; int xc
+
+#define COMPYDP \
+	u64 ym; int ye; int ys; int yc
+
+#define COMPZDP \
+	u64 zm; int ze; int zs; int zc
+
+#define EXPLODEDP(v, vc, vs, ve, vm)					\
+{									\
+	vm = DPMANT(v);							\
+	vs = DPSIGN(v);							\
+	ve = DPBEXP(v);							\
+	if (ve == DP_EMAX+1+DP_EBIAS) {					\
+		if (vm == 0)						\
+			vc = IEEE754_CLASS_INF;				\
+		else if (ieee754_csr.nan2008 ^ !(vm & DP_MBIT(DP_FBITS - 1))) \
+			vc = IEEE754_CLASS_QNAN;			\
+		else							\
+			vc = IEEE754_CLASS_SNAN;			\
+	} else if (ve == DP_EMIN-1+DP_EBIAS) {				\
+		if (vm) {						\
+			ve = DP_EMIN;					\
+			vc = IEEE754_CLASS_DNORM;			\
+		} else							\
+			vc = IEEE754_CLASS_ZERO;			\
+	} else {							\
+		ve -= DP_EBIAS;						\
+		vm |= DP_HIDDEN_BIT;					\
+		vc = IEEE754_CLASS_NORM;				\
+	}								\
+}
+#define EXPLODEXDP EXPLODEDP(x, xc, xs, xe, xm)
+#define EXPLODEYDP EXPLODEDP(y, yc, ys, ye, ym)
+#define EXPLODEZDP EXPLODEDP(z, zc, zs, ze, zm)
+
+#define FLUSHDP(v, vc, vs, ve, vm)					\
+	if (vc==IEEE754_CLASS_DNORM) {					\
+		if (ieee754_csr.nod) {					\
+			ieee754_setcx(IEEE754_INEXACT);			\
+			vc = IEEE754_CLASS_ZERO;			\
+			ve = DP_EMIN-1+DP_EBIAS;			\
+			vm = 0;						\
+			v = ieee754dp_zero(vs);				\
+		}							\
+	}
+
+#define FLUSHSP(v, vc, vs, ve, vm)					\
+	if (vc==IEEE754_CLASS_DNORM) {					\
+		if (ieee754_csr.nod) {					\
+			ieee754_setcx(IEEE754_INEXACT);			\
+			vc = IEEE754_CLASS_ZERO;			\
+			ve = SP_EMIN-1+SP_EBIAS;			\
+			vm = 0;						\
+			v = ieee754sp_zero(vs);				\
+		}							\
+	}
+
+#define FLUSHXDP FLUSHDP(x, xc, xs, xe, xm)
+#define FLUSHYDP FLUSHDP(y, yc, ys, ye, ym)
+#define FLUSHZDP FLUSHDP(z, zc, zs, ze, zm)
+#define FLUSHXSP FLUSHSP(x, xc, xs, xe, xm)
+#define FLUSHYSP FLUSHSP(y, yc, ys, ye, ym)
+#define FLUSHZSP FLUSHSP(z, zc, zs, ze, zm)
+
+#endif /* __IEEE754INT_H  */

linux-7.2.3/arch/mips/math-emu/ieee754sp.c

@@ -0,0 +1,196 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include <linux/compiler.h>
+
+#include "ieee754sp.h"
+
+int ieee754sp_class(union ieee754sp x)
+{
+	COMPXSP;
+	EXPLODEXSP;
+	return xc;
+}
+
+static inline int ieee754sp_isnan(union ieee754sp x)
+{
+	return ieee754_class_nan(ieee754sp_class(x));
+}
+
+static inline int ieee754sp_issnan(union ieee754sp x)
+{
+	int qbit;
+
+	assert(ieee754sp_isnan(x));
+	qbit = (SPMANT(x) & SP_MBIT(SP_FBITS - 1)) == SP_MBIT(SP_FBITS - 1);
+	return ieee754_csr.nan2008 ^ qbit;
+}
+
+
+/*
+ * Raise the Invalid Operation IEEE 754 exception
+ * and convert the signaling NaN supplied to a quiet NaN.
+ */
+union ieee754sp __cold ieee754sp_nanxcpt(union ieee754sp r)
+{
+	assert(ieee754sp_issnan(r));
+
+	ieee754_setcx(IEEE754_INVALID_OPERATION);
+	if (ieee754_csr.nan2008) {
+		SPMANT(r) |= SP_MBIT(SP_FBITS - 1);
+	} else {
+		SPMANT(r) &= ~SP_MBIT(SP_FBITS - 1);
+		if (!ieee754sp_isnan(r))
+			SPMANT(r) |= SP_MBIT(SP_FBITS - 2);
+	}
+
+	return r;
+}
+
+static unsigned int ieee754sp_get_rounding(int sn, unsigned int xm)
+{
+	/* inexact must round of 3 bits
+	 */
+	if (xm & (SP_MBIT(3) - 1)) {
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RZ:
+			break;
+		case FPU_CSR_RN:
+			xm += 0x3 + ((xm >> 3) & 1);
+			/* xm += (xm&0x8)?0x4:0x3 */
+			break;
+		case FPU_CSR_RU:	/* toward +Infinity */
+			if (!sn)	/* ?? */
+				xm += 0x8;
+			break;
+		case FPU_CSR_RD:	/* toward -Infinity */
+			if (sn) /* ?? */
+				xm += 0x8;
+			break;
+		}
+	}
+	return xm;
+}
+
+
+/* generate a normal/denormal number with over,under handling
+ * sn is sign
+ * xe is an unbiased exponent
+ * xm is 3bit extended precision value.
+ */
+union ieee754sp ieee754sp_format(int sn, int xe, unsigned int xm)
+{
+	assert(xm);		/* we don't gen exact zeros (probably should) */
+
+	assert((xm >> (SP_FBITS + 1 + 3)) == 0);	/* no excess */
+	assert(xm & (SP_HIDDEN_BIT << 3));
+
+	if (xe < SP_EMIN) {
+		/* strip lower bits */
+		int es = SP_EMIN - xe;
+
+		if (ieee754_csr.nod) {
+			ieee754_setcx(IEEE754_UNDERFLOW);
+			ieee754_setcx(IEEE754_INEXACT);
+
+			switch(ieee754_csr.rm) {
+			case FPU_CSR_RN:
+			case FPU_CSR_RZ:
+				return ieee754sp_zero(sn);
+			case FPU_CSR_RU:      /* toward +Infinity */
+				if (sn == 0)
+					return ieee754sp_min(0);
+				else
+					return ieee754sp_zero(1);
+			case FPU_CSR_RD:      /* toward -Infinity */
+				if (sn == 0)
+					return ieee754sp_zero(0);
+				else
+					return ieee754sp_min(1);
+			}
+		}
+
+		if (xe == SP_EMIN - 1 &&
+		    ieee754sp_get_rounding(sn, xm) >> (SP_FBITS + 1 + 3))
+		{
+			/* Not tiny after rounding */
+			ieee754_setcx(IEEE754_INEXACT);
+			xm = ieee754sp_get_rounding(sn, xm);
+			xm >>= 1;
+			/* Clear grs bits */
+			xm &= ~(SP_MBIT(3) - 1);
+			xe++;
+		} else {
+			/* sticky right shift es bits
+			 */
+			xm = XSPSRS(xm, es);
+			xe += es;
+			assert((xm & (SP_HIDDEN_BIT << 3)) == 0);
+			assert(xe == SP_EMIN);
+		}
+	}
+	if (xm & (SP_MBIT(3) - 1)) {
+		ieee754_setcx(IEEE754_INEXACT);
+		if ((xm & (SP_HIDDEN_BIT << 3)) == 0) {
+			ieee754_setcx(IEEE754_UNDERFLOW);
+		}
+
+		/* inexact must round of 3 bits
+		 */
+		xm = ieee754sp_get_rounding(sn, xm);
+		/* adjust exponent for rounding add overflowing
+		 */
+		if (xm >> (SP_FBITS + 1 + 3)) {
+			/* add causes mantissa overflow */
+			xm >>= 1;
+			xe++;
+		}
+	}
+	/* strip grs bits */
+	xm >>= 3;
+
+	assert((xm >> (SP_FBITS + 1)) == 0);	/* no excess */
+	assert(xe >= SP_EMIN);
+
+	if (xe > SP_EMAX) {
+		ieee754_setcx(IEEE754_OVERFLOW);
+		ieee754_setcx(IEEE754_INEXACT);
+		/* -O can be table indexed by (rm,sn) */
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RN:
+			return ieee754sp_inf(sn);
+		case FPU_CSR_RZ:
+			return ieee754sp_max(sn);
+		case FPU_CSR_RU:	/* toward +Infinity */
+			if (sn == 0)
+				return ieee754sp_inf(0);
+			else
+				return ieee754sp_max(1);
+		case FPU_CSR_RD:	/* toward -Infinity */
+			if (sn == 0)
+				return ieee754sp_max(0);
+			else
+				return ieee754sp_inf(1);
+		}
+	}
+	/* gen norm/denorm/zero */
+
+	if ((xm & SP_HIDDEN_BIT) == 0) {
+		/* we underflow (tiny/zero) */
+		assert(xe == SP_EMIN);
+		if (ieee754_csr.mx & IEEE754_UNDERFLOW)
+			ieee754_setcx(IEEE754_UNDERFLOW);
+		return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm);
+	} else {
+		assert((xm >> (SP_FBITS + 1)) == 0);	/* no excess */
+		assert(xm & SP_HIDDEN_BIT);
+
+		return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
+	}
+}

linux-7.2.3/arch/mips/math-emu/ieee754sp.h

@@ -0,0 +1,77 @@
+/* SPDX-License-Identifier: GPL-2.0-only */
+/*
+ * IEEE754 floating point
+ * double precision internal header file
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include <linux/compiler.h>
+
+#include "ieee754int.h"
+
+#define assert(expr) ((void)0)
+
+#define SP_EBIAS	127
+#define SP_EMIN		(-126)
+#define SP_EMAX		127
+#define SP_FBITS	23
+#define SP_MBITS	23
+
+#define SP_MBIT(x)	((u32)1 << (x))
+#define SP_HIDDEN_BIT	SP_MBIT(SP_FBITS)
+#define SP_SIGN_BIT	SP_MBIT(31)
+
+#define SPSIGN(sp)	(sp.sign)
+#define SPBEXP(sp)	(sp.bexp)
+#define SPMANT(sp)	(sp.mant)
+
+static inline int ieee754sp_finite(union ieee754sp x)
+{
+	return SPBEXP(x) != SP_EMAX + 1 + SP_EBIAS;
+}
+
+/* 64 bit right shift with rounding */
+#define XSPSRS64(v, rs)						\
+	(((rs) >= 64) ? ((v) != 0) : ((v) >> (rs)) | ((v) << (64-(rs)) != 0))
+
+/* 3bit extended single precision sticky right shift */
+#define XSPSRS(v, rs)						\
+	((rs > (SP_FBITS+3))?1:((v) >> (rs)) | ((v) << (32-(rs)) != 0))
+
+#define XSPSRS1(m) \
+	((m >> 1) | (m & 1))
+
+#define SPXSRSX1() \
+	(xe++, (xm = XSPSRS1(xm)))
+
+#define SPXSRSY1() \
+	(ye++, (ym = XSPSRS1(ym)))
+
+/* convert denormal to normalized with extended exponent */
+#define SPDNORMx(m,e) \
+	while ((m >> SP_FBITS) == 0) { m <<= 1; e--; }
+#define SPDNORMX	SPDNORMx(xm, xe)
+#define SPDNORMY	SPDNORMx(ym, ye)
+#define SPDNORMZ	SPDNORMx(zm, ze)
+
+static inline union ieee754sp buildsp(int s, int bx, unsigned int m)
+{
+	union ieee754sp r;
+
+	assert((s) == 0 || (s) == 1);
+	assert((bx) >= SP_EMIN - 1 + SP_EBIAS
+	       && (bx) <= SP_EMAX + 1 + SP_EBIAS);
+	assert(((m) >> SP_FBITS) == 0);
+
+	r.sign = s;
+	r.bexp = bx;
+	r.mant = m;
+
+	return r;
+}
+
+extern union ieee754sp __cold ieee754sp_nanxcpt(union ieee754sp);
+extern union ieee754sp ieee754sp_format(int, int, unsigned);

linux-7.2.3/arch/mips/math-emu/me-debugfs.c

@@ -0,0 +1,353 @@
+// SPDX-License-Identifier: GPL-2.0
+#include <linux/cpumask.h>
+#include <linux/debugfs.h>
+#include <linux/fs.h>
+#include <linux/init.h>
+#include <linux/percpu.h>
+#include <linux/types.h>
+#include <asm/debug.h>
+#include <asm/fpu_emulator.h>
+#include <asm/local.h>
+
+DEFINE_PER_CPU(struct mips_fpu_emulator_stats, fpuemustats);
+
+static int fpuemu_stat_get(void *data, u64 *val)
+{
+	int cpu;
+	unsigned long sum = 0;
+
+	for_each_online_cpu(cpu) {
+		struct mips_fpu_emulator_stats *ps;
+		local_t *pv;
+
+		ps = &per_cpu(fpuemustats, cpu);
+		pv = (void *)ps + (unsigned long)data;
+		sum += local_read(pv);
+	}
+	*val = sum;
+	return 0;
+}
+DEFINE_SIMPLE_ATTRIBUTE(fops_fpuemu_stat, fpuemu_stat_get, NULL, "%llu\n");
+
+/*
+ * Used to obtain names for a debugfs instruction counter, given field name
+ * in fpuemustats structure. For example, for input "cmp_sueq_d", the output
+ * would be "cmp.sueq.d". This is needed since dots are not allowed to be
+ * used in structure field names, and are, on the other hand, desired to be
+ * used in debugfs item names to be clearly associated to corresponding
+ * MIPS FPU instructions.
+ */
+static void adjust_instruction_counter_name(char *out_name, char *in_name, size_t len)
+{
+	int i = 0;
+
+	strscpy(out_name, in_name, len);
+	while (in_name[i] != '\0') {
+		if (out_name[i] == '_')
+			out_name[i] = '.';
+		i++;
+	}
+}
+
+static int fpuemustats_clear_show(struct seq_file *s, void *unused)
+{
+	__this_cpu_write((fpuemustats).emulated, 0);
+	__this_cpu_write((fpuemustats).loads, 0);
+	__this_cpu_write((fpuemustats).stores, 0);
+	__this_cpu_write((fpuemustats).branches, 0);
+	__this_cpu_write((fpuemustats).cp1ops, 0);
+	__this_cpu_write((fpuemustats).cp1xops, 0);
+	__this_cpu_write((fpuemustats).errors, 0);
+	__this_cpu_write((fpuemustats).ieee754_inexact, 0);
+	__this_cpu_write((fpuemustats).ieee754_underflow, 0);
+	__this_cpu_write((fpuemustats).ieee754_overflow, 0);
+	__this_cpu_write((fpuemustats).ieee754_zerodiv, 0);
+	__this_cpu_write((fpuemustats).ieee754_invalidop, 0);
+	__this_cpu_write((fpuemustats).ds_emul, 0);
+
+	__this_cpu_write((fpuemustats).abs_s, 0);
+	__this_cpu_write((fpuemustats).abs_d, 0);
+	__this_cpu_write((fpuemustats).add_s, 0);
+	__this_cpu_write((fpuemustats).add_d, 0);
+	__this_cpu_write((fpuemustats).bc1eqz, 0);
+	__this_cpu_write((fpuemustats).bc1nez, 0);
+	__this_cpu_write((fpuemustats).ceil_w_s, 0);
+	__this_cpu_write((fpuemustats).ceil_w_d, 0);
+	__this_cpu_write((fpuemustats).ceil_l_s, 0);
+	__this_cpu_write((fpuemustats).ceil_l_d, 0);
+	__this_cpu_write((fpuemustats).class_s, 0);
+	__this_cpu_write((fpuemustats).class_d, 0);
+	__this_cpu_write((fpuemustats).cmp_af_s, 0);
+	__this_cpu_write((fpuemustats).cmp_af_d, 0);
+	__this_cpu_write((fpuemustats).cmp_eq_s, 0);
+	__this_cpu_write((fpuemustats).cmp_eq_d, 0);
+	__this_cpu_write((fpuemustats).cmp_le_s, 0);
+	__this_cpu_write((fpuemustats).cmp_le_d, 0);
+	__this_cpu_write((fpuemustats).cmp_lt_s, 0);
+	__this_cpu_write((fpuemustats).cmp_lt_d, 0);
+	__this_cpu_write((fpuemustats).cmp_ne_s, 0);
+	__this_cpu_write((fpuemustats).cmp_ne_d, 0);
+	__this_cpu_write((fpuemustats).cmp_or_s, 0);
+	__this_cpu_write((fpuemustats).cmp_or_d, 0);
+	__this_cpu_write((fpuemustats).cmp_ueq_s, 0);
+	__this_cpu_write((fpuemustats).cmp_ueq_d, 0);
+	__this_cpu_write((fpuemustats).cmp_ule_s, 0);
+	__this_cpu_write((fpuemustats).cmp_ule_d, 0);
+	__this_cpu_write((fpuemustats).cmp_ult_s, 0);
+	__this_cpu_write((fpuemustats).cmp_ult_d, 0);
+	__this_cpu_write((fpuemustats).cmp_un_s, 0);
+	__this_cpu_write((fpuemustats).cmp_un_d, 0);
+	__this_cpu_write((fpuemustats).cmp_une_s, 0);
+	__this_cpu_write((fpuemustats).cmp_une_d, 0);
+	__this_cpu_write((fpuemustats).cmp_saf_s, 0);
+	__this_cpu_write((fpuemustats).cmp_saf_d, 0);
+	__this_cpu_write((fpuemustats).cmp_seq_s, 0);
+	__this_cpu_write((fpuemustats).cmp_seq_d, 0);
+	__this_cpu_write((fpuemustats).cmp_sle_s, 0);
+	__this_cpu_write((fpuemustats).cmp_sle_d, 0);
+	__this_cpu_write((fpuemustats).cmp_slt_s, 0);
+	__this_cpu_write((fpuemustats).cmp_slt_d, 0);
+	__this_cpu_write((fpuemustats).cmp_sne_s, 0);
+	__this_cpu_write((fpuemustats).cmp_sne_d, 0);
+	__this_cpu_write((fpuemustats).cmp_sor_s, 0);
+	__this_cpu_write((fpuemustats).cmp_sor_d, 0);
+	__this_cpu_write((fpuemustats).cmp_sueq_s, 0);
+	__this_cpu_write((fpuemustats).cmp_sueq_d, 0);
+	__this_cpu_write((fpuemustats).cmp_sule_s, 0);
+	__this_cpu_write((fpuemustats).cmp_sule_d, 0);
+	__this_cpu_write((fpuemustats).cmp_sult_s, 0);
+	__this_cpu_write((fpuemustats).cmp_sult_d, 0);
+	__this_cpu_write((fpuemustats).cmp_sun_s, 0);
+	__this_cpu_write((fpuemustats).cmp_sun_d, 0);
+	__this_cpu_write((fpuemustats).cmp_sune_s, 0);
+	__this_cpu_write((fpuemustats).cmp_sune_d, 0);
+	__this_cpu_write((fpuemustats).cvt_d_l, 0);
+	__this_cpu_write((fpuemustats).cvt_d_s, 0);
+	__this_cpu_write((fpuemustats).cvt_d_w, 0);
+	__this_cpu_write((fpuemustats).cvt_l_s, 0);
+	__this_cpu_write((fpuemustats).cvt_l_d, 0);
+	__this_cpu_write((fpuemustats).cvt_s_d, 0);
+	__this_cpu_write((fpuemustats).cvt_s_l, 0);
+	__this_cpu_write((fpuemustats).cvt_s_w, 0);
+	__this_cpu_write((fpuemustats).cvt_w_s, 0);
+	__this_cpu_write((fpuemustats).cvt_w_d, 0);
+	__this_cpu_write((fpuemustats).div_s, 0);
+	__this_cpu_write((fpuemustats).div_d, 0);
+	__this_cpu_write((fpuemustats).floor_w_s, 0);
+	__this_cpu_write((fpuemustats).floor_w_d, 0);
+	__this_cpu_write((fpuemustats).floor_l_s, 0);
+	__this_cpu_write((fpuemustats).floor_l_d, 0);
+	__this_cpu_write((fpuemustats).maddf_s, 0);
+	__this_cpu_write((fpuemustats).maddf_d, 0);
+	__this_cpu_write((fpuemustats).max_s, 0);
+	__this_cpu_write((fpuemustats).max_d, 0);
+	__this_cpu_write((fpuemustats).maxa_s, 0);
+	__this_cpu_write((fpuemustats).maxa_d, 0);
+	__this_cpu_write((fpuemustats).min_s, 0);
+	__this_cpu_write((fpuemustats).min_d, 0);
+	__this_cpu_write((fpuemustats).mina_s, 0);
+	__this_cpu_write((fpuemustats).mina_d, 0);
+	__this_cpu_write((fpuemustats).mov_s, 0);
+	__this_cpu_write((fpuemustats).mov_d, 0);
+	__this_cpu_write((fpuemustats).msubf_s, 0);
+	__this_cpu_write((fpuemustats).msubf_d, 0);
+	__this_cpu_write((fpuemustats).mul_s, 0);
+	__this_cpu_write((fpuemustats).mul_d, 0);
+	__this_cpu_write((fpuemustats).neg_s, 0);
+	__this_cpu_write((fpuemustats).neg_d, 0);
+	__this_cpu_write((fpuemustats).recip_s, 0);
+	__this_cpu_write((fpuemustats).recip_d, 0);
+	__this_cpu_write((fpuemustats).rint_s, 0);
+	__this_cpu_write((fpuemustats).rint_d, 0);
+	__this_cpu_write((fpuemustats).round_w_s, 0);
+	__this_cpu_write((fpuemustats).round_w_d, 0);
+	__this_cpu_write((fpuemustats).round_l_s, 0);
+	__this_cpu_write((fpuemustats).round_l_d, 0);
+	__this_cpu_write((fpuemustats).rsqrt_s, 0);
+	__this_cpu_write((fpuemustats).rsqrt_d, 0);
+	__this_cpu_write((fpuemustats).sel_s, 0);
+	__this_cpu_write((fpuemustats).sel_d, 0);
+	__this_cpu_write((fpuemustats).seleqz_s, 0);
+	__this_cpu_write((fpuemustats).seleqz_d, 0);
+	__this_cpu_write((fpuemustats).selnez_s, 0);
+	__this_cpu_write((fpuemustats).selnez_d, 0);
+	__this_cpu_write((fpuemustats).sqrt_s, 0);
+	__this_cpu_write((fpuemustats).sqrt_d, 0);
+	__this_cpu_write((fpuemustats).sub_s, 0);
+	__this_cpu_write((fpuemustats).sub_d, 0);
+	__this_cpu_write((fpuemustats).trunc_w_s, 0);
+	__this_cpu_write((fpuemustats).trunc_w_d, 0);
+	__this_cpu_write((fpuemustats).trunc_l_s, 0);
+	__this_cpu_write((fpuemustats).trunc_l_d, 0);
+
+	return 0;
+}
+
+DEFINE_SHOW_ATTRIBUTE(fpuemustats_clear);
+
+static int __init debugfs_fpuemu(void)
+{
+	struct dentry *fpuemu_debugfs_base_dir;
+	struct dentry *fpuemu_debugfs_inst_dir;
+	char name[32];
+
+	fpuemu_debugfs_base_dir = debugfs_create_dir("fpuemustats",
+						     mips_debugfs_dir);
+
+	debugfs_create_file("fpuemustats_clear", 0444, mips_debugfs_dir, NULL,
+			    &fpuemustats_clear_fops);
+
+#define FPU_EMU_STAT_OFFSET(m)						\
+	offsetof(struct mips_fpu_emulator_stats, m)
+
+#define FPU_STAT_CREATE(m)						\
+do {									\
+	debugfs_create_file(#m, 0444, fpuemu_debugfs_base_dir,		\
+				(void *)FPU_EMU_STAT_OFFSET(m),		\
+				&fops_fpuemu_stat);			\
+} while (0)
+
+	FPU_STAT_CREATE(emulated);
+	FPU_STAT_CREATE(loads);
+	FPU_STAT_CREATE(stores);
+	FPU_STAT_CREATE(branches);
+	FPU_STAT_CREATE(cp1ops);
+	FPU_STAT_CREATE(cp1xops);
+	FPU_STAT_CREATE(errors);
+	FPU_STAT_CREATE(ieee754_inexact);
+	FPU_STAT_CREATE(ieee754_underflow);
+	FPU_STAT_CREATE(ieee754_overflow);
+	FPU_STAT_CREATE(ieee754_zerodiv);
+	FPU_STAT_CREATE(ieee754_invalidop);
+	FPU_STAT_CREATE(ds_emul);
+
+	fpuemu_debugfs_inst_dir = debugfs_create_dir("instructions",
+						     fpuemu_debugfs_base_dir);
+
+#define FPU_STAT_CREATE_EX(m)						\
+do {									\
+	adjust_instruction_counter_name(name, #m, sizeof(name));			\
+									\
+	debugfs_create_file(name, 0444, fpuemu_debugfs_inst_dir,	\
+				(void *)FPU_EMU_STAT_OFFSET(m),		\
+				&fops_fpuemu_stat);			\
+} while (0)
+
+	FPU_STAT_CREATE_EX(abs_s);
+	FPU_STAT_CREATE_EX(abs_d);
+	FPU_STAT_CREATE_EX(add_s);
+	FPU_STAT_CREATE_EX(add_d);
+	FPU_STAT_CREATE_EX(bc1eqz);
+	FPU_STAT_CREATE_EX(bc1nez);
+	FPU_STAT_CREATE_EX(ceil_w_s);
+	FPU_STAT_CREATE_EX(ceil_w_d);
+	FPU_STAT_CREATE_EX(ceil_l_s);
+	FPU_STAT_CREATE_EX(ceil_l_d);
+	FPU_STAT_CREATE_EX(class_s);
+	FPU_STAT_CREATE_EX(class_d);
+	FPU_STAT_CREATE_EX(cmp_af_s);
+	FPU_STAT_CREATE_EX(cmp_af_d);
+	FPU_STAT_CREATE_EX(cmp_eq_s);
+	FPU_STAT_CREATE_EX(cmp_eq_d);
+	FPU_STAT_CREATE_EX(cmp_le_s);
+	FPU_STAT_CREATE_EX(cmp_le_d);
+	FPU_STAT_CREATE_EX(cmp_lt_s);
+	FPU_STAT_CREATE_EX(cmp_lt_d);
+	FPU_STAT_CREATE_EX(cmp_ne_s);
+	FPU_STAT_CREATE_EX(cmp_ne_d);
+	FPU_STAT_CREATE_EX(cmp_or_s);
+	FPU_STAT_CREATE_EX(cmp_or_d);
+	FPU_STAT_CREATE_EX(cmp_ueq_s);
+	FPU_STAT_CREATE_EX(cmp_ueq_d);
+	FPU_STAT_CREATE_EX(cmp_ule_s);
+	FPU_STAT_CREATE_EX(cmp_ule_d);
+	FPU_STAT_CREATE_EX(cmp_ult_s);
+	FPU_STAT_CREATE_EX(cmp_ult_d);
+	FPU_STAT_CREATE_EX(cmp_un_s);
+	FPU_STAT_CREATE_EX(cmp_un_d);
+	FPU_STAT_CREATE_EX(cmp_une_s);
+	FPU_STAT_CREATE_EX(cmp_une_d);
+	FPU_STAT_CREATE_EX(cmp_saf_s);
+	FPU_STAT_CREATE_EX(cmp_saf_d);
+	FPU_STAT_CREATE_EX(cmp_seq_s);
+	FPU_STAT_CREATE_EX(cmp_seq_d);
+	FPU_STAT_CREATE_EX(cmp_sle_s);
+	FPU_STAT_CREATE_EX(cmp_sle_d);
+	FPU_STAT_CREATE_EX(cmp_slt_s);
+	FPU_STAT_CREATE_EX(cmp_slt_d);
+	FPU_STAT_CREATE_EX(cmp_sne_s);
+	FPU_STAT_CREATE_EX(cmp_sne_d);
+	FPU_STAT_CREATE_EX(cmp_sor_s);
+	FPU_STAT_CREATE_EX(cmp_sor_d);
+	FPU_STAT_CREATE_EX(cmp_sueq_s);
+	FPU_STAT_CREATE_EX(cmp_sueq_d);
+	FPU_STAT_CREATE_EX(cmp_sule_s);
+	FPU_STAT_CREATE_EX(cmp_sule_d);
+	FPU_STAT_CREATE_EX(cmp_sult_s);
+	FPU_STAT_CREATE_EX(cmp_sult_d);
+	FPU_STAT_CREATE_EX(cmp_sun_s);
+	FPU_STAT_CREATE_EX(cmp_sun_d);
+	FPU_STAT_CREATE_EX(cmp_sune_s);
+	FPU_STAT_CREATE_EX(cmp_sune_d);
+	FPU_STAT_CREATE_EX(cvt_d_l);
+	FPU_STAT_CREATE_EX(cvt_d_s);
+	FPU_STAT_CREATE_EX(cvt_d_w);
+	FPU_STAT_CREATE_EX(cvt_l_s);
+	FPU_STAT_CREATE_EX(cvt_l_d);
+	FPU_STAT_CREATE_EX(cvt_s_d);
+	FPU_STAT_CREATE_EX(cvt_s_l);
+	FPU_STAT_CREATE_EX(cvt_s_w);
+	FPU_STAT_CREATE_EX(cvt_w_s);
+	FPU_STAT_CREATE_EX(cvt_w_d);
+	FPU_STAT_CREATE_EX(div_s);
+	FPU_STAT_CREATE_EX(div_d);
+	FPU_STAT_CREATE_EX(floor_w_s);
+	FPU_STAT_CREATE_EX(floor_w_d);
+	FPU_STAT_CREATE_EX(floor_l_s);
+	FPU_STAT_CREATE_EX(floor_l_d);
+	FPU_STAT_CREATE_EX(maddf_s);
+	FPU_STAT_CREATE_EX(maddf_d);
+	FPU_STAT_CREATE_EX(max_s);
+	FPU_STAT_CREATE_EX(max_d);
+	FPU_STAT_CREATE_EX(maxa_s);
+	FPU_STAT_CREATE_EX(maxa_d);
+	FPU_STAT_CREATE_EX(min_s);
+	FPU_STAT_CREATE_EX(min_d);
+	FPU_STAT_CREATE_EX(mina_s);
+	FPU_STAT_CREATE_EX(mina_d);
+	FPU_STAT_CREATE_EX(mov_s);
+	FPU_STAT_CREATE_EX(mov_d);
+	FPU_STAT_CREATE_EX(msubf_s);
+	FPU_STAT_CREATE_EX(msubf_d);
+	FPU_STAT_CREATE_EX(mul_s);
+	FPU_STAT_CREATE_EX(mul_d);
+	FPU_STAT_CREATE_EX(neg_s);
+	FPU_STAT_CREATE_EX(neg_d);
+	FPU_STAT_CREATE_EX(recip_s);
+	FPU_STAT_CREATE_EX(recip_d);
+	FPU_STAT_CREATE_EX(rint_s);
+	FPU_STAT_CREATE_EX(rint_d);
+	FPU_STAT_CREATE_EX(round_w_s);
+	FPU_STAT_CREATE_EX(round_w_d);
+	FPU_STAT_CREATE_EX(round_l_s);
+	FPU_STAT_CREATE_EX(round_l_d);
+	FPU_STAT_CREATE_EX(rsqrt_s);
+	FPU_STAT_CREATE_EX(rsqrt_d);
+	FPU_STAT_CREATE_EX(sel_s);
+	FPU_STAT_CREATE_EX(sel_d);
+	FPU_STAT_CREATE_EX(seleqz_s);
+	FPU_STAT_CREATE_EX(seleqz_d);
+	FPU_STAT_CREATE_EX(selnez_s);
+	FPU_STAT_CREATE_EX(selnez_d);
+	FPU_STAT_CREATE_EX(sqrt_s);
+	FPU_STAT_CREATE_EX(sqrt_d);
+	FPU_STAT_CREATE_EX(sub_s);
+	FPU_STAT_CREATE_EX(sub_d);
+	FPU_STAT_CREATE_EX(trunc_w_s);
+	FPU_STAT_CREATE_EX(trunc_w_d);
+	FPU_STAT_CREATE_EX(trunc_l_s);
+	FPU_STAT_CREATE_EX(trunc_l_d);
+
+	return 0;
+}
+arch_initcall(debugfs_fpuemu);

linux-7.2.3/arch/mips/math-emu/sp_2008class.c

@@ -0,0 +1,52 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * IEEE754 floating point arithmetic
+ * single precision: CLASS.f
+ * FPR[fd] = class(FPR[fs])
+ *
+ * MIPS floating point support
+ * Copyright (C) 2015 Imagination Technologies, Ltd.
+ * Author: Markos Chandras <markos.chandras@imgtec.com>
+ */
+
+#include "ieee754sp.h"
+
+int ieee754sp_2008class(union ieee754sp x)
+{
+	COMPXSP;
+
+	EXPLODEXSP;
+
+	/*
+	 * 10 bit mask as follows:
+	 *
+	 * bit0 = SNAN
+	 * bit1 = QNAN
+	 * bit2 = -INF
+	 * bit3 = -NORM
+	 * bit4 = -DNORM
+	 * bit5 = -ZERO
+	 * bit6 = INF
+	 * bit7 = NORM
+	 * bit8 = DNORM
+	 * bit9 = ZERO
+	 */
+
+	switch(xc) {
+	case IEEE754_CLASS_SNAN:
+		return 0x01;
+	case IEEE754_CLASS_QNAN:
+		return 0x02;
+	case IEEE754_CLASS_INF:
+		return 0x04 << (xs ? 0 : 4);
+	case IEEE754_CLASS_NORM:
+		return 0x08 << (xs ? 0 : 4);
+	case IEEE754_CLASS_DNORM:
+		return 0x10 << (xs ? 0 : 4);
+	case IEEE754_CLASS_ZERO:
+		return 0x20 << (xs ? 0 : 4);
+	default:
+		pr_err("Unknown class: %d\n", xc);
+		return 0;
+	}
+}

linux-7.2.3/arch/mips/math-emu/sp_add.c

@@ -0,0 +1,164 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_add(union ieee754sp x, union ieee754sp y)
+{
+	int s;
+
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+
+	ieee754_clearcx();
+
+	FLUSHXSP;
+	FLUSHYSP;
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754sp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754sp_nanxcpt(x);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return x;
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		if (xs == ys)
+			return x;
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754sp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+		return x;
+
+	/*
+	 * Zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		if (xs == ys)
+			return x;
+		else
+			return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD);
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		SPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		break;
+	}
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+	/*
+	 * Provide guard, round and stick bit space.
+	 */
+	xm <<= 3;
+	ym <<= 3;
+
+	if (xe > ye) {
+		/*
+		 * Have to shift y fraction right to align.
+		 */
+		s = xe - ye;
+		ym = XSPSRS(ym, s);
+		ye += s;
+	} else if (ye > xe) {
+		/*
+		 * Have to shift x fraction right to align.
+		 */
+		s = ye - xe;
+		xm = XSPSRS(xm, s);
+		xe += s;
+	}
+	assert(xe == ye);
+	assert(xe <= SP_EMAX);
+
+	if (xs == ys) {
+		/*
+		 * Generate 28 bit result of adding two 27 bit numbers
+		 * leaving result in xm, xs and xe.
+		 */
+		xm = xm + ym;
+
+		if (xm >> (SP_FBITS + 1 + 3)) { /* carry out */
+			SPXSRSX1();
+		}
+	} else {
+		if (xm >= ym) {
+			xm = xm - ym;
+		} else {
+			xm = ym - xm;
+			xs = ys;
+		}
+		if (xm == 0)
+			return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD);
+
+		/*
+		 * Normalize in extended single precision
+		 */
+		while ((xm >> (SP_FBITS + 3)) == 0) {
+			xm <<= 1;
+			xe--;
+		}
+	}
+
+	return ieee754sp_format(xs, xe, xm);
+}

linux-7.2.3/arch/mips/math-emu/sp_cmp.c

@@ -0,0 +1,47 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+int ieee754sp_cmp(union ieee754sp x, union ieee754sp y, int cmp, int sig)
+{
+	int vx;
+	int vy;
+
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+	FLUSHXSP;
+	FLUSHYSP;
+	ieee754_clearcx();	/* Even clear inexact flag here */
+
+	if (ieee754_class_nan(xc) || ieee754_class_nan(yc)) {
+		if (sig ||
+		    xc == IEEE754_CLASS_SNAN || yc == IEEE754_CLASS_SNAN)
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return (cmp & IEEE754_CUN) != 0;
+	} else {
+		vx = x.bits;
+		vy = y.bits;
+
+		if (vx < 0)
+			vx = -vx ^ SP_SIGN_BIT;
+		if (vy < 0)
+			vy = -vy ^ SP_SIGN_BIT;
+
+		if (vx < vy)
+			return (cmp & IEEE754_CLT) != 0;
+		else if (vx == vy)
+			return (cmp & IEEE754_CEQ) != 0;
+		else
+			return (cmp & IEEE754_CGT) != 0;
+	}
+}

linux-7.2.3/arch/mips/math-emu/sp_div.c

@@ -0,0 +1,142 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_div(union ieee754sp x, union ieee754sp y)
+{
+	unsigned int rm;
+	int re;
+	unsigned int bm;
+
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+
+	ieee754_clearcx();
+
+	FLUSHXSP;
+	FLUSHYSP;
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754sp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754sp_nanxcpt(x);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return x;
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754sp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+		return ieee754sp_zero(xs ^ ys);
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+		return ieee754sp_inf(xs ^ ys);
+
+	/*
+	 * Zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754sp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		ieee754_setcx(IEEE754_ZERO_DIVIDE);
+		return ieee754sp_inf(xs ^ ys);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return ieee754sp_zero(xs == ys ? 0 : 1);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		SPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		break;
+	}
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+	/* provide rounding space */
+	xm <<= 3;
+	ym <<= 3;
+
+	/* now the dirty work */
+
+	rm = 0;
+	re = xe - ye;
+
+	for (bm = SP_MBIT(SP_FBITS + 2); bm; bm >>= 1) {
+		if (xm >= ym) {
+			xm -= ym;
+			rm |= bm;
+			if (xm == 0)
+				break;
+		}
+		xm <<= 1;
+	}
+
+	rm <<= 1;
+	if (xm)
+		rm |= 1;	/* have remainder, set sticky */
+
+	assert(rm);
+
+	/* normalise rm to rounding precision ?
+	 */
+	while ((rm >> (SP_FBITS + 3)) == 0) {
+		rm <<= 1;
+		re--;
+	}
+
+	return ieee754sp_format(xs == ys ? 0 : 1, re, rm);
+}

linux-7.2.3/arch/mips/math-emu/sp_fdp.c

@@ -0,0 +1,73 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+#include "ieee754dp.h"
+
+static inline union ieee754sp ieee754sp_nan_fdp(int xs, u64 xm)
+{
+	return buildsp(xs, SP_EMAX + 1 + SP_EBIAS,
+		       xm >> (DP_FBITS - SP_FBITS));
+}
+
+union ieee754sp ieee754sp_fdp(union ieee754dp x)
+{
+	union ieee754sp y;
+	u32 rm;
+
+	COMPXDP;
+	COMPYSP;
+
+	EXPLODEXDP;
+
+	ieee754_clearcx();
+
+	FLUSHXDP;
+
+	switch (xc) {
+	case IEEE754_CLASS_SNAN:
+		x = ieee754dp_nanxcpt(x);
+		EXPLODEXDP;
+		fallthrough;
+	case IEEE754_CLASS_QNAN:
+		y = ieee754sp_nan_fdp(xs, xm);
+		if (!ieee754_csr.nan2008) {
+			EXPLODEYSP;
+			if (!ieee754_class_nan(yc))
+				y = ieee754sp_indef();
+		}
+		return y;
+
+	case IEEE754_CLASS_INF:
+		return ieee754sp_inf(xs);
+
+	case IEEE754_CLASS_ZERO:
+		return ieee754sp_zero(xs);
+
+	case IEEE754_CLASS_DNORM:
+		/* can't possibly be sp representable */
+		ieee754_setcx(IEEE754_UNDERFLOW);
+		ieee754_setcx(IEEE754_INEXACT);
+		if ((ieee754_csr.rm == FPU_CSR_RU && !xs) ||
+				(ieee754_csr.rm == FPU_CSR_RD && xs))
+			return ieee754sp_mind(xs);
+		return ieee754sp_zero(xs);
+
+	case IEEE754_CLASS_NORM:
+		break;
+	}
+
+	/*
+	 * Convert from DP_FBITS to SP_FBITS+3 with sticky right shift.
+	 */
+	rm = (xm >> (DP_FBITS - (SP_FBITS + 3))) |
+	     ((xm << (64 - (DP_FBITS - (SP_FBITS + 3)))) != 0);
+
+	return ieee754sp_format(xs, xe, rm);
+}

linux-7.2.3/arch/mips/math-emu/sp_fint.c

@@ -0,0 +1,53 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_fint(int x)
+{
+	unsigned int xm;
+	int xe;
+	int xs;
+
+	ieee754_clearcx();
+
+	if (x == 0)
+		return ieee754sp_zero(0);
+	if (x == 1 || x == -1)
+		return ieee754sp_one(x < 0);
+	if (x == 10 || x == -10)
+		return ieee754sp_ten(x < 0);
+
+	xs = (x < 0);
+	if (xs) {
+		if (x == (1 << 31))
+			xm = ((unsigned) 1 << 31);	/* max neg can't be safely negated */
+		else
+			xm = -x;
+	} else {
+		xm = x;
+	}
+	xe = SP_FBITS + 3;
+
+	if (xm >> (SP_FBITS + 1 + 3)) {
+		/* shunt out overflow bits
+		 */
+		while (xm >> (SP_FBITS + 1 + 3)) {
+			SPXSRSX1();
+		}
+	} else {
+		/* normalize in grs extended single precision
+		 */
+		while ((xm >> (SP_FBITS + 3)) == 0) {
+			xm <<= 1;
+			xe--;
+		}
+	}
+	return ieee754sp_format(xs, xe, xm);
+}

linux-7.2.3/arch/mips/math-emu/sp_flong.c

@@ -0,0 +1,52 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_flong(s64 x)
+{
+	u64 xm;		/* <--- need 64-bit mantissa temp */
+	int xe;
+	int xs;
+
+	ieee754_clearcx();
+
+	if (x == 0)
+		return ieee754sp_zero(0);
+	if (x == 1 || x == -1)
+		return ieee754sp_one(x < 0);
+	if (x == 10 || x == -10)
+		return ieee754sp_ten(x < 0);
+
+	xs = (x < 0);
+	if (xs) {
+		if (x == (1ULL << 63))
+			xm = (1ULL << 63);	/* max neg can't be safely negated */
+		else
+			xm = -x;
+	} else {
+		xm = x;
+	}
+	xe = SP_FBITS + 3;
+
+	if (xm >> (SP_FBITS + 1 + 3)) {
+		/* shunt out overflow bits
+		 */
+		while (xm >> (SP_FBITS + 1 + 3)) {
+			SPXSRSX1();
+		}
+	} else {
+		/* normalize in grs extended single precision */
+		while ((xm >> (SP_FBITS + 3)) == 0) {
+			xm <<= 1;
+			xe--;
+		}
+	}
+	return ieee754sp_format(xs, xe, xm);
+}

linux-7.2.3/arch/mips/math-emu/sp_fmax.c

@@ -0,0 +1,252 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * IEEE754 floating point arithmetic
+ * single precision: MAX{,A}.f
+ * MAX : Scalar Floating-Point Maximum
+ * MAXA: Scalar Floating-Point argument with Maximum Absolute Value
+ *
+ * MAX.S : FPR[fd] = maxNum(FPR[fs],FPR[ft])
+ * MAXA.S: FPR[fd] = maxNumMag(FPR[fs],FPR[ft])
+ *
+ * MIPS floating point support
+ * Copyright (C) 2015 Imagination Technologies, Ltd.
+ * Author: Markos Chandras <markos.chandras@imgtec.com>
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_fmax(union ieee754sp x, union ieee754sp y)
+{
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+
+	FLUSHXSP;
+	FLUSHYSP;
+
+	ieee754_clearcx();
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754sp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754sp_nanxcpt(x);
+
+	/*
+	 * Quiet NaN handling
+	 */
+
+	/*
+	 *    The case of both inputs quiet NaNs
+	 */
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+		return x;
+
+	/*
+	 *    The cases of exactly one input quiet NaN (numbers
+	 *    are here preferred as returned values to NaNs)
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return y;
+
+	/*
+	 * Infinity and zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return xs ? y : x;
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return ys ? x : y;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		return ieee754sp_zero(xs & ys);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		SPDNORMX;
+	}
+
+	/* Finally get to do some computation */
+
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+	/* Compare signs */
+	if (xs > ys)
+		return y;
+	else if (xs < ys)
+		return x;
+
+	/* Signs of inputs are equal, let's compare exponents */
+	if (xs == 0) {
+		/* Inputs are both positive */
+		if (xe > ye)
+			return x;
+		else if (xe < ye)
+			return y;
+	} else {
+		/* Inputs are both negative */
+		if (xe > ye)
+			return y;
+		else if (xe < ye)
+			return x;
+	}
+
+	/* Signs and exponents of inputs are equal, let's compare mantissas */
+	if (xs == 0) {
+		/* Inputs are both positive, with equal signs and exponents */
+		if (xm <= ym)
+			return y;
+		return x;
+	}
+	/* Inputs are both negative, with equal signs and exponents */
+	if (xm <= ym)
+		return x;
+	return y;
+}
+
+union ieee754sp ieee754sp_fmaxa(union ieee754sp x, union ieee754sp y)
+{
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+
+	FLUSHXSP;
+	FLUSHYSP;
+
+	ieee754_clearcx();
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754sp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754sp_nanxcpt(x);
+
+	/*
+	 * Quiet NaN handling
+	 */
+
+	/*
+	 *    The case of both inputs quiet NaNs
+	 */
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+		return x;
+
+	/*
+	 *    The cases of exactly one input quiet NaN (numbers
+	 *    are here preferred as returned values to NaNs)
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return y;
+
+	/*
+	 * Infinity and zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		return ieee754sp_inf(xs & ys);
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		return ieee754sp_zero(xs & ys);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		SPDNORMX;
+	}
+
+	/* Finally get to do some computation */
+
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+	/* Compare exponent */
+	if (xe > ye)
+		return x;
+	else if (xe < ye)
+		return y;
+
+	/* Compare mantissa */
+	if (xm < ym)
+		return y;
+	else if (xm > ym)
+		return x;
+	else if (xs == 0)
+		return x;
+	return y;
+}

linux-7.2.3/arch/mips/math-emu/sp_fmin.c

@@ -0,0 +1,252 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * IEEE754 floating point arithmetic
+ * single precision: MIN{,A}.f
+ * MIN : Scalar Floating-Point Minimum
+ * MINA: Scalar Floating-Point argument with Minimum Absolute Value
+ *
+ * MIN.S : FPR[fd] = minNum(FPR[fs],FPR[ft])
+ * MINA.S: FPR[fd] = maxNumMag(FPR[fs],FPR[ft])
+ *
+ * MIPS floating point support
+ * Copyright (C) 2015 Imagination Technologies, Ltd.
+ * Author: Markos Chandras <markos.chandras@imgtec.com>
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_fmin(union ieee754sp x, union ieee754sp y)
+{
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+
+	FLUSHXSP;
+	FLUSHYSP;
+
+	ieee754_clearcx();
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754sp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754sp_nanxcpt(x);
+
+	/*
+	 * Quiet NaN handling
+	 */
+
+	/*
+	 *    The case of both inputs quiet NaNs
+	 */
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+		return x;
+
+	/*
+	 *    The cases of exactly one input quiet NaN (numbers
+	 *    are here preferred as returned values to NaNs)
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return y;
+
+	/*
+	 * Infinity and zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return xs ? x : y;
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return ys ? y : x;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		return ieee754sp_zero(xs | ys);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		SPDNORMX;
+	}
+
+	/* Finally get to do some computation */
+
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+	/* Compare signs */
+	if (xs > ys)
+		return x;
+	else if (xs < ys)
+		return y;
+
+	/* Signs of inputs are the same, let's compare exponents */
+	if (xs == 0) {
+		/* Inputs are both positive */
+		if (xe > ye)
+			return y;
+		else if (xe < ye)
+			return x;
+	} else {
+		/* Inputs are both negative */
+		if (xe > ye)
+			return x;
+		else if (xe < ye)
+			return y;
+	}
+
+	/* Signs and exponents of inputs are equal, let's compare mantissas */
+	if (xs == 0) {
+		/* Inputs are both positive, with equal signs and exponents */
+		if (xm <= ym)
+			return x;
+		return y;
+	}
+	/* Inputs are both negative, with equal signs and exponents */
+	if (xm <= ym)
+		return y;
+	return x;
+}
+
+union ieee754sp ieee754sp_fmina(union ieee754sp x, union ieee754sp y)
+{
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+
+	FLUSHXSP;
+	FLUSHYSP;
+
+	ieee754_clearcx();
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754sp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754sp_nanxcpt(x);
+
+	/*
+	 * Quiet NaN handling
+	 */
+
+	/*
+	 *    The case of both inputs quiet NaNs
+	 */
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+		return x;
+
+	/*
+	 *    The cases of exactly one input quiet NaN (numbers
+	 *    are here preferred as returned values to NaNs)
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return y;
+
+	/*
+	 * Infinity and zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		return ieee754sp_inf(xs | ys);
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		return ieee754sp_zero(xs | ys);
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		SPDNORMX;
+	}
+
+	/* Finally get to do some computation */
+
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+	/* Compare exponent */
+	if (xe > ye)
+		return y;
+	else if (xe < ye)
+		return x;
+
+	/* Compare mantissa */
+	if (xm < ym)
+		return x;
+	else if (xm > ym)
+		return y;
+	else if (xs == 1)
+		return x;
+	return y;
+}

linux-7.2.3/arch/mips/math-emu/sp_maddf.c

@@ -0,0 +1,278 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/*
+ * IEEE754 floating point arithmetic
+ * single precision: MADDF.f (Fused Multiply Add)
+ * MADDF.fmt: FPR[fd] = FPR[fd] + (FPR[fs] x FPR[ft])
+ *
+ * MIPS floating point support
+ * Copyright (C) 2015 Imagination Technologies, Ltd.
+ * Author: Markos Chandras <markos.chandras@imgtec.com>
+ */
+
+#include "ieee754sp.h"
+
+
+static union ieee754sp _sp_maddf(union ieee754sp z, union ieee754sp x,
+				 union ieee754sp y, enum maddf_flags flags)
+{
+	int re;
+	int rs;
+	unsigned int rm;
+	u64 rm64;
+	u64 zm64;
+	int s;
+
+	COMPXSP;
+	COMPYSP;
+	COMPZSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+	EXPLODEZSP;
+
+	FLUSHXSP;
+	FLUSHYSP;
+	FLUSHZSP;
+
+	ieee754_clearcx();
+
+	rs = xs ^ ys;
+	if (flags & MADDF_NEGATE_PRODUCT)
+		rs ^= 1;
+	if (flags & MADDF_NEGATE_ADDITION)
+		zs ^= 1;
+
+	/*
+	 * Handle the cases when at least one of x, y or z is a NaN.
+	 * Order of precedence is sNaN, qNaN and z, x, y.
+	 */
+	if (zc == IEEE754_CLASS_SNAN)
+		return ieee754sp_nanxcpt(z);
+	if (xc == IEEE754_CLASS_SNAN)
+		return ieee754sp_nanxcpt(x);
+	if (yc == IEEE754_CLASS_SNAN)
+		return ieee754sp_nanxcpt(y);
+	if (zc == IEEE754_CLASS_QNAN)
+		return z;
+	if (xc == IEEE754_CLASS_QNAN)
+		return x;
+	if (yc == IEEE754_CLASS_QNAN)
+		return y;
+
+	if (zc == IEEE754_CLASS_DNORM)
+		SPDNORMZ;
+	/* ZERO z cases are handled separately below */
+
+	switch (CLPAIR(xc, yc)) {
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754sp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		if ((zc == IEEE754_CLASS_INF) && (zs != rs)) {
+			/*
+			 * Cases of addition of infinities with opposite signs
+			 * or subtraction of infinities with same signs.
+			 */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754sp_indef();
+		}
+		/*
+		 * z is here either not an infinity, or an infinity having the
+		 * same sign as product (x*y). The result must be an infinity,
+		 * and its sign is determined only by the sign of product (x*y).
+		 */
+		return ieee754sp_inf(rs);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		if (zc == IEEE754_CLASS_INF)
+			return ieee754sp_inf(zs);
+		if (zc == IEEE754_CLASS_ZERO) {
+			/* Handle cases +0 + (-0) and similar ones. */
+			if (zs == rs)
+				/*
+				 * Cases of addition of zeros of equal signs
+				 * or subtraction of zeroes of opposite signs.
+				 * The sign of the resulting zero is in any
+				 * such case determined only by the sign of z.
+				 */
+				return z;
+
+			return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD);
+		}
+		/* x*y is here 0, and z is not 0, so just return z */
+		return z;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		if (zc == IEEE754_CLASS_INF)
+			return ieee754sp_inf(zs);
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		if (zc == IEEE754_CLASS_INF)
+			return ieee754sp_inf(zs);
+		SPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		if (zc == IEEE754_CLASS_INF)
+			return ieee754sp_inf(zs);
+		/* continue to real computations */
+	}
+
+	/* Finally get to do some computation */
+
+	/*
+	 * Do the multiplication bit first
+	 *
+	 * rm = xm * ym, re = xe + ye basically
+	 *
+	 * At this point xm and ym should have been normalized.
+	 */
+
+	/* rm = xm * ym, re = xe+ye basically */
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+	re = xe + ye;
+
+	/* Multiple 24 bit xm and ym to give 48 bit results */
+	rm64 = (uint64_t)xm * ym;
+
+	/* Shunt to top of word */
+	rm64 = rm64 << 16;
+
+	/* Put explicit bit at bit 62 if necessary */
+	if ((int64_t) rm64 < 0) {
+		rm64 = rm64 >> 1;
+		re++;
+	}
+
+	assert(rm64 & (1 << 62));
+
+	if (zc == IEEE754_CLASS_ZERO) {
+		/*
+		 * Move explicit bit from bit 62 to bit 26 since the
+		 * ieee754sp_format code expects the mantissa to be
+		 * 27 bits wide (24 + 3 rounding bits).
+		 */
+		rm = XSPSRS64(rm64, (62 - 26));
+		return ieee754sp_format(rs, re, rm);
+	}
+
+	/* Move explicit bit from bit 23 to bit 62 */
+	zm64 = (uint64_t)zm << (62 - 23);
+	assert(zm64 & (1 << 62));
+
+	/* Make the exponents the same */
+	if (ze > re) {
+		/*
+		 * Have to shift r fraction right to align.
+		 */
+		s = ze - re;
+		rm64 = XSPSRS64(rm64, s);
+		re += s;
+	} else if (re > ze) {
+		/*
+		 * Have to shift z fraction right to align.
+		 */
+		s = re - ze;
+		zm64 = XSPSRS64(zm64, s);
+		ze += s;
+	}
+	assert(ze == re);
+	assert(ze <= SP_EMAX);
+
+	/* Do the addition */
+	if (zs == rs) {
+		/*
+		 * Generate 64 bit result by adding two 63 bit numbers
+		 * leaving result in zm64, zs and ze.
+		 */
+		zm64 = zm64 + rm64;
+		if ((int64_t)zm64 < 0) {	/* carry out */
+			zm64 = XSPSRS1(zm64);
+			ze++;
+		}
+	} else {
+		if (zm64 >= rm64) {
+			zm64 = zm64 - rm64;
+		} else {
+			zm64 = rm64 - zm64;
+			zs = rs;
+		}
+		if (zm64 == 0)
+			return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD);
+
+		/*
+		 * Put explicit bit at bit 62 if necessary.
+		 */
+		while ((zm64 >> 62) == 0) {
+			zm64 <<= 1;
+			ze--;
+		}
+	}
+
+	/*
+	 * Move explicit bit from bit 62 to bit 26 since the
+	 * ieee754sp_format code expects the mantissa to be
+	 * 27 bits wide (24 + 3 rounding bits).
+	 */
+	zm = XSPSRS64(zm64, (62 - 26));
+
+	return ieee754sp_format(zs, ze, zm);
+}
+
+union ieee754sp ieee754sp_maddf(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y)
+{
+	return _sp_maddf(z, x, y, 0);
+}
+
+union ieee754sp ieee754sp_msubf(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y)
+{
+	return _sp_maddf(z, x, y, MADDF_NEGATE_PRODUCT);
+}
+
+union ieee754sp ieee754sp_madd(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y)
+{
+	return _sp_maddf(z, x, y, 0);
+}
+
+union ieee754sp ieee754sp_msub(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y)
+{
+	return _sp_maddf(z, x, y, MADDF_NEGATE_ADDITION);
+}
+
+union ieee754sp ieee754sp_nmadd(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y)
+{
+	return _sp_maddf(z, x, y, MADDF_NEGATE_PRODUCT|MADDF_NEGATE_ADDITION);
+}
+
+union ieee754sp ieee754sp_nmsub(union ieee754sp z, union ieee754sp x,
+				union ieee754sp y)
+{
+	return _sp_maddf(z, x, y, MADDF_NEGATE_PRODUCT);
+}

linux-7.2.3/arch/mips/math-emu/sp_mul.c

@@ -0,0 +1,154 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_mul(union ieee754sp x, union ieee754sp y)
+{
+	int re;
+	int rs;
+	unsigned int rm;
+	unsigned short lxm;
+	unsigned short hxm;
+	unsigned short lym;
+	unsigned short hym;
+	unsigned int lrm;
+	unsigned int hrm;
+	unsigned int t;
+	unsigned int at;
+
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+
+	ieee754_clearcx();
+
+	FLUSHXSP;
+	FLUSHYSP;
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754sp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754sp_nanxcpt(x);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return x;
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754sp_indef();
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		return ieee754sp_inf(xs ^ ys);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return ieee754sp_zero(xs ^ ys);
+
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		SPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		break;
+	}
+	/* rm = xm * ym, re = xe+ye basically */
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+	re = xe + ye;
+	rs = xs ^ ys;
+
+	/* shunt to top of word */
+	xm <<= 32 - (SP_FBITS + 1);
+	ym <<= 32 - (SP_FBITS + 1);
+
+	/*
+	 * Multiply 32 bits xm, ym to give high 32 bits rm with stickness.
+	 */
+	lxm = xm & 0xffff;
+	hxm = xm >> 16;
+	lym = ym & 0xffff;
+	hym = ym >> 16;
+
+	lrm = lxm * lym;	/* 16 * 16 => 32 */
+	hrm = hxm * hym;	/* 16 * 16 => 32 */
+
+	t = lxm * hym; /* 16 * 16 => 32 */
+	at = lrm + (t << 16);
+	hrm += at < lrm;
+	lrm = at;
+	hrm = hrm + (t >> 16);
+
+	t = hxm * lym; /* 16 * 16 => 32 */
+	at = lrm + (t << 16);
+	hrm += at < lrm;
+	lrm = at;
+	hrm = hrm + (t >> 16);
+
+	rm = hrm | (lrm != 0);
+
+	/*
+	 * Sticky shift down to normal rounding precision.
+	 */
+	if ((int) rm < 0) {
+		rm = (rm >> (32 - (SP_FBITS + 1 + 3))) |
+		    ((rm << (SP_FBITS + 1 + 3)) != 0);
+		re++;
+	} else {
+		rm = (rm >> (32 - (SP_FBITS + 1 + 3 + 1))) |
+		     ((rm << (SP_FBITS + 1 + 3 + 1)) != 0);
+	}
+	assert(rm & (SP_HIDDEN_BIT << 3));
+
+	return ieee754sp_format(rs, re, rm);
+}

linux-7.2.3/arch/mips/math-emu/sp_rint.c

@@ -0,0 +1,79 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ * Copyright (C) 2017 Imagination Technologies, Ltd.
+ * Author: Aleksandar Markovic <aleksandar.markovic@imgtec.com>
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_rint(union ieee754sp x)
+{
+	union ieee754sp ret;
+	u32 residue;
+	int sticky;
+	int round;
+	int odd;
+
+	COMPXDP;		/* <-- DP needed for 64-bit mantissa tmp */
+
+	ieee754_clearcx();
+
+	EXPLODEXSP;
+	FLUSHXSP;
+
+	if (xc == IEEE754_CLASS_SNAN)
+		return ieee754sp_nanxcpt(x);
+
+	if ((xc == IEEE754_CLASS_QNAN) ||
+	    (xc == IEEE754_CLASS_INF) ||
+	    (xc == IEEE754_CLASS_ZERO))
+		return x;
+
+	if (xe >= SP_FBITS)
+		return x;
+
+	if (xe < -1) {
+		residue = xm;
+		round = 0;
+		sticky = residue != 0;
+		xm = 0;
+	} else {
+		residue = xm << (xe + 1);
+		residue <<= 31 - SP_FBITS;
+		round = (residue >> 31) != 0;
+		sticky = (residue << 1) != 0;
+		xm >>= SP_FBITS - xe;
+	}
+
+	odd = (xm & 0x1) != 0x0;
+
+	switch (ieee754_csr.rm) {
+	case FPU_CSR_RN:	/* toward nearest */
+		if (round && (sticky || odd))
+			xm++;
+		break;
+	case FPU_CSR_RZ:	/* toward zero */
+		break;
+	case FPU_CSR_RU:	/* toward +infinity */
+		if ((round || sticky) && !xs)
+			xm++;
+		break;
+	case FPU_CSR_RD:	/* toward -infinity */
+		if ((round || sticky) && xs)
+			xm++;
+		break;
+	}
+
+	if (round || sticky)
+		ieee754_setcx(IEEE754_INEXACT);
+
+	ret = ieee754sp_flong(xm);
+	SPSIGN(ret) = xs;
+
+	return ret;
+}

linux-7.2.3/arch/mips/math-emu/sp_simple.c

@@ -0,0 +1,49 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_neg(union ieee754sp x)
+{
+	union ieee754sp y;
+
+	if (ieee754_csr.abs2008) {
+		y = x;
+		SPSIGN(y) = !SPSIGN(x);
+	} else {
+		unsigned int oldrm;
+
+		oldrm = ieee754_csr.rm;
+		ieee754_csr.rm = FPU_CSR_RD;
+		y = ieee754sp_sub(ieee754sp_zero(0), x);
+		ieee754_csr.rm = oldrm;
+	}
+	return y;
+}
+
+union ieee754sp ieee754sp_abs(union ieee754sp x)
+{
+	union ieee754sp y;
+
+	if (ieee754_csr.abs2008) {
+		y = x;
+		SPSIGN(y) = 0;
+	} else {
+		unsigned int oldrm;
+
+		oldrm = ieee754_csr.rm;
+		ieee754_csr.rm = FPU_CSR_RD;
+		if (SPSIGN(x))
+			y = ieee754sp_sub(ieee754sp_zero(0), x);
+		else
+			y = ieee754sp_add(ieee754sp_zero(0), x);
+		ieee754_csr.rm = oldrm;
+	}
+	return y;
+}

linux-7.2.3/arch/mips/math-emu/sp_sqrt.c

@@ -0,0 +1,103 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision square root
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_sqrt(union ieee754sp x)
+{
+	int ix, s, q, m, t, i;
+	unsigned int r;
+	COMPXSP;
+
+	/* take care of Inf and NaN */
+
+	EXPLODEXSP;
+	ieee754_clearcx();
+	FLUSHXSP;
+
+	/* x == INF or NAN? */
+	switch (xc) {
+	case IEEE754_CLASS_SNAN:
+		return ieee754sp_nanxcpt(x);
+
+	case IEEE754_CLASS_QNAN:
+		/* sqrt(Nan) = Nan */
+		return x;
+
+	case IEEE754_CLASS_ZERO:
+		/* sqrt(0) = 0 */
+		return x;
+
+	case IEEE754_CLASS_INF:
+		if (xs) {
+			/* sqrt(-Inf) = Nan */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754sp_indef();
+		}
+		/* sqrt(+Inf) = Inf */
+		return x;
+
+	case IEEE754_CLASS_DNORM:
+	case IEEE754_CLASS_NORM:
+		if (xs) {
+			/* sqrt(-x) = Nan */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754sp_indef();
+		}
+		break;
+	}
+
+	ix = x.bits;
+
+	/* normalize x */
+	m = (ix >> 23);
+	if (m == 0) {		/* subnormal x */
+		for (i = 0; (ix & 0x00800000) == 0; i++)
+			ix <<= 1;
+		m -= i - 1;
+	}
+	m -= 127;		/* unbias exponent */
+	ix = (ix & 0x007fffff) | 0x00800000;
+	if (m & 1)		/* odd m, double x to make it even */
+		ix += ix;
+	m >>= 1;		/* m = [m/2] */
+
+	/* generate sqrt(x) bit by bit */
+	ix += ix;
+	s = 0;
+	q = 0;			/* q = sqrt(x) */
+	r = 0x01000000;		/* r = moving bit from right to left */
+
+	while (r != 0) {
+		t = s + r;
+		if (t <= ix) {
+			s = t + r;
+			ix -= t;
+			q += r;
+		}
+		ix += ix;
+		r >>= 1;
+	}
+
+	if (ix != 0) {
+		ieee754_setcx(IEEE754_INEXACT);
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RU:
+			q += 2;
+			break;
+		case FPU_CSR_RN:
+			q += (q & 1);
+			break;
+		}
+	}
+	ix = (q >> 1) + 0x3f000000;
+	ix += (m << 23);
+	x.bits = ix;
+	return x;
+}

linux-7.2.3/arch/mips/math-emu/sp_sub.c

@@ -0,0 +1,168 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+union ieee754sp ieee754sp_sub(union ieee754sp x, union ieee754sp y)
+{
+	int s;
+
+	COMPXSP;
+	COMPYSP;
+
+	EXPLODEXSP;
+	EXPLODEYSP;
+
+	ieee754_clearcx();
+
+	FLUSHXSP;
+	FLUSHYSP;
+
+	switch (CLPAIR(xc, yc)) {
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
+		return ieee754sp_nanxcpt(y);
+
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
+		return ieee754sp_nanxcpt(x);
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
+	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
+		return x;
+
+
+	/*
+	 * Infinity handling
+	 */
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
+		if (xs != ys)
+			return x;
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754sp_indef();
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
+		return ieee754sp_inf(ys ^ 1);
+
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
+		return x;
+
+	/*
+	 * Zero handling
+	 */
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
+		if (xs != ys)
+			return x;
+		else
+			return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD);
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
+		return x;
+
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
+	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
+		/* quick fix up */
+		SPSIGN(y) ^= 1;
+		return y;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
+		SPDNORMX;
+		fallthrough;
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
+		SPDNORMY;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
+		SPDNORMX;
+		break;
+
+	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
+		break;
+	}
+	/* flip sign of y and handle as add */
+	ys ^= 1;
+
+	assert(xm & SP_HIDDEN_BIT);
+	assert(ym & SP_HIDDEN_BIT);
+
+
+	/* provide guard,round and stick bit space */
+	xm <<= 3;
+	ym <<= 3;
+
+	if (xe > ye) {
+		/*
+		 * have to shift y fraction right to align
+		 */
+		s = xe - ye;
+		ym = XSPSRS(ym, s);
+		ye += s;
+	} else if (ye > xe) {
+		/*
+		 * have to shift x fraction right to align
+		 */
+		s = ye - xe;
+		xm = XSPSRS(xm, s);
+		xe += s;
+	}
+	assert(xe == ye);
+	assert(xe <= SP_EMAX);
+
+	if (xs == ys) {
+		/* generate 28 bit result of adding two 27 bit numbers
+		 */
+		xm = xm + ym;
+
+		if (xm >> (SP_FBITS + 1 + 3)) { /* carry out */
+			SPXSRSX1();	/* shift preserving sticky */
+		}
+	} else {
+		if (xm >= ym) {
+			xm = xm - ym;
+		} else {
+			xm = ym - xm;
+			xs = ys;
+		}
+		if (xm == 0) {
+			if (ieee754_csr.rm == FPU_CSR_RD)
+				return ieee754sp_zero(1);	/* round negative inf. => sign = -1 */
+			else
+				return ieee754sp_zero(0);	/* other round modes   => sign = 1 */
+		}
+		/* normalize to rounding precision
+		 */
+		while ((xm >> (SP_FBITS + 3)) == 0) {
+			xm <<= 1;
+			xe--;
+		}
+	}
+
+	return ieee754sp_format(xs, xe, xm);
+}

linux-7.2.3/arch/mips/math-emu/sp_tint.c

@@ -0,0 +1,100 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+int ieee754sp_tint(union ieee754sp x)
+{
+	u32 residue;
+	int round;
+	int sticky;
+	int odd;
+
+	COMPXSP;
+
+	ieee754_clearcx();
+
+	EXPLODEXSP;
+	FLUSHXSP;
+
+	switch (xc) {
+	case IEEE754_CLASS_SNAN:
+	case IEEE754_CLASS_QNAN:
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754si_indef();
+
+	case IEEE754_CLASS_INF:
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754si_overflow(xs);
+
+	case IEEE754_CLASS_ZERO:
+		return 0;
+
+	case IEEE754_CLASS_DNORM:
+	case IEEE754_CLASS_NORM:
+		break;
+	}
+	if (xe >= 31) {
+		/* look for valid corner case */
+		if (xe == 31 && xs && xm == SP_HIDDEN_BIT)
+			return -0x80000000;
+		/* Set invalid. We will only use overflow for floating
+		   point overflow */
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754si_overflow(xs);
+	}
+	/* oh gawd */
+	if (xe > SP_FBITS) {
+		xm <<= xe - SP_FBITS;
+	} else {
+		if (xe < -1) {
+			residue = xm;
+			round = 0;
+			sticky = residue != 0;
+			xm = 0;
+		} else {
+			/* Shifting a u32 32 times does not work,
+			* so we do it in two steps. Be aware that xe
+			* may be -1 */
+			residue = xm << (xe + 1);
+			residue <<= 31 - SP_FBITS;
+			round = (residue >> 31) != 0;
+			sticky = (residue << 1) != 0;
+			xm >>= SP_FBITS - xe;
+		}
+		odd = (xm & 0x1) != 0x0;
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RN:
+			if (round && (sticky || odd))
+				xm++;
+			break;
+		case FPU_CSR_RZ:
+			break;
+		case FPU_CSR_RU:	/* toward +Infinity */
+			if ((round || sticky) && !xs)
+				xm++;
+			break;
+		case FPU_CSR_RD:	/* toward -Infinity */
+			if ((round || sticky) && xs)
+				xm++;
+			break;
+		}
+		if ((xm >> 31) != 0) {
+			/* This can happen after rounding */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754si_overflow(xs);
+		}
+		if (round || sticky)
+			ieee754_setcx(IEEE754_INEXACT);
+	}
+	if (xs)
+		return -xm;
+	else
+		return xm;
+}

linux-7.2.3/arch/mips/math-emu/sp_tlong.c

@@ -0,0 +1,96 @@
+// SPDX-License-Identifier: GPL-2.0-only
+/* IEEE754 floating point arithmetic
+ * single precision
+ */
+/*
+ * MIPS floating point support
+ * Copyright (C) 1994-2000 Algorithmics Ltd.
+ */
+
+#include "ieee754sp.h"
+
+s64 ieee754sp_tlong(union ieee754sp x)
+{
+	u32 residue;
+	int round;
+	int sticky;
+	int odd;
+
+	COMPXDP;		/* <-- need 64-bit mantissa tmp */
+
+	ieee754_clearcx();
+
+	EXPLODEXSP;
+	FLUSHXSP;
+
+	switch (xc) {
+	case IEEE754_CLASS_SNAN:
+	case IEEE754_CLASS_QNAN:
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754di_indef();
+
+	case IEEE754_CLASS_INF:
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754di_overflow(xs);
+
+	case IEEE754_CLASS_ZERO:
+		return 0;
+
+	case IEEE754_CLASS_DNORM:
+	case IEEE754_CLASS_NORM:
+		break;
+	}
+	if (xe >= 63) {
+		/* look for valid corner case */
+		if (xe == 63 && xs && xm == SP_HIDDEN_BIT)
+			return -0x8000000000000000LL;
+		/* Set invalid. We will only use overflow for floating
+		   point overflow */
+		ieee754_setcx(IEEE754_INVALID_OPERATION);
+		return ieee754di_overflow(xs);
+	}
+	/* oh gawd */
+	if (xe > SP_FBITS) {
+		xm <<= xe - SP_FBITS;
+	} else if (xe < SP_FBITS) {
+		if (xe < -1) {
+			residue = xm;
+			round = 0;
+			sticky = residue != 0;
+			xm = 0;
+		} else {
+			residue = xm << (32 - SP_FBITS + xe);
+			round = (residue >> 31) != 0;
+			sticky = (residue << 1) != 0;
+			xm >>= SP_FBITS - xe;
+		}
+		odd = (xm & 0x1) != 0x0;
+		switch (ieee754_csr.rm) {
+		case FPU_CSR_RN:
+			if (round && (sticky || odd))
+				xm++;
+			break;
+		case FPU_CSR_RZ:
+			break;
+		case FPU_CSR_RU:	/* toward +Infinity */
+			if ((round || sticky) && !xs)
+				xm++;
+			break;
+		case FPU_CSR_RD:	/* toward -Infinity */
+			if ((round || sticky) && xs)
+				xm++;
+			break;
+		}
+		if ((xm >> 63) != 0) {
+			/* This can happen after rounding */
+			ieee754_setcx(IEEE754_INVALID_OPERATION);
+			return ieee754di_overflow(xs);
+		}
+		if (round || sticky)
+			ieee754_setcx(IEEE754_INEXACT);
+	}
+	if (xs)
+		return -xm;
+	else
+		return xm;
+}