From e09b41010ba33a20a87472ee821fa407a5b8da36 Mon Sep 17 00:00:00 2001 From: José Pekkarinen Date: Mon, 11 Apr 2016 10:41:07 +0300 Subject: These changes are the raw update to linux-4.4.6-rt14. Kernel sources are taken from kernel.org, and rt patch from the rt wiki download page. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit During the rebasing, the following patch collided: Force tick interrupt and get rid of softirq magic(I70131fb85). Collisions have been removed because its logic was found on the source already. Change-Id: I7f57a4081d9deaa0d9ccfc41a6c8daccdee3b769 Signed-off-by: José Pekkarinen --- kernel/arch/mips/math-emu/Makefile | 4 +- kernel/arch/mips/math-emu/cp1emu.c | 380 ++++++++++++++++++++++++++++++- kernel/arch/mips/math-emu/dp_2008class.c | 55 +++++ kernel/arch/mips/math-emu/dp_fmax.c | 213 +++++++++++++++++ kernel/arch/mips/math-emu/dp_fmin.c | 213 +++++++++++++++++ kernel/arch/mips/math-emu/dp_maddf.c | 265 +++++++++++++++++++++ kernel/arch/mips/math-emu/dp_msubf.c | 269 ++++++++++++++++++++++ kernel/arch/mips/math-emu/dsemul.c | 1 - kernel/arch/mips/math-emu/ieee754.h | 19 ++ kernel/arch/mips/math-emu/ieee754int.h | 8 +- kernel/arch/mips/math-emu/me-debugfs.c | 4 +- kernel/arch/mips/math-emu/sp_2008class.c | 55 +++++ kernel/arch/mips/math-emu/sp_fmax.c | 213 +++++++++++++++++ kernel/arch/mips/math-emu/sp_fmin.c | 213 +++++++++++++++++ kernel/arch/mips/math-emu/sp_maddf.c | 255 +++++++++++++++++++++ kernel/arch/mips/math-emu/sp_msubf.c | 258 +++++++++++++++++++++ 16 files changed, 2404 insertions(+), 21 deletions(-) create mode 100644 kernel/arch/mips/math-emu/dp_2008class.c create mode 100644 kernel/arch/mips/math-emu/dp_fmax.c create mode 100644 kernel/arch/mips/math-emu/dp_fmin.c create mode 100644 kernel/arch/mips/math-emu/dp_maddf.c create mode 100644 kernel/arch/mips/math-emu/dp_msubf.c create mode 100644 kernel/arch/mips/math-emu/sp_2008class.c create mode 100644 kernel/arch/mips/math-emu/sp_fmax.c create mode 100644 kernel/arch/mips/math-emu/sp_fmin.c create mode 100644 kernel/arch/mips/math-emu/sp_maddf.c create mode 100644 kernel/arch/mips/math-emu/sp_msubf.c (limited to 'kernel/arch/mips/math-emu') diff --git a/kernel/arch/mips/math-emu/Makefile b/kernel/arch/mips/math-emu/Makefile index 2e5f96275..a19641d3a 100644 --- a/kernel/arch/mips/math-emu/Makefile +++ b/kernel/arch/mips/math-emu/Makefile @@ -4,9 +4,9 @@ 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_tint.o dp_fint.o dp_maddf.o dp_msubf.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_tint.o sp_fint.o sp_maddf.o sp_msubf.o sp_2008class.o sp_fmin.o sp_fmax.o \ dsemul.o lib-y += ieee754d.o \ diff --git a/kernel/arch/mips/math-emu/cp1emu.c b/kernel/arch/mips/math-emu/cp1emu.c index 2b95e34fa..32f0e19a0 100644 --- a/kernel/arch/mips/math-emu/cp1emu.c +++ b/kernel/arch/mips/math-emu/cp1emu.c @@ -551,7 +551,7 @@ static int isBranchInstr(struct pt_regs *regs, struct mm_decoded_insn dec_insn, dec_insn.next_pc_inc; return 1; case blezl_op: - if (NO_R6EMU) + if (!insn.i_format.rt && NO_R6EMU) break; case blez_op: @@ -588,7 +588,7 @@ static int isBranchInstr(struct pt_regs *regs, struct mm_decoded_insn dec_insn, dec_insn.next_pc_inc; return 1; case bgtzl_op: - if (NO_R6EMU) + if (!insn.i_format.rt && NO_R6EMU) break; case bgtz_op: /* @@ -1394,6 +1394,14 @@ static const unsigned char cmptab[8] = { IEEE754_CLT | IEEE754_CEQ | IEEE754_CUN, /* cmp_ule (sig) cmp_ngt */ }; +static const unsigned char negative_cmptab[8] = { + 0, /* Reserved */ + IEEE754_CLT | IEEE754_CGT | IEEE754_CEQ, + IEEE754_CLT | IEEE754_CGT | IEEE754_CUN, + IEEE754_CLT | IEEE754_CGT, + /* Reserved */ +}; + /* * Additional MIPS4 instructions @@ -1735,6 +1743,126 @@ static int fpu_emu(struct pt_regs *xcp, struct mips_fpu_struct *ctx, SPFROMREG(rv.s, MIPSInst_FS(ir)); break; + case fseleqz_op: + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(rv.s, MIPSInst_FT(ir)); + if (rv.w & 0x1) + rv.w = 0; + else + SPFROMREG(rv.s, MIPSInst_FS(ir)); + break; + + case fselnez_op: + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(rv.s, MIPSInst_FT(ir)); + if (rv.w & 0x1) + SPFROMREG(rv.s, MIPSInst_FS(ir)); + else + rv.w = 0; + break; + + case fmaddf_op: { + union ieee754sp ft, fs, fd; + + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(ft, MIPSInst_FT(ir)); + SPFROMREG(fs, MIPSInst_FS(ir)); + SPFROMREG(fd, MIPSInst_FD(ir)); + rv.s = ieee754sp_maddf(fd, fs, ft); + break; + } + + case fmsubf_op: { + union ieee754sp ft, fs, fd; + + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(ft, MIPSInst_FT(ir)); + SPFROMREG(fs, MIPSInst_FS(ir)); + SPFROMREG(fd, MIPSInst_FD(ir)); + rv.s = ieee754sp_msubf(fd, fs, ft); + break; + } + + case frint_op: { + union ieee754sp fs; + + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.l = ieee754sp_tlong(fs); + rv.s = ieee754sp_flong(rv.l); + goto copcsr; + } + + case fclass_op: { + union ieee754sp fs; + + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.w = ieee754sp_2008class(fs); + rfmt = w_fmt; + break; + } + + case fmin_op: { + union ieee754sp fs, ft; + + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(ft, MIPSInst_FT(ir)); + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.s = ieee754sp_fmin(fs, ft); + break; + } + + case fmina_op: { + union ieee754sp fs, ft; + + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(ft, MIPSInst_FT(ir)); + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.s = ieee754sp_fmina(fs, ft); + break; + } + + case fmax_op: { + union ieee754sp fs, ft; + + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(ft, MIPSInst_FT(ir)); + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.s = ieee754sp_fmax(fs, ft); + break; + } + + case fmaxa_op: { + union ieee754sp fs, ft; + + if (!cpu_has_mips_r6) + return SIGILL; + + SPFROMREG(ft, MIPSInst_FT(ir)); + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.s = ieee754sp_fmaxa(fs, ft); + break; + } + case fabs_op: handler.u = ieee754sp_abs; goto scopuop; @@ -1838,7 +1966,7 @@ copcsr: goto copcsr; default: - if (MIPSInst_FUNC(ir) >= fcmp_op) { + if (!NO_R6EMU && MIPSInst_FUNC(ir) >= fcmp_op) { unsigned cmpop = MIPSInst_FUNC(ir) - fcmp_op; union ieee754sp fs, ft; @@ -1932,6 +2060,127 @@ copcsr: return 0; DPFROMREG(rv.d, MIPSInst_FS(ir)); break; + + case fseleqz_op: + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(rv.d, MIPSInst_FT(ir)); + if (rv.l & 0x1) + rv.l = 0; + else + DPFROMREG(rv.d, MIPSInst_FS(ir)); + break; + + case fselnez_op: + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(rv.d, MIPSInst_FT(ir)); + if (rv.l & 0x1) + DPFROMREG(rv.d, MIPSInst_FS(ir)); + else + rv.l = 0; + break; + + case fmaddf_op: { + union ieee754dp ft, fs, fd; + + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(ft, MIPSInst_FT(ir)); + DPFROMREG(fs, MIPSInst_FS(ir)); + DPFROMREG(fd, MIPSInst_FD(ir)); + rv.d = ieee754dp_maddf(fd, fs, ft); + break; + } + + case fmsubf_op: { + union ieee754dp ft, fs, fd; + + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(ft, MIPSInst_FT(ir)); + DPFROMREG(fs, MIPSInst_FS(ir)); + DPFROMREG(fd, MIPSInst_FD(ir)); + rv.d = ieee754dp_msubf(fd, fs, ft); + break; + } + + case frint_op: { + union ieee754dp fs; + + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.l = ieee754dp_tlong(fs); + rv.d = ieee754dp_flong(rv.l); + goto copcsr; + } + + case fclass_op: { + union ieee754dp fs; + + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.w = ieee754dp_2008class(fs); + rfmt = w_fmt; + break; + } + + case fmin_op: { + union ieee754dp fs, ft; + + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(ft, MIPSInst_FT(ir)); + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.d = ieee754dp_fmin(fs, ft); + break; + } + + case fmina_op: { + union ieee754dp fs, ft; + + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(ft, MIPSInst_FT(ir)); + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.d = ieee754dp_fmina(fs, ft); + break; + } + + case fmax_op: { + union ieee754dp fs, ft; + + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(ft, MIPSInst_FT(ir)); + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.d = ieee754dp_fmax(fs, ft); + break; + } + + case fmaxa_op: { + union ieee754dp fs, ft; + + if (!cpu_has_mips_r6) + return SIGILL; + + DPFROMREG(ft, MIPSInst_FT(ir)); + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.d = ieee754dp_fmaxa(fs, ft); + break; + } + case fabs_op: handler.u = ieee754dp_abs; goto dcopuop; @@ -2015,7 +2264,7 @@ dcopuop: goto copcsr; default: - if (MIPSInst_FUNC(ir) >= fcmp_op) { + if (!NO_R6EMU && MIPSInst_FUNC(ir) >= fcmp_op) { unsigned cmpop = MIPSInst_FUNC(ir) - fcmp_op; union ieee754dp fs, ft; @@ -2039,8 +2288,11 @@ dcopuop: break; } break; + } + + case w_fmt: { + union ieee754dp fs; - case w_fmt: switch (MIPSInst_FUNC(ir)) { case fcvts_op: /* convert word to single precision real */ @@ -2054,10 +2306,65 @@ dcopuop: rv.d = ieee754dp_fint(fs.bits); rfmt = d_fmt; goto copcsr; - default: - return SIGILL; + default: { + /* Emulating the new CMP.condn.fmt R6 instruction */ +#define CMPOP_MASK 0x7 +#define SIGN_BIT (0x1 << 3) +#define PREDICATE_BIT (0x1 << 4) + + int cmpop = MIPSInst_FUNC(ir) & CMPOP_MASK; + int sig = MIPSInst_FUNC(ir) & SIGN_BIT; + union ieee754sp fs, ft; + + /* This is an R6 only instruction */ + if (!cpu_has_mips_r6 || + (MIPSInst_FUNC(ir) & 0x20)) + return SIGILL; + + /* fmt is w_fmt for single precision so fix it */ + rfmt = s_fmt; + /* default to false */ + rv.w = 0; + + /* CMP.condn.S */ + SPFROMREG(fs, MIPSInst_FS(ir)); + SPFROMREG(ft, MIPSInst_FT(ir)); + + /* positive predicates */ + if (!(MIPSInst_FUNC(ir) & PREDICATE_BIT)) { + if (ieee754sp_cmp(fs, ft, cmptab[cmpop], + sig)) + rv.w = -1; /* true, all 1s */ + if ((sig) && + ieee754_cxtest(IEEE754_INVALID_OPERATION)) + rcsr = FPU_CSR_INV_X | FPU_CSR_INV_S; + else + goto copcsr; + } else { + /* negative predicates */ + switch (cmpop) { + case 1: + case 2: + case 3: + if (ieee754sp_cmp(fs, ft, + negative_cmptab[cmpop], + sig)) + rv.w = -1; /* true, all 1s */ + if (sig && + ieee754_cxtest(IEEE754_INVALID_OPERATION)) + rcsr = FPU_CSR_INV_X | FPU_CSR_INV_S; + else + goto copcsr; + break; + default: + /* Reserved R6 ops */ + pr_err("Reserved MIPS R6 CMP.condn.S operation\n"); + return SIGILL; + } + } + break; + } } - break; } case l_fmt: @@ -2078,11 +2385,60 @@ dcopuop: rv.d = ieee754dp_flong(bits); rfmt = d_fmt; goto copcsr; - default: - return SIGILL; - } - break; + default: { + /* Emulating the new CMP.condn.fmt R6 instruction */ + int cmpop = MIPSInst_FUNC(ir) & CMPOP_MASK; + int sig = MIPSInst_FUNC(ir) & SIGN_BIT; + union ieee754dp fs, ft; + + if (!cpu_has_mips_r6 || + (MIPSInst_FUNC(ir) & 0x20)) + return SIGILL; + + /* fmt is l_fmt for double precision so fix it */ + rfmt = d_fmt; + /* default to false */ + rv.l = 0; + + /* CMP.condn.D */ + DPFROMREG(fs, MIPSInst_FS(ir)); + DPFROMREG(ft, MIPSInst_FT(ir)); + /* positive predicates */ + if (!(MIPSInst_FUNC(ir) & PREDICATE_BIT)) { + if (ieee754dp_cmp(fs, ft, + cmptab[cmpop], sig)) + rv.l = -1LL; /* true, all 1s */ + if (sig && + ieee754_cxtest(IEEE754_INVALID_OPERATION)) + rcsr = FPU_CSR_INV_X | FPU_CSR_INV_S; + else + goto copcsr; + } else { + /* negative predicates */ + switch (cmpop) { + case 1: + case 2: + case 3: + if (ieee754dp_cmp(fs, ft, + negative_cmptab[cmpop], + sig)) + rv.l = -1LL; /* true, all 1s */ + if (sig && + ieee754_cxtest(IEEE754_INVALID_OPERATION)) + rcsr = FPU_CSR_INV_X | FPU_CSR_INV_S; + else + goto copcsr; + break; + default: + /* Reserved R6 ops */ + pr_err("Reserved MIPS R6 CMP.condn.D operation\n"); + return SIGILL; + } + } + break; + } + } default: return SIGILL; } diff --git a/kernel/arch/mips/math-emu/dp_2008class.c b/kernel/arch/mips/math-emu/dp_2008class.c new file mode 100644 index 000000000..9dc39fc48 --- /dev/null +++ b/kernel/arch/mips/math-emu/dp_2008class.c @@ -0,0 +1,55 @@ +/* + * 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 + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#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; + } +} diff --git a/kernel/arch/mips/math-emu/dp_fmax.c b/kernel/arch/mips/math-emu/dp_fmax.c new file mode 100644 index 000000000..fd71b8daa --- /dev/null +++ b/kernel/arch/mips/math-emu/dp_fmax.c @@ -0,0 +1,213 @@ +/* + * 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 + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#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); + + /* numbers are preferred 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_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 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): + if (xs == ys) + return x; + return ieee754dp_zero(1); + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + DPDNORMX; + + 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; + + /* Compare exponent */ + if (xe > ye) + return x; + else if (xe < ye) + return y; + + /* Compare mantissa */ + if (xm <= ym) + return y; + return x; +} + +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); + + /* numbers are preferred 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_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 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 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 y; + + case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): + if (xs == ys) + return x; + return ieee754dp_zero(1); + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + DPDNORMX; + + 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; + return x; +} diff --git a/kernel/arch/mips/math-emu/dp_fmin.c b/kernel/arch/mips/math-emu/dp_fmin.c new file mode 100644 index 000000000..c1072b0df --- /dev/null +++ b/kernel/arch/mips/math-emu/dp_fmin.c @@ -0,0 +1,213 @@ +/* + * 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 + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#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); + + /* numbers are preferred 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_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 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): + if (xs == ys) + return x; + return ieee754dp_zero(1); + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + DPDNORMX; + + 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; + + /* Compare exponent */ + if (xe > ye) + return y; + else if (xe < ye) + return x; + + /* Compare mantissa */ + if (xm <= ym) + return x; + return y; +} + +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); + + /* numbers are preferred 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_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 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 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 y; + + case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): + if (xs == ys) + return x; + return ieee754dp_zero(1); + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + DPDNORMX; + + 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; + return y; +} diff --git a/kernel/arch/mips/math-emu/dp_maddf.c b/kernel/arch/mips/math-emu/dp_maddf.c new file mode 100644 index 000000000..119eda9fa --- /dev/null +++ b/kernel/arch/mips/math-emu/dp_maddf.c @@ -0,0 +1,265 @@ +/* + * 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 + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#include "ieee754dp.h" + +union ieee754dp ieee754dp_maddf(union ieee754dp z, union ieee754dp x, + union ieee754dp y) +{ + int re; + int rs; + u64 rm; + unsigned lxm; + unsigned hxm; + unsigned lym; + unsigned hym; + u64 lrm; + u64 hrm; + u64 t; + u64 at; + int s; + + COMPXDP; + COMPYDP; + + u64 zm; int ze; int zs __maybe_unused; int zc; + + EXPLODEXDP; + EXPLODEYDP; + EXPLODEDP(z, zc, zs, ze, zm) + + FLUSHXDP; + FLUSHYDP; + FLUSHDP(z, zc, zs, ze, zm); + + ieee754_clearcx(); + + switch (zc) { + case IEEE754_CLASS_SNAN: + ieee754_setcx(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(z); + case IEEE754_CLASS_DNORM: + DPDNORMx(zm, ze); + /* QNAN is handled separately below */ + } + + 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): + if (zc == IEEE754_CLASS_QNAN) + return z; + 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_QNAN) + return z; + 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): + if (zc == IEEE754_CLASS_INF) + return ieee754dp_inf(zs); + /* Multiplication is 0 so just return z */ + return z; + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + DPDNORMX; + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754dp_inf(zs); + DPDNORMY; + break; + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754dp_inf(zs); + DPDNORMX; + break; + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754dp_inf(zs); + /* fall through 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; + rs = xs ^ ys; + + /* shunt to top of word */ + xm <<= 64 - (DP_FBITS + 1); + ym <<= 64 - (DP_FBITS + 1); + + /* + * Multiply 32 bits xm, ym to give high 32 bits rm with stickness. + */ + + /* 32 * 32 => 64 */ +#define DPXMULT(x, y) ((u64)(x) * (u64)y) + + 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)); + + /* And now the addition */ + assert(zm & DP_HIDDEN_BIT); + + /* + * Provide guard,round and stick bit space. + */ + zm <<= 3; + + if (ze > re) { + /* + * Have to shift y fraction right to align. + */ + s = ze - re; + rm = XDPSRS(rm, s); + re += s; + } else if (re > ze) { + /* + * Have to shift x fraction right to align. + */ + s = re - ze; + zm = XDPSRS(zm, s); + ze += s; + } + assert(ze == re); + assert(ze <= DP_EMAX); + + if (zs == rs) { + /* + * Generate 28 bit result of adding two 27 bit numbers + * leaving result in xm, xs and xe. + */ + zm = zm + rm; + + if (zm >> (DP_FBITS + 1 + 3)) { /* carry out */ + zm = XDPSRS1(zm); + ze++; + } + } else { + if (zm >= rm) { + zm = zm - rm; + } else { + zm = rm - zm; + zs = rs; + } + if (zm == 0) + return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD); + + /* + * Normalize to rounding precision. + */ + while ((zm >> (DP_FBITS + 3)) == 0) { + zm <<= 1; + ze--; + } + } + + return ieee754dp_format(zs, ze, zm); +} diff --git a/kernel/arch/mips/math-emu/dp_msubf.c b/kernel/arch/mips/math-emu/dp_msubf.c new file mode 100644 index 000000000..12241262f --- /dev/null +++ b/kernel/arch/mips/math-emu/dp_msubf.c @@ -0,0 +1,269 @@ +/* + * IEEE754 floating point arithmetic + * double precision: MSUB.f (Fused Multiply Subtract) + * MSUBF.fmt: FPR[fd] = FPR[fd] - (FPR[fs] x FPR[ft]) + * + * MIPS floating point support + * Copyright (C) 2015 Imagination Technologies, Ltd. + * Author: Markos Chandras + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#include "ieee754dp.h" + +union ieee754dp ieee754dp_msubf(union ieee754dp z, union ieee754dp x, + union ieee754dp y) +{ + int re; + int rs; + u64 rm; + unsigned lxm; + unsigned hxm; + unsigned lym; + unsigned hym; + u64 lrm; + u64 hrm; + u64 t; + u64 at; + int s; + + COMPXDP; + COMPYDP; + + u64 zm; int ze; int zs __maybe_unused; int zc; + + EXPLODEXDP; + EXPLODEYDP; + EXPLODEDP(z, zc, zs, ze, zm) + + FLUSHXDP; + FLUSHYDP; + FLUSHDP(z, zc, zs, ze, zm); + + ieee754_clearcx(); + + switch (zc) { + case IEEE754_CLASS_SNAN: + ieee754_setcx(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(z); + case IEEE754_CLASS_DNORM: + DPDNORMx(zm, ze); + /* QNAN is handled separately below */ + } + + 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): + if (zc == IEEE754_CLASS_QNAN) + return z; + 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_QNAN) + return z; + 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): + if (zc == IEEE754_CLASS_INF) + return ieee754dp_inf(zs); + /* Multiplication is 0 so just return z */ + return z; + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + DPDNORMX; + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754dp_inf(zs); + DPDNORMY; + break; + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754dp_inf(zs); + DPDNORMX; + break; + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754dp_inf(zs); + /* fall through 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; + rs = xs ^ ys; + + /* shunt to top of word */ + xm <<= 64 - (DP_FBITS + 1); + ym <<= 64 - (DP_FBITS + 1); + + /* + * Multiply 32 bits xm, ym to give high 32 bits rm with stickness. + */ + + /* 32 * 32 => 64 */ +#define DPXMULT(x, y) ((u64)(x) * (u64)y) + + 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)); + + /* And now the subtraction */ + + /* flip sign of r and handle as add */ + rs ^= 1; + + assert(zm & DP_HIDDEN_BIT); + + /* + * Provide guard,round and stick bit space. + */ + zm <<= 3; + + if (ze > re) { + /* + * Have to shift y fraction right to align. + */ + s = ze - re; + rm = XDPSRS(rm, s); + re += s; + } else if (re > ze) { + /* + * Have to shift x fraction right to align. + */ + s = re - ze; + zm = XDPSRS(zm, s); + ze += s; + } + assert(ze == re); + assert(ze <= DP_EMAX); + + if (zs == rs) { + /* + * Generate 28 bit result of adding two 27 bit numbers + * leaving result in xm, xs and xe. + */ + zm = zm + rm; + + if (zm >> (DP_FBITS + 1 + 3)) { /* carry out */ + zm = XDPSRS1(zm); + ze++; + } + } else { + if (zm >= rm) { + zm = zm - rm; + } else { + zm = rm - zm; + zs = rs; + } + if (zm == 0) + return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD); + + /* + * Normalize to rounding precision. + */ + while ((zm >> (DP_FBITS + 3)) == 0) { + zm <<= 1; + ze--; + } + } + + return ieee754dp_format(zs, ze, zm); +} diff --git a/kernel/arch/mips/math-emu/dsemul.c b/kernel/arch/mips/math-emu/dsemul.c index e0b5cc27d..cbb36c14b 100644 --- a/kernel/arch/mips/math-emu/dsemul.c +++ b/kernel/arch/mips/math-emu/dsemul.c @@ -33,7 +33,6 @@ struct emuframe { int mips_dsemul(struct pt_regs *regs, mips_instruction ir, unsigned long cpc) { - extern asmlinkage void handle_dsemulret(void); struct emuframe __user *fr; int err; diff --git a/kernel/arch/mips/math-emu/ieee754.h b/kernel/arch/mips/math-emu/ieee754.h index a5ca108ce..df9472071 100644 --- a/kernel/arch/mips/math-emu/ieee754.h +++ b/kernel/arch/mips/math-emu/ieee754.h @@ -75,6 +75,16 @@ 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); +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) */ @@ -99,6 +109,15 @@ 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); +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 diff --git a/kernel/arch/mips/math-emu/ieee754int.h b/kernel/arch/mips/math-emu/ieee754int.h index 05389d5e3..6383e2c5c 100644 --- a/kernel/arch/mips/math-emu/ieee754int.h +++ b/kernel/arch/mips/math-emu/ieee754int.h @@ -65,8 +65,8 @@ static inline int ieee754_class_nan(int xc) vc = IEEE754_CLASS_INF; \ else if (vm & SP_MBIT(SP_FBITS-1)) \ vc = IEEE754_CLASS_SNAN; \ - else \ - vc = IEEE754_CLASS_QNAN; \ + else \ + vc = IEEE754_CLASS_QNAN; \ } else if (ve == SP_EMIN-1+SP_EBIAS) { \ if (vm) { \ ve = SP_EMIN; \ @@ -105,8 +105,8 @@ static inline int ieee754_class_nan(int xc) if (vm) { \ ve = DP_EMIN; \ vc = IEEE754_CLASS_DNORM; \ - } else \ - vc = IEEE754_CLASS_ZERO; \ + } else \ + vc = IEEE754_CLASS_ZERO; \ } else { \ ve -= DP_EBIAS; \ vm |= DP_HIDDEN_BIT; \ diff --git a/kernel/arch/mips/math-emu/me-debugfs.c b/kernel/arch/mips/math-emu/me-debugfs.c index f308e0f05..be650ed7d 100644 --- a/kernel/arch/mips/math-emu/me-debugfs.c +++ b/kernel/arch/mips/math-emu/me-debugfs.c @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -27,7 +28,6 @@ static int fpuemu_stat_get(void *data, u64 *val) } DEFINE_SIMPLE_ATTRIBUTE(fops_fpuemu_stat, fpuemu_stat_get, NULL, "%llu\n"); -extern struct dentry *mips_debugfs_dir; static int __init debugfs_fpuemu(void) { struct dentry *d, *dir; @@ -65,4 +65,4 @@ do { \ return 0; } -__initcall(debugfs_fpuemu); +arch_initcall(debugfs_fpuemu); diff --git a/kernel/arch/mips/math-emu/sp_2008class.c b/kernel/arch/mips/math-emu/sp_2008class.c new file mode 100644 index 000000000..ff62606a1 --- /dev/null +++ b/kernel/arch/mips/math-emu/sp_2008class.c @@ -0,0 +1,55 @@ +/* + * 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 + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#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; + } +} diff --git a/kernel/arch/mips/math-emu/sp_fmax.c b/kernel/arch/mips/math-emu/sp_fmax.c new file mode 100644 index 000000000..4d000844e --- /dev/null +++ b/kernel/arch/mips/math-emu/sp_fmax.c @@ -0,0 +1,213 @@ +/* + * 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 + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#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); + + /* numbers are preferred 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_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 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): + if (xs == ys) + return x; + return ieee754sp_zero(1); + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + SPDNORMX; + + 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; + + /* Compare exponent */ + if (xe > ye) + return x; + else if (xe < ye) + return y; + + /* Compare mantissa */ + if (xm <= ym) + return y; + return x; +} + +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); + + /* numbers are preferred 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_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 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 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 y; + + case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): + if (xs == ys) + return x; + return ieee754sp_zero(1); + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + SPDNORMX; + + 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; + return x; +} diff --git a/kernel/arch/mips/math-emu/sp_fmin.c b/kernel/arch/mips/math-emu/sp_fmin.c new file mode 100644 index 000000000..4eb1bb9e9 --- /dev/null +++ b/kernel/arch/mips/math-emu/sp_fmin.c @@ -0,0 +1,213 @@ +/* + * 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 + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#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); + + /* numbers are preferred 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_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 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): + if (xs == ys) + return x; + return ieee754sp_zero(1); + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + SPDNORMX; + + 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; + + /* Compare exponent */ + if (xe > ye) + return y; + else if (xe < ye) + return x; + + /* Compare mantissa */ + if (xm <= ym) + return x; + return y; +} + +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); + + /* numbers are preferred 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_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 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 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 y; + + case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): + if (xs == ys) + return x; + return ieee754sp_zero(1); + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + SPDNORMX; + + 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; + return y; +} diff --git a/kernel/arch/mips/math-emu/sp_maddf.c b/kernel/arch/mips/math-emu/sp_maddf.c new file mode 100644 index 000000000..dd1dd83e3 --- /dev/null +++ b/kernel/arch/mips/math-emu/sp_maddf.c @@ -0,0 +1,255 @@ +/* + * 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 + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#include "ieee754sp.h" + +union ieee754sp ieee754sp_maddf(union ieee754sp z, union ieee754sp x, + union ieee754sp y) +{ + int re; + int rs; + unsigned rm; + unsigned short lxm; + unsigned short hxm; + unsigned short lym; + unsigned short hym; + unsigned lrm; + unsigned hrm; + unsigned t; + unsigned at; + int s; + + COMPXSP; + COMPYSP; + u32 zm; int ze; int zs __maybe_unused; int zc; + + EXPLODEXSP; + EXPLODEYSP; + EXPLODESP(z, zc, zs, ze, zm) + + FLUSHXSP; + FLUSHYSP; + FLUSHSP(z, zc, zs, ze, zm); + + ieee754_clearcx(); + + switch (zc) { + case IEEE754_CLASS_SNAN: + ieee754_setcx(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(z); + case IEEE754_CLASS_DNORM: + SPDNORMx(zm, ze); + /* QNAN is handled separately below */ + } + + 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): + if (zc == IEEE754_CLASS_QNAN) + return z; + 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_QNAN) + return z; + 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): + if (zc == IEEE754_CLASS_INF) + return ieee754sp_inf(zs); + /* Multiplication is 0 so just return z */ + return z; + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + SPDNORMX; + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754sp_inf(zs); + SPDNORMY; + break; + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754sp_inf(zs); + SPDNORMX; + break; + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754sp_inf(zs); + /* fall through 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; + 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)); + + /* And now the addition */ + + assert(zm & SP_HIDDEN_BIT); + + /* + * Provide guard,round and stick bit space. + */ + zm <<= 3; + + if (ze > re) { + /* + * Have to shift y fraction right to align. + */ + s = ze - re; + SPXSRSYn(s); + } else if (re > ze) { + /* + * Have to shift x fraction right to align. + */ + s = re - ze; + SPXSRSYn(s); + } + assert(ze == re); + assert(ze <= SP_EMAX); + + if (zs == rs) { + /* + * Generate 28 bit result of adding two 27 bit numbers + * leaving result in zm, zs and ze. + */ + zm = zm + rm; + + if (zm >> (SP_FBITS + 1 + 3)) { /* carry out */ + SPXSRSX1(); + } + } else { + if (zm >= rm) { + zm = zm - rm; + } else { + zm = rm - zm; + zs = rs; + } + if (zm == 0) + return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD); + + /* + * Normalize in extended single precision + */ + while ((zm >> (SP_MBITS + 3)) == 0) { + zm <<= 1; + ze--; + } + + } + return ieee754sp_format(zs, ze, zm); +} diff --git a/kernel/arch/mips/math-emu/sp_msubf.c b/kernel/arch/mips/math-emu/sp_msubf.c new file mode 100644 index 000000000..81c38b980 --- /dev/null +++ b/kernel/arch/mips/math-emu/sp_msubf.c @@ -0,0 +1,258 @@ +/* + * IEEE754 floating point arithmetic + * single precision: MSUB.f (Fused Multiply Subtract) + * MSUBF.fmt: FPR[fd] = FPR[fd] - (FPR[fs] x FPR[ft]) + * + * MIPS floating point support + * Copyright (C) 2015 Imagination Technologies, Ltd. + * Author: Markos Chandras + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License as published by the + * Free Software Foundation; version 2 of the License. + */ + +#include "ieee754sp.h" + +union ieee754sp ieee754sp_msubf(union ieee754sp z, union ieee754sp x, + union ieee754sp y) +{ + int re; + int rs; + unsigned rm; + unsigned short lxm; + unsigned short hxm; + unsigned short lym; + unsigned short hym; + unsigned lrm; + unsigned hrm; + unsigned t; + unsigned at; + int s; + + COMPXSP; + COMPYSP; + u32 zm; int ze; int zs __maybe_unused; int zc; + + EXPLODEXSP; + EXPLODEYSP; + EXPLODESP(z, zc, zs, ze, zm) + + FLUSHXSP; + FLUSHYSP; + FLUSHSP(z, zc, zs, ze, zm); + + ieee754_clearcx(); + + switch (zc) { + case IEEE754_CLASS_SNAN: + ieee754_setcx(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(z); + case IEEE754_CLASS_DNORM: + SPDNORMx(zm, ze); + /* QNAN is handled separately below */ + } + + 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): + if (zc == IEEE754_CLASS_QNAN) + return z; + 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_QNAN) + return z; + 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): + if (zc == IEEE754_CLASS_INF) + return ieee754sp_inf(zs); + /* Multiplication is 0 so just return z */ + return z; + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): + SPDNORMX; + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754sp_inf(zs); + SPDNORMY; + break; + + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754sp_inf(zs); + SPDNORMX; + break; + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): + if (zc == IEEE754_CLASS_QNAN) + return z; + else if (zc == IEEE754_CLASS_INF) + return ieee754sp_inf(zs); + /* fall through to real compuation */ + } + + /* 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; + 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)); + + /* And now the subtraction */ + + /* Flip sign of r and handle as add */ + rs ^= 1; + + assert(zm & SP_HIDDEN_BIT); + + /* + * Provide guard,round and stick bit space. + */ + zm <<= 3; + + if (ze > re) { + /* + * Have to shift y fraction right to align. + */ + s = ze - re; + SPXSRSYn(s); + } else if (re > ze) { + /* + * Have to shift x fraction right to align. + */ + s = re - ze; + SPXSRSYn(s); + } + assert(ze == re); + assert(ze <= SP_EMAX); + + if (zs == rs) { + /* + * Generate 28 bit result of adding two 27 bit numbers + * leaving result in zm, zs and ze. + */ + zm = zm + rm; + + if (zm >> (SP_FBITS + 1 + 3)) { /* carry out */ + SPXSRSX1(); /* shift preserving sticky */ + } + } else { + if (zm >= rm) { + zm = zm - rm; + } else { + zm = rm - zm; + zs = rs; + } + if (zm == 0) + return ieee754sp_zero(ieee754_csr.rm == FPU_CSR_RD); + + /* + * Normalize in extended single precision + */ + while ((zm >> (SP_MBITS + 3)) == 0) { + zm <<= 1; + ze--; + } + + } + return ieee754sp_format(zs, ze, zm); +} -- cgit 1.2.3-korg