分割大笔款项;处理无穷大的总和-解网

问：

我正在运行一种生成对的算法，其中和都是正 s。它们不是独立的; 取决于。算法完成后，我需要将 all 的总和除以 all 的总和。(t[i], x[i])t[i]x[i]double(t[i], x[i])(t[i-1], x[i-1])ax[i]bt[i]

这是故意的，因为它会趋向于无穷大，因为这种计算背后的数学告诉我，在这种情况下，趋向于我实际想要计算的值。binfinityia/b

但是，当我尝试实现此算法时，这是有问题的。我怎样才能防止在某个时候实际上等于无穷大？我想我需要做的是在算法完成之前以某种方式除以的除法。例如，如果我们能写成，那么我可以除以，设置为零，进一步求和，最后除以。但是，我们再次需要确保它不会在某个时候得到无穷大，并且我没有简单的方法来知道我应该如何提出因式分解。babb = b[1] * b[2]ab[1]b[2]b[2]b[2]b = b[1] * b[2]

另一种解决方案可能是某种“巨大的浮点”类型。我愿意接受任何建议，但我需要提到的是，性能在我的应用程序中至关重要。

C++ 数学浮点 C++20 除法

假设（t[i]， x[i]）取决于（t[i−1]， x[i−1]），那么最终结果 sum（t[i]）/sum（x[i]）纯粹是初始值（t[0]， x[0]）的函数。那么，是否有必要计算总和，或者可以通过分析求解呢？还是减少？您是否想要任何任意依赖关系的解决方案（x[i]， t[i]） = f（（x[i−1]， t[i−1]）或只是一个特定的 f？对于起始值还是仅针对特定值？您需要多少个术语？t[i] 和 x[i] 是整数吗？您需要结果的准确性如何？举一些例子。

3赞 lastchance 8/21/2023

我们能看到这些系列的定义吗？

答：

2赞 chux - Reinstate Monica 8/20/2023 #1

将 2 个大值相加可能会溢出到无穷大。

因此，求和必须在加法之前测试大小。

检查任一加法操作数是否接近无穷大，并根据需要缩放值。

// Pseudo code
power_of_2_max = DBL_MAX_EXP - 1; // Max power of 2 before addition
scale_power_of_2 = 0;
sum = 0.0;
for (i=0; i<n; i++) {
  xi = x[i] / pow(2, scale_power_of_2);
  xi_power_of_2 = get_power_of_2(xi);
  if (xi_power_of_2 > power_of_2_max) {
    xi /= 2;
    sum /= 2;
    scale_power_of_2++;
  }
  sum_power_of_2 = get_power_of_2(sum);
  if (sum_power_of_2 > power_of_2_max) {
    xi /= 2;
    sum /= 2;
    scale_power_of_2++;
  }
  sum += xi;
}

// True sum is sum * pow(2, scale_power_of_2)

更深

将许多术语相加会导致不精确。@Eric Postpischil
这是否重要取决于未说明的 OP 目标。

缓解措施：

使用范围更广/更精确的浮点类型进行求和。
对值进行排序，并将值从最小到最大求和。

下面是未经测试的代码，以更好地说明。关键是总和和按 2 的幂缩小，这通常是一个精确的操作，比使用 .x[i]log()

由于的缩放比例至少为 1 in 和，因此加法不会溢出到无穷大。x[i]sum.d += ldexp(x[i], -sum.expo)sum.d <= DBL_MAX/2

#include <assert.h>
#include <float.h>
#include <math.h>

static_assert(FLT_RADIX != 10, "TBD code for base 10");

typedef struct {
  double d;
  int expo;
} double_extra_range;

// Return the sum of the array x[] in .d normalized to about 1.0
// with the extended exponent in .expo
// Designed to work with positive numbers.
double_extra_range sum_double(size_t n, const double *x) {
  double_extra_range sum = {.d = 0.0, .expo = 1};
  for (size_t i = 0; i < n; i++) {
    if (sum.d > DBL_MAX / 2) {
      sum.d /= 2;  
      sum.expo++;
    }
    sum.d += ldexp(x[i], -sum.expo);
  }
  // Normalize to about 1.0
  int p;
  sum.d = frexp(sum.d, &p);
  sum.expo += p;
  return sum;
}

double divide_large_sum(size_t n, const double *t, const double *x) {
  double_extra_range t_sum = sum_double(n, t);
  double_extra_range x_sum = sum_double(n, x);
  double quotinet = x_sum.d / t_sum.d;
  return ldexp(quotinet, x_sum.expo - t_sum.expo);
}

@TobySpeight，我想OP已经尝试过了。不幸的是，有时与 10 字节相同或仅与 10 字节相同，而不是慷慨的 16 字节。OP 具有“性能至关重要”，因此可能不符合该标准。AFAIK，其嵌入式并导致较慢的吞吐量。long doubledoublelong doublelong double

上一个：何时使用浮点运算计算 ceiling-integer-log2 失败？

下一个：Python ROUND_HALF_EVEN的工作原理

分割大笔款项;处理无穷大的总和

Division of large sums; deal with summation up to infinity

评论

评论