ack subdir for combined ack/gcc library build

lib/ack/Makefile

@@ -0,0 +1,12 @@
+
+SUBDIRS = \
+	float \
+	fphook \
+	i386 \
+	libm2 \
+	libp \
+	liby \
+	math \
+	rts \
+
+include ../Makefile.inc

lib/ack/Makefile.ack

@@ -0,0 +1,11 @@
+
+all:
+	cd float && make
+	cd fphook && make
+	cd `arch` && make
+	cd libm2 && make
+	cd libp && make
+	cd liby && make
+	cd math && make
+	cd rts && make
+

lib/ack/float/FP.compile

@@ -0,0 +1,19 @@
+#!/bin/sh
+#							Author: Kees J. Bot
+# Compile one soft FP source file.
+# (These files shouldn't be optimized normally, but the 16-bit C compiler
+# only optimizes scratch register allocation a bit with -O.  To the 32-bit
+# compiler -O is a no-op.)
+
+case $#:$1 in
+1:*.fc)	;;
+*)	echo "$0: $1: not a FC file" >&2; exit 1
+esac
+
+base="`basename "$1" .fc`"
+trap 'rm -f tmp.c tmp.s"; exit 1' 2
+
+cp "$1" tmp.c &&
+cc -O -I. -D_MINIX -D_POSIX_SOURCE -S tmp.c &&
+sed -f FP.script tmp.s > "$base.s" &&
+rm tmp.c tmp.s

lib/ack/float/FP.script

@@ -0,0 +1,39 @@
+s/_adf4/.adf4/
+s/_adf8/.adf8/
+s/_cff4/.cff4/
+s/_cff8/.cff8/
+s/_cfi/.cfi/
+s/_cfu/.cfu/
+s/_cif4/.cif4/
+s/_cif8/.cif8/
+s/_cmf4/.cmf4/
+s/_cmf8/.cmf8/
+s/_cuf4/.cuf4/
+s/_cuf8/.cuf8/
+s/_dvf4/.dvf4/
+s/_dvf8/.dvf8/
+s/_fef4/.fef4/
+s/_fef8/.fef8/
+s/_fif4/.fif4/
+s/_fif8/.fif8/
+s/_mlf4/.mlf4/
+s/_mlf8/.mlf8/
+s/_ngf4/.ngf4/
+s/_ngf8/.ngf8/
+s/_sbf4/.sbf4/
+s/_sbf8/.sbf8/
+s/_zrf4/.zrf4/
+s/_zrf8/.zrf8/
+s/_add_ext/.add_ext/
+s/_div_ext/.div_ext/
+s/_mul_ext/.mul_ext/
+s/_nrm_ext/.nrm_ext/
+s/_sft_ext/.sft_ext/
+s/_sub_ext/.sub_ext/
+s/_zrf_ext/.zrf_ext/
+s/_compact/.compact/
+s/_extend/.extend/
+s/_b64_add/.b64_add/
+s/_b64_sft/.b64_sft/
+s/_b64_rsft/.b64_rsft/
+s/_b64_lsft/.b64_lsft/

lib/ack/float/FP_bias.h

@@ -0,0 +1,28 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	include file for floating point package
+*/
+
+		/*	FLOAT FORMAT EXPONENT BIAS	*/
+
+#define	SGL_BIAS	 127	/* excess  128 notation used	*/
+#define	DBL_BIAS	1023	/* excess 1024 notation used	*/
+#define	EXT_BIAS	   0	/* 2s-complement notation used	*/
+				/* this is possible because the	*/
+				/* sign is in a seperate word	*/
+		
+		/*	VARIOUS MAX AND MIN VALUES	*/
+		/*	1) FOR THE DIFFERENT FORMATS	*/
+
+#define	SGL_MAX		   254	/*	standard definition	*/
+#define	SGL_MIN		     1	/*	standard definition	*/
+#define	DBL_MAX		  2046	/*	standard definition	*/
+#define	DBL_MIN		     1	/*	standard definition	*/
+#define EXT_MAX		 16383	/*	standard minimum	*/
+#define EXT_MIN		-16382	/*	standard minimum	*/

lib/ack/float/FP_shift.h

@@ -0,0 +1,49 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	include file for floating point package
+*/
+
+# define	CARRYBIT	0x80000000L
+# define	NORMBIT		0x80000000L
+# define	EXP_STORE	16
+
+
+				/* parameters for Single Precision */
+#define SGL_EXPSHIFT	7
+#define SGL_M1LEFT	8
+#define SGL_ZERO	0xffffff80L
+#define SGL_EXACT	0xff
+#define SGL_RUNPACK	SGL_M1LEFT
+
+#define SGL_ROUNDUP	0x80
+#define	SGL_CARRYOUT	0x01000000L
+#define	SGL_MASK	0x007fffffL
+
+				/* parameters for Double Precision */
+				/* used in extend.c */
+
+#define DBL_EXPSHIFT	4
+
+#define DBL_M1LEFT	11
+
+#define	DBL_RPACK	(32-DBL_M1LEFT)
+#define	DBL_LPACK	DBL_M1LEFT
+
+				/* used in compact.c */
+
+#define DBL_ZERO	0xfffffd00L
+
+#define DBL_EXACT	0x7ff
+
+#define DBL_RUNPACK	DBL_M1LEFT
+#define DBL_LUNPACK	(32-DBL_RUNPACK)
+
+#define DBL_ROUNDUP	0x400
+#define	DBL_CARRYOUT	0x00200000L
+#define	DBL_MASK	0x000fffffL

lib/ack/float/FP_trap.h

@@ -0,0 +1,22 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	include file for floating point package
+*/
+
+			/*	EM TRAPS	*/
+
+#define	EIOVFL	3	/* Integer  Overflow		*/
+#define	EFOVFL	4	/* Floating Overflow		*/
+#define	EFUNFL	5	/* Floating Underflow		*/
+#define	EIDIVZ	6	/* Integer  Divide by 0		*/
+#define	EFDIVZ	7	/* Floating Divide by 0.0	*/
+#define	EIUND	8	/* Integer  Undefined Number	*/
+#define	EFUND	9	/* Floating Undefined Number	*/
+#define	ECONV	10	/* Conversion Error		*/
+# define trap(x) _fptrp(x)

lib/ack/float/FP_types.h

