Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for FreeBSD x86 reals #7689

Merged
merged 4 commits into from
Nov 13, 2020
Merged

Add support for FreeBSD x86 reals #7689

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
c955773
std.complex: Relax unittests for FreeBSD x86 53-bit precision reals
ibuclaw Nov 9, 2020
51c457d
std.conv: Fix unittests for FreeBSD x86 53-bit precision reals
ibuclaw Nov 9, 2020
a7ba764
std.internal.math.gammafunction: Define MAXGAMMA, MAXLOG, MINLOG for …
ibuclaw Nov 9, 2020
ab8221e
std.math: Add support for FreeBSD x86 53-bit precision reals
ibuclaw Nov 9, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 20 additions & 15 deletions std/complex.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -901,6 +901,17 @@ Complex!(CommonType!(T, U)) fromPolar(T, U)(const T modulus, const U argument)
assert(approxEqual(z.im, 1.0L, real.epsilon));
}

version (StdUnittest)
{
// Helper function for comparing two Complex numbers.
int ceqrel(T)(const Complex!T x, const Complex!T y) @safe pure nothrow @nogc
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

private

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah you just moved it, nvm

{
import std.math : feqrel;
const r = feqrel(x.re, y.re);
const i = feqrel(x.im, y.im);
return r < i ? r : i;
}
}

