정수 제수를 일정하게 사용하여 효율적인 부동 소수점 나눗셈

정수 제수를 일정하게 사용하여 효율적인 부동 소수점 나눗셈

최근 컴파일러가 부동소수점 나눗셈을 부동소수점 곱셈으로 대체할 수 있는지 여부에 대한 질문에서 이 질문을 하게 되었습니다.

부호 변환 후의 결과가 실제의 분할 연산과 비트 단위로 동일해야 한다는 엄격한 요건에서는, 2진수 IEEE-754 산술의 경우, 2승의 제수가 가능한 것을 알 수 없다.제수의 역수가 표현 가능한 한, 제수의 역수에 곱하면 나눗셈과 동일한 결과가 나옵니다.를 들면, 「」에 의한 곱셈, 「」에 의한 등입니다.0.5할 수 2.0.

그런 다음 이러한 대체가 어떤 다른 약수가 작동하는지 궁금할 수 있습니다. 이러한 약수는 분할을 대체하지만 훨씬 더 빠르게 실행되며 비트와 동일한 결과를 제공한다고 가정합니다.특히 일반 곱셈 외에 퓨전 다중 덧셈 연산을 허용합니다.코멘트에서 나는 다음과 같은 관련 문서를 지적했다.

니콜라스 브리세바레, 장 미셸 뮐러, 사우라브 쿠마르 레이나.제수가 미리 알려진 경우 올바르게 반올림된 부동 소수점 나눗셈을 가속합니다.컴퓨터에 관한 IEEE 거래, 제53권, 제8호, 2004년 8월, 페이지 1069-1072.

이 논문의 저자들이 주장하는 기법은 z = 1 / y_hl, z = fma(-y_h, z, 1) / y와_l 같이 제수 y의 역수를 미리 계산한다. 나중에_h, 나눗셈 q = x / y는 q = fma(z_h, x, z_l * x)로 계산된다.이 논문은 이 알고리즘이 작동하기 위해 제수 y가 충족해야 하는 다양한 조건을 도출합니다.쉽게 관찰할 수 있듯이 이 알고리즘은 머리와 꼬리의 부호가 다를 때 무한과 0에 문제가 있다.더 중요한 것은, 배당이 매우 작은 x에 대해 정확한 결과를 도출하지 못할 것이다. 왜냐하면 몫 꼬리(z * x)의_l 연산이 언더플로우이기 때문이다.

이 백서는 또한 Peter Markstein이 IBM에 있을 때 개척한 대체 FMA 기반 분할 알고리즘에 대해서도 언급하고 있습니다.관련 참고 자료는 다음과 같습니다.

P. W. 마크스타인IBM RISC System/6000 프로세서의 기본 함수 계산.IBM Journal of Research & Development, Vol. 34, No.1, 1990년 1월, 페이지 111-119

마크스타인의 알고리즘에서는 먼저 역수 rc를 계산하고, 여기서 초기 몫 q = x * rc가 형성된다.그런 다음, 나머지 나눗셈은 r = fma(-y, q, x)로 FMA를 사용하여 정확하게 계산하고, 개선되고 더 정확한 지수를 q = fma(r, rc, q)로 계산한다.

이 알고리즘에는 제로 또는 무한대(적절한 조건부 실행으로 쉽게 해결 가능)의 문제도 있지만 IEEE-754를 사용한 철저한 테스트입니다.float데이터는 이 많은 작은 정수 중 많은 제수 y에 대한 모든 가능한 배당 x에 대해 올바른 몫이 전달된다는 것을 보여줍니다.이 C코드에 의해 구현됩니다.

/* precompute reciprocal */
rc = 1.0f / y;

/* compute quotient q=x/y */
q = x * rc;
if ((x != 0) && (!isinf(x))) {
    r = fmaf (-y, q, x);
    q = fmaf (r, rc, q);
}

대부분의 프로세서 아키텍처에서는 이것이 예측, 조건부 이동 또는 선택 유형의 명령을 사용하여 분기 없는 일련의 명령으로 변환됩니다.: 에 대하여3.0f , . . . . . . . .nvcc 7의 머신 를 생성합니다.CUDA 7.5는 GPU의 머신 코드를 합니다.

    LDG.E R5, [R2];                        // load x
    FSETP.NEU.AND P0, PT, |R5|, +INF , PT; // pred0 = fabsf(x) != INF
    FMUL32I R2, R5, 0.3333333432674408;    // q = x * (1.0f/3.0f)
    FSETP.NEU.AND P0, PT, R5, RZ, P0;      // pred0 = (x != 0.0f) && (fabsf(x) != INF)
    FMA R5, R2, -3, R5;                    // r = fmaf (q, -3.0f, x);
    MOV R4, R2                             // q
