在 for(int i=0;...) { for(int j=0; ...){summation}} 嵌套 for 循环中求和势能不起作用

Question

这是重现我得到的错误的简单代码：

#include <math.h> 
#include <iostream>
//#include <omp.h>
//handling Not a number exception:
#include <fenv.h>
#include <signal.h>
#include "unistd.h"

void handler(int sig)
{
  printf("Floating Point Exception\n");
  exit(0);
}
#define EKCOR
const float alpha=200.0/137;
const int N=4096;//4096;//8192;//16384;
const float md=940;
const float Ep=0.1f;
float E1;
int STEP=1;
struct float3
{
  float x, y, z;
};
float3 Pi;
struct Particle
{
  float x;
  float y;
  float z;
  float t;
  float vx;
  float vy;
  float vz;
  float m;
};
Particle p[N] __attribute__((aligned(64)));
inline float3 RandomDirection()
{
  float number1 = static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
  float z   = 2.0*number1 - 1.;  
  float rho = sqrtf((1.+z)*(1.-z));
  float number2 = static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
  float phi = M_PI*2.0*number2;
  float3 result={rho*cosf(phi), rho*sinf(phi), z};
  return result;
}
void function()
{
  float K=0.0;
  Pi={0.0, 0.0, 0.0};
  double Px=0.0;
  double Py=0.0;
  double Pz=0.0;
  float P3=0.0;
  float P4=0.0;
  //#1
  for(int i=0; i<N; ++i)
  {
    Px+=p[i].vx*p[i].m;
    Py+=p[i].vy*p[i].m;
    Pz+=p[i].vz*p[i].m;
    float Penergy=0.0;
  #pragma novector
    for(int j=0; j<N; ++j) if(i!=j)
    {
      float rdist1=sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x)+(p[i].y-p[j].y)*(p[i].y-p[j].y)+(p[i].z-p[j].z)*(p[i].z-p[j].z));
      Penergy+=alpha/rdist1;
      P4+=alpha/rdist1;
    }
    P3+=Penergy;
    float v2=p[i].vx*p[i].vx+p[i].vy*p[i].vy+p[i].vz*p[i].vz;
    K+=p[i].m*v2/2;
  }
  P4/=2;
  Pi.x=Px;
  Pi.y=Py;
  Pi.z=Pz;
  P3/=2;
  float E2=K+P3;
  float r=(E1-P3)/K;
  std::cout<<"r="<<r<<",E1="<<E1<<",P3="<<P3<<",K="<<K<<std::endl;
  float rc=sqrt(r);
  std::cout<<"E2="<<K+P3<<",K="<<K<<",P3="<<P3<<",P4="<<P4<<",Px="<<Pi.x<<",Py="<<Pi.y<<",Pz="<<Pi.z<<std::endl;
}
void init()
{
  const double pi=3.1415926536;   
  float RADIUS=pow(50.0*N,1.0/3.0);
  Pi={0.0, 0.0, 0.0};
  double Px=0.0, Py=0.0, Pz=0.0;
#pragma omp single
  for(int i=0; i<N; ++i)
  {
    float DISTANCE=0.0f;
    if(i>0)
    {
      while(DISTANCE<=1.0f)
      {
        float theta=pi*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
        float phi=2*pi*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
        float rr=RADIUS*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));       
        p[i].x =rr*sin(theta)*cos(phi);
        p[i].y =rr*sin(theta)*sin(phi);
        p[i].z =rr*cos(theta);
        DISTANCE=10000.0f;
      #pragma simd reduction(min:DISTANCE)     
        for(int j=0; j<i; ++j)
        {
          float dij=sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x)+(p[i].y-p[j].y)*(p[i].y-p[j].y)+(p[i].z-p[j].z)*(p[i].z-p[j].z));
          if(dij<DISTANCE) DISTANCE=dij;
        }
      }
    }
    else
    {
      float theta=pi*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
      float phi=2*pi*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));
      float rr=RADIUS*static_cast <float> (rand())/(static_cast <float> (RAND_MAX));       
      p[i].x =rr*sin(theta)*cos(phi);
      p[i].y =rr*sin(theta)*sin(phi);
      p[i].z =rr*cos(theta);
    }
    float modv=sqrt(2.0*Ep/md);
    float3 res=RandomDirection();
    float3 v;
    v.x=modv*res.x;
    v.y=modv*res.y;
    v.z=modv*res.z; 
    p[i].vx =v.x;
    p[i].vy =v.y;
    p[i].vz =v.z;
    p[i].m=md;
    Px+=p[i].vx*p[i].m;
    Py+=p[i].vy*p[i].m;
    Pz+=p[i].vz*p[i].m;   
  }
  Px/=N;
  Py/=N;
  Pz/=N;
