Fiquei surpreso ao descobrir que é difícil encontrar informações sobre esse problema em russo, como se poucas pessoas se importassem com o fato de que as bibliotecas matemáticas usadas em compiladores modernos às vezes não fornecem um resultado arredondado corretamente. Estou preocupado com esta situação, já que estou apenas trabalhando no desenvolvimento dessas bibliotecas matemáticas. Na literatura estrangeira, este problema está bem coberto, então decidi apresentá-lo em russo de uma forma científica popular, contando com fontes ocidentais e ainda um pouco de experiência pessoal.

---

Amigos, para sua comodidade, o artigo também está disponível em formato de apresentação de vídeo (quase 34 minutos), este formato é mais adequado para aqueles leitores que têm dificuldade em construir na cabeça as imagens matemáticas necessárias, visto que há muito material ilustrativo na apresentação. As informações do vídeo são totalmente idênticas ao conteúdo do artigo. Aja conforme sua conveniência.

---

Repito que este não é um artigo científico, mas sim um artigo de ciência popular, após a leitura do qual você se familiarizará brevemente com ele.

• Funções elementares transcendentais (exp, sin, log, cosh e outras) que funcionam com aritmética de ponto flutuante são arredondadas incorretamente, às vezes cometem um erro no último bit.
• A razão dos erros nem sempre reside na preguiça ou na baixa qualificação dos desenvolvedores, mas em uma circunstância fundamental, que a ciência moderna ainda não foi capaz de superar.
• «», - .
• , , , , exp2(x) pow(2.0, x).

Para entender o conteúdo deste artigo, você precisa estar familiarizado com o formato de ponto flutuante IEEE-754. Será suficiente se você pelo menos entender que, por exemplo, isto é: 0x400921FB54442D18 - número pi em formato de precisão dupla (binário64, ou duplo), ou seja, você acabou de entender o que quero dizer com esse registro; Não preciso ser capaz de fazer essas transformações instantaneamente. E vou lembrá-lo sobre os modos de arredondamento neste artigo, esta é uma parte importante da história. Também é desejável saber inglês para "programadores", porque haverá termos e citações da literatura ocidental, mas você pode se dar bem com um tradutor online.



Exemplos primeiro, para que você entenda imediatamente qual é o assunto da conversa. Agora darei o código em C ++, mas se esta não for a sua linguagem, tenho certeza de que ainda entenderá facilmente o que está escrito. Por favor, olhe este código:

#include <stdio.h>
#include <cmath>

int main() {
  float x = 0.00296957581304013729095458984375f;  // ,  .
  float z;
  z = exp2f(x);  // z = 2**x  .
  printf ("%.8f\n", z);  //      8   .
  z = powf(2.0f, x);  // z = 2**x  
  printf ("%.8f\n", z);  //   .
  return 0;
}

O número x é intencionalmente escrito com tal número de dígitos significativos para que seja exatamente representável no tipo float, ou seja, para que o compilador o converta em um código binário sem arredondamento. Afinal, você sabe muito bem que alguns compiladores não conseguem arredondar sem erros (se não souber, indique nos comentários, vou escrever um artigo separado com exemplos). Em seguida, no programa, precisamos calcular 2 ^x , mas vamos fazer de duas maneiras: a função exp2f (x) e a exponenciação explícita de dois, powf (2.0f, x). O resultado, claro, será diferente, porque eu disse acima que as funções elementares não podem funcionar corretamente em todos os casos, e selecionei especialmente um exemplo para mostrar isso. Aqui está o resultado:

1.00206053
1.00206041

Quatro compiladores me deram estes valores: Microsoft C ++ (19.00.23026), Intel C ++ 15.0, GCC (6.3.0) e Clang (3.7.0). Eles diferem em um bit menos significativo. Aqui está o código hexadecimal para esses números:

0x3F804385  // 
0x3F804384  // 