@P0 FFMA R4, R5, c[0x2][0x0], R2;          // if (pred0) q = fmaf (r, (1.0f/3.0f), q)
    ST.E [R6], R4;                         // store q

나의 실험에서는, 다음의 작은 C 테스트 프로그램을 작성했습니다.정수 나눗셈을 순서대로 하고, 각각의 코드 시퀀스를 적절한 나눗셈에 대해 철저하게 테스트합니다.이 완전 테스트에 합격한 제수의 리스트가 표시됩니다.부분 출력은 다음과 같습니다.

PASS: 1, 2, 3, 4, 5, 7, 8, 9, 11, 13, 15, 16, 17, 19, 21, 23, 25, 27, 29, 31, 32, 33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 64, 65, 67, 69,

치환 알고리즘을 최적화로서 컴파일러에 짜넣기 위해서는, 상기의 코드 변환을 안전하게 적용할 수 있는 제수의 화이트리스트는 실용적이지 않습니다.에서는, 「」( 「1의 가능한 , 하고 있는 을 알 수 .x y2번물론 증거가 아닌 일화적인 증거죠

위의 코드 시퀀스로의 나눗셈 변환이 안전한지 여부를 결정할 수 있는 수학적 조건은 무엇입니까?답변은 모든 부동소수점 연산이 기본 반올림 모드인 "가장 가까운 반올림 또는 짝수"에서 수행된다고 가정할 수 있습니다.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

int main (void)
{
    float r, q, x, y, rc;
    volatile union {
        float f;
        unsigned int i;
    } arg, res, ref;
    int err;

    y = 1.0f;
    printf ("PASS: ");
    while (1) {
        /* precompute reciprocal */
        rc = 1.0f / y;

        arg.i = 0x80000000;
        err = 0;
        do {
            /* do the division, fast */
            x = arg.f;
            q = x * rc;
            if ((x != 0) && (!isinf(x))) {
                r = fmaf (-y, q, x);
                q = fmaf (r, rc, q);
            }
            res.f = q;
            /* compute the reference, slowly */
            ref.f = x / y;

            if (res.i != ref.i) {
                err = 1;
                break;
            }
            arg.i--;
        } while (arg.i != 0x80000000);

        if (!err) printf ("%g, ", y);
        y += 1.0f;
    }
    return EXIT_SUCCESS;
}

세 번째 재기동하겠습니다.우리는 가속화를 시도하고 있다.

    q = x / y

서 ''는y이고, 「정수」입니다.q,x , , , , 입니다.y는 모두 IEEE 754-2008 binary32 부동소수점 값입니다.아래,fmaf(a,b,c) 곱셈을 a * b + c을 사용합니다.

순진한 알고리즘은 미리 계산된 역수를 통해,

    C = 1.0f / y

실행 시 (빠른 속도로) 곱셈으로 다음을 충족합니다.

    q = x * C

Brisebarre-Muller-Raina 가속은 사전 계산된 두 개의 상수를 사용합니다.

    zh = 1.0f / y
    zl = -fmaf(zh, y, -1.0f) / y

따라서 런타임에 하나의 곱셈과 하나의 융합된 곱셈 덧셈으로 충분합니다.

    q = fmaf(x, zh, x * zl)

마크스타인 알고리즘은 순진한 접근방식을 두 개의 퓨전된 다중 애드와 결합하며, 순진한 접근방식이 사전 계산을 통해 최소 유의한 위치에서 1 단위 내에서 결과를 산출할 경우 올바른 결과를 산출한다.

    C1 = 1.0f / y
    C2 = -y

따라서 점괘는 다음과 같이 근사할 수 있습니다.

    t1 = x * C1
    t2 = fmaf(C1, t1, x)
    q  = fmaf(C2, t2, t1)

---

은 두 가지 .y이데올로기 때문에를 들어,7, 30의 , , 수 7, , 14, ,, 28 15, 28 30에 대해 를 얻을 수 .x.

의 Brisebarre-Muller-Raina의 거의 합니다.y, 적은 , , , , , , , , "x 결과(모든 의 0.)를 산출하다x 에 .y를 참조해 주세요.

Brisebarre-Muller-Raina 기사는 순진한 접근법의 최대 오차는 ±1.5 ULP임을 보여준다.

의 Markstein의 한다.y 정수인 y(Markstein 접근법에 대해 실패한 홀수 제수는 발견되지 않았습니다.)

---

Markstein 접근법에 대해서는 1부터 19700까지의 약수를 분석했다(여기서 원시 데이터).

의 수, 값 수)를 그립니다.xMarkstein 접근법이 해당 제수에 대해 실패함)에서는 단순한 패턴이 발생하는 것을 볼 수 있습니다.

