Replace a proliferation of buggy MD implementations of modf() with a

working MI one.  The MI one only needs to be overridden on machines
with non-IEEE754 arithmetic.  (The last supported one was the VAX.)
It can also be overridden if someone comes up with a faster one that
actually passes the regression tests -- but this is harder than it sounds.
This commit is contained in:
das 2011-10-21 06:40:36 +00:00
parent 35a6cd6492
commit 9373f7b9c4
28 changed files with 149 additions and 987 deletions

View File

@ -26,7 +26,6 @@ FBSD_1.0 {
__infinity;
__nan;
makecontext;
modf;
rfork_thread;
setjmp;
longjmp;

View File

@ -2,7 +2,7 @@
# $FreeBSD$
SRCS+= _setjmp.S _set_tp.c rfork_thread.S setjmp.S sigsetjmp.S \
fabs.S modf.S \
fabs.S \
infinity.c ldexp.c makecontext.c signalcontext.c \
flt_rounds.c fpgetmask.c fpsetmask.c fpgetprec.c fpsetprec.c \
fpgetround.c fpsetround.c fpgetsticky.c

View File

@ -1,91 +0,0 @@
/*-
* Copyright (c) 1990 The Regents of the University of California.
* All rights reserved.
*
* This code is derived from software contributed to Berkeley by
* Sean Eric Fagan.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 4. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* from: @(#)modf.s 5.5 (Berkeley) 3/18/91
*/
#include <machine/asm.h>
#if defined(LIBC_SCCS)
RCSID("$NetBSD: modf.S,v 1.5 1997/07/16 14:37:18 christos Exp $")
#endif
__FBSDID("$FreeBSD$");
/*
* modf(value, iptr): return fractional part of value, and stores the
* integral part into iptr (a pointer to double).
*
* Written by Sean Eric Fagan (sef@kithrup.COM)
* Sun Mar 11 20:27:30 PST 1990
*/
/* With CHOP mode on, frndint behaves as TRUNC does. Useful. */
ENTRY(modf)
/*
* Set chop mode.
*/
fnstcw -12(%rsp)
movw -12(%rsp),%dx
orw $3072,%dx
movw %dx,-16(%rsp)
fldcw -16(%rsp)
/*
* Get integral part.
*/
movsd %xmm0,-24(%rsp)
fldl -24(%rsp)
frndint
fstpl -8(%rsp)
/*
* Restore control word.
*/
fldcw -12(%rsp)
/*
* Store integral part.
*/
movsd -8(%rsp),%xmm0
movsd %xmm0,(%rdi)
/*
* Get fractional part and return it.
*/
fldl -24(%rsp)
fsubl -8(%rsp)
fstpl -8(%rsp)
movsd -8(%rsp),%xmm0
ret
END(modf)
.section .note.GNU-stack,"",%progbits

View File

@ -19,7 +19,6 @@ FBSD_1.0 {
__infinity;
__nan;
makecontext;
modf;
setjmp;
longjmp;
sigsetjmp;

View File

@ -2,5 +2,5 @@
# $FreeBSD$
SRCS+= _ctx_start.S _setjmp.S _set_tp.c alloca.S fabs.c \
infinity.c ldexp.c makecontext.c modf.c \
infinity.c ldexp.c makecontext.c \
setjmp.S signalcontext.c sigsetjmp.S divsi3.S

View File

@ -1,107 +0,0 @@
/*
* Copyright (c) 1994, 1995 Carnegie-Mellon University.
* All rights reserved.
*
* Author: Chris G. Demetriou
*
* Permission to use, copy, modify and distribute this software and
* its documentation is hereby granted, provided that both the copyright
* notice and this permission notice appear in all copies of the
* software, derivative works or modified versions, and any portions
* thereof, and that both notices appear in supporting documentation.
*
* CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS"
* CONDITION. CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND
* FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
*
* Carnegie Mellon requests users of this software to return to
*
* Software Distribution Coordinator or Software.Distribution@CS.CMU.EDU
* School of Computer Science
* Carnegie Mellon University
* Pittsburgh PA 15213-3890
*
* any improvements or extensions that they make and grant Carnegie the
* rights to redistribute these changes.
*
* $NetBSD: modf.c,v 1.1 1995/02/10 17:50:25 cgd Exp $
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <sys/types.h>
#include <errno.h>
#include <math.h>
#include <machine/ieee.h>
/*
* double modf(double val, double *iptr)
* returns: f and i such that |f| < 1.0, (f + i) = val, and
* sign(f) == sign(i) == sign(val).
*
* Beware signedness when doing subtraction, and also operand size!
*/
double
modf(val, iptr)
double val, *iptr;
{
union doub {
double v;
struct ieee_double s;
} u, v;
u_int64_t frac;
/*
* If input is Inf or NaN, return it and leave i alone.
*/
u.v = val;
if (u.s.dbl_exp == DBL_EXP_INFNAN)
return (u.v);
/*
* If input can't have a fractional part, return
* (appropriately signed) zero, and make i be the input.
*/
if ((int)u.s.dbl_exp - DBL_EXP_BIAS > DBL_FRACBITS - 1) {
*iptr = u.v;
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
return (v.v);
}
/*
* If |input| < 1.0, return it, and set i to the appropriately
* signed zero.
*/
if (u.s.dbl_exp < DBL_EXP_BIAS) {
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
*iptr = v.v;
return (u.v);
}
/*
* There can be a fractional part of the input.
* If you look at the math involved for a few seconds, it's
* plain to see that the integral part is the input, with the
* low (DBL_FRACBITS - (exponent - DBL_EXP_BIAS)) bits zeroed,
* the fractional part is the part with the rest of the
* bits zeroed. Just zeroing the high bits to get the
* fractional part would yield a fraction in need of
* normalization. Therefore, we take the easy way out, and
* just use subtraction to get the fractional part.
*/
v.v = u.v;
/* Zero the low bits of the fraction, the sleazy way. */
frac = ((u_int64_t)v.s.dbl_frach << 32) + v.s.dbl_fracl;
frac >>= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
frac <<= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
v.s.dbl_fracl = frac & 0xffffffff;
v.s.dbl_frach = frac >> 32;
*iptr = v.v;
u.v -= v.v;
u.s.dbl_sign = v.s.dbl_sign;
return (u.v);
}

View File

@ -35,6 +35,8 @@ SRCS+= __getosreldate.c __xuname.c \
usleep.c utime.c utxdb.c valloc.c vis.c wait.c wait3.c waitpid.c \
wordexp.c
MISRCS+=modf.c
CANCELPOINTS_SRCS=sem.c sem_new.c
.for src in ${CANCELPOINTS_SRCS}
SRCS+=cancelpoints_${src}

View File

@ -213,6 +213,7 @@ FBSD_1.0 {
ldexp;
lockf;
lrand48;
modf;
mrand48;
nftw;
nice;

138
lib/libc/gen/modf.c Normal file
View File

@ -0,0 +1,138 @@
/* @(#)s_modf.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* modf(double x, double *iptr)
* return fraction part of x, and return x's integral part in *iptr.
* Method:
* Bit twiddling.
*
* Exception:
* No exception.
*/
#include <sys/types.h>
#include <machine/endian.h>
#include <math.h>
/* Bit fiddling routines copied from msun/src/math_private.h,v 1.15 */
#if BYTE_ORDER == BIG_ENDIAN
typedef union
{
double value;
struct
{
u_int32_t msw;
u_int32_t lsw;
} parts;
} ieee_double_shape_type;
#endif
#if BYTE_ORDER == LITTLE_ENDIAN
typedef union
{
double value;
struct
{
u_int32_t lsw;
u_int32_t msw;
} parts;
} ieee_double_shape_type;
#endif
/* 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)
/* Set a double from two 32 bit ints. */
#define INSERT_WORDS(d,ix0,ix1) \
do { \
ieee_double_shape_type iw_u; \
iw_u.parts.msw = (ix0); \
iw_u.parts.lsw = (ix1); \
(d) = iw_u.value; \
} while (0)
static const double one = 1.0;
double
modf(double x, double *iptr)
{
int32_t i0,i1,j0;
u_int32_t i;
EXTRACT_WORDS(i0,i1,x);
j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */
if(j0<20) { /* integer part in high x */
if(j0<0) { /* |x|<1 */
INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */
return x;
} else {
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) { /* x is integral */
u_int32_t high;
*iptr = x;
GET_HIGH_WORD(high,x);
INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
return x;
} else {
INSERT_WORDS(*iptr,i0&(~i),0);
return x - *iptr;
}
}
} else if (j0>51) { /* no fraction part */
u_int32_t high;
if (j0 == 0x400) { /* inf/NaN */
*iptr = x;
return 0.0 / x;
}
*iptr = x*one;
GET_HIGH_WORD(high,x);
INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
return x;
} else { /* fraction part in low x */
i = ((u_int32_t)(0xffffffff))>>(j0-20);
if((i1&i)==0) { /* x is integral */
u_int32_t high;
*iptr = x;
GET_HIGH_WORD(high,x);
INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
return x;
} else {
INSERT_WORDS(*iptr,i0,i1&(~i));
return x - *iptr;
}
}
}

View File

@ -20,7 +20,6 @@ FBSD_1.0 {
__nan;
__infinity;
makecontext;
modf;
rfork_thread;
setjmp;
longjmp;

View File

@ -2,5 +2,5 @@
# $FreeBSD$
SRCS+= _ctx_start.S _setjmp.S _set_tp.c fabs.S \
flt_rounds.c infinity.c ldexp.c makecontext.c modf.S \
flt_rounds.c infinity.c ldexp.c makecontext.c \
rfork_thread.S setjmp.S signalcontext.c sigsetjmp.S

View File

@ -1,87 +0,0 @@
/*-
* Copyright (c) 1990 The Regents of the University of California.
* All rights reserved.
*
* This code is derived from software contributed to Berkeley by
* Sean Eric Fagan.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 4. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#if defined(LIBC_SCCS) && !defined(lint)
.asciz "@(#)modf.s 5.5 (Berkeley) 3/18/91"
#endif /* LIBC_SCCS and not lint */
#include <machine/asm.h>
__FBSDID("$FreeBSD$");
/*
* modf(value, iptr): return fractional part of value, and stores the
* integral part into iptr (a pointer to double).
*
* Written by Sean Eric Fagan (sef@kithrup.COM)
* Sun Mar 11 20:27:30 PST 1990
*/
/* With CHOP mode on, frndint behaves as TRUNC does. Useful. */
ENTRY(modf)
pushl %ebp
movl %esp,%ebp
/* Check for Inf/NaN */
movl 12(%ebp),%eax
andl $0x7fffffff,%eax
cmpl $0x7ff00000,%eax
jae 1f
/* Finite value */
subl $16,%esp
fnstcw -12(%ebp)
movw -12(%ebp),%dx
orw $3072,%dx
movw %dx,-16(%ebp)
fldcw -16(%ebp)
fldl 8(%ebp)
frndint
fstpl -8(%ebp)
fldcw -12(%ebp)
movl 16(%ebp),%eax
movl -8(%ebp),%edx
movl -4(%ebp),%ecx
movl %edx,(%eax)
movl %ecx,4(%eax)
fldl 8(%ebp)
fsubl -8(%ebp)
leave
ret
/* Inf/NaN handling */
1: fldl 8(%ebp)
movl 16(%ebp),%edx
fstl (%edx)
fldz
fdivp /* return +/- 0 for +/- Inf, NaN for NaN */
leave
ret
END(modf)
.section .note.GNU-stack,"",%progbits

View File

@ -23,7 +23,6 @@ FBSD_1.0 {
__infinity;
__nan;
makecontext;
modf;
setjmp;
longjmp;
sigsetjmp;

View File

@ -3,7 +3,7 @@
SRCS+= __divdf3.S __divdi3.S __divsf3.S __divsi3.S __moddi3.S __modsi3.S \
__udivdi3.S __udivsi3.S __umoddi3.S __umodsi3.S _mcount.S _set_tp.c \
_setjmp.S fabs.S flt_rounds.c fpgetmask.c fpgetround.c fpsetmask.c \
fpsetround.c infinity.c ldexp.c makecontext.c modf.c setjmp.S \
fpsetround.c infinity.c ldexp.c makecontext.c setjmp.S \
signalcontext.c sigsetjmp.S
# The following may go away if function _Unwind_FindTableEntry()

View File

@ -1,106 +0,0 @@
/* $NetBSD: modf.c,v 1.1 1995/02/10 17:50:25 cgd Exp $ */
/*
* Copyright (c) 1994, 1995 Carnegie-Mellon University.
* All rights reserved.
*
* Author: Chris G. Demetriou
*
* Permission to use, copy, modify and distribute this software and
* its documentation is hereby granted, provided that both the copyright
* notice and this permission notice appear in all copies of the
* software, derivative works or modified versions, and any portions
* thereof, and that both notices appear in supporting documentation.
*
* CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS"
* CONDITION. CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND
* FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
*
* Carnegie Mellon requests users of this software to return to
*
* Software Distribution Coordinator or Software.Distribution@CS.CMU.EDU
* School of Computer Science
* Carnegie Mellon University
* Pittsburgh PA 15213-3890
*
* any improvements or extensions that they make and grant Carnegie the
* rights to redistribute these changes.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <sys/types.h>
#include <machine/ieee.h>
#include <math.h>
/*
* double modf(double val, double *iptr)
* returns: f and i such that |f| < 1.0, (f + i) = val, and
* sign(f) == sign(i) == sign(val).
*
* Beware signedness when doing subtraction, and also operand size!
*/
double
modf(val, iptr)
double val, *iptr;
{
union doub {
double v;
struct ieee_double s;
} u, v;
u_int64_t frac;
/*
* If input is Inf or NaN, return it and leave i alone.
*/
u.v = val;
if (u.s.dbl_exp == DBL_EXP_INFNAN)
return (u.v);
/*
* If input can't have a fractional part, return
* (appropriately signed) zero, and make i be the input.
*/
if ((int)u.s.dbl_exp - DBL_EXP_BIAS > DBL_FRACBITS - 1) {
*iptr = u.v;
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
return (v.v);
}
/*
* If |input| < 1.0, return it, and set i to the appropriately
* signed zero.
*/
if (u.s.dbl_exp < DBL_EXP_BIAS) {
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
*iptr = v.v;
return (u.v);
}
/*
* There can be a fractional part of the input.
* If you look at the math involved for a few seconds, it's
* plain to see that the integral part is the input, with the
* low (DBL_FRACBITS - (exponent - DBL_EXP_BIAS)) bits zeroed,
* the fractional part is the part with the rest of the
* bits zeroed. Just zeroing the high bits to get the
* fractional part would yield a fraction in need of
* normalization. Therefore, we take the easy way out, and
* just use subtraction to get the fractional part.
*/
v.v = u.v;
/* Zero the low bits of the fraction, the sleazy way. */
frac = ((u_int64_t)v.s.dbl_frach << 32) + v.s.dbl_fracl;
frac >>= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
frac <<= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
v.s.dbl_fracl = frac & 0xffffffff;
v.s.dbl_frach = frac >> 32;
*iptr = v.v;
u.v -= v.v;
u.s.dbl_sign = v.s.dbl_sign;
return (u.v);
}

View File

@ -18,7 +18,6 @@ FBSD_1.0 {
__infinity;
__nan;
makecontext;
modf;
setjmp;
longjmp;
sigsetjmp;

View File

@ -1,7 +1,7 @@
# $NetBSD: Makefile.inc,v 1.27 2005/10/07 17:16:40 tsutsui Exp $
# $FreeBSD$
SRCS+= infinity.c fabs.c ldexp.c modf.c
SRCS+= infinity.c fabs.c ldexp.c
# SRCS+= flt_rounds.c fpgetmask.c fpgetround.c fpgetsticky.c fpsetmask.c \
# fpsetround.c fpsetsticky.c

View File

@ -1,82 +0,0 @@
/* $NetBSD: modf.S,v 1.10 2003/08/07 16:42:15 agc Exp $ */
/*-
* Copyright (c) 1991, 1993, 1995
* The Regents of the University of California. All rights reserved.
*
* This code is derived from software contributed to Berkeley by
* Ralph Campbell.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <machine/asm.h>
__FBSDID("$FreeBSD$");
#if defined(LIBC_SCCS) && !defined(lint)
ASMSTR("from: @(#)modf.s 8.1 (Berkeley) 6/4/93")
ASMSTR("$NetBSD: modf.S,v 1.10 2003/08/07 16:42:15 agc Exp $")
#endif /* LIBC_SCCS and not lint */
#ifdef __ABICALLS__
.abicalls
#endif
/*
* double modf(val, iptr)
* double val, *iptr;
* returns: xxx and n (in *iptr) where val == n.xxx
*/
LEAF(modf)
#ifdef __ABICALLS__
.set noreorder
.cpload t9
.set reorder
#endif
cfc1 t0, $31 # get the control register
li.d $f2, 4503599627370496e0 # f2 <- 2^52
or t1, t0, 0x3 # set rounding mode to round to zero
xor t1, t1, 0x2 # (i.e., 01)
ctc1 t1, $31
mov.d $f0, $f12 # f0 <- f12
abs.d $f4, $f12 # f4 <- |f12|
c.olt.d $f4, $f2 # f4 ? < f2
bc1f 1f # leave f0 alone if Nan, infinity
# or >=2^52
c.eq.d $f12,$f4 # was f12 positive ?
add.d $f4,$f2,$f4 # round off to integer
bc1f 2f # No -> will have to negate result
sub.d $f0,$f4,$f2 # Remove fudge factor
j 1f # integer fraction got
2:
sub.d $f0,$f2,$f4 # Remove fudge factor and negate
1:
ctc1 t0, $31 # restore old rounding mode
s.d $f0, 0(a2) # save the integer part
sub.d $f0, $f12, $f0 # subtract val - integer part
j ra
END(modf)

View File

@ -1,107 +0,0 @@
/*
* Copyright (c) 1994, 1995 Carnegie-Mellon University.
* All rights reserved.
*
* Author: Chris G. Demetriou
*
* Permission to use, copy, modify and distribute this software and
* its documentation is hereby granted, provided that both the copyright
* notice and this permission notice appear in all copies of the
* software, derivative works or modified versions, and any portions
* thereof, and that both notices appear in supporting documentation.
*
* CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS"
* CONDITION. CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND
* FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
*
* Carnegie Mellon requests users of this software to return to
*
* Software Distribution Coordinator or Software.Distribution@CS.CMU.EDU
* School of Computer Science
* Carnegie Mellon University
* Pittsburgh PA 15213-3890
*
* any improvements or extensions that they make and grant Carnegie the
* rights to redistribute these changes.
*
* $NetBSD: modf.c,v 1.1 1995/02/10 17:50:25 cgd Exp $
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <sys/types.h>
#include <errno.h>
#include <math.h>
#include <machine/ieee.h>
/*
* double modf(double val, double *iptr)
* returns: f and i such that |f| < 1.0, (f + i) = val, and
* sign(f) == sign(i) == sign(val).
*
* Beware signedness when doing subtraction, and also operand size!
*/
double
modf(val, iptr)
double val, *iptr;
{
union doub {
double v;
struct ieee_double s;
} u, v;
u_int64_t frac;
/*
* If input is Inf or NaN, return it and leave i alone.
*/
u.v = val;
if (u.s.dbl_exp == DBL_EXP_INFNAN)
return (u.v);
/*
* If input can't have a fractional part, return
* (appropriately signed) zero, and make i be the input.
*/
if ((int)u.s.dbl_exp - DBL_EXP_BIAS > DBL_FRACBITS - 1) {
*iptr = u.v;
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
return (v.v);
}
/*
* If |input| < 1.0, return it, and set i to the appropriately
* signed zero.
*/
if (u.s.dbl_exp < DBL_EXP_BIAS) {
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
*iptr = v.v;
return (u.v);
}
/*
* There can be a fractional part of the input.
* If you look at the math involved for a few seconds, it's
* plain to see that the integral part is the input, with the
* low (DBL_FRACBITS - (exponent - DBL_EXP_BIAS)) bits zeroed,
* the fractional part is the part with the rest of the
* bits zeroed. Just zeroing the high bits to get the
* fractional part would yield a fraction in need of
* normalization. Therefore, we take the easy way out, and
* just use subtraction to get the fractional part.
*/
v.v = u.v;
/* Zero the low bits of the fraction, the sleazy way. */
frac = ((u_int64_t)v.s.dbl_frach << 32) + v.s.dbl_fracl;
frac >>= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
frac <<= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
v.s.dbl_fracl = frac & 0xffffffff;
v.s.dbl_frach = frac >> 32;
*iptr = v.v;
u.v -= v.v;
u.s.dbl_sign = v.s.dbl_sign;
return (u.v);
}

View File

@ -24,7 +24,6 @@ FBSD_1.0 {
__infinity;
__nan;
makecontext;
modf;
setjmp;
longjmp;
sigsetjmp;

View File

@ -2,7 +2,7 @@
SRCS += _ctx_start.S fabs.S flt_rounds.c fpgetmask.c fpgetround.c \
fpgetsticky.c fpsetmask.c fpsetround.c \
infinity.c ldexp.c makecontext.c modf.c _setjmp.S \
infinity.c ldexp.c makecontext.c _setjmp.S \
setjmp.S sigsetjmp.S signalcontext.c syncicache.c \
_set_tp.c

View File

@ -1,107 +0,0 @@
/*
* Copyright (c) 1994, 1995 Carnegie-Mellon University.
* All rights reserved.
*
* Author: Chris G. Demetriou
*
* Permission to use, copy, modify and distribute this software and
* its documentation is hereby granted, provided that both the copyright
* notice and this permission notice appear in all copies of the
* software, derivative works or modified versions, and any portions
* thereof, and that both notices appear in supporting documentation.
*
* CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS"
* CONDITION. CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND
* FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
*
* Carnegie Mellon requests users of this software to return to
*
* Software Distribution Coordinator or Software.Distribution@CS.CMU.EDU
* School of Computer Science
* Carnegie Mellon University
* Pittsburgh PA 15213-3890
*
* any improvements or extensions that they make and grant Carnegie the
* rights to redistribute these changes.
*
* $NetBSD: modf.c,v 1.1 1995/02/10 17:50:25 cgd Exp $
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <sys/types.h>
#include <machine/ieee.h>
#include <errno.h>
#include <math.h>
/*
* double modf(double val, double *iptr)
* returns: f and i such that |f| < 1.0, (f + i) = val, and
* sign(f) == sign(i) == sign(val).
*
* Beware signedness when doing subtraction, and also operand size!
*/
double
modf(val, iptr)
double val, *iptr;
{
union doub {
double v;
struct ieee_double s;
} u, v;
u_int64_t frac;
/*
* If input is Inf or NaN, return it and leave i alone.
*/
u.v = val;
if (u.s.dbl_exp == DBL_EXP_INFNAN)
return (u.v);
/*
* If input can't have a fractional part, return
* (appropriately signed) zero, and make i be the input.
*/
if ((int)u.s.dbl_exp - DBL_EXP_BIAS > DBL_FRACBITS - 1) {
*iptr = u.v;
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
return (v.v);
}
/*
* If |input| < 1.0, return it, and set i to the appropriately
* signed zero.
*/
if (u.s.dbl_exp < DBL_EXP_BIAS) {
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
*iptr = v.v;
return (u.v);
}
/*
* There can be a fractional part of the input.
* If you look at the math involved for a few seconds, it's
* plain to see that the integral part is the input, with the
* low (DBL_FRACBITS - (exponent - DBL_EXP_BIAS)) bits zeroed,
* the fractional part is the part with the rest of the
* bits zeroed. Just zeroing the high bits to get the
* fractional part would yield a fraction in need of
* normalization. Therefore, we take the easy way out, and
* just use subtraction to get the fractional part.
*/
v.v = u.v;
/* Zero the low bits of the fraction, the sleazy way. */
frac = ((u_int64_t)v.s.dbl_frach << 32) + v.s.dbl_fracl;
frac >>= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
frac <<= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
v.s.dbl_fracl = frac & 0xffffffff;
v.s.dbl_frach = frac >> 32;
*iptr = v.v;
u.v -= v.v;
u.s.dbl_sign = v.s.dbl_sign;
return (u.v);
}

View File

@ -24,7 +24,6 @@ FBSD_1.0 {
__infinity;
__nan;
makecontext;
modf;
setjmp;
longjmp;
sigsetjmp;

View File

@ -2,7 +2,7 @@
SRCS += _ctx_start.S fabs.S flt_rounds.c fpgetmask.c fpgetround.c \
fpgetsticky.c fpsetmask.c fpsetround.c \
infinity.c ldexp.c makecontext.c modf.c _setjmp.S \
infinity.c ldexp.c makecontext.c _setjmp.S \
setjmp.S sigsetjmp.S signalcontext.c syncicache.c \
_set_tp.c

View File

@ -1,107 +0,0 @@
/*
* Copyright (c) 1994, 1995 Carnegie-Mellon University.
* All rights reserved.
*
* Author: Chris G. Demetriou
*
* Permission to use, copy, modify and distribute this software and
* its documentation is hereby granted, provided that both the copyright
* notice and this permission notice appear in all copies of the
* software, derivative works or modified versions, and any portions
* thereof, and that both notices appear in supporting documentation.
*
* CARNEGIE MELLON ALLOWS FREE USE OF THIS SOFTWARE IN ITS "AS IS"
* CONDITION. CARNEGIE MELLON DISCLAIMS ANY LIABILITY OF ANY KIND
* FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE USE OF THIS SOFTWARE.
*
* Carnegie Mellon requests users of this software to return to
*
* Software Distribution Coordinator or Software.Distribution@CS.CMU.EDU
* School of Computer Science
* Carnegie Mellon University
* Pittsburgh PA 15213-3890
*
* any improvements or extensions that they make and grant Carnegie the
* rights to redistribute these changes.
*
* $NetBSD: modf.c,v 1.1 1995/02/10 17:50:25 cgd Exp $
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <sys/types.h>
#include <machine/ieee.h>
#include <errno.h>
#include <math.h>
/*
* double modf(double val, double *iptr)
* returns: f and i such that |f| < 1.0, (f + i) = val, and
* sign(f) == sign(i) == sign(val).
*
* Beware signedness when doing subtraction, and also operand size!
*/
double
modf(val, iptr)
double val, *iptr;
{
union doub {
double v;
struct ieee_double s;
} u, v;
u_int64_t frac;
/*
* If input is Inf or NaN, return it and leave i alone.
*/
u.v = val;
if (u.s.dbl_exp == DBL_EXP_INFNAN)
return (u.v);
/*
* If input can't have a fractional part, return
* (appropriately signed) zero, and make i be the input.
*/
if ((int)u.s.dbl_exp - DBL_EXP_BIAS > DBL_FRACBITS - 1) {
*iptr = u.v;
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
return (v.v);
}
/*
* If |input| < 1.0, return it, and set i to the appropriately
* signed zero.
*/
if (u.s.dbl_exp < DBL_EXP_BIAS) {
v.v = 0.0;
v.s.dbl_sign = u.s.dbl_sign;
*iptr = v.v;
return (u.v);
}
/*
* There can be a fractional part of the input.
* If you look at the math involved for a few seconds, it's
* plain to see that the integral part is the input, with the
* low (DBL_FRACBITS - (exponent - DBL_EXP_BIAS)) bits zeroed,
* the fractional part is the part with the rest of the
* bits zeroed. Just zeroing the high bits to get the
* fractional part would yield a fraction in need of
* normalization. Therefore, we take the easy way out, and
* just use subtraction to get the fractional part.
*/
v.v = u.v;
/* Zero the low bits of the fraction, the sleazy way. */
frac = ((u_int64_t)v.s.dbl_frach << 32) + v.s.dbl_fracl;
frac >>= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
frac <<= DBL_FRACBITS - (u.s.dbl_exp - DBL_EXP_BIAS);
v.s.dbl_fracl = frac & 0xffffffff;
v.s.dbl_frach = frac >> 32;
*iptr = v.v;
u.v -= v.v;
u.s.dbl_sign = v.s.dbl_sign;
return (u.v);
}

View File

@ -24,7 +24,6 @@ FBSD_1.0 {
__infinity;
__nan;
makecontext;
modf;
setjmp;
longjmp;
sigsetjmp;

View File

@ -2,5 +2,5 @@
SRCS+= _ctx_start.S _setjmp.S fabs.S fixunsdfsi.S flt_rounds.c fpgetmask.c \
fpgetround.c fpgetsticky.c fpsetmask.c fpsetround.c \
infinity.c ldexp.c makecontext.c modf.S \
infinity.c ldexp.c makecontext.c \
signalcontext.c setjmp.S sigsetjmp.S _set_tp.c

View File

@ -1,177 +0,0 @@
/*
* Copyright (c) 1992, 1993
* The Regents of the University of California. All rights reserved.
*
* This software was developed by the Computer Systems Engineering group
* at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
* contributed to Berkeley.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 4. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* from: Header: modf.s,v 1.3 92/06/20 00:00:54 torek Exp
*/
#if defined(LIBC_SCCS) && !defined(lint)
.asciz "@(#)modf.s 8.1 (Berkeley) 6/4/93"
#if 0
.asciz "$NetBSD: modf.S,v 1.2 2000/07/23 07:12:22 eeh Exp $"
#endif
#endif /* LIBC_SCCS and not lint */
#include <machine/asm.h>
__FBSDID("$FreeBSD$");
#include "assym.s"
/*
* double modf(double val, double *iptr)
*
* Returns the fractional part of `val', storing the integer part of
* `val' in *iptr. Both *iptr and the return value have the same sign
* as `val'.
*
* Method:
*
* We use the fpu's normalization hardware to compute the integer portion
* of the double precision argument. Sun IEEE double precision numbers
* have 52 bits of mantissa, 11 bits of exponent, and one bit of sign,
* with the sign occupying bit 31 of word 0, and the exponent bits 30:20
* of word 0. Thus, values >= 2^52 are by definition integers.
*
* If we take a value that is in the range [+0..2^52) and add 2^52, all
* of the fractional bits fall out and all of the integer bits are summed
* with 2^52. If we then subtract 2^52, we get those integer bits back.
* This must be done with rounding set to `towards 0' or `towards -inf'.
* `Toward -inf' fails when the value is 0 (we get -0 back)....
*
* Note that this method will work anywhere, but is machine dependent in
* various aspects.
*
* Stack usage:
* 4@[%fp + SPOFF - 4] saved %fsr
* 4@[%fp + SPOFF - 8] new %fsr with rounding set to `towards 0'
* 8@[%fp + SPOFF - 16] space for moving between %i and %f registers
* Register usage:
* %f0:f1 double val;
* %l0 scratch
* %l1 sign bit (0x80000000)
* %i1 double *iptr;
* %f2:f3 `magic number' 2^52, in fpu registers
* %f4:f5 double v, in fpu registers
* %f6:f7 double temp.
*/
.align 8
.Lmagic:
.word 0x43300000 ! sign = 0, exponent = 52 + 1023, mantissa = 0
.word 0 ! (i.e., .double 0r4503599627370496e+00)
.L0:
.word 0 ! 0.0
.word 0
ENTRY(modf)
save %sp, -CCFSZ - 16, %sp
PIC_PROLOGUE(%l6, %l7)
/*
* First, compute v = abs(val)
*/
fabsd %f0, %f4 ! %f4:f5 = v
fcmped %fcc1, %f0, %f4 ! %fcc1 = (val == abs(val))
SET(.Lmagic, %l7, %l0)
ldd [%l0], %f2
/*
* Is %f4:f5 >= %f2:f3 ? If so, it is all integer bits.
* It is probably less, though.
*/
fcmped %f4, %f2
fbuge .Lbig ! if >= (or unordered), go out
nop
/*
* v < 2^52, so add 2^52, then subtract 2^52, but do it all
* with rounding set towards zero. We leave any enabled
* traps enabled, but change the rounding mode. This might
* not be so good. Oh well....
*/
st %fsr, [%fp + SPOFF - 4] ! %l5 = current FSR mode
set FSR_RD_MASK, %l3 ! %l3 = rounding direction mask
ld [%fp + SPOFF - 4], %l5
set FSR_RD_RD_Z, %l4
andn %l5, %l3, %l6
or %l6, %l4, %l6 ! round towards zero, please
and %l5, %l3, %l5 ! save original rounding mode
st %l6, [%fp + SPOFF - 8]
ld [%fp + SPOFF - 8], %fsr
faddd %f4, %f2, %f4 ! %f4:f5 += 2^52
fsubd %f4, %f2, %f4 ! %f4:f5 -= 2^52
/*
* Restore %fsr, but leave exceptions accrued.
*/
st %fsr, [%fp + SPOFF - 4]
ld [%fp + SPOFF - 4], %l6
andn %l6, %l3, %l6 ! %l6 = %fsr & ~FSR_RD_MASK;
or %l5, %l6, %l5 ! %l5 |= %l6;
st %l5, [%fp + SPOFF - 4]
ld [%fp + SPOFF - 4], %fsr ! restore %fsr, leaving accrued stuff
/*
* Now insert the original sign in %f4:f5.
* %fcc1 should still have the results of (val == abs(val))
* from above, so we use a conditional move on %fcc1 to:
*
* %f4 = (val == abs(val)) ? %f4 : -%f4
*
*/
fnegd %f4, %f6
fmovdnz %fcc1, %f6, %f4
1:
/*
* The value in %f4:f5 is now the integer portion of the original
* argument. We need to store this in *ival (%i1), subtract it
* from the original value argument (%d0), and return the result.
*/
std %f4, [%i1] ! *ival = %f4:f5;
fsubd %f0, %f4, %f0 ! %f0:f1 -= %f4:f5;
ret
restore
.Lbig:
/*
* We get here if the original comparison of %f4:f5 (v) to
* %f2:f3 (2^52) came out `greater or unordered'. In this
* case the integer part is the original value, and the
* fractional part is 0.
*/
SET(.L0, %l7, %l0)
std %f0, [%i1] ! *ival = val;
ldd [%l0], %f0 ! return 0.0;
ret
restore
END(modf)