Fixed version of nextafter[f]/c for 80386 math lib v2.0

[Glenn Geers]

Hi,
	nextafter[f]() was overflowing the 80387 register set if called repeatedly.
A fixed version is attached to this posting.
					Glenn


#!/bin/sh
# This is a shell archive (shar 3.47)
# made 12/17/1990 19:29 UTC by glenn at trantor
# Source directory /tmp
#
# existing files will NOT be overwritten unless -c is specified
#
# This shar contains:
# length  mode       name
# ------ ---------- ------------------------------------------
#   1743 -rw-r--r-- nextafter.c
#   1735 -rw-r--r-- nextafterf.c
#
# ============= nextafter.c ==============
if test -f 'nextafter.c' -a X"$1" != X"-c"; then
	echo 'x - skipping nextafter.c (File already exists)'
else
echo 'x - extracting nextafter.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'nextafter.c' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
** A mix of C and assembler - well I've got the functions so I might 
** as well use them!
**
*/
X
#include "fpumath.h"
X
asm(".align 4");
asm(".Lulp:");
asm(".double 1.1102230246251565e-16");
X
asm(".align 4");
asm(".Lulpup:");
asm(".double 2.2204460492503131e-16");
X
double
nextafter(x, y)
double x, y;
{
X	asm("subl $8, %esp");
X
X	if (isnan(x) || isnan(y))
X		return(quiet_nan(1.0));
X
X	if (isinf(x))
X		if (y > x)
X			return(-max_normal());
X		else
X			if (y < x)
X				return(max_normal());
X
X	if (x == 0.0)
X		if (y > 0.0)
X			return(min_subnormal());
X		else
X			return(-min_subnormal());
X
X	if (isnormal(x)) {
X		if ((x == min_normal()) && y < x)
X			return(max_subnormal());
X
X		if ((x == max_normal()) && y > x)
X			return(infinity());
X
X		if ((x == -max_normal()) && y < x)
X			return(-infinity());
X
X		asm("movl 12(%ebp), %eax");
X		asm("andl $0x7ff00000, %eax");
X		asm("movl %eax, -12(%ebp)");
X		asm("movl $0x0, -16(%ebp)");
X		asm("fincstp");
X
X		if (fabs(x) <= 2.0 && y < x)
X			asm("fldl .Lulp");
X		else
X			asm("fldl .Lulpup");
X
X		asm("fldl -16(%ebp)");
X		asm("fmulp");
X
X		if (y > x) {
X			asm("fldl 8(%ebp)");
X			asm("faddp");
X			asm("leave");
X			asm("ret");
X		}
X		if (y < x) {
X			asm("fldl 8(%ebp)");
X			asm("fsubp");
X			asm("leave");
X			asm("ret");
X		}
X	}
X	else
X	if (issubnormal(x)) {
X		if ((x == max_subnormal()) && y > x)
X			return(min_normal());
X
X		if ((x == -max_subnormal()) && y < x)
X			return(-min_normal());
X
X		if (y > x) 
X			return(x + min_subnormal());
X
X		if (y < x)
X			return(x - min_subnormal());
X	}
X	
X	return(x);
}
SHAR_EOF
chmod 0644 nextafter.c ||
echo 'restore of nextafter.c failed'
Wc_c="`wc -c < 'nextafter.c'`"
test 1743 -eq "$Wc_c" ||
	echo 'nextafter.c: original size 1743, current size' "$Wc_c"
fi
# ============= nextafterf.c ==============
if test -f 'nextafterf.c' -a X"$1" != X"-c"; then
	echo 'x - skipping nextafterf.c (File already exists)'
else
echo 'x - extracting nextafterf.c (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'nextafterf.c' &&
/*
** This file is part of the alternative 80386 math library and is
** covered by the GNU General Public license with my modification
** as noted in the README file that accompanied this file.
**
** Copyright 1990 G. Geers
**
** A mix of C and assembler - well I've got the functions so I might 
** as well use them!
**
*/
X
#include "fpumath.h"
X
asm(".align 4");
asm(".Lulp:");
asm(".double 5.9604644775390625e-08");
X
asm(".align 4");
asm(".Lulpup:");
asm(".double 1.1920928955078125e-07");
X
float
nextafterf(float x, float y)
{
X	asm("subl $8, %esp");
X
X	if (isnanf(x) || isnanf(y))
X		return(quiet_nanf(1.0));
X
X	if (isinff(x))
X		if (y > x)
X			return(-max_normalf());
X		else
X		if (y < x)
X			return(max_normalf());
X
X	if (x == 0.0) {
X		if (y > 0.0)
X			return(min_subnormalf());
X		else
X			return(-min_subnormalf());
X	}
X
X	if (isnormalf(x)) {
X		if ((x == min_normalf()) && y < x)
X			return(max_subnormalf());
X
X		if ((x == max_normalf()) && y > x)
X			return(infinityf());
X
X		if ((x == -max_normalf()) && y < x)
X			return(-infinityf());
X
X		asm("movl 8(%ebp), %eax");
X		asm("andl $0x7f800000, %eax");
X		asm("movl %eax, -8(%ebp)");
X		asm("fincstp");
X
X		if (fabsf(x) <= 2.0 && y < x) 
X			asm("fldl .Lulp");
X		else
X			asm("fldl .Lulpup");
X
X		asm("flds -8(%ebp)");
X		asm("fmulp");
X
X		if (y > x) {
X			asm("flds 8(%ebp)");
X			asm("faddp");
X			asm("leave");
X			asm("ret");
X		}
X		if (y < x) {
X			asm("flds 8(%ebp)");
X			asm("fsubp");
X			asm("leave");
X			asm("ret");
X		}
X	}
X	else
X	if (issubnormalf(x)) {
X		if ((x == max_subnormalf()) && y > x)
X			return(min_normalf());
X
X		if ((x == -max_subnormalf()) && y < x)
X			return(-min_normalf());
X
X		if (y > x) 
X			return(x + min_subnormalf());
X
X		if (y < x)
X			return(x - min_subnormalf());
X	}
X	
X	return(x);
}
SHAR_EOF
chmod 0644 nextafterf.c ||
echo 'restore of nextafterf.c failed'
Wc_c="`wc -c < 'nextafterf.c'`"
test 1735 -eq "$Wc_c" ||
	echo 'nextafterf.c: original size 1735, current size' "$Wc_c"
fi
exit 0
--
Glenn Geers                       | "So when it's over, we're back to people.
Department of Theoretical Physics |  Just to prove that human touch can have
The University of Sydney          |  no equal."
Sydney NSW 2006 Australia         |  - Basia Trzetrzelewska, 'Prime Time TV'