mirrors/phobos
1
0
Fork 0
mirror of https://github.com/dlang/phobos.git synced 2025-04-29 14:40:30 +03:00
Code Issues Projects Releases Packages Wiki Activity Actions 1

std.math: Add support for FreeBSD x86 53-bit precision reals

Browse source
This commit is contained in:
Iain Buclaw 2020-11-09 14:24:32 +01:00
parent a7ba7644ff
commit ab8221efb0
2 changed files with 40 additions and 17 deletions
Download patch file Download diff file Expand all files Collapse all files

3
std/complex.d
View file

@ -936,7 +936,8 @@ Complex!T sin(T)(Complex!T z) @safe pure nothrow @nogc
@safe pure nothrow unittest @safe pure nothrow unittest
{ {
assert(ceqrel(sin(complex(2.0L, 0)), std.math.sin(2.0L)) >= real.mant_dig - 1); static import std.math;
assert(ceqrel(sin(complex(2.0L, 0)), complex(std.math.sin(2.0L))) >= real.mant_dig - 1);
} }
/// ditto /// ditto

54
std/math.d
View file

@ -1715,7 +1715,8 @@ private T expImpl(T)(T x) @safe pure nothrow @nogc
enum T OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024) enum T OF = 7.09782712893383996732E2; // ln((1-2^-53) * 2^1024)
enum T UF = -7.451332191019412076235E2; // ln(2^-1075) enum T UF = -7.451332191019412076235E2; // ln(2^-1075)
} }
else static if (F.realFormat == RealFormat.ieeeExtended) else static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
// Coefficients for exp(x) // Coefficients for exp(x)
static immutable T[3] P = [ static immutable T[3] P = [
@ -1837,7 +1838,8 @@ private T expImpl(T)(T x) @safe pure nothrow @nogc
[-0x1p+30L, 0 ], // far underflow [-0x1p+30L, 0 ], // far underflow
]; ];
} }
else static if (realFormat == RealFormat.ieeeExtended) else static if (realFormat == RealFormat.ieeeExtended ||
realFormat == RealFormat.ieeeExtended53)
{ {
static immutable T[2][] exptestpoints = static immutable T[2][] exptestpoints =
[ // x exp(x) [ // x exp(x)
@ -2777,7 +2779,8 @@ if (isFloatingPoint!T)
alias F = floatTraits!T; alias F = floatTraits!T;
ex = vu[F.EXPPOS_SHORT] & F.EXPMASK; ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
static if (F.realFormat == RealFormat.ieeeExtended) static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
if (ex) if (ex)
{ // If exponent is non-zero { // If exponent is non-zero
@ -3120,7 +3123,8 @@ if (isFloatingPoint!T)
y.rv = x; y.rv = x;
int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK; int ex = y.vu[F.EXPPOS_SHORT] & F.EXPMASK;
static if (F.realFormat == RealFormat.ieeeExtended) static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
if (ex) if (ex)
{ {
@ -3393,6 +3397,7 @@ float ldexp(float n, int exp) @safe pure nothrow @nogc { return core.math.ldex
@safe pure nothrow @nogc unittest @safe pure nothrow @nogc unittest
{ {
static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended || static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
floatTraits!(real).realFormat == RealFormat.ieeeQuadruple) floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
{ {
assert(ldexp(1.0L, -16384) == 0x1p-16384L); assert(ldexp(1.0L, -16384) == 0x1p-16384L);
@ -4776,12 +4781,16 @@ long lrint(real x) @trusted pure nothrow @nogc
return sign ? -result : result; return sign ? -result : result;
} }
else static if (F.realFormat == RealFormat.ieeeExtended) else static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
long result; long result;
// Rounding limit when casting from real(80-bit) to ulong. // Rounding limit when casting from real(80-bit) to ulong.
enum real OF = 9.22337203685477580800E18L; static if (F.realFormat == RealFormat.ieeeExtended)
enum real OF = 9.22337203685477580800E18L;
else
enum real OF = 4.50359962737049600000E15L;
ushort* vu = cast(ushort*)(&x); ushort* vu = cast(ushort*)(&x);
uint* vi = cast(uint*)(&x); uint* vi = cast(uint*)(&x);
@ -5998,7 +6007,8 @@ if (isFloatingPoint!(X))
// At least one bit among the least significant 52 bits should be set. // At least one bit among the least significant 52 bits should be set.
return (p & 0x7FFF_FFFF_FFFF_FFFF) > 0x7FF0_0000_0000_0000; return (p & 0x7FFF_FFFF_FFFF_FFFF) > 0x7FF0_0000_0000_0000;
} }
else static if (F.realFormat == RealFormat.ieeeExtended) else static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
const ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]; const ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
const ulong ps = *cast(ulong *)&x; const ulong ps = *cast(ulong *)&x;
@ -6222,7 +6232,8 @@ bool isSubnormal(X)(X x) @trusted pure nothrow @nogc
return (e == 0 && return (e == 0 &&
((ps[MANTISSA_LSB]|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF)) != 0)); ((ps[MANTISSA_LSB]|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF)) != 0));
} }
else static if (F.realFormat == RealFormat.ieeeExtended) else static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
ushort* pe = cast(ushort *)&x; ushort* pe = cast(ushort *)&x;
long* ps = cast(long *)&x; long* ps = cast(long *)&x;
@ -6283,7 +6294,8 @@ if (isFloatingPoint!(X))
return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF) return ((*cast(ulong *)&x) & 0x7FFF_FFFF_FFFF_FFFF)
== 0x7FF0_0000_0000_0000; == 0x7FF0_0000_0000_0000;
} }
else static if (F.realFormat == RealFormat.ieeeExtended) else static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
const ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]); const ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]);
const ulong ps = *cast(ulong *)&x; const ulong ps = *cast(ulong *)&x;
@ -6644,7 +6656,8 @@ if (isFloatingPoint!F || isIntegral!F)
real NaN(ulong payload) @trusted pure nothrow @nogc real NaN(ulong payload) @trusted pure nothrow @nogc
{ {
alias F = floatTraits!(real); alias F = floatTraits!(real);
static if (F.realFormat == RealFormat.ieeeExtended) static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
// real80 (in x86 real format, the implied bit is actually // real80 (in x86 real format, the implied bit is actually
// not implied but a real bit which is stored in the real) // not implied but a real bit which is stored in the real)
@ -6866,11 +6879,14 @@ real nextUp(real x) @trusted pure nothrow @nogc
} }
return x; return x;
} }
else static if (F.realFormat == RealFormat.ieeeExtended) else static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
// For 80-bit reals, the "implied bit" is a nuisance... // For 80-bit reals, the "implied bit" is a nuisance...
ushort *pe = cast(ushort *)&x; ushort *pe = cast(ushort *)&x;
ulong *ps = cast(ulong *)&x; ulong *ps = cast(ulong *)&x;
// EPSILON is 1 for 64-bit, and 2048 for 53-bit precision reals.
enum ulong EPSILON = 2UL ^^ (64 - real.mant_dig);
if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK) if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK)
{ {
@ -6881,7 +6897,7 @@ real nextUp(real x) @trusted pure nothrow @nogc
if (pe[F.EXPPOS_SHORT] & 0x8000) if (pe[F.EXPPOS_SHORT] & 0x8000)
{ {
// Negative number -- need to decrease the significand // Negative number -- need to decrease the significand
--*ps; *ps -= EPSILON;
// Need to mask with 0x7FFF... so subnormals are treated correctly. // Need to mask with 0x7FFF... so subnormals are treated correctly.
if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF) if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF)
{ {
@ -6906,7 +6922,7 @@ real nextUp(real x) @trusted pure nothrow @nogc
{ {
// Positive number -- need to increase the significand. // Positive number -- need to increase the significand.
// Works automatically for positive zero. // Works automatically for positive zero.
++*ps; *ps += EPSILON;
if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0) if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0)
{ {
// change in exponent // change in exponent
@ -8069,6 +8085,7 @@ if (isFloatingPoint!(X))
static if (F.realFormat == RealFormat.ieeeSingle static if (F.realFormat == RealFormat.ieeeSingle
|| F.realFormat == RealFormat.ieeeDouble || F.realFormat == RealFormat.ieeeDouble
|| F.realFormat == RealFormat.ieeeExtended || F.realFormat == RealFormat.ieeeExtended
|| F.realFormat == RealFormat.ieeeExtended53
|| F.realFormat == RealFormat.ieeeQuadruple) || F.realFormat == RealFormat.ieeeQuadruple)
{ {
if (x == y) if (x == y)
@ -9685,7 +9702,8 @@ do
alias F = floatTraits!(T); alias F = floatTraits!(T);
T u; T u;
static if (F.realFormat == RealFormat.ieeeExtended) static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
// There's slight additional complexity because they are actually // There's slight additional complexity because they are actually
// 79-bit reals... // 79-bit reals...
@ -10003,7 +10021,8 @@ T floorImpl(T)(const T x) @trusted pure nothrow @nogc
else else
int pos = 3; int pos = 3;
} }
else static if (F.realFormat == RealFormat.ieeeExtended) else static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{ {
int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff; int exp = (y.vu[F.EXPPOS_SHORT] & 0x7fff) - 0x3fff;
@ -10050,7 +10069,10 @@ T floorImpl(T)(const T x) @trusted pure nothrow @nogc
} }
else else
{ {
exp = (T.mant_dig - 1) - exp; static if (F.realFormat == RealFormat.ieeeExtended53)
exp = (T.mant_dig + 11 - 1) - exp; // mant_dig is really 64
else
exp = (T.mant_dig - 1) - exp;
// Zero 16 bits at a time. // Zero 16 bits at a time.
while (exp >= 16) while (exp >= 16)