FSIN,
FSINCOS,
FSQRT,
- FSTD, FSTS, FSTPX, FSTPD, FSTPS,
+ FSTD, FSTS, FSTP, FSTPX, FSTPD, FSTPS,
FSTCW,
FSTENV,
FSTSW,
{ FSTCW, "fnstcw" },
{ FSTD, "fstl" },
{ FSTENV, "fnstenv" },
+ { FSTP, "fstp" },
{ FSTPD, "fstpl" },
{ FSTPS, "fstps" },
{ FSTPX, "fstpt" },
{ "fstcw", FSTCW, WORD },
{ "fstd", FSTD, WORD },
{ "fstenv", FSTENV, WORD },
+ { "fstp", FSTP, WORD },
{ "fstpd", FSTPD, WORD },
{ "fstps", FSTPS, WORD },
{ "fstpx", FSTPX, WORD },
--- /dev/null
+#ifndef _FENV_H
+#define _FENV_H
+
+#include <stdint.h>
+
+#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
#include <ansi.h>
#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) );
#ifdef _POSIX_SOURCE /* STD-C? */
#include <mathconst.h>
+
+#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 */
int sign = 1;
int currexp;
- if (__IsNan(fl)) {
+ if (isnan(fl)) {
errno = EDOM;
return fl;
}
SUBDIRS="\
int64 \
+ math \
misc \
rts \
string"
--- /dev/null
+# 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
--- /dev/null
+#include <assert.h>
+#include <math.h>
+
+#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);
+}
--- /dev/null
+#include <errno.h>
+#include <fenv.h>
+#include <math.h>
+
+#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);
+}
--- /dev/null
+#include <assert.h>
+#include <fenv.h>
+
+#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;
+}
--- /dev/null
+#include <errno.h>
+#include <fenv.h>
+
+#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;
+}
--- /dev/null
+#include <errno.h>
+#include <fenv.h>
+
+#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;
+}
--- /dev/null
+#ifndef __FPU_CW__
+#define __FPU_CW__
+
+#include <stdint.h>
+
+/*
+ * 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__) */
--- /dev/null
+! 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
--- /dev/null
+#ifndef __FPU_RNDINT__
+#define __FPU_RNDINT__
+
+void fpu_rndint(double *value);
+void fpu_remainder(double *x, double y);
+
+#endif /* !defined(__FPU_RNDINT__) */
--- /dev/null
+! 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
--- /dev/null
+#ifndef __FPU_SW__
+#define __FPU_SW__
+
+#include <stdint.h>
+
+#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__) */
+
--- /dev/null
+! 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
+
atan.c \
atan2.c \
ceil.c \
+ compare.c \
exp.c \
fabs.c \
floor.c \
--- /dev/null
+#include <assert.h>
+#include <math.h>
+
+/* 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;
+}
#include <math.h>
double
-__huge_val(void)
+__infinity(void)
{
#if (CHIP == INTEL)
static unsigned char ieee_infinity[] = {
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
+}
+
*/
/* $Header$ */
+#define __IsNan isnan
+
/* some constants (Hart & Cheney) */
#define M_PI 3.14159265358979323846264338327950288
#define M_2PI 6.28318530717958647692528676655900576
#include <math.h>
#include <float.h>
#include <errno.h>
+#include "localmath.h"
#define NITER 5
--- /dev/null
+.TH FEHOLDEXCEPT 3 "December 22, 2009"
+.UC 4
+.SH NAME
+feholdexcept \- disable floating point exceptions
+.SH SYNOPSIS
+.nf
+.ft B
+#include <fenv.h>
+
+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.
+
+
--- /dev/null
+.TH FESETROUND 3 "December 18, 2009"
+.UC 4
+.SH NAME
+fegetround, fesetround \- floating point rounding mode control
+.SH SYNOPSIS
+.nf
+.ft B
+#include <fenv.h>
+
+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)
--- /dev/null
+.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 <math.h>
+
+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.
--- /dev/null
+.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 <math.h>
+
+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.
--- /dev/null
+.TH NEARBYINT 3 "December 18, 2009"
+.UC 4
+.SH NAME
+ceil, floor, nearbyint, trunc \- floating point rounding
+.SH SYNOPSIS
+.nf
+.ft B
+#include <math.h>
+
+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)
--- /dev/null
+.TH REMAINDER 3 "December 18, 2009"
+.UC 4
+.SH NAME
+remainder \- floating point remainder after division
+.SH SYNOPSIS
+.nf
+.ft B
+#include <math.h>
+
+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.
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
test45: test45.c test45.h
test45-gcc: test45.c test45.h
test46: test46.c
+test47: test47.c
# 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
--- /dev/null
+#include <assert.h>
+#include <fenv.h>
+#include <math.h>
+#include <signal.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <sys/sigcontext.h>
+
+#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;
+}