Por favor, lembre-se deste exemplo, é nele que veremos a essência do problema um pouco mais tarde, mas por agora, para que você tenha uma impressão mais clara, consulte exemplos para o tipo de dados de precisão dupla (duplo, binário64) com algumas outras funções elementares. Eu apresento os resultados na tabela. As respostas corretas (quando disponíveis) têm * no final.







Espero que não fique com a impressão de que fiz alguns testes completamente únicos que você dificilmente consegue encontrar. Nesse caso, vamos preparar de joelhos uma enumeração completa de todos os argumentos fracionários possíveis para a função 2 ^x para o tipo de dados float. É claro que estamos interessados ​​apenas em valores de x entre 0 e 1, pois outros argumentos produzirão um resultado que difere apenas no valor do campo expoente e não são de interesse. Você mesmo entende:

$$

$$2^x = 2^{[x]}\cdot2^{\{x\}}.$$

Tendo escrito tal programa (o texto oculto estará abaixo), verifiquei a função exp2f e quantos valores errôneos ela produz no intervalo x de 0 a 1.







É claro a partir do programa abaixo que o número de argumentos x testados foi 197612997. Acontece que, por exemplo, o Microsoft C ++ calcula incorretamente a função 2 ^x para quase um por cento deles. Não fiquem felizes, queridos fãs do GCC e do Clang, é que esta função está implementada corretamente nestes compiladores, mas cheia de erros em outros.





Não vou aborrecer o leitor com esses exemplos, o principal aqui foi mostrar que as implementações modernas de funções transcendentais podem arredondar o último bit incorretamente e diferentes compiladores cometem erros em diferentes lugares, mas nenhum deles funcionará corretamente. A propósito, o padrão IEEE-754 permite esse erro no último bit (sobre o qual falarei a seguir), mas ainda me parece estranho: ok double, esse é um tipo de dado grande, mas float pode ser verificado pela força bruta! Foi tão difícil de fazer? Não é nada difícil e já mostrei um exemplo.



Nosso código de enumeração contém uma função "autoescrita" de cálculo correto 2 ^xusando um polinômio aproximado de 10º grau, e foi escrito em poucos minutos, porque tais polinômios são derivados automaticamente, por exemplo, no sistema de álgebra computacional Maple. Basta definir uma condição para que o polinômio forneça 54 bits de precisão (para esta função, 2 ^x ). Por que 54? Mas você logo descobrirá, logo depois que eu lhe contar a essência do problema e dizer por que, em princípio, agora é impossível criar funções transcendentais rápidas e corretas para o tipo de dados de precisão quádrupla (binário128), embora já haja tentativas de atacar esse problema em teoria.

Arredondamento padrão e o problema com ele

Se você não estiver imerso no desenvolvimento de bibliotecas matemáticas, não há nada de errado em esquecer a regra de arredondamento padrão para números de ponto flutuante de acordo com o padrão IEEE-754. Portanto, vou lembrá-lo disso. Se você se lembra de tudo bem, dê uma olhada pelo menos no final desta seção de qualquer maneira, você terá uma surpresa: eu vou lhe mostrar uma situação em que arredondar um número pode ser muito difícil.



Você pode facilmente lembrar o que é “arredondar para cima” (para mais infinito), “arredondar para baixo” (para menos infinito) ou “arredondar para zero” pelo nome (se houver alguma coisa, existe a Wikipedia) As principais dificuldades para os programadores surgem com o arredondamento "para o mais próximo, mas no caso de distância igual do mais próximo - para aquele com o último dígito par". Sim, é assim que se traduz este modo de arredondamento, que a literatura ocidental denomina resumidamente: "Arredondar os laços mais próximos ao par".



Este modo de arredondamento é usado por padrão e funciona da seguinte maneira. Se, como resultado dos cálculos, o comprimento da mantissa for maior do que o tipo de dados resultante pode acomodar, o arredondamento é executado para o mais próximo de dois valores possíveis. No entanto, pode surgir uma situação em que o número original acabou ficando exatamente no meio entre os dois mais próximos, então é escolhido o resultado para o qual o último bit (após o arredondamento) acaba sendo par, ou seja, igual a zero. Considere quatro exemplos em que você precisa arredondar para dois bits após o ponto decimal binário:

1. Arredondar 1,00 1 001. O terceiro bit depois da vírgula é 1, mas depois existe outro 6º bit, que é 1, o que significa que o arredondamento será para cima, porque o número original está mais próximo de 1,01 do que de 1,00.
2. 1,001000. , 1,00 1,01, .
3. 1,011000. 1,01 1,10. , .
4. 1,010111. , 1,01, 1,10.

A partir desses exemplos, pode parecer que tudo é simples, mas não é. O fato é que às vezes não podemos dizer com certeza se estamos realmente no meio entre dois valores. Veja um exemplo. Vamos novamente arredondar para dois bits após o ponto decimal:



1,00 1 0000000000000000000000000000000000001



Agora é óbvio para você que o arredondamento deve ser para cima, ou seja, para o número 1,01. No entanto, você está olhando para um número com 40 bits após o ponto decimal. E se o seu algoritmo não pudesse fornecer 40 bits de precisão e atingir apenas 30 bits? Em seguida, ele dará outro número:



1,00 1 000000000000000000000000000



Sem saber que na 40ª posição (que o algoritmo não é capaz de calcular) haverá uma posição querida, você arredonda esse número e obtém 1,00, o que está errado. Você arredondou a última parte incorretamente - esse é o assunto da nossa discussão. Do exposto, verifica-se que, para obter apenas o segundo bit correto, você deve calcular a função em até 40 bits! Uau! E se a "locomotiva" dos zeros for ainda mais longa? É sobre isso que falaremos na próxima seção.



A propósito, este é o erro que muitos compiladores cometem ao converter a notação decimal de um número de ponto flutuante para o formato binário resultante. Se o número decimal original no código do programa estiver muito próximo do meio entre dois valores binários representáveis ​​com precisão, ele não será arredondado corretamente. Mas este não é o tema deste artigo, mas um motivo para uma história separada.

A essência do problema de arredondar a última parte significativa

O problema se manifesta por dois motivos. O primeiro é a rejeição deliberada de cálculos demorados em favor da velocidade. Nesse caso, contanto que a precisão especificada seja observada, e quais bits haverá na resposta é uma questão secundária. A segunda razão é o dilema do criador de mesa, que é o assunto principal de nossa conversa. Vamos considerar os dois motivos com mais detalhes.

Primeira razão

Você, é claro, entende que o cálculo de funções transcendentais é implementado por alguns métodos aproximados, por exemplo, pelo método de aproximação de polinômios ou mesmo (raramente) por expansão em série. Para fazer os cálculos o mais rápido possível, os desenvolvedores concordam em realizar o mínimo possível de iterações do método numérico (ou tomar um polinômio do menor grau possível), desde que o algoritmo permita um erro que não exceda a metade do valor do último bit da mantissa. Na literatura, isso é escrito como 0,5ulp (ulp = unidade em último lugar ).



Por exemplo, se estamos falando de um número x do tipo float no intervalo (0,5; 1), o valor ulp = 2 ^-23 . No intervalo (1; 2) ulp = 2 ^-22 . Em outras palavras, se x está no intervalo (0; 1), então 2 ^xestará no intervalo (1,2), e para garantir uma precisão de 0,5ulp, você precisa, grosso modo, escolher EPS = 2 ^-23 (como iremos denotar a constante "épsilon", nas pessoas comuns chamadas "erro", ou "precisão", para quem como quiser, não encontre falhas, por favor).



Para cálculos aplicados, isso é suficiente, mas o fato de os últimos bits não coincidirem com o resultado absoluto não é importante para quase 100% dos programadores, pois para eles não importa quais serão os bits, mas qual será a precisão.