#pragma novector
  for(int i=0; i<N; ++i)
  {
    p[i].vx-=Px/p[i].m;
    p[i].vy-=Py/p[i].m;
    p[i].vz-=Pz/p[i].m;
  }
  Px=0.0, Py=0.0, Pz=0.0;
  float K1=0.0;
  float P1=0.0;
  float P2=0.0;
  //#2
#pragma novector
  for(int i=0; i<N; ++i)
  {
    Px+=p[i].vx*p[i].m;
    Py+=p[i].vy*p[i].m;
    Pz+=p[i].vz*p[i].m;
    K1+=p[i].vx*p[i].vx+p[i].vy*p[i].vy+p[i].vz*p[i].vz;
    float pp=0.0;
  #pragma novector
    for(int j=0; j<N; ++j) if(i!=j)
    {
       float rd=sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x)+(p[i].y-p[j].y)*(p[i].y-p[j].y)+(p[i].z-p[j].z)*(p[i].z-p[j].z));
       P1+=alpha/rd;
       pp+=alpha/rd;
    }
    P2+=pp;
  }
  Pi.x=Px;
  Pi.y=Py;
  Pi.z=Pz;
  K1*=md/2;
  P1/=2;
  P2/=2;
  E1=K1+P1;
  std::cout<<"INIT Px="<<Pi.x<<" Py="<<Pi.y<<" Pz="<<Pi.z<<" K1="<<K1<<" P1="<<P1<<" P2="<<P2<<" E1="<<E1<<std::endl;
}

int
main(int argc, char **argv)
{
  //handling Not a number exception:
  feenableexcept(FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW | FE_UNDERFLOW);
  signal(SIGFPE, handler);
  //
  init();
  function();
 }

在 N

N=256 P1=3492.48 P2=3492.5 P3=3492.5 P4=3492.48

但在 N=1024 和 N>1024 时，差异变得越来越大：

N=1024 P1=34968.6 P2=34969.7 P3=34969.7 P4=34968.6
N=2048 P1=114493 P2=114482 P3=114482 P4=114493
N=4096 P1=357880 P2=362032 r=-9.14142

这里程序崩溃了，因为 r=-9.14142 和 sqrt(r) 抛出浮点异常。

我的操作系统是 Fedora 23，处理器是 Intel Core i7-3770，我使用了编译器 gcc 版本 5.3.1 和 intel c++ 编译器 icpc 版本 17.0.1，如果是 必要的。即使不使用 OpenMP，他们都给出了错误。

故障描述在代码下方。 我的问题是：

1. 为什么在 N>=1024 时 P1 与 P2 不同，而 P3 与 P4 的差异如此之大（可以使用 Intel (icpc) 或 gcc (g++) 编译器不带参数进行编译）？该程序在 1 个线程中运行。 它们的值必须相同！

2. 我需要编写代码，以便嵌套的 for 循环 #1 和 #2 并行使用

   #pragma omp parallel for reduction(+:P) for(int i=0; i(p[i].x-p[j].x)+(p[i].y-p[j].y)(p[i].y-p[j].y)+( p[i].z-p[j].z)*(p[i].z-p[j].z)); PP+=阿尔法/r； } P+=PP; } P/=2；

   并使用所有优化标志（我使用集合 -DCMAKE_CXX_FLAGS="-march=native -mtune=native -ipo16 -fp-model fast=2 -O3 -qopt-report=5 -mcmodel=large" 用于英特尔编译器）。 我做不到（即使只有“-O0”）。可能是因为 1) 错误，它给了我一个 错误的值。

Answer 1

您可能还对 randomascii 上的 why floating point arithmetic don't usually do what you think it should do 系列感兴趣。这只是一篇文章的引述，该文章探讨了为什么计算机在浮点（类数学）计算中不精确的问题：

浮点数学不精确。 像0.1这样的简单值无法用二进制浮点数精确表示，浮点数的精度有限意味着运算顺序的细微变化 strong> 或中间体的精度可以改变结果。这意味着比较两个浮点数以查看它们是否相等通常不是您想要的。

(...)