_{(출처 : nominal-animal.net )}

이 그림에는 수평축과 수직축이 모두 로그로 표시됩니다.홀수 제수에 대한 점은 없습니다. 테스트한 모든 홀수 제수에 대해 올바른 결과를 얻을 수 있기 때문입니다.

x축을 제수의 비트 역방향(역순으로 16진수, 즉 0b11101101 → 0b101101, 데이터)으로 변경하면 매우 명확한 패턴을 얻을 수 있습니다.
_{(출처 : nominal-animal.net )}

집합의 곡선이 .4194304/x플로트의 플로트를 모두 할 때는 두 8388608/x ★★★★★★★★★★★★★★★★★」2097152/x전체 오류 패턴을 완전히 묶습니다.

을 사용하면, 을 할 수 있습니다.rev(y)의 y , , , 「 」8388608/rev(y)이 아닌 Markstein의 수( 플로트의입니다.y (또는16777216/rev(x)상한을 설정합니다.)

2016-02-28 추가: 임의의 정수(이진수 32) 제수가 주어졌을 때 Markstein 접근방식을 사용하여 에러 건수의 근사치를 구했다.다음은 의사 코드와 같습니다.

function markstein_failure_estimate(divisor):
    if (divisor is zero)
        return no estimate
    if (divisor is not an integer)
        return no estimate

    if (divisor is negative)
        negate divisor

    # Consider, for avoiding underflow cases,
    if (divisor is very large, say 1e+30 or larger)
        return no estimate - do as division

    while (divisor > 16777216)
        divisor = divisor / 2

    if (divisor is a power of two)
        return 0

    if (divisor is odd)
        return 0

    while (divisor is not odd)
        divisor = divisor / 2

    # Use return (1 + 83833608 / divisor) / 2
    # if only nonnegative finite float divisors are counted!
    return 1 + 8388608 / divisor

테스트한 Markstein 고장 사례에서 오차 추정치는 ±1 이내가 됩니다(8388608보다 큰 약수는 아직 테스트하지 않았습니다).마지막 구분은 잘못된 0이 보고되지 않도록 해야 하지만, 저는 장담할 수 없습니다.언더플로우 문제가 있는 매우 큰 약수(0x1p100, 1e+30, 그 이상 크기)는 고려하지 않습니다.어쨌든 이러한 약수는 가속에서 제외됩니다.

예비 테스트에서, 그 추정치는 믿을 수 없을 정도로 정확해 보인다.모든 점이 그림에서 정확히 일치하기 때문에 추정치와 20000의 실제 오차를 비교하는 그림을 그리지 않았습니다(이 범위 내에서 추정치가 정확하거나 너무 큽니다).기본적으로 추정치는 이 답변의 첫 번째 그림을 정확하게 재현합니다.

---

Markstein 접근법의 실패 패턴은 규칙적이고 매우 흥미롭습니다.이 방법은 2개의 제수와 홀수 정수 제수의 모든 거듭제곱에 대해 작동합니다.

16777216보다 큰 제수의 경우 16777216보다 작은 값을 얻기 위해 2의 최소제곱으로 나눈 제수와 동일한 오차가 나타납니다.예를 들어 0x1.3cdfa4p+23 및 0x1.3cdfa4p+41, 0x1.d8874p+32, 0x1.cf84f8p+23 및 0x1.cf84f8p+34, 0x1.1e4a7fp+23 및 0x1.437p+fp+inside입니다.

테스트 벤치에 오류가 없다고 가정하면, 마크스타인 접근방식은 16777216보다 큰 약수(단, 1e+30보다 작은 약수)를 2의 최소 거듭제곱으로 나누면 16777216보다 작은 약수가 되고, 그 몫이 홀수인 경우 또한 작용한다는 것을 의미한다.

질문입니다.Y할 수 입니다.x / Y의 FMA에 대해 FMA를 하여 보다 을 할 수 .x은 정적 하여 값의 x는 변환된 코드가 원래 분할과 다른 값이 발생하지 않는다는 것을 알고 일반적으로 건전한 변환을 적용할 수 있습니다.

부동소수점 계산의 문제에 잘 적응하는 부동소수점 값 집합의 표현을 사용하면 함수의 시작부터 시작하는 순방향 분석에서도 유용한 정보를 생성할 수 있습니다.예:

float f(float z) {
  float x = 1.0f + z;
  float r = x / Y;
  return r;
}

를 전제로, 의 함수 「」에서는 「」(*)」(「Round-to-Nearest 모드)를 사용합니다.x는 NaN(입력값이 NaN인 경우), +0.0f 또는 2보다^-24 큰 수치만 사용할 수 있으며 -0.0f 또는 2보다^-24 0에 가까운 값은 사용할 수 없습니다.이는 상수의 많은 값에 대한 질문에 표시된 두 가지 형태 중 하나로 변환하는 것을 정당화한다.Y.

으로 사용하지 한 많은 하고 C 가 이미 (*)#pragma STDC FENV_ACCESS ON

---

A는 다음 정보를 예측하는 스태틱 분석을 전송합니다.x상기의 표현은, 식이 다음의 튜플로서 취할 수 있는 부동 소수점 값의 세트의 표현에 근거할 수 있습니다.

• 가능한 NaN 값의 집합에 대한 표현(NaN의 동작이 충분히 지정되지 않았기 때문에 선택은 부울만 사용하는 것이다.true즉, 일부 NaN이 존재할 수 있습니다.falseNaN이 존재하지 않음을 나타냅니다.)
• 각각 +inf, -inf, +0.0, -0.0의 존재를 나타내는 4개의 부울플래그
• 음의 유한 부동 소수점 값의 포함 간격
• 양의 유한 부동 소수점 값의 포함 간격.

이 방법을 따르려면 C 프로그램에서 발생할 수 있는 모든 부동소수점 연산을 스태틱아나라이저로 이해해야 합니다.예를 들어, 덧셈은 U와 V의 집합 사이에 있으며, 이 값을 처리하는 데 사용됩니다.+는 다음과 같이 구현될 수 있습니다.

• 피연산자 중 하나에 NaN이 존재하거나 피연산자가 반대 부호의 무한대일 수 있는 경우 결과에 NaN이 존재한다.
• 0이 U 값과 V 값을 더한 결과일 수 없는 경우 표준 간격 산술을 사용합니다.결과의 상한은 U에서 가장 큰 값과 V에서 가장 큰 값을 반올림에서 가장 가까운 값으로 더하기 위해 구해지므로 이러한 한계를 반올림에서 가장 가까운 값으로 계산해야 한다.
• 0이 U의 양의 값과 V의 음의 값을 더한 결과일 경우, M이 V에 존재하도록 M을 U에서 가장 작은 양의 값으로 한다.
  • succ(M)가 U에 있으면 이 값의 쌍은 결과의 양수 값에 succ(M) - M을 기여한다.
  • V에 -succ(M)가 있으면 이 값의 쌍은 결과의 음수 값에 음수 값 M - succ(M)을 기여한다.
  • pred(M)가 U에 있는 경우, 이 값의 쌍은 결과의 음수 값에 pred(M) - M을 기여합니다.
  • -pred(M)가 V에 있는 경우, 이 값의 쌍은 결과의 양의 값에 M - pred(M) 값을 기여합니다.
• 0이 음수 값 U와 양수 값 V를 더한 결과일 경우 동일한 작업을 수행한다.

인식: 상기 아이디어는 "부동소수점 가산 및 감산 제약 개선", Bruno Marre & Claude Michel에서 차용됩니다.

---

예: 함수의 컴파일f이하에 나타냅니다.

float f(float z, float t) {
  float x = 1.0f + z;
  if (x + t == 0.0f) {
    float r = x / 6.0f;
    return r;
  }
  return 0.0f;
}

질문의 접근방식은 기능상의 분담을 변환하는 것을 거부한다.f6은 나눗셈을 무조건 변환할 수 있는 값 중 하나가 아니기 때문에 대체 형식으로 변환합니다.대신, 제가 제안하는 것은 함수의 시작부터 간단한 값 분석을 적용하는 것입니다. 이 경우, 이 경우,x유한 플로트이거나+0.0f또는 최소^-24 2개의 규모, 그리고 이 정보를 사용하여 브리세바레 외 연구진의 변환을 적용하기 위해 다음과 같은 지식에 확신을 갖는다.x * C2는 언더플로하지 않습니다.

명확하게 하기 위해 다음과 같은 알고리즘을 사용하여 나눗셈을 보다 단순한 것으로 변환할지 여부를 결정할 것을 제안합니다.

1. 이Y알고리즘에 따라 Brisebarre 등의 방법을 사용하여 변환할 수 있는 값 중 하나?
2. 그들의 방법의 C1과 C2는 같은 부호를 가지고 있는가, 아니면 배당이 무한할 가능성을 배제할 수 있는가?
3. C1과 C2가 같은 부호를 가지거나x0의 두 표현 중 하나만 선택하시겠습니까?C1과 C2의 부호가 다른 경우 및x0의 표현은 1개뿐입니다.FMA 기반 계산의 부호를 만지작거리고(**) 다음 경우에 올바른 0을 생성하도록 하는 것을 잊지 마십시오.x0 입니다.
4. 배당금의 규모가 다음과 같은 가능성을 배제할 수 있을 만큼 충분히 크다는 것을 보장할 수 있는가?x * C2언더플로우?

네 가지 질문에 대한 대답이 "예"인 경우, 나눗셈은 컴파일 중인 함수의 맥락에서 곱셈과 FMA로 변환할 수 있습니다.위에서 설명한 정적 분석은 질문 2., 3. 및 4에 답하는 데 도움이 됩니다.

(**) '부호를 만지작거린다'는 것은 x가 두 개의 부호 0 중 하나일 경우에만 결과를 정확하게 도출하기 위해 필요한 경우 FMA(C1, x, (-C2)*x) 대신 -FMA(-C1, x, (-C2)*x)를 사용하는 것을 의미한다.

@Pascal의 답변은 마음에 들지만 최적화에서는 완벽한 솔루션보다는 단순하고 잘 이해된 변환 서브셋을 사용하는 것이 좋습니다.

현재 및 일반적인 모든 과거 부동소수점 형식에는 이진수 가수라는 한 가지 공통점이 있었습니다.

따라서 모든 분수는 형태의 유리수였다.

x / 2ⁿ

이는 프로그램의 상수(및 모든 가능한 기본 10분수)와 대조되며, 이는 형식의 유리수이다.

x / (2ⁿ * 5^m )

따라서 하나의 최적화는 단순히 m == 0에 대해 입력과 역수를 테스트하는 것이다. 왜냐하면 이러한 숫자는 정확히 FP 형식으로 표현되며 이 숫자를 사용한 연산은 형식 내에서 정확한 숫자를 생성해야 하기 때문이다.

예를 들어, (10진수 2자리) 범위 내에서.01로.0.99다음과 같은 수로 나누거나 곱하는 것이 최적화됩니다.

.25 .50 .75

그리고 다른 건 다 그렇지 않아. (먼저 테스트해봐, 하하)

부동소수점 나눗셈의 결과는 다음과 같습니다.

• 기호의 깃발
• 중요 인물
• 지수
• 깃발 세트(조기, 언더플로, 부정확 등) - 참조fenv())