Para aqueles que não entendem, darei um exemplo no sistema numérico decimal. Aqui estão dois números: 1.999999 e 2.0. Digamos que o primeiro é o que o programador recebeu e o segundo é o padrão do que deveria ser obtido se tivéssemos possibilidades ilimitadas. A diferença entre eles é de apenas um milionésimo, ou seja, a resposta é calculada com um erro de EPS = 10 ^-6 . No entanto, não há um único número correto nesta resposta. É ruim? Não, do ponto de vista do programa aplicativo, isso é roxo, o programador arredondará a resposta, digamos, para duas casas decimais e obterá 2,00 (por exemplo, era sobre moeda, $ 2,00), ele não precisa de mais, mas o fato de que ele coloquei EPS = 10 ^-6 no meu programa , então bem feito, tirou uma margem para o erro de cálculos intermediários e resolveu o problema corretamente.



Em outras palavras, não se confunda: a precisão e o número de bits (ou dígitos) corretos são duas coisas diferentes. Aqueles que precisam de precisão (quase 100% dos programadores), o problema discutido não diz respeito a eles. Qualquer pessoa que precise de uma sequência de bits para corresponder a uma referência corretamente arredondada está muito preocupada com este problema, por exemplo, desenvolvedores de bibliotecas de funções elementares. No entanto, é útil para todos saberem sobre isso para desenvolvimento geral.



Deixe-me lembrá-lo de que essa foi a primeira direção do problema: as últimas partes da resposta podem estar erradas porque essa é uma solução intencional. O principal é manter a precisão de 0,5ulp (ou superior). Portanto, o algoritmo numérico é selecionado apenas a partir desta condição, se funcionar de forma extremamente rápida. Ao mesmo tempo, o padrão permite a implementação de funções elementares sem o arredondamento correto do último bit. Cito [1, seção 12.1] (inglês):

A versão de 1985 do Padrão IEEE 754 para Aritmética de Ponto Flutuante não especifica nada a respeito da função elementar. Isso porque se acreditou por anos que funções arredondadas corretamente seriam muito lentas, pelo menos para alguns argumentos de entrada. A situação mudou desde então e a versão 2008 da norma recomenda (mas não exige) que algumas funções sejam arredondadas corretamente.

As funções a seguir são recomendadas, mas não precisam ser arredondadas corretamente:

A segunda razão

Finalmente chegamos ao tópico da conversa: Table Maker's Dilemma (abreviado como TMD). Não consegui traduzir adequadamente seu nome para o russo, ele foi introduzido por William Kahan (fundador do IEEE-754) no artigo [2]. Talvez se você ler o artigo, você entenderá porque o nome é exatamente esse. Em suma, a essência do dilema é que precisamos obter um arredondamento absolutamente preciso da função z = f (x), como se tivéssemos um registro de bits infinitos do resultado z perfeitamente calculado à nossa disposição. Mas está claro para todos que não podemos obter uma sequência infinita. Quantos bits levar então? Acima, mostrei um exemplo em que precisamos ver 40 bits do resultado para obter pelo menos 2 bits corretos após o arredondamento. E a essência do problema TMD é que não sabemos com antecedência, até quantos bits calcular o valor de z para obter tantos bits corretos após o arredondamento quantos forem necessários. E se houver cem ou mil? Não sabemos com antecedência!



Por exemplo, como eu disse, para a função 2 ^x , para o tipo de dados float, onde a parte fracionária da mantissa tem apenas 23 bits, precisamos realizar o cálculo com uma precisão de 2 ^-54 para que o arredondamento ocorra corretamente para todos os argumentos x possíveis, sem exceção. Não é difícil obter esta estimativa por pesquisa exaustiva, mas para a maioria das outras funções, especialmente para os tipos double ou long double (coloque "classe" se você souber o que é), tais estimativas são desconhecidas .



Já vamos entender porque isso está acontecendo. Deliberadamente dei o primeiro exemplo neste artigo com o tipo de dado float e pedi que você se lembrasse dele, porque neste tipo existem apenas 32 bits e será mais fácil olhar para ele, em outros tipos de dados a situação é semelhante.