以下是可能出现的不精确性的一个示例：

float f = 0.1f;
float sum;
sum = 0;

for (int i = 0; i < 10; ++i)
    sum += f;
float product = f * 10;
printf("sum = %1.15f, mul = %1.15f, mul2 = %1.15f\n",
        sum, product, f * 10);

此代码尝试以三种不同的方式计算“一”：重复加法和两种轻微的乘法变体。自然我们会得到三个不同的结果，其中只有一个是 1.0：

sum=1.000000119209290, mul=1.000000000000000,  mul2=1.000000014901161

(...)

以下是 0.1、float(0.1) 和 double(0.1) 的确切值：

==================================================== ========================= |号码 |价值 | |------------|------------------------------------ ----------------------| | 0.1 | 0.1（当然）| |浮动 0.1 | 0.100000001490116119384765625 | |双0.1 | 0.10000000000000000055511151231257827021181583404541015625 | ==================================================== =========================

解决了这个问题，我们来看看上面代码的结果：

1. sum = 1.000000119209290：这个计算从一个四舍五入的值开始，然后将其相加十次，每次相加都可能进行四舍五入，因此有很大的误差空间。最终结果不是 1.0，也不是10 * 浮点数（0.1）。但是它是 1.0 以上的下一个可表示的浮点数，因此非常接近。
2. mul = 1.000000000000000：这个计算从一个四舍五入的值开始，然后乘以十，因此出现错误的机会更少。事实证明，从 0.1 到 float(0.1) 的转换向上取整，但乘法在这种情况下，十点恰好是向下舍入，有时两轮是正确的。 因此，我们会因为错误的原因而得到正确的答案。或者可能是错误的答案，因为它实际上不是 float(0.1) 的十倍！
3. mul2 = 1.000000014901161：此计算从一个舍入值开始，然后将 double 精度乘以十，从而避免任何后续舍入错误。所以我们得到了一个不同的正确答案——10 * float(0.1) 的确切值（可以存储在 double 中，但不能存储在 float 中） .

所以，答案之一是不正确的，但它只有一个 float 的距离。答案二是正确的（但不准确），而答案三是完全正确的（但似乎是错误的）。

强调和标记是我的。 randomascii 上的帖子甚至为这个不精确性问题提出了一些可能的解决方案，但他们并没有解决问题（他们只是将不精确性转移到浮点数线的不同部分）。

因此，在处理浮点运算时，您永远不会得到精确的结果。但是您可以采取一些措施来提高计算的精度：

1. 增加浮点中的有效位数。 C++ 的 floats 有 21 个有效位（大约 7 个有效数字），doubles 有 52 个有效位 （大约 ~17 个有效数字）
2. 减少涉及的计算次数（因此4.0*c 比c+c+c+c 更精确）
3. 尽量保证你会以完全相同的顺序进行完全相同的计算（只有这样你才能==/!=这两个值并得到一个合理的结果）

因此，例如，如果您将代码 floats（7 位精度）更改为 doubles（17 位精度），您将看到您的结果变得更加准确并且 显示更多数字。如果您尝试在代码中使用并行化，您的计算可能（也可能不会，取决于实现）在不同的线程/内核上以不同的顺序发生，从而导致所涉及的每个数字的浮点精度大不相同。

作为一个例子，这里是使用double而不是floats的randomascii代码：

  double f = 0.1;
  double sum;
  sum = 0;

  for (int i = 0; i < 10; ++i)
      sum += f;
  double product = f * 10;
  printf("sum = %1.15f, mul = %1.15f, mul2 = %1.15f\n",
          sum, product, f * 10);

哪些输出：

  sum = 1.000000000000000, mul = 1.000000000000000, mul2 = 1.000000000000000

这似乎是正确的，但是当您将 printf 的精度从 1.15f 提高到 1.17f 时：

  sum = 0.99999999999999989, mul = 1.00000000000000000, mul2 = 1.00000000000000000

再次，您可以看到不精确性已经蔓延。sum 执行了 10 次 + 操作，而 mul 和 mul2 分别执行了一次操作 *，这就是为什么 sum 不精确度大于其他两个的不精确。

如果 17 位精度对您来说还不够，那么您可能会对 C++ 的任意精度解决方案感兴趣。

Definition of BigNum from wikipedia:

在计算机科学中，任意精度算术，也称为 bignum 算术、多精度算术或有时是无限精度算术，表示仅对精度位数有限的数字执行计算由主机系统的可用内存决定。

