d/freebsd-dev
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update scalbn* functions to the musl versions

Browse Source 
The only diff compared to musl is a minor change to scalbnl() to replace
musl's union ldshape with union IEEEl2bits.
This fixes the scalbn tests on non-x86 (since x86 has an assembly version
that is used instead).

Musl commit messages:
commit 8c44a060243f04283ca68dad199aab90336141db
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Mon Apr 3 02:38:13 2017 +0200

    fix scalbn when result is in the subnormal range

    in nearest rounding mode scalbn could introduce double rounding error
    when an intermediate value and the final result were both in the
    subnormal range e.g.

      scalbn(0x1.7ffffffffffffp-1, -1073)

    returned 0x1p-1073 instead of 0x1p-1074, because the intermediate
    computation got rounded to 0x1.8p-1023.

    with the fix an intermediate value can only be in the subnormal range
    if the final result is 0 which is correct even after double rounding.
    (there still can be two roundings so signals may be raised twice, but
    that's only observable with trapping exceptions which is not supported.)

commit 2eaed464e2080d8321d3903b71086a1ecfc4ee4a
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Wed Sep 4 15:52:54 2013 +0000

    math: use float_t and double_t in scalbnf and scalbn

    remove STRICT_ASSIGN (c99 semantics is assumed) and use the conventional
    union to prepare the scaling factor (so libm.h is no longer needed)

commit 1b77b9072f374bd26eb0574b83a0d5f18d75ec60
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Thu Aug 15 10:07:46 2013 +0000

    math: minor scalbn*.c simplification

commit c4359e01303da2755fe7e8033826b132eb3659b1
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Tue Nov 13 10:55:35 2012 +0100

    math: excess precision fix modf, modff, scalbn, scalbnf

    old code was correct only if the result was stored (without the
    excess precision) or musl was compiled with -ffloat-store.
    now we use STRICT_ASSIGN to work around the issue.
    (see note 160 in c11 section 6.8.6.4)

commit 666271c105e4137bdfa195e217799d74143370d4
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Tue Nov 13 10:30:40 2012 +0100

    math: fix scalbn and scalbnf on overflow/underflow

    old code was correct only if the result was stored (without the
    excess precision) or musl was compiled with -ffloat-store.
    (see note 160 in n1570.pdf section 6.8.6.4)

commit 8051e08e10d2b739fcfcbc6bc7466e8d77fa49f1
Author: nsz <nsz@port70.net>
Date:   Mon Mar 19 10:54:07 2012 +0100

    simplify scalbn*.c implementations

    The old scalbn.c was wrong and slow, the new one is just slow.
    (scalbn(0x1p+1023,-2097) should give 0x1p-1074, but the old code gave 0)

Reviewed By:	dim
Differential Revision: https://reviews.freebsd.org/D28872
This commit is contained in:
Alex Richardson 2021-02-24 16:56:53 +00:00
parent 98202829d1
commit 760b2ffc55
3 changed files with 90 additions and 166 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

86
lib/msun/src/s_scalbn.c
View File

@ -1,63 +1,35 @@
/* @(#)s_scalbn.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 <math.h>
#include <stdint.h>
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
/*
* scalbn (double x, int n)
* scalbn(x,n) returns x* 2**n computed by exponent
* manipulation rather than by actually performing an
* exponentiation or a multiplication.
*/
#include <float.h>
#include "math.h"
#include "math_private.h"
static const double
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
huge = 1.0e+300,
tiny = 1.0e-300;
double
scalbn (double x, int n)
double scalbn(double x, int n)
{
int32_t k,hx,lx;
EXTRACT_WORDS(hx,lx,x);
k = (hx&0x7ff00000)>>20; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
GET_HIGH_WORD(hx,x);
k = ((hx&0x7ff00000)>>20) - 54;
if (n< -50000) return tiny*x; /*underflow*/
}
if (k==0x7ff) return x+x; /* NaN or Inf */
k = k+n;
if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */
if (k > 0) /* normal result */
{SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
if (k <= -54) {
if (n > 50000) /* in case integer overflow in n+k */
return huge*copysign(huge,x); /*overflow*/
else
return tiny*copysign(tiny,x); /*underflow*/
union {double f; uint64_t i;} u;
double_t y = x;
if (n > 1023) {
y *= 0x1p1023;
n -= 1023;
if (n > 1023) {
y *= 0x1p1023;
n -= 1023;
if (n > 1023)
n = 1023;
}
} else if (n < -1022) {
/* make sure final n < -53 to avoid double
rounding in the subnormal range */
y *= 0x1p-1022 * 0x1p53;
n += 1022 - 53;
if (n < -1022) {
y *= 0x1p-1022 * 0x1p53;
n += 1022 - 53;
if (n < -1022)
n = -1022;
}
}
k += 54; /* subnormal result */
SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
return x*twom54;
u.i = (uint64_t)(0x3ff+n)<<52;
x = y * u.f;
return x;
}
#if (LDBL_MANT_DIG == 53)

78
lib/msun/src/s_scalbnf.c
View File

@ -1,57 +1,33 @@
/* s_scalbnf.c -- float version of s_scalbn.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
#include <math.h>
#include <stdint.h>
/*
* ====================================================
* 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$");
#include "math.h"
#include "math_private.h"
static const float
two25 = 3.355443200e+07, /* 0x4c000000 */
twom25 = 2.9802322388e-08, /* 0x33000000 */
huge = 1.0e+30,
tiny = 1.0e-30;
float
scalbnf (float x, int n)
float scalbnf(float x, int n)
{
int32_t k,ix;
GET_FLOAT_WORD(ix,x);
k = (ix&0x7f800000)>>23; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((ix&0x7fffffff)==0) return x; /* +-0 */
x *= two25;
GET_FLOAT_WORD(ix,x);
k = ((ix&0x7f800000)>>23) - 25;
if (n< -50000) return tiny*x; /*underflow*/
}
if (k==0xff) return x+x; /* NaN or Inf */
k = k+n;
if (k > 0xfe) return huge*copysignf(huge,x); /* overflow */
if (k > 0) /* normal result */
{SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
if (k <= -25) {
if (n > 50000) /* in case integer overflow in n+k */
return huge*copysignf(huge,x); /*overflow*/
else
return tiny*copysignf(tiny,x); /*underflow*/
union {float f; uint32_t i;} u;
float_t y = x;
if (n > 127) {
y *= 0x1p127f;
n -= 127;
if (n > 127) {
y *= 0x1p127f;
n -= 127;
if (n > 127)
n = 127;
}
} else if (n < -126) {
y *= 0x1p-126f * 0x1p24f;
n += 126 - 24;
if (n < -126) {
y *= 0x1p-126f * 0x1p24f;
n += 126 - 24;
if (n < -126)
n = -126;
}
}
k += 25; /* subnormal result */
SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
return x*twom25;
u.i = (uint32_t)(0x7f+n)<<23;
x = y * u.f;
return x;
}
__strong_reference(scalbnf, ldexpf);

92
lib/msun/src/s_scalbnl.c
View File

@ -1,18 +1,7 @@
/* @(#)s_scalbn.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$");
#include <math.h>
#include <float.h>
#include "math_private.h"
#include "fpmath.h"
/*
* scalbnl (long double x, int n)
* scalbnl(x,n) returns x* 2**n computed by exponent
@ -20,52 +9,39 @@ __FBSDID("$FreeBSD$");
* exponentiation or a multiplication.
*/
/*
* We assume that a long double has a 15-bit exponent. On systems
* where long double is the same as double, scalbnl() is an alias
* for scalbn(), so we don't use this routine.
*/
#include <float.h>
#include <math.h>
#include "fpmath.h"
#if LDBL_MAX_EXP != 0x4000
#error "Unsupported long double format"
#endif
static const long double
huge = 0x1p16000L,
tiny = 0x1p-16000L;
long double
scalbnl (long double x, int n)
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double scalbnl(long double x, int n)
{
return scalbn(x, n);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
long double scalbnl(long double x, int n)
{
union IEEEl2bits u;
int k;
u.e = x;
k = u.bits.exp; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((u.bits.manh|u.bits.manl)==0) return x; /* +-0 */
u.e *= 0x1p+128;
k = u.bits.exp - 128;
if (n< -50000) return tiny*x; /*underflow*/
}
if (k==0x7fff) return x+x; /* NaN or Inf */
k = k+n;
if (k >= 0x7fff) return huge*copysignl(huge,x); /* overflow */
if (k > 0) /* normal result */
{u.bits.exp = k; return u.e;}
if (k <= -128) {
if (n > 50000) /* in case integer overflow in n+k */
return huge*copysign(huge,x); /*overflow*/
else
return tiny*copysign(tiny,x); /*underflow*/
if (n > 16383) {
x *= 0x1p16383L;
n -= 16383;
if (n > 16383) {
x *= 0x1p16383L;
n -= 16383;
if (n > 16383)
n = 16383;
}
} else if (n < -16382) {
x *= 0x1p-16382L * 0x1p113L;
n += 16382 - 113;
if (n < -16382) {
x *= 0x1p-16382L * 0x1p113L;
n += 16382 - 113;
if (n < -16382)
n = -16382;
}
}
k += 128; /* subnormal result */
u.bits.exp = k;
return u.e*0x1p-128;
u.e = 1.0;
u.xbits.expsign = 0x3fff + n;
return x * u.e;
}
#endif
__strong_reference(scalbnl, ldexpl);