Extend internal math library from GNU libc to accomodate older systems

with glibc 2.2.X or simply no C99 capable C library. Fix vrfiz instruction
to really truncate on float values.
This commit is contained in:
gbeauche 2005-06-22 12:32:32 +00:00
parent 28a167b3bd
commit f226f793dd
10 changed files with 580 additions and 71 deletions

View File

@ -347,7 +347,8 @@ AC_CHECK_FUNCS(sigaction signal)
AC_CHECK_FUNCS(mmap mprotect munmap)
AC_CHECK_FUNCS(vm_allocate vm_deallocate vm_protect)
AC_CHECK_FUNCS(posix_memalign memalign valloc)
AC_CHECK_FUNCS(exp2f log2f exp2 log2 trunc)
AC_CHECK_FUNCS(exp2f log2f exp2 log2)
AC_CHECK_FUNCS(floorf roundf ceilf truncf floor round ceil trunc)
dnl Darwin seems to define mach_task_self() instead of task_self().
AC_CHECK_FUNCS(mach_task_self task_self)
@ -999,6 +1000,146 @@ EOF
fi
AC_MSG_RESULT($WANT_ADDRESSING_MODE)
dnl Utility macro used by next two tests.
dnl AC_EXAMINE_OBJECT(C source code,
dnl commands examining object file,
dnl [commands to run if compile failed]):
dnl
dnl Compile the source code to an object file; then convert it into a
dnl printable representation. All unprintable characters and
dnl asterisks (*) are replaced by dots (.). All white space is
dnl deleted. Newlines (ASCII 0x10) in the input are preserved in the
dnl output, but runs of newlines are compressed to a single newline.
dnl Finally, line breaks are forcibly inserted so that no line is
dnl longer than 80 columns and the file ends with a newline. The
dnl result of all this processing is in the file conftest.dmp, which
dnl may be examined by the commands in the second argument.
dnl
AC_DEFUN([gcc_AC_EXAMINE_OBJECT],
[AC_LANG_SAVE
AC_LANG_C
dnl Next bit cribbed from AC_TRY_COMPILE.
cat > conftest.$ac_ext <<EOF
[#line __oline__ "configure"
#include "confdefs.h"
$1
]EOF
if AC_TRY_EVAL(ac_compile); then
od -c conftest.o |
sed ['s/^[0-7]*[ ]*/ /
s/\*/./g
s/ \\n/*/g
s/ [0-9][0-9][0-9]/./g
s/ \\[^ ]/./g'] |
tr -d '
' | tr -s '*' '
' | fold | sed '$a\
' > conftest.dmp
$2
ifelse($3, , , else
$3
)dnl
fi
rm -rf conftest*
AC_LANG_RESTORE])
dnl Floating point format probe.
dnl The basic concept is the same as the above: grep the object
dnl file for an interesting string. We have to watch out for
dnl rounding changing the values in the object, however; this is
dnl handled by ignoring the least significant byte of the float.
dnl
dnl Does not know about VAX G-float or C4x idiosyncratic format.
dnl It does know about PDP-10 idiosyncratic format, but this is
dnl not presently supported by GCC. S/390 "binary floating point"
dnl is in fact IEEE (but maybe we should have that in EBCDIC as well
dnl as ASCII?)
dnl
AC_DEFUN([gcc_AC_C_FLOAT_FORMAT],
[AC_CACHE_CHECK(floating point format, ac_cv_c_float_format,
[gcc_AC_EXAMINE_OBJECT(
[/* This will not work unless sizeof(double) == 8. */
extern char sizeof_double_must_be_8 [sizeof(double) == 8 ? 1 : -1];
/* This structure must have no internal padding. */
struct possibility {
char prefix[8];
double candidate;
char postfix[8];
};
#define C(cand) { "\nformat:", cand, ":tamrof\n" }
struct possibility table [] =
{
C( 3.25724264705901305206e+01), /* @@IEEEFP - IEEE 754 */
C( 3.53802595280598432000e+18), /* D__float - VAX */
C( 5.32201830133125317057e-19), /* D.PDP-10 - PDP-10 - the dot is 0x13a */
C( 1.77977764695171661377e+10), /* IBMHEXFP - s/390 format, ascii */
C(-5.22995989424860458374e+10) /* IBMHEXFP - s/390 format, EBCDIC */
};],
[if grep 'format:.@IEEEF.:tamrof' conftest.dmp >/dev/null 2>&1; then
ac_cv_c_float_format='IEEE (big-endian)'
elif grep 'format:.I@@PFE.:tamrof' conftest.dmp >/dev/null 2>&1; then
ac_cv_c_float_format='IEEE (big-endian)'
elif grep 'format:.FEEEI@.:tamrof' conftest.dmp >/dev/null 2>&1; then
ac_cv_c_float_format='IEEE (little-endian)'
elif grep 'format:.EFP@@I.:tamrof' conftest.dmp >/dev/null 2>&1; then
ac_cv_c_float_format='IEEE (little-endian)'
elif grep 'format:.__floa.:tamrof' conftest.dmp >/dev/null 2>&1; then
ac_cv_c_float_format='VAX D-float'
elif grep 'format:..PDP-1.:tamrof' conftest.dmp >/dev/null 2>&1; then
ac_cv_c_float_format='PDP-10'
elif grep 'format:.BMHEXF.:tamrof' conftest.dmp >/dev/null 2>&1; then
ac_cv_c_float_format='IBM 370 hex'
else
AC_MSG_ERROR(Unknown floating point format)
fi],
[AC_MSG_ERROR(compile failed)])
])
# IEEE is the default format. If the float endianness isn't the same
# as the integer endianness, we have to set FLOAT_WORDS_BIG_ENDIAN
# (which is a tristate: yes, no, default). This is only an issue with
# IEEE; the other formats are only supported by a few machines each,
# all with the same endianness.
format=IEEE_FLOAT_FORMAT
fbigend=
case $ac_cv_c_float_format in
'IEEE (big-endian)' )
if test $ac_cv_c_bigendian = no; then
fbigend=1
fi
;;
'IEEE (little-endian)' )
if test $ac_cv_c_bigendian = yes; then
fbigend=0
fi
;;
'VAX D-float' )
format=VAX_FLOAT_FORMAT
;;
'PDP-10' )
format=PDP10_FLOAT_FORMAT
;;
'IBM 370 hex' )
format=IBM_FLOAT_FORMAT
;;
esac
AC_DEFINE_UNQUOTED(HOST_FLOAT_FORMAT, $format,
[Define to the floating point format of the host machine.])
if test -n "$fbigend"; then
AC_DEFINE_UNQUOTED(HOST_FLOAT_WORDS_BIG_ENDIAN, $fbigend,
[Define to 1 if the host machine stores floating point numbers in
memory with the word containing the sign bit at the lowest address,
or to 0 if it does it the other way around.
This macro should not be defined if the ordering is the same as for
multi-word integers.])
fi
])
dnl Check for host float format
gcc_AC_C_FLOAT_FORMAT
dnl Platform specific binary postprocessor
AC_PATH_PROG(BLESS, "true")
if [[ "x$ac_cv_pagezero_hack" = "xyes" ]]; then
@ -1057,6 +1198,7 @@ dnl CPU emulator sources
if [[ "x$EMULATED_PPC" = "xyes" ]]; then
CPUSRCS="\
../kpx_cpu/src/mathlib/ieeefp.cpp \
../kpx_cpu/src/mathlib/mathlib.cpp \
../kpx_cpu/src/cpu/ppc/ppc-cpu.cpp \
../kpx_cpu/src/cpu/ppc/ppc-decode.cpp \
../kpx_cpu/src/cpu/ppc/ppc-execute.cpp \