Começamos com o número x = 0,00296957581304013729095458984375, este é um número exatamente representável no tipo de dado float, ou seja, é escrito de forma que possa ser convertido para o sistema float binário sem arredondamento. Calculamos 2 ^x , e se tivéssemos uma calculadora com precisão infinita, então deveríamos obter (para que você possa verificar, o cálculo é feito no sistema online WolframAlpha ):



1.0020604729652405753669743044108123031635398201893943954577320057 ...



Vamos traduzir esse número para o binário, digamos que 64 bits sejam suficientes:



1.00000000100001110000100 1 000000000000000000000000000001101111101



O bit de arredondamento (24º bit após o ponto decimal) é sublinhado. Pergunta: para onde arredondar? Para cima ou para baixo? Claramente, você sabe disso porque vê pedaços suficientes e pode tomar uma decisão. Mas olhe com cuidado ...



Após o bit de arredondamento, temos 29 zeros. Isso significa que estamos muito, muito próximos do meio entre os dois números mais próximos e basta apenas nos movermos um pouco para baixo, pois a direção do arredondamento mudará. Mas a questão é: onde será essa mudança? O algoritmo numérico pode sequencialmente, passo a passo, aproximar o valor exato de diferentes lados, e até que passemos todos esses 29 zeros e alcancemos uma precisão que exceda o valor do último zero nesta "locomotiva", não saberemos a direção do arredondamento ... E se, de fato, a resposta correta fosse:



1.00000000100001110000100 0 1111111111111111111111111111111?



Então, o arredondamento será reduzido.



Não sabemos disso até que nossa precisão alcance 54º bit após a vírgula decimal. Quando o 54º bit for conhecido exatamente, saberemos exatamente de qual dos dois números mais próximos estamos realmente mais próximos. Esses números são chamados de pontos mais difíceis de arredondar [1, seção 12.3] (pontos críticos para arredondamento), e o número 54 é chamado de dureza para arredondar e é denotado pela letra m no livro citado.



A complexidade do arredondamento (m) é o número de bits que é o mínimo necessário para garantir que para todos os argumentos de alguma função f (x) e para um intervalo pré-selecionado, a função f (x) seja arredondada corretamente para o último bit (para diferentes modos de arredondamento, pode haver diferentes valor m). Em outras palavras, para o tipo de dados float e para o argumento x do intervalo (0; 1) para o modo de arredondamento "mais próximo par", a complexidade do arredondamento é m = 54. Isso significa que, para absolutamente todos os x do intervalo (0; 1), podemos colocar no algoritmo a mesma precisão ESP = 2 ^-54 , e todos os resultados serão arredondados corretamente para 23 bits após o ponto decimal binário.



Na verdade, alguns algoritmos são capazes de fornecer um resultado exato e com base em 53 e até 52 bits, a força bruta mostra isso, mas teoricamente, 54 é necessário. Se não fosse pela possibilidade de desencadear a força bruta, não seríamos capazes e economize alguns bits, como fiz no programa de força bruta acima. Peguei um polinômio com um grau menor do que deveria, mas ainda funciona, só porque tive sorte.



Assim, independentemente do modo de arredondamento, temos duas situações possíveis: ou uma “locomotiva a vapor” de zeros surge na área do arredondamento, ou uma “locomotiva a vapor” de uns. A tarefa do algoritmo correto para calcular a função transcendental f (x) é refinar o valor desta função até que a precisão exceda o valor do último bit desta "locomotiva", e até que se torne precisamente claro que, como resultado de flutuações subsequentes do algoritmo numérico para calcular f (x) zeros não se transformarão em uns, ou vice-versa. Assim que tudo se estabilizar e o algoritmo atingir uma precisão que ultrapassa os limites da “locomotiva a vapor”, podemos terminar como se tivéssemos um número infinito de bits. E esse arredondamento ficará com o último bit correto. Mas como isso pode ser alcançado?

"Muletas"

Conforme mencionado, o principal problema é fazer com que o algoritmo supere a locomotiva de zeros ou uns que vem imediatamente após o bit de arredondamento. Quando a locomotiva é vencida e a vemos como um todo, isso equivale ao fato de que esses zeros ou uns já foram calculados com exatidão , e já sabemos exatamente em que direção o arredondamento ocorrerá agora. Mas se não sabemos o comprimento da locomotiva, como podemos projetar um algoritmo?

