x86/kprobes: Fix crash when probing CS CALL instructionsgrafted

When using eBPF to probe CS CALL instructions within a function,
a crash can be triggered.

The eBPF tool probes offset 257 of the __hrtimer_run_queues()
function:

<__hrtimer_run_queues+249>:  nopl   0x0(%rax,%rax,1)
<__hrtimer_run_queues+254>:  mov    %r14,%rdi
<__hrtimer_run_queues+257>:  cs call <__x86_indirect_thunk_r12>
<__hrtimer_run_queues+263>:  mov    %eax,%r12d
<__hrtimer_run_queues+266>:  xchg   %ax,%ax
<__hrtimer_run_queues+268>:  mov    %r13,%rdi

Which triggers this crash:

  BUG: unable to handle page fault for address: 00000000000f41c9
  #PF: supervisor write access in kernel mode
  #PF: error_code(0x0002) - not-present page
  PGD 0 P4D 0
  Oops: 0002 [#1] SMP NOPTI
  CPU: 1 PID: 0 Comm: swapper/1 Kdump: loaded Tainted: P
  RIP: 0010:__hrtimer_run_queues+0x106/0x230

Note that __hrtimer_run_queues+0x106 is __hrtimer_run_queues+262, which is
at the 6th byte of the above CS CALL instruction. Since the CS CALL
instruction occupies 6 bytes, the exception occurred in the middle of that
call instruction.

The root cause is that when using eBPF tools to probe in the middle of a
function, a kprobe with INT3 is used as the underlying implementation.

During single-step emulation of the original CALL instruction,
int3_emulate_call() assumes that the probed CALL instruction is 5 bytes
long. However, the actual CS-prefixed CALL instruction occupies 6 bytes,
so it constructs an incorrect exception return address. When the CPU
returns from the kprobe handler, the next instruction to be executed is at
the address of the last byte of that CS CALL instruction. Coincidentally,
starting from that address, the CPU fetches and decodes a completely
different instruction, which ultimately triggers a kernel crash.

Fix the issue by using the actual instruction length obtained from
the instruction decoder when constructing the exception return
address, rather than relying on the hardcoded CALL_INSN_SIZE macro.

[ mingo: Refined the changelog ]

Fixes: 6256e668b7af ("x86/kprobes: Use int3 instead of debug trap for single-step")
Suggested-by: Masami Hiramatsu (Google) <mhiramat@kernel.org>
Signed-off-by: Jinke Han <jinkehan@didiglobal.com>
Signed-off-by: Ingo Molnar <mingo@kernel.org>
Reviewed-by: Masami Hiramatsu (Google) <mhiramat@kernel.org>
Acked-by: Yafang Shao <laoar.shao@gmail.com>
Acked-by: Borislav Petkov <bp@alien8.de>
Cc: Peter Zijlstra <peterz@infradead.org>
Link: https://patch.msgid.link/20260908073742.GA10517@didi-ThinkCentre-M920t-N000

24 files changed, 7775 insertions, 0 deletions

diff --git a/arch/arm/nwfpe/ARM-gcc.h b/arch/arm/nwfpe/ARM-gcc.h
new file mode 100644
index 000000000..74389f5b0
--- /dev/null
+++ b/arch/arm/nwfpe/ARM-gcc.h
@@ -0,0 +1,121 @@
+/* SPDX-License-Identifier: GPL-2.0 */
+/*
+-------------------------------------------------------------------------------
+The macro `BITS64' can be defined to indicate that 64-bit integer types are
+supported by the compiler.
+-------------------------------------------------------------------------------
+*/
+#define BITS64
+
+/*
+-------------------------------------------------------------------------------
+Each of the following `typedef's defines the most convenient type that holds
+integers of at least as many bits as specified.  For example, `uint8' should
+be the most convenient type that can hold unsigned integers of as many as
+8 bits.  The `flag' type must be able to hold either a 0 or 1.  For most
+implementations of C, `flag', `uint8', and `int8' should all be `typedef'ed
+to the same as `int'.
+-------------------------------------------------------------------------------
+*/
+typedef char flag;
+typedef unsigned char uint8;
+typedef signed char int8;
+typedef int uint16;
+typedef int int16;
+typedef unsigned int uint32;
+typedef signed int int32;
+#ifdef BITS64
+typedef unsigned long long int bits64;
+typedef signed long long int sbits64;
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Each of the following `typedef's defines a type that holds integers
+of _exactly_ the number of bits specified.  For instance, for most
+implementation of C, `bits16' and `sbits16' should be `typedef'ed to
+`unsigned short int' and `signed short int' (or `short int'), respectively.
+-------------------------------------------------------------------------------
+*/
+typedef unsigned char bits8;
+typedef signed char sbits8;
+typedef unsigned short int bits16;
+typedef signed short int sbits16;
+typedef unsigned int bits32;
+typedef signed int sbits32;
+#ifdef BITS64
+typedef unsigned long long int uint64;
+typedef signed long long int int64;
+#endif
+
+#ifdef BITS64
+/*
+-------------------------------------------------------------------------------
+The `LIT64' macro takes as its argument a textual integer literal and if
+necessary ``marks'' the literal as having a 64-bit integer type.  For
+example, the Gnu C Compiler (`gcc') requires that 64-bit literals be
+appended with the letters `LL' standing for `long long', which is `gcc's
+name for the 64-bit integer type.  Some compilers may allow `LIT64' to be
+defined as the identity macro:  `#define LIT64( a ) a'.
+-------------------------------------------------------------------------------
+*/
+#define LIT64( a ) a##LL
+#endif
+
+/*
+-------------------------------------------------------------------------------
+The macro `INLINE' can be used before functions that should be inlined.  If
+a compiler does not support explicit inlining, this macro should be defined
+to be `static'.
+-------------------------------------------------------------------------------
+*/
+#define INLINE static inline
+
+
+/* For use as a GCC soft-float library we need some special function names. */
+
+#ifdef __LIBFLOAT__
+
+/* Some 32-bit ops can be mapped straight across by just changing the name. */
+#define float32_add			__addsf3
+#define float32_sub			__subsf3
+#define float32_mul			__mulsf3
+#define float32_div			__divsf3
+#define int32_to_float32		__floatsisf
+#define float32_to_int32_round_to_zero	__fixsfsi
+#define float32_to_uint32_round_to_zero	__fixunssfsi
+
+/* These ones go through the glue code.  To avoid namespace pollution
+   we rename the internal functions too.  */
+#define float32_eq			___float32_eq
+#define float32_le			___float32_le
+#define float32_lt			___float32_lt
+
+/* All the 64-bit ops have to go through the glue, so we pull the same
+   trick.  */
+#define float64_add			___float64_add
+#define float64_sub			___float64_sub
+#define float64_mul			___float64_mul
+#define float64_div			___float64_div
+#define int32_to_float64		___int32_to_float64
+#define float64_to_int32_round_to_zero	___float64_to_int32_round_to_zero
+#define float64_to_uint32_round_to_zero	___float64_to_uint32_round_to_zero
+#define float64_to_float32		___float64_to_float32
+#define float32_to_float64		___float32_to_float64
+#define float64_eq			___float64_eq
+#define float64_le			___float64_le
+#define float64_lt			___float64_lt
+
+#if 0
+#define float64_add			__adddf3
+#define float64_sub			__subdf3
+#define float64_mul			__muldf3
+#define float64_div			__divdf3
+#define int32_to_float64		__floatsidf
+#define float64_to_int32_round_to_zero	__fixdfsi
+#define float64_to_uint32_round_to_zero	__fixunsdfsi
+#define float64_to_float32		__truncdfsf2
+#define float32_to_float64		__extendsfdf2
+#endif
+
+#endif
diff --git a/arch/arm/nwfpe/ChangeLog b/arch/arm/nwfpe/ChangeLog
new file mode 100644
index 000000000..fa8028b1e
--- /dev/null
+++ b/arch/arm/nwfpe/ChangeLog
@@ -0,0 +1,91 @@
+2003-03-22  Ralph Siemsen <ralphs@netwinder.org>
+	* Reformat all but softfloat files to get a consistent coding style.
+	  Used "indent -kr -i8 -ts8 -sob -l132 -ss" and a few manual fixups.
+	* Removed dead code and fixed function protypes to match definitions.
+	* Consolidated use of (opcode && MASK_ARITHMETIC_OPCODE) >> 20.
+	* Make 80-bit precision a compile-time option. (1%)
+	* Only initialize FPE state once in repeat-FP situations. (6%)
+
+2002-01-19  Russell King <rmk@arm.linux.org.uk>
+
+	* fpa11.h - Add documentation
+	          - remove userRegisters pointer from this structure.
+	          - add new method to obtain integer register values.
+	* softfloat.c - Remove float128
+	* softfloat.h - Remove float128
+	* softfloat-specialize - Remove float128
+
+	* The FPA11 structure is not a kernel-specific data structure.
+	  It is used by users of ptrace to examine the values of the
+	  floating point registers.  Therefore, any changes to the
+	  FPA11 structure (size or position of elements contained
+	  within) have to be well thought out.
+
+	* Since 128-bit float requires the FPA11 structure to change
+	  size, it has been removed.  128-bit float is currently unused,
+	  and needs various things to be re-worked so that we won't
+	  overflow the available space in the task structure.
+
+	* The changes are designed to break any patch that goes on top
+	  of this code, so that the authors properly review their changes.
+
+1999-08-19  Scott Bambrough  <scottb@netwinder.org>
+
+	* fpmodule.c - Changed version number to 0.95
+	* fpa11.h - modified FPA11, FPREG structures
+	* fpa11.c - Changes due to FPA11, FPREG structure alterations.
+	* fpa11_cpdo.c - Changes due to FPA11, FPREG structure alterations.
+	* fpa11_cpdt.c - Changes due to FPA11, FPREG structure alterations.
+	* fpa11_cprt.c - Changes due to FPA11, FPREG structure alterations.
+	* single_cpdo.c - Changes due to FPA11, FPREG structure alterations.
+	* double_cpdo.c - Changes due to FPA11, FPREG structure alterations.
+	* extended_cpdo.c - Changes due to FPA11, FPREG structure alterations.
+
+	* I discovered several bugs.  First and worst is that the kernel
+	  passes in a pointer to the FPE's state area.	This is defined
+	  as a struct user_fp (see user.h).  This pointer was cast to a
+	  FPA11*.  Unfortunately FPA11 and user_fp are of different sizes;
+	  user_fp is smaller.  This meant that the FPE scribbled on things
+	  below its area, which is bad, as the area is in the thread_struct
+	  embedded in the process task structure.  Thus we were scribbling
+	  over one of the most important structures in the entire OS.
+
+	* user_fp and FPA11 have now been harmonized.  Most of the changes
+	  in the above code were dereferencing problems due to moving the
+	  register type out of FPREG, and getting rid of the union variable
+	  fpvalue.
+
+	* Second I noticed resetFPA11 was not always being called for a
+	  task.  This should happen on the first floating point exception
+	  that occurs.	It is controlled by init_flag in FPA11.  The
+	  comment in the code beside init_flag state the kernel guarantees
+	  this to be zero.  Not so.  I found that the kernel recycles task
+	  structures, and that recycled ones may not have init_flag zeroed.
+	  I couldn't even find anything that guarantees it is zeroed when
+	  when the task structure is initially allocated.  In any case
+	  I now initialize the entire FPE state in the thread structure to
+	  zero when allocated and recycled.  See alloc_task_struct() and
+	  flush_thread() in arch/arm/process.c.  The change to
+	  alloc_task_struct() may not be necessary, but I left it in for
+	  completeness (better safe than sorry).
+
+1998-11-23  Scott Bambrough  <scottb@netwinder.org>
+
+	* README.FPE - fix typo in description of lfm/sfm instructions
+	* NOTES - Added file to describe known bugs/problems
+	* fpmodule.c - Changed version number to 0.94
+
+1998-11-20  Scott Bambrough  <scottb@netwinder.org>
+
+	* README.FPE - fix description of URD, NRM instructions
+	* TODO - remove URD, NRM instructions from TODO list
+	* single_cpdo.c - implement URD, NRM
+	* double_cpdo.c - implement URD, NRM
+	* extended_cpdo.c - implement URD, NRM
+
+1998-11-19  Scott Bambrough  <scottb@netwinder.org>
+
+	* ChangeLog - Added this file to track changes made.
+	* fpa11.c - added code to initialize register types to typeNone
+	* fpa11_cpdt.c - fixed bug in storeExtended (typeExtended changed to
+	  typeDouble in switch statement)
diff --git a/arch/arm/nwfpe/Makefile b/arch/arm/nwfpe/Makefile
new file mode 100644
index 000000000..2aec85ab1
--- /dev/null
+++ b/arch/arm/nwfpe/Makefile
@@ -0,0 +1,19 @@
+# SPDX-License-Identifier: GPL-2.0
+#
+# Copyright (C) 1998, 1999, 2001 Philip Blundell
+#
+
+obj-$(CONFIG_FPE_NWFPE)		+= nwfpe.o
+
+nwfpe-y				+= fpa11.o fpa11_cpdo.o fpa11_cpdt.o \
+				   fpa11_cprt.o fpmodule.o fpopcode.o \
+				   softfloat.o single_cpdo.o double_cpdo.o \
+				   entry.o
+
+nwfpe-$(CONFIG_FPE_NWFPE_XP)	+= extended_cpdo.o
+
+# Try really hard to avoid generating calls to __aeabi_uldivmod() from
+# float64_rem() due to loop elision.
+ifdef CONFIG_CC_IS_CLANG
+CFLAGS_softfloat.o	+= -mllvm -replexitval=never
+endif
diff --git a/arch/arm/nwfpe/double_cpdo.c b/arch/arm/nwfpe/double_cpdo.c
new file mode 100644
index 000000000..344bb7c46
--- /dev/null
+++ b/arch/arm/nwfpe/double_cpdo.c
@@ -0,0 +1,155 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+#include "softfloat.h"
+#include "fpopcode.h"
+
+union float64_components {
+	float64 f64;
+	unsigned int i[2];
+};
+
+float64 float64_exp(float64 Fm);
+float64 float64_ln(float64 Fm);
+float64 float64_sin(float64 rFm);
+float64 float64_cos(float64 rFm);
+float64 float64_arcsin(float64 rFm);
+float64 float64_arctan(float64 rFm);
+float64 float64_log(float64 rFm);
+float64 float64_tan(float64 rFm);
+float64 float64_arccos(float64 rFm);
+float64 float64_pow(float64 rFn, float64 rFm);
+float64 float64_pol(float64 rFn, float64 rFm);
+
+static float64 float64_rsf(struct roundingData *roundData, float64 rFn, float64 rFm)
+{
+	return float64_sub(roundData, rFm, rFn);
+}
+
+static float64 float64_rdv(struct roundingData *roundData, float64 rFn, float64 rFm)
+{
+	return float64_div(roundData, rFm, rFn);
+}
+
+static float64 (*const dyadic_double[16])(struct roundingData*, float64 rFn, float64 rFm) = {
+	[ADF_CODE >> 20] = float64_add,
+	[MUF_CODE >> 20] = float64_mul,
+	[SUF_CODE >> 20] = float64_sub,
+	[RSF_CODE >> 20] = float64_rsf,
+	[DVF_CODE >> 20] = float64_div,
+	[RDF_CODE >> 20] = float64_rdv,
+	[RMF_CODE >> 20] = float64_rem,
+
+	/* strictly, these opcodes should not be implemented */
+	[FML_CODE >> 20] = float64_mul,
+	[FDV_CODE >> 20] = float64_div,
+	[FRD_CODE >> 20] = float64_rdv,
+};
+
+static float64 float64_mvf(struct roundingData *roundData,float64 rFm)
+{
+	return rFm;
+}
+
+static float64 float64_mnf(struct roundingData *roundData,float64 rFm)
+{
+	union float64_components u;
+
+	u.f64 = rFm;
+#ifdef __ARMEB__
+	u.i[0] ^= 0x80000000;
+#else
+	u.i[1] ^= 0x80000000;
+#endif
+
+	return u.f64;
+}
+
+static float64 float64_abs(struct roundingData *roundData,float64 rFm)
+{
+	union float64_components u;
+
+	u.f64 = rFm;
+#ifdef __ARMEB__
+	u.i[0] &= 0x7fffffff;
+#else
+	u.i[1] &= 0x7fffffff;
+#endif
+
+	return u.f64;
+}
+
+static float64 (*const monadic_double[16])(struct roundingData *, float64 rFm) = {
+	[MVF_CODE >> 20] = float64_mvf,
+	[MNF_CODE >> 20] = float64_mnf,
+	[ABS_CODE >> 20] = float64_abs,
+	[RND_CODE >> 20] = float64_round_to_int,
+	[URD_CODE >> 20] = float64_round_to_int,
+	[SQT_CODE >> 20] = float64_sqrt,
+	[NRM_CODE >> 20] = float64_mvf,
+};
+
+unsigned int DoubleCPDO(struct roundingData *roundData, const unsigned int opcode, FPREG * rFd)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	float64 rFm;
+	unsigned int Fm, opc_mask_shift;
+
+	Fm = getFm(opcode);
+	if (CONSTANT_FM(opcode)) {
+		rFm = getDoubleConstant(Fm);
+	} else {
+		switch (fpa11->fType[Fm]) {
+		case typeSingle:
+			rFm = float32_to_float64(fpa11->fpreg[Fm].fSingle);
+			break;
+
+		case typeDouble:
+			rFm = fpa11->fpreg[Fm].fDouble;
+			break;
+
+		default:
+			return 0;
+		}
+	}
+
+	opc_mask_shift = (opcode & MASK_ARITHMETIC_OPCODE) >> 20;
+	if (!MONADIC_INSTRUCTION(opcode)) {
+		unsigned int Fn = getFn(opcode);
+		float64 rFn;
+
+		switch (fpa11->fType[Fn]) {
+		case typeSingle:
+			rFn = float32_to_float64(fpa11->fpreg[Fn].fSingle);
+			break;
+
+		case typeDouble:
+			rFn = fpa11->fpreg[Fn].fDouble;
+			break;
+
+		default:
+			return 0;
+		}
+
+		if (dyadic_double[opc_mask_shift]) {
+			rFd->fDouble = dyadic_double[opc_mask_shift](roundData, rFn, rFm);
+		} else {
+			return 0;
+		}
+	} else {
+		if (monadic_double[opc_mask_shift]) {
+			rFd->fDouble = monadic_double[opc_mask_shift](roundData, rFm);
+		} else {
+			return 0;
+		}
+	}
+
+	return 1;
+}
diff --git a/arch/arm/nwfpe/entry.S b/arch/arm/nwfpe/entry.S
new file mode 100644
index 000000000..354d297a1
--- /dev/null
+++ b/arch/arm/nwfpe/entry.S
@@ -0,0 +1,190 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998
+    (c) 1998, 1999 Philip Blundell
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+#include <linux/linkage.h>
+#include <asm/assembler.h>
+#include <asm/opcodes.h>
+
+/* This is the kernel's entry point into the floating point emulator.
+It is called from the kernel with code similar to this:
+
+	sub	r4, r5, #4
+	ldrt	r0, [r4]			@ r0  = instruction
+	adrsvc	al, r9, ret_from_exception	@ r9  = normal FP return
+	adrsvc	al, lr, fpundefinstr		@ lr  = undefined instr return
+
+	get_current_task r10
+	mov	r8, #1
+	strb	r8, [r10, #TSK_USED_MATH]	@ set current->used_math
+	add	r10, r10, #TSS_FPESAVE		@ r10 = workspace
+	ldr	r4, .LC2
+	ldr	pc, [r4]			@ Call FP emulator entry point
+
+The kernel expects the emulator to return via one of two possible
+points of return it passes to the emulator.  The emulator, if
+successful in its emulation, jumps to ret_from_exception (passed in
+r9) and the kernel takes care of returning control from the trap to
+the user code.  If the emulator is unable to emulate the instruction,
+it returns via _fpundefinstr (passed via lr) and the kernel halts the
+user program with a core dump.
+
+On entry to the emulator r10 points to an area of private FP workspace
+reserved in the thread structure for this process.  This is where the
+emulator saves its registers across calls.  The first word of this area
+is used as a flag to detect the first time a process uses floating point,
+so that the emulator startup cost can be avoided for tasks that don't
+want it.
+
+This routine does three things:
+
+1) The kernel has created a struct pt_regs on the stack and saved the
+user registers into it.  See /usr/include/asm/proc/ptrace.h for details.
+
+2) It calls EmulateAll to emulate a floating point instruction.
+EmulateAll returns 1 if the emulation was successful, or 0 if not.
+
+3) If an instruction has been emulated successfully, it looks ahead at
+the next instruction.  If it is a floating point instruction, it
+executes the instruction, without returning to user space.  In this
+way it repeatedly looks ahead and executes floating point instructions
+until it encounters a non floating point instruction, at which time it
+returns via _fpreturn.
+
+This is done to reduce the effect of the trap overhead on each
+floating point instructions.  GCC attempts to group floating point
+instructions to allow the emulator to spread the cost of the trap over
+several floating point instructions.  */
+
+#include <asm/asm-offsets.h>
+
+	.globl	nwfpe_enter
+nwfpe_enter:
+	mov	r4, lr			@ save the failure-return addresses
+	mov	sl, sp			@ we access the registers via 'sl'
+
+	ldr	r5, [sp, #S_PC]		@ get contents of PC;
+	mov	r6, r0			@ save the opcode
+emulate:
+	ldr	r1, [sp, #S_PSR]	@ fetch the PSR
+	bl	arm_check_condition	@ check the condition
+	cmp	r0, #ARM_OPCODE_CONDTEST_PASS	@ condition passed?
+
+	@ if condition code failed to match, next insn
+	bne	next			@ get the next instruction;
+
+	mov	r0, r6			@ prepare for EmulateAll()
+	bl	EmulateAll		@ emulate the instruction
+	cmp	r0, #0			@ was emulation successful
+	reteq	r4			@ no, return failure
+
+next:
+	uaccess_enable r3
+.Lx1:	ldrt	r6, [r5], #4		@ get the next instruction and
+					@ increment PC
+	uaccess_disable r3
+	and	r2, r6, #0x0F000000	@ test for FP insns
+	teq	r2, #0x0C000000
+	teqne	r2, #0x0D000000
+	teqne	r2, #0x0E000000
+	retne	r9			@ return ok if not a fp insn
+
+	str	r5, [sp, #S_PC]		@ update PC copy in regs
+
+	mov	r0, r6			@ save a copy
+	b	emulate			@ check condition and emulate
+
+	@ We need to be prepared for the instructions at .Lx1 and .Lx2 
+	@ to fault.  Emit the appropriate exception gunk to fix things up.
+	@ ??? For some reason, faults can happen at .Lx2 even with a
+	@ plain LDR instruction.  Weird, but it seems harmless.
+	.pushsection .text.fixup,"ax"
+	.align	2
+.Lrep:	str     r4, [sp, #S_PC]		@ retry current instruction
+.Lfix:	ret	r9			@ let the user eat segfaults
+	.popsection
+
+	.pushsection __ex_table,"a"
+	.align	3
+	.long	.Lx1, .Lfix
+	.popsection
+
+	@
+	@ Check whether the instruction is a co-processor instruction.
+	@ If yes, we need to call the relevant co-processor handler.
+	@ Only FPE instructions are dispatched here, everything else
+	@ is handled by undef hooks.
+	@
+	@ Emulators may wish to make use of the following registers:
+	@  r4  = PC value to resume execution after successful emulation
+	@  r9  = normal "successful" return address
+	@  lr  = unrecognised instruction return address
+	@ IRQs enabled, FIQs enabled.
+	@
+ENTRY(call_fpe)
+	mov	r2, r4
+	sub	r4, r4, #4			@ ARM instruction at user PC - 4
+USERL(	.Lrep,	ldrt r0, [r4])			@ load opcode from user space
+ARM_BE8(rev	r0, r0)				@ little endian instruction
+
+	uaccess_disable ip
+
+	get_thread_info r10			@ get current thread
+	tst	r0, #0x08000000			@ only CDP/CPRT/LDC/STC have bit 27
+	reteq	lr
+	and	r8, r0, #0x00000f00		@ mask out CP number
+#ifdef CONFIG_IWMMXT
+	@ Test if we need to give access to iWMMXt coprocessors
+	ldr	r5, [r10, #TI_FLAGS]
+	rsbs	r7, r8, #(1 << 8)		@ CP 0 or 1 only
+	movscs	r7, r5, lsr #(TIF_USING_IWMMXT + 1)
+	movcs	r0, sp				@ pass struct pt_regs
+	bcs	iwmmxt_task_enable
+#endif
+	add	pc, pc, r8, lsr #6
+	nop
+
+	ret	lr				@ CP#0
+	b	do_fpe				@ CP#1 (FPE)
+	b	do_fpe				@ CP#2 (FPE)
+	ret	lr				@ CP#3
+	ret	lr				@ CP#4
+	ret	lr				@ CP#5
+	ret	lr				@ CP#6
+	ret	lr				@ CP#7
+	ret	lr				@ CP#8
+	ret	lr				@ CP#9
+	ret	lr				@ CP#10 (VFP)
+	ret	lr				@ CP#11 (VFP)
+	ret	lr				@ CP#12
+	ret	lr				@ CP#13
+	ret	lr				@ CP#14 (Debug)
+	ret	lr				@ CP#15 (Control)
+
+do_fpe:
+	add	r10, r10, #TI_FPSTATE		@ r10 = workspace
+	ldr_va	pc, fp_enter, tmp=r4		@ Call FP module USR entry point
+
+	@
+	@ The FP module is called with these registers set:
+	@  r0  = instruction
+	@  r2  = PC+4
+	@  r9  = normal "successful" return address
+	@  r10 = FP workspace
+	@  lr  = unrecognised FP instruction return address
+	@
+
+	.pushsection .data
+	.align	2
+ENTRY(fp_enter)
+	.word	no_fp
+	.popsection
+
+no_fp:
+	ret	lr
+ENDPROC(no_fp)
diff --git a/arch/arm/nwfpe/extended_cpdo.c b/arch/arm/nwfpe/extended_cpdo.c
new file mode 100644
index 000000000..529b047ba
--- /dev/null
+++ b/arch/arm/nwfpe/extended_cpdo.c
@@ -0,0 +1,142 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+#include "softfloat.h"
+#include "fpopcode.h"
+
+floatx80 floatx80_exp(floatx80 Fm);
+floatx80 floatx80_ln(floatx80 Fm);
+floatx80 floatx80_sin(floatx80 rFm);
+floatx80 floatx80_cos(floatx80 rFm);
+floatx80 floatx80_arcsin(floatx80 rFm);
+floatx80 floatx80_arctan(floatx80 rFm);
+floatx80 floatx80_log(floatx80 rFm);
+floatx80 floatx80_tan(floatx80 rFm);
+floatx80 floatx80_arccos(floatx80 rFm);
+floatx80 floatx80_pow(floatx80 rFn, floatx80 rFm);
+floatx80 floatx80_pol(floatx80 rFn, floatx80 rFm);
+
+static floatx80 floatx80_rsf(struct roundingData *roundData, floatx80 rFn, floatx80 rFm)
+{
+	return floatx80_sub(roundData, rFm, rFn);
+}
+
+static floatx80 floatx80_rdv(struct roundingData *roundData, floatx80 rFn, floatx80 rFm)
+{
+	return floatx80_div(roundData, rFm, rFn);
+}
+
+static floatx80 (*const dyadic_extended[16])(struct roundingData*, floatx80 rFn, floatx80 rFm) = {
+	[ADF_CODE >> 20] = floatx80_add,
+	[MUF_CODE >> 20] = floatx80_mul,
+	[SUF_CODE >> 20] = floatx80_sub,
+	[RSF_CODE >> 20] = floatx80_rsf,
+	[DVF_CODE >> 20] = floatx80_div,
+	[RDF_CODE >> 20] = floatx80_rdv,
+	[RMF_CODE >> 20] = floatx80_rem,
+
+	/* strictly, these opcodes should not be implemented */
+	[FML_CODE >> 20] = floatx80_mul,
+	[FDV_CODE >> 20] = floatx80_div,
+	[FRD_CODE >> 20] = floatx80_rdv,
+};
+
+static floatx80 floatx80_mvf(struct roundingData *roundData, floatx80 rFm)
+{
+	return rFm;
+}
+
+static floatx80 floatx80_mnf(struct roundingData *roundData, floatx80 rFm)
+{
+	rFm.high ^= 0x8000;
+	return rFm;
+}
+
+static floatx80 floatx80_abs(struct roundingData *roundData, floatx80 rFm)
+{
+	rFm.high &= 0x7fff;
+	return rFm;
+}
+
+static floatx80 (*const monadic_extended[16])(struct roundingData*, floatx80 rFm) = {
+	[MVF_CODE >> 20] = floatx80_mvf,
+	[MNF_CODE >> 20] = floatx80_mnf,
+	[ABS_CODE >> 20] = floatx80_abs,
+	[RND_CODE >> 20] = floatx80_round_to_int,
+	[URD_CODE >> 20] = floatx80_round_to_int,
+	[SQT_CODE >> 20] = floatx80_sqrt,
+	[NRM_CODE >> 20] = floatx80_mvf,
+};
+
+unsigned int ExtendedCPDO(struct roundingData *roundData, const unsigned int opcode, FPREG * rFd)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	floatx80 rFm;
+	unsigned int Fm, opc_mask_shift;
+
+	Fm = getFm(opcode);
+	if (CONSTANT_FM(opcode)) {
+		rFm = getExtendedConstant(Fm);
+	} else {
+		switch (fpa11->fType[Fm]) {
+		case typeSingle:
+			rFm = float32_to_floatx80(fpa11->fpreg[Fm].fSingle);
+			break;
+
+		case typeDouble:
+			rFm = float64_to_floatx80(fpa11->fpreg[Fm].fDouble);
+			break;
+
+		case typeExtended:
+			rFm = fpa11->fpreg[Fm].fExtended;
+			break;
+
+		default:
+			return 0;
+		}
+	}
+
+	opc_mask_shift = (opcode & MASK_ARITHMETIC_OPCODE) >> 20;
+	if (!MONADIC_INSTRUCTION(opcode)) {
+		unsigned int Fn = getFn(opcode);
+		floatx80 rFn;
+
+		switch (fpa11->fType[Fn]) {
+		case typeSingle:
+			rFn = float32_to_floatx80(fpa11->fpreg[Fn].fSingle);
+			break;
+
+		case typeDouble:
+			rFn = float64_to_floatx80(fpa11->fpreg[Fn].fDouble);
+			break;
+
+		case typeExtended:
+			rFn = fpa11->fpreg[Fn].fExtended;
+			break;
+
+		default:
+			return 0;
+		}
+
+		if (dyadic_extended[opc_mask_shift]) {
+			rFd->fExtended = dyadic_extended[opc_mask_shift](roundData, rFn, rFm);
+		} else {
+			return 0;
+		}
+	} else {
+		if (monadic_extended[opc_mask_shift]) {
+			rFd->fExtended = monadic_extended[opc_mask_shift](roundData, rFm);
+		} else {
+			return 0;
+		}
+	}
+
+	return 1;
+}
diff --git a/arch/arm/nwfpe/fpa11.c b/arch/arm/nwfpe/fpa11.c
new file mode 100644
index 000000000..2379e2cc7
--- /dev/null
+++ b/arch/arm/nwfpe/fpa11.c
@@ -0,0 +1,116 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+    (c) Philip Blundell, 2001
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+#include "fpopcode.h"
+
+#include "fpmodule.h"
+#include "fpmodule.inl"
+
+#include <linux/compiler.h>
+#include <linux/string.h>
+
+/* Reset the FPA11 chip.  Called to initialize and reset the emulator. */
+static void resetFPA11(void)
+{
+	int i;
+	FPA11 *fpa11 = GET_FPA11();
+
+	/* initialize the register type array */
+	for (i = 0; i <= 7; i++) {
+		fpa11->fType[i] = typeNone;
+	}
+
+	/* FPSR: set system id to FP_EMULATOR, set AC, clear all other bits */
+	fpa11->fpsr = FP_EMULATOR | BIT_AC;
+}
+
+int8 SetRoundingMode(const unsigned int opcode)
+{
+	switch (opcode & MASK_ROUNDING_MODE) {
+	default:
+	case ROUND_TO_NEAREST:
+		return float_round_nearest_even;
+
+	case ROUND_TO_PLUS_INFINITY:
+		return float_round_up;
+
+	case ROUND_TO_MINUS_INFINITY:
+		return float_round_down;
+
+	case ROUND_TO_ZERO:
+		return float_round_to_zero;
+	}
+}
+
+int8 SetRoundingPrecision(const unsigned int opcode)
+{
+#ifdef CONFIG_FPE_NWFPE_XP
+	switch (opcode & MASK_ROUNDING_PRECISION) {
+	case ROUND_SINGLE:
+		return 32;
+
+	case ROUND_DOUBLE:
+		return 64;
+
+	case ROUND_EXTENDED:
+		return 80;
+
+	default:
+		return 80;
+	}
+#endif
+	return 80;
+}
+
+void nwfpe_init_fpa(union fp_state *fp)
+{
+	FPA11 *fpa11 = (FPA11 *)fp;
+#ifdef NWFPE_DEBUG
+	printk("NWFPE: setting up state.\n");
+#endif
+ 	memset(fpa11, 0, sizeof(FPA11));
+	resetFPA11();
+	fpa11->initflag = 1;
+}
+
+/* Emulate the instruction in the opcode. */
+unsigned int EmulateAll(unsigned int opcode)
+{
+	unsigned int code;
+
+#ifdef NWFPE_DEBUG
+	printk("NWFPE: emulating opcode %08x\n", opcode);
+#endif
+	code = opcode & 0x00000f00;
+	if (code == 0x00000100 || code == 0x00000200) {
+		/* For coprocessor 1 or 2 (FPA11) */
+		code = opcode & 0x0e000000;
+		if (code == 0x0e000000) {
+			if (opcode & 0x00000010) {
+				/* Emulate conversion opcodes. */
+				/* Emulate register transfer opcodes. */
+				/* Emulate comparison opcodes. */
+				return EmulateCPRT(opcode);
+			} else {
+				/* Emulate monadic arithmetic opcodes. */
+				/* Emulate dyadic arithmetic opcodes. */
+				return EmulateCPDO(opcode);
+			}
+		} else if (code == 0x0c000000) {
+			/* Emulate load/store opcodes. */
+			/* Emulate load/store multiple opcodes. */
+			return EmulateCPDT(opcode);
+		}
+	}
+
+	/* Invalid instruction detected.  Return FALSE. */
+	return 0;
+}
diff --git a/arch/arm/nwfpe/fpa11.h b/arch/arm/nwfpe/fpa11.h
new file mode 100644
index 000000000..d51195d41
--- /dev/null
+++ b/arch/arm/nwfpe/fpa11.h
@@ -0,0 +1,109 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.com, 1998-1999
+    
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#ifndef __FPA11_H__
+#define __FPA11_H__
+
+#define GET_FPA11() ((FPA11 *)(&current_thread_info()->fpstate))
+
+/*
+ * The processes registers are always at the very top of the 8K
+ * stack+task struct.  Use the same method as 'current' uses to
+ * reach them.
+ */
+#define GET_USERREG() ((struct pt_regs *)(THREAD_START_SP + (unsigned long)current_thread_info()) - 1)
+
+#include <linux/thread_info.h>
+
+/* includes */
+#include "fpsr.h"		/* FP control and status register definitions */
+#include "milieu.h"
+
+struct roundingData {
+    int8 mode;
+    int8 precision;
+    signed char exception;
+};
+
+#include "softfloat.h"
+
+#define		typeNone		0x00
+#define		typeSingle		0x01
+#define		typeDouble		0x02
+#define		typeExtended		0x03
+
+/*
+ * This must be no more and no less than 12 bytes.
+ */
+typedef union tagFPREG {
+	float32 fSingle;
+	float64 fDouble;
+#ifdef CONFIG_FPE_NWFPE_XP
+	floatx80 fExtended;
+#else
+	u32 padding[3];
+#endif
+} __attribute__ ((packed,aligned(4))) FPREG;
+
+/*
+ * FPA11 device model.
+ *
+ * This structure is exported to user space.  Do not re-order.
+ * Only add new stuff to the end, and do not change the size of
+ * any element.  Elements of this structure are used by user
+ * space, and must match struct user_fp in <asm/user.h>.
+ * We include the byte offsets below for documentation purposes.
+ *
+ * The size of this structure and FPREG are checked by fpmodule.c
+ * on initialisation.  If the rules have been broken, NWFPE will
+ * not initialise.
+ */
+typedef struct tagFPA11 {
+/*   0 */ FPREG fpreg[8];	/* 8 floating point registers */
+/*  96 */ FPSR fpsr;		/* floating point status register */
+/* 100 */ FPCR fpcr;		/* floating point control register */
+/* 104 */ unsigned char fType[8];	/* type of floating point value held in
+					   floating point registers.  One of
+					   none, single, double or extended. */
+/* 112 */ int initflag;		/* this is special.  The kernel guarantees
+				   to set it to 0 when a thread is launched,
+				   so we can use it to detect whether this
+				   instance of the emulator needs to be
+				   initialised. */
+} __attribute__ ((packed,aligned(4))) FPA11;
+
+extern int8 SetRoundingMode(const unsigned int);
+extern int8 SetRoundingPrecision(const unsigned int);
+extern void nwfpe_init_fpa(union fp_state *fp);
+
+extern unsigned int EmulateAll(unsigned int opcode);
+
+extern unsigned int EmulateCPDT(const unsigned int opcode);
+extern unsigned int EmulateCPDO(const unsigned int opcode);
+extern unsigned int EmulateCPRT(const unsigned int opcode);
+
+/* fpa11_cpdt.c */
+extern unsigned int PerformLDF(const unsigned int opcode);
+extern unsigned int PerformSTF(const unsigned int opcode);
+extern unsigned int PerformLFM(const unsigned int opcode);
+extern unsigned int PerformSFM(const unsigned int opcode);
+
+/* single_cpdo.c */
+
+extern unsigned int SingleCPDO(struct roundingData *roundData,
+			       const unsigned int opcode, FPREG * rFd);
+/* double_cpdo.c */
+extern unsigned int DoubleCPDO(struct roundingData *roundData,
+			       const unsigned int opcode, FPREG * rFd);
+
+/* extneded_cpdo.c */
+extern unsigned int ExtendedCPDO(struct roundingData *roundData,
+				 const unsigned int opcode, FPREG * rFd);
+
+#endif
diff --git a/arch/arm/nwfpe/fpa11.inl b/arch/arm/nwfpe/fpa11.inl
new file mode 100644
index 000000000..ab8d68262
--- /dev/null
+++ b/arch/arm/nwfpe/fpa11.inl
@@ -0,0 +1,51 @@
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that 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., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#include "fpa11.h"
+
+/* Read and write floating point status register */
+static inline unsigned int readFPSR(void)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	return (fpa11->fpsr);
+}
+
+static inline void writeFPSR(FPSR reg)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	/* the sysid byte in the status register is readonly */
+	fpa11->fpsr = (fpa11->fpsr & MASK_SYSID) | (reg & ~MASK_SYSID);
+}
+
+/* Read and write floating point control register */
+static inline FPCR readFPCR(void)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	/* clear SB, AB and DA bits before returning FPCR */
+	return (fpa11->fpcr & ~MASK_RFC);
+}
+
+static inline void writeFPCR(FPCR reg)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	fpa11->fpcr &= ~MASK_WFC;		/* clear SB, AB and DA bits */
+	fpa11->fpcr |= (reg & MASK_WFC);	/* write SB, AB and DA bits */
+}
diff --git a/arch/arm/nwfpe/fpa11_cpdo.c b/arch/arm/nwfpe/fpa11_cpdo.c
new file mode 100644
index 000000000..631333fc5
--- /dev/null
+++ b/arch/arm/nwfpe/fpa11_cpdo.c
@@ -0,0 +1,125 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+    (c) Philip Blundell, 2001
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+#include "fpopcode.h"
+
+unsigned int SingleCPDO(struct roundingData *roundData, const unsigned int opcode, FPREG * rFd);
+unsigned int DoubleCPDO(struct roundingData *roundData, const unsigned int opcode, FPREG * rFd);
+unsigned int ExtendedCPDO(struct roundingData *roundData, const unsigned int opcode, FPREG * rFd);
+
+unsigned int EmulateCPDO(const unsigned int opcode)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	FPREG *rFd;
+	unsigned int nType, nDest, nRc;
+	struct roundingData roundData;
+
+	/* Get the destination size.  If not valid let Linux perform
+	   an invalid instruction trap. */
+	nDest = getDestinationSize(opcode);
+	if (typeNone == nDest)
+		return 0;
+
+	roundData.mode = SetRoundingMode(opcode);
+	roundData.precision = SetRoundingPrecision(opcode);
+	roundData.exception = 0;
+
+	/* Compare the size of the operands in Fn and Fm.
+	   Choose the largest size and perform operations in that size,
+	   in order to make use of all the precision of the operands.
+	   If Fm is a constant, we just grab a constant of a size
+	   matching the size of the operand in Fn. */
+	if (MONADIC_INSTRUCTION(opcode))
+		nType = nDest;
+	else
+		nType = fpa11->fType[getFn(opcode)];
+
+	if (!CONSTANT_FM(opcode)) {
+		register unsigned int Fm = getFm(opcode);
+		if (nType < fpa11->fType[Fm]) {
+			nType = fpa11->fType[Fm];
+		}
+	}
+
+	rFd = &fpa11->fpreg[getFd(opcode)];
+
+	switch (nType) {
+	case typeSingle:
+		nRc = SingleCPDO(&roundData, opcode, rFd);
+		break;
+	case typeDouble:
+		nRc = DoubleCPDO(&roundData, opcode, rFd);
+		break;
+#ifdef CONFIG_FPE_NWFPE_XP
+	case typeExtended:
+		nRc = ExtendedCPDO(&roundData, opcode, rFd);
+		break;
+#endif
+	default:
+		nRc = 0;
+	}
+
+	/* The CPDO functions used to always set the destination type
+	   to be the same as their working size. */
+
+	if (nRc != 0) {
+		/* If the operation succeeded, check to see if the result in the
+		   destination register is the correct size.  If not force it
+		   to be. */
+
+		fpa11->fType[getFd(opcode)] = nDest;
+
+#ifdef CONFIG_FPE_NWFPE_XP
+		if (nDest != nType) {
+			switch (nDest) {
+			case typeSingle:
+				{
+					if (typeDouble == nType)
+						rFd->fSingle = float64_to_float32(&roundData, rFd->fDouble);
+					else
+						rFd->fSingle = floatx80_to_float32(&roundData, rFd->fExtended);
+				}
+				break;
+
+			case typeDouble:
+				{
+					if (typeSingle == nType)
+						rFd->fDouble = float32_to_float64(rFd->fSingle);
+					else
+						rFd->fDouble = floatx80_to_float64(&roundData, rFd->fExtended);
+				}
+				break;
+
+			case typeExtended:
+				{
+					if (typeSingle == nType)
+						rFd->fExtended = float32_to_floatx80(rFd->fSingle);
+					else
+						rFd->fExtended = float64_to_floatx80(rFd->fDouble);
+				}
+				break;
+			}
+		}
+#else
+		if (nDest != nType) {
+			if (nDest == typeSingle)
+				rFd->fSingle = float64_to_float32(&roundData, rFd->fDouble);
+			else
+				rFd->fDouble = float32_to_float64(rFd->fSingle);
+		}
+#endif
+	}
+
+	if (roundData.exception)
+		float_raise(roundData.exception);
+
+	return nRc;
+}
diff --git a/arch/arm/nwfpe/fpa11_cpdt.c b/arch/arm/nwfpe/fpa11_cpdt.c
new file mode 100644
index 000000000..cd4c57fd4
--- /dev/null
+++ b/arch/arm/nwfpe/fpa11_cpdt.c
@@ -0,0 +1,395 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.com, 1998-1999
+    (c) Philip Blundell, 1998, 2001
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+#include "softfloat.h"
+#include "fpopcode.h"
+#include "fpmodule.h"
+#include "fpmodule.inl"
+
+#include <linux/uaccess.h>
+
+static inline void loadSingle(const unsigned int Fn, const unsigned int __user *pMem)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	fpa11->fType[Fn] = typeSingle;
+	get_user(fpa11->fpreg[Fn].fSingle, pMem);
+}
+
+static inline void loadDouble(const unsigned int Fn, const unsigned int __user *pMem)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	unsigned int *p;
+	p = (unsigned int *) &fpa11->fpreg[Fn].fDouble;
+	fpa11->fType[Fn] = typeDouble;
+#ifdef __ARMEB__
+	get_user(p[0], &pMem[0]);	/* sign & exponent */
+	get_user(p[1], &pMem[1]);
+#else
+	get_user(p[0], &pMem[1]);
+	get_user(p[1], &pMem[0]);	/* sign & exponent */
+#endif
+}
+
+#ifdef CONFIG_FPE_NWFPE_XP
+static inline void loadExtended(const unsigned int Fn, const unsigned int __user *pMem)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	unsigned int *p;
+	p = (unsigned int *) &fpa11->fpreg[Fn].fExtended;
+	fpa11->fType[Fn] = typeExtended;
+	get_user(p[0], &pMem[0]);	/* sign & exponent */
+#ifdef __ARMEB__
+	get_user(p[1], &pMem[1]);	/* ms bits */
+	get_user(p[2], &pMem[2]);	/* ls bits */
+#else
+	get_user(p[1], &pMem[2]);	/* ls bits */
+	get_user(p[2], &pMem[1]);	/* ms bits */
+#endif
+}
+#endif
+
+static inline void loadMultiple(const unsigned int Fn, const unsigned int __user *pMem)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	register unsigned int *p;
+	unsigned long x;
+
+	p = (unsigned int *) &(fpa11->fpreg[Fn]);
+	get_user(x, &pMem[0]);
+	fpa11->fType[Fn] = (x >> 14) & 0x00000003;
+
+	switch (fpa11->fType[Fn]) {
+	case typeSingle:
+	case typeDouble:
+		{
+			get_user(p[0], &pMem[2]);	/* Single */
+			get_user(p[1], &pMem[1]);	/* double msw */
+			p[2] = 0;			/* empty */
+		}
+		break;
+
+#ifdef CONFIG_FPE_NWFPE_XP
+	case typeExtended:
+		{
+			get_user(p[1], &pMem[2]);
+			get_user(p[2], &pMem[1]);	/* msw */
+			p[0] = (x & 0x80003fff);
+		}
+		break;
+#endif
+	}
+}
+
+static inline void storeSingle(struct roundingData *roundData, const unsigned int Fn, unsigned int __user *pMem)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	union {
+		float32 f;
+		unsigned int i[1];
+	} val;
+
+	switch (fpa11->fType[Fn]) {
+	case typeDouble:
+		val.f = float64_to_float32(roundData, fpa11->fpreg[Fn].fDouble);
+		break;
+
+#ifdef CONFIG_FPE_NWFPE_XP
+	case typeExtended:
+		val.f = floatx80_to_float32(roundData, fpa11->fpreg[Fn].fExtended);
+		break;
+#endif
+
+	default:
+		val.f = fpa11->fpreg[Fn].fSingle;
+	}
+
+	put_user(val.i[0], pMem);
+}
+
+static inline void storeDouble(struct roundingData *roundData, const unsigned int Fn, unsigned int __user *pMem)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	union {
+		float64 f;
+		unsigned int i[2];
+	} val;
+
+	switch (fpa11->fType[Fn]) {
+	case typeSingle:
+		val.f = float32_to_float64(fpa11->fpreg[Fn].fSingle);
+		break;
+
+#ifdef CONFIG_FPE_NWFPE_XP
+	case typeExtended:
+		val.f = floatx80_to_float64(roundData, fpa11->fpreg[Fn].fExtended);
+		break;
+#endif
+
+	default:
+		val.f = fpa11->fpreg[Fn].fDouble;
+	}
+
+#ifdef __ARMEB__
+	put_user(val.i[0], &pMem[0]);	/* msw */
+	put_user(val.i[1], &pMem[1]);	/* lsw */
+#else
+	put_user(val.i[1], &pMem[0]);	/* msw */
+	put_user(val.i[0], &pMem[1]);	/* lsw */
+#endif
+}
+
+#ifdef CONFIG_FPE_NWFPE_XP
+static inline void storeExtended(const unsigned int Fn, unsigned int __user *pMem)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	union {
+		floatx80 f;
+		unsigned int i[3];
+	} val;
+
+	switch (fpa11->fType[Fn]) {
+	case typeSingle:
+		val.f = float32_to_floatx80(fpa11->fpreg[Fn].fSingle);
+		break;
+
+	case typeDouble:
+		val.f = float64_to_floatx80(fpa11->fpreg[Fn].fDouble);
+		break;
+
+	default:
+		val.f = fpa11->fpreg[Fn].fExtended;
+	}
+
+	put_user(val.i[0], &pMem[0]);	/* sign & exp */
+#ifdef __ARMEB__
+	put_user(val.i[1], &pMem[1]);	/* msw */
+	put_user(val.i[2], &pMem[2]);
+#else
+	put_user(val.i[1], &pMem[2]);
+	put_user(val.i[2], &pMem[1]);	/* msw */
+#endif
+}
+#endif
+
+static inline void storeMultiple(const unsigned int Fn, unsigned int __user *pMem)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	register unsigned int nType, *p;
+
+	p = (unsigned int *) &(fpa11->fpreg[Fn]);
+	nType = fpa11->fType[Fn];
+
+	switch (nType) {
+	case typeSingle:
+	case typeDouble:
+		{
+			put_user(p[0], &pMem[2]);	/* single */
+			put_user(p[1], &pMem[1]);	/* double msw */
+			put_user(nType << 14, &pMem[0]);
+		}
+		break;
+
+#ifdef CONFIG_FPE_NWFPE_XP
+	case typeExtended:
+		{
+			put_user(p[2], &pMem[1]);	/* msw */
+			put_user(p[1], &pMem[2]);
+			put_user((p[0] & 0x80003fff) | (nType << 14), &pMem[0]);
+		}
+		break;
+#endif
+	}
+}
+
+unsigned int PerformLDF(const unsigned int opcode)
+{
+	unsigned int __user *pBase, *pAddress, *pFinal;
+	unsigned int nRc = 1, write_back = WRITE_BACK(opcode);
+
+	pBase = (unsigned int __user *) readRegister(getRn(opcode));
+	if (REG_PC == getRn(opcode)) {
+		pBase += 2;
+		write_back = 0;
+	}
+
+	pFinal = pBase;
+	if (BIT_UP_SET(opcode))
+		pFinal += getOffset(opcode);
+	else
+		pFinal -= getOffset(opcode);
+
+	if (PREINDEXED(opcode))
+		pAddress = pFinal;
+	else
+		pAddress = pBase;
+
+	switch (opcode & MASK_TRANSFER_LENGTH) {
+	case TRANSFER_SINGLE:
+		loadSingle(getFd(opcode), pAddress);
+		break;
+	case TRANSFER_DOUBLE:
+		loadDouble(getFd(opcode), pAddress);
+		break;
+#ifdef CONFIG_FPE_NWFPE_XP
+	case TRANSFER_EXTENDED:
+		loadExtended(getFd(opcode), pAddress);
+		break;
+#endif
+	default:
+		nRc = 0;
+	}
+
+	if (write_back)
+		writeRegister(getRn(opcode), (unsigned long) pFinal);
+	return nRc;
+}
+
+unsigned int PerformSTF(const unsigned int opcode)
+{
+	unsigned int __user *pBase, *pAddress, *pFinal;
+	unsigned int nRc = 1, write_back = WRITE_BACK(opcode);
+	struct roundingData roundData;
+
+	roundData.mode = SetRoundingMode(opcode);
+	roundData.precision = SetRoundingPrecision(opcode);
+	roundData.exception = 0;
+
+	pBase = (unsigned int __user *) readRegister(getRn(opcode));
+	if (REG_PC == getRn(opcode)) {
+		pBase += 2;
+		write_back = 0;
+	}
+
+	pFinal = pBase;
+	if (BIT_UP_SET(opcode))
+		pFinal += getOffset(opcode);
+	else
+		pFinal -= getOffset(opcode);
+
+	if (PREINDEXED(opcode))
+		pAddress = pFinal;
+	else
+		pAddress = pBase;
+
+	switch (opcode & MASK_TRANSFER_LENGTH) {
+	case TRANSFER_SINGLE:
+		storeSingle(&roundData, getFd(opcode), pAddress);
+		break;
+	case TRANSFER_DOUBLE:
+		storeDouble(&roundData, getFd(opcode), pAddress);
+		break;
+#ifdef CONFIG_FPE_NWFPE_XP
+	case TRANSFER_EXTENDED:
+		storeExtended(getFd(opcode), pAddress);
+		break;
+#endif
+	default:
+		nRc = 0;
+	}
+
+	if (roundData.exception)
+		float_raise(roundData.exception);
+
+	if (write_back)
+		writeRegister(getRn(opcode), (unsigned long) pFinal);
+	return nRc;
+}
+
+unsigned int PerformLFM(const unsigned int opcode)
+{
+	unsigned int __user *pBase, *pAddress, *pFinal;
+	unsigned int i, Fd, write_back = WRITE_BACK(opcode);
+
+	pBase = (unsigned int __user *) readRegister(getRn(opcode));
+	if (REG_PC == getRn(opcode)) {
+		pBase += 2;
+		write_back = 0;
+	}
+
+	pFinal = pBase;
+	if (BIT_UP_SET(opcode))
+		pFinal += getOffset(opcode);
+	else
+		pFinal -= getOffset(opcode);
+
+	if (PREINDEXED(opcode))
+		pAddress = pFinal;
+	else
+		pAddress = pBase;
+
+	Fd = getFd(opcode);
+	for (i = getRegisterCount(opcode); i > 0; i--) {
+		loadMultiple(Fd, pAddress);
+		pAddress += 3;
+		Fd++;
+		if (Fd == 8)
+			Fd = 0;
+	}
+
+	if (write_back)
+		writeRegister(getRn(opcode), (unsigned long) pFinal);
+	return 1;
+}
+
+unsigned int PerformSFM(const unsigned int opcode)
+{
+	unsigned int __user *pBase, *pAddress, *pFinal;
+	unsigned int i, Fd, write_back = WRITE_BACK(opcode);
+
+	pBase = (unsigned int __user *) readRegister(getRn(opcode));
+	if (REG_PC == getRn(opcode)) {
+		pBase += 2;
+		write_back = 0;
+	}
+
+	pFinal = pBase;
+	if (BIT_UP_SET(opcode))
+		pFinal += getOffset(opcode);
+	else
+		pFinal -= getOffset(opcode);
+
+	if (PREINDEXED(opcode))
+		pAddress = pFinal;
+	else
+		pAddress = pBase;
+
+	Fd = getFd(opcode);
+	for (i = getRegisterCount(opcode); i > 0; i--) {
+		storeMultiple(Fd, pAddress);
+		pAddress += 3;
+		Fd++;
+		if (Fd == 8)
+			Fd = 0;
+	}
+
+	if (write_back)
+		writeRegister(getRn(opcode), (unsigned long) pFinal);
+	return 1;
+}
+
+unsigned int EmulateCPDT(const unsigned int opcode)
+{
+	unsigned int nRc = 0;
+
+	if (LDF_OP(opcode)) {
+		nRc = PerformLDF(opcode);
+	} else if (LFM_OP(opcode)) {
+		nRc = PerformLFM(opcode);
+	} else if (STF_OP(opcode)) {
+		nRc = PerformSTF(opcode);
+	} else if (SFM_OP(opcode)) {
+		nRc = PerformSFM(opcode);
+	} else {
+		nRc = 0;
+	}
+
+	return nRc;
+}
diff --git a/arch/arm/nwfpe/fpa11_cprt.c b/arch/arm/nwfpe/fpa11_cprt.c
new file mode 100644
index 000000000..a8c142e45
--- /dev/null
+++ b/arch/arm/nwfpe/fpa11_cprt.c
@@ -0,0 +1,361 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+    (c) Philip Blundell, 1999, 2001
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+#include "fpopcode.h"
+#include "fpa11.inl"
+#include "fpmodule.h"
+#include "fpmodule.inl"
+#include "softfloat.h"
+
+unsigned int PerformFLT(const unsigned int opcode);
+unsigned int PerformFIX(const unsigned int opcode);
+
+static unsigned int PerformComparison(const unsigned int opcode);
+
+unsigned int EmulateCPRT(const unsigned int opcode)
+{
+
+	if (opcode & 0x800000) {
+		/* This is some variant of a comparison (PerformComparison
+		   will sort out which one).  Since most of the other CPRT
+		   instructions are oddball cases of some sort or other it
+		   makes sense to pull this out into a fast path.  */
+		return PerformComparison(opcode);
+	}
+
+	/* Hint to GCC that we'd like a jump table rather than a load of CMPs */
+	switch ((opcode & 0x700000) >> 20) {
+	case FLT_CODE >> 20:
+		return PerformFLT(opcode);
+		break;
+	case FIX_CODE >> 20:
+		return PerformFIX(opcode);
+		break;
+
+	case WFS_CODE >> 20:
+		writeFPSR(readRegister(getRd(opcode)));
+		break;
+	case RFS_CODE >> 20:
+		writeRegister(getRd(opcode), readFPSR());
+		break;
+
+	default:
+		return 0;
+	}
+
+	return 1;
+}
+
+unsigned int PerformFLT(const unsigned int opcode)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	struct roundingData roundData;
+
+	roundData.mode = SetRoundingMode(opcode);
+	roundData.precision = SetRoundingPrecision(opcode);
+	roundData.exception = 0;
+
+	switch (opcode & MASK_ROUNDING_PRECISION) {
+	case ROUND_SINGLE:
+		{
+			fpa11->fType[getFn(opcode)] = typeSingle;
+			fpa11->fpreg[getFn(opcode)].fSingle = int32_to_float32(&roundData, readRegister(getRd(opcode)));
+		}
+		break;
+
+	case ROUND_DOUBLE:
+		{
+			fpa11->fType[getFn(opcode)] = typeDouble;
+			fpa11->fpreg[getFn(opcode)].fDouble = int32_to_float64(readRegister(getRd(opcode)));
+		}
+		break;
+
+#ifdef CONFIG_FPE_NWFPE_XP
+	case ROUND_EXTENDED:
+		{
+			fpa11->fType[getFn(opcode)] = typeExtended;
+			fpa11->fpreg[getFn(opcode)].fExtended = int32_to_floatx80(readRegister(getRd(opcode)));
+		}
+		break;
+#endif
+
+	default:
+		return 0;
+	}
+
+	if (roundData.exception)
+		float_raise(roundData.exception);
+
+	return 1;
+}
+
+unsigned int PerformFIX(const unsigned int opcode)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	unsigned int Fn = getFm(opcode);
+	struct roundingData roundData;
+
+	roundData.mode = SetRoundingMode(opcode);
+	roundData.precision = SetRoundingPrecision(opcode);
+	roundData.exception = 0;
+
+	switch (fpa11->fType[Fn]) {
+	case typeSingle:
+		{
+			writeRegister(getRd(opcode), float32_to_int32(&roundData, fpa11->fpreg[Fn].fSingle));
+		}
+		break;
+
+	case typeDouble:
+		{
+			writeRegister(getRd(opcode), float64_to_int32(&roundData, fpa11->fpreg[Fn].fDouble));
+		}
+		break;
+
+#ifdef CONFIG_FPE_NWFPE_XP
+	case typeExtended:
+		{
+			writeRegister(getRd(opcode), floatx80_to_int32(&roundData, fpa11->fpreg[Fn].fExtended));
+		}
+		break;
+#endif
+
+	default:
+		return 0;
+	}
+
+	if (roundData.exception)
+		float_raise(roundData.exception);
+
+	return 1;
+}
+
+/* This instruction sets the flags N, Z, C, V in the FPSR. */
+static unsigned int PerformComparison(const unsigned int opcode)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	unsigned int Fn = getFn(opcode), Fm = getFm(opcode);
+	int e_flag = opcode & 0x400000;	/* 1 if CxFE */
+	int n_flag = opcode & 0x200000;	/* 1 if CNxx */
+	unsigned int flags = 0;
+
+#ifdef CONFIG_FPE_NWFPE_XP
+	floatx80 rFn, rFm;
+
+	/* Check for unordered condition and convert all operands to 80-bit
+	   format.
+	   ?? Might be some mileage in avoiding this conversion if possible.
+	   Eg, if both operands are 32-bit, detect this and do a 32-bit
+	   comparison (cheaper than an 80-bit one).  */
+	switch (fpa11->fType[Fn]) {
+	case typeSingle:
+		//printk("single.\n");
+		if (float32_is_nan(fpa11->fpreg[Fn].fSingle))
+			goto unordered;
+		rFn = float32_to_floatx80(fpa11->fpreg[Fn].fSingle);
+		break;
+
+	case typeDouble:
+		//printk("double.\n");
+		if (float64_is_nan(fpa11->fpreg[Fn].fDouble))
+			goto unordered;
+		rFn = float64_to_floatx80(fpa11->fpreg[Fn].fDouble);
+		break;
+
+	case typeExtended:
+		//printk("extended.\n");
+		if (floatx80_is_nan(fpa11->fpreg[Fn].fExtended))
+			goto unordered;
+		rFn = fpa11->fpreg[Fn].fExtended;
+		break;
+
+	default:
+		return 0;
+	}
+
+	if (CONSTANT_FM(opcode)) {
+		//printk("Fm is a constant: #%d.\n",Fm);
+		rFm = getExtendedConstant(Fm);
+		if (floatx80_is_nan(rFm))
+			goto unordered;
+	} else {
+		//printk("Fm = r%d which contains a ",Fm);
+		switch (fpa11->fType[Fm]) {
+		case typeSingle:
+			//printk("single.\n");
+			if (float32_is_nan(fpa11->fpreg[Fm].fSingle))
+				goto unordered;
+			rFm = float32_to_floatx80(fpa11->fpreg[Fm].fSingle);
+			break;
+
+		case typeDouble:
+			//printk("double.\n");
+			if (float64_is_nan(fpa11->fpreg[Fm].fDouble))
+				goto unordered;
+			rFm = float64_to_floatx80(fpa11->fpreg[Fm].fDouble);
+			break;
+
+		case typeExtended:
+			//printk("extended.\n");
+			if (floatx80_is_nan(fpa11->fpreg[Fm].fExtended))
+				goto unordered;
+			rFm = fpa11->fpreg[Fm].fExtended;
+			break;
+
+		default:
+			return 0;
+		}
+	}
+
+	if (n_flag)
+		rFm.high ^= 0x8000;
+
+	/* test for less than condition */
+	if (floatx80_lt(rFn, rFm))
+		flags |= CC_NEGATIVE;
+
+	/* test for equal condition */
+	if (floatx80_eq(rFn, rFm))
+		flags |= CC_ZERO;
+
+	/* test for greater than or equal condition */
+	if (floatx80_lt(rFm, rFn))
+		flags |= CC_CARRY;
+
+#else
+	if (CONSTANT_FM(opcode)) {
+		/* Fm is a constant.  Do the comparison in whatever precision
+		   Fn happens to be stored in.  */
+		if (fpa11->fType[Fn] == typeSingle) {
+			float32 rFm = getSingleConstant(Fm);
+			float32 rFn = fpa11->fpreg[Fn].fSingle;
+
+			if (float32_is_nan(rFn))
+				goto unordered;
+
+			if (n_flag)
+				rFm ^= 0x80000000;
+
+			/* test for less than condition */
+			if (float32_lt_nocheck(rFn, rFm))
+				flags |= CC_NEGATIVE;
+
+			/* test for equal condition */
+			if (float32_eq_nocheck(rFn, rFm))
+				flags |= CC_ZERO;
+
+			/* test for greater than or equal condition */
+			if (float32_lt_nocheck(rFm, rFn))
+				flags |= CC_CARRY;
+		} else {
+			float64 rFm = getDoubleConstant(Fm);
+			float64 rFn = fpa11->fpreg[Fn].fDouble;
+
+			if (float64_is_nan(rFn))
+				goto unordered;
+
+			if (n_flag)
+				rFm ^= 0x8000000000000000ULL;
+
+			/* test for less than condition */
+			if (float64_lt_nocheck(rFn, rFm))
+				flags |= CC_NEGATIVE;
+
+			/* test for equal condition */
+			if (float64_eq_nocheck(rFn, rFm))
+				flags |= CC_ZERO;
+
+			/* test for greater than or equal condition */
+			if (float64_lt_nocheck(rFm, rFn))
+				flags |= CC_CARRY;
+		}
+	} else {
+		/* Both operands are in registers.  */
+		if (fpa11->fType[Fn] == typeSingle
+		    && fpa11->fType[Fm] == typeSingle) {
+			float32 rFm = fpa11->fpreg[Fm].fSingle;
+			float32 rFn = fpa11->fpreg[Fn].fSingle;
+
+			if (float32_is_nan(rFn)
+			    || float32_is_nan(rFm))
+				goto unordered;
+
+			if (n_flag)
+				rFm ^= 0x80000000;
+
+			/* test for less than condition */
+			if (float32_lt_nocheck(rFn, rFm))
+				flags |= CC_NEGATIVE;
+
+			/* test for equal condition */
+			if (float32_eq_nocheck(rFn, rFm))
+				flags |= CC_ZERO;
+
+			/* test for greater than or equal condition */
+			if (float32_lt_nocheck(rFm, rFn))
+				flags |= CC_CARRY;
+		} else {
+			/* Promote 32-bit operand to 64 bits.  */
+			float64 rFm, rFn;
+
+			rFm = (fpa11->fType[Fm] == typeSingle) ?
+			    float32_to_float64(fpa11->fpreg[Fm].fSingle)
+			    : fpa11->fpreg[Fm].fDouble;
+
+			rFn = (fpa11->fType[Fn] == typeSingle) ?
+			    float32_to_float64(fpa11->fpreg[Fn].fSingle)
+			    : fpa11->fpreg[Fn].fDouble;
+
+			if (float64_is_nan(rFn)
+			    || float64_is_nan(rFm))
+				goto unordered;
+
+			if (n_flag)
+				rFm ^= 0x8000000000000000ULL;
+
+			/* test for less than condition */
+			if (float64_lt_nocheck(rFn, rFm))
+				flags |= CC_NEGATIVE;
+
+			/* test for equal condition */
+			if (float64_eq_nocheck(rFn, rFm))
+				flags |= CC_ZERO;
+
+			/* test for greater than or equal condition */
+			if (float64_lt_nocheck(rFm, rFn))
+				flags |= CC_CARRY;
+		}
+	}
+
+#endif
+
+	writeConditionCodes(flags);
+
+	return 1;
+
+      unordered:
+	/* ?? The FPA data sheet is pretty vague about this, in particular
+	   about whether the non-E comparisons can ever raise exceptions.
+	   This implementation is based on a combination of what it says in
+	   the data sheet, observation of how the Acorn emulator actually
+	   behaves (and how programs expect it to) and guesswork.  */
+	flags |= CC_OVERFLOW;
+	flags &= ~(CC_ZERO | CC_NEGATIVE);
+
+	if (BIT_AC & readFPSR())
+		flags |= CC_CARRY;
+
+	if (e_flag)
+		float_raise(float_flag_invalid);
+
+	writeConditionCodes(flags);
+	return 1;
+}
diff --git a/arch/arm/nwfpe/fpmodule.c b/arch/arm/nwfpe/fpmodule.c
new file mode 100644
index 000000000..1122ed459
--- /dev/null
+++ b/arch/arm/nwfpe/fpmodule.c
@@ -0,0 +1,178 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.com, 1998-1999
+    (c) Philip Blundell, 1998-1999
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+
+#include <linux/module.h>
+#include <linux/moduleparam.h>
+
+/* XXX */
+#include <linux/errno.h>
+#include <linux/types.h>
+#include <linux/kernel.h>
+#include <linux/signal.h>
+#include <linux/sched/signal.h>
+#include <linux/init.h>
+
+#include <asm/thread_notify.h>
+
+#include "softfloat.h"
+#include "fpopcode.h"
+#include "fpmodule.h"
+#include "fpa11.inl"
+
+/* kernel symbols required for signal handling */
+#ifdef CONFIG_FPE_NWFPE_XP
+#define NWFPE_BITS "extended"
+#else
+#define NWFPE_BITS "double"
+#endif
+
+#ifdef MODULE
+void fp_send_sig(unsigned long sig, struct task_struct *p, int priv);
+#else
+#define fp_send_sig	send_sig
+#define kern_fp_enter	fp_enter
+
+extern char fpe_type[];
+#endif
+
+static int nwfpe_notify(struct notifier_block *self, unsigned long cmd, void *v)
+{
+	struct thread_info *thread = v;
+
+	if (cmd == THREAD_NOTIFY_FLUSH)
+		nwfpe_init_fpa(&thread->fpstate);
+
+	return NOTIFY_DONE;
+}
+
+static struct notifier_block nwfpe_notifier_block = {
+	.notifier_call = nwfpe_notify,
+};
+
+/* kernel function prototypes required */
+void fp_setup(void);
+
+/* external declarations for saved kernel symbols */
+extern void (*kern_fp_enter)(void);
+
+/* Original value of fp_enter from kernel before patched by fpe_init. */
+static void (*orig_fp_enter)(void);
+
+/* forward declarations */
+extern void nwfpe_enter(void);
+
+static int __init fpe_init(void)
+{
+	if (sizeof(FPA11) > sizeof(union fp_state)) {
+		pr_err("nwfpe: bad structure size\n");
+		return -EINVAL;
+	}
+
+	if (sizeof(FPREG) != 12) {
+		pr_err("nwfpe: bad register size\n");
+		return -EINVAL;
+	}
+	if (fpe_type[0] && strcmp(fpe_type, "nwfpe"))
+		return 0;
+
+	/* Display title, version and copyright information. */
+	pr_info("NetWinder Floating Point Emulator V0.97 ("
+	        NWFPE_BITS " precision)\n");
+
+	thread_register_notifier(&nwfpe_notifier_block);
+
+	/* Save pointer to the old FP handler and then patch ourselves in */
+	orig_fp_enter = kern_fp_enter;
+	kern_fp_enter = nwfpe_enter;
+
+	return 0;
+}
+
+static void __exit fpe_exit(void)
+{
+	thread_unregister_notifier(&nwfpe_notifier_block);
+	/* Restore the values we saved earlier. */
+	kern_fp_enter = orig_fp_enter;
+}
+
+/*
+ScottB:  November 4, 1998
+
+Moved this function out of softfloat-specialize into fpmodule.c.
+This effectively isolates all the changes required for integrating with the
+Linux kernel into fpmodule.c.  Porting to NetBSD should only require modifying
+fpmodule.c to integrate with the NetBSD kernel (I hope!).
+
+[1/1/99: Not quite true any more unfortunately.  There is Linux-specific
+code to access data in user space in some other source files at the 
+moment (grep for get_user / put_user calls).  --philb]
+
+This function is called by the SoftFloat routines to raise a floating
+point exception.  We check the trap enable byte in the FPSR, and raise
+a SIGFPE exception if necessary.  If not the relevant bits in the 
+cumulative exceptions flag byte are set and we return.
+*/
+
+#ifdef CONFIG_DEBUG_USER
+/* By default, ignore inexact errors as there are far too many of them to log */
+static int debug = ~BIT_IXC;
+#endif
+
+void float_raise(signed char flags)
+{
+	register unsigned int fpsr, cumulativeTraps;
+
+#ifdef CONFIG_DEBUG_USER
+	if (flags & debug)
+ 		printk(KERN_DEBUG
+		       "NWFPE: %s[%d] takes exception %08x at %ps from %08lx\n",
+		       current->comm, current->pid, flags,
+		       __builtin_return_address(0), GET_USERREG()->ARM_pc);
+#endif
+
+	/* Read fpsr and initialize the cumulativeTraps.  */
+	fpsr = readFPSR();
+	cumulativeTraps = 0;
+
+	/* For each type of exception, the cumulative trap exception bit is only
+	   set if the corresponding trap enable bit is not set.  */
+	if ((!(fpsr & BIT_IXE)) && (flags & BIT_IXC))
+		cumulativeTraps |= BIT_IXC;
+	if ((!(fpsr & BIT_UFE)) && (flags & BIT_UFC))
+		cumulativeTraps |= BIT_UFC;
+	if ((!(fpsr & BIT_OFE)) && (flags & BIT_OFC))
+		cumulativeTraps |= BIT_OFC;
+	if ((!(fpsr & BIT_DZE)) && (flags & BIT_DZC))
+		cumulativeTraps |= BIT_DZC;
+	if ((!(fpsr & BIT_IOE)) && (flags & BIT_IOC))
+		cumulativeTraps |= BIT_IOC;
+
+	/* Set the cumulative exceptions flags.  */
+	if (cumulativeTraps)
+		writeFPSR(fpsr | cumulativeTraps);
+
+	/* Raise an exception if necessary.  */
+	if (fpsr & (flags << 16))
+		fp_send_sig(SIGFPE, current, 1);
+}
+
+module_init(fpe_init);
+module_exit(fpe_exit);
+
+MODULE_AUTHOR("Scott Bambrough <scottb@rebel.com>");
+MODULE_DESCRIPTION("NWFPE floating point emulator (" NWFPE_BITS " precision)");
+MODULE_LICENSE("GPL");
+
+#ifdef CONFIG_DEBUG_USER
+module_param(debug, int, 0644);
+#endif
diff --git a/arch/arm/nwfpe/fpmodule.h b/arch/arm/nwfpe/fpmodule.h
new file mode 100644
index 000000000..77a3a60a0
--- /dev/null
+++ b/arch/arm/nwfpe/fpmodule.h
@@ -0,0 +1,35 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.com, 1998-1999
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+
+    27/03/03 Ian Molton Clean up CONFIG_CPU
+*/
+
+#ifndef __FPMODULE_H__
+#define __FPMODULE_H__
+
+#define REG_ORIG_R0	17
+#define REG_CPSR	16
+#define REG_PC		15
+#define REG_LR		14
+#define REG_SP		13
+#define REG_IP		12
+#define REG_FP		11
+#define REG_R10		10
+#define REG_R9		9
+#define REG_R9		9
+#define REG_R8		8
+#define REG_R7		7
+#define REG_R6		6
+#define REG_R5		5
+#define REG_R4		4
+#define REG_R3		3
+#define REG_R2		2
+#define REG_R1		1
+#define REG_R0		0
+
+#endif
diff --git a/arch/arm/nwfpe/fpmodule.inl b/arch/arm/nwfpe/fpmodule.inl
new file mode 100644
index 000000000..2c39ad408
--- /dev/null
+++ b/arch/arm/nwfpe/fpmodule.inl
@@ -0,0 +1,74 @@
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that 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., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+static inline unsigned long readRegister(const unsigned int nReg)
+{
+	/* Note: The CPU thinks it has dealt with the current instruction.
+	   As a result the program counter has been advanced to the next
+	   instruction, and points 4 bytes beyond the actual instruction
+	   that caused the invalid instruction trap to occur.  We adjust
+	   for this in this routine.  LDF/STF instructions with Rn = PC
+	   depend on the PC being correct, as they use PC+8 in their
+	   address calculations. */
+	struct pt_regs *regs = GET_USERREG();
+	unsigned int val = regs->uregs[nReg];
+	if (REG_PC == nReg)
+		val -= 4;
+	return val;
+}
+
+static inline void
+writeRegister(const unsigned int nReg, const unsigned long val)
+{
+	struct pt_regs *regs = GET_USERREG();
+	regs->uregs[nReg] = val;
+}
+
+static inline unsigned long readCPSR(void)
+{
+	return (readRegister(REG_CPSR));
+}
+
+static inline void writeCPSR(const unsigned long val)
+{
+	writeRegister(REG_CPSR, val);
+}
+
+static inline unsigned long readConditionCodes(void)
+{
+#ifdef __FPEM_TEST__
+	return (0);
+#else
+	return (readCPSR() & CC_MASK);
+#endif
+}
+
+static inline void writeConditionCodes(const unsigned long val)
+{
+	struct pt_regs *regs = GET_USERREG();
+	unsigned long rval;
+	/*
+	 * Operate directly on userRegisters since
+	 * the CPSR may be the PC register itself.
+	 */
+	rval = regs->ARM_cpsr & ~CC_MASK;
+	regs->ARM_cpsr = rval | (val & CC_MASK);
+}
diff --git a/arch/arm/nwfpe/fpopcode.c b/arch/arm/nwfpe/fpopcode.c
new file mode 100644
index 000000000..c680bf1bb
--- /dev/null
+++ b/arch/arm/nwfpe/fpopcode.c
@@ -0,0 +1,51 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+#include "softfloat.h"
+#include "fpopcode.h"
+#include "fpsr.h"
+#include "fpmodule.h"
+#include "fpmodule.inl"
+
+#ifdef CONFIG_FPE_NWFPE_XP
+const floatx80 floatx80Constant[] = {
+	{ .high = 0x0000, .low = 0x0000000000000000ULL},/* extended 0.0 */
+	{ .high = 0x3fff, .low = 0x8000000000000000ULL},/* extended 1.0 */
+	{ .high = 0x4000, .low = 0x8000000000000000ULL},/* extended 2.0 */
+	{ .high = 0x4000, .low = 0xc000000000000000ULL},/* extended 3.0 */
+	{ .high = 0x4001, .low = 0x8000000000000000ULL},/* extended 4.0 */
+	{ .high = 0x4001, .low = 0xa000000000000000ULL},/* extended 5.0 */
+	{ .high = 0x3ffe, .low = 0x8000000000000000ULL},/* extended 0.5 */
+	{ .high = 0x4002, .low = 0xa000000000000000ULL},/* extended 10.0 */
+};
+#endif
+
+const float64 float64Constant[] = {
+	0x0000000000000000ULL,	/* double 0.0 */
+	0x3ff0000000000000ULL,	/* double 1.0 */
+	0x4000000000000000ULL,	/* double 2.0 */
+	0x4008000000000000ULL,	/* double 3.0 */
+	0x4010000000000000ULL,	/* double 4.0 */
+	0x4014000000000000ULL,	/* double 5.0 */
+	0x3fe0000000000000ULL,	/* double 0.5 */
+	0x4024000000000000ULL	/* double 10.0 */
+};
+
+const float32 float32Constant[] = {
+	0x00000000,		/* single 0.0 */
+	0x3f800000,		/* single 1.0 */
+	0x40000000,		/* single 2.0 */
+	0x40400000,		/* single 3.0 */
+	0x40800000,		/* single 4.0 */
+	0x40a00000,		/* single 5.0 */
+	0x3f000000,		/* single 0.5 */
+	0x41200000		/* single 10.0 */
+};
+
diff --git a/arch/arm/nwfpe/fpopcode.h b/arch/arm/nwfpe/fpopcode.h
new file mode 100644
index 000000000..5729f1e1b
--- /dev/null
+++ b/arch/arm/nwfpe/fpopcode.h
@@ -0,0 +1,469 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+    (c) Philip Blundell, 2001
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#ifndef __FPOPCODE_H__
+#define __FPOPCODE_H__
+
+
+/*
+ARM Floating Point Instruction Classes
+| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 
+|c o n d|1 1 0 P|U|u|W|L|   Rn  |v|  Fd |0|0|0|1|  o f f s e t  | CPDT
+|c o n d|1 1 0 P|U|w|W|L|   Rn  |x|  Fd |0|0|1|0|  o f f s e t  | CPDT (copro 2)
+| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 
+|c o n d|1 1 1 0|a|b|c|d|e|  Fn |j|  Fd |0|0|0|1|f|g|h|0|i|  Fm | CPDO
+|c o n d|1 1 1 0|a|b|c|L|e|  Fn |   Rd  |0|0|0|1|f|g|h|1|i|  Fm | CPRT
+|c o n d|1 1 1 0|a|b|c|1|e|  Fn |1|1|1|1|0|0|0|1|f|g|h|1|i|  Fm | comparisons
+| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 
+
+CPDT		data transfer instructions
+		LDF, STF, LFM (copro 2), SFM (copro 2)
+		
+CPDO		dyadic arithmetic instructions
+		ADF, MUF, SUF, RSF, DVF, RDF,
+		POW, RPW, RMF, FML, FDV, FRD, POL
+
+CPDO		monadic arithmetic instructions
+		MVF, MNF, ABS, RND, SQT, LOG, LGN, EXP,
+		SIN, COS, TAN, ASN, ACS, ATN, URD, NRM
+		
+CPRT		joint arithmetic/data transfer instructions
+		FIX (arithmetic followed by load/store)
+		FLT (load/store followed by arithmetic)
+		CMF, CNF CMFE, CNFE (comparisons)
+		WFS, RFS (write/read floating point status register)
+		WFC, RFC (write/read floating point control register)
+
+cond		condition codes
+P		pre/post index bit: 0 = postindex, 1 = preindex
+U		up/down bit: 0 = stack grows down, 1 = stack grows up
+W		write back bit: 1 = update base register (Rn)
+L		load/store bit: 0 = store, 1 = load
+Rn		base register
+Rd		destination/source register		
+Fd		floating point destination register
+Fn		floating point source register
+Fm		floating point source register or floating point constant
+
+uv		transfer length (TABLE 1)
+wx		register count (TABLE 2)
+abcd		arithmetic opcode (TABLES 3 & 4)
+ef		destination size (rounding precision) (TABLE 5)
+gh		rounding mode (TABLE 6)
+j		dyadic/monadic bit: 0 = dyadic, 1 = monadic
+i 		constant bit: 1 = constant (TABLE 6)
+*/
+
+/*
+TABLE 1
++-------------------------+---+---+---------+---------+
+|  Precision              | u | v | FPSR.EP | length  |
++-------------------------+---+---+---------+---------+
+| Single                  | 0 | 0 |    x    | 1 words |
+| Double                  | 1 | 1 |    x    | 2 words |
+| Extended                | 1 | 1 |    x    | 3 words |
+| Packed decimal          | 1 | 1 |    0    | 3 words |
+| Expanded packed decimal | 1 | 1 |    1    | 4 words |
++-------------------------+---+---+---------+---------+
+Note: x = don't care
+*/
+
+/*
+TABLE 2
++---+---+---------------------------------+
+| w | x | Number of registers to transfer |
++---+---+---------------------------------+
+| 0 | 1 |  1                              |
+| 1 | 0 |  2                              |
+| 1 | 1 |  3                              |
+| 0 | 0 |  4                              |
++---+---+---------------------------------+
+*/
+
+/*
+TABLE 3: Dyadic Floating Point Opcodes
++---+---+---+---+----------+-----------------------+-----------------------+
+| a | b | c | d | Mnemonic | Description           | Operation             |
++---+---+---+---+----------+-----------------------+-----------------------+
+| 0 | 0 | 0 | 0 | ADF      | Add                   | Fd := Fn + Fm         |
+| 0 | 0 | 0 | 1 | MUF      | Multiply              | Fd := Fn * Fm         |
+| 0 | 0 | 1 | 0 | SUF      | Subtract              | Fd := Fn - Fm         |
+| 0 | 0 | 1 | 1 | RSF      | Reverse subtract      | Fd := Fm - Fn         |
+| 0 | 1 | 0 | 0 | DVF      | Divide                | Fd := Fn / Fm         |
+| 0 | 1 | 0 | 1 | RDF      | Reverse divide        | Fd := Fm / Fn         |
+| 0 | 1 | 1 | 0 | POW      | Power                 | Fd := Fn ^ Fm         |
+| 0 | 1 | 1 | 1 | RPW      | Reverse power         | Fd := Fm ^ Fn         |
+| 1 | 0 | 0 | 0 | RMF      | Remainder             | Fd := IEEE rem(Fn/Fm) |
+| 1 | 0 | 0 | 1 | FML      | Fast Multiply         | Fd := Fn * Fm         |
+| 1 | 0 | 1 | 0 | FDV      | Fast Divide           | Fd := Fn / Fm         |
+| 1 | 0 | 1 | 1 | FRD      | Fast reverse divide   | Fd := Fm / Fn         |
+| 1 | 1 | 0 | 0 | POL      | Polar angle (ArcTan2) | Fd := arctan2(Fn,Fm)  |
+| 1 | 1 | 0 | 1 |          | undefined instruction | trap                  |
+| 1 | 1 | 1 | 0 |          | undefined instruction | trap                  |
+| 1 | 1 | 1 | 1 |          | undefined instruction | trap                  |
++---+---+---+---+----------+-----------------------+-----------------------+
+Note: POW, RPW, POL are deprecated, and are available for backwards
+      compatibility only.
+*/
+
+/*
+TABLE 4: Monadic Floating Point Opcodes
++---+---+---+---+----------+-----------------------+-----------------------+
+| a | b | c | d | Mnemonic | Description           | Operation             |
++---+---+---+---+----------+-----------------------+-----------------------+
+| 0 | 0 | 0 | 0 | MVF      | Move                  | Fd := Fm              |
+| 0 | 0 | 0 | 1 | MNF      | Move negated          | Fd := - Fm            |
+| 0 | 0 | 1 | 0 | ABS      | Absolute value        | Fd := abs(Fm)         |
+| 0 | 0 | 1 | 1 | RND      | Round to integer      | Fd := int(Fm)         |
+| 0 | 1 | 0 | 0 | SQT      | Square root           | Fd := sqrt(Fm)        |
+| 0 | 1 | 0 | 1 | LOG      | Log base 10           | Fd := log10(Fm)       |
+| 0 | 1 | 1 | 0 | LGN      | Log base e            | Fd := ln(Fm)          |
+| 0 | 1 | 1 | 1 | EXP      | Exponent              | Fd := e ^ Fm          |
+| 1 | 0 | 0 | 0 | SIN      | Sine                  | Fd := sin(Fm)         |
+| 1 | 0 | 0 | 1 | COS      | Cosine                | Fd := cos(Fm)         |
+| 1 | 0 | 1 | 0 | TAN      | Tangent               | Fd := tan(Fm)         |
+| 1 | 0 | 1 | 1 | ASN      | Arc Sine              | Fd := arcsin(Fm)      |
+| 1 | 1 | 0 | 0 | ACS      | Arc Cosine            | Fd := arccos(Fm)      |
+| 1 | 1 | 0 | 1 | ATN      | Arc Tangent           | Fd := arctan(Fm)      |
+| 1 | 1 | 1 | 0 | URD      | Unnormalized round    | Fd := int(Fm)         |
+| 1 | 1 | 1 | 1 | NRM      | Normalize             | Fd := norm(Fm)        |
++---+---+---+---+----------+-----------------------+-----------------------+
+Note: LOG, LGN, EXP, SIN, COS, TAN, ASN, ACS, ATN are deprecated, and are
+      available for backwards compatibility only.
+*/
+
+/*
+TABLE 5
++-------------------------+---+---+
+|  Rounding Precision     | e | f |
++-------------------------+---+---+
+| IEEE Single precision   | 0 | 0 |
+| IEEE Double precision   | 0 | 1 |
+| IEEE Extended precision | 1 | 0 |
+| undefined (trap)        | 1 | 1 |
++-------------------------+---+---+
+*/
+
+/*
+TABLE 5
++---------------------------------+---+---+
+|  Rounding Mode                  | g | h |
++---------------------------------+---+---+
+| Round to nearest (default)      | 0 | 0 |
+| Round toward plus infinity      | 0 | 1 |
+| Round toward negative infinity  | 1 | 0 |
+| Round toward zero               | 1 | 1 |
++---------------------------------+---+---+
+*/
+
+/*
+===
+=== Definitions for load and store instructions
+===
+*/
+
+/* bit masks */
+#define BIT_PREINDEX	0x01000000
+#define BIT_UP		0x00800000
+#define BIT_WRITE_BACK	0x00200000
+#define BIT_LOAD	0x00100000
+
+/* masks for load/store */
+#define MASK_CPDT		0x0c000000	/* data processing opcode */
+#define MASK_OFFSET		0x000000ff
+#define MASK_TRANSFER_LENGTH	0x00408000
+#define MASK_REGISTER_COUNT	MASK_TRANSFER_LENGTH
+#define MASK_COPROCESSOR	0x00000f00
+
+/* Tests for transfer length */
+#define TRANSFER_SINGLE		0x00000000
+#define TRANSFER_DOUBLE		0x00008000
+#define TRANSFER_EXTENDED	0x00400000
+#define TRANSFER_PACKED		MASK_TRANSFER_LENGTH
+
+/* Get the coprocessor number from the opcode. */
+#define getCoprocessorNumber(opcode)	((opcode & MASK_COPROCESSOR) >> 8)
+
+/* Get the offset from the opcode. */
+#define getOffset(opcode)		(opcode & MASK_OFFSET)
+
+/* Tests for specific data transfer load/store opcodes. */
+#define TEST_OPCODE(opcode,mask)	(((opcode) & (mask)) == (mask))
+
+#define LOAD_OP(opcode)   TEST_OPCODE((opcode),MASK_CPDT | BIT_LOAD)
+#define STORE_OP(opcode)  ((opcode & (MASK_CPDT | BIT_LOAD)) == MASK_CPDT)
+
+#define LDF_OP(opcode)	(LOAD_OP(opcode) && (getCoprocessorNumber(opcode) == 1))
+#define LFM_OP(opcode)	(LOAD_OP(opcode) && (getCoprocessorNumber(opcode) == 2))
+#define STF_OP(opcode)	(STORE_OP(opcode) && (getCoprocessorNumber(opcode) == 1))
+#define SFM_OP(opcode)	(STORE_OP(opcode) && (getCoprocessorNumber(opcode) == 2))
+
+#define PREINDEXED(opcode)		((opcode & BIT_PREINDEX) != 0)
+#define POSTINDEXED(opcode)		((opcode & BIT_PREINDEX) == 0)
+#define BIT_UP_SET(opcode)		((opcode & BIT_UP) != 0)
+#define BIT_UP_CLEAR(opcode)		((opcode & BIT_DOWN) == 0)
+#define WRITE_BACK(opcode)		((opcode & BIT_WRITE_BACK) != 0)
+#define LOAD(opcode)			((opcode & BIT_LOAD) != 0)
+#define STORE(opcode)			((opcode & BIT_LOAD) == 0)
+
+/*
+===
+=== Definitions for arithmetic instructions
+===
+*/
+/* bit masks */
+#define BIT_MONADIC	0x00008000
+#define BIT_CONSTANT	0x00000008
+
+#define CONSTANT_FM(opcode)		((opcode & BIT_CONSTANT) != 0)
+#define MONADIC_INSTRUCTION(opcode)	((opcode & BIT_MONADIC) != 0)
+
+/* instruction identification masks */
+#define MASK_CPDO		0x0e000000	/* arithmetic opcode */
+#define MASK_ARITHMETIC_OPCODE	0x00f08000
+#define MASK_DESTINATION_SIZE	0x00080080
+
+/* dyadic arithmetic opcodes. */
+#define ADF_CODE	0x00000000
+#define MUF_CODE	0x00100000
+#define SUF_CODE	0x00200000
+#define RSF_CODE	0x00300000
+#define DVF_CODE	0x00400000
+#define RDF_CODE	0x00500000
+#define POW_CODE	0x00600000
+#define RPW_CODE	0x00700000
+#define RMF_CODE	0x00800000
+#define FML_CODE	0x00900000
+#define FDV_CODE	0x00a00000
+#define FRD_CODE	0x00b00000
+#define POL_CODE	0x00c00000
+/* 0x00d00000 is an invalid dyadic arithmetic opcode */
+/* 0x00e00000 is an invalid dyadic arithmetic opcode */
+/* 0x00f00000 is an invalid dyadic arithmetic opcode */
+
+/* monadic arithmetic opcodes. */
+#define MVF_CODE	0x00008000
+#define MNF_CODE	0x00108000
+#define ABS_CODE	0x00208000
+#define RND_CODE	0x00308000
+#define SQT_CODE	0x00408000
+#define LOG_CODE	0x00508000
+#define LGN_CODE	0x00608000
+#define EXP_CODE	0x00708000
+#define SIN_CODE	0x00808000
+#define COS_CODE	0x00908000
+#define TAN_CODE	0x00a08000
+#define ASN_CODE	0x00b08000
+#define ACS_CODE	0x00c08000
+#define ATN_CODE	0x00d08000
+#define URD_CODE	0x00e08000
+#define NRM_CODE	0x00f08000
+
+/*
+===
+=== Definitions for register transfer and comparison instructions
+===
+*/
+
+#define MASK_CPRT		0x0e000010	/* register transfer opcode */
+#define MASK_CPRT_CODE		0x00f00000
+#define FLT_CODE		0x00000000
+#define FIX_CODE		0x00100000
+#define WFS_CODE		0x00200000
+#define RFS_CODE		0x00300000
+#define WFC_CODE		0x00400000
+#define RFC_CODE		0x00500000
+#define CMF_CODE		0x00900000
+#define CNF_CODE		0x00b00000
+#define CMFE_CODE		0x00d00000
+#define CNFE_CODE		0x00f00000
+
+/*
+===
+=== Common definitions
+===
+*/
+
+/* register masks */
+#define MASK_Rd		0x0000f000
+#define MASK_Rn		0x000f0000
+#define MASK_Fd		0x00007000
+#define MASK_Fm		0x00000007
+#define MASK_Fn		0x00070000
+
+/* condition code masks */
+#define CC_MASK		0xf0000000
+#define CC_NEGATIVE	0x80000000
+#define CC_ZERO		0x40000000
+#define CC_CARRY	0x20000000
+#define CC_OVERFLOW	0x10000000
+#define CC_EQ		0x00000000
+#define CC_NE		0x10000000
+#define CC_CS		0x20000000
+#define CC_HS		CC_CS
+#define CC_CC		0x30000000
+#define CC_LO		CC_CC
+#define CC_MI		0x40000000
+#define CC_PL		0x50000000
+#define CC_VS		0x60000000
+#define CC_VC		0x70000000
+#define CC_HI		0x80000000
+#define CC_LS		0x90000000
+#define CC_GE		0xa0000000
+#define CC_LT		0xb0000000
+#define CC_GT		0xc0000000
+#define CC_LE		0xd0000000
+#define CC_AL		0xe0000000
+#define CC_NV		0xf0000000
+
+/* rounding masks/values */
+#define MASK_ROUNDING_MODE	0x00000060
+#define ROUND_TO_NEAREST	0x00000000
+#define ROUND_TO_PLUS_INFINITY	0x00000020
+#define ROUND_TO_MINUS_INFINITY	0x00000040
+#define ROUND_TO_ZERO		0x00000060
+
+#define MASK_ROUNDING_PRECISION	0x00080080
+#define ROUND_SINGLE		0x00000000
+#define ROUND_DOUBLE		0x00000080
+#define ROUND_EXTENDED		0x00080000
+
+/* Get the condition code from the opcode. */
+#define getCondition(opcode)		(opcode >> 28)
+
+/* Get the source register from the opcode. */
+#define getRn(opcode)			((opcode & MASK_Rn) >> 16)
+
+/* Get the destination floating point register from the opcode. */
+#define getFd(opcode)			((opcode & MASK_Fd) >> 12)
+
+/* Get the first source floating point register from the opcode. */
+#define getFn(opcode)		((opcode & MASK_Fn) >> 16)
+
+/* Get the second source floating point register from the opcode. */
+#define getFm(opcode)		(opcode & MASK_Fm)
+
+/* Get the destination register from the opcode. */
+#define getRd(opcode)		((opcode & MASK_Rd) >> 12)
+
+/* Get the rounding mode from the opcode. */
+#define getRoundingMode(opcode)		((opcode & MASK_ROUNDING_MODE) >> 5)
+
+#ifdef CONFIG_FPE_NWFPE_XP
+static inline floatx80 __pure getExtendedConstant(const unsigned int nIndex)
+{
+	extern const floatx80 floatx80Constant[];
+	return floatx80Constant[nIndex];
+}
+#endif
+
+static inline float64 __pure getDoubleConstant(const unsigned int nIndex)
+{
+	extern const float64 float64Constant[];
+	return float64Constant[nIndex];
+}
+
+static inline float32 __pure getSingleConstant(const unsigned int nIndex)
+{
+	extern const float32 float32Constant[];
+	return float32Constant[nIndex];
+}
+
+static inline unsigned int getTransferLength(const unsigned int opcode)
+{
+	unsigned int nRc;
+
+	switch (opcode & MASK_TRANSFER_LENGTH) {
+	case 0x00000000:
+		nRc = 1;
+		break;		/* single precision */
+	case 0x00008000:
+		nRc = 2;
+		break;		/* double precision */
+	case 0x00400000:
+		nRc = 3;
+		break;		/* extended precision */
+	default:
+		nRc = 0;
+	}
+
+	return (nRc);
+}
+
+static inline unsigned int getRegisterCount(const unsigned int opcode)
+{
+	unsigned int nRc;
+
+	switch (opcode & MASK_REGISTER_COUNT) {
+	case 0x00000000:
+		nRc = 4;
+		break;
+	case 0x00008000:
+		nRc = 1;
+		break;
+	case 0x00400000:
+		nRc = 2;
+		break;
+	case 0x00408000:
+		nRc = 3;
+		break;
+	default:
+		nRc = 0;
+	}
+
+	return (nRc);
+}
+
+static inline unsigned int getRoundingPrecision(const unsigned int opcode)
+{
+	unsigned int nRc;
+
+	switch (opcode & MASK_ROUNDING_PRECISION) {
+	case 0x00000000:
+		nRc = 1;
+		break;
+	case 0x00000080:
+		nRc = 2;
+		break;
+	case 0x00080000:
+		nRc = 3;
+		break;
+	default:
+		nRc = 0;
+	}
+
+	return (nRc);
+}
+
+static inline unsigned int getDestinationSize(const unsigned int opcode)
+{
+	unsigned int nRc;
+
+	switch (opcode & MASK_DESTINATION_SIZE) {
+	case 0x00000000:
+		nRc = typeSingle;
+		break;
+	case 0x00000080:
+		nRc = typeDouble;
+		break;
+	case 0x00080000:
+		nRc = typeExtended;
+		break;
+	default:
+		nRc = typeNone;
+	}
+
+	return (nRc);
+}
+
+extern const float64 float64Constant[];
+extern const float32 float32Constant[];
+
+#endif
diff --git a/arch/arm/nwfpe/fpsr.h b/arch/arm/nwfpe/fpsr.h
new file mode 100644
index 000000000..62c996fdd
--- /dev/null
+++ b/arch/arm/nwfpe/fpsr.h
@@ -0,0 +1,96 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.com, 1998-1999
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#ifndef __FPSR_H__
+#define __FPSR_H__
+
+/*
+The FPSR is a 32 bit register consisting of 4 parts, each exactly
+one byte.
+
+	SYSTEM ID
+	EXCEPTION TRAP ENABLE BYTE
+	SYSTEM CONTROL BYTE
+	CUMULATIVE EXCEPTION FLAGS BYTE
+
+The FPCR is a 32 bit register consisting of bit flags.
+*/
+
+/* SYSTEM ID
+------------
+Note: the system id byte is read only  */
+
+typedef unsigned int FPSR;	/* type for floating point status register */
+typedef unsigned int FPCR;	/* type for floating point control register */
+
+#define MASK_SYSID		0xff000000
+#define BIT_HARDWARE		0x80000000
+#define FP_EMULATOR		0x01000000	/* System ID for emulator */
+#define FP_ACCELERATOR		0x81000000	/* System ID for FPA11 */
+
+/* EXCEPTION TRAP ENABLE BYTE
+----------------------------- */
+
+#define MASK_TRAP_ENABLE	0x00ff0000
+#define MASK_TRAP_ENABLE_STRICT	0x001f0000
+#define BIT_IXE		0x00100000	/* inexact exception enable */
+#define BIT_UFE		0x00080000	/* underflow exception enable */
+#define BIT_OFE		0x00040000	/* overflow exception enable */
+#define BIT_DZE		0x00020000	/* divide by zero exception enable */
+#define BIT_IOE		0x00010000	/* invalid operation exception enable */
+
+/* SYSTEM CONTROL BYTE
+---------------------- */
+
+#define MASK_SYSTEM_CONTROL	0x0000ff00
+#define MASK_TRAP_STRICT	0x00001f00
+
+#define BIT_AC	0x00001000	/* use alternative C-flag definition
+				   for compares */
+#define BIT_EP	0x00000800	/* use expanded packed decimal format */
+#define BIT_SO	0x00000400	/* select synchronous operation of FPA */
+#define BIT_NE	0x00000200	/* NaN exception bit */
+#define BIT_ND	0x00000100	/* no denormalized numbers bit */
+
+/* CUMULATIVE EXCEPTION FLAGS BYTE
+---------------------------------- */
+
+#define MASK_EXCEPTION_FLAGS		0x000000ff
+#define MASK_EXCEPTION_FLAGS_STRICT	0x0000001f
+
+#define BIT_IXC		0x00000010	/* inexact exception flag */
+#define BIT_UFC		0x00000008	/* underflow exception flag */
+#define BIT_OFC		0x00000004	/* overfloat exception flag */
+#define BIT_DZC		0x00000002	/* divide by zero exception flag */
+#define BIT_IOC		0x00000001	/* invalid operation exception flag */
+
+/* Floating Point Control Register
+----------------------------------*/
+
+#define BIT_RU		0x80000000	/* rounded up bit */
+#define BIT_IE		0x10000000	/* inexact bit */
+#define BIT_MO		0x08000000	/* mantissa overflow bit */
+#define BIT_EO		0x04000000	/* exponent overflow bit */
+#define BIT_SB		0x00000800	/* store bounce */
+#define BIT_AB		0x00000400	/* arithmetic bounce */
+#define BIT_RE		0x00000200	/* rounding exception */
+#define BIT_DA		0x00000100	/* disable FPA */
+
+#define MASK_OP		0x00f08010	/* AU operation code */
+#define MASK_PR		0x00080080	/* AU precision */
+#define MASK_S1		0x00070000	/* AU source register 1 */
+#define MASK_S2		0x00000007	/* AU source register 2 */
+#define MASK_DS		0x00007000	/* AU destination register */
+#define MASK_RM		0x00000060	/* AU rounding mode */
+#define MASK_ALU	0x9cfff2ff	/* only ALU can write these bits */
+#define MASK_RESET	0x00000d00	/* bits set on reset, all others cleared */
+#define MASK_WFC	MASK_RESET
+#define MASK_RFC	~MASK_RESET
+
+#endif
diff --git a/arch/arm/nwfpe/milieu.h b/arch/arm/nwfpe/milieu.h
new file mode 100644
index 000000000..09a4f2dde
--- /dev/null
+++ b/arch/arm/nwfpe/milieu.h
@@ -0,0 +1,48 @@
+
+/*
+===============================================================================
+
+This C header file is part of the SoftFloat IEC/IEEE Floating-point
+Arithmetic Package, Release 2.
+
+Written by John R. Hauser.  This work was made possible in part by the
+International Computer Science Institute, located at Suite 600, 1947 Center
+Street, Berkeley, California 94704.  Funding was partially provided by the
+National Science Foundation under grant MIP-9311980.  The original version
+of this code was written as part of a project to build a fixed-point vector
+processor in collaboration with the University of California at Berkeley,
+overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
+is available through the Web page
+http://www.jhauser.us/arithmetic/SoftFloat-2b/SoftFloat-source.txt
+
+THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
+has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
+TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
+PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
+AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
+
+Derivative works are acceptable, even for commercial purposes, so long as
+(1) they include prominent notice that the work is derivative, and (2) they
+include prominent notice akin to these three paragraphs for those parts of
+this code that are retained.
+
+===============================================================================
+*/
+
+/*
+-------------------------------------------------------------------------------
+Include common integer types and flags.
+-------------------------------------------------------------------------------
+*/
+#include "ARM-gcc.h"
+
+/*
+-------------------------------------------------------------------------------
+Symbolic Boolean literals.
+-------------------------------------------------------------------------------
+*/
+enum {
+    FALSE = 0,
+    TRUE  = 1
+};
+
diff --git a/arch/arm/nwfpe/single_cpdo.c b/arch/arm/nwfpe/single_cpdo.c
new file mode 100644
index 000000000..4e8be20f6
--- /dev/null
+++ b/arch/arm/nwfpe/single_cpdo.c
@@ -0,0 +1,112 @@
+// SPDX-License-Identifier: GPL-2.0-or-later
+/*
+    NetWinder Floating Point Emulator
+    (c) Rebel.COM, 1998,1999
+    (c) Philip Blundell, 2001
+
+    Direct questions, comments to Scott Bambrough <scottb@netwinder.org>
+
+*/
+
+#include "fpa11.h"
+#include "softfloat.h"
+#include "fpopcode.h"
+
+float32 float32_exp(float32 Fm);
+float32 float32_ln(float32 Fm);
+float32 float32_sin(float32 rFm);
+float32 float32_cos(float32 rFm);
+float32 float32_arcsin(float32 rFm);
+float32 float32_arctan(float32 rFm);
+float32 float32_log(float32 rFm);
+float32 float32_tan(float32 rFm);
+float32 float32_arccos(float32 rFm);
+float32 float32_pow(float32 rFn, float32 rFm);
+float32 float32_pol(float32 rFn, float32 rFm);
+
+static float32 float32_rsf(struct roundingData *roundData, float32 rFn, float32 rFm)
+{
+	return float32_sub(roundData, rFm, rFn);
+}
+
+static float32 float32_rdv(struct roundingData *roundData, float32 rFn, float32 rFm)
+{
+	return float32_div(roundData, rFm, rFn);
+}
+
+static float32 (*const dyadic_single[16])(struct roundingData *, float32 rFn, float32 rFm) = {
+	[ADF_CODE >> 20] = float32_add,
+	[MUF_CODE >> 20] = float32_mul,
+	[SUF_CODE >> 20] = float32_sub,
+	[RSF_CODE >> 20] = float32_rsf,
+	[DVF_CODE >> 20] = float32_div,
+	[RDF_CODE >> 20] = float32_rdv,
+	[RMF_CODE >> 20] = float32_rem,
+
+	[FML_CODE >> 20] = float32_mul,
+	[FDV_CODE >> 20] = float32_div,
+	[FRD_CODE >> 20] = float32_rdv,
+};
+
+static float32 float32_mvf(struct roundingData *roundData, float32 rFm)
+{
+	return rFm;
+}
+
+static float32 float32_mnf(struct roundingData *roundData, float32 rFm)
+{
+	return rFm ^ 0x80000000;
+}
+
+static float32 float32_abs(struct roundingData *roundData, float32 rFm)
+{
+	return rFm & 0x7fffffff;
+}
+
+static float32 (*const monadic_single[16])(struct roundingData*, float32 rFm) = {
+	[MVF_CODE >> 20] = float32_mvf,
+	[MNF_CODE >> 20] = float32_mnf,
+	[ABS_CODE >> 20] = float32_abs,
+	[RND_CODE >> 20] = float32_round_to_int,
+	[URD_CODE >> 20] = float32_round_to_int,
+	[SQT_CODE >> 20] = float32_sqrt,
+	[NRM_CODE >> 20] = float32_mvf,
+};
+
+unsigned int SingleCPDO(struct roundingData *roundData, const unsigned int opcode, FPREG * rFd)
+{
+	FPA11 *fpa11 = GET_FPA11();
+	float32 rFm;
+	unsigned int Fm, opc_mask_shift;
+
+	Fm = getFm(opcode);
+	if (CONSTANT_FM(opcode)) {
+		rFm = getSingleConstant(Fm);
+	} else if (fpa11->fType[Fm] == typeSingle) {
+		rFm = fpa11->fpreg[Fm].fSingle;
+	} else {
+		return 0;
+	}
+
+	opc_mask_shift = (opcode & MASK_ARITHMETIC_OPCODE) >> 20;
+	if (!MONADIC_INSTRUCTION(opcode)) {
+		unsigned int Fn = getFn(opcode);
+		float32 rFn;
+
+		if (fpa11->fType[Fn] == typeSingle &&
+		    dyadic_single[opc_mask_shift]) {
+			rFn = fpa11->fpreg[Fn].fSingle;
+			rFd->fSingle = dyadic_single[opc_mask_shift](roundData, rFn, rFm);
+		} else {
+			return 0;
+		}
+	} else {
+		if (monadic_single[opc_mask_shift]) {
+			rFd->fSingle = monadic_single[opc_mask_shift](roundData, rFm);
+		} else {
+			return 0;
+		}
+	}
+
+	return 1;
+}
diff --git a/arch/arm/nwfpe/softfloat-macros b/arch/arm/nwfpe/softfloat-macros
new file mode 100644
index 000000000..cf2a61731
--- /dev/null
+++ b/arch/arm/nwfpe/softfloat-macros
@@ -0,0 +1,754 @@
+
+/*
+===============================================================================
+
+This C source fragment is part of the SoftFloat IEC/IEEE Floating-point
+Arithmetic Package, Release 2.
+
+Written by John R. Hauser.  This work was made possible in part by the
+International Computer Science Institute, located at Suite 600, 1947 Center
+Street, Berkeley, California 94704.  Funding was partially provided by the
+National Science Foundation under grant MIP-9311980.  The original version
+of this code was written as part of a project to build a fixed-point vector
+processor in collaboration with the University of California at Berkeley,
+overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
+is available through the web page
+http://www.jhauser.us/arithmetic/SoftFloat-2b/SoftFloat-source.txt
+
+THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
+has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
+TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
+PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
+AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
+
+Derivative works are acceptable, even for commercial purposes, so long as
+(1) they include prominent notice that the work is derivative, and (2) they
+include prominent notice akin to these three paragraphs for those parts of
+this code that are retained.
+
+===============================================================================
+*/
+
+/*
+-------------------------------------------------------------------------------
+Shifts `a' right by the number of bits given in `count'.  If any nonzero
+bits are shifted off, they are ``jammed'' into the least significant bit of
+the result by setting the least significant bit to 1.  The value of `count'
+can be arbitrarily large; in particular, if `count' is greater than 32, the
+result will be either 0 or 1, depending on whether `a' is zero or nonzero.
+The result is stored in the location pointed to by `zPtr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void shift32RightJamming( bits32 a, int16 count, bits32 *zPtr )
+{
+    bits32 z;
+    if ( count == 0 ) {
+        z = a;
+    }
+    else if ( count < 32 ) {
+        z = ( a>>count ) | ( ( a<<( ( - count ) & 31 ) ) != 0 );
+    }
+    else {
+        z = ( a != 0 );
+    }
+    *zPtr = z;
+}
+
+/*
+-------------------------------------------------------------------------------
+Shifts `a' right by the number of bits given in `count'.  If any nonzero
+bits are shifted off, they are ``jammed'' into the least significant bit of
+the result by setting the least significant bit to 1.  The value of `count'
+can be arbitrarily large; in particular, if `count' is greater than 64, the
+result will be either 0 or 1, depending on whether `a' is zero or nonzero.
+The result is stored in the location pointed to by `zPtr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void shift64RightJamming( bits64 a, int16 count, bits64 *zPtr )
+{
+    bits64 z;
+
+ __asm__("@shift64RightJamming -- start");   
+    if ( count == 0 ) {
+        z = a;
+    }
+    else if ( count < 64 ) {
+        z = ( a>>count ) | ( ( a<<( ( - count ) & 63 ) ) != 0 );
+    }
+    else {
+        z = ( a != 0 );
+    }
+ __asm__("@shift64RightJamming -- end");   
+    *zPtr = z;
+}
+
+/*
+-------------------------------------------------------------------------------
+Shifts the 128-bit value formed by concatenating `a0' and `a1' right by 64
+_plus_ the number of bits given in `count'.  The shifted result is at most
+64 nonzero bits; this is stored at the location pointed to by `z0Ptr'.  The
+bits shifted off form a second 64-bit result as follows:  The _last_ bit
+shifted off is the most-significant bit of the extra result, and the other
+63 bits of the extra result are all zero if and only if _all_but_the_last_
+bits shifted off were all zero.  This extra result is stored in the location
+pointed to by `z1Ptr'.  The value of `count' can be arbitrarily large.
+    (This routine makes more sense if `a0' and `a1' are considered to form a
+fixed-point value with binary point between `a0' and `a1'.  This fixed-point
+value is shifted right by the number of bits given in `count', and the
+integer part of the result is returned at the location pointed to by
+`z0Ptr'.  The fractional part of the result may be slightly corrupted as
+described above, and is returned at the location pointed to by `z1Ptr'.)
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ shift64ExtraRightJamming(
+     bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
+{
+    bits64 z0, z1;
+    int8 negCount = ( - count ) & 63;
+
+    if ( count == 0 ) {
+        z1 = a1;
+        z0 = a0;
+    }
+    else if ( count < 64 ) {
+        z1 = ( a0<<negCount ) | ( a1 != 0 );
+        z0 = a0>>count;
+    }
+    else {
+        if ( count == 64 ) {
+            z1 = a0 | ( a1 != 0 );
+        }
+        else {
+            z1 = ( ( a0 | a1 ) != 0 );
+        }
+        z0 = 0;
+    }
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Shifts the 128-bit value formed by concatenating `a0' and `a1' right by the
+number of bits given in `count'.  Any bits shifted off are lost.  The value
+of `count' can be arbitrarily large; in particular, if `count' is greater
+than 128, the result will be 0.  The result is broken into two 64-bit pieces
+which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ shift128Right(
+     bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
+{
+    bits64 z0, z1;
+    int8 negCount = ( - count ) & 63;
+
+    if ( count == 0 ) {
+        z1 = a1;
+        z0 = a0;
+    }
+    else if ( count < 64 ) {
+        z1 = ( a0<<negCount ) | ( a1>>count );
+        z0 = a0>>count;
+    }
+    else {
+        z1 = ( count < 64 ) ? ( a0>>( count & 63 ) ) : 0;
+        z0 = 0;
+    }
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Shifts the 128-bit value formed by concatenating `a0' and `a1' right by the
+number of bits given in `count'.  If any nonzero bits are shifted off, they
+are ``jammed'' into the least significant bit of the result by setting the
+least significant bit to 1.  The value of `count' can be arbitrarily large;
+in particular, if `count' is greater than 128, the result will be either 0
+or 1, depending on whether the concatenation of `a0' and `a1' is zero or
+nonzero.  The result is broken into two 64-bit pieces which are stored at
+the locations pointed to by `z0Ptr' and `z1Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ shift128RightJamming(
+     bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
+{
+    bits64 z0, z1;
+    int8 negCount = ( - count ) & 63;
+
+    if ( count == 0 ) {
+        z1 = a1;
+        z0 = a0;
+    }
+    else if ( count < 64 ) {
+        z1 = ( a0<<negCount ) | ( a1>>count ) | ( ( a1<<negCount ) != 0 );
+        z0 = a0>>count;
+    }
+    else {
+        if ( count == 64 ) {
+            z1 = a0 | ( a1 != 0 );
+        }
+        else if ( count < 128 ) {
+            z1 = ( a0>>( count & 63 ) ) | ( ( ( a0<<negCount ) | a1 ) != 0 );
+        }
+        else {
+            z1 = ( ( a0 | a1 ) != 0 );
+        }
+        z0 = 0;
+    }
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Shifts the 192-bit value formed by concatenating `a0', `a1', and `a2' right
+by 64 _plus_ the number of bits given in `count'.  The shifted result is
+at most 128 nonzero bits; these are broken into two 64-bit pieces which are
+stored at the locations pointed to by `z0Ptr' and `z1Ptr'.  The bits shifted
+off form a third 64-bit result as follows:  The _last_ bit shifted off is
+the most-significant bit of the extra result, and the other 63 bits of the
+extra result are all zero if and only if _all_but_the_last_ bits shifted off
+were all zero.  This extra result is stored in the location pointed to by
+`z2Ptr'.  The value of `count' can be arbitrarily large.
+    (This routine makes more sense if `a0', `a1', and `a2' are considered
+to form a fixed-point value with binary point between `a1' and `a2'.  This
+fixed-point value is shifted right by the number of bits given in `count',
+and the integer part of the result is returned at the locations pointed to
+by `z0Ptr' and `z1Ptr'.  The fractional part of the result may be slightly
+corrupted as described above, and is returned at the location pointed to by
+`z2Ptr'.)
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ shift128ExtraRightJamming(
+     bits64 a0,
+     bits64 a1,
+     bits64 a2,
+     int16 count,
+     bits64 *z0Ptr,
+     bits64 *z1Ptr,
+     bits64 *z2Ptr
+ )
+{
+    bits64 z0, z1, z2;
+    int8 negCount = ( - count ) & 63;
+
+    if ( count == 0 ) {
+        z2 = a2;
+        z1 = a1;
+        z0 = a0;
+    }
+    else {
+        if ( count < 64 ) {
+            z2 = a1<<negCount;
+            z1 = ( a0<<negCount ) | ( a1>>count );
+            z0 = a0>>count;
+        }
+        else {
+            if ( count == 64 ) {
+                z2 = a1;
+                z1 = a0;
+            }
+            else {
+                a2 |= a1;
+                if ( count < 128 ) {
+                    z2 = a0<<negCount;
+                    z1 = a0>>( count & 63 );
+                }
+                else {
+                    z2 = ( count == 128 ) ? a0 : ( a0 != 0 );
+                    z1 = 0;
+                }
+            }
+            z0 = 0;
+        }
+        z2 |= ( a2 != 0 );
+    }
+    *z2Ptr = z2;
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Shifts the 128-bit value formed by concatenating `a0' and `a1' left by the
+number of bits given in `count'.  Any bits shifted off are lost.  The value
+of `count' must be less than 64.  The result is broken into two 64-bit
+pieces which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ shortShift128Left(
+     bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
+{
+
+    *z1Ptr = a1<<count;
+    *z0Ptr =
+        ( count == 0 ) ? a0 : ( a0<<count ) | ( a1>>( ( - count ) & 63 ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Shifts the 192-bit value formed by concatenating `a0', `a1', and `a2' left
+by the number of bits given in `count'.  Any bits shifted off are lost.
+The value of `count' must be less than 64.  The result is broken into three
+64-bit pieces which are stored at the locations pointed to by `z0Ptr',
+`z1Ptr', and `z2Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ shortShift192Left(
+     bits64 a0,
+     bits64 a1,
+     bits64 a2,
+     int16 count,
+     bits64 *z0Ptr,
+     bits64 *z1Ptr,
+     bits64 *z2Ptr
+ )
+{
+    bits64 z0, z1, z2;
+    int8 negCount;
+
+    z2 = a2<<count;
+    z1 = a1<<count;
+    z0 = a0<<count;
+    if ( 0 < count ) {
+        negCount = ( ( - count ) & 63 );
+        z1 |= a2>>negCount;
+        z0 |= a1>>negCount;
+    }
+    *z2Ptr = z2;
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Adds the 128-bit value formed by concatenating `a0' and `a1' to the 128-bit
+value formed by concatenating `b0' and `b1'.  Addition is modulo 2^128, so
+any carry out is lost.  The result is broken into two 64-bit pieces which
+are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ add128(
+     bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 *z0Ptr, bits64 *z1Ptr )
+{
+    bits64 z1;
+
+    z1 = a1 + b1;
+    *z1Ptr = z1;
+    *z0Ptr = a0 + b0 + ( z1 < a1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Adds the 192-bit value formed by concatenating `a0', `a1', and `a2' to the
+192-bit value formed by concatenating `b0', `b1', and `b2'.  Addition is
+modulo 2^192, so any carry out is lost.  The result is broken into three
+64-bit pieces which are stored at the locations pointed to by `z0Ptr',
+`z1Ptr', and `z2Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ add192(
+     bits64 a0,
+     bits64 a1,
+     bits64 a2,
+     bits64 b0,
+     bits64 b1,
+     bits64 b2,
+     bits64 *z0Ptr,
+     bits64 *z1Ptr,
+     bits64 *z2Ptr
+ )
+{
+    bits64 z0, z1, z2;
+    int8 carry0, carry1;
+
+    z2 = a2 + b2;
+    carry1 = ( z2 < a2 );
+    z1 = a1 + b1;
+    carry0 = ( z1 < a1 );
+    z0 = a0 + b0;
+    z1 += carry1;
+    z0 += ( z1 < carry1 );
+    z0 += carry0;
+    *z2Ptr = z2;
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Subtracts the 128-bit value formed by concatenating `b0' and `b1' from the
+128-bit value formed by concatenating `a0' and `a1'.  Subtraction is modulo
+2^128, so any borrow out (carry out) is lost.  The result is broken into two
+64-bit pieces which are stored at the locations pointed to by `z0Ptr' and
+`z1Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ sub128(
+     bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 *z0Ptr, bits64 *z1Ptr )
+{
+
+    *z1Ptr = a1 - b1;
+    *z0Ptr = a0 - b0 - ( a1 < b1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Subtracts the 192-bit value formed by concatenating `b0', `b1', and `b2'
+from the 192-bit value formed by concatenating `a0', `a1', and `a2'.
+Subtraction is modulo 2^192, so any borrow out (carry out) is lost.  The
+result is broken into three 64-bit pieces which are stored at the locations
+pointed to by `z0Ptr', `z1Ptr', and `z2Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ sub192(
+     bits64 a0,
+     bits64 a1,
+     bits64 a2,
+     bits64 b0,
+     bits64 b1,
+     bits64 b2,
+     bits64 *z0Ptr,
+     bits64 *z1Ptr,
+     bits64 *z2Ptr
+ )
+{
+    bits64 z0, z1, z2;
+    int8 borrow0, borrow1;
+
+    z2 = a2 - b2;
+    borrow1 = ( a2 < b2 );
+    z1 = a1 - b1;
+    borrow0 = ( a1 < b1 );
+    z0 = a0 - b0;
+    z0 -= ( z1 < borrow1 );
+    z1 -= borrow1;
+    z0 -= borrow0;
+    *z2Ptr = z2;
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Multiplies `a' by `b' to obtain a 128-bit product.  The product is broken
+into two 64-bit pieces which are stored at the locations pointed to by
+`z0Ptr' and `z1Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void mul64To128( bits64 a, bits64 b, bits64 *z0Ptr, bits64 *z1Ptr )
+{
+    bits32 aHigh, aLow, bHigh, bLow;
+    bits64 z0, zMiddleA, zMiddleB, z1;
+
+    aLow = a;
+    aHigh = a>>32;
+    bLow = b;
+    bHigh = b>>32;
+    z1 = ( (bits64) aLow ) * bLow;
+    zMiddleA = ( (bits64) aLow ) * bHigh;
+    zMiddleB = ( (bits64) aHigh ) * bLow;
+    z0 = ( (bits64) aHigh ) * bHigh;
+    zMiddleA += zMiddleB;
+    z0 += ( ( (bits64) ( zMiddleA < zMiddleB ) )<<32 ) + ( zMiddleA>>32 );
+    zMiddleA <<= 32;
+    z1 += zMiddleA;
+    z0 += ( z1 < zMiddleA );
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Multiplies the 128-bit value formed by concatenating `a0' and `a1' by `b' to
+obtain a 192-bit product.  The product is broken into three 64-bit pieces
+which are stored at the locations pointed to by `z0Ptr', `z1Ptr', and
+`z2Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ mul128By64To192(
+     bits64 a0,
+     bits64 a1,
+     bits64 b,
+     bits64 *z0Ptr,
+     bits64 *z1Ptr,
+     bits64 *z2Ptr
+ )
+{
+    bits64 z0, z1, z2, more1;
+
+    mul64To128( a1, b, &z1, &z2 );
+    mul64To128( a0, b, &z0, &more1 );
+    add128( z0, more1, 0, z1, &z0, &z1 );
+    *z2Ptr = z2;
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Multiplies the 128-bit value formed by concatenating `a0' and `a1' to the
+128-bit value formed by concatenating `b0' and `b1' to obtain a 256-bit
+product.  The product is broken into four 64-bit pieces which are stored at
+the locations pointed to by `z0Ptr', `z1Ptr', `z2Ptr', and `z3Ptr'.
+-------------------------------------------------------------------------------
+*/
+INLINE void
+ mul128To256(
+     bits64 a0,
+     bits64 a1,
+     bits64 b0,
+     bits64 b1,
+     bits64 *z0Ptr,
+     bits64 *z1Ptr,
+     bits64 *z2Ptr,
+     bits64 *z3Ptr
+ )
+{
+    bits64 z0, z1, z2, z3;
+    bits64 more1, more2;
+
+    mul64To128( a1, b1, &z2, &z3 );
+    mul64To128( a1, b0, &z1, &more2 );
+    add128( z1, more2, 0, z2, &z1, &z2 );
+    mul64To128( a0, b0, &z0, &more1 );
+    add128( z0, more1, 0, z1, &z0, &z1 );
+    mul64To128( a0, b1, &more1, &more2 );
+    add128( more1, more2, 0, z2, &more1, &z2 );
+    add128( z0, z1, 0, more1, &z0, &z1 );
+    *z3Ptr = z3;
+    *z2Ptr = z2;
+    *z1Ptr = z1;
+    *z0Ptr = z0;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns an approximation to the 64-bit integer quotient obtained by dividing
+`b' into the 128-bit value formed by concatenating `a0' and `a1'.  The
+divisor `b' must be at least 2^63.  If q is the exact quotient truncated
+toward zero, the approximation returned lies between q and q + 2 inclusive.
+If the exact quotient q is larger than 64 bits, the maximum positive 64-bit
+unsigned integer is returned.
+-------------------------------------------------------------------------------
+*/
+static bits64 estimateDiv128To64( bits64 a0, bits64 a1, bits64 b )
+{
+    bits64 b0, b1;
+    bits64 rem0, rem1, term0, term1;
+    bits64 z;
+    if ( b <= a0 ) return LIT64( 0xFFFFFFFFFFFFFFFF );
+    b0 = b>>32;  /* hence b0 is 32 bits wide now */
+    if ( b0<<32 <= a0 ) {
+        z = LIT64( 0xFFFFFFFF00000000 );
+    }  else {
+        z = a0;
+        do_div( z, b0 );
+        z <<= 32;
+    }
+    mul64To128( b, z, &term0, &term1 );
+    sub128( a0, a1, term0, term1, &rem0, &rem1 );
+    while ( ( (sbits64) rem0 ) < 0 ) {
+        z -= LIT64( 0x100000000 );
+        b1 = b<<32;
+        add128( rem0, rem1, b0, b1, &rem0, &rem1 );
+    }
+    rem0 = ( rem0<<32 ) | ( rem1>>32 );
+    if ( b0<<32 <= rem0 ) {
+        z |= 0xFFFFFFFF;
+    } else {
+        do_div( rem0, b0 );
+        z |= rem0;
+    }
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns an approximation to the square root of the 32-bit significand given
+by `a'.  Considered as an integer, `a' must be at least 2^31.  If bit 0 of
+`aExp' (the least significant bit) is 1, the integer returned approximates
+2^31*sqrt(`a'/2^31), where `a' is considered an integer.  If bit 0 of `aExp'
+is 0, the integer returned approximates 2^31*sqrt(`a'/2^30).  In either
+case, the approximation returned lies strictly within +/-2 of the exact
+value.
+-------------------------------------------------------------------------------
+*/
+static bits32 estimateSqrt32( int16 aExp, bits32 a )
+{
+    static const bits16 sqrtOddAdjustments[] = {
+        0x0004, 0x0022, 0x005D, 0x00B1, 0x011D, 0x019F, 0x0236, 0x02E0,
+        0x039C, 0x0468, 0x0545, 0x0631, 0x072B, 0x0832, 0x0946, 0x0A67
+    };
+    static const bits16 sqrtEvenAdjustments[] = {
+        0x0A2D, 0x08AF, 0x075A, 0x0629, 0x051A, 0x0429, 0x0356, 0x029E,
+        0x0200, 0x0179, 0x0109, 0x00AF, 0x0068, 0x0034, 0x0012, 0x0002
+    };
+    int8 index;
+    bits32 z;
+    bits64 A;
+
+    index = ( a>>27 ) & 15;
+    if ( aExp & 1 ) {
+        z = 0x4000 + ( a>>17 ) - sqrtOddAdjustments[ index ];
+        z = ( ( a / z )<<14 ) + ( z<<15 );
+        a >>= 1;
+    }
+    else {
+        z = 0x8000 + ( a>>17 ) - sqrtEvenAdjustments[ index ];
+        z = a / z + z;
+        z = ( 0x20000 <= z ) ? 0xFFFF8000 : ( z<<15 );
+        if ( z <= a ) return (bits32) ( ( (sbits32) a )>>1 );
+    }
+    A = ( (bits64) a )<<31;
+    do_div( A, z );
+    return ( (bits32) A ) + ( z>>1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the number of leading 0 bits before the most-significant 1 bit
+of `a'.  If `a' is zero, 32 is returned.
+-------------------------------------------------------------------------------
+*/
+static int8 countLeadingZeros32( bits32 a )
+{
+    static const int8 countLeadingZerosHigh[] = {
+        8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
+        3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
+        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
+        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
+        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+    };
+    int8 shiftCount;
+
+    shiftCount = 0;
+    if ( a < 0x10000 ) {
+        shiftCount += 16;
+        a <<= 16;
+    }
+    if ( a < 0x1000000 ) {
+        shiftCount += 8;
+        a <<= 8;
+    }
+    shiftCount += countLeadingZerosHigh[ a>>24 ];
+    return shiftCount;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the number of leading 0 bits before the most-significant 1 bit
+of `a'.  If `a' is zero, 64 is returned.
+-------------------------------------------------------------------------------
+*/
+static int8 countLeadingZeros64( bits64 a )
+{
+    int8 shiftCount;
+
+    shiftCount = 0;
+    if ( a < ( (bits64) 1 )<<32 ) {
+        shiftCount += 32;
+    }
+    else {
+        a >>= 32;
+    }
+    shiftCount += countLeadingZeros32( a );
+    return shiftCount;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the 128-bit value formed by concatenating `a0' and `a1'
+is equal to the 128-bit value formed by concatenating `b0' and `b1'.
+Otherwise, returns 0.
+-------------------------------------------------------------------------------
+*/
+INLINE flag eq128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
+{
+
+    return ( a0 == b0 ) && ( a1 == b1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is less
+than or equal to the 128-bit value formed by concatenating `b0' and `b1'.
+Otherwise, returns 0.
+-------------------------------------------------------------------------------
+*/
+INLINE flag le128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
+{
+
+    return ( a0 < b0 ) || ( ( a0 == b0 ) && ( a1 <= b1 ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is less
+than the 128-bit value formed by concatenating `b0' and `b1'.  Otherwise,
+returns 0.
+-------------------------------------------------------------------------------
+*/
+INLINE flag lt128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
+{
+
+    return ( a0 < b0 ) || ( ( a0 == b0 ) && ( a1 < b1 ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is
+not equal to the 128-bit value formed by concatenating `b0' and `b1'.
+Otherwise, returns 0.
+-------------------------------------------------------------------------------
+*/
+INLINE flag ne128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
+{
+
+    return ( a0 != b0 ) || ( a1 != b1 );
+
+}
+
diff --git a/arch/arm/nwfpe/softfloat-specialize b/arch/arm/nwfpe/softfloat-specialize
new file mode 100644
index 000000000..679a0269d
--- /dev/null
+++ b/arch/arm/nwfpe/softfloat-specialize
@@ -0,0 +1,367 @@
+
+/*
+===============================================================================
+
+This C source fragment is part of the SoftFloat IEC/IEEE Floating-point
+Arithmetic Package, Release 2.
+
+Written by John R. Hauser.  This work was made possible in part by the
+International Computer Science Institute, located at Suite 600, 1947 Center
+Street, Berkeley, California 94704.  Funding was partially provided by the
+National Science Foundation under grant MIP-9311980.  The original version
+of this code was written as part of a project to build a fixed-point vector
+processor in collaboration with the University of California at Berkeley,
+overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
+is available through the Web page
+http://www.jhauser.us/arithmetic/SoftFloat-2b/SoftFloat-source.txt
+
+THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
+has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
+TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
+PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
+AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
+
+Derivative works are acceptable, even for commercial purposes, so long as
+(1) they include prominent notice that the work is derivative, and (2) they
+include prominent notice akin to these three paragraphs for those parts of
+this code that are retained.
+
+===============================================================================
+*/
+
+/*
+-------------------------------------------------------------------------------
+Underflow tininess-detection mode, statically initialized to default value.
+(The declaration in `softfloat.h' must match the `int8' type here.)
+-------------------------------------------------------------------------------
+*/
+int8 float_detect_tininess = float_tininess_after_rounding;
+
+/*
+-------------------------------------------------------------------------------
+Raises the exceptions specified by `flags'.  Floating-point traps can be
+defined here if desired.  It is currently not possible for such a trap to
+substitute a result value.  If traps are not implemented, this routine
+should be simply `float_exception_flags |= flags;'.
+
+ScottB:  November 4, 1998
+Moved this function out of softfloat-specialize into fpmodule.c.
+This effectively isolates all the changes required for integrating with the
+Linux kernel into fpmodule.c.  Porting to NetBSD should only require modifying
+fpmodule.c to integrate with the NetBSD kernel (I hope!).
+-------------------------------------------------------------------------------
+void float_raise( int8 flags )
+{
+    float_exception_flags |= flags;
+}
+*/
+
+/*
+-------------------------------------------------------------------------------
+Internal canonical NaN format.
+-------------------------------------------------------------------------------
+*/
+typedef struct {
+    flag sign;
+    bits64 high, low;
+} commonNaNT;
+
+/*
+-------------------------------------------------------------------------------
+The pattern for a default generated single-precision NaN.
+-------------------------------------------------------------------------------
+*/
+#define float32_default_nan 0xFFFFFFFF
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the single-precision floating-point value `a' is a NaN;
+otherwise returns 0.
+-------------------------------------------------------------------------------
+*/
+flag float32_is_nan( float32 a )
+{
+
+    return ( 0xFF000000 < (bits32) ( a<<1 ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the single-precision floating-point value `a' is a signaling
+NaN; otherwise returns 0.
+-------------------------------------------------------------------------------
+*/
+flag float32_is_signaling_nan( float32 a )
+{
+
+    return ( ( ( a>>22 ) & 0x1FF ) == 0x1FE ) && ( a & 0x003FFFFF );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the single-precision floating-point NaN
+`a' to the canonical NaN format.  If `a' is a signaling NaN, the invalid
+exception is raised.
+-------------------------------------------------------------------------------
+*/
+static commonNaNT float32ToCommonNaN( float32 a )
+{
+    commonNaNT z;
+
+    if ( float32_is_signaling_nan( a ) ) float_raise( float_flag_invalid );
+    z.sign = a>>31;
+    z.low = 0;
+    z.high = ( (bits64) a )<<41;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the canonical NaN `a' to the single-
+precision floating-point format.
+-------------------------------------------------------------------------------
+*/
+static float32 commonNaNToFloat32( commonNaNT a )
+{
+
+    return ( ( (bits32) a.sign )<<31 ) | 0x7FC00000 | ( a.high>>41 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes two single-precision floating-point values `a' and `b', one of which
+is a NaN, and returns the appropriate NaN result.  If either `a' or `b' is a
+signaling NaN, the invalid exception is raised.
+-------------------------------------------------------------------------------
+*/
+static float32 propagateFloat32NaN( float32 a, float32 b )
+{
+    flag aIsNaN, aIsSignalingNaN, bIsNaN, bIsSignalingNaN;
+
+    aIsNaN = float32_is_nan( a );
+    aIsSignalingNaN = float32_is_signaling_nan( a );
+    bIsNaN = float32_is_nan( b );
+    bIsSignalingNaN = float32_is_signaling_nan( b );
+    a |= 0x00400000;
+    b |= 0x00400000;
+    if ( aIsSignalingNaN | bIsSignalingNaN ) float_raise( float_flag_invalid );
+    if ( aIsNaN ) {
+        return ( aIsSignalingNaN & bIsNaN ) ? b : a;
+    }
+    else {
+        return b;
+    }
+
+}
+
+/*
+-------------------------------------------------------------------------------
+The pattern for a default generated double-precision NaN.
+-------------------------------------------------------------------------------
+*/
+#define float64_default_nan LIT64( 0xFFFFFFFFFFFFFFFF )
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the double-precision floating-point value `a' is a NaN;
+otherwise returns 0.
+-------------------------------------------------------------------------------
+*/
+flag float64_is_nan( float64 a )
+{
+
+    return ( LIT64( 0xFFE0000000000000 ) < (bits64) ( a<<1 ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the double-precision floating-point value `a' is a signaling
+NaN; otherwise returns 0.
+-------------------------------------------------------------------------------
+*/
+flag float64_is_signaling_nan( float64 a )
+{
+
+    return
+           ( ( ( a>>51 ) & 0xFFF ) == 0xFFE )
+        && ( a & LIT64( 0x0007FFFFFFFFFFFF ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the double-precision floating-point NaN
+`a' to the canonical NaN format.  If `a' is a signaling NaN, the invalid
+exception is raised.
+-------------------------------------------------------------------------------
+*/
+static commonNaNT float64ToCommonNaN( float64 a )
+{
+    commonNaNT z;
+
+    if ( float64_is_signaling_nan( a ) ) float_raise( float_flag_invalid );
+    z.sign = a>>63;
+    z.low = 0;
+    z.high = a<<12;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the canonical NaN `a' to the double-
+precision floating-point format.
+-------------------------------------------------------------------------------
+*/
+static float64 commonNaNToFloat64( commonNaNT a )
+{
+
+    return
+          ( ( (bits64) a.sign )<<63 )
+        | LIT64( 0x7FF8000000000000 )
+        | ( a.high>>12 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes two double-precision floating-point values `a' and `b', one of which
+is a NaN, and returns the appropriate NaN result.  If either `a' or `b' is a
+signaling NaN, the invalid exception is raised.
+-------------------------------------------------------------------------------
+*/
+static float64 propagateFloat64NaN( float64 a, float64 b )
+{
+    flag aIsNaN, aIsSignalingNaN, bIsNaN, bIsSignalingNaN;
+
+    aIsNaN = float64_is_nan( a );
+    aIsSignalingNaN = float64_is_signaling_nan( a );
+    bIsNaN = float64_is_nan( b );
+    bIsSignalingNaN = float64_is_signaling_nan( b );
+    a |= LIT64( 0x0008000000000000 );
+    b |= LIT64( 0x0008000000000000 );
+    if ( aIsSignalingNaN | bIsSignalingNaN ) float_raise( float_flag_invalid );
+    if ( aIsNaN ) {
+        return ( aIsSignalingNaN & bIsNaN ) ? b : a;
+    }
+    else {
+        return b;
+    }
+
+}
+
+#ifdef FLOATX80
+
+/*
+-------------------------------------------------------------------------------
+The pattern for a default generated extended double-precision NaN.  The
+`high' and `low' values hold the most- and least-significant bits,
+respectively.
+-------------------------------------------------------------------------------
+*/
+#define floatx80_default_nan_high 0xFFFF
+#define floatx80_default_nan_low  LIT64( 0xFFFFFFFFFFFFFFFF )
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the extended double-precision floating-point value `a' is a
+NaN; otherwise returns 0.
+-------------------------------------------------------------------------------
+*/
+flag floatx80_is_nan( floatx80 a )
+{
+
+    return ( ( a.high & 0x7FFF ) == 0x7FFF ) && (bits64) ( a.low<<1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the extended double-precision floating-point value `a' is a
+signaling NaN; otherwise returns 0.
+-------------------------------------------------------------------------------
+*/
+flag floatx80_is_signaling_nan( floatx80 a )
+{
+    //register int lr;
+    bits64 aLow;
+
+    //__asm__("mov %0, lr" : : "g" (lr));
+    //fp_printk("floatx80_is_signalling_nan() called from 0x%08x\n",lr);
+    aLow = a.low & ~ LIT64( 0x4000000000000000 );
+    return
+           ( ( a.high & 0x7FFF ) == 0x7FFF )
+        && (bits64) ( aLow<<1 )
+        && ( a.low == aLow );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the extended double-precision floating-
+point NaN `a' to the canonical NaN format.  If `a' is a signaling NaN, the
+invalid exception is raised.
+-------------------------------------------------------------------------------
+*/
+static commonNaNT floatx80ToCommonNaN( floatx80 a )
+{
+    commonNaNT z;
+
+    if ( floatx80_is_signaling_nan( a ) ) float_raise( float_flag_invalid );
+    z.sign = a.high>>15;
+    z.low = 0;
+    z.high = a.low<<1;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the canonical NaN `a' to the extended
+double-precision floating-point format.
+-------------------------------------------------------------------------------
+*/
+static floatx80 commonNaNToFloatx80( commonNaNT a )
+{
+    floatx80 z;
+
+    z.low = LIT64( 0xC000000000000000 ) | ( a.high>>1 );
+    z.high = ( ( (bits16) a.sign )<<15 ) | 0x7FFF;
+    z.__padding = 0;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes two extended double-precision floating-point values `a' and `b', one
+of which is a NaN, and returns the appropriate NaN result.  If either `a' or
+`b' is a signaling NaN, the invalid exception is raised.
+-------------------------------------------------------------------------------
+*/
+static floatx80 propagateFloatx80NaN( floatx80 a, floatx80 b )
+{
+    flag aIsNaN, aIsSignalingNaN, bIsNaN, bIsSignalingNaN;
+
+    aIsNaN = floatx80_is_nan( a );
+    aIsSignalingNaN = floatx80_is_signaling_nan( a );
+    bIsNaN = floatx80_is_nan( b );
+    bIsSignalingNaN = floatx80_is_signaling_nan( b );
+    a.low |= LIT64( 0xC000000000000000 );
+    b.low |= LIT64( 0xC000000000000000 );
+    if ( aIsSignalingNaN | bIsSignalingNaN ) float_raise( float_flag_invalid );
+    if ( aIsNaN ) {
+        return ( aIsSignalingNaN & bIsNaN ) ? b : a;
+    }
+    else {
+        return b;
+    }
+
+}
+
+#endif
diff --git a/arch/arm/nwfpe/softfloat.c b/arch/arm/nwfpe/softfloat.c
new file mode 100644
index 000000000..ffa6b4387
--- /dev/null
+++ b/arch/arm/nwfpe/softfloat.c
@@ -0,0 +1,3435 @@
+/*
+===============================================================================
+
+This C source file is part of the SoftFloat IEC/IEEE Floating-point
+Arithmetic Package, Release 2.
+
+Written by John R. Hauser.  This work was made possible in part by the
+International Computer Science Institute, located at Suite 600, 1947 Center
+Street, Berkeley, California 94704.  Funding was partially provided by the
+National Science Foundation under grant MIP-9311980.  The original version
+of this code was written as part of a project to build a fixed-point vector
+processor in collaboration with the University of California at Berkeley,
+overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
+is available through the web page
+http://www.jhauser.us/arithmetic/SoftFloat-2b/SoftFloat-source.txt
+
+THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
+has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
+TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
+PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
+AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
+
+Derivative works are acceptable, even for commercial purposes, so long as
+(1) they include prominent notice that the work is derivative, and (2) they
+include prominent notice akin to these three paragraphs for those parts of
+this code that are retained.
+
+===============================================================================
+*/
+
+#include <asm/div64.h>
+
+#include "fpa11.h"
+//#include "milieu.h"
+//#include "softfloat.h"
+
+/*
+-------------------------------------------------------------------------------
+Primitive arithmetic functions, including multi-word arithmetic, and
+division and square root approximations.  (Can be specialized to target if
+desired.)
+-------------------------------------------------------------------------------
+*/
+#include "softfloat-macros"
+
+/*
+-------------------------------------------------------------------------------
+Functions and definitions to determine:  (1) whether tininess for underflow
+is detected before or after rounding by default, (2) what (if anything)
+happens when exceptions are raised, (3) how signaling NaNs are distinguished
+from quiet NaNs, (4) the default generated quiet NaNs, and (5) how NaNs
+are propagated from function inputs to output.  These details are target-
+specific.
+-------------------------------------------------------------------------------
+*/
+#include "softfloat-specialize"
+
+/*
+-------------------------------------------------------------------------------
+Takes a 64-bit fixed-point value `absZ' with binary point between bits 6
+and 7, and returns the properly rounded 32-bit integer corresponding to the
+input.  If `zSign' is nonzero, the input is negated before being converted
+to an integer.  Bit 63 of `absZ' must be zero.  Ordinarily, the fixed-point
+input is simply rounded to an integer, with the inexact exception raised if
+the input cannot be represented exactly as an integer.  If the fixed-point
+input is too large, however, the invalid exception is raised and the largest
+positive or negative integer is returned.
+-------------------------------------------------------------------------------
+*/
+static int32 roundAndPackInt32( struct roundingData *roundData, flag zSign, bits64 absZ )
+{
+    int8 roundingMode;
+    flag roundNearestEven;
+    int8 roundIncrement, roundBits;
+    int32 z;
+
+    roundingMode = roundData->mode;
+    roundNearestEven = ( roundingMode == float_round_nearest_even );
+    roundIncrement = 0x40;
+    if ( ! roundNearestEven ) {
+        if ( roundingMode == float_round_to_zero ) {
+            roundIncrement = 0;
+        }
+        else {
+            roundIncrement = 0x7F;
+            if ( zSign ) {
+                if ( roundingMode == float_round_up ) roundIncrement = 0;
+            }
+            else {
+                if ( roundingMode == float_round_down ) roundIncrement = 0;
+            }
+        }
+    }
+    roundBits = absZ & 0x7F;
+    absZ = ( absZ + roundIncrement )>>7;
+    absZ &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
+    z = absZ;
+    if ( zSign ) z = - z;
+    if ( ( absZ>>32 ) || ( z && ( ( z < 0 ) ^ zSign ) ) ) {
+        roundData->exception |= float_flag_invalid;
+        return zSign ? 0x80000000 : 0x7FFFFFFF;
+    }
+    if ( roundBits ) roundData->exception |= float_flag_inexact;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the fraction bits of the single-precision floating-point value `a'.
+-------------------------------------------------------------------------------
+*/
+INLINE bits32 extractFloat32Frac( float32 a )
+{
+
+    return a & 0x007FFFFF;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the exponent bits of the single-precision floating-point value `a'.
+-------------------------------------------------------------------------------
+*/
+INLINE int16 extractFloat32Exp( float32 a )
+{
+
+    return ( a>>23 ) & 0xFF;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the sign bit of the single-precision floating-point value `a'.
+-------------------------------------------------------------------------------
+*/
+#if 0	/* in softfloat.h */
+INLINE flag extractFloat32Sign( float32 a )
+{
+
+    return a>>31;
+
+}
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Normalizes the subnormal single-precision floating-point value represented
+by the denormalized significand `aSig'.  The normalized exponent and
+significand are stored at the locations pointed to by `zExpPtr' and
+`zSigPtr', respectively.
+-------------------------------------------------------------------------------
+*/
+static void
+ normalizeFloat32Subnormal( bits32 aSig, int16 *zExpPtr, bits32 *zSigPtr )
+{
+    int8 shiftCount;
+
+    shiftCount = countLeadingZeros32( aSig ) - 8;
+    *zSigPtr = aSig<<shiftCount;
+    *zExpPtr = 1 - shiftCount;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
+single-precision floating-point value, returning the result.  After being
+shifted into the proper positions, the three fields are simply added
+together to form the result.  This means that any integer portion of `zSig'
+will be added into the exponent.  Since a properly normalized significand
+will have an integer portion equal to 1, the `zExp' input should be 1 less
+than the desired result exponent whenever `zSig' is a complete, normalized
+significand.
+-------------------------------------------------------------------------------
+*/
+INLINE float32 packFloat32( flag zSign, int16 zExp, bits32 zSig )
+{
+#if 0
+   float32 f;
+   __asm__("@ packFloat32				\n\
+   	    mov %0, %1, asl #31				\n\
+   	    orr %0, %2, asl #23				\n\
+   	    orr %0, %3"
+   	    : /* no outputs */
+   	    : "g" (f), "g" (zSign), "g" (zExp), "g" (zSig)
+   	    : "cc");
+   return f;
+#else
+    return ( ( (bits32) zSign )<<31 ) + ( ( (bits32) zExp )<<23 ) + zSig;
+#endif 
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes an abstract floating-point value having sign `zSign', exponent `zExp',
+and significand `zSig', and returns the proper single-precision floating-
+point value corresponding to the abstract input.  Ordinarily, the abstract
+value is simply rounded and packed into the single-precision format, with
+the inexact exception raised if the abstract input cannot be represented
+exactly.  If the abstract value is too large, however, the overflow and
+inexact exceptions are raised and an infinity or maximal finite value is
+returned.  If the abstract value is too small, the input value is rounded to
+a subnormal number, and the underflow and inexact exceptions are raised if
+the abstract input cannot be represented exactly as a subnormal single-
+precision floating-point number.
+    The input significand `zSig' has its binary point between bits 30
+and 29, which is 7 bits to the left of the usual location.  This shifted
+significand must be normalized or smaller.  If `zSig' is not normalized,
+`zExp' must be 0; in that case, the result returned is a subnormal number,
+and it must not require rounding.  In the usual case that `zSig' is
+normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
+The handling of underflow and overflow follows the IEC/IEEE Standard for
+Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static float32 roundAndPackFloat32( struct roundingData *roundData, flag zSign, int16 zExp, bits32 zSig )
+{
+    int8 roundingMode;
+    flag roundNearestEven;
+    int8 roundIncrement, roundBits;
+    flag isTiny;
+
+    roundingMode = roundData->mode;
+    roundNearestEven = ( roundingMode == float_round_nearest_even );
+    roundIncrement = 0x40;
+    if ( ! roundNearestEven ) {
+        if ( roundingMode == float_round_to_zero ) {
+            roundIncrement = 0;
+        }
+        else {
+            roundIncrement = 0x7F;
+            if ( zSign ) {
+                if ( roundingMode == float_round_up ) roundIncrement = 0;
+            }
+            else {
+                if ( roundingMode == float_round_down ) roundIncrement = 0;
+            }
+        }
+    }
+    roundBits = zSig & 0x7F;
+    if ( 0xFD <= (bits16) zExp ) {
+        if (    ( 0xFD < zExp )
+             || (    ( zExp == 0xFD )
+                  && ( (sbits32) ( zSig + roundIncrement ) < 0 ) )
+           ) {
+            roundData->exception |= float_flag_overflow | float_flag_inexact;
+            return packFloat32( zSign, 0xFF, 0 ) - ( roundIncrement == 0 );
+        }
+        if ( zExp < 0 ) {
+            isTiny =
+                   ( float_detect_tininess == float_tininess_before_rounding )
+                || ( zExp < -1 )
+                || ( zSig + roundIncrement < 0x80000000 );
+            shift32RightJamming( zSig, - zExp, &zSig );
+            zExp = 0;
+            roundBits = zSig & 0x7F;
+            if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow;
+        }
+    }
+    if ( roundBits ) roundData->exception |= float_flag_inexact;
+    zSig = ( zSig + roundIncrement )>>7;
+    zSig &= ~ ( ( ( roundBits ^ 0x40 ) == 0 ) & roundNearestEven );
+    if ( zSig == 0 ) zExp = 0;
+    return packFloat32( zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes an abstract floating-point value having sign `zSign', exponent `zExp',
+and significand `zSig', and returns the proper single-precision floating-
+point value corresponding to the abstract input.  This routine is just like
+`roundAndPackFloat32' except that `zSig' does not have to be normalized in
+any way.  In all cases, `zExp' must be 1 less than the ``true'' floating-
+point exponent.
+-------------------------------------------------------------------------------
+*/
+static float32
+ normalizeRoundAndPackFloat32( struct roundingData *roundData, flag zSign, int16 zExp, bits32 zSig )
+{
+    int8 shiftCount;
+
+    shiftCount = countLeadingZeros32( zSig ) - 1;
+    return roundAndPackFloat32( roundData, zSign, zExp - shiftCount, zSig<<shiftCount );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the fraction bits of the double-precision floating-point value `a'.
+-------------------------------------------------------------------------------
+*/
+INLINE bits64 extractFloat64Frac( float64 a )
+{
+
+    return a & LIT64( 0x000FFFFFFFFFFFFF );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the exponent bits of the double-precision floating-point value `a'.
+-------------------------------------------------------------------------------
+*/
+INLINE int16 extractFloat64Exp( float64 a )
+{
+
+    return ( a>>52 ) & 0x7FF;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the sign bit of the double-precision floating-point value `a'.
+-------------------------------------------------------------------------------
+*/
+#if 0	/* in softfloat.h */
+INLINE flag extractFloat64Sign( float64 a )
+{
+
+    return a>>63;
+
+}
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Normalizes the subnormal double-precision floating-point value represented
+by the denormalized significand `aSig'.  The normalized exponent and
+significand are stored at the locations pointed to by `zExpPtr' and
+`zSigPtr', respectively.
+-------------------------------------------------------------------------------
+*/
+static void
+ normalizeFloat64Subnormal( bits64 aSig, int16 *zExpPtr, bits64 *zSigPtr )
+{
+    int8 shiftCount;
+
+    shiftCount = countLeadingZeros64( aSig ) - 11;
+    *zSigPtr = aSig<<shiftCount;
+    *zExpPtr = 1 - shiftCount;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Packs the sign `zSign', exponent `zExp', and significand `zSig' into a
+double-precision floating-point value, returning the result.  After being
+shifted into the proper positions, the three fields are simply added
+together to form the result.  This means that any integer portion of `zSig'
+will be added into the exponent.  Since a properly normalized significand
+will have an integer portion equal to 1, the `zExp' input should be 1 less
+than the desired result exponent whenever `zSig' is a complete, normalized
+significand.
+-------------------------------------------------------------------------------
+*/
+INLINE float64 packFloat64( flag zSign, int16 zExp, bits64 zSig )
+{
+
+    return ( ( (bits64) zSign )<<63 ) + ( ( (bits64) zExp )<<52 ) + zSig;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes an abstract floating-point value having sign `zSign', exponent `zExp',
+and significand `zSig', and returns the proper double-precision floating-
+point value corresponding to the abstract input.  Ordinarily, the abstract
+value is simply rounded and packed into the double-precision format, with
+the inexact exception raised if the abstract input cannot be represented
+exactly.  If the abstract value is too large, however, the overflow and
+inexact exceptions are raised and an infinity or maximal finite value is
+returned.  If the abstract value is too small, the input value is rounded to
+a subnormal number, and the underflow and inexact exceptions are raised if
+the abstract input cannot be represented exactly as a subnormal double-
+precision floating-point number.
+    The input significand `zSig' has its binary point between bits 62
+and 61, which is 10 bits to the left of the usual location.  This shifted
+significand must be normalized or smaller.  If `zSig' is not normalized,
+`zExp' must be 0; in that case, the result returned is a subnormal number,
+and it must not require rounding.  In the usual case that `zSig' is
+normalized, `zExp' must be 1 less than the ``true'' floating-point exponent.
+The handling of underflow and overflow follows the IEC/IEEE Standard for
+Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static float64 roundAndPackFloat64( struct roundingData *roundData, flag zSign, int16 zExp, bits64 zSig )
+{
+    int8 roundingMode;
+    flag roundNearestEven;
+    int16 roundIncrement, roundBits;
+    flag isTiny;
+
+    roundingMode = roundData->mode;
+    roundNearestEven = ( roundingMode == float_round_nearest_even );
+    roundIncrement = 0x200;
+    if ( ! roundNearestEven ) {
+        if ( roundingMode == float_round_to_zero ) {
+            roundIncrement = 0;
+        }
+        else {
+            roundIncrement = 0x3FF;
+            if ( zSign ) {
+                if ( roundingMode == float_round_up ) roundIncrement = 0;
+            }
+            else {
+                if ( roundingMode == float_round_down ) roundIncrement = 0;
+            }
+        }
+    }
+    roundBits = zSig & 0x3FF;
+    if ( 0x7FD <= (bits16) zExp ) {
+        if (    ( 0x7FD < zExp )
+             || (    ( zExp == 0x7FD )
+                  && ( (sbits64) ( zSig + roundIncrement ) < 0 ) )
+           ) {
+            //register int lr = __builtin_return_address(0);
+            //printk("roundAndPackFloat64 called from 0x%08x\n",lr);
+            roundData->exception |= float_flag_overflow | float_flag_inexact;
+            return packFloat64( zSign, 0x7FF, 0 ) - ( roundIncrement == 0 );
+        }
+        if ( zExp < 0 ) {
+            isTiny =
+                   ( float_detect_tininess == float_tininess_before_rounding )
+                || ( zExp < -1 )
+                || ( zSig + roundIncrement < LIT64( 0x8000000000000000 ) );
+            shift64RightJamming( zSig, - zExp, &zSig );
+            zExp = 0;
+            roundBits = zSig & 0x3FF;
+            if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow;
+        }
+    }
+    if ( roundBits ) roundData->exception |= float_flag_inexact;
+    zSig = ( zSig + roundIncrement )>>10;
+    zSig &= ~ ( ( ( roundBits ^ 0x200 ) == 0 ) & roundNearestEven );
+    if ( zSig == 0 ) zExp = 0;
+    return packFloat64( zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes an abstract floating-point value having sign `zSign', exponent `zExp',
+and significand `zSig', and returns the proper double-precision floating-
+point value corresponding to the abstract input.  This routine is just like
+`roundAndPackFloat64' except that `zSig' does not have to be normalized in
+any way.  In all cases, `zExp' must be 1 less than the ``true'' floating-
+point exponent.
+-------------------------------------------------------------------------------
+*/
+static float64
+ normalizeRoundAndPackFloat64( struct roundingData *roundData, flag zSign, int16 zExp, bits64 zSig )
+{
+    int8 shiftCount;
+
+    shiftCount = countLeadingZeros64( zSig ) - 1;
+    return roundAndPackFloat64( roundData, zSign, zExp - shiftCount, zSig<<shiftCount );
+
+}
+
+#ifdef FLOATX80
+
+/*
+-------------------------------------------------------------------------------
+Returns the fraction bits of the extended double-precision floating-point
+value `a'.
+-------------------------------------------------------------------------------
+*/
+INLINE bits64 extractFloatx80Frac( floatx80 a )
+{
+
+    return a.low;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the exponent bits of the extended double-precision floating-point
+value `a'.
+-------------------------------------------------------------------------------
+*/
+INLINE int32 extractFloatx80Exp( floatx80 a )
+{
+
+    return a.high & 0x7FFF;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the sign bit of the extended double-precision floating-point value
+`a'.
+-------------------------------------------------------------------------------
+*/
+INLINE flag extractFloatx80Sign( floatx80 a )
+{
+
+    return a.high>>15;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Normalizes the subnormal extended double-precision floating-point value
+represented by the denormalized significand `aSig'.  The normalized exponent
+and significand are stored at the locations pointed to by `zExpPtr' and
+`zSigPtr', respectively.
+-------------------------------------------------------------------------------
+*/
+static void
+ normalizeFloatx80Subnormal( bits64 aSig, int32 *zExpPtr, bits64 *zSigPtr )
+{
+    int8 shiftCount;
+
+    shiftCount = countLeadingZeros64( aSig );
+    *zSigPtr = aSig<<shiftCount;
+    *zExpPtr = 1 - shiftCount;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Packs the sign `zSign', exponent `zExp', and significand `zSig' into an
+extended double-precision floating-point value, returning the result.
+-------------------------------------------------------------------------------
+*/
+INLINE floatx80 packFloatx80( flag zSign, int32 zExp, bits64 zSig )
+{
+    floatx80 z;
+
+    z.low = zSig;
+    z.high = ( ( (bits16) zSign )<<15 ) + zExp;
+    z.__padding = 0;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes an abstract floating-point value having sign `zSign', exponent `zExp',
+and extended significand formed by the concatenation of `zSig0' and `zSig1',
+and returns the proper extended double-precision floating-point value
+corresponding to the abstract input.  Ordinarily, the abstract value is
+rounded and packed into the extended double-precision format, with the
+inexact exception raised if the abstract input cannot be represented
+exactly.  If the abstract value is too large, however, the overflow and
+inexact exceptions are raised and an infinity or maximal finite value is
+returned.  If the abstract value is too small, the input value is rounded to
+a subnormal number, and the underflow and inexact exceptions are raised if
+the abstract input cannot be represented exactly as a subnormal extended
+double-precision floating-point number.
+    If `roundingPrecision' is 32 or 64, the result is rounded to the same
+number of bits as single or double precision, respectively.  Otherwise, the
+result is rounded to the full precision of the extended double-precision
+format.
+    The input significand must be normalized or smaller.  If the input
+significand is not normalized, `zExp' must be 0; in that case, the result
+returned is a subnormal number, and it must not require rounding.  The
+handling of underflow and overflow follows the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static floatx80
+ roundAndPackFloatx80(
+     struct roundingData *roundData, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
+ )
+{
+    int8 roundingMode, roundingPrecision;
+    flag roundNearestEven, increment, isTiny;
+    int64 roundIncrement, roundMask, roundBits;
+
+    roundingMode = roundData->mode;
+    roundingPrecision = roundData->precision;
+    roundNearestEven = ( roundingMode == float_round_nearest_even );
+    if ( roundingPrecision == 80 ) goto precision80;
+    if ( roundingPrecision == 64 ) {
+        roundIncrement = LIT64( 0x0000000000000400 );
+        roundMask = LIT64( 0x00000000000007FF );
+    }
+    else if ( roundingPrecision == 32 ) {
+        roundIncrement = LIT64( 0x0000008000000000 );
+        roundMask = LIT64( 0x000000FFFFFFFFFF );
+    }
+    else {
+        goto precision80;
+    }
+    zSig0 |= ( zSig1 != 0 );
+    if ( ! roundNearestEven ) {
+        if ( roundingMode == float_round_to_zero ) {
+            roundIncrement = 0;
+        }
+        else {
+            roundIncrement = roundMask;
+            if ( zSign ) {
+                if ( roundingMode == float_round_up ) roundIncrement = 0;
+            }
+            else {
+                if ( roundingMode == float_round_down ) roundIncrement = 0;
+            }
+        }
+    }
+    roundBits = zSig0 & roundMask;
+    if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
+        if (    ( 0x7FFE < zExp )
+             || ( ( zExp == 0x7FFE ) && ( zSig0 + roundIncrement < zSig0 ) )
+           ) {
+            goto overflow;
+        }
+        if ( zExp <= 0 ) {
+            isTiny =
+                   ( float_detect_tininess == float_tininess_before_rounding )
+                || ( zExp < 0 )
+                || ( zSig0 <= zSig0 + roundIncrement );
+            shift64RightJamming( zSig0, 1 - zExp, &zSig0 );
+            zExp = 0;
+            roundBits = zSig0 & roundMask;
+            if ( isTiny && roundBits ) roundData->exception |= float_flag_underflow;
+            if ( roundBits ) roundData->exception |= float_flag_inexact;
+            zSig0 += roundIncrement;
+            if ( (sbits64) zSig0 < 0 ) zExp = 1;
+            roundIncrement = roundMask + 1;
+            if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
+                roundMask |= roundIncrement;
+            }
+            zSig0 &= ~ roundMask;
+            return packFloatx80( zSign, zExp, zSig0 );
+        }
+    }
+    if ( roundBits ) roundData->exception |= float_flag_inexact;
+    zSig0 += roundIncrement;
+    if ( zSig0 < roundIncrement ) {
+        ++zExp;
+        zSig0 = LIT64( 0x8000000000000000 );
+    }
+    roundIncrement = roundMask + 1;
+    if ( roundNearestEven && ( roundBits<<1 == roundIncrement ) ) {
+        roundMask |= roundIncrement;
+    }
+    zSig0 &= ~ roundMask;
+    if ( zSig0 == 0 ) zExp = 0;
+    return packFloatx80( zSign, zExp, zSig0 );
+ precision80:
+    increment = ( (sbits64) zSig1 < 0 );
+    if ( ! roundNearestEven ) {
+        if ( roundingMode == float_round_to_zero ) {
+            increment = 0;
+        }
+        else {
+            if ( zSign ) {
+                increment = ( roundingMode == float_round_down ) && zSig1;
+            }
+            else {
+                increment = ( roundingMode == float_round_up ) && zSig1;
+            }
+        }
+    }
+    if ( 0x7FFD <= (bits32) ( zExp - 1 ) ) {
+        if (    ( 0x7FFE < zExp )
+             || (    ( zExp == 0x7FFE )
+                  && ( zSig0 == LIT64( 0xFFFFFFFFFFFFFFFF ) )
+                  && increment
+                )
+           ) {
+            roundMask = 0;
+ overflow:
+            roundData->exception |= float_flag_overflow | float_flag_inexact;
+            if (    ( roundingMode == float_round_to_zero )
+                 || ( zSign && ( roundingMode == float_round_up ) )
+                 || ( ! zSign && ( roundingMode == float_round_down ) )
+               ) {
+                return packFloatx80( zSign, 0x7FFE, ~ roundMask );
+            }
+            return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
+        }
+        if ( zExp <= 0 ) {
+            isTiny =
+                   ( float_detect_tininess == float_tininess_before_rounding )
+                || ( zExp < 0 )
+                || ! increment
+                || ( zSig0 < LIT64( 0xFFFFFFFFFFFFFFFF ) );
+            shift64ExtraRightJamming( zSig0, zSig1, 1 - zExp, &zSig0, &zSig1 );
+            zExp = 0;
+            if ( isTiny && zSig1 ) roundData->exception |= float_flag_underflow;
+            if ( zSig1 ) roundData->exception |= float_flag_inexact;
+            if ( roundNearestEven ) {
+                increment = ( (sbits64) zSig1 < 0 );
+            }
+            else {
+                if ( zSign ) {
+                    increment = ( roundingMode == float_round_down ) && zSig1;
+                }
+                else {
+                    increment = ( roundingMode == float_round_up ) && zSig1;
+                }
+            }
+            if ( increment ) {
+                ++zSig0;
+                zSig0 &= ~ ( ( zSig1 + zSig1 == 0 ) & roundNearestEven );
+                if ( (sbits64) zSig0 < 0 ) zExp = 1;
+            }
+            return packFloatx80( zSign, zExp, zSig0 );
+        }
+    }
+    if ( zSig1 ) roundData->exception |= float_flag_inexact;
+    if ( increment ) {
+        ++zSig0;
+        if ( zSig0 == 0 ) {
+            ++zExp;
+            zSig0 = LIT64( 0x8000000000000000 );
+        }
+        else {
+            zSig0 &= ~ ( ( zSig1 + zSig1 == 0 ) & roundNearestEven );
+        }
+    }
+    else {
+        if ( zSig0 == 0 ) zExp = 0;
+    }
+    
+    return packFloatx80( zSign, zExp, zSig0 );
+}
+
+/*
+-------------------------------------------------------------------------------
+Takes an abstract floating-point value having sign `zSign', exponent
+`zExp', and significand formed by the concatenation of `zSig0' and `zSig1',
+and returns the proper extended double-precision floating-point value
+corresponding to the abstract input.  This routine is just like
+`roundAndPackFloatx80' except that the input significand does not have to be
+normalized.
+-------------------------------------------------------------------------------
+*/
+static floatx80
+ normalizeRoundAndPackFloatx80(
+     struct roundingData *roundData, flag zSign, int32 zExp, bits64 zSig0, bits64 zSig1
+ )
+{
+    int8 shiftCount;
+
+    if ( zSig0 == 0 ) {
+        zSig0 = zSig1;
+        zSig1 = 0;
+        zExp -= 64;
+    }
+    shiftCount = countLeadingZeros64( zSig0 );
+    shortShift128Left( zSig0, zSig1, shiftCount, &zSig0, &zSig1 );
+    zExp -= shiftCount;
+    return
+        roundAndPackFloatx80( roundData, zSign, zExp, zSig0, zSig1 );
+
+}
+
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the 32-bit two's complement integer `a' to
+the single-precision floating-point format.  The conversion is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 int32_to_float32(struct roundingData *roundData, int32 a)
+{
+    flag zSign;
+
+    if ( a == 0 ) return 0;
+    if ( a == 0x80000000 ) return packFloat32( 1, 0x9E, 0 );
+    zSign = ( a < 0 );
+    return normalizeRoundAndPackFloat32( roundData, zSign, 0x9C, zSign ? - a : a );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the 32-bit two's complement integer `a' to
+the double-precision floating-point format.  The conversion is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 int32_to_float64( int32 a )
+{
+    flag aSign;
+    uint32 absA;
+    int8 shiftCount;
+    bits64 zSig;
+
+    if ( a == 0 ) return 0;
+    aSign = ( a < 0 );
+    absA = aSign ? - a : a;
+    shiftCount = countLeadingZeros32( absA ) + 21;
+    zSig = absA;
+    return packFloat64( aSign, 0x432 - shiftCount, zSig<<shiftCount );
+
+}
+
+#ifdef FLOATX80
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the 32-bit two's complement integer `a'
+to the extended double-precision floating-point format.  The conversion
+is performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 int32_to_floatx80( int32 a )
+{
+    flag zSign;
+    uint32 absA;
+    int8 shiftCount;
+    bits64 zSig;
+
+    if ( a == 0 ) return packFloatx80( 0, 0, 0 );
+    zSign = ( a < 0 );
+    absA = zSign ? - a : a;
+    shiftCount = countLeadingZeros32( absA ) + 32;
+    zSig = absA;
+    return packFloatx80( zSign, 0x403E - shiftCount, zSig<<shiftCount );
+
+}
+
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the single-precision floating-point value
+`a' to the 32-bit two's complement integer format.  The conversion is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic---which means in particular that the conversion is rounded
+according to the current rounding mode.  If `a' is a NaN, the largest
+positive integer is returned.  Otherwise, if the conversion overflows, the
+largest integer with the same sign as `a' is returned.
+-------------------------------------------------------------------------------
+*/
+int32 float32_to_int32( struct roundingData *roundData, float32 a )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    bits32 aSig;
+    bits64 zSig;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
+    if ( aExp ) aSig |= 0x00800000;
+    shiftCount = 0xAF - aExp;
+    zSig = aSig;
+    zSig <<= 32;
+    if ( 0 < shiftCount ) shift64RightJamming( zSig, shiftCount, &zSig );
+    return roundAndPackInt32( roundData, aSign, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the single-precision floating-point value
+`a' to the 32-bit two's complement integer format.  The conversion is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic, except that the conversion is always rounded toward zero.  If
+`a' is a NaN, the largest positive integer is returned.  Otherwise, if the
+conversion overflows, the largest integer with the same sign as `a' is
+returned.
+-------------------------------------------------------------------------------
+*/
+int32 float32_to_int32_round_to_zero( float32 a )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    bits32 aSig;
+    int32 z;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    shiftCount = aExp - 0x9E;
+    if ( 0 <= shiftCount ) {
+        if ( a == 0xCF000000 ) return 0x80000000;
+        float_raise( float_flag_invalid );
+        if ( ! aSign || ( ( aExp == 0xFF ) && aSig ) ) return 0x7FFFFFFF;
+        return 0x80000000;
+    }
+    else if ( aExp <= 0x7E ) {
+        if ( aExp | aSig ) float_raise( float_flag_inexact );
+        return 0;
+    }
+    aSig = ( aSig | 0x00800000 )<<8;
+    z = aSig>>( - shiftCount );
+    if ( (bits32) ( aSig<<( shiftCount & 31 ) ) ) {
+        float_raise( float_flag_inexact );
+    }
+    return aSign ? - z : z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the single-precision floating-point value
+`a' to the double-precision floating-point format.  The conversion is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 float32_to_float64( float32 a )
+{
+    flag aSign;
+    int16 aExp;
+    bits32 aSig;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    if ( aExp == 0xFF ) {
+        if ( aSig ) return commonNaNToFloat64( float32ToCommonNaN( a ) );
+        return packFloat64( aSign, 0x7FF, 0 );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat64( aSign, 0, 0 );
+        normalizeFloat32Subnormal( aSig, &aExp, &aSig );
+        --aExp;
+    }
+    return packFloat64( aSign, aExp + 0x380, ( (bits64) aSig )<<29 );
+
+}
+
+#ifdef FLOATX80
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the single-precision floating-point value
+`a' to the extended double-precision floating-point format.  The conversion
+is performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 float32_to_floatx80( float32 a )
+{
+    flag aSign;
+    int16 aExp;
+    bits32 aSig;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    if ( aExp == 0xFF ) {
+        if ( aSig ) return commonNaNToFloatx80( float32ToCommonNaN( a ) );
+        return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
+        normalizeFloat32Subnormal( aSig, &aExp, &aSig );
+    }
+    aSig |= 0x00800000;
+    return packFloatx80( aSign, aExp + 0x3F80, ( (bits64) aSig )<<40 );
+
+}
+
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Rounds the single-precision floating-point value `a' to an integer, and
+returns the result as a single-precision floating-point value.  The
+operation is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 float32_round_to_int( struct roundingData *roundData, float32 a )
+{
+    flag aSign;
+    int16 aExp;
+    bits32 lastBitMask, roundBitsMask;
+    int8 roundingMode;
+    float32 z;
+
+    aExp = extractFloat32Exp( a );
+    if ( 0x96 <= aExp ) {
+        if ( ( aExp == 0xFF ) && extractFloat32Frac( a ) ) {
+            return propagateFloat32NaN( a, a );
+        }
+        return a;
+    }
+    roundingMode = roundData->mode;
+    if ( aExp <= 0x7E ) {
+        if ( (bits32) ( a<<1 ) == 0 ) return a;
+        roundData->exception |= float_flag_inexact;
+        aSign = extractFloat32Sign( a );
+        switch ( roundingMode ) {
+         case float_round_nearest_even:
+            if ( ( aExp == 0x7E ) && extractFloat32Frac( a ) ) {
+                return packFloat32( aSign, 0x7F, 0 );
+            }
+            break;
+         case float_round_down:
+            return aSign ? 0xBF800000 : 0;
+         case float_round_up:
+            return aSign ? 0x80000000 : 0x3F800000;
+        }
+        return packFloat32( aSign, 0, 0 );
+    }
+    lastBitMask = 1;
+    lastBitMask <<= 0x96 - aExp;
+    roundBitsMask = lastBitMask - 1;
+    z = a;
+    if ( roundingMode == float_round_nearest_even ) {
+        z += lastBitMask>>1;
+        if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
+    }
+    else if ( roundingMode != float_round_to_zero ) {
+        if ( extractFloat32Sign( z ) ^ ( roundingMode == float_round_up ) ) {
+            z += roundBitsMask;
+        }
+    }
+    z &= ~ roundBitsMask;
+    if ( z != a ) roundData->exception |= float_flag_inexact;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of adding the absolute values of the single-precision
+floating-point values `a' and `b'.  If `zSign' is true, the sum is negated
+before being returned.  `zSign' is ignored if the result is a NaN.  The
+addition is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static float32 addFloat32Sigs( struct roundingData *roundData, float32 a, float32 b, flag zSign )
+{
+    int16 aExp, bExp, zExp;
+    bits32 aSig, bSig, zSig;
+    int16 expDiff;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    bSig = extractFloat32Frac( b );
+    bExp = extractFloat32Exp( b );
+    expDiff = aExp - bExp;
+    aSig <<= 6;
+    bSig <<= 6;
+    if ( 0 < expDiff ) {
+        if ( aExp == 0xFF ) {
+            if ( aSig ) return propagateFloat32NaN( a, b );
+            return a;
+        }
+        if ( bExp == 0 ) {
+            --expDiff;
+        }
+        else {
+            bSig |= 0x20000000;
+        }
+        shift32RightJamming( bSig, expDiff, &bSig );
+        zExp = aExp;
+    }
+    else if ( expDiff < 0 ) {
+        if ( bExp == 0xFF ) {
+            if ( bSig ) return propagateFloat32NaN( a, b );
+            return packFloat32( zSign, 0xFF, 0 );
+        }
+        if ( aExp == 0 ) {
+            ++expDiff;
+        }
+        else {
+            aSig |= 0x20000000;
+        }
+        shift32RightJamming( aSig, - expDiff, &aSig );
+        zExp = bExp;
+    }
+    else {
+        if ( aExp == 0xFF ) {
+            if ( aSig | bSig ) return propagateFloat32NaN( a, b );
+            return a;
+        }
+        if ( aExp == 0 ) return packFloat32( zSign, 0, ( aSig + bSig )>>6 );
+        zSig = 0x40000000 + aSig + bSig;
+        zExp = aExp;
+        goto roundAndPack;
+    }
+    aSig |= 0x20000000;
+    zSig = ( aSig + bSig )<<1;
+    --zExp;
+    if ( (sbits32) zSig < 0 ) {
+        zSig = aSig + bSig;
+        ++zExp;
+    }
+ roundAndPack:
+    return roundAndPackFloat32( roundData, zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of subtracting the absolute values of the single-
+precision floating-point values `a' and `b'.  If `zSign' is true, the
+difference is negated before being returned.  `zSign' is ignored if the
+result is a NaN.  The subtraction is performed according to the IEC/IEEE
+Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static float32 subFloat32Sigs( struct roundingData *roundData, float32 a, float32 b, flag zSign )
+{
+    int16 aExp, bExp, zExp;
+    bits32 aSig, bSig, zSig;
+    int16 expDiff;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    bSig = extractFloat32Frac( b );
+    bExp = extractFloat32Exp( b );
+    expDiff = aExp - bExp;
+    aSig <<= 7;
+    bSig <<= 7;
+    if ( 0 < expDiff ) goto aExpBigger;
+    if ( expDiff < 0 ) goto bExpBigger;
+    if ( aExp == 0xFF ) {
+        if ( aSig | bSig ) return propagateFloat32NaN( a, b );
+        roundData->exception |= float_flag_invalid;
+        return float32_default_nan;
+    }
+    if ( aExp == 0 ) {
+        aExp = 1;
+        bExp = 1;
+    }
+    if ( bSig < aSig ) goto aBigger;
+    if ( aSig < bSig ) goto bBigger;
+    return packFloat32( roundData->mode == float_round_down, 0, 0 );
+ bExpBigger:
+    if ( bExp == 0xFF ) {
+        if ( bSig ) return propagateFloat32NaN( a, b );
+        return packFloat32( zSign ^ 1, 0xFF, 0 );
+    }
+    if ( aExp == 0 ) {
+        ++expDiff;
+    }
+    else {
+        aSig |= 0x40000000;
+    }
+    shift32RightJamming( aSig, - expDiff, &aSig );
+    bSig |= 0x40000000;
+ bBigger:
+    zSig = bSig - aSig;
+    zExp = bExp;
+    zSign ^= 1;
+    goto normalizeRoundAndPack;
+ aExpBigger:
+    if ( aExp == 0xFF ) {
+        if ( aSig ) return propagateFloat32NaN( a, b );
+        return a;
+    }
+    if ( bExp == 0 ) {
+        --expDiff;
+    }
+    else {
+        bSig |= 0x40000000;
+    }
+    shift32RightJamming( bSig, expDiff, &bSig );
+    aSig |= 0x40000000;
+ aBigger:
+    zSig = aSig - bSig;
+    zExp = aExp;
+ normalizeRoundAndPack:
+    --zExp;
+    return normalizeRoundAndPackFloat32( roundData, zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of adding the single-precision floating-point values `a'
+and `b'.  The operation is performed according to the IEC/IEEE Standard for
+Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 float32_add( struct roundingData *roundData, float32 a, float32 b )
+{
+    flag aSign, bSign;
+
+    aSign = extractFloat32Sign( a );
+    bSign = extractFloat32Sign( b );
+    if ( aSign == bSign ) {
+        return addFloat32Sigs( roundData, a, b, aSign );
+    }
+    else {
+        return subFloat32Sigs( roundData, a, b, aSign );
+    }
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of subtracting the single-precision floating-point values
+`a' and `b'.  The operation is performed according to the IEC/IEEE Standard
+for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 float32_sub( struct roundingData *roundData, float32 a, float32 b )
+{
+    flag aSign, bSign;
+
+    aSign = extractFloat32Sign( a );
+    bSign = extractFloat32Sign( b );
+    if ( aSign == bSign ) {
+        return subFloat32Sigs( roundData, a, b, aSign );
+    }
+    else {
+        return addFloat32Sigs( roundData, a, b, aSign );
+    }
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of multiplying the single-precision floating-point values
+`a' and `b'.  The operation is performed according to the IEC/IEEE Standard
+for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 float32_mul( struct roundingData *roundData, float32 a, float32 b )
+{
+    flag aSign, bSign, zSign;
+    int16 aExp, bExp, zExp;
+    bits32 aSig, bSig;
+    bits64 zSig64;
+    bits32 zSig;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    bSig = extractFloat32Frac( b );
+    bExp = extractFloat32Exp( b );
+    bSign = extractFloat32Sign( b );
+    zSign = aSign ^ bSign;
+    if ( aExp == 0xFF ) {
+        if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
+            return propagateFloat32NaN( a, b );
+        }
+        if ( ( bExp | bSig ) == 0 ) {
+            roundData->exception |= float_flag_invalid;
+            return float32_default_nan;
+        }
+        return packFloat32( zSign, 0xFF, 0 );
+    }
+    if ( bExp == 0xFF ) {
+        if ( bSig ) return propagateFloat32NaN( a, b );
+        if ( ( aExp | aSig ) == 0 ) {
+            roundData->exception |= float_flag_invalid;
+            return float32_default_nan;
+        }
+        return packFloat32( zSign, 0xFF, 0 );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
+        normalizeFloat32Subnormal( aSig, &aExp, &aSig );
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) return packFloat32( zSign, 0, 0 );
+        normalizeFloat32Subnormal( bSig, &bExp, &bSig );
+    }
+    zExp = aExp + bExp - 0x7F;
+    aSig = ( aSig | 0x00800000 )<<7;
+    bSig = ( bSig | 0x00800000 )<<8;
+    shift64RightJamming( ( (bits64) aSig ) * bSig, 32, &zSig64 );
+    zSig = zSig64;
+    if ( 0 <= (sbits32) ( zSig<<1 ) ) {
+        zSig <<= 1;
+        --zExp;
+    }
+    return roundAndPackFloat32( roundData, zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of dividing the single-precision floating-point value `a'
+by the corresponding value `b'.  The operation is performed according to the
+IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 float32_div( struct roundingData *roundData, float32 a, float32 b )
+{
+    flag aSign, bSign, zSign;
+    int16 aExp, bExp, zExp;
+    bits32 aSig, bSig, zSig;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    bSig = extractFloat32Frac( b );
+    bExp = extractFloat32Exp( b );
+    bSign = extractFloat32Sign( b );
+    zSign = aSign ^ bSign;
+    if ( aExp == 0xFF ) {
+        if ( aSig ) return propagateFloat32NaN( a, b );
+        if ( bExp == 0xFF ) {
+            if ( bSig ) return propagateFloat32NaN( a, b );
+            roundData->exception |= float_flag_invalid;
+            return float32_default_nan;
+        }
+        return packFloat32( zSign, 0xFF, 0 );
+    }
+    if ( bExp == 0xFF ) {
+        if ( bSig ) return propagateFloat32NaN( a, b );
+        return packFloat32( zSign, 0, 0 );
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) {
+            if ( ( aExp | aSig ) == 0 ) {
+                roundData->exception |= float_flag_invalid;
+                return float32_default_nan;
+            }
+            roundData->exception |= float_flag_divbyzero;
+            return packFloat32( zSign, 0xFF, 0 );
+        }
+        normalizeFloat32Subnormal( bSig, &bExp, &bSig );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat32( zSign, 0, 0 );
+        normalizeFloat32Subnormal( aSig, &aExp, &aSig );
+    }
+    zExp = aExp - bExp + 0x7D;
+    aSig = ( aSig | 0x00800000 )<<7;
+    bSig = ( bSig | 0x00800000 )<<8;
+    if ( bSig <= ( aSig + aSig ) ) {
+        aSig >>= 1;
+        ++zExp;
+    }
+    {
+        bits64 tmp = ( (bits64) aSig )<<32;
+        do_div( tmp, bSig );
+        zSig = tmp;
+    }
+    if ( ( zSig & 0x3F ) == 0 ) {
+        zSig |= ( ( (bits64) bSig ) * zSig != ( (bits64) aSig )<<32 );
+    }
+    return roundAndPackFloat32( roundData, zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the remainder of the single-precision floating-point value `a'
+with respect to the corresponding value `b'.  The operation is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 float32_rem( struct roundingData *roundData, float32 a, float32 b )
+{
+    flag aSign, bSign, zSign;
+    int16 aExp, bExp, expDiff;
+    bits32 aSig, bSig;
+    bits32 q;
+    bits64 aSig64, bSig64, q64;
+    bits32 alternateASig;
+    sbits32 sigMean;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    bSig = extractFloat32Frac( b );
+    bExp = extractFloat32Exp( b );
+    bSign = extractFloat32Sign( b );
+    if ( aExp == 0xFF ) {
+        if ( aSig || ( ( bExp == 0xFF ) && bSig ) ) {
+            return propagateFloat32NaN( a, b );
+        }
+        roundData->exception |= float_flag_invalid;
+        return float32_default_nan;
+    }
+    if ( bExp == 0xFF ) {
+        if ( bSig ) return propagateFloat32NaN( a, b );
+        return a;
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) {
+            roundData->exception |= float_flag_invalid;
+            return float32_default_nan;
+        }
+        normalizeFloat32Subnormal( bSig, &bExp, &bSig );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return a;
+        normalizeFloat32Subnormal( aSig, &aExp, &aSig );
+    }
+    expDiff = aExp - bExp;
+    aSig |= 0x00800000;
+    bSig |= 0x00800000;
+    if ( expDiff < 32 ) {
+        aSig <<= 8;
+        bSig <<= 8;
+        if ( expDiff < 0 ) {
+            if ( expDiff < -1 ) return a;
+            aSig >>= 1;
+        }
+        q = ( bSig <= aSig );
+        if ( q ) aSig -= bSig;
+        if ( 0 < expDiff ) {
+            bits64 tmp = ( (bits64) aSig )<<32;
+            do_div( tmp, bSig );
+            q = tmp;
+            q >>= 32 - expDiff;
+            bSig >>= 2;
+            aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
+        }
+        else {
+            aSig >>= 2;
+            bSig >>= 2;
+        }
+    }
+    else {
+        if ( bSig <= aSig ) aSig -= bSig;
+        aSig64 = ( (bits64) aSig )<<40;
+        bSig64 = ( (bits64) bSig )<<40;
+        expDiff -= 64;
+        while ( 0 < expDiff ) {
+            q64 = estimateDiv128To64( aSig64, 0, bSig64 );
+            q64 = ( 2 < q64 ) ? q64 - 2 : 0;
+            aSig64 = - ( ( bSig * q64 )<<38 );
+            expDiff -= 62;
+        }
+        expDiff += 64;
+        q64 = estimateDiv128To64( aSig64, 0, bSig64 );
+        q64 = ( 2 < q64 ) ? q64 - 2 : 0;
+        q = q64>>( 64 - expDiff );
+        bSig <<= 6;
+        aSig = ( ( aSig64>>33 )<<( expDiff - 1 ) ) - bSig * q;
+    }
+    do {
+        alternateASig = aSig;
+        ++q;
+        aSig -= bSig;
+    } while ( 0 <= (sbits32) aSig );
+    sigMean = aSig + alternateASig;
+    if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
+        aSig = alternateASig;
+    }
+    zSign = ( (sbits32) aSig < 0 );
+    if ( zSign ) aSig = - aSig;
+    return normalizeRoundAndPackFloat32( roundData, aSign ^ zSign, bExp, aSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the square root of the single-precision floating-point value `a'.
+The operation is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 float32_sqrt( struct roundingData *roundData, float32 a )
+{
+    flag aSign;
+    int16 aExp, zExp;
+    bits32 aSig, zSig;
+    bits64 rem, term;
+
+    aSig = extractFloat32Frac( a );
+    aExp = extractFloat32Exp( a );
+    aSign = extractFloat32Sign( a );
+    if ( aExp == 0xFF ) {
+        if ( aSig ) return propagateFloat32NaN( a, 0 );
+        if ( ! aSign ) return a;
+        roundData->exception |= float_flag_invalid;
+        return float32_default_nan;
+    }
+    if ( aSign ) {
+        if ( ( aExp | aSig ) == 0 ) return a;
+        roundData->exception |= float_flag_invalid;
+        return float32_default_nan;
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return 0;
+        normalizeFloat32Subnormal( aSig, &aExp, &aSig );
+    }
+    zExp = ( ( aExp - 0x7F )>>1 ) + 0x7E;
+    aSig = ( aSig | 0x00800000 )<<8;
+    zSig = estimateSqrt32( aExp, aSig ) + 2;
+    if ( ( zSig & 0x7F ) <= 5 ) {
+        if ( zSig < 2 ) {
+            zSig = 0xFFFFFFFF;
+        }
+        else {
+            aSig >>= aExp & 1;
+            term = ( (bits64) zSig ) * zSig;
+            rem = ( ( (bits64) aSig )<<32 ) - term;
+            while ( (sbits64) rem < 0 ) {
+                --zSig;
+                rem += ( ( (bits64) zSig )<<1 ) | 1;
+            }
+            zSig |= ( rem != 0 );
+        }
+    }
+    shift32RightJamming( zSig, 1, &zSig );
+    return roundAndPackFloat32( roundData, 0, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the single-precision floating-point value `a' is equal to the
+corresponding value `b', and 0 otherwise.  The comparison is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float32_eq( float32 a, float32 b )
+{
+
+    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
+         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
+       ) {
+        if ( float32_is_signaling_nan( a ) || float32_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid );
+        }
+        return 0;
+    }
+    return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the single-precision floating-point value `a' is less than or
+equal to the corresponding value `b', and 0 otherwise.  The comparison is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float32_le( float32 a, float32 b )
+{
+    flag aSign, bSign;
+
+    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
+         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    aSign = extractFloat32Sign( a );
+    bSign = extractFloat32Sign( b );
+    if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
+    return ( a == b ) || ( aSign ^ ( a < b ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the single-precision floating-point value `a' is less than
+the corresponding value `b', and 0 otherwise.  The comparison is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float32_lt( float32 a, float32 b )
+{
+    flag aSign, bSign;
+
+    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
+         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    aSign = extractFloat32Sign( a );
+    bSign = extractFloat32Sign( b );
+    if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
+    return ( a != b ) && ( aSign ^ ( a < b ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the single-precision floating-point value `a' is equal to the
+corresponding value `b', and 0 otherwise.  The invalid exception is raised
+if either operand is a NaN.  Otherwise, the comparison is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float32_eq_signaling( float32 a, float32 b )
+{
+
+    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
+         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    return ( a == b ) || ( (bits32) ( ( a | b )<<1 ) == 0 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the single-precision floating-point value `a' is less than or
+equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not
+cause an exception.  Otherwise, the comparison is performed according to the
+IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float32_le_quiet( float32 a, float32 b )
+{
+    flag aSign, bSign;
+    //int16 aExp, bExp;
+
+    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
+         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
+       ) {
+        /* Do nothing, even if NaN as we're quiet */
+        return 0;
+    }
+    aSign = extractFloat32Sign( a );
+    bSign = extractFloat32Sign( b );
+    if ( aSign != bSign ) return aSign || ( (bits32) ( ( a | b )<<1 ) == 0 );
+    return ( a == b ) || ( aSign ^ ( a < b ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the single-precision floating-point value `a' is less than
+the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
+exception.  Otherwise, the comparison is performed according to the IEC/IEEE
+Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float32_lt_quiet( float32 a, float32 b )
+{
+    flag aSign, bSign;
+
+    if (    ( ( extractFloat32Exp( a ) == 0xFF ) && extractFloat32Frac( a ) )
+         || ( ( extractFloat32Exp( b ) == 0xFF ) && extractFloat32Frac( b ) )
+       ) {
+        /* Do nothing, even if NaN as we're quiet */
+        return 0;
+    }
+    aSign = extractFloat32Sign( a );
+    bSign = extractFloat32Sign( b );
+    if ( aSign != bSign ) return aSign && ( (bits32) ( ( a | b )<<1 ) != 0 );
+    return ( a != b ) && ( aSign ^ ( a < b ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the double-precision floating-point value
+`a' to the 32-bit two's complement integer format.  The conversion is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic---which means in particular that the conversion is rounded
+according to the current rounding mode.  If `a' is a NaN, the largest
+positive integer is returned.  Otherwise, if the conversion overflows, the
+largest integer with the same sign as `a' is returned.
+-------------------------------------------------------------------------------
+*/
+int32 float64_to_int32( struct roundingData *roundData, float64 a )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    bits64 aSig;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
+    if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
+    shiftCount = 0x42C - aExp;
+    if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig );
+    return roundAndPackInt32( roundData, aSign, aSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the double-precision floating-point value
+`a' to the 32-bit two's complement integer format.  The conversion is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic, except that the conversion is always rounded toward zero.  If
+`a' is a NaN, the largest positive integer is returned.  Otherwise, if the
+conversion overflows, the largest integer with the same sign as `a' is
+returned.
+-------------------------------------------------------------------------------
+*/
+int32 float64_to_int32_round_to_zero( float64 a )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    bits64 aSig, savedASig;
+    int32 z;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    shiftCount = 0x433 - aExp;
+    if ( shiftCount < 21 ) {
+        if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
+        goto invalid;
+    }
+    else if ( 52 < shiftCount ) {
+        if ( aExp || aSig ) float_raise( float_flag_inexact );
+        return 0;
+    }
+    aSig |= LIT64( 0x0010000000000000 );
+    savedASig = aSig;
+    aSig >>= shiftCount;
+    z = aSig;
+    if ( aSign ) z = - z;
+    if ( ( z < 0 ) ^ aSign ) {
+ invalid:
+        float_raise( float_flag_invalid );
+        return aSign ? 0x80000000 : 0x7FFFFFFF;
+    }
+    if ( ( aSig<<shiftCount ) != savedASig ) {
+        float_raise( float_flag_inexact );
+    }
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the double-precision floating-point value
+`a' to the 32-bit two's complement unsigned integer format.  The conversion
+is performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic---which means in particular that the conversion is rounded
+according to the current rounding mode.  If `a' is a NaN, the largest
+positive integer is returned.  Otherwise, if the conversion overflows, the
+largest positive integer is returned.
+-------------------------------------------------------------------------------
+*/
+int32 float64_to_uint32( struct roundingData *roundData, float64 a )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    bits64 aSig;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = 0; //extractFloat64Sign( a );
+    //if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
+    if ( aExp ) aSig |= LIT64( 0x0010000000000000 );
+    shiftCount = 0x42C - aExp;
+    if ( 0 < shiftCount ) shift64RightJamming( aSig, shiftCount, &aSig );
+    return roundAndPackInt32( roundData, aSign, aSig );
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the double-precision floating-point value
+`a' to the 32-bit two's complement integer format.  The conversion is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic, except that the conversion is always rounded toward zero.  If
+`a' is a NaN, the largest positive integer is returned.  Otherwise, if the
+conversion overflows, the largest positive integer is returned.
+-------------------------------------------------------------------------------
+*/
+int32 float64_to_uint32_round_to_zero( float64 a )
+{
+    flag aSign;
+    int16 aExp, shiftCount;
+    bits64 aSig, savedASig;
+    int32 z;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    shiftCount = 0x433 - aExp;
+    if ( shiftCount < 21 ) {
+        if ( ( aExp == 0x7FF ) && aSig ) aSign = 0;
+        goto invalid;
+    }
+    else if ( 52 < shiftCount ) {
+        if ( aExp || aSig ) float_raise( float_flag_inexact );
+        return 0;
+    }
+    aSig |= LIT64( 0x0010000000000000 );
+    savedASig = aSig;
+    aSig >>= shiftCount;
+    z = aSig;
+    if ( aSign ) z = - z;
+    if ( ( z < 0 ) ^ aSign ) {
+ invalid:
+        float_raise( float_flag_invalid );
+        return aSign ? 0x80000000 : 0x7FFFFFFF;
+    }
+    if ( ( aSig<<shiftCount ) != savedASig ) {
+        float_raise( float_flag_inexact );
+    }
+    return z;
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the double-precision floating-point value
+`a' to the single-precision floating-point format.  The conversion is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 float64_to_float32( struct roundingData *roundData, float64 a )
+{
+    flag aSign;
+    int16 aExp;
+    bits64 aSig;
+    bits32 zSig;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    if ( aExp == 0x7FF ) {
+        if ( aSig ) return commonNaNToFloat32( float64ToCommonNaN( a ) );
+        return packFloat32( aSign, 0xFF, 0 );
+    }
+    shift64RightJamming( aSig, 22, &aSig );
+    zSig = aSig;
+    if ( aExp || zSig ) {
+        zSig |= 0x40000000;
+        aExp -= 0x381;
+    }
+    return roundAndPackFloat32( roundData, aSign, aExp, zSig );
+
+}
+
+#ifdef FLOATX80
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the double-precision floating-point value
+`a' to the extended double-precision floating-point format.  The conversion
+is performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 float64_to_floatx80( float64 a )
+{
+    flag aSign;
+    int16 aExp;
+    bits64 aSig;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    if ( aExp == 0x7FF ) {
+        if ( aSig ) return commonNaNToFloatx80( float64ToCommonNaN( a ) );
+        return packFloatx80( aSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloatx80( aSign, 0, 0 );
+        normalizeFloat64Subnormal( aSig, &aExp, &aSig );
+    }
+    return
+        packFloatx80(
+            aSign, aExp + 0x3C00, ( aSig | LIT64( 0x0010000000000000 ) )<<11 );
+
+}
+
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Rounds the double-precision floating-point value `a' to an integer, and
+returns the result as a double-precision floating-point value.  The
+operation is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 float64_round_to_int( struct roundingData *roundData, float64 a )
+{
+    flag aSign;
+    int16 aExp;
+    bits64 lastBitMask, roundBitsMask;
+    int8 roundingMode;
+    float64 z;
+
+    aExp = extractFloat64Exp( a );
+    if ( 0x433 <= aExp ) {
+        if ( ( aExp == 0x7FF ) && extractFloat64Frac( a ) ) {
+            return propagateFloat64NaN( a, a );
+        }
+        return a;
+    }
+    if ( aExp <= 0x3FE ) {
+        if ( (bits64) ( a<<1 ) == 0 ) return a;
+        roundData->exception |= float_flag_inexact;
+        aSign = extractFloat64Sign( a );
+        switch ( roundData->mode ) {
+         case float_round_nearest_even:
+            if ( ( aExp == 0x3FE ) && extractFloat64Frac( a ) ) {
+                return packFloat64( aSign, 0x3FF, 0 );
+            }
+            break;
+         case float_round_down:
+            return aSign ? LIT64( 0xBFF0000000000000 ) : 0;
+         case float_round_up:
+            return
+            aSign ? LIT64( 0x8000000000000000 ) : LIT64( 0x3FF0000000000000 );
+        }
+        return packFloat64( aSign, 0, 0 );
+    }
+    lastBitMask = 1;
+    lastBitMask <<= 0x433 - aExp;
+    roundBitsMask = lastBitMask - 1;
+    z = a;
+    roundingMode = roundData->mode;
+    if ( roundingMode == float_round_nearest_even ) {
+        z += lastBitMask>>1;
+        if ( ( z & roundBitsMask ) == 0 ) z &= ~ lastBitMask;
+    }
+    else if ( roundingMode != float_round_to_zero ) {
+        if ( extractFloat64Sign( z ) ^ ( roundingMode == float_round_up ) ) {
+            z += roundBitsMask;
+        }
+    }
+    z &= ~ roundBitsMask;
+    if ( z != a ) roundData->exception |= float_flag_inexact;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of adding the absolute values of the double-precision
+floating-point values `a' and `b'.  If `zSign' is true, the sum is negated
+before being returned.  `zSign' is ignored if the result is a NaN.  The
+addition is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static float64 addFloat64Sigs( struct roundingData *roundData, float64 a, float64 b, flag zSign )
+{
+    int16 aExp, bExp, zExp;
+    bits64 aSig, bSig, zSig;
+    int16 expDiff;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    bSig = extractFloat64Frac( b );
+    bExp = extractFloat64Exp( b );
+    expDiff = aExp - bExp;
+    aSig <<= 9;
+    bSig <<= 9;
+    if ( 0 < expDiff ) {
+        if ( aExp == 0x7FF ) {
+            if ( aSig ) return propagateFloat64NaN( a, b );
+            return a;
+        }
+        if ( bExp == 0 ) {
+            --expDiff;
+        }
+        else {
+            bSig |= LIT64( 0x2000000000000000 );
+        }
+        shift64RightJamming( bSig, expDiff, &bSig );
+        zExp = aExp;
+    }
+    else if ( expDiff < 0 ) {
+        if ( bExp == 0x7FF ) {
+            if ( bSig ) return propagateFloat64NaN( a, b );
+            return packFloat64( zSign, 0x7FF, 0 );
+        }
+        if ( aExp == 0 ) {
+            ++expDiff;
+        }
+        else {
+            aSig |= LIT64( 0x2000000000000000 );
+        }
+        shift64RightJamming( aSig, - expDiff, &aSig );
+        zExp = bExp;
+    }
+    else {
+        if ( aExp == 0x7FF ) {
+            if ( aSig | bSig ) return propagateFloat64NaN( a, b );
+            return a;
+        }
+        if ( aExp == 0 ) return packFloat64( zSign, 0, ( aSig + bSig )>>9 );
+        zSig = LIT64( 0x4000000000000000 ) + aSig + bSig;
+        zExp = aExp;
+        goto roundAndPack;
+    }
+    aSig |= LIT64( 0x2000000000000000 );
+    zSig = ( aSig + bSig )<<1;
+    --zExp;
+    if ( (sbits64) zSig < 0 ) {
+        zSig = aSig + bSig;
+        ++zExp;
+    }
+ roundAndPack:
+    return roundAndPackFloat64( roundData, zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of subtracting the absolute values of the double-
+precision floating-point values `a' and `b'.  If `zSign' is true, the
+difference is negated before being returned.  `zSign' is ignored if the
+result is a NaN.  The subtraction is performed according to the IEC/IEEE
+Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static float64 subFloat64Sigs( struct roundingData *roundData, float64 a, float64 b, flag zSign )
+{
+    int16 aExp, bExp, zExp;
+    bits64 aSig, bSig, zSig;
+    int16 expDiff;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    bSig = extractFloat64Frac( b );
+    bExp = extractFloat64Exp( b );
+    expDiff = aExp - bExp;
+    aSig <<= 10;
+    bSig <<= 10;
+    if ( 0 < expDiff ) goto aExpBigger;
+    if ( expDiff < 0 ) goto bExpBigger;
+    if ( aExp == 0x7FF ) {
+        if ( aSig | bSig ) return propagateFloat64NaN( a, b );
+        roundData->exception |= float_flag_invalid;
+        return float64_default_nan;
+    }
+    if ( aExp == 0 ) {
+        aExp = 1;
+        bExp = 1;
+    }
+    if ( bSig < aSig ) goto aBigger;
+    if ( aSig < bSig ) goto bBigger;
+    return packFloat64( roundData->mode == float_round_down, 0, 0 );
+ bExpBigger:
+    if ( bExp == 0x7FF ) {
+        if ( bSig ) return propagateFloat64NaN( a, b );
+        return packFloat64( zSign ^ 1, 0x7FF, 0 );
+    }
+    if ( aExp == 0 ) {
+        ++expDiff;
+    }
+    else {
+        aSig |= LIT64( 0x4000000000000000 );
+    }
+    shift64RightJamming( aSig, - expDiff, &aSig );
+    bSig |= LIT64( 0x4000000000000000 );
+ bBigger:
+    zSig = bSig - aSig;
+    zExp = bExp;
+    zSign ^= 1;
+    goto normalizeRoundAndPack;
+ aExpBigger:
+    if ( aExp == 0x7FF ) {
+        if ( aSig ) return propagateFloat64NaN( a, b );
+        return a;
+    }
+    if ( bExp == 0 ) {
+        --expDiff;
+    }
+    else {
+        bSig |= LIT64( 0x4000000000000000 );
+    }
+    shift64RightJamming( bSig, expDiff, &bSig );
+    aSig |= LIT64( 0x4000000000000000 );
+ aBigger:
+    zSig = aSig - bSig;
+    zExp = aExp;
+ normalizeRoundAndPack:
+    --zExp;
+    return normalizeRoundAndPackFloat64( roundData, zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of adding the double-precision floating-point values `a'
+and `b'.  The operation is performed according to the IEC/IEEE Standard for
+Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 float64_add( struct roundingData *roundData, float64 a, float64 b )
+{
+    flag aSign, bSign;
+
+    aSign = extractFloat64Sign( a );
+    bSign = extractFloat64Sign( b );
+    if ( aSign == bSign ) {
+        return addFloat64Sigs( roundData, a, b, aSign );
+    }
+    else {
+        return subFloat64Sigs( roundData, a, b, aSign );
+    }
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of subtracting the double-precision floating-point values
+`a' and `b'.  The operation is performed according to the IEC/IEEE Standard
+for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 float64_sub( struct roundingData *roundData, float64 a, float64 b )
+{
+    flag aSign, bSign;
+
+    aSign = extractFloat64Sign( a );
+    bSign = extractFloat64Sign( b );
+    if ( aSign == bSign ) {
+        return subFloat64Sigs( roundData, a, b, aSign );
+    }
+    else {
+        return addFloat64Sigs( roundData, a, b, aSign );
+    }
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of multiplying the double-precision floating-point values
+`a' and `b'.  The operation is performed according to the IEC/IEEE Standard
+for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 float64_mul( struct roundingData *roundData, float64 a, float64 b )
+{
+    flag aSign, bSign, zSign;
+    int16 aExp, bExp, zExp;
+    bits64 aSig, bSig, zSig0, zSig1;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    bSig = extractFloat64Frac( b );
+    bExp = extractFloat64Exp( b );
+    bSign = extractFloat64Sign( b );
+    zSign = aSign ^ bSign;
+    if ( aExp == 0x7FF ) {
+        if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
+            return propagateFloat64NaN( a, b );
+        }
+        if ( ( bExp | bSig ) == 0 ) {
+            roundData->exception |= float_flag_invalid;
+            return float64_default_nan;
+        }
+        return packFloat64( zSign, 0x7FF, 0 );
+    }
+    if ( bExp == 0x7FF ) {
+        if ( bSig ) return propagateFloat64NaN( a, b );
+        if ( ( aExp | aSig ) == 0 ) {
+            roundData->exception |= float_flag_invalid;
+            return float64_default_nan;
+        }
+        return packFloat64( zSign, 0x7FF, 0 );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
+        normalizeFloat64Subnormal( aSig, &aExp, &aSig );
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) return packFloat64( zSign, 0, 0 );
+        normalizeFloat64Subnormal( bSig, &bExp, &bSig );
+    }
+    zExp = aExp + bExp - 0x3FF;
+    aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
+    bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
+    mul64To128( aSig, bSig, &zSig0, &zSig1 );
+    zSig0 |= ( zSig1 != 0 );
+    if ( 0 <= (sbits64) ( zSig0<<1 ) ) {
+        zSig0 <<= 1;
+        --zExp;
+    }
+    return roundAndPackFloat64( roundData, zSign, zExp, zSig0 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of dividing the double-precision floating-point value `a'
+by the corresponding value `b'.  The operation is performed according to
+the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 float64_div( struct roundingData *roundData, float64 a, float64 b )
+{
+    flag aSign, bSign, zSign;
+    int16 aExp, bExp, zExp;
+    bits64 aSig, bSig, zSig;
+    bits64 rem0, rem1;
+    bits64 term0, term1;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    bSig = extractFloat64Frac( b );
+    bExp = extractFloat64Exp( b );
+    bSign = extractFloat64Sign( b );
+    zSign = aSign ^ bSign;
+    if ( aExp == 0x7FF ) {
+        if ( aSig ) return propagateFloat64NaN( a, b );
+        if ( bExp == 0x7FF ) {
+            if ( bSig ) return propagateFloat64NaN( a, b );
+            roundData->exception |= float_flag_invalid;
+            return float64_default_nan;
+        }
+        return packFloat64( zSign, 0x7FF, 0 );
+    }
+    if ( bExp == 0x7FF ) {
+        if ( bSig ) return propagateFloat64NaN( a, b );
+        return packFloat64( zSign, 0, 0 );
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) {
+            if ( ( aExp | aSig ) == 0 ) {
+                roundData->exception |= float_flag_invalid;
+                return float64_default_nan;
+            }
+            roundData->exception |= float_flag_divbyzero;
+            return packFloat64( zSign, 0x7FF, 0 );
+        }
+        normalizeFloat64Subnormal( bSig, &bExp, &bSig );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloat64( zSign, 0, 0 );
+        normalizeFloat64Subnormal( aSig, &aExp, &aSig );
+    }
+    zExp = aExp - bExp + 0x3FD;
+    aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<10;
+    bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
+    if ( bSig <= ( aSig + aSig ) ) {
+        aSig >>= 1;
+        ++zExp;
+    }
+    zSig = estimateDiv128To64( aSig, 0, bSig );
+    if ( ( zSig & 0x1FF ) <= 2 ) {
+        mul64To128( bSig, zSig, &term0, &term1 );
+        sub128( aSig, 0, term0, term1, &rem0, &rem1 );
+        while ( (sbits64) rem0 < 0 ) {
+            --zSig;
+            add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
+        }
+        zSig |= ( rem1 != 0 );
+    }
+    return roundAndPackFloat64( roundData, zSign, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the remainder of the double-precision floating-point value `a'
+with respect to the corresponding value `b'.  The operation is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 float64_rem( struct roundingData *roundData, float64 a, float64 b )
+{
+    flag aSign, bSign, zSign;
+    int16 aExp, bExp, expDiff;
+    bits64 aSig, bSig;
+    bits64 q, alternateASig;
+    sbits64 sigMean;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    bSig = extractFloat64Frac( b );
+    bExp = extractFloat64Exp( b );
+    bSign = extractFloat64Sign( b );
+    if ( aExp == 0x7FF ) {
+        if ( aSig || ( ( bExp == 0x7FF ) && bSig ) ) {
+            return propagateFloat64NaN( a, b );
+        }
+        roundData->exception |= float_flag_invalid;
+        return float64_default_nan;
+    }
+    if ( bExp == 0x7FF ) {
+        if ( bSig ) return propagateFloat64NaN( a, b );
+        return a;
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) {
+            roundData->exception |= float_flag_invalid;
+            return float64_default_nan;
+        }
+        normalizeFloat64Subnormal( bSig, &bExp, &bSig );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return a;
+        normalizeFloat64Subnormal( aSig, &aExp, &aSig );
+    }
+    expDiff = aExp - bExp;
+    aSig = ( aSig | LIT64( 0x0010000000000000 ) )<<11;
+    bSig = ( bSig | LIT64( 0x0010000000000000 ) )<<11;
+    if ( expDiff < 0 ) {
+        if ( expDiff < -1 ) return a;
+        aSig >>= 1;
+    }
+    q = ( bSig <= aSig );
+    if ( q ) aSig -= bSig;
+    expDiff -= 64;
+    while ( 0 < expDiff ) {
+        q = estimateDiv128To64( aSig, 0, bSig );
+        q = ( 2 < q ) ? q - 2 : 0;
+        aSig = - ( ( bSig>>2 ) * q );
+        expDiff -= 62;
+    }
+    expDiff += 64;
+    if ( 0 < expDiff ) {
+        q = estimateDiv128To64( aSig, 0, bSig );
+        q = ( 2 < q ) ? q - 2 : 0;
+        q >>= 64 - expDiff;
+        bSig >>= 2;
+        aSig = ( ( aSig>>1 )<<( expDiff - 1 ) ) - bSig * q;
+    }
+    else {
+        aSig >>= 2;
+        bSig >>= 2;
+    }
+    do {
+        alternateASig = aSig;
+        ++q;
+        aSig -= bSig;
+    } while ( 0 <= (sbits64) aSig );
+    sigMean = aSig + alternateASig;
+    if ( ( sigMean < 0 ) || ( ( sigMean == 0 ) && ( q & 1 ) ) ) {
+        aSig = alternateASig;
+    }
+    zSign = ( (sbits64) aSig < 0 );
+    if ( zSign ) aSig = - aSig;
+    return normalizeRoundAndPackFloat64( roundData, aSign ^ zSign, bExp, aSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the square root of the double-precision floating-point value `a'.
+The operation is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 float64_sqrt( struct roundingData *roundData, float64 a )
+{
+    flag aSign;
+    int16 aExp, zExp;
+    bits64 aSig, zSig;
+    bits64 rem0, rem1, term0, term1; //, shiftedRem;
+    //float64 z;
+
+    aSig = extractFloat64Frac( a );
+    aExp = extractFloat64Exp( a );
+    aSign = extractFloat64Sign( a );
+    if ( aExp == 0x7FF ) {
+        if ( aSig ) return propagateFloat64NaN( a, a );
+        if ( ! aSign ) return a;
+        roundData->exception |= float_flag_invalid;
+        return float64_default_nan;
+    }
+    if ( aSign ) {
+        if ( ( aExp | aSig ) == 0 ) return a;
+        roundData->exception |= float_flag_invalid;
+        return float64_default_nan;
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return 0;
+        normalizeFloat64Subnormal( aSig, &aExp, &aSig );
+    }
+    zExp = ( ( aExp - 0x3FF )>>1 ) + 0x3FE;
+    aSig |= LIT64( 0x0010000000000000 );
+    zSig = estimateSqrt32( aExp, aSig>>21 );
+    zSig <<= 31;
+    aSig <<= 9 - ( aExp & 1 );
+    zSig = estimateDiv128To64( aSig, 0, zSig ) + zSig + 2;
+    if ( ( zSig & 0x3FF ) <= 5 ) {
+        if ( zSig < 2 ) {
+            zSig = LIT64( 0xFFFFFFFFFFFFFFFF );
+        }
+        else {
+            aSig <<= 2;
+            mul64To128( zSig, zSig, &term0, &term1 );
+            sub128( aSig, 0, term0, term1, &rem0, &rem1 );
+            while ( (sbits64) rem0 < 0 ) {
+                --zSig;
+                shortShift128Left( 0, zSig, 1, &term0, &term1 );
+                term1 |= 1;
+                add128( rem0, rem1, term0, term1, &rem0, &rem1 );
+            }
+            zSig |= ( ( rem0 | rem1 ) != 0 );
+        }
+    }
+    shift64RightJamming( zSig, 1, &zSig );
+    return roundAndPackFloat64( roundData, 0, zExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the double-precision floating-point value `a' is equal to the
+corresponding value `b', and 0 otherwise.  The comparison is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float64_eq( float64 a, float64 b )
+{
+
+    if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
+         || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
+       ) {
+        if ( float64_is_signaling_nan( a ) || float64_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid );
+        }
+        return 0;
+    }
+    return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the double-precision floating-point value `a' is less than or
+equal to the corresponding value `b', and 0 otherwise.  The comparison is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float64_le( float64 a, float64 b )
+{
+    flag aSign, bSign;
+
+    if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
+         || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    aSign = extractFloat64Sign( a );
+    bSign = extractFloat64Sign( b );
+    if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
+    return ( a == b ) || ( aSign ^ ( a < b ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the double-precision floating-point value `a' is less than
+the corresponding value `b', and 0 otherwise.  The comparison is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float64_lt( float64 a, float64 b )
+{
+    flag aSign, bSign;
+
+    if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
+         || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    aSign = extractFloat64Sign( a );
+    bSign = extractFloat64Sign( b );
+    if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
+    return ( a != b ) && ( aSign ^ ( a < b ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the double-precision floating-point value `a' is equal to the
+corresponding value `b', and 0 otherwise.  The invalid exception is raised
+if either operand is a NaN.  Otherwise, the comparison is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float64_eq_signaling( float64 a, float64 b )
+{
+
+    if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
+         || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    return ( a == b ) || ( (bits64) ( ( a | b )<<1 ) == 0 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the double-precision floating-point value `a' is less than or
+equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs do not
+cause an exception.  Otherwise, the comparison is performed according to the
+IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float64_le_quiet( float64 a, float64 b )
+{
+    flag aSign, bSign;
+    //int16 aExp, bExp;
+
+    if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
+         || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
+       ) {
+        /* Do nothing, even if NaN as we're quiet */
+        return 0;
+    }
+    aSign = extractFloat64Sign( a );
+    bSign = extractFloat64Sign( b );
+    if ( aSign != bSign ) return aSign || ( (bits64) ( ( a | b )<<1 ) == 0 );
+    return ( a == b ) || ( aSign ^ ( a < b ) );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the double-precision floating-point value `a' is less than
+the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause an
+exception.  Otherwise, the comparison is performed according to the IEC/IEEE
+Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag float64_lt_quiet( float64 a, float64 b )
+{
+    flag aSign, bSign;
+
+    if (    ( ( extractFloat64Exp( a ) == 0x7FF ) && extractFloat64Frac( a ) )
+         || ( ( extractFloat64Exp( b ) == 0x7FF ) && extractFloat64Frac( b ) )
+       ) {
+        /* Do nothing, even if NaN as we're quiet */
+        return 0;
+    }
+    aSign = extractFloat64Sign( a );
+    bSign = extractFloat64Sign( b );
+    if ( aSign != bSign ) return aSign && ( (bits64) ( ( a | b )<<1 ) != 0 );
+    return ( a != b ) && ( aSign ^ ( a < b ) );
+
+}
+
+#ifdef FLOATX80
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the extended double-precision floating-
+point value `a' to the 32-bit two's complement integer format.  The
+conversion is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic---which means in particular that the conversion
+is rounded according to the current rounding mode.  If `a' is a NaN, the
+largest positive integer is returned.  Otherwise, if the conversion
+overflows, the largest integer with the same sign as `a' is returned.
+-------------------------------------------------------------------------------
+*/
+int32 floatx80_to_int32( struct roundingData *roundData, floatx80 a )
+{
+    flag aSign;
+    int32 aExp, shiftCount;
+    bits64 aSig;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+    if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
+    shiftCount = 0x4037 - aExp;
+    if ( shiftCount <= 0 ) shiftCount = 1;
+    shift64RightJamming( aSig, shiftCount, &aSig );
+    return roundAndPackInt32( roundData, aSign, aSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the extended double-precision floating-
+point value `a' to the 32-bit two's complement integer format.  The
+conversion is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic, except that the conversion is always rounded
+toward zero.  If `a' is a NaN, the largest positive integer is returned.
+Otherwise, if the conversion overflows, the largest integer with the same
+sign as `a' is returned.
+-------------------------------------------------------------------------------
+*/
+int32 floatx80_to_int32_round_to_zero( floatx80 a )
+{
+    flag aSign;
+    int32 aExp, shiftCount;
+    bits64 aSig, savedASig;
+    int32 z;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+    shiftCount = 0x403E - aExp;
+    if ( shiftCount < 32 ) {
+        if ( ( aExp == 0x7FFF ) && (bits64) ( aSig<<1 ) ) aSign = 0;
+        goto invalid;
+    }
+    else if ( 63 < shiftCount ) {
+        if ( aExp || aSig ) float_raise( float_flag_inexact );
+        return 0;
+    }
+    savedASig = aSig;
+    aSig >>= shiftCount;
+    z = aSig;
+    if ( aSign ) z = - z;
+    if ( ( z < 0 ) ^ aSign ) {
+ invalid:
+        float_raise( float_flag_invalid );
+        return aSign ? 0x80000000 : 0x7FFFFFFF;
+    }
+    if ( ( aSig<<shiftCount ) != savedASig ) {
+        float_raise( float_flag_inexact );
+    }
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the extended double-precision floating-
+point value `a' to the single-precision floating-point format.  The
+conversion is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float32 floatx80_to_float32( struct roundingData *roundData, floatx80 a )
+{
+    flag aSign;
+    int32 aExp;
+    bits64 aSig;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+    if ( aExp == 0x7FFF ) {
+        if ( (bits64) ( aSig<<1 ) ) {
+            return commonNaNToFloat32( floatx80ToCommonNaN( a ) );
+        }
+        return packFloat32( aSign, 0xFF, 0 );
+    }
+    shift64RightJamming( aSig, 33, &aSig );
+    if ( aExp || aSig ) aExp -= 0x3F81;
+    return roundAndPackFloat32( roundData, aSign, aExp, aSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of converting the extended double-precision floating-
+point value `a' to the double-precision floating-point format.  The
+conversion is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+float64 floatx80_to_float64( struct roundingData *roundData, floatx80 a )
+{
+    flag aSign;
+    int32 aExp;
+    bits64 aSig, zSig;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+    if ( aExp == 0x7FFF ) {
+        if ( (bits64) ( aSig<<1 ) ) {
+            return commonNaNToFloat64( floatx80ToCommonNaN( a ) );
+        }
+        return packFloat64( aSign, 0x7FF, 0 );
+    }
+    shift64RightJamming( aSig, 1, &zSig );
+    if ( aExp || aSig ) aExp -= 0x3C01;
+    return roundAndPackFloat64( roundData, aSign, aExp, zSig );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Rounds the extended double-precision floating-point value `a' to an integer,
+and returns the result as an extended quadruple-precision floating-point
+value.  The operation is performed according to the IEC/IEEE Standard for
+Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 floatx80_round_to_int( struct roundingData *roundData, floatx80 a )
+{
+    flag aSign;
+    int32 aExp;
+    bits64 lastBitMask, roundBitsMask;
+    int8 roundingMode;
+    floatx80 z;
+
+    aExp = extractFloatx80Exp( a );
+    if ( 0x403E <= aExp ) {
+        if ( ( aExp == 0x7FFF ) && (bits64) ( extractFloatx80Frac( a )<<1 ) ) {
+            return propagateFloatx80NaN( a, a );
+        }
+        return a;
+    }
+    if ( aExp <= 0x3FFE ) {
+        if (    ( aExp == 0 )
+             && ( (bits64) ( extractFloatx80Frac( a )<<1 ) == 0 ) ) {
+            return a;
+        }
+        roundData->exception |= float_flag_inexact;
+        aSign = extractFloatx80Sign( a );
+        switch ( roundData->mode ) {
+         case float_round_nearest_even:
+            if ( ( aExp == 0x3FFE ) && (bits64) ( extractFloatx80Frac( a )<<1 )
+               ) {
+                return
+                    packFloatx80( aSign, 0x3FFF, LIT64( 0x8000000000000000 ) );
+            }
+            break;
+         case float_round_down:
+            return
+                  aSign ?
+                      packFloatx80( 1, 0x3FFF, LIT64( 0x8000000000000000 ) )
+                : packFloatx80( 0, 0, 0 );
+         case float_round_up:
+            return
+                  aSign ? packFloatx80( 1, 0, 0 )
+                : packFloatx80( 0, 0x3FFF, LIT64( 0x8000000000000000 ) );
+        }
+        return packFloatx80( aSign, 0, 0 );
+    }
+    lastBitMask = 1;
+    lastBitMask <<= 0x403E - aExp;
+    roundBitsMask = lastBitMask - 1;
+    z = a;
+    roundingMode = roundData->mode;
+    if ( roundingMode == float_round_nearest_even ) {
+        z.low += lastBitMask>>1;
+        if ( ( z.low & roundBitsMask ) == 0 ) z.low &= ~ lastBitMask;
+    }
+    else if ( roundingMode != float_round_to_zero ) {
+        if ( extractFloatx80Sign( z ) ^ ( roundingMode == float_round_up ) ) {
+            z.low += roundBitsMask;
+        }
+    }
+    z.low &= ~ roundBitsMask;
+    if ( z.low == 0 ) {
+        ++z.high;
+        z.low = LIT64( 0x8000000000000000 );
+    }
+    if ( z.low != a.low ) roundData->exception |= float_flag_inexact;
+    return z;
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of adding the absolute values of the extended double-
+precision floating-point values `a' and `b'.  If `zSign' is true, the sum is
+negated before being returned.  `zSign' is ignored if the result is a NaN.
+The addition is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static floatx80 addFloatx80Sigs( struct roundingData *roundData, floatx80 a, floatx80 b, flag zSign )
+{
+    int32 aExp, bExp, zExp;
+    bits64 aSig, bSig, zSig0, zSig1;
+    int32 expDiff;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    bSig = extractFloatx80Frac( b );
+    bExp = extractFloatx80Exp( b );
+    expDiff = aExp - bExp;
+    if ( 0 < expDiff ) {
+        if ( aExp == 0x7FFF ) {
+            if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
+            return a;
+        }
+        if ( bExp == 0 ) --expDiff;
+        shift64ExtraRightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
+        zExp = aExp;
+    }
+    else if ( expDiff < 0 ) {
+        if ( bExp == 0x7FFF ) {
+            if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
+            return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
+        }
+        if ( aExp == 0 ) ++expDiff;
+        shift64ExtraRightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
+        zExp = bExp;
+    }
+    else {
+        if ( aExp == 0x7FFF ) {
+            if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
+                return propagateFloatx80NaN( a, b );
+            }
+            return a;
+        }
+        zSig1 = 0;
+        zSig0 = aSig + bSig;
+        if ( aExp == 0 ) {
+            normalizeFloatx80Subnormal( zSig0, &zExp, &zSig0 );
+            goto roundAndPack;
+        }
+        zExp = aExp;
+        goto shiftRight1;
+    }
+    
+    zSig0 = aSig + bSig;
+
+    if ( (sbits64) zSig0 < 0 ) goto roundAndPack; 
+ shiftRight1:
+    shift64ExtraRightJamming( zSig0, zSig1, 1, &zSig0, &zSig1 );
+    zSig0 |= LIT64( 0x8000000000000000 );
+    ++zExp;
+ roundAndPack:
+    return
+        roundAndPackFloatx80(
+            roundData, zSign, zExp, zSig0, zSig1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of subtracting the absolute values of the extended
+double-precision floating-point values `a' and `b'.  If `zSign' is true,
+the difference is negated before being returned.  `zSign' is ignored if the
+result is a NaN.  The subtraction is performed according to the IEC/IEEE
+Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+static floatx80 subFloatx80Sigs( struct roundingData *roundData, floatx80 a, floatx80 b, flag zSign )
+{
+    int32 aExp, bExp, zExp;
+    bits64 aSig, bSig, zSig0, zSig1;
+    int32 expDiff;
+    floatx80 z;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    bSig = extractFloatx80Frac( b );
+    bExp = extractFloatx80Exp( b );
+    expDiff = aExp - bExp;
+    if ( 0 < expDiff ) goto aExpBigger;
+    if ( expDiff < 0 ) goto bExpBigger;
+    if ( aExp == 0x7FFF ) {
+        if ( (bits64) ( ( aSig | bSig )<<1 ) ) {
+            return propagateFloatx80NaN( a, b );
+        }
+        roundData->exception |= float_flag_invalid;
+        z.low = floatx80_default_nan_low;
+        z.high = floatx80_default_nan_high;
+        z.__padding = 0;
+        return z;
+    }
+    if ( aExp == 0 ) {
+        aExp = 1;
+        bExp = 1;
+    }
+    zSig1 = 0;
+    if ( bSig < aSig ) goto aBigger;
+    if ( aSig < bSig ) goto bBigger;
+    return packFloatx80( roundData->mode == float_round_down, 0, 0 );
+ bExpBigger:
+    if ( bExp == 0x7FFF ) {
+        if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
+        return packFloatx80( zSign ^ 1, 0x7FFF, LIT64( 0x8000000000000000 ) );
+    }
+    if ( aExp == 0 ) ++expDiff;
+    shift128RightJamming( aSig, 0, - expDiff, &aSig, &zSig1 );
+ bBigger:
+    sub128( bSig, 0, aSig, zSig1, &zSig0, &zSig1 );
+    zExp = bExp;
+    zSign ^= 1;
+    goto normalizeRoundAndPack;
+ aExpBigger:
+    if ( aExp == 0x7FFF ) {
+        if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
+        return a;
+    }
+    if ( bExp == 0 ) --expDiff;
+    shift128RightJamming( bSig, 0, expDiff, &bSig, &zSig1 );
+ aBigger:
+    sub128( aSig, 0, bSig, zSig1, &zSig0, &zSig1 );
+    zExp = aExp;
+ normalizeRoundAndPack:
+    return
+        normalizeRoundAndPackFloatx80(
+            roundData, zSign, zExp, zSig0, zSig1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of adding the extended double-precision floating-point
+values `a' and `b'.  The operation is performed according to the IEC/IEEE
+Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 floatx80_add( struct roundingData *roundData, floatx80 a, floatx80 b )
+{
+    flag aSign, bSign;
+    
+    aSign = extractFloatx80Sign( a );
+    bSign = extractFloatx80Sign( b );
+    if ( aSign == bSign ) {
+        return addFloatx80Sigs( roundData, a, b, aSign );
+    }
+    else {
+        return subFloatx80Sigs( roundData, a, b, aSign );
+    }
+    
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of subtracting the extended double-precision floating-
+point values `a' and `b'.  The operation is performed according to the
+IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 floatx80_sub( struct roundingData *roundData, floatx80 a, floatx80 b )
+{
+    flag aSign, bSign;
+
+    aSign = extractFloatx80Sign( a );
+    bSign = extractFloatx80Sign( b );
+    if ( aSign == bSign ) {
+        return subFloatx80Sigs( roundData, a, b, aSign );
+    }
+    else {
+        return addFloatx80Sigs( roundData, a, b, aSign );
+    }
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of multiplying the extended double-precision floating-
+point values `a' and `b'.  The operation is performed according to the
+IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 floatx80_mul( struct roundingData *roundData, floatx80 a, floatx80 b )
+{
+    flag aSign, bSign, zSign;
+    int32 aExp, bExp, zExp;
+    bits64 aSig, bSig, zSig0, zSig1;
+    floatx80 z;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+    bSig = extractFloatx80Frac( b );
+    bExp = extractFloatx80Exp( b );
+    bSign = extractFloatx80Sign( b );
+    zSign = aSign ^ bSign;
+    if ( aExp == 0x7FFF ) {
+        if (    (bits64) ( aSig<<1 )
+             || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
+            return propagateFloatx80NaN( a, b );
+        }
+        if ( ( bExp | bSig ) == 0 ) goto invalid;
+        return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
+    }
+    if ( bExp == 0x7FFF ) {
+        if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
+        if ( ( aExp | aSig ) == 0 ) {
+ invalid:
+            roundData->exception |= float_flag_invalid;
+            z.low = floatx80_default_nan_low;
+            z.high = floatx80_default_nan_high;
+            z.__padding = 0;
+            return z;
+        }
+        return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
+        normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) return packFloatx80( zSign, 0, 0 );
+        normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
+    }
+    zExp = aExp + bExp - 0x3FFE;
+    mul64To128( aSig, bSig, &zSig0, &zSig1 );
+    if ( 0 < (sbits64) zSig0 ) {
+        shortShift128Left( zSig0, zSig1, 1, &zSig0, &zSig1 );
+        --zExp;
+    }
+    return
+        roundAndPackFloatx80(
+            roundData, zSign, zExp, zSig0, zSig1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the result of dividing the extended double-precision floating-point
+value `a' by the corresponding value `b'.  The operation is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 floatx80_div( struct roundingData *roundData, floatx80 a, floatx80 b )
+{
+    flag aSign, bSign, zSign;
+    int32 aExp, bExp, zExp;
+    bits64 aSig, bSig, zSig0, zSig1;
+    bits64 rem0, rem1, rem2, term0, term1, term2;
+    floatx80 z;
+
+    aSig = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+    bSig = extractFloatx80Frac( b );
+    bExp = extractFloatx80Exp( b );
+    bSign = extractFloatx80Sign( b );
+    zSign = aSign ^ bSign;
+    if ( aExp == 0x7FFF ) {
+        if ( (bits64) ( aSig<<1 ) ) return propagateFloatx80NaN( a, b );
+        if ( bExp == 0x7FFF ) {
+            if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
+            goto invalid;
+        }
+        return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
+    }
+    if ( bExp == 0x7FFF ) {
+        if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
+        return packFloatx80( zSign, 0, 0 );
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) {
+            if ( ( aExp | aSig ) == 0 ) {
+ invalid:
+                roundData->exception |= float_flag_invalid;
+                z.low = floatx80_default_nan_low;
+                z.high = floatx80_default_nan_high;
+                z.__padding = 0;
+                return z;
+            }
+            roundData->exception |= float_flag_divbyzero;
+            return packFloatx80( zSign, 0x7FFF, LIT64( 0x8000000000000000 ) );
+        }
+        normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
+    }
+    if ( aExp == 0 ) {
+        if ( aSig == 0 ) return packFloatx80( zSign, 0, 0 );
+        normalizeFloatx80Subnormal( aSig, &aExp, &aSig );
+    }
+    zExp = aExp - bExp + 0x3FFE;
+    rem1 = 0;
+    if ( bSig <= aSig ) {
+        shift128Right( aSig, 0, 1, &aSig, &rem1 );
+        ++zExp;
+    }
+    zSig0 = estimateDiv128To64( aSig, rem1, bSig );
+    mul64To128( bSig, zSig0, &term0, &term1 );
+    sub128( aSig, rem1, term0, term1, &rem0, &rem1 );
+    while ( (sbits64) rem0 < 0 ) {
+        --zSig0;
+        add128( rem0, rem1, 0, bSig, &rem0, &rem1 );
+    }
+    zSig1 = estimateDiv128To64( rem1, 0, bSig );
+    if ( (bits64) ( zSig1<<1 ) <= 8 ) {
+        mul64To128( bSig, zSig1, &term1, &term2 );
+        sub128( rem1, 0, term1, term2, &rem1, &rem2 );
+        while ( (sbits64) rem1 < 0 ) {
+            --zSig1;
+            add128( rem1, rem2, 0, bSig, &rem1, &rem2 );
+        }
+        zSig1 |= ( ( rem1 | rem2 ) != 0 );
+    }
+    return
+        roundAndPackFloatx80(
+            roundData, zSign, zExp, zSig0, zSig1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the remainder of the extended double-precision floating-point value
+`a' with respect to the corresponding value `b'.  The operation is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 floatx80_rem( struct roundingData *roundData, floatx80 a, floatx80 b )
+{
+    flag aSign, bSign, zSign;
+    int32 aExp, bExp, expDiff;
+    bits64 aSig0, aSig1, bSig;
+    bits64 q, term0, term1, alternateASig0, alternateASig1;
+    floatx80 z;
+
+    aSig0 = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+    bSig = extractFloatx80Frac( b );
+    bExp = extractFloatx80Exp( b );
+    bSign = extractFloatx80Sign( b );
+    if ( aExp == 0x7FFF ) {
+        if (    (bits64) ( aSig0<<1 )
+             || ( ( bExp == 0x7FFF ) && (bits64) ( bSig<<1 ) ) ) {
+            return propagateFloatx80NaN( a, b );
+        }
+        goto invalid;
+    }
+    if ( bExp == 0x7FFF ) {
+        if ( (bits64) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b );
+        return a;
+    }
+    if ( bExp == 0 ) {
+        if ( bSig == 0 ) {
+ invalid:
+            roundData->exception |= float_flag_invalid;
+            z.low = floatx80_default_nan_low;
+            z.high = floatx80_default_nan_high;
+            z.__padding = 0;
+            return z;
+        }
+        normalizeFloatx80Subnormal( bSig, &bExp, &bSig );
+    }
+    if ( aExp == 0 ) {
+        if ( (bits64) ( aSig0<<1 ) == 0 ) return a;
+        normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
+    }
+    bSig |= LIT64( 0x8000000000000000 );
+    zSign = aSign;
+    expDiff = aExp - bExp;
+    aSig1 = 0;
+    if ( expDiff < 0 ) {
+        if ( expDiff < -1 ) return a;
+        shift128Right( aSig0, 0, 1, &aSig0, &aSig1 );
+        expDiff = 0;
+    }
+    q = ( bSig <= aSig0 );
+    if ( q ) aSig0 -= bSig;
+    expDiff -= 64;
+    while ( 0 < expDiff ) {
+        q = estimateDiv128To64( aSig0, aSig1, bSig );
+        q = ( 2 < q ) ? q - 2 : 0;
+        mul64To128( bSig, q, &term0, &term1 );
+        sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
+        shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 );
+        expDiff -= 62;
+    }
+    expDiff += 64;
+    if ( 0 < expDiff ) {
+        q = estimateDiv128To64( aSig0, aSig1, bSig );
+        q = ( 2 < q ) ? q - 2 : 0;
+        q >>= 64 - expDiff;
+        mul64To128( bSig, q<<( 64 - expDiff ), &term0, &term1 );
+        sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
+        shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 );
+        while ( le128( term0, term1, aSig0, aSig1 ) ) {
+            ++q;
+            sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 );
+        }
+    }
+    else {
+        term1 = 0;
+        term0 = bSig;
+    }
+    sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 );
+    if (    lt128( alternateASig0, alternateASig1, aSig0, aSig1 )
+         || (    eq128( alternateASig0, alternateASig1, aSig0, aSig1 )
+              && ( q & 1 ) )
+       ) {
+        aSig0 = alternateASig0;
+        aSig1 = alternateASig1;
+        zSign = ! zSign;
+    }
+
+    return
+        normalizeRoundAndPackFloatx80(
+            roundData, zSign, bExp + expDiff, aSig0, aSig1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns the square root of the extended double-precision floating-point
+value `a'.  The operation is performed according to the IEC/IEEE Standard
+for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+floatx80 floatx80_sqrt( struct roundingData *roundData, floatx80 a )
+{
+    flag aSign;
+    int32 aExp, zExp;
+    bits64 aSig0, aSig1, zSig0, zSig1;
+    bits64 rem0, rem1, rem2, rem3, term0, term1, term2, term3;
+    bits64 shiftedRem0, shiftedRem1;
+    floatx80 z;
+
+    aSig0 = extractFloatx80Frac( a );
+    aExp = extractFloatx80Exp( a );
+    aSign = extractFloatx80Sign( a );
+    if ( aExp == 0x7FFF ) {
+        if ( (bits64) ( aSig0<<1 ) ) return propagateFloatx80NaN( a, a );
+        if ( ! aSign ) return a;
+        goto invalid;
+    }
+    if ( aSign ) {
+        if ( ( aExp | aSig0 ) == 0 ) return a;
+ invalid:
+        roundData->exception |= float_flag_invalid;
+        z.low = floatx80_default_nan_low;
+        z.high = floatx80_default_nan_high;
+        z.__padding = 0;
+        return z;
+    }
+    if ( aExp == 0 ) {
+        if ( aSig0 == 0 ) return packFloatx80( 0, 0, 0 );
+        normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 );
+    }
+    zExp = ( ( aExp - 0x3FFF )>>1 ) + 0x3FFF;
+    zSig0 = estimateSqrt32( aExp, aSig0>>32 );
+    zSig0 <<= 31;
+    aSig1 = 0;
+    shift128Right( aSig0, 0, ( aExp & 1 ) + 2, &aSig0, &aSig1 );
+    zSig0 = estimateDiv128To64( aSig0, aSig1, zSig0 ) + zSig0 + 4;
+    if ( 0 <= (sbits64) zSig0 ) zSig0 = LIT64( 0xFFFFFFFFFFFFFFFF );
+    shortShift128Left( aSig0, aSig1, 2, &aSig0, &aSig1 );
+    mul64To128( zSig0, zSig0, &term0, &term1 );
+    sub128( aSig0, aSig1, term0, term1, &rem0, &rem1 );
+    while ( (sbits64) rem0 < 0 ) {
+        --zSig0;
+        shortShift128Left( 0, zSig0, 1, &term0, &term1 );
+        term1 |= 1;
+        add128( rem0, rem1, term0, term1, &rem0, &rem1 );
+    }
+    shortShift128Left( rem0, rem1, 63, &shiftedRem0, &shiftedRem1 );
+    zSig1 = estimateDiv128To64( shiftedRem0, shiftedRem1, zSig0 );
+    if ( (bits64) ( zSig1<<1 ) <= 10 ) {
+        if ( zSig1 == 0 ) zSig1 = 1;
+        mul64To128( zSig0, zSig1, &term1, &term2 );
+        shortShift128Left( term1, term2, 1, &term1, &term2 );
+        sub128( rem1, 0, term1, term2, &rem1, &rem2 );
+        mul64To128( zSig1, zSig1, &term2, &term3 );
+        sub192( rem1, rem2, 0, 0, term2, term3, &rem1, &rem2, &rem3 );
+        while ( (sbits64) rem1 < 0 ) {
+            --zSig1;
+            shortShift192Left( 0, zSig0, zSig1, 1, &term1, &term2, &term3 );
+            term3 |= 1;
+            add192(
+                rem1, rem2, rem3, term1, term2, term3, &rem1, &rem2, &rem3 );
+        }
+        zSig1 |= ( ( rem1 | rem2 | rem3 ) != 0 );
+    }
+    return
+        roundAndPackFloatx80(
+            roundData, 0, zExp, zSig0, zSig1 );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the extended double-precision floating-point value `a' is
+equal to the corresponding value `b', and 0 otherwise.  The comparison is
+performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag floatx80_eq( floatx80 a, floatx80 b )
+{
+
+    if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+       ) {
+        if (    floatx80_is_signaling_nan( a )
+             || floatx80_is_signaling_nan( b ) ) {
+            float_raise( float_flag_invalid );
+        }
+        return 0;
+    }
+    return
+           ( a.low == b.low )
+        && (    ( a.high == b.high )
+             || (    ( a.low == 0 )
+                  && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
+           );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the extended double-precision floating-point value `a' is
+less than or equal to the corresponding value `b', and 0 otherwise.  The
+comparison is performed according to the IEC/IEEE Standard for Binary
+Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag floatx80_le( floatx80 a, floatx80 b )
+{
+    flag aSign, bSign;
+
+    if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    aSign = extractFloatx80Sign( a );
+    bSign = extractFloatx80Sign( b );
+    if ( aSign != bSign ) {
+        return
+               aSign
+            || (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+                 == 0 );
+    }
+    return
+          aSign ? le128( b.high, b.low, a.high, a.low )
+        : le128( a.high, a.low, b.high, b.low );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the extended double-precision floating-point value `a' is
+less than the corresponding value `b', and 0 otherwise.  The comparison
+is performed according to the IEC/IEEE Standard for Binary Floating-point
+Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag floatx80_lt( floatx80 a, floatx80 b )
+{
+    flag aSign, bSign;
+
+    if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    aSign = extractFloatx80Sign( a );
+    bSign = extractFloatx80Sign( b );
+    if ( aSign != bSign ) {
+        return
+               aSign
+            && (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+                 != 0 );
+    }
+    return
+          aSign ? lt128( b.high, b.low, a.high, a.low )
+        : lt128( a.high, a.low, b.high, b.low );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the extended double-precision floating-point value `a' is equal
+to the corresponding value `b', and 0 otherwise.  The invalid exception is
+raised if either operand is a NaN.  Otherwise, the comparison is performed
+according to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag floatx80_eq_signaling( floatx80 a, floatx80 b )
+{
+
+    if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+       ) {
+        float_raise( float_flag_invalid );
+        return 0;
+    }
+    return
+           ( a.low == b.low )
+        && (    ( a.high == b.high )
+             || (    ( a.low == 0 )
+                  && ( (bits16) ( ( a.high | b.high )<<1 ) == 0 ) )
+           );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the extended double-precision floating-point value `a' is less
+than or equal to the corresponding value `b', and 0 otherwise.  Quiet NaNs
+do not cause an exception.  Otherwise, the comparison is performed according
+to the IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag floatx80_le_quiet( floatx80 a, floatx80 b )
+{
+    flag aSign, bSign;
+
+    if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+       ) {
+        /* Do nothing, even if NaN as we're quiet */
+        return 0;
+    }
+    aSign = extractFloatx80Sign( a );
+    bSign = extractFloatx80Sign( b );
+    if ( aSign != bSign ) {
+        return
+               aSign
+            || (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+                 == 0 );
+    }
+    return
+          aSign ? le128( b.high, b.low, a.high, a.low )
+        : le128( a.high, a.low, b.high, b.low );
+
+}
+
+/*
+-------------------------------------------------------------------------------
+Returns 1 if the extended double-precision floating-point value `a' is less
+than the corresponding value `b', and 0 otherwise.  Quiet NaNs do not cause
+an exception.  Otherwise, the comparison is performed according to the
+IEC/IEEE Standard for Binary Floating-point Arithmetic.
+-------------------------------------------------------------------------------
+*/
+flag floatx80_lt_quiet( floatx80 a, floatx80 b )
+{
+    flag aSign, bSign;
+
+    if (    (    ( extractFloatx80Exp( a ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( a )<<1 ) )
+         || (    ( extractFloatx80Exp( b ) == 0x7FFF )
+              && (bits64) ( extractFloatx80Frac( b )<<1 ) )
+       ) {
+        /* Do nothing, even if NaN as we're quiet */
+        return 0;
+    }
+    aSign = extractFloatx80Sign( a );
+    bSign = extractFloatx80Sign( b );
+    if ( aSign != bSign ) {
+        return
+               aSign
+            && (    ( ( (bits16) ( ( a.high | b.high )<<1 ) ) | a.low | b.low )
+                 != 0 );
+    }
+    return
+          aSign ? lt128( b.high, b.low, a.high, a.low )
+        : lt128( a.high, a.low, b.high, b.low );
+
+}
+
+#endif
+
diff --git a/arch/arm/nwfpe/softfloat.h b/arch/arm/nwfpe/softfloat.h
new file mode 100644
index 000000000..df4d243a2
--- /dev/null
+++ b/arch/arm/nwfpe/softfloat.h
@@ -0,0 +1,281 @@
+
+/*
+===============================================================================
+
+This C header file is part of the SoftFloat IEC/IEEE Floating-point
+Arithmetic Package, Release 2.
+
+Written by John R. Hauser.  This work was made possible in part by the
+International Computer Science Institute, located at Suite 600, 1947 Center
+Street, Berkeley, California 94704.  Funding was partially provided by the
+National Science Foundation under grant MIP-9311980.  The original version
+of this code was written as part of a project to build a fixed-point vector
+processor in collaboration with the University of California at Berkeley,
+overseen by Profs. Nelson Morgan and John Wawrzynek.  More information
+is available through the Web page
+http://www.jhauser.us/arithmetic/SoftFloat-2b/SoftFloat-source.txt
+
+THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort
+has been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT
+TIMES RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO
+PERSONS AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ANY
+AND ALL LOSSES, COSTS, OR OTHER PROBLEMS ARISING FROM ITS USE.
+
+Derivative works are acceptable, even for commercial purposes, so long as
+(1) they include prominent notice that the work is derivative, and (2) they
+include prominent notice akin to these three paragraphs for those parts of
+this code that are retained.
+
+===============================================================================
+*/
+
+#ifndef __SOFTFLOAT_H__
+#define __SOFTFLOAT_H__
+
+
+/*
+-------------------------------------------------------------------------------
+The macro `FLOATX80' must be defined to enable the extended double-precision
+floating-point format `floatx80'.  If this macro is not defined, the
+`floatx80' type will not be defined, and none of the functions that either
+input or output the `floatx80' type will be defined.
+-------------------------------------------------------------------------------
+*/
+#ifdef CONFIG_FPE_NWFPE_XP
+#define FLOATX80
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE floating-point types.
+-------------------------------------------------------------------------------
+*/
+typedef u32 float32;
+typedef u64 float64;
+typedef struct {
+#ifdef __ARMEB__
+    u16 __padding;
+    u16 high;
+#else
+    u16 high;
+    u16 __padding;
+#endif
+    u64 low;
+}  __attribute__ ((packed,aligned(4))) floatx80;
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE floating-point underflow tininess-detection mode.
+-------------------------------------------------------------------------------
+*/
+extern signed char float_detect_tininess;
+enum {
+    float_tininess_after_rounding  = 0,
+    float_tininess_before_rounding = 1
+};
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE floating-point rounding mode.
+-------------------------------------------------------------------------------
+*/
+//extern int8 float_rounding_mode;
+enum {
+    float_round_nearest_even = 0,
+    float_round_to_zero      = 1,
+    float_round_down         = 2,
+    float_round_up           = 3
+};
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE floating-point exception flags.
+-------------------------------------------------------------------------------
+enum {
+    float_flag_inexact   =  1,
+    float_flag_underflow =  2,
+    float_flag_overflow  =  4,
+    float_flag_divbyzero =  8,
+    float_flag_invalid   = 16
+};
+
+ScottB: November 4, 1998
+Changed the enumeration to match the bit order in the FPA11.
+*/
+
+enum {
+    float_flag_invalid   =  1,
+    float_flag_divbyzero =  2,
+    float_flag_overflow  =  4,
+    float_flag_underflow =  8,
+    float_flag_inexact   = 16
+};
+
+/*
+-------------------------------------------------------------------------------
+Routine to raise any or all of the software IEC/IEEE floating-point
+exception flags.
+-------------------------------------------------------------------------------
+*/
+void float_raise( signed char );
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE integer-to-floating-point conversion routines.
+-------------------------------------------------------------------------------
+*/
+float32 int32_to_float32( struct roundingData *, signed int );
+float64 int32_to_float64( signed int );
+#ifdef FLOATX80
+floatx80 int32_to_floatx80( signed int );
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE single-precision conversion routines.
+-------------------------------------------------------------------------------
+*/
+signed int float32_to_int32( struct roundingData *, float32 );
+signed int float32_to_int32_round_to_zero( float32 );
+float64 float32_to_float64( float32 );
+#ifdef FLOATX80
+floatx80 float32_to_floatx80( float32 );
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE single-precision operations.
+-------------------------------------------------------------------------------
+*/
+float32 float32_round_to_int( struct roundingData*, float32 );
+float32 float32_add( struct roundingData *, float32, float32 );
+float32 float32_sub( struct roundingData *, float32, float32 );
+float32 float32_mul( struct roundingData *, float32, float32 );
+float32 float32_div( struct roundingData *, float32, float32 );
+float32 float32_rem( struct roundingData *, float32, float32 );
+float32 float32_sqrt( struct roundingData*, float32 );
+char float32_eq( float32, float32 );
+char float32_le( float32, float32 );
+char float32_lt( float32, float32 );
+char float32_eq_signaling( float32, float32 );
+char float32_le_quiet( float32, float32 );
+char float32_lt_quiet( float32, float32 );
+char float32_is_signaling_nan( float32 );
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE double-precision conversion routines.
+-------------------------------------------------------------------------------
+*/
+signed int float64_to_int32( struct roundingData *, float64 );
+signed int float64_to_int32_round_to_zero( float64 );
+float32 float64_to_float32( struct roundingData *, float64 );
+#ifdef FLOATX80
+floatx80 float64_to_floatx80( float64 );
+#endif
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE double-precision operations.
+-------------------------------------------------------------------------------
+*/
+float64 float64_round_to_int( struct roundingData *, float64 );
+float64 float64_add( struct roundingData *, float64, float64 );
+float64 float64_sub( struct roundingData *, float64, float64 );
+float64 float64_mul( struct roundingData *, float64, float64 );
+float64 float64_div( struct roundingData *, float64, float64 );
+float64 float64_rem( struct roundingData *, float64, float64 );
+float64 float64_sqrt( struct roundingData *, float64 );
+char float64_eq( float64, float64 );
+char float64_le( float64, float64 );
+char float64_lt( float64, float64 );
+char float64_eq_signaling( float64, float64 );
+char float64_le_quiet( float64, float64 );
+char float64_lt_quiet( float64, float64 );
+char float64_is_signaling_nan( float64 );
+
+#ifdef FLOATX80
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE extended double-precision conversion routines.
+-------------------------------------------------------------------------------
+*/
+signed int floatx80_to_int32( struct roundingData *, floatx80 );
+signed int floatx80_to_int32_round_to_zero( floatx80 );
+float32 floatx80_to_float32( struct roundingData *, floatx80 );
+float64 floatx80_to_float64( struct roundingData *, floatx80 );
+
+/*
+-------------------------------------------------------------------------------
+Software IEC/IEEE extended double-precision operations.
+-------------------------------------------------------------------------------
+*/
+floatx80 floatx80_round_to_int( struct roundingData *, floatx80 );
+floatx80 floatx80_add( struct roundingData *, floatx80, floatx80 );
+floatx80 floatx80_sub( struct roundingData *, floatx80, floatx80 );
+floatx80 floatx80_mul( struct roundingData *, floatx80, floatx80 );
+floatx80 floatx80_div( struct roundingData *, floatx80, floatx80 );
+floatx80 floatx80_rem( struct roundingData *, floatx80, floatx80 );
+floatx80 floatx80_sqrt( struct roundingData *, floatx80 );
+char floatx80_eq( floatx80, floatx80 );
+char floatx80_le( floatx80, floatx80 );
+char floatx80_lt( floatx80, floatx80 );
+char floatx80_eq_signaling( floatx80, floatx80 );
+char floatx80_le_quiet( floatx80, floatx80 );
+char floatx80_lt_quiet( floatx80, floatx80 );
+char floatx80_is_signaling_nan( floatx80 );
+
+extern flag floatx80_is_nan(floatx80);
+
+#endif
+
+static inline flag extractFloat32Sign(float32 a)
+{
+	return a >> 31;
+}
+
+static inline flag float32_eq_nocheck(float32 a, float32 b)
+{
+	return (a == b) || ((bits32) ((a | b) << 1) == 0);
+}
+
+static inline flag float32_lt_nocheck(float32 a, float32 b)
+{
+	flag aSign, bSign;
+
+	aSign = extractFloat32Sign(a);
+	bSign = extractFloat32Sign(b);
+	if (aSign != bSign)
+		return aSign && ((bits32) ((a | b) << 1) != 0);
+	return (a != b) && (aSign ^ (a < b));
+}
+
+static inline flag extractFloat64Sign(float64 a)
+{
+	return a >> 63;
+}
+
+static inline flag float64_eq_nocheck(float64 a, float64 b)
+{
+	return (a == b) || ((bits64) ((a | b) << 1) == 0);
+}
+
+static inline flag float64_lt_nocheck(float64 a, float64 b)
+{
+	flag aSign, bSign;
+
+	aSign = extractFloat64Sign(a);
+	bSign = extractFloat64Sign(b);
+	if (aSign != bSign)
+		return aSign && ((bits64) ((a | b) << 1) != 0);
+	return (a != b) && (aSign ^ (a < b));
+}
+
+extern flag float32_is_nan( float32 a );
+extern flag float64_is_nan( float64 a );
+
+extern int32 float64_to_uint32( struct roundingData *roundData, float64 a );
+extern int32 float64_to_uint32_round_to_zero( float64 a );
+
+#endif