c - number - ieee floating point




향상된 정확도로 효율적으로(a-K)/(a+K) 컴퓨팅 (4)

예를 들어 수학 함수에 대한 인수 감소와 같은 다양한 상황에서 a 는 양의 변수 인수이고 K 는 상수 인 (a - K) / (a + K) 를 계산해야합니다. 대부분의 경우, K 는 2의 제곱으로, 제 작업과 관련된 유스 케이스입니다. 나는 직접적 분할로 달성 할 수있는 것보다 더 정확하게이 지수를 계산하는 효율적인 방법을 찾고있다. 이 작업은 현재 모든 주요 CPU 및 GPU 아키텍처에서 제공되며 fma()fmaf() 함수를 통해 C / C ++에서 사용할 수 있으므로 FMA (fused multiply-add)에 대한 하드웨어 지원을 가정 할 수 있습니다.

탐험을 쉽게하기 위해 나는 float 산수를 실험하고 있습니다. double 연산에 대한 접근 방식을 포팅 할 계획이기 때문에 인수와 결과 모두 네이티브 정밀도보다 높은 연산을 사용할 수는 없습니다. 내 최고의 솔루션은 지금까지 :

 /* Compute q = (a - K) / (a + K) with improved accuracy. Variant 1 */
 m = a - K;
 p = a + K;
 r = 1.0f / p;
 q = m * r;
 t = fmaf (q, -2.0f*K, m);
 e = fmaf (q, -m, t);
 q = fmaf (r, e, q);

[K/2, 4.23*K] 간격의 인수 a 의 경우, 위의 코드는 K 가 2의 제곱이고, 모든 입력에 대해 거의 정확하게 반올림 된 몫을 계산합니다 (최대 오류는 0.5 ulps에 매우 가깝습니다). 중간 결과에는 오버플로 또는 언더 플로우가 발생하지 않습니다. 2의 거듭 제곱이 아닌 K 경우이 코드는 나눗셈을 기반으로 한 순진 알고리즘보다 훨씬 정확합니다. 성능면에서이 코드는 부동 소수점 역수가 부동 소수점 분할보다 빠르게 계산 될 수있는 플랫폼에서의 순진 접근보다 빠를 수 있습니다.

나는 K = 2 n 일 때 다음과 같은 관찰을한다 : 작업 간격의 상한선이 8*K , 16*K 증가 할 때, 최대 오차는 점진적으로 증가하고 순진한 계산의 최대 오차는 아래에서부터 천천히 근사하기 시작한다 . 불행하게도, 간격의 하한에 대해서도 같은 것은 나타나지 않습니다. 하한값이 0.25*K 떨어지면, 위의 개선 된 방법의 최대 오차는 순진 기법의 최대 오차와 동일하다.

순진한 방법과 위의 코드 시퀀스에 비해 더 큰 최대 오차 ( ulp 대 수학적 결과로 측정)를 얻을 수있는 q = (a-K) / (a ​​+ K)를 계산하는 방법이 있습니까? , 특히 하한이 0.5*K 미만인 간격에 대해? 효율성은 중요하지만 위의 코드에서 사용 된 것보다 몇 가지 작업이 더 용인 될 수 있습니다.

아래의 한 가지 대답에서 두 개의 피연산자, 즉 머리 꼬리 쌍 q:qlo 의 평가되지 않은 합으로 몫을 반환하여 정확도를 향상시킬 수 있다고 지적했습니다 q:qlo 즉 잘 알려진 double float 및 더블 double 포맷. 위의 코드에서 이것은 qlo = r * e 마지막 줄을 변경하는 것을 의미합니다.

이 접근법은 분명히 유용하며, pow() 에서 사용하기위한 확장 된 대수에 대한 사용을 이미 고려했습니다. 그러나 향상된 계산이 더 정확한 지수를 제공하는 간격의 원하는 확대에 근본적으로 도움이되지 않습니다. 내가 바라는 특정한 경우에, 기본 근사 간격을 좁게 유지하기 위해 K=2 (단 정밀도의 경우) 또는 K=4 (배정 밀도의 경우)를 사용하고 싶습니다. 간격은 대략 [0.28 ]. 내가 직면하고있는 실용적인 문제는 <0.25 * K의 주장에 대해 개선 된 부문의 정확성이 순진한 방법보다 실질적으로 좋지 않다는 것입니다.


a가 K에 비해 크다면 (aK) / (a ​​+ K) = 1 - 2K / (a ​​+ K)는 좋은 근사치를 제공합니다. a가 K에 비해 작 으면 2a / (a ​​+ K) - 1은 좋은 근사치를 제공합니다. K / 2 ≤ a ≤ 2K이면 aK는 정확한 연산이므로 나누기를 수행하면 적절한 결과를 얻을 수 있습니다.


나는 정말로 대답을하지 못했다. (적절한 부동 소수점 오차 분석은 매우 지루하다.) 그러나 몇 가지 관찰들 :

  • 빠른 상호 지침 (예 : RCPSS )은 사업부만큼 정확하지 않으므로 이들을 사용하면 정확성이 저하 될 수 있습니다.
  • m 은 정확히 ε [0.5 × K b , 2 1 + n × K b ], K b 는 K 아래의 2의 거듭 제곱 (또는 K가 2의 거듭 제곱 인 경우 K 자체)이고, n은 (K가 2의 거듭 제곱이면 n = 23).
  • 이는 Dekker (1971)div2 알고리즘의 단순화 된 형식과 유사합니다. 범위를 확장하려면 (특히 하한), 아마도 더 많은 수정 조항을 통합해야 할 것입니다 (즉, m 을 2 float s를 사용하거나 double 사용하십시오).

