2

我正在用 C# 编写一个 BigDecimal 类。我已经成功实现了 +、- 和 * 运算符。但我想不出一种方法来计算 2 BigDecimals 的除法。使用这 3 个运算符实现除法的最快方法是什么？还是有更好的方法来做到这一点？（同时考虑开发时间和算法速度）

目标 1：我希望结果是另一个具有固定精度的 BigDecimal（应该是可变的）

目标 2：正如您所提到的，BigDecimal 的目的不是固定精度。那么我怎样才能达到无限精度呢？

另一个问题：BigRational在此线程中使用 Microsoft BCL 的类进行任意精度算术然后使用 Christopher Currens 的扩展方法是否更好（关于速度和灵活性） ： C# 中是否有 BigFloat 类？获得十进制表示而不是编写一个新类？

4

3 回答 3

10

首先，我假设“大十进制”是指您代表一个理性，其中分母被限制为十的任何幂。

您需要认真考虑您希望两位小数的除法输出是什么。小数在加法、减法和乘法上是封闭的，但在除法上不是封闭的。也就是说，任何两个小数相乘产生第三个：(7 / 10) * (9 / 100) 得到 63 / 1000，这是另一个小数。但是将这两个小数相除，你会得到一个分母中没有十次方的有理数。

要回答您实际提出的问题：就像乘法可以通过循环中的加法构建一样，除法可以通过循环中的减法构建。要将 23 除以 7，请说：

  • 7 < 23 吗？是的。
  • 从 23 中减去 7 得到 16。
  • 7 < 16 吗？是的。
  • 从 16 中减去 7 得到 9
  • 7 < 9 吗？是的。
  • 从 9 中减去 7 得到 2。
  • 7 < 2？不，我们有我们的第一个数字。我们做了三个减法，所以第一个数字是 3。
  • 将 2 乘以 10 得到 20
  • 7 < 20 吗？是的。
  • 从 20 中减去 7 得到 13。
  • 7 < 13 吗？是的。
  • 从 15 中减去 7 得到 6。
  • 7 < 6？不，我们做了两次减法和一次乘以 10，所以下一个数字是 2。
  • 将 6 乘以 10 得到 60
  • 7 < 60 吗？是的...
  • ...
  • 我们做了八次减法，所以下一个数字是 8...
  • ... 等等

你知道为此目的更快的算法吗？

当然，有很多更快的除法算法。这是一个：Goldschmidt 算法。

首先，我希望很清楚，如果您尝试计算，X / D那么您可以先计算1 / D然后将其乘以X. 此外，让我们假设 WOLOG D 严格介于 0 和 1 之间。

如果不是呢？如果 D 是负数，则将它和 X 反转；如果 D 为零，则给出错误；如果 D 为 1，则答案为 X；如果 D 大于 1，则将它和 X 都除以 10，这对您来说应该很容易，因为您的系统是十进制的。继续应用这些规则，直到你的 D 介于 0 和 1 之间。（作为附加优化：当 D 非常小时，算法最慢，所以如果 D 小于 0.1，则将 X 和 D 乘以 10，直到 D 大于或等于 0.1。）

好的，所以我们的问题是我们有一个介于 0 和 1 之间的数字 D，我们希望计算1 / D。可能最简单的事情就是做一个例子。假设我们正在尝试计算1 / 0.7. 正确答案是1.42857142857...

首先从 2 中减去 0.7 得到 1.3。现在将分数的两个部分都乘以 1.3：

(1 / 0.7) * (1.3 / 1.3) = 1.3 / 0.91

伟大的。我们现在已经计算1 / 0.7到一位数的准确性。

现在再做一次。从 2 中减去 0.91 得到 1.09。将分数的两个部分乘以 1.09：

(1.3 / 0.91) * (1.09 / 1.09) = 1.417 / 0.9919

太好了，现在我们有两个正确的数字。现在再做一次。从 2 中减去 0.9919 得到 1.0081。将顶部和底部乘以 1.0081：

(1.417 / 0.9919) * (1.0081 / 1.0081) = 1.4284777 / 0.99993439

嘿，现在我们有四个正确的数字。看看情况如何？分母的每一步都越来越接近 1，因此分子越来越接近1 / 0.7。