View File

@ -410,6 +410,14 @@ typedef struct timespec tm_time_t;
typedef struct timeval tm_time_t;
#endif
/* Define codes for all the float formats that we know of.
* Though we only handle IEEE format. */
#define UNKNOWN_FLOAT_FORMAT 0
#define IEEE_FLOAT_FORMAT 1
#define VAX_FLOAT_FORMAT 2
#define IBM_FLOAT_FORMAT 3
#define C4X_FLOAT_FORMAT 4
// High-precision timing
#if defined(HAVE_PTHREADS) && defined(HAVE_CLOCK_NANOSLEEP) && 0
#define PRECISE_TIMING 1

View File

@ -23,6 +23,7 @@
#include <math.h>
#include "mathlib/ieeefp.hpp"
#include "mathlib/mathlib.hpp"
/**
* Define an unary/binary/trinary operation
@ -128,25 +129,6 @@ DEFINE_ALIAS_OP(xor_64, xor, uint64);
// Floating-point basic operations
#ifndef HAVE_EXP2F
#ifdef HAVE_EXP2
#define exp2f(x) (float)exp2(x)
#else
#define exp2f(x) powf(2.0, (x))
#endif
#endif
#ifndef HAVE_LOG2F
#ifdef HAVE_LOG2
#define log2f(x) (float)log2(x)
#else
#ifndef M_LN2
#define M_LN2 logf(2.0)
#endif
#define log2f(x) logf(x) / M_LN2
#endif
#endif
DEFINE_OP1(fnop, double, x);
DEFINE_OP1(fabs, double, fabs(x));
DEFINE_OP2(fadd, double, x + y);
@ -180,7 +162,7 @@ DEFINE_OP1(frsqrt, float, 1 / sqrt(x));
DEFINE_OP1(frim, float, floorf(x));
DEFINE_OP1(frin, float, roundf(x));
DEFINE_OP1(frip, float, ceilf(x));
DEFINE_OP1(friz, float, trunc(x));
DEFINE_OP1(friz, float, truncf(x));
// Misc operations used in AltiVec instructions

View File

@ -1,5 +1,6 @@
/*
* ieeefp-i386.cpp - Access to FPU environment, x86 specific code
* Code largely derived from GNU libc
*
* Kheperix (C) 2003-2005 Gwenole Beauchesne
*
@ -108,21 +109,3 @@ int fesetround(int round)
return 0;
}
// Truncate double value
#ifndef HAVE_TRUNC
#define HAVE_TRUNC
double trunc(double x)
{
volatile unsigned short int cw;
volatile unsigned short int cwtmp;
double value;
__asm__ __volatile__("fnstcw %0" : "=m" (cw));
cwtmp = (cw & 0xf3ff) | 0x0c00; /* toward zero */
__asm__ __volatile__("fldcw %0" : : "m" (cwtmp));
__asm__ __volatile__("frndint" : "=t" (value) : "0" (x));
__asm__ __volatile__("fldcw %0" : : "m" (cw));
return value;
}
#endif

