hook fma to the build

This commit is contained in:
christos 2013-02-11 01:29:58 +00:00
parent 3a4e9cb991
commit e3f9b211a9
5 changed files with 72 additions and 22 deletions

View File

@ -1,4 +1,4 @@
# $NetBSD: Makefile,v 1.128 2013/02/09 22:33:13 christos Exp $
# $NetBSD: Makefile,v 1.129 2013/02/11 01:29:58 christos Exp $
#
# @(#)Makefile 5.1beta 93/09/24
#
@ -67,7 +67,7 @@ COMMON_SRCS+= fenv.c
.PATH: ${.CURDIR}/arch/i387
COMMON_SRCS+= fenv.c s_nextafterl.c s_nexttoward.c s_nexttowardf.c \
s_nearbyint.c s_rintl.c # s_nan.c
s_nearbyint.c s_rintl.c s_fma.c s_fmaf.c s_fmal.c # s_nan.c
ARCH_SRCS = e_acos.S e_asin.S e_atan2.S e_exp.S e_expf.S e_fmod.S e_log.S \
e_logf.S e_log10.S e_log10f.S e_log2.S e_log2f.S e_remainder.S \
e_remainderf.S e_scalb.S e_scalbf.S e_sqrt.S e_sqrtf.S s_atan.S \

View File

@ -11,7 +11,7 @@
/*
* from: @(#)fdlibm.h 5.1 93/09/24
* $NetBSD: math_private.h,v 1.17 2012/05/05 17:54:14 christos Exp $
* $NetBSD: math_private.h,v 1.18 2013/02/11 01:29:58 christos Exp $
*/
#ifndef _MATH_PRIVATE_H_
@ -48,6 +48,9 @@ typedef union
u_int32_t msw;
u_int32_t lsw;
} parts;
struct {
u_int64_t w;
} xparts;
} ieee_double_shape_type;
#endif
@ -63,6 +66,9 @@ typedef union
u_int32_t lsw;
u_int32_t msw;
} parts;
struct {
u_int64_t w;
} xparts;
} ieee_double_shape_type;
#endif
@ -77,6 +83,15 @@ do { \
(ix1) = ew_u.parts.lsw; \
} while (/*CONSTCOND*/0)
/* Get a 64-bit int from a double. */
#define EXTRACT_WORD64(ix,d) \
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix) = ew_u.xparts.w; \
} while (/*CONSTCOND*/0)
/* Get the more significant 32 bit int from a double. */
#define GET_HIGH_WORD(i,d) \
@ -105,6 +120,15 @@ do { \
(d) = iw_u.value; \
} while (/*CONSTCOND*/0)
/* Set a double from a 64-bit int. */
#define INSERT_WORD64(d,ix) \
do { \
ieee_double_shape_type iw_u; \
iw_u.xparts.w = (ix); \
(d) = iw_u.value; \
} while (/*CONSTCOND*/0)
/* Set the more significant 32 bits of a double from an int. */
#define SET_HIGH_WORD(d,v) \

View File

@ -1,3 +1,5 @@
/* $NetBSD: s_fma.c,v 1.2 2013/02/11 01:29:58 christos Exp $ */
/*-
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
@ -25,8 +27,13 @@
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_fma.c,v 1.8 2011/10/21 06:30:43 das Exp $");
#else
__RCSID("$NetBSD: s_fma.c,v 1.2 2013/02/11 01:29:58 christos Exp $");
#endif
#include <machine/ieee.h>
#include <fenv.h>
#include <float.h>
#include <math.h>
@ -117,7 +124,7 @@ add_and_denormalize(double a, double b, int scale)
if (sum.lo != 0) {
EXTRACT_WORD64(hibits, sum.hi);
bits_lost = -((int)(hibits >> 52) & 0x7ff) - scale + 1;
if (bits_lost != 1 ^ (int)(hibits & 1)) {
if ((bits_lost != 1) ^ (int)(hibits & 1)) {
/* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
EXTRACT_WORD64(lobits, sum.lo);
hibits += 1 - (((hibits ^ lobits) >> 62) & 2);
@ -216,17 +223,17 @@ fma(double x, double y, double z)
case FE_TONEAREST:
return (z);
case FE_TOWARDZERO:
if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0))
return (z);
else
return (nextafter(z, 0));
case FE_DOWNWARD:
if (x > 0.0 ^ y < 0.0)
if ((x > 0.0) ^ (y < 0.0))
return (z);
else
return (nextafter(z, -INFINITY));
default: /* FE_UPWARD */
if (x > 0.0 ^ y < 0.0)
if ((x > 0.0) ^ (y < 0.0))
return (nextafter(z, INFINITY));
else
return (z);
@ -258,8 +265,10 @@ fma(double x, double y, double z)
* the correct sign.
*/
fesetround(oround);
{
volatile double vzs = zs; /* XXX gcc CSE bug workaround */
return (xy.hi + vzs + ldexp(xy.lo, spread));
}
}
if (oround != FE_TONEAREST) {
@ -279,6 +288,6 @@ fma(double x, double y, double z)
return (add_and_denormalize(r.hi, adj, spread));
}
#if (LDBL_MANT_DIG == 53)
#ifndef EXT_FRACBITS
__weak_reference(fma, fmal);
#endif

