Como C computa sin () e outras funções matemáticas?

Eu tenho estudado as desmontagens do .NET e o código-fonte do GCC, mas parece que não consigo encontrar em lugar algum a implementação real do sin() e outras funções matemáticas … elas sempre parecem estar fazendo referência a outra coisa.

Alguém pode me ajudar a encontrá-los? Eu sinto que é improvável que todo o hardware que C irá rodar suporte funções trigonométricas no hardware, então deve haver um algoritmo de software em algum lugar , certo?


Estou ciente de várias maneiras pelas quais as funções podem ser calculadas e escrevi minhas próprias rotinas para calcular funções usando a série taylor para diversão. Estou curioso sobre como as linguagens de produção reais fazem isso, já que todas as minhas implementações são sempre várias ordens de magnitude mais lentas, embora eu ache que meus algoritmos são bem inteligentes (obviamente, não são).

No GNU libm, a implementação do sin é dependente do sistema. Portanto, você pode encontrar a implementação, para cada plataforma, em algum lugar no subdiretório apropriado de sysdeps .

Um diretório inclui uma implementação em C, contribuída pela IBM. Desde outubro de 2011, esse é o código que realmente é executado quando você chama sin() em um sistema Linux x86-64 típico. É aparentemente mais rápido que a fsin assembly do fsin . Código fonte: sysdeps / ieee754 / dbl-64 / s_sin.c , procure por __sin (double x) .

Esse código é muito complexo. Nenhum algoritmo de software é tão rápido quanto possível e também precisa em todo o intervalo de valores de x , portanto, a biblioteca implementa muitos algoritmos diferentes e seu primeiro trabalho é examinar xe decidir qual algoritmo usar. Em algumas regiões, usa o que parece ser a série familiar de Taylor. Vários algoritmos calculam primeiro um resultado rápido, então, se isso não for preciso, descartá-lo e recorrer a um algoritmo mais lento.

Versões de 32 bits mais antigas do GCC / glibc usavam a instrução fsin , que é surpreendentemente imprecisa para algumas inputs. Há um post fascinante no blog ilustrando isso com apenas 2 linhas de código .

A implementação do sin pelo fdlibm em C puro é muito mais simples que a da glibc e é bem comentada. Código fonte: fdlibm / s_sin.c e fdlibm / k_sin.c

OK kiddies, tempo para os profissionais …. Esta é uma das minhas maiores queixas com engenheiros de software inexperientes. Eles vêm calculando funções transcendentais a partir do zero (usando a série de Taylor) como se ninguém tivesse feito esses cálculos antes em suas vidas. Não é verdade. Este é um problema bem definido e tem sido abordado milhares de vezes por engenheiros de software e hardware muito inteligentes e tem uma solução bem definida. Basicamente, a maioria das funções transcendentais usa polinômios de Chebyshev para calculá-los. Quanto a quais polinômios são usados ​​depende das circunstâncias. Primeiro, a bíblia sobre esse assunto é um livro chamado “Computer Approximations”, de Hart e Cheney. Nesse livro, você pode decidir se tem um sumdor de hardware, multiplicador, divisor etc. e decidir quais operações são mais rápidas. Por exemplo, se você tiver um divisor realmente rápido, o caminho mais rápido para calcular o seno pode ser P1 (x) / P2 (x) onde P1, P2 são polinômios de Chebyshev. Sem o divisor rápido, pode ser apenas P (x), onde P tem muito mais termos que P1 ou P2 …. então seria mais lento. Então, o primeiro passo é determinar seu hardware e o que ele pode fazer. Então você escolhe a combinação apropriada de polinômios de Chebyshev (geralmente é da forma cos (ax) = aP (x) para cosseno por exemplo, novamente onde P é um polinômio de Chebyshev). Então você decide qual precisão decimal você quer. Por exemplo, se você quer 7 dígitos de precisão, você olha para cima na tabela apropriada no livro que mencionei, e ele lhe dará (para precisão = 7,33) um número N = 4 e um número polinomial 3502. N é a ordem do polinomial (por isso é p4.x ^ 4 + p3.x ^ 3 + p2.x ^ 2 + p1.x + p0), porque N = 4. Então você procura o valor real dos valores de p4, p3, p2, p1, p0 na parte de trás do livro abaixo de 3502 (eles estarão em ponto flutuante). Então você implementa seu algoritmo no software no formato: (((p4.x + p3) .x + p2) .x + p1) .x + p0 …. e é assim que você calcula o co-seno para 7 casas decimais coloca nesse hardware.

