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

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 128 ms, we should be able to hit 200 ms with purpose-built code that gets another 11 bits of precision.

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. Target range was [UInt64.MaxValue - 1million, <-thatValue 9999]. I also tested [0, 9999] to make sure it worked at both ends of the spectrum, but didn't time it.

The following values are all done in debug mode. I probably should have done release mode tests, but I didn't think about it until after I'd already deleted most of the code.

Numbers in parenthesis are what's left after removing the PRNG time.

78ms (0): the random generator with no scaling.
128ms (50): using double conversion, with reduced precision.
635ms (557): this C Code, converted to C#.
950ms (900): double128 library. Requires a paid license for commercial use.
1575ms (1497): this UInt128 library.
2000ms: the same UInt128 library when I target x64.

I couldn't get these to give proper results, but got some timing information.

205ms (127): this MulDiv64 library, linked to the main application with DllImport.
1000ms (932): QPFloat, compiled to x64, created a MulDiv64 function in the C code.
1315 ms (1237): this Java code.
NA: 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.

CodePudding user response：

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;
}

CodePudding user response：

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;