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

Answer 2

I don't know of a spigot algorithm for binary or hexadecimal sqrt(1/2) but that doesn't mean one doesn't exist. It isn't my field but I am interested in PLSQ which can handle similar approximation problems.

The BPP algorithm you recall for pi is something of a peculiarity. In that it is possible to compute isolated digits of the hexadecimal representation of pi correctly using a very specific formula obtained by sophisticated computer assisted high precision experimental mathematics. It was a big surprise when it was discovered in 1995. Later papers on PLSQ and LLL still use its derivation as an example of the power of newer faster integer relations and lattice reduction algorithms.

There is however the option of Newton-Raphson available in a form that doubles the number of significant digits each time you use it and can work in scaled integer arbitrary precision form if you are careful about it.

Wanting a precise sqrt(1/2) is particularly helpful here. Multiplies and divides by 2 are trivial in base 2.

I think there may be a way forward for you using a scaled integer version of an old Newton-Raphson based method for computing 1/sqrt(y) where y in this case is 2. The reciprocal square root algorithm works as follows and dates to before there was FP hardware support for doing it quickly. Some cbrts still use this method which neatly avoids doing a slow division. Divisions remain pretty slow today when compared to all of the other binary operations. Below is the simple derivation of rsqrt()

f(x)  = y - 1/x^2
f'(x) = 2/x^3

Let the starting guess be x0 then NR iteration gives

x1 = x0 - f(x0)/f'(x0)

x1 = x0 - (y - 1/x0^2)*x0^3/2  = x0 - x0*(y*x0^2-1)/2

Then set y = 2 (your special case) when the NR formula for simplifies further

x1 = x0 - x0*(x0*x0-1/2)  (#1)

x1 = 1 - 1*(1 - 1/2)    = 1/2
x2 = 1/2 - 1/2*(1/4-1/2) = 5/8
x3 = 3/4 - 5/8*(25/64-1/2) = 355/512

and with this you can double the number of accurate bits every iteration. The computation cost is relatively modest 2 * and 2 - computed at twice the length of the previous result.

But be careful of the last few bits at any fixed depth since although the algorithm is self correcting there could be a carry from the next stage that affects several bits in the least significant part of the present one.

I know that (#1) can be rewritten as but it is less accurate in that form (at least in normal FP work)

 x1 = x0*(3/2-x0*x0)  (#1a)

You could also do it with Halley's method but that requires a non-trivial division at full working precision and so is bound to be slower here.

I would be inclined to precompute these values in an environment that allows easy access to fast arbitrary precision floating point like Julia or Python and import them as needed at runtime in the C# program