这比减法方法收敛得快得多。

你明白它为什么有效吗？

由于 D 介于 0 和 1 之间，因此存在一个数 E 使得 D = 1 - E，并且 E 也在 0 和 1 之间。

当我们将 D 乘以 (2 - D) 时，我们将 (1 - E) 乘以 (1 + E)，得到 1 - E 2。

由于 0 < E < 1，显然 E 2小于 E 并且也在 0 和 1 之间，这意味着 1 - E 2更接近于 1。事实上它更接近 1。通过多次重复这个过程，我们很快就接近 1。实际上，我们在这里所做的是将每一步正确数字的数量大致增加一倍。显然，这比减法方法要好得多，减法方法在每一步都给出一个额外的数字。

继续这样做，直到获得所需的准确性。由于您在每个步骤中将准确数字的数量大致翻了一番，因此您应该能够很快达到可接受的准确度。由于我们一开始就已经安排 D 大于等于 0.1，所以 E 永远不会大于 0.9；反复平方 0.9 会让你很快降到一个非常小的数字。

于 2013-08-29T15:05:20.047 回答
1

Eric Lippert 的回答看起来很像“以最慢的方式进行长除法”。它是准确的，但它并不快。假设您有一种用双精度值逼近 BigDecimal 的方法，并且由于您已经在您的BD类中实现了 +、- 和 *，您可以（大约）执行以下操作：

BD operator/(BD a, BD b) {
  BD r, result=0, residual;
  double aa, bb, rr;
  residual = a;
  bb = b;
  while (residual > desired) {
    rr = (double)residual / bb;
    r = rr; // assuming you have a conversion from double to bigdecimal with loss of precision
    result += r;
    residual = a - (result * b);
  }
  return result;
}

这将进行逐次逼近，但一次有多个数字（基本上您使用 BD 算术找到错误，然后使用双算术将其除掉）。我认为这应该是一种相对直接和有效的方法。

我已经使用floatand实现了上述的一个版本double- 只是为了向自己证明这个原理是可行的。使用简单的 C 框架，只需两次迭代即可将精度降低到double除法级别。我希望你能明白。

#include <stdio.h>

double divide(double a, double b);
int main(void) {
  double a, b, result;
  float fa, fb, fr;
  a = 123.5;
  b = 234.6;
  fa = a;
  fb = b;
  fr = fa / fb;
  printf("using float: %f\n", fr);
  result = divide(a, b);
  printf("using double: %lf\n", result);
  printf("difference: %le\n", result - fr);
}

double divide(double a, double b) {
      double r, result=0, residual;
      float aa, bb, rr;
      residual = a;
      bb = b;
      while (residual > 1e-8) {
        rr = (float)residual / bb;
        r = rr; // assuming you have a conversion from double to bigdecimal with loss of precision
        result += r;
        residual = a - (result * b);
        printf("residual is now %le\n", residual);
      }
      return result;
    }

输出：

using float: 0.526428
residual is now 8.881092e-06
residual is now 5.684342e-14
using double: 0.526428
difference: 3.785632e-08
于 2013-09-05T14:44:06.673 回答
1
  • 您可以查看 Java BigDecimal实现以获取一些想法
  • 对于计算 a/b，您可以使用二分搜索算法来查找 b * c = a 的 c（您应该运行算法，直到达到所需的精度）。
  • 你也可以看看这个BigFloat类，它有有趣的指数形式的实现
  • 对于无限精度，您可以将您存储BigDecimal为二的有理分数BigInteger。

有理分数表示的代数：

(x1/x2) + (y1/y2) = (x1*y2+x2*y1)/(x2*y2)    
(x1/x2) - (y1/y2) = (x1*y2-x2*y1)/(x2*y2)    
(x1/x2) * (y1/y2) = (x1*y1)/(x2*y2)
(x1/x2) / (y1/y2) = (x1*y2)/(x2*y1)

双精度（固定精度）的实现示例：

public double Divide(double a, double b, double eps)
{
    double l = 0, r = a;
    while (r - l > eps)
    {
        double m = (l + r) / 2;
        if (m * b < a)
            l = m;
        else
            r = m;
    }
    return (l + r) / 2;
}
于 2013-08-29T14:14:34.653 回答