A primeira "muleta"

Pode parecer ao leitor que a resposta é óbvia: pegue a aritmética com precisão infinita e coloque um número deliberadamente excessivo de bits e, se não for suficiente, coloque outra e recalcule. Em geral, está correto. Isso é feito quando a velocidade e os recursos do computador não desempenham um papel especial. Essa abordagem tem um nome: estratégia multinível de Ziv [1, seção 12.3]. Sua essência é extremamente simples. O algoritmo deve suportar cálculos em vários níveis: um cálculo preliminar rápido (na maioria dos casos, acaba sendo final), um cálculo mais lento, porém mais preciso (salva na maioria dos casos críticos), ainda mais lento, mas um cálculo ainda mais preciso (quando é absolutamente "ruim "Tive que fazer isso) e assim por diante.



Na esmagadora maioria dos casos, é suficiente levar a precisão ligeiramente superior a 0,5ulp, mas se aparecer uma "locomotiva", então a aumentamos. Enquanto a "locomotiva a vapor" permanecer, aumentamos a precisão até que fique estritamente claro que outras flutuações do método numérico não afetarão esta "locomotiva a vapor". Assim, por exemplo, em nosso caso, se alcançamos ESP = 2 ^-54 , então na 54ª posição aparece uma unidade, que, por assim dizer, "protege" a locomotiva de zeros e garante que nenhuma subtração de um valor maior ou igual a ^2-53 e os zeros não se transformarão em uns, arrastando o bit de arredondamento para zero com ele.



Foi uma apresentação científica popular, mesmo assim com o teste de arredondamento de Ziv, onde é mostrado a rapidez, em uma etapa, para verificar se alcançamos a precisão desejada, pode ser lida em [1, Capítulo 12], ou em [3, Seção 10,5].



O problema com essa abordagem é óbvio. É necessário projetar um algoritmo para o cálculo de cada função transcendental f (x) de modo que no decorrer da peça seja possível aumentar a precisão dos cálculos. Para implementação de software, isso ainda não é tão assustador, por exemplo, o método de Newton permite, grosso modo, dobrar o número de bits exatos após a vírgula decimal em cada iteração. Você pode dobrar até que se torne "suficiente", embora seja um processo bastante demorado, devo admitir que o método de Newton nem sempre se justifica, pois requer o cálculo da função inversa f ^-1(x), que em alguns casos pode não ser mais simples do que calcular o próprio f (x). Para implementação de hardware, a "estratégia Ziva" é completamente inadequada. O algoritmo, conectado ao processador, deve executar uma série de ações com o número de bits já predefinido, e isso é bastante problemático de implementar se não soubermos esse número com antecedência. Faça um balanço? Quanto?



A abordagem probabilística para resolver o problema [1, Seção 12.6] nos permite estimar o valor de m (lembre-se, esse é o número de bits, que é suficiente para o arredondamento correto). Acontece que o comprimento da "locomotiva" no sentido probabilístico é ligeiramente maior do que o comprimento da mantissa do número. Assim, na maioria dos casos, será suficiente levar m um pouco mais do que o dobro do valor da mantissa, e apenas em casos muito raros será necessário levar ainda mais. Cito os autores deste trabalho: "deduzimos que, na prática, m deve ser ligeiramente maior que 2p" (eles têm p - o comprimento da mantissa junto com a parte inteira, ou seja, p = 24 para float). Mais adiante no texto, eles mostram que a probabilidade de erro com tal estratégia é próxima de zero, mas ainda positiva, e isso é confirmado por experimentos.



No entanto, ainda existem casos em que o valor de m deve ser considerado ainda mais, e o pior caso não é conhecido com antecedência. Existem estimativas teóricas do pior caso [1, seção 12.7.2], mas elas geram milhões de bits impensáveis, o que não é bom. Aqui está uma tabela do trabalho citado (isto é para a função exp (x) no intervalo de -ln (2) a ln (2)):

