Merge tag 'trace-v7.3-rc2' of git://git.kernel.org/pub/scm/linux/kernel/git/trace/linux-tracegrafted

Pull tracing fixes from Steven Rostedt:

 - Don't destroy user event fields when removal fails

   User event fields are destroyed before the event is removed from
   visibility. But that can fail leaving the still visible event with no
   fields. Move the destroying of the fields to after the event is
   successfully removed from visibility.

 - Initialize function graph state is fork before calling
   copy_exec_state()

   For non-CLONE_VM forks, copy_exec_state() allocates a new
   task_exec_state. If that allocation fails, ftrace_graph_exit_task()
   will free the tasks ret_stack pointer. Since that pointer is still
   using the parent's ret_stack, it mistakenly frees the parent's
   pointer too.

   Call ftrace_graph_init() on the task first which will NULL out the
   new tasks's ret_stack and if the copy fails, it will not free
   anything.

 - Remove FGRAPH_MAX_INDEX

   The macro FGRAPH_MAX_INDEX was added but never used. Remove it.

 - Save ent_size in function graph printing of nested functions

   The function graph tracer needs to look at the next event to see if
   the next event is the return of the current function entry. If it is,
   it prints a single line:

	ktime_get();

   Otherwise it prints it like a nested function:

	tick_nohz_irq_exit() {
	    ktime_get();
	    kcpustat_irq_exit();
	}

   In order to look at the next event, it must save the current event so
   that it has the information to print from it. It saves the event in
   the iterator descriptor called "ent". What it doesn't save is the
   ent_size of the event which is now used to know if the function graph
   arguments are to be printed. The peek doesn't save the size so the
   size used happens to be that of the size of the last event that was
   seen.

   Save the entry event size in the iterator descriptor so that the
   correct size is used.

 - Fix several errors with freeing data in the histogram code

   The histogram code had a lot of leaked or or incorrect accounting
   when failures happen. Correct them.

 - Fix histogram regression of .percent and .graph modifiers

   Up until 6.3 histogram values could have "percent" or "graph"
   modifiers that changed how they were printed. But a change that added
   restricting histograms values from being strings, stack traces and
   other modifiers inadvertently prevented them from using the percent
   and graph modifiers, which were legal use cases for values.

   Put back the percent and graph modifiers.

 - Fix various typos in the comments

 - Set the trace_clock before initializing a histogram with clock
   argument

   The histogram API allows the user to specific which trace clock to
   use via a "clock=" string. The histogram is set up first before the
   clock is checked. If the passed in clock is not valid, it exits
   without fully fixing up the histogram leaving it on the list and a
   use-after-free can trigger.

   Update the clock argument first and if it fails then exit gracefully
   before the histogram trigger is placed on any lists.

 - Restore :mod: trailer after parsing in ftrace_set_clr_event

   The function ftrace_set_clr_event() modifies the parse string and
   needs to put it back to what was passed in. It searches for ":mod:"
   via a strsep() but fails to put back the first ':' in the string.

   Add back the ':' in the passed in string.

 - Take trace_array reference when opening a tracer options file

   The options files are dynamically created and some tracers add their
   own options. When a tracer adds their own list of options, the
   trace_array holding them has an array to hold the list of options for
   each tracer. This array increases in size via a krealloc(), and the
   new entry gets a newly allocated array to hold the options of the new
   tracer being added.

   The element in each entry of the tracer's option array holds a
   pointer back to the trace_array, a pointer to the tracer it is
   associated to, a pointer to the flags of the option.

   The issue is that these arrays are freed when the trace_array is
   freed when its instance it represents is removed from the instances
   directory. There's a race that an open of one of these options files
   can happen when the instance is being removed.

   Add a new helper function to be called by the open function of the
   options file to iterate all existing trace_arrays under a lock and
   find the one that has the given option element in one of it's tracer
   arrays. If found, then update the associated trace_array's reference
   counter to keep it from being freed. If not found, have the open call
   return -ENODEV.

 - Disable interrupts when acquiring the lock in rb_wake_up_waiters()

   The function rb_wake_up_waiters() assumes it will be called in
   interrupt context and does not disable irqs when taking
   cpu_buffer->reader_lock, which can be called in hard interrupt
   context. The issue is in PREEMPT_RT, this function is called in
   thread context leaving this lock open to a deadlock.

   Take the lock with interrupts disabled.

 - Use rcu_assign_pointer() for tmp_ops filter hash

   The tmp_ops used in update_ftrace_direct_mod() assigns its
   filter_hash field directly, but that field is annotated as __rcu and
   sparse complains. Assign it with rcu_assign_pointer()

 - Fix use-after-free in enable_trigger_private_data_free()

   The trace_event_call is accessed through the event_trigger_data's
   trace_event_file pointer to put the trace_event_call on freeing. The
   issue is that the trace_event_file data may have been freed already
   causing a use-after-free. Add a field to the event_trigger_data that
   points directly to the trace_event_call so that it can decrement its
   reference directly without needing to go through the
   trace_event_file.

 - Fix accounting of buffer data remote headers

   trace_buffer_desc_size() and trace_remote_alloc_buffer() undercount
   the number of pages is needed for the asked for size as it doesn't
   take into account the meta data on each page. Add a helper function
   to do the calculation properly and use that in these functions.

 - Catch nr_page_va overflow in ring_buffer_desc sizing

   The number of pages per remote ring buffer is capped by
   ring_buffer_desc::nr_page_va (32 bits). A buffer_size large enough to
   overflow that field would silently allocate a descriptor smaller than
   what was asked for.

 - Do not resize the subbuf order if any per_cpu buffer is disabled

   The mmapping of ring buffers disables resizing the subbuffers, but it
   is done per-cpu whereas the subbuf size change is done for all the
   per_cpu buffers under the buffer->mutex. It could change the size of
   some while the mapping is happening on others. Have the resize of the
   subbuf order check all the per_cpu buffers under the lock to see if
   any of them is disabled before starting and causing an inconsistency
   between buffers that are being mapped.

* tag 'trace-v7.3-rc2' of git://git.kernel.org/pub/scm/linux/kernel/git/trace/linux-trace: (25 commits)
  ring-buffer: Check resize_disabled before publishing the new subbuf order
  tracing/remotes: Catch nr_page_va overflow in ring_buffer_desc sizing
  tracing/remotes: Account for ring buffer page header in size calculation
  tracing: Don't dereference trace_event_file in deferred trigger free
  ftrace: Use rcu_assign_pointer() for tmp_ops filter hash
  ring-buffer: Acquire the lock with irqsave in rb_wake_up_waiters()
  tracing: Take trace_array reference when opening a tracer options file
  tracing: Fix ring_buffer_read_page_size() kernel-doc
  tracing: Restore :mod: trailer after parsing in ftrace_set_clr_event()
  tracing: Fix memory corruption from a "STACKTRACE" histogram key
  tracing: Fix memory corruption from the histogram stacktrace modifier
  tracing: Undo the registration when enabling the histogram trigger fails
  tracing: Take the reference before publishing the named histogram trigger
  tracing: Set the trace clock before registering the histogram trigger
  tracing: Fix typo "preceeded" in comment
  tracing: Fix typo "availabe" in comment
  tracing: Let histogram values keep the percent and graph modifiers
  tracing: Keep the entry count when the histogram stats allocation fails
  tracing: Free histogram the field rejected for a bad modifier
  tracing: Free histogram the var ref when its initialization fails
  ...

9 files changed, 3328 insertions, 0 deletions

