How can I quickly and accurately multiply a 64-bit integer by a 64-bit fraction in C#?

Question

There are a lot of similar questions asked on SO, but I've yet to find one that works and is easily portable to C#. Most involve C++ or similar, and the (presumably) working answers rely on either embedded assembly or native C/C++ functions that don't exist in C#. Several functions work for part of the range, but fail at other parts. I found one working answer I was able to port to C#, but it was very slow (turns out it's decently-fast when I compile to x64 instead of x86, so I posted it as the answer to beat).

Problem

I need a function that allows me to multiply any 64-bit integer by a fraction 0 to 1 (or -1 to 1) that is derived from two 64-bit integers. Ideally, the answer would work for both Int64 and UInt64, but it's probably not hard to make one work from the other.

In my case, I have a random 64-bit Int64/UInt64 (using the xoshiro256p algorithm, though that's likely irrelevant). I want to scale that number to any arbitrary range in the type's allowed values. For example, I might want to scale Int64 to the range [1000, 35000]. This is, conceptually, easy enough:

UInt64 minVal = 1000;
UInt64 maxVal = 35000;
UInt64 maxInt = UInt64.MaxValue;

UInt64 randInt = NextUInt64(); // Random value between 0 and maxInt.
UInt64 diff = maxVal - minVal + 1;
UInt64 scaledInt = randInt * diff / maxInt; // This line can overflow.
return scaledInt + minVal;

As noted by many other people, and the comment above, the problem is that randInt * diff can potentially overflow.

On paper, I could simply store that intermediate result in a 128-bit integer, then store the result of the division in the 64-bit output. But 128-bit math isn't native to 64-bit systems, and I'd rather avoid arbitrary-precision libraries since I'll be making lots of calls to this function and efficiency will be notable.

I could multiply by a double to get 53 bits of precision, which is fine for what I'm currently doing, but I'd rather come up with a proper solution.

I could create a C++ library with one of the ASM solutions and call that library, but I'd like something that's pure C#.

Requirements

• Needs to be pure C#.
• Needs to work for any set of inputs such that randInt * diff / maxInt is in the range [0, maxInt] (and each value itself is in the same range).
• Shouldn't require an external library.
• Needs to be +-1 from the mathematically-correct answer.
• Needs to be reasonably quick. Maybe I'm just asking for miracles, but I feel like if doubles can do 5-10 ms, we should be able to hit 20 ms with purpose-built code that gets another 11 bits of precision.
• Ideally works relatively well in both release and debug modes. My code has about a 3:1 ratio, so I'd think we could get debug under 5-ish times the release time.

My Testing

I've tested the following solutions for relative performance. Each test ran 1 million iterations of my random number generator, scaling using various methods. I started by generating random numbers and putting them in lists (one for signed, one for unsigned). Then I ran through each list and scaled it into a second list.

I initially had a bunch of tests in debug mode. It mostly didn't matter (we're testing relative performance), but the Int128/UInt128 libraries fared much better in release mode.

Numbers in parenthesis are the debug time. I include them here because I still want decent performance while debugging. The Int128 library, for example, is great for release mode, but terrible for debug. It might be useful to use something that has a better balance until you're ready for final release. Because I'm testing a million samples, the time in milliseconds is also the time in nanoseconds per operation (all million UInt64s get generated in 33 ms, so each one is generated in 33 ns).

Source code for my testing can be found here, on GitGub.

• 86 ms (267): the Int64 random generator.
• 33 ms (80): the UInt64 random generator.
• 4 ms (5): using double conversion to Int64, with reduced precision.
• 8 ms (10): again for UInt64.
• 76 ms (197): this C Code for Int64, converted to C# (exact code in my answer below).
• 72 ms (187): again for UInt64.
• 54 ms (1458): this UInt128 library, for Int64.
• 40 ms (1476): again for UInt64.
• 1446 ms (1455): double128 library for Int64. Requires a paid license for commercial use.
• 1374 ms (1397): again for UInt64.

I couldn't get these to give proper results.

• this MulDiv64 library, linked to the main application with DllImport.
• QPFloat, compiled to x64, created a MulDiv64 function in the C++ code.
• this Java code.
• the MFllMulDiv function from the Microsoft Media Foundation library. I tried to test it, but couldn't figure out how to get VS to link into my C++ project properly.

Similar Questions

Most accurate way to do a combined multiply-and-divide operation in 64-bit?

• Answers by phuclv, Soonts, Mysticial, and 500 - Internal Server Error involve external libraries, assembly, or MSVC-specific functions.
• Answers by timos, AnT, Alexey Frunze, and Michael Burr don't actually answer anything.
• Answers by Serge Rogatch and Pubby aren't precise.
• Answer by AProgrammer works, but is very slow (and I have no idea how it works) -- I ended up using it anyways and getting decent results in x64 compilation.

How can I descale x by n/d, when x*n overflows?

• The only answer, by Abhay Aravinda, isn't real code, I wasn't sure how to implement the last section, and the comments suggest it can overflow for large values anyways.

