标签: numerical-computing

为什么乘法的顺序会影响结果？请考虑以下代码

a=47.215419672114173
b=-0.45000000000000007
c=-0.91006620964286644
result1=a*b*c
temp=b*c
result2=a*temp
result1==result2

我们都知道,RESULT1 应该是平等的,以RESULT2,但是我们得到:

result1==result2 #FALSE!

差异很小

result1-result2 #3.552713678800501e-15

但是,对于特定应用程序,此错误可能会放大,以便执行相同计算的两个程序的输出(一个使用result1,另一个使用result2)可能完全不同.

为什么会如此以及如何在数字/科学应用中解决这些问题呢？

谢谢!

UPDATE

很好的答案,但我还是很怀念的原因,为什么在顺序的乘法问题,如

temp2=a*b
result3=temp2*c
result1==result3 #True

所以似乎编译器/解释器将*b*c视为(a*b)*c

我已经完成了从Fortran到C++的端口,但发现了COMPLEX类型的一些差异.请考虑以下代码:

PROGRAM CMPLX
    COMPLEX*16 c
    REAL*8 a
    c = (1.23456789, 3.45678901)
    a = AIMAG(1.0 / c)
    WRITE (*, *) a
END

而C++:

#include <complex>
#include <iostream>
#include <iomanip>

int main()
{
    std::complex<double> c(1.23456789, 3.45678901);
    double a = (1.0 / c).imag();

    std::cout << std::setprecision(15) << " " << a << std::endl;
}

使用clang ++或g ++编译C++版本,我得到输出:-0.256561150444368编译Fortran版本但是给了我:-0.25656115049876993

我的意思是,这两种语言都不符合IEEE 754吗？如果我在Octave(Matlab)中运行以下代码:

octave:1> c=1.23456789+ 3.45678901i
c =  1.2346 + 3.4568i
octave:2> c
c =  1.2346 + 3.4568i
octave:3> output_precision(15)
octave:4> c
c =  1.23456789000000e+00 + 3.45678901000000e+00i …

我发现了这个:http://en.cppreference.com/w/cpp/numeric/math/isinf但它似乎检查正无穷大.我只想检查一个值是否恰好等于负无穷大,或者换句话说是log(0)

谢谢你的回答!根据下面的响应,这里有一些代码显示了什么有效.

#include <iostream>
#include <cmath>
#include <math.h>
using namespace std;
int main()
{
    double c = std::log(0.0);
    auto result = c == - INFINITY;
    cout << result << endl;
    return 0;
}

我使用c ++来计算各种类型的特殊函数(例如Lambert函数,用于评估反演的迭代方法等).在许多情况下,直接使用尾数和指数有明显更好的方法.

我找到了许多答案如何提取尾数和指数部分,但是所有这些都只是"学术案例,计算速度不是很快",对我来说有点无用(我用尾数和指数操作的动机是改进计算速度).有时我需要调用一些特定功能大约十亿次(非常昂贵的计算),所以每个节省的计算工作都很好.并且使用将尾数作为double的"frexp"不太合适.

我的问题是(对于带有IEEE 754浮点的c ++编译器):

1)如何读取浮点数/双精度的尾数的具体位？

2)如何将整个尾数读入float/double的整数/字节？

3)与指数1),2)相同的问题.

4)与写入1),2),3)相同的问题.

如果我直接使用尾数或指数,我的动机是更快的计算.我想必须有一个非常简单的解决方案.

在MATLAB中,如果我这样做,realmax - 1000000 == realmax我得到一个逻辑1(真)作为答案.这是为什么？

我正在尝试编写一段代码，用于在周期性边界条件下以二阶数值计算向量场的旋度。然而，我制作的算法非常慢，我想知道是否有人知道任何替代算法。

给出更具体的上下文：我使用 3xAxBxC numpy 数组作为向量场，其中第一个轴指的是笛卡尔方向 (x,y,z)，A、B、C 指的是该笛卡尔方向上的 bin 数量（即决议）。例如，我可能有一个向量场 F = np.zeros((3,64,64,64))，其中 Fx = F[0] 本身就是一个 64x64x64 笛卡尔晶格。到目前为止，我的解决方案是使用 3 点中心差分模板来计算导数，并使用嵌套循环来使用模块化算术迭代所有不同的维度，以强制执行周期性边界条件（例如见下文）。然而，随着我的分辨率增加（A、B、C 的大小），这开始需要很长时间（最多 2 分钟，如果我为模拟执行数百次，这就会增加 - 这只是一个小部分）更大的算法）。我想知道是否有人知道执行此操作的替代方法？

