Initial revision

git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@8496 138bc75d-0d04-0410-961f-82ee72b054a4

9 files changed, 3892 insertions, 0 deletions

diff --git a/gcc/config/m68k/fpgnulib.c b/gcc/config/m68k/fpgnulib.c
new file mode 100644
index 00000000000..bd54058e77e
--- /dev/null
+++ b/gcc/config/m68k/fpgnulib.c
@@ -0,0 +1,439 @@
+/* This is a stripped down version of floatlib.c.  It supplies only those
+   functions which exist in libgcc, but for which there is not assembly
+   language versions in m68k/lb1sf68.asm.
+
+   It also includes simplistic support for extended floats (by working in
+   double precision).  You must compile this file again with -DEXTFLOAT
+   to get this support.  */
+
+/*
+** gnulib support for software floating point.
+** Copyright (C) 1991 by Pipeline Associates, Inc.  All rights reserved.
+** Permission is granted to do *anything* you want with this file,
+** commercial or otherwise, provided this message remains intact.  So there!
+** I would appreciate receiving any updates/patches/changes that anyone
+** makes, and am willing to be the repository for said changes (am I
+** making a big mistake?).
+**
+** Pat Wood
+** Pipeline Associates, Inc.
+** pipeline!phw@motown.com or
+** sun!pipeline!phw or
+** uunet!motown!pipeline!phw
+**
+** 05/01/91 -- V1.0 -- first release to gcc mailing lists
+** 05/04/91 -- V1.1 -- added float and double prototypes and return values
+**                  -- fixed problems with adding and subtracting zero
+**                  -- fixed rounding in truncdfsf2
+**                  -- fixed SWAP define and tested on 386
+*/
+
+/*
+** The following are routines that replace the gnulib soft floating point
+** routines that are called automatically when -msoft-float is selected.
+** The support single and double precision IEEE format, with provisions
+** for byte-swapped machines (tested on 386).  Some of the double-precision
+** routines work at full precision, but most of the hard ones simply punt
+** and call the single precision routines, producing a loss of accuracy.
+** long long support is not assumed or included.
+** Overall accuracy is close to IEEE (actually 68882) for single-precision
+** arithmetic.  I think there may still be a 1 in 1000 chance of a bit
+** being rounded the wrong way during a multiply.  I'm not fussy enough to
+** bother with it, but if anyone is, knock yourself out.
+**
+** Efficiency has only been addressed where it was obvious that something
+** would make a big difference.  Anyone who wants to do this right for
+** best speed should go in and rewrite in assembler.
+**
+** I have tested this only on a 68030 workstation and 386/ix integrated
+** in with -msoft-float.
+*/
+
+/* the following deal with IEEE single-precision numbers */
+#define EXCESS		126L
+#define SIGNBIT		0x80000000L
+#define HIDDEN		(1L << 23L)
+#define SIGN(fp)	((fp) & SIGNBIT)
+#define EXP(fp)		(((fp) >> 23L) & 0xFF)
+#define MANT(fp)	(((fp) & 0x7FFFFFL) | HIDDEN)
+#define PACK(s,e,m)	((s) | ((e) << 23L) | (m))
+
+/* the following deal with IEEE double-precision numbers */
+#define EXCESSD		1022
+#define HIDDEND		(1L << 20L)
+#define EXPDBITS	11
+#define EXPDMASK	0x7FF
+#define EXPD(fp)	(((fp.l.upper) >> 20L) & 0x7FFL)
+#define SIGND(fp)	((fp.l.upper) & SIGNBIT)
+#define MANTD(fp)	(((((fp.l.upper) & 0xFFFFF) | HIDDEND) << 10) | \
+				(fp.l.lower >> 22))
+#define MANTDMASK	0xFFFFF /* mask of upper part */
+
+/* the following deal with IEEE extended-precision numbers */
+#define EXCESSX		16382
+#define HIDDENX		(1L << 31L)
+#define EXPXBITS	15
+#define EXPXMASK	0x7FFF
+#define EXPX(fp)	(((fp.l.upper) >> 16) & EXPXMASK)
+#define SIGNX(fp)	((fp.l.upper) & SIGNBIT)
+#define MANTXMASK	0x7FFFFFFF /* mask of upper part */
+
+union double_long 
+{
+  double d;
+  struct {
+      long upper;
+      unsigned long lower;
+    } l;
+};
+
+union float_long {
+  float f;
+  long l;
+};
+
+union long_double_long
+{
+  long double ld;
+  struct
+    {
+      long upper;
+      unsigned long middle;
+      unsigned long lower;
+    } l;
+};
+
+#ifndef EXTFLOAT
+
+/* convert int to double */
+double
+__floatsidf (int a1)
+{
+  long sign = 0, exp = 31 + EXCESSD;
+  union double_long dl;
+
+  if (!a1)
+    {
+      dl.l.upper = dl.l.lower = 0;
+      return dl.d;
+    }
+
+  if (a1 < 0)
+    {
+      sign = SIGNBIT;
+      a1 = -a1;
+      if (a1 < 0)
+	{
+	  dl.l.upper = SIGNBIT | ((32 + EXCESSD) << 20L);
+	  dl.l.lower = 0;
+	  return dl.d;
+        }
+    }
+
+  while (a1 < 0x1000000)
+    {
+      a1 <<= 4;
+      exp -= 4;
+    }
+
+  while (a1 < 0x40000000)
+    {
+      a1 <<= 1;
+      exp--;
+    }
+
+  /* pack up and go home */
+  dl.l.upper = sign;
+  dl.l.upper |= exp << 20L;
+  dl.l.upper |= (a1 >> 10L) & ~HIDDEND;
+  dl.l.lower = a1 << 22L;
+
+  return dl.d;
+}
+
+/* convert int to float */
+float
+__floatsisf (int l)
+{
+  double foo = __floatsidf (l);
+  return foo;
+}
+
+/* convert float to double */
+double
+__extendsfdf2 (float a1)
+{
+  register union float_long fl1;
+  register union double_long dl;
+  register long exp;
+
+  fl1.f = a1;
+
+  if (!fl1.l)
+    {
+      dl.l.upper = dl.l.lower = 0;
+      return dl.d;
+    }
+
+  dl.l.upper = SIGN (fl1.l);
+  exp = EXP (fl1.l) - EXCESS + EXCESSD;
+  dl.l.upper |= exp << 20;
+  dl.l.upper |= (MANT (fl1.l) & ~HIDDEN) >> 3;
+  dl.l.lower = MANT (fl1.l) << 29;
+	
+  return dl.d;
+}
+
+/* convert double to float */
+float
+__truncdfsf2 (double a1)
+{
+  register long exp;
+  register long mant;
+  register union float_long fl;
+  register union double_long dl1;
+
+  dl1.d = a1;
+
+  if (!dl1.l.upper && !dl1.l.lower)
+    return 0;
+
+  exp = EXPD (dl1) - EXCESSD + EXCESS;
+
+  /* shift double mantissa 6 bits so we can round */
+  mant = MANTD (dl1) >> 6;
+
+  /* now round and shift down */
+  mant += 1;
+  mant >>= 1;
+
+  /* did the round overflow? */
+  if (mant & 0xFF000000)
+    {
+      mant >>= 1;
+      exp++;
+    }
+
+  mant &= ~HIDDEN;
+
+  /* pack up and go home */
+  fl.l = PACK (SIGND (dl1), exp, mant);
+  return (fl.f);
+}
+
+/* convert double to int */
+int
+__fixdfsi (double a1)
+{
+  register union double_long dl1;
+  register long exp;
+  register long l;
+
+  dl1.d = a1;
+
+  if (!dl1.l.upper && !dl1.l.lower) 
+    return 0;
+
+  exp = EXPD (dl1) - EXCESSD - 31;
+  l = MANTD (dl1);
+
+  if (exp > 0) 
+    {
+      /* Return largest integer.  */
+      return SIGND (dl1) ? 0x80000000 : 0x7fffffff;
+    }
+
+  /* shift down until exp = 0 or l = 0 */
+  if (exp < 0 && exp > -32 && l) 
+    l >>= -exp;
+
+  return (SIGND (dl1) ? -l : l);
+}
+
+/* convert float to int */
+int
+__fixsfsi (float a1)
+{
+  double foo = a1;
+  return __fixdfsi (foo);
+}
+
+#else /* EXTFLOAT */
+
+/* Primitive extended precision floating point support.
+
+   We assume all numbers are normalized, don't do any rounding, etc.  */
+
+/* Prototypes for the above in case we use them.  */
+double __floatsidf (int);
+float __floatsisf (int);
+double __extendsfdf2 (float);
+float __truncdfsf2 (double);
+int __fixdfsi (double);
+int __fixsfsi (float);
+
+/* convert double to long double */
+long double
+__extenddfxf2 (double d)
+{
+  register union double_long dl;
+  register union long_double_long ldl;
+  register long exp;
+
+  dl.d = d;
+  /*printf ("dfxf in: %g\n", d);*/
+
+  if (!dl.l.upper && !dl.l.lower)
+    return 0;
+
+  ldl.l.upper = SIGND (dl);
+  exp = EXPD (dl) - EXCESSD + EXCESSX;
+  ldl.l.upper |= exp << 16;
+  ldl.l.middle = HIDDENX;
+  /* 31-20: # mantissa bits in ldl.l.middle - # mantissa bits in dl.l.upper */
+  ldl.l.middle |= (dl.l.upper & MANTDMASK) << (31 - 20);
+  /* 1+20: explicit-integer-bit + # mantissa bits in dl.l.upper */
+  ldl.l.middle |= dl.l.lower >> (1 + 20);
+  /* 32 - 21: # bits of dl.l.lower in ldl.l.middle */
+  ldl.l.lower = dl.l.lower << (32 - 21);
+
+  /*printf ("dfxf out: %s\n", dumpxf (ldl.ld));*/
+  return ldl.ld;
+}
+
+/* convert long double to double */
+double
+__truncxfdf2 (long double ld)
+{
+  register long exp;
+  register union double_long dl;
+  register union long_double_long ldl;
+
+  ldl.ld = ld;
+  /*printf ("xfdf in: %s\n", dumpxf (ld));*/
+
+  if (!ldl.l.upper && !ldl.l.middle && !ldl.l.lower)
+    return 0;
+
+  exp = EXPX (ldl) - EXCESSX + EXCESSD;
+  /* ??? quick and dirty: keep `exp' sane */
+  if (exp >= EXPDMASK)
+    exp = EXPDMASK - 1;
+  dl.l.upper = SIGNX (ldl);
+  dl.l.upper |= exp << (32 - (EXPDBITS + 1));
+  /* +1-1: add one for sign bit, but take one off for explicit-integer-bit */
+  dl.l.upper |= (ldl.l.middle & MANTXMASK) >> (EXPDBITS + 1 - 1);
+  dl.l.lower = (ldl.l.middle & MANTXMASK) << (32 - (EXPDBITS + 1 - 1));
+  dl.l.lower |= ldl.l.lower >> (EXPDBITS + 1 - 1);
+
+  /*printf ("xfdf out: %g\n", dl.d);*/
+  return dl.d;
+}
+
+/* convert a float to a long double */
+long double
+__extendsfxf2 (float f)
+{
+  long double foo = __extenddfxf2 (__extendsfdf2 (f));
+  return foo;
+}
+
+/* convert a long double to a float */
+float
+__truncxfsf2 (long double ld)
+{
+  float foo = __truncdfsf2 (__truncxfdf2 (ld));
+  return foo;
+}
+
+/* convert an int to a long double */
+long double
+__floatsixf (int l)
+{
+  double foo = __floatsidf (l);
+  return foo;
+}
+
+/* convert a long double to an int */
+int
+__fixxfsi (long double ld)
+{
+  int foo = __fixdfsi ((double) ld);
+  return foo;
+}
+
+/* The remaining provide crude math support by working in double precision.  */
+
+long double
+__addxf3 (long double x1, long double x2)
+{
+  return (double) x1 + (double) x2;
+}
+
+long double
+__subxf3 (long double x1, long double x2)
+{
+  return (double) x1 - (double) x2;
+}
+
+long double
+__mulxf3 (long double x1, long double x2)
+{
+  return (double) x1 * (double) x2;
+}
+
+long double
+__divxf3 (long double x1, long double x2)
+{
+  return (double) x1 / (double) x2;
+}
+
+long double
+__negxf2 (long double x1)
+{
+  return - (double) x1;
+}
+
+long
+__cmpxf2 (long double x1, long double x2)
+{
+  return __cmpdf2 ((double) x1, (double) x2);
+}
+
+long
+__eqxf2 (long double x1, long double x2)
+{
+  return __cmpdf2 ((double) x1, (double) x2);
+}
+
+long
+__nexf2 (long double x1, long double x2)
+{
+  return __cmpdf2 ((double) x1, (double) x2);
+}
+
+long
+__ltxf2 (long double x1, long double x2)
+{
+  return __cmpdf2 ((double) x1, (double) x2);
+}
+
+long
+__lexf2 (long double x1, long double x2)
+{
+  return __cmpdf2 ((double) x1, (double) x2);
+}
+
+long
+__gtxf2 (long double x1, long double x2)
+{
+  return __cmpdf2 ((double) x1, (double) x2);
+}
+
+long
+__gexf2 (long double x1, long double x2)
+{
+  return __cmpdf2 ((double) x1, (double) x2);
+}
+
+#endif /* EXTFLOAT */
diff --git a/gcc/config/m68k/lb1sf68.asm b/gcc/config/m68k/lb1sf68.asm
new file mode 100644
index 00000000000..6e9124d165e
--- /dev/null
+++ b/gcc/config/m68k/lb1sf68.asm
@@ -0,0 +1,3191 @@
+/* libgcc1 routines for 68000 w/o floating-point hardware. */
+/* Copyright (C) 1994 Free Software Foundation, Inc.
+
+This file is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file with other programs, and to distribute
+those programs without any restriction coming from the use of this
+file.  (The General Public License restrictions do apply in other
+respects; for example, they cover modification of the file, and
+distribution when not linked into another program.)
+
+This file is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; see the file COPYING.  If not, write to
+the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
+
+/* As a special exception, if you link this library with files
+   compiled with GCC to produce an executable, this does not cause
+   the resulting executable to be covered by the GNU General Public License.
+   This exception does not however invalidate any other reasons why
+   the executable file might be covered by the GNU General Public License.  */
+
+/* Use this one for any 680x0; assumes no floating point hardware.
+   The trailing " '" appearing on some lines is for ANSI preprocessors.  Yuk.
+   Some of this code comes from MINIX, via the folks at ericsson.
+   D. V. Henkel-Wallace (gumby@cygnus.com) Fete Bastille, 1992
+*/
+
+/* These are predefined by new versions of GNU cpp.  */
+
+#ifndef __USER_LABEL_PREFIX__
+#define __USER_LABEL_PREFIX__ _
+#endif
+
+#ifndef __REGISTER_PREFIX__
+#define __REGISTER_PREFIX__
+#endif
+
+/* ANSI concatenation macros.  */
+
+#define CONCAT1(a, b) CONCAT2(a, b)
+#define CONCAT2(a, b) a ## b
+
+/* Use the right prefix for global labels.  */
+
+#define SYM(x) CONCAT1 (__USER_LABEL_PREFIX__, x)
+
+/* Use the right prefix for registers.  */
+
+#define REG(x) CONCAT1 (__REGISTER_PREFIX__, x)
+
+#define d0 REG (d0)
+#define d1 REG (d1)
+#define d2 REG (d2)
+#define d3 REG (d3)
+#define d4 REG (d4)
+#define d5 REG (d5)
+#define d6 REG (d6)
+#define d7 REG (d7)
+#define a0 REG (a0)
+#define a1 REG (a1)
+#define a2 REG (a2)
+#define a3 REG (a3)
+#define a4 REG (a4)
+#define a5 REG (a5)
+#define a6 REG (a6)
+#define fp REG (fp)
+#define sp REG (sp)
+
+#ifdef L_floatex
+
+| This is an attempt at a decent floating point (single, double and 
+| extended double) code for the GNU C compiler. It should be easy to
+| adapt to other compilers (but beware of the local labels!).
+
+| Starting date: 21 October, 1990
+
+| It is convenient to introduce the notation (s,e,f) for a floating point
+| number, where s=sign, e=exponent, f=fraction. We will call a floating
+| point number fpn to abbreviate, independently of the precision.
+| Let MAX_EXP be in each case the maximum exponent (255 for floats, 1023 
+| for doubles and 16383 for long doubles). We then have the following 
+| different cases:
+|  1. Normalized fpns have 0 < e < MAX_EXP. They correspond to 
+|     (-1)^s x 1.f x 2^(e-bias-1).
+|  2. Denormalized fpns have e=0. They correspond to numbers of the form
+|     (-1)^s x 0.f x 2^(-bias).
+|  3. +/-INFINITY have e=MAX_EXP, f=0.
+|  4. Quiet NaN (Not a Number) have all bits set.
+|  5. Signaling NaN (Not a Number) have s=0, e=MAX_EXP, f=1.
+
+|=============================================================================
+|                                  exceptions
+|=============================================================================
+
+| This is the floating point condition code register (_fpCCR):
+|
+| struct {
+|   short _exception_bits;	
+|   short _trap_enable_bits;	
+|   short _sticky_bits;
+|   short _rounding_mode;
+|   short _format;
+|   short _last_operation;
+|   union {
+|     float sf;
+|     double df;
+|   } _operand1;
+|   union {
+|     float sf;
+|     double df;
+|   } _operand2;
+| } _fpCCR;
+
+	.data
+	.even
+
+	.globl	SYM (_fpCCR)
+	
+SYM (_fpCCR):
+__exception_bits:
+	.word	0
+__trap_enable_bits:
+	.word	0
+__sticky_bits:
+	.word	0
+__rounding_mode:
+	.word	ROUND_TO_NEAREST
+__format:
+	.word	NIL
+__last_operation:
+	.word	NOOP
+__operand1:
+	.long	0
+	.long	0
+__operand2:
+	.long 	0
+	.long	0
+
+| Offsets:
+EBITS  = __exception_bits - SYM (_fpCCR)
+TRAPE  = __trap_enable_bits - SYM (_fpCCR)
+STICK  = __sticky_bits - SYM (_fpCCR)
+ROUND  = __rounding_mode - SYM (_fpCCR)
+FORMT  = __format - SYM (_fpCCR)
+LASTO  = __last_operation - SYM (_fpCCR)
+OPER1  = __operand1 - SYM (_fpCCR)
+OPER2  = __operand2 - SYM (_fpCCR)
+
+| The following exception types are supported:
+INEXACT_RESULT 		= 0x0001
+UNDERFLOW 		= 0x0002
+OVERFLOW 		= 0x0004
+DIVIDE_BY_ZERO 		= 0x0008
+INVALID_OPERATION 	= 0x0010
+
+| The allowed rounding modes are:
+UNKNOWN           = -1
+ROUND_TO_NEAREST  = 0 | round result to nearest representable value
+ROUND_TO_ZERO     = 1 | round result towards zero
+ROUND_TO_PLUS     = 2 | round result towards plus infinity
+ROUND_TO_MINUS    = 3 | round result towards minus infinity
+
+| The allowed values of format are:
+NIL          = 0
+SINGLE_FLOAT = 1
+DOUBLE_FLOAT = 2
+LONG_FLOAT   = 3
+
+| The allowed values for the last operation are:
+NOOP         = 0
+ADD          = 1
+MULTIPLY     = 2
+DIVIDE       = 3
+NEGATE       = 4
+COMPARE      = 5
+EXTENDSFDF   = 6
+TRUNCDFSF    = 7
+
+|=============================================================================
+|                           __clear_sticky_bits
+|=============================================================================
+
+| The sticky bits are normally not cleared (thus the name), whereas the 
+| exception type and exception value reflect the last computation. 
+| This routine is provided to clear them (you can also write to _fpCCR,
+| since it is globally visible).
+
+	.globl  SYM (__clear_sticky_bit)
+
+	.text
+	.even
+
+| void __clear_sticky_bits(void);
+SYM (__clear_sticky_bit):		
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@(STICK)
+	rts
+
+|=============================================================================
+|                           $_exception_handler
+|=============================================================================
+
+	.globl  $_exception_handler
+
+	.text
+	.even
+
+| This is the common exit point if an exception occurs.
+| NOTE: it is NOT callable from C!
+| It expects the exception type in d7, the format (SINGLE_FLOAT,
+| DOUBLE_FLOAT or LONG_FLOAT) in d6, and the last operation code in d5.
+| It sets the corresponding exception and sticky bits, and the format. 
+| Depending on the format if fills the corresponding slots for the 
+| operands which produced the exception (all this information is provided
+| so if you write your own exception handlers you have enough information
+| to deal with the problem).
+| Then checks to see if the corresponding exception is trap-enabled, 
+| in which case it pushes the address of _fpCCR and traps through 
+| trap FPTRAP (15 for the moment).
+
+FPTRAP = 15
+
+$_exception_handler:
+	lea	SYM (_fpCCR),a0
+	movew	d7,a0@(EBITS)	| set __exception_bits
+	orw	d7,a0@(STICK)	| and __sticky_bits
+	movew	d6,a0@(FORMT)	| and __format
+	movew	d5,a0@(LASTO)	| and __last_operation
+
+| Now put the operands in place:
+	cmpw	#SINGLE_FLOAT,d6
+	beq	1f
+	movel	a6@(8),a0@(OPER1)
+	movel	a6@(12),a0@(OPER1+4)
+	movel	a6@(16),a0@(OPER2)
+	movel	a6@(20),a0@(OPER2+4)
+	bra	2f
+1:	movel	a6@(8),a0@(OPER1)
+	movel	a6@(12),a0@(OPER2)
+2:
+| And check whether the exception is trap-enabled:
+	andw	a0@(TRAPE),d7	| is exception trap-enabled?
+	beq	1f		| no, exit
+	pea	SYM (_fpCCR)	| yes, push address of _fpCCR
+	trap	#FPTRAP		| and trap
+1:	moveml	sp@+,d2-d7	| restore data registers
+	unlk	a6		| and return
+	rts
+#endif /* L_floatex */
+
+#ifdef  L_mulsi3
+	.text
+	.proc
+|#PROC# 04
+	.globl	SYM (__mulsi3)
+SYM (__mulsi3):
+|#PROLOGUE# 0
+	link	a6,#0
+	addl	#-LF14,sp
+	moveml	#LS14,sp@
+|#PROLOGUE# 1
+	movew	a6@(0x8), d0	/* x0 -> d0 */
+	muluw	a6@(0xe), d0	/* x0*y1 */
+	movew	a6@(0xa), d1	/* x1 -> d1 */
+	muluw	a6@(0xc), d1	/* x1*y0 */
+	addw	d1, d0
+	lsll	#8, d0
+	lsll	#8, d0
+	movew	a6@(0xa), d1	/* x1 -> d1 */
+	muluw	a6@(0xe), d1	/* x1*y1 */
+	addl	d1, d0
+	jra	LE14
+LE14:
+|#PROLOGUE# 2
+	moveml	sp@, #LS14
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+	LF14 = 4
+	LS14 = 0x0002		/* d1 will be saved and restored */
+	LFF14 = 0
+	LSS14 = 0x0
+	LV14 = 0
+#endif /* L_mulsi3 */
+
+#ifdef  L_udivsi3
+	.text
+	.proc
+|#PROC# 04
+	.globl	SYM (__udivsi3)
+SYM (__udivsi3):
+|#PROLOGUE# 0
+	link	a6,#0
+	addl	#-LF14,sp
+	moveml	#LS14,sp@
+|#PROLOGUE# 1
+	movel	a6@(0xc), d0	/* d0 = divisor */
+	movel	a6@(0x8), d1	/* d1 = dividend */
+	movel	d1, d3
+
+
+	cmpl	#0x10000, d0	/* divisor >= 2 ^ 16 ?   */
+	bge	l4		/* then try next algorithm */
+	movel	d1, d2
+	lsrl	#8, d2		/* get high dividend */
+	lsrl	#8, d2
+	divu	d0, d2          /* high quotient in lower word */
+	movew	d2, d1		/* save high quotient */
+	swap	d1
+	movew	d3, d2		/* get low dividend + high rest */
+	divu	d0, d2		/* low quotient */
+	movew	d2, d1
+	jra	l5
+
+l4:	movel	d0, d2		/* use d2 as divisor backup */
+l4a:	lsrl	#1, d0		/* shift divisor */
+	lsrl	#1, d1		/* shift dividend */
+	cmpl	#0x10000, d0	/* still divisor >= 2 ^ 16 ?  */
+	bge	l4a
+	divu	d0, d1		/* now we have 16bit divisor => compute remainder */
+	andl	#0xffff, d1
+	movel	d1, sp@-	/* multiply divisor with */
+	movel	d2, sp@-	/* remainder             */
+	jbsr	SYM (__umulsi3)	/* and                   */
+	addql	#8, sp
+	cmpl	d0, d3		/* compare the result with the dividend */
+	bge	l5		/* if dividend >= result => nofix */
+	subql	#1, d1
+
+l5:	movel	d1, d0	
+
+l6:	jra	LE14
+LE14:
+|#PROLOGUE# 2
+	moveml	sp@, #LS14
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+	LF14 = 16
+	LS14 = 0x000e		/* d1-d3 will be saved and restored */
+	LFF14 = 0
+	LSS14 = 0x0
+	LV14 = 0
+#endif /* L_udivsi3 */
+
+#ifdef L_umulsi3
+	.text
+	.proc
+|#PROC# 04
+	.globl	SYM (__umulsi3)
+SYM (__umulsi3):
+|#PROLOGUE# 0
+	link	a6,#0
+	addl	#-LF14,sp
+	moveml	#LS14,sp@
+|#PROLOGUE# 1
+	movew	a6@(0x8), d0	/* x0 -> d0 */
+	muluw	a6@(0xe), d0	/* x0*y1 */
+	movew	a6@(0xa), d1	/* x1 -> d1 */
+	muluw	a6@(0xc), d1	/* x1*y0 */
+	addw	d1, d0
+	lsll	#8, d0
+	lsll	#8, d0
+	movew	a6@(0xa), d1	/* x1 -> d1 */
+	muluw	a6@(0xe), d1	/* x1*y1 */
+	addl	d1, d0
+	jra	LE15
+LE15:
+|#PROLOGUE# 2
+	moveml	sp@, #LS14
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+	LF14 = 4
+	LS14 = 0x0002		/* d1 will be saved and restored */
+	LFF14 = 0
+	LSS14 = 0x0
+	LV14 = 0
+#endif /* L_umulsi3 */
+
+#ifdef  L_divsi3
+	.text
+	.proc
+|#PROC# 04
+	.globl	SYM (__divsi3)
+SYM (__divsi3):
+|#PROLOGUE# 0
+	link	a6,#0
+	addl	#-LF14,sp
+	moveml	#LS14,sp@
+|#PROLOGUE# 1
+	moveb	#1, d4		/* sign of result stored in d4 (=1 or =-1) */
+	movel	a6@(0xc), d0	/* d0 = divisor */
+	bpl	l1
+	negl	d0
+	negb	d4		/* change sign because divisor <0  */
+l1:	movel	a6@(0x8), d1	/* d1 = dividend */
+	bpl	l2
+	negl	d1
+	negb	d4
+l2:	movel	d1, d3
+
+
+	cmpl	#0x10000, d0	/* divisor >= 2 ^ 16 ?   */
+	bge	l4		/* then try next algorithm */
+	movel	d1, d2
+	lsrl	#8, d2		/* get high dividend */
+	lsrl	#8, d2
+	divu	d0, d2          /* high quotient in lower word */
+	movew	d2, d1		/* save high quotient */
+	swap	d1
+	movew	d3, d2		/* get low dividend + high rest */
+	divu	d0, d2		/* low quotient */
+	movew	d2, d1
+	jra	l5
+
+l4:	movel	d0, d2		/* use d2 as divisor backup */
+l4a:	lsrl	#1, d0		/* shift divisor */
+	lsrl	#1, d1		/* shift dividend */
+	cmpl	#0x10000, d0	/* still divisor >= 2 ^ 16 ?  */
+	bge	l4a
+	divu	d0, d1		/* now we have 16bit divisor => compute remainder */
+	andl	#0xffff, d1
+	movel	d1, sp@-	/* multiply divisor with */
+	movel	d2, sp@-	/* remainder             */
+	jbsr	SYM (__umulsi3)	/* and                   */
+	addql	#8, sp
+	cmpl	d0, d3		/* compare the result with the dividend */
+	bge	l5		/* if dividend >= result => nofix */
+	subql	#1, d1
+
+l5:	movel	d1, d0	
+	tstb	d4
+	bpl	l6
+	negl	d0
+
+l6:	jra	LE14
+LE14:
+|#PROLOGUE# 2
+	moveml	sp@, #LS14
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+	LF14 = 16
+	LS14 = 0x001e		/* d1-d4 will be saved and restored */
+	LFF14 = 8
+	LSS14 = 0x0
+	LV14 = 8
+#endif /* L_divsi3 */
+
+#ifdef  L_umodsi3
+	.text
+	.proc
+|#PROC# 04
+	.globl	SYM (__umodsi3)
+SYM (__umodsi3):
+|#PROLOGUE# 0
+	link	a6,#0
+	addl	#-LF14,sp
+	moveml	#LS14,sp@
+|#PROLOGUE# 1
+	movel	a6@(0xc),d1	/* divisor */
+	movel	a6@(0x8),d2	/* dividend */
+	movel	d1, sp@-
+	movel	d2, sp@-
+	jbsr	SYM (__udivsi3)	/* d0 = a/b */
+	addql	#8, sp
+	movel	d0, sp@-
+	movel	d1, sp@-
+	jbsr	SYM (__umulsi3)	/* d0 = (a/b)*b */
+	addql	#8, sp
+	negl	d0
+	addl	d2, d0		/* d0 = a - (a/b)*b */
+	jra	LE14
+LE14:
+|#PROLOGUE# 2
+	moveml	sp@, #LS14
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+	LF14 = 8
+	LS14 = 0x006		/* d1-d2 will be saved and restored */
+	LFF14 = 0
+	LSS14 = 0x0
+	LV14 = 0
+#endif /* L_umodsi3 */
+
+#ifdef  L_modsi3
+	.text
+	.proc
+|#PROC# 04
+	.globl	SYM (__modsi3)
+SYM (__modsi3):
+|#PROLOGUE# 0
+	link	a6,#0
+	addl	#-LF14,sp
+	moveml	#LS14,sp@
+|#PROLOGUE# 1
+	movel	a6@(0xc),d1	/* divisor */
+	movel	a6@(0x8),d2	/* dividend */
+	movel	d1, sp@-
+	movel	d2, sp@-
+	jbsr	SYM (__divsi3)	/* d0 = a/b */
+	addql	#8, sp
+	movel	d0, sp@-
+	movel	d1, sp@-
+	jbsr	SYM (__mulsi3)	/* d0 = (a/b)*b */
+	addql	#8, sp
+	negl	d0
+	addl	d2, d0		/* d0 = a - (a/b)*b */
+	jra	LE14
+LE14:
+|#PROLOGUE# 2
+	moveml	sp@, #LS14
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+	LF14 = 8
+	LS14 = 0x006		/* d1-d2 will be saved and restored */
+	LFF14 = 0
+	LSS14 = 0x0
+	LV14 = 0
+#endif /* L_modsi3 */
+
+#ifdef  L_lshrsi3
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	4
+	LS18	=	128
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__lshrsi3)
+SYM (__lshrsi3):
+|#PROLOGUE# 0
+	link	a6,#-4
+|#PROLOGUE# 1
+	movl	a6@(8),d0
+	movw	a6@(14),d1
+	lsrl	d1,d0
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_lshrsi3 */
+
+#ifdef  L_lshlsi3
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	4
+	LS18	=	128
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__lshlsi3)
+SYM (__lshlsi3):
+|#PROLOGUE# 0
+	link	a6,#-4
+|#PROLOGUE# 1
+	movl	a6@(8),d0
+	movw	a6@(14),d1
+	lsll	d1,d0
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_lshlsi3 */
+
+#ifdef  L_ashrsi3
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	4
+	LS18	=	128
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__ashrsi3)
+SYM (__ashrsi3):
+|#PROLOGUE# 0
+	link	a6,#-4
+|#PROLOGUE# 1
+	movl	a6@(8),d0
+	movw	a6@(14),d1
+	asrl	d1,d0
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_ashrsi3 */
+
+#ifdef  L_ashlsi3
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	4
+	LS18	=	128
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__ashlsi3)
+SYM (__ashlsi3):
+|#PROLOGUE# 0
+	link	a6,#-4
+|#PROLOGUE# 1
+	movl	a6@(8),d0
+	movw	a6@(14),d1
+	asll	d1,d0
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_ashlsi3 */
+
+#ifdef  L_double
+
+	.globl	SYM (_fpCCR)
+	.globl  $_exception_handler
+
+QUIET_NaN      = 0xffffffff
+
+D_MAX_EXP      = 0x07ff
+D_BIAS         = 1022
+DBL_MAX_EXP    = D_MAX_EXP - D_BIAS
+DBL_MIN_EXP    = 1 - D_BIAS
+DBL_MANT_DIG   = 53
+
+INEXACT_RESULT 		= 0x0001
+UNDERFLOW 		= 0x0002
+OVERFLOW 		= 0x0004
+DIVIDE_BY_ZERO 		= 0x0008
+INVALID_OPERATION 	= 0x0010
+
+DOUBLE_FLOAT = 2
+
+NOOP         = 0
+ADD          = 1
+MULTIPLY     = 2
+DIVIDE       = 3
+NEGATE       = 4
+COMPARE      = 5
+EXTENDSFDF   = 6
+TRUNCDFSF    = 7
+
+UNKNOWN           = -1
+ROUND_TO_NEAREST  = 0 | round result to nearest representable value
+ROUND_TO_ZERO     = 1 | round result towards zero
+ROUND_TO_PLUS     = 2 | round result towards plus infinity
+ROUND_TO_MINUS    = 3 | round result towards minus infinity
+
+| Entry points:
+
+	.globl SYM (__adddf3)
+	.globl SYM (__subdf3)
+	.globl SYM (__muldf3)
+	.globl SYM (__divdf3)
+	.globl SYM (__negdf2)
+	.globl SYM (__cmpdf2)
+
+	.text
+	.even
+
+| These are common routines to return and signal exceptions.	
+
+Ld$den:
+| Return and signal a denormalized number
+	orl	d7,d0
+	movew	#UNDERFLOW,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#DOUBLE_FLOAT,d6
+	jmp	$_exception_handler
+
+Ld$infty:
+Ld$overflow:
+| Return a properly signed INFINITY and set the exception flags 
+	movel	#0x7ff00000,d0
+	movel	#0,d1
+	orl	d7,d0
+	movew	#OVERFLOW,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#DOUBLE_FLOAT,d6
+	jmp	$_exception_handler
+
+Ld$underflow:
+| Return 0 and set the exception flags 
+	movel	#0,d0
+	movel	d0,d1
+	movew	#UNDERFLOW,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#DOUBLE_FLOAT,d6
+	jmp	$_exception_handler
+
+Ld$inop:
+| Return a quiet NaN and set the exception flags
+	movel	#QUIET_NaN,d0
+	movel	d0,d1
+	movew	#INVALID_OPERATION,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#DOUBLE_FLOAT,d6
+	jmp	$_exception_handler
+
+Ld$div$0:
+| Return a properly signed INFINITY and set the exception flags
+	movel	#0x7ff00000,d0
+	movel	#0,d1
+	orl	d7,d0
+	movew	#DIVIDE_BY_ZERO,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#DOUBLE_FLOAT,d6
+	jmp	$_exception_handler
+
+|=============================================================================
+|=============================================================================
+|                         double precision routines
+|=============================================================================
+|=============================================================================
+
+| A double precision floating point number (double) has the format:
+|
+| struct _double {
+|  unsigned int sign      : 1;  /* sign bit */ 
+|  unsigned int exponent  : 11; /* exponent, shifted by 126 */
+|  unsigned int fraction  : 52; /* fraction */
+| } double;
+| 
+| Thus sizeof(double) = 8 (64 bits). 
+|
+| All the routines are callable from C programs, and return the result 
+| in the register pair d0-d1. They also preserve all registers except 
+| d0-d1 and a0-a1.
+
+|=============================================================================
+|                              __subdf3
+|=============================================================================
+
+| double __subdf3(double, double);
+SYM (__subdf3):
+	bchg	#31,sp@(12)	| change sign of second operand
+				| and fall through, so we always add
+|=============================================================================
+|                              __adddf3
+|=============================================================================
+
+| double __adddf3(double, double);
+SYM (__adddf3):
+	link	a6,#0		| everything will be done in registers
+	moveml	d2-d7,sp@-	| save all data registers and a2 (but d0-d1)
+	movel	a6@(8),d0	| get first operand
+	movel	a6@(12),d1	| 
+	movel	a6@(16),d2	| get second operand
+	movel	a6@(20),d3	| 
+
+	movel	d0,d7		| get d0's sign bit in d7 '
+	addl	d1,d1		| check and clear sign bit of a, and gain one
+	addxl	d0,d0		| bit of extra precision
+	beq	Ladddf$b	| if zero return second operand
+
+	movel	d2,d6		| save sign in d6 
+	addl	d3,d3		| get rid of sign bit and gain one bit of
+	addxl	d2,d2		| extra precision
+	beq	Ladddf$a	| if zero return first operand
+
+	andl	#0x80000000,d7	| isolate a's sign bit '
+        swap	d6		| and also b's sign bit '
+	andw	#0x8000,d6	|
+	orw	d6,d7		| and combine them into d7, so that a's sign '
+				| bit is in the high word and b's is in the '
+				| low word, so d6 is free to be used
+	movel	d7,a0		| now save d7 into a0, so d7 is free to
+                		| be used also
+
+| Get the exponents and check for denormalized and/or infinity.
+
+	movel	#0x001fffff,d6	| mask for the fraction
+	movel	#0x00200000,d7	| mask to put hidden bit back
+
+	movel	d0,d4		| 
+	andl	d6,d0		| get fraction in d0
+	notl	d6		| make d6 into mask for the exponent
+	andl	d6,d4		| get exponent in d4
+	beq	Ladddf$a$den	| branch if a is denormalized
+	cmpl	d6,d4		| check for INFINITY or NaN
+	beq	Ladddf$nf       | 
+	orl	d7,d0		| and put hidden bit back
+Ladddf$1:
+	swap	d4		| shift right exponent so that it starts
+	lsrw	#5,d4		| in bit 0 and not bit 20
+| Now we have a's exponent in d4 and fraction in d0-d1 '
+	movel	d2,d5		| save b to get exponent
+	andl	d6,d5		| get exponent in d5
+	beq	Ladddf$b$den	| branch if b is denormalized
+	cmpl	d6,d5		| check for INFINITY or NaN
+	beq	Ladddf$nf
+	notl	d6		| make d6 into mask for the fraction again
+	andl	d6,d2		| and get fraction in d2
+	orl	d7,d2		| and put hidden bit back
+Ladddf$2:
+	swap	d5		| shift right exponent so that it starts
+	lsrw	#5,d5		| in bit 0 and not bit 20
+
+| Now we have b's exponent in d5 and fraction in d2-d3. '
+
+| The situation now is as follows: the signs are combined in a0, the 
+| numbers are in d0-d1 (a) and d2-d3 (b), and the exponents in d4 (a)
+| and d5 (b). To do the rounding correctly we need to keep all the
+| bits until the end, so we need to use d0-d1-d2-d3 for the first number
+| and d4-d5-d6-d7 for the second. To do this we store (temporarily) the
+| exponents in a2-a3.
+
+	moveml	a2-a3,sp@-	| save the address registers
+
+	movel	d4,a2		| save the exponents
+	movel	d5,a3		| 
+
+	movel	#0,d7		| and move the numbers around
+	movel	d7,d6		|
+	movel	d3,d5		|
+	movel	d2,d4		|
+	movel	d7,d3		|
+	movel	d7,d2		|
+
+| Here we shift the numbers until the exponents are the same, and put 
+| the largest exponent in a2.
+	exg	d4,a2		| get exponents back
+	exg	d5,a3		|
+	cmpw	d4,d5		| compare the exponents
+	beq	Ladddf$3	| if equal don't shift '
+	bhi	9f		| branch if second exponent is higher
+
+| Here we have a's exponent larger than b's, so we have to shift b. We do 
+| this by using as counter d2:
+1:	movew	d4,d2		| move largest exponent to d2
+	subw	d5,d2		| and substract second exponent
+	exg	d4,a2		| get back the longs we saved
+	exg	d5,a3		|
+| if difference is too large we don't shift (actually, we can just exit) '
+	cmpw	#DBL_MANT_DIG+2,d2
+	bge	Ladddf$b$small
+	cmpw	#32,d2		| if difference >= 32, shift by longs
+	bge	5f
+2:	cmpw	#16,d2		| if difference >= 16, shift by words	
+	bge	6f
+	bra	3f		| enter dbra loop
+
+4:	lsrl	#1,d4
+	roxrl	#1,d5
+	roxrl	#1,d6
+	roxrl	#1,d7
+3:	dbra	d2,4b
+	movel	#0,d2
+	movel	d2,d3	
+	bra	Ladddf$4
+5:
+	movel	d6,d7
+	movel	d5,d6
+	movel	d4,d5
+	movel	#0,d4
+	subw	#32,d2
+	bra	2b
+6:
+	movew	d6,d7
+	swap	d7
+	movew	d5,d6
+	swap	d6
+	movew	d4,d5
+	swap	d5
+	movew	#0,d4
+	swap	d4
+	subw	#16,d2
+	bra	3b
+	
+9:	exg	d4,d5
+	movew	d4,d6
+	subw	d5,d6		| keep d5 (largest exponent) in d4
+	exg	d4,a2
+	exg	d5,a3
+| if difference is too large we don't shift (actually, we can just exit) '
+	cmpw	#DBL_MANT_DIG+2,d6
+	bge	Ladddf$a$small
+	cmpw	#32,d6		| if difference >= 32, shift by longs
+	bge	5f
+2:	cmpw	#16,d6		| if difference >= 16, shift by words	
+	bge	6f
+	bra	3f		| enter dbra loop
+
+4:	lsrl	#1,d0
+	roxrl	#1,d1
+	roxrl	#1,d2
+	roxrl	#1,d3
+3:	dbra	d6,4b
+	movel	#0,d7
+	movel	d7,d6
+	bra	Ladddf$4
+5:
+	movel	d2,d3
+	movel	d1,d2
+	movel	d0,d1
+	movel	#0,d0
+	subw	#32,d6
+	bra	2b
+6:
+	movew	d2,d3
+	swap	d3
+	movew	d1,d2
+	swap	d2
+	movew	d0,d1
+	swap	d1
+	movew	#0,d0
+	swap	d0
+	subw	#16,d6
+	bra	3b
+Ladddf$3:
+	exg	d4,a2	
+	exg	d5,a3
+Ladddf$4:	
+| Now we have the numbers in d0--d3 and d4--d7, the exponent in a2, and
+| the signs in a4.
+
+| Here we have to decide whether to add or substract the numbers:
+	exg	d7,a0		| get the signs 
+	exg	d6,a3		| a3 is free to be used
+	movel	d7,d6		|
+	movew	#0,d7	        | get a's sign in d7 '
+	swap	d6              |
+	movew	#0,d6           | and b's sign in d6 '
+	eorl	d7,d6		| compare the signs
+	bmi	Lsubdf$0	| if the signs are different we have 
+				| to substract
+	exg	d7,a0		| else we add the numbers
+	exg	d6,a3		|
+	addl	d7,d3		|
+	addxl	d6,d2		|
+	addxl	d5,d1		| 
+	addxl	d4,d0           |
+
+	movel	a2,d4		| return exponent to d4
+	movel	a0,d7		| 
+	andl	#0x80000000,d7	| d7 now has the sign
+
+	moveml	sp@+,a2-a3	
+
+| Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
+| the case of denormalized numbers in the rounding routine itself).
+| As in the addition (not in the substraction!) we could have set 
+| one more bit we check this:
+	btst	#DBL_MANT_DIG+1,d0	
+	beq	1f
+	lsrl	#1,d0
+	roxrl	#1,d1
+	roxrl	#1,d2
+	roxrl	#1,d3
+	addw	#1,d4
+1:
+	lea	Ladddf$5,a0	| to return from rounding routine
+	lea	SYM (_fpCCR),a1	| check the rounding mode
+	movew	a1@(6),d6	| rounding mode in d6
+	beq	Lround$to$nearest
+	cmpw	#ROUND_TO_PLUS,d6
+	bhi	Lround$to$minus
+	blt	Lround$to$zero
+	bra	Lround$to$plus
+Ladddf$5:
+| Put back the exponent and check for overflow
+	cmpw	#0x7ff,d4	| is the exponent big?
+	bge	1f
+	bclr	#DBL_MANT_DIG-1,d0
+	lslw	#4,d4		| put exponent back into position
+	swap	d0		| 
+	orw	d4,d0		|
+	swap	d0		|
+	bra	Ladddf$ret
+1:
+	movew	#ADD,d5
+	bra	Ld$overflow
+
+Lsubdf$0:
+| Here we do the substraction.
+	exg	d7,a0		| put sign back in a0
+	exg	d6,a3		|
+	subl	d7,d3		|
+	subxl	d6,d2		|
+	subxl	d5,d1		|
+	subxl	d4,d0		|
+	beq	Ladddf$ret$1	| if zero just exit
+	bpl	1f		| if positive skip the following
+	exg	d7,a0		|
+	bchg	#31,d7		| change sign bit in d7
+	exg	d7,a0		|
+	negl	d3		|
+	negxl	d2		|
+	negxl	d1              | and negate result
+	negxl	d0              |
+1:	
+	movel	a2,d4		| return exponent to d4
+	movel	a0,d7
+	andl	#0x80000000,d7	| isolate sign bit
+	moveml	sp@+,a2-a3	|
+
+| Before rounding normalize so bit #DBL_MANT_DIG is set (we will consider
+| the case of denormalized numbers in the rounding routine itself).
+| As in the addition (not in the substraction!) we could have set 
+| one more bit we check this:
+	btst	#DBL_MANT_DIG+1,d0	
+	beq	1f
+	lsrl	#1,d0
+	roxrl	#1,d1
+	roxrl	#1,d2
+	roxrl	#1,d3
+	addw	#1,d4
+1:
+	lea	Lsubdf$1,a0	| to return from rounding routine
+	lea	SYM (_fpCCR),a1	| check the rounding mode
+	movew	a1@(6),d6	| rounding mode in d6
+	beq	Lround$to$nearest
+	cmpw	#ROUND_TO_PLUS,d6
+	bhi	Lround$to$minus
+	blt	Lround$to$zero
+	bra	Lround$to$plus
+Lsubdf$1:
+| Put back the exponent and sign (we don't have overflow). '
+	bclr	#DBL_MANT_DIG-1,d0	
+	lslw	#4,d4		| put exponent back into position
+	swap	d0		| 
+	orw	d4,d0		|
+	swap	d0		|
+	bra	Ladddf$ret
+
+| If one of the numbers was too small (difference of exponents >= 
+| DBL_MANT_DIG+1) we return the other (and now we don't have to '
+| check for finiteness or zero).
+Ladddf$a$small:
+	moveml	sp@+,a2-a3	
+	movel	a6@(16),d0
+	movel	a6@(20),d1
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7	| restore data registers
+	unlk	a6		| and return
+	rts
+
+Ladddf$b$small:
+	moveml	sp@+,a2-a3	
+	movel	a6@(8),d0
+	movel	a6@(12),d1
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7	| restore data registers
+	unlk	a6		| and return
+	rts
+
+Ladddf$a$den:
+	movel	d7,d4		| d7 contains 0x00200000
+	bra	Ladddf$1
+
+Ladddf$b$den:
+	movel	d7,d5           | d7 contains 0x00200000
+	notl	d6
+	bra	Ladddf$2
+
+Ladddf$b:
+| Return b (if a is zero)
+	movel	d2,d0
+	movel	d3,d1
+	bra	1f
+Ladddf$a:
+	movel	a6@(8),d0
+	movel	a6@(12),d1
+1:
+	movew	#ADD,d5
+| Check for NaN and +/-INFINITY.
+	movel	d0,d7         	|
+	andl	#0x80000000,d7  |
+	bclr	#31,d0		|
+	cmpl	#0x7ff00000,d0	|
+	bge	2f		|
+	movel	d0,d0           | check for zero, since we don't  '
+	bne	Ladddf$ret	| want to return -0 by mistake
+	bclr	#31,d7		|
+	bra	Ladddf$ret	|
+2:
+	andl	#0x000fffff,d0	| check for NaN (nonzero fraction)
+	orl	d1,d0		|
+	bne	Ld$inop         |
+	bra	Ld$infty	|
+	
+Ladddf$ret$1:
+	moveml	sp@+,a2-a3	| restore regs and exit
+
+Ladddf$ret:
+| Normal exit.
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	orl	d7,d0		| put sign bit back
+	moveml	sp@+,d2-d7
+	unlk	a6
+	rts
+
+Ladddf$ret$den:
+| Return a denormalized number.
+	lsrl	#1,d0		| shift right once more
+	roxrl	#1,d1           |
+	bra	Ladddf$ret
+
+Ladddf$nf:
+	movew	#ADD,d5
+| This could be faster but it is not worth the effort, since it is not
+| executed very often. We sacrifice speed for clarity here.
+	movel	a6@(8),d0	| get the numbers back (remember that we
+	movel	a6@(12),d1	| did some processing already)
+	movel	a6@(16),d2	| 
+	movel	a6@(20),d3	| 
+	movel	#0x7ff00000,d4	| useful constant (INFINITY)
+	movel	d0,d7		| save sign bits
+	movel	d2,d6		| 
+	bclr	#31,d0		| clear sign bits
+	bclr	#31,d2		| 
+| We know that one of them is either NaN of +/-INFINITY
+| Check for NaN (if either one is NaN return NaN)
+	cmpl	d4,d0		| check first a (d0)
+	bhi	Ld$inop		| if d0 > 0x7ff00000 or equal and
+	bne	2f
+	tstl	d1		| d1 > 0, a is NaN
+	bne	Ld$inop		| 
+2:	cmpl	d4,d2		| check now b (d1)
+	bhi	Ld$inop		| 
+	bne	3f
+	tstl	d3		| 
+	bne	Ld$inop		| 
+3:
+| Now comes the check for +/-INFINITY. We know that both are (maybe not
+| finite) numbers, but we have to check if both are infinite whether we
+| are adding or substracting them.
+	eorl	d7,d6		| to check sign bits
+	bmi	1f
+	andl	#0x80000000,d7	| get (common) sign bit
+	bra	Ld$infty
+1:
+| We know one (or both) are infinite, so we test for equality between the
+| two numbers (if they are equal they have to be infinite both, so we
+| return NaN).
+	cmpl	d2,d0		| are both infinite?
+	bne	1f		| if d0 <> d2 they are not equal
+	cmpl	d3,d1		| if d0 == d2 test d3 and d1
+	beq	Ld$inop		| if equal return NaN
+1:	
+	andl	#0x80000000,d7	| get a's sign bit '
+	cmpl	d4,d0		| test now for infinity
+	beq	Ld$infty	| if a is INFINITY return with this sign
+	bchg	#31,d7		| else we know b is INFINITY and has
+	bra	Ld$infty	| the opposite sign
+
+|=============================================================================
+|                              __muldf3
+|=============================================================================
+
+| double __muldf3(double, double);
+SYM (__muldf3):
+	link	a6,#0
+	moveml	d2-d7,sp@-
+	movel	a6@(8),d0	| get a into d0-d1
+	movel	a6@(12),d1	| 
+	movel	a6@(16),d2	| and b into d2-d3
+	movel	a6@(20),d3	|
+	movel	d0,d7		| d7 will hold the sign of the product
+	eorl	d2,d7		|
+	andl	#0x80000000,d7	|
+	movel	d7,a0		| save sign bit into a0 
+	movel	#0x7ff00000,d7	| useful constant (+INFINITY)
+	movel	d7,d6		| another (mask for fraction)
+	notl	d6		| 
+	bclr	#31,d0		| get rid of a's sign bit '
+	movel	d0,d4		| 
+	orl	d1,d4		| 
+	beq	Lmuldf$a$0	| branch if a is zero
+	movel	d0,d4		| 
+	bclr	#31,d2		| get rid of b's sign bit '
+	movel	d2,d5		| 
+	orl	d3,d5		| 
+	beq	Lmuldf$b$0	| branch if b is zero
+	movel	d2,d5		| 
+	cmpl	d7,d0		| is a big?
+	bhi	Lmuldf$inop	| if a is NaN return NaN
+	beq	Lmuldf$a$nf	| we still have to check d1 and b ...
+	cmpl	d7,d2		| now compare b with INFINITY
+	bhi	Lmuldf$inop	| is b NaN?
+	beq	Lmuldf$b$nf 	| we still have to check d3 ...
+| Here we have both numbers finite and nonzero (and with no sign bit).
+| Now we get the exponents into d4 and d5.
+	andl	d7,d4		| isolate exponent in d4
+	beq	Lmuldf$a$den	| if exponent is zero we have a denormalized
+	andl	d6,d0		| isolate fraction
+	orl	#0x00100000,d0	| and put hidden bit back
+	swap	d4		| I like exponents in the first byte
+	lsrw	#4,d4		| 
+Lmuldf$1:			
+	andl	d7,d5		|
+	beq	Lmuldf$b$den	|
+	andl	d6,d2		|
+	orl	#0x00100000,d2	| and put hidden bit back
+	swap	d5		|
+	lsrw	#4,d5		|
+Lmuldf$2:			|
+	addw	d5,d4		| add exponents
+	subw	#D_BIAS+1,d4	| and substract bias (plus one)
+
+| We are now ready to do the multiplication. The situation is as follows:
+| both a and b have bit 52 ( bit 20 of d0 and d2) set (even if they were 
+| denormalized to start with!), which means that in the product bit 104 
+| (which will correspond to bit 8 of the fourth long) is set.
+
+| Here we have to do the product.
+| To do it we have to juggle the registers back and forth, as there are not
+| enough to keep everything in them. So we use the address registers to keep
+| some intermediate data.
+
+	moveml	a2-a3,sp@-	| save a2 and a3 for temporary use
+	movel	#0,a2		| a2 is a null register
+	movel	d4,a3		| and a3 will preserve the exponent
+
+| First, shift d2-d3 so bit 20 becomes bit 31:
+	rorl	#5,d2		| rotate d2 5 places right
+	swap	d2		| and swap it
+	rorl	#5,d3		| do the same thing with d3
+	swap	d3		|
+	movew	d3,d6		| get the rightmost 11 bits of d3
+	andw	#0x07ff,d6	|
+	orw	d6,d2		| and put them into d2
+	andw	#0xf800,d3	| clear those bits in d3
+
+	movel	d2,d6		| move b into d6-d7
+	movel	d3,d7           | move a into d4-d5
+	movel	d0,d4           | and clear d0-d1-d2-d3 (to put result)
+	movel	d1,d5           |
+	movel	#0,d3           |
+	movel	d3,d2           |
+	movel	d3,d1           |
+	movel	d3,d0	        |
+
+| We use a1 as counter:	
+	movel	#DBL_MANT_DIG-1,a1		
+	exg	d7,a1
+
+1:	exg	d7,a1		| put counter back in a1
+	addl	d3,d3		| shift sum once left
+	addxl	d2,d2           |
+	addxl	d1,d1           |
+	addxl	d0,d0           |
+	addl	d7,d7		|
+	addxl	d6,d6		|
+	bcc	2f		| if bit clear skip the following
+	exg	d7,a2		|
+	addl	d5,d3		| else add a to the sum
+	addxl	d4,d2		|
+	addxl	d7,d1		|
+	addxl	d7,d0		|
+	exg	d7,a2		| 
+2:	exg	d7,a1		| put counter in d7
+	dbf	d7,1b		| decrement and branch
+
+	movel	a3,d4		| restore exponent
+	moveml	sp@+,a2-a3
+
+| Now we have the product in d0-d1-d2-d3, with bit 8 of d0 set. The 
+| first thing to do now is to normalize it so bit 8 becomes bit 
+| DBL_MANT_DIG-32 (to do the rounding); later we will shift right.
+	swap	d0
+	swap	d1
+	movew	d1,d0
+	swap	d2
+	movew	d2,d1
+	swap	d3
+	movew	d3,d2
+	movew	#0,d3
+	lsrl	#1,d0
+	roxrl	#1,d1
+	roxrl	#1,d2
+	roxrl	#1,d3
+	lsrl	#1,d0
+	roxrl	#1,d1
+	roxrl	#1,d2
+	roxrl	#1,d3
+	lsrl	#1,d0
+	roxrl	#1,d1
+	roxrl	#1,d2
+	roxrl	#1,d3
+	
+| Now round, check for over- and underflow, and exit.
+	movel	a0,d7		| get sign bit back into d7
+	movew	#MULTIPLY,d5
+
+	btst	#DBL_MANT_DIG+1-32,d0
+	beq	Lround$exit
+	lsrl	#1,d0
+	roxrl	#1,d1
+	addw	#1,d4
+	bra	Lround$exit
+
+Lmuldf$inop:
+	movew	#MULTIPLY,d5
+	bra	Ld$inop
+
+Lmuldf$b$nf:
+	movew	#MULTIPLY,d5
+	movel	a0,d7		| get sign bit back into d7
+	tstl	d3		| we know d2 == 0x7ff00000, so check d3
+	bne	Ld$inop		| if d3 <> 0 b is NaN
+	bra	Ld$overflow	| else we have overflow (since a is finite)
+
+Lmuldf$a$nf:
+	movew	#MULTIPLY,d5
+	movel	a0,d7		| get sign bit back into d7
+	tstl	d1		| we know d0 == 0x7ff00000, so check d1
+	bne	Ld$inop		| if d1 <> 0 a is NaN
+	bra	Ld$overflow	| else signal overflow
+
+| If either number is zero return zero, unless the other is +/-INFINITY or
+| NaN, in which case we return NaN.
+Lmuldf$b$0:
+	movew	#MULTIPLY,d5
+	exg	d2,d0		| put b (==0) into d0-d1
+	exg	d3,d1		| and a (with sign bit cleared) into d2-d3
+	bra	1f
+Lmuldf$a$0:
+	movel	a6@(16),d2	| put b into d2-d3 again
+	movel	a6@(20),d3	|
+	bclr	#31,d2		| clear sign bit
+1:	cmpl	#0x7ff00000,d2	| check for non-finiteness
+	bge	Ld$inop		| in case NaN or +/-INFINITY return NaN
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7
+	unlk	a6
+	rts
+
+| If a number is denormalized we put an exponent of 1 but do not put the 
+| hidden bit back into the fraction; instead we shift left until bit 21
+| (the hidden bit) is set, adjusting the exponent accordingly. We do this
+| to ensure that the product of the fractions is close to 1.
+Lmuldf$a$den:
+	movel	#1,d4
+	andl	d6,d0
+1:	addl	d1,d1           | shift a left until bit 20 is set
+	addxl	d0,d0		|
+	subw	#1,d4		| and adjust exponent
+	btst	#20,d0          |
+	bne	Lmuldf$1        |
+	bra	1b
+
+Lmuldf$b$den:
+	movel	#1,d5
+	andl	d6,d2
+1:	addl	d3,d3		| shift b left until bit 20 is set
+	addxl	d2,d2		|
+	subw	#1,d5		| and adjust exponent
+	btst	#20,d2		|
+	bne	Lmuldf$2	|
+	bra	1b
+
+
+|=============================================================================
+|                              __divdf3
+|=============================================================================
+
+| double __divdf3(double, double);
+SYM (__divdf3):
+	link	a6,#0
+	moveml	d2-d7,sp@-
+	movel	a6@(8),d0	| get a into d0-d1
+	movel	a6@(12),d1	| 
+	movel	a6@(16),d2	| and b into d2-d3
+	movel	a6@(20),d3	|
+	movel	d0,d7		| d7 will hold the sign of the result
+	eorl	d2,d7		|
+	andl	#0x80000000,d7	| 
+	movel	d7,a0		| save sign into a0
+	movel	#0x7ff00000,d7	| useful constant (+INFINITY)
+	movel	d7,d6		| another (mask for fraction)
+	notl	d6		|
+	bclr	#31,d0		| get rid of a's sign bit '
+	movel	d0,d4		|
+	orl	d1,d4		|
+	beq	Ldivdf$a$0	| branch if a is zero
+	movel	d0,d4		|
+	bclr	#31,d2		| get rid of b's sign bit '
+	movel	d2,d5		|
+	orl	d3,d5		|
+	beq	Ldivdf$b$0	| branch if b is zero
+	movel	d2,d5
+	cmpl	d7,d0		| is a big?
+	bhi	Ldivdf$inop	| if a is NaN return NaN
+	beq	Ldivdf$a$nf	| if d0 == 0x7ff00000 we check d1
+	cmpl	d7,d2		| now compare b with INFINITY 
+	bhi	Ldivdf$inop	| if b is NaN return NaN
+	beq	Ldivdf$b$nf	| if d2 == 0x7ff00000 we check d3
+| Here we have both numbers finite and nonzero (and with no sign bit).
+| Now we get the exponents into d4 and d5 and normalize the numbers to
+| ensure that the ratio of the fractions is around 1. We do this by
+| making sure that both numbers have bit #DBL_MANT_DIG-32-1 (hidden bit)
+| set, even if they were denormalized to start with.
+| Thus, the result will satisfy: 2 > result > 1/2.
+	andl	d7,d4		| and isolate exponent in d4
+	beq	Ldivdf$a$den	| if exponent is zero we have a denormalized
+	andl	d6,d0		| and isolate fraction
+	orl	#0x00100000,d0	| and put hidden bit back
+	swap	d4		| I like exponents in the first byte
+	lsrw	#4,d4		| 
+Ldivdf$1:			| 
+	andl	d7,d5		|
+	beq	Ldivdf$b$den	|
+	andl	d6,d2		|
+	orl	#0x00100000,d2	|
+	swap	d5		|
+	lsrw	#4,d5		|
+Ldivdf$2:			|
+	subw	d5,d4		| substract exponents
+	addw	#D_BIAS,d4	| and add bias
+
+| We are now ready to do the division. We have prepared things in such a way
+| that the ratio of the fractions will be less than 2 but greater than 1/2.
+| At this point the registers in use are:
+| d0-d1	hold a (first operand, bit DBL_MANT_DIG-32=0, bit 
+| DBL_MANT_DIG-1-32=1)
+| d2-d3	hold b (second operand, bit DBL_MANT_DIG-32=1)
+| d4	holds the difference of the exponents, corrected by the bias
+| a0	holds the sign of the ratio
+
+| To do the rounding correctly we need to keep information about the
+| nonsignificant bits. One way to do this would be to do the division
+| using four registers; another is to use two registers (as originally
+| I did), but use a sticky bit to preserve information about the 
+| fractional part. Note that we can keep that info in a1, which is not
+| used.
+	movel	#0,d6		| d6-d7 will hold the result
+	movel	d6,d7		| 
+	movel	#0,a1		| and a1 will hold the sticky bit
+
+	movel	#DBL_MANT_DIG-32+1,d5	
+	
+1:	cmpl	d0,d2		| is a < b?
+	bhi	3f		| if b > a skip the following
+	beq	4f		| if d0==d2 check d1 and d3
+2:	subl	d3,d1		| 
+	subxl	d2,d0		| a <-- a - b
+	bset	d5,d6		| set the corresponding bit in d6
+3:	addl	d1,d1		| shift a by 1
+	addxl	d0,d0		|
+	dbra	d5,1b		| and branch back
+	bra	5f			
+4:	cmpl	d1,d3		| here d0==d2, so check d1 and d3
+	bhi	3b		| if d1 > d2 skip the substraction
+	bra	2b		| else go do it
+5:
+| Here we have to start setting the bits in the second long.
+	movel	#31,d5		| again d5 is counter
+
+1:	cmpl	d0,d2		| is a < b?
+	bhi	3f		| if b > a skip the following
+	beq	4f		| if d0==d2 check d1 and d3
+2:	subl	d3,d1		| 
+	subxl	d2,d0		| a <-- a - b
+	bset	d5,d7		| set the corresponding bit in d7
+3:	addl	d1,d1		| shift a by 1
+	addxl	d0,d0		|
+	dbra	d5,1b		| and branch back
+	bra	5f			
+4:	cmpl	d1,d3		| here d0==d2, so check d1 and d3
+	bhi	3b		| if d1 > d2 skip the substraction
+	bra	2b		| else go do it
+5:
+| Now go ahead checking until we hit a one, which we store in d2.
+	movel	#DBL_MANT_DIG,d5
+1:	cmpl	d2,d0		| is a < b?
+	bhi	4f		| if b < a, exit
+	beq	3f		| if d0==d2 check d1 and d3
+2:	addl	d1,d1		| shift a by 1
+	addxl	d0,d0		|
+	dbra	d5,1b		| and branch back
+	movel	#0,d2		| here no sticky bit was found
+	movel	d2,d3
+	bra	5f			
+3:	cmpl	d1,d3		| here d0==d2, so check d1 and d3
+	bhi	2b		| if d1 > d2 go back
+4:
+| Here put the sticky bit in d2-d3 (in the position which actually corresponds
+| to it; if you don't do this the algorithm loses in some cases). '
+	movel	#0,d2
+	movel	d2,d3
+	subw	#DBL_MANT_DIG,d5
+	addw	#63,d5
+	cmpw	#31,d5
+	bhi	2f
+1:	bset	d5,d3
+	bra	5f
+	subw	#32,d5
+2:	bset	d5,d2
+5:
+| Finally we are finished! Move the longs in the address registers to
+| their final destination:
+	movel	d6,d0
+	movel	d7,d1
+	movel	#0,d3
+
+| Here we have finished the division, with the result in d0-d1-d2-d3, with
+| 2^21 <= d6 < 2^23. Thus bit 23 is not set, but bit 22 could be set.
+| If it is not, then definitely bit 21 is set. Normalize so bit 22 is
+| not set:
+	btst	#DBL_MANT_DIG-32+1,d0
+	beq	1f
+	lsrl	#1,d0
+	roxrl	#1,d1
+	roxrl	#1,d2
+	roxrl	#1,d3
+	addw	#1,d4
+1:
+| Now round, check for over- and underflow, and exit.
+	movel	a0,d7		| restore sign bit to d7
+	movew	#DIVIDE,d5
+	bra	Lround$exit
+
+Ldivdf$inop:
+	movew	#DIVIDE,d5
+	bra	Ld$inop
+
+Ldivdf$a$0:
+| If a is zero check to see whether b is zero also. In that case return
+| NaN; then check if b is NaN, and return NaN also in that case. Else
+| return zero.
+	movew	#DIVIDE,d5
+	bclr	#31,d2		|
+	movel	d2,d4		| 
+	orl	d3,d4		| 
+	beq	Ld$inop		| if b is also zero return NaN
+	cmpl	#0x7ff00000,d2	| check for NaN
+	bhi	Ld$inop		| 
+	blt	1f		|
+	tstl	d3		|
+	bne	Ld$inop		|
+1:	movel	#0,d0		| else return zero
+	movel	d0,d1		| 
+	lea	SYM (_fpCCR),a0	| clear exception flags
+	movew	#0,a0@		|
+	moveml	sp@+,d2-d7	| 
+	unlk	a6		| 
+	rts			| 	
+
+Ldivdf$b$0:
+	movew	#DIVIDE,d5
+| If we got here a is not zero. Check if a is NaN; in that case return NaN,
+| else return +/-INFINITY. Remember that a is in d0 with the sign bit 
+| cleared already.
+	movel	a0,d7		| put a's sign bit back in d7 '
+	cmpl	#0x7ff00000,d0	| compare d0 with INFINITY
+	bhi	Ld$inop		| if larger it is NaN
+	tstl	d1		| 
+	bne	Ld$inop		| 
+	bra	Ld$div$0	| else signal DIVIDE_BY_ZERO
+
+Ldivdf$b$nf:
+	movew	#DIVIDE,d5
+| If d2 == 0x7ff00000 we have to check d3.
+	tstl	d3		|
+	bne	Ld$inop		| if d3 <> 0, b is NaN
+	bra	Ld$underflow	| else b is +/-INFINITY, so signal underflow
+
+Ldivdf$a$nf:
+	movew	#DIVIDE,d5
+| If d0 == 0x7ff00000 we have to check d1.
+	tstl	d1		|
+	bne	Ld$inop		| if d1 <> 0, a is NaN
+| If a is INFINITY we have to check b
+	cmpl	d7,d2		| compare b with INFINITY 
+	bge	Ld$inop		| if b is NaN or INFINITY return NaN
+	tstl	d3		|
+	bne	Ld$inop		| 
+	bra	Ld$overflow	| else return overflow
+
+| If a number is denormalized we put an exponent of 1 but do not put the 
+| bit back into the fraction.
+Ldivdf$a$den:
+	movel	#1,d4
+	andl	d6,d0
+1:	addl	d1,d1		| shift a left until bit 20 is set
+	addxl	d0,d0
+	subw	#1,d4		| and adjust exponent
+	btst	#DBL_MANT_DIG-32-1,d0
+	bne	Ldivdf$1
+	bra	1b
+
+Ldivdf$b$den:
+	movel	#1,d5
+	andl	d6,d2
+1:	addl	d3,d3		| shift b left until bit 20 is set
+	addxl	d2,d2
+	subw	#1,d5		| and adjust exponent
+	btst	#DBL_MANT_DIG-32-1,d2
+	bne	Ldivdf$2
+	bra	1b
+
+Lround$exit:
+| This is a common exit point for __muldf3 and __divdf3. When they enter
+| this point the sign of the result is in d7, the result in d0-d1, normalized
+| so that 2^21 <= d0 < 2^22, and the exponent is in the lower byte of d4.
+
+| First check for underlow in the exponent:
+	cmpw	#-DBL_MANT_DIG-1,d4		
+	blt	Ld$underflow	
+| It could happen that the exponent is less than 1, in which case the 
+| number is denormalized. In this case we shift right and adjust the 
+| exponent until it becomes 1 or the fraction is zero (in the latter case 
+| we signal underflow and return zero).
+	movel	d7,a0		|
+	movel	#0,d6		| use d6-d7 to collect bits flushed right
+	movel	d6,d7		| use d6-d7 to collect bits flushed right
+	cmpw	#1,d4		| if the exponent is less than 1 we 
+	bge	2f		| have to shift right (denormalize)
+1:	addw	#1,d4		| adjust the exponent
+	lsrl	#1,d0		| shift right once 
+	roxrl	#1,d1		|
+	roxrl	#1,d2		|
+	roxrl	#1,d3		|
+	roxrl	#1,d6		| 
+	roxrl	#1,d7		|
+	cmpw	#1,d4		| is the exponent 1 already?
+	beq	2f		| if not loop back
+	bra	1b              |
+	bra	Ld$underflow	| safety check, shouldn't execute '
+2:	orl	d6,d2		| this is a trick so we don't lose  '
+	orl	d7,d3		| the bits which were flushed right
+	movel	a0,d7		| get back sign bit into d7
+| Now call the rounding routine (which takes care of denormalized numbers):
+	lea	Lround$0,a0	| to return from rounding routine
+	lea	SYM (_fpCCR),a1	| check the rounding mode
+	movew	a1@(6),d6	| rounding mode in d6
+	beq	Lround$to$nearest
+	cmpw	#ROUND_TO_PLUS,d6
+	bhi	Lround$to$minus
+	blt	Lround$to$zero
+	bra	Lround$to$plus
+Lround$0:
+| Here we have a correctly rounded result (either normalized or denormalized).
+
+| Here we should have either a normalized number or a denormalized one, and
+| the exponent is necessarily larger or equal to 1 (so we don't have to  '
+| check again for underflow!). We have to check for overflow or for a 
+| denormalized number (which also signals underflow).
+| Check for overflow (i.e., exponent >= 0x7ff).
+	cmpw	#0x07ff,d4
+	bge	Ld$overflow
+| Now check for a denormalized number (exponent==0):
+	movew	d4,d4
+	beq	Ld$den
+1:
+| Put back the exponents and sign and return.
+	lslw	#4,d4		| exponent back to fourth byte
+	bclr	#DBL_MANT_DIG-32-1,d0
+	swap	d0		| and put back exponent
+	orw	d4,d0		| 
+	swap	d0		|
+	orl	d7,d0		| and sign also
+
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7
+	unlk	a6
+	rts
+
+|=============================================================================
+|                              __negdf2
+|=============================================================================
+
+| double __negdf2(double, double);
+SYM (__negdf2):
+	link	a6,#0
+	moveml	d2-d7,sp@-
+	movew	#NEGATE,d5
+	movel	a6@(8),d0	| get number to negate in d0-d1
+	movel	a6@(12),d1	|
+	bchg	#31,d0		| negate
+	movel	d0,d2		| make a positive copy (for the tests)
+	bclr	#31,d2		|
+	movel	d2,d4		| check for zero
+	orl	d1,d4		|
+	beq	2f		| if zero (either sign) return +zero
+	cmpl	#0x7ff00000,d2	| compare to +INFINITY
+	blt	1f		| if finite, return
+	bhi	Ld$inop		| if larger (fraction not zero) is NaN
+	tstl	d1		| if d2 == 0x7ff00000 check d1
+	bne	Ld$inop		|
+	movel	d0,d7		| else get sign and return INFINITY
+	andl	#0x80000000,d7
+	bra	Ld$infty		
+1:	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7
+	unlk	a6
+	rts
+2:	bclr	#31,d0
+	bra	1b
+
+|=============================================================================
+|                              __cmpdf2
+|=============================================================================
+
+GREATER =  1
+LESS    = -1
+EQUAL   =  0
+
+| int __cmpdf2(double, double);
+SYM (__cmpdf2):
+	link	a6,#0
+	moveml	d2-d7,sp@- 	| save registers
+	movew	#COMPARE,d5
+	movel	a6@(8),d0	| get first operand
+	movel	a6@(12),d1	|
+	movel	a6@(16),d2	| get second operand
+	movel	a6@(20),d3	|
+| First check if a and/or b are (+/-) zero and in that case clear
+| the sign bit.
+	movel	d0,d6		| copy signs into d6 (a) and d7(b)
+	bclr	#31,d0		| and clear signs in d0 and d2
+	movel	d2,d7		|
+	bclr	#31,d2		|
+	cmpl	#0x7fff0000,d0	| check for a == NaN
+	bhi	Ld$inop		| if d0 > 0x7ff00000, a is NaN
+	beq	Lcmpdf$a$nf	| if equal can be INFINITY, so check d1
+	movel	d0,d4		| copy into d4 to test for zero
+	orl	d1,d4		|
+	beq	Lcmpdf$a$0	|
+Lcmpdf$0:
+	cmpl	#0x7fff0000,d2	| check for b == NaN
+	bhi	Ld$inop		| if d2 > 0x7ff00000, b is NaN
+	beq	Lcmpdf$b$nf	| if equal can be INFINITY, so check d3
+	movel	d2,d4		|
+	orl	d3,d4		|
+	beq	Lcmpdf$b$0	|
+Lcmpdf$1:
+| Check the signs
+	eorl	d6,d7
+	bpl	1f
+| If the signs are not equal check if a >= 0
+	tstl	d6
+	bpl	Lcmpdf$a$gt$b	| if (a >= 0 && b < 0) => a > b
+	bmi	Lcmpdf$b$gt$a	| if (a < 0 && b >= 0) => a < b
+1:
+| If the signs are equal check for < 0
+	tstl	d6
+	bpl	1f
+| If both are negative exchange them
+	exg	d0,d2
+	exg	d1,d3
+1:
+| Now that they are positive we just compare them as longs (does this also
+| work for denormalized numbers?).
+	cmpl	d0,d2
+	bhi	Lcmpdf$b$gt$a	| |b| > |a|
+	bne	Lcmpdf$a$gt$b	| |b| < |a|
+| If we got here d0 == d2, so we compare d1 and d3.
+	cmpl	d1,d3
+	bhi	Lcmpdf$b$gt$a	| |b| > |a|
+	bne	Lcmpdf$a$gt$b	| |b| < |a|
+| If we got here a == b.
+	movel	#EQUAL,d0
+	moveml	sp@+,d2-d7 	| put back the registers
+	unlk	a6
+	rts
+Lcmpdf$a$gt$b:
+	movel	#GREATER,d0
+	moveml	sp@+,d2-d7 	| put back the registers
+	unlk	a6
+	rts
+Lcmpdf$b$gt$a:
+	movel	#LESS,d0
+	moveml	sp@+,d2-d7 	| put back the registers
+	unlk	a6
+	rts
+
+Lcmpdf$a$0:	
+	bclr	#31,d6
+	bra	Lcmpdf$0
+Lcmpdf$b$0:
+	bclr	#31,d7
+	bra	Lcmpdf$1
+
+Lcmpdf$a$nf:
+	tstl	d1
+	bne	Ld$inop
+	bra	Lcmpdf$0
+
+Lcmpdf$b$nf:
+	tstl	d3
+	bne	Ld$inop
+	bra	Lcmpdf$1
+
+|=============================================================================
+|                           rounding routines
+|=============================================================================
+
+| The rounding routines expect the number to be normalized in registers
+| d0-d1-d2-d3, with the exponent in register d4. They assume that the 
+| exponent is larger or equal to 1. They return a properly normalized number
+| if possible, and a denormalized number otherwise. The exponent is returned
+| in d4.
+
+Lround$to$nearest:
+| We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
+| Here we assume that the exponent is not too small (this should be checked
+| before entering the rounding routine), but the number could be denormalized.
+
+| Check for denormalized numbers:
+1:	btst	#DBL_MANT_DIG-32,d0
+	bne	2f		| if set the number is normalized
+| Normalize shifting left until bit #DBL_MANT_DIG-32 is set or the exponent 
+| is one (remember that a denormalized number corresponds to an 
+| exponent of -D_BIAS+1).
+	cmpw	#1,d4		| remember that the exponent is at least one
+ 	beq	2f		| an exponent of one means denormalized
+	addl	d3,d3		| else shift and adjust the exponent
+	addxl	d2,d2		|
+	addxl	d1,d1		|
+	addxl	d0,d0		|
+	dbra	d4,1b		|
+2:
+| Now round: we do it as follows: after the shifting we can write the
+| fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
+| If delta < 1, do nothing. If delta > 1, add 1 to f. 
+| If delta == 1, we make sure the rounded number will be even (odd?) 
+| (after shifting).
+	btst	#0,d1		| is delta < 1?
+	beq	2f		| if so, do not do anything
+	orl	d2,d3		| is delta == 1?
+	bne	1f		| if so round to even
+	movel	d1,d3		| 
+	andl	#2,d3		| bit 1 is the last significant bit
+	movel	#0,d2		|
+	addl	d3,d1		|
+	addxl	d2,d0		|
+	bra	2f		| 
+1:	movel	#1,d3		| else add 1 
+	movel	#0,d2		|
+	addl	d3,d1		|
+	addxl	d2,d0
+| Shift right once (because we used bit #DBL_MANT_DIG-32!).
+2:	lsrl	#1,d0
+	roxrl	#1,d1		
+
+| Now check again bit #DBL_MANT_DIG-32 (rounding could have produced a
+| 'fraction overflow' ...).
+	btst	#DBL_MANT_DIG-32,d0	
+	beq	1f
+	lsrl	#1,d0
+	roxrl	#1,d1
+	addw	#1,d4
+1:
+| If bit #DBL_MANT_DIG-32-1 is clear we have a denormalized number, so we 
+| have to put the exponent to zero and return a denormalized number.
+	btst	#DBL_MANT_DIG-32-1,d0
+	beq	1f
+	jmp	a0@
+1:	movel	#0,d4
+	jmp	a0@
+
+Lround$to$zero:
+Lround$to$plus:
+Lround$to$minus:
+	jmp	a0@
+#endif /* L_double */
+
+#ifdef  L_float
+
+	.globl	SYM (_fpCCR)
+	.globl  $_exception_handler
+
+QUIET_NaN    = 0xffffffff
+SIGNL_NaN    = 0x7f800001
+INFINITY     = 0x7f800000
+
+F_MAX_EXP      = 0xff
+F_BIAS         = 126
+FLT_MAX_EXP    = F_MAX_EXP - F_BIAS
+FLT_MIN_EXP    = 1 - F_BIAS
+FLT_MANT_DIG   = 24
+
+INEXACT_RESULT 		= 0x0001
+UNDERFLOW 		= 0x0002
+OVERFLOW 		= 0x0004
+DIVIDE_BY_ZERO 		= 0x0008
+INVALID_OPERATION 	= 0x0010
+
+SINGLE_FLOAT = 1
+
+NOOP         = 0
+ADD          = 1
+MULTIPLY     = 2
+DIVIDE       = 3
+NEGATE       = 4
+COMPARE      = 5
+EXTENDSFDF   = 6
+TRUNCDFSF    = 7
+
+UNKNOWN           = -1
+ROUND_TO_NEAREST  = 0 | round result to nearest representable value
+ROUND_TO_ZERO     = 1 | round result towards zero
+ROUND_TO_PLUS     = 2 | round result towards plus infinity
+ROUND_TO_MINUS    = 3 | round result towards minus infinity
+
+| Entry points:
+
+	.globl SYM (__addsf3)
+	.globl SYM (__subsf3)
+	.globl SYM (__mulsf3)
+	.globl SYM (__divsf3)
+	.globl SYM (__negsf2)
+	.globl SYM (__cmpsf2)
+
+| These are common routines to return and signal exceptions.	
+
+	.text
+	.even
+
+Lf$den:
+| Return and signal a denormalized number
+	orl	d7,d0
+	movew	#UNDERFLOW,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#SINGLE_FLOAT,d6
+	jmp	$_exception_handler
+
+Lf$infty:
+Lf$overflow:
+| Return a properly signed INFINITY and set the exception flags 
+	movel	#INFINITY,d0
+	orl	d7,d0
+	movew	#OVERFLOW,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#SINGLE_FLOAT,d6
+	jmp	$_exception_handler
+
+Lf$underflow:
+| Return 0 and set the exception flags 
+	movel	#0,d0
+	movew	#UNDERFLOW,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#SINGLE_FLOAT,d6
+	jmp	$_exception_handler
+
+Lf$inop:
+| Return a quiet NaN and set the exception flags
+	movel	#QUIET_NaN,d0
+	movew	#INVALID_OPERATION,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#SINGLE_FLOAT,d6
+	jmp	$_exception_handler
+
+Lf$div$0:
+| Return a properly signed INFINITY and set the exception flags
+	movel	#INFINITY,d0
+	orl	d7,d0
+	movew	#DIVIDE_BY_ZERO,d7
+	orw	#INEXACT_RESULT,d7
+	movew	#SINGLE_FLOAT,d6
+	jmp	$_exception_handler
+
+|=============================================================================
+|=============================================================================
+|                         single precision routines
+|=============================================================================
+|=============================================================================
+
+| A single precision floating point number (float) has the format:
+|
+| struct _float {
+|  unsigned int sign      : 1;  /* sign bit */ 
+|  unsigned int exponent  : 8;  /* exponent, shifted by 126 */
+|  unsigned int fraction  : 23; /* fraction */
+| } float;
+| 
+| Thus sizeof(float) = 4 (32 bits). 
+|
+| All the routines are callable from C programs, and return the result 
+| in the single register d0. They also preserve all registers except 
+| d0-d1 and a0-a1.
+
+|=============================================================================
+|                              __subsf3
+|=============================================================================
+
+| float __subsf3(float, float);
+SYM (__subsf3):
+	bchg	#31,sp@(8)	| change sign of second operand
+				| and fall through
+|=============================================================================
+|                              __addsf3
+|=============================================================================
+
+| float __addsf3(float, float);
+SYM (__addsf3):
+	link	a6,#0		| everything will be done in registers
+	moveml	d2-d7,sp@-	| save all data registers but d0-d1
+	movel	a6@(8),d0	| get first operand
+	movel	a6@(12),d1	| get second operand
+	movel	d0,d6		| get d0's sign bit '
+	addl	d0,d0		| check and clear sign bit of a
+	beq	Laddsf$b	| if zero return second operand
+	movel	d1,d7		| save b's sign bit '
+	addl	d1,d1		| get rid of sign bit
+	beq	Laddsf$a	| if zero return first operand
+
+	movel	d6,a0		| save signs in address registers
+	movel	d7,a1		| so we can use d6 and d7
+
+| Get the exponents and check for denormalized and/or infinity.
+
+	movel	#0x00ffffff,d4	| mask to get fraction
+	movel	#0x01000000,d5	| mask to put hidden bit back
+
+	movel	d0,d6		| save a to get exponent
+	andl	d4,d0		| get fraction in d0
+	notl 	d4		| make d4 into a mask for the exponent
+	andl	d4,d6		| get exponent in d6
+	beq	Laddsf$a$den	| branch if a is denormalized
+	cmpl	d4,d6		| check for INFINITY or NaN
+	beq	Laddsf$nf
+	swap	d6		| put exponent into first word
+	orl	d5,d0		| and put hidden bit back
+Laddsf$1:
+| Now we have a's exponent in d6 (second byte) and the mantissa in d0. '
+	movel	d1,d7		| get exponent in d7
+	andl	d4,d7		| 
+	beq	Laddsf$b$den	| branch if b is denormalized
+	cmpl	d4,d7		| check for INFINITY or NaN
+	beq	Laddsf$nf
+	swap	d7		| put exponent into first word
+	notl 	d4		| make d4 into a mask for the fraction
+	andl	d4,d1		| get fraction in d1
+	orl	d5,d1		| and put hidden bit back
+Laddsf$2:
+| Now we have b's exponent in d7 (second byte) and the mantissa in d1. '
+
+| Note that the hidden bit corresponds to bit #FLT_MANT_DIG-1, and we 
+| shifted right once, so bit #FLT_MANT_DIG is set (so we have one extra
+| bit).
+
+	movel	d1,d2		| move b to d2, since we want to use
+				| two registers to do the sum
+	movel	#0,d1		| and clear the new ones
+	movel	d1,d3		|
+
+| Here we shift the numbers in registers d0 and d1 so the exponents are the
+| same, and put the largest exponent in d6. Note that we are using two
+| registers for each number (see the discussion by D. Knuth in "Seminumerical 
+| Algorithms").
+	cmpw	d6,d7		| compare exponents
+	beq	Laddsf$3	| if equal don't shift '
+	bhi	5f		| branch if second exponent largest
+1:
+	subl	d6,d7		| keep the largest exponent
+	negl	d7
+	lsrw	#8,d7		| put difference in lower byte
+| if difference is too large we don't shift (actually, we can just exit) '
+	cmpw	#FLT_MANT_DIG+2,d7		
+	bge	Laddsf$b$small
+	cmpw	#16,d7		| if difference >= 16 swap
+	bge	4f
+2:
+	subw	#1,d7
+3:	lsrl	#1,d2		| shift right second operand
+	roxrl	#1,d3
+	dbra	d7,3b
+	bra	Laddsf$3
+4:
+	movew	d2,d3
+	swap	d3
+	movew	d3,d2
+	swap	d2
+	subw	#16,d7
+	bra	2b
+5:
+	exg	d6,d7		| exchange the exponents
+	subl	d6,d7		| keep the largest exponent
+	negl	d7		|
+	lsrw	#8,d7		| put difference in lower byte
+| if difference is too large we don't shift (and exit!) '
+	cmpw	#FLT_MANT_DIG+2,d7		
+	bge	Laddsf$a$small
+	cmpw	#16,d7		| if difference >= 16 swap
+	bge	8f
+6:
+	subw	#1,d7
+7:	lsrl	#1,d0		| shift right first operand
+	roxrl	#1,d1
+	dbra	d7,7b
+	bra	Laddsf$3
+8:
+	movew	d0,d1
+	swap	d1
+	movew	d1,d0
+	swap	d0
+	subw	#16,d7
+	bra	6b
+
+| Now we have a in d0-d1, b in d2-d3, and the largest exponent in d6 (the
+| signs are stored in a0 and a1).
+
+Laddsf$3:
+| Here we have to decide whether to add or substract the numbers
+	exg	d6,a0		| get signs back
+	exg	d7,a1		| and save the exponents
+	eorl	d6,d7		| combine sign bits
+	bmi	Lsubsf$0	| if negative a and b have opposite 
+				| sign so we actually substract the
+				| numbers
+
+| Here we have both positive or both negative
+	exg	d6,a0		| now we have the exponent in d6
+	movel	a0,d7		| and sign in d7
+	andl	#0x80000000,d7	|
+| Here we do the addition.
+	addl	d3,d1
+	addxl	d2,d0
+| Note: now we have d2, d3, d4 and d5 to play with! 
+
+| Put the exponent, in the first byte, in d2, to use the "standard" rounding
+| routines:
+	movel	d6,d2
+	lsrw	#8,d2
+
+| Before rounding normalize so bit #FLT_MANT_DIG is set (we will consider
+| the case of denormalized numbers in the rounding routine itself).
+| As in the addition (not in the substraction!) we could have set 
+| one more bit we check this:
+	btst	#FLT_MANT_DIG+1,d0	
+	beq	1f
+	lsrl	#1,d0
+	roxrl	#1,d1
+	addl	#1,d2
+1:
+	lea	Laddsf$4,a0	| to return from rounding routine
+	lea	SYM (_fpCCR),a1	| check the rounding mode
+	movew	a1@(6),d6	| rounding mode in d6
+	beq	Lround$to$nearest
+	cmpw	#ROUND_TO_PLUS,d6
+	bhi	Lround$to$minus
+	blt	Lround$to$zero
+	bra	Lround$to$plus
+Laddsf$4:
+| Put back the exponent, but check for overflow.
+	cmpw	#0xff,d2
+	bhi	1f
+	bclr	#FLT_MANT_DIG-1,d0
+	lslw	#7,d2
+	swap	d2
+	orl	d2,d0
+	bra	Laddsf$ret
+1:
+	movew	#ADD,d5
+	bra	Lf$overflow
+
+Lsubsf$0:
+| We are here if a > 0 and b < 0 (sign bits cleared).
+| Here we do the substraction.
+	movel	d6,d7		| put sign in d7
+	andl	#0x80000000,d7	|
+
+	subl	d3,d1		| result in d0-d1
+	subxl	d2,d0		|
+	beq	Laddsf$ret	| if zero just exit
+	bpl	1f		| if positive skip the following
+	bchg	#31,d7		| change sign bit in d7
+	negl	d1
+	negxl	d0
+1:
+	exg	d2,a0		| now we have the exponent in d2
+	lsrw	#8,d2		| put it in the first byte
+
+| Now d0-d1 is positive and the sign bit is in d7.
+
+| Note that we do not have to normalize, since in the substraction bit
+| #FLT_MANT_DIG+1 is never set, and denormalized numbers are handled by
+| the rounding routines themselves.
+	lea	Lsubsf$1,a0	| to return from rounding routine
+	lea	SYM (_fpCCR),a1	| check the rounding mode
+	movew	a1@(6),d6	| rounding mode in d6
+	beq	Lround$to$nearest
+	cmpw	#ROUND_TO_PLUS,d6
+	bhi	Lround$to$minus
+	blt	Lround$to$zero
+	bra	Lround$to$plus
+Lsubsf$1:
+| Put back the exponent (we can't have overflow!). '
+	bclr	#FLT_MANT_DIG-1,d0
+	lslw	#7,d2
+	swap	d2
+	orl	d2,d0
+	bra	Laddsf$ret
+
+| If one of the numbers was too small (difference of exponents >= 
+| FLT_MANT_DIG+2) we return the other (and now we don't have to '
+| check for finiteness or zero).
+Laddsf$a$small:
+	movel	a6@(12),d0
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7	| restore data registers
+	unlk	a6		| and return
+	rts
+
+Laddsf$b$small:
+	movel	a6@(8),d0
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7	| restore data registers
+	unlk	a6		| and return
+	rts
+
+| If the numbers are denormalized remember to put exponent equal to 1.
+
+Laddsf$a$den:
+	movel	d5,d6		| d5 contains 0x01000000
+	swap	d6
+	bra	Laddsf$1
+
+Laddsf$b$den:
+	movel	d5,d7
+	swap	d7
+	notl 	d4		| make d4 into a mask for the fraction
+				| (this was not executed after the jump)
+	bra	Laddsf$2
+
+| The rest is mainly code for the different results which can be 
+| returned (checking always for +/-INFINITY and NaN).
+
+Laddsf$b:
+| Return b (if a is zero).
+	movel	a6@(12),d0
+	bra	1f
+Laddsf$a:
+| Return a (if b is zero).
+	movel	a6@(8),d0
+1:
+	movew	#ADD,d5
+| We have to check for NaN and +/-infty.
+	movel	d0,d7
+	andl	#0x80000000,d7	| put sign in d7
+	bclr	#31,d0		| clear sign
+	cmpl	#INFINITY,d0	| check for infty or NaN
+	bge	2f
+	movel	d0,d0		| check for zero (we do this because we don't '
+	bne	Laddsf$ret	| want to return -0 by mistake
+	bclr	#31,d7		| if zero be sure to clear sign
+	bra	Laddsf$ret	| if everything OK just return
+2:
+| The value to be returned is either +/-infty or NaN
+	andl	#0x007fffff,d0	| check for NaN
+	bne	Lf$inop		| if mantissa not zero is NaN
+	bra	Lf$infty
+
+Laddsf$ret:
+| Normal exit (a and b nonzero, result is not NaN nor +/-infty).
+| We have to clear the exception flags (just the exception type).
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	orl	d7,d0		| put sign bit
+	moveml	sp@+,d2-d7	| restore data registers
+	unlk	a6		| and return
+	rts
+
+Laddsf$ret$den:
+| Return a denormalized number (for addition we don't signal underflow) '
+	lsrl	#1,d0		| remember to shift right back once
+	bra	Laddsf$ret	| and return
+
+| Note: when adding two floats of the same sign if either one is 
+| NaN we return NaN without regard to whether the other is finite or 
+| not. When substracting them (i.e., when adding two numbers of 
+| opposite signs) things are more complicated: if both are INFINITY 
+| we return NaN, if only one is INFINITY and the other is NaN we return
+| NaN, but if it is finite we return INFINITY with the corresponding sign.
+
+Laddsf$nf:
+	movew	#ADD,d5
+| This could be faster but it is not worth the effort, since it is not
+| executed very often. We sacrifice speed for clarity here.
+	movel	a6@(8),d0	| get the numbers back (remember that we
+	movel	a6@(12),d1	| did some processing already)
+	movel	#INFINITY,d4	| useful constant (INFINITY)
+	movel	d0,d2		| save sign bits
+	movel	d1,d3
+	bclr	#31,d0		| clear sign bits
+	bclr	#31,d1
+| We know that one of them is either NaN of +/-INFINITY
+| Check for NaN (if either one is NaN return NaN)
+	cmpl	d4,d0		| check first a (d0)
+	bhi	Lf$inop		
+	cmpl	d4,d1		| check now b (d1)
+	bhi	Lf$inop		
+| Now comes the check for +/-INFINITY. We know that both are (maybe not
+| finite) numbers, but we have to check if both are infinite whether we
+| are adding or substracting them.
+	eorl	d3,d2		| to check sign bits
+	bmi	1f
+	movel	d0,d7
+	andl	#0x80000000,d7	| get (common) sign bit
+	bra	Lf$infty
+1:
+| We know one (or both) are infinite, so we test for equality between the
+| two numbers (if they are equal they have to be infinite both, so we
+| return NaN).
+	cmpl	d1,d0		| are both infinite?
+	beq	Lf$inop		| if so return NaN
+
+	movel	d0,d7
+	andl	#0x80000000,d7	| get a's sign bit '
+	cmpl	d4,d0		| test now for infinity
+	beq	Lf$infty	| if a is INFINITY return with this sign
+	bchg	#31,d7		| else we know b is INFINITY and has
+	bra	Lf$infty	| the opposite sign
+
+|=============================================================================
+|                             __mulsf3
+|=============================================================================
+
+| float __mulsf3(float, float);
+SYM (__mulsf3):
+	link	a6,#0
+	moveml	d2-d7,sp@-
+	movel	a6@(8),d0	| get a into d0
+	movel	a6@(12),d1	| and b into d1
+	movel	d0,d7		| d7 will hold the sign of the product
+	eorl	d1,d7		|
+	andl	#0x80000000,d7	| 
+	movel	#INFINITY,d6	| useful constant (+INFINITY)
+	movel	d6,d5		| another (mask for fraction)
+	notl	d5		|
+	movel	#0x00800000,d4	| this is to put hidden bit back
+	bclr	#31,d0		| get rid of a's sign bit '
+	movel	d0,d2		|
+	beq	Lmulsf$a$0	| branch if a is zero
+	bclr	#31,d1		| get rid of b's sign bit '
+	movel	d1,d3		|
+	beq	Lmulsf$b$0	| branch if b is zero
+	cmpl	d6,d0		| is a big?
+	bhi	Lmulsf$inop	| if a is NaN return NaN
+	beq	Lmulsf$inf	| if a is INFINITY we have to check b
+	cmpl	d6,d1		| now compare b with INFINITY
+	bhi	Lmulsf$inop	| is b NaN?
+	beq	Lmulsf$overflow | is b INFINITY?
+| Here we have both numbers finite and nonzero (and with no sign bit).
+| Now we get the exponents into d2 and d3.
+	andl	d6,d2		| and isolate exponent in d2
+	beq	Lmulsf$a$den	| if exponent is zero we have a denormalized
+	andl	d5,d0		| and isolate fraction
+	orl	d4,d0		| and put hidden bit back
+	swap	d2		| I like exponents in the first byte
+	lsrw	#7,d2		| 
+Lmulsf$1:			| number
+	andl	d6,d3		|
+	beq	Lmulsf$b$den	|
+	andl	d5,d1		|
+	orl	d4,d1		|
+	swap	d3		|
+	lsrw	#7,d3		|
+Lmulsf$2:			|
+	addw	d3,d2		| add exponents
+	subw	#F_BIAS+1,d2	| and substract bias (plus one)
+
+| We are now ready to do the multiplication. The situation is as follows:
+| both a and b have bit FLT_MANT_DIG-1 set (even if they were 
+| denormalized to start with!), which means that in the product 
+| bit 2*(FLT_MANT_DIG-1) (that is, bit 2*FLT_MANT_DIG-2-32 of the 
+| high long) is set. 
+
+| To do the multiplication let us move the number a little bit around ...
+	movel	d1,d6		| second operand in d6
+	movel	d0,d5		| first operand in d4-d5
+	movel	#0,d4
+	movel	d4,d1		| the sums will go in d0-d1
+	movel	d4,d0
+
+| now bit FLT_MANT_DIG-1 becomes bit 31:
+	lsll	#31-FLT_MANT_DIG+1,d6		
+
+| Start the loop (we loop #FLT_MANT_DIG times):
+	movew	#FLT_MANT_DIG-1,d3	
+1:	addl	d1,d1		| shift sum 
+	addxl	d0,d0
+	lsll	#1,d6		| get bit bn
+	bcc	2f		| if not set skip sum
+	addl	d5,d1		| add a
+	addxl	d4,d0
+2:	dbf	d3,1b		| loop back
+
+| Now we have the product in d0-d1, with bit (FLT_MANT_DIG - 1) + FLT_MANT_DIG
+| (mod 32) of d0 set. The first thing to do now is to normalize it so bit 
+| FLT_MANT_DIG is set (to do the rounding).
+	rorl	#6,d1
+	swap	d1
+	movew	d1,d3
+	andw	#0x03ff,d3
+	andw	#0xfd00,d1
+	lsll	#8,d0
+	addl	d0,d0
+	addl	d0,d0
+	orw	d3,d0
+
+	movew	#MULTIPLY,d5
+	
+	btst	#FLT_MANT_DIG+1,d0
+	beq	Lround$exit
+	lsrl	#1,d0
+	roxrl	#1,d1
+	addw	#1,d2
+	bra	Lround$exit
+
+Lmulsf$inop:
+	movew	#MULTIPLY,d5
+	bra	Lf$inop
+
+Lmulsf$overflow:
+	movew	#MULTIPLY,d5
+	bra	Lf$overflow
+
+Lmulsf$inf:
+	movew	#MULTIPLY,d5
+| If either is NaN return NaN; else both are (maybe infinite) numbers, so
+| return INFINITY with the correct sign (which is in d7).
+	cmpl	d6,d1		| is b NaN?
+	bhi	Lf$inop		| if so return NaN
+	bra	Lf$overflow	| else return +/-INFINITY
+
+| If either number is zero return zero, unless the other is +/-INFINITY, 
+| or NaN, in which case we return NaN.
+Lmulsf$b$0:
+| Here d1 (==b) is zero.
+	movel	d1,d0		| put b into d0 (just a zero)
+	movel	a6@(8),d1	| get a again to check for non-finiteness
+	bra	1f
+Lmulsf$a$0:
+	movel	a6@(12),d1	| get b again to check for non-finiteness
+1:	bclr	#31,d1		| clear sign bit 
+	cmpl	#INFINITY,d1	| and check for a large exponent
+	bge	Lf$inop		| if b is +/-INFINITY or NaN return NaN
+	lea	SYM (_fpCCR),a0	| else return zero
+	movew	#0,a0@		| 
+	moveml	sp@+,d2-d7	| 
+	unlk	a6		| 
+	rts			| 
+
+| If a number is denormalized we put an exponent of 1 but do not put the 
+| hidden bit back into the fraction; instead we shift left until bit 23
+| (the hidden bit) is set, adjusting the exponent accordingly. We do this
+| to ensure that the product of the fractions is close to 1.
+Lmulsf$a$den:
+	movel	#1,d2
+	andl	d5,d0
+1:	addl	d0,d0		| shift a left (until bit 23 is set)
+	subw	#1,d2		| and adjust exponent
+	btst	#FLT_MANT_DIG-1,d0
+	bne	Lmulsf$1	|
+	bra	1b		| else loop back
+
+Lmulsf$b$den:
+	movel	#1,d3
+	andl	d5,d1
+1:	addl	d1,d1		| shift b left until bit 23 is set
+	subw	#1,d3		| and adjust exponent
+	btst	#FLT_MANT_DIG-1,d1
+	bne	Lmulsf$2	|
+	bra	1b		| else loop back
+
+|=============================================================================
+|                             __divsf3
+|=============================================================================
+
+| float __divsf3(float, float);
+SYM (__divsf3):
+	link	a6,#0
+	moveml	d2-d7,sp@-
+	movel	a6@(8),d0	| get a into d0
+	movel	a6@(12),d1	| and b into d1
+	movel	d0,d7		| d7 will hold the sign of the result
+	eorl	d1,d7		|
+	andl	#0x80000000,d7	| 
+	movel	#INFINITY,d6	| useful constant (+INFINITY)
+	movel	d6,d5		| another (mask for fraction)
+	notl	d5		|
+	movel	#0x00800000,d4	| this is to put hidden bit back
+	bclr	#31,d0		| get rid of a's sign bit '
+	movel	d0,d2		|
+	beq	Ldivsf$a$0	| branch if a is zero
+	bclr	#31,d1		| get rid of b's sign bit '
+	movel	d1,d3		|
+	beq	Ldivsf$b$0	| branch if b is zero
+	cmpl	d6,d0		| is a big?
+	bhi	Ldivsf$inop	| if a is NaN return NaN
+	beq	Ldivsf$inf	| if a is INIFINITY we have to check b
+	cmpl	d6,d1		| now compare b with INFINITY 
+	bhi	Ldivsf$inop	| if b is NaN return NaN
+	beq	Ldivsf$underflow
+| Here we have both numbers finite and nonzero (and with no sign bit).
+| Now we get the exponents into d2 and d3 and normalize the numbers to
+| ensure that the ratio of the fractions is close to 1. We do this by
+| making sure that bit #FLT_MANT_DIG-1 (hidden bit) is set.
+	andl	d6,d2		| and isolate exponent in d2
+	beq	Ldivsf$a$den	| if exponent is zero we have a denormalized
+	andl	d5,d0		| and isolate fraction
+	orl	d4,d0		| and put hidden bit back
+	swap	d2		| I like exponents in the first byte
+	lsrw	#7,d2		| 
+Ldivsf$1:			| 
+	andl	d6,d3		|
+	beq	Ldivsf$b$den	|
+	andl	d5,d1		|
+	orl	d4,d1		|
+	swap	d3		|
+	lsrw	#7,d3		|
+Ldivsf$2:			|
+	subw	d3,d2		| substract exponents
+ 	addw	#F_BIAS,d2	| and add bias
+ 
+| We are now ready to do the division. We have prepared things in such a way
+| that the ratio of the fractions will be less than 2 but greater than 1/2.
+| At this point the registers in use are:
+| d0	holds a (first operand, bit FLT_MANT_DIG=0, bit FLT_MANT_DIG-1=1)
+| d1	holds b (second operand, bit FLT_MANT_DIG=1)
+| d2	holds the difference of the exponents, corrected by the bias
+| d7	holds the sign of the ratio
+| d4, d5, d6 hold some constants
+	movel	d7,a0		| d6-d7 will hold the ratio of the fractions
+	movel	#0,d6		| 
+	movel	d6,d7
+
+	movew	#FLT_MANT_DIG+1,d3
+1:	cmpl	d0,d1		| is a < b?
+	bhi	2f		|
+	bset	d3,d6		| set a bit in d6
+	subl	d1,d0		| if a >= b  a <-- a-b
+	beq	3f		| if a is zero, exit
+2:	addl	d0,d0		| multiply a by 2
+	dbra	d3,1b
+
+| Now we keep going to set the sticky bit ...
+	movew	#FLT_MANT_DIG,d3
+1:	cmpl	d0,d1
+	ble	2f
+	addl	d0,d0
+	dbra	d3,1b
+	movel	#0,d1
+	bra	3f
+2:	movel	#0,d1
+	subw	#FLT_MANT_DIG,d3
+	addw	#31,d3
+	bset	d3,d1
+3:
+	movel	d6,d0		| put the ratio in d0-d1
+	movel	a0,d7		| get sign back
+
+| Because of the normalization we did before we are guaranteed that 
+| d0 is smaller than 2^26 but larger than 2^24. Thus bit 26 is not set,
+| bit 25 could be set, and if it is not set then bit 24 is necessarily set.
+	btst	#FLT_MANT_DIG+1,d0		
+	beq	1f              | if it is not set, then bit 24 is set
+	lsrl	#1,d0           |
+	addw	#1,d2           |
+1:
+| Now round, check for over- and underflow, and exit.
+	movew	#DIVIDE,d5
+	bra	Lround$exit
+
+Ldivsf$inop:
+	movew	#DIVIDE,d5
+	bra	Lf$inop
+
+Ldivsf$overflow:
+	movew	#DIVIDE,d5
+	bra	Lf$overflow
+
+Ldivsf$underflow:
+	movew	#DIVIDE,d5
+	bra	Lf$underflow
+
+Ldivsf$a$0:
+	movew	#DIVIDE,d5
+| If a is zero check to see whether b is zero also. In that case return
+| NaN; then check if b is NaN, and return NaN also in that case. Else
+| return zero.
+	andl	#0x7fffffff,d1	| clear sign bit and test b
+	beq	Lf$inop		| if b is also zero return NaN
+	cmpl	#INFINITY,d1	| check for NaN
+	bhi	Lf$inop		| 
+	movel	#0,d0		| else return zero
+	lea	SYM (_fpCCR),a0	|
+	movew	#0,a0@		|
+	moveml	sp@+,d2-d7	| 
+	unlk	a6		| 
+	rts			| 
+	
+Ldivsf$b$0:
+	movew	#DIVIDE,d5
+| If we got here a is not zero. Check if a is NaN; in that case return NaN,
+| else return +/-INFINITY. Remember that a is in d0 with the sign bit 
+| cleared already.
+	cmpl	#INFINITY,d0	| compare d0 with INFINITY
+	bhi	Lf$inop		| if larger it is NaN
+	bra	Lf$div$0	| else signal DIVIDE_BY_ZERO
+
+Ldivsf$inf:
+	movew	#DIVIDE,d5
+| If a is INFINITY we have to check b
+	cmpl	#INFINITY,d1	| compare b with INFINITY 
+	bge	Lf$inop		| if b is NaN or INFINITY return NaN
+	bra	Lf$overflow	| else return overflow
+
+| If a number is denormalized we put an exponent of 1 but do not put the 
+| bit back into the fraction.
+Ldivsf$a$den:
+	movel	#1,d2
+	andl	d5,d0
+1:	addl	d0,d0		| shift a left until bit FLT_MANT_DIG-1 is set
+	subw	#1,d2		| and adjust exponent
+	btst	#FLT_MANT_DIG-1,d0
+	bne	Ldivsf$1
+	bra	1b
+
+Ldivsf$b$den:
+	movel	#1,d3
+	andl	d5,d1
+1:	addl	d1,d1		| shift b left until bit FLT_MANT_DIG is set
+	subw	#1,d3		| and adjust exponent
+	btst	#FLT_MANT_DIG-1,d1
+	bne	Ldivsf$2
+	bra	1b
+
+Lround$exit:
+| This is a common exit point for __mulsf3 and __divsf3. 
+
+| First check for underlow in the exponent:
+	cmpw	#-FLT_MANT_DIG-1,d2		
+	blt	Lf$underflow	
+| It could happen that the exponent is less than 1, in which case the 
+| number is denormalized. In this case we shift right and adjust the 
+| exponent until it becomes 1 or the fraction is zero (in the latter case 
+| we signal underflow and return zero).
+	movel	#0,d6		| d6 is used temporarily
+	cmpw	#1,d2		| if the exponent is less than 1 we 
+	bge	2f		| have to shift right (denormalize)
+1:	addw	#1,d2		| adjust the exponent
+	lsrl	#1,d0		| shift right once 
+	roxrl	#1,d1		|
+	roxrl	#1,d6		| d6 collect bits we would lose otherwise
+	cmpw	#1,d2		| is the exponent 1 already?
+	beq	2f		| if not loop back
+	bra	1b              |
+	bra	Lf$underflow	| safety check, shouldn't execute '
+2:	orl	d6,d1		| this is a trick so we don't lose  '
+				| the extra bits which were flushed right
+| Now call the rounding routine (which takes care of denormalized numbers):
+	lea	Lround$0,a0	| to return from rounding routine
+	lea	SYM (_fpCCR),a1	| check the rounding mode
+	movew	a1@(6),d6	| rounding mode in d6
+	beq	Lround$to$nearest
+	cmpw	#ROUND_TO_PLUS,d6
+	bhi	Lround$to$minus
+	blt	Lround$to$zero
+	bra	Lround$to$plus
+Lround$0:
+| Here we have a correctly rounded result (either normalized or denormalized).
+
+| Here we should have either a normalized number or a denormalized one, and
+| the exponent is necessarily larger or equal to 1 (so we don't have to  '
+| check again for underflow!). We have to check for overflow or for a 
+| denormalized number (which also signals underflow).
+| Check for overflow (i.e., exponent >= 255).
+	cmpw	#0x00ff,d2
+	bge	Lf$overflow
+| Now check for a denormalized number (exponent==0).
+	movew	d2,d2
+	beq	Lf$den
+1:
+| Put back the exponents and sign and return.
+	lslw	#7,d2		| exponent back to fourth byte
+	bclr	#FLT_MANT_DIG-1,d0
+	swap	d0		| and put back exponent
+	orw	d2,d0		| 
+	swap	d0		|
+	orl	d7,d0		| and sign also
+
+	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7
+	unlk	a6
+	rts
+
+|=============================================================================
+|                             __negsf2
+|=============================================================================
+
+| This is trivial and could be shorter if we didn't bother checking for NaN '
+| and +/-INFINITY.
+
+| float __negsf2(float);
+SYM (__negsf2):
+	link	a6,#0
+	moveml	d2-d7,sp@-
+	movew	#NEGATE,d5
+	movel	a6@(8),d0	| get number to negate in d0
+	bchg	#31,d0		| negate
+	movel	d0,d1		| make a positive copy
+	bclr	#31,d1		|
+	tstl	d1		| check for zero
+	beq	2f		| if zero (either sign) return +zero
+	cmpl	#INFINITY,d1	| compare to +INFINITY
+	blt	1f		|
+	bhi	Lf$inop		| if larger (fraction not zero) is NaN
+	movel	d0,d7		| else get sign and return INFINITY
+	andl	#0x80000000,d7
+	bra	Lf$infty		
+1:	lea	SYM (_fpCCR),a0
+	movew	#0,a0@
+	moveml	sp@+,d2-d7
+	unlk	a6
+	rts
+2:	bclr	#31,d0
+	bra	1b
+
+|=============================================================================
+|                             __cmpsf2
+|=============================================================================
+
+GREATER =  1
+LESS    = -1
+EQUAL   =  0
+
+| int __cmpsf2(float, float);
+SYM (__cmpsf2):
+	link	a6,#0
+	moveml	d2-d7,sp@- 	| save registers
+	movew	#COMPARE,d5
+	movel	a6@(8),d0	| get first operand
+	movel	a6@(12),d1	| get second operand
+| Check if either is NaN, and in that case return garbage and signal
+| INVALID_OPERATION. Check also if either is zero, and clear the signs
+| if necessary.
+	movel	d0,d6
+	andl	#0x7fffffff,d0
+	beq	Lcmpsf$a$0
+	cmpl	#0x7f800000,d0
+	bhi	Lf$inop
+Lcmpsf$1:
+	movel	d1,d7
+	andl	#0x7fffffff,d1
+	beq	Lcmpsf$b$0
+	cmpl	#0x7f800000,d1
+	bhi	Lf$inop
+Lcmpsf$2:
+| Check the signs
+	eorl	d6,d7
+	bpl	1f
+| If the signs are not equal check if a >= 0
+	tstl	d6
+	bpl	Lcmpsf$a$gt$b	| if (a >= 0 && b < 0) => a > b
+	bmi	Lcmpsf$b$gt$a	| if (a < 0 && b >= 0) => a < b
+1:
+| If the signs are equal check for < 0
+	tstl	d6
+	bpl	1f
+| If both are negative exchange them
+	exg	d0,d1
+1:
+| Now that they are positive we just compare them as longs (does this also
+| work for denormalized numbers?).
+	cmpl	d0,d1
+	bhi	Lcmpsf$b$gt$a	| |b| > |a|
+	bne	Lcmpsf$a$gt$b	| |b| < |a|
+| If we got here a == b.
+	movel	#EQUAL,d0
+	moveml	sp@+,d2-d7 	| put back the registers
+	unlk	a6
+	rts
+Lcmpsf$a$gt$b:
+	movel	#GREATER,d0
+	moveml	sp@+,d2-d7 	| put back the registers
+	unlk	a6
+	rts
+Lcmpsf$b$gt$a:
+	movel	#LESS,d0
+	moveml	sp@+,d2-d7 	| put back the registers
+	unlk	a6
+	rts
+
+Lcmpsf$a$0:	
+	bclr	#31,d6
+	bra	Lcmpsf$1
+Lcmpsf$b$0:
+	bclr	#31,d7
+	bra	Lcmpsf$2
+
+|=============================================================================
+|                           rounding routines
+|=============================================================================
+
+| The rounding routines expect the number to be normalized in registers
+| d0-d1, with the exponent in register d2. They assume that the 
+| exponent is larger or equal to 1. They return a properly normalized number
+| if possible, and a denormalized number otherwise. The exponent is returned
+| in d2.
+
+Lround$to$nearest:
+| We now normalize as suggested by D. Knuth ("Seminumerical Algorithms"):
+| Here we assume that the exponent is not too small (this should be checked
+| before entering the rounding routine), but the number could be denormalized.
+
+| Check for denormalized numbers:
+1:	btst	#FLT_MANT_DIG,d0
+	bne	2f		| if set the number is normalized
+| Normalize shifting left until bit #FLT_MANT_DIG is set or the exponent 
+| is one (remember that a denormalized number corresponds to an 
+| exponent of -F_BIAS+1).
+	cmpw	#1,d2		| remember that the exponent is at least one
+ 	beq	2f		| an exponent of one means denormalized
+	addl	d1,d1		| else shift and adjust the exponent
+	addxl	d0,d0		|
+	dbra	d2,1b		|
+2:
+| Now round: we do it as follows: after the shifting we can write the
+| fraction part as f + delta, where 1 < f < 2^25, and 0 <= delta <= 2.
+| If delta < 1, do nothing. If delta > 1, add 1 to f. 
+| If delta == 1, we make sure the rounded number will be even (odd?) 
+| (after shifting).
+	btst	#0,d0		| is delta < 1?
+	beq	2f		| if so, do not do anything
+	tstl	d1		| is delta == 1?
+	bne	1f		| if so round to even
+	movel	d0,d1		| 
+	andl	#2,d1		| bit 1 is the last significant bit
+	addl	d1,d0		| 
+	bra	2f		| 
+1:	movel	#1,d1		| else add 1 
+	addl	d1,d0		|
+| Shift right once (because we used bit #FLT_MANT_DIG!).
+2:	lsrl	#1,d0		
+| Now check again bit #FLT_MANT_DIG (rounding could have produced a
+| 'fraction overflow' ...).
+	btst	#FLT_MANT_DIG,d0	
+	beq	1f
+	lsrl	#1,d0
+	addw	#1,d2
+1:
+| If bit #FLT_MANT_DIG-1 is clear we have a denormalized number, so we 
+| have to put the exponent to zero and return a denormalized number.
+	btst	#FLT_MANT_DIG-1,d0
+	beq	1f
+	jmp	a0@
+1:	movel	#0,d2
+	jmp	a0@
+
+Lround$to$zero:
+Lround$to$plus:
+Lround$to$minus:
+	jmp	a0@
+#endif /* L_float */
+
+| gcc expects the routines __eqdf2, __nedf2, __gtdf2, __gedf2,
+| __ledf2, __ltdf2 to all return the same value as a direct call to
+| __cmpdf2 would.  In this implementation, each of these routines
+| simply calls __cmpdf2.  It would be more efficient to give the
+| __cmpdf2 routine several names, but separating them out will make it
+| easier to write efficient versions of these routines someday.
+
+#ifdef  L_eqdf2
+LL0:
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	4
+	LS18	=	128
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__eqdf2)
+SYM (__eqdf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(20),sp@-
+	movl	a6@(16),sp@-
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpdf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_eqdf2 */
+
+#ifdef  L_nedf2
+LL0:
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__nedf2)
+SYM (__nedf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(20),sp@-
+	movl	a6@(16),sp@-
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpdf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_nedf2 */
+
+#ifdef  L_gtdf2
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__gtdf2)
+SYM (__gtdf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(20),sp@-
+	movl	a6@(16),sp@-
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpdf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_gtdf2 */
+
+#ifdef  L_gedf2
+LL0:
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__gedf2)
+SYM (__gedf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(20),sp@-
+	movl	a6@(16),sp@-
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpdf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_gedf2 */
+
+#ifdef  L_ltdf2
+LL0:
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__ltdf2)
+SYM (__ltdf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(20),sp@-
+	movl	a6@(16),sp@-
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpdf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_ltdf2 */
+
+#ifdef  L_ledf2
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__ledf2)
+SYM (__ledf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(20),sp@-
+	movl	a6@(16),sp@-
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpdf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_ledf2 */
+
+| The comments above about __eqdf2, et. al., also apply to __eqsf2,
+| et. al., except that the latter call __cmpsf2 rather than __cmpdf2.
+
+#ifdef  L_eqsf2
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	4
+	LS18	=	128
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__eqsf2)
+SYM (__eqsf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpsf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_eqsf2 */
+
+#ifdef  L_nesf2
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__nesf2)
+SYM (__nesf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpsf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_nesf2 */
+
+#ifdef  L_gtsf2
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__gtsf2)
+SYM (__gtsf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpsf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_gtsf2 */
+
+#ifdef  L_gesf2
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__gesf2)
+SYM (__gesf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpsf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_gesf2 */
+
+#ifdef  L_ltsf2
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__ltsf2)
+SYM (__ltsf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpsf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_ltsf2 */
+
+#ifdef  L_lesf2
+	.text
+	.proc
+|#PROC# 04
+	LF18	=	8
+	LS18	=	132
+	LFF18	=	0
+	LSS18	=	0
+	LV18	=	0
+	.text
+	.globl	SYM (__lesf2)
+SYM (__lesf2):
+|#PROLOGUE# 0
+	link	a6,#0
+|#PROLOGUE# 1
+	movl	a6@(12),sp@-
+	movl	a6@(8),sp@-
+	jbsr	SYM (__cmpsf2)
+|#PROLOGUE# 2
+	unlk	a6
+|#PROLOGUE# 3
+	rts
+#endif /* L_lesf2 */
+
diff --git a/gcc/config/m68k/m68k-aout.h b/gcc/config/m68k/m68k-aout.h
new file mode 100644
index 00000000000..8e3b66181f2
--- /dev/null
+++ b/gcc/config/m68k/m68k-aout.h
@@ -0,0 +1,30 @@
+/* Definitions of target machine for GNU compiler.  "naked" 68020,
+   a.out object files and debugging, version.
+   Copyright (C) 1994 Free Software Foundation, Inc.
+
+This file is part of GNU CC.
+
+GNU CC is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+
+GNU CC is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GNU CC; see the file COPYING.  If not, write to
+the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
+
+/* This comment is here to see if it will keep Sun's cpp from dying.  */
+
+#include "m68k/m68k-none.h"
+#include "m68k/m68kemb.h"
+#include "aoutos.h"
+
+#define DBX_DEBUGGING_INFO
+#undef SDB_DEBUGGING_INFO
+
+/* end of m68k-aout.h */
diff --git a/gcc/config/m68k/m68k-none.h b/gcc/config/m68k/m68k-none.h
new file mode 100644
index 00000000000..5e584681ce0
--- /dev/null
+++ b/gcc/config/m68k/m68k-none.h
@@ -0,0 +1,94 @@
+/* Definitions of target machine for GNU compiler.  "naked" 68020.
+   Copyright (C) 1994 Free Software Foundation, Inc.
+
+This file is part of GNU CC.
+
+GNU CC is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+
+GNU CC is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GNU CC; see the file COPYING.  If not, write to
+the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
+
+#include "m68k/m68k.h"
+
+/* See m68k.h.  7 means 68020 with 68881.  */
+
+#ifndef TARGET_DEFAULT
+#define TARGET_DEFAULT 7
+#endif
+
+/* Names to predefine in the preprocessor for this target machine.
+   Always define mc68000.  Other definitions depend on switches given
+   to the compiler:
+
+   -m68000: define nothing else
+   default, -m68020, -mc68020: define mc68020
+   -m68030: define mc68030
+   -m68040: define mc68040
+   -m68020-40: define mc68020 mc68030 mc68040
+   -m68302: define mc68302
+   -m68332: define mc68332
+   */
+
+#ifndef CPP_PREDEFINES
+#define CPP_PREDEFINES "-Dmc68000"
+#endif
+
+#ifndef CPP_SPEC
+
+#if TARGET_DEFAULT & 02
+
+/* -m68881 is the default */
+#define CPP_SPEC \
+"%{!mc68000:%{!m68000:%{!m68332:%{!msoft-float:%{mfpa:-D__HAVE_FPA__ }%{!mfpa:-D__HAVE_68881__ }}}}}\
+%{!ansi:%{m68010:-Dmc68010 }%{m68020:-Dmc68020 }%{mc68020:-Dmc68020 }%{m68030:-Dmc68030 }%{m68040:-Dmc68040 }%{m68020-40:-Dmc68020 -Dmc68030 -Dmc68040 }%{m68302:-Dmc68302 }%{m68332:-Dmc68332 }%{!mc68000:%{!m68000:%{!m68010:%{!mc68020:%{!m68020:%{!m68030:%{!m68040:%{!m68020-40:%{!m68302:%{!m68332:-Dmc68020 }}}}}}}}}}}\
+%{m68010:-D__mc68010__ -D__mc68010 }%{m68020:-D__mc68020__ -D__mc68020 }%{mc68020:-D__mc68020__ -D__mc68020 }%{m68030:-D__mc68030__ -D__mc68030 }%{m68040:-D__mc68040__ -D__mc68040 }%{m68020-40:-D__mc68020__ -D__mc68030__ -D__mc68040__ -D__mc68020 -D__mc68030 -D__mc68040 }%{m68302:-D__mc68302__ -D__mc68302 }%{m68332:-D__mc68332__ -D__mc68332 }%{!mc68000:%{!m68000:%{!m68010:%{!mc68020:%{!m68020:%{!m68030:%{!m68040:%{!m68020-40:%{!m68302:%{!m68332:-D__mc68020__ -D__mc68020 }}}}}}}}}}"
+
+#else
+#if TARGET_DEFAULT & 0100
+
+/* -mfpa is the default */
+#define CPP_SPEC \
+"%{!msoft-float:%{m68881:-D__HAVE_68881__ }%{!m68881:-D__HAVE_FPA__ }}\
+%{!ansi:%{m68010:-Dmc68010 }%{m68020:-Dmc68020 }%{mc68020:-Dmc68020 }%{m68030:-Dmc68030 }%{m68040:-Dmc68040 }%{m68020-40:-Dmc68020 -Dmc68030 -Dmc68040 }%{m68302:-Dmc68302 }%{m68332:-Dmc68332 }%{!mc68000:%{!m68000:%{!m68010:%{!mc68020:%{!m68020:%{!m68030:%{!m68040:%{!m68020-40:%{!m68302:%{!m68332:-Dmc68020 }}}}}}}}}}}\
+%{m68010:-D__mc68010__ -D__mc68010 }%{m68020:-D__mc68020__ -D__mc68020 }%{mc68020:-D__mc68020__ -D__mc68020 }%{m68030:-D__mc68030__ -D__mc68030 }%{m68040:-D__mc68040__ -D__mc68040 }%{m68020-40:-D__mc68020__ -D__mc68030__ -D__mc68040__ -D__mc68020 -D__mc68030 -D__mc68040 }%{m68302:-D__mc68302__ -D__mc68302 }%{m68332:-D__mc68332__ -D__mc68332 }%{!mc68000:%{!m68000:%{!m68010:%{!mc68020:%{!m68020:%{!m68030:%{!m68040:%{!m68020-40:%{!m68302:%{!m68332:-D__mc68020__ -D__mc68020 }}}}}}}}}}"
+
+#else
+
+/* -msoft-float is the default */
+#define CPP_SPEC \
+"%{m68881:-D__HAVE_68881__ }%{mfpa:-D__HAVE_FPA__ }\
+%{!ansi:%{m68010:-Dmc68010 }%{m68020:-Dmc68020 }%{mc68020:-Dmc68020 }%{m68030:-Dmc68030 }%{m68040:-Dmc68040 }%{m68020-40:-Dmc68020 -Dmc68030 -Dmc68040 }%{m68302:-Dmc68302 }%{m68332:-Dmc68332 }%{!mc68000:%{!m68000:%{!m68010:%{!mc68020:%{!m68020:%{!m68030:%{!m68040:%{!m68020-40:%{!m68302:%{!m68332:-Dmc68020 }}}}}}}}}}}\
+%{m68010:-D__mc68010__ -D__mc68010 }%{m68020:-D__mc68020__ -D__mc68020 }%{mc68020:-D__mc68020__ -D__mc68020 }%{m68030:-D__mc68030__ -D__mc68030 }%{m68040:-D__mc68040__ -D__mc68040 }%{m68020-40:-D__mc68020__ -D__mc68030__ -D__mc68040__ -D__mc68020 -D__mc68030 -D__mc68040 }%{m68302:-D__mc68302__ -D__mc68302 }%{m68332:-D__mc68332__ -D__mc68332 }%{!mc68000:%{!m68000:%{!m68010:%{!mc68020:%{!m68020:%{!m68030:%{!m68040:%{!m68020-40:%{!m68302:%{!m68332:-D__mc68020__ -D__mc68020 }}}}}}}}}}"
+
+#endif
+#endif
+
+#endif
+
+/* Pass flags to gas indicating which type of processor we have.  */
+
+#ifndef ASM_SPEC
+
+#define ASM_SPEC \
+"%{m68851}%{mno-68851}%{m68881}%{mno-68881}%{msoft-float:-mno-68881 }\
+%{m68000}%{mc68000}%{m68010}%{m68020}%{mc68020}%{m68030}%{m68040}%{m68020-40:-mc68040}%{m68302}%{m68332}%{!m68000:%{!mc68000:%{!m68010:%{!mc68020:%{!m68020:%{!m68030:%{!m68040:%{!m68020-40:%{!m68302:%{!m68332:-mc68020}}}}}}}}}}"
+
+#endif
+
+#ifndef CC1_SPEC
+
+#define CC1_SPEC \
+ "%{m68000:%{!m68881:-msoft-float }}%{m68302:-m68000}%{m68332:-m68020 -mnobitfield %{!m68881:-msoft-float}}%{!m68000:%{!mc68000:%{!m68010:%{!mc68020:%{!m68020:%{!m68030:%{!m68040:%{!m68020-40:%{!m68302:%{!m68332:-m68020}}}}}}}}}}"
+
+#endif
+
+/* end of m68k-none.h */
diff --git a/gcc/config/m68k/m68kemb.h b/gcc/config/m68k/m68kemb.h
new file mode 100644
index 00000000000..7b7521ab2ee
--- /dev/null
+++ b/gcc/config/m68k/m68kemb.h
@@ -0,0 +1,39 @@
+/* Definitions of target machine for GNU compiler.  "embedded" 68XXX.
+   This is meant to be included after m68k.h.
+   Copyright (C) 1994 Free Software Foundation, Inc.  */
+
+#define PTRDIFF_TYPE "long int"
+#define SIZE_TYPE "long unsigned int"
+
+/* In order for bitfields to work on a 68000, or with -mnobitfield, we must
+   define either PCC_BITFIELD_TYPE_MATTERS or STRUCTURE_SIZE_BOUNDARY.
+   Defining STRUCTURE_SIZE_BOUNDARY results in structure packing problems,
+   so we define PCC_BITFIELD_TYPE_MATTERS.  */
+#define PCC_BITFIELD_TYPE_MATTERS 1
+
+/* Undef PCC_STATIC_STRUCT_RETURN so that we get a re-entrant calling
+   convention.  */
+#undef PCC_STATIC_STRUCT_RETURN
+
+/* Don't default to pcc-struct-return, so that we can return small structures
+   and unions in registers, which is slightly more efficient.  */
+#define DEFAULT_PCC_STRUCT_RETURN 0
+
+/* Return floating point values in a fp register.  This make fp code a
+   little bit faster.  It also makes -msoft-float code incompatible with
+   -m68881 code, so people have to be careful not to mix the two.  */
+#undef FUNCTION_VALUE
+#define FUNCTION_VALUE(VALTYPE,FUNC) LIBCALL_VALUE (TYPE_MODE (VALTYPE))
+
+#undef LIBCALL_VALUE
+#define LIBCALL_VALUE(MODE)                                                \
+ gen_rtx (REG, (MODE),                                                     \
+          ((TARGET_68881                                                   \
+            && ((MODE) == SFmode || (MODE) == DFmode || (MODE) == XFmode)) \
+           ? 16 : 0))
+
+#undef FUNCTION_VALUE_REGNO_P
+#define FUNCTION_VALUE_REGNO_P(N) ((N) == 0 || (TARGET_68881 && (N) == 16))
+
+#undef NEEDS_UNTYPED_CALL
+#define NEEDS_UNTYPED_CALL 1
diff --git a/gcc/config/m68k/t-m68kbare b/gcc/config/m68k/t-m68kbare
new file mode 100644
index 00000000000..732ad4d1632
--- /dev/null
+++ b/gcc/config/m68k/t-m68kbare
@@ -0,0 +1,23 @@
+CROSS_LIBGCC1 = libgcc1-asm.a
+LIB1ASMSRC = m68k/lb1sf68.asm
+LIB1ASMFUNCS = _mulsi3 _udivsi3 _umulsi3 _divsi3 _umodsi3 _modsi3 \
+   _double _float _floatex \
+   _eqdf2 _nedf2 _gtdf2 _gedf2 _ltdf2 _ledf2 \
+   _eqsf2 _nesf2 _gtsf2 _gesf2 _ltsf2 _lesf2
+
+# These are really part of libgcc1, but this will cause them to be
+# built correctly, so...
+LIB2FUNCS_EXTRA = fpgnulib.c xfgnulib.c
+
+fpgnulib.c: $(srcdir)/config/m68k/fpgnulib.c
+	cp $(srcdir)/config/m68k/fpgnulib.c fpgnulib.c
+xfgnulib.c: $(srcdir)/config/m68k/fpgnulib.c
+	echo '#define EXTFLOAT' > xfgnulib.c
+	cat $(srcdir)/config/m68k/fpgnulib.c >> xfgnulib.c
+
+MULTILIB_OPTIONS=m68000/m68020 m68881/msoft-float
+MULTILIB_DIRNAMES=
+MULTILIB_MATCHES=m68020=m68040 m68000=mc68000 m68020=mc68020
+
+LIBGCC = stmp-multilib
+INSTALL_LIBGCC = install-multilib
diff --git a/gcc/config/m88k/m88k-aout.h b/gcc/config/m88k/m88k-aout.h
new file mode 100644
index 00000000000..e3cc87bca63
--- /dev/null
+++ b/gcc/config/m88k/m88k-aout.h
@@ -0,0 +1,31 @@
+/* Definitions for "naked" Motorola 88k using a.out object format files
+   and stabs debugging info.
+
+   Copyright (C) 1994 Free Software Foundation, Inc.
+
+This file is part of GNU CC.
+
+GNU CC is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+
+GNU CC is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GNU CC; see the file COPYING.  If not, write to
+the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
+
+#undef SDB_DEBUGGING_INFO
+#define DBX_DEBUGGING_INFO
+
+#include "m88k/m88k.h"
+#include "aoutos.h"
+
+#undef CPP_PREDEFINES
+#define CPP_PREDEFINES "-Dm88000 -Dm88k"
+
+/* end of m88k-aout.h */
diff --git a/gcc/config/m88k/m88k-coff.h b/gcc/config/m88k/m88k-coff.h
new file mode 100644
index 00000000000..b86582a22b7
--- /dev/null
+++ b/gcc/config/m88k/m88k-coff.h
@@ -0,0 +1,33 @@
+/* Definitions for "naked" Motorola 88k using coff object format files
+   and coff debugging info.
+
+   Copyright (C) 1994 Free Software Foundation, Inc.
+
+This file is part of GNU CC.
+
+GNU CC is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+
+GNU CC is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GNU CC; see the file COPYING.  If not, write to
+the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
+
+#include "m88k/m88k.h"
+
+#undef CPP_PREDEFINES
+#define CPP_PREDEFINES "-Dm88000 -Dm88k"
+
+#undef DBX_DEBUGGING_INFO
+#define SDB_DEBUGGING_INFO
+
+#undef PREFERRED_DEBUGGING_TYPE
+#define PREFERRED_DEBUGGING_TYPE SDB_DEBUG
+
+/* end of m88k-coff.h */
diff --git a/gcc/config/m88k/t-bug b/gcc/config/m88k/t-bug
new file mode 100644
index 00000000000..a5e71ddf3fc
--- /dev/null
+++ b/gcc/config/m88k/t-bug
@@ -0,0 +1,12 @@
+# Specify how to create the *.asm files
+
+MOVE_ASM = moveHI15x.asm moveQI16x.asm moveSI46x.asm moveSI64n.asm \
+	   moveHI48x.asm moveSI45x.asm moveSI47x.asm moveSI96x.asm \
+	   moveDI96x.asm
+
+$(MOVE_ASM): $(srcdir)/config/m88k/m88k-move.sh
+	$(srcdir)/config/m88k/m88k-move.sh
+
+LIB2FUNCS_EXTRA = $(MOVE_ASM)
+LIBGCC1 = libgcc1.null
+CROSS_LIBGCC1 = libgcc1.null