我怎样能在C#中快速准确地将一个64位整数乘以一个64位分数?

3
有很多类似的问题在SO上问过,但我还没有找到一个适用于C#且易于移植的解决方案。大多数涉及C++或类似语言,而(可能)有效的答案依赖于嵌入式汇编或本地C/C++函数,这些函数在C#中不存在。一些功能可以在部分范围内工作,但在其他部分失败。我找到了一个可行的答案,并将其移植到了C#,但速度非常慢(原来当我编译为x64时速度还不错,所以我把它发表为答案)。

问题

在我的情况下,我有一个随机的64位Int64/UInt64(使用xoshiro256p算法,尽管这可能与问题无关)。我想将该数字缩放到类型允许范围内的任意范围。例如,我可能希望将Int64缩放到范围[1000, 35000]。从概念上讲,这很容易:

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;
正如许多其他人和上面的评论所指出的那样,问题在于randInt * diff 可能会溢出。 在纸面上,我可以将中间结果存储在128位整数中,然后将除法的结果存储在64位输出中。但是,128位数学对64位系统来说不是本地支持的,而且我宁愿避免使用任意精度库,因为我将要频繁调用此函数,效率将会受到影响。 我可以通过乘以双精度数来获得53位精度,这对于我目前的工作来说足够了,但我宁愿找到一个适当的解决方案。 我可以创建一个带有ASM解决方案的C++库并调用该库,但我想要一个纯C#的解决方案。

要求

  • 需要是纯C#。
  • 需要适用于任何一组输入,使得randInt * diff / maxInt在范围[0,maxInt]内(且每个值本身也在同一范围内)。
  • 不应要求外部库。
  • 需要与数学上正确的答案相差+-1。
  • 需要相当快。也许我只是在寻求奇迹,但我觉得如果双精度可以做到5-10毫秒,我们应该能够用专门构建的代码达到20毫秒,并获得另外11位精度。
  • 最好在发布和调试模式下都能相对良好地工作。我的代码大约有3:1的比率,所以我认为我们可以将调试时间缩短至不到5倍的发布时间。

我的测试

我已经测试了以下解决方案的相对性能。每个测试运行了100万次我的随机数生成器,使用各种方法进行缩放。我首先生成随机数并将它们放入列表中(一个用于有符号数,一个用于无符号数)。然后我遍历每个列表并将其缩放到第二个列表中。

我最初在调试模式下进行了许多测试。这大多数情况下并不重要(我们测试相对性能),但是Int128 / UInt128库在发布模式下表现得好得多。

括号中的数字为调试时间。我在这里包含它们,因为我仍然希望在调试时获得良好的性能。例如,Int128库非常适合发布模式,但调试时很糟糕。在你准备好进行最终发布之前,使用某些具有更好平衡的东西可能会很有用。因为我正在测试一百万个样本,所以毫秒时间也是每个操作的纳秒时间(所有一百万个UInt64均在33ms内生成,因此每个均在33ns内生成)。

我的测试源代码可以在这里,在GitGub上找到。

  • 86毫秒(267):Int64随机数生成器。
  • 33毫秒(80):UInt64随机数生成器。
  • 4毫秒(5):使用双精度转换为Int64,精度降低。
  • 8毫秒(10):再次为UInt64。
  • 76毫秒(197):this C Code适用于Int64,转换为C#(完整代码在我的下面答案中)。
  • 72毫秒(187):再次适用于UInt64。
  • 54毫秒(1458):this UInt128 library适用于Int64。
  • 40毫秒(1476):再次适用于UInt64。
  • 1446毫秒(1455):double128库适用于Int64。商业使用需要付费许可证。
  • 1374毫秒(1397):再次适用于UInt64。

我无法使它们给出正确的结果。

  • 使用DllImport与主应用程序链接的this MulDiv64库
  • 编译为x64的QPFloat,在C++代码中创建了MulDiv64函数。
  • this Java代码
  • 来自Microsoft Media Foundation库的MFllMulDiv函数。我试图测试它,但无法找出如何使VS正确链接到我的C++项目。

类似问题

最准确的方法进行64位组合乘除运算是什么?

  • Phuclv、Soonts、Mysticial和500 - Internal Server Error的回答涉及外部库、汇编或MSVC特定函数。
  • Timos、AnT、Alexey Frunze和Michael Burr的回答实际上没有回答任何问题。
  • Serge Rogatch和Pubby的回答不够精确。
  • AProgrammer的回答有效,但速度非常慢(我不知道它是如何工作的)--我最终仍然使用它,并在x64编译中获得了不错的结果。

当x * n溢出时,如何将x按n / d缩小?

  • 唯一的答案是Abhay Aravinda提供的,但它不是真正的代码,我不确定如何实现最后一部分,而且评论表明它对于大值仍可能发生溢出。

