Cálculo robusto da média de dois números em ponto flutuante?

Vamos x, yser dois números de ponto flutuante. Qual é o caminho certo para calcular a média deles?

A maneira ingênua (x+y)/2pode resultar em estouros quando xe ysão muito grandes. Acho que 0.5 * x + 0.5 * ytalvez seja melhor, mas envolve duas multiplicações (o que talvez seja ineficiente), e não tenho certeza se é bom o suficiente. Existe uma maneira melhor?

Outra idéia com a qual estou brincando é (y/2)(1 + x/y)se x<=y. Mais uma vez, não tenho certeza de como analisar isso e provar que ele atende aos meus requisitos.

Além disso, preciso de uma garantia de que a média calculada será >= min(x,y)e <= max(x,y). Como apontado na resposta de Don Hatch , talvez uma maneira melhor de fazer essa pergunta seja: o que é uma implementação da média de dois números que sempre fornece o resultado mais preciso possível? Ou seja, se xe ysão números de ponto flutuante, como calcular o número de ponto flutuante mais próximo (x+y)/2? Nesse caso, a média calculada é automaticamente >= min(x,y)e <= max(x,y). Veja a resposta de Don Hatch para detalhes.

Nota: Minha prioridade é uma precisão robusta. Eficiência é dispensável. No entanto, se houver muitos algoritmos robustos e precisos, eu escolheria o mais eficiente.

Penso que a precisão e a estabilidade dos algoritmos numéricos de Higham abordam como se pode analisar esses tipos de problemas. Veja o Capítulo 2, especialmente o exercício 2.8.

Nesta resposta, gostaria de salientar algo que não é realmente abordado no livro de Higham (não parece ser muito conhecido, por falar nisso). Se você estiver interessado em provar propriedades de algoritmos numéricos simples como esses, poderá usar o poder dos modernos solucionadores de SMT ( Teorias do Módulo de Satisfação ), como o z3 , usando um pacote como o sbv em Haskell. Isso é um pouco mais fácil do que usar lápis e papel.

Suponha que me seja dado e gostaria de saber se z = ( x + y ) / 2 satisfaz x ≤ z ≤ y . O seguinte código Haskell$0\leq x\leq y$$z=(x+y)/2$$x\leq z\leq y$

import Data.SBV

test1 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test1 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ 0 .<= x &&& x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

test2 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test2 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

vai me deixar fazer isso automaticamente . Aqui test1 funestá a proposição de que para todos os flutuadores finitos x , y com 0 ≤ x ≤ y .$x \leq \mathit{fun}(x,y) \leq y$$x,y$$0\leq x\leq y$

λ> prove $ test1 (\x y -> (x + y) / 2)
Falsifiable. Counter-example:
  x = 2.3089316e36 :: Float
  y = 3.379786e38 :: Float

Transborda. Suponha que agora eu use sua outra fórmula: $z=x/2+y/2$

λ> prove $ test1 (\x y -> x/2 + y/2)
Falsifiable. Counter-example:
  x = 2.3509886e-38 :: Float
  y = 2.3509886e-38 :: Float

Não funciona (devido a underflow gradual: , que pode ser unintuitive devido a toda base-2 estar aritmética).$(x/2)\times2 \neq x$

Agora tente :$z=x + (y-x)/2$

λ> prove $ test1 (\x y -> x + (y-x)/2)
Q.E.D.

Trabalho! A Q.E.D.é uma prova de que a test1propriedade possui para todos os carros alegóricos, conforme definido acima.

E o mesmo, mas restrito a (em vez de 0 ≤ x ≤ y )?$x\leq y$$0\leq x\leq y$

λ> prove $ test2 (\x y -> x + (y-x)/2)
Falsifiable. Counter-example:
  x = -3.1300826e34 :: Float
  y = 3.402721e38 :: Float

Ok, então, se estourar, que tal z = x + ( y / 2 - x / 2 ) ?$y-x$$z = x + (y/2-x/2)$