Observe que a maioria das implementações de hardware de operações transcendentais em uma FPU geralmente envolve algum microcódigo e operações como essa (depende do hardware). Os polinômios de Chebyshev são usados ​​para a maioria dos transcendentais, mas não para todos. Por exemplo, a raiz quadrada é mais rápida para usar uma iteração dupla do método Newton Raphson usando uma tabela de pesquisa primeiro. Mais uma vez, o livro “Computer Approximations” (Aproximações de Computador) lhe dirá isso.

Se você planeja implantar essas funções, recomendo a qualquer pessoa que obtenha uma cópia desse livro. É realmente a bíblia para esses tipos de algoritmos. Observe que existem vários meios alternativos para calcular esses valores, como cordics, etc, mas eles tendem a ser melhores para algoritmos específicos, nos quais você precisa apenas de baixa precisão. Para garantir a precisão toda vez, os polinômios chebyshev são o caminho a percorrer. Como eu disse, problema bem definido. Foi resolvido por 50 anos agora ….. e é assim que é feito.

Agora, dito isto, há técnicas pelas quais os polinômios de Chebyshev podem ser usados ​​para obter um único resultado de precisão com um polinômio de baixo grau (como o exemplo do cosseno acima). Em seguida, existem outras técnicas para interpolar entre valores para aumentar a precisão sem ter que ir para um polinômio muito maior, como “Método de tabelas precisas de Gal”. Esta última técnica é o que o post referente à literatura do ACM está se referindo. Mas, em última análise, os polinômios de Chebyshev são usados ​​para chegar a 90% do caminho até lá.

Apreciar.

Funções como seno e cosseno são implementadas em microcódigo dentro de microprocessadores. Chips da Intel, por exemplo, possuem instruções de assembly para estes. O compilador AC gerará código que chama essas instruções de assembly. (Por outro lado, um compilador Java não. O Java avalia as funções trigonométricas no software, e não no hardware, e por isso é executado muito mais lentamente.)

Chips não usam a série de Taylor para computar funções trigonométricas, pelo menos não inteiramente. Primeiro de tudo eles usam CORDIC , mas eles também podem usar uma pequena série de Taylor para polir o resultado de CORDIC ou para casos especiais, como computação senoidal com alta precisão relativa para ângulos muito pequenos. Para mais explicações, veja esta resposta do StackOverflow .

Sim, existem algoritmos de software para calcular o sin também. Basicamente, calcular esse tipo de coisa com um computador digital geralmente é feito usando methods numéricos como aproximar a série de Taylor representando a function.

Os methods numéricos podem aproximar as funções a uma quantidade arbitrária de precisão e, como a quantidade de precisão que você tem em um número flutuante é finita, elas atendem bem a essas tarefas.

É uma questão complexa. A CPU da família x86 da Intel tem uma implementação de hardware da function sin() , mas é parte da FPU x87 e não é mais usada no modo de 64 bits (onde os registros SSE2 são usados). Nesse modo, uma implementação de software é usada.

Existem várias implementações desse tipo por aí. Um está no fdlibm e é usado em Java. Tanto quanto sei, a implementação do glibc contém partes do fdlibm e outras partes contribuídas pela IBM.

Implementações de software de funções transcendentais, como sin() tipicamente usam aproximações por polinômios, freqüentemente obtidas da série de Taylor.

use a série taylor e tente encontrar a relação entre os termos da série para que você não calcule as coisas de novo e de novo

Aqui está um exemplo para cosinus:

 double cosinus(double x,double prec) { double t , s ; int p; p = 0; s = 1.0; t = 1.0; while(fabs(t/s) > prec) { p++; t = (-t * x * x) / ((2 * p - 1) * (2 * p)); s += t; } return s;} 

usando isso, podemos obter o novo termo da sum usando o já usado (evitamos o fatorial e x ^ 2p)

explicação http://img514.imageshack.us/img514/1959/82098830.jpg

Os polinômios de Chebyshev, como mencionado em outra resposta, são os polinômios onde a maior diferença entre a function e o polinômio é a menor possível. Esse é um excelente começo.

Em alguns casos, o erro máximo não é o que você está interessado, mas o erro relativo máximo. Por exemplo, para a function seno, o erro próximo de x = 0 deve ser muito menor que para valores maiores; você quer um pequeno erro relativo . Então, você calcularia o polinômio de Chebyshev para sin x / x e multiplicaria esse polinômio por x.

