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

Unify hypot and abs of complex values

Browse source 
This ports the fixes from abs(T)(ComplexT) regarding NaN and skipping
potentially inaccuracy inducing mathematical no-ops when one arg is 0 to
hypot. Both functions do the same thing and should be deduplicated.
hypot also has some logic regarding under and overflows, and while I
don't fully understand it, it should probably not be removed for complex
numbers.
This commit is contained in:
Inkrementator 2025-03-17 21:19:44 +01:00
parent e1d9d0208b
commit 91aa5feead
2 changed files with 25 additions and 25 deletions
Download patch file Download diff file Expand all files Collapse all files

20
std/complex.d
View file

@ -767,24 +767,8 @@ if (is(T R == Complex!R))
*/
T abs(T)(Complex!T z) @safe pure nothrow @nogc
{
import std.math.algebraic : sqrt;
import std.math.traits : isInfinity, isNaN;
import std.math : fabs;
import std.algorithm.comparison : max;
if (isNaN(z.re) || isNaN(z.im))
return T.nan;
if (isInfinity(z.re) || isInfinity(z.im))
return T.infinity;
const T maxabs = max(fabs(z.re), fabs(z.im));
if (maxabs == 0.0)
return maxabs;
const T qx = z.re / maxabs;
const T qy = z.im / maxabs;
return maxabs * sqrt(qx * qx + qy * qy);
import std.math.algebraic : hypot;
return hypot(z.re, z.im);
}
///

30
std/math/algebraic.d
View file

@ -294,10 +294,11 @@ real cbrt(real x) @trusted pure nothrow @nogc
* hypot(x, -y) are equivalent.
*
* $(TABLE_SV
* $(TR $(TH x) $(TH y) $(TH hypot(x, y)) $(TH invalid?))
* $(TR $(TD x) $(TD $(PLUSMN)0.0) $(TD |x|) $(TD no))
* $(TR $(TD $(PLUSMNINF)) $(TD y) $(TD +$(INFIN)) $(TD no))
* $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD +$(INFIN)) $(TD no))
* $(TR $(TH x) $(TH y) $(TH hypot(x, y)) $(TH invalid?))
* $(TR $(TD x) $(TD $(PLUSMN)0.0) $(TD |x|) $(TD no))
* $(TR $(TD $(PLUSMNINF)) $(TD y) $(TD +$(INFIN)) $(TD no))
* $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD +$(INFIN)) $(TD no))
* $(TR $(TD $(NAN)) $(TD y != $(PLUSMNINF)) $(TD $(NAN)) $(TD no))
* )
*/
T hypot(T)(const T x, const T y) @safe pure nothrow @nogc
@ -308,13 +309,14 @@ if (isFloatingPoint!T)
// If both are huge, avoid overflow by scaling by 2^^-N.
// If both are tiny, avoid underflow by scaling by 2^^N.
import core.math : fabs, sqrt;
import std.math.traits : floatTraits, RealFormat;
import std.math.traits : floatTraits, RealFormat, isNaN;
import std.algorithm.comparison : max;
alias F = floatTraits!T;
T u = fabs(x);
T v = fabs(y);
if (!(u >= v)) // check for NaN as well.
if (!(u >= v))
{
v = u;
u = fabs(y);
@ -322,6 +324,15 @@ if (isFloatingPoint!T)
if (v == T.infinity) return v; // hypot(nan, inf) == inf
}
if (u.isNaN || v.isNaN)
return T.nan;
const maxabs = max(u,v);
if (v == 0.0)
{
return u;
}
static if (F.realFormat == RealFormat.ieeeSingle)
{
enum SQRTMIN = 0x1p-60f;
@ -375,18 +386,23 @@ if (isFloatingPoint!T)
}
// both are in the normal range
return ratio * sqrt(u*u + v*v);
return ratio * sqrt(u^^2 + v^^2);
}
///
@safe unittest
{
import std.math.operations : feqrel;
import std.math.traits : isNaN;
assert(hypot(1.0, 1.0).feqrel(1.4142) > 16);
assert(hypot(3.0, 4.0).feqrel(5.0) > 16);
assert(hypot(real.infinity, 1.0L) == real.infinity);
assert(hypot(1.0L, real.infinity) == real.infinity);
assert(hypot(real.infinity, real.nan) == real.infinity);
assert(hypot(real.nan, real.infinity) == real.infinity);
assert(hypot(real.nan, 1.0L).isNaN);
assert(hypot(1.0L, real.nan).isNaN);
}
@safe unittest