Há algum tempo eu escrevi um código que tentava calcular o sem usar as funções da biblioteca. Ontem, eu estava revisando o código antigo e tentei torná-lo o mais rápido possível (e correto). Aqui está minha tentativa até agora:
const double ee = exp(1);
double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 )
n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(1 - x) = -x - x**2/2 - x**3/3... */
n = 1 - n;
now = term = n;
for ( i = 1 ; ; ){
lgVal -= now;
term *= n;
now = term / ++i;
if ( now < 1e-17 ) break;
}
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Aqui, estou tentando encontrar para que acabe de n e, em seguida, adiciono o valor do logaritmo de , que é menor que 1. Neste ponto, a expansão do de Taylor pode ser usada sem se preocupar.n log(1-x)
Recentemente, interessei-me pela análise numérica, e é por isso que não consigo deixar de fazer a pergunta: quanto mais rápido esse segmento de código pode ser executado na prática, enquanto está correto o suficiente? Preciso mudar para outros métodos, por exemplo, usando a fração continuada, como esta ?
A função fornecida com a biblioteca padrão C é quase 5,1 vezes mais rápida que esta implementação.
ATUALIZAÇÃO 1 : Usando a série arctan hiperbólica mencionada na Wikipedia , o cálculo parece quase 2,2 vezes mais lento que a função de log de biblioteca padrão C. No entanto, não verifiquei extensivamente o desempenho e, para números maiores, minha implementação atual parece REALMENTE lenta. Quero verificar minha implementação quanto ao limite de erros e ao tempo médio para obter uma ampla variedade de números, se eu puder gerenciar. Aqui está o meu segundo esforço.
double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
for ( i = 3 ; ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Qualquer sugestão ou crítica é apreciada.
ATUALIZAÇÃO 2: Com base nas sugestões feitas abaixo, adicionei algumas mudanças incrementais aqui, que são cerca de 2,5 vezes mais lentas que a implementação da biblioteca padrão. No entanto, eu testei apenas para números inteiros desta vez, para números maiores, o tempo de execução aumentaria. Por enquanto. Ainda não conheço técnicas para gerar números duplos aleatórios , por isso ainda não está totalmente aferido. Para tornar o código mais robusto, adicionei correções para casos de canto. O erro médio para os testes que fiz é de cerca de . ≤ 1 e 308 4 e - 15
double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n == 0 ) return -1./0.; /* -inf */
if ( n < 0 ) return 0./0.; /* NaN*/
if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
/* the cutoff iteration is 650, as over e**650, term multiplication would
overflow. For larger numbers, the loop dominates the arctanh approximation
loop (with having 13-15 iterations on average for tested numbers so far */
for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
if ( lgVal == 650 ){
n /= term;
for ( term = 1 ; term < n ; term *= ee, lgVal++ );
}
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
/* limiting the iteration for worst case scenario, maximum 24 iteration */
for ( i = 3 ; i < 50 ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
fonte
frexppara obter a mantissa e o expoente do número , de modo que .A resposta de Kirill já abordou um grande número de questões relevantes. Gostaria de expandir alguns deles com base na experiência prática em design de bibliotecas de matemática. Uma observação adiantada: os projetistas de bibliotecas de matemática tendem a usar todas as otimizações algorítmicas publicadas, bem como muitas otimizações específicas de máquinas, nem todas serão publicadas. O código freqüentemente será escrito em linguagem assembly, em vez de usar código compilado. Portanto, é improvável que uma implementação direta e compilada atinja mais de 75% do desempenho de uma implementação existente de biblioteca matemática de alta qualidade, assumindo conjuntos de recursos idênticos (precisão, tratamento de casos especiais, relatório de erros, suporte ao modo de arredondamento).
Em termos de precisão, para funções elementares (por exemplo , , ) e funções especiais simples (por exemplo , , ), o erro ulp substituiu o erro relativo como a métrica de erro relevante. Um objetivo de projeto comum para funções elementares é um erro máximo de menos de 1 ulp, resultando em uma função fielmente arredondada. Uma função fielmente arredondada retorna o número de ponto flutuante mais próximo do resultado infinitamente preciso ou um número de ponto flutuante imediatamente adjacente.exp log erfc Γ
A precisão é normalmente avaliada por comparação com uma referência de alta precisão (de terceiros). As funções de precisão única de argumento único podem ser facilmente testadas exaustivamente; outras funções requerem teste com vetores de teste aleatórios (direcionados). Claramente, não se pode calcular resultados de referência infinitamente precisos, mas a pesquisa no dilema do criador de tabelas sugere que, para muitas funções simples, basta calcular uma referência com uma precisão de cerca de três vezes a precisão do objetivo. Veja por exemplo:
Vincent Lefèvre, Jean-Michel Muller, "Piores casos de arredondamento correto das funções elementares em dupla precisão". Em Proceedings 15th IEEE Symposium on Computer Arithmetic , 2001, 1111-118). (pré-impressão online)
Em termos de desempenho, é preciso distinguir entre otimizar a latência (importante quando se olha o tempo de execução das operações dependentes) e otimizar o throughput (relevante ao considerar o tempo de execução das operações independentes). Nos últimos vinte anos, a proliferação de técnicas de paralelização de hardware, como paralelismo em nível de instrução (por exemplo, superescalar, processadores fora de ordem), paralelismo em nível de dados (por exemplo, instruções SIMD) e paralelismo em nível de encadeamento (por exemplo, hiperencadeamento, processadores multinúcleo) levou a uma ênfase na taxa de transferência computacional como a métrica mais relevante.
Aproximações principais para funções elementares são quase exclusivamente aproximações mínimas . Devido ao custo relativamente alto da divisão, aproximações polinomiais de minimax. Ferramentas como Mathematica ou Maple possuem recursos internos para gerá-los; também existem ferramentas especializadas como o Sollya . Para o logaritmo, as opções básicas de aproximação do núcleo, após a redução do argumento para valores próximos à unidade, são e , em que é uma aproximação polinomial de minimax. Em termos de desempenho, o primeiro é geralmente preferido para implementações de precisão única (consulte esta respostal O g ( x ) = 2 ⋅ um t um n h ( ( x - 1 ) / ( x + 1 ) ) = p ( ( ( x - 1 ) / ( x + 1 ) ) 2 ) plog(1+x)=p(x) log(x)=2⋅atanh((x−1)/(x+1))=p(((x−1)/(x+1))2) p para um exemplo trabalhado), enquanto o último é preferido para implementações de precisão dupla.
A operação de adição múltipla ( FMA ), introduzida pela IBM há 25 anos e agora disponível em todas as principais arquiteturas de processadores, é um componente essencial das implementações modernas de bibliotecas de matemática. Ele fornece redução de erro de arredondamento, oferece proteção limitada contra cancelamento subtrativo e simplifica enormemente a aritmética duplo-duplo .
A233
C99implementação exemplar de precisão dupla IEEE-754log()abaixo demonstra o uso de FMA (expostoC99como afma()função matemática padrão), juntamente com o uso muito limitado de técnicas de duplo duplo para melhorar a precisão de produtos com constantes transcendentais. São fornecidas duas aproximações principais diferentes, ambas com resultados fielmente arredondados, conforme demonstrado pelo teste com vetores de teste aleatórios. As aproximações mínimas utilizadas foram calculadas com minhas próprias ferramentas, com base no algoritmo Remez .fonte