import numpy as np

F =np.array([np.ones([128,128,128]),2*np.ones([128,128,128]),
               3*np.ones([128,128,128])])


VxF =np.array([np.zeros([128,128,128]),np.zeros([128,128,128]),
               np.zeros([128,128,128])])

for i in range(0,128):
     for j in range(0,128):
          for k in range(0,128):
           VxF[0][i,j,k] = 0.5*((F[2][i,(j+1)%128,k]-
                    F[2][i,j-1,k])-(F[1][i,j,(k+1)%128]-F[1][i,j,k-1]))
           VxF[1][i,j,k] = 0.5*((F[0][i,j,(k+1)%128]-
                    F[0][i,j,k-1])-(F[2][(i+1)%128,j,k]-F[2][i-1,j,k]))
           VxF[2][i,j,k] = 0.5*((F[1][(i+1)%128,j,k]-
                    F[1][i-1,j,k])-(F[0][i,(j+1)%128,k]-F[0][i,j-1,k]))

重申一下，我正在寻找一种算法，该算法可以比我现有的算法更快地在给定周期性边界条件的情况下将向量场数组的旋度计算为二阶。也许没有什么可以做到这一点，但我只想在继续花时间运行这个算法之前检查一下。感谢。大家提前！

我正在使用以下算法进行双精度除法，并尝试使其在浮点软件模拟中正确舍入。设a为被除数，b为除数。

所有操作均在 Q2.62 中执行。

倒数的初始近似值为。

b/2是b的有效数，加上其隐含位，并右移一位。接下来，当写成a或b时，它的意思是a或b的有效数加上其隐含位。

这近似为0x17504f333f9de6( 0x5D413CCCFE779800Q2.62 中的)。

之后，用 Newton-Raphson 迭代来近似倒数：

对于倒数r有 6 次这样的迭代。商q是通过将r乘以a（的有效数）来计算的。

商的附加调整步骤：

最终的舍入结果为：

if a <= (a - q * b/2):
      result = final_biased_exponent | q
else
      result = final_biased_exponent | adjusted_q

除以下两种情况外，此方法可以正常工作：a）结果低于正常值或 b）a和b均低于正常值。在这些情况下，它不会正确舍入，并且结果会偏离 1 位（与 x86 结果相比）。（数字a和b被归一化，并且当a或b中的任何一个被归一化时，指数也会相应地缩放。） …

IT++是一个流行的 GPL 库（建议是 LGPL，但目前没有迹象）用于科学计算，尤其是在电信领域。这是同事推荐给我的。

在我的 Linux/Windows 系统上安装这个库看起来很难。和安装指南似乎很难。

除了花时间在这个安装上之外，我还有很多事情要做。任何人都可以通过一些快速安装来节省我的时间，如果可能的话，一些你好，世界！输入说明？

我需要sqrt(1 + (x/2)^2) + x/2数值计算，对于 positive x。对于非常大的 值，使用此表达式直接失败x。我怎样才能重写它以获得更准确的评估？

由于0.0 / 0.0在数学上未定义，IEEE-754 浮点标准合理地定义了NaN其结果。现在，因为与 不同infinity，NaN不是一个明确定义的值，而是一组值，因此是否0.0 / 0.0是一个明确定义的常量的问题也是合理的。

值得一提的是x / 0.0，infinity如果x != 0.0.

IEEE-754 浮点标准中是否有0.0 / 0.0明确定义的常量值？NaN换句话说，它是否具有明确定义的位模式？

我的代码应该计算序列的第 100 个元素$x_0=1 ; x_i=\dfrac{x_{i-1}+1}{x_{i-1}+2}, i=1,2, \ldots$

我写了迭代和递归函数，但结果不相等。是因为丢失了小数吗？

这是我的驱动程序代码。文件中的数据是 i=100。

int main()
{
    int i;

    ifstream f ("data.txt");
    f >> i;

    double x_0= 1.00;

    double x_100 = l(x_0, i);

    ofstream g ("results.txt", ios::app);
    g <<"\n100th element (by looping): " << x_100;

    x_100 = r(x_0);
    g <<"\n100th element (by recursion): " << x_100;

 
   return 0;
}

l()是迭代函数， r()是递归函数

double l(double x, int i)
{
    for (int j = 0; j<i ; j++){
            x = (x + 1)/(x+2);
    } …