Existe alguma maneira de obter o arredondamento correto com a instrução i387 fsqrt?

Question

Existe alguma maneira de obter o arredondamento correto com a instrução i387 fsqrt?...

...além de alterar o modo de precisão na palavra de controle x87 - eu sei que isso é possível, mas não é uma solução razoável porque tem problemas desagradáveis ​​​​do tipo reentrada, onde o modo de precisão estará errado se a operação sqrt for interrompida.

O problema com o qual estou lidando é o seguinte:o x87 fsqrt opcode executa uma operação de raiz quadrada arredondada corretamente (de acordo com IEEE 754) na precisão dos registros fpu, que presumo ser uma precisão estendida (80 bits).No entanto, quero usá-lo para implementar funções eficientes de raiz quadrada de precisão simples e dupla com os resultados arredondados corretamente (de acordo com o modo de arredondamento atual).Como o resultado tem excesso de precisão, a segunda etapa é converter o resultado em arredondamentos de precisão simples ou dupla novamente, possivelmente deixando um resultado não arredondado corretamente.

Com algumas operações é possível contornar isso com vieses.Por exemplo, posso evitar o excesso de precisão nos resultados da adição adicionando uma tendência na forma de uma potência de dois que força os 52 bits significativos de um valor de precisão dupla para os últimos 52 bits da mantissa de precisão estendida de 63 bits. .Mas não vejo nenhuma maneira óbvia de fazer esse truque com raiz quadrada.

Alguma ideia inteligente?

(Também marcado como C porque o aplicativo pretendido é a implementação do C sqrt e sqrtf funções.)

Answer 1

Primeiro, vamos tirar o óbvio do caminho:você deveria usar SSE em vez de x87.A ESS sqrtss e sqrtsd as instruções fazem exatamente o que você deseja, são suportadas em todos os sistemas x86 modernos e também são significativamente mais rápidas.

Agora, se você insiste em usar x87, começarei com uma boa notícia:você não precisa fazer nada para flutuar.Você precisa 2p + 2 bits para calcular uma raiz quadrada corretamente arredondada em um formato de ponto flutuante de p bits.Porque 80 > 2*24 + 2, o arredondamento adicional para precisão simples sempre arredondará corretamente e você terá uma raiz quadrada arredondada corretamente.

Agora a má notícia: 80 < 2*53 + 2, então não tive tanta sorte com precisão dupla.Posso sugerir várias soluções alternativas;aqui está uma boa e fácil que veio da minha cabeça.

1. deixar y = round_to_double(x87_square_root(x));
2. use um produto Dekker (head-tail) para calcular a e b de tal modo que y*y = a + b exatamente.
3. calcular o resíduo r = x - a - b.
4. if (r == 0) return y
5. if (r > 0), deixar y1 = y + 1 ulp, e calcular a1, b1 st. y1*y1 = a1 + b1.Comparar r1 = x - a1 - b1 para r, e retorne também y ou y1, dependendo de qual tem o resíduo menor (ou aquele com zero bit de ordem inferior, se os resíduos forem iguais em magnitude).
6. if (r < 0), faça a mesma coisa para y1 = y - 1 ulp.

Este procedimento trata apenas do modo de arredondamento padrão;entretanto, nos modos de arredondamento direcionado, simplesmente arredondar para o formato de destino faz a coisa certa.

Answer 2

OK, acho que tenho uma solução melhor:

1. Calcular y=sqrt(x) em precisão estendida (fsqrt).
2. Se os últimos 11 bits não forem 0x400, simplesmente converta para precisão dupla e retorne.
3. Adicionar 0x100-(fpu_status_word&0x200) à palavra baixa da representação de precisão estendida.
4. Converta para precisão dupla e retorne.

O passo 3 é baseado no fato de que o bit C1 (0x200) da palavra de status é 1 se e somente se fsqrto resultado foi arredondado.Isto é válido porque, devido ao teste na etapa 2, x não era um quadrado perfeito;se fosse um quadrado perfeito, y não teria bits além da precisão dupla.

Pode ser mais rápido executar a etapa 3 com um ponto flutuante condicional operando em vez de trabalhar na representação e recarregamento do bit.

Aqui está o código (parece funcionar em todos os casos):

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

Pode não ser o que você deseja, pois não aproveita o 387 fsqrt instrução, mas há uma surpreendentemente eficiente sqrtf(float) em glibc implementado com aritmética inteira de 32 bits.Ele até lida com NaNs, Infs e subnormais corretamente - pode ser possível eliminar algumas dessas verificações com instruções x87 reais/sinalizadores de palavras de controle FP.ver: glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

O dbl-64/e_sqrt.c o código não é tão amigável.É difícil dizer à primeira vista quais suposições estão sendo feitas.Curiosamente, o i386 da biblioteca sqrt[f|l] implementações basta ligar fsqrt, mas carregue o valor de forma diferente. flds para SP, fldl para DP.