From 915950ede281243b2c7a5400980ef16681cf3ab4 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Wed, 3 Apr 2002 10:26:12 +0000 Subject: run dos2unix on these files --- test/crypt/Makefile | 146 +- test/math/Makefile | 182 +- test/math/econst.c | 192 +- test/math/eexp.c | 154 +- test/math/ehead.h | 84 +- test/math/elog.c | 184 +- test/math/eparanoi.c | 7100 +++++++++++++++++++++---------------------- test/math/epow.c | 430 +-- test/math/etanh.c | 104 +- test/math/etodec.c | 362 +-- test/math/ieee.c | 8238 +++++++++++++++++++++++++------------------------- test/math/ieetst.c | 1700 +++++------ test/math/ieetst.doc | 264 +- test/math/mconf.h | 216 +- test/math/mtherr.c | 192 +- 15 files changed, 9774 insertions(+), 9774 deletions(-) (limited to 'test') diff --git a/test/crypt/Makefile b/test/crypt/Makefile index fe743aa45..a965466f3 100644 --- a/test/crypt/Makefile +++ b/test/crypt/Makefile @@ -1,73 +1,73 @@ -# Makefile for uClibc -# -# Copyright (C) 2002 Erik Andersen -# -# This program is free software; you can redistribute it and/or modify it under -# the terms of the GNU Library General Public License as published by the Free -# Software Foundation; either version 2 of the License, or (at your option) any -# later version. -# -# This program is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more -# details. -# -# You should have received a copy of the GNU Library General Public License -# along with this program; if not, write to the Free Software Foundation, Inc., -# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - - - -TESTDIR=../ -include $(TESTDIR)/Rules.mak - -TARGETS=diff md5c-test -EXTRA_LIBS=-lcrypt - -all: $(TARGETS) - -crypt: crypt.c Makefile $(TESTDIR)/Config $(TESTDIR)/Rules.mak $(CC) - -@ echo "-------" - -@ echo " " - -@ echo "Compiling vs uClibc: " - -@ echo " " - $(CC) $(CFLAGS) -c $< -o $@.o - $(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS) - $(STRIPTOOL) -x -R .note -R .comment $@ - -./$@ < crypt.input > $@.out 2>&1 - -@ echo " " - -crypt_glibc: crypt.c Makefile $(TESTDIR)/Config $(TESTDIR)/Rules.mak $(CC) - -@ echo "-------" - -@ echo " " - -@ echo "Compiling vs uClibc: " - -@ echo " " - $(HOST_CC) $(GLIBC_CFLAGS) -c $< -o $@.o - $(HOST_CC) $(GLIBC_LDFLAGS) $@.o -o $@ $(EXTRA_LIBS) - $(STRIPTOOL) -x -R .note -R .comment $@ - -./$@ < crypt.input > $@.out 2>&1 - -@ echo " " - -diff: crypt_glibc crypt - -@ echo "-------" - -@ echo " " - -@ echo "Diffing output: " - -@ echo " " - -diff -u crypt_glibc.out crypt.out - -@ echo " " - -md5c-test: md5c-test.c Makefile $(TESTDIR)/Config $(TESTDIR)/Rules.mak $(CC) - -@ echo "-------" - -@ echo " " - -@ echo "Compiling vs uClibc: " - -@ echo " " - $(CC) $(CFLAGS) -c $< -o $@.o - $(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS) - $(STRIPTOOL) -x -R .note -R .comment $@ - -./$@ - -@ echo " " - -clean: - rm -f *.[oa] *~ core crypt_glibc crypt crypt_glibc.out crypt.out md5c-test - - +# Makefile for uClibc +# +# Copyright (C) 2002 Erik Andersen +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU Library General Public License as published by the Free +# Software Foundation; either version 2 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more +# details. +# +# You should have received a copy of the GNU Library General Public License +# along with this program; if not, write to the Free Software Foundation, Inc., +# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + +TESTDIR=../ +include $(TESTDIR)/Rules.mak + +TARGETS=crypt md5c-test +EXTRA_LIBS=-lcrypt + +all: $(TARGETS) + +crypt: crypt.c Makefile $(TESTDIR)/Config $(TESTDIR)/Rules.mak $(CC) + -@ echo "-------" + -@ echo " " + -@ echo "Compiling vs uClibc: " + -@ echo " " + $(CC) $(CFLAGS) -c $< -o $@.o + $(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS) + $(STRIPTOOL) -x -R .note -R .comment $@ + -./$@ < crypt.input #> $@.out 2>&1 + -@ echo " " + +crypt_glibc: crypt.c Makefile $(TESTDIR)/Config $(TESTDIR)/Rules.mak $(CC) + -@ echo "-------" + -@ echo " " + -@ echo "Compiling vs uClibc: " + -@ echo " " + $(HOST_CC) $(GLIBC_CFLAGS) -c $< -o $@.o + $(HOST_CC) $(GLIBC_LDFLAGS) $@.o -o $@ $(EXTRA_LIBS) + $(STRIPTOOL) -x -R .note -R .comment $@ + -./$@ < crypt.input > $@.out 2>&1 + -@ echo " " + +diff: crypt_glibc crypt + -@ echo "-------" + -@ echo " " + -@ echo "Diffing output: " + -@ echo " " + -diff -u crypt_glibc.out crypt.out + -@ echo " " + +md5c-test: md5c-test.c Makefile $(TESTDIR)/Config $(TESTDIR)/Rules.mak $(CC) + -@ echo "-------" + -@ echo " " + -@ echo "Compiling vs uClibc: " + -@ echo " " + $(CC) $(CFLAGS) -c $< -o $@.o + $(CC) $(LDFLAGS) $@.o -o $@ $(EXTRA_LIBS) + $(STRIPTOOL) -x -R .note -R .comment $@ + -./$@ + -@ echo " " + +clean: + rm -f *.[oa] *~ core crypt_glibc crypt crypt_glibc.out crypt.out md5c-test + + diff --git a/test/math/Makefile b/test/math/Makefile index be7f261de..4d0bc5ee6 100644 --- a/test/math/Makefile +++ b/test/math/Makefile @@ -1,91 +1,91 @@ -# Makefile for uClibc -# -# Copyright (C) 2000,2001 Erik Andersen -# -# This program is free software; you can redistribute it and/or modify it under -# the terms of the GNU Library General Public License as published by the Free -# Software Foundation; either version 2 of the License, or (at your option) any -# later version. -# -# This program is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more -# details. -# -# You should have received a copy of the GNU Library General Public License -# along with this program; if not, write to the Free Software Foundation, Inc., -# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - - - -# Unix makefile for ieetst, eparanoi. -# Set LARGEMEM 1 in qcalc.h for 32-bit memory addresses. -# Define computer type and/or endianness in mconf.h. -# -# Configure eparanoi.c for desired arithmetic test; -# also define appropriate version of setprec.o, or use a stub that -# does no FPU setup. To test native arithmetic, eparanoi uses -# the system libraries only; compile simply by `cc eparanoi.c -lm'. -# - -TESTDIR=../ -include $(TESTDIR)/Rules.mak - - -#CC = gcc -#CFLAGS= -O -INCS= mconf.h ehead.h -OBJS = ieee.o econst.o eexp.o elog.o epow.o etanh.o etodec.o mtherr.o #setprec.o -TARGETS=ieetst eparanoi - -all: $(TARGETS) - -ieetst: ieetst.o $(OBJS) drand.o $(INCS) - $(CC) -o ieetst ieetst.o $(OBJS) drand.o -lc -lm - -eparanoi: eparanoi.o $(OBJS) $(INCS) - $(CC) -o eparanoi eparanoi.o $(OBJS) -lc -lm - -#setprec.o: setprec.387 -# as -o setprec.o setprec.387 - -#setprec.o: setprec.688 -# as -o setprec.o setprec.688 - -ieee.o: ieee.c $(INCS) - $(CC) $(CFLAGS) -c ieee.c - -econst.o: econst.c $(INCS) - $(CC) $(CFLAGS) -c econst.c - -elog.o: elog.c $(INCS) - $(CC) $(CFLAGS) -c elog.c - -eexp.o: eexp.c $(INCS) - $(CC) $(CFLAGS) -c eexp.c - -etanh.o: etanh.c $(INCS) - $(CC) $(CFLAGS) -c etanh.c - -epow.o: epow.c $(INCS) - $(CC) $(CFLAGS) -c epow.c - -mtherr.o: mtherr.c $(INCS) - $(CC) $(CFLAGS) -c mtherr.c - -ieetst.o: ieetst.c $(INCS) - $(CC) $(CFLAGS) -c ieetst.c - -drand.o: drand.c $(INCS) - $(CC) $(CFLAGS) -c drand.c - -etodec.o: etodec.c $(INCS) - $(CC) $(CFLAGS) -c etodec.c - -eparanoi.o: eparanoi.c $(INCS) - $(CC) $(CFLAGS) -c eparanoi.c - -clean: - rm -f *.[oa] *~ core $(TARGETS) - - +# Makefile for uClibc +# +# Copyright (C) 2000,2001 Erik Andersen +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU Library General Public License as published by the Free +# Software Foundation; either version 2 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more +# details. +# +# You should have received a copy of the GNU Library General Public License +# along with this program; if not, write to the Free Software Foundation, Inc., +# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + + +# Unix makefile for ieetst, eparanoi. +# Set LARGEMEM 1 in qcalc.h for 32-bit memory addresses. +# Define computer type and/or endianness in mconf.h. +# +# Configure eparanoi.c for desired arithmetic test; +# also define appropriate version of setprec.o, or use a stub that +# does no FPU setup. To test native arithmetic, eparanoi uses +# the system libraries only; compile simply by `cc eparanoi.c -lm'. +# + +TESTDIR=../ +include $(TESTDIR)/Rules.mak + + +#CC = gcc +#CFLAGS= -O +INCS= mconf.h ehead.h +OBJS = ieee.o econst.o eexp.o elog.o epow.o etanh.o etodec.o mtherr.o #setprec.o +TARGETS=ieetst eparanoi + +all: $(TARGETS) + +ieetst: ieetst.o $(OBJS) drand.o $(INCS) + $(CC) -o ieetst ieetst.o $(OBJS) drand.o -lc -lm + +eparanoi: eparanoi.o $(OBJS) $(INCS) + $(CC) -o eparanoi eparanoi.o $(OBJS) -lc -lm + +#setprec.o: setprec.387 +# as -o setprec.o setprec.387 + +#setprec.o: setprec.688 +# as -o setprec.o setprec.688 + +ieee.o: ieee.c $(INCS) + $(CC) $(CFLAGS) -c ieee.c + +econst.o: econst.c $(INCS) + $(CC) $(CFLAGS) -c econst.c + +elog.o: elog.c $(INCS) + $(CC) $(CFLAGS) -c elog.c + +eexp.o: eexp.c $(INCS) + $(CC) $(CFLAGS) -c eexp.c + +etanh.o: etanh.c $(INCS) + $(CC) $(CFLAGS) -c etanh.c + +epow.o: epow.c $(INCS) + $(CC) $(CFLAGS) -c epow.c + +mtherr.o: mtherr.c $(INCS) + $(CC) $(CFLAGS) -c mtherr.c + +ieetst.o: ieetst.c $(INCS) + $(CC) $(CFLAGS) -c ieetst.c + +drand.o: drand.c $(INCS) + $(CC) $(CFLAGS) -c drand.c + +etodec.o: etodec.c $(INCS) + $(CC) $(CFLAGS) -c etodec.c + +eparanoi.o: eparanoi.c $(INCS) + $(CC) $(CFLAGS) -c eparanoi.c + +clean: + rm -f *.[oa] *~ core $(TARGETS) + + diff --git a/test/math/econst.c b/test/math/econst.c index 21523e568..cfddbe3e2 100644 --- a/test/math/econst.c +++ b/test/math/econst.c @@ -1,96 +1,96 @@ -/* econst.c */ -/* e type constants used by high precision check routines */ - -#include "ehead.h" - - -#if NE == 10 -/* 0.0 */ -unsigned short ezero[NE] = - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,}; - -/* 5.0E-1 */ -unsigned short ehalf[NE] = - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,}; - -/* 1.0E0 */ -unsigned short eone[NE] = - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,}; - -/* 2.0E0 */ -unsigned short etwo[NE] = - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,}; - -/* 3.2E1 */ -unsigned short e32[NE] = - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,}; - -/* 6.93147180559945309417232121458176568075500134360255E-1 */ -unsigned short elog2[NE] = - {0x40f3, 0xf6af, 0x03f2, 0xb398, - 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,}; - -/* 1.41421356237309504880168872420969807856967187537695E0 */ -unsigned short esqrt2[NE] = - {0x1d6f, 0xbe9f, 0x754a, 0x89b3, - 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,}; - -/* 3.14159265358979323846264338327950288419716939937511E0 */ -unsigned short epi[NE] = - {0x2902, 0x1cd1, 0x80dc, 0x628b, - 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,}; - -/* 5.7721566490153286060651209008240243104215933593992E-1 */ -unsigned short eeul[NE] = { -0xd1be,0xc7a4,0076660,0063743,0111704,0x3ffe,}; - -#else - -/* 0.0 */ -unsigned short ezero[NE] = { -0, 0000000,0000000,0000000,0000000,0000000,}; -/* 5.0E-1 */ -unsigned short ehalf[NE] = { -0, 0000000,0000000,0000000,0100000,0x3ffe,}; -/* 1.0E0 */ -unsigned short eone[NE] = { -0, 0000000,0000000,0000000,0100000,0x3fff,}; -/* 2.0E0 */ -unsigned short etwo[NE] = { -0, 0000000,0000000,0000000,0100000,0040000,}; -/* 3.2E1 */ -unsigned short e32[NE] = { -0, 0000000,0000000,0000000,0100000,0040004,}; -/* 6.93147180559945309417232121458176568075500134360255E-1 */ -unsigned short elog2[NE] = { -0xc9e4,0x79ab,0150717,0013767,0130562,0x3ffe,}; -/* 1.41421356237309504880168872420969807856967187537695E0 */ -unsigned short esqrt2[NE] = { -0x597e,0x6484,0174736,0171463,0132404,0x3fff,}; -/* 2/sqrt(PI) = - * 1.12837916709551257389615890312154517168810125865800E0 */ -unsigned short eoneopi[NE] = { -0x71d5,0x688d,0012333,0135202,0110156,0x3fff,}; -/* 3.14159265358979323846264338327950288419716939937511E0 */ -unsigned short epi[NE] = { -0xc4c6,0xc234,0020550,0155242,0144417,0040000,}; -/* 5.7721566490153286060651209008240243104215933593992E-1 */ -unsigned short eeul[NE] = { -0xd1be,0xc7a4,0076660,0063743,0111704,0x3ffe,}; -#endif -extern unsigned short ezero[]; -extern unsigned short ehalf[]; -extern unsigned short eone[]; -extern unsigned short etwo[]; -extern unsigned short e32[]; -extern unsigned short elog2[]; -extern unsigned short esqrt2[]; -extern unsigned short eoneopi[]; -extern unsigned short epi[]; -extern unsigned short eeul[]; - +/* econst.c */ +/* e type constants used by high precision check routines */ + +#include "ehead.h" + + +#if NE == 10 +/* 0.0 */ +unsigned short ezero[NE] = + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,}; + +/* 5.0E-1 */ +unsigned short ehalf[NE] = + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,}; + +/* 1.0E0 */ +unsigned short eone[NE] = + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,}; + +/* 2.0E0 */ +unsigned short etwo[NE] = + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,}; + +/* 3.2E1 */ +unsigned short e32[NE] = + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,}; + +/* 6.93147180559945309417232121458176568075500134360255E-1 */ +unsigned short elog2[NE] = + {0x40f3, 0xf6af, 0x03f2, 0xb398, + 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,}; + +/* 1.41421356237309504880168872420969807856967187537695E0 */ +unsigned short esqrt2[NE] = + {0x1d6f, 0xbe9f, 0x754a, 0x89b3, + 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,}; + +/* 3.14159265358979323846264338327950288419716939937511E0 */ +unsigned short epi[NE] = + {0x2902, 0x1cd1, 0x80dc, 0x628b, + 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,}; + +/* 5.7721566490153286060651209008240243104215933593992E-1 */ +unsigned short eeul[NE] = { +0xd1be,0xc7a4,0076660,0063743,0111704,0x3ffe,}; + +#else + +/* 0.0 */ +unsigned short ezero[NE] = { +0, 0000000,0000000,0000000,0000000,0000000,}; +/* 5.0E-1 */ +unsigned short ehalf[NE] = { +0, 0000000,0000000,0000000,0100000,0x3ffe,}; +/* 1.0E0 */ +unsigned short eone[NE] = { +0, 0000000,0000000,0000000,0100000,0x3fff,}; +/* 2.0E0 */ +unsigned short etwo[NE] = { +0, 0000000,0000000,0000000,0100000,0040000,}; +/* 3.2E1 */ +unsigned short e32[NE] = { +0, 0000000,0000000,0000000,0100000,0040004,}; +/* 6.93147180559945309417232121458176568075500134360255E-1 */ +unsigned short elog2[NE] = { +0xc9e4,0x79ab,0150717,0013767,0130562,0x3ffe,}; +/* 1.41421356237309504880168872420969807856967187537695E0 */ +unsigned short esqrt2[NE] = { +0x597e,0x6484,0174736,0171463,0132404,0x3fff,}; +/* 2/sqrt(PI) = + * 1.12837916709551257389615890312154517168810125865800E0 */ +unsigned short eoneopi[NE] = { +0x71d5,0x688d,0012333,0135202,0110156,0x3fff,}; +/* 3.14159265358979323846264338327950288419716939937511E0 */ +unsigned short epi[NE] = { +0xc4c6,0xc234,0020550,0155242,0144417,0040000,}; +/* 5.7721566490153286060651209008240243104215933593992E-1 */ +unsigned short eeul[NE] = { +0xd1be,0xc7a4,0076660,0063743,0111704,0x3ffe,}; +#endif +extern unsigned short ezero[]; +extern unsigned short ehalf[]; +extern unsigned short eone[]; +extern unsigned short etwo[]; +extern unsigned short e32[]; +extern unsigned short elog2[]; +extern unsigned short esqrt2[]; +extern unsigned short eoneopi[]; +extern unsigned short epi[]; +extern unsigned short eeul[]; + diff --git a/test/math/eexp.c b/test/math/eexp.c index dd2a64e1c..14ea9899d 100644 --- a/test/math/eexp.c +++ b/test/math/eexp.c @@ -1,77 +1,77 @@ -/* xexp.c */ -/* exponential function check routine */ -/* by Stephen L. Moshier. */ - - -#include "ehead.h" - -/* -extern int powinited; -extern short maxposint[], maxnegint[]; -*/ - -void eexp( x, y ) -unsigned short *x, *y; -{ -unsigned short num[NE], den[NE], x2[NE]; -long i; -unsigned short sign, expchk; - -/* range reduction theory: x = i + f, 0<=f<1; - * e**x = e**i * e**f - * e**i = 2**(i/log 2). - * Let i/log2 = i1 + f1, 0<=f1<1. - * Then e**i = 2**i1 * 2**f1, so - * e**x = 2**i1 * e**(log 2 * f1) * e**f. - */ -/* -if( powinited == 0 ) - initpow(); -*/ -if( ecmp(x, ezero) == 0 ) - { - emov( eone, y ); - return; - } -emov(x, x2); -expchk = x2[NE-1]; -sign = expchk & 0x8000; -x2[NE-1] &= 0x7fff; - -/* Test for excessively large argument */ -expchk &= 0x7fff; -if( expchk > (EXONE + 15) ) - { - eclear( y ); - if( sign == 0 ) - einfin( y ); - return; - } - -eifrac( x2, &i, num ); /* x = i + f */ - -if( i != 0 ) - { - ltoe( &i, den ); /* floating point i */ - ediv( elog2, den, den ); /* i/log 2 */ - eifrac( den, &i, den ); /* i/log 2 = i1 + f1 */ - emul( elog2, den, den ); /* log 2 * f1 */ - eadd( den, num, x2 ); /* log 2 * f1 + f */ - } - -/*x2[NE-1] -= 1;*/ -eldexp( x2, -1L, x2 ); /* divide by 2 */ -etanh( x2, x2 ); /* tanh( x/2 ) */ -eadd( x2, eone, num ); /* 1 + tanh */ -eneg( x2 ); -eadd( x2, eone, den ); /* 1 - tanh */ -ediv( den, num, y ); /* (1 + tanh)/(1 - tanh) */ - -/*y[NE-1] += i;*/ -if( sign ) - { - ediv( y, eone, y ); - i = -i; - } -eldexp( y, i, y ); /* multiply by 2**i */ -} +/* xexp.c */ +/* exponential function check routine */ +/* by Stephen L. Moshier. */ + + +#include "ehead.h" + +/* +extern int powinited; +extern short maxposint[], maxnegint[]; +*/ + +void eexp( x, y ) +unsigned short *x, *y; +{ +unsigned short num[NE], den[NE], x2[NE]; +long i; +unsigned short sign, expchk; + +/* range reduction theory: x = i + f, 0<=f<1; + * e**x = e**i * e**f + * e**i = 2**(i/log 2). + * Let i/log2 = i1 + f1, 0<=f1<1. + * Then e**i = 2**i1 * 2**f1, so + * e**x = 2**i1 * e**(log 2 * f1) * e**f. + */ +/* +if( powinited == 0 ) + initpow(); +*/ +if( ecmp(x, ezero) == 0 ) + { + emov( eone, y ); + return; + } +emov(x, x2); +expchk = x2[NE-1]; +sign = expchk & 0x8000; +x2[NE-1] &= 0x7fff; + +/* Test for excessively large argument */ +expchk &= 0x7fff; +if( expchk > (EXONE + 15) ) + { + eclear( y ); + if( sign == 0 ) + einfin( y ); + return; + } + +eifrac( x2, &i, num ); /* x = i + f */ + +if( i != 0 ) + { + ltoe( &i, den ); /* floating point i */ + ediv( elog2, den, den ); /* i/log 2 */ + eifrac( den, &i, den ); /* i/log 2 = i1 + f1 */ + emul( elog2, den, den ); /* log 2 * f1 */ + eadd( den, num, x2 ); /* log 2 * f1 + f */ + } + +/*x2[NE-1] -= 1;*/ +eldexp( x2, -1L, x2 ); /* divide by 2 */ +etanh( x2, x2 ); /* tanh( x/2 ) */ +eadd( x2, eone, num ); /* 1 + tanh */ +eneg( x2 ); +eadd( x2, eone, den ); /* 1 - tanh */ +ediv( den, num, y ); /* (1 + tanh)/(1 - tanh) */ + +/*y[NE-1] += i;*/ +if( sign ) + { + ediv( y, eone, y ); + i = -i; + } +eldexp( y, i, y ); /* multiply by 2**i */ +} diff --git a/test/math/ehead.h b/test/math/ehead.h index 746b60df7..24c95ce05 100644 --- a/test/math/ehead.h +++ b/test/math/ehead.h @@ -1,42 +1,42 @@ - -/* Include file for extended precision arithmetic programs. - */ - -/* Number of 16 bit words in external x type format */ -#define NE 6 - -/* Number of 16 bit words in internal format */ -#define NI (NE+3) - -/* Array offset to exponent */ -#define E 1 - -/* Array offset to high guard word */ -#define M 2 - -/* Number of bits of precision */ -#define NBITS ((NI-4)*16) - -/* Maximum number of decimal digits in ASCII conversion - * = NBITS*log10(2) - */ -#define NDEC (NBITS*8/27) - -/* The exponent of 1.0 */ -#define EXONE (0x3fff) - -void eadd(), esub(), emul(), ediv(); -int ecmp(), enormlz(), eshift(); -void eshup1(), eshup8(), eshup6(), eshdn1(), eshdn8(), eshdn6(); -void eabs(), eneg(), emov(), eclear(), einfin(), efloor(); -void eldexp(), efrexp(), eifrac(), ltoe(); -void esqrt(), elog(), eexp(), etanh(), epow(); -void asctoe(), asctoe24(), asctoe53(), asctoe64(); -void etoasc(), e24toasc(), e53toasc(), e64toasc(); -void etoe64(), etoe53(), etoe24(), e64toe(), e53toe(), e24toe(); -void mtherr(); -extern unsigned short ezero[], ehalf[], eone[], etwo[]; -extern unsigned short elog2[], esqrt2[]; - - -/* by Stephen L. Moshier. */ + +/* Include file for extended precision arithmetic programs. + */ + +/* Number of 16 bit words in external x type format */ +#define NE 6 + +/* Number of 16 bit words in internal format */ +#define NI (NE+3) + +/* Array offset to exponent */ +#define E 1 + +/* Array offset to high guard word */ +#define M 2 + +/* Number of bits of precision */ +#define NBITS ((NI-4)*16) + +/* Maximum number of decimal digits in ASCII conversion + * = NBITS*log10(2) + */ +#define NDEC (NBITS*8/27) + +/* The exponent of 1.0 */ +#define EXONE (0x3fff) + +void eadd(), esub(), emul(), ediv(); +int ecmp(), enormlz(), eshift(); +void eshup1(), eshup8(), eshup6(), eshdn1(), eshdn8(), eshdn6(); +void eabs(), eneg(), emov(), eclear(), einfin(), efloor(); +void eldexp(), efrexp(), eifrac(), ltoe(); +void esqrt(), elog(), eexp(), etanh(), epow(); +void asctoe(), asctoe24(), asctoe53(), asctoe64(); +void etoasc(), e24toasc(), e53toasc(), e64toasc(); +void etoe64(), etoe53(), etoe24(), e64toe(), e53toe(), e24toe(); +void mtherr(); +extern unsigned short ezero[], ehalf[], eone[], etwo[]; +extern unsigned short elog2[], esqrt2[]; + + +/* by Stephen L. Moshier. */ diff --git a/test/math/elog.c b/test/math/elog.c index 8afc1e5b4..bc517b197 100644 --- a/test/math/elog.c +++ b/test/math/elog.c @@ -1,92 +1,92 @@ -/* xlog.c */ -/* natural logarithm */ -/* by Stephen L. Moshier. */ - -#include "mconf.h" -#include "ehead.h" - - - -void elog( x, y ) -unsigned short *x, *y; -{ -unsigned short xx[NE], z[NE], a[NE], b[NE], t[NE], qj[NE]; -long ex; -int fex; - - -if( x[NE-1] & (unsigned short )0x8000 ) - { - eclear(y); - mtherr( "elog", DOMAIN ); - return; - } -if( ecmp( x, ezero ) == 0 ) - { - einfin( y ); - eneg(y); - mtherr( "elog", SING ); - return; - } -if( ecmp( x, eone ) == 0 ) - { - eclear( y ); - return; - } - -/* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */ -efrexp( x, &fex, xx ); -/* -emov(x, xx ); -ex = xx[NX-1] & 0x7fff; -ex -= 0x3ffe; -xx[NX-1] = 0x3ffe; -*/ - -/* Adjust range to 1/sqrt(2), sqrt(2) */ -esqrt2[NE-1] -= 1; -if( ecmp( xx, esqrt2 ) < 0 ) - { - fex -= 1; - emul( xx, etwo, xx ); - } -esqrt2[NE-1] += 1; - -esub( eone, xx, a ); -if( a[NE-1] == 0 ) - { - eclear( y ); - goto logdon; - } -eadd( eone, xx, b ); -ediv( b, a, y ); /* store (x-1)/(x+1) in y */ - -emul( y, y, z ); - -emov( eone, a ); -emov( eone, b ); -emov( eone, qj ); -do - { - eadd( etwo, qj, qj ); /* 2 * i + 1 */ - emul( z, a, a ); - ediv( qj, a, t ); - eadd( t, b, b ); - } -while( ((b[NE-1] & 0x7fff) - (t[NE-1] & 0x7fff)) < NBITS ); - - -emul( b, y, y ); -emul( y, etwo, y ); - -logdon: - -/* now add log of 2**ex */ -if( fex != 0 ) - { - ex = fex; - ltoe( &ex, b ); - emul( elog2, b, b ); - eadd( b, y, y ); - } -} +/* xlog.c */ +/* natural logarithm */ +/* by Stephen L. Moshier. */ + +#include "mconf.h" +#include "ehead.h" + + + +void elog( x, y ) +unsigned short *x, *y; +{ +unsigned short xx[NE], z[NE], a[NE], b[NE], t[NE], qj[NE]; +long ex; +int fex; + + +if( x[NE-1] & (unsigned short )0x8000 ) + { + eclear(y); + mtherr( "elog", DOMAIN ); + return; + } +if( ecmp( x, ezero ) == 0 ) + { + einfin( y ); + eneg(y); + mtherr( "elog", SING ); + return; + } +if( ecmp( x, eone ) == 0 ) + { + eclear( y ); + return; + } + +/* range reduction: log x = log( 2**ex * m ) = ex * log2 + log m */ +efrexp( x, &fex, xx ); +/* +emov(x, xx ); +ex = xx[NX-1] & 0x7fff; +ex -= 0x3ffe; +xx[NX-1] = 0x3ffe; +*/ + +/* Adjust range to 1/sqrt(2), sqrt(2) */ +esqrt2[NE-1] -= 1; +if( ecmp( xx, esqrt2 ) < 0 ) + { + fex -= 1; + emul( xx, etwo, xx ); + } +esqrt2[NE-1] += 1; + +esub( eone, xx, a ); +if( a[NE-1] == 0 ) + { + eclear( y ); + goto logdon; + } +eadd( eone, xx, b ); +ediv( b, a, y ); /* store (x-1)/(x+1) in y */ + +emul( y, y, z ); + +emov( eone, a ); +emov( eone, b ); +emov( eone, qj ); +do + { + eadd( etwo, qj, qj ); /* 2 * i + 1 */ + emul( z, a, a ); + ediv( qj, a, t ); + eadd( t, b, b ); + } +while( ((b[NE-1] & 0x7fff) - (t[NE-1] & 0x7fff)) < NBITS ); + + +emul( b, y, y ); +emul( y, etwo, y ); + +logdon: + +/* now add log of 2**ex */ +if( fex != 0 ) + { + ex = fex; + ltoe( &ex, b ); + emul( elog2, b, b ); + eadd( b, y, y ); + } +} diff --git a/test/math/eparanoi.c b/test/math/eparanoi.c index 0e479f9f5..84cab73f8 100644 --- a/test/math/eparanoi.c +++ b/test/math/eparanoi.c @@ -1,3550 +1,3550 @@ -/* paranoia.c arithmetic tester - * - * This is an implementation of the PARANOIA program. It substitutes - * subroutine calls for ALL floating point arithmetic operations. - * This permits you to substitute your own experimental versions of - * arithmetic routines. It also defeats compiler optimizations, - * so for native arithmetic you can be pretty sure you are testing - * the arithmetic and not the compiler. - * - * This version of PARANOIA omits the display of division by zero. - * It also omits the test for extra precise subexpressions, since - * they cannot occur in this context. Otherwise it includes all the - * tests of the 27 Jan 86 distribution, plus a few additional tests. - * Commentary has been reduced to a minimum in order to make the program - * smaller. - * - * The original PARANOIA program, written by W. Kahan, C version - * by Thos Sumner and David Gay, can be downloaded free from the - * Internet NETLIB. An MSDOS disk can be obtained for $15 from: - * Richard Karpinski - * 6521 Raymond Street - * Oakland, CA 94609 - * - * Steve Moshier, 28 Oct 88 - * last rev: 23 May 92 - */ - -#define DEBUG 0 - -/* To use the native arithmetic of the computer, define NATIVE - * to be 1. To use your own supplied arithmetic routines, NATIVE is 0. - */ -#define NATIVE 0 - -/* gcc real.c interface */ -#define L128DOUBLE 0 - -#include - - - - -/* Data structure of floating point number. If NATIVE was - * selected above, you can define LDOUBLE 1 to test 80-bit long double - * precision or define it 0 to test 64-bit double precision. -*/ -#define LDOUBLE 0 -#if NATIVE - -#define NE 1 -#if LDOUBLE -#define FSIZE long double -#define FLOAT(x) FSIZE x[NE] -static FSIZE eone[NE] = {1.0L}; /* The constant 1.0 */ -#define ZSQRT sqrtl -#define ZLOG logl -#define ZFLOOR floorl -#define ZPOW powl -long double sqrtl(), logl(), floorl(), powl(); -#define FSETUP einit -#else /* not LDOUBLE */ -#define FSIZE double -#define FLOAT(x) FSIZE x[NE] -static FSIZE eone[NE] = {1.0}; /* The constant 1.0 */ -#define ZSQRT sqrt -#define ZLOG log -#define ZFLOOR floor -#define ZPOW pow -double sqrt(), log(), floor(), pow(); -/* Coprocessor initialization, - * defeat underflow trap or what have you. - * This is required mainly on i386 and 68K processors. - */ -#define FSETUP dprec -#endif /* double, not LDOUBLE */ - -#else /* not NATIVE */ - -/* Setup for extended double type. - * Put NE = 10 for real.c operating with TFmode support (16-byte reals) - * Put NE = 6 for real.c operating with XFmode support (10- or 12-byte reals) - * The value of NE must agree with that in ehead.h, if ieee.c is used. - */ -#define NE 6 -#define FSIZE unsigned short -#define FLOAT(x) unsigned short x[NE] -extern unsigned short eone[]; -#define FSETUP einit - -/* default for FSETUP */ -/* -einit() -{} -*/ - -error(s) -char *s; -{ -printf( "error: %s\n", s ); -} - -#endif /* not NATIVE */ - - - -#if L128DOUBLE -/* real.c interface */ - -#undef FSETUP -#define FSETUP efsetup - -FLOAT(enone); - -#define ONE enone - -/* Use emov to convert from widest type to widest type, ... */ -/* -#define ENTOE emov -#define ETOEN emov -*/ - -/* ... else choose e24toe, e53toe, etc. */ -#define ENTOE e64toe -#define ETOEN etoe64 -#define NNBITS 64 - -#define NIBITS ((NE-1)*16) -extern int rndprc; - -efsetup() -{ -rndprc = NNBITS; -ETOEN(eone, enone); -} - -add(a,b,c) -FLOAT(a); -FLOAT(b); -FLOAT(c); -{ -unsigned short aa[10], bb[10], cc[10]; - -ENTOE(a,aa); -ENTOE(b,bb); -eadd(aa,bb,cc); -ETOEN(cc,c); -} - -sub(a,b,c) -FLOAT(a); -FLOAT(b); -FLOAT(c); -{ -unsigned short aa[10], bb[10], cc[10]; - -ENTOE(a,aa); -ENTOE(b,bb); -esub(aa,bb,cc); -ETOEN(cc,c); -} - -mul(a,b,c) -FLOAT(a); -FLOAT(b); -FLOAT(c); -{ -unsigned short aa[10], bb[10], cc[10]; - -ENTOE(a,aa); -ENTOE(b,bb); -emul(aa,bb,cc); -ETOEN(cc,c); -} - -div(a,b,c) -FLOAT(a); -FLOAT(b); -FLOAT(c); -{ -unsigned short aa[10], bb[10], cc[10]; - -ENTOE(a,aa); -ENTOE(b,bb); -ediv(aa,bb,cc); -ETOEN(cc,c); -} - -int cmp(a,b) -FLOAT(a); -FLOAT(b); -{ -unsigned short aa[10], bb[10]; -int c; -int ecmp(); - -ENTOE(a,aa); -ENTOE(b,bb); -c = ecmp(aa,bb); -return(c); -} - -mov(a,b) -FLOAT(a); -FLOAT(b); -{ -int i; - -for( i=0; ib, 0 if a==b, -1 if a (*b) ) - return( 1 ); -if( (*a) < (*b) ) - return( -1 ); -if( (*a) != (*b) ) - goto cmpf; -if( (*a) == (*b) ) - return( 0 ); -cmpf: -printf( "Compare fails\n" ); -return(0); -} - - -emov( a, b ) -FSIZE *a, *b; -{ -*b = *a; -} - -eneg( a ) -FSIZE *a; -{ -*a = -(*a); -} - -eclear(a) -FSIZE *a; -{ -*a = 0.0; -} - -eabs(x) -FSIZE *x; -{ -if( (*x) < 0.0 ) - *x = -(*x); -} - -efloor(x,y) -FSIZE *x, *y; -{ - -*y = (FSIZE )ZFLOOR( *x ); -} - -elog(x,y) -FSIZE *x, *y; -{ - -*y = (FSIZE )ZLOG( *x ); -} - -epow(x,y,z) -FSIZE *x, *y, *z; -{ - -*z = (FSIZE )ZPOW( *x, *y ); -} - -esqrt(x,y) -FSIZE *x, *y; -{ - -*y = (FSIZE )ZSQRT( *x ); -} - - -eifrac(x,i,f) -FSIZE *x; -long *i; -FSIZE *f; -{ -FSIZE y; - -y = (FSIZE )ZFLOOR( *x ); -if( y < 0.0 ) - { - *f = y - *x; - *i = -y; - } -else - { - *f = *x - y; - *i = y; - } -} - - -ltoe(i,x) -long *i; -FSIZE *x; -{ -*x = *i; -} - - -etoasc(a,str,n) -FSIZE *a; -char *str; -int n; -{ -double x; - -x = (double )(*a); -sprintf( str, " %.17e ", x ); -} - -/* default for FSETUP */ -einit() -{} - -#endif /* NATIVE */ - - - - -FLOAT(Radix); -FLOAT(BInvrse); -FLOAT(RadixD2); -FLOAT(BMinusU2); -/*Small floating point constants.*/ -FLOAT(Zero); -FLOAT(Half); -FLOAT(One); -FLOAT(Two); -FLOAT(Three); -FLOAT(Four); -FLOAT(Five); -FLOAT(Six); -FLOAT(Eight); -FLOAT(Nine); -FLOAT(Ten); -FLOAT(TwentySeven); -FLOAT(ThirtyTwo); -FLOAT(TwoForty); -FLOAT(MinusOne ); -FLOAT(OneAndHalf); - -/*Integer constants*/ -int NoTrials = 20; /*Number of tests for commutativity. */ -#define False 0 -#define True 1 - -/* Definitions for declared types - Guard == (Yes, No); - Rounding == (Chopped, Rounded, Other); - Message == packed array [1..40] of char; - Class == (Flaw, Defect, Serious, Failure); - */ -#define Yes 1 -#define No 0 -#define Chopped 2 -#define Rounded 1 -#define Other 0 -#define Flaw 3 -#define Defect 2 -#define Serious 1 -#define Failure 0 - -typedef int Guard, Rounding, Class; -typedef char Message; - -/* Declarations of Variables */ -FLOAT(AInvrse); -FLOAT(A1); -FLOAT(C); -FLOAT(CInvrse); -FLOAT(D); -FLOAT(FourD); -FLOAT(E0); -FLOAT(E1); -FLOAT(Exp2); -FLOAT(E3); -FLOAT(MinSqEr); -FLOAT(SqEr); -FLOAT(MaxSqEr); -FLOAT(E9); -FLOAT(Third); -FLOAT(F6); -FLOAT(F9); -FLOAT(H); -FLOAT(HInvrse); -FLOAT(StickyBit); -FLOAT(J); -FLOAT(MyZero); -FLOAT(Precision); -FLOAT(Q); -FLOAT(Q9); -FLOAT(R); -FLOAT(Random9); -FLOAT(T); -FLOAT(Underflow); -FLOAT(S); -FLOAT(OneUlp); -FLOAT(UfThold); -FLOAT(U1); -FLOAT(U2); -FLOAT(V); -FLOAT(V0); -FLOAT(V9); -FLOAT(W); -FLOAT(X); -FLOAT(X1); -FLOAT(X2); -FLOAT(X8); -FLOAT(Random1); -FLOAT(Y); -FLOAT(YY1); -FLOAT(Y2); -FLOAT(Random2); -FLOAT(Z); -FLOAT(PseudoZero); -FLOAT(Z1); -FLOAT(Z2); -FLOAT(Z9); -static FLOAT(t); -FLOAT(t2); -FLOAT(Sqarg); -int ErrCnt[4]; -int fpecount; -int Milestone; -int PageNo; -int I, M, N, N1, stkflg; -Guard GMult, GDiv, GAddSub; -Rounding RMult, RDiv, RAddSub, RSqrt; -int Break, Done, NotMonot, Monot, Anomaly, IEEE; -int SqRWrng, UfNGrad; -int k, k2; -int Indx; -char ch[8]; - -long lngint, lng2; /* intermediate for conversion between int and FLOAT */ - -/* Computed constants. */ -/*U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */ -/*U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */ - - -show( x ) -short x[]; -{ -int i; -char s[80]; - -/* Number of 16-bit groups to display */ -#if NATIVE -#if LDOUBLE -#define NPRT (sizeof( long double )/2) -#else -#define NPRT (sizeof( double )/2) -#endif -#else -#define NPRT NE -#endif - -TOASC( x, s, 70 ); -printf( "%s\n", s ); -for( i=0; i -#endif -#include -jmp_buf ovfl_buf; -/*typedef int (*Sig_type)();*/ -typedef void (*Sig_type)(); -Sig_type sigsave; - -/* Floating point exception receiver */ -void sigfpe() -{ -fpecount++; -printf( "\n* * * FLOATING-POINT ERROR * * *\n" ); -/* reinitialize the floating point unit */ -FSETUP(); -fflush(stdout); -if( sigsave ) - { -#ifndef NOSIGNAL - signal( SIGFPE, sigsave ); -#endif - sigsave = 0; - longjmp( ovfl_buf, 1 ); - } -abort(); -} - - -main() -{ - -/* Do coprocessor or other initializations */ -FSETUP(); - -printf( - "This version of paranoia omits test for extra precise subexpressions\n" ); -printf( "and includes a few additional tests.\n" ); - -clear(Zero); -printf( "0 = " ); -show( Zero ); -mov( ONE, One); -printf( "1 = " ); -show( One ); -add( One, One, Two ); -printf( "1+1 = " ); -show( Two ); -add( Two, One, Three ); -add( Three, One, Four ); -add( Four, One, Five ); -add( Five, One, Six ); -add( Four, Four, Eight ); -mul( Three, Three, Nine ); -add( Nine, One, Ten ); -mul( Nine, Three, TwentySeven ); -mul( Four, Eight, ThirtyTwo ); -mul( Four, Five, t ); -mul( t, Three, t ); -mul( t, Four, TwoForty ); -mov( One, MinusOne ); -neg( MinusOne ); -div( Two, One, Half ); -add( One, Half, OneAndHalf ); -ErrCnt[Failure] = 0; -ErrCnt[Serious] = 0; -ErrCnt[Defect] = 0; -ErrCnt[Flaw] = 0; -PageNo = 1; -#ifndef NOSIGNAL -signal( SIGFPE, sigfpe ); -#endif -printf("Program is now RUNNING tests on small integers:\n"); - -add( Zero, Zero, t ); -if( cmp( t, Zero ) != 0) - { - printf( "0+0 != 0\n" ); - ErrCnt[Failure] += 1; - } -sub( One, One, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "1-1 != 0\n" ); - ErrCnt[Failure] += 1; - } -if( cmp( One, Zero ) <= 0 ) - { - printf( "1 <= 0\n" ); - ErrCnt[Failure] += 1; - } -add( One, One, t ); -if( cmp( t, Two ) != 0 ) - { - printf( "1+1 != 2\n" ); - ErrCnt[Failure] += 1; - } -mov( Zero, Z ); -neg( Z ); -FLOOR( Z, t ); -if( cmp(t,Zero) != 0 ) - { - ErrCnt[Serious] += 1; - printf( "FLOOR(-0) should equal 0, is = " ); - show( t ); - } -if( cmp(Z, Zero) != 0) - { - ErrCnt[Failure] += 1; - printf("Comparison alleges that -0.0 is Non-zero!\n"); - } -else - { - div( TwoForty, One, U1 ); /* U1 = 0.001 */ - mov( One, Radix ); - TstPtUf(); - } -add( Two, One, t ); -if( cmp( t, Three ) != 0 ) - { - printf( "2+1 != 3\n" ); - ErrCnt[Failure] += 1; - } -add( Three, One, t ); -if( cmp( t, Four ) != 0 ) - { - printf( "3+1 != 4\n" ); - ErrCnt[Failure] += 1; - } -mov( Two, t ); -neg( t ); -mul( Two, t, t ); -add( Four, t, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "4+2*(-2) != 0\n" ); - ErrCnt[Failure] += 1; - } -sub( Three, Four, t ); -sub( One, t, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "4-3-1 != 0\n" ); - ErrCnt[Failure] += 1; - } - sub( One, Zero, t ); -if( cmp( t, MinusOne ) != 0 ) - { - printf( "-1 != 0-1\n" ); - ErrCnt[Failure] += 1; - } -add( One, MinusOne, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "1+(-1) != 0\n" ); - ErrCnt[Failure] += 1; - } -mov( One, t ); -FABS( t ); -add( MinusOne, t, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "-1+abs(1) != 0\n" ); - ErrCnt[Failure] += 1; - } -mul( MinusOne, MinusOne, t ); -add( MinusOne, t, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "-1+(-1)*(-1) != 0\n" ); - ErrCnt[Failure] += 1; - } -add( Half, MinusOne, t ); -add( Half, t, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "1/2 + (-1) + 1/2 != 0\n" ); - ErrCnt[Failure] += 1; - } -Milestone = 10; -mul( Three, Three, t ); -if( cmp( t, Nine ) != 0 ) - { - printf( "3*3 != 9\n" ); - ErrCnt[Failure] += 1; - } -mul( Nine, Three, t ); -if( cmp( t, TwentySeven ) != 0 ) - { - printf( "3*9 != 27\n" ); - ErrCnt[Failure] += 1; - } -add( Four, Four, t ); -if( cmp( t, Eight ) != 0 ) - { - printf( "4+4 != 8\n" ); - ErrCnt[Failure] += 1; - } -mul( Eight, Four, t ); -if( cmp( t, ThirtyTwo ) != 0 ) - { - printf( "8*4 != 32\n" ); - ErrCnt[Failure] += 1; - } -sub( TwentySeven, ThirtyTwo, t ); -sub( Four, t, t ); -sub( One, t, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "32-27-4-1 != 0\n" ); - ErrCnt[Failure] += 1; - } -add( Four, One, t ); -if( cmp( t, Five ) != 0 ) - { - printf( "4+1 != 5\n" ); - ErrCnt[Failure] += 1; - } -mul( Four, Five, t ); -mul( Three, t, t ); -mul( Four, t, t ); -if( cmp( t, TwoForty ) != 0 ) - { - printf( "4*5*3*4 != 240\n" ); - ErrCnt[Failure] += 1; - } -div( Three, TwoForty, t ); -mul( Four, Four, t2 ); -mul( Five, t2, t2 ); -sub( t2, t2, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "240/3 - 4*4*5 != 0\n" ); - ErrCnt[Failure] += 1; - } -div( Four, TwoForty, t ); -mul( Five, Three, t2 ); -mul( Four, t2, t2 ); -sub( t2, t, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "240/4 - 5*3*4 != 0\n" ); - ErrCnt[Failure] += 1; - } -div( Five, TwoForty, t ); -mul( Four, Three, t2 ); -mul( Four, t2, t2 ); -sub( t2, t, t ); -if( cmp( t, Zero ) != 0 ) - { - printf( "240/5 - 4*3*4 != 0\n" ); - ErrCnt[Failure] += 1; - } -if(ErrCnt[Failure] == 0) - { -printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n\n"); - } -printf("Searching for Radix and Precision.\n"); -mov( One, W ); -do - { - add( W, W, W ); - add( W, One, Y ); - sub( W, Y, Z ); - sub( One, Z, Y ); - mov( Y, t ); - FABS(t); - add( MinusOne, t, t ); - k = cmp( t, Zero ); - } -while( k < 0 ); -/*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/ -mov( Zero, Precision ); -mov( One, Y ); -do - { - add( W, Y, Radix ); - add( Y, Y, Y ); - sub( W, Radix, Radix ); - k = cmp( Radix, Zero ); - } -while( k == 0); - -if( cmp(Radix, Two) < 0 ) - mov( One, Radix ); -printf("Radix = " ); -show( Radix ); -if( cmp(Radix, One) != 0) - { - mov( One, W ); - do - { - add( One, Precision, Precision ); - mul( W, Radix, W ); - add( W, One, Y ); - sub( W, Y, t ); - k = cmp( t, One ); - } - while( k == 0 ); - } -/* now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1 */ -div( W, One, U1 ); -mul( Radix, U1, U2 ); -printf( "Closest relative separation found is U 1 = " ); -show( U1 ); -printf( "Recalculating radix and precision." ); - -/*save old values*/ -mov( Radix, E0 ); -mov( U1, E1 ); -mov( U2, E9 ); -mov( Precision, E3 ); - -div( Three, Four, X ); -sub( One, X, Third ); -sub( Third, Half, F6 ); -add( F6, F6, X ); -sub( Third, X, X ); -FABS( X ); -if( cmp(X, U2) < 0 ) - mov( U2, X ); - -/*... now X = (unknown no.) ulps of 1+...*/ -do - { - mov( X, U2 ); -/* Y = Half * U2 + ThirtyTwo * U2 * U2; */ - mul( ThirtyTwo, U2, t ); - mul( t, U2, t ); - mul( Half, U2, Y ); - add( t, Y, Y ); - add( One, Y, Y ); - sub( One, Y, X ); - k = cmp( U2, X ); - k2 = cmp( X, Zero ); - } -while ( ! ((k <= 0) || (k2 <= 0))); - -/*... now U2 == 1 ulp of 1 + ... */ -div( Three, Two, X ); -sub( Half, X, F6 ); -add( F6, F6, Third ); -sub( Half, Third, X ); -add( F6, X, X ); -FABS( X ); -if( cmp(X, U1) < 0 ) - mov( U1, X ); - -/*... now X == (unknown no.) ulps of 1 -... */ -do - { - mov( X, U1 ); - /* Y = Half * U1 + ThirtyTwo * U1 * U1;*/ - mul( ThirtyTwo, U1, t ); - mul( U1, t, t ); - mul( Half, U1, Y ); - add( t, Y, Y ); - sub( Y, Half, Y ); - add( Half, Y, X ); - sub( X, Half, Y ); - add( Half, Y, X ); - k = cmp( U1, X ); - k2 = cmp( X, Zero ); - } while ( ! ((k <= 0) || (k2 <= 0))); -/*... now U1 == 1 ulp of 1 - ... */ -if( cmp( U1, E1 ) == 0 ) - printf("confirms closest relative separation U1 .\n"); -else - { - printf("gets better closest relative separation U1 = " ); - show( U1 ); - } -div( U1, One, W ); -sub( U1, Half, F9 ); -add( F9, Half, F9 ); -div( U1, U2, t ); -div( TwoForty, One, t2 ); -add( t2, t, t ); -FLOOR( t, Radix ); -if( cmp(Radix, E0) == 0 ) - printf("Radix confirmed.\n"); -else - { - printf("MYSTERY: recalculated Radix = " ); - show( Radix ); - mov( E0, Radix ); - } -add( Eight, Eight, t ); -if( cmp( Radix, t ) > 0 ) - { - printf( "Radix is too big: roundoff problems\n" ); - ErrCnt[Defect] += 1; - } -k = 1; -if( cmp( Radix, Two ) == 0 ) - k = 0; -if( cmp( Radix, Ten ) == 0 ) - k = 0; -if( cmp( Radix, One ) == 0 ) - k = 0; -if( k != 0 ) - { - printf( "Radix is not as good as 2 or 10\n" ); - ErrCnt[Flaw] += 1; - } -/*=============================================*/ -Milestone = 20; -/*=============================================*/ -sub( Half, F9, t ); -if( cmp( t, Half ) >= 0 ) - { - printf( "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?\n" ); - ErrCnt[Failure] += 1; - } -mov( F9, X ); -I = 1; -sub( Half, X, Y ); -sub( Half, Y, Z ); -if( (cmp( X, One ) == 0) && (cmp( Z, Zero) != 0) ) - { - printf( "Comparison is fuzzy ,X=1 but X-1/2-1/2 != 0\n" ); - ErrCnt[Failure] += 1; - } -add( One, U2, X ); -I = 0; -/*=============================================*/ -Milestone = 25; -/*=============================================*/ -/*... BMinusU2 = nextafter(Radix, 0) */ - -sub( One, Radix, BMinusU2 ); -sub( U2, BMinusU2, t ); -add( One, t, BMinusU2 ); -/* Purify Integers */ -if( cmp(Radix,One) != 0 ) - { -/*X = - TwoForty * LOG(U1) / LOG(Radix);*/ - LOG( U1, X ); - LOG( Radix, t ); - div( t, X, X ); - mul( TwoForty, X, X ); - neg( X ); - - add( Half, X, Y ); - FLOOR( Y, Y ); - sub( Y, X, t ); - FABS( t ); - mul( Four, t, t ); - if( cmp( t, One ) < 0 ) - mov( Y, X ); - div( TwoForty, X, Precision ); - add( Half, Precision, Y ); - FLOOR( Y, Y ); - sub( Y, Precision, t ); - FABS( t ); - mul( TwoForty, t, t ); - if( cmp( t, Half ) < 0 ) - mov( Y, Precision ); - } -FLOOR( Precision, t ); -if( (cmp( Precision, t ) != 0) || (cmp( Radix, One ) == 0) ) - { - printf("Precision cannot be characterized by an Integer number\n"); - printf("of significant digits but, by itself, this is a minor flaw.\n"); - } -if( cmp(Radix, One) == 0 ) - printf("logarithmic encoding has precision characterized solely by U1.\n"); -else - { - printf("The number of significant digits of the Radix is " ); - show( Precision ); - } -mul( U2, Nine, t ); -mul( Nine, t, t ); -mul( TwoForty, t, t ); -if( cmp( t, One ) >= 0 ) - { - printf( "Precision worse than 5 decimal figures\n" ); - ErrCnt[Serious] += 1; - } -/*=============================================*/ -Milestone = 30; -/*=============================================*/ -/* Test for extra-precise subepressions has been deleted. */ -Milestone = 35; -/*=============================================*/ -if( cmp(Radix,Two) >= 0 ) - { - mul( Radix, Radix, t ); - div( t, W, X ); - add( X, One, Y ); - sub( X, Y, Z ); - add( Z, U2, T ); - sub( Z, T, X ); - if( cmp( X, U2 ) != 0 ) - { - printf( "Subtraction is not normalized X=Y,X+Z != Y+Z!\n" ); - ErrCnt[Failure] += 1; - } - if( cmp(X,U2) == 0 ) - printf("Subtraction appears to be normalized, as it should be."); - } - -printf("\nChecking for guard digit in *, /, and -.\n"); -mul( F9, One, Y ); -mul( One, F9, Z ); -sub( Half, F9, X ); -sub( Half, Y, Y ); -sub( X, Y, Y ); -sub( Half, Z, Z ); -sub( X, Z, Z ); -add( One, U2, X ); -mul( X, Radix, T ); -mul( Radix, X, R ); -sub( Radix, T, X ); -mul( Radix, U2, t ); -sub( t, X, X ); -sub( Radix, R, T ); -mul( Radix, U2, t ); -sub( t, T, T ); -sub( One, Radix, t ); -mul( t, X, X ); -sub( One, Radix, t ); -mul( t, T, T ); - -k = cmp(X,Zero); -k |= cmp(Y,Zero); -k |= cmp(Z,Zero); -k |= cmp(T,Zero); -if( k == 0 ) - GMult = Yes; -else - { - GMult = No; - ErrCnt[Serious] += 1; - printf( "* lacks a Guard Digit, so 1*X != X\n" ); - } -mul( Radix, U2, Z ); -add( One, Z, X ); -add( X, Z, Y ); -mul( X, X, t ); -sub( t, Y, Y ); -FABS( Y ); -sub( U2, Y, Y ); -sub( U2, One, X ); -sub( U2, X, Z ); -mul( X, X, t ); -sub( t, Z, Z ); -FABS( Z ); -sub( U1, Z, Z ); -if( (cmp(Y,Zero) > 0) || (cmp(Z,Zero) > 0) ) - { - ErrCnt[Failure] += 1; - printf( "* gets too many final digits wrong.\n" ); - } -sub( U2, One, Y ); -add( One, U2, X ); -div( Y, One, Z ); -sub( X, Z, Y ); -div( Three, One, X ); -div( Nine, Three, Z ); -sub( Z, X, X ); -div( TwentySeven, Nine, T ); -sub( T, Z, Z ); -k = cmp( X, Zero ); -k |= cmp( Y, Zero ); -k |= cmp( Z, Zero ); -if( k ) - { - ErrCnt[Defect] += 1; -printf( "Division lacks a Guard Digit, so error can exceed 1 ulp\n" ); -printf( "or 1/3 and 3/9 and 9/27 may disagree\n" ); - } -div( One, F9, Y ); -sub( Half, F9, X ); -sub( Half, Y, Y ); -sub( X, Y, Y ); -add( One, U2, X ); -div( One, X, T ); -sub( X, T, X ); -k = cmp( X, Zero ); -k |= cmp( Y, Zero ); -k |= cmp( Z, Zero ); -if( k == 0 ) - GDiv = Yes; -else - { - GDiv = No; - ErrCnt[Serious] += 1; - printf( "Division lacks a Guard Digit, so X/1 != X\n" ); - } -add( One, U2, X ); -div( X, One, X ); -sub( Half, X, Y ); -sub( Half, Y, Y ); -if( cmp(Y,Zero) >= 0 ) - { - ErrCnt[Serious] += 1; - printf( "Computed value of 1/1.000..1 >= 1\n" ); - } -sub( U2, One, X ); -mul( Radix, U2, Y ); -add( One, Y, Y ); -mul( X, Radix, Z ); -mul( Y, Radix, T ); -div( Radix, Z, R ); -div( Radix, T, StickyBit ); -sub( X, R, X ); -sub( Y, StickyBit, Y ); -k = cmp( X, Zero ); -k |= cmp( Y, Zero ); -if( k ) - { - ErrCnt[Failure] += 1; - printf( "* and/or / gets too many last digits wrong\n" ); - } -sub( U1, One, Y ); -sub( F9, One, X ); -sub( Y, One, Y ); -sub( U2, Radix, T ); -sub( BMinusU2, Radix, Z ); -sub( T, Radix, T ); -k = cmp( X, U1 ); -k |= cmp( Y, U1 ); -k |= cmp( Z, U2 ); -k |= cmp( T, U2 ); -if( k == 0 ) - GAddSub = Yes; -else - { - GAddSub = No; - ErrCnt[Serious] += 1; - printf( "- lacks Guard Digit, so cancellation is obscured\n" ); - } -sub( One, F9, t ); -if( (cmp(F9,One) != 0) && (cmp(t,Zero) >= 0) ) - { - ErrCnt[Serious] += 1; - printf("comparison alleges (1-U1) < 1 although\n"); - printf(" subtration yields (1-U1) - 1 = 0 , thereby vitiating\n"); - printf(" such precautions against division by zero as\n"); - printf(" ... if (X == 1.0) {.....} else {.../(X-1.0)...}\n"); - } -if (GMult == Yes && GDiv == Yes && GAddSub == Yes) - printf(" *, /, and - appear to have guard digits, as they should.\n"); -/*=============================================*/ -Milestone = 40; -/*=============================================*/ -printf("Checking rounding on multiply, divide and add/subtract.\n"); -RMult = Other; -RDiv = Other; -RAddSub = Other; -div( Two, Radix, RadixD2 ); -mov( Two, A1 ); -Done = False; -do - { - mov( Radix, AInvrse ); - do - { - mov( AInvrse, X ); - div( A1, AInvrse, AInvrse ); - FLOOR( AInvrse, t ); - k = cmp( t, AInvrse ); - } - while( ! (k != 0 ) ); - k = cmp( X, One ); - k2 = cmp( A1, Three ); - Done = (k == 0) || (k2 > 0); - if(! Done) - add( Nine, One, A1 ); - } -while( ! (Done)); -if( cmp(X, One) == 0 ) - mov( Radix, A1 ); -div( A1, One, AInvrse ); -mov( A1, X ); -mov( AInvrse, Y ); -Done = False; -do - { - mul( X, Y, Z ); - sub( Half, Z, Z ); - if( cmp( Z, Half ) != 0 ) - { - ErrCnt[Failure] += 1; - printf( "X * (1/X) differs from 1\n" ); - } - k = cmp( X, Radix ); - Done = (k == 0); - mov( Radix, X ); - div( X, One, Y ); - } -while( ! (Done)); - -add( One, U2, Y2 ); -sub( U2, One, YY1 ); -sub( U2, OneAndHalf, X ); -add( OneAndHalf, U2, Y ); -sub( U2, X, Z ); -mul( Z, Y2, Z ); -mul( Y, YY1, T ); -sub( X, Z, Z ); -sub( X, T, T ); -mul( X, Y2, X ); -add( Y, U2, Y ); -mul( Y, YY1, Y ); -sub( OneAndHalf, X, X ); -sub( OneAndHalf, Y, Y ); -k = cmp( X, Zero ); -k |= cmp( Y, Zero ); -k |= cmp( Z, Zero ); -if( cmp( T, Zero ) > 0 ) - k = 1; -if( k == 0 ) - { - add( OneAndHalf, U2, X ); - mul( X, Y2, X ); - sub( U2, OneAndHalf, Y ); - sub( U2, Y, Y ); - add( OneAndHalf, U2, Z ); - add( U2, Z, Z ); - sub( U2, OneAndHalf, T ); - mul( T, YY1, T ); - add( Z, U2, t ); - sub( t, X, X ); - mul( Y, YY1, StickyBit ); - mul( Z, Y2, S ); - sub( Y, T, T ); - sub( Y, U2, Y ); - add( StickyBit, Y, Y ); -/* Z = S - (Z + U2 + U2); */ - add( Z, U2, t ); - add( t, U2, t ); - sub( t, S, Z ); - add( Y2, U2, t ); - mul( t, YY1, StickyBit ); - mul( Y2, YY1, YY1 ); - sub( Y2, StickyBit, StickyBit ); - sub( Half, YY1, YY1 ); - k = cmp( X, Zero ); - k |= cmp( Y, Zero ); - k |= cmp( Z, Zero ); - k |= cmp( T, Zero ); - k |= cmp( StickyBit, Zero ); - k |= cmp( YY1, Half ); - if( k == 0 ) - { - RMult = Rounded; - printf("Multiplication appears to round correctly.\n"); - } - else - { - add( X, U2, t ); - k = cmp( t, Zero ); - if( cmp( Y, Zero ) >= 0 ) - k |= 1; - add( Z, U2, t ); - k |= cmp( t, Zero ); - if( cmp( T, Zero ) >= 0 ) - k |= 1; - add( StickyBit, U2, t ); - k |= cmp( t, Zero ); - if( cmp(YY1, Half) >= 0 ) - k |= 1; - if( k == 0 ) - { - printf("Multiplication appears to chop.\n"); - } - else - { - printf("* is neither chopped nor correctly rounded.\n"); - } - if( (RMult == Rounded) && (GMult == No) ) - printf("Multiplication has inconsistent result"); - } - } -else - printf("* is neither chopped nor correctly rounded.\n"); - -/*=============================================*/ -Milestone = 45; -/*=============================================*/ -add( One, U2, Y2 ); -sub( U2, One, YY1 ); -add( OneAndHalf, U2, Z ); -add( Z, U2, Z ); -div( Y2, Z, X ); -sub( U2, OneAndHalf, T ); -sub( U2, T, T ); -sub( U2, T, Y ); -div( YY1, Y, Y ); -add( Z, U2, Z ); -div( Y2, Z, Z ); -sub( OneAndHalf, X, X ); -sub( T, Y, Y ); -div( YY1, T, T ); -add( OneAndHalf, U2, t ); -sub( t, Z, Z ); -sub( OneAndHalf, U2, t ); -add( t, T, T ); -k = 0; -if( cmp( X, Zero ) > 0 ) - k = 1; -if( cmp( Y, Zero ) > 0 ) - k = 1; -if( cmp( Z, Zero ) > 0 ) - k = 1; -if( cmp( T, Zero ) > 0 ) - k = 1; -if( k == 0 ) - { - div( Y2, OneAndHalf, X ); - sub( U2, OneAndHalf, Y ); - add( U2, OneAndHalf, Z ); - sub( Y, X, X ); - div( YY1, OneAndHalf, T ); - div( YY1, Y, Y ); - add( Z, U2, t ); - sub( t, T, T ); - sub( Z, Y, Y ); - div( Y2, Z, Z ); - add( Y2, U2, YY1 ); - div( Y2, YY1, YY1 ); - sub( OneAndHalf, Z, Z ); - sub( Y2, YY1, Y2 ); - sub( U1, F9, YY1 ); - div( F9, YY1, YY1 ); - k = cmp( X, Zero ); - k |= cmp( Y, Zero ); - k |= cmp( Z, Zero ); - k |= cmp( T, Zero ); - k |= cmp( Y2, Zero ); - sub( Half, YY1, t ); - sub( Half, F9, t2 ); - k |= cmp( t, t2 ); - if( k == 0 ) - { - RDiv = Rounded; - printf("Division appears to round correctly.\n"); - if(GDiv == No) - printf("Division test inconsistent\n"); - } - else - { - k = 0; - if( cmp( X, Zero ) >= 0 ) - k = 1; - if( cmp( Y, Zero ) >= 0 ) - k = 1; - if( cmp( Z, Zero ) >= 0 ) - k = 1; - if( cmp( T, Zero ) >= 0 ) - k = 1; - if( cmp( Y, Zero ) >= 0 ) - k = 1; - sub( Half, YY1, t ); - sub( Half, F9, t2 ); - if( cmp( t, t2 ) >= 0 ) - k = 1; - if( k == 0 ) - { - RDiv = Chopped; - printf("Division appears to chop.\n"); - } - } - } -if(RDiv == Other) - printf("/ is neither chopped nor correctly rounded.\n"); -div( Radix, One, BInvrse ); -mul( BInvrse, Radix, t ); -sub( Half, t, t ); -if( cmp( t, Half ) != 0 ) - { - ErrCnt[Failure] += 1; - printf( "Radix * ( 1 / Radix ) differs from 1\n" ); - } - -Milestone = 50; -/*=============================================*/ -add( F9, U1, t ); -sub( Half, t, t ); -k = cmp( t, Half ); -add( BMinusU2, U2, t ); -sub( One, t, t ); -sub( One, Radix, t2 ); -k |= cmp( t, t2 ); -if( k != 0 ) - { - ErrCnt[Failure] += 1; - printf( "Incomplete carry-propagation in Addition\n" ); - } -mul( U1, U1, X ); -sub( X, One, X ); -sub( U2, One, Y ); -mul( U2, Y, Y ); -add( One, Y, Y ); -sub( Half, F9, Z ); -sub( Half, X, X ); -sub( Z, X, X ); -sub( One, Y, Y ); -if( (cmp(X,Zero) == 0) && (cmp(Y,Zero) == 0) ) - { - RAddSub = Chopped; - printf("Add/Subtract appears to be chopped.\n"); - } -if(GAddSub == Yes) - { - add( Half, U2, X ); - mul( X, U2, X ); - sub( U2, Half, Y ); - mul( Y, U2, Y ); - add( One, X, X ); - add( One, Y, Y ); - add( One, U2, t ); - sub( X, t, X ); - sub( Y, One, Y ); - k = cmp(X,Zero); - if( k ) - printf( "1+U2-[u2(1/2+U2)+1] != 0\n" ); - k2 = cmp(Y,Zero); - if( k2 ) - printf( "1-[U2(1/2-U2)+1] != 0\n" ); - k |= k2; - if( k == 0 ) - { - add( Half, U2, X ); - mul( X, U1, X ); - sub( U2, Half, Y ); - mul( Y, U1, Y ); - sub( X, One, X ); - sub( Y, One, Y ); - sub( X, F9, X ); - sub( Y, One, Y ); - k = cmp(X,Zero); - if( k ) - printf( "F9-[1-U1(1/2+U2)] != 0\n" ); - k2 = cmp(Y,Zero); - if( k2 ) - printf( "1-[1-U1(1/2-U2)] != 0\n" ); - k |= k2; - if( k == 0 ) - { - RAddSub = Rounded; - printf("Addition/Subtraction appears to round correctly.\n"); - if(GAddSub == No) - printf( "Add/Subtract test inconsistent\n"); - } - else - { - printf("Addition/Subtraction neither rounds nor chops.\n"); - } - } - else - printf("Addition/Subtraction neither rounds nor chops.\n"); - } -else - printf("Addition/Subtraction neither rounds nor chops.\n"); - -mov( One, S ); -add( One, Half, X ); -mul( Half, X, X ); -add( One, X, X ); -add( One, U2, Y ); -mul( Y, Half, Y ); -sub( Y, X, Z ); -sub( X, Y, T ); -add( Z, T, StickyBit ); -if( cmp(StickyBit, Zero) != 0 ) - { - mov( Zero, S ); - ErrCnt[Flaw] += 1; - printf( "(X - Y) + (Y - X) is non zero!\n" ); - } -mov( Zero, StickyBit ); -FLOOR( RadixD2, t ); -k2 = cmp( t, RadixD2 ); -k = 1; -if( (GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes) - && (RMult == Rounded) && (RDiv == Rounded) - && (RAddSub == Rounded) && (k2 == 0) ) - { - printf("Checking for sticky bit.\n"); - k = 0; - add( Half, U1, X ); - mul( X, U2, X ); - mul( Half, U2, Y ); - add( One, Y, Z ); - add( One, X, T ); - sub( One, Z, t ); - sub( One, T, t2 ); - if( cmp(t,Zero) > 0 ) - { - k = 1; - printf( "[1+(1/2)U2]-1 > 0\n" ); - } - if( cmp(t2,U2) < 0 ) - { - k = 1; - printf( "[1+U2(1/2+U1)]-1 < U2\n" ); - } - add( T, Y, Z ); - sub( X, Z, Y ); - sub( T, Z, t ); - sub( T, Y, t2 ); - if( cmp(t,U2) < 0 ) - { - k = 1; - printf( "[[1+U2(1/2+U1)]+(1/2)U2]-[1+U2(1/2+U1)] < U2\n" ); - } - if( cmp(t2,Zero) != 0 ) - { - k = 1; - printf( "(1/2)U2-[1+U2(1/2+U1)] != 0\n" ); - } - add( Half, U1, X ); - mul( X, U1, X ); - mul( Half, U1, Y ); - sub( Y, One, Z ); - sub( X, One, T ); - sub( One, Z, t ); - sub( F9, T, t2 ); - if( cmp(t,Zero) != 0 ) - { - k = 1; - printf( "(1-(1/2)U1)-1 != 0\n" ); - } - if( cmp(t2,Zero) != 0 ) - { - k = 1; - printf( "[1-U1(1/2+U1)]-F9 != 0\n" ); - } - sub( U1, Half, Z ); - mul( Z, U1, Z ); - sub( Z, F9, T ); - sub( Y, F9, Q ); - sub( F9, T, t ); - if( cmp( t, Zero ) != 0 ) - { - k = 1; - printf( "[F9-U1(1/2-U1)]-F9 != 0\n" ); - } - sub( U1, F9, t ); - sub( Q, t, t ); - if( cmp( t, Zero ) != 0 ) - { - k = 1; - printf( "(F9-U1)-(F9-(1/2)U1) != 0\n" ); - } - add( One, U2, Z ); - mul( Z, OneAndHalf, Z ); - add( OneAndHalf, U2, T ); - sub( Z, T, T ); - add( U2, T, T ); - div( Radix, Half, X ); - add( One, X, X ); - mul( Radix, U2, Y ); - add( One, Y, Y ); - mul( X, Y, Z ); - if( cmp( T, Zero ) != 0 ) - { - k = 1; - printf( "(3/2+U2)-3/2(1+U2)+U2 != 0\n" ); - } - mul( Radix, U2, t ); - add( X, t, t ); - sub( Z, t, t ); - if( cmp( t, Zero ) != 0 ) - { - k = 1; - printf( "(1+1/2Radix)+Radix*U2-[1+1/(2Radix)][1+Radix*U2] != 0\n" ); - } - if( cmp(Radix, Two) != 0 ) - { - add( Two, U2, X ); - div( Two, X, Y ); - sub( One, Y, t ); - if( cmp( t, Zero) != 0 ) - k = 1; - } - } -if( k == 0 ) - { - printf("Sticky bit apparently used correctly.\n"); - mov( One, StickyBit ); - } -else - { - printf("Sticky bit used incorrectly or not at all.\n"); - } - -if( GMult == No || GDiv == No || GAddSub == No || - RMult == Other || RDiv == Other || RAddSub == Other) - { - ErrCnt[Flaw] += 1; - printf("lack(s) of guard digits or failure(s) to correctly round or chop\n"); -printf( "(noted above) count as one flaw in the final tally below\n" ); - } -/*=============================================*/ -Milestone = 60; -/*=============================================*/ -printf("\n"); -printf("Does Multiplication commute? "); -printf("Testing on %d random pairs.\n", NoTrials); -SQRT( Three, Random9 ); -mov( Third, Random1 ); -I = 1; -do - { - Random(); - mov( Random1, X ); - Random(); - mov( Random1, Y ); - mul( Y, X, Z9 ); - mul( X, Y, Z ); - sub( Z9, Z, Z9 ); - I = I + 1; - } -while ( ! ((I > NoTrials) || (cmp(Z9,Zero) != 0))); -if(I == NoTrials) - { - div( Three, Half, t ); - add( One, t, Random1 ); - add( U2, U1, t ); - add( t, One, Random2 ); - mul( Random1, Random2, Z ); - mul( Random2, Random1, Y ); -/* Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half / - * Three) * ((U2 + U1) + One); - */ - div( Three, Half, t2 ); - add( One, t2, t2 ); - add( U2, U1, t ); - add( t, One, t ); - mul( t2, t, Z9 ); - mul( t2, t, t ); - sub( t, Z9, Z9 ); - } -if(! ((I == NoTrials) || (cmp(Z9,Zero) == 0))) - { - ErrCnt[Defect] += 1; - printf( "X * Y == Y * X trial fails.\n"); - } -else - { - printf(" No failures found in %d integer pairs.\n", NoTrials); - } -/*=============================================*/ -Milestone = 70; -/*=============================================*/ -sqtest(); -Milestone = 90; -pow1test(); - -Milestone = 110; - -printf("Seeking Underflow thresholds UfThold and E0.\n"); -mov( U1, D ); -FLOOR( Precision, t ); -if( cmp(Precision, t) != 0 ) - { - mov( BInvrse, D ); - mov( Precision, X ); - do - { - mul( D, BInvrse, D ); - sub( One, X, X ); - } - while( cmp(X, Zero) > 0 ); - } -mov( One, Y ); -mov( D, Z ); -/* ... D is power of 1/Radix < 1. */ -sigsave = sigfpe; -if( setjmp(ovfl_buf) ) - goto under0; -do - { - mov( Y, C ); - mov( Z, Y ); - mul( Y, Y, Z ); - add( Z, Z, t ); - } -while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) ); - -under0: -sigsave = 0; - -mov( C, Y ); -mul( Y, D, Z ); -sigsave = sigfpe; -if( setjmp(ovfl_buf) ) - goto under1; -do - { - mov( Y, C ); - mov( Z, Y ); - mul( Y, D, Z ); - add( Z, Z, t ); - } -while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) ); - -under1: -sigsave = 0; - -if( cmp(Radix,Two) < 0 ) - mov( Two, HInvrse ); -else - mov( Radix, HInvrse ); -div( HInvrse, One, H ); -/* ... 1/HInvrse == H == Min(1/Radix, 1/2) */ -div( C, One, CInvrse ); -mov( C, E0 ); -mul( E0, H, Z ); -/* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */ -sigsave = sigfpe; -if( setjmp(ovfl_buf) ) - goto under2; -do - { - mov( E0, Y ); - mov( Z, E0 ); - mul( E0, H, Z ); - add( Z, Z, t ); - } -while( (cmp(E0,Z) > 0) && (cmp(t,Z) > 0) ); - -under2: -sigsave = 0; - -mov( E0, UfThold ); -mov( Zero, E1 ); -mov( Zero, Q ); -mov( U2, E9 ); -add( One, E9, S ); -mul( C, S, D ); -if( cmp(D,C) <= 0 ) - { - mul( Radix, U2, E9 ); - add( One, E9, S ); - mul( C, S, D ); - if( cmp(D, C) <= 0 ) - { - ErrCnt[Failure] += 1; - printf( "multiplication gets too many last digits wrong.\n" ); - mov( E0, Underflow ); - mov( Zero, YY1 ); - mov( Z, PseudoZero ); - } - } -else - { - mov( D, Underflow ); - mul( Underflow, H, PseudoZero ); - mov( Zero, UfThold ); - do - { - mov( Underflow, YY1 ); - mov( PseudoZero, Underflow ); - add( E1, E1, t ); - if( cmp(t, E1) <= 0) - { - mul( Underflow, HInvrse, Y2 ); - sub( Y2, YY1, E1 ); - FABS( E1 ); - mov( YY1, Q ); - if( (cmp( UfThold, Zero ) == 0) - && (cmp(YY1, Y2) != 0) ) - mov( YY1, UfThold ); - } - mul( PseudoZero, H, PseudoZero ); - add( PseudoZero, PseudoZero, t ); - } - while( (cmp(Underflow, PseudoZero) > 0) - && (cmp(t, PseudoZero) > 0) ); - } -/* Comment line 4530 .. 4560 */ -if( cmp(PseudoZero, Zero) != 0 ) - { - printf("\n"); - mov(PseudoZero, Z ); -/* ... Test Pseudo