快速方法,在不使用浮点数或溢出的情况下将整数乘以适当的分数

  • Taron和chux - Reinstate Monica的答案是近似值或仅适用于MSVC。
  • R.. GitHub STOP HELPING ICE的答案只是使用64位数学,因为该问题涉及Int32的乘法。

(a * b) / c MulDiv以及处理中间乘法溢出的方法

  • Jeff Penfold的答案对我没有用(我认为在从Java转换到C#时缺少一些逻辑运算符),而且速度非常慢。
  • greybeard的答案看起来不错,但我不确定如何将其翻译成C#。
  • tohoho和Dave的答案会导致溢出。
  • David Eisenstat的答案需要BigInt库。

如何在C++中将64位整数乘以分数并最小化误差?

  • 所有答案在不同情况下都会溢出。

1
禁用优化进行测量对于这样的代码几乎没有用处。使用Benchmark.Net可能是获取准确值的好方法。您可能希望缩放计时值,以便它们表示单个操作的时间,即将ms更改为ns。 - JonasH
2个回答

阿里云服务器只需要99元/年,新老用户同享,点击查看详情
1
“但是128位数学并不是64位系统的本地功能。” 虽然这基本上是正确的,但是有一个不错的方法可以获得两个64位整数的完整128位乘积:Math.BigMul(适用于.NET 5及更高版本)。 x64有一个相应的具有128位输入的除法,这样一对完全乘以宽除法就可以实现这个“按适当分数缩放整数”的操作(限制是分数不能大于1,否则可能会导致溢出)。但是,C#没有访问宽除法的功能,即使有,它在大多数硬件上也不会非常有效。 但是您也可以直接使用BigMul,因为除数应该从一开始就是2的64次方(而不是2的64次方-1),并且BigMul会自动除以2的64次方。 所以代码变成了:(未经测试)
ulong ignore;
ulong scaled = Math.BigMul(randInt, diff, out ignore);
return scaled + minVal;

对于旧版本的.NET,获取产品的高64位可以像这样完成:
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;
}
不幸的是,这不如 Math.BigMul, 但至少还没有除法。

看起来很有前途,但是我无法获取此版本的BigMul。我安装了所有的Framework 4.8东西,升级了VS2019到最新版本,但是我得到的只是32位的BigMul。BigMul页面说4.8支持它,所以我可能做错了什么。我尝试了一个全新的解决方案,但毫无效果。我还没有尝试过VS2022。 - MichaelS
@MichaelS 请查看上面的文档,Math.BigMul仅适用于.NET Core 2.0而不是.NET 4.8。 - phuclv

0

我通过告诉编译器不使用AnyCpu设置,使用AProgrammer's C code,将时间降至约250毫秒。

在发布模式下,PRNG大约需要5毫秒(我有点怀疑;当我尝试运行PRNG时,我认为它被优化掉了),总时间降至约77毫秒。

我仍然不确定它是如何工作的,但链接的答案说该代码对于十进制支持有一些冗余操作。如果我知道足够的工作原理,我认为可以通过优化掉十进制支持来进一步减少时间。

Int64(有符号)速度稍慢(发布版78 vs 77毫秒,调试版慢约20毫秒),但基本上相同的速度。如果min=Int64.MinValue和max=Int64.MaxValue,则会失败并每次返回min,但对我能想到的其他组合都有效。

有符号数学对于直接缩放不太有用。我只是做了一些适用于我的用例的东西。因此,我进行了一些转换,似乎适用于一般的有符号情况,但可能还可以进行一些优化。

无符号缩放算法,转换为C#。

/// <summary>
/// Returns an accurate, 64-bit result from value * multiplier / divisor without overflow.
/// From https://dev59.com/lWoy5IYBdhLWcg3wKq4S#8757419
/// </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使用无符号版本进行有符号转换。在我的用例中速度较慢(调试时290毫秒对260毫秒,发布时95毫秒对81毫秒),但适用于一般情况。不适用于Int64.MinValue(引发异常:“否定补码数的最小值是无效的。”)。

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 函数返回一个介于最小值和最大值之间(包括最小值和最大值)的无符号 UInt64。缩放直接从早期函数进行。

/// <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 函数返回一个在最小值和最大值之间的有符号 Int64。虽然不容易转换为一般缩放,但在这种情况下比上述方法更快。

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

1
即使与本机.NET类相比,C++代码也很慢,因为它未能使用SIMD操作,例如System.Numerics.Vectors命名空间中或通过System.Runtime.Intrinsics命名空间提供的硬件操作。 - Panagiotis Kanavos
1
此外,调试时间是没有意义的。即使平均执行1M次也无效,因为它总是包括预热、噪声,并且未能考虑内存分配和垃圾回收成本。要获得真正有意义的结果,请使用BenchmarkDotNet。 - Panagiotis Kanavos

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,