--- /dev/null
+FP.script
+FP_bias.h
+FP_shift.h
+FP_trap.h
+FP_types.h
+Makefile
+add_ext.c
+adder.c
+adder.h
+adf4.c
+adf8.c
+cff4.c
+cff8.c
+cfi.c
+cfu.c
+cif4.c
+cif8.c
+cmf4.c
+cmf8.c
+compact.c
+cuf4.c
+cuf8.c
+div_ext.c
+dvf4.c
+dvf8.c
+extend.c
+fef4.c
+fef8.c
+fif4.c
+fif8.c
+fptrp.e
+get_put.h
+mlf4.c
+mlf8.c
+mul_ext.c
+ngf4.c
+ngf8.c
+nrm_ext.c
+prt_dbl.c
+prt_ext.c
+sbf4.c
+sbf8.c
+sft_ext.c
+shifter.c
+sub_ext.c
+zrf4.c
+zrf8.c
+zrf_ext.c
--- /dev/null
+g/_adf4/s//.adf4/g
+g/_adf8/s//.adf8/g
+g/_cff4/s//.cff4/g
+g/_cff8/s//.cff8/g
+g/_cfi/s//.cfi/g
+g/_cfu/s//.cfu/g
+g/_cif4/s//.cif4/g
+g/_cif8/s//.cif8/g
+g/_cmf4/s//.cmf4/g
+g/_cmf8/s//.cmf8/g
+g/_cuf4/s//.cuf4/g
+g/_cuf8/s//.cuf8/g
+g/_dvf4/s//.dvf4/g
+g/_dvf8/s//.dvf8/g
+g/_fef4/s//.fef4/g
+g/_fef8/s//.fef8/g
+g/_fif4/s//.fif4/g
+g/_fif8/s//.fif8/g
+g/_mlf4/s//.mlf4/g
+g/_mlf8/s//.mlf8/g
+g/_ngf4/s//.ngf4/g
+g/_ngf8/s//.ngf8/g
+g/_sbf4/s//.sbf4/g
+g/_sbf8/s//.sbf8/g
+g/_zrf4/s//.zrf4/g
+g/_zrf8/s//.zrf8/g
+g/_add_ext/s//.add_ext/g
+g/_div_ext/s//.div_ext/g
+g/_mul_ext/s//.mul_ext/g
+g/_nrm_ext/s//.nrm_ext/g
+g/_prt_ext/s//.prt_ext/g
+g/_sft_ext/s//.sft_ext/g
+g/_sub_ext/s//.sub_ext/g
+g/_zrf_ext/s//.zrf_ext/g
+g/_compact/s//.compact/g
+g/_extend/s//.extend/g
+g/_load4/s//.load4/g
+g/_store4/s//.store4/g
+g/_b32_add/s//.b32_add/g
+g/_b64_add/s//.b64_add/g
+w
+q
--- /dev/null
+/*
+ include file for floating point package
+*/
+
+#ifdef IEEEFORMAT
+ /* 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 255 /* standard definition */
+#define SGL_MIN 0 /* standard definition */
+#define DBL_MAX 2047 /* standard definition */
+#define DBL_MIN 0 /* standard definition */
+#define EXT_MAX 16384 /* standard minimum */
+#define EXT_MIN -16383 /* standard minimum */
+#else
+
+ /* FLOAT FORMAT EXPONENT BIAS */
+
+#define SGL_BIAS 127 /* excess 128 notation used */
+#define DBL_BIAS 127 /* excess 128 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 255 /* standard definition */
+#define SGL_MIN 0 /* standard definition */
+#define DBL_MAX 255 /* standard definition */
+#define DBL_MIN 0 /* standard definition */
+#define EXT_MAX 16384 /* standard minimum */
+#define EXT_MIN -16383 /* standard minimum */
+#endif
--- /dev/null
+/*
+ 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 */
+#ifndef IEEEFORMAT
+
+#define DBL_EXPSHIFT SGL_EXPSHIFT
+#define DBL_M1LEFT SGL_M1LEFT
+
+#define DBL_LPACK DBL_RUNPACK
+#define DBL_RPACK DBL_LUNPACK
+
+#define DBL_ZERO SGL_ZERO
+#define DBL_EXACT SGL_EXACT
+
+#define DBL_RUNPACK DBL_M1LEFT
+#define DBL_LUNPACK 32-DBL_M1LEFT
+
+#define DBL_ROUNDUP SGL_ROUNDUP
+#define DBL_CARRYOUT SGL_CARRYOUT
+#define DBL_MASK SGL_MASK
+
+#else
+ /* 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
+
+#endif IEEEFORMAT
--- /dev/null
+/*
+ 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)
--- /dev/null
+/********************************************************/
+/*
+ include file for floating point package
+*/
+/********************************************************/
+/*
+ THESE STRUCTURES ARE USED TO ADDRESS THE INDIVIDUAL
+ PARTS OF THE FLOATING NUMBER REPRESENTATIONS.
+
+ THREE STRUCTURES ARE DEFINED:
+ SINGLE: single precision floating format
+ DOUBLE: double precision floating format
+ EXTEND: double precision extended format
+*/
+/********************************************************/
+typedef unsigned long _float;
+
+typedef union {
+ double _dbl;
+ unsigned long __double[2];
+} _double;
+
+typedef union { /* address parts of float */
+ /* field to extract exponent */
+ short sgl_exp[sizeof(float)/sizeof(short)];
+ /* same as fract below */
+ _float f[sizeof(float)/sizeof(_float)];
+ unsigned long fract;
+} SINGLE;
+
+typedef union {
+ /* field to extract exponent */
+ short dbl_exp[sizeof(double)/sizeof(short)];
+ /* for normal syntax use */
+ _double _f8[sizeof(double)/sizeof(_double)];
+ /* a float in a double - for returns */
+ _float _f4[sizeof(double)/sizeof(_float)];
+ /* to be deleted eventually */
+ struct { /* address parts of float */
+ SINGLE p1; /* part one */
+ unsigned long p2; /* part two */
+ } _s;
+} DOUBLE;
+
+typedef struct { /* expanded float format */
+ short sign;
+ short exp;
+ unsigned long m1;
+ unsigned long m2; /* includes guard byte */
+} EXTEND;
+#ifdef PRT_EXT
+#include <stdio.h>
+#endif
--- /dev/null
+EMHOME=../../..
+SUF=s
+MACH=m68k4
+ASAR=arch
+CFLAGS=-O
+# must use -r option of make so that default rules
+# are not loaded
+#
+# $Header$
+#
+# various flags that can be used during compilation
+# define DEBUG
+# define PRT_ADD
+# define PRT_ALL
+# define PRT_DBL
+# define PRT_ENTRY
+# define PRT_EXIT
+# define PRT_EXT
+# define PRT_EXT2
+# define PRT_LONG
+# define PRT_RNDMSG
+# define PRT_STDERR
+# define PRT_TRAP
+#
+# DFLAGS=-DPRT_ADD -DPRT_ALL -DPRT_DBL -DPRT_ENTRY -DPRT_EXIT -DPRT_EXT -DPRT_EXT2 -DPRT_LONG -DPRT_RNDMSG -DPRT_STDERR -DPRT_TRAP
+DFLAGS=
+EMFLAGS= -L -LIB -I. $(DFLAGS) $(CFLAGS)
+# AS=ack -m$(MACH) -c.$(SUF)
+# CC=ack -m$(MACH) -c.s
+# CCFLAGS=$(EMFLAGS)
+CDIR=$(EMHOME)/mach/proto/fp
+
+LIST = cff4.$(SUF) cff8.$(SUF) cfu.$(SUF) cmf4.$(SUF) cmf8.$(SUF)\
+ cuf4.$(SUF) cuf8.$(SUF)\
+ dvf4.$(SUF) dvf8.$(SUF) fef4.$(SUF) fef8.$(SUF)\
+ fif4.$(SUF) fif8.$(SUF)\
+ cfi.$(SUF) cif4.$(SUF) cif8.$(SUF) mlf4.$(SUF) mlf8.$(SUF)\
+ ngf4.$(SUF)\
+ ngf8.$(SUF) sbf4.$(SUF) sbf8.$(SUF) adf4.$(SUF) adf8.$(SUF)\
+ zrf4.$(SUF) zrf8.$(SUF)\
+ extend.$(SUF) compact.$(SUF)\
+ add_ext.$(SUF) div_ext.$(SUF) mul_ext.$(SUF) nrm_ext.$(SUF)\
+ sft_ext.$(SUF) sub_ext.$(SUF) zrf_ext.$(SUF)\
+ adder.$(SUF) shifter.$(SUF)\
+ prt_dbl.$(SUF)\
+ fptrp.$(SUF) prt_ext.$(SUF) # debugging
+SLIST = cff4.s cff8.s cfu.s cmf4.s cmf8.s\
+ cuf4.s cuf8.s\
+ dvf4.s dvf8.s fef4.s fef8.s\
+ fif4.s fif8.s\
+ cfi.s cif4.s cif8.s mlf4.s mlf8.s\
+ ngf4.s\
+ ngf8.s sbf4.s sbf8.s adf4.s adf8.s\
+ zrf4.s zrf8.s\
+ extend.s compact.s\
+ add_ext.s div_ext.s mul_ext.s nrm_ext.s\
+ sft_ext.s sub_ext.s zrf_ext.s\
+ adder.s shifter.s\
+ prt_dbl.s\
+ fptrp.s prt_ext.s # debugging
+
+SRC = FP_bias.h FP_shift.h FP_trap.h FP_types.h adder.h get_put.h\
+ cff4.c cff8.c cfu.c cmf4.c cmf8.c\
+ cuf4.c cuf8.c\
+ dvf4.c dvf8.c fef4.c fef8.c\
+ fif4.c fif8.c\
+ cfi.c cif4.c cif8.c mlf4.c mlf8.c\
+ ngf4.c\
+ ngf8.c sbf4.c sbf8.c adf4.c adf8.c\
+ zrf4.c zrf8.c\
+ extend.c compact.c\
+ add_ext.c div_ext.c mul_ext.c nrm_ext.c\
+ sft_ext.c sub_ext.c zrf_ext.c\
+ adder.c shifter.c\
+ prt_dbl.c\
+ fptrp.e prt_ext.c
+
+all: FP_$(MACH).a
+
+install: tail_fp
+
+tail_fp: FP_$(MACH).a
+ ../../install FP_$(MACH).a tail_fp
+
+clean:
+ rm -f $(LIST) FP_$(MACH).a
+ rm -f $(SLIST)
+
+opr:
+ make pr | opr
+
+pr:
+ @pr Makefile FP.script $(SRC)
+
+FP_$(MACH).a: $(LIST)
+ $(ASAR) rv $@ $?
+
+fptrp.$(SUF): $(CDIR)/fptrp.e
+ ack -m$(MACH) -L -LIB -c $(CDIR)/fptrp.e
+
+extend.$(SUF) compact.$(SUF): byte_order.h $(CDIR)/get_put.h
+
+cff4.$(SUF): $(CDIR)/cff4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cff4.c
+ ed - cff4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cff4.s
+
+cff8.$(SUF): $(CDIR)/cff8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cff8.c
+ ed - cff8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cff8.s
+
+cfu.$(SUF): $(CDIR)/cfu.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cfu.c
+ ed - cfu.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cfu.s
+
+cmf4.$(SUF): $(CDIR)/cmf4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cmf4.c
+ ed - cmf4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cmf4.s
+
+cmf8.$(SUF): $(CDIR)/cmf8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cmf8.c
+ ed - cmf8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cmf8.s
+
+cuf4.$(SUF): $(CDIR)/cuf4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cuf4.c
+ ed - cuf4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cuf4.s
+
+cuf8.$(SUF): $(CDIR)/cuf8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cuf8.c
+ ed - cuf8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cuf8.s
+
+dvf4.$(SUF): $(CDIR)/dvf4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/dvf4.c
+ ed - dvf4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) dvf4.s
+
+dvf8.$(SUF): $(CDIR)/dvf8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/dvf8.c
+ ed - dvf8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) dvf8.s
+
+fef4.$(SUF): $(CDIR)/fef4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/fef4.c
+ ed - fef4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) fef4.s
+
+fef8.$(SUF): $(CDIR)/fef8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/fef8.c
+ ed - fef8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) fef8.s
+
+fif4.$(SUF): $(CDIR)/fif4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/fif4.c
+ ed - fif4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) fif4.s
+
+fif8.$(SUF): $(CDIR)/fif8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/fif8.c
+ ed - fif8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) fif8.s
+
+cfi.$(SUF): $(CDIR)/cfi.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cfi.c
+ ed - cfi.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cfi.s
+
+cif4.$(SUF): $(CDIR)/cif4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cif4.c
+ ed - cif4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cif4.s
+
+cif8.$(SUF): $(CDIR)/cif8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/cif8.c
+ ed - cif8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) cif8.s
+
+mlf4.$(SUF): $(CDIR)/mlf4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/mlf4.c
+ ed - mlf4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) mlf4.s
+
+mlf8.$(SUF): $(CDIR)/mlf8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/mlf8.c
+ ed - mlf8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) mlf8.s
+
+ngf4.$(SUF): $(CDIR)/ngf4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/ngf4.c
+ ed - ngf4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) ngf4.s
+
+ngf8.$(SUF): $(CDIR)/ngf8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/ngf8.c
+ ed - ngf8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) ngf8.s
+
+sbf4.$(SUF): $(CDIR)/sbf4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/sbf4.c
+ ed - sbf4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) sbf4.s
+
+sbf8.$(SUF): $(CDIR)/sbf8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/sbf8.c
+ ed - sbf8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) sbf8.s
+
+adf4.$(SUF): $(CDIR)/adf4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/adf4.c
+ ed - adf4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) adf4.s
+
+adf8.$(SUF): $(CDIR)/adf8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/adf8.c
+ ed - adf8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) adf8.s
+
+zrf4.$(SUF): $(CDIR)/zrf4.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/zrf4.c
+ ed - zrf4.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) zrf4.s
+
+zrf8.$(SUF): $(CDIR)/zrf8.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/zrf8.c
+ ed - zrf8.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) zrf8.s
+
+extend.$(SUF): $(CDIR)/extend.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/extend.c
+ ed - extend.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) extend.s
+
+compact.$(SUF): $(CDIR)/compact.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/compact.c
+ ed - compact.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) compact.s
+
+add_ext.$(SUF): $(CDIR)/add_ext.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/add_ext.c
+ ed - add_ext.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) add_ext.s
+
+div_ext.$(SUF): $(CDIR)/div_ext.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/div_ext.c
+ ed - div_ext.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) div_ext.s
+
+mul_ext.$(SUF): $(CDIR)/mul_ext.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/mul_ext.c
+ ed - mul_ext.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) mul_ext.s
+
+nrm_ext.$(SUF): $(CDIR)/nrm_ext.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/nrm_ext.c
+ ed - nrm_ext.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) nrm_ext.s
+
+sft_ext.$(SUF): $(CDIR)/sft_ext.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/sft_ext.c
+ ed - sft_ext.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) sft_ext.s
+
+sub_ext.$(SUF): $(CDIR)/sub_ext.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/sub_ext.c
+ ed - sub_ext.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) sub_ext.s
+
+zrf_ext.$(SUF): $(CDIR)/zrf_ext.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/zrf_ext.c
+ ed - zrf_ext.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) zrf_ext.s
+
+adder.$(SUF): $(CDIR)/adder.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/adder.c
+ ed - adder.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) adder.s
+
+shifter.$(SUF): $(CDIR)/shifter.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/shifter.c
+ ed - shifter.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) shifter.s
+
+prt_dbl.$(SUF): $(CDIR)/prt_dbl.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/prt_dbl.c
+ ed - prt_dbl.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) prt_dbl.s
+
+prt_ext.$(SUF): $(CDIR)/prt_ext.c
+ ack -c.s -m$(MACH) $(EMFLAGS) $(CDIR)/prt_ext.c
+ ed - prt_ext.s <$(CDIR)/FP.script
+ ack -c -m$(MACH) $(EMFLAGS) prt_ext.s
--- /dev/null
+/*
+#define PRT_EXT
+ ADD TWO EXTENDED FORMAT NUMBERS
+*/
+
+#include "FP_types.h"
+
+add_ext(e1,e2)
+register EXTEND *e1,*e2;
+{
+#ifdef PRT_EXT
+ prt_ext("before ADD #1:",e1);
+ prt_ext("before ADD #2:",e2);
+#endif PRT_EXT
+ if (b64_add(&e1->m1,&e2->m1)) { /* addition carry */
+ b64_sft(&e1->m1,1); /* shift mantissa one bit RIGHT */
+ e1->m1 |= 0x80000000L; /* set max bit */
+ e1->exp++; /* increase the exponent */
+ }
+#ifdef PRT_EXT
+ prt_ext("AFTER ADD :",e1);
+#endif
+ nrm_ext(e1);
+#ifdef PRT_EXT
+ prt_ext("AFTER NRM :",e1);
+#endif
+}
--- /dev/null
+/*
+ * these are the routines the routines to do 32 and 64-bit addition
+ */
+
+# ifdef DEBUG
+# include <stdio.h>
+# endif
+
+# include "adder.h"
+# define UNKNOWN -1
+# define TRUE 1
+# define FALSE 0
+# define MAXBIT 0x80000000L
+
+ /*
+ * add 64 bits
+ */
+b64_add(e1,e2)
+ /*
+ * pointers to 64 bit 'registers'
+ */
+register B64 *e1,*e2;
+{
+ register short overflow;
+ short carry;
+
+ /* add higher pair of 32 bits */
+ overflow = b32_add(&e1->h_32,&e2->h_32);
+
+ /* add lower pair of 32 bits */
+ carry = b32_add(&e1->l_32,&e2->l_32);
+# ifdef 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 */
+ else
+ return(overflow); /* return status from higher add */
+}
+
+ /*
+ * add 32 bits (unsigned longs)
+ * and return the carry status
+ */
+
+b32_add(e1,e2)
+register unsigned long *e1,*e2;
+{
+ register short carry;
+
+ if (*e1 & *e2 & MAXBIT) /* both max_bits are set */
+ carry = TRUE; /* so there is a carry */
+ else
+ carry = ((*e1 | *e2) & MAXBIT)
+ /* only one is set - might be a carry */
+ ? UNKNOWN
+ /* both are clear - no carry */
+ : FALSE;
+# ifdef DEBUG
+ fflush(stdout);
+ printf("\t\t\t\t\tb32_add: overflow before add(%d) test(%d)\n",
+ carry,(*e1&MAXBIT)?FALSE:TRUE);
+ printf("%08X\n%08X\n",*e1,*e2);
+# endif
+
+ *e1 += *e2;
+# ifdef DEBUG
+ printf("%08X\n",*e1);
+ fflush(stdout);
+# endif
+ if (carry != UNKNOWN)
+ return(carry);
+ else
+ /*
+ * if maxbit in answer is set there is no carry
+ * return the NAND of this bit
+ */
+ return((*e1&MAXBIT)?FALSE:TRUE);
+}
--- /dev/null
+/*
+ * include file for 32 & 64 bit addition
+ */
+
+typedef struct {
+ unsigned long h_32; /* higher 32 bits of 64 */
+ unsigned long l_32; /* lower 32 bits of 64 */
+} B64;
--- /dev/null
+/*
+ ADD TWO FLOATS - SINGLE
+*/
+
+#include "FP_types.h"
+
+_float
+adf4(s2,s1)
+_float s1,s2;
+{
+ EXTEND e1,e2;
+ int swap = 0;
+
+ extend((_double *)&s1,&e1,sizeof(SINGLE));
+ extend((_double *)&s2,&e2,sizeof(SINGLE));
+ /* if signs differ do subtraction */
+ /* result in e1 */
+ if (e1.sign ^ e2.sign) {
+#ifdef DEBUG
+ printf("\t\t\tADF calls SUBTRACT\n");
+#endif DEBUG
+ /* set sign of e1 to sign of largest */
+ swap = (e2.exp > e1.exp) ? 1 : (e2.exp < e1.exp) ? 0 :
+ (e2.m1 > e1.m1 ) ? 1 : 0;
+ /* adjust mantissas to equal powers */
+ sft_ext(&e1,&e2);
+ /* subtract the extended formats */
+ if (swap) {
+#ifdef DEBUG
+ printf("\t\t\t\tSWAP\n");
+#endif DEBUG
+ sub_ext(&e2,&e1);
+ e1 = e2;
+ }
+ else
+ sub_ext(&e1,&e2);
+ }
+ else {
+#ifdef DEBUG
+ printf("\t\t\tADF calls ADDITION\n");
+#endif DEBUG
+ /* adjust mantissas to equal powers */
+ sft_ext(&e1,&e2);
+ add_ext(&e1,&e2);
+ }
+ compact(&e1,(_double *)&s1,sizeof(SINGLE));
+ return(s1);
+}
--- /dev/null
+/*
+ ADD TWO FLOATS - DOUBLE
+*/
+
+#include "FP_types.h"
+
+_double
+adf8(s2,s1)
+_double s1,s2;
+{
+ EXTEND e1,e2;
+ short swap;
+ extend(&s1,&e1,sizeof(_double));
+ extend(&s2,&e2,sizeof(_double));
+#ifdef PRT_EXT
+ prt_ext("ADF8: e1",&e1);
+ prt_ext("ADF8: e2",&e2);
+#endif
+ /* adjust mantissas to equal powers */
+ if (e1.sign ^ e2.sign) { /* signs are different */
+ /* determine which is largest number */
+ swap = (e2.exp > e1.exp) ? 1 : (e2.exp < e1.exp) ? 0 :
+ (e2.m1 > e1.m1 ) ? 1 : (e2.m1 < e1.m1 ) ? 0 :
+ (e2.m2 > e1.m2 ) ? 1 : 0;
+ /* adjust mantissas to equal powers */
+ sft_ext(&e1,&e2);
+ /* subtract the extended formats */
+ if (swap) { /* &e2 is the largest number */
+#ifdef PRT_EXT
+ fprintf(stderr,"ADF8: swaps and subtracts extended\n");
+#endif
+ sub_ext(&e2,&e1);
+ e1 = e2;
+ }
+ else {
+ sub_ext(&e1,&e2);
+ }
+ }
+ else {
+ /* adjust mantissas to equal powers */
+ sft_ext(&e1,&e2);
+ add_ext(&e1,&e2);
+ }
+#ifdef PRT_EXT
+ prt_ext("ADF8: e1 result",&e1);
+#endif
+ compact(&e1,&s1,sizeof(_double));
+ return(s1);
+}
--- /dev/null
+/*
+ CONVERT DOUBLE TO FLOAT
+
+ 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"
+
+cff4(src)
+_double src; /* the source itself - THIS TIME it's DOUBLE */
+{
+ EXTEND buf;
+
+ extend(&src,&buf,8); /* no matter what */
+#ifdef PRT_EXT
+ prt_ext("CFF4() entry:",&buf);
+ fprintf(stderr,"ds(%d),ss(%d),src(%08X%08X)\n",8,4,src.__double[0],
+ src.__double[1]);
+#endif PRT_EXT
+ compact(&buf,(_double *) &(src.__double[1]),4);
+#ifdef PRT_EXT
+ fprintf(stderr,"CFF4() exit : %08X\n",src.__double[1]);
+#endif PRT_EXT
+}
--- /dev/null
+/*
+ CONVERT FLOAT TO DOUBLE
+
+ 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"
+
+cff8(src)
+_float src; /* the space on the stack is for a double - see cg/table */
+{
+ EXTEND buf;
+
+ extend((_double *) &src,&buf,4); /* no matter what */
+#ifdef PRT_EXT
+ prt_ext("CFF8() entry:",&buf);
+ fprintf(stderr,"ds(%d),ss(%d),src(%08X)\n",4,8,src);
+#endif
+ compact(&buf,&src,8);
+#ifdef DEBUG
+ fprintf(stderr,"CFF8() exit : %08X",src.__double[0]);
+ fprintf(stderr,"%08X\n",src.__double[1]);
+#endif DEBUG
+}
--- /dev/null
+/*
+ CONVERT FLOAT TO UNSIGNED
+
+ 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
+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 newint, max_exp;
+
+ extend(&src,&buf,ss); /* get extended format */
+#ifdef PRT_EXT
+ prt_ext("CFI() entry:",&buf);
+#endif PRT_EXT
+ buf.exp--; /* additional bias correction */
+ if (buf.exp < 1) { /* no conversion needed */
+ src.__double[ss == 8] = 0L;
+ return(0L);
+ }
+ max_exp = (ds << 3) - 1; /* signed numbers */
+ /* have more limited max_exp */
+ if (buf.exp > max_exp) {
+#ifdef PRT_EXT
+ prt_ext("CFI() INT OVERFLOW", &buf);
+#endif PRT_EXT
+ trap(EIOVFL); /* integer overflow */
+ buf.exp %= max_exp; /* truncate */
+ }
+ new = buf.m1 >> (32-buf.exp);
+ if (buf.sign)
+ new = -new;
+done:
+ src.__double[ss == 8] = new;
+ return(new);
+}
--- /dev/null
+/*
+ CONVERT FLOAT TO UNSIGNED
+
+ 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,&buf,ss); /* get extended format */
+#ifdef PRT_EXT
+ prt_ext("CFU() entry:",&buf);
+#endif PRT_EXT
+ buf.exp--; /* additional bias correction */
+ if (buf.exp < 1) { /* no conversion needed */
+ src.__double[ss == 8] = 0L;
+ return(0L);
+ }
+ max_exp = (ds << 3);
+ if (buf.exp > max_exp) {
+#ifdef PRT_EXT
+ prt_ext("CFU() INT OVERFLOW",&buf);
+#endif PRT_EXT
+ trap(EIOVFL); /* integer overflow */
+ buf.exp %= max_exp;
+ }
+ new = buf.m1 >> (32-buf.exp);
+done:
+ src.__double[ss == 8] = new;
+ return(new);
+}
--- /dev/null
+/*
+ CONVERT INTEGER TO FLOAT
+
+ 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"
+
+_float
+cif4(ss,src)
+int ss; /* source size */
+long src; /* largest possible integer to convert */
+{
+ EXTEND buf;
+ short *ipt;
+ long i_src;
+ _float *result;
+
+ zrf_ext(&buf);
+ if (ss == sizeof(long)) {
+ buf.exp = 33;
+ i_src = src;
+ result = (_float *) &src;
+ }
+ else {
+ ipt = (short *) &src;
+ i_src = (long) *ipt;
+ buf.exp = 17;
+ result = (_float *) &ss;
+ }
+#ifdef PRT_STDERR
+ fprintf(stderr,"CIF4(ds(%d),ss(%d),src(%D))\n\n",4,ss,i_src);
+#endif
+ if (i_src == 0) {
+ *result = (_float) 0L;
+ return(0L);
+ }
+ /* 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 */
+#ifdef PRT_STDERR
+ fprintf(stderr,"CIF() buf.exp after nrm_ext() == %d\n\n",buf.exp);
+#endif
+ compact(&buf,(_double *) result,4); /* put on stack */
+ return(*result);
+}
--- /dev/null
+/*
+ CONVERT INTEGER TO FLOAT
+
+ 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"
+
+_double
+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 *) &ss; /* always */
+ zrf_ext(&buf);
+ if (ss == sizeof(long)) {
+ buf.exp = 33;
+ i_src = src;
+ }
+ else {
+ ipt = (short *) &src;
+ i_src = (long) *ipt;
+ buf.exp = 17;
+ }
+#ifdef PRT_STDERR
+ fprintf(stderr,"CIF8(ds(%d),ss(%d),src(%D))\n\n",8,ss,i_src);
+#endif
+ if (i_src == 0) {
+ zrf8(result);
+ return(*result);
+ }
+ /* 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);
+#ifdef PRT_STDERR
+ fprintf(stderr,"CIF() buf.exp after nrm_ext() == %d\n\n",buf.exp);
+#endif
+ compact(&buf,result,8);
+ return(*result);
+}
--- /dev/null
+/*
+ COMPARE DOUBLES
+*/
+
+#include "FP_types.h"
+#include "get_put.h"
+
+short
+cmf4(f1,f2)
+_float 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)
+ return ((sign1 > 0) ? -1 : 1);
+
+ return (sign1 * ((l1 < l2) ? 1 : -1));
+}
--- /dev/null
+/*
+ COMPARE DOUBLES
+*/
+
+#include "FP_types.h"
+#include "get_put.h"
+
+short
+cmf8(d1,d2)
+_double d1,d2;
+{
+#define SIGN(x) (((x) < 0) ? -1 : 1)
+ /*
+ * return ((d1 < d2) ? 1 : (d1 > d2) ? -1 : 0))
+ */
+ long l1,l2;
+ unsigned short *s1,*s2;
+ int sign1,sign2;
+
+ l1 = get4((char *)&d1);
+ l2 = get4((char *)&d2);
+ sign1 = SIGN(l1);
+ sign2 = SIGN(l2);
+ if (sign1 != sign2)
+ return ((sign1 > 0) ? -1 : 1);
+ if (l1 != l2) { /* we can decide here */
+ s1 = (unsigned short *) &l1;
+ s2 = (unsigned short *) &l2;
+ /* set both signs positive */
+ *s1 &= ~0x8000;
+ *s2 &= ~0x8000;
+ /* we already know they aren't equal so */
+ return (sign1 * ((l1 < l2) ? 1 : -1));
+ }
+ else { /* decide in 2nd half */
+ l1 = get4(((char *)&d1 + 4));
+ l2 = get4(((char *)&d2 + 4));
+ if (l1 == l2)
+ return(0);
+ else {
+ s1 = (unsigned short *) &l1;
+ s2 = (unsigned short *) &l2;
+ if (*s1 == *s2) {
+ s1++;
+ s2++;
+ }
+ return (sign1 * ((*s1 < *s2) ? 1 : -1));
+ }
+ }
+}
--- /dev/null
+/*
+#define PRT_EXIT
+#define PRT_TRAP
+#define PRT_ENTRY
+ 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"
+
+compact(f,to,size)
+EXTEND *f;
+_double *to;
+int size;
+{
+ DOUBLE *DBL;
+ int error = 0;
+ SINGLE *SGL;
+ int exact;
+
+#ifdef PRT_ENTRY
+ prt_ext("enter compact:",f);
+#endif PRT_ENTRY
+ if (size == sizeof(_double))
+/********************************************************/
+/*
+ COMPACT EXTENDED INTO DOUBLE
+*/
+/********************************************************/
+ {
+ if ((f->m1|(f->m2 & DBL_ZERO)) == 0L) {
+ zrf8(to);
+ goto leave;
+ }
+ f->exp += DBL_BIAS; /* restore proper bias */
+ if (f->exp > DBL_MAX) {
+dbl_over: trap(EFOVFL);
+#ifdef PRT_TRAP
+ prt_ext("FCOMPACT DBL OVERFLOW",f);
+#endif PRT_TRAP
+ f->exp = DBL_MAX;
+ f->m1 = f->m2 = 0L;
+ if (error++)
+ return;
+ }
+ else if (f->exp < DBL_MIN) {
+#ifdef PRT_TRAP
+ prt_ext("FCOMPACT DBL UNDERFLOW",f);
+#endif PRT_TRAP
+ trap(EFUNFL);
+ f->exp = DBL_MIN;
+ f->m1 = f->m2 = 0L;
+ if (error++)
+ return;
+ }
+
+ /* local CAST conversion */
+ DBL = (DOUBLE *) to;
+ /* check if answer is exact */
+ /* (the last 11 bits are zero (0) */
+
+ exact = ((f->m2 & DBL_EXACT) == 0) ? 1 : 0;
+
+ /* because of special format shift only 10 bits */
+ /* bit shift mantissa 10 bits */
+
+ /* first align within words, then do store operation */
+
+ DBL->_s.p1.fract = f->m1 >> DBL_RUNPACK; /* plus 22 == 32 */
+ DBL->_s.p2 = f->m2 >> DBL_RUNPACK; /* plus 22 == 32 */
+ DBL->_s.p2 |= (f->m1 << DBL_LUNPACK); /* plus 10 == 32 */
+
+ /* if not exact then round to nearest */
+
+ if (!exact) {
+ /* INEXACT(); */
+ if (f->m2 & DBL_ROUNDUP) {
+ DBL->_s.p2++; /* rounding up */
+ if (DBL->_s.p2 == 0L) { /* carry out */
+#ifdef PRT_RNDMSG
+ write(2,"rounding up lsb\n",16);
+#endif PRT_RNDMSG
+ DBL->_s.p1.fract++;
+ if (DBL->_s.p1.fract & DBL_CARRYOUT) { /* carry out */
+#ifdef PRT_RNDMSG
+ write(2,"shift due to rounding\n",22);
+#endif PRT_RNDMSG
+ if (DBL->_s.p1.fract & 01)
+ DBL->_s.p2 = CARRYBIT;
+ DBL->_s.p1.fract >>= 1;
+ f->exp++;
+ }
+ }
+ }
+ }
+ /* check for overflow */
+ if (f->exp >= DBL_MAX)
+ goto dbl_over;
+
+ /* STORE EXPONENT: */
+
+ /* 1) clear leading bits (B4-B15) */
+ DBL->_s.p1.fract &= DBL_MASK;
+
+ /* 2) shift and store exponent */
+ f->exp <<= DBL_EXPSHIFT;
+ DBL->_s.p1.fract |= ((long) f->exp << EXP_STORE);
+ }
+ else
+/********************************************************/
+/*
+ COMPACT EXTENDED INTO FLOAT
+*/
+/********************************************************/
+ {
+ /* local CAST conversion */
+ SGL = (SINGLE *) to;
+ if ((f->m1 & SGL_ZERO) == 0L) {
+ SGL->fract = 0L;
+ goto leave;
+ }
+ f->exp += SGL_BIAS; /* restore bias */
+ if (f->exp > SGL_MAX) {
+sgl_over: trap(EFOVFL);
+#ifdef PRT_TRAP
+ prt_ext("FCOMPACT FLOAT OVERFLOW",f);
+#endif PRT_TRAP
+ f->exp = SGL_MAX;
+ f->m1 = f->m2 = 0L;
+ if (error++)
+ return;
+ }
+ else if (f->exp < SGL_MIN) {
+#ifdef PRT_TRAP
+ prt_ext("FCOMPACT FLOAT UNDERFLOW",f);
+#endif PRT_TRAP
+ trap(EFUNFL);
+ f->exp = SGL_MIN;
+ f->m1 = f->m2 = 0L;
+ if (error++)
+ return;
+ }
+ /* check if the answer is exact */
+ /* the last 40 bits are zero */
+ /* check first last bits of mantissa 1 */
+ exact = ((f->m1 & SGL_EXACT) == 0) ? 1 : 0;
+
+ /* check last 32 bits in mantissa 2 */
+ if (exact) /* first part masks to zero */
+ exact = (f->m2 == 0L) ? 1 : 0;
+
+ /* shift mantissa and store */
+ SGL->fract = (f->m1 >> SGL_RUNPACK);
+
+ /* check for rounding to nearest */
+ if (!exact) {
+ /* INEXACT(); */
+ if (f->m1 & SGL_ROUNDUP) {
+ SGL->fract++;
+#ifdef PRT_RNDMSG
+ write(2,"rounding up lsb\n",16);
+#endif PRT_RNDMSG
+ /* check normal */
+ if (SGL->fract & SGL_CARRYOUT) {
+ SGL->fract >>= 1;
+ f->exp++;
+ }
+ }
+ }
+ if (f->exp >= SGL_MAX)
+ goto sgl_over;
+
+ /* STORE EXPONENT */
+ /* 1) clear leading bit of fraction */
+ SGL->fract &= SGL_MASK; /* B23-B31 are 0 */
+
+ /* 2) shift and store exponent */
+ f->exp <<= SGL_EXPSHIFT;
+ SGL->fract |= ((long) f->exp << EXP_STORE);
+ }
+
+/********************************************************/
+/*
+ STORE SIGN BIT
+*/
+/********************************************************/
+ if (f->sign) {
+ SGL = (SINGLE *) to; /* make sure */
+ SGL->fract |= CARRYBIT;
+ }
+/********************************************************/
+/*
+ STORE MANTISSA
+/*
+/********************************************************/
+
+ if (size == sizeof(_double)) {
+ put4(DBL->_s.p1.fract, (char *) &DBL->_s.p1.fract);
+ put4(DBL->_s.p2, (char *) &DBL->_s.p2);
+ }
+ else
+ put4(SGL->fract, (char *) &SGL->fract);
+
+leave:
+#ifdef PRT_EXIT
+ prt_ext("exit compact:",f);
+ prt_dbl((DOUBLE *) to,size); getchar();
+#endif PRT_EXIT
+ ; /* end of print statement or null statement */
+}
--- /dev/null
+/*
+ CONVERT INTEGER TO FLOAT
+
+ 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"
+
+cuf4(ss,src)
+int ss; /* source size */
+long src; /* largest possible integer to convert */
+{
+ EXTEND buf;
+ short *ipt;
+ _float *result;
+ long i_src;
+
+ zrf_ext(&buf);
+ if (ss == sizeof(long)) {
+ buf.exp = 33;
+ i_src = src;
+ result = (_float *) &src;
+ }
+ else {
+ ipt = (short *) &src;
+ i_src = (long) *ipt;
+ buf.exp = 17;
+ result = (_float *) &ss;
+ }
+#ifdef PRT_STDERR
+ fprintf(stderr,"CUF4(ds(%d),ss(%d),src(%D))\n\n",4,ss,i_src);
+#endif
+ if (i_src == 0) {
+ src = 0L;
+ }
+ /* 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);
+#ifdef PRT_STDERR
+ fprintf(stderr,"CUF() buf.exp after nrm_ext() == %d\n\n",buf.exp);
+#endif
+ compact(&buf,(_double *) result,4);
+}
--- /dev/null
+/*
+ CONVERT INTEGER TO FLOAT
+
+ 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"
+
+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 = 33;
+ i_src = src;
+ }
+ else {
+ ipt = (short *) &src;
+ i_src = (long) *ipt;
+ buf.exp = 17;
+ }
+#ifdef PRT_STDERR
+ fprintf(stderr,"CUF8(ds(%d),ss(%d),src(%D))\n\n",8,ss,i_src);
+#endif
+ if (i_src == 0) {
+ zrf8(&src);
+ 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);
+#ifdef PRT_STDERR
+ fprintf(stderr,"CUF() buf.exp after nrm_ext() == %d\n\n",buf.exp);
+#endif
+ compact(&buf,(_double *) &ss,8);
+}
--- /dev/null
+/*
+#define PRT_EXT
+#define PRT_ALL
+ 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.
+ It is based on the partial products method
+ and makes no use possible machine instructions
+ to divide (hardware dividers). It is intended
+ that it be rewritten to do so, but expedieancy
+ requires that something be written NOW - and
+ this is it.
+*/
+/********************************************************/
+
+div_ext(e1,e2)
+EXTEND *e1,*e2;
+{
+ short count;
+ short error = 0;
+ unsigned long result[2];
+ register unsigned long *lp;
+
+#ifdef PRT_EXT
+ fprintf("stderr:start div_ext:\n");
+ prt_ext("dividend:",e1);
+ prt_ext("divisor :",e2);
+#endif
+ if ((e1->m1 | e1->m2) == 0) { /* 0 / anything == 0 */
+ e1->exp = 0; /* make sure */
+ return;
+ }
+ /*
+ * numbers are right shifted one bit to make sure
+ * that m1 is quaranteed to be larger if its
+ * maximum bit is set
+ */
+ b64_sft(&e1->m1,1); /* 64 bit shift right */
+ b64_sft(&e2->m1,1); /* 64 bit shift right */
+ e1->exp++;
+ e2->exp++;
+ /* check for underflow, divide by zero, etc */
+ e1->sign ^= e2->sign;
+ e1->exp -= e2->exp;
+ e1->exp += 2; /* bias correction */
+ if (e1->exp < EXT_MIN) {
+ error++;
+#ifdef PRT_EXT
+ prt_ext("DIV_EXT UNDERFLOW",e1);
+#endif PRT_EXT
+ trap(EFUNFL); /* underflow */
+ e1->exp = EXT_MIN;
+ e1->m1 = e1->m2 = 0L;
+ }
+ if ((e2->m1 | e2->m2) == 0) {
+ error++;
+#ifdef PRT_EXT
+ prt_ext("DIV_EXT DIV 0.0",e2);
+#endif PRT_EXT
+ trap(EFDIVZ);
+ e1->m1 = e1->m2 = 0L;
+ e1->exp = EXT_MAX;
+ }
+ if (error)
+ return;
+
+ /* do division of mantissas */
+ /* uses partial product method */
+ /* init control variables */
+
+ count = 64;
+ lp = result; /* result[0] == high word */
+ /* result[0] == low word */
+ *lp++ = 0L; /* high word */
+ *lp-- = 0L; /* low word */
+
+ /* partial product division loop */
+
+ while (count--) {
+ /* first left shift result 1 bit */
+ /* this is ALWAYS done */
+
+ b64_sft(result,-1);
+
+ /* compare dividend and divisor */
+ /* if dividend >= divisor add a bit */
+ /* and subtract divisior from dividend */
+#ifdef PRT_ALL
+ prt_ext("dividend:",e1);
+ prt_ext("divisor :",e2);
+#endif
+
+ if ( (e1->m1 < e2->m1) ||
+ ((e1->m1 == e2->m1) && (e1->m2 < e2->m2) ))
+ ; /* null statement */
+ /* i.e., don't add or subtract */
+ else {
+ result[1]++; /* ADD */
+ if (e2->m2 > e1->m2)
+ e1->m1 -= 1; /* carry in */
+ e1->m1 -= e2->m1; /* do SUBTRACTION */
+ e1->m2 -= e2->m2; /* SUBTRACTION */
+#ifdef PRT_ALL
+ prt_ext("result :",e1);
+#endif
+ }
+#ifdef PRT_ALL
+ fprintf(stderr,"div_ext %d %08X%08X\n\n",64-count,
+ result[0],result[1]);
+ fflush(stderr);
+#endif
+
+ /* 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_sft(&e1->m1,-1);
+ continue;
+ }
+ else
+ break; /* leave loop */
+ } /* end of divide by subtraction loop */
+
+ /* DISPLAY RESULTS FOR DEBUGGING */
+#ifdef PRT_ALL
+ prt_ext("dividend:",e1);
+ prt_ext("divisor :",e2);
+ fprintf(stderr,"div_ext %d %08X%08X\n",64-count,
+ result[0],result[1]);
+#endif
+
+ if (count > 0) {
+ lp = result;
+ 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[1] */
+ if (*lp) {
+ result[0] |= (*lp >> 32-count);
+ *lp <<= count;
+ }
+ }
+ /*
+ if (error)
+ INEXACT();
+ */
+ e1->m1 = result[0];
+ e1->m2 = result[1];
+#ifdef PRT_EXT
+ prt_ext("result :",e1);
+#endif
+ nrm_ext(e1);
+#ifdef PRT_EXT
+ prt_ext("after nrm:",e1);
+ /*sleep(4);*/
+#endif
+}
--- /dev/null
+/*
+ DIVIDE TWO FLOATS - SINGLE Precision
+*/
+
+#include "FP_types.h"
+
+dvf4(s2,s1)
+_float s1,s2;
+{
+ EXTEND e1,e2;
+
+ extend((_double *)&s1,&e1,sizeof(_float));
+ extend((_double *)&s2,&e2,sizeof(_float));
+
+ /* do a divide */
+ div_ext(&e1,&e2);
+ compact(&e1,(_double *)&s1,sizeof(_float));
+}
--- /dev/null
+/*
+ DIVIDE TWO FLOATS - DOUBLE Precision
+*/
+
+#include "FP_types.h"
+
+dvf8(s2,s1)
+_double s1,s2;
+{
+ EXTEND e1,e2;
+
+ extend(&s1,&e1,sizeof(_double));
+ extend(&s2,&e2,sizeof(_double));
+
+ /* do a divide */
+ div_ext(&e1,&e2);
+ compact(&e1,&s1,sizeof(_double));
+}
--- /dev/null
+/*
+#define PRT_EXIT
+#define PRT_ENTRY
+#define PRT_DBL
+ 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"
+/********************************************************/
+
+extend(from,to,size)
+_double *from;
+EXTEND *to;
+int size;
+{
+ DOUBLE *f;
+ register char *cpt1,*cpt2;
+ unsigned long tmp;
+ int leadbit = 0;
+
+#ifdef PRT_ENTRY
+ write(2,"entry extend: ",14);
+#ifdef PRT_DBL
+ prt_dbl(from,size);
+#else
+ write(2,"\n",1);
+#endif PRT_DBL
+#endif PRT_ENTRY
+ f = (DOUBLE *) from; /* local cast conversion */
+ if (f->_s.p1.fract == 0L) {
+ if (size == sizeof(SINGLE)) {
+zero: zrf_ext(to);
+ goto ready;
+ }
+ else if (f->_s.p2 == 0L)
+ goto zero;
+ }
+/* there is a number to convert so lets get started */
+/* first extract the exponent; its always in the first two bytes */
+
+ cpt1 = (char *) from;
+ to->exp = uget2(cpt1);
+ 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 */
+
+ to->m1 = get4(cpt1);
+
+ if (size == sizeof(DOUBLE)) {
+ to->m1 <<= DBL_M1LEFT; /* shift */
+ to->exp -= DBL_BIAS; /* remove bias */
+ cpt1 += 4;
+ tmp = get4(cpt1);
+ to->m1 |= (tmp>>DBL_RPACK); /* plus 10 == 32 */
+ to->m2 = (tmp<<DBL_LPACK); /* plus 22 == 32 */
+ }
+ else { /* size == sizeof(SINGLE) */
+ to->m1 <<= SGL_M1LEFT; /* shift */
+ 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 */
+ready:
+#ifdef PRT_EXIT
+ prt_ext("exit extend:",to)
+#endif PRT_EXIT
+ ; /* end of print statement or null statement */
+}
--- /dev/null
+/*
+ SEPERATE INTO EXPONENT AND FRACTION
+*/
+
+#include "FP_types.h"
+
+struct fef4_returns {
+ short e;
+ _float f;
+};
+
+fef4(s1)
+_float s1;
+{
+ struct fef4_returns *r = (struct fef4_returns *) &s1;
+ EXTEND buf;
+
+ extend((_double *) &s1,&buf,sizeof(_float));
+ r->e = buf.exp-1;
+ buf.exp = 1;
+ compact(&buf,(_double *) &r->f,sizeof(_float));
+}
--- /dev/null
+/*
+ SEPERATE DOUBLE INTO EXPONENT AND FRACTION
+*/
+
+#include "FP_types.h"
+
+struct fef8_returns {
+ short e;
+ _double f;
+};
+
+fef8(s1)
+_double s1;
+{
+ EXTEND buf;
+ struct fef8_returns *r = (struct fef8_returns *) &s1;
+
+#ifdef DEBUG
+ printf("FEF8(): ");
+#endif DEBUG
+ extend(&s1,&buf,sizeof(_double));
+ r->e = buf.exp - 1;
+ buf.exp = 1;
+ compact(&buf,&r->f,sizeof(_double));
+#ifdef DEBUG
+ printf("exponent = %3d fraction = 0x%08X%08X: ",
+ r->f.__double[0],r->f.__double[1]);
+ printf("FEF8()\n");
+#endif DEBUG
+}
--- /dev/null
+/*
+ MULTIPLY AND DISMEMBER PARTS
+*/
+
+#include "FP_types.h"
+#include "FP_shift.h"
+
+_float mlf4();
+_float sbf4();
+
+fif4(x,y)
+_float x,y;
+{
+ EXTEND e;
+
+ y = mlf4(x,y);
+ extend((_double *)&y,&e,sizeof(SINGLE));
+ e.exp--; /* additional bias correction */
+ if (e.exp < 1) {
+ x = 0;
+ return;
+ }
+ if (e.exp > 31 - SGL_M1LEFT) {
+ x = y;
+ y = 0;
+ return;
+ }
+ b64_sft(&e.m1, 64 - e.exp);
+ b64_sft(&e.m1, e.exp - 64); /* "loose" low order bits */
+ e.exp++;
+ compact(&e,(_double *) &x, sizeof(SINGLE));
+ y = sbf4(x, y);
+}
--- /dev/null
+/*
+ MULTIPLY AND DISMEMBER PARTS
+*/
+
+#include "FP_types.h"
+#include "FP_shift.h"
+
+_double mlf8();
+_double sbf8();
+
+fif8(x,y)
+_double x,y;
+{
+ EXTEND e;
+
+ y = mlf8(x,y);
+ extend((_double *)&y,&e,sizeof(DOUBLE));
+ e.exp--; /* additional bias correction */
+ if (e.exp < 1) {
+ x.__double[0] = 0;
+ x.__double[1] = 0;
+ return;
+ }
+ if (e.exp > 63 - DBL_M1LEFT) {
+ x.__double[0] = y.__double[0];
+ x.__double[1] = y.__double[1];
+ y.__double[0] = 0;
+ y.__double[1] = 0;
+ return;
+ }
+ b64_sft(&e.m1, 64 - e.exp);
+ b64_sft(&e.m1, e.exp - 64); /* "loose" low order bits */
+ e.exp++;
+ compact(&e, &x, sizeof(DOUBLE));
+ y = sbf8(x, y);
+}
--- /dev/null
+#
+
+ mes 2,EM_WSIZE,EM_PSIZE
+
+#define TRAP 0
+
+; _fptrp is called with one parameter:
+; - trap number (TRAP)
+
+ exp $_fptrp
+ pro $_fptrp,0
+ lol TRAP
+ trp
+ ret 0
+ end ?
--- /dev/null
+#include <byte_order.h>
+
+#if CHAR_UNSIGNED
+#define Xchar(ch) (ch)
+#else
+#define Xchar(ch) ((ch) & 0377)
+#endif
+
+#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
--- /dev/null
+/*
+ * Multiply Single Precesion Float
+ */
+
+#include "FP_types.h"
+
+_float
+mlf4(s2,s1)
+_float s1,s2;
+{
+ EXTEND e1,e2;
+
+ extend((_double *)&s1,&e1,sizeof(_float));
+ extend((_double *)&s2,&e2,sizeof(_float));
+ /* do a multiply */
+ mul_ext(&e1,&e2);
+ compact(&e1,(_double *)&s1,sizeof(_float));
+ return(s1);
+}
--- /dev/null
+/*
+ * Multiply Single Precesion Float
+ */
+
+#include "FP_types.h"
+
+_double
+mlf8(s2,s1)
+_double s1,s2;
+{
+ EXTEND e1,e2;
+
+ extend(&s1,&e1,sizeof(_double));
+ extend(&s2,&e2,sizeof(_double));
+ /* do a multiply */
+ mul_ext(&e1,&e2);
+ compact(&e1,&s1,sizeof(_double));
+ return(s1);
+}
--- /dev/null
+/*
+ ROUTINE TO MULTIPLY TWO EXTENDED FORMAT NUMBERS
+*/
+
+# include "adder.h"
+# include "FP_bias.h"
+# include "FP_trap.h"
+# include "FP_types.h"
+
+mul_ext(e1,e2)
+EXTEND *e1,*e2;
+{
+ register short k,i,j; /* loop control */
+ long unsigned *reg[7];
+ long unsigned tmp[4];
+ short unsigned mp[4]; /* multiplier */
+ short unsigned mc[4]; /* multipcand */
+ B64 low64,tmp64; /* 64 bit storage */
+
+#ifdef PRT_EXT
+ prt_ext("before MUL_EXT() e1:",e1);
+ prt_ext("before MUL_EXT() e2:",e2);
+#endif
+ /* first save the sign (XOR) */
+
+ e1->sign ^= e2->sign;
+
+ /********************************************************/
+ /* INCREASE EXPONENT BY ONE (1) */
+ /* */
+ /* the nature of the multiplication algorithm used */
+ /* results in an exponent that is small by an additive */
+ /* factor of one (1); */
+ /* if the maximum bit is set it will not be subtracted */
+ /* during normalization -> this is correct and can be */
+ /* expected often with normalized numbers */
+ /* HOWEVER, it is also possible that unnormalized */
+ /* numbers are used. Rather than shifting here */
+ /* always(!) (unless L bit is set) I chose to */
+ /* increase the exponent by one - a simple (FAST) */
+ /* process - and to decrease it later during */
+ /* normalization. */
+ /* */
+ /********************************************************/
+ /* The effects of bias (as used here) */
+ /* and the multiplication algorithm used cancel */
+ /* so these statements are commented out */
+ /* August 1985 - if changing the Leading Bit (or NORMBIT) */
+ /* this problem with the multiplication algorithm no longer */
+ /* exists - bias must be subtracted now */
+ /* */
+ /* e1->exp++; */
+ /********************************************************/
+
+ /* next add the exponents */
+
+ e1->exp += e2->exp;
+ e1->exp -= 1; /* correction for bias */
+
+ /* check for overflow */
+ if (e1->exp >= EXT_MAX) {
+#ifdef PRT_EXT
+ prt_ext("EXT_MUL OVERFLOW",e1);
+#endif
+ trap(EFOVFL);
+ /* if caught */
+ /* return signed infinity */
+ e1->exp = EXT_MAX;
+infinity: e1->m1 = e1->m2 =0L;
+#ifdef PRT_EXT
+ prt_ext("after MUL_EXT() e1:",e1);
+#endif
+ return;
+ }
+ /* check for underflow */
+ if (e1->exp < EXT_MIN) {
+#ifdef PRT_EXT
+ prt_ext("EXT_MUL UNDERFLOW",e1);
+#endif
+ trap(EFUNFL);
+ e1->exp = EXT_MIN;
+ goto infinity;
+ }
+
+ /* 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;
+# ifdef DEBUG
+ for(i=0;i<4;i++)
+ printf("%04x",mp[i]);
+ putchar('\r');
+ putchar('\n');
+ for(i=0;i<4;i++)
+ printf("%04x",mc[i]);
+ putchar('\r');
+ putchar('\n');
+# endif
+ /*
+ * assign pointers
+ */
+ reg[0] = &e1->m1; /* the answer goes here */
+ reg[1] = &tmp[1];
+ reg[2] = &e1->m2; /* and here */
+ reg[3] = &tmp[2];
+ reg[4] = &low64.h_32;
+ reg[5] = &tmp[3];
+ reg[6] = &low64.l_32;
+
+ /*
+ * zero registers
+ */
+ for(i=7;i--;)
+ *reg[i] = 0;
+
+ /*
+ * fill registers with their components
+ */
+ for(i=4;i--;) if (mp[i])
+ for(j=4;j--;) if (mc[j]) {
+ k = i+j;
+ tmp[0] = (long)mp[i] * (long)mc[j];
+# ifdef PRT_EXT2
+ printf("%04x * %04x == %08X ",mp[i],mc[j],tmp[0]);
+ printf("index == %d ",k);
+ printf("register before add == %08X\n",*reg[k]);
+ fflush(stdout);
+# endif
+#ifdef PRT_ADD
+ printf("REGISTERS-----\n");
+ printf("%08X %08X %08X %08X\n0000%04x %04x%04x %04x%04x %04x0000\n",
+ *reg[0],*reg[2],*reg[4],*reg[6],
+ (short)(*reg[1]>>16),(short)(*reg[1]),(short)(*reg[3]>>16),
+ (short)(*reg[3]),(short)(*reg[5]>>16),(short)(*reg[5]));
+# endif
+ if (b32_add(reg[k],tmp)) {
+ for(tmp[0] = 0x10000L;k>0;)
+ if (b32_add(reg[--k],tmp) == 0)
+ break;
+#ifdef PRT_ADD
+ printf("CARRY---------\n");
+ printf("%08X %08X %08X %08X\n0000%04x %04x%04x %04x%04x %04x0000\n",
+ *reg[0],*reg[2],*reg[4],*reg[6],
+ (short)(*reg[1]>>16),(short)(*reg[1]),(short)(*reg[3]>>16),
+ (short)(*reg[3]),(short)(*reg[5]>>16),(short)(*reg[5]));
+#endif
+ }
+ }
+
+ /*
+ * combine the registers to a total
+ */
+#ifdef PRT_ADD
+ printf("%08X %08X %08X %08X\n0000%04x %04x%04x %04x%04x %04x0000\n",
+ *reg[0],*reg[2],*reg[4],*reg[6],
+ (short)(*reg[1]>>16),(short)(*reg[1]),(short)(*reg[3]>>16),
+ (short)(*reg[3]),(short)(*reg[5]>>16),(short)(*reg[5]));
+# endif
+ tmp64.h_32 = (*reg[1]>>16);
+ tmp64.l_32 = (*reg[1]<<16) + (*reg[3]>>16);
+# ifdef PRT_ALL
+ printf("%08X%08X tmp64\n",tmp64.h_32,tmp64.l_32);
+ fflush(stdout);
+ printf("%08X%08X e1->m1\n",e1->m1,e1->m2);
+ fflush(stdout);
+# endif
+ b64_add((B64 *)&e1->m1,&tmp64);
+# ifdef PRT_ALL
+ printf("b64_add:\n");
+ printf("%08X%08X e1->m1\n",e1->m1,e1->m2);
+ fflush(stdout);
+# endif
+ tmp64.l_32 = *reg[5]<<16;
+ tmp64.h_32 = (*reg[5]>>16) + (*reg[3]<<16);
+ if (b64_add(&low64,&tmp64))
+ if (++e1->m2 == 0)
+ e1->m1++;
+
+# ifdef PRT_ADD
+ printf("%08X %08X %08X %08X\n",e1->m1,e1->m2,low64.h_32,low64.l_32);
+ fflush(stdout);
+#endif
+#ifdef PRT_EXT
+ prt_ext("after MUL_EXT() e1:",e1);
+#endif PRT_EXT
+ nrm_ext(e1);
+#ifdef PRT_EXT
+ prt_ext("after NRM_EXT() e1:",e1);
+ sleep(4);
+#endif PRT_EXT
+}
--- /dev/null
+/*
+ NEGATE A FLOATING POINT
+*/
+/********************************************************/
+/*
+ Assumes exponent is located in bytes 0 & 1
+*/
+/********************************************************/
+
+#include "FP_types.h"
+
+ngf4(f)
+_float f;
+{
+ char unsigned *p;
+
+ p = (char unsigned *) &f;
+ *p ^= 0x80;
+}
+
--- /dev/null
+/*
+ NEGATE A FLOATING POINT
+*/
+/********************************************************/
+/*
+ Assumes exponent is located in bytes 0 & 1
+*/
+/********************************************************/
+
+#include "FP_types.h"
+
+ngf8(f)
+_double f;
+{
+ unsigned char *p;
+
+ p = (unsigned char *) &f;
+ *p ^= 0x80;
+}
+
--- /dev/null
+/********************************************************/
+/*
+ NORMALIZE an EXTENDED FORMAT NUMBER
+*/
+/********************************************************/
+
+#include "FP_shift.h"
+#include "FP_types.h"
+
+nrm_ext(e1)
+EXTEND *e1;
+{
+ register unsigned long *mant_1;
+ register unsigned long *mant_2;
+
+ /* local CAST conversion */
+#ifdef PRT_EXT
+ prt_ext("before NRM_EXT() e1:",e1);
+#endif PRT_EXT
+ mant_1 = (unsigned long *) &e1->m1;
+ /*
+ THIS RESULTS IN A BAD CODE !!!!
+ ANOTHER BUG IN EM CODE MAYBE????
+ mant_2 = (unsigned long *) &e1->m2;
+ */
+ /* statement that works */
+ mant_2 = mant_1 + 1;
+
+ /* 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 ((*mant_1 | *mant_2) == 0L)
+ return;
+
+ /* if top word is zero mov low word */
+ /* to top word, adjust exponent value */
+ if (*mant_1 == 0L) {
+ *mant_1++ = e1->m2;
+ *mant_1-- = 0L;
+ e1->exp -= 32;
+ }
+#ifdef OLD
+ /* check that e1->m1 is not too large */
+ if (*mant_1 & CARRYBIT) { /* carry occured */
+ e1->exp++; /* increase exponent */
+ *mant_2 >>= 1; /* right shift mantissa */
+ if ((short) *mant_1 & 01)
+ *mant_2 |= CARRYBIT;
+ *mant_1 >>= 1;
+ }
+#endif
+
+
+ while ((*mant_1 & NORMBIT) == 0) {
+ e1->exp--;
+ *mant_1 <<= 1;
+ if ((*mant_2 & CARRYBIT) == 0)
+ ; /* empty statement */
+ else {
+ *mant_1 += 1;
+ }
+ *mant_2 <<= 1;
+ }
+#ifdef PRT_EXT
+ prt_ext("after NRM_EXT() e1:",e1);
+#endif PRT_EXT
+}
--- /dev/null
+# include "FP_types.h"
+
+prt_dbl(dbl,size)
+DOUBLE *dbl;
+int size;
+{
+#ifdef PRT_DBL
+ unsigned long *l;
+
+ fprintf(stderr,"PRT_DBL SIZE = %d ",size);
+ fprintf(stderr,"_s.p1.fract = 0x%08X ",dbl->_s.p1.fract);
+ if (size == 8)
+ fprintf(stderr,"_s.p2 = 0x%08X",dbl->_s.p2);
+ l = (unsigned long *) dbl;
+#ifdef PRT_LONG
+ fprintf(stderr,"\nl[0] = 0x%08X ",*l++);
+ if (size == 8)
+ fprintf(stderr,"l[1] = 0x%08X",*l);
+#endif PRT_LONG
+ putc('\r',stderr);
+ putc('\n',stderr);
+ fflush(stderr);
+#endif
+}
+
--- /dev/null
+/********************************************************/
+/*
+ PRINT EXTENDED FORMAT AND MESSAGE
+ DEBUG ROUTINE
+*/
+/********************************************************/
+
+#include "FP_types.h"
+
+prt_ext(m,e)
+char *m;
+EXTEND *e;
+{
+#ifdef PRT_EXT
+ fprintf(stderr,"%s ",m);
+ fprintf(stderr,"%c",(e->sign) ? '-' : '+');
+ fprintf(stderr,"m1:0x%08X m2:0x%08X ^ %03d 0x%x\n",
+ e->m1,e->m2,e->exp,e->exp);
+ fprintf(stderr,"hit any key\n\r");
+ fflush(stderr);
+ getchar();
+#endif
+}
--- /dev/null
+/*
+ SUBTRACT TWO FLOATS - SINGLE Precision
+*/
+
+#include "FP_types.h"
+
+extern _float adf4();
+
+_float
+sbf4(s2,s1)
+_float s1,s2;
+{
+ /* changing the sign directly */
+ /* is faster than the code: */
+ /* s2 = -s2 */
+ char unsigned *p;
+
+ p = (char unsigned *) &s2;
+ *p ^= 0x80; /* change sign of s2 */
+ s1 = adf4(s2,s1);
+ return(s1); /* add and return result */
+}
+
--- /dev/null
+/*
+ SUBTRACT TWO FLOATS - DOUBLE Precision
+*/
+
+#include "FP_types.h"
+
+extern _double adf8();
+
+_double
+sbf8(s2,s1)
+_double s1,s2;
+{
+ /* changing the sign directly */
+ /* is faster than the code line */
+ /* s2 = -s2; */
+ char unsigned *p; /* sufficient to access sign bit */
+
+#ifdef PRT_EXT
+ fprintf(stderr,"SBF8 ():\n");
+#endif
+ p = (char unsigned *) &s2;
+ *p ^= 0x80; /* change sign of s2 */
+ s1 = adf8(s2,s1); /* add and return result */
+ return(s1);
+}
+
--- /dev/null
+/*
+#define PRT_EXT
+ SHIFT TWO EXTENDED NUMBERS INTO PROPER
+ ALIGNMENT FOR ADDITION (exponents are equal)
+*/
+
+#include "FP_types.h"
+
+sft_ext(e1,e2)
+EXTEND *e1,*e2;
+{
+ register EXTEND *s;
+ register short diff;
+ long tmp;
+
+#ifdef PRT_EXT
+ prt_ext("enter sft_ext e1:",e1);
+ prt_ext("enter sft_ext e2:",e2);
+#endif PRT_EXT
+ 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;
+
+ if (diff > 63) { /* no relative value */
+ s->m1 = 0L;
+ s->m2 = 0L;
+ return;
+ }
+
+ if (diff > 32) {
+ diff -= 32;
+ s->m2 = s->m1;
+ s->m1 = 0L;
+ }
+ if (diff) {
+ if (s->m1) {
+ tmp = s->m1;
+ tmp <<= (32-diff);
+ s->m1 >>= diff;
+ }
+ else
+ tmp = 0L;
+ s->m2 >>= diff;
+ s->m2 |= tmp;
+ }
+#ifdef PRT_EXT
+ prt_ext("exit sft_ext e1:",e1);
+ prt_ext("exit sft_ext e2:",e2);
+#endif PRT_EXT
+}
--- /dev/null
+# include "adder.h"
+
+b64_sft(e1,n)
+B64 *e1;
+short n;
+{
+ if (n>0) do { /* RIGHT shift n bits */
+ e1->l_32 >>= 1; /* shift 64 bits */
+ if (e1->h_32 & 1)
+ e1->l_32 |= 0x80000000L;
+ e1->h_32 >>= 1;
+ } while (--n);
+ else /* LEFT shift n bits */
+ while (n++) {
+ e1->h_32 <<= 1; /* shift 64 bits */
+ if (e1->l_32 & 0x80000000L)
+ e1->h_32 |= 1;
+ e1->l_32 <<= 1;
+ }
+}
--- /dev/null
+/*
+#define PRT_EXT
+ SUBTRACT EXTENDED FORMAT
+*/
+ /* assumes that e1 >= e2 on entry */
+ /* no test is made to check this */
+ /* so make sure yourself */
+
+#include "FP_types.h"
+sub_ext(e1,e2)
+EXTEND *e1,*e2;
+{
+#ifdef PRT_EXT
+ prt_ext("before SUB_EXT() e1:",e1);
+ prt_ext("before SUB_EXT() e2:",e2);
+#endif PRT_EXT
+ if (e2->m2 > e1->m2)
+ e1->m1 -= 1; /* carry in */
+ e1->m1 -= e2->m1;
+ e1->m2 -= e2->m2;
+#ifdef PRT_EXT
+ prt_ext("after SUB_EXT() e1:",e1);
+#endif PRT_EXT
+ nrm_ext(e1);
+#ifdef PRT_EXT
+ prt_ext("after NRM_EXT() e1:",e1);
+#endif PRT_EXT
+}
--- /dev/null
+/*
+ return a zero float
+*/
+
+#include "FP_types.h"
+
+zrf4(l)
+long *l;
+{
+ *l = 0L;
+}
--- /dev/null
+/*
+ return a zero double
+*/
+
+#include "FP_types.h"
+
+zrf8(z)
+_double *z;
+{
+
+ z->__double[0] = 0L;
+ z->__double[1] = 0L;
+}
--- /dev/null
+/*
+ ZERO and return EXTEND FORMAT FLOAT
+*/
+
+#include "FP_types.h"
+
+zrf_ext(e)
+EXTEND *e;
+{
+ register short *ipt;
+ register short i;
+ register short zero = 0;
+
+ /* local CAST conversion */
+ ipt = (short *) e;
+
+ i = sizeof(EXTEND)/sizeof(short);
+ while (i--)
+ *ipt++ = zero;
+}
public / ack.git / commitdiff