i387 fsqrt 명령으로 올바른 반올림을 얻을 수 있는 방법이 있습니까?

StackOverflow https://stackoverflow.com//questions/9678224

문제

i387 fsqrt 명령으로 올바른 반올림을 얻을 수 있는 방법이 있습니까?...

...정밀 모드를 변경하는 것 외에도 x87 제어 단어에서 - 가능하다는 것을 알고 있지만 sqrt 작업이 중단되면 정밀도 모드가 잘못되는 불쾌한 재진입 유형 문제가 있기 때문에 합리적인 해결책이 아닙니다.

제가 다루고 있는 문제는 다음과 같습니다.x87 fsqrt opcode는 fpu 레지스터의 정밀도에서 올바르게 반올림된(IEEE 754에 따라) 제곱근 연산을 수행합니다. 이는 확장된(80비트) 정밀도라고 가정하겠습니다.그러나 이를 사용하여 결과를 올바르게 반올림하여(현재 반올림 모드에 따라) 효율적인 단일 및 이중 정밀도 제곱근 함수를 구현하고 싶습니다.결과의 정밀도가 너무 높기 때문에 결과를 단정밀도 또는 배정밀도 반올림으로 다시 변환하는 두 번째 단계에서는 올바르지 않게 반올림된 결과가 남을 수 있습니다.

일부 작업에서는 편향을 사용하여 이 문제를 해결할 수 있습니다.예를 들어, 배정밀도 값의 52개 유효 비트를 63비트 확장 정밀도 가수의 마지막 52비트에 강제로 적용하는 2의 거듭제곱 형태로 바이어스를 추가하여 덧셈 결과에서 과도한 정밀도를 방지할 수 있습니다. .그러나 제곱근을 사용하여 그러한 트릭을 수행하는 확실한 방법은 없습니다.

기발한 아이디어가 있나요?

(의도된 애플리케이션이 C의 구현이기 때문에 C 태그도 지정되었습니다. sqrt 그리고 sqrtf 기능.)

도움이 되었습니까?

해결책

먼저, 명백한 내용을 제거해 보겠습니다.x87 대신 SSE를 사용해야 합니다.SSE sqrtss 그리고 sqrtsd 명령어는 원하는 작업을 정확하게 수행하고 모든 최신 x86 시스템에서 지원되며 속도도 훨씬 빠릅니다.

이제 x87 사용을 고집하신다면 좋은 소식부터 시작하겠습니다.float에 대해서는 아무것도 할 필요가 없습니다.당신은 필요 2p + 2 비트를 사용하여 p비트 부동 소수점 형식으로 올바르게 반올림된 제곱근을 계산합니다.왜냐하면 80 > 2*24 + 2, 단정밀도에 대한 추가 반올림은 항상 올바르게 반올림되며 올바르게 반올림된 제곱근을 갖게 됩니다.

이제 나쁜 소식은 다음과 같습니다. 80 < 2*53 + 2, 따라서 이중 정밀도에는 그런 행운이 없습니다.몇 가지 해결 방법을 제안해 드릴 수 있습니다.여기 내 머리 꼭대기에 좋은 쉬운 것이 있습니다.

  1. 허락하다 y = round_to_double(x87_square_root(x));
  2. Dekker(헤드-테일) 제품을 사용하여 계산 a 그리고 b 그렇게 y*y = a + b 정확히.
  3. 잔차를 계산하다 r = x - a - b.
  4. if (r == 0) return y
  5. if (r > 0), 허락하다 y1 = y + 1 ulp, 그리고 계산 a1, b1 성. y1*y1 = a1 + b1.비교하다 r1 = x - a1 - b1 에게 r, 다음 중 하나를 반환합니다. y 또는 y1, 어느 것이 더 작은 잔차를 갖는지에 따라(또는 잔차의 크기가 동일한 경우 하위 비트가 0인 것).
  6. if (r < 0), 동일한 작업을 수행합니다. y1 = y - 1 ulp.

이 절차는 기본 반올림 모드만 처리합니다.그러나 방향성 반올림 모드에서는 단순히 대상 형식으로 반올림하는 것이 올바른 일입니다.

다른 팁

좋습니다. 더 나은 해결책이 있는 것 같습니다.

  1. 컴퓨팅 y=sqrt(x) 확장된 정밀도(fsqrt).
  2. 마지막 11비트가 아닌 경우 0x400, 단순히 배정밀도로 변환하고 반환하면 됩니다.
  3. 추가하다 0x100-(fpu_status_word&0x200) 확장된 정밀도 표현의 낮은 단어로.
  4. 배정밀도로 변환하고 반환합니다.

3단계는 다음과 같은 경우에만 상태어의 C1 비트(0x200)가 1이라는 사실을 기반으로 합니다. fsqrt의 결과는 반올림되었습니다.이는 2단계의 테스트로 인해 유효합니다. x 완전한 정사각형은 아니었습니다.만약 완전제곱수라면, y 배정밀도 이상의 비트는 없습니다.

비트 표현 작업을 하고 다시 로드하는 것보다 조건부 부동 소수점 연산을 사용하여 3단계를 수행하는 것이 더 빠를 수도 있습니다.

코드는 다음과 같습니다(모든 경우에 작동하는 것 같습니다).

sqrt:
    fldl 4(%esp)
    fsqrt
    fstsw %ax
    sub $12,%esp
    fld %st(0)
    fstpt (%esp)
    mov (%esp),%ecx
    and $0x7ff,%ecx
    cmp $0x400,%ecx
    jnz 1f
    and $0x200,%eax
    sub $0x100,%eax
    sub %eax,(%esp)
    fstp %st(0)
    fldt (%esp)
1:  add $12,%esp
    fstpl 4(%esp)
    fldl 4(%esp)
    ret

387의 장점을 활용하지 않기 때문에 원하는 것이 아닐 수도 있습니다. fsqrt 지시사항이 있지만 놀랍도록 효율적입니다. sqrtf(float) ~에 glibc 32비트 정수 연산으로 구현되었습니다.NaN, Infs, 비정규도 올바르게 처리합니다. 실제 x87 명령어/FP 제어 단어 플래그를 사용하여 이러한 검사 중 일부를 제거하는 것이 가능할 수도 있습니다.보다: glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

그만큼 dbl-64/e_sqrt.c 코드가 그다지 친숙하지 않습니다.한 눈에 어떤 가정이 만들어지고 있는지 알기는 어렵습니다.신기하게도 도서관의 i386은 sqrt[f|l] 구현은 그냥 호출 fsqrt, 이지만 값을 다르게 로드합니다. flds SP의 경우, fldl DP를 위해.

라이센스 : CC-BY-SA ~와 함께 속성
제휴하지 않습니다 StackOverflow
scroll top