왜 두 가지 유사한 부동 소수점 계산이 두 가지 다른 결과를 제공합니까?

Question

다음 코드는 고유 벡터 을 컨테이너로 사용하거나 또는 간단한 C- 배열을 사용하여 동일한 계산을 실현합니다. 그것은 닫히지 만 비트 대 비트가 아닌 동일한 결과를 생성합니다.

최종 수학 연산은 x * alpha + y * beta입니다.

#include <Eigen/Eigen> 

int main() 
{ 
    Eigen::VectorXd x(2); 
    double* y = new double[2]; 
    long long int a = 4603016991731078785; 
    double ga = *(double*)(&a); 
    long long int b = -4617595986472363966; 
    double gb = *(double*)(&b); 
    long long int x0 = 451; 
    long long int x1 = -9223372036854775100; 
    x[0] = *(double*)(&x0); 
    y[0] = *(double*)(&x0); 
    x[1] = *(double*)(&x1); 
    y[1] = *(double*)(&x1); 
    double r = ga*x[0] + gb*x[1]; 
    double s = ga*y[0] + gb*y[1]; 
} 

왜 그렇습니까?

MSVC 및 gcc (64 비트 OS)를 사용하면 결과가 달라집니다.

Answer 1

이것은 한 계산이 80 비트 정밀도로 FPU (부동 소수점 단위) 내에서 완전하게 수행되는 반면 다른 계산은 부분적으로 64 비트 정밀도 (double 크기)를 사용하기 때문일 수 있습니다. 이것은 Eigen을 사용하지 않고도 시연 할 수 있습니다. 다음과 같은 프로그램 봐 : 당신은 최적화없이 컴파일하면

int main() 
{ 
    // Load ga, gb, y[0], y[1] as in original program 
    double* y = new double[2]; 
    long long int a = 4603016991731078785; 
    double ga = *(double*)(&a); 
    long long int b = -4617595986472363966; 
    double gb = *(double*)(&b); 
    long long int x0 = 451; 
    long long int x1 = -9223372036854775100; 
    y[0] = *(double*)(&x0); 
    y[1] = *(double*)(&x1); 

    // Compute s as in original program 
    double s = ga*y[0] + gb*y[1]; 

    // Same computation, but in steps 
    double r1 = ga*y[0]; 
    double r2 = gb*y[1]; 
    double r = r1+r2; 
} 

, 당신은 R s는 (적어도, 내 컴퓨터에서 봤어요) 다른 값을 볼 수 있습니다. 어셈블리 코드를 살펴보면 첫 번째 계산에서 ga, y [0], gb 및 y [1]의 값이 FPU에로드되고 계산 ga * y [0] + gb * y [1]가 완료되면 결과가 메모리에 저장됩니다. FPU는 모든 계산을 80 비트로 수행하지만 결과가 메모리에 저장되면 숫자가 반올림되어 이중 변수의 64 비트 내에 맞습니다.

두 번째 계산은 다르게 진행됩니다. 먼저, ga와 y [0]은 FPU에로드되고 곱해진 다음 64 비트 숫자로 반올림되어 메모리에 저장됩니다. 그런 다음, gb와 y [1]이 FPU에로드되고 곱해진 다음 64 비트 숫자로 반올림되어 메모리에 저장됩니다. 마지막으로 r1과 r2가 FPU에로드되고 추가되어 64 비트 숫자로 반올림되어 메모리에 저장됩니다. 이번에는 컴퓨터가 중간 결과를 반올림하므로 차이가 발생합니다.

비정규 숫자로 작업하기 때문에 반올림은 상당히 큰 효과가 있습니다.

자, 여기서 확실하지 않은 부분이 있습니다. (그리고 이것이 여러분의 질문이라면 사과드립니다.) 이것이 x가 Eigen 컨테이너 인 원래 프로그램과 어떤 관련이 있습니까? 여기에서 계산은 다음과 같습니다 : Eigen의 함수가 x [0]을 얻기 위해 호출되면, 그 함수의 결과와 ga가 FPU에로드되고 곱 해져서 임시 메모리 위치 (64 비트이므로이 값은 반올림). 그런 다음 gb와 x [1]은 FPU에로드되고 곱 해져서 임시 메모리 위치에 저장된 중간 결과에 추가되고 마지막으로 x에 저장됩니다. 따라서 원래 프로그램에서 r을 계산할 때 ga * x [0]의 결과는 64 비트로 반올림됩니다. 아마도 이것은 부동 소수점 스택이 함수 호출에서 보존되지 않기 때문일 수 있습니다.