Em seguida, você tem que descobrir como avaliar o polinômio. Você quer avaliá-lo de tal forma que os valores intermediários sejam pequenos e, portanto, os erros de arredondamento sejam pequenos. Caso contrário, os erros de arredondamento podem se tornar muito maiores do que erros no polinômio. E com funções como a function seno, se você é descuidado, então pode ser possível que o resultado que você calculou para o sin x seja maior que o resultado para o sin y mesmo quando x

Por exemplo, sen x = x – x ^ 3/6 + x ^ 5/120 – x ^ 7/5040 … Se você calcular ingenuamente sen x = x * (1 – x ^ 2/6 + x ^ 4 / 120 – x ^ 6/5040 …), então essa function entre parênteses está diminuindo, e acontecerá que se y for o próximo número maior para x, então às vezes sen y será menor que sin x. Em vez disso, calcule sin x = x – x ^ 3 * (1/6 – x ^ 2/120 + x ^ 4/5040 …) onde isso não pode acontecer.

Ao calcular os polinômios de Chebyshev, geralmente é necessário arredondar os coeficientes para precisão dupla, por exemplo. Mas enquanto um polinômio de Chebyshev é ótimo, o polinômio de Chebyshev com coeficientes arredondados para precisão dupla não é o polinômio ideal com coeficientes de precisão dupla!

Por exemplo, para sin (x), onde você precisa de coeficientes para x, x ^ 3, x ^ 5, x ^ 7 etc., faça o seguinte: Calcule a melhor aproximação de sin x com um polinômio (ax + bx ^ 3 + cx ^ 5 + dx ^ 7) com precisão maior que o dobro, então arredonde para precisão dupla, dando A. A diferença entre a e A seria bem grande. Agora calcule a melhor aproximação de (sen x – Ax) com um polinômio (bx ^ 3 + cx ^ 5 + dx ^ 7). Você obtém coeficientes diferentes, porque eles se adaptam à diferença entre a e A. Rodada b para precisão dupla B. Então aproximada (sin x – Ax – Bx ^ 3) com um polinômio cx ^ 5 + dx ^ 7 e assim por diante. Você obterá um polinômio que é quase tão bom quanto o polinômio Chebyshev original, mas muito melhor do que Chebyshev arredondado para precisão dupla.

Em seguida, você deve levar em conta os erros de arredondamento na escolha do polinômio. Você encontrou um polinômio com erro mínimo no polinômio ignorando o erro de arredondamento, mas deseja otimizar o erro polinomial mais arredondamento. Depois de ter o polinômio Chebyshev, você pode calcular limites para o erro de arredondamento. Diga que f (x) é sua function, P (x) é o polinômio e E (x) é o erro de arredondamento. Você não quer otimizar | f (x) – P (x) |, você quer otimizar | f (x) – P (x) +/- E (x) | Você obterá um polinômio ligeiramente diferente que tenta manter os erros polinomiais abaixo do local onde o erro de arredondamento é grande e relaxa os erros do polinômio um pouco quando o erro de arredondamento é pequeno.

Tudo isso fará você arredondar erros de no máximo 0,55 vezes o último bit, onde +, -, *, / têm erros de arredondamento de no máximo 0,50 vezes o último bit.

Para o sin especificamente, usar a expansão de Taylor lhe daria:

sin (x): = x – x ^ 3/3! + x ^ 5/5! – x ^ 7/7! + … (1)

você continuaria adicionando termos até que a diferença entre eles fosse menor do que um nível de tolerância aceito ou apenas por uma quantidade finita de etapas (mais rápido, mas menos preciso). Um exemplo seria algo como:

 float sin(float x) { float res=0, pow=x, fact=1; for(int i=0; i<5; ++i) { res+=pow/fact; pow*=-1*x*x; fact*=(2*(i+1))*(2*(i+1)+1); } return res; } 

Nota: (1) funciona por causa da aproximação sin (x) = x para pequenos ângulos. Para ângulos maiores, você precisa calcular mais e mais termos para obter resultados aceitáveis. Você pode usar um argumento while e continuar com uma certa precisão:

 double sin (double x){ int i = 1; double cur = x; double acc = 1; double fact= 1; double pow = x; while (fabs(acc) > .00000001 && i < 100){ fact *= ((2*i)*(2*i+1)); pow *= -1 * x*x; acc = pow / fact; cur += acc; i++; } return cur; } 

A implementação real das funções da biblioteca depende do compilador específico e / ou do provedor da biblioteca. Quer seja feito em hardware ou software, seja uma expansão de Taylor ou não, etc., irá variar.

Eu percebo que absolutamente não ajuda.

Em relação à function trigonométrica como sin() , cos() , tan() não houve menção, após 5 anos, de um aspecto importante das funções trigonométricas de alta qualidade: Redução do alcance .