@@ -0,0 +1,113 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/********************************************************/
+/*
+	Type definitions for C Floating Point Package
+	include file for floating point package
+*/
+/********************************************************/
+/*
+	THESE STRUCTURES ARE USED TO ADDRESS THE INDIVIDUAL
+	PARTS OF THE FLOATING POINT NUMBER REPRESENTATIONS.
+
+	THREE STRUCTURES ARE DEFINED:
+		SINGLE:	single precision floating format
+		DOUBLE:	double precision floating format
+		EXTEND:	double precision extended format
+*/
+/********************************************************/
+
+#ifndef __FPTYPES
+#define __FPTYPES
+
+typedef	struct	{
+	unsigned long	h_32;	/* higher 32 bits of 64 */
+	unsigned long	l_32;	/* lower  32 bits of 64 */
+}	B64;
+
+typedef	unsigned long	SINGLE;
+
+typedef	struct	{
+	unsigned long	d[2];
+}	DOUBLE;
+
+typedef	struct	{	/* expanded float format	*/
+	short	sign;
+	short	exp;
+	B64	mantissa;
+#define m1 mantissa.h_32
+#define m2 mantissa.l_32
+} EXTEND;
+
+struct	fef4_returns {
+	int	e;
+	SINGLE	f;
+};
+
+struct	fef8_returns {
+	int	e;
+	DOUBLE	f;
+};
+
+struct fif4_returns {
+	SINGLE ipart;
+	SINGLE fpart;
+};
+
+struct fif8_returns {
+	DOUBLE ipart;
+	DOUBLE fpart;
+};
+
+#if __STDC__
+#define _PROTOTYPE(function, params)	function params
+#else
+#define _PROTOTYPE(function, params)	function()
+#endif
+_PROTOTYPE( void add_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void mul_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void div_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void sub_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void sft_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void nrm_ext, (EXTEND *e1));
+_PROTOTYPE( void zrf_ext, (EXTEND *e1));
+_PROTOTYPE( void extend, (unsigned long *from, EXTEND *to, int size));
+_PROTOTYPE( void compact, (EXTEND *from, unsigned long *to, int size));
+_PROTOTYPE( void _fptrp, (int));
+_PROTOTYPE( void adf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( void adf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( void sbf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( void sbf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( void dvf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( void dvf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( void mlf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( void mlf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( void ngf4, (SINGLE f));
+_PROTOTYPE( void ngf8, (DOUBLE f));
+_PROTOTYPE( void zrf4, (SINGLE *l));
+_PROTOTYPE( void zrf8, (DOUBLE *z));
+_PROTOTYPE( void cff4, (DOUBLE src));
+_PROTOTYPE( void cff8, (SINGLE src));
+_PROTOTYPE( void cif4, (int ss, long src));
+_PROTOTYPE( void cif8, (int ss, long src));
+_PROTOTYPE( void cuf4, (int ss, long src));
+_PROTOTYPE( void cuf8, (int ss, long src));
+_PROTOTYPE( long cfu, (int ds, int ss, DOUBLE src));
+_PROTOTYPE( long cfi, (int ds, int ss, DOUBLE src));
+_PROTOTYPE( int cmf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( int cmf8, (DOUBLE d1, DOUBLE d2));
+_PROTOTYPE( void fef4, (struct fef4_returns *r, SINGLE s1));
+_PROTOTYPE( void fef8, (struct fef8_returns *r, DOUBLE s1));
+_PROTOTYPE( void fif4, (struct fif4_returns *p, SINGLE x, SINGLE y));
+_PROTOTYPE( void fif8, (struct fif8_returns *p, DOUBLE x, DOUBLE y));
+
+_PROTOTYPE( void b64_sft, (B64 *, int));
+_PROTOTYPE( void b64_lsft, (B64 *));
+_PROTOTYPE( void b64_rsft, (B64 *));
+_PROTOTYPE( int b64_add, (B64 *, B64 *));
+#endif

lib/ack/float/Makefile

@@ -0,0 +1,61 @@
+# Makefile for lib/float.
+
+CC1	= /bin/sh ./FP.compile
+
+LIBRARIES = libfp
+
+libfp_OBJECTS	= \
+	add_ext.o \
+	adder.o \
+	adf4.o \
+	adf8.o \
+	cff4.o \
+	cff8.o \
+	cfi.o \
+	cfu.o \
+	cif4.o \
+	cif8.o \
+	cmf4.o \
+	cmf8.o \
+	compact.o \
+	cuf4.o \
+	cuf8.o \
+	div_ext.o \
+	dvf4.o \
+	dvf8.o \
+	extend.o \
+	fef4.o \
+	fef8.o \
+	fif4.o \
+	fif8.o \
+	fptrp.o \
+	mlf4.o \
+	mlf8.o \
+	mul_ext.o \
+	ngf4.o \
+	ngf8.o \
+	nrm_ext.o \
+	sbf4.o \
+	sbf8.o \
+	sft_ext.o \
+	shifter.o \
+	sub_ext.o \
+	zrf4.o \
+	zrf8.o \
+	zrf_ext.o \
+
+include ../../Makefile.ack.inc
+
+#extra commands to convert the c files to the correct assembler files
+
+%.s: %.fc
+	/bin/sh ./FP.compile $<
+
+#1. make a assembler file of the c file
+#%.fs: %.fc
+#	-cp $< $(<:.fc=.c) && cc -O -I. -D_MINIX -D_POSIX_SOURCE -S $(<:.fc=.c) && cp $(<:.fc=.s) $(<:.fc=.fs)
+#	@rm $(<:.fc=.c) $(<:.fc=.s)
+
+#2. modify the assembler file
+#%.s: %.fs
+#	sed -f FP.script $< > $@

lib/ack/float/add_ext.fc

@@ -0,0 +1,56 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ADD TWO EXTENDED FORMAT NUMBERS
+*/
+
+#include "FP_types.h"
+
+void
+add_ext(e1,e2)
+register EXTEND	*e1,*e2;
+{
+	if ((e2->m1 | e2->m2) == 0L) {
+		return;
+	}
+	if ((e1->m1 | e1->m2) == 0L) {
+		*e1 = *e2;
+		return;
+	}
+	sft_ext(e1, e2);	/* adjust mantissas to equal powers */
+	if (e1->sign != e2->sign) {
+		/* e1 + e2 = e1 - (-e2) */
+		if (e2->m1 > e1->m1 ||
+                    (e2->m1 == e1->m1 && e2->m2 > e1->m2)) {
+                	/*      abs(e2) > abs(e1) */
+			EXTEND x;
+
+			x = *e1;
+			*e1 = *e2;
+                	if (x.m2 > e1->m2) {
+                        	e1->m1 -= 1;    /* carry in */
+                	}
+                	e1->m1 -= x.m1;
+                	e1->m2 -= x.m2;
+        	}
+        	else {
+                	if (e2->m2 > e1->m2)
+                        	e1->m1 -= 1;    /* carry in */
+                	e1->m1 -= e2->m1;
+                	e1->m2 -= e2->m2;
+        	}
+	}
+	else {
+		if (b64_add(&e1->mantissa,&e2->mantissa)) {	/* addition carry */
+			b64_rsft(&e1->mantissa);	/* shift mantissa one bit RIGHT */
+			e1->m1 |= 0x80000000L;	/* set max bit	*/
+			e1->exp++;		/* increase the exponent */
+		}
+	}
+	nrm_ext(e1);
+}

lib/ack/float/adder.fc

@@ -0,0 +1,50 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+ *	these are the routines the routines to do 32 and  64-bit addition
+ */
+
+# ifdef	EXT_DEBUG
+# include <stdio.h>
+# endif
+
+# include "FP_types.h"
+# define	UNKNOWN -1
+# define	TRUE	 1
+# define	FALSE	 0
+# define	MAXBIT	0x80000000L
+
+	/*
+	 *	add 64 bits
+	 */
+int
+b64_add(e1,e2)
+		/*
+		 * pointers to 64 bit 'registers'
+		 */
+register	B64	*e1,*e2;
+{
+		register	int	overflow;
+				int	carry;
+
+			/* add higher pair of 32 bits */
+	overflow = ((unsigned long) 0xFFFFFFFF - e1->h_32 < e2->h_32);
+	e1->h_32 += e2->h_32;
+
+			/* add lower pair of 32 bits */
+	carry = ((unsigned long) 0xFFFFFFFF - e1->l_32 < e2->l_32);
+	e1->l_32 += e2->l_32;
+# ifdef	EXT_DEBUG
+	printf("\t\t\t\t\tb64_add: overflow (%d); internal carry(%d)\n",
+					overflow,carry);
+	fflush(stdout);
+# endif
+	if ((carry) && (++e1->h_32 == 0))
+		return(TRUE);		/* had a 64 bit overflow */
+	return(overflow);		/* return status from higher add */
+}

lib/ack/float/adder.h

@@ -0,0 +1,15 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+ *	include file for 32 & 64 bit addition
+ */
+
+typedef	struct	B64 {
+	unsigned long	h_32;	/* higher 32 bits of 64 */
+	unsigned long	l_32;	/* lower  32 bits of 64 */
+}	B64;

lib/ack/float/adf4.fc

@@ -0,0 +1,32 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ADD TWO FLOATS - SINGLE (ADF 4)
+*/
+
+#include	"FP_types.h"
+
+void
+adf4(s2,s1)
+SINGLE	s1,s2;
+{
+	EXTEND	e1,e2;
+	int	swap = 0;
+
+	if (s1 == (SINGLE) 0) {
+		s1 = s2;
+		return;
+	}
+	if (s2 == (SINGLE) 0) {
+		return;
+	}
+	extend(&s1,&e1,sizeof(SINGLE));
+	extend(&s2,&e2,sizeof(SINGLE));
+	add_ext(&e1,&e2);
+	compact(&e1,&s1,sizeof(SINGLE));
+}

lib/ack/float/adf8.fc

@@ -0,0 +1,32 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ADD TWO FLOATS - DOUBLE (ADF 8)
+*/
+
+#include	"FP_types.h"
+
+void
+adf8(s2,s1)
+DOUBLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	if (s1.d[0] == 0 && s1.d[1] == 0) {
+		s1 = s2;
+		return;
+	}
+	if (s2.d[0] == 0 && s2.d[1] == 0) {
+		return;
+	}
+
+	extend(&s1.d[0],&e1,sizeof(DOUBLE));
+	extend(&s2.d[0],&e2,sizeof(DOUBLE));
+	add_ext(&e1,&e2);
+	compact(&e1,&s1.d[0],sizeof(DOUBLE));
+}

lib/ack/float/byte_order.h

@@ -0,0 +1,6 @@
+#define CHAR_UNSIGNED	0
+#define MSB_AT_LOW_ADDRESS	0
+#define MSW_AT_LOW_ADDRESS	0
+#define FL_MSB_AT_LOW_ADDRESS	0
+#define FL_MSW_AT_LOW_ADDRESS	0
+#define FL_MSL_AT_LOW_ADDRESS	0

lib/ack/float/cff4.fc

@@ -0,0 +1,28 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		CONVERT DOUBLE TO SINGLE (CFF 8 4)
+
+	This routine works quite simply. A floating point
+	of size 08 is converted to extended format.
+	This extended variable is converted back to
+	a floating point of size 04.
+
+*/
+
+#include	"FP_types.h"
+
+void
+cff4(src)
+DOUBLE	src;	/* the source itself -	THIS TIME it's DOUBLE */
+{
+	EXTEND	buf;
+
+	extend(&src.d[0],&buf,sizeof(DOUBLE));	/* no matter what */
+	compact(&buf,&(src.d[1]),sizeof(SINGLE));
+}

lib/ack/float/cff8.fc

@@ -0,0 +1,28 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		CONVERT SINGLE TO DOUBLE (CFF 4 8)
+
+	This routine works quite simply. A floating point
+	of size 04 is converted to extended format.
+	This extended variable is converted back to
+	a floating point of size 08.
+
+*/
+
+#include "FP_types.h"
+
+void
+cff8(src)
+SINGLE	src;
+{
+	EXTEND	buf;
+
+	extend(&src,&buf,sizeof(SINGLE));	/* no matter what */
+	compact(&buf, &src,sizeof(DOUBLE));
+}

lib/ack/float/cfi.fc

@@ -0,0 +1,52 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		CONVERT FLOAT TO SIGNED (CFI m n)
+
+		N.B. The caller must know what it is getting.
+		     A LONG is always returned. If it is an
+		     integer the high byte is cleared first.
+*/
+
+#include "FP_trap.h"
+#include "FP_types.h"
+#include "FP_shift.h"
+
+long
+cfi(ds,ss,src)
+int	ds;	/* destination size (2 or 4) */
+int	ss;	/* source size	    (4 or 8) */
+DOUBLE	src;	/* assume worst case */
+{
+	EXTEND	buf;
+	long	new;
+	short	max_exp;
+
+	extend(&src.d[0],&buf,ss);	/* get extended format */
+	if (buf.exp < 0) {	/* no conversion needed */
+		src.d[ss == 8] = 0L;
+		return(0L);
+	}
+	max_exp = (ds << 3) - 2;	/* signed numbers */
+				/* have more limited max_exp */
+	if (buf.exp > max_exp) {
+		if (buf.exp == max_exp+1 && buf.sign && buf.m1 == NORMBIT &&
+		    buf.m2 == 0L) {
+		}
+		else {
+			trap(EIOVFL);	/* integer overflow	*/
+			buf.exp %= max_exp; /* truncate	*/
+		}
+	}
+	new = buf.m1 >> (31-buf.exp);
+	if (buf.sign)
+		new = -new;
+done:
+	src.d[ss == 8] = new;
+	return(new);
+}

lib/ack/float/cfu.fc

@@ -0,0 +1,43 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		CONVERT FLOAT TO UNSIGNED (CFU m n)
+
+		N.B. The caller must know what it is getting.
+		     A LONG is always returned. If it is an
+		     integer the high byte is cleared first.
+*/
+
+#include "FP_trap.h"
+#include "FP_types.h"
+
+long
+cfu(ds,ss,src)
+int	ds;	/* destination size (2 or 4) */
+int	ss;	/* source size	    (4 or 8) */
+DOUBLE	src;	/* assume worst case */
+{
+	EXTEND	buf;
+	long	new;
+	short	newint, max_exp;
+
+	extend(&src.d[0],&buf,ss);	/* get extended format	*/
+	if (buf.exp < 0) {	/* no conversion needed	*/
+		src.d[ss == 8] = 0L;
+		return(0L);
+	}
+	max_exp = (ds << 3) - 1;
+	if (buf.exp > max_exp) {
+		trap(EIOVFL);	/* integer overflow	*/
+		buf.exp %= max_exp;
+	}
+	new = buf.m1 >> (31-buf.exp);
+done:
+	src.d[ss == 8] = new;
+	return(new);
+}

lib/ack/float/cif4.fc

@@ -0,0 +1,56 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERT INTEGER TO SINGLE (CIF n 4)
+
+	THIS ROUTINE WORKS BY FILLING AN EXTENDED
+	WITH THE INTEGER VALUE IN EXTENDED FORMAT
+	AND USES COMPACT() TO PUT IT INTO THE PROPER
+	FLOATING POINT PRECISION.
+*/
+
+#include "FP_types.h"
+
+void
+cif4(ss,src)
+int	ss;	/* source size */
+long	src;	/* largest possible integer to convert */
+{
+	EXTEND	buf;
+	short	*ipt;
+	long	i_src;
+	SINGLE	*result;
+
+	zrf_ext(&buf);
+	if (ss == sizeof(long))	{
+		buf.exp = 31;
+		i_src = src;
+		result = (SINGLE *) &src;
+	}
+	else	{
+		ipt = (short *) &src;
+		i_src = (long) *ipt;
+		buf.exp = 15;
+		result = (SINGLE *) &ss;
+	}
+	if (i_src == 0)	{
+		*result = (SINGLE) 0L;
+		return;
+	}
+			/* ESTABLISHED THAT src != 0	*/
+			/* adjust exponent field	*/
+	buf.sign = (i_src < 0) ? 0x8000 : 0;
+			/* clear sign bit of integer	*/
+			/* move to mantissa field	*/
+	buf.m1 = (i_src < 0) ? -i_src : i_src;
+			/* adjust mantissa field	*/
+	if (ss != sizeof(long))
+		buf.m1 <<= 16;
+	nrm_ext(&buf);		/* adjust mantissa field	*/
+	compact(&buf, result,sizeof(SINGLE));	/* put on stack */
+}

lib/ack/float/cif8.fc

@@ -0,0 +1,55 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERT INTEGER TO FLOAT (CIF n 8)
+
+	THIS ROUTINE WORKS BY FILLING AN EXTENDED
+	WITH THE INTEGER VALUE IN EXTENDED FORMAT
+	AND USES COMPACT() TO PUT IT INTO THE PROPER
+	FLOATING POINT PRECISION.
+*/
+
+#include "FP_types.h"
+
+void
+cif8(ss,src)
+int	ss;	/* source size */
+long	src;	/* largest possible integer to convert */
+{
+	EXTEND	buf;
+	DOUBLE	*result;	/* for return value */
+	short	*ipt;
+	long	i_src;
+
+	result = (DOUBLE *) ((void *) &ss);	/* always */
+	zrf_ext(&buf);
+	if (ss == sizeof(long))	{
+		buf.exp = 31;
+		i_src = src;
+	}
+	else	{
+		ipt = (short *) &src;
+		i_src = (long) *ipt;
+		buf.exp = 15;
+	}
+	if (i_src == 0)	{
+		zrf8(result);
+		return;
+	}
+			/* ESTABLISHED THAT src != 0	*/
+			/* adjust exponent field	*/
+	buf.sign = (i_src < 0) ? 0x8000 : 0;
+			/* clear sign bit of integer	*/
+			/* move to mantissa field	*/
+	buf.m1 = (i_src < 0) ? -i_src : i_src;
+			/* adjust mantissa field	*/
+	if (ss != sizeof(long))
+		buf.m1 <<= 16;
+	nrm_ext(&buf);
+	compact(&buf,&result->d[0],8);
+}

lib/ack/float/cmf4.fc

@@ -0,0 +1,40 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	COMPARE	SINGLES (CMF 4)
+*/
+
+#include	"FP_types.h"
+#include	"get_put.h"
+
+int
+cmf4(f1,f2)
+SINGLE	f1,f2;
+{
+		/*
+		 * return ((f1 < f2) ? 1 : (f1 - f2))
+		 */
+#define SIGN(x) (((x) < 0) ? -1 : 1)
+        int	sign1,sign2;
+	long	l1,l2;
+
+	l1 = get4((char *) &f1);
+	l2 = get4((char *) &f2);
+
+	if (l1 == l2) return 0;
+
+        sign1 = SIGN(l1);
+        sign2 = SIGN(l2);
+        if (sign1 != sign2) {
+		if ((l1 & 0x7fffffff) == 0 &&
+		    (l2 & 0x7fffffff) == 0) return 0;
+                return ((sign1 > 0) ? -1 : 1);
+	}
+
+	return (sign1 * ((l1 < l2) ? 1 : -1));
+}

lib/ack/float/cmf8.fc

@@ -0,0 +1,61 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	COMPARE	DOUBLES (CMF 8)
+*/
+
+#include	"FP_types.h"
+#include	"get_put.h"
+
+int
+cmf8(d1,d2)
+DOUBLE	d1,d2;
+{
+#define	SIGN(x)	(((x) < 0) ? -1 : 1)
+		/*
+		 * return ((d1 < d2) ? 1 : (d1 > d2) ? -1 : 0))
+		 */
+	long	l1,l2;
+	int	sign1,sign2;
+	int	rv;
+
+#if FL_MSL_AT_LOW_ADDRESS
+	l1 = get4((char *)&d1);
+	l2 = get4((char *)&d2);
+#else
+	l1 = get4(((char *)&d1+4));
+	l2 = get4(((char *)&d2+4));
+#endif
+	sign1 = SIGN(l1);
+	sign2 = SIGN(l2);
+	if (sign1 != sign2) {
+		l1 &= 0x7fffffff;
+		l2 &= 0x7fffffff;
+		if (l1 != 0 || l2 != 0) {
+			return ((sign1 > 0) ? -1 : 1);
+		}
+	}
+	if (l1 != l2)	{	/* we can decide here */
+		rv = l1 < l2 ? 1 : -1;
+	}
+	else	{ 		/* decide in 2nd half */
+		unsigned long u1, u2;
+#if FL_MSL_AT_LOW_ADDRESS
+		u1 = get4(((char *)&d1 + 4));
+		u2 = get4(((char *)&d2 + 4));
+#else
+		u1 = get4((char *)&d1);
+		u2 = get4((char *)&d2);
+#endif
+		if (u1 == u2)
+			return(0);
+		if (u1 < u2) rv = 1;
+		else rv = -1;
+	}
+	return sign1 * rv;
+}

lib/ack/float/compact.fc

@@ -0,0 +1,202 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	COMPACT EXTEND FORMAT INTO FLOAT OF PROPER SIZE
+*/
+
+# include "FP_bias.h"
+# include "FP_shift.h"
+# include "FP_trap.h"
+# include "FP_types.h"
+# include "get_put.h"
+
+void
+compact(f,to,size)
+EXTEND	*f;
+unsigned long	*to;
+int	size;
+{
+	int	error = 0;
+
+	if (size == sizeof(DOUBLE)) {
+	/*
+	 * COMPACT EXTENDED INTO DOUBLE
+	 */
+		DOUBLE *DBL = (DOUBLE *) (void *) to;
+
+		if ((f->m1|(f->m2 & DBL_ZERO)) == 0L)	{
+			zrf8(DBL);
+			return;
+		}
+		f->exp += DBL_BIAS;	/* restore proper bias	*/
+		if (f->exp > DBL_MAX)	{
+dbl_over:			trap(EFOVFL);
+			f->exp = DBL_MAX+1;
+			f->m1 = 0;
+			f->m2 = 0;
+			if (error++)
+				return;
+		}
+		else if (f->exp < DBL_MIN)	{
+			b64_rsft(&(f->mantissa));
+			if (f->exp < 0) {
+				b64_sft(&(f->mantissa), -f->exp);
+				f->exp = 0;
+			}
+			/* underflow ??? */
+		}
+			
+		/* local CAST conversion		*/
+
+		/* because of special format shift only 10 bits */
+		/* bit shift mantissa 10 bits		*/
+
+		/* first align within words, then do store operation */
+
+		DBL->d[0] = f->m1 >> DBL_RUNPACK;   /* plus 22 == 32 */
+		DBL->d[1] = f->m2 >> DBL_RUNPACK;   /* plus 22 == 32 */
+		DBL->d[1] |= (f->m1 << DBL_LUNPACK); /* plus 10 == 32 */
+
+		/* if not exact then round to nearest	*/
+		/* on a tie, round to even */
+
+#ifdef EXCEPTION_INEXACT
+		if ((f->m2 & DBL_EXACT) != 0) {
+		    INEXACT();
+#endif
+		    if (((f->m2 & DBL_EXACT) > DBL_ROUNDUP)
+			|| ((f->m2 & DBL_EXACT) == DBL_ROUNDUP
+			    && (f->m2 & (DBL_ROUNDUP << 1)))) {
+			DBL->d[1]++;	/* rounding up	*/
+			if (DBL->d[1] == 0L) { /* carry out	*/
+			    DBL->d[0]++;
+
+			    if (f->exp == 0 && (DBL->d[0] & ~DBL_MASK)) {
+					f->exp++;
+				}
+			    if (DBL->d[0] & DBL_CARRYOUT) { /* carry out */
+				if (DBL->d[0] & 01)
+				    DBL->d[1] = CARRYBIT;
+				DBL->d[0] >>= 1;
+				f->exp++;
+			    }
+			}
+			/*	check for overflow			*/
+			if (f->exp > DBL_MAX)
+		    		goto dbl_over;
+		    }
+#ifdef EXCEPTION_INEXACT
+		}
+#endif
+
+		/*
+		 * STORE EXPONENT AND SIGN:
+		 *
+		 * 1) clear leading bits (B4-B15)
+		 * 2) shift and store exponent
+		 */
+
+		DBL->d[0] &= DBL_MASK;
+		DBL->d[0] |= 
+			((long) (f->exp << DBL_EXPSHIFT) << EXP_STORE);
+		if (f->sign)
+			DBL->d[0] |= CARRYBIT;
+
+		/*
+		 * STORE MANTISSA
+		 */
+
+#if FL_MSL_AT_LOW_ADDRESS
+		put4(DBL->d[0], (char *) &DBL->d[0]);
+		put4(DBL->d[1], (char *) &DBL->d[1]);
+#else
+		{ unsigned long l;
+		  put4(DBL->d[1], (char *) &l);
+		  put4(DBL->d[0], (char *) &DBL->d[1]);
+		  DBL->d[0] = l;
+		}
+#endif
+	}
+	else {
+		/*
+		 * COMPACT EXTENDED INTO FLOAT
+		 */
+		SINGLE	*SGL;
+
+		/* local CAST conversion		*/
+		SGL = (SINGLE *) (void *) to;
+		if ((f->m1 & SGL_ZERO) == 0L)	{
+			*SGL = 0L;
+			return;
+		}
+		f->exp += SGL_BIAS;	/* restore bias	*/
+		if (f->exp > SGL_MAX)	{
+sgl_over:			trap(EFOVFL);
+			f->exp = SGL_MAX+1;
+			f->m1 = 0L;
+			f->m2 = 0L;
+			if (error++)
+				return;
+		}
+		else if (f->exp < SGL_MIN)	{
+			b64_rsft(&(f->mantissa));
+			if (f->exp < 0) {
+				b64_sft(&(f->mantissa), -f->exp);
+				f->exp = 0;
+			}
+			/* underflow ??? */
+		}
+
+		/* shift mantissa and store	*/
+		*SGL = (f->m1 >> SGL_RUNPACK);
+
+		/* check for rounding to nearest	*/
+		/* on a tie, round to even		*/
+#ifdef EXCEPTION_INEXACT
+		if (f->m2 != 0 ||
+		    (f->m1 & SGL_EXACT) != 0L) {
+			INEXACT();
+#endif
+		        if (((f->m1 & SGL_EXACT) > SGL_ROUNDUP)
+			    || ((f->m1 & SGL_EXACT) == SGL_ROUNDUP
+			        && (f->m1 & (SGL_ROUNDUP << 1)))) {
+				(*SGL)++;
+				if (f->exp == 0 && (*SGL & ~SGL_MASK)) {
+					f->exp++;
+				}
+			/* check normal */
+				if (*SGL & SGL_CARRYOUT)	{
+					*SGL >>= 1;
+					f->exp++;
+				}
+				if (f->exp > SGL_MAX)
+					goto sgl_over;
+			}
+#ifdef EXCEPTION_INEXACT
+		}
+#endif
+
+		/*
+		 * STORE EXPONENT AND SIGN:
+		 *
+		 * 1) clear leading bit of fraction
+		 * 2) shift and store exponent
+		 */
+
+		*SGL &= SGL_MASK; /* B23-B31 are 0 */
+		*SGL |= ((long) (f->exp << SGL_EXPSHIFT) << EXP_STORE);
+		if (f->sign)
+			*SGL |= CARRYBIT;
+
+		/*
+		 * STORE MANTISSA
+		 */
+
+		put4(*SGL, (char *) &SGL);
+	}
+}

lib/ack/float/cuf4.fc

@@ -0,0 +1,57 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERT INTEGER TO SINGLE (CUF n 4)
+
+	THIS ROUTINE WORKS BY FILLING AN EXTENDED
+	WITH THE INTEGER VALUE IN EXTENDED FORMAT
+	AND USES COMPACT() TO PUT IT INTO THE PROPER
+	FLOATING POINT PRECISION.
+*/
+
+#include "FP_types.h"
+
+void
+cuf4(ss,src)
+int	ss;	/* source size */
+long	src;	/* largest possible integer to convert */
+{
+	EXTEND	buf;
+	short	*ipt;
+	SINGLE	*result;
+	long	i_src;
+
+	zrf_ext(&buf);
+	if (ss == sizeof(long))	{
+		buf.exp = 31;
+		i_src = src;
+		result = (SINGLE *) &src;
+	}
+	else	{
+		ipt = (short *) &src;
+		i_src = (long) *ipt;
+		buf.exp = 15;
+		result = (SINGLE *) ((void *) &ss);
+	}
+	if (i_src == 0)	{
+		*result = (SINGLE) 0L;
+		return;
+	}
+			/* ESTABLISHED THAT src != 0	*/
+
+			/* adjust exponent field	*/
+	if (ss != sizeof(long))
+		i_src <<= 16;
+
+			/* move to mantissa field	*/
+	buf.m1 = i_src;
+
+			/* adjust mantissa field	*/
+	nrm_ext(&buf);
+	compact(&buf,result,4);
+}

lib/ack/float/cuf8.fc

@@ -0,0 +1,54 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERT INTEGER TO FLOAT (CUF n 8)
+
+	THIS ROUTINE WORKS BY FILLING AN EXTENDED
+	WITH THE INTEGER VALUE IN EXTENDED FORMAT
+	AND USES COMPACT() TO PUT IT INTO THE PROPER
+	FLOATING POINT PRECISION.
+*/
+
+#include "FP_types.h"
+
+void
+cuf8(ss,src)
+int	ss;	/* source size */
+long	src;	/* largest possible integer to convert */
+{
+	EXTEND	buf;
+	short	*ipt;
+	long	i_src;
+
+	zrf_ext(&buf);
+	if (ss == sizeof(long))	{
+		buf.exp = 31;
+		i_src = src;
+	}
+	else	{
+		ipt = (short *) &src;
+		i_src = (long) *ipt;
+		buf.exp = 15;
+	}
+	if (i_src == 0)	{
+		zrf8((DOUBLE *)((void *)&ss));
+		return;
+	}
+			/* ESTABLISHED THAT src != 0	*/
+
+			/* adjust exponent field	*/
+	if (ss != sizeof(long))
+		i_src <<= 16;
+
+			/* move to mantissa field	*/
+	buf.m1 = i_src;
+
+			/* adjust mantissa field	*/
+	nrm_ext(&buf);
+	compact(&buf,(unsigned long *) (void *)&ss,8);
+}

lib/ack/float/div_ext.fc

@@ -0,0 +1,266 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	DIVIDE EXTENDED FORMAT
+*/
+
+#include "FP_bias.h"
+#include "FP_trap.h"
+#include "FP_types.h"
+
+/*
+	November 15, 1984
+
+	This is a routine to do the work.
+	There are two versions: 
+	One is based on the partial products method
+	and makes no use possible machine instructions
+	to divide (hardware dividers).
+	The other is used when USE_DIVIDE is defined. It is much faster on
+	machines with fast 4 byte operations.
+*/
+/********************************************************/
+
+void
+div_ext(e1,e2)
+EXTEND	*e1,*e2;
+{
+	short	error = 0;
+	B64		result;
+	register	unsigned long	*lp;
+#ifndef USE_DIVIDE
+	short	count;
+#else
+	unsigned short u[9], v[5];
+	register int j;
+	register unsigned short *u_p = u;
+	int maxv = 4;
+#endif
+
+	if ((e2->m1 | e2->m2) == 0) {
+                /*
+                 * Exception 8.2 - Divide by zero
+                 */
+		trap(EFDIVZ);
+		e1->m1 = e1->m2 = 0L;
+		e1->exp = EXT_MAX;
+		return;
+	}
+	if ((e1->m1 | e1->m2) == 0) {	/* 0 / anything == 0 */
+		e1->exp = 0;	/* make sure */
+		return;
+	}
+#ifndef USE_DIVIDE
+	/*
+	 * numbers are right shifted one bit to make sure
+	 * that m1 is quaranteed to be larger if its
+	 * maximum bit is set
+	 */
+	b64_rsft(&e1->mantissa);	/* 64 bit shift right */
+	b64_rsft(&e2->mantissa);	/* 64 bit shift right */
+	e1->exp++;
+	e2->exp++;
+#endif
+	/*	check for underflow, divide by zero, etc	*/
+	e1->sign ^= e2->sign;
+	e1->exp -= e2->exp;
+
+#ifndef USE_DIVIDE
+		/* do division of mantissas	*/
+		/* uses partial product method	*/
+		/* init control variables	*/
+
+	count = 64;
+	result.h_32 = 0L;
+	result.l_32 = 0L;
+
+		/* partial product division loop */
+
+	while (count--)	{
+		/* first left shift result 1 bit	*/
+		/* this is ALWAYS done			*/
+
+		b64_lsft(&result);
+
+		/* compare dividend and divisor		*/
+		/* if dividend >= divisor add a bit	*/
+		/* and subtract divisior from dividend	*/
+
+		if ( (e1->m1 < e2->m1) ||
+			((e1->m1 == e2->m1) && (e1->m2 < e2->m2) ))
+			;	/* null statement */
+				/* i.e., don't add or subtract */
+		else	{
+			result.l_32++;	/* ADD	*/
+			if (e2->m2 > e1->m2)
+				e1->m1 -= 1;	/* carry in */
+			e1->m1 -= e2->m1;	/* do SUBTRACTION */
+			e1->m2 -= e2->m2;	/*    SUBTRACTION */
+		}
+
+		/*	shift dividend left one bit OR	*/
+		/*	IF it equals ZERO we can break out	*/
+		/*	of the loop, but still must shift	*/
+		/*	the quotient the remaining count bits	*/
+		/* NB	save the results of this test in error	*/
+		/*	if not zero, then the result is inexact. */
+		/* 	this would be reported in IEEE standard	*/
+
+		/*	lp points to dividend			*/
+		lp = &e1->m1;
+
+		error = ((*lp | *(lp+1)) != 0L) ? 1 : 0;
+		if (error)	{	/* more work */
+			/*	assume max bit == 0 (see above)	*/
+			b64_lsft(&e1->mantissa);
+			continue;
+		}
+		else
+			break;	/* leave loop	*/
+	}	/* end of divide by subtraction loop	*/
+
+	if (count > 0)	{
+		lp = &result.h_32;
+		if (count > 31) {	/* move to higher word */
+			*lp = *(lp+1);
+			count -= 32;
+			*(lp+1) = 0L;	/* clear low word	*/
+		}
+		if (*lp)
+			*lp <<= count;	/* shift rest of way	*/
+		lp++;	/*  == &result.l_32	*/
+		if (*lp) {
+			result.h_32 |= (*lp >> 32-count);
+			*lp <<= count;
+		}
+	}
+#else /* USE_DIVIDE */
+
+	u[4] = (e1->m2 & 1) << 15;
+	b64_rsft(&(e1->mantissa));
+	u[0] = e1->m1 >> 16;
+	u[1] = e1->m1;
+	u[2] = e1->m2 >> 16;
+	u[3] = e1->m2;
+	u[5] = 0; u[6] = 0; u[7] = 0;
+	v[1] = e2->m1 >> 16;
+	v[2] = e2->m1;
+	v[3] = e2->m2 >> 16;
+	v[4] = e2->m2;
+	while (! v[maxv]) maxv--;
+	result.h_32 = 0;
+	result.l_32 = 0;
+	lp = &result.h_32;
+
+	/*
+	 * Use an algorithm of Knuth (The art of programming, Seminumerical
+	 * algorithms), to divide u by v. u and v are both seen as numbers
+	 * with base 65536. 
+	 */
+	for (j = 0; j <= 3; j++, u_p++) {
+		unsigned long q_est, temp;
+
+		if (j == 2) lp++;
+		if (u_p[0] == 0 && u_p[1] < v[1]) continue;
+		temp = ((unsigned long)u_p[0] << 16) + u_p[1];
+		if (u_p[0] >= v[1]) {
+			q_est = 0x0000FFFFL;
+		}
+		else {
+			q_est = temp / v[1];
+		}
+		temp -= q_est * v[1];
+		while (temp < 0x10000 && v[2]*q_est > ((temp<<16)+u_p[2])) {
+			q_est--;
+			temp += v[1];
+		}
+		/*	Now, according to Knuth, we have an estimate of the
+			quotient, that is either correct or one too big, but
+			almost always correct.
+		*/
+		if (q_est != 0)  {
+			int i;
+			unsigned long k = 0;
+			int borrow = 0;
+
+			for (i = maxv; i > 0; i--) {
+				unsigned long tmp = q_est * v[i] + k + borrow;
+				unsigned short md = tmp;
+
+				borrow = (md > u_p[i]);
+				u_p[i] -= md;
+				k = tmp >> 16;
+			}
+			k += borrow;
+			borrow = u_p[0] < k;
+			u_p[0] -= k;
+
+			if (borrow) {
+				/* So, this does not happen often; the estimate
+				   was one too big; correct this
+				*/
+				*lp |= (j & 1) ? (q_est - 1) : ((q_est-1)<<16);
+				borrow = 0;
+				for (i = maxv; i > 0; i--) {
+					unsigned long tmp 
+					    = v[i]+(unsigned long)u_p[i]+borrow;
+					
+					u_p[i] = tmp;
+					borrow = tmp >> 16;
+				}
+				u_p[0] += borrow;
+			}
+			else *lp |= (j & 1) ? q_est : (q_est<<16);
+		}
+	}
+#ifdef	EXCEPTION_INEXACT
+	u_p = &u[0];
+	for (j = 7; j >= 0; j--) {
+		if (*u_p++) {
+			error = 1;
+			break;
+		}
+	}
+#endif
+#endif
+
+#ifdef  EXCEPTION_INEXACT
+        if (error)      {
+                /*
+                 * report here exception 8.5 - Inexact
+                 * from Draft 8.0 of IEEE P754:
+                 * In the absence of an invalid operation exception,
+                 * if the rounded result of an operation is not exact or if
+                 * it overflows without a trap, then the inexact exception
+                 * shall be assigned. The rounded or overflowed result
+                 * shall be delivered to the destination.
+                 */
+                INEXACT();
+#endif
+	e1->mantissa = result;
+
+	nrm_ext(e1);
+	if (e1->exp < EXT_MIN)	{
+		/*
+		 * Exception 8.4 - Underflow
+		 */
+		trap(EFUNFL);	/* underflow */
+		e1->exp = EXT_MIN;
+		e1->m1 = e1->m2 = 0L;
+		return;
+	}
+	if (e1->exp >= EXT_MAX) {
+                /*
+                 * Exception 8.3 - Overflow
+                 */
+                trap(EFOVFL);   /* overflow */
+                e1->exp = EXT_MAX;
+                e1->m1 = e1->m2 = 0L;
+                return;
+        }
+}

lib/ack/float/dvf4.fc

@@ -0,0 +1,26 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	DIVIDE TWO SINGLES - SINGLE Precision (dvf 4)
+*/
+
+#include	"FP_types.h"
+
+void
+dvf4(s2,s1)
+SINGLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	extend(&s1,&e1,sizeof(SINGLE));
+	extend(&s2,&e2,sizeof(SINGLE));
+
+		/* do a divide */
+	div_ext(&e1,&e2);
+	compact(&e1,&s1,sizeof(SINGLE));
+}

lib/ack/float/dvf8.fc

@@ -0,0 +1,26 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	DIVIDE TWO FLOATS - DOUBLE Precision (DVF 8)
+*/
+
+#include	"FP_types.h"
+
+void
+dvf8(s2,s1)
+DOUBLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	extend(&s1.d[0],&e1,sizeof(DOUBLE));
+	extend(&s2.d[0],&e2,sizeof(DOUBLE));
+
+		/* do a divide */
+	div_ext(&e1,&e2);
+	compact(&e1,&s1.d[0],sizeof(DOUBLE));
+}

lib/ack/float/extend.fc

@@ -0,0 +1,111 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERTS FLOATING POINT TO EXTENDED FORMAT
+
+	Two sizes of FLOATING Point are known:
+		SINGLE and DOUBLE
+*/
+/********************************************************/
+/*
+	It is not required to normalize in extended
+	format, but it has been chosen to do so.
+	Extended Format is as follows (at exit):
+
+->sign	S000 0000 | 0000 0000		<SIGN>
+->exp	0EEE EEEE | EEEE EEEE		<EXPONENT>
+->m1	LFFF FFFF | FFFF FFFF		<L.Fraction>
+	FFFF FFFF | FFFF FFFF		<Fraction>
+->m2	FFFF FFFF | FFFF FFFF		<Fraction>
+	FFFF F000 | 0000 0000		<Fraction>
+*/
+/********************************************************/
+
+#include "FP_bias.h"
+#include "FP_shift.h"
+#include "FP_types.h"
+#include "get_put.h"
+/********************************************************/
+
+void
+extend(from,to,size)
+unsigned long	*from;
+EXTEND	*to;
+int	size;
+{
+	register char *cpt1;
+	unsigned long	tmp;
+	int	leadbit = 0;
+
+	cpt1 = (char *) from;
+
+#if FL_MSL_AT_LOW_ADDRESS
+#if FL_MSW_AT_LOW_ADDRESS
+	to->exp = uget2(cpt1);
+#else
+	to->exp = uget2(cpt1+2);
+#endif
+#else
+#if FL_MSW_AT_LOW_ADDRESS
+	to->exp = uget2(cpt1+(size == sizeof(DOUBLE) ? 4 : 0));
+#else
+	to->exp = uget2(cpt1+(size == sizeof(DOUBLE) ? 6 : 2));
+#endif
+#endif
+	to->sign = (to->exp & 0x8000);	/* set sign bit */
+	to->exp ^= to->sign;
+	if (size == sizeof(DOUBLE))
+		to->exp >>= DBL_EXPSHIFT;
+	else
+		to->exp >>= SGL_EXPSHIFT;
+	if (to->exp > 0)
+		leadbit++;	/* will set Lead bit later	*/
+	else to->exp++;
+
+	if (size == sizeof(DOUBLE))	{
+#if FL_MSL_AT_LOW_ADDRESS
+		to->m1 = get4(cpt1);
+		cpt1 += 4;
+		tmp = get4(cpt1);
+#else
+		tmp = get4(cpt1);
+		cpt1 += 4;
+		to->m1 = get4(cpt1);
+#endif
+		if (to->exp == 1 && to->m1 == 0 && tmp == 0) {
+			to->exp = 0;
+			to->sign = 0;
+			to->m1 = 0;
+			to->m2 = 0;
+			return;
+		}
+		to->m1 <<= DBL_M1LEFT;		/* shift	*/
+		to->exp -= DBL_BIAS;		/* remove bias	*/
+		to->m1 |= (tmp>>DBL_RPACK);	/* plus 10 == 32	*/
+		to->m2 = (tmp<<DBL_LPACK);	/* plus 22 == 32	*/
+	}
+	else	{	/* size == sizeof(SINGLE)		*/
+		to->m1 = get4(cpt1);
+		to->m1  <<= SGL_M1LEFT;	/* shift	*/
+		if (to->exp == 1 && to->m1 == 0) {
+			to->exp = 0;
+			to->sign = 0;
+			to->m1 = 0;
+			to->m2 = 0;
+			return;
+		}
+		to->exp -= SGL_BIAS;		/* remove bias	*/
+		to->m2 = 0L;
+	}
+
+	to->m1 |= NORMBIT;				/* set bit L	*/
+	if (leadbit == 0) {		/* set or clear Leading Bit	*/
+		to->m1 &= ~NORMBIT;			/* clear bit L	*/
+		nrm_ext(to);				/* and normalize */
+	}
+}

lib/ack/float/fef4.fc

@@ -0,0 +1,33 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SEPERATE INTO EXPONENT AND FRACTION (FEF 4)
+*/
+
+#include	"FP_types.h"
+
+void
+fef4(r,s1)
+SINGLE	s1;
+struct fef4_returns	*r;
+{
+	EXTEND	buf;
+	register struct fef4_returns	*p = r;	/* make copy; r might refer
+						   to itself (see table)
+						*/
+
+	extend(&s1,&buf,sizeof(SINGLE));
+	if (buf.exp == 0 && buf.m1 == 0 && buf.m2 == 0) {
+		p->e = 0;
+	}
+	else {
+		p->e = buf.exp+1;
+		buf.exp = -1;
+	}
+	compact(&buf,&p->f,sizeof(SINGLE));
+}

lib/ack/float/fef8.fc

@@ -0,0 +1,33 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SEPERATE DOUBLE INTO EXPONENT AND FRACTION (FEF 8)
+*/
+
+#include	"FP_types.h"
+
+void
+fef8(r, s1)
+DOUBLE	s1;
+struct fef8_returns *r;
+{
+	EXTEND	buf;
+	register struct fef8_returns *p = r;	/* make copy, r might refer
+						   to itself (see table)
+						*/
+
+	extend(&s1.d[0],&buf,sizeof(DOUBLE));
+	if (buf.exp == 0 && buf.m1 == 0 && buf.m2 == 0) {
+		p->e = 0;
+	}
+	else {
+		p->e = buf.exp + 1;
+		buf.exp = -1;
+	}
+	compact(&buf,&p->f.d[0],sizeof(DOUBLE));
+}

lib/ack/float/fif4.fc

@@ -0,0 +1,46 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	MULTIPLY AND DISMEMBER PARTS (FIF 4)
+*/
+
+#include "FP_types.h"
+#include "FP_shift.h"
+
+void
+fif4(p,x,y)
+SINGLE	x,y;
+struct fif4_returns *p;
+{
+
+	EXTEND	e1,e2;
+
+	extend(&y,&e1,sizeof(SINGLE));
+	extend(&x,&e2,sizeof(SINGLE));
+		/* do a multiply */
+	mul_ext(&e1,&e2);
+	e2 = e1;
+	compact(&e2,&y,sizeof(SINGLE));
+	if (e1.exp < 0) {
+		p->ipart = 0;
+		p->fpart = y;
+		return;
+	}
+	if (e1.exp > 30 - SGL_M1LEFT) {
+		p->ipart = y;
+		p->fpart = 0;
+		return;
+	}
+	b64_sft(&e1.mantissa, 63 - e1.exp);
+	b64_sft(&e1.mantissa, e1.exp - 63);	/* "loose" low order bits */
+	compact(&e1,&(p->ipart),sizeof(SINGLE));
+	extend(&(p->ipart), &e2, sizeof(SINGLE));
+	extend(&y, &e1, sizeof(SINGLE));
+	sub_ext(&e1, &e2);
+	compact(&e1, &(p->fpart), sizeof(SINGLE));
+}

lib/ack/float/fif8.fc

@@ -0,0 +1,48 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	MULTIPLY AND DISMEMBER PARTS (FIF 8)
+*/
+
+#include "FP_types.h"
+#include "FP_shift.h"
+
+void
+fif8(p,x,y)
+DOUBLE	x,y;
+struct fif8_returns *p;
+{
+
+	EXTEND	e1,e2;
+
+	extend(&y.d[0],&e1,sizeof(DOUBLE));
+	extend(&x.d[0],&e2,sizeof(DOUBLE));
+		/* do a multiply */
+	mul_ext(&e1,&e2);
+	e2 = e1;
+	compact(&e2, &y.d[0], sizeof(DOUBLE));
+	if (e1.exp < 0) {
+		p->ipart.d[0] = 0;
+		p->ipart.d[1] = 0;
+		p->fpart = y;
+		return;
+	}
+	if (e1.exp > 62 - DBL_M1LEFT) {
+		p->ipart = y;
+		p->fpart.d[0] = 0;
+		p->fpart.d[1] = 0;
+		return;
+	}
+	b64_sft(&e1.mantissa, 63 - e1.exp);
+	b64_sft(&e1.mantissa, e1.exp - 63);	/* "loose" low order bits */
+	compact(&e1, &(p->ipart.d[0]), sizeof(DOUBLE));
+	extend(&(p->ipart.d[0]), &e2, sizeof(DOUBLE));
+	extend(&y.d[0], &e1, sizeof(DOUBLE));
+	sub_ext(&e1, &e2);
+	compact(&e1, &(p->fpart.d[0]), sizeof(DOUBLE));
+}

lib/ack/float/fptrp.s

@@ -0,0 +1,19 @@
+#
+.sect .text; .sect .rom; .sect .data; .sect .bss
+.define __fptrp
+.sect .text
+__fptrp:
+#if __i386
+	push	ebp
+	mov	ebp, esp
+	mov	eax, 8(bp)
+	call	.Xtrp
+	leave
+	ret
+#else /* i86 */
+	push	bp
+	mov	bp, sp
+	mov	ax, 4(bp)
+	call	.Xtrp
+	jmp	.cret
+#endif

lib/ack/float/get_put.h

@@ -0,0 +1,41 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+#include <byte_order.h>
+
+#if CHAR_UNSIGNED
+#define Xchar(ch)	(ch)
+#else
+#define Xchar(ch)	((ch) & 0377)
+#endif
+
+#define BYTES_REVERSED (MSB_AT_LOW_ADDRESS != FL_MSB_AT_LOW_ADDRESS)
+#define WORDS_REVERSED (MSW_AT_LOW_ADDRESS != FL_MSW_AT_LOW_ADDRESS)
+#define LONGS_REVERSED (FL_MSL_AT_LOW_ADDRESS)
+
+#if BYTES_REVERSED
+#define uget2(c)	(Xchar((c)[1]) | ((unsigned) Xchar((c)[0]) << 8))
+#define Xput2(i, c)	(((c)[1] = (i)), ((c)[0] = (i) >> 8))
+#define put2(i, c)	{ register int j = (i); Xput2(j, c); }
+#else
+#define uget2(c)	(* ((unsigned short *) (c)))
+#define Xput2(i, c)	(* ((short *) (c)) = (i))
+#define put2(i, c)	Xput2(i, c)
+#endif
+
+#define get2(c)		((short) uget2(c))
+
+#if WORDS_REVERSED || BYTES_REVERSED
+#define get4(c)		(uget2((c)+2) | ((long) uget2(c) << 16))
+#define put4(l, c)	{ register long x=(l); \
+			  Xput2((int)x,(c)+2); \
+			  Xput2((int)(x>>16),(c)); \
+			}
+#else
+#define get4(c)		(* ((long *) (c)))
+#define put4(l, c)	(* ((long *) (c)) = (l))
+#endif

lib/ack/float/mlf4.fc

@@ -0,0 +1,25 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+ * Multiply Single Precesion Float (MLF 4)
+ */
+
+#include	"FP_types.h"
+
+void
+mlf4(s2,s1)
+SINGLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	extend(&s1,&e1,sizeof(SINGLE));
+	extend(&s2,&e2,sizeof(SINGLE));
+		/* do a multiply */
+	mul_ext(&e1,&e2);
+	compact(&e1,&s1,sizeof(SINGLE));
+}

lib/ack/float/mlf8.fc

@@ -0,0 +1,25 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+ * Multiply Double Precision Float (MLF 8)
+ */
+
+#include	"FP_types.h"
+
+void
+mlf8(s2,s1)
+DOUBLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	extend(&s1.d[0],&e1,sizeof(DOUBLE));
+	extend(&s2.d[0],&e2,sizeof(DOUBLE));
+		/* do a multiply */
+	mul_ext(&e1,&e2);
+	compact(&e1,&s1.d[0],sizeof(DOUBLE));
+}

lib/ack/float/mul_ext.fc

@@ -0,0 +1,98 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ROUTINE TO MULTIPLY TWO EXTENDED FORMAT NUMBERS
+*/
+
+# include "FP_bias.h"
+# include "FP_trap.h"
+# include "FP_types.h"
+# include "FP_shift.h"
+
+void
+mul_ext(e1,e2)
+EXTEND	*e1,*e2;
+{
+	register int	i,j;		/* loop control	*/
+	unsigned short	mp[4];		/* multiplier */
+	unsigned short	mc[4];		/* multipcand */
+	unsigned short	result[8];	/* result */
+	register unsigned short *pres;
+
+	/* first save the sign (XOR)			*/
+	e1->sign ^= e2->sign;
+
+	/* compute new exponent */
+	e1->exp += e2->exp + 1;
+	/* 128 bit multiply of mantissas			*/
+
+		/* assign unknown long formats		*/
+		/* to known unsigned word formats	*/
+	mp[0] = e1->m1 >> 16;
+	mp[1] = (unsigned short) e1->m1;
+	mp[2] = e1->m2 >> 16;
+	mp[3] = (unsigned short) e1->m2;
+	mc[0] = e2->m1 >> 16;
+	mc[1] = (unsigned short) e2->m1;
+	mc[2] = e2->m2 >> 16;
+	mc[3] = (unsigned short) e2->m2;
+	for (i = 8; i--;) {
+		result[i] = 0;
+	}
+	/*
+	 *	fill registers with their components
+	 */
+	for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
+		unsigned short k = 0;
+		unsigned long mpi = mp[i];
+		for(j=4;j--;) {
+			unsigned long tmp = (unsigned long)pres[j] + k;
+			if (mc[j]) tmp += mpi * mc[j];
+			pres[j] = tmp;
+			k = tmp >> 16;
+		}
+		pres[-1] = k;
+	}
+        if (! (result[0] & 0x8000)) {
+                e1->exp--;
+                for (i = 0; i <= 3; i++) {
+                        result[i] <<= 1;
+                        if (result[i+1]&0x8000) result[i] |= 1;
+                }
+                result[4] <<= 1;
+        }
+
+	/*
+	 *	combine the registers to a total
+	 */
+	e1->m1 = ((unsigned long)(result[0]) << 16) + result[1];
+	e1->m2 = ((unsigned long)(result[2]) << 16) + result[3];
+	if (result[4] & 0x8000) {
+		if (++e1->m2 == 0)
+			if (++e1->m1 == 0) {
+				e1->m1 = NORMBIT;
+				e1->exp++;
+			}
+	}
+
+					/* check for overflow	*/
+	if (e1->exp >= EXT_MAX)	{
+		trap(EFOVFL);
+			/* if caught 			*/
+			/* return signed infinity	*/
+		e1->exp = EXT_MAX;
+infinity:	e1->m1 = e1->m2 =0L;
+		return;
+	}
+				/* check for underflow	*/
+	if (e1->exp < EXT_MIN)	{
+		trap(EFUNFL);
+		e1->exp = EXT_MIN;
+		goto infinity;
+	}
+}

lib/ack/float/ngf4.fc

@@ -0,0 +1,27 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		NEGATE A FLOATING POINT (NGF 4)
+*/
+/********************************************************/
+
+#include "FP_types.h"
+#include "get_put.h"
+
+#define OFF ((FL_MSW_AT_LOW_ADDRESS ? 0 : 2) + (FL_MSB_AT_LOW_ADDRESS ? 0 : 1))
+void
+ngf4(f)
+SINGLE	f;
+{
+	unsigned char *p;
+
+	if (f != (SINGLE) 0) {
+		p = (unsigned char *) &f + OFF;
+		*p ^= 0x80;
+	}
+}

lib/ack/float/ngf8.fc

@@ -0,0 +1,28 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		NEGATE A FLOATING POINT (NGF 8)
+*/
+/********************************************************/
+
+#include "FP_types.h"
+#include "get_put.h"
+
+#define OFF ((FL_MSL_AT_LOW_ADDRESS ? 0 : 4) + (FL_MSW_AT_LOW_ADDRESS ? 0 : 2) + (FL_MSB_AT_LOW_ADDRESS ? 0 : 1))
+
+void
+ngf8(f)
+DOUBLE	f;
+{
+	unsigned char	*p;
+
+	if (f.d[0] != 0 || f.d[1] != 0) {
+		p = (unsigned char *) &f + OFF;
+		*p ^= 0x80;
+	}
+}

lib/ack/float/nrm_ext.fc

@@ -0,0 +1,50 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/********************************************************/
+/*
+	NORMALIZE an EXTENDED FORMAT NUMBER
+*/
+/********************************************************/
+
+#include "FP_shift.h"
+#include "FP_types.h"
+
+void
+nrm_ext(e1)
+EXTEND	*e1;
+{
+		/* we assume that the mantissa != 0	*/
+		/* if it is then just return		*/
+		/* to let it be a problem elsewhere	*/
+		/* THAT IS, The exponent is not set to	*/
+		/* zero. If we don't test here an	*/
+		/* infinite loop is generated when	*/
+		/* mantissa is zero			*/
+
+	if ((e1->m1 | e1->m2) == 0L)
+		return;
+
+		/* if top word is zero mov low word	*/
+		/* to top word, adjust exponent value	*/
+	if (e1->m1 == 0L)	{
+		e1->m1 = e1->m2;
+		e1->m2 = 0L;
+		e1->exp -= 32;
+	}
+	if ((e1->m1 & NORMBIT) == 0) {
+		unsigned long l = ((unsigned long)NORMBIT >> 1);
+		int cnt = -1;
+
+		while (! (l & e1->m1)) {
+			l >>= 1;
+			cnt--;
+		}
+		e1->exp += cnt;
+		b64_sft(&(e1->mantissa), cnt);
+	}
+}

lib/ack/float/sbf4.fc

@@ -0,0 +1,27 @@
+/*
+ (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SUBTRACT TWO FLOATS - SINGLE Precision (SBF 4)
+*/
+
+#include	"FP_types.h"
+
+void
+sbf4(s2,s1)
+SINGLE	s1,s2;
+{
+	EXTEND e1,e2;
+
+	if (s2 == (SINGLE) 0) {
+		return;
+	}
+	extend(&s1,&e1,sizeof(SINGLE));
+	extend(&s2,&e2,sizeof(SINGLE));
+	sub_ext(&e1,&e2);
+	compact(&e1,&s1,sizeof(SINGLE));
+}

lib/ack/float/sbf8.fc

@@ -0,0 +1,27 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SUBTRACT TWO FLOATS - DOUBLE Precision (SBF 8)
+*/
+
+#include	"FP_types.h"
+
+void
+sbf8(s2,s1)
+DOUBLE	s1,s2;
+{
+	EXTEND e1, e2;
+
+	if (s2.d[0] == 0 && s2.d[1] == 0) {
+		return;
+	}
+	extend(&s1.d[0],&e1,sizeof(DOUBLE));
+	extend(&s2.d[0],&e2,sizeof(DOUBLE));
+	sub_ext(&e1,&e2);
+	compact(&e1,&s1.d[0],sizeof(DOUBLE));
+}

lib/ack/float/sft_ext.fc

@@ -0,0 +1,39 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SHIFT TWO EXTENDED NUMBERS INTO PROPER
+	ALIGNMENT FOR ADDITION (exponents are equal)
+	Numbers should not be zero on entry.
+*/
+
+#include "FP_types.h"
+
+void
+sft_ext(e1,e2)
+EXTEND	*e1,*e2;
+{
+	register	EXTEND	*s;
+	register	int	diff;
+
+	diff = e1->exp - e2->exp;
+
+	if (!diff)
+		return;	/* exponents are equal	*/
+
+	if (diff < 0)	{ /* e2 is larger	*/
+			/* shift e1		*/
+		diff = -diff;
+		s = e1;
+	}
+	else		/* e1 is larger		*/
+			/* shift e2		*/
+		s = e2;
+
+	s->exp += diff;
+	b64_sft(&(s->mantissa), diff);
+}

lib/ack/float/shifter.fc

@@ -0,0 +1,75 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+# include "FP_types.h"
+
+void
+b64_sft(e1,n)
+B64	*e1;
+int	n;
+{
+	if (n > 0) {
+		if (n > 63) {
+			e1->l_32 = 0;
+			e1->h_32 = 0;
+			return;
+		}
+		if (n >= 32) {
+			e1->l_32 = e1->h_32;
+			e1->h_32 = 0;
+			n -= 32;
+		}
+		if (n > 0) {
+			e1->l_32 >>= n;
+			if (e1->h_32 != 0) {
+				e1->l_32 |= (e1->h_32 << (32 - n));
+				e1->h_32 >>= n;
+			}
+		}
+		return;
+	}
+	n = -n;
+	if (n > 0) {
+		if (n > 63) {
+			e1->l_32 = 0;
+			e1->h_32 = 0;
+			return;
+		}
+		if (n >= 32) {
+			e1->h_32 = e1->l_32;
+			e1->l_32 = 0;
+			n -= 32;
+		}
+		if (n > 0) {
+			e1->h_32 <<= n;
+			if (e1->l_32 != 0) {
+				e1->h_32 |= (e1->l_32 >> (32 - n));
+				e1->l_32 <<= n;
+			}
+		}
+	}
+}
+
+void
+b64_lsft(e1)
+B64	*e1;
+{
+	/*	shift left 1 bit */
+	e1->h_32 <<= 1;
+	if (e1->l_32 & 0x80000000L) e1->h_32 |= 1;
+	e1->l_32 <<= 1;
+}
+
+void
+b64_rsft(e1)
+B64	*e1;
+{
+	/*	shift right 1 bit */
+	e1->l_32 >>= 1;
+	if (e1->h_32 & 1) e1->l_32 |= 0x80000000L;
+	e1->h_32 >>= 1;
+}

lib/ack/float/sub_ext.fc

@@ -0,0 +1,53 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SUBTRACT 2 EXTENDED FORMAT NUMBERS
+*/
+
+#include "FP_types.h"
+
+void
+sub_ext(e1,e2)
+EXTEND	*e1,*e2;
+{
+	if ((e2->m1 | e2->m2) == 0L) {
+		return;
+	}
+	if ((e1->m1 | e1->m2) == 0L) {
+		*e1 = *e2;
+		e1->sign = e2->sign ? 0 : 1;
+		return;
+	}
+	sft_ext(e1, e2);
+	if (e1->sign != e2->sign) {
+		/* e1 - e2 = e1 + (-e2) */
+		if (b64_add(&e1->mantissa,&e2->mantissa)) { /* addition carry */
+                	b64_rsft(&e1->mantissa);      /* shift mantissa one bit RIGHT */
+                	e1->m1 |= 0x80000000L;  /* set max bit  */
+                	e1->exp++;              /* increase the exponent */
+        	}
+	}
+        else if (e2->m1 > e1->m1 ||
+                 (e2->m1 == e1->m1 && e2->m2 > e1->m2)) {
+		/*	abs(e2) > abs(e1) */
+		if (e1->m2 > e2->m2) {
+			e2->m1 -= 1;	/* carry in */
+		}
+		e2->m1 -= e1->m1;
+		e2->m2 -= e1->m2;
+		*e1 = *e2;
+		e1->sign = e2->sign ? 0 : 1;
+	}
+	else {
+		if (e2->m2 > e1->m2)
+			e1->m1 -= 1;	/* carry in */
+		e1->m1 -= e2->m1;
+		e1->m2 -= e2->m2;
+	}
+	nrm_ext(e1);
+}

lib/ack/float/zrf4.fc

@@ -0,0 +1,19 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	return a zero float (ZRF 4)
+*/
+
+#include "FP_types.h"
+
+void
+zrf4(l)
+SINGLE	*l;
+{
+	*l = 0L;
+}

lib/ack/float/zrf8.fc

@@ -0,0 +1,21 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	return a zero double (ZRF 8)
+*/
+
+#include "FP_types.h"
+
+void
+zrf8(z)
+DOUBLE	*z;
+{
+
+	z->d[0] = 0L;
+	z->d[1] = 0L;
+}

lib/ack/float/zrf_ext.fc

@@ -0,0 +1,22 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ZERO and return EXTEND FORMAT FLOAT
+*/
+
+#include "FP_types.h"
+
+void
+zrf_ext(e)
+EXTEND	*e;
+{
+	e->m1 = 0;
+	e->m2 = 0;
+	e->exp = 0;
+	e->sign = 0;
+}

lib/ack/fphook/FP.compile

@@ -0,0 +1,19 @@
+#!/bin/sh
+#							Author: Kees J. Bot
+# Compile one soft FP source file.
+# (These files shouldn't be optimized normally, but the 16-bit C compiler
+# only optimizes scratch register allocation a bit with -O.  To the 32-bit
+# compiler -O is a no-op.)
+
+case $#:$1 in
+1:*.fc)	;;
+*)	echo "$0: $1: not a FC file" >&2; exit 1
+esac
+
+base="`basename "$1" .fc`"
+trap 'rm -f tmp.c tmp.s"; exit 1' 2
+
+cp "$1" tmp.c &&
+cc -O -I. -D_MINIX -D_POSIX_SOURCE -S tmp.c &&
+sed -f FP.script tmp.s > "$base.s" &&
+rm tmp.c tmp.s

lib/ack/fphook/FP.script

@@ -0,0 +1,39 @@
+s/_adf4/.adf4/
+s/_adf8/.adf8/
+s/_cff4/.cff4/
+s/_cff8/.cff8/
+s/_cfi/.cfi/
+s/_cfu/.cfu/
+s/_cif4/.cif4/
+s/_cif8/.cif8/
+s/_cmf4/.cmf4/
+s/_cmf8/.cmf8/
+s/_cuf4/.cuf4/
+s/_cuf8/.cuf8/
+s/_dvf4/.dvf4/
+s/_dvf8/.dvf8/
+s/_fef4/.fef4/
+s/_fef8/.fef8/
+s/_fif4/.fif4/
+s/_fif8/.fif8/
+s/_mlf4/.mlf4/
+s/_mlf8/.mlf8/
+s/_ngf4/.ngf4/
+s/_ngf8/.ngf8/
+s/_sbf4/.sbf4/
+s/_sbf8/.sbf8/
+s/_zrf4/.zrf4/
+s/_zrf8/.zrf8/
+s/_add_ext/.add_ext/
+s/_div_ext/.div_ext/
+s/_mul_ext/.mul_ext/
+s/_nrm_ext/.nrm_ext/
+s/_sft_ext/.sft_ext/
+s/_sub_ext/.sub_ext/
+s/_zrf_ext/.zrf_ext/
+s/_compact/.compact/
+s/_extend/.extend/
+s/_b64_add/.b64_add/
+s/_b64_sft/.b64_sft/
+s/_b64_rsft/.b64_rsft/
+s/_b64_lsft/.b64_lsft/

lib/ack/fphook/FP_bias.h

@@ -0,0 +1,28 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	include file for floating point package
+*/
+
+		/*	FLOAT FORMAT EXPONENT BIAS	*/
+
+#define	SGL_BIAS	 127	/* excess  128 notation used	*/
+#define	DBL_BIAS	1023	/* excess 1024 notation used	*/
+#define	EXT_BIAS	   0	/* 2s-complement notation used	*/
+				/* this is possible because the	*/
+				/* sign is in a seperate word	*/
+		
+		/*	VARIOUS MAX AND MIN VALUES	*/
+		/*	1) FOR THE DIFFERENT FORMATS	*/
+
+#define	SGL_MAX		   254	/*	standard definition	*/
+#define	SGL_MIN		     1	/*	standard definition	*/
+#define	DBL_MAX		  2046	/*	standard definition	*/
+#define	DBL_MIN		     1	/*	standard definition	*/
+#define EXT_MAX		 16383	/*	standard minimum	*/
+#define EXT_MIN		-16382	/*	standard minimum	*/

lib/ack/fphook/FP_shift.h

@@ -0,0 +1,49 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	include file for floating point package
+*/
+
+# define	CARRYBIT	0x80000000L
+# define	NORMBIT		0x80000000L
+# define	EXP_STORE	16
+
+
+				/* parameters for Single Precision */
+#define SGL_EXPSHIFT	7
+#define SGL_M1LEFT	8
+#define SGL_ZERO	0xffffff80L
+#define SGL_EXACT	0xff
+#define SGL_RUNPACK	SGL_M1LEFT
+
+#define SGL_ROUNDUP	0x80
+#define	SGL_CARRYOUT	0x01000000L
+#define	SGL_MASK	0x007fffffL
+
+				/* parameters for Double Precision */
+				/* used in extend.c */
+
+#define DBL_EXPSHIFT	4
+
+#define DBL_M1LEFT	11
+
+#define	DBL_RPACK	(32-DBL_M1LEFT)
+#define	DBL_LPACK	DBL_M1LEFT
+
+				/* used in compact.c */
+
+#define DBL_ZERO	0xfffffd00L
+
+#define DBL_EXACT	0x7ff
+
+#define DBL_RUNPACK	DBL_M1LEFT
+#define DBL_LUNPACK	(32-DBL_RUNPACK)
+
+#define DBL_ROUNDUP	0x400
+#define	DBL_CARRYOUT	0x00200000L
+#define	DBL_MASK	0x000fffffL

lib/ack/fphook/FP_trap.h

@@ -0,0 +1,22 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	include file for floating point package
+*/
+
+			/*	EM TRAPS	*/
+
+#define	EIOVFL	3	/* Integer  Overflow		*/
+#define	EFOVFL	4	/* Floating Overflow		*/
+#define	EFUNFL	5	/* Floating Underflow		*/
+#define	EIDIVZ	6	/* Integer  Divide by 0		*/
+#define	EFDIVZ	7	/* Floating Divide by 0.0	*/
+#define	EIUND	8	/* Integer  Undefined Number	*/
+#define	EFUND	9	/* Floating Undefined Number	*/
+#define	ECONV	10	/* Conversion Error		*/
+# define trap(x) _fptrp(x)

lib/ack/fphook/FP_types.h

@@ -0,0 +1,113 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/********************************************************/
+/*
+	Type definitions for C Floating Point Package
+	include file for floating point package
+*/
+/********************************************************/
+/*
+	THESE STRUCTURES ARE USED TO ADDRESS THE INDIVIDUAL
+	PARTS OF THE FLOATING POINT NUMBER REPRESENTATIONS.
+
+	THREE STRUCTURES ARE DEFINED:
+		SINGLE:	single precision floating format
+		DOUBLE:	double precision floating format
+		EXTEND:	double precision extended format
+*/
+/********************************************************/
+
+#ifndef __FPTYPES
+#define __FPTYPES
+
+typedef	struct	{
+	unsigned long	h_32;	/* higher 32 bits of 64 */
+	unsigned long	l_32;	/* lower  32 bits of 64 */
+}	B64;
+
+typedef	unsigned long	SINGLE;
+
+typedef	struct	{
+	unsigned long	d[2];
+}	DOUBLE;
+
+typedef	struct	{	/* expanded float format	*/
+	short	sign;
+	short	exp;
+	B64	mantissa;
+#define m1 mantissa.h_32
+#define m2 mantissa.l_32
+} EXTEND;
+
+struct	fef4_returns {
+	int	e;
+	SINGLE	f;
+};
+
+struct	fef8_returns {
+	int	e;
+	DOUBLE	f;
+};
+
+struct fif4_returns {
+	SINGLE ipart;
+	SINGLE fpart;
+};
+
+struct fif8_returns {
+	DOUBLE ipart;
+	DOUBLE fpart;
+};
+
+#if __STDC__
+#define _PROTOTYPE(function, params)	function params
+#else
+#define _PROTOTYPE(function, params)	function()
+#endif
+_PROTOTYPE( void add_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void mul_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void div_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void sub_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void sft_ext, (EXTEND *e1, EXTEND *e2));
+_PROTOTYPE( void nrm_ext, (EXTEND *e1));
+_PROTOTYPE( void zrf_ext, (EXTEND *e1));
+_PROTOTYPE( void extend, (unsigned long *from, EXTEND *to, int size));
+_PROTOTYPE( void compact, (EXTEND *from, unsigned long *to, int size));
+_PROTOTYPE( void _fptrp, (int));
+_PROTOTYPE( void adf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( void adf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( void sbf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( void sbf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( void dvf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( void dvf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( void mlf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( void mlf8, (DOUBLE s2, DOUBLE s1));
+_PROTOTYPE( void ngf4, (SINGLE f));
+_PROTOTYPE( void ngf8, (DOUBLE f));
+_PROTOTYPE( void zrf4, (SINGLE *l));
+_PROTOTYPE( void zrf8, (DOUBLE *z));
+_PROTOTYPE( void cff4, (DOUBLE src));
+_PROTOTYPE( void cff8, (SINGLE src));
+_PROTOTYPE( void cif4, (int ss, long src));
+_PROTOTYPE( void cif8, (int ss, long src));
+_PROTOTYPE( void cuf4, (int ss, long src));
+_PROTOTYPE( void cuf8, (int ss, long src));
+_PROTOTYPE( long cfu, (int ds, int ss, DOUBLE src));
+_PROTOTYPE( long cfi, (int ds, int ss, DOUBLE src));
+_PROTOTYPE( int cmf4, (SINGLE s2, SINGLE s1));
+_PROTOTYPE( int cmf8, (DOUBLE d1, DOUBLE d2));
+_PROTOTYPE( void fef4, (struct fef4_returns *r, SINGLE s1));
+_PROTOTYPE( void fef8, (struct fef8_returns *r, DOUBLE s1));
+_PROTOTYPE( void fif4, (struct fif4_returns *p, SINGLE x, SINGLE y));
+_PROTOTYPE( void fif8, (struct fif8_returns *p, DOUBLE x, DOUBLE y));
+
+_PROTOTYPE( void b64_sft, (B64 *, int));
+_PROTOTYPE( void b64_lsft, (B64 *));
+_PROTOTYPE( void b64_rsft, (B64 *));
+_PROTOTYPE( int b64_add, (B64 *, B64 *));
+#endif

lib/ack/fphook/Makefile

@@ -0,0 +1,19 @@
+# Makefile for lib/ack/fphook.
+
+# The ACK ANSI C compiler has an nice trick to reduce the size of programs
+# that do not use floating point.  If a program uses floating point then the
+# compiler generates an external reference to the label '_fp_hook'.  This makes
+# the loader bring in the floating point printing and conversion routines
+# '_f_print' and 'strtod' from the library libd.a.  Otherwise two dummy
+# routines are found in libc.a.  (The printf and scanf need floating point
+# for the %f formats, whether you use them or not.)
+
+CFLAGS	= -O -D_MINIX -D_POSIX_SOURCE -I..
+
+LIBRARIES = libd libc
+
+libd_OBJECTS = fphook.o
+
+libc_OBJECTS = fltpr.o strtod.o
+
+include ../../Makefile.ack.inc

lib/ack/fphook/Makefile.ack.conv

@@ -0,0 +1,32 @@
+# Makefile for lib/i86/fphook.
+
+# The ACK ANSI C compiler has an nice trick to reduce the size of programs
+# that do not use floating point.  If a program uses floating point then the
+# compiler generates an external reference to the label '_fp_hook'.  This makes
+# the loader bring in the floating point printing and conversion routines
+# '_f_print' and 'strtod' from the library libd.a.  Otherwise two dummy
+# routines are found in libc.a.  (The printf and scanf need floating point
+# for the %f formats, whether you use them or not.)
+
+CFLAGS	= -O -D_MINIX -D_POSIX_SOURCE
+CC1	= $(CC) $(CFLAGS) -c -I$(SRCDIR)/ack
+
+LIBD	= ../../libd.a
+LIBC	= ../../libc.a
+
+all:	$(LIBD) $(LIBC)
+
+$(LIBD):	fphook.c
+	$(CC1) fphook.c
+	aal cr $@ fphook.o
+	rm fphook.o
+
+$(LIBC):	$(LIBC)(fltpr.o) $(LIBC)(strtod.o)
+	aal cr $@ *.o
+	rm *.o
+
+$(LIBC)(fltpr.o): fltpr.c
+	$(CC1) fltpr.c
+
+$(LIBC)(strtod.o): strtod.c
+	$(CC1) strtod.c

lib/ack/fphook/add_ext.fc

@@ -0,0 +1,56 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ADD TWO EXTENDED FORMAT NUMBERS
+*/
+
+#include "FP_types.h"
+
+void
+add_ext(e1,e2)
+register EXTEND	*e1,*e2;
+{
+	if ((e2->m1 | e2->m2) == 0L) {
+		return;
+	}
+	if ((e1->m1 | e1->m2) == 0L) {
+		*e1 = *e2;
+		return;
+	}
+	sft_ext(e1, e2);	/* adjust mantissas to equal powers */
+	if (e1->sign != e2->sign) {
+		/* e1 + e2 = e1 - (-e2) */
+		if (e2->m1 > e1->m1 ||
+                    (e2->m1 == e1->m1 && e2->m2 > e1->m2)) {
+                	/*      abs(e2) > abs(e1) */
+			EXTEND x;
+
+			x = *e1;
+			*e1 = *e2;
+                	if (x.m2 > e1->m2) {
+                        	e1->m1 -= 1;    /* carry in */
+                	}
+                	e1->m1 -= x.m1;
+                	e1->m2 -= x.m2;
+        	}
+        	else {
+                	if (e2->m2 > e1->m2)
+                        	e1->m1 -= 1;    /* carry in */
+                	e1->m1 -= e2->m1;
+                	e1->m2 -= e2->m2;
+        	}
+	}
+	else {
+		if (b64_add(&e1->mantissa,&e2->mantissa)) {	/* addition carry */
+			b64_rsft(&e1->mantissa);	/* shift mantissa one bit RIGHT */
+			e1->m1 |= 0x80000000L;	/* set max bit	*/
+			e1->exp++;		/* increase the exponent */
+		}
+	}
+	nrm_ext(e1);
+}

lib/ack/fphook/adder.fc

@@ -0,0 +1,50 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+ *	these are the routines the routines to do 32 and  64-bit addition
+ */
+
+# ifdef	EXT_DEBUG
+# include <stdio.h>
+# endif
+
+# include "FP_types.h"
+# define	UNKNOWN -1
+# define	TRUE	 1
+# define	FALSE	 0
+# define	MAXBIT	0x80000000L
+
+	/*
+	 *	add 64 bits
+	 */
+int
+b64_add(e1,e2)
+		/*
+		 * pointers to 64 bit 'registers'
+		 */
+register	B64	*e1,*e2;
+{
+		register	int	overflow;
+				int	carry;
+
+			/* add higher pair of 32 bits */
+	overflow = ((unsigned long) 0xFFFFFFFF - e1->h_32 < e2->h_32);
+	e1->h_32 += e2->h_32;
+
+			/* add lower pair of 32 bits */
+	carry = ((unsigned long) 0xFFFFFFFF - e1->l_32 < e2->l_32);
+	e1->l_32 += e2->l_32;
+# ifdef	EXT_DEBUG
+	printf("\t\t\t\t\tb64_add: overflow (%d); internal carry(%d)\n",
+					overflow,carry);
+	fflush(stdout);
+# endif
+	if ((carry) && (++e1->h_32 == 0))
+		return(TRUE);		/* had a 64 bit overflow */
+	return(overflow);		/* return status from higher add */
+}

lib/ack/fphook/adder.h

@@ -0,0 +1,15 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+ *	include file for 32 & 64 bit addition
+ */
+
+typedef	struct	B64 {
+	unsigned long	h_32;	/* higher 32 bits of 64 */
+	unsigned long	l_32;	/* lower  32 bits of 64 */
+}	B64;

lib/ack/fphook/adf4.fc

@@ -0,0 +1,32 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ADD TWO FLOATS - SINGLE (ADF 4)
+*/
+
+#include	"FP_types.h"
+
+void
+adf4(s2,s1)
+SINGLE	s1,s2;
+{
+	EXTEND	e1,e2;
+	int	swap = 0;
+
+	if (s1 == (SINGLE) 0) {
+		s1 = s2;
+		return;
+	}
+	if (s2 == (SINGLE) 0) {
+		return;
+	}
+	extend(&s1,&e1,sizeof(SINGLE));
+	extend(&s2,&e2,sizeof(SINGLE));
+	add_ext(&e1,&e2);
+	compact(&e1,&s1,sizeof(SINGLE));
+}

lib/ack/fphook/adf8.fc

@@ -0,0 +1,32 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ADD TWO FLOATS - DOUBLE (ADF 8)
+*/
+
+#include	"FP_types.h"
+
+void
+adf8(s2,s1)
+DOUBLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	if (s1.d[0] == 0 && s1.d[1] == 0) {
+		s1 = s2;
+		return;
+	}
+	if (s2.d[0] == 0 && s2.d[1] == 0) {
+		return;
+	}
+
+	extend(&s1.d[0],&e1,sizeof(DOUBLE));
+	extend(&s2.d[0],&e2,sizeof(DOUBLE));
+	add_ext(&e1,&e2);
+	compact(&e1,&s1.d[0],sizeof(DOUBLE));
+}

lib/ack/fphook/byte_order.h

@@ -0,0 +1,6 @@
+#define CHAR_UNSIGNED	0
+#define MSB_AT_LOW_ADDRESS	0
+#define MSW_AT_LOW_ADDRESS	0
+#define FL_MSB_AT_LOW_ADDRESS	0
+#define FL_MSW_AT_LOW_ADDRESS	0
+#define FL_MSL_AT_LOW_ADDRESS	0

lib/ack/fphook/cff4.fc

@@ -0,0 +1,28 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		CONVERT DOUBLE TO SINGLE (CFF 8 4)
+
+	This routine works quite simply. A floating point
+	of size 08 is converted to extended format.
+	This extended variable is converted back to
+	a floating point of size 04.
+
+*/
+
+#include	"FP_types.h"
+
+void
+cff4(src)
+DOUBLE	src;	/* the source itself -	THIS TIME it's DOUBLE */
+{
+	EXTEND	buf;
+
+	extend(&src.d[0],&buf,sizeof(DOUBLE));	/* no matter what */
+	compact(&buf,&(src.d[1]),sizeof(SINGLE));
+}

lib/ack/fphook/cff8.fc

@@ -0,0 +1,28 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		CONVERT SINGLE TO DOUBLE (CFF 4 8)
+
+	This routine works quite simply. A floating point
+	of size 04 is converted to extended format.
+	This extended variable is converted back to
+	a floating point of size 08.
+
+*/
+
+#include "FP_types.h"
+
+void
+cff8(src)
+SINGLE	src;
+{
+	EXTEND	buf;
+
+	extend(&src,&buf,sizeof(SINGLE));	/* no matter what */
+	compact(&buf, &src,sizeof(DOUBLE));
+}

lib/ack/fphook/cfi.fc

@@ -0,0 +1,52 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		CONVERT FLOAT TO SIGNED (CFI m n)
+
+		N.B. The caller must know what it is getting.
+		     A LONG is always returned. If it is an
+		     integer the high byte is cleared first.
+*/
+
+#include "FP_trap.h"
+#include "FP_types.h"
+#include "FP_shift.h"
+
+long
+cfi(ds,ss,src)
+int	ds;	/* destination size (2 or 4) */
+int	ss;	/* source size	    (4 or 8) */
+DOUBLE	src;	/* assume worst case */
+{
+	EXTEND	buf;
+	long	new;
+	short	max_exp;
+
+	extend(&src.d[0],&buf,ss);	/* get extended format */
+	if (buf.exp < 0) {	/* no conversion needed */
+		src.d[ss == 8] = 0L;
+		return(0L);
+	}
+	max_exp = (ds << 3) - 2;	/* signed numbers */
+				/* have more limited max_exp */
+	if (buf.exp > max_exp) {
+		if (buf.exp == max_exp+1 && buf.sign && buf.m1 == NORMBIT &&
+		    buf.m2 == 0L) {
+		}
+		else {
+			trap(EIOVFL);	/* integer overflow	*/
+			buf.exp %= max_exp; /* truncate	*/
+		}
+	}
+	new = buf.m1 >> (31-buf.exp);
+	if (buf.sign)
+		new = -new;
+done:
+	src.d[ss == 8] = new;
+	return(new);
+}

lib/ack/fphook/cfu.fc

@@ -0,0 +1,43 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		CONVERT FLOAT TO UNSIGNED (CFU m n)
+
+		N.B. The caller must know what it is getting.
+		     A LONG is always returned. If it is an
+		     integer the high byte is cleared first.
+*/
+
+#include "FP_trap.h"
+#include "FP_types.h"
+
+long
+cfu(ds,ss,src)
+int	ds;	/* destination size (2 or 4) */
+int	ss;	/* source size	    (4 or 8) */
+DOUBLE	src;	/* assume worst case */
+{
+	EXTEND	buf;
+	long	new;
+	short	newint, max_exp;
+
+	extend(&src.d[0],&buf,ss);	/* get extended format	*/
+	if (buf.exp < 0) {	/* no conversion needed	*/
+		src.d[ss == 8] = 0L;
+		return(0L);
+	}
+	max_exp = (ds << 3) - 1;
+	if (buf.exp > max_exp) {
+		trap(EIOVFL);	/* integer overflow	*/
+		buf.exp %= max_exp;
+	}
+	new = buf.m1 >> (31-buf.exp);
+done:
+	src.d[ss == 8] = new;
+	return(new);
+}

lib/ack/fphook/cif4.fc

@@ -0,0 +1,56 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERT INTEGER TO SINGLE (CIF n 4)
+
+	THIS ROUTINE WORKS BY FILLING AN EXTENDED
+	WITH THE INTEGER VALUE IN EXTENDED FORMAT
+	AND USES COMPACT() TO PUT IT INTO THE PROPER
+	FLOATING POINT PRECISION.
+*/
+
+#include "FP_types.h"
+
+void
+cif4(ss,src)
+int	ss;	/* source size */
+long	src;	/* largest possible integer to convert */
+{
+	EXTEND	buf;
+	short	*ipt;
+	long	i_src;
+	SINGLE	*result;
+
+	zrf_ext(&buf);
+	if (ss == sizeof(long))	{
+		buf.exp = 31;
+		i_src = src;
+		result = (SINGLE *) &src;
+	}
+	else	{
+		ipt = (short *) &src;
+		i_src = (long) *ipt;
+		buf.exp = 15;
+		result = (SINGLE *) &ss;
+	}
+	if (i_src == 0)	{
+		*result = (SINGLE) 0L;
+		return;
+	}
+			/* ESTABLISHED THAT src != 0	*/
+			/* adjust exponent field	*/
+	buf.sign = (i_src < 0) ? 0x8000 : 0;
+			/* clear sign bit of integer	*/
+			/* move to mantissa field	*/
+	buf.m1 = (i_src < 0) ? -i_src : i_src;
+			/* adjust mantissa field	*/
+	if (ss != sizeof(long))
+		buf.m1 <<= 16;
+	nrm_ext(&buf);		/* adjust mantissa field	*/
+	compact(&buf, result,sizeof(SINGLE));	/* put on stack */
+}

lib/ack/fphook/cif8.fc

@@ -0,0 +1,55 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERT INTEGER TO FLOAT (CIF n 8)
+
+	THIS ROUTINE WORKS BY FILLING AN EXTENDED
+	WITH THE INTEGER VALUE IN EXTENDED FORMAT
+	AND USES COMPACT() TO PUT IT INTO THE PROPER
+	FLOATING POINT PRECISION.
+*/
+
+#include "FP_types.h"
+
+void
+cif8(ss,src)
+int	ss;	/* source size */
+long	src;	/* largest possible integer to convert */
+{
+	EXTEND	buf;
+	DOUBLE	*result;	/* for return value */
+	short	*ipt;
+	long	i_src;
+
+	result = (DOUBLE *) ((void *) &ss);	/* always */
+	zrf_ext(&buf);
+	if (ss == sizeof(long))	{
+		buf.exp = 31;
+		i_src = src;
+	}
+	else	{
+		ipt = (short *) &src;
+		i_src = (long) *ipt;
+		buf.exp = 15;
+	}
+	if (i_src == 0)	{
+		zrf8(result);
+		return;
+	}
+			/* ESTABLISHED THAT src != 0	*/
+			/* adjust exponent field	*/
+	buf.sign = (i_src < 0) ? 0x8000 : 0;
+			/* clear sign bit of integer	*/
+			/* move to mantissa field	*/
+	buf.m1 = (i_src < 0) ? -i_src : i_src;
+			/* adjust mantissa field	*/
+	if (ss != sizeof(long))
+		buf.m1 <<= 16;
+	nrm_ext(&buf);
+	compact(&buf,&result->d[0],8);
+}

lib/ack/fphook/cmf4.fc

@@ -0,0 +1,40 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	COMPARE	SINGLES (CMF 4)
+*/
+
+#include	"FP_types.h"
+#include	"get_put.h"
+
+int
+cmf4(f1,f2)
+SINGLE	f1,f2;
+{
+		/*
+		 * return ((f1 < f2) ? 1 : (f1 - f2))
+		 */
+#define SIGN(x) (((x) < 0) ? -1 : 1)
+        int	sign1,sign2;
+	long	l1,l2;
+
+	l1 = get4((char *) &f1);
+	l2 = get4((char *) &f2);
+
+	if (l1 == l2) return 0;
+
+        sign1 = SIGN(l1);
+        sign2 = SIGN(l2);
+        if (sign1 != sign2) {
+		if ((l1 & 0x7fffffff) == 0 &&
+		    (l2 & 0x7fffffff) == 0) return 0;
+                return ((sign1 > 0) ? -1 : 1);
+	}
+
+	return (sign1 * ((l1 < l2) ? 1 : -1));
+}

lib/ack/fphook/cmf8.fc

@@ -0,0 +1,61 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	COMPARE	DOUBLES (CMF 8)
+*/
+
+#include	"FP_types.h"
+#include	"get_put.h"
+
+int
+cmf8(d1,d2)
+DOUBLE	d1,d2;
+{
+#define	SIGN(x)	(((x) < 0) ? -1 : 1)
+		/*
+		 * return ((d1 < d2) ? 1 : (d1 > d2) ? -1 : 0))
+		 */
+	long	l1,l2;
+	int	sign1,sign2;
+	int	rv;
+
+#if FL_MSL_AT_LOW_ADDRESS
+	l1 = get4((char *)&d1);
+	l2 = get4((char *)&d2);
+#else
+	l1 = get4(((char *)&d1+4));
+	l2 = get4(((char *)&d2+4));
+#endif
+	sign1 = SIGN(l1);
+	sign2 = SIGN(l2);
+	if (sign1 != sign2) {
+		l1 &= 0x7fffffff;
+		l2 &= 0x7fffffff;
+		if (l1 != 0 || l2 != 0) {
+			return ((sign1 > 0) ? -1 : 1);
+		}
+	}
+	if (l1 != l2)	{	/* we can decide here */
+		rv = l1 < l2 ? 1 : -1;
+	}
+	else	{ 		/* decide in 2nd half */
+		unsigned long u1, u2;
+#if FL_MSL_AT_LOW_ADDRESS
+		u1 = get4(((char *)&d1 + 4));
+		u2 = get4(((char *)&d2 + 4));
+#else
+		u1 = get4((char *)&d1);
+		u2 = get4((char *)&d2);
+#endif
+		if (u1 == u2)
+			return(0);
+		if (u1 < u2) rv = 1;
+		else rv = -1;
+	}
+	return sign1 * rv;
+}

lib/ack/fphook/compact.fc

@@ -0,0 +1,202 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	COMPACT EXTEND FORMAT INTO FLOAT OF PROPER SIZE
+*/
+
+# include "FP_bias.h"
+# include "FP_shift.h"
+# include "FP_trap.h"
+# include "FP_types.h"
+# include "get_put.h"
+
+void
+compact(f,to,size)
+EXTEND	*f;
+unsigned long	*to;
+int	size;
+{
+	int	error = 0;
+
+	if (size == sizeof(DOUBLE)) {
+	/*
+	 * COMPACT EXTENDED INTO DOUBLE
+	 */
+		DOUBLE *DBL = (DOUBLE *) (void *) to;
+
+		if ((f->m1|(f->m2 & DBL_ZERO)) == 0L)	{
+			zrf8(DBL);
+			return;
+		}
+		f->exp += DBL_BIAS;	/* restore proper bias	*/
+		if (f->exp > DBL_MAX)	{
+dbl_over:			trap(EFOVFL);
+			f->exp = DBL_MAX+1;
+			f->m1 = 0;
+			f->m2 = 0;
+			if (error++)
+				return;
+		}
+		else if (f->exp < DBL_MIN)	{
+			b64_rsft(&(f->mantissa));
+			if (f->exp < 0) {
+				b64_sft(&(f->mantissa), -f->exp);
+				f->exp = 0;
+			}
+			/* underflow ??? */
+		}
+			
+		/* local CAST conversion		*/
+
+		/* because of special format shift only 10 bits */
+		/* bit shift mantissa 10 bits		*/
+
+		/* first align within words, then do store operation */
+
+		DBL->d[0] = f->m1 >> DBL_RUNPACK;   /* plus 22 == 32 */
+		DBL->d[1] = f->m2 >> DBL_RUNPACK;   /* plus 22 == 32 */
+		DBL->d[1] |= (f->m1 << DBL_LUNPACK); /* plus 10 == 32 */
+
+		/* if not exact then round to nearest	*/
+		/* on a tie, round to even */
+
+#ifdef EXCEPTION_INEXACT
+		if ((f->m2 & DBL_EXACT) != 0) {
+		    INEXACT();
+#endif
+		    if (((f->m2 & DBL_EXACT) > DBL_ROUNDUP)
+			|| ((f->m2 & DBL_EXACT) == DBL_ROUNDUP
+			    && (f->m2 & (DBL_ROUNDUP << 1)))) {
+			DBL->d[1]++;	/* rounding up	*/
+			if (DBL->d[1] == 0L) { /* carry out	*/
+			    DBL->d[0]++;
+
+			    if (f->exp == 0 && (DBL->d[0] & ~DBL_MASK)) {
+					f->exp++;
+				}
+			    if (DBL->d[0] & DBL_CARRYOUT) { /* carry out */
+				if (DBL->d[0] & 01)
+				    DBL->d[1] = CARRYBIT;
+				DBL->d[0] >>= 1;
+				f->exp++;
+			    }
+			}
+			/*	check for overflow			*/
+			if (f->exp > DBL_MAX)
+		    		goto dbl_over;
+		    }
+#ifdef EXCEPTION_INEXACT
+		}
+#endif
+
+		/*
+		 * STORE EXPONENT AND SIGN:
+		 *
+		 * 1) clear leading bits (B4-B15)
+		 * 2) shift and store exponent
+		 */
+
+		DBL->d[0] &= DBL_MASK;
+		DBL->d[0] |= 
+			((long) (f->exp << DBL_EXPSHIFT) << EXP_STORE);
+		if (f->sign)
+			DBL->d[0] |= CARRYBIT;
+
+		/*
+		 * STORE MANTISSA
+		 */
+
+#if FL_MSL_AT_LOW_ADDRESS
+		put4(DBL->d[0], (char *) &DBL->d[0]);
+		put4(DBL->d[1], (char *) &DBL->d[1]);
+#else
+		{ unsigned long l;
+		  put4(DBL->d[1], (char *) &l);
+		  put4(DBL->d[0], (char *) &DBL->d[1]);
+		  DBL->d[0] = l;
+		}
+#endif
+	}
+	else {
+		/*
+		 * COMPACT EXTENDED INTO FLOAT
+		 */
+		SINGLE	*SGL;
+
+		/* local CAST conversion		*/
+		SGL = (SINGLE *) (void *) to;
+		if ((f->m1 & SGL_ZERO) == 0L)	{
+			*SGL = 0L;
+			return;
+		}
+		f->exp += SGL_BIAS;	/* restore bias	*/
+		if (f->exp > SGL_MAX)	{
+sgl_over:			trap(EFOVFL);
+			f->exp = SGL_MAX+1;
+			f->m1 = 0L;
+			f->m2 = 0L;
+			if (error++)
+				return;
+		}
+		else if (f->exp < SGL_MIN)	{
+			b64_rsft(&(f->mantissa));
+			if (f->exp < 0) {
+				b64_sft(&(f->mantissa), -f->exp);
+				f->exp = 0;
+			}
+			/* underflow ??? */
+		}
+
+		/* shift mantissa and store	*/
+		*SGL = (f->m1 >> SGL_RUNPACK);
+
+		/* check for rounding to nearest	*/
+		/* on a tie, round to even		*/
+#ifdef EXCEPTION_INEXACT
+		if (f->m2 != 0 ||
+		    (f->m1 & SGL_EXACT) != 0L) {
+			INEXACT();
+#endif
+		        if (((f->m1 & SGL_EXACT) > SGL_ROUNDUP)
+			    || ((f->m1 & SGL_EXACT) == SGL_ROUNDUP
+			        && (f->m1 & (SGL_ROUNDUP << 1)))) {
+				(*SGL)++;
+				if (f->exp == 0 && (*SGL & ~SGL_MASK)) {
+					f->exp++;
+				}
+			/* check normal */
+				if (*SGL & SGL_CARRYOUT)	{
+					*SGL >>= 1;
+					f->exp++;
+				}
+				if (f->exp > SGL_MAX)
+					goto sgl_over;
+			}
+#ifdef EXCEPTION_INEXACT
+		}
+#endif
+
+		/*
+		 * STORE EXPONENT AND SIGN:
+		 *
+		 * 1) clear leading bit of fraction
+		 * 2) shift and store exponent
+		 */
+
+		*SGL &= SGL_MASK; /* B23-B31 are 0 */
+		*SGL |= ((long) (f->exp << SGL_EXPSHIFT) << EXP_STORE);
+		if (f->sign)
+			*SGL |= CARRYBIT;
+
+		/*
+		 * STORE MANTISSA
+		 */
+
+		put4(*SGL, (char *) &SGL);
+	}
+}

lib/ack/fphook/cuf4.fc

@@ -0,0 +1,57 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERT INTEGER TO SINGLE (CUF n 4)
+
+	THIS ROUTINE WORKS BY FILLING AN EXTENDED
+	WITH THE INTEGER VALUE IN EXTENDED FORMAT
+	AND USES COMPACT() TO PUT IT INTO THE PROPER
+	FLOATING POINT PRECISION.
+*/
+
+#include "FP_types.h"
+
+void
+cuf4(ss,src)
+int	ss;	/* source size */
+long	src;	/* largest possible integer to convert */
+{
+	EXTEND	buf;
+	short	*ipt;
+	SINGLE	*result;
+	long	i_src;
+
+	zrf_ext(&buf);
+	if (ss == sizeof(long))	{
+		buf.exp = 31;
+		i_src = src;
+		result = (SINGLE *) &src;
+	}
+	else	{
+		ipt = (short *) &src;
+		i_src = (long) *ipt;
+		buf.exp = 15;
+		result = (SINGLE *) ((void *) &ss);
+	}
+	if (i_src == 0)	{
+		*result = (SINGLE) 0L;
+		return;
+	}
+			/* ESTABLISHED THAT src != 0	*/
+
+			/* adjust exponent field	*/
+	if (ss != sizeof(long))
+		i_src <<= 16;
+
+			/* move to mantissa field	*/
+	buf.m1 = i_src;
+
+			/* adjust mantissa field	*/
+	nrm_ext(&buf);
+	compact(&buf,result,4);
+}

lib/ack/fphook/cuf8.fc

@@ -0,0 +1,54 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERT INTEGER TO FLOAT (CUF n 8)
+
+	THIS ROUTINE WORKS BY FILLING AN EXTENDED
+	WITH THE INTEGER VALUE IN EXTENDED FORMAT
+	AND USES COMPACT() TO PUT IT INTO THE PROPER
+	FLOATING POINT PRECISION.
+*/
+
+#include "FP_types.h"
+
+void
+cuf8(ss,src)
+int	ss;	/* source size */
+long	src;	/* largest possible integer to convert */
+{
+	EXTEND	buf;
+	short	*ipt;
+	long	i_src;
+
+	zrf_ext(&buf);
+	if (ss == sizeof(long))	{
+		buf.exp = 31;
+		i_src = src;
+	}
+	else	{
+		ipt = (short *) &src;
+		i_src = (long) *ipt;
+		buf.exp = 15;
+	}
+	if (i_src == 0)	{
+		zrf8((DOUBLE *)((void *)&ss));
+		return;
+	}
+			/* ESTABLISHED THAT src != 0	*/
+
+			/* adjust exponent field	*/
+	if (ss != sizeof(long))
+		i_src <<= 16;
+
+			/* move to mantissa field	*/
+	buf.m1 = i_src;
+
+			/* adjust mantissa field	*/
+	nrm_ext(&buf);
+	compact(&buf,(unsigned long *) (void *)&ss,8);
+}

lib/ack/fphook/div_ext.fc

@@ -0,0 +1,266 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	DIVIDE EXTENDED FORMAT
+*/
+
+#include "FP_bias.h"
+#include "FP_trap.h"
+#include "FP_types.h"
+
+/*
+	November 15, 1984
+
+	This is a routine to do the work.
+	There are two versions: 
+	One is based on the partial products method
+	and makes no use possible machine instructions
+	to divide (hardware dividers).
+	The other is used when USE_DIVIDE is defined. It is much faster on
+	machines with fast 4 byte operations.
+*/
+/********************************************************/
+
+void
+div_ext(e1,e2)
+EXTEND	*e1,*e2;
+{
+	short	error = 0;
+	B64		result;
+	register	unsigned long	*lp;
+#ifndef USE_DIVIDE
+	short	count;
+#else
+	unsigned short u[9], v[5];
+	register int j;
+	register unsigned short *u_p = u;
+	int maxv = 4;
+#endif
+
+	if ((e2->m1 | e2->m2) == 0) {
+                /*
+                 * Exception 8.2 - Divide by zero
+                 */
+		trap(EFDIVZ);
+		e1->m1 = e1->m2 = 0L;
+		e1->exp = EXT_MAX;
+		return;
+	}
+	if ((e1->m1 | e1->m2) == 0) {	/* 0 / anything == 0 */
+		e1->exp = 0;	/* make sure */
+		return;
+	}
+#ifndef USE_DIVIDE
+	/*
+	 * numbers are right shifted one bit to make sure
+	 * that m1 is quaranteed to be larger if its
+	 * maximum bit is set
+	 */
+	b64_rsft(&e1->mantissa);	/* 64 bit shift right */
+	b64_rsft(&e2->mantissa);	/* 64 bit shift right */
+	e1->exp++;
+	e2->exp++;
+#endif
+	/*	check for underflow, divide by zero, etc	*/
+	e1->sign ^= e2->sign;
+	e1->exp -= e2->exp;
+
+#ifndef USE_DIVIDE
+		/* do division of mantissas	*/
+		/* uses partial product method	*/
+		/* init control variables	*/
+
+	count = 64;
+	result.h_32 = 0L;
+	result.l_32 = 0L;
+
+		/* partial product division loop */
+
+	while (count--)	{
+		/* first left shift result 1 bit	*/
+		/* this is ALWAYS done			*/
+
+		b64_lsft(&result);
+
+		/* compare dividend and divisor		*/
+		/* if dividend >= divisor add a bit	*/
+		/* and subtract divisior from dividend	*/
+
+		if ( (e1->m1 < e2->m1) ||
+			((e1->m1 == e2->m1) && (e1->m2 < e2->m2) ))
+			;	/* null statement */
+				/* i.e., don't add or subtract */
+		else	{
+			result.l_32++;	/* ADD	*/
+			if (e2->m2 > e1->m2)
+				e1->m1 -= 1;	/* carry in */
+			e1->m1 -= e2->m1;	/* do SUBTRACTION */
+			e1->m2 -= e2->m2;	/*    SUBTRACTION */
+		}
+
+		/*	shift dividend left one bit OR	*/
+		/*	IF it equals ZERO we can break out	*/
+		/*	of the loop, but still must shift	*/
+		/*	the quotient the remaining count bits	*/
+		/* NB	save the results of this test in error	*/
+		/*	if not zero, then the result is inexact. */
+		/* 	this would be reported in IEEE standard	*/
+
+		/*	lp points to dividend			*/
+		lp = &e1->m1;
+
+		error = ((*lp | *(lp+1)) != 0L) ? 1 : 0;
+		if (error)	{	/* more work */
+			/*	assume max bit == 0 (see above)	*/
+			b64_lsft(&e1->mantissa);
+			continue;
+		}
+		else
+			break;	/* leave loop	*/
+	}	/* end of divide by subtraction loop	*/
+
+	if (count > 0)	{
+		lp = &result.h_32;
+		if (count > 31) {	/* move to higher word */
+			*lp = *(lp+1);
+			count -= 32;
+			*(lp+1) = 0L;	/* clear low word	*/
+		}
+		if (*lp)
+			*lp <<= count;	/* shift rest of way	*/
+		lp++;	/*  == &result.l_32	*/
+		if (*lp) {
+			result.h_32 |= (*lp >> 32-count);
+			*lp <<= count;
+		}
+	}
+#else /* USE_DIVIDE */
+
+	u[4] = (e1->m2 & 1) << 15;
+	b64_rsft(&(e1->mantissa));
+	u[0] = e1->m1 >> 16;
+	u[1] = e1->m1;
+	u[2] = e1->m2 >> 16;
+	u[3] = e1->m2;
+	u[5] = 0; u[6] = 0; u[7] = 0;
+	v[1] = e2->m1 >> 16;
+	v[2] = e2->m1;
+	v[3] = e2->m2 >> 16;
+	v[4] = e2->m2;
+	while (! v[maxv]) maxv--;
+	result.h_32 = 0;
+	result.l_32 = 0;
+	lp = &result.h_32;
+
+	/*
+	 * Use an algorithm of Knuth (The art of programming, Seminumerical
+	 * algorithms), to divide u by v. u and v are both seen as numbers
+	 * with base 65536. 
+	 */
+	for (j = 0; j <= 3; j++, u_p++) {
+		unsigned long q_est, temp;
+
+		if (j == 2) lp++;
+		if (u_p[0] == 0 && u_p[1] < v[1]) continue;
+		temp = ((unsigned long)u_p[0] << 16) + u_p[1];
+		if (u_p[0] >= v[1]) {
+			q_est = 0x0000FFFFL;
+		}
+		else {
+			q_est = temp / v[1];
+		}
+		temp -= q_est * v[1];
+		while (temp < 0x10000 && v[2]*q_est > ((temp<<16)+u_p[2])) {
+			q_est--;
+			temp += v[1];
+		}
+		/*	Now, according to Knuth, we have an estimate of the
+			quotient, that is either correct or one too big, but
+			almost always correct.
+		*/
+		if (q_est != 0)  {
+			int i;
+			unsigned long k = 0;
+			int borrow = 0;
+
+			for (i = maxv; i > 0; i--) {
+				unsigned long tmp = q_est * v[i] + k + borrow;
+				unsigned short md = tmp;
+
+				borrow = (md > u_p[i]);
+				u_p[i] -= md;
+				k = tmp >> 16;
+			}
+			k += borrow;
+			borrow = u_p[0] < k;
+			u_p[0] -= k;
+
+			if (borrow) {
+				/* So, this does not happen often; the estimate
+				   was one too big; correct this
+				*/
+				*lp |= (j & 1) ? (q_est - 1) : ((q_est-1)<<16);
+				borrow = 0;
+				for (i = maxv; i > 0; i--) {
+					unsigned long tmp 
+					    = v[i]+(unsigned long)u_p[i]+borrow;
+					
+					u_p[i] = tmp;
+					borrow = tmp >> 16;
+				}
+				u_p[0] += borrow;
+			}
+			else *lp |= (j & 1) ? q_est : (q_est<<16);
+		}
+	}
+#ifdef	EXCEPTION_INEXACT
+	u_p = &u[0];
+	for (j = 7; j >= 0; j--) {
+		if (*u_p++) {
+			error = 1;
+			break;
+		}
+	}
+#endif
+#endif
+
+#ifdef  EXCEPTION_INEXACT
+        if (error)      {
+                /*
+                 * report here exception 8.5 - Inexact
+                 * from Draft 8.0 of IEEE P754:
+                 * In the absence of an invalid operation exception,
+                 * if the rounded result of an operation is not exact or if
+                 * it overflows without a trap, then the inexact exception
+                 * shall be assigned. The rounded or overflowed result
+                 * shall be delivered to the destination.
+                 */
+                INEXACT();
+#endif
+	e1->mantissa = result;
+
+	nrm_ext(e1);
+	if (e1->exp < EXT_MIN)	{
+		/*
+		 * Exception 8.4 - Underflow
+		 */
+		trap(EFUNFL);	/* underflow */
+		e1->exp = EXT_MIN;
+		e1->m1 = e1->m2 = 0L;
+		return;
+	}
+	if (e1->exp >= EXT_MAX) {
+                /*
+                 * Exception 8.3 - Overflow
+                 */
+                trap(EFOVFL);   /* overflow */
+                e1->exp = EXT_MAX;
+                e1->m1 = e1->m2 = 0L;
+                return;
+        }
+}

lib/ack/fphook/dvf4.fc

@@ -0,0 +1,26 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	DIVIDE TWO SINGLES - SINGLE Precision (dvf 4)
+*/
+
+#include	"FP_types.h"
+
+void
+dvf4(s2,s1)
+SINGLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	extend(&s1,&e1,sizeof(SINGLE));
+	extend(&s2,&e2,sizeof(SINGLE));
+
+		/* do a divide */
+	div_ext(&e1,&e2);
+	compact(&e1,&s1,sizeof(SINGLE));
+}

lib/ack/fphook/dvf8.fc

@@ -0,0 +1,26 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	DIVIDE TWO FLOATS - DOUBLE Precision (DVF 8)
+*/
+
+#include	"FP_types.h"
+
+void
+dvf8(s2,s1)
+DOUBLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	extend(&s1.d[0],&e1,sizeof(DOUBLE));
+	extend(&s2.d[0],&e2,sizeof(DOUBLE));
+
+		/* do a divide */
+	div_ext(&e1,&e2);
+	compact(&e1,&s1.d[0],sizeof(DOUBLE));
+}

lib/ack/fphook/extend.fc

@@ -0,0 +1,111 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	CONVERTS FLOATING POINT TO EXTENDED FORMAT
+
+	Two sizes of FLOATING Point are known:
+		SINGLE and DOUBLE
+*/
+/********************************************************/
+/*
+	It is not required to normalize in extended
+	format, but it has been chosen to do so.
+	Extended Format is as follows (at exit):
+
+->sign	S000 0000 | 0000 0000		<SIGN>
+->exp	0EEE EEEE | EEEE EEEE		<EXPONENT>
+->m1	LFFF FFFF | FFFF FFFF		<L.Fraction>
+	FFFF FFFF | FFFF FFFF		<Fraction>
+->m2	FFFF FFFF | FFFF FFFF		<Fraction>
+	FFFF F000 | 0000 0000		<Fraction>
+*/
+/********************************************************/
+
+#include "FP_bias.h"
+#include "FP_shift.h"
+#include "FP_types.h"
+#include "get_put.h"
+/********************************************************/
+
+void
+extend(from,to,size)
+unsigned long	*from;
+EXTEND	*to;
+int	size;
+{
+	register char *cpt1;
+	unsigned long	tmp;
+	int	leadbit = 0;
+
+	cpt1 = (char *) from;
+
+#if FL_MSL_AT_LOW_ADDRESS
+#if FL_MSW_AT_LOW_ADDRESS
+	to->exp = uget2(cpt1);
+#else
+	to->exp = uget2(cpt1+2);
+#endif
+#else
+#if FL_MSW_AT_LOW_ADDRESS
+	to->exp = uget2(cpt1+(size == sizeof(DOUBLE) ? 4 : 0));
+#else
+	to->exp = uget2(cpt1+(size == sizeof(DOUBLE) ? 6 : 2));
+#endif
+#endif
+	to->sign = (to->exp & 0x8000);	/* set sign bit */
+	to->exp ^= to->sign;
+	if (size == sizeof(DOUBLE))
+		to->exp >>= DBL_EXPSHIFT;
+	else
+		to->exp >>= SGL_EXPSHIFT;
+	if (to->exp > 0)
+		leadbit++;	/* will set Lead bit later	*/
+	else to->exp++;
+
+	if (size == sizeof(DOUBLE))	{
+#if FL_MSL_AT_LOW_ADDRESS
+		to->m1 = get4(cpt1);
+		cpt1 += 4;
+		tmp = get4(cpt1);
+#else
+		tmp = get4(cpt1);
+		cpt1 += 4;
+		to->m1 = get4(cpt1);
+#endif
+		if (to->exp == 1 && to->m1 == 0 && tmp == 0) {
+			to->exp = 0;
+			to->sign = 0;
+			to->m1 = 0;
+			to->m2 = 0;
+			return;
+		}
+		to->m1 <<= DBL_M1LEFT;		/* shift	*/
+		to->exp -= DBL_BIAS;		/* remove bias	*/
+		to->m1 |= (tmp>>DBL_RPACK);	/* plus 10 == 32	*/
+		to->m2 = (tmp<<DBL_LPACK);	/* plus 22 == 32	*/
+	}
+	else	{	/* size == sizeof(SINGLE)		*/
+		to->m1 = get4(cpt1);
+		to->m1  <<= SGL_M1LEFT;	/* shift	*/
+		if (to->exp == 1 && to->m1 == 0) {
+			to->exp = 0;
+			to->sign = 0;
+			to->m1 = 0;
+			to->m2 = 0;
+			return;
+		}
+		to->exp -= SGL_BIAS;		/* remove bias	*/
+		to->m2 = 0L;
+	}
+
+	to->m1 |= NORMBIT;				/* set bit L	*/
+	if (leadbit == 0) {		/* set or clear Leading Bit	*/
+		to->m1 &= ~NORMBIT;			/* clear bit L	*/
+		nrm_ext(to);				/* and normalize */
+	}
+}

lib/ack/fphook/fef4.fc

@@ -0,0 +1,33 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SEPERATE INTO EXPONENT AND FRACTION (FEF 4)
+*/
+
+#include	"FP_types.h"
+
+void
+fef4(r,s1)
+SINGLE	s1;
+struct fef4_returns	*r;
+{
+	EXTEND	buf;
+	register struct fef4_returns	*p = r;	/* make copy; r might refer
+						   to itself (see table)
+						*/
+
+	extend(&s1,&buf,sizeof(SINGLE));
+	if (buf.exp == 0 && buf.m1 == 0 && buf.m2 == 0) {
+		p->e = 0;
+	}
+	else {
+		p->e = buf.exp+1;
+		buf.exp = -1;
+	}
+	compact(&buf,&p->f,sizeof(SINGLE));
+}

lib/ack/fphook/fef8.fc

@@ -0,0 +1,33 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SEPERATE DOUBLE INTO EXPONENT AND FRACTION (FEF 8)
+*/
+
+#include	"FP_types.h"
+
+void
+fef8(r, s1)
+DOUBLE	s1;
+struct fef8_returns *r;
+{
+	EXTEND	buf;
+	register struct fef8_returns *p = r;	/* make copy, r might refer
+						   to itself (see table)
+						*/
+
+	extend(&s1.d[0],&buf,sizeof(DOUBLE));
+	if (buf.exp == 0 && buf.m1 == 0 && buf.m2 == 0) {
+		p->e = 0;
+	}
+	else {
+		p->e = buf.exp + 1;
+		buf.exp = -1;
+	}
+	compact(&buf,&p->f.d[0],sizeof(DOUBLE));
+}

lib/ack/fphook/fif4.fc

@@ -0,0 +1,46 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	MULTIPLY AND DISMEMBER PARTS (FIF 4)
+*/
+
+#include "FP_types.h"
+#include "FP_shift.h"
+
+void
+fif4(p,x,y)
+SINGLE	x,y;
+struct fif4_returns *p;
+{
+
+	EXTEND	e1,e2;
+
+	extend(&y,&e1,sizeof(SINGLE));
+	extend(&x,&e2,sizeof(SINGLE));
+		/* do a multiply */
+	mul_ext(&e1,&e2);
+	e2 = e1;
+	compact(&e2,&y,sizeof(SINGLE));
+	if (e1.exp < 0) {
+		p->ipart = 0;
+		p->fpart = y;
+		return;
+	}
+	if (e1.exp > 30 - SGL_M1LEFT) {
+		p->ipart = y;
+		p->fpart = 0;
+		return;
+	}
+	b64_sft(&e1.mantissa, 63 - e1.exp);
+	b64_sft(&e1.mantissa, e1.exp - 63);	/* "loose" low order bits */
+	compact(&e1,&(p->ipart),sizeof(SINGLE));
+	extend(&(p->ipart), &e2, sizeof(SINGLE));
+	extend(&y, &e1, sizeof(SINGLE));
+	sub_ext(&e1, &e2);
+	compact(&e1, &(p->fpart), sizeof(SINGLE));
+}

lib/ack/fphook/fif8.fc

@@ -0,0 +1,48 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	MULTIPLY AND DISMEMBER PARTS (FIF 8)
+*/
+
+#include "FP_types.h"
+#include "FP_shift.h"
+
+void
+fif8(p,x,y)
+DOUBLE	x,y;
+struct fif8_returns *p;
+{
+
+	EXTEND	e1,e2;
+
+	extend(&y.d[0],&e1,sizeof(DOUBLE));
+	extend(&x.d[0],&e2,sizeof(DOUBLE));
+		/* do a multiply */
+	mul_ext(&e1,&e2);
+	e2 = e1;
+	compact(&e2, &y.d[0], sizeof(DOUBLE));
+	if (e1.exp < 0) {
+		p->ipart.d[0] = 0;
+		p->ipart.d[1] = 0;
+		p->fpart = y;
+		return;
+	}
+	if (e1.exp > 62 - DBL_M1LEFT) {
+		p->ipart = y;
+		p->fpart.d[0] = 0;
+		p->fpart.d[1] = 0;
+		return;
+	}
+	b64_sft(&e1.mantissa, 63 - e1.exp);
+	b64_sft(&e1.mantissa, e1.exp - 63);	/* "loose" low order bits */
+	compact(&e1, &(p->ipart.d[0]), sizeof(DOUBLE));
+	extend(&(p->ipart.d[0]), &e2, sizeof(DOUBLE));
+	extend(&y.d[0], &e1, sizeof(DOUBLE));
+	sub_ext(&e1, &e2);
+	compact(&e1, &(p->fpart.d[0]), sizeof(DOUBLE));
+}

lib/ack/fphook/fltpr.c

@@ -0,0 +1,12 @@
+#include	<stdio.h>
+#include	<stdlib.h>
+#include	"../stdio/loc_incl.h"
+
+int _fp_hook = 1;
+
+char *
+_f_print(va_list *ap, int flags, char *s, char c, int precision)
+{
+	fprintf(stderr,"cannot print floating point\n");
+	exit(EXIT_FAILURE);
+}

lib/ack/fphook/fphook.c

@@ -0,0 +1,195 @@
+/*
+ * fltpr.c - print floating point numbers
+ */
+/* $Header$ */
+
+#ifndef	NOFLOAT
+#include	<string.h>
+#include	<stdarg.h>
+#include	"../stdio/loc_incl.h"
+int _fp_hook = 1;
+
+static char *
+_pfloat(long double r, register char *s, int n, int flags)
+{
+	register char *s1;
+	int sign, dp;
+	register int i;
+
+	s1 = _fcvt(r, n, &dp, &sign);
+	if (sign)
+		*s++ = '-';
+	else if (flags & FL_SIGN)
+		*s++ = '+';
+	else if (flags & FL_SPACE)
+		*s++ = ' ';
+
+	if (dp<=0)
+		*s++ = '0';
+	for (i=dp; i>0; i--)
+		if (*s1) *s++ = *s1++;
+		else *s++ = '0';
+	if (((i=n) > 0) || (flags & FL_ALT))
+		*s++ = '.';
+	while (++dp <= 0) {
+		if (--i<0)
+			break;
+		*s++ = '0';
+	}
+	while (--i >= 0)
+		if (*s1) *s++ = *s1++;
+		else *s++ = '0';
+	return s;
+}
+
+static char *
+_pscien(long double r, register char *s, int n, int flags)
+{
+	int sign, dp; 
+	register char *s1;
+
+	s1 = _ecvt(r, n + 1, &dp, &sign);
+	if (sign)
+		*s++ = '-';
+	else if (flags & FL_SIGN)
+		*s++ = '+';
+	else if (flags & FL_SPACE)
+		*s++ = ' ';
+
+	*s++ = *s1++;
+	if ((n > 0) || (flags & FL_ALT))
+		*s++ = '.';
+	while (--n >= 0)
+		if (*s1) *s++ = *s1++;
+		else *s++ = '0';
+	*s++ = 'e';
+	if ( r != 0 ) --dp ;
+	if ( dp<0 ) {
+		*s++ = '-' ; dp= -dp ;
+	} else {
+		*s++ = '+' ;
+	}
+	if (dp >= 100) {
+		*s++ = '0' + (dp / 100);
+		dp %= 100;
+	}
+	*s++ = '0' + (dp/10);
+	*s++ = '0' + (dp%10);
+	return s;
+}
+
+#define	NDIGINEXP(exp)		(((exp) >= 100 || (exp) <= -100) ? 3 : 2)
+#define	LOW_EXP			-4
+#define	USE_EXP(exp, ndigits)	(((exp) < LOW_EXP + 1) || (exp >= ndigits + 1))
+
+static char *
+_gcvt(long double value, int ndigit, char *s, int flags)
+{
+	int sign, dp;
+	register char *s1, *s2;
+	register int i;
+	register int nndigit = ndigit;
+
+	s1 = _ecvt(value, ndigit, &dp, &sign);
+	s2 = s;
+	if (sign) *s2++ = '-';
+	else if (flags & FL_SIGN)
+		*s2++ = '+';
+	else if (flags & FL_SPACE)
+		*s2++ = ' ';
+
+	if (!(flags & FL_ALT))
+		for (i = nndigit - 1; i > 0 && s1[i] == '0'; i--)
+			nndigit--;
+
+	if (USE_EXP(dp,ndigit))	{
+		/* Use E format */
+		dp--;
+		*s2++ = *s1++;
+		if ((nndigit > 1) || (flags & FL_ALT)) *s2++ = '.';
+		while (--nndigit > 0) *s2++ = *s1++;
+		*s2++ = 'e';
+		if (dp < 0) {
+			*s2++ = '-';
+			dp = -dp;
+		}
+		else	 *s2++ = '+';
+		s2 += NDIGINEXP(dp);
+		*s2 = 0;
+		for (i = NDIGINEXP(dp); i > 0; i--) {
+			*--s2 = dp % 10 + '0';
+			dp /= 10;
+		}
+		return s;
+	}
+	/* Use f format */
+	if (dp <= 0) {
+		if (*s1 != '0')	{
+			/* otherwise the whole number is 0 */
+			*s2++ = '0';
+			*s2++ = '.';
+		}
+		while (dp < 0) {
+			dp++;
+			*s2++ = '0';
+		}
+	}
+	for (i = 1; i <= nndigit; i++) {
+		*s2++ = *s1++;
+		if (i == dp) *s2++ = '.';
+	}
+	if (i <= dp) {
+		while (i++ <= dp) *s2++ = '0';
+		*s2++ = '.';
+	}
+	if ((s2[-1]=='.') && !(flags & FL_ALT)) s2--;
+	*s2 = '\0';
+	return s;
+}
+
+char *
+_f_print(va_list *ap, int flags, char *s, char c, int precision)
+{
+	register char *old_s = s;
+	long double ld_val;
+
+	if (flags & FL_LONGDOUBLE) ld_val = va_arg(*ap, long double);
+	else ld_val = (long double) va_arg(*ap, double);
+
+	switch(c) {
+	case 'f':
+		s = _pfloat(ld_val, s, precision, flags);
+		break;
+	case 'e':
+	case 'E':
+		s = _pscien(ld_val, s, precision , flags);
+		break;
+	case 'g':
+	case 'G':
+		s = _gcvt(ld_val, precision, s, flags);
+		s += strlen(s);
+		break;
+	}
+	if ( c == 'E' || c == 'G') {
+		while (*old_s && *old_s != 'e') old_s++;
+		if (*old_s == 'e') *old_s = 'E';
+	}
+	return s;
+}
+#endif	/* NOFLOAT */
+/* $Header$ */
+
+#include <stdlib.h>
+#include "../ansi/ext_fmt.h"
+
+void _str_ext_cvt(const char *s, char **ss, struct EXTEND *e);
+double _ext_dbl_cvt(struct EXTEND *e);
+
+double
+strtod(const char *p, char **pp)
+{
+	struct EXTEND e;
+
+	_str_ext_cvt(p, pp, &e);
+	return _ext_dbl_cvt(&e);
+}

lib/ack/fphook/fptrp.s

@@ -0,0 +1,19 @@
+#
+.sect .text; .sect .rom; .sect .data; .sect .bss
+.define __fptrp
+.sect .text
+__fptrp:
+#if __i386
+	push	ebp
+	mov	ebp, esp
+	mov	eax, 8(bp)
+	call	.Xtrp
+	leave
+	ret
+#else /* i86 */
+	push	bp
+	mov	bp, sp
+	mov	ax, 4(bp)
+	call	.Xtrp
+	jmp	.cret
+#endif

lib/ack/fphook/get_put.h

@@ -0,0 +1,41 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+#include <byte_order.h>
+
+#if CHAR_UNSIGNED
+#define Xchar(ch)	(ch)
+#else
+#define Xchar(ch)	((ch) & 0377)
+#endif
+
+#define BYTES_REVERSED (MSB_AT_LOW_ADDRESS != FL_MSB_AT_LOW_ADDRESS)
+#define WORDS_REVERSED (MSW_AT_LOW_ADDRESS != FL_MSW_AT_LOW_ADDRESS)
+#define LONGS_REVERSED (FL_MSL_AT_LOW_ADDRESS)
+
+#if BYTES_REVERSED
+#define uget2(c)	(Xchar((c)[1]) | ((unsigned) Xchar((c)[0]) << 8))
+#define Xput2(i, c)	(((c)[1] = (i)), ((c)[0] = (i) >> 8))
+#define put2(i, c)	{ register int j = (i); Xput2(j, c); }
+#else
+#define uget2(c)	(* ((unsigned short *) (c)))
+#define Xput2(i, c)	(* ((short *) (c)) = (i))
+#define put2(i, c)	Xput2(i, c)
+#endif
+
+#define get2(c)		((short) uget2(c))
+
+#if WORDS_REVERSED || BYTES_REVERSED
+#define get4(c)		(uget2((c)+2) | ((long) uget2(c) << 16))
+#define put4(l, c)	{ register long x=(l); \
+			  Xput2((int)x,(c)+2); \
+			  Xput2((int)(x>>16),(c)); \
+			}
+#else
+#define get4(c)		(* ((long *) (c)))
+#define put4(l, c)	(* ((long *) (c)) = (l))
+#endif

lib/ack/fphook/mlf4.fc

@@ -0,0 +1,25 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+ * Multiply Single Precesion Float (MLF 4)
+ */
+
+#include	"FP_types.h"
+
+void
+mlf4(s2,s1)
+SINGLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	extend(&s1,&e1,sizeof(SINGLE));
+	extend(&s2,&e2,sizeof(SINGLE));
+		/* do a multiply */
+	mul_ext(&e1,&e2);
+	compact(&e1,&s1,sizeof(SINGLE));
+}

lib/ack/fphook/mlf8.fc

@@ -0,0 +1,25 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+ * Multiply Double Precision Float (MLF 8)
+ */
+
+#include	"FP_types.h"
+
+void
+mlf8(s2,s1)
+DOUBLE	s1,s2;
+{
+	EXTEND	e1,e2;
+
+	extend(&s1.d[0],&e1,sizeof(DOUBLE));
+	extend(&s2.d[0],&e2,sizeof(DOUBLE));
+		/* do a multiply */
+	mul_ext(&e1,&e2);
+	compact(&e1,&s1.d[0],sizeof(DOUBLE));
+}

lib/ack/fphook/mul_ext.fc

@@ -0,0 +1,98 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	ROUTINE TO MULTIPLY TWO EXTENDED FORMAT NUMBERS
+*/
+
+# include "FP_bias.h"
+# include "FP_trap.h"
+# include "FP_types.h"
+# include "FP_shift.h"
+
+void
+mul_ext(e1,e2)
+EXTEND	*e1,*e2;
+{
+	register int	i,j;		/* loop control	*/
+	unsigned short	mp[4];		/* multiplier */
+	unsigned short	mc[4];		/* multipcand */
+	unsigned short	result[8];	/* result */
+	register unsigned short *pres;
+
+	/* first save the sign (XOR)			*/
+	e1->sign ^= e2->sign;
+
+	/* compute new exponent */
+	e1->exp += e2->exp + 1;
+	/* 128 bit multiply of mantissas			*/
+
+		/* assign unknown long formats		*/
+		/* to known unsigned word formats	*/
+	mp[0] = e1->m1 >> 16;
+	mp[1] = (unsigned short) e1->m1;
+	mp[2] = e1->m2 >> 16;
+	mp[3] = (unsigned short) e1->m2;
+	mc[0] = e2->m1 >> 16;
+	mc[1] = (unsigned short) e2->m1;
+	mc[2] = e2->m2 >> 16;
+	mc[3] = (unsigned short) e2->m2;
+	for (i = 8; i--;) {
+		result[i] = 0;
+	}
+	/*
+	 *	fill registers with their components
+	 */
+	for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
+		unsigned short k = 0;
+		unsigned long mpi = mp[i];
+		for(j=4;j--;) {
+			unsigned long tmp = (unsigned long)pres[j] + k;
+			if (mc[j]) tmp += mpi * mc[j];
+			pres[j] = tmp;
+			k = tmp >> 16;
+		}
+		pres[-1] = k;
+	}
+        if (! (result[0] & 0x8000)) {
+                e1->exp--;
+                for (i = 0; i <= 3; i++) {
+                        result[i] <<= 1;
+                        if (result[i+1]&0x8000) result[i] |= 1;
+                }
+                result[4] <<= 1;
+        }
+
+	/*
+	 *	combine the registers to a total
+	 */
+	e1->m1 = ((unsigned long)(result[0]) << 16) + result[1];
+	e1->m2 = ((unsigned long)(result[2]) << 16) + result[3];
+	if (result[4] & 0x8000) {
+		if (++e1->m2 == 0)
+			if (++e1->m1 == 0) {
+				e1->m1 = NORMBIT;
+				e1->exp++;
+			}
+	}
+
+					/* check for overflow	*/
+	if (e1->exp >= EXT_MAX)	{
+		trap(EFOVFL);
+			/* if caught 			*/
+			/* return signed infinity	*/
+		e1->exp = EXT_MAX;
+infinity:	e1->m1 = e1->m2 =0L;
+		return;
+	}
+				/* check for underflow	*/
+	if (e1->exp < EXT_MIN)	{
+		trap(EFUNFL);
+		e1->exp = EXT_MIN;
+		goto infinity;
+	}
+}

lib/ack/fphook/ngf4.fc

@@ -0,0 +1,27 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		NEGATE A FLOATING POINT (NGF 4)
+*/
+/********************************************************/
+
+#include "FP_types.h"
+#include "get_put.h"
+
+#define OFF ((FL_MSW_AT_LOW_ADDRESS ? 0 : 2) + (FL_MSB_AT_LOW_ADDRESS ? 0 : 1))
+void
+ngf4(f)
+SINGLE	f;
+{
+	unsigned char *p;
+
+	if (f != (SINGLE) 0) {
+		p = (unsigned char *) &f + OFF;
+		*p ^= 0x80;
+	}
+}

lib/ack/fphook/ngf8.fc

@@ -0,0 +1,28 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+		NEGATE A FLOATING POINT (NGF 8)
+*/
+/********************************************************/
+
+#include "FP_types.h"
+#include "get_put.h"
+
+#define OFF ((FL_MSL_AT_LOW_ADDRESS ? 0 : 4) + (FL_MSW_AT_LOW_ADDRESS ? 0 : 2) + (FL_MSB_AT_LOW_ADDRESS ? 0 : 1))
+
+void
+ngf8(f)
+DOUBLE	f;
+{
+	unsigned char	*p;
+
+	if (f.d[0] != 0 || f.d[1] != 0) {
+		p = (unsigned char *) &f + OFF;
+		*p ^= 0x80;
+	}
+}

lib/ack/fphook/nrm_ext.fc

@@ -0,0 +1,50 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/********************************************************/
+/*
+	NORMALIZE an EXTENDED FORMAT NUMBER
+*/
+/********************************************************/
+
+#include "FP_shift.h"
+#include "FP_types.h"
+
+void
+nrm_ext(e1)
+EXTEND	*e1;
+{
+		/* we assume that the mantissa != 0	*/
+		/* if it is then just return		*/
+		/* to let it be a problem elsewhere	*/
+		/* THAT IS, The exponent is not set to	*/
+		/* zero. If we don't test here an	*/
+		/* infinite loop is generated when	*/
+		/* mantissa is zero			*/
+
+	if ((e1->m1 | e1->m2) == 0L)
+		return;
+
+		/* if top word is zero mov low word	*/
+		/* to top word, adjust exponent value	*/
+	if (e1->m1 == 0L)	{
+		e1->m1 = e1->m2;
+		e1->m2 = 0L;
+		e1->exp -= 32;
+	}
+	if ((e1->m1 & NORMBIT) == 0) {
+		unsigned long l = ((unsigned long)NORMBIT >> 1);
+		int cnt = -1;
+
+		while (! (l & e1->m1)) {
+			l >>= 1;
+			cnt--;
+		}
+		e1->exp += cnt;
+		b64_sft(&(e1->mantissa), cnt);
+	}
+}

lib/ack/fphook/sbf4.fc

@@ -0,0 +1,27 @@
+/*
+ (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SUBTRACT TWO FLOATS - SINGLE Precision (SBF 4)
+*/
+
+#include	"FP_types.h"
+
+void
+sbf4(s2,s1)
+SINGLE	s1,s2;
+{
+	EXTEND e1,e2;
+
+	if (s2 == (SINGLE) 0) {
+		return;
+	}
+	extend(&s1,&e1,sizeof(SINGLE));
+	extend(&s2,&e2,sizeof(SINGLE));
+	sub_ext(&e1,&e2);
+	compact(&e1,&s1,sizeof(SINGLE));
+}

lib/ack/fphook/sbf8.fc

@@ -0,0 +1,27 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SUBTRACT TWO FLOATS - DOUBLE Precision (SBF 8)
+*/
+
+#include	"FP_types.h"
+
+void
+sbf8(s2,s1)
+DOUBLE	s1,s2;
+{
+	EXTEND e1, e2;
+
+	if (s2.d[0] == 0 && s2.d[1] == 0) {
+		return;
+	}
+	extend(&s1.d[0],&e1,sizeof(DOUBLE));
+	extend(&s2.d[0],&e2,sizeof(DOUBLE));
+	sub_ext(&e1,&e2);
+	compact(&e1,&s1.d[0],sizeof(DOUBLE));
+}

lib/ack/fphook/sft_ext.fc

@@ -0,0 +1,39 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SHIFT TWO EXTENDED NUMBERS INTO PROPER
+	ALIGNMENT FOR ADDITION (exponents are equal)
+	Numbers should not be zero on entry.
+*/
+
+#include "FP_types.h"
+
+void
+sft_ext(e1,e2)
+EXTEND	*e1,*e2;
+{
+	register	EXTEND	*s;
+	register	int	diff;
+
+	diff = e1->exp - e2->exp;
+
+	if (!diff)
+		return;	/* exponents are equal	*/
+
+	if (diff < 0)	{ /* e2 is larger	*/
+			/* shift e1		*/
+		diff = -diff;
+		s = e1;
+	}
+	else		/* e1 is larger		*/
+			/* shift e2		*/
+		s = e2;
+
+	s->exp += diff;
+	b64_sft(&(s->mantissa), diff);
+}

lib/ack/fphook/shifter.fc

@@ -0,0 +1,75 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+# include "FP_types.h"
+
+void
+b64_sft(e1,n)
+B64	*e1;
+int	n;
+{
+	if (n > 0) {
+		if (n > 63) {
+			e1->l_32 = 0;
+			e1->h_32 = 0;
+			return;
+		}
+		if (n >= 32) {
+			e1->l_32 = e1->h_32;
+			e1->h_32 = 0;
+			n -= 32;
+		}
+		if (n > 0) {
+			e1->l_32 >>= n;
+			if (e1->h_32 != 0) {
+				e1->l_32 |= (e1->h_32 << (32 - n));
+				e1->h_32 >>= n;
+			}
+		}
+		return;
+	}
+	n = -n;
+	if (n > 0) {
+		if (n > 63) {
+			e1->l_32 = 0;
+			e1->h_32 = 0;
+			return;
+		}
+		if (n >= 32) {
+			e1->h_32 = e1->l_32;
+			e1->l_32 = 0;
+			n -= 32;
+		}
+		if (n > 0) {
+			e1->h_32 <<= n;
+			if (e1->l_32 != 0) {
+				e1->h_32 |= (e1->l_32 >> (32 - n));
+				e1->l_32 <<= n;
+			}
+		}
+	}
+}
+
+void
+b64_lsft(e1)
+B64	*e1;
+{
+	/*	shift left 1 bit */
+	e1->h_32 <<= 1;
+	if (e1->l_32 & 0x80000000L) e1->h_32 |= 1;
+	e1->l_32 <<= 1;
+}
+
+void
+b64_rsft(e1)
+B64	*e1;
+{
+	/*	shift right 1 bit */
+	e1->l_32 >>= 1;
+	if (e1->h_32 & 1) e1->l_32 |= 0x80000000L;
+	e1->h_32 >>= 1;
+}

lib/ack/fphook/strtod.c

@@ -0,0 +1,9 @@
+#include	<stdio.h>
+#include	<stdlib.h>
+
+double
+strtod(const char *p, char **pp)
+{
+	fprintf(stderr,"cannot print floating point\n");
+	exit(EXIT_FAILURE);
+}

lib/ack/fphook/sub_ext.fc

@@ -0,0 +1,53 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	SUBTRACT 2 EXTENDED FORMAT NUMBERS
+*/
+
+#include "FP_types.h"
+
+void
+sub_ext(e1,e2)
+EXTEND	*e1,*e2;
+{
+	if ((e2->m1 | e2->m2) == 0L) {
+		return;
+	}
+	if ((e1->m1 | e1->m2) == 0L) {
+		*e1 = *e2;
+		e1->sign = e2->sign ? 0 : 1;
+		return;
+	}
+	sft_ext(e1, e2);
+	if (e1->sign != e2->sign) {
+		/* e1 - e2 = e1 + (-e2) */
+		if (b64_add(&e1->mantissa,&e2->mantissa)) { /* addition carry */
+                	b64_rsft(&e1->mantissa);      /* shift mantissa one bit RIGHT */
+                	e1->m1 |= 0x80000000L;  /* set max bit  */
+                	e1->exp++;              /* increase the exponent */
+        	}
+	}
+        else if (e2->m1 > e1->m1 ||
+                 (e2->m1 == e1->m1 && e2->m2 > e1->m2)) {
+		/*	abs(e2) > abs(e1) */
+		if (e1->m2 > e2->m2) {
+			e2->m1 -= 1;	/* carry in */
+		}
+		e2->m1 -= e1->m1;
+		e2->m2 -= e1->m2;
+		*e1 = *e2;
+		e1->sign = e2->sign ? 0 : 1;
+	}
+	else {
+		if (e2->m2 > e1->m2)
+			e1->m1 -= 1;	/* carry in */
+		e1->m1 -= e2->m1;
+		e1->m2 -= e2->m2;
+	}
+	nrm_ext(e1);
+}

lib/ack/fphook/zrf4.fc

@@ -0,0 +1,19 @@
+/*
+  (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
+  See the copyright notice in the ACK home directory, in the file "Copyright".
+*/
+
+/* $Header$ */
+
+/*
+	return a zero float (ZRF 4)
+*/
+
+#include "FP_types.h"
+
+void
+zrf4(l)
+SINGLE	*l;
+{
+	*l = 0L;
+}