glenn@qed.physics.su.OZ.AU (Glenn Geers) (12/18/90)
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@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'