Um passo inicial em qualquer uma dessas funções é reduzir o ângulo, em radianos, para um intervalo de 2 * π. Mas π é irracional, então reduções simples como x = remainder(x, 2*M_PI) introduzem erro, pois M_PI , ou máquina pi, é uma aproximação de π. Então, como fazer x = remainder(x, 2*π) ?

As primeiras bibliotecas usavam a precisão ampliada ou a programação criada para dar resultados de qualidade, mas ainda em um intervalo limitado de double . Quando um grande valor foi solicitado como sin(pow(2,30)) , os resultados foram sem sentido ou 0.0 e talvez com um sinalizador de erro definido como algo como perda total de precisão de PLOSS ou PLOSS perda parcial de precisão.

Uma boa redução de intervalo de grandes valores para um intervalo como -π a π é um problema desafiador que rivaliza com os desafios da function trigonométrica básica, como o próprio sin() .

Um bom relatório é a redução de argumentos para grandes argumentos: Good to the last bit (1992). Aborda bem a questão: discute a necessidade e como as coisas estavam em várias plataformas (SPARC, PC, HP, 30+ outras) e fornece um algoritmo de solução que fornece resultados de qualidade para todos os double de -DBL_MAX a DBL_MAX .


Se os argumentos originais estiverem em graus, ainda que possam ser de grande valor, use fmod() primeiro para maior precisão. Um bom fmod() não apresentará erros e, portanto, fornecerá uma excelente redução de faixa.

 // sin(degrees2radians(x)) sin(degrees2radians(fmod(x, 360.0))); // -360.0 < = fmod(x,360) <= +360.0 

Várias identidades trigonométricas e remquo() oferecem ainda mais melhorias. Amostra: sind ()

Eles são normalmente implementados em software e não usarão as chamadas de hardware correspondentes (ou seja, na maioria dos casos). No entanto, como Jason apontou, estas são específicas da implementação.

Observe que essas rotinas de software não fazem parte das fonts do compilador, mas serão encontradas na biblioteca correspoding, como o clib, ou glibc para o compilador GNU. Veja http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions

Se você quer um maior controle, você deve avaliar cuidadosamente o que você precisa exatamente. Alguns dos methods típicos são a interpolação de tabelas de consulta, a chamada de assembly (que geralmente é lenta) ou outros esquemas de aproximação, como o Newton-Raphson, para raízes quadradas.

Se você quer uma implementação em software, não em hardware, o lugar para procurar uma resposta definitiva para essa questão é o Capítulo 5 das Receitas Numéricas . Minha cópia está em uma checkbox, então não posso dar detalhes, mas a versão curta (se bem me lembro disso) é que você toma tan(theta/2) como sua operação primitiva e computa os outros a partir daí. A computação é feita com uma aproximação de série, mas é algo que converge muito mais rapidamente que uma série de Taylor.

Desculpe, não posso me lembrar mais sem colocar minha mão no livro.

Como muitas pessoas apontaram, é dependente da implementação. Mas, tanto quanto eu entendi sua pergunta, você estava interessado em uma implementação de software real de funções matemáticas, mas simplesmente não conseguiu encontrar um. Se este for o caso, então você está aqui:

  • Faça o download do código fonte da glibc em http://ftp.gnu.org/gnu/glibc/
  • Veja o arquivo dosincos.c localizado na dosincos.c descompactada glibc root \ sysdeps \ ieee754 \ dbl-64
  • Da mesma forma você pode encontrar implementações do resto da biblioteca de matemática, basta procurar o arquivo com o nome apropriado

Você também pode dar uma olhada nos arquivos com a extensão .tbl , seus conteúdos nada mais são do que grandes tabelas de valores pré-computados de diferentes funções em uma forma binária. É por isso que a implementação é tão rápida: em vez de computar todos os coeficientes de qualquer série que eles usem, eles simplesmente fazem uma pesquisa rápida, que é muito mais rápida. BTW, eles usam série Tailor para calcular seno e cosseno.

Eu espero que isso ajude.

Vou tentar responder para o caso de sin() em um programa C, compilado com o compilador C do GCC em um processador x86 atual (digamos um Intel Core 2 Duo).

Na linguagem C, a Biblioteca Padrão C inclui funções matemáticas comuns, não incluídas na própria linguagem (por exemplo, pow , sin e cos para power, sine e cosine respectivamente). Os headers estão incluídos em math.h.