View File

@ -55,23 +55,4 @@ enum {
#endif
// 7.12.14 Comparison macros
#if defined(__GNUC__)
#ifndef isless
#define isless(x, y) \
({ register char __result; \
__asm__ ("fucompp; fnstsw; testb $0x45, %%ah; setz %%al" \
: "=a" (__result) : "u" (x), "t" (y) : "cc", "st", "st(1)"); \
__result; })
#endif
#ifndef isgreater
#define isgreater(x, y) \
({ register char __result; \
__asm__ ("fucompp; fnstsw; testb $0x45, %%ah; setz %%al" \
: "=a" (__result) : "u" (y), "t" (x) : "cc", "st", "st(1)"); \
__result; })
#endif
#endif
#endif /* IEEEFP_I386_H */

View File

@ -36,17 +36,4 @@ extern "C" int fesetround(int);
#endif /* FENV_H */
// Rounding
#ifndef HAVE_TRUNC
extern "C" double trunc(double);
#endif
// Comparison macros
#ifndef isless
#define isless(x, y) ((x) < (y))
#endif
#ifndef isgreater
#define isgreater(x, y) ((x) > (y))
#endif
#endif /* IEEEFP_H */

View File

@ -0,0 +1,38 @@
/*
* mathlib-i386.cpp - Math library wrapper, x86 specific code
* Code largely derived from GNU libc
*
* Kheperix (C) 2003-2005 Gwenole Beauchesne
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU 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 General Public License for more details.
*
* You should have received a copy of the GNU 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
*/
// 7.12.9.8 The trunc functions
#ifndef HAVE_TRUNC
#define HAVE_TRUNC
double trunc(double x)
{
volatile unsigned short int cw;
volatile unsigned short int cwtmp;
double value;
__asm__ __volatile__("fnstcw %0" : "=m" (cw));
cwtmp = (cw & 0xf3ff) | 0x0c00; /* toward zero */
__asm__ __volatile__("fldcw %0" : : "m" (cwtmp));
__asm__ __volatile__("frndint" : "=t" (value) : "0" (x));
__asm__ __volatile__("fldcw %0" : : "m" (cw));
return value;
}
#endif

View File

@ -0,0 +1,44 @@
/*
* mathlib-i386.hpp - Math library wrapper, x86 specific code
* Code largely derived from GNU libc
*
* Kheperix (C) 2003-2005 Gwenole Beauchesne
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU 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 General Public License for more details.
*
* You should have received a copy of the GNU 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
*/
#ifndef MATHLIB_I386_H
#define MATHLIB_I386_H
// 7.12.14 Comparison macros
#if defined(__GNUC__)
#ifndef isless
#define isless(x, y) \
({ register char __result; \
__asm__ ("fucompp; fnstsw; testb $0x45, %%ah; setz %%al" \
: "=a" (__result) : "u" (x), "t" (y) : "cc", "st", "st(1)"); \
__result; })
#endif
#ifndef isgreater
#define isgreater(x, y) \
({ register char __result; \
__asm__ ("fucompp; fnstsw; testb $0x45, %%ah; setz %%al" \
: "=a" (__result) : "u" (y), "t" (x) : "cc", "st", "st(1)"); \
__result; })
#endif
#endif
#endif /* MATHLIB_I386_H */

View File

@ -0,0 +1,182 @@
/*
* mathlib.cpp - Math library wrapper
* Code largely derived from GNU libc
*
* Kheperix (C) 2003-2005 Gwenole Beauchesne
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU 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 General Public License for more details.
*
* You should have received a copy of the GNU 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
*/
#include "sysdeps.h"
#include "mathlib/mathlib.hpp"
#include <stdio.h>
#include <stdlib.h>
#if defined(HOST_FLOAT_WORDS_BIG_ENDIAN) || defined(WORDS_BIGENDIAN)
#define FLOAT_WORD_ORDER_BIG_ENDIAN
#endif
/**
* Representation of an IEEE 754 float
**/
union ieee_float_shape_type {
float value;
uint32 word;
};
// Get a 32 bit int from a float
#define GET_FLOAT_WORD(i,d) \
do { \
ieee_float_shape_type gf_u; \
gf_u.value = (d); \
(i) = gf_u.word; \
} while (0)
/**
* Representation of an IEEE 754 double
**/
union ieee_double_shape_type {
double value;
struct {
#ifdef FLOAT_WORD_ORDER_BIG_ENDIAN
uint32 msw;
uint32 lsw;
#else
uint32 lsw;
uint32 msw;
#endif
} parts;
};
// Get two 32-bit ints from a double
#define EXTRACT_WORDS(ix0,ix1,d) \
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
// Get the more significant 32 bit int from a double
#define GET_HIGH_WORD(i,d) \
do { \
ieee_double_shape_type gh_u; \
gh_u.value = (d); \
(i) = gh_u.parts.msw; \
} while (0)
// Get the less significant 32 bit int from a double
#define GET_LOW_WORD(i,d) \
do { \
ieee_double_shape_type gl_u; \
gl_u.value = (d); \
(i) = gl_u.parts.lsw; \
} while (0)
/**
* Arch-dependent optimizations
**/
#if defined(__i386__)
#include "mathlib/mathlib-i386.cpp"
#endif
/**
* Helper functions
**/
static void unimplemented(const char *function)
{
fprintf(stderr, "MATHLIB: unimplemented function '%s', aborting execution\n", function);
abort();
}
/**
* 7.12.3.1 The fpclassify macro
**/
int mathlib_fpclassifyf (float x)
{
uint32 wx;
int retval = FP_NORMAL;
GET_FLOAT_WORD (wx, x);
wx &= 0x7fffffff;
if (wx == 0)
retval = FP_ZERO;
else if (wx < 0x800000)
retval = FP_SUBNORMAL;
else if (wx >= 0x7f800000)
retval = wx > 0x7f800000 ? FP_NAN : FP_INFINITE;
return retval;
}
int mathlib_fpclassify (double x)
{
uint32 hx, lx;
int retval = FP_NORMAL;
EXTRACT_WORDS (hx, lx, x);
lx |= hx & 0xfffff;
hx &= 0x7ff00000;
if ((hx | lx) == 0)
retval = FP_ZERO;
else if (hx == 0)
retval = FP_SUBNORMAL;
else if (hx == 0x7ff00000)
retval = lx != 0 ? FP_NAN : FP_INFINITE;
return retval;
}
int mathlib_fpclassifyl(long double x)
{
unimplemented("fpclassifyl");
}
/**
* 7.12.3.6 The signbit macro
**/
int mathlib_signbitf (float x)
{
int32 hx;
GET_FLOAT_WORD (hx, x);
return hx & 0x80000000;
}
int mathlib_signbit (double x)
{
int32 hx;
GET_HIGH_WORD (hx, x);
return hx & 0x80000000;
}
int mathlib_signbitl(long double x)
{
unimplemented("signbitl");
}

View File

@ -0,0 +1,162 @@
/*
* mathlib.hpp - Math library wrapper
*
* Kheperix (C) 2003-2005 Gwenole Beauchesne
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU 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 General Public License for more details.
*
* You should have received a copy of the GNU 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
*/
#ifndef MATHLIB_H
#define MATHLIB_H
#include <math.h>
// 7.12 Mathematics <math.h> [#6]
#ifndef FP_NAN
enum {
FP_NAN,
# define FP_NAN FP_NAN
FP_INFINITE,
# define FP_INFINITE FP_INFINITE
FP_ZERO,
# define FP_ZERO FP_ZERO
FP_SUBNORMAL,
# define FP_SUBNORMAL FP_SUBNORMAL
FP_NORMAL
# define FP_NORMAL FP_NORMAL
};
#endif
// Arch-dependent definitions
#if defined(__i386__)
#include "mathlib/mathlib-i386.hpp"
#endif
// 7.12.6.2 The exp2 functions
#ifdef HAVE_EXP2F
extern "C" float exp2f(float x);
#else
#ifdef HAVE_EXP2
extern "C" double exp2(double x);
#define exp2f(x) (float)exp2(x)
#else
#define exp2f(x) powf(2.0, (x))
#endif
#endif
// 7.12.6.10 The log2 functions
#ifdef HAVE_LOG2F
extern "C" float log2f(float x);
#else
#ifdef HAVE_LOG2
extern "C" double log2(double x);
#define log2f(x) (float)log2(x)
#else
#ifndef M_LN2
#define M_LN2 logf(2.0)
#endif
#define log2f(x) logf(x) / M_LN2
#endif
#endif
// 7.12.9.1 The ceil functions
#ifdef HAVE_CEILF
extern "C" float ceilf(float x);
#else
#ifdef HAVE_CEIL
extern "C" double ceil(double x);
#define ceilf(x) (float)ceil(x)
#endif
#endif
// 7.12.9.2 The floor functions
#ifdef HAVE_FLOORF
extern "C" float floorf(float x);
#else
#ifdef HAVE_FLOOR
extern "C" double floor(double x);
#define floorf(x) (float)floor(x)
#endif
#endif
// 7.12.9.6 The round functions
#ifdef HAVE_ROUNDF
extern "C" float roundf(float x);
#else
#ifdef HAVE_ROUND
extern "C" double round(double x);
#define roundf(x) (float)round(x)
#endif
#endif
// 7.12.9.8 The trunc functions
#ifdef HAVE_TRUNCF
extern "C" float truncf(float x);
#else
#ifdef HAVE_TRUNC
extern "C" double trunc(double x);
#define truncf(x) (float)trunc(x)
#endif
#endif
// 7.12.3.1 The fpclassify macro
#ifndef fpclassify
#ifndef mathlib_fpclassifyf
extern int mathlib_fpclassifyf(float x);
#endif
#ifndef mathlib_fpclassify
extern int mathlib_fpclassify(double x);
#endif
#ifndef mathlib_fpclassifyl
extern int mathlib_fpclassifyl(long double x);
#endif
#define fpclassify(x) \
(sizeof (x) == sizeof (float) \
? mathlib_fpclassifyf (x) \
: sizeof (x) == sizeof (double) \
? mathlib_fpclassify (x) : mathlib_fpclassifyl (x))
#endif
// 7.12.3.6 The signbit macro
#ifndef signbit
#ifndef mathlib_signbitf
extern int mathlib_signbitf(float x);
#endif
#ifndef mathlib_signbit
extern int mathlib_signbit(double x);
#endif
#ifndef mathlib_signbitl
extern int mathlib_signbitl(long double x);
#endif
#define signbit(x) \
(sizeof (x) == sizeof (float) \
? mathlib_signbitf (x) \
: sizeof (x) == sizeof (double) \
? mathlib_signbit (x) : mathlib_signbitl (x))
#endif
// 7.12.14.1 The isgreater macro
// FIXME: this is wrong for unordered values
#ifndef isgreater
#define isgreater(x, y) ((x) > (y))
#endif
// 7.12.14.3 The isless macro
// FIXME: this is wrong for unordered values
#ifndef isless
#define isless(x, y) ((x) < (y))
#endif
#endif /* MATHLIB_H */