문제는 (a + K) 의 덧셈입니다. (a + K) 의 모든 정밀도 손실은 나누기로 확대됩니다. 문제는 부서 자체가 아닙니다.

aK 의 지수가 같으면 (거의) 정밀도가 손실되지 않으며 지수의 절대 차가 significand 크기보다 큰 경우 (a + K) == a ( a 가 더 큰 경우) 또는 (a + K) == K ( K 가 더 큰 경우).

이것을 막을 방법이 없습니다. 유효 숫자 크기를 늘리면 (예 : 80x86에서 80 비트 "확장 double"사용) "정확한 결과 범위"를 약간 넓힐 수 있습니다. 이유를 이해하려면 smallest + largest 것을 고려하십시오 (여기서 smallest 것은 32 비트 부동 소수점 숫자가 될 수있는 가장 작은 양의 비정규 값입니다). 이 경우 (32 비트 부동 소수점의 경우) 정밀도 손실을 완전히 피하기 위해 결과에 대해 유효 크기 약 260 비트가 필요합니다. (예 : temp = 1/(a + K); result = a * temp - K / temp; temp = 1/(a + K); result = a * temp - K / temp; (a + K) 문제가 여전히 똑같기 때문에 많은 도움이되지 않습니다 (그러나 (a - K) 에서 비슷한 문제는 피할 것입니다). 또한 result = anything / p + anything_error/p_error 를 할 수 없습니다. 왜냐하면 부서가 그렇게 작동하지 않기 때문입니다.

32 비트 부동 소수점에 적합 할 수 a 모든 가능한 양수 값에 대해 0.5 ulps에 가깝게 생각할 수있는 3 가지 대안 만 있습니다. 아무도 받아들이 기 어려울 것입니다.

첫 번째 대안은 32 비트 부동 소수점의 경우 약 2 GiB가되는 (모든 트릭을 사용하여) 64 비트의 부동 소수점에 대해 완전히 미친, a 모든 값에 대해 룩업 테이블 ( "큰 실수" 비트 부동 소수점). 물론, a 의 가능한 값의 범위가 "32 비트 float에 들어갈 수있는 임의의 양수 값"보다 작은 경우에는 조회 테이블의 크기가 줄어 듭니다.

두 번째 대안은 실행시 계산 (그리고 32 비트 부동 소수점으로 변환 / 변환)에 다른 것을 사용하는 것입니다 ( "큰 실수").

세 번째 대안은 "무언가"(나는 그것이 무엇인지 불리지는 모르지만 값이 비싸다)와 관련이있다. 반올림 모드를 "양의 무한대로 반올림"으로 설정하고 temp1 = (a + K); if(a < K) temp2 = (a - K); 계산합니다 temp1 = (a + K); if(a < K) temp2 = (a - K); temp1 = (a + K); if(a < K) temp2 = (a - K); "음의 무한대로 반올림"으로 전환하고 if(a >= K) temp2 = (a - K); lower_bound = temp2 / temp1; if(a >= K) temp2 = (a - K); lower_bound = temp2 / temp1; . 다음으로 a_lower = a 하고 가능한 한 가장 작은 양만큼 a_lower 를 줄이고 "lower_bound"계산을 반복하고 lower_bound에 대해 다른 값을 얻을 때까지 계속 수행 한 다음 a_lower 의 이전 값으로 a_lower . 그 후 upper_bounda_upper ( a 의 원래 값으로 시작)를 결정하기 위해 근본적으로 같은 (그러나 반대 반올림 모드와 감소하지 않는 증가) 작업을 수행합니다. 마지막으로 a_range = a_upper - a_lower; result = upper_bound * (a_upper - a) / a_range + lower_bound * (a - a_lower) / a_range; 와 같이 보간합니다 a_range = a_upper - a_lower; result = upper_bound * (a_upper - a) / a_range + lower_bound * (a - a_lower) / a_range; a_range = a_upper - a_lower; result = upper_bound * (a_upper - a) / a_range + lower_bound * (a - a_lower) / a_range; . 초기 상한과 하한을 계산하고이 값이 같으면이 모든 것을 건너 뛰기를 원할 것입니다. 또한이 모든 것이 "이론적으로, 완전히 테스트되지 않았다"는 경고를 받으면서 아마 어딘가에서 그것을 골랐다.

주로 내가 말하는 것은 (제 의견으로는) 당신이 포기하고 0.5 ulp에 가깝게하기 위해 할 수있는 일이 없다는 것을 받아 들여야한다는 것입니다. 죄송합니다.. :)


오류를 모델링하는 다른 변수를 반환하기 위해 API를 완화 할 수 있다면 솔루션은 훨씬 간단 해집니다.

float foo(float a, float k, float *res)
{
    float ret=(a-k)/(a+k);
    *res = fmaf(-ret,a+k,a-k)/(a+k);
    return ret;
}

이 솔루션은 나눗셈의 절단 오류 만 처리하지만 a+kak 의 정밀도 손실은 처리하지 않습니다.

이러한 오류를 처리하기 위해 고정 소수점을 사용하기 위해 이중 정밀도 또는 둔감도를 사용해야한다고 생각합니다.

입력에서 0이 아닌 최하위 비트를 인위적으로 생성하도록 테스트 코드가 업데이트됩니다.

테스트 코드

https://ideone.com/bHxAg8





floating-accuracy