C'è un modo per avere il corretto arrotondamento con il i387 fsqrt istruzione?

Question

C'è un modo per avere il corretto arrotondamento con il i387 fsqrt istruzione?...

...oltre a cambiare la modalità di precisione nel 87 x parola di controllo - so che è possibile, ma non è una soluzione ragionevole, perché è brutto rientranza-tipo di problemi in cui la modalità di precisione sarà sbagliato se il sqrt operazione viene interrotta.

Il problema che sto per trattare è il seguente:il 87 x fsqrt opcode esegue correttamente arrotondato (per lo standard IEEE 754) radice quadrata di funzionamento nella precisione delle fpu registri, che si assume è estesa (80 bit) di precisione.Tuttavia, voglio utilizzare per attuare efficienti precisione singola e doppia radice quadrata con funzioni correttamente i risultati arrotondati (secondo le attuali modalità di arrotondamento).Dal momento che il risultato è un eccesso di precisione, la seconda fase di conversione il risultato a precisione singola o doppia, tornate di nuovo, lasciando eventualmente un non-correttamente-risultato arrotondato.

Con alcune operazioni è possibile risolvere il problema con i pregiudizi.Per esempio, posso evitare un eccesso di precisione nei risultati (oltre che con l'aggiunta di un bias nella forma di una potenza di due, che costringe il 52 bit più significativi di un doppio valore di precisione nelle ultime 52 bit 63-bit extended precisione mantissa.Ma io non vedo alcun modo ovvio per fare un trucco con la radice quadrata.

Eventuali idee intelligenti?

(Anche tagged C, in quanto la sua applicazione è l'attuazione del C sqrt e sqrtf funzioni.)

Answer 1

Iniziamo prima di tutto con l'evidente fuori strada:si dovrebbe essere utilizzando SSE invece di 87 x.SSE sqrtss e sqrtsd istruzioni per fare esattamente quello che vuoi, sono supportati su tutti i moderni sistemi x86, e sono significativamente più veloce.

Ora, se ti ostini a usare 87 x, io inizio con la buona notizia:non c'è bisogno di fare nulla per il flottante.Hai bisogno 2p + 2 bit per calcolare correttamente un quadrato con angoli arrotondati-radice p-bit floating-point formato.Perché 80 > 2*24 + 2, l'ulteriore arrotondamento a singola precisione, sempre correttamente, e si dispone di un arrotondate, radice quadrata.

Ora la cattiva notizia: 80 < 2*53 + 2, in modo che nessuna tale fortuna per la doppia precisione.Posso suggerire diverse soluzioni alternative;ecco una bella e facile largo della parte superiore della mia testa.

1. lasciate y = round_to_double(x87_square_root(x));
2. utilizzare un Dekker (testa-coda) prodotto da calcolare a e b tale che y*y = a + b esattamente.
3. calcolare il residuo r = x - a - b.
4. if (r == 0) return y
5. if (r > 0), lasciare y1 = y + 1 ulp, e calcolare a1, b1 s.t. y1*y1 = a1 + b1.Confrontare r1 = x - a1 - b1 per r, e restituire y o y1, a seconda di quale è il più piccolo residuo (o quella con zero bit di ordine inferiore, se i residui sono uguali in modulo).
6. if (r < 0), fare la stessa cosa per y1 = y - 1 ulp.

Questa procedura gestisce solo la modalità predefinita;tuttavia, in diretto modalità di arrotondamento, semplicemente arrotondamento per il formato di destinazione è la cosa giusta.

Answer 2

OK, penso di avere una soluzione migliore:

.
1. Compute y=sqrt(x) in precisione estesa (fsqrt).
2. Se gli ultimi 11 bit non sono 0x400, semplicemente converti con doppia precisione e ritorno.
3. Aggiungi 0x100-(fpu_status_word&0x200) alla parola bassa della rappresentazione di precisione estesa.
4. Converti in doppia precisione e ritorno.

Step 3 si basa sul fatto che il bit C1 (0x200) della parola di stato è 1 se e solo se il risultato di fsqrt è stato arrotondato.Questo è valido perché, a causa del test al punto 2, x non era un quadrato perfetto;Se fosse un quadrato perfetto, y non avrebbe bit oltre a doppia precisione.

Potrebbe essere più veloce eseguire il passaggio 3 con un punto flottante condizionale che funziona piuttosto che lavorare sulla rappresentazione di bit e la ricarica.

Ecco il codice (sembra funzionare in tutti i casi):

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

Potrebbe non essere quello che vuoi, poiché non approfitta delle 387 istruzioni fsqrt, ma c'è un sqrtf(float) di generazione sorprendentemente efficiente in ibc implementato con un numero intero a 32 bit Aritmetico.Gestisce anche Nans, Infs, Subnormals correttamente - Potrebbe essere possibile eliminare alcuni di questi controlli con le istruzioni Real X87 / Bandiere della parola di controllo FP.Vedi: glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

Il codice dbl-64/e_sqrt.c non è così amichevole.È difficile dire quali ipotesi vengono fatte a colpo d'occhio.Curiosamente, le implementazioni sqrt[f|l] I386 della biblioteca hanno appena chiamato fsqrt, ma carica il valore in modo diverso.flds per SP, fldl per DP.