Skip to content

Commit

Permalink
Exposing IRootFunctions.Hypot and IRootFunctions.Root (#71010)
Browse files Browse the repository at this point in the history 
* Exposing IRootFunctions.Hypot and IRootFunctions.Root

* Adding tests for IRootFunctions.Hypot and IRootFunctions.Root
  • Loading branch information
@tannergooding
tannergooding authored Jun 20, 2022
1 parent 19fcc31 commit 715c71a
Show file tree
Hide file tree
Showing 12 changed files with 1,404 additions and 69 deletions.
29 changes: 28 additions & 1 deletion THIRD-PARTY-NOTICES.TXT
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1108,4 +1108,31 @@ Copyright (C) 1999 Lucent Technologies
Excerpted from 'The Practice of Programming
by Brian W. Kernighan and Rob Pike

You may use this code for any purpose, as long as you leave the copyright notice and book citation attached.
You may use this code for any purpose, as long as you leave the copyright notice and book citation attached.

License notice for amd/aocl-libm-ose
-------------------------------

Copyright (C) 2008-2020 Advanced Micro Devices, Inc. All rights reserved.

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 copyright holder 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 COPYRIGHT HOLDERS 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 COPYRIGHT HOLDER 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.
290 changes: 284 additions & 6 deletions src/libraries/System.Private.CoreLib/src/System/Double.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
===========================================================*/

using System.Buffers.Binary;
using System.Diagnostics;
using System.Diagnostics.CodeAnalysis;
using System.Globalization;
using System.Numerics;
Expand Down Expand Up @@ -106,6 +107,9 @@ public readonly struct Double
internal const ulong MinTrailingSignificand = 0x0000_0000_0000_0000;
internal const ulong MaxTrailingSignificand = 0x000F_FFFF_FFFF_FFFF;

internal const int TrailingSignificandLength = 52;
internal const int SignificandLength = TrailingSignificandLength + 1;

internal ushort BiasedExponent
{
get
Expand Down Expand Up @@ -917,7 +921,7 @@ bool IFloatingPoint<double>.TryWriteSignificandLittleEndian(Span<byte> destinati
public static double Clamp(double value, double min, double max) => Math.Clamp(value, min, max);

/// <inheritdoc cref="INumber{TSelf}.CopySign(TSelf, TSelf)" />
public static double CopySign(double x, double y) => Math.CopySign(x, y);
public static double CopySign(double value, double sign) => Math.CopySign(value, sign);

/// <inheritdoc cref="INumber{TSelf}.Max(TSelf, TSelf)" />
public static double Max(double x, double y) => Math.Max(x, y);
Expand Down Expand Up @@ -1353,15 +1357,289 @@ private static bool TryConvertTo<TOther>(double value, [NotNullWhen(true)] out T
/// <inheritdoc cref="IRootFunctions{TSelf}.Cbrt(TSelf)" />
public static double Cbrt(double x) => Math.Cbrt(x);

// /// <inheritdoc cref="IRootFunctions{TSelf}.Hypot(TSelf, TSelf)" />
// public static double Hypot(double x, double y) => Math.Hypot(x, y);
/// <inheritdoc cref="IRootFunctions{TSelf}.Hypot(TSelf, TSelf)" />
public static double Hypot(double x, double y)
{
// This code is based on `hypot` from amd/aocl-libm-ose
// Copyright (C) 2008-2020 Advanced Micro Devices, Inc. All rights reserved.
//
// Licensed under the BSD 3-Clause "New" or "Revised" License
// See THIRD-PARTY-NOTICES.TXT for the full license text

double result;

if (IsFinite(x) && IsFinite(y))
{
double ax = Abs(x);
double ay = Abs(y);

if (ax == 0.0f)
{
result = ay;
}
else if (ay == 0.0f)
{
result = ax;
}
else
{
ulong xBits = BitConverter.DoubleToUInt64Bits(ax);
ulong yBits = BitConverter.DoubleToUInt64Bits(ay);

uint xExp = (uint)((xBits >> BiasedExponentShift) & ShiftedExponentMask);
uint yExp = (uint)((yBits >> BiasedExponentShift) & ShiftedExponentMask);

int expDiff = (int)(xExp - yExp);
double expFix = 1.0;

if ((expDiff <= (SignificandLength + 1)) && (expDiff >= (-SignificandLength - 1)))
{
if ((xExp > (ExponentBias + 500)) || (yExp > (ExponentBias + 500)))
{
// To prevent overflow, scale down by 2^+600
expFix = 4.149515568880993E+180;

xBits -= 0x2580000000000000;
yBits -= 0x2580000000000000;
}
else if ((xExp < (ExponentBias - 500)) || (yExp < (ExponentBias - 500)))
{
// To prevent underflow, scale up by 2^-600
expFix = 2.409919865102884E-181;

xBits += 0x2580000000000000;
yBits += 0x2580000000000000;

// For subnormal values, do an additional fixing up changing the
// adjustment to scale up by 2^601 instead and then subtract a
// correction of 2^601 to account for the implicit bit.

if (xExp == 0) // x is subnormal
{
xBits += 0x0010000000000000;

ax = BitConverter.UInt64BitsToDouble(xBits);
ax -= 9.232978617785736E-128;

xBits = BitConverter.DoubleToUInt64Bits(ax);
}

if (yExp == 0) // y is subnormal
{
yBits += 0x0010000000000000;

ay = BitConverter.UInt64BitsToDouble(yBits);
ay -= 9.232978617785736E-128;

yBits = BitConverter.DoubleToUInt64Bits(ay);
}
}

ax = BitConverter.UInt64BitsToDouble(xBits);
ay = BitConverter.UInt64BitsToDouble(yBits);

if (ax < ay)
{
// Sort so ax is greater than ay
double tmp = ax;

ax = ay;
ay = tmp;

ulong tmpBits = xBits;

xBits = yBits;
yBits = tmpBits;
}

Debug.Assert(ax >= ay);

// Split ax and ay into a head and tail portion

double xHead = BitConverter.UInt64BitsToDouble(xBits & 0xFFFF_FFFF_F800_0000);
double yHead = BitConverter.UInt64BitsToDouble(yBits & 0xFFFF_FFFF_F800_0000);

double xTail = ax - xHead;
double yTail = ay - yHead;

// Compute (x * x) + (y * y) with extra precision
//
// This includes taking into account expFix which may
// cause an underflow or overflow, but if it does that
// will still be the correct result.

double xx = ax * ax;
double yy = ay * ay;

double rHead = xx + yy;
double rTail = (xx - rHead) + yy;

rTail += (xHead * xHead) - xx;
rTail += 2 * xHead * xTail;
rTail += xTail * xTail;

if (expDiff == 0)
{
// We only need to do extra accounting when ax and ay have equal exponents

rTail += (yHead * yHead) - yy;
rTail += 2 * yHead * yTail;
rTail += yTail * yTail;
}

result = Sqrt(rHead + rTail) * expFix;
}
else
{
// x or y is insignificant compared to the other
result = x + y;
}
}
}
else if (IsInfinity(x) || IsInfinity(y))
{
// IEEE 754 requires that we return +Infinity
// even if one of the inputs is NaN

result = PositiveInfinity;
}
else
{
// IEEE 754 requires that we return NaN
// if either input is NaN and neither is Infinity

result = NaN;
}

return result;
}

/// <inheritdoc cref="IRootFunctions{TSelf}.Root(TSelf, int)" />
public static double Root(double x, int n)
{
double result;

if (n > 0)
{
if (n == 2)
{
result = (x != 0.0) ? Sqrt(x) : 0.0;
}
else if (n == 3)
{
result = Cbrt(x);
}
else
{
result = PositiveN(x, n);
}
}
else if (n < 0)
{
result = NegativeN(x, n);
}
else
{
Debug.Assert(n == 0);
result = NaN;
}

return result;

static double PositiveN(double x, int n)
{
double result;

if (IsFinite(x))
{
if (x != 0)
{
if ((x > 0) || IsOddInteger(n))
{
result = Pow(Abs(x), 1.0 / n);
result = CopySign(result, x);
}
else
{
result = NaN;
}
}
else if (IsEvenInteger(n))
{
result = 0.0;
}
else
{
result = CopySign(0.0, x);
}
}
else if (IsNaN(x))
{
result = NaN;
}
else if (x > 0)
{
Debug.Assert(IsPositiveInfinity(x));
result = PositiveInfinity;
}
else
{
Debug.Assert(IsNegativeInfinity(x));
result = int.IsOddInteger(n) ? NegativeInfinity : NaN;
}

return result;
}

static double NegativeN(double x, int n)
{
double result;

if (IsFinite(x))
{
if (x != 0)
{
if ((x > 0) || IsOddInteger(n))
{
result = Pow(Abs(x), 1.0 / n);
result = CopySign(result, x);
}
else
{
result = NaN;
}
}
else if (IsEvenInteger(n))
{
result = PositiveInfinity;
}
else
{
result = CopySign(PositiveInfinity, x);
}
}
else if (IsNaN(x))
{
result = NaN;
}
else if (x > 0)
{
Debug.Assert(IsPositiveInfinity(x));
result = 0.0;
}
else
{
Debug.Assert(IsNegativeInfinity(x));
result = int.IsOddInteger(n) ? -0.0 : NaN;
}

return result;
}
}

/// <inheritdoc cref="IRootFunctions{TSelf}.Sqrt(TSelf)" />
public static double Sqrt(double x) => Math.Sqrt(x);

// /// <inheritdoc cref="IRootFunctions{TSelf}.Root(TSelf, TSelf)" />
// public static double Root(double x, double n) => Math.Root(x, n);

//
// ISignedNumber
//
Expand Down
12 changes: 6 additions & 6 deletions src/libraries/System.Private.CoreLib/src/System/Half.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1359,7 +1359,7 @@ bool IFloatingPoint<Half>.TryWriteSignificandLittleEndian(Span<byte> destination
public static Half Clamp(Half value, Half min, Half max) => (Half)Math.Clamp((float)value, (float)min, (float)max);

/// <inheritdoc cref="INumber{TSelf}.CopySign(TSelf, TSelf)" />
public static Half CopySign(Half x, Half y) => (Half)MathF.CopySign((float)x, (float)y);
public static Half CopySign(Half value, Half sign) => (Half)MathF.CopySign((float)value, (float)sign);

/// <inheritdoc cref="INumber{TSelf}.Max(TSelf, TSelf)" />
public static Half Max(Half x, Half y) => (Half)MathF.Max((float)x, (float)y);
Expand Down Expand Up @@ -1788,15 +1788,15 @@ private static bool TryConvertTo<TOther>(Half value, [NotNullWhen(true)] out TOt
/// <inheritdoc cref="IRootFunctions{TSelf}.Cbrt(TSelf)" />
public static Half Cbrt(Half x) => (Half)MathF.Cbrt((float)x);

// /// <inheritdoc cref="IRootFunctions{TSelf}.Hypot(TSelf, TSelf)" />
// public static Half Hypot(Half x, Half y) => (Half)MathF.Hypot((float)x, (float)y);
/// <inheritdoc cref="IRootFunctions{TSelf}.Hypot(TSelf, TSelf)" />
public static Half Hypot(Half x, Half y) => (Half)float.Hypot((float)x, (float)y);

/// <inheritdoc cref="IRootFunctions{TSelf}.Root(TSelf, int)" />
public static Half Root(Half x, int n) => (Half)float.Root((float)x, n);

/// <inheritdoc cref="IRootFunctions{TSelf}.Sqrt(TSelf)" />
public static Half Sqrt(Half x) => (Half)MathF.Sqrt((float)x);

// /// <inheritdoc cref="IRootFunctions{TSelf}.Root(TSelf, TSelf)" />
// public static Half Root(Half x, Half n) => (Half)MathF.Root((float)x, (float)n);

//
// ISignedNumber
//
Expand Down
Loading

0 comments on commit 715c71a

Please sign in to comment.