Fast method to multiply integer by proper fraction without floats or overflow

• Answers by Taron and chux - Reinstate Monica are approximations or MSVC-specific.
• Answer by R.. GitHub STOP HELPING ICE just uses 64-bit math since that question is about multiplying Int32.

(a * b) / c MulDiv and dealing with overflow from intermediate multiplication

• Answer by Jeff Penfold didn't work for me (I think I'm missing something in the logical operators converting from Java to C#), and it was very slow.
• Answer by greybeard looks nice, but I wasn't sure how to translate it to C#.
• Answers by tohoho and Dave overflow.
• Answer by David Eisenstat requires BigInt libraries.

How to multiply a 64 bit integer by a fraction in C++ while minimizing error?

• All the answers overflow in different circumstances.

Answer 1

But 128-bit math isn't native to 64-bit systems

While that is mostly true, there is a decent way to get the full 128-bit product of two 64-bit integers: Math.BigMul (for .NET 5 and later)

x64 has a corresponding division with a 128-bit input, and such a pair of full-multiply followed by a wide-division would implement this "scale integer by a proper fraction" operation (with the limitation that the fraction must not be greater than 1, otherwise an overflow could result). However, C# doesn't have access to wide division, and even if it did, it wouldn't be very efficient on most hardware.

But you can just use BigMul directly too, because the divisor should really be 2⁶⁴ to begin with (not 2⁶⁴ - 1), and BigMul automatically divides by 2⁶⁴.

So the code becomes: (not tested)

ulong ignore;
ulong scaled = Math.BigMul(randInt, diff, out ignore);
return scaled + minVal;

---

For older versions of .NET, getting the high 64 bits of the product could be done like this:

static ulong High64BitsOfProduct(ulong a, ulong b)
{
    // decompose into 32bit blocks (in ulong to avoid casts later)
    ulong al = (uint)a;
    ulong ah = a >> 32;
    ulong bl = (uint)b;
    ulong bh = b >> 32;
    // low times low and high times high
    ulong l = al * bl;
    ulong h = ah * bh;
    // cross terms
    ulong x1 = al * bh;
    ulong x2 = ah * bl;
    // carry from low half of product into high half
    ulong carry = ((l >> 32) + (uint)x1 + (uint)x2) >> 32;
    // add up all the parts
    return h + (x1 >> 32) + (x2 >> 32) + carry;
}

Unfortunately that's not as good as Math.BigMul, but at least there is still no division.

Answer 2

I was able to get down to about 250 ms using AProgrammer's C code by telling the compiler to NOT prefer 32-bit code using the AnyCpu setup.