(...)

任意精度用于算术速度不是限制因素的应用，或需要精确结果具有非常大的数字需要的应用>.

再次强调我的。

Here's a related answer suggesting a BigNum library for C++:

GNU 多精度算术库可以满足您的需求http://gmplib.org/

这是使用 GMP 实现的先前代码（使用 64 位精度或大约 21 位有效数字）：

 // Compile like: g++ question.cpp -o out.tmp -lgmpxx -lgmp
 #include <stdio.h>
 #include <gmpxx.h>

 int main(){
      mpf_class f("0.1", 64);
      mpf_class sum("0", 64);

      for (int i = 0; i < 10; ++i)
          sum += f;
      mpf_class product = f * 10;
      printf("sum = %1.17f, mul = %1.17f, mul2 = %1.17f\n",
             sum.get_d(), product.get_d(), ((mpf_class) (f * 10)).get_d());
 }

哪些输出：

  sum = 0.99999999999999989, mul = 0.99999999999999989, mul2 = 0.99999999999999989

这是以 64 位精度进行计算，然后舍入到 51 位（C++ 的 double）并打印出来的结果。

但是，您可以直接从 GMP 打印值：

 // Compile like: g++ question.cpp -o out.tmp -lgmpxx -lgmp
 #include <stdio.h>
 #include <gmpxx.h>
 #include <string>

 int main(){
      mpf_class f("0.1", 64);
      mpf_class sum("0", 64);

      for (int i = 0; i < 10; ++i)
          sum += f;
      mpf_class product = f * 10;
      long exp = 10;
      int base = 10;
      int digits = 21;
      printf("sum = %s, mul = %s, mul2 = %s\n",
             sum.get_str(exp, base, digits).c_str(),
             product.get_str(exp, base, digits).c_str(),
             ((mpf_class) (f * 10)).get_str(exp, base, digits).c_str());
 }

哪些输出：

      sum = 1, mul = 1, mul2 = 1

这是比double 表示更精确的结果。您可以检查 GMP C++ 接口here 和here。 但是请注意，任意精度库通常比内置的 floats 或 doubles 慢。 好处是，为了提高精度，您只需更改 mpf_class variable(expression, precision);行。

也不要忘记查看 PaulMcKenzie 的建议Stack Overflow: Is floating point math broken?：

问题：

考虑以下代码：

0.1 + 0.2 == 0.3 -> false

0.1 + 0.2 -> 0.30000000000000004

为什么会出现这些错误？

答案：

二进制浮点数学是这样的。在大多数编程语言中，它基于 IEEE 754 标准。 (...) 问题的症结在于，数字以这种格式表示为整数乘以 2 的幂； 有理数（如 0.1，即 1/10）分母不是 2 的幂的数不能精确表示。

程序中的常量0.2 和0.3 也将近似 到它们的真实值。碰巧 最接近的 double 到 0.2 大于 rational 数字 0.2 但 最接近的 double 到 0.3 小于 rational 数字 0.3。 0.1 和 0.2 的总和最终大于 rational 数字 0.3，因此与代码中的常量不一致。

强调和标记是我的。

Answer 2

你可能需要做更多的分析，但我的第一个猜测是你的求和循环引起了问题。

提高精度损失的三个技巧：

1. 通过增加尺寸对项目进行排序 - 如果它们尚未排序，这可能成本太高。
2. Pairwise summation
3. Kahan summation

Answer 3

请注意，即使在理论上，P1 应该等于 P2，P3 应该等于 P4，但这些都是浮点变量。更重要的是，它们是单精度浮点变量。根据计算的顺序，你肯定会得到不同的结果。由于浮点表示的不精确性，每次计算都会累积错误。

请查看并运行以下代码（tst_float.cpp）：

/* g++ -Wall tst_float.cpp -o tst_float && ./tst_float */

#include <stdio.h>

int main()
{
    int ok;
    int i;
    float x;

    x = 0.0;
    for (i = 0; i < 10; ++i) {
        x += 0.1;
    }

    ok = x == 1.0;

    if (ok) {
        printf("ok!\n");
    } else {
        printf("uh-uh?\n");
    }
    printf("x == %10.9f\n", x);

    return 0;
}

我明白了：

$ g++ -Wall tst_float.cpp -o tst_float && ./tst_float
uh-uh?
x == 1.000000119

总而言之，不要将浮点变量视为具有整数变量的精度。