Segunda "muleta"

Na prática, m não será tão grande. E para determinar o pior caso, uma segunda "muleta" é aplicada, que é chamada de "pré-cálculo exaustivo". Para o tipo de dados float (32 bits), se a função f tiver um argumento (x), então podemos facilmente "executar" todos os valores possíveis de x. O problema surgirá apenas com funções que têm mais de um argumento (entre elas pow (x, y)), para as quais não poderíamos pensar em nada parecido. Depois de verificar todos os valores possíveis de x, calculamos nossa constante m para cada função f (x) e para cada modo de arredondamento. Em seguida, os algoritmos de cálculo que precisam ser implementados em hardware são projetados para fornecer uma precisão de 2 ^m . Então, o arredondamento f (x) é garantido como correto em todos os casos.



Para o tipo duplo (64 bits), a enumeração simples é quase impossível. No entanto, eles estão se resolvendo! Mas como? A resposta é dada em [4]. Vou falar sobre isso brevemente.



O domínio da função f (x) é dividido em segmentos muito pequenos de modo que dentro de cada segmento seja possível substituir f (x) por uma função linear da forma b-ax (os coeficientes aeb são, é claro, diferentes para diferentes segmentos). O tamanho desses segmentos é calculado analiticamente de forma que tal função linear seria de fato quase indistinguível do original em cada segmento.



Então, após algumas operações de escalonamento e deslocamento, chegamos ao seguinte problema: pode uma linha reta b-ax ir "perto o suficiente" de um ponto inteiro?



Acontece que é relativamente fácil responder sim ou não. Ou seja, "sim" - se pontos potencialmente perigosos estão próximos de uma linha reta, e "não" - se nenhum ponto, em princípio, pode se aproximar da linha. A beleza desse método é que a resposta "não" na prática é obtida na esmagadora maioria dos casos, e a resposta "sim", raramente obtida, obriga a percorrer o segmento com uma busca exaustiva para determinar quais pontos específicos se revelaram críticos.



No entanto, a iteração sobre os argumentos de f (x) é reduzida muitas e muitas vezes e torna possível detectar pontos de quebra para números como double (binary64) e long double (80 bits!). Isso é feito em supercomputadores e, claro, placas de vídeo ... em seu tempo livre de mineração. No entanto, ninguém sabe ainda o que fazer com o tipo de dados binary128. Deixe-me lembrá-lo de que a parte fracionária da mantissa desses números é de 112 bits . Portanto, na literatura estrangeira sobre o assunto até agora, pode-se encontrar apenas raciocínios semifilosóficos começando com “esperamos ...” (“esperamos ...”).



Os detalhes do método que permite determinar rapidamente a passagem de uma linha perto de pontos inteiros são inadequados aqui. Para aqueles que desejam aprender o processo, recomendo olhar com mais cuidado para o problema de encontrar a distância entre uma linha reta e Z ² , por exemplo, no artigo [5]. Ele descreve um algoritmo aprimorado, que no decorrer da construção se assemelha ao famoso algoritmo de Euclides para encontrar o maior divisor comum. Darei a mesma imagem de [4] e [5], que mostra a transformação posterior do problema:







Existem tabelas extensas contendo os piores casos de arredondamento em intervalos diferentes para cada função transcendental. Eles são encontrados em [1 seção 12.8.4] e em [3, seção 10.5.3.2], bem como em artigos separados, por exemplo, em [6].



Vou dar alguns exemplos pegando linhas aleatórias de tais tabelas. Enfatizo que estes não são os piores casos para todos os x, mas apenas para alguns pequenos intervalos, consulte a fonte se estiver interessado.







Como ler a tabela? O valor x é especificado em notação dupla de ponto flutuante hexadecimal. Primeiro, como esperado, há um líder, depois 52 bits da parte fracionária da mantissa e a letra P. Essa letra significa "multiplique por dois elevado a uma potência" seguido por um grau. Por exemplo, P-23 significa que a mantissa especificada precisa ser multiplicada por 2 ^-23 .