In release mode, the PRNG takes up about 5 ms (I somewhat doubt this; I think it's being optimized out when I try to just run the PRNG), and the total is down to about 77ms.

I'm still not sure how it works, but the linked answer says the code has some redundant operations for base 10 support. I'm thinking I can reduce the time even further by optimizing out the base 10 support, if I knew how it worked enough to do that.

The Int64 (signed) is a little slower (78 vs 77ms release, about 20ms slower debug), but I'm basically the same speed. It does fail if min=Int64.MinValue and max=Int64.MaxValue, returning min every time, but works for every other combination I could throw at it.

The signed math is less useful for straight scaling. I just made something that worked in my use case. So I made a conversion that seems to work for the general signed case, but it could probably be optimized a bit.

Unsigned scaling algorithm, converted to C#.

/// <summary>
/// Returns an accurate, 64-bit result from value * multiplier / divisor without overflow.
/// From https://stackoverflow.com/a/8757419/5313933
/// </summary>
/// <param name="value">The starting value.</param>
/// <param name="multiplier">The number to multiply by.</param>
/// <param name="divisor">The number to divide by.</param>
/// <returns>The result of value * multiplier / divisor.</returns>
private UInt64 MulDiv64U(UInt64 value, UInt64 multiplier, UInt64 divisor)
{
    UInt64 baseVal = 1UL << 32;
    UInt64 maxdiv = (baseVal - 1) * baseVal + (baseVal - 1);

    // First get the easy thing
    UInt64 res = (value / divisor) * multiplier + (value % divisor) * (multiplier / divisor);
    value %= divisor;
    multiplier %= divisor;
    // Are we done?
    if (value == 0 || multiplier == 0)
        return res;
    // Is it easy to compute what remain to be added?
    if (divisor < baseVal)
        return res + (value * multiplier / divisor);
    // Now 0 < a < c, 0 < b < c, c >= 1ULL
    // Normalize
    UInt64 norm = maxdiv / divisor;
    divisor *= norm;
    value *= norm;
    // split into 2 digits
    UInt64 ah = value / baseVal, al = value % baseVal;
    UInt64 bh = multiplier / baseVal, bl = multiplier % baseVal;
    UInt64 ch = divisor / baseVal, cl = divisor % baseVal;
    // compute the product
    UInt64 p0 = al * bl;
    UInt64 p1 = p0 / baseVal + al * bh;
    p0 %= baseVal;
    UInt64 p2 = p1 / baseVal + ah * bh;
    p1 = (p1 % baseVal) + ah * bl;
    p2 += p1 / baseVal;
    p1 %= baseVal;
    // p2 holds 2 digits, p1 and p0 one

    // first digit is easy, not null only in case of overflow
    UInt64 q2 = p2 / divisor;
    p2 = p2 % divisor;

    // second digit, estimate
    UInt64 q1 = p2 / ch;
    // and now adjust
    UInt64 rhat = p2 % ch;
    // the loop can be unrolled, it will be executed at most twice for
    // even baseVals -- three times for odd one -- due to the normalisation above
    while (q1 >= baseVal || (rhat < baseVal && q1 * cl > rhat * baseVal + p1))
    {
        q1--;
        rhat += ch;
    }
    // subtract 
    p1 = ((p2 % baseVal) * baseVal + p1) - q1 * cl;
    p2 = (p2 / baseVal * baseVal + p1 / baseVal) - q1 * ch;
    p1 = p1 % baseVal + (p2 % baseVal) * baseVal;

    // now p1 hold 2 digits, p0 one and p2 is to be ignored
    UInt64 q0 = p1 / ch;
    rhat = p1 % ch;
    while (q0 >= baseVal || (rhat < baseVal && q0 * cl > rhat * baseVal + p0))
    {
        q0--;
        rhat += ch;
    }
    // we don't need to do the subtraction (needed only to get the remainder,
    // in which case we have to divide it by norm)
    return res + q0 + q1 * baseVal; // + q2 *baseVal*baseVal
}

MulDiv64 uses the unsigned version to get a signed conversion. It's slower in my use case (290ms vs 260ms debug, 95ms vs 81ms release), but works for the general case. Doesn't work for Int64.MinValue (raises an exception: "Negating the minimum value of a twos complement number is invalid.").

public static Int64 MulDiv64(Int64 value, Int64 multiplier, Int64 divisor)
{
    // Get the signs then convert to positive values.
    bool isPositive = true;
    if (value < 0) isPositive = !isPositive;
    UInt64 val = (UInt64)Math.Abs(value);
    if (multiplier < 0) isPositive = !isPositive;
    UInt64 mult = (UInt64)Math.Abs(multiplier);
    if (divisor < 0) isPositive = !isPositive;
    UInt64 div = (UInt64)Math.Abs(divisor);

    // Scaledown.
    UInt64 scaledVal = MulDiv64U(val, mult, div);

    // Convert to signed Int64.
    Int64 result = (Int64)scaledVal;
    if (!isPositive) result *= -1;

    // Finished.
    return result;
}

GetRangeU function returns an unsigned UInt64 between min and max, inclusive. Scaling is straight from the earlier function.

/// <summary>
/// Returns a random unsigned integer between Min and Max, inclusive.
/// </summary>
/// <param name="min">The minimum value that may be returned.</param>
/// <param name="max">The maximum value that may be returned.</param>
/// <returns>The random value selected by the Fates for your application's immediate needs. Or their fickle whims.</returns>
public UInt64 GetRangeU(UInt64 min, UInt64 max)
{
    // Swap inputs if they're in the wrong order.
    if (min > max)
    {
        UInt64 Temp = min;
        min = max;
        max = Temp;
    }

    // Get a random integer.
    UInt64 randInt = NextUInt64();

    // Fraction randInt/MaxValue needs to be strictly less than 1.
    if (randInt == UInt64.MaxValue) randInt = 0;

    // Get the difference between min and max values.
    UInt64 diff = max - min + 1;

    // Scale randInt from the range 0, maxInt to the range 0, diff.
    randInt = MulDiv64U(diff, randInt, UInt64.MaxValue);

    // Add the minimum value and return the result.
    return randInt;// randInt + min;
}

GetRange function returns a signed Int64 between min and max. Not easily convertible to general scaling, but it's faster than the method above in this case.

/// <summary>
/// Returns a random signed integer between Min and Max, inclusive.
/// Returns min if min is Int64.MinValue and max is Int64.MaxValue.
/// </summary>
/// <param name="min">The minimum value that may be returned.</param>
/// <param name="max">The maximum value that may be returned.</param>
/// <returns>The random value selected.</returns>
public Int64 GetRange(Int64 min, Int64 max)
{
    // Swap inputs if they're in the wrong order.
    if (min > max)
    {
        Int64 Temp = min;
        min = max;
        max = Temp;
    }

    // Get a random integer.
    UInt64 randInt = NextUInt64();

    // Fraction randInt/MaxValue needs to be strictly less than 1.
    if (randInt == UInt64.MaxValue) randInt = 0;

    // Get the difference between min and max values.
    UInt64 diff = (UInt64)(max - min) + 1;

    // Scale randInt from the range 0, maxInt to the range 0, diff.
    randInt = MulDiv64U(diff, randInt, UInt64.MaxValue);

    // Convert to signed Int64.
    UInt64 randRem = randInt % 2;
    randInt /= 2;
    Int64 result = min + (Int64)randInt + (Int64)randInt + (Int64)randRem;

    // Finished.
    return result;
}