高效地将双精度数字舍入到给定位数的较低精度的方法

7
在C#中,我想将double类型的数字向下舍入到较低的精度,以便将它们存储在大小不同的桶中,并放入关联数组中。与通常的四舍五入不同,我想要保留指定数量的有效位数进行舍入。因此,大的数字在绝对意义上的变化比小的数字更大,但它们倾向于按相同的比例变化。因此,如果我想将数字舍入到10个二进制位,则找到最重要的十个位,并将所有较低的位清零,可能会添加一个小的数字来进行向上舍入。
我倾向于将"中间值"向上舍入。
如果这是整数类型,以下是可能的算法:
  1. Find: zero-based index of the most significant binary digit set H.
  2. Compute: B = H - P, 
       where P is the number of significant digits of precision to round
       and B is the binary digit to start rounding, where B = 0 is the ones place, 
       B = 1 is the twos place, etc. 
  3. Add: x = x + 2^B 
       This will force a carry if necessary (we round halfway values up).
  4. Zero out: x = x mod 2^(B+1). 
       This clears the B place and all lower digits.
问题在于找到一种高效的方法来找到最高位设定。如果我使用整数,有很酷的位操作技巧可以找到MSB。如果可以避免,我不想调用Round(Log2(x))。该函数将被调用数百万次。
注意:我已阅读了这个SO问题: 什么是将双精度值舍入为(有些)较低精度的好方法? 它适用于C++。我正在使用C#。
更新:
这是代码(修改自答案提供者提供的代码),我正在使用它:
/// <summary>
/// Round numbers to a specified number of significant binary digits.
/// 
/// For example, to 3 places, numbers from zero to seven are unchanged, because they only require 3 binary digits,
/// but larger numbers lose precision:
/// 
///      8    1000 => 1000   8
///      9    1001 => 1010  10
///     10    1010 => 1010  10
///     11    1011 => 1100  12
///     12    1100 => 1100  12
///     13    1101 => 1110  14
///     14    1110 => 1110  14
///     15    1111 =>10000  16
///     16   10000 =>10000  16
///     
/// This is different from rounding in that we are specifying the place where rounding occurs as the distance to the right
/// in binary digits from the highest bit set, not the distance to the left from the zero bit.
/// </summary>
/// <param name="d">Number to be rounded.</param>
/// <param name="digits">Number of binary digits of precision to preserve. </param>
public static double AdjustPrecision(this double d, int digits)
{
    // TODO: Not sure if this will work for both normalized and denormalized doubles. Needs more research.
    var shift = 53 - digits; // IEEE 754 doubles have 53 bits of significand, but one bit is "implied" and not stored.
    ulong significandMask = (0xffffffffffffffffUL >> shift) << shift;
    var local_d = d;
    unsafe
    {
        // double -> fixed point (sorta)
        ulong toLong = *(ulong*)(&local_d);
        // mask off your least-sig bits
        var modLong = toLong & significandMask;
        // fixed point -> float (sorta)
        local_d = *(double*)(&modLong);
    }
    return local_d;
}

更新2:Dekker算法

我从Dekker算法中推导出了这个算法,感谢其他回答者。它会四舍五入到最接近的值,而不是像上面的代码那样截断,并且只使用安全的代码:

private static double[] PowersOfTwoPlusOne;

static NumericalAlgorithms()
{
    PowersOfTwoPlusOne = new double[54];
    for (var i = 0; i < PowersOfTwoPlusOne.Length; i++)
    {
        if (i == 0)
            PowersOfTwoPlusOne[i] = 1; // Special case.
        else
        {
            long two_to_i_plus_one = (1L << i) + 1L;
            PowersOfTwoPlusOne[i] = (double)two_to_i_plus_one;
        }
    }
}

public static double AdjustPrecisionSafely(this double d, int digits)
{
    double t = d * PowersOfTwoPlusOne[53 - digits];
    double adjusted = t - (t - d);
    return adjusted;
}

更新2:时间

我进行了一项测试,发现Dekker算法比TWICE快了一倍以上!

测试中的调用次数:100,000,000
不安全时间 = 1.922(秒)
安全时间 = 0.799(秒)


1
阅读了如何将双精度值舍入到(稍微)较低的精度?之后,您在编写解决方案方面有什么想法。为什么不尝试一下所读内容,如果需要重构或优化,那么请在这里报告结果和/或错误。除非您已经有代码,否则请将其与原始问题一起发布。 - MethodMan
你可能需要具体说明实现你所链接的SO帖子时遇到了什么问题。 - Ani
当然。另一个帖子是关于C++的,并调用了ldexp和frexp函数,对于这些函数我不知道在C#中有什么相似的函数。 - Paul Chernoch
Scalb(y, n)和Logb(x)不就是scalb ilogb的等效物吗?如果是,它们具有frexp和ldexp的功能,只是API略有不同... - aka.nice
2个回答