Além disso, imagine que a função f (x) é calculada com precisão infinita e os primeiros 53 bits são cortados dela (sem arredondamento!). São esses 53 bits (um deles até a vírgula) que são indicados na coluna f (x). Os bits subsequentes são indicados na última coluna. O sinal de "grau" da sequência de bits na última coluna significa o número de repetições de bits, ou seja, por exemplo, 10 ⁵³1011 significa que primeiro há um bit igual a 1, depois 53 zeros e depois 1011. Em seguida, três pontos, o que significa que, em geral, não precisamos do resto dos bits.



Além disso, é uma questão de tecnologia - conhecemos os piores casos para cada intervalo de uma função tomada separadamente e podemos escolher para esse intervalo tal aproximação de modo que cubra o pior caso com sua precisão. Com apenas anos de computação de supercomputador, é possível criar implementações de hardware rápidas e precisas de funções elementares. A questão é pequena: falta ensinar pelo menos os desenvolvedores de compiladores a usar essas tabelas.

Por que isso é necessário?

Ótima pergunta! Afinal, eu falei repetidamente acima que quase 100% dos programadores não precisam saber uma função elementar para dentro do último bit arredondado corretamente (muitas vezes eles nem precisam da metade dos bits), por que os cientistas dirigem supercomputadores e compilam tabelas para resolver um problema "inútil"?



Primeiro, o desafio é fundamental. É bastante interessante não obter um arredondamento exato para obter um arredondamento preciso, mas, em princípio, para entender como esse interessante problema poderia ser resolvido, que segredos da matemática computacional sua solução nos revelará? Como esses segredos podem ser usados ​​em outras tarefas? Ciências fundamentais - elas são assim, você pode fazer algum tipo de "bobagem" por décadas, e então cem anos depois, graças a essa "bobagem", ocorre um avanço científico em alguma outra área.



Em segundo lugar, a questão da portabilidade do código. Se uma função pode lidar com os últimos bits do resultado da maneira que quiser, isso significa que em diferentes plataformas e em diferentes compiladores, resultados ligeiramente diferentes podem ser obtidos (mesmo se dentro do erro especificado). Em alguns casos, isso não é importante, mas em alguns pode ser significativo, especialmente quando o programa contém um erro que aparece em uma plataforma, mas não aparece em outra plataforma precisamente por causa dos bits diferentes do resultado. Mas por que estou descrevendo a você a conhecida dor de cabeça associada a diferentes comportamentos de programas? Você sabe de tudo isso sem mim. Seria ótimo ter um sistema matemático que funcionasse exatamente da mesma forma em todas as plataformas, não importa o quão compilado esteja. Isso é o que você precisa fazer corretamente terminar a última parte.

Lista de fontes

[1] Jean-Michel Muller, "Elementary Functions: Algorithms and Implementation", 2016



[2] William Kahan, " A Logarithm Too Clever by Half ", 2004



[3] Jean-Michel Muller, "Handbook of floating-point arithmetic" , 2018



[4] Vincent Lefèvre, Jean-Michel Muller, “Toward Correctly Rounded Transcendentals”, IEEE TRANSACTIONS ON COMPUTERS, VOL. 47, NO. 11, NOVEMBRO 1998. pp. 1235-1243



[5] Vincent Lefèvre. “Novos resultados sobre a distância entre um segmento e Z ² ”. Aplicação ao Arredondamento Exato. 17º Simpósio IEEE em Aritmética Computacional - Arith'17, junho de 2005, Cape Cod, MA,

Estados Unidos. pp. 68-75



[6] Vincent Lefèvre, Jean-Michel Muller, “Piores Casos para o Arredondamento Correto das Funções Elementares em Dupla Precisão”, Rapport de recherche (INSTITUT NA TIONAL DE RECHERCHE EN INFORMA TIQUE ET EN AUTOMA TIQUE) n˚4044 - novembro de 2000 - 19 páginas.