Есть ли способ получить правильное округление с помощью инструкции i387 fsqrt?

Question

Есть ли способ получить правильное округление с помощью инструкции i387 fsqrt?...

...кроме изменения режима точности в управляющем слове x87 - я знаю, что это возможно, но это неразумное решение, поскольку оно имеет неприятные проблемы типа повторного входа, когда режим точности будет неправильным, если операция sqrt будет прервана.

Проблема, с которой я имею дело, заключается в следующем:х87 fsqrt код операции выполняет правильно округленную (согласно IEEE 754) операцию квадратного корня с точностью регистров fpu, что, как я предполагаю, является расширенной (80-битной) точностью.Однако я хочу использовать его для реализации эффективных функций квадратного корня одинарной и двойной точности с правильным округлением результатов (в соответствии с текущим режимом округления).Поскольку результат имеет избыточную точность, на втором этапе снова преобразуется результат в округления с одинарной или двойной точностью, что может привести к получению неправильно округленного результата.

С помощью некоторых операций это можно обойти с помощью смещений.Например, я могу избежать избыточной точности результатов сложения, добавив смещение в виде степени двойки, которая принудительно помещает 52 значащих бита значения двойной точности в последние 52 бита 63-битной мантиссы повышенной точности. .Но я не вижу очевидного способа проделать такой трюк с квадратным корнем.

Есть умные идеи?

(Также помечен как C, потому что предполагаемое приложение — это реализация C. sqrt и sqrtf функции.)

Answer 1

Для начала давайте разберемся с очевидным:вам следует использовать SSE вместо x87.ССЕ 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. используйте произведение Деккера (голова-хвост) для вычисления 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, в зависимости от того, какой из них имеет меньший остаток (или тот, у которого младший бит равен нулю, если остатки равны по величине).
6. if (r < 0), сделайте то же самое для y1 = y - 1 ulp.

Эта процедура обрабатывает только режим округления по умолчанию;однако в режимах направленного округления достаточно простого округления до целевого формата.

Answer 2

Хорошо, я думаю, у меня есть лучшее решение:

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

Answer 3

Возможно, это не то, что вам нужно, поскольку оно не использует преимущества 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 для СП, fldl для ДП.