From 6dc5d42798b6342a1e58cc17fe424ef1936cf7d9 Mon Sep 17 00:00:00 2001 From: Erik van der Kouwe Date: Thu, 24 Dec 2009 20:22:41 +0000 Subject: [PATCH] Floating point support functions --- commands/i386/asmconv/asm86.h | 2 +- commands/i386/asmconv/emit_gnu.c | 1 + commands/i386/asmconv/parse_ack.c | 1 + include/fenv.h | 21 +++ include/math.h | 30 +++- lib/ack/math/ldexp.c | 2 +- lib/i386/Makefile.in | 1 + lib/i386/math/Makefile.in | 17 ++ lib/i386/math/arch_compare.c | 118 +++++++++++++ lib/i386/math/arch_round.c | 57 +++++++ lib/i386/math/fegetround.c | 23 +++ lib/i386/math/feholdexcept.c | 17 ++ lib/i386/math/fesetround.c | 27 +++ lib/i386/math/fpu_cw.h | 33 ++++ lib/i386/math/fpu_cw.s | 20 +++ lib/i386/math/fpu_round.h | 7 + lib/i386/math/fpu_round.s | 36 ++++ lib/i386/math/fpu_sw.h | 26 +++ lib/i386/math/fpu_sw.s | 38 +++++ lib/math/Makefile.in | 1 + lib/math/compare.c | 39 +++++ lib/math/hugeval.c | 17 +- lib/math/localmath.h | 2 + lib/math/sqrt.c | 1 + man/man3/feholdexcept.3 | 28 +++ man/man3/fesetround.3 | 34 ++++ man/man3/fpclassify.3 | 25 +++ man/man3/islessgreater.3 | 26 +++ man/man3/nearbyint.3 | 28 +++ man/man3/remainder.3 | 20 +++ test/Makefile | 3 +- test/run | 2 +- test/test47.c | 275 ++++++++++++++++++++++++++++++ 33 files changed, 970 insertions(+), 8 deletions(-) create mode 100644 include/fenv.h create mode 100644 lib/i386/math/Makefile.in create mode 100644 lib/i386/math/arch_compare.c create mode 100644 lib/i386/math/arch_round.c create mode 100644 lib/i386/math/fegetround.c create mode 100644 lib/i386/math/feholdexcept.c create mode 100644 lib/i386/math/fesetround.c create mode 100644 lib/i386/math/fpu_cw.h create mode 100644 lib/i386/math/fpu_cw.s create mode 100644 lib/i386/math/fpu_round.h create mode 100644 lib/i386/math/fpu_round.s create mode 100644 lib/i386/math/fpu_sw.h create mode 100644 lib/i386/math/fpu_sw.s create mode 100644 lib/math/compare.c create mode 100755 man/man3/feholdexcept.3 create mode 100644 man/man3/fesetround.3 create mode 100644 man/man3/fpclassify.3 create mode 100644 man/man3/islessgreater.3 create mode 100644 man/man3/nearbyint.3 create mode 100644 man/man3/remainder.3 create mode 100644 test/test47.c diff --git a/commands/i386/asmconv/asm86.h b/commands/i386/asmconv/asm86.h index a85abc4af..7c1876fe1 100644 --- a/commands/i386/asmconv/asm86.h +++ b/commands/i386/asmconv/asm86.h @@ -89,7 +89,7 @@ typedef enum opcode { /* 80486 opcodes, from the i486 reference manual. FSIN, FSINCOS, FSQRT, - FSTD, FSTS, FSTPX, FSTPD, FSTPS, + FSTD, FSTS, FSTP, FSTPX, FSTPD, FSTPS, FSTCW, FSTENV, FSTSW, diff --git a/commands/i386/asmconv/emit_gnu.c b/commands/i386/asmconv/emit_gnu.c index 4ec8ec867..54726be4c 100644 --- a/commands/i386/asmconv/emit_gnu.c +++ b/commands/i386/asmconv/emit_gnu.c @@ -151,6 +151,7 @@ static mnemonic_t mnemtab[] = { { FSTCW, "fnstcw" }, { FSTD, "fstl" }, { FSTENV, "fnstenv" }, + { FSTP, "fstp" }, { FSTPD, "fstpl" }, { FSTPS, "fstps" }, { FSTPX, "fstpt" }, diff --git a/commands/i386/asmconv/parse_ack.c b/commands/i386/asmconv/parse_ack.c index f5518684c..b0c17142e 100644 --- a/commands/i386/asmconv/parse_ack.c +++ b/commands/i386/asmconv/parse_ack.c @@ -160,6 +160,7 @@ static mnemonic_t mnemtab[] = { /* This array is sorted. */ { "fstcw", FSTCW, WORD }, { "fstd", FSTD, WORD }, { "fstenv", FSTENV, WORD }, + { "fstp", FSTP, WORD }, { "fstpd", FSTPD, WORD }, { "fstps", FSTPS, WORD }, { "fstpx", FSTPX, WORD }, diff --git a/include/fenv.h b/include/fenv.h new file mode 100644 index 000000000..34a8ded40 --- /dev/null +++ b/include/fenv.h @@ -0,0 +1,21 @@ +#ifndef _FENV_H +#define _FENV_H + +#include + +#define FE_TONEAREST 1 +#define FE_DOWNWARD 2 +#define FE_UPWARD 3 +#define FE_TOWARDZERO 4 + +typedef struct +{ + u16_t cw; + u16_t sw; +} fenv_t; + +int feholdexcept(fenv_t *envp); +int fegetround(void); +int fesetround(int round); + +#endif diff --git a/include/math.h b/include/math.h index 96d6e61c1..8efaab15c 100644 --- a/include/math.h +++ b/include/math.h @@ -7,11 +7,13 @@ #include #endif -#define HUGE_VAL (__huge_val()) /* may be infinity */ +#define INFINITY (__infinity()) +#define NAN (__qnan()) +#define HUGE_VAL INFINITY /* Function Prototypes. */ -_PROTOTYPE( double __huge_val, (void) ); -_PROTOTYPE( int __IsNan, (double _x) ); +_PROTOTYPE( double __infinity, (void) ); +_PROTOTYPE( double __qnan, (void) ); _PROTOTYPE( double acos, (double _x) ); _PROTOTYPE( double asin, (double _x) ); @@ -40,6 +42,28 @@ _PROTOTYPE( double hypot, (double _x, double _y) ); #ifdef _POSIX_SOURCE /* STD-C? */ #include + +#define FP_INFINITE 1 +#define FP_NAN 2 +#define FP_NORMAL 3 +#define FP_SUBNORMAL 4 +#define FP_ZERO 5 + +_PROTOTYPE( int fpclassify, (double x) ); +_PROTOTYPE( int isfinite, (double x) ); +_PROTOTYPE( int isinf, (double x) ); +_PROTOTYPE( int isnan, (double x) ); +_PROTOTYPE( int isnormal, (double x) ); +_PROTOTYPE( int signbit, (double x) ); +_PROTOTYPE( int isgreater, (double x, double y) ); +_PROTOTYPE( int isgreaterequal, (double x, double y) ); +_PROTOTYPE( int isless, (double x, double y) ); +_PROTOTYPE( int islessequal, (double x, double y) ); +_PROTOTYPE( int islessgreater, (double x, double y) ); +_PROTOTYPE( int isunordered, (double x, double y) ); +_PROTOTYPE( double nearbyint, (double x) ); +_PROTOTYPE( double remainder, (double x, double y) ); +_PROTOTYPE( double trunc, (double x) ); #endif #endif /* _MATH_H */ diff --git a/lib/ack/math/ldexp.c b/lib/ack/math/ldexp.c index 501dac452..a892dcadf 100644 --- a/lib/ack/math/ldexp.c +++ b/lib/ack/math/ldexp.c @@ -14,7 +14,7 @@ ldexp(double fl, int exp) int sign = 1; int currexp; - if (__IsNan(fl)) { + if (isnan(fl)) { errno = EDOM; return fl; } diff --git a/lib/i386/Makefile.in b/lib/i386/Makefile.in index 4a8ed4e64..f1e741226 100644 --- a/lib/i386/Makefile.in +++ b/lib/i386/Makefile.in @@ -2,6 +2,7 @@ SUBDIRS="\ int64 \ + math \ misc \ rts \ string" diff --git a/lib/i386/math/Makefile.in b/lib/i386/math/Makefile.in new file mode 100644 index 000000000..c955d2ae1 --- /dev/null +++ b/lib/i386/math/Makefile.in @@ -0,0 +1,17 @@ +# Makefile for lib/i386/math. + +CFLAGS="-O -D_MINIX -D_POSIX_SOURCE" + +LIBRARIES=libc + +libc_FILES=" \ + arch_compare.c \ + arch_round.c \ + fpu_cw.s \ + fpu_sw.s \ + fpu_round.s \ + fegetround.c \ + feholdexcept.c \ + fesetround.c" + +TYPE=both diff --git a/lib/i386/math/arch_compare.c b/lib/i386/math/arch_compare.c new file mode 100644 index 000000000..5127db3a9 --- /dev/null +++ b/lib/i386/math/arch_compare.c @@ -0,0 +1,118 @@ +#include +#include + +#include "fpu_cw.h" +#include "fpu_sw.h" + +#define FPUSW_FLAG_MASK \ + (FPUSW_CONDITION_C0 | FPUSW_CONDITION_C2 | FPUSW_CONDITION_C3) + +#define DOUBLE_NORMAL_MIN 2.2250738585072013830902327173324e-308 /* 2^-1022 */ + +int fpclassify(double x) +{ + /* use status word returned by fpu_xam to determine type */ + switch (fpu_xam(x) & FPUSW_FLAG_MASK) + { + case FPUSW_CONDITION_C0: + return FP_NAN; + + case FPUSW_CONDITION_C2: + /* + * unfortunately, fxam always operates on long doubles + * regardless of the precision setting. This means some + * subnormals are incorrectly classified as normals, + * since they can be normalized using the additional + * exponent bits available. However, if we already know + * that the number is normal as a long double, finding + * out whether it would be subnormal as a double is just + * a simple comparison. + */ + if (-DOUBLE_NORMAL_MIN < x && x < DOUBLE_NORMAL_MIN) + return FP_SUBNORMAL; + else + return FP_NORMAL; + + case FPUSW_CONDITION_C0 | FPUSW_CONDITION_C2: + return FP_INFINITE; + + case FPUSW_CONDITION_C3: + return FP_ZERO; + + case FPUSW_CONDITION_C3 | FPUSW_CONDITION_C2: + return FP_SUBNORMAL; + } + + /* we don't expect any other case: unsupported, emtpy or reserved */ + assert(0); + return -1; +} + +int signbit(double x) +{ + u16_t sw; + + /* examine and use status word to determine sign */ + sw = fpu_xam(x); + return (sw & FPUSW_CONDITION_C1) ? 1 : 0; +} + +#define FPUSW_GREATER 0 +#define FPUSW_LESS FPUSW_CONDITION_C0 +#define FPUSW_EQUAL FPUSW_CONDITION_C3 +#define FPUSW_UNORDERED \ + (FPUSW_CONDITION_C0 | FPUSW_CONDITION_C2 | FPUSW_CONDITION_C3) + +static int fpcompare(double x, double y, u16_t sw1, u16_t sw2) +{ + u16_t sw; + + /* compare and check sanity */ + sw = fpu_compare(x, y) & FPUSW_FLAG_MASK; + switch (sw) + { + case FPUSW_GREATER: + case FPUSW_LESS: + case FPUSW_EQUAL: + case FPUSW_UNORDERED: + break; + + default: + /* other values are not possible (see IA32 Dev Man) */ + assert(0); + return -1; + } + + /* test whether FPUSW equals either sw1 or sw2 */ + return sw == sw1 || sw == sw2; +} + +int isgreater(double x, double y) +{ + return fpcompare(x, y, FPUSW_GREATER, -1); +} + +int isgreaterequal(double x, double y) +{ + return fpcompare(x, y, FPUSW_GREATER, FPUSW_EQUAL); +} + +int isless(double x, double y) +{ + return fpcompare(x, y, FPUSW_LESS, -1); +} + +int islessequal(double x, double y) +{ + return fpcompare(x, y, FPUSW_LESS, FPUSW_EQUAL); +} + +int islessgreater(double x, double y) +{ + return fpcompare(x, y, FPUSW_LESS, FPUSW_GREATER); +} + +int isunordered(double x, double y) +{ + return fpcompare(x, y, FPUSW_UNORDERED, -1); +} diff --git a/lib/i386/math/arch_round.c b/lib/i386/math/arch_round.c new file mode 100644 index 000000000..a239f3636 --- /dev/null +++ b/lib/i386/math/arch_round.c @@ -0,0 +1,57 @@ +#include +#include +#include + +#include "fpu_cw.h" +#include "fpu_round.h" + +static double rndint(double x, u16_t cw_bits, u16_t cw_mask) +{ + u16_t cw; + + /* set FPUCW to the right value */ + cw = fpu_cw_get(); + fpu_cw_set((cw & cw_mask) | cw_bits); + + /* perform the round */ + fpu_rndint(&x); + + /* restore FPUCW */ + fpu_cw_set(cw); + return x; +} + +double nearbyint(double x) +{ + /* round, disabling floating point precision error */ + return rndint(x, FPUCW_EXCEPTION_MASK_PM, ~0); +} + +double remainder(double x, double y) +{ + int xclass, yclass; + + /* check arguments */ + xclass = fpclassify(x); + yclass = fpclassify(y); + if (xclass == FP_NAN || yclass == FP_NAN) + return NAN; + + if (xclass == FP_INFINITE || yclass == FP_ZERO) + { + errno = EDOM; + return NAN; + } + + /* call the assembly implementation */ + fpu_remainder(&x, y); + return x; +} + +double trunc(double x) +{ + /* round in truncate mode, disabling floating point precision error */ + return rndint(x, + FPUCW_EXCEPTION_MASK_PM | FPUCW_ROUNDING_CONTROL_TRUNC, + ~FPUCW_ROUNDING_CONTROL); +} diff --git a/lib/i386/math/fegetround.c b/lib/i386/math/fegetround.c new file mode 100644 index 000000000..f9225ad41 --- /dev/null +++ b/lib/i386/math/fegetround.c @@ -0,0 +1,23 @@ +#include +#include + +#include "fpu_cw.h" + +int fegetround(void) +{ + u16_t cw; + + /* read and categorize FPUCW */ + cw = fpu_cw_get(); + switch (cw & FPUCW_ROUNDING_CONTROL) + { + case FPUCW_ROUNDING_CONTROL_NEAREST: return FE_TONEAREST; + case FPUCW_ROUNDING_CONTROL_DOWN: return FE_DOWNWARD; + case FPUCW_ROUNDING_CONTROL_UP: return FE_UPWARD; + case FPUCW_ROUNDING_CONTROL_TRUNC: return FE_TOWARDZERO; + } + + /* each case has been handled, otherwise the constants are wrong */ + assert(0); + return -1; +} diff --git a/lib/i386/math/feholdexcept.c b/lib/i386/math/feholdexcept.c new file mode 100644 index 000000000..22184a44d --- /dev/null +++ b/lib/i386/math/feholdexcept.c @@ -0,0 +1,17 @@ +#include +#include + +#include "fpu_cw.h" +#include "fpu_sw.h" + +int feholdexcept(fenv_t *envp) +{ + /* read FPUCW and FPUSW */ + envp->cw = fpu_cw_get(); + envp->sw = fpu_sw_get(); + + /* update FPUCW to block exceptions */ + fpu_cw_set(envp->cw | FPUCW_EXCEPTION_MASK); + + return 0; +} diff --git a/lib/i386/math/fesetround.c b/lib/i386/math/fesetround.c new file mode 100644 index 000000000..131d64446 --- /dev/null +++ b/lib/i386/math/fesetround.c @@ -0,0 +1,27 @@ +#include +#include + +#include "fpu_cw.h" + +int fesetround(int round) +{ + u16_t cw; + + /* read and update FPUCW */ + cw = fpu_cw_get() & ~FPUCW_ROUNDING_CONTROL; + switch (round) + { + case FE_TONEAREST: cw |= FPUCW_ROUNDING_CONTROL_NEAREST; break; + case FE_DOWNWARD: cw |= FPUCW_ROUNDING_CONTROL_DOWN; break; + case FE_UPWARD: cw |= FPUCW_ROUNDING_CONTROL_UP; break; + case FE_TOWARDZERO: cw |= FPUCW_ROUNDING_CONTROL_TRUNC; break; + + default: + errno = EINVAL; + return -1; + } + + /* set FPUCW to the updated value */ + fpu_cw_set(cw); + return 0; +} diff --git a/lib/i386/math/fpu_cw.h b/lib/i386/math/fpu_cw.h new file mode 100644 index 000000000..f981cab94 --- /dev/null +++ b/lib/i386/math/fpu_cw.h @@ -0,0 +1,33 @@ +#ifndef __FPU_CW__ +#define __FPU_CW__ + +#include + +/* + * see section 8.1.5 "x87 FPU Control Word" in "Intel 64 and IA-32 Architectures + * Software Developer's Manual Volume 1 Basic Architecture" + */ +#define FPUCW_EXCEPTION_MASK 0x003f +#define FPUCW_EXCEPTION_MASK_IM 0x0001 +#define FPUCW_EXCEPTION_MASK_DM 0x0002 +#define FPUCW_EXCEPTION_MASK_ZM 0x0004 +#define FPUCW_EXCEPTION_MASK_OM 0x0008 +#define FPUCW_EXCEPTION_MASK_UM 0x0010 +#define FPUCW_EXCEPTION_MASK_PM 0x0020 + +#define FPUCW_PRECISION_CONTROL 0x0300 +#define FPUCW_PRECISION_CONTROL_SINGLE 0x0000 +#define FPUCW_PRECISION_CONTROL_DOUBLE 0x0200 +#define FPUCW_PRECISION_CONTROL_XDOUBLE 0x0300 + +#define FPUCW_ROUNDING_CONTROL 0x0c00 +#define FPUCW_ROUNDING_CONTROL_NEAREST 0x0000 +#define FPUCW_ROUNDING_CONTROL_DOWN 0x0400 +#define FPUCW_ROUNDING_CONTROL_UP 0x0800 +#define FPUCW_ROUNDING_CONTROL_TRUNC 0x0c00 + +/* get and set FPU control word */ +u16_t fpu_cw_get(void); +void fpu_cw_set(u16_t fpu_cw); + +#endif /* !defined(__FPU_CW__) */ diff --git a/lib/i386/math/fpu_cw.s b/lib/i386/math/fpu_cw.s new file mode 100644 index 000000000..f637497ee --- /dev/null +++ b/lib/i386/math/fpu_cw.s @@ -0,0 +1,20 @@ +! fpu_cw_get() - get FPU control word Author: Erik van der Kouwe +! fpu_cw_set() - set FPU control word 9 Dec 2009 +.sect .text +.define _fpu_cw_get +.define _fpu_cw_set + +! u16_t fpu_cw_get(void) +_fpu_cw_get: + ! clear unused bits just to be sure + xor eax, eax + push eax + fstcw (esp) + pop eax + ret + +! void fpu_cw_set(u16_t fpu_cw) +_fpu_cw_set: + ! load control word from parameter + fldcw 4(esp) + ret diff --git a/lib/i386/math/fpu_round.h b/lib/i386/math/fpu_round.h new file mode 100644 index 000000000..d5b482fa5 --- /dev/null +++ b/lib/i386/math/fpu_round.h @@ -0,0 +1,7 @@ +#ifndef __FPU_RNDINT__ +#define __FPU_RNDINT__ + +void fpu_rndint(double *value); +void fpu_remainder(double *x, double y); + +#endif /* !defined(__FPU_RNDINT__) */ diff --git a/lib/i386/math/fpu_round.s b/lib/i386/math/fpu_round.s new file mode 100644 index 000000000..404f2fd82 --- /dev/null +++ b/lib/i386/math/fpu_round.s @@ -0,0 +1,36 @@ +! fpu_rndint() - round integer Author: Erik van der Kouwe +! 17 Dec 2009 +.sect .text +.define _fpu_rndint +.define _fpu_remainder + +! void fpu_rndint(double *value) +_fpu_rndint: + ! move the value onto the floating point stack + mov eax, 4(esp) + fldd (eax) + + ! round it (beware of precision exception!) + frndint + + ! store the result + fstpd (eax) + ret + +! void fpu_remainder(double *x, double y) +_fpu_remainder: + ! move the values onto the floating point stack + fldd 8(esp) + mov edx, 4(esp) + fldd (edx) + + ! compute remainder, multiple iterations may be needed +1: fprem1 + .data1 0xdf, 0xe0 ! fnstsw ax + sahf + jp 1b + + ! store the result and pop the divisor + fstpd (edx) + fstp st + ret diff --git a/lib/i386/math/fpu_sw.h b/lib/i386/math/fpu_sw.h new file mode 100644 index 000000000..e66889ddd --- /dev/null +++ b/lib/i386/math/fpu_sw.h @@ -0,0 +1,26 @@ +#ifndef __FPU_SW__ +#define __FPU_SW__ + +#include + +#define FPUSW_EXCEPTION_IE 0x0001 +#define FPUSW_EXCEPTION_DE 0x0002 +#define FPUSW_EXCEPTION_ZE 0x0004 +#define FPUSW_EXCEPTION_OE 0x0008 +#define FPUSW_EXCEPTION_UE 0x0010 +#define FPUSW_EXCEPTION_PE 0x0020 +#define FPUSW_STACK_FAULT 0x0040 +#define FPUSW_ERROR_SUMMARY 0x0080 +#define FPUSW_CONDITION_C0 0x0100 +#define FPUSW_CONDITION_C1 0x0200 +#define FPUSW_CONDITION_C2 0x0400 +#define FPUSW_CONDITION_C3 0x4000 +#define FPUSW_BUSY 0x8000 + +u16_t fpu_compare(double x, double y); +u16_t fpu_sw_get(void); +void fpu_sw_set(u16_t value); +u16_t fpu_xam(double value); + +#endif /* !defined(__FPU_SW__) */ + diff --git a/lib/i386/math/fpu_sw.s b/lib/i386/math/fpu_sw.s new file mode 100644 index 000000000..3d6f30da5 --- /dev/null +++ b/lib/i386/math/fpu_sw.s @@ -0,0 +1,38 @@ +! fpu_compare() - compare doubles Author: Erik van der Kouwe +! fpu_sw_get() - get FPU status 17 Dec 2009 +! fpu_xam() - examine double +.sect .text +.define _fpu_compare +.define _fpu_sw_get +.define _fpu_xam + +! u16_t fpu_compare(double x, double y) +_fpu_compare: + ! move the values onto the floating point stack + fldd 12(esp) + fldd 4(esp) + + ! compare values and return status word + fcompp + jmp _fpu_sw_get + +! u16_t fpu_sw_get(void) +_fpu_sw_get: + ! clear unused high-order word and get status word + xor eax, eax + .data1 0xdf, 0xe0 ! fnstsw ax + ret + +! u16_t fpu_xam(double value) +_fpu_xam: + ! move the value onto the floating point stack + fldd 4(esp) + + ! examine value and get status word + fxam + call _fpu_sw_get + + ! pop the value + fstp st + ret + diff --git a/lib/math/Makefile.in b/lib/math/Makefile.in index 61f3f999e..434ca242e 100644 --- a/lib/math/Makefile.in +++ b/lib/math/Makefile.in @@ -9,6 +9,7 @@ libc_FILES=" \ atan.c \ atan2.c \ ceil.c \ + compare.c \ exp.c \ fabs.c \ floor.c \ diff --git a/lib/math/compare.c b/lib/math/compare.c new file mode 100644 index 000000000..e3dd91d63 --- /dev/null +++ b/lib/math/compare.c @@ -0,0 +1,39 @@ +#include +#include + +/* functions missing here are architecture-specific and are in i386/float */ + +int isfinite(double x) +{ + /* return value based on classification */ + switch (fpclassify(x)) + { + case FP_INFINITE: + case FP_NAN: + return 0; + + case FP_NORMAL: + case FP_SUBNORMAL: + case FP_ZERO: + return 1; + } + + /* if we get here, fpclassify is buggy */ + assert(0); + return -1; +} + +int isinf(double x) +{ + return fpclassify(x) == FP_INFINITE; +} + +int isnan(double x) +{ + return fpclassify(x) == FP_NAN; +} + +int isnormal(double x) +{ + return fpclassify(x) == FP_NORMAL; +} diff --git a/lib/math/hugeval.c b/lib/math/hugeval.c index 3f4bf13e0..cb9952581 100644 --- a/lib/math/hugeval.c +++ b/lib/math/hugeval.c @@ -9,7 +9,7 @@ #include double -__huge_val(void) +__infinity(void) { #if (CHIP == INTEL) static unsigned char ieee_infinity[] = { @@ -21,3 +21,18 @@ __huge_val(void) return 1.0e+1000; /* This will generate a warning */ #endif } + +double +__qnan(void) +{ +#if (CHIP == INTEL) + static unsigned char ieee_qnan[] = { + 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x7f }; + + assert(sizeof(double) == sizeof(ieee_qnan)); + return *(double *) ieee_qnan; +#else +#error QNaN not defined on this architecture +#endif +} + diff --git a/lib/math/localmath.h b/lib/math/localmath.h index 0b2aaea83..29432af59 100644 --- a/lib/math/localmath.h +++ b/lib/math/localmath.h @@ -3,6 +3,8 @@ */ /* $Header$ */ +#define __IsNan isnan + /* some constants (Hart & Cheney) */ #define M_PI 3.14159265358979323846264338327950288 #define M_2PI 6.28318530717958647692528676655900576 diff --git a/lib/math/sqrt.c b/lib/math/sqrt.c index e02638812..2f2402e39 100644 --- a/lib/math/sqrt.c +++ b/lib/math/sqrt.c @@ -9,6 +9,7 @@ #include #include #include +#include "localmath.h" #define NITER 5 diff --git a/man/man3/feholdexcept.3 b/man/man3/feholdexcept.3 new file mode 100755 index 000000000..9e7257c69 --- /dev/null +++ b/man/man3/feholdexcept.3 @@ -0,0 +1,28 @@ +.TH FEHOLDEXCEPT 3 "December 22, 2009" +.UC 4 +.SH NAME +feholdexcept \- disable floating point exceptions +.SH SYNOPSIS +.nf +.ft B +#include + +int feholdexcept(fenv_t *\fIenvp\fP); +.fi +.SH DESCRIPTION +The feholdexcept function disables floating point exceptions. The floating +point environment prior to disabling exceptions is stored in the struct +pointed to by \fIenvp\fP. +.SH "RETURN VALUE" +On x86, this function never fails and the return value is always 0. A portable +program should, however, consider the possibility that this function may fail +and will return -1 to indicate this. +.SH "KNOWN ISSUES" +POSIX requires the function to clear exception flags. However, as none of the +other functions dealing with FPU flags and FPU state have been implemented this +has no effect on MINIX. A full implementation may need to store more FPU state +as there is no instruction to change the status word independent of the rest +of the FPU state. Do not depend on the fields of fenv_t as they may change even +on future versions of MINIX/x86. + + diff --git a/man/man3/fesetround.3 b/man/man3/fesetround.3 new file mode 100644 index 000000000..97c9109d2 --- /dev/null +++ b/man/man3/fesetround.3 @@ -0,0 +1,34 @@ +.TH FESETROUND 3 "December 18, 2009" +.UC 4 +.SH NAME +fegetround, fesetround \- floating point rounding mode control +.SH SYNOPSIS +.nf +.ft B +#include + +int fegetround(void); +int fesetround(int \fIround\fP); +.fi +.SH DESCRIPTION +These functions control the floating point rounding mode. One can use +fegetround to determine the rounding mode currently in effect and fesetround +the change the rounding mode. Four rounding modes are recognized: FE_TONEAREST, +FE_DOWNWARD, FE_UPWARD and FE_TOWARDZERO. When rounding a value \fIx\fP to an +integer \fIn\fP, the following approaches can be taken: FE_TONEAREST selects +the \fIn\fP for which abs(fIx\fP - \fIn\fP) is minimal, preferring an even +\fIn\fP if there are two possibilities; FE_DOWNWARD selects the largest \fIn\fP +for which \fIn\fP <= fIx\fP; FE_UPWARD selects the smallest \fIn\fP for which +\fIn\fP >= fIx\fP; and FE_TOWARDZERO selects \fIn\fP with the largest +abs(\fIn\fP) for which abs(\fIn\fP) <= abs(\fIx\fP). +.SH "RETURN VALUE" +fegetround returns the current rounding mode. fesetround returns 0 if it was +succesful in adjusting the rounding mode and -1 otherwise. The only possible +reason for failure is an unsupported value for \fIround\fP, which is signalled +by setting errno to EINVAL. +.SH DIAGNOSTICS +fegetround always succeeds. The only possible reason for failure of fesetround +is an unsupported value for \fIround\fP, which is signalled by setting errno to +EINVAL. +.SH "SEE ALSO" +nearbyint(3) diff --git a/man/man3/fpclassify.3 b/man/man3/fpclassify.3 new file mode 100644 index 000000000..9f1fd3f05 --- /dev/null +++ b/man/man3/fpclassify.3 @@ -0,0 +1,25 @@ +.TH FPCLASSIFY 3 "December 18, 2009" +.UC 4 +.SH NAME +fpclassify, isfinite, isinf, isnan, isnormal, signbit \- classify floating point numbers +.SH SYNOPSIS +.nf +.ft B +#include + +int fpclassify(double \fIx\fP): +int isfinite(double \fIx\fP); +int isinf(double \fIx\fP); +int isnan(double \fIx\fP); +int isnormal(double \fIx\fP); +int signbit(double \fIx\fP); +.fi +.SH DESCRIPTION +These functions provide information about the specified floating point number +\fIx\fP. fpclassify returns one of the values FP_INFINITE, FP_NAN, FP_NORMAL, +FP_SUBNORMAL and FP_ZERO depending on the type of number provided. The isinf, +isinf, isnan and isnormal test for specific number classes, returning a +non-zero value is and only if the specified number belongs to the class +specified by the function name. The signbit function returns a non-zero value +if and only if the sign bit is set, which for non-NaN values (including zero) +means that the number is negative. diff --git a/man/man3/islessgreater.3 b/man/man3/islessgreater.3 new file mode 100644 index 000000000..96cba3326 --- /dev/null +++ b/man/man3/islessgreater.3 @@ -0,0 +1,26 @@ +.TH ISLESSGREATER 3 "December 18, 2009" +.UC 4 +.SH NAME +isgreater, isgreaterequal, isless, islessequal, islessgreater, isunordered \- order floating point numbers +.SH SYNOPSIS +.nf +.ft B +#include + +int isgreater(double \fIx\fP, double \fIy\fP); +int isgreaterequal(double \fIx\fP, double \fIy\fP); +int isless(double \fIx\fP, double \fIy\fP); +int islessequal(double \fIx\fP, double \fIy\fP); +int islessgreater(double \fIx\fP, double \fIy\fP); +int isunordered(double \fIx\fP, double \fIy\fP); +.fi +.SH DESCRIPTION +These functions compare the specified floating point numbers \fIx\fP and +\fIy\fP. They differ from the regular comparison operators in that they +recognize and additional 'unordered' relationship between the numbers if one of +them is NaN. As a consequence, these functions do not raise floating point +expceptions and can be safely used with NaN values. +.SH "RETURN VALUE" +The functions return a non-zero value if and only if the relationship specified +in the name holds between \fIx\fP and \fIy\fP. All functions except for +'isunordered' return zero if either of the numbers is NaN. diff --git a/man/man3/nearbyint.3 b/man/man3/nearbyint.3 new file mode 100644 index 000000000..c91303257 --- /dev/null +++ b/man/man3/nearbyint.3 @@ -0,0 +1,28 @@ +.TH NEARBYINT 3 "December 18, 2009" +.UC 4 +.SH NAME +ceil, floor, nearbyint, trunc \- floating point rounding +.SH SYNOPSIS +.nf +.ft B +#include + +double ceil(double \fIx\fP); +double floor(double \fIx\fP); +double nearbyint(double \fIx\fP); +double trunc(double \fIx\fP); +.fi +.SH DESCRIPTION +These functions round the specified floating point number to a nearby integer. +For nearbyint, the rounding mode is determined by the value set using the +fesetround function. The other functions are not influenced by the rounding +mode. The ceil function rounds upwards, selecting the smallest integer that is +larger than or equal to \fIx\fP. The floor function rounds downwards, selecting +the largest integer that is smaller than or equal to \fIx\fP. The trunc function +rounds towards zero, selecting the largest integer with the largest absolute +value of which the absolute value is smaller than or equal to the absolute +value of \fIx\fP. +.SH "RETURN VALUE" +The functions return an integer close to \fIx\fP. +.SH "SEE ALSO" +fesetround(3) diff --git a/man/man3/remainder.3 b/man/man3/remainder.3 new file mode 100644 index 000000000..d84944f0a --- /dev/null +++ b/man/man3/remainder.3 @@ -0,0 +1,20 @@ +.TH REMAINDER 3 "December 18, 2009" +.UC 4 +.SH NAME +remainder \- floating point remainder after division +.SH SYNOPSIS +.nf +.ft B +#include + +double remainder(double \fIx\fP, double \fIy\fP); +.fi +.SH DESCRIPTION +This function returns the remainder of a division of \fIx\fP by \fIy\fP. More +formally, it computes which integer \fIn\fP is closest to the fraction +\fIx\fP/\fIy\fP and returns \fIx\fP - \fIn\fP*\fIy\fP. An even value for +\fIn\fP is preferred over an odd one if both are equally far away (that is, +banker's rounding is applied). If \fIx\fP is infinite or \fIy\fP is zero, +a NaN value is returned. +.SH "RETURN VALUE" +The function returns the remainder of a division of \fIx\fP by \fIy\fP. diff --git a/test/Makefile b/test/Makefile index a275ae6eb..b74829891 100644 --- a/test/Makefile +++ b/test/Makefile @@ -10,7 +10,7 @@ OBJ= test1 test2 test3 test4 test5 test6 test7 test8 test9 \ test21 test22 test23 test25 test26 test27 test28 test29 \ test30 test31 test32 test34 test35 test36 test37 test38 \ test39 t10a t11a t11b test40 t40a t40b t40c t40d t40e t40f test41 \ - test42 test44 test45 + test42 test44 test45 test47 BIGOBJ= test20 test24 ROOTOBJ= test11 test33 test43 test46 @@ -94,4 +94,5 @@ test44: test44.c test45: test45.c test45.h test45-gcc: test45.c test45.h test46: test46.c +test47: test47.c diff --git a/test/run b/test/run index 1dc5ee68b..1455f3b75 100755 --- a/test/run +++ b/test/run @@ -19,7 +19,7 @@ echo " " # Run all the tests, keeping track of who failed. for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 \ 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 \ - 41 42 43 44 45 45-gcc 46 sh1.sh sh2.sh + 41 42 43 44 45 45-gcc 46 47 sh1.sh sh2.sh do if [ -x ./test$i ] then diff --git a/test/test47.c b/test/test47.c new file mode 100644 index 000000000..f75a04df0 --- /dev/null +++ b/test/test47.c @@ -0,0 +1,275 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#define MAX_ERROR 4 +static int errct; + +/* maximum allowed FP difference for our tests */ +#define EPSILON 0.00000000023283064365386962890625 /* 2^(-32) */ + +static void quit(void) +{ + if (errct == 0) + { + printf("ok\n"); + exit(0); + } + else + { + printf("%d errors\n", errct); + exit(1); + } +} + +#define ERR(x, y) e(__LINE__, (x), (y)) + +static void e(int n, double x, double y) +{ + printf("Line %d, x=%.20g, y=%.20g\n", n, x, y); + if (errct++ > MAX_ERROR) + { + printf("Too many errors; test aborted\n"); + exit(1); + } +} + +static void signal_handler(int signum) +{ + struct sigframe *sigframe; + + /* report signal */ + sigframe = (struct sigframe *) ((char *) &signum - + (char *) &((struct sigframe *) NULL)->sf_signo); + printf("Signal %d at 0x%x\n", signum, sigframe->sf_scp->sc_regs.pc); + + /* count as error */ + ERR(0, 0); + fflush(stdout); + + /* handle signa again next time */ + signal(signum, signal_handler); +} + +static void test_fpclassify(double value, int class, int test_sign) +{ + /* test fpclassify */ + if (fpclassify(value) != class) ERR(value, 0); + if (test_sign) + { + if (fpclassify(-value) != class) ERR(-value, 0); + + /* test signbit */ + if (signbit(value)) ERR(value, 0); + if (!signbit(-value)) ERR(-value, 0); + } +} + +/* Maximum normal double: (2 - 2^(-53)) * 2^1023 */ +#define NORMAL_DOUBLE_MAX 1.797693134862315708145274237317e+308 + +/* Minimum normal double: 2^(-1022) */ +#define NORMAL_DOUBLE_MIN 2.2250738585072013830902327173324e-308 + +/* Maximum subnormal double: (1 - 2^(-53)) * 2^(-1022) */ +#define SUBNORMAL_DOUBLE_MAX 2.2250738585072008890245868760859e-308 + +/* Minimum subnormal double: 2^(-52) * 2^(-1023) */ +#define SUBNORMAL_DOUBLE_MIN 4.9406564584124654417656879286822e-324 + +static void test_fpclassify_values(void) +{ + double d; + char negzero[] = { 0, 0, 0, 0, 0, 0, 0, 0x80 }; + + /* test some corner cases for fpclassify and signbit */ + test_fpclassify(INFINITY, FP_INFINITE, 1); + test_fpclassify(NAN, FP_NAN, 0); + test_fpclassify(0, FP_ZERO, 0); + test_fpclassify(1, FP_NORMAL, 1); + test_fpclassify(NORMAL_DOUBLE_MAX, FP_NORMAL, 1); + test_fpclassify(NORMAL_DOUBLE_MIN, FP_NORMAL, 1); + test_fpclassify(SUBNORMAL_DOUBLE_MAX, FP_SUBNORMAL, 1); + test_fpclassify(SUBNORMAL_DOUBLE_MIN, FP_SUBNORMAL, 1); + + /* + * unfortunately the minus operator does not change the sign of zero, + * so a special case is needed to test it + */ + assert(sizeof(negzero) == sizeof(double)); + test_fpclassify(*(double *) negzero, FP_ZERO, 0); + if (!signbit(*(double *) negzero)) ERR(0, 0); + + /* test other small numbers for fpclassify and signbit */ + d = 1; + while (d >= NORMAL_DOUBLE_MIN) + { + test_fpclassify(d, FP_NORMAL, 1); + d /= 10; + } + while (d >= SUBNORMAL_DOUBLE_MIN) + { + test_fpclassify(d, FP_SUBNORMAL, 1); + d /= 10; + } + test_fpclassify(d, FP_ZERO, 0); + + /* test other large numbers for fpclassify and signbit */ + d = 1; + while (d <= NORMAL_DOUBLE_MAX / 10) + { + test_fpclassify(d, FP_NORMAL, 1); + d *= 10; + } +} + +/* expected rounding: up, down or identical */ +#define ROUND_EQ 1 +#define ROUND_DN 2 +#define ROUND_UP 3 + +static void test_round_value_mode_func(double value, int mode, double (*func)(double), int exp) +{ + int mode_old; + double rounded; + + /* update and check rounding mode */ + mode_old = fegetround(); + fesetround(mode); + if (fegetround() != mode) ERR(0, 0); + + /* perform rounding */ + rounded = func(value); + + /* check direction of rounding */ + switch (exp) + { + case ROUND_EQ: if (rounded != value) ERR(value, rounded); break; + case ROUND_DN: if (rounded >= value) ERR(value, rounded); break; + case ROUND_UP: if (rounded <= value) ERR(value, rounded); break; + default: assert(0); + } + + /* check whether the number is sufficiently close */ + if (fabs(value - rounded) >= 1) ERR(value, rounded); + + /* check whether the number is integer */ + if (remainder(rounded, 1)) ERR(value, rounded); + + /* re-check and restore rounding mode */ + if (fegetround() != mode) ERR(0, 0); + fesetround(mode_old); +} + +static void test_round_value_mode(double value, int mode, int exp_nearbyint, + int exp_ceil, int exp_floor, int exp_trunc) +{ + /* test both nearbyint and trunc */ + test_round_value_mode_func(value, mode, nearbyint, exp_nearbyint); + test_round_value_mode_func(value, mode, ceil, exp_ceil); + test_round_value_mode_func(value, mode, floor, exp_floor); + test_round_value_mode_func(value, mode, trunc, exp_trunc); +} + +static void test_round_value(double value, int exp_down, int exp_near, int exp_up, int exp_trunc) +{ + /* test each rounding mode */ + test_round_value_mode(value, FE_DOWNWARD, exp_down, exp_up, exp_down, exp_trunc); + test_round_value_mode(value, FE_TONEAREST, exp_near, exp_up, exp_down, exp_trunc); + test_round_value_mode(value, FE_UPWARD, exp_up, exp_up, exp_down, exp_trunc); + test_round_value_mode(value, FE_TOWARDZERO, exp_trunc, exp_up, exp_down, exp_trunc); +} + +static void test_round_values(void) +{ + int i; + + /* test various values */ + for (i = -100000; i < 100000; i++) + { + test_round_value(i + 0.00, ROUND_EQ, ROUND_EQ, ROUND_EQ, ROUND_EQ); + test_round_value(i + 0.25, ROUND_DN, ROUND_DN, ROUND_UP, (i >= 0) ? ROUND_DN : ROUND_UP); + test_round_value(i + 0.50, ROUND_DN, (i % 2) ? ROUND_UP : ROUND_DN, ROUND_UP, (i >= 0) ? ROUND_DN : ROUND_UP); + test_round_value(i + 0.75, ROUND_DN, ROUND_UP, ROUND_UP, (i >= 0) ? ROUND_DN : ROUND_UP); + } +} + +static void test_remainder_value(double x, double y) +{ + int mode_old; + double r1, r2, z; + + assert(y != 0); + + /* compute remainder using the function */ + r1 = remainder(x, y); + + /* compute remainder using alternative approach */ + mode_old = fegetround(); + fesetround(FE_TONEAREST); + r2 = x - nearbyint(x / y) * y; + fesetround(mode_old); + + /* Compare results */ + if (fabs(r1 - r2) > EPSILON && fabs(r1 + r2) > EPSILON) + { + printf("%.20g != %.20g\n", r1, r2); + ERR(x, y); + } +} + +static void test_remainder_values_y(double x) +{ + int i, j; + + /* try various rational and transcendental values for y (except zero) */ + for (i = -50; i <= 50; i++) + if (i != 0) + for (j = 1; j < 10; j++) + { + test_remainder_value(x, (double) i / (double) j); + test_remainder_value(x, i * M_E + j * M_PI); + } +} + +static void test_remainder_values_xy(void) +{ + int i, j; + + /* try various rational and transcendental values for x */ + for (i = -50; i <= 50; i++) + for (j = 1; j < 10; j++) + { + test_remainder_values_y((double) i / (double) j); + test_remainder_values_y(i * M_E + j * M_PI); + } +} + +int main(int argc, char **argv) +{ + fenv_t fenv; + int i; + + printf("Test 47 "); + fflush(stdout); + + /* no FPU errors, please */ + if (feholdexcept(&fenv) < 0) ERR(0, 0); + + /* some signals count as errors */ + for (i = 0; i < _NSIG; i++) + if (i != SIGINT && i != SIGTERM && i != SIGKILL) + signal(i, signal_handler); + + /* test various floating point support functions */ + test_fpclassify_values(); + test_remainder_values_xy(); + test_round_values(); + quit(); + return -1; +}