Има ли някакъв начин да получите правилно закръгляване с инструкцията 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. използвайте продукт на 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 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-битова целочислена аритметика. Той дори обработва NaNs, Infs, 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.