View File

@ -1,3 +1,5 @@
/* $NetBSD: s_fmaf.c,v 1.2 2013/02/11 01:29:58 christos Exp $ */
/*-
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
@ -25,7 +27,11 @@
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_fmaf.c,v 1.3 2011/10/15 04:16:58 das Exp $");
#else
__RCSID("$NetBSD: s_fmaf.c,v 1.2 2013/02/11 01:29:58 christos Exp $");
#endif
#include <fenv.h>

View File

@ -1,3 +1,5 @@
/* $NetBSD: s_fmal.c,v 1.2 2013/02/11 01:29:58 christos Exp $ */
/*-
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
@ -25,14 +27,20 @@
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_fmal.c,v 1.7 2011/10/21 06:30:43 das Exp $");
#else
__RCSID("$NetBSD: s_fmal.c,v 1.2 2013/02/11 01:29:58 christos Exp $");
#endif
#include <machine/ieee.h>
#include <fenv.h>
#include <float.h>
#include <math.h>
#include "fpmath.h"
#include "math_private.h"
#ifdef EXT_EXPBITS
/*
* A struct dd represents a floating-point number with twice the precision
* of a long double. We maintain the invariant that "hi" stores the high-order
@ -75,12 +83,12 @@ static inline long double
add_adjusted(long double a, long double b)
{
struct dd sum;
union IEEEl2bits u;
union ieee_ext_u u;
sum = dd_add(a, b);
if (sum.lo != 0) {
u.e = sum.hi;
if ((u.bits.manl & 1) == 0)
u.extu_ld = sum.hi;
if ((u.extu_ext.ext_fracl & 1) == 0)
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return (sum.hi);
@ -96,7 +104,7 @@ add_and_denormalize(long double a, long double b, int scale)
{
struct dd sum;
int bits_lost;
union IEEEl2bits u;
union ieee_ext_u u;
sum = dd_add(a, b);
@ -111,12 +119,12 @@ add_and_denormalize(long double a, long double b, int scale)
* break the ties manually.
*/
if (sum.lo != 0) {
u.e = sum.hi;
bits_lost = -u.bits.exp - scale + 1;
if (bits_lost != 1 ^ (int)(u.bits.manl & 1))
u.extu_ld = sum.hi;
bits_lost = -u.extu_ext.ext_exp - scale + 1;
if ((bits_lost != 1) ^ (int)(u.extu_ext.ext_fracl & 1))
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return (ldexp(sum.hi, scale));
return (ldexp((double)sum.hi, scale));
}
/*
@ -204,18 +212,18 @@ fmal(long double x, long double y, long double z)
case FE_TONEAREST:
return (z);
case FE_TOWARDZERO:
if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0))
return (z);
else
return (nextafterl(z, 0));
case FE_DOWNWARD:
if (x > 0.0 ^ y < 0.0)
if ((x > 0.0) ^ (y < 0.0))
return (z);
else
return (nextafterl(z, -INFINITY));
return (nextafterl(z, (long double)-INFINITY));
default: /* FE_UPWARD */
if (x > 0.0 ^ y < 0.0)
return (nextafterl(z, INFINITY));
if ((x > 0.0) ^ (y < 0.0))
return (nextafterl(z, (long double)INFINITY));
else
return (z);
}
@ -246,8 +254,10 @@ fmal(long double x, long double y, long double z)
* the correct sign.
*/
fesetround(oround);
{
volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
return (xy.hi + vzs + ldexpl(xy.lo, spread));
}
}
if (oround != FE_TONEAREST) {
@ -266,3 +276,4 @@ fmal(long double x, long double y, long double z)
else
return (add_and_denormalize(r.hi, adj, spread));
}
#endif