diff options
Diffstat (limited to 'arch/mips/math-emu')
47 files changed, 6787 insertions, 0 deletions
diff --git a/arch/mips/math-emu/Makefile b/arch/mips/math-emu/Makefile new file mode 100644 index 0000000..121a848 --- /dev/null +++ b/arch/mips/math-emu/Makefile @@ -0,0 +1,11 @@ +# +# Makefile for the Linux/MIPS kernel FPU emulation. +# + +obj-y := cp1emu.o ieee754m.o ieee754d.o ieee754dp.o ieee754sp.o ieee754.o \ + ieee754xcpt.o dp_frexp.o dp_modf.o dp_div.o dp_mul.o dp_sub.o \ + dp_add.o dp_fsp.o dp_cmp.o dp_logb.o dp_scalb.o dp_simple.o \ + dp_tint.o dp_fint.o dp_tlong.o dp_flong.o sp_frexp.o sp_modf.o \ + sp_div.o sp_mul.o sp_sub.o sp_add.o sp_fdp.o sp_cmp.o sp_logb.o \ + sp_scalb.o sp_simple.o sp_tint.o sp_fint.o sp_tlong.o sp_flong.o \ + dp_sqrt.o sp_sqrt.o kernel_linkage.o dsemul.o diff --git a/arch/mips/math-emu/cp1emu.c b/arch/mips/math-emu/cp1emu.c new file mode 100644 index 0000000..20a552b --- /dev/null +++ b/arch/mips/math-emu/cp1emu.c @@ -0,0 +1,1322 @@ +/* + * cp1emu.c: a MIPS coprocessor 1 (fpu) instruction emulator + * + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * Kevin D. Kissell, kevink@mips.com and Carsten Langgaard, carstenl@mips.com + * Copyright (C) 2000 MIPS Technologies, Inc. + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * A complete emulator for MIPS coprocessor 1 instructions. This is + * required for #float(switch) or #float(trap), where it catches all + * COP1 instructions via the "CoProcessor Unusable" exception. + * + * More surprisingly it is also required for #float(ieee), to help out + * the hardware fpu at the boundaries of the IEEE-754 representation + * (denormalised values, infinities, underflow, etc). It is made + * quite nasty because emulation of some non-COP1 instructions is + * required, e.g. in branch delay slots. + * + * Note if you know that you won't have an fpu, then you'll get much + * better performance by compiling with -msoft-float! + */ +#include <linux/sched.h> + +#include <asm/inst.h> +#include <asm/bootinfo.h> +#include <asm/cpu.h> +#include <asm/cpu-features.h> +#include <asm/processor.h> +#include <asm/ptrace.h> +#include <asm/signal.h> +#include <asm/mipsregs.h> +#include <asm/fpu_emulator.h> +#include <asm/uaccess.h> +#include <asm/branch.h> + +#include "ieee754.h" +#include "dsemul.h" + +/* Strap kernel emulator for full MIPS IV emulation */ + +#ifdef __mips +#undef __mips +#endif +#define __mips 4 + +/* Function which emulates a floating point instruction. */ + +static int fpu_emu(struct pt_regs *, struct mips_fpu_soft_struct *, + mips_instruction); + +#if __mips >= 4 && __mips != 32 +static int fpux_emu(struct pt_regs *, + struct mips_fpu_soft_struct *, mips_instruction); +#endif + +/* Further private data for which no space exists in mips_fpu_soft_struct */ + +struct mips_fpu_emulator_private fpuemuprivate; + +/* Control registers */ + +#define FPCREG_RID 0 /* $0 = revision id */ +#define FPCREG_CSR 31 /* $31 = csr */ + +/* Convert Mips rounding mode (0..3) to IEEE library modes. */ +static const unsigned char ieee_rm[4] = { + IEEE754_RN, IEEE754_RZ, IEEE754_RU, IEEE754_RD +}; + +#if __mips >= 4 +/* convert condition code register number to csr bit */ +static const unsigned int fpucondbit[8] = { + FPU_CSR_COND0, + FPU_CSR_COND1, + FPU_CSR_COND2, + FPU_CSR_COND3, + FPU_CSR_COND4, + FPU_CSR_COND5, + FPU_CSR_COND6, + FPU_CSR_COND7 +}; +#endif + + +/* + * Redundant with logic already in kernel/branch.c, + * embedded in compute_return_epc. At some point, + * a single subroutine should be used across both + * modules. + */ +static int isBranchInstr(mips_instruction * i) +{ + switch (MIPSInst_OPCODE(*i)) { + case spec_op: + switch (MIPSInst_FUNC(*i)) { + case jalr_op: + case jr_op: + return 1; + } + break; + + case bcond_op: + switch (MIPSInst_RT(*i)) { + case bltz_op: + case bgez_op: + case bltzl_op: + case bgezl_op: + case bltzal_op: + case bgezal_op: + case bltzall_op: + case bgezall_op: + return 1; + } + break; + + case j_op: + case jal_op: + case jalx_op: + case beq_op: + case bne_op: + case blez_op: + case bgtz_op: + case beql_op: + case bnel_op: + case blezl_op: + case bgtzl_op: + return 1; + + case cop0_op: + case cop1_op: + case cop2_op: + case cop1x_op: + if (MIPSInst_RS(*i) == bc_op) + return 1; + break; + } + + return 0; +} + +/* + * In the Linux kernel, we support selection of FPR format on the + * basis of the Status.FR bit. This does imply that, if a full 32 + * FPRs are desired, there needs to be a flip-flop that can be written + * to one at that bit position. In any case, O32 MIPS ABI uses + * only the even FPRs (Status.FR = 0). + */ + +#define CP0_STATUS_FR_SUPPORT + +#ifdef CP0_STATUS_FR_SUPPORT +#define FR_BIT ST0_FR +#else +#define FR_BIT 0 +#endif + +#define SIFROMREG(si,x) ((si) = \ + (xcp->cp0_status & FR_BIT) || !(x & 1) ? \ + (int)ctx->fpr[x] : \ + (int)(ctx->fpr[x & ~1] >> 32 )) +#define SITOREG(si,x) (ctx->fpr[x & ~((xcp->cp0_status & FR_BIT) == 0)] = \ + (xcp->cp0_status & FR_BIT) || !(x & 1) ? \ + ctx->fpr[x & ~1] >> 32 << 32 | (u32)(si) : \ + ctx->fpr[x & ~1] << 32 >> 32 | (u64)(si) << 32) + +#define DIFROMREG(di,x) ((di) = \ + ctx->fpr[x & ~((xcp->cp0_status & FR_BIT) == 0)]) +#define DITOREG(di,x) (ctx->fpr[x & ~((xcp->cp0_status & FR_BIT) == 0)] \ + = (di)) + +#define SPFROMREG(sp,x) SIFROMREG((sp).bits,x) +#define SPTOREG(sp,x) SITOREG((sp).bits,x) +#define DPFROMREG(dp,x) DIFROMREG((dp).bits,x) +#define DPTOREG(dp,x) DITOREG((dp).bits,x) + +/* + * Emulate the single floating point instruction pointed at by EPC. + * Two instructions if the instruction is in a branch delay slot. + */ + +static int cop1Emulate(struct pt_regs *xcp, struct mips_fpu_soft_struct *ctx) +{ + mips_instruction ir; + vaddr_t emulpc, contpc; + unsigned int cond; + + if (get_user(ir, (mips_instruction *) xcp->cp0_epc)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + + /* XXX NEC Vr54xx bug workaround */ + if ((xcp->cp0_cause & CAUSEF_BD) && !isBranchInstr(&ir)) + xcp->cp0_cause &= ~CAUSEF_BD; + + if (xcp->cp0_cause & CAUSEF_BD) { + /* + * The instruction to be emulated is in a branch delay slot + * which means that we have to emulate the branch instruction + * BEFORE we do the cop1 instruction. + * + * This branch could be a COP1 branch, but in that case we + * would have had a trap for that instruction, and would not + * come through this route. + * + * Linux MIPS branch emulator operates on context, updating the + * cp0_epc. + */ + emulpc = REG_TO_VA(xcp->cp0_epc + 4); /* Snapshot emulation target */ + + if (__compute_return_epc(xcp)) { +#ifdef CP1DBG + printk("failed to emulate branch at %p\n", + REG_TO_VA(xcp->cp0_epc)); +#endif + return SIGILL; + } + if (get_user(ir, (mips_instruction *) emulpc)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + /* __compute_return_epc() will have updated cp0_epc */ + contpc = REG_TO_VA xcp->cp0_epc; + /* In order not to confuse ptrace() et al, tweak context */ + xcp->cp0_epc = VA_TO_REG emulpc - 4; + } + else { + emulpc = REG_TO_VA xcp->cp0_epc; + contpc = REG_TO_VA(xcp->cp0_epc + 4); + } + + emul: + fpuemuprivate.stats.emulated++; + switch (MIPSInst_OPCODE(ir)) { +#ifndef SINGLE_ONLY_FPU + case ldc1_op:{ + u64 *va = REG_TO_VA(xcp->regs[MIPSInst_RS(ir)] + + MIPSInst_SIMM(ir)); + u64 val; + + fpuemuprivate.stats.loads++; + if (get_user(val, va)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + DITOREG(val, MIPSInst_RT(ir)); + break; + } + + case sdc1_op:{ + u64 *va = REG_TO_VA(xcp->regs[MIPSInst_RS(ir)] + + MIPSInst_SIMM(ir)); + u64 val; + + fpuemuprivate.stats.stores++; + DIFROMREG(val, MIPSInst_RT(ir)); + if (put_user(val, va)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + break; + } +#endif + + case lwc1_op:{ + u32 *va = REG_TO_VA(xcp->regs[MIPSInst_RS(ir)] + + MIPSInst_SIMM(ir)); + u32 val; + + fpuemuprivate.stats.loads++; + if (get_user(val, va)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } +#ifdef SINGLE_ONLY_FPU + if (MIPSInst_RT(ir) & 1) { + /* illegal register in single-float mode */ + return SIGILL; + } +#endif + SITOREG(val, MIPSInst_RT(ir)); + break; + } + + case swc1_op:{ + u32 *va = REG_TO_VA(xcp->regs[MIPSInst_RS(ir)] + + MIPSInst_SIMM(ir)); + u32 val; + + fpuemuprivate.stats.stores++; +#ifdef SINGLE_ONLY_FPU + if (MIPSInst_RT(ir) & 1) { + /* illegal register in single-float mode */ + return SIGILL; + } +#endif + SIFROMREG(val, MIPSInst_RT(ir)); + if (put_user(val, va)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + break; + } + + case cop1_op: + switch (MIPSInst_RS(ir)) { + +#if __mips64 && !defined(SINGLE_ONLY_FPU) + case dmfc_op: + /* copregister fs -> gpr[rt] */ + if (MIPSInst_RT(ir) != 0) { + DIFROMREG(xcp->regs[MIPSInst_RT(ir)], + MIPSInst_RD(ir)); + } + break; + + case dmtc_op: + /* copregister fs <- rt */ + DITOREG(xcp->regs[MIPSInst_RT(ir)], MIPSInst_RD(ir)); + break; +#endif + + case mfc_op: + /* copregister rd -> gpr[rt] */ +#ifdef SINGLE_ONLY_FPU + if (MIPSInst_RD(ir) & 1) { + /* illegal register in single-float mode */ + return SIGILL; + } +#endif + if (MIPSInst_RT(ir) != 0) { + SIFROMREG(xcp->regs[MIPSInst_RT(ir)], + MIPSInst_RD(ir)); + } + break; + + case mtc_op: + /* copregister rd <- rt */ +#ifdef SINGLE_ONLY_FPU + if (MIPSInst_RD(ir) & 1) { + /* illegal register in single-float mode */ + return SIGILL; + } +#endif + SITOREG(xcp->regs[MIPSInst_RT(ir)], MIPSInst_RD(ir)); + break; + + case cfc_op:{ + /* cop control register rd -> gpr[rt] */ + u32 value; + + if (ir == CP1UNDEF) { + return do_dsemulret(xcp); + } + if (MIPSInst_RD(ir) == FPCREG_CSR) { + value = ctx->fcr31; +#ifdef CSRTRACE + printk("%p gpr[%d]<-csr=%08x\n", + REG_TO_VA(xcp->cp0_epc), + MIPSInst_RT(ir), value); +#endif + } + else if (MIPSInst_RD(ir) == FPCREG_RID) + value = 0; + else + value = 0; + if (MIPSInst_RT(ir)) + xcp->regs[MIPSInst_RT(ir)] = value; + break; + } + + case ctc_op:{ + /* copregister rd <- rt */ + u32 value; + + if (MIPSInst_RT(ir) == 0) + value = 0; + else + value = xcp->regs[MIPSInst_RT(ir)]; + + /* we only have one writable control reg + */ + if (MIPSInst_RD(ir) == FPCREG_CSR) { +#ifdef CSRTRACE + printk("%p gpr[%d]->csr=%08x\n", + REG_TO_VA(xcp->cp0_epc), + MIPSInst_RT(ir), value); +#endif + ctx->fcr31 = value; + /* copy new rounding mode and + flush bit to ieee library state! */ + ieee754_csr.nod = (ctx->fcr31 & 0x1000000) != 0; + ieee754_csr.rm = ieee_rm[value & 0x3]; + } + if ((ctx->fcr31 >> 5) & ctx->fcr31 & FPU_CSR_ALL_E) { + return SIGFPE; + } + break; + } + + case bc_op:{ + int likely = 0; + + if (xcp->cp0_cause & CAUSEF_BD) + return SIGILL; + +#if __mips >= 4 + cond = ctx->fcr31 & fpucondbit[MIPSInst_RT(ir) >> 2]; +#else + cond = ctx->fcr31 & FPU_CSR_COND; +#endif + switch (MIPSInst_RT(ir) & 3) { + case bcfl_op: + likely = 1; + case bcf_op: + cond = !cond; + break; + case bctl_op: + likely = 1; + case bct_op: + break; + default: + /* thats an illegal instruction */ + return SIGILL; + } + + xcp->cp0_cause |= CAUSEF_BD; + if (cond) { + /* branch taken: emulate dslot + * instruction + */ + xcp->cp0_epc += 4; + contpc = REG_TO_VA + (xcp->cp0_epc + + (MIPSInst_SIMM(ir) << 2)); + + if (get_user(ir, (mips_instruction *) + REG_TO_VA xcp->cp0_epc)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + + switch (MIPSInst_OPCODE(ir)) { + case lwc1_op: + case swc1_op: +#if (__mips >= 2 || __mips64) && !defined(SINGLE_ONLY_FPU) + case ldc1_op: + case sdc1_op: +#endif + case cop1_op: +#if __mips >= 4 && __mips != 32 + case cop1x_op: +#endif + /* its one of ours */ + goto emul; +#if __mips >= 4 + case spec_op: + if (MIPSInst_FUNC(ir) == movc_op) + goto emul; + break; +#endif + } + + /* + * Single step the non-cp1 + * instruction in the dslot + */ + return mips_dsemul(xcp, ir, VA_TO_REG contpc); + } + else { + /* branch not taken */ + if (likely) { + /* + * branch likely nullifies + * dslot if not taken + */ + xcp->cp0_epc += 4; + contpc += 4; + /* + * else continue & execute + * dslot as normal insn + */ + } + } + break; + } + + default: + if (!(MIPSInst_RS(ir) & 0x10)) + return SIGILL; + { + int sig; + + /* a real fpu computation instruction */ + if ((sig = fpu_emu(xcp, ctx, ir))) + return sig; + } + } + break; + +#if __mips >= 4 && __mips != 32 + case cop1x_op:{ + int sig; + + if ((sig = fpux_emu(xcp, ctx, ir))) + return sig; + break; + } +#endif + +#if __mips >= 4 + case spec_op: + if (MIPSInst_FUNC(ir) != movc_op) + return SIGILL; + cond = fpucondbit[MIPSInst_RT(ir) >> 2]; + if (((ctx->fcr31 & cond) != 0) == ((MIPSInst_RT(ir) & 1) != 0)) + xcp->regs[MIPSInst_RD(ir)] = + xcp->regs[MIPSInst_RS(ir)]; + break; +#endif + + default: + return SIGILL; + } + + /* we did it !! */ + xcp->cp0_epc = VA_TO_REG(contpc); + xcp->cp0_cause &= ~CAUSEF_BD; + return 0; +} + +/* + * Conversion table from MIPS compare ops 48-63 + * cond = ieee754dp_cmp(x,y,IEEE754_UN,sig); + */ +static const unsigned char cmptab[8] = { + 0, /* cmp_0 (sig) cmp_sf */ + IEEE754_CUN, /* cmp_un (sig) cmp_ngle */ + IEEE754_CEQ, /* cmp_eq (sig) cmp_seq */ + IEEE754_CEQ | IEEE754_CUN, /* cmp_ueq (sig) cmp_ngl */ + IEEE754_CLT, /* cmp_olt (sig) cmp_lt */ + IEEE754_CLT | IEEE754_CUN, /* cmp_ult (sig) cmp_nge */ + IEEE754_CLT | IEEE754_CEQ, /* cmp_ole (sig) cmp_le */ + IEEE754_CLT | IEEE754_CEQ | IEEE754_CUN, /* cmp_ule (sig) cmp_ngt */ +}; + + +#if __mips >= 4 && __mips != 32 + +/* + * Additional MIPS4 instructions + */ + +#define DEF3OP(name, p, f1, f2, f3) \ +static ieee754##p fpemu_##p##_##name (ieee754##p r, ieee754##p s, \ + ieee754##p t) \ +{ \ + struct ieee754_csr ieee754_csr_save; \ + s = f1 (s, t); \ + ieee754_csr_save = ieee754_csr; \ + s = f2 (s, r); \ + ieee754_csr_save.cx |= ieee754_csr.cx; \ + ieee754_csr_save.sx |= ieee754_csr.sx; \ + s = f3 (s); \ + ieee754_csr.cx |= ieee754_csr_save.cx; \ + ieee754_csr.sx |= ieee754_csr_save.sx; \ + return s; \ +} + +static ieee754dp fpemu_dp_recip(ieee754dp d) +{ + return ieee754dp_div(ieee754dp_one(0), d); +} + +static ieee754dp fpemu_dp_rsqrt(ieee754dp d) +{ + return ieee754dp_div(ieee754dp_one(0), ieee754dp_sqrt(d)); +} + +static ieee754sp fpemu_sp_recip(ieee754sp s) +{ + return ieee754sp_div(ieee754sp_one(0), s); +} + +static ieee754sp fpemu_sp_rsqrt(ieee754sp s) +{ + return ieee754sp_div(ieee754sp_one(0), ieee754sp_sqrt(s)); +} + +DEF3OP(madd, sp, ieee754sp_mul, ieee754sp_add,); +DEF3OP(msub, sp, ieee754sp_mul, ieee754sp_sub,); +DEF3OP(nmadd, sp, ieee754sp_mul, ieee754sp_add, ieee754sp_neg); +DEF3OP(nmsub, sp, ieee754sp_mul, ieee754sp_sub, ieee754sp_neg); +DEF3OP(madd, dp, ieee754dp_mul, ieee754dp_add,); +DEF3OP(msub, dp, ieee754dp_mul, ieee754dp_sub,); +DEF3OP(nmadd, dp, ieee754dp_mul, ieee754dp_add, ieee754dp_neg); +DEF3OP(nmsub, dp, ieee754dp_mul, ieee754dp_sub, ieee754dp_neg); + +static int fpux_emu(struct pt_regs *xcp, struct mips_fpu_soft_struct *ctx, + mips_instruction ir) +{ + unsigned rcsr = 0; /* resulting csr */ + + fpuemuprivate.stats.cp1xops++; + + switch (MIPSInst_FMA_FFMT(ir)) { + case s_fmt:{ /* 0 */ + + ieee754sp(*handler) (ieee754sp, ieee754sp, ieee754sp); + ieee754sp fd, fr, fs, ft; + u32 *va; + u32 val; + + switch (MIPSInst_FUNC(ir)) { + case lwxc1_op: + va = REG_TO_VA(xcp->regs[MIPSInst_FR(ir)] + + xcp->regs[MIPSInst_FT(ir)]); + + fpuemuprivate.stats.loads++; + if (get_user(val, va)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } +#ifdef SINGLE_ONLY_FPU + if (MIPSInst_FD(ir) & 1) { + /* illegal register in single-float + * mode + */ + return SIGILL; + } +#endif + SITOREG(val, MIPSInst_FD(ir)); + break; + + case swxc1_op: + va = REG_TO_VA(xcp->regs[MIPSInst_FR(ir)] + + xcp->regs[MIPSInst_FT(ir)]); + + fpuemuprivate.stats.stores++; +#ifdef SINGLE_ONLY_FPU + if (MIPSInst_FS(ir) & 1) { + /* illegal register in single-float + * mode + */ + return SIGILL; + } +#endif + + SIFROMREG(val, MIPSInst_FS(ir)); + if (put_user(val, va)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + break; + + case madd_s_op: + handler = fpemu_sp_madd; + goto scoptop; + case msub_s_op: + handler = fpemu_sp_msub; + goto scoptop; + case nmadd_s_op: + handler = fpemu_sp_nmadd; + goto scoptop; + case nmsub_s_op: + handler = fpemu_sp_nmsub; + goto scoptop; + + scoptop: + SPFROMREG(fr, MIPSInst_FR(ir)); + SPFROMREG(fs, MIPSInst_FS(ir)); + SPFROMREG(ft, MIPSInst_FT(ir)); + fd = (*handler) (fr, fs, ft); + SPTOREG(fd, MIPSInst_FD(ir)); + + copcsr: + if (ieee754_cxtest(IEEE754_INEXACT)) + rcsr |= FPU_CSR_INE_X | FPU_CSR_INE_S; + if (ieee754_cxtest(IEEE754_UNDERFLOW)) + rcsr |= FPU_CSR_UDF_X | FPU_CSR_UDF_S; + if (ieee754_cxtest(IEEE754_OVERFLOW)) + rcsr |= FPU_CSR_OVF_X | FPU_CSR_OVF_S; + if (ieee754_cxtest(IEEE754_INVALID_OPERATION)) + rcsr |= FPU_CSR_INV_X | FPU_CSR_INV_S; + + ctx->fcr31 = (ctx->fcr31 & ~FPU_CSR_ALL_X) | rcsr; + if (ieee754_csr.nod) + ctx->fcr31 |= 0x1000000; + if ((ctx->fcr31 >> 5) & ctx->fcr31 & FPU_CSR_ALL_E) { + /*printk ("SIGFPE: fpu csr = %08x\n", + ctx->fcr31); */ + return SIGFPE; + } + + break; + + default: + return SIGILL; + } + break; + } + +#ifndef SINGLE_ONLY_FPU + case d_fmt:{ /* 1 */ + ieee754dp(*handler) (ieee754dp, ieee754dp, ieee754dp); + ieee754dp fd, fr, fs, ft; + u64 *va; + u64 val; + + switch (MIPSInst_FUNC(ir)) { + case ldxc1_op: + va = REG_TO_VA(xcp->regs[MIPSInst_FR(ir)] + + xcp->regs[MIPSInst_FT(ir)]); + + fpuemuprivate.stats.loads++; + if (get_user(val, va)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + DITOREG(val, MIPSInst_FD(ir)); + break; + + case sdxc1_op: + va = REG_TO_VA(xcp->regs[MIPSInst_FR(ir)] + + xcp->regs[MIPSInst_FT(ir)]); + + fpuemuprivate.stats.stores++; + DIFROMREG(val, MIPSInst_FS(ir)); + if (put_user(val, va)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + break; + + case madd_d_op: + handler = fpemu_dp_madd; + goto dcoptop; + case msub_d_op: + handler = fpemu_dp_msub; + goto dcoptop; + case nmadd_d_op: + handler = fpemu_dp_nmadd; + goto dcoptop; + case nmsub_d_op: + handler = fpemu_dp_nmsub; + goto dcoptop; + + dcoptop: + DPFROMREG(fr, MIPSInst_FR(ir)); + DPFROMREG(fs, MIPSInst_FS(ir)); + DPFROMREG(ft, MIPSInst_FT(ir)); + fd = (*handler) (fr, fs, ft); + DPTOREG(fd, MIPSInst_FD(ir)); + goto copcsr; + + default: + return SIGILL; + } + break; + } +#endif + + case 0x7: /* 7 */ + if (MIPSInst_FUNC(ir) != pfetch_op) { + return SIGILL; + } + /* ignore prefx operation */ + break; + + default: + return SIGILL; + } + + return 0; +} +#endif + + + +/* + * Emulate a single COP1 arithmetic instruction. + */ +static int fpu_emu(struct pt_regs *xcp, struct mips_fpu_soft_struct *ctx, + mips_instruction ir) +{ + int rfmt; /* resulting format */ + unsigned rcsr = 0; /* resulting csr */ + unsigned cond; + union { + ieee754dp d; + ieee754sp s; + int w; +#if __mips64 + s64 l; +#endif + } rv; /* resulting value */ + + fpuemuprivate.stats.cp1ops++; + switch (rfmt = (MIPSInst_FFMT(ir) & 0xf)) { + case s_fmt:{ /* 0 */ + union { + ieee754sp(*b) (ieee754sp, ieee754sp); + ieee754sp(*u) (ieee754sp); + } handler; + + switch (MIPSInst_FUNC(ir)) { + /* binary ops */ + case fadd_op: + handler.b = ieee754sp_add; + goto scopbop; + case fsub_op: + handler.b = ieee754sp_sub; + goto scopbop; + case fmul_op: + handler.b = ieee754sp_mul; + goto scopbop; + case fdiv_op: + handler.b = ieee754sp_div; + goto scopbop; + + /* unary ops */ +#if __mips >= 2 || __mips64 + case fsqrt_op: + handler.u = ieee754sp_sqrt; + goto scopuop; +#endif +#if __mips >= 4 && __mips != 32 + case frsqrt_op: + handler.u = fpemu_sp_rsqrt; + goto scopuop; + case frecip_op: + handler.u = fpemu_sp_recip; + goto scopuop; +#endif +#if __mips >= 4 + case fmovc_op: + cond = fpucondbit[MIPSInst_FT(ir) >> 2]; + if (((ctx->fcr31 & cond) != 0) != + ((MIPSInst_FT(ir) & 1) != 0)) + return 0; + SPFROMREG(rv.s, MIPSInst_FS(ir)); + break; + case fmovz_op: + if (xcp->regs[MIPSInst_FT(ir)] != 0) + return 0; + SPFROMREG(rv.s, MIPSInst_FS(ir)); + break; + case fmovn_op: + if (xcp->regs[MIPSInst_FT(ir)] == 0) + return 0; + SPFROMREG(rv.s, MIPSInst_FS(ir)); + break; +#endif + case fabs_op: + handler.u = ieee754sp_abs; + goto scopuop; + case fneg_op: + handler.u = ieee754sp_neg; + goto scopuop; + case fmov_op: + /* an easy one */ + SPFROMREG(rv.s, MIPSInst_FS(ir)); + goto copcsr; + + /* binary op on handler */ + scopbop: + { + ieee754sp fs, ft; + + SPFROMREG(fs, MIPSInst_FS(ir)); + SPFROMREG(ft, MIPSInst_FT(ir)); + + rv.s = (*handler.b) (fs, ft); + goto copcsr; + } + scopuop: + { + ieee754sp fs; + + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.s = (*handler.u) (fs); + goto copcsr; + } + copcsr: + if (ieee754_cxtest(IEEE754_INEXACT)) + rcsr |= FPU_CSR_INE_X | FPU_CSR_INE_S; + if (ieee754_cxtest(IEEE754_UNDERFLOW)) + rcsr |= FPU_CSR_UDF_X | FPU_CSR_UDF_S; + if (ieee754_cxtest(IEEE754_OVERFLOW)) + rcsr |= FPU_CSR_OVF_X | FPU_CSR_OVF_S; + if (ieee754_cxtest(IEEE754_ZERO_DIVIDE)) + rcsr |= FPU_CSR_DIV_X | FPU_CSR_DIV_S; + if (ieee754_cxtest(IEEE754_INVALID_OPERATION)) + rcsr |= FPU_CSR_INV_X | FPU_CSR_INV_S; + break; + + /* unary conv ops */ + case fcvts_op: + return SIGILL; /* not defined */ + case fcvtd_op:{ +#ifdef SINGLE_ONLY_FPU + return SIGILL; /* not defined */ +#else + ieee754sp fs; + + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.d = ieee754dp_fsp(fs); + rfmt = d_fmt; + goto copcsr; + } +#endif + case fcvtw_op:{ + ieee754sp fs; + + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.w = ieee754sp_tint(fs); + rfmt = w_fmt; + goto copcsr; + } + +#if __mips >= 2 || __mips64 + case fround_op: + case ftrunc_op: + case fceil_op: + case ffloor_op:{ + unsigned int oldrm = ieee754_csr.rm; + ieee754sp fs; + + SPFROMREG(fs, MIPSInst_FS(ir)); + ieee754_csr.rm = ieee_rm[MIPSInst_FUNC(ir) & 0x3]; + rv.w = ieee754sp_tint(fs); + ieee754_csr.rm = oldrm; + rfmt = w_fmt; + goto copcsr; + } +#endif /* __mips >= 2 */ + +#if __mips64 && !defined(SINGLE_ONLY_FPU) + case fcvtl_op:{ + ieee754sp fs; + + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.l = ieee754sp_tlong(fs); + rfmt = l_fmt; + goto copcsr; + } + + case froundl_op: + case ftruncl_op: + case fceill_op: + case ffloorl_op:{ + unsigned int oldrm = ieee754_csr.rm; + ieee754sp fs; + + SPFROMREG(fs, MIPSInst_FS(ir)); + ieee754_csr.rm = ieee_rm[MIPSInst_FUNC(ir) & 0x3]; + rv.l = ieee754sp_tlong(fs); + ieee754_csr.rm = oldrm; + rfmt = l_fmt; + goto copcsr; + } +#endif /* __mips64 && !fpu(single) */ + + default: + if (MIPSInst_FUNC(ir) >= fcmp_op) { + unsigned cmpop = MIPSInst_FUNC(ir) - fcmp_op; + ieee754sp fs, ft; + + SPFROMREG(fs, MIPSInst_FS(ir)); + SPFROMREG(ft, MIPSInst_FT(ir)); + rv.w = ieee754sp_cmp(fs, ft, + cmptab[cmpop & 0x7], cmpop & 0x8); + rfmt = -1; + if ((cmpop & 0x8) && ieee754_cxtest + (IEEE754_INVALID_OPERATION)) + rcsr = FPU_CSR_INV_X | FPU_CSR_INV_S; + else + goto copcsr; + + } + else { + return SIGILL; + } + break; + } + break; + } + +#ifndef SINGLE_ONLY_FPU + case d_fmt:{ + union { + ieee754dp(*b) (ieee754dp, ieee754dp); + ieee754dp(*u) (ieee754dp); + } handler; + + switch (MIPSInst_FUNC(ir)) { + /* binary ops */ + case fadd_op: + handler.b = ieee754dp_add; + goto dcopbop; + case fsub_op: + handler.b = ieee754dp_sub; + goto dcopbop; + case fmul_op: + handler.b = ieee754dp_mul; + goto dcopbop; + case fdiv_op: + handler.b = ieee754dp_div; + goto dcopbop; + + /* unary ops */ +#if __mips >= 2 || __mips64 + case fsqrt_op: + handler.u = ieee754dp_sqrt; + goto dcopuop; +#endif +#if __mips >= 4 && __mips != 32 + case frsqrt_op: + handler.u = fpemu_dp_rsqrt; + goto dcopuop; + case frecip_op: + handler.u = fpemu_dp_recip; + goto dcopuop; +#endif +#if __mips >= 4 + case fmovc_op: + cond = fpucondbit[MIPSInst_FT(ir) >> 2]; + if (((ctx->fcr31 & cond) != 0) != + ((MIPSInst_FT(ir) & 1) != 0)) + return 0; + DPFROMREG(rv.d, MIPSInst_FS(ir)); + break; + case fmovz_op: + if (xcp->regs[MIPSInst_FT(ir)] != 0) + return 0; + DPFROMREG(rv.d, MIPSInst_FS(ir)); + break; + case fmovn_op: + if (xcp->regs[MIPSInst_FT(ir)] == 0) + return 0; + DPFROMREG(rv.d, MIPSInst_FS(ir)); + break; +#endif + case fabs_op: + handler.u = ieee754dp_abs; + goto dcopuop; + + case fneg_op: + handler.u = ieee754dp_neg; + goto dcopuop; + + case fmov_op: + /* an easy one */ + DPFROMREG(rv.d, MIPSInst_FS(ir)); + goto copcsr; + + /* binary op on handler */ + dcopbop:{ + ieee754dp fs, ft; + + DPFROMREG(fs, MIPSInst_FS(ir)); + DPFROMREG(ft, MIPSInst_FT(ir)); + + rv.d = (*handler.b) (fs, ft); + goto copcsr; + } + dcopuop:{ + ieee754dp fs; + + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.d = (*handler.u) (fs); + goto copcsr; + } + + /* unary conv ops */ + case fcvts_op:{ + ieee754dp fs; + + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.s = ieee754sp_fdp(fs); + rfmt = s_fmt; + goto copcsr; + } + case fcvtd_op: + return SIGILL; /* not defined */ + + case fcvtw_op:{ + ieee754dp fs; + + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.w = ieee754dp_tint(fs); /* wrong */ + rfmt = w_fmt; + goto copcsr; + } + +#if __mips >= 2 || __mips64 + case fround_op: + case ftrunc_op: + case fceil_op: + case ffloor_op:{ + unsigned int oldrm = ieee754_csr.rm; + ieee754dp fs; + + DPFROMREG(fs, MIPSInst_FS(ir)); + ieee754_csr.rm = ieee_rm[MIPSInst_FUNC(ir) & 0x3]; + rv.w = ieee754dp_tint(fs); + ieee754_csr.rm = oldrm; + rfmt = w_fmt; + goto copcsr; + } +#endif + +#if __mips64 && !defined(SINGLE_ONLY_FPU) + case fcvtl_op:{ + ieee754dp fs; + + DPFROMREG(fs, MIPSInst_FS(ir)); + rv.l = ieee754dp_tlong(fs); + rfmt = l_fmt; + goto copcsr; + } + + case froundl_op: + case ftruncl_op: + case fceill_op: + case ffloorl_op:{ + unsigned int oldrm = ieee754_csr.rm; + ieee754dp fs; + + DPFROMREG(fs, MIPSInst_FS(ir)); + ieee754_csr.rm = ieee_rm[MIPSInst_FUNC(ir) & 0x3]; + rv.l = ieee754dp_tlong(fs); + ieee754_csr.rm = oldrm; + rfmt = l_fmt; + goto copcsr; + } +#endif /* __mips >= 3 && !fpu(single) */ + + default: + if (MIPSInst_FUNC(ir) >= fcmp_op) { + unsigned cmpop = MIPSInst_FUNC(ir) - fcmp_op; + ieee754dp fs, ft; + + DPFROMREG(fs, MIPSInst_FS(ir)); + DPFROMREG(ft, MIPSInst_FT(ir)); + rv.w = ieee754dp_cmp(fs, ft, + cmptab[cmpop & 0x7], cmpop & 0x8); + rfmt = -1; + if ((cmpop & 0x8) + && + ieee754_cxtest + (IEEE754_INVALID_OPERATION)) + rcsr = FPU_CSR_INV_X | FPU_CSR_INV_S; + else + goto copcsr; + + } + else { + return SIGILL; + } + break; + } + break; + } +#endif /* ifndef SINGLE_ONLY_FPU */ + + case w_fmt:{ + ieee754sp fs; + + switch (MIPSInst_FUNC(ir)) { + case fcvts_op: + /* convert word to single precision real */ + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.s = ieee754sp_fint(fs.bits); + rfmt = s_fmt; + goto copcsr; +#ifndef SINGLE_ONLY_FPU + case fcvtd_op: + /* convert word to double precision real */ + SPFROMREG(fs, MIPSInst_FS(ir)); + rv.d = ieee754dp_fint(fs.bits); + rfmt = d_fmt; + goto copcsr; +#endif + default: + return SIGILL; + } + break; + } + +#if __mips64 && !defined(SINGLE_ONLY_FPU) + case l_fmt:{ + switch (MIPSInst_FUNC(ir)) { + case fcvts_op: + /* convert long to single precision real */ + rv.s = ieee754sp_flong(ctx->fpr[MIPSInst_FS(ir)]); + rfmt = s_fmt; + goto copcsr; + case fcvtd_op: + /* convert long to double precision real */ + rv.d = ieee754dp_flong(ctx->fpr[MIPSInst_FS(ir)]); + rfmt = d_fmt; + goto copcsr; + default: + return SIGILL; + } + break; + } +#endif + + default: + return SIGILL; + } + + /* + * Update the fpu CSR register for this operation. + * If an exception is required, generate a tidy SIGFPE exception, + * without updating the result register. + * Note: cause exception bits do not accumulate, they are rewritten + * for each op; only the flag/sticky bits accumulate. + */ + ctx->fcr31 = (ctx->fcr31 & ~FPU_CSR_ALL_X) | rcsr; + if ((ctx->fcr31 >> 5) & ctx->fcr31 & FPU_CSR_ALL_E) { + /*printk ("SIGFPE: fpu csr = %08x\n",ctx->fcr31); */ + return SIGFPE; + } + + /* + * Now we can safely write the result back to the register file. + */ + switch (rfmt) { + case -1:{ +#if __mips >= 4 + cond = fpucondbit[MIPSInst_FD(ir) >> 2]; +#else + cond = FPU_CSR_COND; +#endif + if (rv.w) + ctx->fcr31 |= cond; + else + ctx->fcr31 &= ~cond; + break; + } +#ifndef SINGLE_ONLY_FPU + case d_fmt: + DPTOREG(rv.d, MIPSInst_FD(ir)); + break; +#endif + case s_fmt: + SPTOREG(rv.s, MIPSInst_FD(ir)); + break; + case w_fmt: + SITOREG(rv.w, MIPSInst_FD(ir)); + break; +#if __mips64 && !defined(SINGLE_ONLY_FPU) + case l_fmt: + DITOREG(rv.l, MIPSInst_FD(ir)); + break; +#endif + default: + return SIGILL; + } + + return 0; +} + +int fpu_emulator_cop1Handler(int xcptno, struct pt_regs *xcp, + struct mips_fpu_soft_struct *ctx) +{ + gpreg_t oldepc, prevepc; + mips_instruction insn; + int sig = 0; + + oldepc = xcp->cp0_epc; + do { + prevepc = xcp->cp0_epc; + + if (get_user(insn, (mips_instruction *) xcp->cp0_epc)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + if (insn == 0) + xcp->cp0_epc += 4; /* skip nops */ + else { + /* Update ieee754_csr. Only relevant if we have a + h/w FPU */ + ieee754_csr.nod = (ctx->fcr31 & 0x1000000) != 0; + ieee754_csr.rm = ieee_rm[ctx->fcr31 & 0x3]; + ieee754_csr.cx = (ctx->fcr31 >> 12) & 0x1f; + sig = cop1Emulate(xcp, ctx); + } + + if (cpu_has_fpu) + break; + if (sig) + break; + + cond_resched(); + } while (xcp->cp0_epc > prevepc); + + /* SIGILL indicates a non-fpu instruction */ + if (sig == SIGILL && xcp->cp0_epc != oldepc) + /* but if epc has advanced, then ignore it */ + sig = 0; + + return sig; +} diff --git a/arch/mips/math-emu/dp_add.c b/arch/mips/math-emu/dp_add.c new file mode 100644 index 0000000..bcf73bb --- /dev/null +++ b/arch/mips/math-emu/dp_add.c @@ -0,0 +1,183 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + * + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_add(ieee754dp x, ieee754dp y) +{ + COMPXDP; + COMPYDP; + + EXPLODEXDP; + EXPLODEYDP; + + CLEARCX; + + FLUSHXDP; + FLUSHYDP; + + switch (CLPAIR(xc, yc)) { + case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): + case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): + case CLPAIR(IEEE754_CLASS_SNAN, 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): + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "add", x, y); + + 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; + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_xcpt(ieee754dp_indef(), "add", x, y); + + 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 == + IEEE754_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; + + /* FALL THROUGH */ + + 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 + */ + int s = xe - ye; + ym = XDPSRS(ym, s); + ye += s; + } else if (ye > xe) { + /* have to shift x fraction right to align + */ + int 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,xe + */ + xm = xm + ym; + xe = xe; + xs = xs; + + if (xm >> (DP_MBITS + 1 + 3)) { /* carry out */ + xm = XDPSRS1(xm); + xe++; + } + } else { + if (xm >= ym) { + xm = xm - ym; + xe = xe; + xs = xs; + } else { + xm = ym - xm; + xe = xe; + xs = ys; + } + if (xm == 0) + return ieee754dp_zero(ieee754_csr.rm == + IEEE754_RD); + + /* normalize to rounding precision */ + while ((xm >> (DP_MBITS + 3)) == 0) { + xm <<= 1; + xe--; + } + + } + DPNORMRET2(xs, xe, xm, "add", x, y); +} diff --git a/arch/mips/math-emu/dp_cmp.c b/arch/mips/math-emu/dp_cmp.c new file mode 100644 index 0000000..8ab4f32 --- /dev/null +++ b/arch/mips/math-emu/dp_cmp.c @@ -0,0 +1,67 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +int ieee754dp_cmp(ieee754dp x, ieee754dp y, int cmp, int sig) +{ + COMPXDP; + COMPYDP; + + EXPLODEXDP; + EXPLODEYDP; + FLUSHXDP; + FLUSHYDP; + CLEARCX; /* Even clear inexact flag here */ + + if (ieee754dp_isnan(x) || ieee754dp_isnan(y)) { + if (sig || xc == IEEE754_CLASS_SNAN || yc == IEEE754_CLASS_SNAN) + SETCX(IEEE754_INVALID_OPERATION); + if (cmp & IEEE754_CUN) + return 1; + if (cmp & (IEEE754_CLT | IEEE754_CGT)) { + if (sig && SETANDTESTCX(IEEE754_INVALID_OPERATION)) + return ieee754si_xcpt(0, "fcmpf", x); + } + return 0; + } else { + s64 vx = x.bits; + s64 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; + } +} diff --git a/arch/mips/math-emu/dp_div.c b/arch/mips/math-emu/dp_div.c new file mode 100644 index 0000000..6acedce --- /dev/null +++ b/arch/mips/math-emu/dp_div.c @@ -0,0 +1,157 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_div(ieee754dp x, ieee754dp y) +{ + COMPXDP; + COMPYDP; + + EXPLODEXDP; + EXPLODEYDP; + + CLEARCX; + + FLUSHXDP; + FLUSHYDP; + + switch (CLPAIR(xc, yc)) { + case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): + case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): + case CLPAIR(IEEE754_CLASS_SNAN, 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): + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "div", x, y); + + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_xcpt(ieee754dp_indef(), "div", x, y); + + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_xcpt(ieee754dp_indef(), "div", x, y); + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO): + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO): + SETCX(IEEE754_ZERO_DIVIDE); + return ieee754dp_xcpt(ieee754dp_inf(xs ^ ys), "div", x, y); + + 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; + + 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 */ + + u64 rm = 0; + int re = xe - ye; + u64 bm; + + for (bm = DP_MBIT(DP_MBITS + 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_MBITS + 3)) == 0) { + rm <<= 1; + re--; + } + + DPNORMRET2(xs == ys ? 0 : 1, re, rm, "div", x, y); + } +} diff --git a/arch/mips/math-emu/dp_fint.c b/arch/mips/math-emu/dp_fint.c new file mode 100644 index 0000000..0065dea --- /dev/null +++ b/arch/mips/math-emu/dp_fint.c @@ -0,0 +1,80 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_fint(int x) +{ + COMPXDP; + + CLEARCX; + + xc = ( 0 ? xc : xc ); + + 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; + } + +#if 1 + /* normalize - result can never be inexact or overflow */ + xe = DP_MBITS; + while ((xm >> DP_MBITS) == 0) { + xm <<= 1; + xe--; + } + return builddp(xs, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); +#else + /* normalize */ + xe = DP_MBITS + 3; + while ((xm >> (DP_MBITS + 3)) == 0) { + xm <<= 1; + xe--; + } + DPNORMRET1(xs, xe, xm, "fint", x); +#endif +} + +ieee754dp ieee754dp_funs(unsigned int u) +{ + if ((int) u < 0) + return ieee754dp_add(ieee754dp_1e31(), + ieee754dp_fint(u & ~(1 << 31))); + return ieee754dp_fint(u); +} diff --git a/arch/mips/math-emu/dp_flong.c b/arch/mips/math-emu/dp_flong.c new file mode 100644 index 0000000..cb105b1 --- /dev/null +++ b/arch/mips/math-emu/dp_flong.c @@ -0,0 +1,78 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_flong(s64 x) +{ + COMPXDP; + + CLEARCX; + + xc = ( 0 ? xc : xc ); + + 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_MBITS + 3; + if (xm >> (DP_MBITS + 1 + 3)) { + /* shunt out overflow bits */ + while (xm >> (DP_MBITS + 1 + 3)) { + XDPSRSX1(); + } + } else { + /* normalize in grs extended double precision */ + while ((xm >> (DP_MBITS + 3)) == 0) { + xm <<= 1; + xe--; + } + } + DPNORMRET1(xs, xe, xm, "dp_flong", x); +} + +ieee754dp ieee754dp_fulong(u64 u) +{ + if ((s64) u < 0) + return ieee754dp_add(ieee754dp_1e63(), + ieee754dp_flong(u & ~(1ULL << 63))); + return ieee754dp_flong(u); +} diff --git a/arch/mips/math-emu/dp_frexp.c b/arch/mips/math-emu/dp_frexp.c new file mode 100644 index 0000000..e650cb1 --- /dev/null +++ b/arch/mips/math-emu/dp_frexp.c @@ -0,0 +1,53 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +/* close to ieeep754dp_logb +*/ +ieee754dp ieee754dp_frexp(ieee754dp x, int *eptr) +{ + COMPXDP; + CLEARCX; + EXPLODEXDP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + case IEEE754_CLASS_ZERO: + *eptr = 0; + return x; + case IEEE754_CLASS_DNORM: + DPDNORMX; + break; + case IEEE754_CLASS_NORM: + break; + } + *eptr = xe + 1; + return builddp(xs, -1 + DP_EBIAS, xm & ~DP_HIDDEN_BIT); +} diff --git a/arch/mips/math-emu/dp_fsp.c b/arch/mips/math-emu/dp_fsp.c new file mode 100644 index 0000000..494d19a --- /dev/null +++ b/arch/mips/math-emu/dp_fsp.c @@ -0,0 +1,74 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_fsp(ieee754sp x) +{ + COMPXSP; + + EXPLODEXSP; + + CLEARCX; + + FLUSHXSP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "fsp"); + case IEEE754_CLASS_QNAN: + return ieee754dp_nanxcpt(builddp(xs, + DP_EMAX + 1 + DP_EBIAS, + ((u64) xm + << (DP_MBITS - + SP_MBITS))), "fsp", + x); + case IEEE754_CLASS_INF: + return ieee754dp_inf(xs); + case IEEE754_CLASS_ZERO: + return ieee754dp_zero(xs); + case IEEE754_CLASS_DNORM: + /* normalize */ + while ((xm >> SP_MBITS) == 0) { + xm <<= 1; + xe--; + } + break; + case IEEE754_CLASS_NORM: + break; + } + + /* CANT possibly overflow,underflow, or need rounding + */ + + /* drop the hidden bit */ + xm &= ~SP_HIDDEN_BIT; + + return builddp(xs, xe + DP_EBIAS, + (u64) xm << (DP_MBITS - SP_MBITS)); +} diff --git a/arch/mips/math-emu/dp_logb.c b/arch/mips/math-emu/dp_logb.c new file mode 100644 index 0000000..6033886 --- /dev/null +++ b/arch/mips/math-emu/dp_logb.c @@ -0,0 +1,54 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_logb(ieee754dp x) +{ + COMPXDP; + + CLEARCX; + + EXPLODEXDP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + return ieee754dp_nanxcpt(x, "logb", x); + case IEEE754_CLASS_QNAN: + return x; + case IEEE754_CLASS_INF: + return ieee754dp_inf(0); + case IEEE754_CLASS_ZERO: + return ieee754dp_inf(1); + case IEEE754_CLASS_DNORM: + DPDNORMX; + break; + case IEEE754_CLASS_NORM: + break; + } + return ieee754dp_fint(xe); +} diff --git a/arch/mips/math-emu/dp_modf.c b/arch/mips/math-emu/dp_modf.c new file mode 100644 index 0000000..25861a4 --- /dev/null +++ b/arch/mips/math-emu/dp_modf.c @@ -0,0 +1,80 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +/* modf function is always exact for a finite number +*/ +ieee754dp ieee754dp_modf(ieee754dp x, ieee754dp * ip) +{ + COMPXDP; + + CLEARCX; + + EXPLODEXDP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + case IEEE754_CLASS_ZERO: + *ip = x; + return x; + case IEEE754_CLASS_DNORM: + /* far to small */ + *ip = ieee754dp_zero(xs); + return x; + case IEEE754_CLASS_NORM: + break; + } + if (xe < 0) { + *ip = ieee754dp_zero(xs); + return x; + } + if (xe >= DP_MBITS) { + *ip = x; + return ieee754dp_zero(xs); + } + /* generate ipart mantissa by clearing bottom bits + */ + *ip = builddp(xs, xe + DP_EBIAS, + ((xm >> (DP_MBITS - xe)) << (DP_MBITS - xe)) & + ~DP_HIDDEN_BIT); + + /* generate fpart mantissa by clearing top bits + * and normalizing (must be able to normalize) + */ + xm = (xm << (64 - (DP_MBITS - xe))) >> (64 - (DP_MBITS - xe)); + if (xm == 0) + return ieee754dp_zero(xs); + + while ((xm >> DP_MBITS) == 0) { + xm <<= 1; + xe--; + } + return builddp(xs, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); +} diff --git a/arch/mips/math-emu/dp_mul.c b/arch/mips/math-emu/dp_mul.c new file mode 100644 index 0000000..f237390 --- /dev/null +++ b/arch/mips/math-emu/dp_mul.c @@ -0,0 +1,177 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y) +{ + COMPXDP; + COMPYDP; + + EXPLODEXDP; + EXPLODEYDP; + + CLEARCX; + + FLUSHXDP; + FLUSHYDP; + + switch (CLPAIR(xc, yc)) { + case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): + case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): + case CLPAIR(IEEE754_CLASS_SNAN, 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): + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "mul", x, y); + + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_xcpt(ieee754dp_indef(), "mul", x, y); + + 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; + + 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 basicly */ + assert(xm & DP_HIDDEN_BIT); + assert(ym & DP_HIDDEN_BIT); + { + int re = xe + ye; + int rs = xs ^ ys; + u64 rm; + + /* shunt to top of word */ + xm <<= 64 - (DP_MBITS + 1); + ym <<= 64 - (DP_MBITS + 1); + + /* multiply 32bits xm,ym to give high 32bits rm with stickness + */ + + /* 32 * 32 => 64 */ +#define DPXMULT(x,y) ((u64)(x) * (u64)y) + + { + unsigned lxm = xm; + unsigned hxm = xm >> 32; + unsigned lym = ym; + unsigned hym = ym >> 32; + u64 lrm; + u64 hrm; + + lrm = DPXMULT(lxm, lym); + hrm = DPXMULT(hxm, hym); + + { + u64 t = DPXMULT(lxm, hym); + { + u64 at = + lrm + (t << 32); + hrm += at < lrm; + lrm = at; + } + hrm = hrm + (t >> 32); + } + + { + u64 t = DPXMULT(hxm, lym); + { + u64 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_MBITS + 1 + 3))) | + ((rm << (DP_MBITS + 1 + 3)) != 0); + re++; + } else { + rm = + (rm >> (64 - (DP_MBITS + 1 + 3 + 1))) | + ((rm << (DP_MBITS + 1 + 3 + 1)) != 0); + } + assert(rm & (DP_HIDDEN_BIT << 3)); + DPNORMRET2(rs, re, rm, "mul", x, y); + } +} diff --git a/arch/mips/math-emu/dp_scalb.c b/arch/mips/math-emu/dp_scalb.c new file mode 100644 index 0000000..b84e633 --- /dev/null +++ b/arch/mips/math-emu/dp_scalb.c @@ -0,0 +1,58 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_scalb(ieee754dp x, int n) +{ + COMPXDP; + + CLEARCX; + + EXPLODEXDP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + return ieee754dp_nanxcpt(x, "scalb", x, n); + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + case IEEE754_CLASS_ZERO: + return x; + case IEEE754_CLASS_DNORM: + DPDNORMX; + break; + case IEEE754_CLASS_NORM: + break; + } + DPNORMRET2(xs, xe + n, xm << 3, "scalb", x, n); +} + + +ieee754dp ieee754dp_ldexp(ieee754dp x, int n) +{ + return ieee754dp_scalb(x, n); +} diff --git a/arch/mips/math-emu/dp_simple.c b/arch/mips/math-emu/dp_simple.c new file mode 100644 index 0000000..495c1ac --- /dev/null +++ b/arch/mips/math-emu/dp_simple.c @@ -0,0 +1,84 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +int ieee754dp_finite(ieee754dp x) +{ + return DPBEXP(x) != DP_EMAX + 1 + DP_EBIAS; +} + +ieee754dp ieee754dp_copysign(ieee754dp x, ieee754dp y) +{ + CLEARCX; + DPSIGN(x) = DPSIGN(y); + return x; +} + + +ieee754dp ieee754dp_neg(ieee754dp x) +{ + COMPXDP; + + EXPLODEXDP; + CLEARCX; + FLUSHXDP; + + if (xc == IEEE754_CLASS_SNAN) { + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "neg"); + } + + if (ieee754dp_isnan(x)) /* but not infinity */ + return ieee754dp_nanxcpt(x, "neg", x); + + /* quick fix up */ + DPSIGN(x) ^= 1; + return x; +} + + +ieee754dp ieee754dp_abs(ieee754dp x) +{ + COMPXDP; + + EXPLODEXDP; + CLEARCX; + FLUSHXDP; + + if (xc == IEEE754_CLASS_SNAN) { + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "neg"); + } + + if (ieee754dp_isnan(x)) /* but not infinity */ + return ieee754dp_nanxcpt(x, "abs", x); + + /* quick fix up */ + DPSIGN(x) = 0; + return x; +} diff --git a/arch/mips/math-emu/dp_sqrt.c b/arch/mips/math-emu/dp_sqrt.c new file mode 100644 index 0000000..c35e871 --- /dev/null +++ b/arch/mips/math-emu/dp_sqrt.c @@ -0,0 +1,165 @@ +/* IEEE754 floating point arithmetic + * double precision square root + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +static const unsigned 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 +}; + +ieee754dp ieee754dp_sqrt(ieee754dp x) +{ + struct ieee754_csr oldcsr; + ieee754dp y, z, t; + unsigned scalx, yh; + COMPXDP; + + EXPLODEXDP; + CLEARCX; + FLUSHXDP; + + /* x == INF or NAN? */ + switch (xc) { + case IEEE754_CLASS_QNAN: + /* sqrt(Nan) = Nan */ + return ieee754dp_nanxcpt(x, "sqrt"); + case IEEE754_CLASS_SNAN: + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt"); + case IEEE754_CLASS_ZERO: + /* sqrt(0) = 0 */ + return x; + case IEEE754_CLASS_INF: + if (xs) { + /* sqrt(-Inf) = Nan */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt"); + } + /* sqrt(+Inf) = Inf */ + return x; + case IEEE754_CLASS_DNORM: + DPDNORMX; + /* fall through */ + case IEEE754_CLASS_NORM: + if (xs) { + /* sqrt(-x) = Nan */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt"); + } + 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 = IEEE754_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; + } + + y = x = builddp(0, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); + + /* 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; */ + z = t = ieee754dp_mul(y, y); + t.parts.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.parts.bexp += 0x001; + y = ieee754dp_add(y, t); + + /* twiddle last bit to force y correctly rounded */ + + /* set RZ, clear INEX flag */ + ieee754_csr.rm = IEEE754_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 IEEE754_RP: + y.bits += 1; + /* drop through */ + case IEEE754_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.parts.bexp += scalx; + + /* restore rounding mode, possibly set inexact */ + ieee754_csr = oldcsr; + + return y; +} diff --git a/arch/mips/math-emu/dp_sub.c b/arch/mips/math-emu/dp_sub.c new file mode 100644 index 0000000..b30c5b1 --- /dev/null +++ b/arch/mips/math-emu/dp_sub.c @@ -0,0 +1,191 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +ieee754dp ieee754dp_sub(ieee754dp x, ieee754dp y) +{ + COMPXDP; + COMPYDP; + + EXPLODEXDP; + EXPLODEYDP; + + CLEARCX; + + FLUSHXDP; + FLUSHYDP; + + switch (CLPAIR(xc, yc)) { + case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): + case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): + case CLPAIR(IEEE754_CLASS_SNAN, 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): + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_nanxcpt(ieee754dp_indef(), "sub", x, y); + + 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; + SETCX(IEEE754_INVALID_OPERATION); + return ieee754dp_xcpt(ieee754dp_indef(), "sub", x, y); + + 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 == + IEEE754_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; + /* FAAL THOROUGH */ + + 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 + */ + int s = xe - ye; + ym = XDPSRS(ym, s); + ye += s; + } else if (ye > xe) { + /* have to shift x fraction right to align + */ + int 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; + xe = xe; + xs = xs; + + if (xm >> (DP_MBITS + 1 + 3)) { /* carry out */ + xm = XDPSRS1(xm); /* shift preserving sticky */ + xe++; + } + } else { + if (xm >= ym) { + xm = xm - ym; + xe = xe; + xs = xs; + } else { + xm = ym - xm; + xe = xe; + xs = ys; + } + if (xm == 0) { + if (ieee754_csr.rm == IEEE754_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_MBITS + 3)) == 0) { + xm <<= 1; + xe--; + } + } + DPNORMRET2(xs, xe, xm, "sub", x, y); +} diff --git a/arch/mips/math-emu/dp_tint.c b/arch/mips/math-emu/dp_tint.c new file mode 100644 index 0000000..77b2b7c --- /dev/null +++ b/arch/mips/math-emu/dp_tint.c @@ -0,0 +1,124 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include <linux/kernel.h> +#include "ieee754dp.h" + +int ieee754dp_tint(ieee754dp x) +{ + COMPXDP; + + CLEARCX; + + EXPLODEXDP; + FLUSHXDP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + SETCX(IEEE754_INVALID_OPERATION); + return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x); + 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 */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x); + } + /* oh gawd */ + if (xe > DP_MBITS) { + xm <<= xe - DP_MBITS; + } else if (xe < DP_MBITS) { + u64 residue; + int round; + int sticky; + int odd; + + if (xe < -1) { + residue = xm; + round = 0; + sticky = residue != 0; + xm = 0; + } + else { + residue = xm << (64 - DP_MBITS + xe); + round = (residue >> 63) != 0; + sticky = (residue << 1) != 0; + xm >>= DP_MBITS - xe; + } + /* Note: At this point upper 32 bits of xm are guaranteed + to be zero */ + odd = (xm & 0x1) != 0x0; + switch (ieee754_csr.rm) { + case IEEE754_RN: + if (round && (sticky || odd)) + xm++; + break; + case IEEE754_RZ: + break; + case IEEE754_RU: /* toward +Infinity */ + if ((round || sticky) && !xs) + xm++; + break; + case IEEE754_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 */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x); + } + if (round || sticky) + SETCX(IEEE754_INEXACT); + } + if (xs) + return -xm; + else + return xm; +} + + +unsigned int ieee754dp_tuns(ieee754dp x) +{ + ieee754dp hb = ieee754dp_1e31(); + + /* what if x < 0 ?? */ + if (ieee754dp_lt(x, hb)) + return (unsigned) ieee754dp_tint(x); + + return (unsigned) ieee754dp_tint(ieee754dp_sub(x, hb)) | + ((unsigned) 1 << 31); +} diff --git a/arch/mips/math-emu/dp_tlong.c b/arch/mips/math-emu/dp_tlong.c new file mode 100644 index 0000000..d71113e --- /dev/null +++ b/arch/mips/math-emu/dp_tlong.c @@ -0,0 +1,127 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +s64 ieee754dp_tlong(ieee754dp x) +{ + COMPXDP; + + CLEARCX; + + EXPLODEXDP; + FLUSHXDP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + SETCX(IEEE754_INVALID_OPERATION); + return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x); + 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 */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x); + } + /* oh gawd */ + if (xe > DP_MBITS) { + xm <<= xe - DP_MBITS; + } else if (xe < DP_MBITS) { + u64 residue; + int round; + int sticky; + int odd; + + 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_MBITS; + round = (residue >> 63) != 0; + sticky = (residue << 1) != 0; + xm >>= DP_MBITS - xe; + } + odd = (xm & 0x1) != 0x0; + switch (ieee754_csr.rm) { + case IEEE754_RN: + if (round && (sticky || odd)) + xm++; + break; + case IEEE754_RZ: + break; + case IEEE754_RU: /* toward +Infinity */ + if ((round || sticky) && !xs) + xm++; + break; + case IEEE754_RD: /* toward -Infinity */ + if ((round || sticky) && xs) + xm++; + break; + } + if ((xm >> 63) != 0) { + /* This can happen after rounding */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x); + } + if (round || sticky) + SETCX(IEEE754_INEXACT); + } + if (xs) + return -xm; + else + return xm; +} + + +u64 ieee754dp_tulong(ieee754dp x) +{ + ieee754dp hb = ieee754dp_1e63(); + + /* what if x < 0 ?? */ + if (ieee754dp_lt(x, hb)) + return (u64) ieee754dp_tlong(x); + + return (u64) ieee754dp_tlong(ieee754dp_sub(x, hb)) | + (1ULL << 63); +} diff --git a/arch/mips/math-emu/dsemul.c b/arch/mips/math-emu/dsemul.c new file mode 100644 index 0000000..aa989c2 --- /dev/null +++ b/arch/mips/math-emu/dsemul.c @@ -0,0 +1,172 @@ +#include <linux/compiler.h> +#include <linux/mm.h> +#include <linux/signal.h> +#include <linux/smp.h> +#include <linux/smp_lock.h> + +#include <asm/asm.h> +#include <asm/bootinfo.h> +#include <asm/byteorder.h> +#include <asm/cpu.h> +#include <asm/inst.h> +#include <asm/processor.h> +#include <asm/uaccess.h> +#include <asm/branch.h> +#include <asm/mipsregs.h> +#include <asm/system.h> +#include <asm/cacheflush.h> + +#include <asm/fpu_emulator.h> + +#include "ieee754.h" +#include "dsemul.h" + +/* Strap kernel emulator for full MIPS IV emulation */ + +#ifdef __mips +#undef __mips +#endif +#define __mips 4 + +extern struct mips_fpu_emulator_private fpuemuprivate; + + +/* + * Emulate the arbritrary instruction ir at xcp->cp0_epc. Required when + * we have to emulate the instruction in a COP1 branch delay slot. Do + * not change cp0_epc due to the instruction + * + * According to the spec: + * 1) it shouldnt be a branch :-) + * 2) it can be a COP instruction :-( + * 3) if we are tring to run a protected memory space we must take + * special care on memory access instructions :-( + */ + +/* + * "Trampoline" return routine to catch exception following + * execution of delay-slot instruction execution. + */ + +struct emuframe { + mips_instruction emul; + mips_instruction badinst; + mips_instruction cookie; + gpreg_t epc; +}; + +int mips_dsemul(struct pt_regs *regs, mips_instruction ir, gpreg_t cpc) +{ + extern asmlinkage void handle_dsemulret(void); + mips_instruction *dsemul_insns; + struct emuframe *fr; + int err; + + if (ir == 0) { /* a nop is easy */ + regs->cp0_epc = cpc; + regs->cp0_cause &= ~CAUSEF_BD; + return 0; + } +#ifdef DSEMUL_TRACE + printk("dsemul %lx %lx\n", regs->cp0_epc, cpc); + +#endif + + /* + * The strategy is to push the instruction onto the user stack + * and put a trap after it which we can catch and jump to + * the required address any alternative apart from full + * instruction emulation!!. + * + * Algorithmics used a system call instruction, and + * borrowed that vector. MIPS/Linux version is a bit + * more heavyweight in the interests of portability and + * multiprocessor support. For Linux we generate a + * an unaligned access and force an address error exception. + * + * For embedded systems (stand-alone) we prefer to use a + * non-existing CP1 instruction. This prevents us from emulating + * branches, but gives us a cleaner interface to the exception + * handler (single entry point). + */ + + /* Ensure that the two instructions are in the same cache line */ + dsemul_insns = (mips_instruction *) REG_TO_VA ((regs->regs[29] - sizeof(struct emuframe)) & ~0x7); + fr = (struct emuframe *) dsemul_insns; + + /* Verify that the stack pointer is not competely insane */ + if (unlikely(!access_ok(VERIFY_WRITE, fr, sizeof(struct emuframe)))) + return SIGBUS; + + err = __put_user(ir, &fr->emul); + err |= __put_user((mips_instruction)BADINST, &fr->badinst); + err |= __put_user((mips_instruction)BD_COOKIE, &fr->cookie); + err |= __put_user(cpc, &fr->epc); + + if (unlikely(err)) { + fpuemuprivate.stats.errors++; + return SIGBUS; + } + + regs->cp0_epc = VA_TO_REG & fr->emul; + + flush_cache_sigtramp((unsigned long)&fr->badinst); + + return SIGILL; /* force out of emulation loop */ +} + +int do_dsemulret(struct pt_regs *xcp) +{ + struct emuframe *fr; + gpreg_t epc; + u32 insn, cookie; + int err = 0; + + fr = (struct emuframe *) (xcp->cp0_epc - sizeof(mips_instruction)); + + /* + * If we can't even access the area, something is very wrong, but we'll + * leave that to the default handling + */ + if (!access_ok(VERIFY_READ, fr, sizeof(struct emuframe))) + return 0; + + /* + * Do some sanity checking on the stackframe: + * + * - Is the instruction pointed to by the EPC an BADINST? + * - Is the following memory word the BD_COOKIE? + */ + err = __get_user(insn, &fr->badinst); + err |= __get_user(cookie, &fr->cookie); + + if (unlikely(err || (insn != BADINST) || (cookie != BD_COOKIE))) { + fpuemuprivate.stats.errors++; + return 0; + } + + /* + * At this point, we are satisfied that it's a BD emulation trap. Yes, + * a user might have deliberately put two malformed and useless + * instructions in a row in his program, in which case he's in for a + * nasty surprise - the next instruction will be treated as a + * continuation address! Alas, this seems to be the only way that we + * can handle signals, recursion, and longjmps() in the context of + * emulating the branch delay instruction. + */ + +#ifdef DSEMUL_TRACE + printk("dsemulret\n"); +#endif + if (__get_user(epc, &fr->epc)) { /* Saved EPC */ + /* This is not a good situation to be in */ + force_sig(SIGBUS, current); + + return 0; + } + + /* Set EPC to return to post-branch instruction */ + xcp->cp0_epc = epc; + + return 1; +} diff --git a/arch/mips/math-emu/dsemul.h b/arch/mips/math-emu/dsemul.h new file mode 100644 index 0000000..dbd85f9 --- /dev/null +++ b/arch/mips/math-emu/dsemul.h @@ -0,0 +1,23 @@ +typedef long gpreg_t; +typedef void *vaddr_t; + +#define REG_TO_VA (vaddr_t) +#define VA_TO_REG (gpreg_t) + +int mips_dsemul(struct pt_regs *regs, mips_instruction ir, gpreg_t cpc); +int do_dsemulret(struct pt_regs *xcp); + +/* Instruction which will always cause an address error */ +#define AdELOAD 0x8c000001 /* lw $0,1($0) */ +/* Instruction which will plainly cause a CP1 exception when FPU is disabled */ +#define CP1UNDEF 0x44400001 /* cfc1 $0,$0 undef */ + +/* Instruction inserted following the badinst to further tag the sequence */ +#define BD_COOKIE 0x0000bd36 /* tne $0,$0 with baggage */ + +/* Setup which instruction to use for trampoline */ +#ifdef STANDALONE_EMULATOR +#define BADINST CP1UNDEF +#else +#define BADINST AdELOAD +#endif diff --git a/arch/mips/math-emu/ieee754.c b/arch/mips/math-emu/ieee754.c new file mode 100644 index 0000000..f0a364a --- /dev/null +++ b/arch/mips/math-emu/ieee754.c @@ -0,0 +1,138 @@ +/* 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. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754int.h" + +#define DP_EBIAS 1023 +#define DP_EMIN (-1022) +#define DP_EMAX 1023 + +#define SP_EBIAS 127 +#define SP_EMIN (-126) +#define SP_EMAX 127 + +/* indexed by class */ +const char *const ieee754_cname[] = { + "Normal", + "Zero", + "Denormal", + "Infinity", + "QNaN", + "SNaN", +}; + +/* the control status register +*/ +struct ieee754_csr ieee754_csr; + +/* special constants +*/ + + +#if (defined(BYTE_ORDER) && BYTE_ORDER == LITTLE_ENDIAN) || defined(__MIPSEL__) +#define SPSTR(s,b,m) {m,b,s} +#define DPSTR(s,b,mh,ml) {ml,mh,b,s} +#endif + +#ifdef __MIPSEB__ +#define SPSTR(s,b,m) {s,b,m} +#define DPSTR(s,b,mh,ml) {s,b,mh,ml} +#endif + +const struct ieee754dp_konst __ieee754dp_spcvals[] = { + DPSTR(0, DP_EMIN - 1 + DP_EBIAS, 0, 0), /* + zero */ + DPSTR(1, DP_EMIN - 1 + DP_EBIAS, 0, 0), /* - zero */ + DPSTR(0, DP_EBIAS, 0, 0), /* + 1.0 */ + DPSTR(1, DP_EBIAS, 0, 0), /* - 1.0 */ + DPSTR(0, 3 + DP_EBIAS, 0x40000, 0), /* + 10.0 */ + DPSTR(1, 3 + DP_EBIAS, 0x40000, 0), /* - 10.0 */ + DPSTR(0, DP_EMAX + 1 + DP_EBIAS, 0, 0), /* + infinity */ + DPSTR(1, DP_EMAX + 1 + DP_EBIAS, 0, 0), /* - infinity */ + DPSTR(0,DP_EMAX+1+DP_EBIAS,0x7FFFF,0xFFFFFFFF), /* + indef quiet Nan */ + DPSTR(0, DP_EMAX + DP_EBIAS, 0xFFFFF, 0xFFFFFFFF), /* + max */ + DPSTR(1, DP_EMAX + DP_EBIAS, 0xFFFFF, 0xFFFFFFFF), /* - max */ + DPSTR(0, DP_EMIN + DP_EBIAS, 0, 0), /* + min normal */ + DPSTR(1, DP_EMIN + DP_EBIAS, 0, 0), /* - min normal */ + DPSTR(0, DP_EMIN - 1 + DP_EBIAS, 0, 1), /* + min denormal */ + DPSTR(1, DP_EMIN - 1 + DP_EBIAS, 0, 1), /* - min denormal */ + DPSTR(0, 31 + DP_EBIAS, 0, 0), /* + 1.0e31 */ + DPSTR(0, 63 + DP_EBIAS, 0, 0), /* + 1.0e63 */ +}; + +const struct ieee754sp_konst __ieee754sp_spcvals[] = { + SPSTR(0, SP_EMIN - 1 + SP_EBIAS, 0), /* + zero */ + SPSTR(1, SP_EMIN - 1 + SP_EBIAS, 0), /* - zero */ + SPSTR(0, SP_EBIAS, 0), /* + 1.0 */ + SPSTR(1, SP_EBIAS, 0), /* - 1.0 */ + SPSTR(0, 3 + SP_EBIAS, 0x200000), /* + 10.0 */ + SPSTR(1, 3 + SP_EBIAS, 0x200000), /* - 10.0 */ + SPSTR(0, SP_EMAX + 1 + SP_EBIAS, 0), /* + infinity */ + SPSTR(1, SP_EMAX + 1 + SP_EBIAS, 0), /* - infinity */ + SPSTR(0,SP_EMAX+1+SP_EBIAS,0x3FFFFF), /* + indef quiet Nan */ + SPSTR(0, SP_EMAX + SP_EBIAS, 0x7FFFFF), /* + max normal */ + SPSTR(1, SP_EMAX + SP_EBIAS, 0x7FFFFF), /* - max normal */ + SPSTR(0, SP_EMIN + SP_EBIAS, 0), /* + min normal */ + SPSTR(1, SP_EMIN + SP_EBIAS, 0), /* - min normal */ + SPSTR(0, SP_EMIN - 1 + SP_EBIAS, 1), /* + min denormal */ + SPSTR(1, SP_EMIN - 1 + SP_EBIAS, 1), /* - min denormal */ + SPSTR(0, 31 + SP_EBIAS, 0), /* + 1.0e31 */ + SPSTR(0, 63 + SP_EBIAS, 0), /* + 1.0e63 */ +}; + + +int ieee754si_xcpt(int r, const char *op, ...) +{ + struct ieee754xctx ax; + + if (!TSTX()) + return r; + ax.op = op; + ax.rt = IEEE754_RT_SI; + ax.rv.si = r; + va_start(ax.ap, op); + ieee754_xcpt(&ax); + return ax.rv.si; +} + +s64 ieee754di_xcpt(s64 r, const char *op, ...) +{ + struct ieee754xctx ax; + + if (!TSTX()) + return r; + ax.op = op; + ax.rt = IEEE754_RT_DI; + ax.rv.di = r; + va_start(ax.ap, op); + ieee754_xcpt(&ax); + return ax.rv.di; +} diff --git a/arch/mips/math-emu/ieee754.h b/arch/mips/math-emu/ieee754.h new file mode 100644 index 0000000..b8772f4 --- /dev/null +++ b/arch/mips/math-emu/ieee754.h @@ -0,0 +1,489 @@ +/* single and double precision fp ops + * missing extended precision. +*/ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + +/************************************************************************** + * 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. + *************************************************************************/ + +#ifdef __KERNEL__ +/* Going from Algorithmics to Linux native environment, add this */ +#include <linux/types.h> + +/* + * Not very pretty, but the Linux kernel's normal va_list definition + * does not allow it to be used as a structure element, as it is here. + */ +#ifndef _STDARG_H +#include <stdarg.h> +#endif + +#else + +/* Note that __KERNEL__ is taken to mean Linux kernel */ + +#if #system(OpenBSD) +#include <machine/types.h> +#endif +#include <machine/endian.h> + +#endif /* __KERNEL__ */ + +#if (defined(BYTE_ORDER) && BYTE_ORDER == LITTLE_ENDIAN) || defined(__MIPSEL__) +struct ieee754dp_konst { + unsigned mantlo:32; + unsigned manthi:20; + unsigned bexp:11; + unsigned sign:1; +}; +struct ieee754sp_konst { + unsigned mant:23; + unsigned bexp:8; + unsigned sign:1; +}; + +typedef union _ieee754dp { + struct ieee754dp_konst oparts; + struct { + u64 mant:52; + unsigned int bexp:11; + unsigned int sign:1; + } parts; + u64 bits; + double d; +} ieee754dp; + +typedef union _ieee754sp { + struct ieee754sp_konst parts; + float f; + u32 bits; +} ieee754sp; +#endif + +#if (defined(BYTE_ORDER) && BYTE_ORDER == BIG_ENDIAN) || defined(__MIPSEB__) +struct ieee754dp_konst { + unsigned sign:1; + unsigned bexp:11; + unsigned manthi:20; + unsigned mantlo:32; +}; +typedef union _ieee754dp { + struct ieee754dp_konst oparts; + struct { + unsigned int sign:1; + unsigned int bexp:11; + u64 mant:52; + } parts; + double d; + u64 bits; +} ieee754dp; + +struct ieee754sp_konst { + unsigned sign:1; + unsigned bexp:8; + unsigned mant:23; +}; + +typedef union _ieee754sp { + struct ieee754sp_konst parts; + float f; + u32 bits; +} ieee754sp; +#endif + +/* + * single precision (often aka float) +*/ +int ieee754sp_finite(ieee754sp x); +int ieee754sp_class(ieee754sp x); + +ieee754sp ieee754sp_abs(ieee754sp x); +ieee754sp ieee754sp_neg(ieee754sp x); +ieee754sp ieee754sp_scalb(ieee754sp x, int); +ieee754sp ieee754sp_logb(ieee754sp x); + +/* x with sign of y */ +ieee754sp ieee754sp_copysign(ieee754sp x, ieee754sp y); + +ieee754sp ieee754sp_add(ieee754sp x, ieee754sp y); +ieee754sp ieee754sp_sub(ieee754sp x, ieee754sp y); +ieee754sp ieee754sp_mul(ieee754sp x, ieee754sp y); +ieee754sp ieee754sp_div(ieee754sp x, ieee754sp y); + +ieee754sp ieee754sp_fint(int x); +ieee754sp ieee754sp_funs(unsigned x); +ieee754sp ieee754sp_flong(s64 x); +ieee754sp ieee754sp_fulong(u64 x); +ieee754sp ieee754sp_fdp(ieee754dp x); + +int ieee754sp_tint(ieee754sp x); +unsigned int ieee754sp_tuns(ieee754sp x); +s64 ieee754sp_tlong(ieee754sp x); +u64 ieee754sp_tulong(ieee754sp x); + +int ieee754sp_cmp(ieee754sp x, ieee754sp y, int cop, int sig); +/* + * basic sp math + */ +ieee754sp ieee754sp_modf(ieee754sp x, ieee754sp * ip); +ieee754sp ieee754sp_frexp(ieee754sp x, int *exp); +ieee754sp ieee754sp_ldexp(ieee754sp x, int exp); + +ieee754sp ieee754sp_ceil(ieee754sp x); +ieee754sp ieee754sp_floor(ieee754sp x); +ieee754sp ieee754sp_trunc(ieee754sp x); + +ieee754sp ieee754sp_sqrt(ieee754sp x); + +/* + * double precision (often aka double) +*/ +int ieee754dp_finite(ieee754dp x); +int ieee754dp_class(ieee754dp x); + +/* x with sign of y */ +ieee754dp ieee754dp_copysign(ieee754dp x, ieee754dp y); + +ieee754dp ieee754dp_add(ieee754dp x, ieee754dp y); +ieee754dp ieee754dp_sub(ieee754dp x, ieee754dp y); +ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y); +ieee754dp ieee754dp_div(ieee754dp x, ieee754dp y); + +ieee754dp ieee754dp_abs(ieee754dp x); +ieee754dp ieee754dp_neg(ieee754dp x); +ieee754dp ieee754dp_scalb(ieee754dp x, int); + +/* return exponent as integer in floating point format + */ +ieee754dp ieee754dp_logb(ieee754dp x); + +ieee754dp ieee754dp_fint(int x); +ieee754dp ieee754dp_funs(unsigned x); +ieee754dp ieee754dp_flong(s64 x); +ieee754dp ieee754dp_fulong(u64 x); +ieee754dp ieee754dp_fsp(ieee754sp x); + +ieee754dp ieee754dp_ceil(ieee754dp x); +ieee754dp ieee754dp_floor(ieee754dp x); +ieee754dp ieee754dp_trunc(ieee754dp x); + +int ieee754dp_tint(ieee754dp x); +unsigned int ieee754dp_tuns(ieee754dp x); +s64 ieee754dp_tlong(ieee754dp x); +u64 ieee754dp_tulong(ieee754dp x); + +int ieee754dp_cmp(ieee754dp x, ieee754dp y, int cop, int sig); +/* + * basic sp math + */ +ieee754dp ieee754dp_modf(ieee754dp x, ieee754dp * ip); +ieee754dp ieee754dp_frexp(ieee754dp x, int *exp); +ieee754dp ieee754dp_ldexp(ieee754dp x, int exp); + +ieee754dp ieee754dp_ceil(ieee754dp x); +ieee754dp ieee754dp_floor(ieee754dp x); +ieee754dp ieee754dp_trunc(ieee754dp x); + +ieee754dp ieee754dp_sqrt(ieee754dp x); + + + +/* 5 types of floating point number +*/ +#define IEEE754_CLASS_NORM 0x00 +#define IEEE754_CLASS_ZERO 0x01 +#define IEEE754_CLASS_DNORM 0x02 +#define IEEE754_CLASS_INF 0x03 +#define IEEE754_CLASS_SNAN 0x04 +#define IEEE754_CLASS_QNAN 0x05 +extern const char *const ieee754_cname[]; + +/* 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 + +/* rounding mode +*/ +#define IEEE754_RN 0 /* round to nearest */ +#define IEEE754_RZ 1 /* round toward zero */ +#define IEEE754_RD 2 /* round toward -Infinity */ +#define IEEE754_RU 3 /* round toward +Infinity */ + +/* other naming */ +#define IEEE754_RM IEEE754_RD +#define IEEE754_RP IEEE754_RU + +/* "normal" comparisons +*/ +static __inline int ieee754sp_eq(ieee754sp x, ieee754sp y) +{ + return ieee754sp_cmp(x, y, IEEE754_CEQ, 0); +} + +static __inline int ieee754sp_ne(ieee754sp x, ieee754sp y) +{ + return ieee754sp_cmp(x, y, + IEEE754_CLT | IEEE754_CGT | IEEE754_CUN, 0); +} + +static __inline int ieee754sp_lt(ieee754sp x, ieee754sp y) +{ + return ieee754sp_cmp(x, y, IEEE754_CLT, 0); +} + +static __inline int ieee754sp_le(ieee754sp x, ieee754sp y) +{ + return ieee754sp_cmp(x, y, IEEE754_CLT | IEEE754_CEQ, 0); +} + +static __inline int ieee754sp_gt(ieee754sp x, ieee754sp y) +{ + return ieee754sp_cmp(x, y, IEEE754_CGT, 0); +} + + +static __inline int ieee754sp_ge(ieee754sp x, ieee754sp y) +{ + return ieee754sp_cmp(x, y, IEEE754_CGT | IEEE754_CEQ, 0); +} + +static __inline int ieee754dp_eq(ieee754dp x, ieee754dp y) +{ + return ieee754dp_cmp(x, y, IEEE754_CEQ, 0); +} + +static __inline int ieee754dp_ne(ieee754dp x, ieee754dp y) +{ + return ieee754dp_cmp(x, y, + IEEE754_CLT | IEEE754_CGT | IEEE754_CUN, 0); +} + +static __inline int ieee754dp_lt(ieee754dp x, ieee754dp y) +{ + return ieee754dp_cmp(x, y, IEEE754_CLT, 0); +} + +static __inline int ieee754dp_le(ieee754dp x, ieee754dp y) +{ + return ieee754dp_cmp(x, y, IEEE754_CLT | IEEE754_CEQ, 0); +} + +static __inline int ieee754dp_gt(ieee754dp x, ieee754dp y) +{ + return ieee754dp_cmp(x, y, IEEE754_CGT, 0); +} + +static __inline int ieee754dp_ge(ieee754dp x, ieee754dp y) +{ + return ieee754dp_cmp(x, y, IEEE754_CGT | IEEE754_CEQ, 0); +} + + +/* like strtod +*/ +ieee754dp ieee754dp_fstr(const char *s, char **endp); +char *ieee754dp_tstr(ieee754dp x, int prec, int fmt, int af); + + +/* the control status register +*/ +struct ieee754_csr { + unsigned pad:13; + unsigned nod:1; /* set 1 for no denormalised numbers */ + unsigned cx:5; /* exceptions this operation */ + unsigned mx:5; /* exception enable mask */ + unsigned sx:5; /* exceptions total */ + unsigned rm:2; /* current rounding mode */ +}; +extern struct ieee754_csr ieee754_csr; + +static __inline unsigned ieee754_getrm(void) +{ + return (ieee754_csr.rm); +} +static __inline unsigned ieee754_setrm(unsigned rm) +{ + return (ieee754_csr.rm = rm); +} + +/* + * get current exceptions + */ +static __inline unsigned ieee754_getcx(void) +{ + return (ieee754_csr.cx); +} + +/* test for current exception condition + */ +static __inline int ieee754_cxtest(unsigned n) +{ + return (ieee754_csr.cx & n); +} + +/* + * get sticky exceptions + */ +static __inline unsigned ieee754_getsx(void) +{ + return (ieee754_csr.sx); +} + +/* clear sticky conditions +*/ +static __inline unsigned ieee754_clrsx(void) +{ + return (ieee754_csr.sx = 0); +} + +/* test for sticky exception condition + */ +static __inline int ieee754_sxtest(unsigned n) +{ + return (ieee754_csr.sx & n); +} + +/* debugging */ +ieee754sp ieee754sp_dump(char *s, ieee754sp x); +ieee754dp ieee754dp_dump(char *s, ieee754dp x); + +#define IEEE754_SPCVAL_PZERO 0 +#define IEEE754_SPCVAL_NZERO 1 +#define IEEE754_SPCVAL_PONE 2 +#define IEEE754_SPCVAL_NONE 3 +#define IEEE754_SPCVAL_PTEN 4 +#define IEEE754_SPCVAL_NTEN 5 +#define IEEE754_SPCVAL_PINFINITY 6 +#define IEEE754_SPCVAL_NINFINITY 7 +#define IEEE754_SPCVAL_INDEF 8 +#define IEEE754_SPCVAL_PMAX 9 /* +max norm */ +#define IEEE754_SPCVAL_NMAX 10 /* -max norm */ +#define IEEE754_SPCVAL_PMIN 11 /* +min norm */ +#define IEEE754_SPCVAL_NMIN 12 /* +min norm */ +#define IEEE754_SPCVAL_PMIND 13 /* +min denorm */ +#define IEEE754_SPCVAL_NMIND 14 /* +min denorm */ +#define IEEE754_SPCVAL_P1E31 15 /* + 1.0e31 */ +#define IEEE754_SPCVAL_P1E63 16 /* + 1.0e63 */ + +extern const struct ieee754dp_konst __ieee754dp_spcvals[]; +extern const struct ieee754sp_konst __ieee754sp_spcvals[]; +#define ieee754dp_spcvals ((const ieee754dp *)__ieee754dp_spcvals) +#define ieee754sp_spcvals ((const 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]) +#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]) +#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 +*/ +#define ieee754si_indef() INT_MAX +#ifdef LONG_LONG_MAX +#define ieee754di_indef() LONG_LONG_MAX +#else +#define ieee754di_indef() ((s64)(~0ULL>>1)) +#endif + +/* IEEE exception context, passed to handler */ +struct ieee754xctx { + const char *op; /* operation name */ + int rt; /* result type */ + union { + ieee754sp sp; /* single precision */ + ieee754dp dp; /* double precision */ +#ifdef IEEE854_XP + ieee754xp xp; /* extended precision */ +#endif + int si; /* standard signed integer (32bits) */ + s64 di; /* extended signed integer (64bits) */ + } rv; /* default result format implied by op */ + va_list ap; +}; + +/* 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 + +extern void ieee754_xcpt(struct ieee754xctx *xcp); + +/* compat */ +#define ieee754dp_fix(x) ieee754dp_tint(x) +#define ieee754sp_fix(x) ieee754sp_tint(x) diff --git a/arch/mips/math-emu/ieee754d.c b/arch/mips/math-emu/ieee754d.c new file mode 100644 index 0000000..7e900f3 --- /dev/null +++ b/arch/mips/math-emu/ieee754d.c @@ -0,0 +1,138 @@ +/* + * Some debug functions + * + * MIPS floating point support + * + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * 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/kernel.h> +#include "ieee754.h" + +#define DP_EBIAS 1023 +#define DP_EMIN (-1022) +#define DP_EMAX 1023 +#define DP_FBITS 52 + +#define SP_EBIAS 127 +#define SP_EMIN (-126) +#define SP_EMAX 127 +#define SP_FBITS 23 + +#define DP_MBIT(x) ((u64)1 << (x)) +#define DP_HIDDEN_BIT DP_MBIT(DP_FBITS) +#define DP_SIGN_BIT DP_MBIT(63) + + +#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.parts.sign) +#define SPBEXP(sp) (sp.parts.bexp) +#define SPMANT(sp) (sp.parts.mant) + +#define DPSIGN(dp) (dp.parts.sign) +#define DPBEXP(dp) (dp.parts.bexp) +#define DPMANT(dp) (dp.parts.mant) + +ieee754dp ieee754dp_dump(char *m, 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; +} + +ieee754sp ieee754sp_dump(char *m, 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; +} + diff --git a/arch/mips/math-emu/ieee754dp.c b/arch/mips/math-emu/ieee754dp.c new file mode 100644 index 0000000..3e214aa --- /dev/null +++ b/arch/mips/math-emu/ieee754dp.c @@ -0,0 +1,243 @@ +/* IEEE754 floating point arithmetic + * double precision: common utilities + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754dp.h" + +int ieee754dp_class(ieee754dp x) +{ + COMPXDP; + EXPLODEXDP; + return xc; +} + +int ieee754dp_isnan(ieee754dp x) +{ + return ieee754dp_class(x) >= IEEE754_CLASS_SNAN; +} + +int ieee754dp_issnan(ieee754dp x) +{ + assert(ieee754dp_isnan(x)); + return ((DPMANT(x) & DP_MBIT(DP_MBITS-1)) == DP_MBIT(DP_MBITS-1)); +} + + +ieee754dp ieee754dp_xcpt(ieee754dp r, const char *op, ...) +{ + struct ieee754xctx ax; + if (!TSTX()) + return r; + + ax.op = op; + ax.rt = IEEE754_RT_DP; + ax.rv.dp = r; + va_start(ax.ap, op); + ieee754_xcpt(&ax); + return ax.rv.dp; +} + +ieee754dp ieee754dp_nanxcpt(ieee754dp r, const char *op, ...) +{ + struct ieee754xctx ax; + + assert(ieee754dp_isnan(r)); + + if (!ieee754dp_issnan(r)) /* QNAN does not cause invalid op !! */ + return r; + + if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) { + /* not enabled convert to a quiet NaN */ + DPMANT(r) &= (~DP_MBIT(DP_MBITS-1)); + if (ieee754dp_isnan(r)) + return r; + else + return ieee754dp_indef(); + } + + ax.op = op; + ax.rt = 0; + ax.rv.dp = r; + va_start(ax.ap, op); + ieee754_xcpt(&ax); + return ax.rv.dp; +} + +ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y) +{ + assert(ieee754dp_isnan(x)); + assert(ieee754dp_isnan(y)); + + if (DPMANT(x) > DPMANT(y)) + return x; + else + return y; +} + + +static u64 get_rounding(int sn, u64 xm) +{ + /* inexact must round of 3 bits + */ + if (xm & (DP_MBIT(3) - 1)) { + switch (ieee754_csr.rm) { + case IEEE754_RZ: + break; + case IEEE754_RN: + xm += 0x3 + ((xm >> 3) & 1); + /* xm += (xm&0x8)?0x4:0x3 */ + break; + case IEEE754_RU: /* toward +Infinity */ + if (!sn) /* ?? */ + xm += 0x8; + break; + case IEEE754_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. + */ +ieee754dp ieee754dp_format(int sn, int xe, u64 xm) +{ + assert(xm); /* we don't gen exact zeros (probably should) */ + + assert((xm >> (DP_MBITS + 1 + 3)) == 0); /* no execess */ + assert(xm & (DP_HIDDEN_BIT << 3)); + + if (xe < DP_EMIN) { + /* strip lower bits */ + int es = DP_EMIN - xe; + + if (ieee754_csr.nod) { + SETCX(IEEE754_UNDERFLOW); + SETCX(IEEE754_INEXACT); + + switch(ieee754_csr.rm) { + case IEEE754_RN: + return ieee754dp_zero(sn); + case IEEE754_RZ: + return ieee754dp_zero(sn); + case IEEE754_RU: /* toward +Infinity */ + if(sn == 0) + return ieee754dp_min(0); + else + return ieee754dp_zero(1); + case IEEE754_RD: /* toward -Infinity */ + if(sn == 0) + return ieee754dp_zero(0); + else + return ieee754dp_min(1); + } + } + + if (xe == DP_EMIN - 1 + && get_rounding(sn, xm) >> (DP_MBITS + 1 + 3)) + { + /* Not tiny after rounding */ + SETCX(IEEE754_INEXACT); + xm = 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)) { + SETCX(IEEE754_INEXACT); + if ((xm & (DP_HIDDEN_BIT << 3)) == 0) { + SETCX(IEEE754_UNDERFLOW); + } + + /* inexact must round of 3 bits + */ + xm = get_rounding(sn, xm); + /* adjust exponent for rounding add overflowing + */ + if (xm >> (DP_MBITS + 3 + 1)) { + /* add causes mantissa overflow */ + xm >>= 1; + xe++; + } + } + /* strip grs bits */ + xm >>= 3; + + assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ + assert(xe >= DP_EMIN); + + if (xe > DP_EMAX) { + SETCX(IEEE754_OVERFLOW); + SETCX(IEEE754_INEXACT); + /* -O can be table indexed by (rm,sn) */ + switch (ieee754_csr.rm) { + case IEEE754_RN: + return ieee754dp_inf(sn); + case IEEE754_RZ: + return ieee754dp_max(sn); + case IEEE754_RU: /* toward +Infinity */ + if (sn == 0) + return ieee754dp_inf(0); + else + return ieee754dp_max(1); + case IEEE754_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) + SETCX(IEEE754_UNDERFLOW); + return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); + } else { + assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ + assert(xm & DP_HIDDEN_BIT); + + return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); + } +} diff --git a/arch/mips/math-emu/ieee754dp.h b/arch/mips/math-emu/ieee754dp.h new file mode 100644 index 0000000..a37370d --- /dev/null +++ b/arch/mips/math-emu/ieee754dp.h @@ -0,0 +1,83 @@ +/* + * IEEE754 floating point + * double precision internal header file + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754int.h" + +#define assert(expr) ((void)0) + +/* 3bit extended double precision sticky right shift */ +#define XDPSRS(v,rs) \ + ((rs > (DP_MBITS+3))?1:((v) >> (rs)) | ((v) << (64-(rs)) != 0)) + +#define XDPSRSX1() \ + (xe++, (xm = (xm >> 1) | (xm & 1))) + +#define XDPSRS1(v) \ + (((v) >> 1) | ((v) & 1)) + +/* convert denormal to normalized with extended exponent */ +#define DPDNORMx(m,e) \ + while( (m >> DP_MBITS) == 0) { m <<= 1; e--; } +#define DPDNORMX DPDNORMx(xm,xe) +#define DPDNORMY DPDNORMx(ym,ye) + +static __inline ieee754dp builddp(int s, int bx, u64 m) +{ + ieee754dp r; + + assert((s) == 0 || (s) == 1); + assert((bx) >= DP_EMIN - 1 + DP_EBIAS + && (bx) <= DP_EMAX + 1 + DP_EBIAS); + assert(((m) >> DP_MBITS) == 0); + + r.parts.sign = s; + r.parts.bexp = bx; + r.parts.mant = m; + return r; +} + +extern int ieee754dp_isnan(ieee754dp); +extern int ieee754dp_issnan(ieee754dp); +extern int ieee754si_xcpt(int, const char *, ...); +extern s64 ieee754di_xcpt(s64, const char *, ...); +extern ieee754dp ieee754dp_xcpt(ieee754dp, const char *, ...); +extern ieee754dp ieee754dp_nanxcpt(ieee754dp, const char *, ...); +extern ieee754dp ieee754dp_bestnan(ieee754dp, ieee754dp); +extern ieee754dp ieee754dp_format(int, int, u64); + + +#define DPNORMRET2(s,e,m,name,a0,a1) \ +{ \ + ieee754dp V = ieee754dp_format(s,e,m); \ + if(TSTX()) \ + return ieee754dp_xcpt(V,name,a0,a1); \ + else \ + return V; \ +} + +#define DPNORMRET1(s,e,m,name,a0) DPNORMRET2(s,e,m,name,a0,a0) diff --git a/arch/mips/math-emu/ieee754int.h b/arch/mips/math-emu/ieee754int.h new file mode 100644 index 0000000..4a5a81d --- /dev/null +++ b/arch/mips/math-emu/ieee754int.h @@ -0,0 +1,165 @@ +/* + * IEEE754 floating point + * common internal header file + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754.h" + +#define DP_EBIAS 1023 +#define DP_EMIN (-1022) +#define DP_EMAX 1023 +#define DP_MBITS 52 + +#define SP_EBIAS 127 +#define SP_EMIN (-126) +#define SP_EMAX 127 +#define SP_MBITS 23 + +#define DP_MBIT(x) ((u64)1 << (x)) +#define DP_HIDDEN_BIT DP_MBIT(DP_MBITS) +#define DP_SIGN_BIT DP_MBIT(63) + +#define SP_MBIT(x) ((u32)1 << (x)) +#define SP_HIDDEN_BIT SP_MBIT(SP_MBITS) +#define SP_SIGN_BIT SP_MBIT(31) + + +#define SPSIGN(sp) (sp.parts.sign) +#define SPBEXP(sp) (sp.parts.bexp) +#define SPMANT(sp) (sp.parts.mant) + +#define DPSIGN(dp) (dp.parts.sign) +#define DPBEXP(dp) (dp.parts.bexp) +#define DPMANT(dp) (dp.parts.mant) + +#define CLPAIR(x,y) ((x)*6+(y)) + +#define CLEARCX \ + (ieee754_csr.cx = 0) + +#define SETCX(x) \ + (ieee754_csr.cx |= (x),ieee754_csr.sx |= (x)) + +#define SETANDTESTCX(x) \ + (SETCX(x),ieee754_csr.mx & (x)) + +#define TSTX() \ + (ieee754_csr.cx & ieee754_csr.mx) + + +#define COMPXSP \ + unsigned xm; int xe; int xs; int xc + +#define COMPYSP \ + unsigned ym; int ye; int ys; int yc + +#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(vm & SP_MBIT(SP_MBITS-1)) \ + vc = IEEE754_CLASS_SNAN;\ + else \ + vc = IEEE754_CLASS_QNAN;\ + } 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 COMPXDP \ +u64 xm; int xe; int xs; int xc + +#define COMPYDP \ +u64 ym; int ye; int ys; int yc + +#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(vm & DP_MBIT(DP_MBITS-1)) \ + vc = IEEE754_CLASS_SNAN;\ + else \ + vc = IEEE754_CLASS_QNAN;\ + } 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 FLUSHDP(v,vc,vs,ve,vm) \ + if(vc==IEEE754_CLASS_DNORM) {\ + if(ieee754_csr.nod) {\ + 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) {\ + 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 FLUSHXSP FLUSHSP(x,xc,xs,xe,xm) +#define FLUSHYSP FLUSHSP(y,yc,ys,ye,ym) diff --git a/arch/mips/math-emu/ieee754m.c b/arch/mips/math-emu/ieee754m.c new file mode 100644 index 0000000..d66896c --- /dev/null +++ b/arch/mips/math-emu/ieee754m.c @@ -0,0 +1,56 @@ +/* + * floor, trunc, ceil + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754.h" + +ieee754dp ieee754dp_floor(ieee754dp x) +{ + ieee754dp i; + + if (ieee754dp_lt(ieee754dp_modf(x, &i), ieee754dp_zero(0))) + return ieee754dp_sub(i, ieee754dp_one(0)); + else + return i; +} + +ieee754dp ieee754dp_ceil(ieee754dp x) +{ + ieee754dp i; + + if (ieee754dp_gt(ieee754dp_modf(x, &i), ieee754dp_zero(0))) + return ieee754dp_add(i, ieee754dp_one(0)); + else + return i; +} + +ieee754dp ieee754dp_trunc(ieee754dp x) +{ + ieee754dp i; + + (void) ieee754dp_modf(x, &i); + return i; +} diff --git a/arch/mips/math-emu/ieee754sp.c b/arch/mips/math-emu/ieee754sp.c new file mode 100644 index 0000000..adda851 --- /dev/null +++ b/arch/mips/math-emu/ieee754sp.c @@ -0,0 +1,243 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +int ieee754sp_class(ieee754sp x) +{ + COMPXSP; + EXPLODEXSP; + return xc; +} + +int ieee754sp_isnan(ieee754sp x) +{ + return ieee754sp_class(x) >= IEEE754_CLASS_SNAN; +} + +int ieee754sp_issnan(ieee754sp x) +{ + assert(ieee754sp_isnan(x)); + return (SPMANT(x) & SP_MBIT(SP_MBITS-1)); +} + + +ieee754sp ieee754sp_xcpt(ieee754sp r, const char *op, ...) +{ + struct ieee754xctx ax; + + if (!TSTX()) + return r; + + ax.op = op; + ax.rt = IEEE754_RT_SP; + ax.rv.sp = r; + va_start(ax.ap, op); + ieee754_xcpt(&ax); + return ax.rv.sp; +} + +ieee754sp ieee754sp_nanxcpt(ieee754sp r, const char *op, ...) +{ + struct ieee754xctx ax; + + assert(ieee754sp_isnan(r)); + + if (!ieee754sp_issnan(r)) /* QNAN does not cause invalid op !! */ + return r; + + if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) { + /* not enabled convert to a quiet NaN */ + SPMANT(r) &= (~SP_MBIT(SP_MBITS-1)); + if (ieee754sp_isnan(r)) + return r; + else + return ieee754sp_indef(); + } + + ax.op = op; + ax.rt = 0; + ax.rv.sp = r; + va_start(ax.ap, op); + ieee754_xcpt(&ax); + return ax.rv.sp; +} + +ieee754sp ieee754sp_bestnan(ieee754sp x, ieee754sp y) +{ + assert(ieee754sp_isnan(x)); + assert(ieee754sp_isnan(y)); + + if (SPMANT(x) > SPMANT(y)) + return x; + else + return y; +} + + +static unsigned get_rounding(int sn, unsigned xm) +{ + /* inexact must round of 3 bits + */ + if (xm & (SP_MBIT(3) - 1)) { + switch (ieee754_csr.rm) { + case IEEE754_RZ: + break; + case IEEE754_RN: + xm += 0x3 + ((xm >> 3) & 1); + /* xm += (xm&0x8)?0x4:0x3 */ + break; + case IEEE754_RU: /* toward +Infinity */ + if (!sn) /* ?? */ + xm += 0x8; + break; + case IEEE754_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. + */ +ieee754sp ieee754sp_format(int sn, int xe, unsigned xm) +{ + assert(xm); /* we don't gen exact zeros (probably should) */ + + assert((xm >> (SP_MBITS + 1 + 3)) == 0); /* no execess */ + assert(xm & (SP_HIDDEN_BIT << 3)); + + if (xe < SP_EMIN) { + /* strip lower bits */ + int es = SP_EMIN - xe; + + if (ieee754_csr.nod) { + SETCX(IEEE754_UNDERFLOW); + SETCX(IEEE754_INEXACT); + + switch(ieee754_csr.rm) { + case IEEE754_RN: + return ieee754sp_zero(sn); + case IEEE754_RZ: + return ieee754sp_zero(sn); + case IEEE754_RU: /* toward +Infinity */ + if(sn == 0) + return ieee754sp_min(0); + else + return ieee754sp_zero(1); + case IEEE754_RD: /* toward -Infinity */ + if(sn == 0) + return ieee754sp_zero(0); + else + return ieee754sp_min(1); + } + } + + if (xe == SP_EMIN - 1 + && get_rounding(sn, xm) >> (SP_MBITS + 1 + 3)) + { + /* Not tiny after rounding */ + SETCX(IEEE754_INEXACT); + xm = get_rounding(sn, xm); + xm >>= 1; + /* Clear grs bits */ + xm &= ~(SP_MBIT(3) - 1); + xe++; + } + else { + /* sticky right shift es bits + */ + SPXSRSXn(es); + assert((xm & (SP_HIDDEN_BIT << 3)) == 0); + assert(xe == SP_EMIN); + } + } + if (xm & (SP_MBIT(3) - 1)) { + SETCX(IEEE754_INEXACT); + if ((xm & (SP_HIDDEN_BIT << 3)) == 0) { + SETCX(IEEE754_UNDERFLOW); + } + + /* inexact must round of 3 bits + */ + xm = get_rounding(sn, xm); + /* adjust exponent for rounding add overflowing + */ + if (xm >> (SP_MBITS + 1 + 3)) { + /* add causes mantissa overflow */ + xm >>= 1; + xe++; + } + } + /* strip grs bits */ + xm >>= 3; + + assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */ + assert(xe >= SP_EMIN); + + if (xe > SP_EMAX) { + SETCX(IEEE754_OVERFLOW); + SETCX(IEEE754_INEXACT); + /* -O can be table indexed by (rm,sn) */ + switch (ieee754_csr.rm) { + case IEEE754_RN: + return ieee754sp_inf(sn); + case IEEE754_RZ: + return ieee754sp_max(sn); + case IEEE754_RU: /* toward +Infinity */ + if (sn == 0) + return ieee754sp_inf(0); + else + return ieee754sp_max(1); + case IEEE754_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) + SETCX(IEEE754_UNDERFLOW); + return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm); + } else { + assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */ + assert(xm & SP_HIDDEN_BIT); + + return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT); + } +} diff --git a/arch/mips/math-emu/ieee754sp.h b/arch/mips/math-emu/ieee754sp.h new file mode 100644 index 0000000..ae82f51 --- /dev/null +++ b/arch/mips/math-emu/ieee754sp.h @@ -0,0 +1,89 @@ +/* + * IEEE754 floating point + * double precision internal header file + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754int.h" + +#define assert(expr) ((void)0) + +/* 3bit extended single precision sticky right shift */ +#define SPXSRSXn(rs) \ + (xe += rs, \ + xm = (rs > (SP_MBITS+3))?1:((xm) >> (rs)) | ((xm) << (32-(rs)) != 0)) + +#define SPXSRSX1() \ + (xe++, (xm = (xm >> 1) | (xm & 1))) + +#define SPXSRSYn(rs) \ + (ye+=rs, \ + ym = (rs > (SP_MBITS+3))?1:((ym) >> (rs)) | ((ym) << (32-(rs)) != 0)) + +#define SPXSRSY1() \ + (ye++, (ym = (ym >> 1) | (ym & 1))) + +/* convert denormal to normalized with extended exponent */ +#define SPDNORMx(m,e) \ + while( (m >> SP_MBITS) == 0) { m <<= 1; e--; } +#define SPDNORMX SPDNORMx(xm,xe) +#define SPDNORMY SPDNORMx(ym,ye) + +static __inline ieee754sp buildsp(int s, int bx, unsigned m) +{ + ieee754sp r; + + assert((s) == 0 || (s) == 1); + assert((bx) >= SP_EMIN - 1 + SP_EBIAS + && (bx) <= SP_EMAX + 1 + SP_EBIAS); + assert(((m) >> SP_MBITS) == 0); + + r.parts.sign = s; + r.parts.bexp = bx; + r.parts.mant = m; + + return r; +} + +extern int ieee754sp_isnan(ieee754sp); +extern int ieee754sp_issnan(ieee754sp); +extern int ieee754si_xcpt(int, const char *, ...); +extern s64 ieee754di_xcpt(s64, const char *, ...); +extern ieee754sp ieee754sp_xcpt(ieee754sp, const char *, ...); +extern ieee754sp ieee754sp_nanxcpt(ieee754sp, const char *, ...); +extern ieee754sp ieee754sp_bestnan(ieee754sp, ieee754sp); +extern ieee754sp ieee754sp_format(int, int, unsigned); + + +#define SPNORMRET2(s,e,m,name,a0,a1) \ +{ \ + ieee754sp V = ieee754sp_format(s,e,m); \ + if(TSTX()) \ + return ieee754sp_xcpt(V,name,a0,a1); \ + else \ + return V; \ +} + +#define SPNORMRET1(s,e,m,name,a0) SPNORMRET2(s,e,m,name,a0,a0) diff --git a/arch/mips/math-emu/ieee754xcpt.c b/arch/mips/math-emu/ieee754xcpt.c new file mode 100644 index 0000000..7d8ef89 --- /dev/null +++ b/arch/mips/math-emu/ieee754xcpt.c @@ -0,0 +1,49 @@ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + +/************************************************************************** + * Nov 7, 2000 + * Added preprocessor hacks to map to Linux kernel diagnostics. + * + * Kevin D. Kissell, kevink@mips.com and Carsten Langgaard, carstenl@mips.com + * Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved. + *************************************************************************/ + +#include <linux/kernel.h> +#include "ieee754.h" + +/* + * Very naff exception handler (you can plug in your own and + * override this). + */ + +static const char *const rtnames[] = { + "sp", "dp", "xp", "si", "di" +}; + +void ieee754_xcpt(struct ieee754xctx *xcp) +{ + printk(KERN_DEBUG "floating point exception in \"%s\", type=%s\n", + xcp->op, rtnames[xcp->rt]); +} + diff --git a/arch/mips/math-emu/kernel_linkage.c b/arch/mips/math-emu/kernel_linkage.c new file mode 100644 index 0000000..04397fec --- /dev/null +++ b/arch/mips/math-emu/kernel_linkage.c @@ -0,0 +1,125 @@ +/* + * Kevin D. Kissell, kevink@mips and Carsten Langgaard, carstenl@mips.com + * Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved. + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * Routines corresponding to Linux kernel FP context + * manipulation primitives for the Algorithmics MIPS + * FPU Emulator + */ +#include <linux/config.h> +#include <linux/sched.h> +#include <asm/processor.h> +#include <asm/signal.h> +#include <asm/uaccess.h> + +#include <asm/fpu_emulator.h> + +extern struct mips_fpu_emulator_private fpuemuprivate; + +#define SIGNALLING_NAN 0x7ff800007ff80000LL + +void fpu_emulator_init_fpu(void) +{ + static int first = 1; + int i; + + if (first) { + first = 0; + printk("Algorithmics/MIPS FPU Emulator v1.5\n"); + } + + current->thread.fpu.soft.fcr31 = 0; + for (i = 0; i < 32; i++) { + current->thread.fpu.soft.fpr[i] = SIGNALLING_NAN; + } +} + + +/* + * Emulator context save/restore to/from a signal context + * presumed to be on the user stack, and therefore accessed + * with appropriate macros from uaccess.h + */ + +int fpu_emulator_save_context(struct sigcontext *sc) +{ + int i; + int err = 0; + + for (i = 0; i < 32; i++) { + err |= + __put_user(current->thread.fpu.soft.fpr[i], + &sc->sc_fpregs[i]); + } + err |= __put_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr); + err |= __put_user(fpuemuprivate.eir, &sc->sc_fpc_eir); + + return err; +} + +int fpu_emulator_restore_context(struct sigcontext *sc) +{ + int i; + int err = 0; + + for (i = 0; i < 32; i++) { + err |= + __get_user(current->thread.fpu.soft.fpr[i], + &sc->sc_fpregs[i]); + } + err |= __get_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr); + err |= __get_user(fpuemuprivate.eir, &sc->sc_fpc_eir); + + return err; +} + +#ifdef CONFIG_MIPS64 +/* + * This is the o32 version + */ + +int fpu_emulator_save_context32(struct sigcontext32 *sc) +{ + int i; + int err = 0; + + for (i = 0; i < 32; i+=2) { + err |= + __put_user(current->thread.fpu.soft.fpr[i], + &sc->sc_fpregs[i]); + } + err |= __put_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr); + err |= __put_user(fpuemuprivate.eir, &sc->sc_fpc_eir); + + return err; +} + +int fpu_emulator_restore_context32(struct sigcontext32 *sc) +{ + int i; + int err = 0; + + for (i = 0; i < 32; i+=2) { + err |= + __get_user(current->thread.fpu.soft.fpr[i], + &sc->sc_fpregs[i]); + } + err |= __get_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr); + err |= __get_user(fpuemuprivate.eir, &sc->sc_fpc_eir); + + return err; +} +#endif diff --git a/arch/mips/math-emu/sp_add.c b/arch/mips/math-emu/sp_add.c new file mode 100644 index 0000000..d8c4211 --- /dev/null +++ b/arch/mips/math-emu/sp_add.c @@ -0,0 +1,177 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_add(ieee754sp x, ieee754sp y) +{ + COMPXSP; + COMPYSP; + + EXPLODEXSP; + EXPLODEYSP; + + CLEARCX; + + FLUSHXSP; + FLUSHYSP; + + switch (CLPAIR(xc, yc)) { + case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): + case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): + case CLPAIR(IEEE754_CLASS_SNAN, 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): + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "add", x, y); + + 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; + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_xcpt(ieee754sp_indef(), "add", x, y); + + 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 == + IEEE754_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; + + 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 + */ + int s = xe - ye; + SPXSRSYn(s); + } else if (ye > xe) { + /* have to shift x fraction right to align + */ + int s = ye - xe; + SPXSRSXn(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,xe + */ + xm = xm + ym; + xe = xe; + xs = xs; + + if (xm >> (SP_MBITS + 1 + 3)) { /* carry out */ + SPXSRSX1(); + } + } else { + if (xm >= ym) { + xm = xm - ym; + xe = xe; + xs = xs; + } else { + xm = ym - xm; + xe = xe; + xs = ys; + } + if (xm == 0) + return ieee754sp_zero(ieee754_csr.rm == + IEEE754_RD); + + /* normalize in extended single precision */ + while ((xm >> (SP_MBITS + 3)) == 0) { + xm <<= 1; + xe--; + } + + } + SPNORMRET2(xs, xe, xm, "add", x, y); +} diff --git a/arch/mips/math-emu/sp_cmp.c b/arch/mips/math-emu/sp_cmp.c new file mode 100644 index 0000000..d3eff6b --- /dev/null +++ b/arch/mips/math-emu/sp_cmp.c @@ -0,0 +1,67 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +int ieee754sp_cmp(ieee754sp x, ieee754sp y, int cmp, int sig) +{ + COMPXSP; + COMPYSP; + + EXPLODEXSP; + EXPLODEYSP; + FLUSHXSP; + FLUSHYSP; + CLEARCX; /* Even clear inexact flag here */ + + if (ieee754sp_isnan(x) || ieee754sp_isnan(y)) { + if (sig || xc == IEEE754_CLASS_SNAN || yc == IEEE754_CLASS_SNAN) + SETCX(IEEE754_INVALID_OPERATION); + if (cmp & IEEE754_CUN) + return 1; + if (cmp & (IEEE754_CLT | IEEE754_CGT)) { + if (sig && SETANDTESTCX(IEEE754_INVALID_OPERATION)) + return ieee754si_xcpt(0, "fcmpf", x); + } + return 0; + } else { + int vx = x.bits; + int 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; + } +} diff --git a/arch/mips/math-emu/sp_div.c b/arch/mips/math-emu/sp_div.c new file mode 100644 index 0000000..2b437fc --- /dev/null +++ b/arch/mips/math-emu/sp_div.c @@ -0,0 +1,157 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_div(ieee754sp x, ieee754sp y) +{ + COMPXSP; + COMPYSP; + + EXPLODEXSP; + EXPLODEYSP; + + CLEARCX; + + FLUSHXSP; + FLUSHYSP; + + switch (CLPAIR(xc, yc)) { + case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): + case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): + case CLPAIR(IEEE754_CLASS_SNAN, 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): + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "div", x, y); + + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_xcpt(ieee754sp_indef(), "div", x, y); + + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_xcpt(ieee754sp_indef(), "div", x, y); + + case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO): + case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO): + SETCX(IEEE754_ZERO_DIVIDE); + return ieee754sp_xcpt(ieee754sp_inf(xs ^ ys), "div", x, y); + + 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; + + 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 */ + + unsigned rm = 0; + int re = xe - ye; + unsigned bm; + + for (bm = SP_MBIT(SP_MBITS + 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_MBITS + 3)) == 0) { + rm <<= 1; + re--; + } + + SPNORMRET2(xs == ys ? 0 : 1, re, rm, "div", x, y); + } +} diff --git a/arch/mips/math-emu/sp_fdp.c b/arch/mips/math-emu/sp_fdp.c new file mode 100644 index 0000000..4093723 --- /dev/null +++ b/arch/mips/math-emu/sp_fdp.c @@ -0,0 +1,77 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_fdp(ieee754dp x) +{ + COMPXDP; + ieee754sp nan; + + EXPLODEXDP; + + CLEARCX; + + FLUSHXDP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "fdp"); + case IEEE754_CLASS_QNAN: + nan = buildsp(xs, SP_EMAX + 1 + SP_EBIAS, (u32) + (xm >> (DP_MBITS - SP_MBITS))); + if (!ieee754sp_isnan(nan)) + nan = ieee754sp_indef(); + return ieee754sp_nanxcpt(nan, "fdp", x); + 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 */ + SETCX(IEEE754_UNDERFLOW); + SETCX(IEEE754_INEXACT); + if ((ieee754_csr.rm == IEEE754_RU && !xs) || + (ieee754_csr.rm == IEEE754_RD && xs)) + return ieee754sp_xcpt(ieee754sp_mind(xs), "fdp", x); + return ieee754sp_xcpt(ieee754sp_zero(xs), "fdp", x); + case IEEE754_CLASS_NORM: + break; + } + + { + u32 rm; + + /* convert from DP_MBITS to SP_MBITS+3 with sticky right shift + */ + rm = (xm >> (DP_MBITS - (SP_MBITS + 3))) | + ((xm << (64 - (DP_MBITS - (SP_MBITS + 3)))) != 0); + + SPNORMRET1(xs, xe, rm, "fdp", x); + } +} diff --git a/arch/mips/math-emu/sp_fint.c b/arch/mips/math-emu/sp_fint.c new file mode 100644 index 0000000..42d9ed4 --- /dev/null +++ b/arch/mips/math-emu/sp_fint.c @@ -0,0 +1,80 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_fint(int x) +{ + COMPXSP; + + CLEARCX; + + xc = ( 0 ? xc : xc ); + + 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_MBITS + 3; + + if (xm >> (SP_MBITS + 1 + 3)) { + /* shunt out overflow bits + */ + while (xm >> (SP_MBITS + 1 + 3)) { + SPXSRSX1(); + } + } else { + /* normalize in grs extended single precision + */ + while ((xm >> (SP_MBITS + 3)) == 0) { + xm <<= 1; + xe--; + } + } + SPNORMRET1(xs, xe, xm, "fint", x); +} + + +ieee754sp ieee754sp_funs(unsigned int u) +{ + if ((int) u < 0) + return ieee754sp_add(ieee754sp_1e31(), + ieee754sp_fint(u & ~(1 << 31))); + return ieee754sp_fint(u); +} diff --git a/arch/mips/math-emu/sp_flong.c b/arch/mips/math-emu/sp_flong.c new file mode 100644 index 0000000..1e26795 --- /dev/null +++ b/arch/mips/math-emu/sp_flong.c @@ -0,0 +1,79 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_flong(s64 x) +{ + COMPXDP; /* <--- need 64-bit mantissa temp */ + + CLEARCX; + + xc = ( 0 ? xc : xc ); + + 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_MBITS + 3; + + if (xm >> (SP_MBITS + 1 + 3)) { + /* shunt out overflow bits + */ + while (xm >> (SP_MBITS + 1 + 3)) { + SPXSRSX1(); + } + } else { + /* normalize in grs extended single precision */ + while ((xm >> (SP_MBITS + 3)) == 0) { + xm <<= 1; + xe--; + } + } + SPNORMRET1(xs, xe, xm, "sp_flong", x); +} + + +ieee754sp ieee754sp_fulong(u64 u) +{ + if ((s64) u < 0) + return ieee754sp_add(ieee754sp_1e63(), + ieee754sp_flong(u & ~(1ULL << 63))); + return ieee754sp_flong(u); +} diff --git a/arch/mips/math-emu/sp_frexp.c b/arch/mips/math-emu/sp_frexp.c new file mode 100644 index 0000000..359c648 --- /dev/null +++ b/arch/mips/math-emu/sp_frexp.c @@ -0,0 +1,53 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +/* close to ieeep754sp_logb +*/ +ieee754sp ieee754sp_frexp(ieee754sp x, int *eptr) +{ + COMPXSP; + CLEARCX; + EXPLODEXSP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + case IEEE754_CLASS_ZERO: + *eptr = 0; + return x; + case IEEE754_CLASS_DNORM: + SPDNORMX; + break; + case IEEE754_CLASS_NORM: + break; + } + *eptr = xe + 1; + return buildsp(xs, -1 + SP_EBIAS, xm & ~SP_HIDDEN_BIT); +} diff --git a/arch/mips/math-emu/sp_logb.c b/arch/mips/math-emu/sp_logb.c new file mode 100644 index 0000000..3c33721 --- /dev/null +++ b/arch/mips/math-emu/sp_logb.c @@ -0,0 +1,54 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_logb(ieee754sp x) +{ + COMPXSP; + + CLEARCX; + + EXPLODEXSP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + return ieee754sp_nanxcpt(x, "logb", x); + case IEEE754_CLASS_QNAN: + return x; + case IEEE754_CLASS_INF: + return ieee754sp_inf(0); + case IEEE754_CLASS_ZERO: + return ieee754sp_inf(1); + case IEEE754_CLASS_DNORM: + SPDNORMX; + break; + case IEEE754_CLASS_NORM: + break; + } + return ieee754sp_fint(xe); +} diff --git a/arch/mips/math-emu/sp_modf.c b/arch/mips/math-emu/sp_modf.c new file mode 100644 index 0000000..4b1dbac --- /dev/null +++ b/arch/mips/math-emu/sp_modf.c @@ -0,0 +1,80 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +/* modf function is always exact for a finite number +*/ +ieee754sp ieee754sp_modf(ieee754sp x, ieee754sp * ip) +{ + COMPXSP; + + CLEARCX; + + EXPLODEXSP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + case IEEE754_CLASS_ZERO: + *ip = x; + return x; + case IEEE754_CLASS_DNORM: + /* far to small */ + *ip = ieee754sp_zero(xs); + return x; + case IEEE754_CLASS_NORM: + break; + } + if (xe < 0) { + *ip = ieee754sp_zero(xs); + return x; + } + if (xe >= SP_MBITS) { + *ip = x; + return ieee754sp_zero(xs); + } + /* generate ipart mantissa by clearing bottom bits + */ + *ip = buildsp(xs, xe + SP_EBIAS, + ((xm >> (SP_MBITS - xe)) << (SP_MBITS - xe)) & + ~SP_HIDDEN_BIT); + + /* generate fpart mantissa by clearing top bits + * and normalizing (must be able to normalize) + */ + xm = (xm << (32 - (SP_MBITS - xe))) >> (32 - (SP_MBITS - xe)); + if (xm == 0) + return ieee754sp_zero(xs); + + while ((xm >> SP_MBITS) == 0) { + xm <<= 1; + xe--; + } + return buildsp(xs, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT); +} diff --git a/arch/mips/math-emu/sp_mul.c b/arch/mips/math-emu/sp_mul.c new file mode 100644 index 0000000..3f070f8 --- /dev/null +++ b/arch/mips/math-emu/sp_mul.c @@ -0,0 +1,171 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_mul(ieee754sp x, ieee754sp y) +{ + COMPXSP; + COMPYSP; + + EXPLODEXSP; + EXPLODEYSP; + + CLEARCX; + + FLUSHXSP; + FLUSHYSP; + + switch (CLPAIR(xc, yc)) { + case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): + case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): + case CLPAIR(IEEE754_CLASS_SNAN, 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): + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "mul", x, y); + + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_xcpt(ieee754sp_indef(), "mul", x, y); + + 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; + + 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 basicly */ + assert(xm & SP_HIDDEN_BIT); + assert(ym & SP_HIDDEN_BIT); + + { + int re = xe + ye; + int rs = xs ^ ys; + unsigned rm; + + /* shunt to top of word */ + xm <<= 32 - (SP_MBITS + 1); + ym <<= 32 - (SP_MBITS + 1); + + /* multiply 32bits xm,ym to give high 32bits rm with stickness + */ + { + unsigned short lxm = xm & 0xffff; + unsigned short hxm = xm >> 16; + unsigned short lym = ym & 0xffff; + unsigned short hym = ym >> 16; + unsigned lrm; + unsigned hrm; + + lrm = lxm * lym; /* 16 * 16 => 32 */ + hrm = hxm * hym; /* 16 * 16 => 32 */ + + { + unsigned t = lxm * hym; /* 16 * 16 => 32 */ + { + unsigned at = lrm + (t << 16); + hrm += at < lrm; + lrm = at; + } + hrm = hrm + (t >> 16); + } + + { + unsigned t = hxm * lym; /* 16 * 16 => 32 */ + { + unsigned 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_MBITS + 1 + 3))) | + ((rm << (SP_MBITS + 1 + 3)) != 0); + re++; + } else { + rm = (rm >> (32 - (SP_MBITS + 1 + 3 + 1))) | + ((rm << (SP_MBITS + 1 + 3 + 1)) != 0); + } + assert(rm & (SP_HIDDEN_BIT << 3)); + + SPNORMRET2(rs, re, rm, "mul", x, y); + } +} diff --git a/arch/mips/math-emu/sp_scalb.c b/arch/mips/math-emu/sp_scalb.c new file mode 100644 index 0000000..44ceb87 --- /dev/null +++ b/arch/mips/math-emu/sp_scalb.c @@ -0,0 +1,58 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_scalb(ieee754sp x, int n) +{ + COMPXSP; + + CLEARCX; + + EXPLODEXSP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + return ieee754sp_nanxcpt(x, "scalb", x, n); + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + case IEEE754_CLASS_ZERO: + return x; + case IEEE754_CLASS_DNORM: + SPDNORMX; + break; + case IEEE754_CLASS_NORM: + break; + } + SPNORMRET2(xs, xe + n, xm << 3, "scalb", x, n); +} + + +ieee754sp ieee754sp_ldexp(ieee754sp x, int n) +{ + return ieee754sp_scalb(x, n); +} diff --git a/arch/mips/math-emu/sp_simple.c b/arch/mips/math-emu/sp_simple.c new file mode 100644 index 0000000..c809830 --- /dev/null +++ b/arch/mips/math-emu/sp_simple.c @@ -0,0 +1,84 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +int ieee754sp_finite(ieee754sp x) +{ + return SPBEXP(x) != SP_EMAX + 1 + SP_EBIAS; +} + +ieee754sp ieee754sp_copysign(ieee754sp x, ieee754sp y) +{ + CLEARCX; + SPSIGN(x) = SPSIGN(y); + return x; +} + + +ieee754sp ieee754sp_neg(ieee754sp x) +{ + COMPXSP; + + EXPLODEXSP; + CLEARCX; + FLUSHXSP; + + if (xc == IEEE754_CLASS_SNAN) { + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "neg"); + } + + if (ieee754sp_isnan(x)) /* but not infinity */ + return ieee754sp_nanxcpt(x, "neg", x); + + /* quick fix up */ + SPSIGN(x) ^= 1; + return x; +} + + +ieee754sp ieee754sp_abs(ieee754sp x) +{ + COMPXSP; + + EXPLODEXSP; + CLEARCX; + FLUSHXSP; + + if (xc == IEEE754_CLASS_SNAN) { + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "abs"); + } + + if (ieee754sp_isnan(x)) /* but not infinity */ + return ieee754sp_nanxcpt(x, "abs", x); + + /* quick fix up */ + SPSIGN(x) = 0; + return x; +} diff --git a/arch/mips/math-emu/sp_sqrt.c b/arch/mips/math-emu/sp_sqrt.c new file mode 100644 index 0000000..8a934b9 --- /dev/null +++ b/arch/mips/math-emu/sp_sqrt.c @@ -0,0 +1,117 @@ +/* IEEE754 floating point arithmetic + * single precision square root + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_sqrt(ieee754sp x) +{ + int ix, s, q, m, t, i; + unsigned int r; + COMPXSP; + + /* take care of Inf and NaN */ + + EXPLODEXSP; + CLEARCX; + FLUSHXSP; + + /* x == INF or NAN? */ + switch (xc) { + case IEEE754_CLASS_QNAN: + /* sqrt(Nan) = Nan */ + return ieee754sp_nanxcpt(x, "sqrt"); + case IEEE754_CLASS_SNAN: + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt"); + case IEEE754_CLASS_ZERO: + /* sqrt(0) = 0 */ + return x; + case IEEE754_CLASS_INF: + if (xs) { + /* sqrt(-Inf) = Nan */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt"); + } + /* sqrt(+Inf) = Inf */ + return x; + case IEEE754_CLASS_DNORM: + case IEEE754_CLASS_NORM: + if (xs) { + /* sqrt(-x) = Nan */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt"); + } + 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; + q = s = 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) { + SETCX(IEEE754_INEXACT); + switch (ieee754_csr.rm) { + case IEEE754_RP: + q += 2; + break; + case IEEE754_RN: + q += (q & 1); + break; + } + } + ix = (q >> 1) + 0x3f000000; + ix += (m << 23); + x.bits = ix; + return x; +} diff --git a/arch/mips/math-emu/sp_sub.c b/arch/mips/math-emu/sp_sub.c new file mode 100644 index 0000000..dbb802c --- /dev/null +++ b/arch/mips/math-emu/sp_sub.c @@ -0,0 +1,184 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +ieee754sp ieee754sp_sub(ieee754sp x, ieee754sp y) +{ + COMPXSP; + COMPYSP; + + EXPLODEXSP; + EXPLODEYSP; + + CLEARCX; + + FLUSHXSP; + FLUSHYSP; + + switch (CLPAIR(xc, yc)) { + case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN): + case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN): + case CLPAIR(IEEE754_CLASS_SNAN, 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): + 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): + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_nanxcpt(ieee754sp_indef(), "sub", x, y); + + 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; + SETCX(IEEE754_INVALID_OPERATION); + return ieee754sp_xcpt(ieee754sp_indef(), "sub", x, y); + + 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 == + IEEE754_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): + SPDNORMX; + + 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 + */ + int s = xe - ye; + SPXSRSYn(s); + } else if (ye > xe) { + /* have to shift x fraction right to align + */ + int s = ye - xe; + SPXSRSXn(s); + } + assert(xe == ye); + assert(xe <= SP_EMAX); + + if (xs == ys) { + /* generate 28 bit result of adding two 27 bit numbers + */ + xm = xm + ym; + xe = xe; + xs = xs; + + if (xm >> (SP_MBITS + 1 + 3)) { /* carry out */ + SPXSRSX1(); /* shift preserving sticky */ + } + } else { + if (xm >= ym) { + xm = xm - ym; + xe = xe; + xs = xs; + } else { + xm = ym - xm; + xe = xe; + xs = ys; + } + if (xm == 0) { + if (ieee754_csr.rm == IEEE754_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_MBITS + 3)) == 0) { + xm <<= 1; + xe--; + } + } + SPNORMRET2(xs, xe, xm, "sub", x, y); +} diff --git a/arch/mips/math-emu/sp_tint.c b/arch/mips/math-emu/sp_tint.c new file mode 100644 index 0000000..1d73d2a --- /dev/null +++ b/arch/mips/math-emu/sp_tint.c @@ -0,0 +1,128 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include <linux/kernel.h> +#include "ieee754sp.h" + +int ieee754sp_tint(ieee754sp x) +{ + COMPXSP; + + CLEARCX; + + EXPLODEXSP; + FLUSHXSP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + SETCX(IEEE754_INVALID_OPERATION); + return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x); + 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 */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x); + } + /* oh gawd */ + if (xe > SP_MBITS) { + xm <<= xe - SP_MBITS; + } else { + u32 residue; + int round; + int sticky; + int odd; + + 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_MBITS; + round = (residue >> 31) != 0; + sticky = (residue << 1) != 0; + xm >>= SP_MBITS - xe; + } + odd = (xm & 0x1) != 0x0; + switch (ieee754_csr.rm) { + case IEEE754_RN: + if (round && (sticky || odd)) + xm++; + break; + case IEEE754_RZ: + break; + case IEEE754_RU: /* toward +Infinity */ + if ((round || sticky) && !xs) + xm++; + break; + case IEEE754_RD: /* toward -Infinity */ + if ((round || sticky) && xs) + xm++; + break; + } + if ((xm >> 31) != 0) { + /* This can happen after rounding */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x); + } + if (round || sticky) + SETCX(IEEE754_INEXACT); + } + if (xs) + return -xm; + else + return xm; +} + + +unsigned int ieee754sp_tuns(ieee754sp x) +{ + ieee754sp hb = ieee754sp_1e31(); + + /* what if x < 0 ?? */ + if (ieee754sp_lt(x, hb)) + return (unsigned) ieee754sp_tint(x); + + return (unsigned) ieee754sp_tint(ieee754sp_sub(x, hb)) | + ((unsigned) 1 << 31); +} diff --git a/arch/mips/math-emu/sp_tlong.c b/arch/mips/math-emu/sp_tlong.c new file mode 100644 index 0000000..4be21aa --- /dev/null +++ b/arch/mips/math-emu/sp_tlong.c @@ -0,0 +1,123 @@ +/* IEEE754 floating point arithmetic + * single precision + */ +/* + * MIPS floating point support + * Copyright (C) 1994-2000 Algorithmics Ltd. + * http://www.algor.co.uk + * + * ######################################################################## + * + * This program is free software; you can distribute it and/or modify it + * under the terms of the GNU General Public License (Version 2) as + * published by the Free Software Foundation. + * + * This program is distributed in the hope it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License along + * with this program; if not, write to the Free Software Foundation, Inc., + * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. + * + * ######################################################################## + */ + + +#include "ieee754sp.h" + +s64 ieee754sp_tlong(ieee754sp x) +{ + COMPXDP; /* <-- need 64-bit mantissa tmp */ + + CLEARCX; + + EXPLODEXSP; + FLUSHXSP; + + switch (xc) { + case IEEE754_CLASS_SNAN: + case IEEE754_CLASS_QNAN: + case IEEE754_CLASS_INF: + SETCX(IEEE754_INVALID_OPERATION); + return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x); + 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 */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x); + } + /* oh gawd */ + if (xe > SP_MBITS) { + xm <<= xe - SP_MBITS; + } else if (xe < SP_MBITS) { + u32 residue; + int round; + int sticky; + int odd; + + if (xe < -1) { + residue = xm; + round = 0; + sticky = residue != 0; + xm = 0; + } + else { + residue = xm << (32 - SP_MBITS + xe); + round = (residue >> 31) != 0; + sticky = (residue << 1) != 0; + xm >>= SP_MBITS - xe; + } + odd = (xm & 0x1) != 0x0; + switch (ieee754_csr.rm) { + case IEEE754_RN: + if (round && (sticky || odd)) + xm++; + break; + case IEEE754_RZ: + break; + case IEEE754_RU: /* toward +Infinity */ + if ((round || sticky) && !xs) + xm++; + break; + case IEEE754_RD: /* toward -Infinity */ + if ((round || sticky) && xs) + xm++; + break; + } + if ((xm >> 63) != 0) { + /* This can happen after rounding */ + SETCX(IEEE754_INVALID_OPERATION); + return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x); + } + if (round || sticky) + SETCX(IEEE754_INEXACT); + } + if (xs) + return -xm; + else + return xm; +} + + +u64 ieee754sp_tulong(ieee754sp x) +{ + ieee754sp hb = ieee754sp_1e63(); + + /* what if x < 0 ?? */ + if (ieee754sp_lt(x, hb)) + return (u64) ieee754sp_tlong(x); + + return (u64) ieee754sp_tlong(ieee754sp_sub(x, hb)) | + (1ULL << 63); +} |
