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

Во-первых, давайте уберем очевидное: вы должны использовать SSE вместо x87. Инструкции 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. используйте произведение Деккера (голова-хвост) для вычисления 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 s.t. y1*y1 = a1 + b1. Сравните r1 = x - a1 - b1 с r и верните либо y, либо y1, в зависимости от того, какой остаток имеет меньший остаток (или тот, у которого нулевой бит младшего разряда, если остатки равны по величине).
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, Inf и субнормальные значения - возможно, можно будет устранить некоторые из этих проверок с помощью реальных инструкций 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.