Agora em um sistema GNU / Linux, essas funções de bibliotecas são fornecidas pela glibc (GNU libc ou GNU C Library). Mas o compilador GCC quer que você vincule à biblioteca de matemática ( libm.so ) usando o -lm compilador -lm para habilitar o uso dessas funções matemáticas. Não sei por que não faz parte da biblioteca C padrão. Estas seriam uma versão de software das funções de ponto flutuante, ou “soft-float”.

Aparte: A razão para ter as funções matemáticas separadas é histórica, e foi meramente destinada a reduzir o tamanho dos programas executáveis ​​em sistemas Unix muito antigos, possivelmente antes que as bibliotecas compartilhadas estivessem disponíveis, tanto quanto eu sei.

Agora o compilador pode otimizar a function padrão da biblioteca C sin() (fornecida pelo libm.so ) para ser substituída por uma chamada para uma instrução nativa para a function sin() libm.so CPU / FPU, que existe como uma instrução FPU ( FSIN para x86 / x87) em processadores mais recentes como a série Core 2 (isso está correto até o i486DX). Isso dependeria de sinalizadores de otimização passados ​​para o compilador gcc. Se o compilador foi instruído a escrever código que seria executado em qualquer processador i386 ou mais recente, ele não faria tal otimização. O -mcpu=486 informaria ao compilador que era seguro fazer tal otimização.

Agora, se o programa executasse a versão do software da function sin (), o faria baseado em um algoritmo CORDIC (Computador de Rotação Coordenada de Coordenadas) ou BKM , ou mais provavelmente um cálculo de tabela ou série de energia que é comumente usado agora para calcular tais funções transcendentais. [Src: http://en.wikipedia.org/wiki/Cordic#Application%5D

Qualquer versão recente (desde 2.9x aprox.) Do gcc também oferece uma versão __builtin_sin() do sin, __builtin_sin() que será usada para replace a chamada padrão para a versão da biblioteca C, como uma otimização.

Tenho certeza de que isso é tão claro quanto a lama, mas espero que você tenha mais informações do que esperava, e muitos pontos para aprender mais.

Não há nada como acertar a fonte e ver como alguém realmente fez isso em uma biblioteca de uso comum; Vamos dar uma olhada em uma implementação da biblioteca C em particular. Eu escolhi o uLibC.

Aqui está a function do pecado:

http://git.uclibc.org/uClibc/tree/libm/s_sin.c

que parece lidar com alguns casos especiais e, em seguida, realiza alguma redução de argumento para mapear a input para o intervalo [-pi / 4, pi / 4], (dividindo o argumento em duas partes, uma grande parte e uma cauda) antes de ligar

http://git.uclibc.org/uClibc/tree/libm/k_sin.c

que então opera nessas duas partes. Se não houver cauda, ​​uma resposta aproximada é gerada usando um polinômio de grau 13. Se houver uma cauda, ​​você recebe uma pequena adição corretiva baseada no princípio de que sin(x+y) = sin(x) + sin'(x')y

Sempre que tal function é avaliada, então, em algum nível, é mais provável que:

  • Uma tabela de valores que é interpolada (para aplicações rápidas e imprecisas – por exemplo, computação gráfica)
  • A avaliação de uma série que converge para o valor desejado – provavelmente não é uma série taylor, mais provavelmente algo baseado em uma quadratura de fantasia como Clenshaw-Curtis.

Se não houver suporte de hardware, o compilador provavelmente usará o último método, emitindo apenas código assembler (sem símbolos de debugging), em vez de usar a biblioteca ac – tornando complicado rastrear o código real no depurador.

Se você quiser olhar para a implementação real do GNU dessas funções em C, confira o último tronco da glibc. Veja a biblioteca C GNU .

Computação seno / cosseno / tangente é realmente muito fácil de fazer através do código usando a série de Taylor. Escrever um você mesmo leva 5 segundos.

Todo o processo pode ser resumido com esta equação aqui: http://sofpt.miximages.com/c/546ecab719ce73dfb34a7496c942972b.png

Aqui estão algumas rotinas que escrevi para C:

 double _pow(double a, double b) { double c = 1; for (int i=0; i 

Don’t use Taylor series. Chebyshev polynomials are both faster and more accurate, as pointed out by a couple of people above. Here is an implementation (originally from the ZX Spectrum ROM): https://albertveli.wordpress.com/2015/01/10/zx-sine/

if you want sin then asm volatile (“fsin” : “=t”(vsin) : “0”(xrads)); if you want cos then asm volatile (“fcos” : “=t”(vcos) : “0”(xrads)); if you want sqrt then asm volatile (“fsqrt” : “=t”(vsqrt) : “0”(value)); so why use inaccurate code when the machine instructions will do.