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

Question

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

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

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

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

Какие-нибудь умные идеи?

(также помечен как C, поскольку предполагаемое приложение - реализация функций C sqrt и sqrtf.)

Answer 1

Во-первых, давайте разберемся с очевидным: вы должны использовать SSE вместо x87.Инструкции SSE sqrtss и sqrtsd выполняют именно то, что вы хотите, поддерживаются во всех современных системах x86 и также значительно быстрее.

Теперь, если вы настаиваете на использовании x87, я начнус хорошими новостями: вам не нужно ничего делать для плавания.Вам нужно 2p + 2 бит для вычисления правильно округленного квадратного корня в формате с плавающей запятой p-бит.Потому что 80 > 2*24 + 2, дополнительное округление до одинарной точности всегда будет правильно округляться, и у вас будет правильно округленный квадратный корень.

Теперь плохие новости: 80 < 2*53 + 2, так что не повезло с двойной точностью.Я могу предложить несколько обходных путей;Вот хороший простой пример:

1. позвольте y = round_to_double(x87_square_root(x));
2. использовать продукт Dekker (head-tail) для вычисления 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, Inf, subnormals - возможно, удастся устранить некоторые из этих проверок с помощью реальных команд 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.