/**
Trigonometric functions on complex numbers.
Expand All @@ -920,9 +931,14 @@ Complex!T sin(T)(Complex!T z) @safe pure nothrow @nogc
{
static import std.math;
assert(sin(complex(0.0)) == 0.0);
assert(sin(complex(2.0L, 0)) == std.math.sin(2.0L));
assert(sin(complex(2.0, 0)) == std.math.sin(2.0));
}

@safe pure nothrow unittest
{
static import std.math;
assert(ceqrel(sin(complex(2.0L, 0)), complex(std.math.sin(2.0L))) >= real.mant_dig - 1);
}

/// ditto
Complex!T cos(T)(Complex!T z) @safe pure nothrow @nogc
Expand All @@ -937,26 +953,15 @@ Complex!T cos(T)(Complex!T z) @safe pure nothrow @nogc
{
static import std.math;
assert(cos(complex(0.0)) == 1.0);
assert(cos(complex(1.3L, 0.0)) == std.math.cos(1.3L));
assert(cos(complex(0.0L, 5.2L)) == std.math.cosh(5.2L));
}

version (StdUnittest)
{
int ceqrel(T)(const Complex!T x, const Complex!T y) @safe pure nothrow @nogc
{
import std.math : feqrel;
const r = feqrel(x.re, y.re);
const i = feqrel(x.im, y.im);
return r < i ? r : i;
}
assert(cos(complex(1.3, 0.0)) == std.math.cos(1.3));
assert(cos(complex(0.0, 5.2)) == std.math.cosh(5.2));
}

@safe pure nothrow unittest
{
static import std.math;
assert(ceqrel(cos(complex(0, 5.2L)), complex(std.math.cosh(5.2L), 0.0L)) >= real.mant_dig - 1);
assert(cos(complex(1.3L)) == std.math.cos(1.3L));
assert(ceqrel(cos(complex(1.3L)), complex(std.math.cos(1.3L))) >= real.mant_dig - 1);
}

/// ditto
Expand Down
26 changes: 20 additions & 6 deletions std/conv.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1750,6 +1750,8 @@ if (!isImplicitlyConvertible!(S, T) && isAssociativeArray!S &&
}
static void testFloatingToIntegral(Floating, Integral)()
{
import std.math : floatTraits, RealFormat;

bool convFails(Source, Target, E)(Source src)
{
try
Expand Down Expand Up @@ -1781,18 +1783,23 @@ if (!isImplicitlyConvertible!(S, T) && isAssociativeArray!S &&
{
a = -a; // -Integral.min not representable as an Integral
assert(convFails!(Floating, Integral, ConvOverflowException)(a)
|| Floating.sizeof <= Integral.sizeof);
|| Floating.sizeof <= Integral.sizeof
|| floatTraits!Floating.realFormat == RealFormat.ieeeExtended53);
}
a = 0.0 + Integral.min;
assert(to!Integral(a) == Integral.min);
--a; // no more representable as an Integral
assert(convFails!(Floating, Integral, ConvOverflowException)(a)
|| Floating.sizeof <= Integral.sizeof);
|| Floating.sizeof <= Integral.sizeof
|| floatTraits!Floating.realFormat == RealFormat.ieeeExtended53);
a = 0.0 + Integral.max;
assert(to!Integral(a) == Integral.max || Floating.sizeof <= Integral.sizeof);
assert(to!Integral(a) == Integral.max
|| Floating.sizeof <= Integral.sizeof
|| floatTraits!Floating.realFormat == RealFormat.ieeeExtended53);
++a; // no more representable as an Integral
assert(convFails!(Floating, Integral, ConvOverflowException)(a)
|| Floating.sizeof <= Integral.sizeof);
|| Floating.sizeof <= Integral.sizeof
|| floatTraits!Floating.realFormat == RealFormat.ieeeExtended53);
// convert a value with a fractional part
a = 3.14;
assert(to!Integral(a) == 3);
Expand Down Expand Up @@ -3278,7 +3285,9 @@ if (isInputRange!Source && isSomeChar!(ElementType!Source) && !is(Source == enum
@system unittest
{
// @system because strtod is not @safe.
static if (real.mant_dig == 53)
import std.math : floatTraits, RealFormat;

static if (floatTraits!real.realFormat == RealFormat.ieeeDouble)
{
import core.stdc.stdlib, std.exception, std.math;

Expand Down Expand Up @@ -3361,7 +3370,8 @@ if (isInputRange!Source && isSomeChar!(ElementType!Source) && !is(Source == enum
{
ushort[8] value;
}
else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended)
else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended ||
floatTraits!real.realFormat == RealFormat.ieeeExtended53)
{
ushort[5] value;
}
Expand All @@ -3383,6 +3393,8 @@ if (isInputRange!Source && isSomeChar!(ElementType!Source) && !is(Source == enum
enum s = "0x1.FFFFFFFFFFFFFFFFFFFFFFFFFFFFp-16382";
else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended)
enum s = "0x1.FFFFFFFFFFFFFFFEp-16382";
else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended53)
enum s = "0x1.FFFFFFFFFFFFFFFEp-16382";
else static if (floatTraits!real.realFormat == RealFormat.ieeeDouble)
enum s = "0x1.FFFFFFFFFFFFFFFEp-1000";
else
Expand All @@ -3400,6 +3412,8 @@ if (isInputRange!Source && isSomeChar!(ElementType!Source) && !is(Source == enum
else
ld1 = strtold(s.ptr, null);
}
else static if (floatTraits!real.realFormat == RealFormat.ieeeExtended53)
ld1 = 0x1.FFFFFFFFFFFFFFFEp-16382L; // strtold rounds to 53 bits.
else
ld1 = strtold(s.ptr, null);

Expand Down
7 changes: 7 additions & 0 deletions std/internal/math/gammafunction.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -253,6 +253,8 @@ static if (floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
enum real MAXGAMMA = 1755.5483429L;
else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
enum real MAXGAMMA = 1755.5483429L;
else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended53)
enum real MAXGAMMA = 1755.5483429L;
else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
enum real MAXGAMMA = 171.6243769L;
else
Expand Down Expand Up @@ -603,6 +605,11 @@ else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended)
enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max)
enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
}
else static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended53)
{
enum real MAXLOG = 0x1.62e42fefa39ef358p+13L; // log(real.max)
enum real MINLOG = -0x1.6436716d5406e6d8p+13L; // log(real.min_normal*real.epsilon) = log(smallest denormal)
}
else static if (floatTraits!(real).realFormat == RealFormat.ieeeDouble)
{
enum real MAXLOG = 0x1.62e42fefa39efp+9L; // log(real.max)
Expand Down
54 changes: 38 additions & 16 deletions std/math.d
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 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)
static immutable T[3] P = [
Expand Down Expand Up @@ -1837,7 +1838,8 @@ private T expImpl(T)(T x) @safe pure nothrow @nogc
[-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 =
[ // x exp(x)
Expand Down Expand Up @@ -2777,7 +2779,8 @@ if (isFloatingPoint!T)
alias F = floatTraits!T;

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 exponent is non-zero
Expand Down Expand Up @@ -3120,7 +3123,8 @@ if (isFloatingPoint!T)
y.rv = x;

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)
{
Expand Down Expand Up @@ -3393,6 +3397,7 @@ float ldexp(float n, int exp) @safe pure nothrow @nogc { return core.math.ldex
@safe pure nothrow @nogc unittest
{
static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended ||
floatTraits!(real).realFormat == RealFormat.ieeeExtended53 ||
floatTraits!(real).realFormat == RealFormat.ieeeQuadruple)
{
assert(ldexp(1.0L, -16384) == 0x1p-16384L);
Expand Down Expand Up @@ -4776,12 +4781,16 @@ long lrint(real x) @trusted pure nothrow @nogc

return sign ? -result : result;
}
else static if (F.realFormat == RealFormat.ieeeExtended)
else static if (F.realFormat == RealFormat.ieeeExtended ||
F.realFormat == RealFormat.ieeeExtended53)
{
long result;

// 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);
uint* vi = cast(uint*)(&x);
Expand Down Expand Up @@ -5998,7 +6007,8 @@ if (isFloatingPoint!(X))
// At least one bit among the least significant 52 bits should be set.
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 ulong ps = *cast(ulong *)&x;
Expand Down Expand Up @@ -6222,7 +6232,8 @@ bool isSubnormal(X)(X x) @trusted pure nothrow @nogc
return (e == 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;
long* ps = cast(long *)&x;
Expand Down Expand Up @@ -6283,7 +6294,8 @@ if (isFloatingPoint!(X))
return ((*cast(ulong *)&x) & 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 = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]);
const ulong ps = *cast(ulong *)&x;
Expand Down Expand Up @@ -6644,7 +6656,8 @@ if (isFloatingPoint!F || isIntegral!F)
real NaN(ulong payload) @trusted pure nothrow @nogc
{
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
// not implied but a real bit which is stored in the real)
Expand Down Expand Up @@ -6866,11 +6879,14 @@ real nextUp(real x) @trusted pure nothrow @nogc
}
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...
ushort *pe = cast(ushort *)&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)
{
Expand All @@ -6881,7 +6897,7 @@ real nextUp(real x) @trusted pure nothrow @nogc
if (pe[F.EXPPOS_SHORT] & 0x8000)
{
// Negative number -- need to decrease the significand
--*ps;
*ps -= EPSILON;
// Need to mask with 0x7FFF... so subnormals are treated correctly.
if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF)
{
Expand All @@ -6906,7 +6922,7 @@ real nextUp(real x) @trusted pure nothrow @nogc
{
// Positive number -- need to increase the significand.
// Works automatically for positive zero.
++*ps;
*ps += EPSILON;
if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0)
{
// change in exponent
Expand Down Expand Up @@ -8069,6 +8085,7 @@ if (isFloatingPoint!(X))
static if (F.realFormat == RealFormat.ieeeSingle
|| F.realFormat == RealFormat.ieeeDouble
|| F.realFormat == RealFormat.ieeeExtended
|| F.realFormat == RealFormat.ieeeExtended53
|| F.realFormat == RealFormat.ieeeQuadruple)
{
if (x == y)
Expand Down Expand Up @@ -9685,7 +9702,8 @@ do

alias F = floatTraits!(T);
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
// 79-bit reals...
Expand Down Expand Up @@ -10003,7 +10021,8 @@ T floorImpl(T)(const T x) @trusted pure nothrow @nogc
else
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;

Expand Down Expand Up @@ -10050,7 +10069,10 @@ T floorImpl(T)(const T x) @trusted pure nothrow @nogc
}
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.
while (exp >= 16)
Expand Down