λ> prove $ test2 (\x y -> x + (y/2 - x/2))
Q.E.D.

$x + (y/2 - x/2)$

$(x+y)/2$

$x \leq x + (y/2-x/2) \leq y$SFloatSDouble

-ffast-math$(x+y)/2$

PPPS Eu me empolguei um pouco olhando apenas para expressões algébricas simples, sem condicionais. A fórmula de Don Hatch é estritamente melhor.

Primeiro, observe que, se você tiver um método que dê uma resposta mais precisa em todos os casos, ele atenderá às condições necessárias. (Note que eu digo uma resposta mais precisa, em vez da resposta mais precisa, pois pode haver dois vencedores.) Prova: Se, ao contrário, você tem uma resposta precisa-como-possível que se não satisfazer a condição exigida, que significa tanto answer<min(x,y)<=max(x,y)(nesse caso, min(x,y)uma resposta melhor, uma contradição) ou min(x,y)<=max(x,y)<answer(nesse caso, max(x,y)uma resposta melhor, uma contradição).

Então, eu acho que isso significa que sua pergunta se resume a encontrar uma resposta mais precisa possível. Supondo aritmética IEEE754, proponho o seguinte:

if max(abs(x),abs(y)) >= 1.:
    return x/2. + y/2.
else:
    return (x+y)/2.

Meu argumento de que isso fornece uma resposta mais precisa é uma análise de caso um tanto tediosa. Aqui vai:

• Caso max(abs(x),abs(y)) >= 1.:

  • Subcasca nem x nem y são desnormalizados: Nesse caso, a resposta computada x/2.+y/2.manipula as mesmas mantissas e, portanto, fornece exatamente a mesma resposta que a computação (x+y)/2produziria se assumíssemos expoentes estendidos para impedir o transbordamento. Essa resposta pode depender do modo de arredondamento, mas, em qualquer caso, é garantida pela IEEE754 como a melhor resposta possível (pelo fato de o computador x+yser garantido como a melhor aproximação para x + y matemático, e a divisão por 2 é exata nesta caso).

  • A subcasca x é desnormalizada (e assim abs(y)>=1):

    answer = x/2. + y/2. = y/2. since abs(x/2.) is so tiny compared to abs(y/2.) = the exact mathematical value of y/2 = a best possible answer.

  • A sub-caixa y é desnormalizada (e assim abs(x)>=1): análoga.

• Caso max(abs(x),abs(y)) < 1.:
  • Subcategoria que o computado x+yé não desnormalizado ou desnormalizado e "par": embora o calculado x+ypossa não ser exato, o IEEE754 garante que a IEEE754 é a melhor aproximação possível da matemática x + y. Nesse caso, a divisão subsequente por 2 na expressão (x+y)/2.é exata, portanto a resposta calculada (x+y)/2.é a melhor aproximação possível à matemática (x + y) / 2.
  • Subcategoria o computado x+yé desnormalizado e "ímpar": Nesse caso, exatamente um de x, y também deve ser desnormalizado e "ímpar", o que significa que o outro de x, y é desnormalizado com o sinal oposto e, portanto, o computado x+yé exatamente o matemático x + y, e, portanto, o (x+y)/2.IEEE754 calculado é garantido como a melhor aproximação possível do matemático (x + y) / 2.

Para os formatos binários de ponto flutuante IEEE-754, exemplificados por binary64computação (precisão dupla), S. Boldo provou formalmente que o algoritmo simples mostrado abaixo fornece a média arredondada corretamente.

Sylvie Boldo, "Verificação formal de programas que computam a média de ponto flutuante". Na Conferência Internacional sobre Métodos Formais de Engenharia , pp. 17-32. Springer, Cham, 2015. ( rascunho online )

$(x+y)/2$$x/2 + y/2$binary64$C \in [2^{-967}, 2^{970}]$$C$ para fornecer o melhor desempenho para um caso de uso específico.

Isso produz o seguinte ISO-C99código exemplar :