처음 3개를 맞히는 것만으로는 충분하지 않습니다(그러나 플래그 세트는 올바르지 않습니다).더 이상의 지식(예를 들어 결과의 어떤 부분이 실제로 중요한지, 배당 가능한 값 등)이 없다면 나눗셈을 상수(및/또는 난해한 FMA 혼란)로 곱셈으로 치환하는 것은 거의 안전하지 않다고 생각합니다.

또, 최신 CPU의 경우, 2개의 FMA로 사업부를 교체하는 것이 항상 개선이라고는 생각하지 않습니다.예를 들어 병목현상이 명령 가져오기/디코딩일 경우 이 "최적화"는 성능을 악화시킵니다.예를 들어 후속 명령이 결과에 의존하지 않는 경우(CPU가 결과를 기다리는 동안 다른 많은 명령을 병렬로 수행할 수 있음) FMA 버전에서는 여러 종속성 정지가 발생하여 성능이 저하될 수 있습니다.세 번째 예로, 모든 레지스터가 사용되고 있는 경우, FMA 버전(추가의 「라이브」변수 필요)은 「스필링」을 증가시켜 퍼포먼스를 악화시킬 가능성이 있습니다.

(다수는 아니지만 많은 경우) 2의 상수 배수로 나누거나 곱하는 것은 추가만으로 실행할 수 있습니다(특히 지수에 시프트 카운트를 추가).

언급URL : https://stackoverflow.com/questions/35527683/efficient-floating-point-division-with-constant-integer-divisors