Calculating the fixed-point representation of (1 - SQRT(0.5)) to arbitrary levels of precision

Question

I have a constant 0.29289321881345247559915563789515..., which can be calculated using the equation (1 - SQRT(0.5)) and then transformed into fixed-point format in order to be used with various sizes of unsigned integers. See the following table for examples:

Word Size	Q Format	Decimal Value	Binary Value
8-bit	Q0.3	5	101
16-bit	Q0.7	75	1001011
32-bit	Q0.15	19195	100101011111011
64-bit	Q0.31	1257966796	1001010111110110000110011001100
128-bit	Q0.63	5402926248376769403	100101011111011000011001100110000000110001000011001101101111011

The problem is that calculating these values becomes infeasible due to the fact that it involves the square root function and numbers that can get arbitrarily large. I remember learning about a "spigot algorithm" that outputs the digits of π using a bounded amount of memory and was hoping that something similar exists for constants such as the one described above. This paper came up during my search, but I have not been able to grok it well enough to translate into code.

How can one accomplish this in a C-like language (preferably C#)? Is there a better way to accomplish the goal of calculating this value for word sizes that are a power of two?

---

Extra Context

I have the following C# snippet that is from another question:

y -= uint.CreateChecked(value: BinaryIntegerConstants<T>.Size) switch {
    8U => (x * ((y * T.CreateChecked(value: 5UL)) >> 4)),
    16U => (x * ((y * T.CreateChecked(value: 75UL)) >> 8)),
    32U => (x * ((y * T.CreateChecked(value: 19195UL)) >> 16)),
    64U => (x * ((y * T.CreateChecked(value: 1257966796UL)) >> 32)),
    128U => (x * ((y * T.CreateChecked(value: 5402926248376769403UL)) >> 64)),
    _ => throw new NotSupportedException(), // TODO: Research a way to calculate the proper constant at runtime.
}

I have manually calculated the constants above by using C# code, but cannot go any further due to the limits of the double type.

Answer 1

I reached out to Daniel Lemire and he was not only gracious enough to help me work out the following solution but wrote a blog post about it too.

public static class BinaryIntegerExtensions
{
    /// <returns>
    /// The result of (1 - SQRT(0.5)) in the fixed-point format Q0.((X >> 1) - 1) where X is the number of bits in type T.
    /// </returns>
    public static T ComputeOneMinusSquareRootOfOneHalf<T>() where T : IBinaryInteger<T> {
        if (typeof(T) == typeof(byte)) { return T.CreateChecked(value: 5); }
        if (typeof(T) == typeof(ushort)) { return T.CreateChecked(value: 75); }

        var @base = (T.One << 8);
        var length = ((T.PopCount(value: T.AllBitsSet) >> 4) + T.One);
        var power = (T.One << 1);
        var result = T.Zero;

        do {
            if (power < ((result << 1) + T.One)) {
                power *= (@base * @base);
                result *= @base;
            }
            else {
                power = ((power - (result << 1)) - T.One);
                result += T.One;
            }
        } while (result < @base.Exponentiate(exponent: length));

        return ((@base.Exponentiate(exponent: (length - T.One)) - ((result / @base) >> 1)) - T.One);
    }
    public static T Exponentiate<T>(this T value, T exponent) where T : IBinaryInteger<T> {
        var result = T.One;

        do {
            if (T.IsOddInteger(value: exponent)) {
                result *= value;
            }

            exponent >>= 1;
            value *= value;
        } while (T.Zero < exponent);

        return result;
    }
}

Example:

BinaryIntegerExtensions.ComputeOneMinusSquareRootOfOneHalf<byte>();    // 5
BinaryIntegerExtensions.ComputeOneMinusSquareRootOfOneHalf<ushort>();  // 75
BinaryIntegerExtensions.ComputeOneMinusSquareRootOfOneHalf<uint>();    // 19195
BinaryIntegerExtensions.ComputeOneMinusSquareRootOfOneHalf<ulong>();   // 1257966796
BinaryIntegerExtensions.ComputeOneMinusSquareRootOfOneHalf<UInt128>(); // 5402926248376769403