diff --git a/include/math-emu/double.h b/include/math-emu/double.h
new file mode 100644
index 000000000..655ccf1c4
--- /dev/null
+++ b/include/math-emu/double.h
@@ -0,0 +1,205 @@
+/* Software floating-point emulation.
+   Definitions for IEEE Double Precision
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef    __MATH_EMU_DOUBLE_H__
+#define    __MATH_EMU_DOUBLE_H__
+
+#if _FP_W_TYPE_SIZE < 32
+#error "Here's a nickel kid.  Go buy yourself a real computer."
+#endif
+
+#if _FP_W_TYPE_SIZE < 64
+#define _FP_FRACTBITS_D		(2 * _FP_W_TYPE_SIZE)
+#else
+#define _FP_FRACTBITS_D		_FP_W_TYPE_SIZE
+#endif
+
+#define _FP_FRACBITS_D		53
+#define _FP_FRACXBITS_D		(_FP_FRACTBITS_D - _FP_FRACBITS_D)
+#define _FP_WFRACBITS_D		(_FP_WORKBITS + _FP_FRACBITS_D)
+#define _FP_WFRACXBITS_D	(_FP_FRACTBITS_D - _FP_WFRACBITS_D)
+#define _FP_EXPBITS_D		11
+#define _FP_EXPBIAS_D		1023
+#define _FP_EXPMAX_D		2047
+
+#define _FP_QNANBIT_D		\
+	((_FP_W_TYPE)1 << (_FP_FRACBITS_D-2) % _FP_W_TYPE_SIZE)
+#define _FP_IMPLBIT_D		\
+	((_FP_W_TYPE)1 << (_FP_FRACBITS_D-1) % _FP_W_TYPE_SIZE)
+#define _FP_OVERFLOW_D		\
+	((_FP_W_TYPE)1 << _FP_WFRACBITS_D % _FP_W_TYPE_SIZE)
+
+#if _FP_W_TYPE_SIZE < 64
+
+union _FP_UNION_D
+{
+  double flt;
+  struct {
+#if __BYTE_ORDER == __BIG_ENDIAN
+    unsigned sign  : 1;
+    unsigned exp   : _FP_EXPBITS_D;
+    unsigned frac1 : _FP_FRACBITS_D - (_FP_IMPLBIT_D != 0) - _FP_W_TYPE_SIZE;
+    unsigned frac0 : _FP_W_TYPE_SIZE;
+#else
+    unsigned frac0 : _FP_W_TYPE_SIZE;
+    unsigned frac1 : _FP_FRACBITS_D - (_FP_IMPLBIT_D != 0) - _FP_W_TYPE_SIZE;
+    unsigned exp   : _FP_EXPBITS_D;
+    unsigned sign  : 1;
+#endif
+  } bits __attribute__((packed));
+};
+
+#define FP_DECL_D(X)		_FP_DECL(2,X)
+#define FP_UNPACK_RAW_D(X,val)	_FP_UNPACK_RAW_2(D,X,val)
+#define FP_UNPACK_RAW_DP(X,val)	_FP_UNPACK_RAW_2_P(D,X,val)
+#define FP_PACK_RAW_D(val,X)	_FP_PACK_RAW_2(D,val,X)
+#define FP_PACK_RAW_DP(val,X)		\
+  do {					\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_2_P(D,val,X);	\
+  } while (0)
+
+#define FP_UNPACK_D(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_2(D,X,val);		\
+    _FP_UNPACK_CANONICAL(D,2,X);	\
+  } while (0)
+
+#define FP_UNPACK_DP(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_2_P(D,X,val);	\
+    _FP_UNPACK_CANONICAL(D,2,X);	\
+  } while (0)
+
+#define FP_PACK_D(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(D,2,X);		\
+    _FP_PACK_RAW_2(D,val,X);		\
+  } while (0)
+
+#define FP_PACK_DP(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(D,2,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_2_P(D,val,X);	\
+  } while (0)
+
+#define FP_ISSIGNAN_D(X)		_FP_ISSIGNAN(D,2,X)
+#define FP_NEG_D(R,X)			_FP_NEG(D,2,R,X)
+#define FP_ADD_D(R,X,Y)			_FP_ADD(D,2,R,X,Y)
+#define FP_SUB_D(R,X,Y)			_FP_SUB(D,2,R,X,Y)
+#define FP_MUL_D(R,X,Y)			_FP_MUL(D,2,R,X,Y)
+#define FP_DIV_D(R,X,Y)			_FP_DIV(D,2,R,X,Y)
+#define FP_SQRT_D(R,X)			_FP_SQRT(D,2,R,X)
+#define _FP_SQRT_MEAT_D(R,S,T,X,Q)	_FP_SQRT_MEAT_2(R,S,T,X,Q)
+
+#define FP_CMP_D(r,X,Y,un)	_FP_CMP(D,2,r,X,Y,un)
+#define FP_CMP_EQ_D(r,X,Y)	_FP_CMP_EQ(D,2,r,X,Y)
+
+#define FP_TO_INT_D(r,X,rsz,rsg)	_FP_TO_INT(D,2,r,X,rsz,rsg)
+#define FP_TO_INT_ROUND_D(r,X,rsz,rsg)	_FP_TO_INT_ROUND(D,2,r,X,rsz,rsg)
+#define FP_FROM_INT_D(X,r,rs,rt)	_FP_FROM_INT(D,2,X,r,rs,rt)
+
+#define _FP_FRAC_HIGH_D(X)	_FP_FRAC_HIGH_2(X)
+#define _FP_FRAC_HIGH_RAW_D(X)	_FP_FRAC_HIGH_2(X)
+
+#else
+
+union _FP_UNION_D
+{
+  double flt;
+  struct {
+#if __BYTE_ORDER == __BIG_ENDIAN
+    unsigned sign : 1;
+    unsigned exp  : _FP_EXPBITS_D;
+    unsigned long frac : _FP_FRACBITS_D - (_FP_IMPLBIT_D != 0);
+#else
+    unsigned long frac : _FP_FRACBITS_D - (_FP_IMPLBIT_D != 0);
+    unsigned exp  : _FP_EXPBITS_D;
+    unsigned sign : 1;
+#endif
+  } bits __attribute__((packed));
+};
+
+#define FP_DECL_D(X)		_FP_DECL(1,X)
+#define FP_UNPACK_RAW_D(X,val)	_FP_UNPACK_RAW_1(D,X,val)
+#define FP_UNPACK_RAW_DP(X,val)	_FP_UNPACK_RAW_1_P(D,X,val)
+#define FP_PACK_RAW_D(val,X)	_FP_PACK_RAW_1(D,val,X)
+#define FP_PACK_RAW_DP(val,X)		\
+  do {					\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_1_P(D,val,X);	\
+  } while (0)
+
+#define FP_UNPACK_D(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_1(D,X,val);		\
+    _FP_UNPACK_CANONICAL(D,1,X);	\
+  } while (0)
+
+#define FP_UNPACK_DP(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_1_P(D,X,val);	\
+    _FP_UNPACK_CANONICAL(D,1,X);	\
+  } while (0)
+
+#define FP_PACK_D(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(D,1,X);		\
+    _FP_PACK_RAW_1(D,val,X);		\
+  } while (0)
+
+#define FP_PACK_DP(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(D,1,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_1_P(D,val,X);	\
+  } while (0)
+
+#define FP_ISSIGNAN_D(X)		_FP_ISSIGNAN(D,1,X)
+#define FP_NEG_D(R,X)			_FP_NEG(D,1,R,X)
+#define FP_ADD_D(R,X,Y)			_FP_ADD(D,1,R,X,Y)
+#define FP_SUB_D(R,X,Y)			_FP_SUB(D,1,R,X,Y)
+#define FP_MUL_D(R,X,Y)			_FP_MUL(D,1,R,X,Y)
+#define FP_DIV_D(R,X,Y)			_FP_DIV(D,1,R,X,Y)
+#define FP_SQRT_D(R,X)			_FP_SQRT(D,1,R,X)
+#define _FP_SQRT_MEAT_D(R,S,T,X,Q)	_FP_SQRT_MEAT_1(R,S,T,X,Q)
+
+/* The implementation of _FP_MUL_D and _FP_DIV_D should be chosen by
+   the target machine.  */
+
+#define FP_CMP_D(r,X,Y,un)	_FP_CMP(D,1,r,X,Y,un)
+#define FP_CMP_EQ_D(r,X,Y)	_FP_CMP_EQ(D,1,r,X,Y)
+
+#define FP_TO_INT_D(r,X,rsz,rsg)	_FP_TO_INT(D,1,r,X,rsz,rsg)
+#define FP_TO_INT_ROUND_D(r,X,rsz,rsg)	_FP_TO_INT_ROUND(D,1,r,X,rsz,rsg)
+#define FP_FROM_INT_D(X,r,rs,rt)	_FP_FROM_INT(D,1,X,r,rs,rt)
+
+#define _FP_FRAC_HIGH_D(X)	_FP_FRAC_HIGH_1(X)
+#define _FP_FRAC_HIGH_RAW_D(X)	_FP_FRAC_HIGH_1(X)
+
+#endif /* W_TYPE_SIZE < 64 */
+
+
+#endif /* __MATH_EMU_DOUBLE_H__ */
diff --git a/include/math-emu/op-1.h b/include/math-emu/op-1.h
new file mode 100644
index 000000000..3be3bb422
--- /dev/null
+++ b/include/math-emu/op-1.h
@@ -0,0 +1,303 @@
+/* Software floating-point emulation.
+   Basic one-word fraction declaration and manipulation.
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef    __MATH_EMU_OP_1_H__
+#define    __MATH_EMU_OP_1_H__
+
+#define _FP_FRAC_DECL_1(X)	_FP_W_TYPE X##_f=0
+#define _FP_FRAC_COPY_1(D,S)	(D##_f = S##_f)
+#define _FP_FRAC_SET_1(X,I)	(X##_f = I)
+#define _FP_FRAC_HIGH_1(X)	(X##_f)
+#define _FP_FRAC_LOW_1(X)	(X##_f)
+#define _FP_FRAC_WORD_1(X,w)	(X##_f)
+
+#define _FP_FRAC_ADDI_1(X,I)	(X##_f += I)
+#define _FP_FRAC_SLL_1(X,N)			\
+  do {						\
+    if (__builtin_constant_p(N) && (N) == 1)	\
+      X##_f += X##_f;				\
+    else					\
+      X##_f <<= (N);				\
+  } while (0)
+#define _FP_FRAC_SRL_1(X,N)	(X##_f >>= N)
+
+/* Right shift with sticky-lsb.  */
+#define _FP_FRAC_SRS_1(X,N,sz)	__FP_FRAC_SRS_1(X##_f, N, sz)
+
+#define __FP_FRAC_SRS_1(X,N,sz)						\
+   (X = (X >> (N) | (__builtin_constant_p(N) && (N) == 1		\
+		     ? X & 1 : (X << (_FP_W_TYPE_SIZE - (N))) != 0)))
+
+#define _FP_FRAC_ADD_1(R,X,Y)	(R##_f = X##_f + Y##_f)
+#define _FP_FRAC_SUB_1(R,X,Y)	(R##_f = X##_f - Y##_f)
+#define _FP_FRAC_DEC_1(X,Y)	(X##_f -= Y##_f)
+#define _FP_FRAC_CLZ_1(z, X)	__FP_CLZ(z, X##_f)
+
+/* Predicates */
+#define _FP_FRAC_NEGP_1(X)	((_FP_WS_TYPE)X##_f < 0)
+#define _FP_FRAC_ZEROP_1(X)	(X##_f == 0)
+#define _FP_FRAC_OVERP_1(fs,X)	(X##_f & _FP_OVERFLOW_##fs)
+#define _FP_FRAC_CLEAR_OVERP_1(fs,X)	(X##_f &= ~_FP_OVERFLOW_##fs)
+#define _FP_FRAC_EQ_1(X, Y)	(X##_f == Y##_f)
+#define _FP_FRAC_GE_1(X, Y)	(X##_f >= Y##_f)
+#define _FP_FRAC_GT_1(X, Y)	(X##_f > Y##_f)
+
+#define _FP_ZEROFRAC_1		0
+#define _FP_MINFRAC_1		1
+#define _FP_MAXFRAC_1		(~(_FP_WS_TYPE)0)
+
+/*
+ * Unpack the raw bits of a native fp value.  Do not classify or
+ * normalize the data.
+ */
+
+#define _FP_UNPACK_RAW_1(fs, X, val)				\
+  do {								\
+    union _FP_UNION_##fs _flo; _flo.flt = (val);		\
+								\
+    X##_f = _flo.bits.frac;					\
+    X##_e = _flo.bits.exp;					\
+    X##_s = _flo.bits.sign;					\
+  } while (0)
+
+#define _FP_UNPACK_RAW_1_P(fs, X, val)				\
+  do {								\
+    union _FP_UNION_##fs *_flo =				\
+      (union _FP_UNION_##fs *)(val);				\
+								\
+    X##_f = _flo->bits.frac;					\
+    X##_e = _flo->bits.exp;					\
+    X##_s = _flo->bits.sign;					\
+  } while (0)
+
+/*
+ * Repack the raw bits of a native fp value.
+ */
+
+#define _FP_PACK_RAW_1(fs, val, X)				\
+  do {								\
+    union _FP_UNION_##fs _flo;					\
+								\
+    _flo.bits.frac = X##_f;					\
+    _flo.bits.exp  = X##_e;					\
+    _flo.bits.sign = X##_s;					\
+								\
+    (val) = _flo.flt;						\
+  } while (0)
+
+#define _FP_PACK_RAW_1_P(fs, val, X)				\
+  do {								\
+    union _FP_UNION_##fs *_flo =				\
+      (union _FP_UNION_##fs *)(val);				\
+								\
+    _flo->bits.frac = X##_f;					\
+    _flo->bits.exp  = X##_e;					\
+    _flo->bits.sign = X##_s;					\
+  } while (0)
+
+
+/*
+ * Multiplication algorithms:
+ */
+
+/* Basic.  Assuming the host word size is >= 2*FRACBITS, we can do the
+   multiplication immediately.  */
+
+#define _FP_MUL_MEAT_1_imm(wfracbits, R, X, Y)				\
+  do {									\
+    R##_f = X##_f * Y##_f;						\
+    /* Normalize since we know where the msb of the multiplicands	\
+       were (bit B), we know that the msb of the of the product is	\
+       at either 2B or 2B-1.  */					\
+    _FP_FRAC_SRS_1(R, wfracbits-1, 2*wfracbits);			\
+  } while (0)
+
+/* Given a 1W * 1W => 2W primitive, do the extended multiplication.  */
+
+#define _FP_MUL_MEAT_1_wide(wfracbits, R, X, Y, doit)			\
+  do {									\
+    _FP_W_TYPE _Z_f0, _Z_f1;						\
+    doit(_Z_f1, _Z_f0, X##_f, Y##_f);					\
+    /* Normalize since we know where the msb of the multiplicands	\
+       were (bit B), we know that the msb of the of the product is	\
+       at either 2B or 2B-1.  */					\
+    _FP_FRAC_SRS_2(_Z, wfracbits-1, 2*wfracbits);			\
+    R##_f = _Z_f0;							\
+  } while (0)
+
+/* Finally, a simple widening multiply algorithm.  What fun!  */
+
+#define _FP_MUL_MEAT_1_hard(wfracbits, R, X, Y)				\
+  do {									\
+    _FP_W_TYPE _xh, _xl, _yh, _yl, _z_f0, _z_f1, _a_f0, _a_f1;		\
+									\
+    /* split the words in half */					\
+    _xh = X##_f >> (_FP_W_TYPE_SIZE/2);					\
+    _xl = X##_f & (((_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE/2)) - 1);		\
+    _yh = Y##_f >> (_FP_W_TYPE_SIZE/2);					\
+    _yl = Y##_f & (((_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE/2)) - 1);		\
+									\
+    /* multiply the pieces */						\
+    _z_f0 = _xl * _yl;							\
+    _a_f0 = _xh * _yl;							\
+    _a_f1 = _xl * _yh;							\
+    _z_f1 = _xh * _yh;							\
+									\
+    /* reassemble into two full words */				\
+    if ((_a_f0 += _a_f1) < _a_f1)					\
+      _z_f1 += (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE/2);			\
+    _a_f1 = _a_f0 >> (_FP_W_TYPE_SIZE/2);				\
+    _a_f0 = _a_f0 << (_FP_W_TYPE_SIZE/2);				\
+    _FP_FRAC_ADD_2(_z, _z, _a);						\
+									\
+    /* normalize */							\
+    _FP_FRAC_SRS_2(_z, wfracbits - 1, 2*wfracbits);			\
+    R##_f = _z_f0;							\
+  } while (0)
+
+
+/*
+ * Division algorithms:
+ */
+
+/* Basic.  Assuming the host word size is >= 2*FRACBITS, we can do the
+   division immediately.  Give this macro either _FP_DIV_HELP_imm for
+   C primitives or _FP_DIV_HELP_ldiv for the ISO function.  Which you
+   choose will depend on what the compiler does with divrem4.  */
+
+#define _FP_DIV_MEAT_1_imm(fs, R, X, Y, doit)		\
+  do {							\
+    _FP_W_TYPE _q, _r;					\
+    X##_f <<= (X##_f < Y##_f				\
+	       ? R##_e--, _FP_WFRACBITS_##fs		\
+	       : _FP_WFRACBITS_##fs - 1);		\
+    doit(_q, _r, X##_f, Y##_f);				\
+    R##_f = _q | (_r != 0);				\
+  } while (0)
+
+/* GCC's longlong.h defines a 2W / 1W => (1W,1W) primitive udiv_qrnnd
+   that may be useful in this situation.  This first is for a primitive
+   that requires normalization, the second for one that does not.  Look
+   for UDIV_NEEDS_NORMALIZATION to tell which your machine needs.  */
+
+#define _FP_DIV_MEAT_1_udiv_norm(fs, R, X, Y)				\
+  do {									\
+    _FP_W_TYPE _nh, _nl, _q, _r, _y;					\
+									\
+    /* Normalize Y -- i.e. make the most significant bit set.  */	\
+    _y = Y##_f << _FP_WFRACXBITS_##fs;					\
+									\
+    /* Shift X op correspondingly high, that is, up one full word.  */	\
+    if (X##_f < Y##_f)							\
+      {									\
+	R##_e--;							\
+	_nl = 0;							\
+	_nh = X##_f;							\
+      }									\
+    else								\
+      {									\
+	_nl = X##_f << (_FP_W_TYPE_SIZE - 1);				\
+	_nh = X##_f >> 1;						\
+      }									\
+    									\
+    udiv_qrnnd(_q, _r, _nh, _nl, _y);					\
+    R##_f = _q | (_r != 0);						\
+  } while (0)
+
+#define _FP_DIV_MEAT_1_udiv(fs, R, X, Y)		\
+  do {							\
+    _FP_W_TYPE _nh, _nl, _q, _r;			\
+    if (X##_f < Y##_f)					\
+      {							\
+	R##_e--;					\
+	_nl = X##_f << _FP_WFRACBITS_##fs;		\
+	_nh = X##_f >> _FP_WFRACXBITS_##fs;		\
+      }							\
+    else						\
+      {							\
+	_nl = X##_f << (_FP_WFRACBITS_##fs - 1);	\
+	_nh = X##_f >> (_FP_WFRACXBITS_##fs + 1);	\
+      }							\
+    udiv_qrnnd(_q, _r, _nh, _nl, Y##_f);		\
+    R##_f = _q | (_r != 0);				\
+  } while (0)
+  
+  
+/*
+ * Square root algorithms:
+ * We have just one right now, maybe Newton approximation
+ * should be added for those machines where division is fast.
+ */
+ 
+#define _FP_SQRT_MEAT_1(R, S, T, X, q)			\
+  do {							\
+    while (q != _FP_WORK_ROUND)				\
+      {							\
+        T##_f = S##_f + q;				\
+        if (T##_f <= X##_f)				\
+          {						\
+            S##_f = T##_f + q;				\
+            X##_f -= T##_f;				\
+            R##_f += q;					\
+          }						\
+        _FP_FRAC_SLL_1(X, 1);				\
+        q >>= 1;					\
+      }							\
+    if (X##_f)						\
+      {							\
+	if (S##_f < X##_f)				\
+	  R##_f |= _FP_WORK_ROUND;			\
+	R##_f |= _FP_WORK_STICKY;			\
+      }							\
+  } while (0)
+
+/*
+ * Assembly/disassembly for converting to/from integral types.  
+ * No shifting or overflow handled here.
+ */
+
+#define _FP_FRAC_ASSEMBLE_1(r, X, rsize)	(r = X##_f)
+#define _FP_FRAC_DISASSEMBLE_1(X, r, rsize)	(X##_f = r)
+
+
+/*
+ * Convert FP values between word sizes
+ */
+
+#define _FP_FRAC_CONV_1_1(dfs, sfs, D, S)				\
+  do {									\
+    D##_f = S##_f;							\
+    if (_FP_WFRACBITS_##sfs > _FP_WFRACBITS_##dfs)			\
+      {									\
+	if (S##_c != FP_CLS_NAN)					\
+	  _FP_FRAC_SRS_1(D, (_FP_WFRACBITS_##sfs-_FP_WFRACBITS_##dfs),	\
+			 _FP_WFRACBITS_##sfs);				\
+	else								\
+	  _FP_FRAC_SRL_1(D, (_FP_WFRACBITS_##sfs-_FP_WFRACBITS_##dfs));	\
+      }									\
+    else								\
+      D##_f <<= _FP_WFRACBITS_##dfs - _FP_WFRACBITS_##sfs;		\
+  } while (0)
+
+#endif /* __MATH_EMU_OP_1_H__ */
diff --git a/include/math-emu/op-2.h b/include/math-emu/op-2.h
new file mode 100644
index 000000000..244522b02
--- /dev/null
+++ b/include/math-emu/op-2.h
@@ -0,0 +1,605 @@
+/* Software floating-point emulation.
+   Basic two-word fraction declaration and manipulation.
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef __MATH_EMU_OP_2_H__
+#define __MATH_EMU_OP_2_H__
+
+#define _FP_FRAC_DECL_2(X)	_FP_W_TYPE X##_f0 = 0, X##_f1 = 0
+#define _FP_FRAC_COPY_2(D,S)	(D##_f0 = S##_f0, D##_f1 = S##_f1)
+#define _FP_FRAC_SET_2(X,I)	__FP_FRAC_SET_2(X, I)
+#define _FP_FRAC_HIGH_2(X)	(X##_f1)
+#define _FP_FRAC_LOW_2(X)	(X##_f0)
+#define _FP_FRAC_WORD_2(X,w)	(X##_f##w)
+#define _FP_FRAC_SLL_2(X, N) (						       \
+	(void) (((N) < _FP_W_TYPE_SIZE)					       \
+	  ? ({								       \
+		if (__builtin_constant_p(N) && (N) == 1) {		       \
+			X##_f1 = X##_f1 + X##_f1 +			       \
+				(((_FP_WS_TYPE) (X##_f0)) < 0);		       \
+			X##_f0 += X##_f0;				       \
+		} else {						       \
+			X##_f1 = X##_f1 << (N) | X##_f0 >>		       \
+						(_FP_W_TYPE_SIZE - (N));       \
+			X##_f0 <<= (N);					       \
+		}							       \
+		0;							       \
+	    })								       \
+	  : ({								       \
+	      X##_f1 = X##_f0 << ((N) - _FP_W_TYPE_SIZE);		       \
+	      X##_f0 = 0;						       \
+	  })))
+
+
+#define _FP_FRAC_SRL_2(X, N) (						       \
+	(void) (((N) < _FP_W_TYPE_SIZE)					       \
+	  ? ({								       \
+	      X##_f0 = X##_f0 >> (N) | X##_f1 << (_FP_W_TYPE_SIZE - (N));      \
+	      X##_f1 >>= (N);						       \
+	    })								       \
+	  : ({								       \
+	      X##_f0 = X##_f1 >> ((N) - _FP_W_TYPE_SIZE);		       \
+	      X##_f1 = 0;						       \
+	    })))
+
+
+/* Right shift with sticky-lsb.  */
+#define _FP_FRAC_SRS_2(X, N, sz) (					       \
+	(void) (((N) < _FP_W_TYPE_SIZE)					       \
+	  ? ({								       \
+	      X##_f0 = (X##_f1 << (_FP_W_TYPE_SIZE - (N)) | X##_f0 >> (N)      \
+			| (__builtin_constant_p(N) && (N) == 1		       \
+			   ? X##_f0 & 1					       \
+			   : (X##_f0 << (_FP_W_TYPE_SIZE - (N))) != 0));       \
+		X##_f1 >>= (N);						       \
+	    })								       \
+	  : ({								       \
+	      X##_f0 = (X##_f1 >> ((N) - _FP_W_TYPE_SIZE)		       \
+			| ((((N) == _FP_W_TYPE_SIZE			       \
+			     ? 0					       \
+			     : (X##_f1 << (2*_FP_W_TYPE_SIZE - (N))))          \
+			    | X##_f0) != 0));				       \
+	      X##_f1 = 0;						       \
+	    })))
+
+#define _FP_FRAC_ADDI_2(X,I)	\
+  __FP_FRAC_ADDI_2(X##_f1, X##_f0, I)
+
+#define _FP_FRAC_ADD_2(R,X,Y)	\
+  __FP_FRAC_ADD_2(R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
+
+#define _FP_FRAC_SUB_2(R,X,Y)	\
+  __FP_FRAC_SUB_2(R##_f1, R##_f0, X##_f1, X##_f0, Y##_f1, Y##_f0)
+
+#define _FP_FRAC_DEC_2(X,Y)	\
+  __FP_FRAC_DEC_2(X##_f1, X##_f0, Y##_f1, Y##_f0)
+
+#define _FP_FRAC_CLZ_2(R,X)	\
+  do {				\
+    if (X##_f1)			\
+      __FP_CLZ(R,X##_f1);	\
+    else 			\
+    {				\
+      __FP_CLZ(R,X##_f0);	\
+      R += _FP_W_TYPE_SIZE;	\
+    }				\
+  } while(0)
+
+/* Predicates */
+#define _FP_FRAC_NEGP_2(X)	((_FP_WS_TYPE)X##_f1 < 0)
+#define _FP_FRAC_ZEROP_2(X)	((X##_f1 | X##_f0) == 0)
+#define _FP_FRAC_OVERP_2(fs,X)	(_FP_FRAC_HIGH_##fs(X) & _FP_OVERFLOW_##fs)
+#define _FP_FRAC_CLEAR_OVERP_2(fs,X)	(_FP_FRAC_HIGH_##fs(X) &= ~_FP_OVERFLOW_##fs)
+#define _FP_FRAC_EQ_2(X, Y)	(X##_f1 == Y##_f1 && X##_f0 == Y##_f0)
+#define _FP_FRAC_GT_2(X, Y)	\
+  (X##_f1 > Y##_f1 || (X##_f1 == Y##_f1 && X##_f0 > Y##_f0))
+#define _FP_FRAC_GE_2(X, Y)	\
+  (X##_f1 > Y##_f1 || (X##_f1 == Y##_f1 && X##_f0 >= Y##_f0))
+
+#define _FP_ZEROFRAC_2		0, 0
+#define _FP_MINFRAC_2		0, 1
+#define _FP_MAXFRAC_2		(~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0)
+
+/*
+ * Internals 
+ */
+
+#define __FP_FRAC_SET_2(X,I1,I0)	(X##_f0 = I0, X##_f1 = I1)
+
+#define __FP_CLZ_2(R, xh, xl)	\
+  do {				\
+    if (xh)			\
+      __FP_CLZ(R,xh);		\
+    else 			\
+    {				\
+      __FP_CLZ(R,xl);		\
+      R += _FP_W_TYPE_SIZE;	\
+    }				\
+  } while(0)
+
+#if 0
+
+#ifndef __FP_FRAC_ADDI_2
+#define __FP_FRAC_ADDI_2(xh, xl, i)	\
+  (xh += ((xl += i) < i))
+#endif
+#ifndef __FP_FRAC_ADD_2
+#define __FP_FRAC_ADD_2(rh, rl, xh, xl, yh, yl)	\
+  (rh = xh + yh + ((rl = xl + yl) < xl))
+#endif
+#ifndef __FP_FRAC_SUB_2
+#define __FP_FRAC_SUB_2(rh, rl, xh, xl, yh, yl)	\
+  (rh = xh - yh - ((rl = xl - yl) > xl))
+#endif
+#ifndef __FP_FRAC_DEC_2
+#define __FP_FRAC_DEC_2(xh, xl, yh, yl)	\
+  do {					\
+    UWtype _t = xl;			\
+    xh -= yh + ((xl -= yl) > _t);	\
+  } while (0)
+#endif
+
+#else
+
+#undef __FP_FRAC_ADDI_2
+#define __FP_FRAC_ADDI_2(xh, xl, i)	add_ssaaaa(xh, xl, xh, xl, 0, i)
+#undef __FP_FRAC_ADD_2
+#define __FP_FRAC_ADD_2			add_ssaaaa
+#undef __FP_FRAC_SUB_2
+#define __FP_FRAC_SUB_2			sub_ddmmss
+#undef __FP_FRAC_DEC_2
+#define __FP_FRAC_DEC_2(xh, xl, yh, yl)	sub_ddmmss(xh, xl, xh, xl, yh, yl)
+
+#endif
+
+/*
+ * Unpack the raw bits of a native fp value.  Do not classify or
+ * normalize the data.
+ */
+
+#define _FP_UNPACK_RAW_2(fs, X, val)			\
+  do {							\
+    union _FP_UNION_##fs _flo; _flo.flt = (val);	\
+							\
+    X##_f0 = _flo.bits.frac0;				\
+    X##_f1 = _flo.bits.frac1;				\
+    X##_e  = _flo.bits.exp;				\
+    X##_s  = _flo.bits.sign;				\
+  } while (0)
+
+#define _FP_UNPACK_RAW_2_P(fs, X, val)			\
+  do {							\
+    union _FP_UNION_##fs *_flo =			\
+      (union _FP_UNION_##fs *)(val);			\
+							\
+    X##_f0 = _flo->bits.frac0;				\
+    X##_f1 = _flo->bits.frac1;				\
+    X##_e  = _flo->bits.exp;				\
+    X##_s  = _flo->bits.sign;				\
+  } while (0)
+
+
+/*
+ * Repack the raw bits of a native fp value.
+ */
+
+#define _FP_PACK_RAW_2(fs, val, X)			\
+  do {							\
+    union _FP_UNION_##fs _flo;				\
+							\
+    _flo.bits.frac0 = X##_f0;				\
+    _flo.bits.frac1 = X##_f1;				\
+    _flo.bits.exp   = X##_e;				\
+    _flo.bits.sign  = X##_s;				\
+							\
+    (val) = _flo.flt;					\
+  } while (0)
+
+#define _FP_PACK_RAW_2_P(fs, val, X)			\
+  do {							\
+    union _FP_UNION_##fs *_flo =			\
+      (union _FP_UNION_##fs *)(val);			\
+							\
+    _flo->bits.frac0 = X##_f0;				\
+    _flo->bits.frac1 = X##_f1;				\
+    _flo->bits.exp   = X##_e;				\
+    _flo->bits.sign  = X##_s;				\
+  } while (0)
+
+
+/*
+ * Multiplication algorithms:
+ */
+
+/* Given a 1W * 1W => 2W primitive, do the extended multiplication.  */
+
+#define _FP_MUL_MEAT_2_wide(wfracbits, R, X, Y, doit)			\
+  do {									\
+    _FP_FRAC_DECL_4(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c);	\
+									\
+    doit(_FP_FRAC_WORD_4(_z,1), _FP_FRAC_WORD_4(_z,0), X##_f0, Y##_f0);	\
+    doit(_b_f1, _b_f0, X##_f0, Y##_f1);					\
+    doit(_c_f1, _c_f0, X##_f1, Y##_f0);					\
+    doit(_FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2), X##_f1, Y##_f1);	\
+									\
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		    _FP_FRAC_WORD_4(_z,1), 0, _b_f1, _b_f0,		\
+		    _FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		    _FP_FRAC_WORD_4(_z,1));				\
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		    _FP_FRAC_WORD_4(_z,1), 0, _c_f1, _c_f0,		\
+		    _FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		    _FP_FRAC_WORD_4(_z,1));				\
+									\
+    /* Normalize since we know where the msb of the multiplicands	\
+       were (bit B), we know that the msb of the of the product is	\
+       at either 2B or 2B-1.  */					\
+    _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits);			\
+    R##_f0 = _FP_FRAC_WORD_4(_z,0);					\
+    R##_f1 = _FP_FRAC_WORD_4(_z,1);					\
+  } while (0)
+
+/* Given a 1W * 1W => 2W primitive, do the extended multiplication.
+   Do only 3 multiplications instead of four. This one is for machines
+   where multiplication is much more expensive than subtraction.  */
+
+#define _FP_MUL_MEAT_2_wide_3mul(wfracbits, R, X, Y, doit)		\
+  do {									\
+    _FP_FRAC_DECL_4(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c);	\
+    _FP_W_TYPE _d;							\
+    int _c1, _c2;							\
+									\
+    _b_f0 = X##_f0 + X##_f1;						\
+    _c1 = _b_f0 < X##_f0;						\
+    _b_f1 = Y##_f0 + Y##_f1;						\
+    _c2 = _b_f1 < Y##_f0;						\
+    doit(_d, _FP_FRAC_WORD_4(_z,0), X##_f0, Y##_f0);			\
+    doit(_FP_FRAC_WORD_4(_z,2), _FP_FRAC_WORD_4(_z,1), _b_f0, _b_f1);	\
+    doit(_c_f1, _c_f0, X##_f1, Y##_f1);					\
+									\
+    _b_f0 &= -_c2;							\
+    _b_f1 &= -_c1;							\
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		    _FP_FRAC_WORD_4(_z,1), (_c1 & _c2), 0, _d,		\
+		    0, _FP_FRAC_WORD_4(_z,2), _FP_FRAC_WORD_4(_z,1));	\
+    __FP_FRAC_ADDI_2(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		     _b_f0);						\
+    __FP_FRAC_ADDI_2(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		     _b_f1);						\
+    __FP_FRAC_DEC_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		    _FP_FRAC_WORD_4(_z,1),				\
+		    0, _d, _FP_FRAC_WORD_4(_z,0));			\
+    __FP_FRAC_DEC_3(_FP_FRAC_WORD_4(_z,3),_FP_FRAC_WORD_4(_z,2),	\
+		    _FP_FRAC_WORD_4(_z,1), 0, _c_f1, _c_f0);		\
+    __FP_FRAC_ADD_2(_FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2),	\
+		    _c_f1, _c_f0,					\
+		    _FP_FRAC_WORD_4(_z,3), _FP_FRAC_WORD_4(_z,2));	\
+									\
+    /* Normalize since we know where the msb of the multiplicands	\
+       were (bit B), we know that the msb of the of the product is	\
+       at either 2B or 2B-1.  */					\
+    _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits);			\
+    R##_f0 = _FP_FRAC_WORD_4(_z,0);					\
+    R##_f1 = _FP_FRAC_WORD_4(_z,1);					\
+  } while (0)
+
+#define _FP_MUL_MEAT_2_gmp(wfracbits, R, X, Y)				\
+  do {									\
+    _FP_FRAC_DECL_4(_z);						\
+    _FP_W_TYPE _x[2], _y[2];						\
+    _x[0] = X##_f0; _x[1] = X##_f1;					\
+    _y[0] = Y##_f0; _y[1] = Y##_f1;					\
+									\
+    mpn_mul_n(_z_f, _x, _y, 2);						\
+									\
+    /* Normalize since we know where the msb of the multiplicands	\
+       were (bit B), we know that the msb of the of the product is	\
+       at either 2B or 2B-1.  */					\
+    _FP_FRAC_SRS_4(_z, wfracbits-1, 2*wfracbits);			\
+    R##_f0 = _z_f[0];							\
+    R##_f1 = _z_f[1];							\
+  } while (0)
+
+/* Do at most 120x120=240 bits multiplication using double floating
+   point multiplication.  This is useful if floating point
+   multiplication has much bigger throughput than integer multiply.
+   It is supposed to work for _FP_W_TYPE_SIZE 64 and wfracbits
+   between 106 and 120 only.  
+   Caller guarantees that X and Y has (1LLL << (wfracbits - 1)) set.
+   SETFETZ is a macro which will disable all FPU exceptions and set rounding
+   towards zero,  RESETFE should optionally reset it back.  */
+
+#define _FP_MUL_MEAT_2_120_240_double(wfracbits, R, X, Y, setfetz, resetfe)	\
+  do {										\
+    static const double _const[] = {						\
+      /* 2^-24 */ 5.9604644775390625e-08,					\
+      /* 2^-48 */ 3.5527136788005009e-15,					\
+      /* 2^-72 */ 2.1175823681357508e-22,					\
+      /* 2^-96 */ 1.2621774483536189e-29,					\
+      /* 2^28 */ 2.68435456e+08,						\
+      /* 2^4 */ 1.600000e+01,							\
+      /* 2^-20 */ 9.5367431640625e-07,						\
+      /* 2^-44 */ 5.6843418860808015e-14,					\
+      /* 2^-68 */ 3.3881317890172014e-21,					\
+      /* 2^-92 */ 2.0194839173657902e-28,					\
+      /* 2^-116 */ 1.2037062152420224e-35};					\
+    double _a240, _b240, _c240, _d240, _e240, _f240, 				\
+	   _g240, _h240, _i240, _j240, _k240;					\
+    union { double d; UDItype i; } _l240, _m240, _n240, _o240,			\
+				   _p240, _q240, _r240, _s240;			\
+    UDItype _t240, _u240, _v240, _w240, _x240, _y240 = 0;			\
+										\
+    if (wfracbits < 106 || wfracbits > 120)					\
+      abort();									\
+										\
+    setfetz;									\
+										\
+    _e240 = (double)(long)(X##_f0 & 0xffffff);					\
+    _j240 = (double)(long)(Y##_f0 & 0xffffff);					\
+    _d240 = (double)(long)((X##_f0 >> 24) & 0xffffff);				\
+    _i240 = (double)(long)((Y##_f0 >> 24) & 0xffffff);				\
+    _c240 = (double)(long)(((X##_f1 << 16) & 0xffffff) | (X##_f0 >> 48));	\
+    _h240 = (double)(long)(((Y##_f1 << 16) & 0xffffff) | (Y##_f0 >> 48));	\
+    _b240 = (double)(long)((X##_f1 >> 8) & 0xffffff);				\
+    _g240 = (double)(long)((Y##_f1 >> 8) & 0xffffff);				\
+    _a240 = (double)(long)(X##_f1 >> 32);					\
+    _f240 = (double)(long)(Y##_f1 >> 32);					\
+    _e240 *= _const[3];								\
+    _j240 *= _const[3];								\
+    _d240 *= _const[2];								\
+    _i240 *= _const[2];								\
+    _c240 *= _const[1];								\
+    _h240 *= _const[1];								\
+    _b240 *= _const[0];								\
+    _g240 *= _const[0];								\
+    _s240.d =							      _e240*_j240;\
+    _r240.d =						_d240*_j240 + _e240*_i240;\
+    _q240.d =				  _c240*_j240 + _d240*_i240 + _e240*_h240;\
+    _p240.d =		    _b240*_j240 + _c240*_i240 + _d240*_h240 + _e240*_g240;\
+    _o240.d = _a240*_j240 + _b240*_i240 + _c240*_h240 + _d240*_g240 + _e240*_f240;\
+    _n240.d = _a240*_i240 + _b240*_h240 + _c240*_g240 + _d240*_f240;		\
+    _m240.d = _a240*_h240 + _b240*_g240 + _c240*_f240;				\
+    _l240.d = _a240*_g240 + _b240*_f240;					\
+    _k240 =   _a240*_f240;							\
+    _r240.d += _s240.d;								\
+    _q240.d += _r240.d;								\
+    _p240.d += _q240.d;								\
+    _o240.d += _p240.d;								\
+    _n240.d += _o240.d;								\
+    _m240.d += _n240.d;								\
+    _l240.d += _m240.d;								\
+    _k240 += _l240.d;								\
+    _s240.d -= ((_const[10]+_s240.d)-_const[10]);				\
+    _r240.d -= ((_const[9]+_r240.d)-_const[9]);					\
+    _q240.d -= ((_const[8]+_q240.d)-_const[8]);					\
+    _p240.d -= ((_const[7]+_p240.d)-_const[7]);					\
+    _o240.d += _const[7];							\
+    _n240.d += _const[6];							\
+    _m240.d += _const[5];							\
+    _l240.d += _const[4];							\
+    if (_s240.d != 0.0) _y240 = 1;						\
+    if (_r240.d != 0.0) _y240 = 1;						\
+    if (_q240.d != 0.0) _y240 = 1;						\
+    if (_p240.d != 0.0) _y240 = 1;						\
+    _t240 = (DItype)_k240;							\
+    _u240 = _l240.i;								\
+    _v240 = _m240.i;								\
+    _w240 = _n240.i;								\
+    _x240 = _o240.i;								\
+    R##_f1 = (_t240 << (128 - (wfracbits - 1)))					\
+	     | ((_u240 & 0xffffff) >> ((wfracbits - 1) - 104));			\
+    R##_f0 = ((_u240 & 0xffffff) << (168 - (wfracbits - 1)))			\
+    	     | ((_v240 & 0xffffff) << (144 - (wfracbits - 1)))			\
+    	     | ((_w240 & 0xffffff) << (120 - (wfracbits - 1)))			\
+    	     | ((_x240 & 0xffffff) >> ((wfracbits - 1) - 96))			\
+    	     | _y240;								\
+    resetfe;									\
+  } while (0)
+
+/*
+ * Division algorithms:
+ */
+
+#define _FP_DIV_MEAT_2_udiv(fs, R, X, Y)				\
+  do {									\
+    _FP_W_TYPE _n_f2, _n_f1, _n_f0, _r_f1, _r_f0, _m_f1, _m_f0;		\
+    if (_FP_FRAC_GT_2(X, Y))						\
+      {									\
+	_n_f2 = X##_f1 >> 1;						\
+	_n_f1 = X##_f1 << (_FP_W_TYPE_SIZE - 1) | X##_f0 >> 1;		\
+	_n_f0 = X##_f0 << (_FP_W_TYPE_SIZE - 1);			\
+      }									\
+    else								\
+      {									\
+	R##_e--;							\
+	_n_f2 = X##_f1;							\
+	_n_f1 = X##_f0;							\
+	_n_f0 = 0;							\
+      }									\
+									\
+    /* Normalize, i.e. make the most significant bit of the 		\
+       denominator set. */						\
+    _FP_FRAC_SLL_2(Y, _FP_WFRACXBITS_##fs);				\
+									\
+    udiv_qrnnd(R##_f1, _r_f1, _n_f2, _n_f1, Y##_f1);			\
+    umul_ppmm(_m_f1, _m_f0, R##_f1, Y##_f0);				\
+    _r_f0 = _n_f0;							\
+    if (_FP_FRAC_GT_2(_m, _r))						\
+      {									\
+	R##_f1--;							\
+	_FP_FRAC_ADD_2(_r, Y, _r);					\
+	if (_FP_FRAC_GE_2(_r, Y) && _FP_FRAC_GT_2(_m, _r))		\
+	  {								\
+	    R##_f1--;							\
+	    _FP_FRAC_ADD_2(_r, Y, _r);					\
+	  }								\
+      }									\
+    _FP_FRAC_DEC_2(_r, _m);						\
+									\
+    if (_r_f1 == Y##_f1)						\
+      {									\
+	/* This is a special case, not an optimization			\
+	   (_r/Y##_f1 would not fit into UWtype).			\
+	   As _r is guaranteed to be < Y,  R##_f0 can be either		\
+	   (UWtype)-1 or (UWtype)-2.  But as we know what kind		\
+	   of bits it is (sticky, guard, round),  we don't care.	\
+	   We also don't care what the reminder is,  because the	\
+	   guard bit will be set anyway.  -jj */			\
+	R##_f0 = -1;							\
+      }									\
+    else								\
+      {									\
+	udiv_qrnnd(R##_f0, _r_f1, _r_f1, _r_f0, Y##_f1);		\
+	umul_ppmm(_m_f1, _m_f0, R##_f0, Y##_f0);			\
+	_r_f0 = 0;							\
+	if (_FP_FRAC_GT_2(_m, _r))					\
+	  {								\
+	    R##_f0--;							\
+	    _FP_FRAC_ADD_2(_r, Y, _r);					\
+	    if (_FP_FRAC_GE_2(_r, Y) && _FP_FRAC_GT_2(_m, _r))		\
+	      {								\
+		R##_f0--;						\
+		_FP_FRAC_ADD_2(_r, Y, _r);				\
+	      }								\
+	  }								\
+	if (!_FP_FRAC_EQ_2(_r, _m))					\
+	  R##_f0 |= _FP_WORK_STICKY;					\
+      }									\
+  } while (0)
+
+
+#define _FP_DIV_MEAT_2_gmp(fs, R, X, Y)					\
+  do {									\
+    _FP_W_TYPE _x[4], _y[2], _z[4];					\
+    _y[0] = Y##_f0; _y[1] = Y##_f1;					\
+    _x[0] = _x[3] = 0;							\
+    if (_FP_FRAC_GT_2(X, Y))						\
+      {									\
+	R##_e++;							\
+	_x[1] = (X##_f0 << (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE) |	\
+		 X##_f1 >> (_FP_W_TYPE_SIZE -				\
+			    (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE)));	\
+	_x[2] = X##_f1 << (_FP_WFRACBITS_##fs-1 - _FP_W_TYPE_SIZE);	\
+      }									\
+    else								\
+      {									\
+	_x[1] = (X##_f0 << (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE) |	\
+		 X##_f1 >> (_FP_W_TYPE_SIZE -				\
+			    (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE)));	\
+	_x[2] = X##_f1 << (_FP_WFRACBITS_##fs - _FP_W_TYPE_SIZE);	\
+      }									\
+									\
+    (void) mpn_divrem (_z, 0, _x, 4, _y, 2);				\
+    R##_f1 = _z[1];							\
+    R##_f0 = _z[0] | ((_x[0] | _x[1]) != 0);				\
+  } while (0)
+
+
+/*
+ * Square root algorithms:
+ * We have just one right now, maybe Newton approximation
+ * should be added for those machines where division is fast.
+ */
+ 
+#define _FP_SQRT_MEAT_2(R, S, T, X, q)			\
+  do {							\
+    while (q)						\
+      {							\
+	T##_f1 = S##_f1 + q;				\
+	if (T##_f1 <= X##_f1)				\
+	  {						\
+	    S##_f1 = T##_f1 + q;			\
+	    X##_f1 -= T##_f1;				\
+	    R##_f1 += q;				\
+	  }						\
+	_FP_FRAC_SLL_2(X, 1);				\
+	q >>= 1;					\
+      }							\
+    q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);		\
+    while (q != _FP_WORK_ROUND)				\
+      {							\
+	T##_f0 = S##_f0 + q;				\
+	T##_f1 = S##_f1;				\
+	if (T##_f1 < X##_f1 || 				\
+	    (T##_f1 == X##_f1 && T##_f0 <= X##_f0))	\
+	  {						\
+	    S##_f0 = T##_f0 + q;			\
+	    S##_f1 += (T##_f0 > S##_f0);		\
+	    _FP_FRAC_DEC_2(X, T);			\
+	    R##_f0 += q;				\
+	  }						\
+	_FP_FRAC_SLL_2(X, 1);				\
+	q >>= 1;					\
+      }							\
+    if (X##_f0 | X##_f1)				\
+      {							\
+	if (S##_f1 < X##_f1 || 				\
+	    (S##_f1 == X##_f1 && S##_f0 < X##_f0))	\
+	  R##_f0 |= _FP_WORK_ROUND;			\
+	R##_f0 |= _FP_WORK_STICKY;			\
+      }							\
+  } while (0)
+
+
+/*
+ * Assembly/disassembly for converting to/from integral types.  
+ * No shifting or overflow handled here.
+ */
+
+#define _FP_FRAC_ASSEMBLE_2(r, X, rsize)	\
+	(void) (((rsize) <= _FP_W_TYPE_SIZE)	\
+		? ({ (r) = X##_f0; })		\
+		: ({				\
+		     (r) = X##_f1;		\
+		     (r) <<= _FP_W_TYPE_SIZE;	\
+		     (r) += X##_f0;		\
+		    }))
+
+#define _FP_FRAC_DISASSEMBLE_2(X, r, rsize)				\
+  do {									\
+    X##_f0 = r;								\
+    X##_f1 = (rsize <= _FP_W_TYPE_SIZE ? 0 : r >> _FP_W_TYPE_SIZE);	\
+  } while (0)
+
+/*
+ * Convert FP values between word sizes
+ */
+
+#define _FP_FRAC_CONV_1_2(dfs, sfs, D, S)				\
+  do {									\
+    if (S##_c != FP_CLS_NAN)						\
+      _FP_FRAC_SRS_2(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs),	\
+		     _FP_WFRACBITS_##sfs);				\
+    else								\
+      _FP_FRAC_SRL_2(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs));	\
+    D##_f = S##_f0;							\
+  } while (0)
+
+#define _FP_FRAC_CONV_2_1(dfs, sfs, D, S)				\
+  do {									\
+    D##_f0 = S##_f;							\
+    D##_f1 = 0;								\
+    _FP_FRAC_SLL_2(D, (_FP_WFRACBITS_##dfs - _FP_WFRACBITS_##sfs));	\
+  } while (0)
+
+#endif
diff --git a/include/math-emu/op-4.h b/include/math-emu/op-4.h
new file mode 100644
index 000000000..ba226f823
--- /dev/null
+++ b/include/math-emu/op-4.h
@@ -0,0 +1,692 @@
+/* Software floating-point emulation.
+   Basic four-word fraction declaration and manipulation.
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef __MATH_EMU_OP_4_H__
+#define __MATH_EMU_OP_4_H__
+
+#define _FP_FRAC_DECL_4(X)	_FP_W_TYPE X##_f[4]
+#define _FP_FRAC_COPY_4(D,S)			\
+  (D##_f[0] = S##_f[0], D##_f[1] = S##_f[1],	\
+   D##_f[2] = S##_f[2], D##_f[3] = S##_f[3])
+#define _FP_FRAC_SET_4(X,I)	__FP_FRAC_SET_4(X, I)
+#define _FP_FRAC_HIGH_4(X)	(X##_f[3])
+#define _FP_FRAC_LOW_4(X)	(X##_f[0])
+#define _FP_FRAC_WORD_4(X,w)	(X##_f[w])
+
+#define _FP_FRAC_SLL_4(X,N)						\
+  do {									\
+    _FP_I_TYPE _up, _down, _skip, _i;					\
+    _skip = (N) / _FP_W_TYPE_SIZE;					\
+    _up = (N) % _FP_W_TYPE_SIZE;					\
+    _down = _FP_W_TYPE_SIZE - _up;					\
+    if (!_up)								\
+      for (_i = 3; _i >= _skip; --_i)					\
+	X##_f[_i] = X##_f[_i-_skip];					\
+    else								\
+      {									\
+	for (_i = 3; _i > _skip; --_i)					\
+	  X##_f[_i] = X##_f[_i-_skip] << _up				\
+		      | X##_f[_i-_skip-1] >> _down;			\
+	X##_f[_i--] = X##_f[0] << _up; 					\
+      }									\
+    for (; _i >= 0; --_i)						\
+      X##_f[_i] = 0;							\
+  } while (0)
+
+/* This one was broken too */
+#define _FP_FRAC_SRL_4(X,N)						\
+  do {									\
+    _FP_I_TYPE _up, _down, _skip, _i;					\
+    _skip = (N) / _FP_W_TYPE_SIZE;					\
+    _down = (N) % _FP_W_TYPE_SIZE;					\
+    _up = _FP_W_TYPE_SIZE - _down;					\
+    if (!_down)								\
+      for (_i = 0; _i <= 3-_skip; ++_i)					\
+	X##_f[_i] = X##_f[_i+_skip];					\
+    else								\
+      {									\
+	for (_i = 0; _i < 3-_skip; ++_i)				\
+	  X##_f[_i] = X##_f[_i+_skip] >> _down				\
+		      | X##_f[_i+_skip+1] << _up;			\
+	X##_f[_i++] = X##_f[3] >> _down;				\
+      }									\
+    for (; _i < 4; ++_i)						\
+      X##_f[_i] = 0;							\
+  } while (0)
+
+
+/* Right shift with sticky-lsb. 
+ * What this actually means is that we do a standard right-shift,
+ * but that if any of the bits that fall off the right hand side
+ * were one then we always set the LSbit.
+ */
+#define _FP_FRAC_SRS_4(X,N,size)					\
+  do {									\
+    _FP_I_TYPE _up, _down, _skip, _i;					\
+    _FP_W_TYPE _s;							\
+    _skip = (N) / _FP_W_TYPE_SIZE;					\
+    _down = (N) % _FP_W_TYPE_SIZE;					\
+    _up = _FP_W_TYPE_SIZE - _down;					\
+    for (_s = _i = 0; _i < _skip; ++_i)					\
+      _s |= X##_f[_i];							\
+    _s |= X##_f[_i] << _up;						\
+/* s is now != 0 if we want to set the LSbit */				\
+    if (!_down)								\
+      for (_i = 0; _i <= 3-_skip; ++_i)					\
+	X##_f[_i] = X##_f[_i+_skip];					\
+    else								\
+      {									\
+	for (_i = 0; _i < 3-_skip; ++_i)				\
+	  X##_f[_i] = X##_f[_i+_skip] >> _down				\
+		      | X##_f[_i+_skip+1] << _up;			\
+	X##_f[_i++] = X##_f[3] >> _down;				\
+      }									\
+    for (; _i < 4; ++_i)						\
+      X##_f[_i] = 0;							\
+    /* don't fix the LSB until the very end when we're sure f[0] is stable */	\
+    X##_f[0] |= (_s != 0);						\
+  } while (0)
+
+#define _FP_FRAC_ADD_4(R,X,Y)						\
+  __FP_FRAC_ADD_4(R##_f[3], R##_f[2], R##_f[1], R##_f[0],		\
+		  X##_f[3], X##_f[2], X##_f[1], X##_f[0],		\
+		  Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
+
+#define _FP_FRAC_SUB_4(R,X,Y)						\
+  __FP_FRAC_SUB_4(R##_f[3], R##_f[2], R##_f[1], R##_f[0],		\
+		  X##_f[3], X##_f[2], X##_f[1], X##_f[0],		\
+		  Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
+
+#define _FP_FRAC_DEC_4(X,Y)						\
+  __FP_FRAC_DEC_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0],		\
+		  Y##_f[3], Y##_f[2], Y##_f[1], Y##_f[0])
+
+#define _FP_FRAC_ADDI_4(X,I)						\
+  __FP_FRAC_ADDI_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0], I)
+
+#define _FP_ZEROFRAC_4  0,0,0,0
+#define _FP_MINFRAC_4   0,0,0,1
+#define _FP_MAXFRAC_4	(~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0), (~(_FP_WS_TYPE)0)
+
+#define _FP_FRAC_ZEROP_4(X)     ((X##_f[0] | X##_f[1] | X##_f[2] | X##_f[3]) == 0)
+#define _FP_FRAC_NEGP_4(X)      ((_FP_WS_TYPE)X##_f[3] < 0)
+#define _FP_FRAC_OVERP_4(fs,X)  (_FP_FRAC_HIGH_##fs(X) & _FP_OVERFLOW_##fs)
+#define _FP_FRAC_CLEAR_OVERP_4(fs,X)  (_FP_FRAC_HIGH_##fs(X) &= ~_FP_OVERFLOW_##fs)
+
+#define _FP_FRAC_EQ_4(X,Y)				\
+ (X##_f[0] == Y##_f[0] && X##_f[1] == Y##_f[1]		\
+  && X##_f[2] == Y##_f[2] && X##_f[3] == Y##_f[3])
+
+#define _FP_FRAC_GT_4(X,Y)				\
+ (X##_f[3] > Y##_f[3] ||				\
+  (X##_f[3] == Y##_f[3] && (X##_f[2] > Y##_f[2] ||	\
+   (X##_f[2] == Y##_f[2] && (X##_f[1] > Y##_f[1] ||	\
+    (X##_f[1] == Y##_f[1] && X##_f[0] > Y##_f[0])	\
+   ))							\
+  ))							\
+ )
+
+#define _FP_FRAC_GE_4(X,Y)				\
+ (X##_f[3] > Y##_f[3] ||				\
+  (X##_f[3] == Y##_f[3] && (X##_f[2] > Y##_f[2] ||	\
+   (X##_f[2] == Y##_f[2] && (X##_f[1] > Y##_f[1] ||	\
+    (X##_f[1] == Y##_f[1] && X##_f[0] >= Y##_f[0])	\
+   ))							\
+  ))							\
+ )
+
+
+#define _FP_FRAC_CLZ_4(R,X)		\
+  do {					\
+    if (X##_f[3])			\
+    {					\
+	__FP_CLZ(R,X##_f[3]);		\
+    }					\
+    else if (X##_f[2])			\
+    {					\
+	__FP_CLZ(R,X##_f[2]);		\
+	R += _FP_W_TYPE_SIZE;		\
+    }					\
+    else if (X##_f[1])			\
+    {					\
+	__FP_CLZ(R,X##_f[2]);		\
+	R += _FP_W_TYPE_SIZE*2;		\
+    }					\
+    else				\
+    {					\
+	__FP_CLZ(R,X##_f[0]);		\
+	R += _FP_W_TYPE_SIZE*3;		\
+    }					\
+  } while(0)
+
+
+#define _FP_UNPACK_RAW_4(fs, X, val)				\
+  do {								\
+    union _FP_UNION_##fs _flo; _flo.flt = (val);		\
+    X##_f[0] = _flo.bits.frac0;					\
+    X##_f[1] = _flo.bits.frac1;					\
+    X##_f[2] = _flo.bits.frac2;					\
+    X##_f[3] = _flo.bits.frac3;					\
+    X##_e  = _flo.bits.exp;					\
+    X##_s  = _flo.bits.sign;					\
+  } while (0)
+
+#define _FP_UNPACK_RAW_4_P(fs, X, val)				\
+  do {								\
+    union _FP_UNION_##fs *_flo =				\
+      (union _FP_UNION_##fs *)(val);				\
+								\
+    X##_f[0] = _flo->bits.frac0;				\
+    X##_f[1] = _flo->bits.frac1;				\
+    X##_f[2] = _flo->bits.frac2;				\
+    X##_f[3] = _flo->bits.frac3;				\
+    X##_e  = _flo->bits.exp;					\
+    X##_s  = _flo->bits.sign;					\
+  } while (0)
+
+#define _FP_PACK_RAW_4(fs, val, X)				\
+  do {								\
+    union _FP_UNION_##fs _flo;					\
+    _flo.bits.frac0 = X##_f[0];					\
+    _flo.bits.frac1 = X##_f[1];					\
+    _flo.bits.frac2 = X##_f[2];					\
+    _flo.bits.frac3 = X##_f[3];					\
+    _flo.bits.exp   = X##_e;					\
+    _flo.bits.sign  = X##_s;					\
+    (val) = _flo.flt;				   		\
+  } while (0)
+
+#define _FP_PACK_RAW_4_P(fs, val, X)				\
+  do {								\
+    union _FP_UNION_##fs *_flo =				\
+      (union _FP_UNION_##fs *)(val);				\
+								\
+    _flo->bits.frac0 = X##_f[0];				\
+    _flo->bits.frac1 = X##_f[1];				\
+    _flo->bits.frac2 = X##_f[2];				\
+    _flo->bits.frac3 = X##_f[3];				\
+    _flo->bits.exp   = X##_e;					\
+    _flo->bits.sign  = X##_s;					\
+  } while (0)
+
+/*
+ * Multiplication algorithms:
+ */
+
+/* Given a 1W * 1W => 2W primitive, do the extended multiplication.  */
+
+#define _FP_MUL_MEAT_4_wide(wfracbits, R, X, Y, doit)			    \
+  do {									    \
+    _FP_FRAC_DECL_8(_z); _FP_FRAC_DECL_2(_b); _FP_FRAC_DECL_2(_c);	    \
+    _FP_FRAC_DECL_2(_d); _FP_FRAC_DECL_2(_e); _FP_FRAC_DECL_2(_f);	    \
+									    \
+    doit(_FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0), X##_f[0], Y##_f[0]); \
+    doit(_b_f1, _b_f0, X##_f[0], Y##_f[1]);				    \
+    doit(_c_f1, _c_f0, X##_f[1], Y##_f[0]);				    \
+    doit(_d_f1, _d_f0, X##_f[1], Y##_f[1]);				    \
+    doit(_e_f1, _e_f0, X##_f[0], Y##_f[2]);				    \
+    doit(_f_f1, _f_f0, X##_f[2], Y##_f[0]);				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),	    \
+		    _FP_FRAC_WORD_8(_z,1), 0,_b_f1,_b_f0,		    \
+		    0,0,_FP_FRAC_WORD_8(_z,1));				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),	    \
+		    _FP_FRAC_WORD_8(_z,1), 0,_c_f1,_c_f0,		    \
+		    _FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2),	    \
+		    _FP_FRAC_WORD_8(_z,1));				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),	    \
+		    _FP_FRAC_WORD_8(_z,2), 0,_d_f1,_d_f0,		    \
+		    0,_FP_FRAC_WORD_8(_z,3),_FP_FRAC_WORD_8(_z,2));	    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),	    \
+		    _FP_FRAC_WORD_8(_z,2), 0,_e_f1,_e_f0,		    \
+		    _FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),	    \
+		    _FP_FRAC_WORD_8(_z,2));				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),	    \
+		    _FP_FRAC_WORD_8(_z,2), 0,_f_f1,_f_f0,		    \
+		    _FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3),	    \
+		    _FP_FRAC_WORD_8(_z,2));				    \
+    doit(_b_f1, _b_f0, X##_f[0], Y##_f[3]);				    \
+    doit(_c_f1, _c_f0, X##_f[3], Y##_f[0]);				    \
+    doit(_d_f1, _d_f0, X##_f[1], Y##_f[2]);				    \
+    doit(_e_f1, _e_f0, X##_f[2], Y##_f[1]);				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),	    \
+		    _FP_FRAC_WORD_8(_z,3), 0,_b_f1,_b_f0,		    \
+		    0,_FP_FRAC_WORD_8(_z,4),_FP_FRAC_WORD_8(_z,3));	    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),	    \
+		    _FP_FRAC_WORD_8(_z,3), 0,_c_f1,_c_f0,		    \
+		    _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),	    \
+		    _FP_FRAC_WORD_8(_z,3));				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),	    \
+		    _FP_FRAC_WORD_8(_z,3), 0,_d_f1,_d_f0,		    \
+		    _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),	    \
+		    _FP_FRAC_WORD_8(_z,3));				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),	    \
+		    _FP_FRAC_WORD_8(_z,3), 0,_e_f1,_e_f0,		    \
+		    _FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4),	    \
+		    _FP_FRAC_WORD_8(_z,3));				    \
+    doit(_b_f1, _b_f0, X##_f[2], Y##_f[2]);				    \
+    doit(_c_f1, _c_f0, X##_f[1], Y##_f[3]);				    \
+    doit(_d_f1, _d_f0, X##_f[3], Y##_f[1]);				    \
+    doit(_e_f1, _e_f0, X##_f[2], Y##_f[3]);				    \
+    doit(_f_f1, _f_f0, X##_f[3], Y##_f[2]);				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),	    \
+		    _FP_FRAC_WORD_8(_z,4), 0,_b_f1,_b_f0,		    \
+		    0,_FP_FRAC_WORD_8(_z,5),_FP_FRAC_WORD_8(_z,4));	    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),	    \
+		    _FP_FRAC_WORD_8(_z,4), 0,_c_f1,_c_f0,		    \
+		    _FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),	    \
+		    _FP_FRAC_WORD_8(_z,4));				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),	    \
+		    _FP_FRAC_WORD_8(_z,4), 0,_d_f1,_d_f0,		    \
+		    _FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5),	    \
+		    _FP_FRAC_WORD_8(_z,4));				    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),	    \
+		    _FP_FRAC_WORD_8(_z,5), 0,_e_f1,_e_f0,		    \
+		    0,_FP_FRAC_WORD_8(_z,6),_FP_FRAC_WORD_8(_z,5));	    \
+    __FP_FRAC_ADD_3(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),	    \
+		    _FP_FRAC_WORD_8(_z,5), 0,_f_f1,_f_f0,		    \
+		    _FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),	    \
+		    _FP_FRAC_WORD_8(_z,5));				    \
+    doit(_b_f1, _b_f0, X##_f[3], Y##_f[3]);				    \
+    __FP_FRAC_ADD_2(_FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6),	    \
+		    _b_f1,_b_f0,					    \
+		    _FP_FRAC_WORD_8(_z,7),_FP_FRAC_WORD_8(_z,6));	    \
+									    \
+    /* Normalize since we know where the msb of the multiplicands	    \
+       were (bit B), we know that the msb of the of the product is	    \
+       at either 2B or 2B-1.  */					    \
+    _FP_FRAC_SRS_8(_z, wfracbits-1, 2*wfracbits);			    \
+    __FP_FRAC_SET_4(R, _FP_FRAC_WORD_8(_z,3), _FP_FRAC_WORD_8(_z,2),	    \
+		    _FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0));	    \
+  } while (0)
+
+#define _FP_MUL_MEAT_4_gmp(wfracbits, R, X, Y)				    \
+  do {									    \
+    _FP_FRAC_DECL_8(_z);						    \
+									    \
+    mpn_mul_n(_z_f, _x_f, _y_f, 4);					    \
+									    \
+    /* Normalize since we know where the msb of the multiplicands	    \
+       were (bit B), we know that the msb of the of the product is	    \
+       at either 2B or 2B-1.  */					    \
+    _FP_FRAC_SRS_8(_z, wfracbits-1, 2*wfracbits);	 		    \
+    __FP_FRAC_SET_4(R, _FP_FRAC_WORD_8(_z,3), _FP_FRAC_WORD_8(_z,2),	    \
+		    _FP_FRAC_WORD_8(_z,1), _FP_FRAC_WORD_8(_z,0));	    \
+  } while (0)
+
+/*
+ * Helper utility for _FP_DIV_MEAT_4_udiv:
+ * pppp = m * nnn
+ */
+#define umul_ppppmnnn(p3,p2,p1,p0,m,n2,n1,n0)				    \
+  do {									    \
+    UWtype _t;								    \
+    umul_ppmm(p1,p0,m,n0);						    \
+    umul_ppmm(p2,_t,m,n1);						    \
+    __FP_FRAC_ADDI_2(p2,p1,_t);						    \
+    umul_ppmm(p3,_t,m,n2);						    \
+    __FP_FRAC_ADDI_2(p3,p2,_t);						    \
+  } while (0)
+
+/*
+ * Division algorithms:
+ */
+
+#define _FP_DIV_MEAT_4_udiv(fs, R, X, Y)				    \
+  do {									    \
+    int _i;								    \
+    _FP_FRAC_DECL_4(_n); _FP_FRAC_DECL_4(_m);				    \
+    _FP_FRAC_SET_4(_n, _FP_ZEROFRAC_4);					    \
+    if (_FP_FRAC_GT_4(X, Y))						    \
+      {									    \
+	_n_f[3] = X##_f[0] << (_FP_W_TYPE_SIZE - 1);			    \
+	_FP_FRAC_SRL_4(X, 1);						    \
+      }									    \
+    else								    \
+      R##_e--;								    \
+									    \
+    /* Normalize, i.e. make the most significant bit of the 		    \
+       denominator set. */						    \
+    _FP_FRAC_SLL_4(Y, _FP_WFRACXBITS_##fs);				    \
+									    \
+    for (_i = 3; ; _i--)						    \
+      {									    \
+        if (X##_f[3] == Y##_f[3])					    \
+          {								    \
+            /* This is a special case, not an optimization		    \
+               (X##_f[3]/Y##_f[3] would not fit into UWtype).		    \
+               As X## is guaranteed to be < Y,  R##_f[_i] can be either	    \
+               (UWtype)-1 or (UWtype)-2.  */				    \
+            R##_f[_i] = -1;						    \
+            if (!_i)							    \
+	      break;							    \
+            __FP_FRAC_SUB_4(X##_f[3], X##_f[2], X##_f[1], X##_f[0],	    \
+			    Y##_f[2], Y##_f[1], Y##_f[0], 0,		    \
+			    X##_f[2], X##_f[1], X##_f[0], _n_f[_i]);	    \
+            _FP_FRAC_SUB_4(X, Y, X);					    \
+            if (X##_f[3] > Y##_f[3])					    \
+              {								    \
+                R##_f[_i] = -2;						    \
+                _FP_FRAC_ADD_4(X, Y, X);				    \
+              }								    \
+          }								    \
+        else								    \
+          {								    \
+            udiv_qrnnd(R##_f[_i], X##_f[3], X##_f[3], X##_f[2], Y##_f[3]);  \
+            umul_ppppmnnn(_m_f[3], _m_f[2], _m_f[1], _m_f[0],		    \
+			  R##_f[_i], Y##_f[2], Y##_f[1], Y##_f[0]);	    \
+            X##_f[2] = X##_f[1];					    \
+            X##_f[1] = X##_f[0];					    \
+            X##_f[0] = _n_f[_i];					    \
+            if (_FP_FRAC_GT_4(_m, X))					    \
+              {								    \
+                R##_f[_i]--;						    \
+                _FP_FRAC_ADD_4(X, Y, X);				    \
+                if (_FP_FRAC_GE_4(X, Y) && _FP_FRAC_GT_4(_m, X))	    \
+                  {							    \
+		    R##_f[_i]--;					    \
+		    _FP_FRAC_ADD_4(X, Y, X);				    \
+                  }							    \
+              }								    \
+            _FP_FRAC_DEC_4(X, _m);					    \
+            if (!_i)							    \
+	      {								    \
+		if (!_FP_FRAC_EQ_4(X, _m))				    \
+		  R##_f[0] |= _FP_WORK_STICKY;				    \
+		break;							    \
+	      }								    \
+          }								    \
+      }									    \
+  } while (0)
+
+
+/*
+ * Square root algorithms:
+ * We have just one right now, maybe Newton approximation
+ * should be added for those machines where division is fast.
+ */
+ 
+#define _FP_SQRT_MEAT_4(R, S, T, X, q)				\
+  do {								\
+    while (q)							\
+      {								\
+	T##_f[3] = S##_f[3] + q;				\
+	if (T##_f[3] <= X##_f[3])				\
+	  {							\
+	    S##_f[3] = T##_f[3] + q;				\
+	    X##_f[3] -= T##_f[3];				\
+	    R##_f[3] += q;					\
+	  }							\
+	_FP_FRAC_SLL_4(X, 1);					\
+	q >>= 1;						\
+      }								\
+    q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);			\
+    while (q)							\
+      {								\
+	T##_f[2] = S##_f[2] + q;				\
+	T##_f[3] = S##_f[3];					\
+	if (T##_f[3] < X##_f[3] || 				\
+	    (T##_f[3] == X##_f[3] && T##_f[2] <= X##_f[2]))	\
+	  {							\
+	    S##_f[2] = T##_f[2] + q;				\
+	    S##_f[3] += (T##_f[2] > S##_f[2]);			\
+	    __FP_FRAC_DEC_2(X##_f[3], X##_f[2],			\
+			    T##_f[3], T##_f[2]);		\
+	    R##_f[2] += q;					\
+	  }							\
+	_FP_FRAC_SLL_4(X, 1);					\
+	q >>= 1;						\
+      }								\
+    q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);			\
+    while (q)							\
+      {								\
+	T##_f[1] = S##_f[1] + q;				\
+	T##_f[2] = S##_f[2];					\
+	T##_f[3] = S##_f[3];					\
+	if (T##_f[3] < X##_f[3] || 				\
+	    (T##_f[3] == X##_f[3] && (T##_f[2] < X##_f[2] ||	\
+	     (T##_f[2] == X##_f[2] && T##_f[1] <= X##_f[1]))))	\
+	  {							\
+	    S##_f[1] = T##_f[1] + q;				\
+	    S##_f[2] += (T##_f[1] > S##_f[1]);			\
+	    S##_f[3] += (T##_f[2] > S##_f[2]);			\
+	    __FP_FRAC_DEC_3(X##_f[3], X##_f[2], X##_f[1],	\
+	    		    T##_f[3], T##_f[2], T##_f[1]);	\
+	    R##_f[1] += q;					\
+	  }							\
+	_FP_FRAC_SLL_4(X, 1);					\
+	q >>= 1;						\
+      }								\
+    q = (_FP_W_TYPE)1 << (_FP_W_TYPE_SIZE - 1);			\
+    while (q != _FP_WORK_ROUND)					\
+      {								\
+	T##_f[0] = S##_f[0] + q;				\
+	T##_f[1] = S##_f[1];					\
+	T##_f[2] = S##_f[2];					\
+	T##_f[3] = S##_f[3];					\
+	if (_FP_FRAC_GE_4(X,T))					\
+	  {							\
+	    S##_f[0] = T##_f[0] + q;				\
+	    S##_f[1] += (T##_f[0] > S##_f[0]);			\
+	    S##_f[2] += (T##_f[1] > S##_f[1]);			\
+	    S##_f[3] += (T##_f[2] > S##_f[2]);			\
+	    _FP_FRAC_DEC_4(X, T);				\
+	    R##_f[0] += q;					\
+	  }							\
+	_FP_FRAC_SLL_4(X, 1);					\
+	q >>= 1;						\
+      }								\
+    if (!_FP_FRAC_ZEROP_4(X))					\
+      {								\
+	if (_FP_FRAC_GT_4(X,S))					\
+	  R##_f[0] |= _FP_WORK_ROUND;				\
+	R##_f[0] |= _FP_WORK_STICKY;				\
+      }								\
+  } while (0)
+
+
+/*
+ * Internals 
+ */
+
+#define __FP_FRAC_SET_4(X,I3,I2,I1,I0)					\
+  (X##_f[3] = I3, X##_f[2] = I2, X##_f[1] = I1, X##_f[0] = I0)
+
+#ifndef __FP_FRAC_ADD_3
+#define __FP_FRAC_ADD_3(r2,r1,r0,x2,x1,x0,y2,y1,y0)		\
+  do {								\
+    int _c1, _c2;							\
+    r0 = x0 + y0;						\
+    _c1 = r0 < x0;						\
+    r1 = x1 + y1;						\
+    _c2 = r1 < x1;						\
+    r1 += _c1;							\
+    _c2 |= r1 < _c1;						\
+    r2 = x2 + y2 + _c2;						\
+  } while (0)
+#endif
+
+#ifndef __FP_FRAC_ADD_4
+#define __FP_FRAC_ADD_4(r3,r2,r1,r0,x3,x2,x1,x0,y3,y2,y1,y0)	\
+  do {								\
+    int _c1, _c2, _c3;						\
+    r0 = x0 + y0;						\
+    _c1 = r0 < x0;						\
+    r1 = x1 + y1;						\
+    _c2 = r1 < x1;						\
+    r1 += _c1;							\
+    _c2 |= r1 < _c1;						\
+    r2 = x2 + y2;						\
+    _c3 = r2 < x2;						\
+    r2 += _c2;							\
+    _c3 |= r2 < _c2;						\
+    r3 = x3 + y3 + _c3;						\
+  } while (0)
+#endif
+
+#ifndef __FP_FRAC_SUB_3
+#define __FP_FRAC_SUB_3(r2,r1,r0,x2,x1,x0,y2,y1,y0)		\
+  do {								\
+    int _c1, _c2;							\
+    r0 = x0 - y0;						\
+    _c1 = r0 > x0;						\
+    r1 = x1 - y1;						\
+    _c2 = r1 > x1;						\
+    r1 -= _c1;							\
+    _c2 |= r1 > _c1;						\
+    r2 = x2 - y2 - _c2;						\
+  } while (0)
+#endif
+
+#ifndef __FP_FRAC_SUB_4
+#define __FP_FRAC_SUB_4(r3,r2,r1,r0,x3,x2,x1,x0,y3,y2,y1,y0)	\
+  do {								\
+    int _c1, _c2, _c3;						\
+    r0 = x0 - y0;						\
+    _c1 = r0 > x0;						\
+    r1 = x1 - y1;						\
+    _c2 = r1 > x1;						\
+    r1 -= _c1;							\
+    _c2 |= r1 > _c1;						\
+    r2 = x2 - y2;						\
+    _c3 = r2 > x2;						\
+    r2 -= _c2;							\
+    _c3 |= r2 > _c2;						\
+    r3 = x3 - y3 - _c3;						\
+  } while (0)
+#endif
+
+#ifndef __FP_FRAC_DEC_3
+#define __FP_FRAC_DEC_3(x2,x1,x0,y2,y1,y0)				\
+  do {									\
+    UWtype _t0, _t1, _t2;						\
+    _t0 = x0, _t1 = x1, _t2 = x2;					\
+    __FP_FRAC_SUB_3 (x2, x1, x0, _t2, _t1, _t0, y2, y1, y0);		\
+  } while (0)
+#endif
+
+#ifndef __FP_FRAC_DEC_4
+#define __FP_FRAC_DEC_4(x3,x2,x1,x0,y3,y2,y1,y0)			\
+  do {									\
+    UWtype _t0, _t1, _t2, _t3;						\
+    _t0 = x0, _t1 = x1, _t2 = x2, _t3 = x3;				\
+    __FP_FRAC_SUB_4 (x3,x2,x1,x0,_t3,_t2,_t1,_t0, y3,y2,y1,y0);		\
+  } while (0)
+#endif
+
+#ifndef __FP_FRAC_ADDI_4
+#define __FP_FRAC_ADDI_4(x3,x2,x1,x0,i)					\
+  do {									\
+    UWtype _t;								\
+    _t = ((x0 += i) < i);						\
+    x1 += _t; _t = (x1 < _t);						\
+    x2 += _t; _t = (x2 < _t);						\
+    x3 += _t;								\
+  } while (0)
+#endif
+
+/* Convert FP values between word sizes. This appears to be more
+ * complicated than I'd have expected it to be, so these might be
+ * wrong... These macros are in any case somewhat bogus because they
+ * use information about what various FRAC_n variables look like 
+ * internally [eg, that 2 word vars are X_f0 and x_f1]. But so do
+ * the ones in op-2.h and op-1.h. 
+ */
+#define _FP_FRAC_CONV_1_4(dfs, sfs, D, S)				\
+   do {									\
+     if (S##_c != FP_CLS_NAN)						\
+       _FP_FRAC_SRS_4(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs),	\
+			  _FP_WFRACBITS_##sfs);				\
+     else								\
+       _FP_FRAC_SRL_4(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs));	\
+     D##_f = S##_f[0];							\
+  } while (0)
+
+#define _FP_FRAC_CONV_2_4(dfs, sfs, D, S)				\
+   do {									\
+     if (S##_c != FP_CLS_NAN)						\
+       _FP_FRAC_SRS_4(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs),	\
+		      _FP_WFRACBITS_##sfs);				\
+     else								\
+       _FP_FRAC_SRL_4(S, (_FP_WFRACBITS_##sfs - _FP_WFRACBITS_##dfs));	\
+     D##_f0 = S##_f[0];							\
+     D##_f1 = S##_f[1];							\
+  } while (0)
+
+/* Assembly/disassembly for converting to/from integral types.  
+ * No shifting or overflow handled here.
+ */
+/* Put the FP value X into r, which is an integer of size rsize. */
+#define _FP_FRAC_ASSEMBLE_4(r, X, rsize)				\
+  do {									\
+    if (rsize <= _FP_W_TYPE_SIZE)					\
+      r = X##_f[0];							\
+    else if (rsize <= 2*_FP_W_TYPE_SIZE)				\
+    {									\
+      r = X##_f[1];							\
+      r <<= _FP_W_TYPE_SIZE;						\
+      r += X##_f[0];							\
+    }									\
+    else								\
+    {									\
+      /* I'm feeling lazy so we deal with int == 3words (implausible)*/	\
+      /* and int == 4words as a single case.			 */	\
+      r = X##_f[3];							\
+      r <<= _FP_W_TYPE_SIZE;						\
+      r += X##_f[2];							\
+      r <<= _FP_W_TYPE_SIZE;						\
+      r += X##_f[1];							\
+      r <<= _FP_W_TYPE_SIZE;						\
+      r += X##_f[0];							\
+    }									\
+  } while (0)
+
+/* "No disassemble Number Five!" */
+/* move an integer of size rsize into X's fractional part. We rely on
+ * the _f[] array consisting of words of size _FP_W_TYPE_SIZE to avoid
+ * having to mask the values we store into it.
+ */
+#define _FP_FRAC_DISASSEMBLE_4(X, r, rsize)				\
+  do {									\
+    X##_f[0] = r;							\
+    X##_f[1] = (rsize <= _FP_W_TYPE_SIZE ? 0 : r >> _FP_W_TYPE_SIZE);	\
+    X##_f[2] = (rsize <= 2*_FP_W_TYPE_SIZE ? 0 : r >> 2*_FP_W_TYPE_SIZE); \
+    X##_f[3] = (rsize <= 3*_FP_W_TYPE_SIZE ? 0 : r >> 3*_FP_W_TYPE_SIZE); \
+  } while (0)
+
+#define _FP_FRAC_CONV_4_1(dfs, sfs, D, S)				\
+   do {									\
+     D##_f[0] = S##_f;							\
+     D##_f[1] = D##_f[2] = D##_f[3] = 0;				\
+     _FP_FRAC_SLL_4(D, (_FP_WFRACBITS_##dfs - _FP_WFRACBITS_##sfs));	\
+   } while (0)
+
+#define _FP_FRAC_CONV_4_2(dfs, sfs, D, S)				\
+   do {									\
+     D##_f[0] = S##_f0;							\
+     D##_f[1] = S##_f1;							\
+     D##_f[2] = D##_f[3] = 0;						\
+     _FP_FRAC_SLL_4(D, (_FP_WFRACBITS_##dfs - _FP_WFRACBITS_##sfs));	\
+   } while (0)
+
+#endif
diff --git a/include/math-emu/op-8.h b/include/math-emu/op-8.h
new file mode 100644
index 000000000..8b8c05eb8
--- /dev/null
+++ b/include/math-emu/op-8.h
@@ -0,0 +1,107 @@
+/* Software floating-point emulation.
+   Basic eight-word fraction declaration and manipulation.
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+                                                         
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef __MATH_EMU_OP_8_H__
+#define __MATH_EMU_OP_8_H__
+
+/* We need just a few things from here for op-4, if we ever need some
+   other macros, they can be added. */
+#define _FP_FRAC_DECL_8(X)	_FP_W_TYPE X##_f[8]
+#define _FP_FRAC_HIGH_8(X)	(X##_f[7])
+#define _FP_FRAC_LOW_8(X)	(X##_f[0])
+#define _FP_FRAC_WORD_8(X,w)	(X##_f[w])
+
+#define _FP_FRAC_SLL_8(X,N)						\
+  do {									\
+    _FP_I_TYPE _up, _down, _skip, _i;					\
+    _skip = (N) / _FP_W_TYPE_SIZE;					\
+    _up = (N) % _FP_W_TYPE_SIZE;					\
+    _down = _FP_W_TYPE_SIZE - _up;					\
+    if (!_up)								\
+      for (_i = 7; _i >= _skip; --_i)					\
+	X##_f[_i] = X##_f[_i-_skip];					\
+    else								\
+      {									\
+	for (_i = 7; _i > _skip; --_i)					\
+	  X##_f[_i] = X##_f[_i-_skip] << _up				\
+		      | X##_f[_i-_skip-1] >> _down;			\
+	X##_f[_i--] = X##_f[0] << _up; 					\
+      }									\
+    for (; _i >= 0; --_i)						\
+      X##_f[_i] = 0;							\
+  } while (0)
+
+#define _FP_FRAC_SRL_8(X,N)						\
+  do {									\
+    _FP_I_TYPE _up, _down, _skip, _i;					\
+    _skip = (N) / _FP_W_TYPE_SIZE;					\
+    _down = (N) % _FP_W_TYPE_SIZE;					\
+    _up = _FP_W_TYPE_SIZE - _down;					\
+    if (!_down)								\
+      for (_i = 0; _i <= 7-_skip; ++_i)					\
+	X##_f[_i] = X##_f[_i+_skip];					\
+    else								\
+      {									\
+	for (_i = 0; _i < 7-_skip; ++_i)				\
+	  X##_f[_i] = X##_f[_i+_skip] >> _down				\
+		      | X##_f[_i+_skip+1] << _up;			\
+	X##_f[_i++] = X##_f[7] >> _down;				\
+      }									\
+    for (; _i < 8; ++_i)						\
+      X##_f[_i] = 0;							\
+  } while (0)
+
+
+/* Right shift with sticky-lsb. 
+ * What this actually means is that we do a standard right-shift,
+ * but that if any of the bits that fall off the right hand side
+ * were one then we always set the LSbit.
+ */
+#define _FP_FRAC_SRS_8(X,N,size)					\
+  do {									\
+    _FP_I_TYPE _up, _down, _skip, _i;					\
+    _FP_W_TYPE _s;							\
+    _skip = (N) / _FP_W_TYPE_SIZE;					\
+    _down = (N) % _FP_W_TYPE_SIZE;					\
+    _up = _FP_W_TYPE_SIZE - _down;					\
+    for (_s = _i = 0; _i < _skip; ++_i)					\
+      _s |= X##_f[_i];							\
+    _s |= X##_f[_i] << _up;						\
+/* s is now != 0 if we want to set the LSbit */				\
+    if (!_down)								\
+      for (_i = 0; _i <= 7-_skip; ++_i)					\
+	X##_f[_i] = X##_f[_i+_skip];					\
+    else								\
+      {									\
+	for (_i = 0; _i < 7-_skip; ++_i)				\
+	  X##_f[_i] = X##_f[_i+_skip] >> _down				\
+		      | X##_f[_i+_skip+1] << _up;			\
+	X##_f[_i++] = X##_f[7] >> _down;				\
+      }									\
+    for (; _i < 8; ++_i)						\
+      X##_f[_i] = 0;							\
+    /* don't fix the LSB until the very end when we're sure f[0] is stable */	\
+    X##_f[0] |= (_s != 0);						\
+  } while (0)
+
+#endif
diff --git a/include/math-emu/op-common.h b/include/math-emu/op-common.h
new file mode 100644
index 000000000..8ce066c03
--- /dev/null
+++ b/include/math-emu/op-common.h
@@ -0,0 +1,885 @@
+/* Software floating-point emulation. Common operations.
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef __MATH_EMU_OP_COMMON_H__
+#define __MATH_EMU_OP_COMMON_H__
+
+#define _FP_DECL(wc, X)			\
+  _FP_I_TYPE X##_c=0, X##_s=0, X##_e=0;	\
+  _FP_FRAC_DECL_##wc(X)
+
+/*
+ * Finish truly unpacking a native fp value by classifying the kind
+ * of fp value and normalizing both the exponent and the fraction.
+ */
+
+#define _FP_UNPACK_CANONICAL(fs, wc, X)					\
+do {									\
+  switch (X##_e)							\
+  {									\
+  default:								\
+    _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_IMPLBIT_##fs;			\
+    _FP_FRAC_SLL_##wc(X, _FP_WORKBITS);					\
+    X##_e -= _FP_EXPBIAS_##fs;						\
+    X##_c = FP_CLS_NORMAL;						\
+    break;								\
+									\
+  case 0:								\
+    if (_FP_FRAC_ZEROP_##wc(X))						\
+      X##_c = FP_CLS_ZERO;						\
+    else								\
+      {									\
+	/* a denormalized number */					\
+	_FP_I_TYPE _shift;						\
+	_FP_FRAC_CLZ_##wc(_shift, X);					\
+	_shift -= _FP_FRACXBITS_##fs;					\
+	_FP_FRAC_SLL_##wc(X, (_shift+_FP_WORKBITS));			\
+	X##_e -= _FP_EXPBIAS_##fs - 1 + _shift;				\
+	X##_c = FP_CLS_NORMAL;						\
+	FP_SET_EXCEPTION(FP_EX_DENORM);					\
+	if (FP_DENORM_ZERO)						\
+	  {								\
+	    FP_SET_EXCEPTION(FP_EX_INEXACT);				\
+	    X##_c = FP_CLS_ZERO;					\
+	  }								\
+      }									\
+    break;								\
+									\
+  case _FP_EXPMAX_##fs:							\
+    if (_FP_FRAC_ZEROP_##wc(X))						\
+      X##_c = FP_CLS_INF;						\
+    else								\
+      {									\
+	X##_c = FP_CLS_NAN;						\
+	/* Check for signaling NaN */					\
+	if (!(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs))		\
+	  FP_SET_EXCEPTION(FP_EX_INVALID | FP_EX_INVALID_SNAN);		\
+      }									\
+    break;								\
+  }									\
+} while (0)
+
+/*
+ * Before packing the bits back into the native fp result, take care
+ * of such mundane things as rounding and overflow.  Also, for some
+ * kinds of fp values, the original parts may not have been fully
+ * extracted -- but that is ok, we can regenerate them now.
+ */
+
+#define _FP_PACK_CANONICAL(fs, wc, X)				\
+do {								\
+  switch (X##_c)						\
+  {								\
+  case FP_CLS_NORMAL:						\
+    X##_e += _FP_EXPBIAS_##fs;					\
+    if (X##_e > 0)						\
+      {								\
+	_FP_ROUND(wc, X);					\
+	if (_FP_FRAC_OVERP_##wc(fs, X))				\
+	  {							\
+	    _FP_FRAC_CLEAR_OVERP_##wc(fs, X);			\
+	    X##_e++;						\
+	  }							\
+	_FP_FRAC_SRL_##wc(X, _FP_WORKBITS);			\
+	if (X##_e >= _FP_EXPMAX_##fs)				\
+	  {							\
+	    /* overflow */					\
+	    switch (FP_ROUNDMODE)				\
+	      {							\
+	      case FP_RND_NEAREST:				\
+		X##_c = FP_CLS_INF;				\
+		break;						\
+	      case FP_RND_PINF:					\
+		if (!X##_s) X##_c = FP_CLS_INF;			\
+		break;						\
+	      case FP_RND_MINF:					\
+		if (X##_s) X##_c = FP_CLS_INF;			\
+		break;						\
+	      }							\
+	    if (X##_c == FP_CLS_INF)				\
+	      {							\
+		/* Overflow to infinity */			\
+		X##_e = _FP_EXPMAX_##fs;			\
+		_FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);	\
+	      }							\
+	    else						\
+	      {							\
+		/* Overflow to maximum normal */		\
+		X##_e = _FP_EXPMAX_##fs - 1;			\
+		_FP_FRAC_SET_##wc(X, _FP_MAXFRAC_##wc);		\
+	      }							\
+	    FP_SET_EXCEPTION(FP_EX_OVERFLOW);			\
+            FP_SET_EXCEPTION(FP_EX_INEXACT);			\
+	  }							\
+      }								\
+    else							\
+      {								\
+	/* we've got a denormalized number */			\
+	X##_e = -X##_e + 1;					\
+	if (X##_e <= _FP_WFRACBITS_##fs)			\
+	  {							\
+	    _FP_FRAC_SRS_##wc(X, X##_e, _FP_WFRACBITS_##fs);	\
+	    if (_FP_FRAC_HIGH_##fs(X)				\
+		& (_FP_OVERFLOW_##fs >> 1))			\
+	      {							\
+	        X##_e = 1;					\
+	        _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);	\
+	      }							\
+	    else						\
+	      {							\
+		_FP_ROUND(wc, X);				\
+		if (_FP_FRAC_HIGH_##fs(X)			\
+		   & (_FP_OVERFLOW_##fs >> 1))			\
+		  {						\
+		    X##_e = 1;					\
+		    _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);	\
+		    FP_SET_EXCEPTION(FP_EX_INEXACT);		\
+		  }						\
+		else						\
+		  {						\
+		    X##_e = 0;					\
+		    _FP_FRAC_SRL_##wc(X, _FP_WORKBITS);		\
+		  }						\
+	      }							\
+	    if ((FP_CUR_EXCEPTIONS & FP_EX_INEXACT) ||		\
+		(FP_TRAPPING_EXCEPTIONS & FP_EX_UNDERFLOW))	\
+		FP_SET_EXCEPTION(FP_EX_UNDERFLOW);		\
+	  }							\
+	else							\
+	  {							\
+	    /* underflow to zero */				\
+	    X##_e = 0;						\
+	    if (!_FP_FRAC_ZEROP_##wc(X))			\
+	      {							\
+	        _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc);		\
+	        _FP_ROUND(wc, X);				\
+	        _FP_FRAC_LOW_##wc(X) >>= (_FP_WORKBITS);	\
+	      }							\
+	    FP_SET_EXCEPTION(FP_EX_UNDERFLOW);			\
+	  }							\
+      }								\
+    break;							\
+								\
+  case FP_CLS_ZERO:						\
+    X##_e = 0;							\
+    _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);			\
+    break;							\
+								\
+  case FP_CLS_INF:						\
+    X##_e = _FP_EXPMAX_##fs;					\
+    _FP_FRAC_SET_##wc(X, _FP_ZEROFRAC_##wc);			\
+    break;							\
+								\
+  case FP_CLS_NAN:						\
+    X##_e = _FP_EXPMAX_##fs;					\
+    if (!_FP_KEEPNANFRACP)					\
+      {								\
+	_FP_FRAC_SET_##wc(X, _FP_NANFRAC_##fs);			\
+	X##_s = _FP_NANSIGN_##fs;				\
+      }								\
+    else							\
+      _FP_FRAC_HIGH_RAW_##fs(X) |= _FP_QNANBIT_##fs;		\
+    break;							\
+  }								\
+} while (0)
+
+/* This one accepts raw argument and not cooked,  returns
+ * 1 if X is a signaling NaN.
+ */
+#define _FP_ISSIGNAN(fs, wc, X)					\
+({								\
+  int __ret = 0;						\
+  if (X##_e == _FP_EXPMAX_##fs)					\
+    {								\
+      if (!_FP_FRAC_ZEROP_##wc(X)				\
+	  && !(_FP_FRAC_HIGH_RAW_##fs(X) & _FP_QNANBIT_##fs))	\
+	__ret = 1;						\
+    }								\
+  __ret;							\
+})
+
+
+
+
+
+/*
+ * Main addition routine.  The input values should be cooked.
+ */
+
+#define _FP_ADD_INTERNAL(fs, wc, R, X, Y, OP)				     \
+do {									     \
+  switch (_FP_CLS_COMBINE(X##_c, Y##_c))				     \
+  {									     \
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL):			     \
+    {									     \
+      /* shift the smaller number so that its exponent matches the larger */ \
+      _FP_I_TYPE diff = X##_e - Y##_e;					     \
+									     \
+      if (diff < 0)							     \
+	{								     \
+	  diff = -diff;							     \
+	  if (diff <= _FP_WFRACBITS_##fs)				     \
+	    _FP_FRAC_SRS_##wc(X, diff, _FP_WFRACBITS_##fs);		     \
+	  else if (!_FP_FRAC_ZEROP_##wc(X))				     \
+	    _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc);			     \
+	  R##_e = Y##_e;						     \
+	}								     \
+      else								     \
+	{								     \
+	  if (diff > 0)							     \
+	    {								     \
+	      if (diff <= _FP_WFRACBITS_##fs)				     \
+	        _FP_FRAC_SRS_##wc(Y, diff, _FP_WFRACBITS_##fs);		     \
+	      else if (!_FP_FRAC_ZEROP_##wc(Y))				     \
+	        _FP_FRAC_SET_##wc(Y, _FP_MINFRAC_##wc);			     \
+	    }								     \
+	  R##_e = X##_e;						     \
+	}								     \
+									     \
+      R##_c = FP_CLS_NORMAL;						     \
+									     \
+      if (X##_s == Y##_s)						     \
+	{								     \
+	  R##_s = X##_s;						     \
+	  _FP_FRAC_ADD_##wc(R, X, Y);					     \
+	  if (_FP_FRAC_OVERP_##wc(fs, R))				     \
+	    {								     \
+	      _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs);		     \
+	      R##_e++;							     \
+	    }								     \
+	}								     \
+      else								     \
+	{								     \
+	  R##_s = X##_s;						     \
+	  _FP_FRAC_SUB_##wc(R, X, Y);					     \
+	  if (_FP_FRAC_ZEROP_##wc(R))					     \
+	    {								     \
+	      /* return an exact zero */				     \
+	      if (FP_ROUNDMODE == FP_RND_MINF)				     \
+		R##_s |= Y##_s;						     \
+	      else							     \
+		R##_s &= Y##_s;						     \
+	      R##_c = FP_CLS_ZERO;					     \
+	    }								     \
+	  else								     \
+	    {								     \
+	      if (_FP_FRAC_NEGP_##wc(R))				     \
+		{							     \
+		  _FP_FRAC_SUB_##wc(R, Y, X);				     \
+		  R##_s = Y##_s;					     \
+		}							     \
+									     \
+	      /* renormalize after subtraction */			     \
+	      _FP_FRAC_CLZ_##wc(diff, R);				     \
+	      diff -= _FP_WFRACXBITS_##fs;				     \
+	      if (diff)							     \
+		{							     \
+		  R##_e -= diff;					     \
+		  _FP_FRAC_SLL_##wc(R, diff);				     \
+		}							     \
+	    }								     \
+	}								     \
+      break;								     \
+    }									     \
+									     \
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN):				     \
+    _FP_CHOOSENAN(fs, wc, R, X, Y, OP);					     \
+    break;								     \
+									     \
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO):			     \
+    R##_e = X##_e;							     \
+	fallthrough;							     \
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL):			     \
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF):				     \
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO):				     \
+    _FP_FRAC_COPY_##wc(R, X);						     \
+    R##_s = X##_s;							     \
+    R##_c = X##_c;							     \
+    break;								     \
+									     \
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL):			     \
+    R##_e = Y##_e;							     \
+	fallthrough;							     \
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN):			     \
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN):				     \
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN):				     \
+    _FP_FRAC_COPY_##wc(R, Y);						     \
+    R##_s = Y##_s;							     \
+    R##_c = Y##_c;							     \
+    break;								     \
+									     \
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF):				     \
+    if (X##_s != Y##_s)							     \
+      {									     \
+	/* +INF + -INF => NAN */					     \
+	_FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);				     \
+	R##_s = _FP_NANSIGN_##fs;					     \
+	R##_c = FP_CLS_NAN;						     \
+	FP_SET_EXCEPTION(FP_EX_INVALID | FP_EX_INVALID_ISI);		     \
+	break;								     \
+      }									     \
+    fallthrough;							     \
+									     \
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL):			     \
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO):				     \
+    R##_s = X##_s;							     \
+    R##_c = FP_CLS_INF;							     \
+    break;								     \
+									     \
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF):			     \
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF):				     \
+    R##_s = Y##_s;							     \
+    R##_c = FP_CLS_INF;							     \
+    break;								     \
+									     \
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO):			     \
+    /* make sure the sign is correct */					     \
+    if (FP_ROUNDMODE == FP_RND_MINF)					     \
+      R##_s = X##_s | Y##_s;						     \
+    else								     \
+      R##_s = X##_s & Y##_s;						     \
+    R##_c = FP_CLS_ZERO;						     \
+    break;								     \
+									     \
+  default:								     \
+    abort();								     \
+  }									     \
+} while (0)
+
+#define _FP_ADD(fs, wc, R, X, Y) _FP_ADD_INTERNAL(fs, wc, R, X, Y, '+')
+#define _FP_SUB(fs, wc, R, X, Y)					     \
+  do {									     \
+    if (Y##_c != FP_CLS_NAN) Y##_s ^= 1;				     \
+    _FP_ADD_INTERNAL(fs, wc, R, X, Y, '-');				     \
+  } while (0)
+
+
+/*
+ * Main negation routine.  FIXME -- when we care about setting exception
+ * bits reliably, this will not do.  We should examine all of the fp classes.
+ */
+
+#define _FP_NEG(fs, wc, R, X)		\
+  do {					\
+    _FP_FRAC_COPY_##wc(R, X);		\
+    R##_c = X##_c;			\
+    R##_e = X##_e;			\
+    R##_s = 1 ^ X##_s;			\
+  } while (0)
+
+
+/*
+ * Main multiplication routine.  The input values should be cooked.
+ */
+
+#define _FP_MUL(fs, wc, R, X, Y)			\
+do {							\
+  R##_s = X##_s ^ Y##_s;				\
+  switch (_FP_CLS_COMBINE(X##_c, Y##_c))		\
+  {							\
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL):	\
+    R##_c = FP_CLS_NORMAL;				\
+    R##_e = X##_e + Y##_e + 1;				\
+							\
+    _FP_MUL_MEAT_##fs(R,X,Y);				\
+							\
+    if (_FP_FRAC_OVERP_##wc(fs, R))			\
+      _FP_FRAC_SRS_##wc(R, 1, _FP_WFRACBITS_##fs);	\
+    else						\
+      R##_e--;						\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN):		\
+    _FP_CHOOSENAN(fs, wc, R, X, Y, '*');		\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL):	\
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF):		\
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO):		\
+    R##_s = X##_s;					\
+	  fallthrough;					\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF):		\
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL):	\
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL):	\
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO):	\
+    _FP_FRAC_COPY_##wc(R, X);				\
+    R##_c = X##_c;					\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN):	\
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN):		\
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN):		\
+    R##_s = Y##_s;					\
+	  fallthrough;					\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF):	\
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO):	\
+    _FP_FRAC_COPY_##wc(R, Y);				\
+    R##_c = Y##_c;					\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO):		\
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF):		\
+    R##_s = _FP_NANSIGN_##fs;				\
+    R##_c = FP_CLS_NAN;					\
+    _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);		\
+    FP_SET_EXCEPTION(FP_EX_INVALID | FP_EX_INVALID_IMZ);\
+    break;						\
+							\
+  default:						\
+    abort();						\
+  }							\
+} while (0)
+
+
+/*
+ * Main division routine.  The input values should be cooked.
+ */
+
+#define _FP_DIV(fs, wc, R, X, Y)			\
+do {							\
+  R##_s = X##_s ^ Y##_s;				\
+  switch (_FP_CLS_COMBINE(X##_c, Y##_c))		\
+  {							\
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NORMAL):	\
+    R##_c = FP_CLS_NORMAL;				\
+    R##_e = X##_e - Y##_e;				\
+							\
+    _FP_DIV_MEAT_##fs(R,X,Y);				\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NAN):		\
+    _FP_CHOOSENAN(fs, wc, R, X, Y, '/');		\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_NORMAL):	\
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_INF):		\
+  case _FP_CLS_COMBINE(FP_CLS_NAN,FP_CLS_ZERO):		\
+    R##_s = X##_s;					\
+    _FP_FRAC_COPY_##wc(R, X);				\
+    R##_c = X##_c;					\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_NAN):	\
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NAN):		\
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NAN):		\
+    R##_s = Y##_s;					\
+    _FP_FRAC_COPY_##wc(R, Y);				\
+    R##_c = Y##_c;					\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_INF):	\
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_INF):		\
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_NORMAL):	\
+    R##_c = FP_CLS_ZERO;				\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_NORMAL,FP_CLS_ZERO):	\
+    FP_SET_EXCEPTION(FP_EX_DIVZERO);			\
+	fallthrough;					\
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_ZERO):		\
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_NORMAL):	\
+    R##_c = FP_CLS_INF;					\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_INF,FP_CLS_INF):		\
+    R##_s = _FP_NANSIGN_##fs;				\
+    R##_c = FP_CLS_NAN;					\
+    _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);		\
+    FP_SET_EXCEPTION(FP_EX_INVALID | FP_EX_INVALID_IDI);\
+    break;						\
+							\
+  case _FP_CLS_COMBINE(FP_CLS_ZERO,FP_CLS_ZERO):	\
+    R##_s = _FP_NANSIGN_##fs;				\
+    R##_c = FP_CLS_NAN;					\
+    _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);		\
+    FP_SET_EXCEPTION(FP_EX_INVALID | FP_EX_INVALID_ZDZ);\
+    break;						\
+							\
+  default:						\
+    abort();						\
+  }							\
+} while (0)
+
+
+/*
+ * Main differential comparison routine.  The inputs should be raw not
+ * cooked.  The return is -1,0,1 for normal values, 2 otherwise.
+ */
+
+#define _FP_CMP(fs, wc, ret, X, Y, un)					\
+  do {									\
+    /* NANs are unordered */						\
+    if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X))		\
+	|| (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y)))	\
+      {									\
+	ret = un;							\
+      }									\
+    else								\
+      {									\
+	int __is_zero_x;						\
+	int __is_zero_y;						\
+									\
+	__is_zero_x = (!X##_e && _FP_FRAC_ZEROP_##wc(X)) ? 1 : 0;	\
+	__is_zero_y = (!Y##_e && _FP_FRAC_ZEROP_##wc(Y)) ? 1 : 0;	\
+									\
+	if (__is_zero_x && __is_zero_y)					\
+		ret = 0;						\
+	else if (__is_zero_x)						\
+		ret = Y##_s ? 1 : -1;					\
+	else if (__is_zero_y)						\
+		ret = X##_s ? -1 : 1;					\
+	else if (X##_s != Y##_s)					\
+	  ret = X##_s ? -1 : 1;						\
+	else if (X##_e > Y##_e)						\
+	  ret = X##_s ? -1 : 1;						\
+	else if (X##_e < Y##_e)						\
+	  ret = X##_s ? 1 : -1;						\
+	else if (_FP_FRAC_GT_##wc(X, Y))				\
+	  ret = X##_s ? -1 : 1;						\
+	else if (_FP_FRAC_GT_##wc(Y, X))				\
+	  ret = X##_s ? 1 : -1;						\
+	else								\
+	  ret = 0;							\
+      }									\
+  } while (0)
+
+
+/* Simplification for strict equality.  */
+
+#define _FP_CMP_EQ(fs, wc, ret, X, Y)					  \
+  do {									  \
+    /* NANs are unordered */						  \
+    if ((X##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(X))		  \
+	|| (Y##_e == _FP_EXPMAX_##fs && !_FP_FRAC_ZEROP_##wc(Y)))	  \
+      {									  \
+	ret = 1;							  \
+      }									  \
+    else								  \
+      {									  \
+	ret = !(X##_e == Y##_e						  \
+		&& _FP_FRAC_EQ_##wc(X, Y)				  \
+		&& (X##_s == Y##_s || !X##_e && _FP_FRAC_ZEROP_##wc(X))); \
+      }									  \
+  } while (0)
+
+/*
+ * Main square root routine.  The input value should be cooked.
+ */
+
+#define _FP_SQRT(fs, wc, R, X)						\
+do {									\
+    _FP_FRAC_DECL_##wc(T); _FP_FRAC_DECL_##wc(S);			\
+    _FP_W_TYPE q;							\
+    switch (X##_c)							\
+    {									\
+    case FP_CLS_NAN:							\
+	_FP_FRAC_COPY_##wc(R, X);					\
+	R##_s = X##_s;							\
+    	R##_c = FP_CLS_NAN;						\
+    	break;								\
+    case FP_CLS_INF:							\
+    	if (X##_s)							\
+    	  {								\
+    	    R##_s = _FP_NANSIGN_##fs;					\
+	    R##_c = FP_CLS_NAN; /* NAN */				\
+	    _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);			\
+	    FP_SET_EXCEPTION(FP_EX_INVALID);				\
+    	  }								\
+    	else								\
+    	  {								\
+    	    R##_s = 0;							\
+    	    R##_c = FP_CLS_INF; /* sqrt(+inf) = +inf */			\
+    	  }								\
+    	break;								\
+    case FP_CLS_ZERO:							\
+	R##_s = X##_s;							\
+	R##_c = FP_CLS_ZERO; /* sqrt(+-0) = +-0 */			\
+	break;								\
+    case FP_CLS_NORMAL:							\
+    	R##_s = 0;							\
+        if (X##_s)							\
+          {								\
+	    R##_c = FP_CLS_NAN; /* sNAN */				\
+	    R##_s = _FP_NANSIGN_##fs;					\
+	    _FP_FRAC_SET_##wc(R, _FP_NANFRAC_##fs);			\
+	    FP_SET_EXCEPTION(FP_EX_INVALID);				\
+	    break;							\
+          }								\
+    	R##_c = FP_CLS_NORMAL;						\
+        if (X##_e & 1)							\
+          _FP_FRAC_SLL_##wc(X, 1);					\
+        R##_e = X##_e >> 1;						\
+        _FP_FRAC_SET_##wc(S, _FP_ZEROFRAC_##wc);			\
+        _FP_FRAC_SET_##wc(R, _FP_ZEROFRAC_##wc);			\
+        q = _FP_OVERFLOW_##fs >> 1;					\
+        _FP_SQRT_MEAT_##wc(R, S, T, X, q);				\
+    }									\
+  } while (0)
+
+/*
+ * Convert from FP to integer
+ */
+
+/* RSIGNED can have following values:
+ * 0:  the number is required to be 0..(2^rsize)-1, if not, NV is set plus
+ *     the result is either 0 or (2^rsize)-1 depending on the sign in such case.
+ * 1:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not, NV is
+ *     set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1 depending
+ *     on the sign in such case.
+ * 2:  the number is required to be -(2^(rsize-1))..(2^(rsize-1))-1, if not, NV is
+ *     set plus the result is truncated to fit into destination.
+ * -1: the number is required to be -(2^(rsize-1))..(2^rsize)-1, if not, NV is
+ *     set plus the result is either -(2^(rsize-1)) or (2^(rsize-1))-1 depending
+ *     on the sign in such case.
+ */
+#define _FP_TO_INT(fs, wc, r, X, rsize, rsigned)				\
+  do {										\
+    switch (X##_c)								\
+      {										\
+      case FP_CLS_NORMAL:							\
+	if (X##_e < 0)								\
+	  {									\
+	    FP_SET_EXCEPTION(FP_EX_INEXACT);					\
+	    fallthrough;							\
+	  case FP_CLS_ZERO:							\
+	    r = 0;								\
+	  }									\
+	else if (X##_e >= rsize - (rsigned > 0 || X##_s)			\
+		 || (!rsigned && X##_s))					\
+	  {	/* overflow */							\
+	    fallthrough;							\
+	  case FP_CLS_NAN:                                                      \
+	  case FP_CLS_INF:							\
+	    if (rsigned == 2)							\
+	      {									\
+		if (X##_c != FP_CLS_NORMAL					\
+		    || X##_e >= rsize - 1 + _FP_WFRACBITS_##fs)			\
+		  r = 0;							\
+		else								\
+		  {								\
+		    _FP_FRAC_SLL_##wc(X, (X##_e - _FP_WFRACBITS_##fs + 1));	\
+		    _FP_FRAC_ASSEMBLE_##wc(r, X, rsize);			\
+		  }								\
+	      }									\
+	    else if (rsigned)							\
+	      {									\
+		r = 1;								\
+		r <<= rsize - 1;						\
+		r -= 1 - X##_s;							\
+	      }									\
+	    else								\
+	      {									\
+		r = 0;								\
+		if (!X##_s)							\
+		  r = ~r;							\
+	      }									\
+	    FP_SET_EXCEPTION(FP_EX_INVALID);					\
+	  }									\
+	else									\
+	  {									\
+	    if (_FP_W_TYPE_SIZE*wc < rsize)					\
+	      {									\
+		_FP_FRAC_ASSEMBLE_##wc(r, X, rsize);				\
+		r <<= X##_e - _FP_WFRACBITS_##fs;				\
+	      }									\
+	    else								\
+	      {									\
+		if (X##_e >= _FP_WFRACBITS_##fs)				\
+		  _FP_FRAC_SLL_##wc(X, (X##_e - _FP_WFRACBITS_##fs + 1));	\
+		else if (X##_e < _FP_WFRACBITS_##fs - 1)			\
+		  {								\
+		    _FP_FRAC_SRS_##wc(X, (_FP_WFRACBITS_##fs - X##_e - 2),	\
+				      _FP_WFRACBITS_##fs);			\
+		    if (_FP_FRAC_LOW_##wc(X) & 1)				\
+		      FP_SET_EXCEPTION(FP_EX_INEXACT);				\
+		    _FP_FRAC_SRL_##wc(X, 1);					\
+		  }								\
+		_FP_FRAC_ASSEMBLE_##wc(r, X, rsize);				\
+	      }									\
+	    if (rsigned && X##_s)						\
+	      r = -r;								\
+	  }									\
+	break;									\
+      }										\
+  } while (0)
+
+#define _FP_TO_INT_ROUND(fs, wc, r, X, rsize, rsigned)				\
+  do {										\
+    r = 0;									\
+    switch (X##_c)								\
+      {										\
+      case FP_CLS_NORMAL:							\
+	if (X##_e >= _FP_FRACBITS_##fs - 1)					\
+	  {									\
+	    if (X##_e < rsize - 1 + _FP_WFRACBITS_##fs)				\
+	      {									\
+		if (X##_e >= _FP_WFRACBITS_##fs - 1)				\
+		  {								\
+		    _FP_FRAC_ASSEMBLE_##wc(r, X, rsize);			\
+		    r <<= X##_e - _FP_WFRACBITS_##fs + 1;			\
+		  }								\
+		else								\
+		  {								\
+		    _FP_FRAC_SRL_##wc(X, _FP_WORKBITS - X##_e			\
+				      + _FP_FRACBITS_##fs - 1);			\
+		    _FP_FRAC_ASSEMBLE_##wc(r, X, rsize);			\
+		  }								\
+	      }									\
+	  }									\
+	else									\
+	  {									\
+	    int _lz0, _lz1;							\
+	    if (X##_e <= -_FP_WORKBITS - 1)					\
+	      _FP_FRAC_SET_##wc(X, _FP_MINFRAC_##wc);				\
+	    else								\
+	      _FP_FRAC_SRS_##wc(X, _FP_FRACBITS_##fs - 1 - X##_e,		\
+				_FP_WFRACBITS_##fs);				\
+	    _FP_FRAC_CLZ_##wc(_lz0, X);						\
+	    _FP_ROUND(wc, X);							\
+	    _FP_FRAC_CLZ_##wc(_lz1, X);						\
+	    if (_lz1 < _lz0)							\
+	      X##_e++; /* For overflow detection.  */				\
+	    _FP_FRAC_SRL_##wc(X, _FP_WORKBITS);					\
+	    _FP_FRAC_ASSEMBLE_##wc(r, X, rsize);				\
+	  }									\
+	if (rsigned && X##_s)							\
+	  r = -r;								\
+	if (X##_e >= rsize - (rsigned > 0 || X##_s)				\
+	    || (!rsigned && X##_s))						\
+	  {	/* overflow */							\
+	    fallthrough;							\
+	  case FP_CLS_NAN:                                                      \
+	  case FP_CLS_INF:							\
+	    if (!rsigned)							\
+	      {									\
+		r = 0;								\
+		if (!X##_s)							\
+		  r = ~r;							\
+	      }									\
+	    else if (rsigned != 2)						\
+	      {									\
+		r = 1;								\
+		r <<= rsize - 1;						\
+		r -= 1 - X##_s;							\
+	      }									\
+	    FP_SET_EXCEPTION(FP_EX_INVALID);					\
+	  }									\
+	break;									\
+      case FP_CLS_ZERO:								\
+        break;									\
+      }										\
+  } while (0)
+
+#define _FP_FROM_INT(fs, wc, X, r, rsize, rtype)			\
+  do {									\
+    if (r)								\
+      {									\
+        unsigned rtype ur_;						\
+	X##_c = FP_CLS_NORMAL;						\
+									\
+	if ((X##_s = (r < 0)))						\
+	  ur_ = (unsigned rtype) -r;					\
+	else								\
+	  ur_ = (unsigned rtype) r;					\
+	(void) (((rsize) <= _FP_W_TYPE_SIZE)				\
+		? ({ __FP_CLZ(X##_e, ur_); })				\
+		: ({							\
+		     __FP_CLZ_2(X##_e, (_FP_W_TYPE)(ur_ >> _FP_W_TYPE_SIZE),  \
+							    (_FP_W_TYPE)ur_); \
+		  }));							\
+	if (rsize < _FP_W_TYPE_SIZE)					\
+		X##_e -= (_FP_W_TYPE_SIZE - rsize);			\
+	X##_e = rsize - X##_e - 1;					\
+									\
+	if (_FP_FRACBITS_##fs < rsize && _FP_WFRACBITS_##fs <= X##_e)	\
+	  __FP_FRAC_SRS_1(ur_, (X##_e - _FP_WFRACBITS_##fs + 1), rsize);\
+	_FP_FRAC_DISASSEMBLE_##wc(X, ur_, rsize);			\
+	if ((_FP_WFRACBITS_##fs - X##_e - 1) > 0)			\
+	  _FP_FRAC_SLL_##wc(X, (_FP_WFRACBITS_##fs - X##_e - 1));	\
+      }									\
+    else								\
+      {									\
+	X##_c = FP_CLS_ZERO, X##_s = 0;					\
+      }									\
+  } while (0)
+
+
+#define FP_CONV(dfs,sfs,dwc,swc,D,S)			\
+  do {							\
+    _FP_FRAC_CONV_##dwc##_##swc(dfs, sfs, D, S);	\
+    D##_e = S##_e;					\
+    D##_c = S##_c;					\
+    D##_s = S##_s;					\
+  } while (0)
+
+/*
+ * Helper primitives.
+ */
+
+/* Count leading zeros in a word.  */
+
+#ifndef __FP_CLZ
+#if _FP_W_TYPE_SIZE < 64
+/* this is just to shut the compiler up about shifts > word length -- PMM 02/1998 */
+#define __FP_CLZ(r, x)				\
+  do {						\
+    _FP_W_TYPE _t = (x);			\
+    r = _FP_W_TYPE_SIZE - 1;			\
+    if (_t > 0xffff) r -= 16;			\
+    if (_t > 0xffff) _t >>= 16;			\
+    if (_t > 0xff) r -= 8;			\
+    if (_t > 0xff) _t >>= 8;			\
+    if (_t & 0xf0) r -= 4;			\
+    if (_t & 0xf0) _t >>= 4;			\
+    if (_t & 0xc) r -= 2;			\
+    if (_t & 0xc) _t >>= 2;			\
+    if (_t & 0x2) r -= 1;			\
+  } while (0)
+#else /* not _FP_W_TYPE_SIZE < 64 */
+#define __FP_CLZ(r, x)				\
+  do {						\
+    _FP_W_TYPE _t = (x);			\
+    r = _FP_W_TYPE_SIZE - 1;			\
+    if (_t > 0xffffffff) r -= 32;		\
+    if (_t > 0xffffffff) _t >>= 32;		\
+    if (_t > 0xffff) r -= 16;			\
+    if (_t > 0xffff) _t >>= 16;			\
+    if (_t > 0xff) r -= 8;			\
+    if (_t > 0xff) _t >>= 8;			\
+    if (_t & 0xf0) r -= 4;			\
+    if (_t & 0xf0) _t >>= 4;			\
+    if (_t & 0xc) r -= 2;			\
+    if (_t & 0xc) _t >>= 2;			\
+    if (_t & 0x2) r -= 1;			\
+  } while (0)
+#endif /* not _FP_W_TYPE_SIZE < 64 */
+#endif /* ndef __FP_CLZ */
+
+#define _FP_DIV_HELP_imm(q, r, n, d)		\
+  do {						\
+    q = n / d, r = n % d;			\
+  } while (0)
+
+#endif /* __MATH_EMU_OP_COMMON_H__ */
diff --git a/include/math-emu/quad.h b/include/math-emu/quad.h
new file mode 100644
index 000000000..616113626
--- /dev/null
+++ b/include/math-emu/quad.h
@@ -0,0 +1,208 @@
+/* Software floating-point emulation.
+   Definitions for IEEE Quad Precision.
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef  __MATH_EMU_QUAD_H__
+#define  __MATH_EMU_QUAD_H__
+
+#if _FP_W_TYPE_SIZE < 32
+#error "Here's a nickel, kid. Go buy yourself a real computer."
+#endif
+
+#if _FP_W_TYPE_SIZE < 64
+#define _FP_FRACTBITS_Q         (4*_FP_W_TYPE_SIZE)
+#else
+#define _FP_FRACTBITS_Q		(2*_FP_W_TYPE_SIZE)
+#endif
+
+#define _FP_FRACBITS_Q		113
+#define _FP_FRACXBITS_Q		(_FP_FRACTBITS_Q - _FP_FRACBITS_Q)
+#define _FP_WFRACBITS_Q		(_FP_WORKBITS + _FP_FRACBITS_Q)
+#define _FP_WFRACXBITS_Q	(_FP_FRACTBITS_Q - _FP_WFRACBITS_Q)
+#define _FP_EXPBITS_Q		15
+#define _FP_EXPBIAS_Q		16383
+#define _FP_EXPMAX_Q		32767
+
+#define _FP_QNANBIT_Q		\
+	((_FP_W_TYPE)1 << (_FP_FRACBITS_Q-2) % _FP_W_TYPE_SIZE)
+#define _FP_IMPLBIT_Q		\
+	((_FP_W_TYPE)1 << (_FP_FRACBITS_Q-1) % _FP_W_TYPE_SIZE)
+#define _FP_OVERFLOW_Q		\
+	((_FP_W_TYPE)1 << (_FP_WFRACBITS_Q % _FP_W_TYPE_SIZE))
+
+#if _FP_W_TYPE_SIZE < 64
+
+union _FP_UNION_Q
+{
+   long double flt;
+   struct 
+   {
+#if __BYTE_ORDER == __BIG_ENDIAN
+      unsigned sign : 1;
+      unsigned exp : _FP_EXPBITS_Q;
+      unsigned long frac3 : _FP_FRACBITS_Q - (_FP_IMPLBIT_Q != 0)-(_FP_W_TYPE_SIZE * 3);
+      unsigned long frac2 : _FP_W_TYPE_SIZE;
+      unsigned long frac1 : _FP_W_TYPE_SIZE;
+      unsigned long frac0 : _FP_W_TYPE_SIZE;
+#else
+      unsigned long frac0 : _FP_W_TYPE_SIZE;
+      unsigned long frac1 : _FP_W_TYPE_SIZE;
+      unsigned long frac2 : _FP_W_TYPE_SIZE;
+      unsigned long frac3 : _FP_FRACBITS_Q - (_FP_IMPLBIT_Q != 0)-(_FP_W_TYPE_SIZE * 3);
+      unsigned exp : _FP_EXPBITS_Q;
+      unsigned sign : 1;
+#endif /* not bigendian */
+   } bits __attribute__((packed));
+};
+
+
+#define FP_DECL_Q(X)		_FP_DECL(4,X)
+#define FP_UNPACK_RAW_Q(X,val)	_FP_UNPACK_RAW_4(Q,X,val)
+#define FP_UNPACK_RAW_QP(X,val)	_FP_UNPACK_RAW_4_P(Q,X,val)
+#define FP_PACK_RAW_Q(val,X)	_FP_PACK_RAW_4(Q,val,X)
+#define FP_PACK_RAW_QP(val,X)		\
+  do {					\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_4_P(Q,val,X);	\
+  } while (0)
+
+#define FP_UNPACK_Q(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_4(Q,X,val);		\
+    _FP_UNPACK_CANONICAL(Q,4,X);	\
+  } while (0)
+
+#define FP_UNPACK_QP(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_4_P(Q,X,val);	\
+    _FP_UNPACK_CANONICAL(Q,4,X);	\
+  } while (0)
+
+#define FP_PACK_Q(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(Q,4,X);		\
+    _FP_PACK_RAW_4(Q,val,X);		\
+  } while (0)
+
+#define FP_PACK_QP(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(Q,4,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_4_P(Q,val,X);	\
+  } while (0)
+
+#define FP_ISSIGNAN_Q(X)		_FP_ISSIGNAN(Q,4,X)
+#define FP_NEG_Q(R,X)			_FP_NEG(Q,4,R,X)
+#define FP_ADD_Q(R,X,Y)			_FP_ADD(Q,4,R,X,Y)
+#define FP_SUB_Q(R,X,Y)			_FP_SUB(Q,4,R,X,Y)
+#define FP_MUL_Q(R,X,Y)			_FP_MUL(Q,4,R,X,Y)
+#define FP_DIV_Q(R,X,Y)			_FP_DIV(Q,4,R,X,Y)
+#define FP_SQRT_Q(R,X)			_FP_SQRT(Q,4,R,X)
+#define _FP_SQRT_MEAT_Q(R,S,T,X,Q)	_FP_SQRT_MEAT_4(R,S,T,X,Q)
+
+#define FP_CMP_Q(r,X,Y,un)	_FP_CMP(Q,4,r,X,Y,un)
+#define FP_CMP_EQ_Q(r,X,Y)	_FP_CMP_EQ(Q,4,r,X,Y)
+
+#define FP_TO_INT_Q(r,X,rsz,rsg)	_FP_TO_INT(Q,4,r,X,rsz,rsg)
+#define FP_TO_INT_ROUND_Q(r,X,rsz,rsg)	_FP_TO_INT_ROUND(Q,4,r,X,rsz,rsg)
+#define FP_FROM_INT_Q(X,r,rs,rt)	_FP_FROM_INT(Q,4,X,r,rs,rt)
+
+#define _FP_FRAC_HIGH_Q(X)	_FP_FRAC_HIGH_4(X)
+#define _FP_FRAC_HIGH_RAW_Q(X)	_FP_FRAC_HIGH_4(X)
+
+#else   /* not _FP_W_TYPE_SIZE < 64 */
+union _FP_UNION_Q
+{
+  long double flt /* __attribute__((mode(TF))) */ ;
+  struct {
+#if __BYTE_ORDER == __BIG_ENDIAN
+    unsigned sign  : 1;
+    unsigned exp   : _FP_EXPBITS_Q;
+    unsigned long frac1 : _FP_FRACBITS_Q-(_FP_IMPLBIT_Q != 0)-_FP_W_TYPE_SIZE;
+    unsigned long frac0 : _FP_W_TYPE_SIZE;
+#else
+    unsigned long frac0 : _FP_W_TYPE_SIZE;
+    unsigned long frac1 : _FP_FRACBITS_Q-(_FP_IMPLBIT_Q != 0)-_FP_W_TYPE_SIZE;
+    unsigned exp   : _FP_EXPBITS_Q;
+    unsigned sign  : 1;
+#endif
+  } bits;
+};
+
+#define FP_DECL_Q(X)		_FP_DECL(2,X)
+#define FP_UNPACK_RAW_Q(X,val)	_FP_UNPACK_RAW_2(Q,X,val)
+#define FP_UNPACK_RAW_QP(X,val)	_FP_UNPACK_RAW_2_P(Q,X,val)
+#define FP_PACK_RAW_Q(val,X)	_FP_PACK_RAW_2(Q,val,X)
+#define FP_PACK_RAW_QP(val,X)		\
+  do {					\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_2_P(Q,val,X);	\
+  } while (0)
+
+#define FP_UNPACK_Q(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_2(Q,X,val);		\
+    _FP_UNPACK_CANONICAL(Q,2,X);	\
+  } while (0)
+
+#define FP_UNPACK_QP(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_2_P(Q,X,val);	\
+    _FP_UNPACK_CANONICAL(Q,2,X);	\
+  } while (0)
+
+#define FP_PACK_Q(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(Q,2,X);		\
+    _FP_PACK_RAW_2(Q,val,X);		\
+  } while (0)
+
+#define FP_PACK_QP(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(Q,2,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_2_P(Q,val,X);	\
+  } while (0)
+
+#define FP_ISSIGNAN_Q(X)		_FP_ISSIGNAN(Q,2,X)
+#define FP_NEG_Q(R,X)			_FP_NEG(Q,2,R,X)
+#define FP_ADD_Q(R,X,Y)			_FP_ADD(Q,2,R,X,Y)
+#define FP_SUB_Q(R,X,Y)			_FP_SUB(Q,2,R,X,Y)
+#define FP_MUL_Q(R,X,Y)			_FP_MUL(Q,2,R,X,Y)
+#define FP_DIV_Q(R,X,Y)			_FP_DIV(Q,2,R,X,Y)
+#define FP_SQRT_Q(R,X)			_FP_SQRT(Q,2,R,X)
+#define _FP_SQRT_MEAT_Q(R,S,T,X,Q)	_FP_SQRT_MEAT_2(R,S,T,X,Q)
+
+#define FP_CMP_Q(r,X,Y,un)	_FP_CMP(Q,2,r,X,Y,un)
+#define FP_CMP_EQ_Q(r,X,Y)	_FP_CMP_EQ(Q,2,r,X,Y)
+
+#define FP_TO_INT_Q(r,X,rsz,rsg)	_FP_TO_INT(Q,2,r,X,rsz,rsg)
+#define FP_TO_INT_ROUND_Q(r,X,rsz,rsg)	_FP_TO_INT_ROUND(Q,2,r,X,rsz,rsg)
+#define FP_FROM_INT_Q(X,r,rs,rt)	_FP_FROM_INT(Q,2,X,r,rs,rt)
+
+#define _FP_FRAC_HIGH_Q(X)	_FP_FRAC_HIGH_2(X)
+#define _FP_FRAC_HIGH_RAW_Q(X)	_FP_FRAC_HIGH_2(X)
+
+#endif /* not _FP_W_TYPE_SIZE < 64 */
+
+#endif /* __MATH_EMU_QUAD_H__ */
diff --git a/include/math-emu/single.h b/include/math-emu/single.h
new file mode 100644
index 000000000..87f90b0f9
--- /dev/null
+++ b/include/math-emu/single.h
@@ -0,0 +1,116 @@
+/* Software floating-point emulation.
+   Definitions for IEEE Single Precision.
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef    __MATH_EMU_SINGLE_H__
+#define    __MATH_EMU_SINGLE_H__
+
+#if _FP_W_TYPE_SIZE < 32
+#error "Here's a nickel kid.  Go buy yourself a real computer."
+#endif
+
+#define _FP_FRACBITS_S		24
+#define _FP_FRACXBITS_S		(_FP_W_TYPE_SIZE - _FP_FRACBITS_S)
+#define _FP_WFRACBITS_S		(_FP_WORKBITS + _FP_FRACBITS_S)
+#define _FP_WFRACXBITS_S	(_FP_W_TYPE_SIZE - _FP_WFRACBITS_S)
+#define _FP_EXPBITS_S		8
+#define _FP_EXPBIAS_S		127
+#define _FP_EXPMAX_S		255
+#define _FP_QNANBIT_S		((_FP_W_TYPE)1 << (_FP_FRACBITS_S-2))
+#define _FP_IMPLBIT_S		((_FP_W_TYPE)1 << (_FP_FRACBITS_S-1))
+#define _FP_OVERFLOW_S		((_FP_W_TYPE)1 << (_FP_WFRACBITS_S))
+
+/* The implementation of _FP_MUL_MEAT_S and _FP_DIV_MEAT_S should be
+   chosen by the target machine.  */
+
+union _FP_UNION_S
+{
+  float flt;
+  struct {
+#if __BYTE_ORDER == __BIG_ENDIAN
+    unsigned sign : 1;
+    unsigned exp  : _FP_EXPBITS_S;
+    unsigned frac : _FP_FRACBITS_S - (_FP_IMPLBIT_S != 0);
+#else
+    unsigned frac : _FP_FRACBITS_S - (_FP_IMPLBIT_S != 0);
+    unsigned exp  : _FP_EXPBITS_S;
+    unsigned sign : 1;
+#endif
+  } bits __attribute__((packed));
+};
+
+#define FP_DECL_S(X)		_FP_DECL(1,X)
+#define FP_UNPACK_RAW_S(X,val)	_FP_UNPACK_RAW_1(S,X,val)
+#define FP_UNPACK_RAW_SP(X,val)	_FP_UNPACK_RAW_1_P(S,X,val)
+#define FP_PACK_RAW_S(val,X)	_FP_PACK_RAW_1(S,val,X)
+#define FP_PACK_RAW_SP(val,X)		\
+  do {					\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_1_P(S,val,X);	\
+  } while (0)
+
+#define FP_UNPACK_S(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_1(S,X,val);		\
+    _FP_UNPACK_CANONICAL(S,1,X);	\
+  } while (0)
+
+#define FP_UNPACK_SP(X,val)		\
+  do {					\
+    _FP_UNPACK_RAW_1_P(S,X,val);	\
+    _FP_UNPACK_CANONICAL(S,1,X);	\
+  } while (0)
+
+#define FP_PACK_S(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(S,1,X);		\
+    _FP_PACK_RAW_1(S,val,X);		\
+  } while (0)
+
+#define FP_PACK_SP(val,X)		\
+  do {					\
+    _FP_PACK_CANONICAL(S,1,X);		\
+    if (!FP_INHIBIT_RESULTS)		\
+      _FP_PACK_RAW_1_P(S,val,X);	\
+  } while (0)
+
+#define FP_ISSIGNAN_S(X)		_FP_ISSIGNAN(S,1,X)
+#define FP_NEG_S(R,X)			_FP_NEG(S,1,R,X)
+#define FP_ADD_S(R,X,Y)			_FP_ADD(S,1,R,X,Y)
+#define FP_SUB_S(R,X,Y)			_FP_SUB(S,1,R,X,Y)
+#define FP_MUL_S(R,X,Y)			_FP_MUL(S,1,R,X,Y)
+#define FP_DIV_S(R,X,Y)			_FP_DIV(S,1,R,X,Y)
+#define FP_SQRT_S(R,X)			_FP_SQRT(S,1,R,X)
+#define _FP_SQRT_MEAT_S(R,S,T,X,Q)	_FP_SQRT_MEAT_1(R,S,T,X,Q)
+
+#define FP_CMP_S(r,X,Y,un)	_FP_CMP(S,1,r,X,Y,un)
+#define FP_CMP_EQ_S(r,X,Y)	_FP_CMP_EQ(S,1,r,X,Y)
+
+#define FP_TO_INT_S(r,X,rsz,rsg)	_FP_TO_INT(S,1,r,X,rsz,rsg)
+#define FP_TO_INT_ROUND_S(r,X,rsz,rsg)	_FP_TO_INT_ROUND(S,1,r,X,rsz,rsg)
+#define FP_FROM_INT_S(X,r,rs,rt)	_FP_FROM_INT(S,1,X,r,rs,rt)
+
+#define _FP_FRAC_HIGH_S(X)	_FP_FRAC_HIGH_1(X)
+#define _FP_FRAC_HIGH_RAW_S(X)	_FP_FRAC_HIGH_1(X)
+
+#endif /* __MATH_EMU_SINGLE_H__ */
diff --git a/include/math-emu/soft-fp.h b/include/math-emu/soft-fp.h
new file mode 100644
index 000000000..5650c1628
--- /dev/null
+++ b/include/math-emu/soft-fp.h
@@ -0,0 +1,207 @@
+/* Software floating-point emulation.
+   Copyright (C) 1997,1998,1999 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Richard Henderson (rth@cygnus.com),
+		  Jakub Jelinek (jj@ultra.linux.cz),
+		  David S. Miller (davem@redhat.com) and
+		  Peter Maydell (pmaydell@chiark.greenend.org.uk).
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License as
+   published by the Free Software Foundation; either version 2 of the
+   License, or (at your option) any later version.
+
+   The GNU C Library 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If
+   not, write to the Free Software Foundation, Inc.,
+   59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
+
+#ifndef __MATH_EMU_SOFT_FP_H__
+#define __MATH_EMU_SOFT_FP_H__
+
+#include <asm/sfp-machine.h>
+
+/* Allow sfp-machine to have its own byte order definitions. */
+#ifndef __BYTE_ORDER
+#include <endian.h>
+#endif
+
+#define _FP_WORKBITS		3
+#define _FP_WORK_LSB		((_FP_W_TYPE)1 << 3)
+#define _FP_WORK_ROUND		((_FP_W_TYPE)1 << 2)
+#define _FP_WORK_GUARD		((_FP_W_TYPE)1 << 1)
+#define _FP_WORK_STICKY		((_FP_W_TYPE)1 << 0)
+
+#ifndef FP_RND_NEAREST
+# define FP_RND_NEAREST		0
+# define FP_RND_ZERO		1
+# define FP_RND_PINF		2
+# define FP_RND_MINF		3
+#ifndef FP_ROUNDMODE
+# define FP_ROUNDMODE		FP_RND_NEAREST
+#endif
+#endif
+
+/* By default don't care about exceptions. */
+#ifndef FP_EX_INVALID
+#define FP_EX_INVALID		0
+#endif
+#ifndef FP_EX_INVALID_SNAN
+#define FP_EX_INVALID_SNAN	0
+#endif
+/* inf - inf */
+#ifndef FP_EX_INVALID_ISI
+#define FP_EX_INVALID_ISI	0
+#endif
+/* inf / inf */
+#ifndef FP_EX_INVALID_IDI
+#define FP_EX_INVALID_IDI	0
+#endif
+/* 0 / 0 */
+#ifndef FP_EX_INVALID_ZDZ
+#define FP_EX_INVALID_ZDZ	0
+#endif
+/* inf * 0 */
+#ifndef FP_EX_INVALID_IMZ
+#define FP_EX_INVALID_IMZ	0
+#endif
+#ifndef FP_EX_OVERFLOW
+#define FP_EX_OVERFLOW		0
+#endif
+#ifndef FP_EX_UNDERFLOW
+#define FP_EX_UNDERFLOW		
+#endif
+#ifndef FP_EX_DIVZERO
+#define FP_EX_DIVZERO		0
+#endif
+#ifndef FP_EX_INEXACT
+#define FP_EX_INEXACT		0
+#endif
+#ifndef FP_EX_DENORM
+#define FP_EX_DENORM		0
+#endif
+
+#ifdef _FP_DECL_EX
+#define FP_DECL_EX					\
+  int _fex = 0;						\
+  _FP_DECL_EX
+#else
+#define FP_DECL_EX int _fex = 0
+#endif
+  
+#ifndef FP_INIT_ROUNDMODE
+#define FP_INIT_ROUNDMODE do {} while (0)
+#endif
+
+#ifndef FP_HANDLE_EXCEPTIONS
+#define FP_HANDLE_EXCEPTIONS do {} while (0)
+#endif
+
+/* By default we never flush denormal input operands to signed zero. */
+#ifndef FP_DENORM_ZERO
+#define FP_DENORM_ZERO 0
+#endif
+
+#ifndef FP_INHIBIT_RESULTS
+/* By default we write the results always.
+ * sfp-machine may override this and e.g.
+ * check if some exceptions are unmasked
+ * and inhibit it in such a case.
+ */
+#define FP_INHIBIT_RESULTS 0
+#endif
+
+#ifndef FP_TRAPPING_EXCEPTIONS
+#define FP_TRAPPING_EXCEPTIONS 0
+#endif
+
+#define FP_SET_EXCEPTION(ex)				\
+  _fex |= (ex)
+  
+#define FP_UNSET_EXCEPTION(ex)				\
+  _fex &= ~(ex)
+
+#define FP_CUR_EXCEPTIONS				\
+  (_fex)
+
+#define FP_CLEAR_EXCEPTIONS				\
+  _fex = 0
+
+#define _FP_ROUND_NEAREST(wc, X)			\
+do {							\
+    if ((_FP_FRAC_LOW_##wc(X) & 15) != _FP_WORK_ROUND)	\
+      _FP_FRAC_ADDI_##wc(X, _FP_WORK_ROUND);		\
+} while (0)
+
+#define _FP_ROUND_ZERO(wc, X)		(void)0
+
+#define _FP_ROUND_PINF(wc, X)				\
+do {							\
+    if (!X##_s && (_FP_FRAC_LOW_##wc(X) & 7))		\
+      _FP_FRAC_ADDI_##wc(X, _FP_WORK_LSB);		\
+} while (0)
+
+#define _FP_ROUND_MINF(wc, X)				\
+do {							\
+    if (X##_s && (_FP_FRAC_LOW_##wc(X) & 7))		\
+      _FP_FRAC_ADDI_##wc(X, _FP_WORK_LSB);		\
+} while (0)
+
+#define _FP_ROUND(wc, X)			\
+do {						\
+	if (_FP_FRAC_LOW_##wc(X) & 7)		\
+	  FP_SET_EXCEPTION(FP_EX_INEXACT);	\
+	switch (FP_ROUNDMODE)			\
+	{					\
+	  case FP_RND_NEAREST:			\
+	    _FP_ROUND_NEAREST(wc,X);		\
+	    break;				\
+	  case FP_RND_ZERO:			\
+	    _FP_ROUND_ZERO(wc,X);		\
+	    break;				\
+	  case FP_RND_PINF:			\
+	    _FP_ROUND_PINF(wc,X);		\
+	    break;				\
+	  case FP_RND_MINF:			\
+	    _FP_ROUND_MINF(wc,X);		\
+	    break;				\
+	}					\
+} while (0)
+
+#define FP_CLS_NORMAL		0
+#define FP_CLS_ZERO		1
+#define FP_CLS_INF		2
+#define FP_CLS_NAN		3
+
+#define _FP_CLS_COMBINE(x,y)	(((x) << 2) | (y))
+
+#include <math-emu/op-1.h>
+#include <math-emu/op-2.h>
+#include <math-emu/op-4.h>
+#include <math-emu/op-8.h>
+#include <math-emu/op-common.h>
+
+/* Sigh.  Silly things longlong.h needs.  */
+#define UWtype		_FP_W_TYPE
+#define W_TYPE_SIZE	_FP_W_TYPE_SIZE
+
+typedef int SItype __attribute__((mode(SI)));
+typedef int DItype __attribute__((mode(DI)));
+typedef unsigned int USItype __attribute__((mode(SI)));
+typedef unsigned int UDItype __attribute__((mode(DI)));
+#if _FP_W_TYPE_SIZE == 32
+typedef unsigned int UHWtype __attribute__((mode(HI)));
+#elif _FP_W_TYPE_SIZE == 64
+typedef USItype UHWtype;
+#endif
+
+#ifndef umul_ppmm
+#include <stdlib/longlong.h>
+#endif
+
+#endif /* __MATH_EMU_SOFT_FP_H__ */