8
Dekker算法将浮点数分成高位和低位两部分。如果尾数中有 s 位(在IEEE 754 64位二进制中为53位),那么*x0会接收高位的-b位,这是您请求的内容,*x1会接收其余的位,可以丢弃。在下面的代码中,Scale应该具有2b的值。如果b在编译时已知,例如常量43,则可以用0x1p43替换Scale。否则,必须以某种方式生成2b
这需要最近舍入模式。IEEE 754算术可以胜任,但其他合理的算术也可能可以。它将四舍五入的情况舍入到偶数,而不是您请求的向上舍入。这是必要的吗?
这假定x * (Scale + 1)不会溢出。操作必须以双精度(而不是更大)进行评估。
void Split(double *x0, double *x1, double x)
{
    double d = x * (Scale + 1);
    double t = d - x;
    *x0 = d - t;
    *x1 = x - *x0;
}

1
它的工作方式如下:首先,我们想要 x*Scale,因为这有一个高位,我们希望在正确的位置进行后续舍入。但实际上,我们想要一个稍微更大的值,使我们能够减去 x 并仍然具有相同的高位。因此,我们计算 x*Scale+x,即 x*(Scale+1)。那就是 d。接下来,观察到 d-(d-x) 如果精确计算,则恰好是 x。但是,d-x 是以浮点数执行的,因此会四舍五入。d 大于 x,因此 x 的低位在有效数字以下。因此,d-(d-x) 给我们带来了舍去了这些低位的 x - Eric Postpischil
问题在于还有另一个Decker算法。https://en.wikipedia.org/wiki/Dekker%27s_algorithm - Z boson
@EricPostpischil,除了溢出问题,我找不到失败的案例。请原谅我,我要吃一块“谦卑的馅饼”。 - chux - Reinstate Monica
@KevinJin:如答案所述,所示代码假定 x * (Scale + 1) 不会溢出。为了处理接近格式上限的情况,您可以简单地缩小比例,例如通过乘以 2 ** -100,应用算法,然后再放大,例如通过乘以 2 ** 100。请注意,如果该数字使其在所选位置上舍入导致其超出有限范围,则正确结果仍可能是无穷大。 - Eric Postpischil
我明白了。可能有一个常数阈值,首先需要将数字缩小,然后再放大。我会看看能否找到这个阈值,以便添加一个条件分支,使得该解决方案适用于整个可表示范围。 - Kevin Jin
显示剩余9条评论

2

有趣...从未听说过这样的需求,但我认为您可以通过一些奇怪而不安全的代码来“完成它”...

void Main()
{
    // how many bits you want "saved"
    var maxBits = 20;

    // create a mask like 0x1111000 where # of 1's == maxBits
    var shift = (sizeof(int) * 8) - maxBits;
    var maxBitsMask = (0xffffffff >> shift) << shift;

    // some floats
    var floats = new []{ 1.04125f, 2.19412347f, 3.1415926f};
    foreach (var f in floats)
    {
        var localf = f;
        unsafe
        {
            // float -> fixed point (sorta)
            int toInt = *(int*)(&localf);
            // mask off your least-sig bits
            var modInt = toInt & maxBitsMask;
            // fixed point -> float (sorta)
            localf = *(float*)(&modInt);
        }
        Console.WriteLine("Was {0}, now {1}", f, localf);
    }
}

并且使用双倍:

void Main()
{
    var maxBits = 50;
    var shift = (sizeof(long) * 8) - maxBits;
    var maxBitsMask = (0xffffffffffffffff >> shift) << shift;
    var doubles = new []{ 1412.04125, 22.19412347, 3.1415926};
    foreach (var d in doubles)
    {
        var local = d;
        unsafe
        {
            var toLong = *(ulong*)(&local);
            var modLong = toLong & maxBitsMask;
            local = *(double*)(&modLong);
        }
        Console.WriteLine("Was {0}, now {1}", d, local);
    }
}

啊哦...我被拒绝了。 :)

为了完整性,这里使用Jeppe的“无害释放”方法:

void Main()
{
    var maxBits = 50;
    var shift = (sizeof(long) * 8) - maxBits;
    var maxBitsMask = (long)((0xffffffffffffffff >> shift) << shift);
    var doubles = new []{ 1412.04125, 22.19412347, 3.1415926};
    foreach (var d in doubles)
    {
        var local = d;
        var asLong = BitConverter.DoubleToInt64Bits(d);
        var modLong = asLong & maxBitsMask;
        local = BitConverter.Int64BitsToDouble(modLong);
        Console.WriteLine("Was {0}, now {1}", d, local);
    }
}

@PaulChernoch 哦,关于归一化/非归一化,一点头绪都没有。 :) - JerKimball
2
不必使用unsafe代码,您可以使用静态方法BitConverter.DoubleToInt64BitsBitConverter.Int64BitsToDoubledouble(System.Double)和long(System.Int64)之间进行转换。然后,在“安全”的上下文中对long值执行按位AND运算。 - Jeppe Stig Nielsen
@jeppe-stig-nielsen 好的点子。鉴于最初的问题,我很想知道它们在性能方面的比较情况。 - JerKimball
当我有时间时,我会调查使用DoubleToInt64Bits或BitConverter。我展示的实现非常快,但我的老板更喜欢我不使用不安全的代码。(在我的一些测试中,我在对单个保险政策执行卷积时调用它3200万次。通常,应用程序将处理数十万个政策。因此,您可以看到,它需要真正快速。) - Paul Chernoch
@PaulChernoch @JeppeStigNielsen 为了完整性添加了“BitConverter”路由示例(因为其他方法更快而被拒绝时,我短暂地感到了悲伤;)) - JerKimball
显示剩余8条评论

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