Algorithms

Algorightm.Math.cs

[StructLayout(LayoutKind.Explicit)]
readonly ref struct FloatToIntBytesConverter
{
    [FieldOffset(0)]
    public readonly float Float;

    [FieldOffset(0)]
    public readonly int Int;

    public FloatToIntBytesConverter(float x) => Float = x;
    public FloatToIntBytesConverter(int x) => Int = x;
}

const int __InvSqrtFloatMagik = 0x5f3759df;

static double SqrtInvFast(float x)
{
    var i = new FloatToIntBytesConverter(x).Int;
    i = __InvSqrtFloatMagik - (i >> 1);
    var y = new FloatToIntBytesConverter(i).Float;
    y *= 1.5f - 0.5f * x * y * y;

    return y;
}

[StructLayout(LayoutKind.Explicit)]
readonly ref struct DoubleToLongBytesConverter
{
    [FieldOffset(0)]
    public readonly double Double;

    [FieldOffset(0)]
    public readonly long Long;

    public DoubleToLongBytesConverter(double x) => Double = x;
    public DoubleToLongBytesConverter(long x) => Long = x;
}

const long __InvSqrtDoubleMagik = 0x5FE6EB50C7B537A9;

static double SqrtInvFast(this double x)
{
    var i = new DoubleToLongBytesConverter(x).Long;
    i = __InvSqrtDoubleMagik - (i >> 1);
    var y = new DoubleToLongBytesConverter(i).Double;
    y *= 1.5 - 0.5 * x * y * y;

    return y;
}

static double Sqrt(double x, double Eps = 1e-7)
{
    switch (x)
    {
        case < 0:
            throw new InvalidOperationException();
        case 0 or 1:
            return x;
    }

    var y = x / 2;
    double prev;

    do
        (prev, y) = (y, 0.5 * (y + x / y)); 
    while 
        (prev - y > Eps);

    return y;
}