double average (double x, double y) 
{
    const double C = 1; /* 0x1p-967 <= C <= 0x1p970 */
    return (C <= fabs (x)) ? (x / 2 + y / 2) : ((x + y) / 2);
}

Em trabalho de acompanhamento recente, S. Boldo e co-autores mostraram como obter os melhores resultados possíveis para os formatos decimais de ponto flutuante IEEE-754, usando as operações FMA (Multiply Add) fundidas e uma conhecida ferramenta de precisão. bloco de construção duplicado (TwoSum):

Sylvie Boldo, Florian Faissole e Vincent Tourneur, "um algoritmo formalmente provado para calcular a média correta dos números decimais de ponto flutuante". No 25º Simpósio IEEE sobre Aritmética Computacional (ARITH 25) , junho de 2018, pp. 69-75. ( rascunho online )

Embora possa não ser supereficiente em termos de desempenho, existe uma maneira muito simples de (1) garantir que nenhum dos números seja maior que xou y(sem estouros) e (2) manter o ponto flutuante tão "preciso" quanto possível. possível (e (3) , como um bônus adicional, mesmo que a subtração esteja sendo usada, nenhum valor será armazenado como números negativos.

float difference = max(x, y) - min(x, y);
return min(x, y) + (difference / 2.0);

De fato, se você realmente deseja obter precisão, nem precisa executar a divisão no local; basta retornar os valores min(x, y)e os differencequais você pode usar para simplificar logicamente ou manipular posteriormente.

Converta para uma precissão mais alta, adicione os valores lá e converta novamente.

Não deve haver excesso na precissão mais alta e, se ambos estiverem na faixa de ponto flutuante válida, o número calculado também deverá estar dentro.

E deve estar entre eles, na pior das hipóteses, apenas metade do número maior, se a precissão não for suficiente.

Teoricamente, x/2pode ser calculado subtraindo 1 da mantissa.

No entanto, a implementação de operações bit a bit como essa não é necessariamente direta, principalmente se você não souber o formato dos seus números de ponto flutuante.

Se você puder fazer isso, toda a operação será reduzida para 3 adições / subtrações, o que deve ser uma melhoria significativa.

Eu estava pensando na mesma linha que @Roland Heath, mas não posso comentar ainda, aqui está a minha opinião:

x/2pode ser calculado subtraindo 1 do expoente (não a mantissa, subtrair 1 da mantissa está subtraindo 2^(value_of_exponent-length_of_mantissa)do valor geral).

Sem restrição do caso geral, vamos assumir x < y. (Se x > y, re-rotule as variáveis. Se x = y, (x+y) / 2é trivial.)

• Transforme (x+y) / 2em x/2 + y/2, que pode ser executado por duas subtrações de número inteiro (por uma do expoente)
  • No entanto, há um limite inferior no expoente, dependendo da sua representação. Se o seu expoente já é mínimo antes de subtrair 1, esse método exigirá tratamento especial de caso. Um expoente mínimo ativado xtornará x/2menor que representável (assumindo que a mantissa seja representada com um líder implícito 1).
  • Em vez de subtrair 1 do expoente de x, mova xa mantissa da direita para uma (e adicione o líder implícito 1, se houver).
  • Subtraia 1 do expoente de y, se não for mínimo. Se for mínimo (y é maior que x, por causa da mantissa), desloque a mantissa para a direita por uma (adicione 1 implícito à esquerda, se houver).
  • Desloque a nova mantissa de xpara a direita de acordo com o expoente de y.
  • Realize adição de número inteiro na mantissa, a menos que a mantissa xtenha sido completamente deslocada. Se os dois expoentes forem mínimos, os principais transbordarão, o que é aceitável, porque esse transbordamento se tornará um líder implícito novamente.
• e uma adição de ponto flutuante.
  • Não consigo pensar em nenhum caso especial aqui; exceto o arredondamento, que também se aplica às mudanças descritas acima.