Como É Que C calcula o sin () e outras funções matemáticas?

Muito bem, meninos, está na hora dos profissionais.... Esta é uma das minhas maiores queixas com engenheiros de software inexperientes. Eles vêm calculando funções transcendentais do zero (usando a série de Taylor) como se nunca ninguém tivesse feito esses cálculos antes em suas vidas. Falso. Este é um problema bem definido e foi 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 usam Polinômios de Chebyshev para calculá-los. A respeito de quais polinômios são usados depende das circunstâncias. Primeiro, a Bíblia sobre este assunto é um livro chamado "aproximações de computador" de Hart e Cheney. Nesse livro, você pode decidir se você tem um adder de hardware, multiplicador, divisor, etc, e decidir quais operações são mais rápidas. por exemplo, se você tivesse um divisor realmente rápido, a maneira mais rápida de calcular seno poderia ser P1(x)/P2(x) onde P1, P2 são polinômios de Chebyshev. Sem o separador rápido, pode ser apenas P (x), onde P tem muito mais termos do que P1 ou P2....so seria mais lento. Então, o primeiro passo é determinar o 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 coseno, por exemplo, novamente ONDE P é um polinômio de Chebyshev). Então você decide que precisão decimal você quer. por exemplo, se você quiser precisão de 7 dígitos, você olhar para cima na tabela apropriada no livro que mencionei, e ele lhe dará (para precision = 7.33) a number N = 4 and a polynomial number 3502. N é a ordem do polinômio (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 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 em software na forma: ((p4.x + p3).x + p2).x + p1).x + p0 ....e é assim que se calcula coseno com 7 casas decimais naquele hardware.

Note que a maior parte do hardware implementações de operações transcendentais em uma FPU geralmente envolvem alguns microcódigos e operações como esta (depende do hardware). Polinômios de Chebyshev são usados para a maioria dos transcendentais, mas não todos. por exemplo, raiz quadrada é mais rápido para usar uma iteração dupla do método Newton raphson usando uma tabela de pesquisa primeiro. Mais uma vez, esse livro "aproximações de computador" vai dizer-lhe isso.

Se planeias implantar estas funções, recomendo a qualquer um que receba uma cópia desse livro. Ele realmente é a Bíblia para estes tipos de algoritmos. Note que existem cachos de meios alternativos para calcular esses valores como cordics, etc, mas estes tendem a ser melhores para algoritmos específicos onde você só precisa de baixa precisão. Para garantir a precisão de cada vez, os polinômios de chebyshev são o caminho a seguir. Como eu disse, problema bem definido. Está resolvido há 50 anos.....e é assim que se faz. Dito isto, há técnicas pelas quais o Chebyshev polinômios podem ser usados para obter um único resultado de precisão com um polinômio de baixo grau (como o exemplo para cosseno acima). Então, 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 referindo-se à literatura ACM está se referindo. Mas, em última análise, os polinômios de Chebyshev são o que são usados para obter 90% do caminho la. Divirtam-se.

Funções como seno e cosseno são implementadas em microcódigo dentro de microprocessadores. Chips Intel, por exemplo, têm instruções de montagem para estes. Um compilador C irá gerar código que chama essas instruções de montagem. (Em contraste, um compilador Java não vai. Java avalia funções trig em software ao invés de hardware, e assim ele corre muito mais devagar.)

Os Chips não usam a série Taylor para calcular funções trig, pelo menos não inteiramente. Em primeiro lugar, eles usam CORDIC , mas eles também podem usar uma série de Taylor curta para polir o resultado do CORDIC ou para casos especiais, como a computação de seno com alta precisão relativa para ângulos muito pequenos. Para mais explicações, veja esta resposta de StackOverflow .

Sim, existem algoritmos de software para calcular sin também. Basicamente, calcular este tipo de coisas com um computador digital é normalmente feito usando métodos numéricos como aproximar a série Taylor representando a função.

Métodos numéricos podem aproximar 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, eles se adequam a essas tarefas muito bem.

É uma questão complexa. CPU Intel-like da família x86 tem uma implementação de hardware da função sin(), mas é parte do X87 FPU e não mais usado no modo 64-bit (onde registros SSE2 são usados em vez). Nesse modo, uma implementação de software é usada. Existem várias implementações deste tipo por aí. Um está em fdlibm e é usado em Java. Tanto quanto sei, a implementação da glibc contém partes de fdlibm, e outras partes contribuíram para IBM.

Implementações de Software de funções transcendentais como sin() normalmente usam aproximações por polinômios, frequentemente obtidas a partir de Séries de Taylor.

Use taylor series e tente encontrar uma relação entre os termos da série para que não calcule as coisas uma e outra vez

Eis um exemplo para cosino:

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 isto podemos obter o novo termo da soma usando o já usado (evitamos o factorial e x^2p)

Anotações 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 função e o polinômio é tão pequena quanto possível. É 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 função sine, o erro próximo de x = 0 deve ser muito menor do que para valores maiores; você quer um pequeno erro relativo. Para que calculasse o Polinomial de Chebyshev para sin x / x, e multiplique esse polinomial por X.

A seguir, tens de descobrir como avaliar o polinômio. Pretende avaliá-la de forma a que os valores intermédios sejam pequenos, pelo que os erros de arredondamento são pequenos. Caso contrário, os erros de arredondamento pode se tornar muito maior do que erros no polinômio. E com funções como a função seno, se você é descuidado, então pode ser possível que o resultado que você calcula para sin x é maior do que o resultado para sin y mesmo quando x

Por exemplo, sin x = x-x^3/6 + x^5 / 120-x^7 / 5040... Se calcular ingenuamente sin x = x * (1 - x^2/6 + x^4/120 - x^6/5040...), então essa função entre parênteses está diminuindo, e irá acontecer que se y é o número maior seguinte A x, então às vezes o pecado y será menor que o 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, normalmente é preciso arredondar os coeficientes para dupla precisão, por exemplo. Mas enquanto um polinômio de Chebyshev é ideal, 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. você faz o seguinte: calcular o a melhor aproximação de sin x com uma precisão polinomial (ax + bx^3 + cx^5 + dx^7) com maior precisão do que a dupla, então arredondando a para precisão dupla, dando A. a diferença entre a e A seria bastante grande. Agora calcule a melhor aproximação de (sin x-Ax) com um polinômio (b x^3 + cx^5 + dx^7). Você obtém coeficientes diferentes, porque eles se adaptam à diferença entre a e A. Round b para dupla precisão B. Então aproximado (sin x - Ax - Bx^3) com um polinômio CX^5 + dx^7 e assim por diante. Você vai obter um polinômio que é quase tão bom quanto o polinômio original de Chebyshev, mas muito melhor do que Chebyshev arredondado para precisão dupla.

A seguir deve ter em conta os erros de arredondamento na escolha do polinômio. Você encontrou um polinômio com um erro mínimo no polinômio ignorando o erro de arredondamento, mas você quer otimizar o erro polinomial mais o erro de arredondamento. Uma vez que você tem o polinômio de Chebyshev, você pode calcular os limites para o erro de arredondamento. Diga f (x) é a sua função, P (x) é o polinômio, 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ê vai obter um polinômio ligeiramente diferente que tenta manter os erros polinomiais para baixo onde o erro de arredondamento é grande, e relaxa os erros polinomiais um pouco onde o erro de arredondamento é pequeno.

Tudo isto lhe dará facilmente erros de arredondamento de, no máximo, 0,55 vezes o último pedaço, onde +, -,*, / têm erros de arredondamento de, no máximo, 0,50 vezes a última parte.

Para sin especificamente, usar a expansão Taylor dar-lhe-ia:

Sin (x): = x-x^3/3! + x^5/5! - X^7/7! + ... (1)

Continuaria a adicionar Termos até que a diferença entre eles fosse inferior a um nível de tolerância aceite ou apenas para uma quantidade finita de passos (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 ângulos pequenos. Para maiores ângulos você precisa calcular mais e mais termos para obter resultados aceitáveis. Você pode usar um argumento de tempo e Continuar para 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 efectiva das funções da biblioteca depende do compilador e/ou do fornecedor específico da biblioteca. Seja feito em hardware ou software, seja uma expansão Taylor ou Não, etc., varia.

Sei que isso não ajuda em nada.

Relativa à função trigonométricasin(), cos(),tan() após 5 anos, não foi mencionado um aspecto importante das funções de trigonometria de alta qualidade: redução da Gama.

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

Bibliotecas antigas utilizadas precisão estendida ou programação trabalhada para dar resultados de qualidade, mas ainda sobre uma gama limitada de double. Quando um grande valor foi solicitado como sin(pow(2,30)), os resultados não tinham sentido ou 0.0 e talvez com uma bandeira de erro definida para algo como TLOSS perda total de precisão ou PLOSS perda parcial de precisão.

Uma boa redução da Gama de valores grandes para um intervalo Como-π Para π é um problema desafiador que rivaliza com os desafios da função trig básica, como sin(), si.

Um bom relatório é argumento redução para argumentos enormes: bom até o último bit (1992). Ele cobre 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 dá resultados de qualidade para Tudo double de -DBL_MAX a DBL_MAX.

---

Se os argumentos originais estão em graus, mas podem ter um grande valor, use fmod() primeiro para uma maior precisão. Uma boa vontade introduza nenhum erro e por isso forneça uma excelente redução de gama.

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

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

Eles são tipicamente implementados em software e não vão usar o hardware correspondente (ou seja, chamadas de aseembly) na maioria dos casos. No entanto, como Jason apontou, estes são implementação específica.

Note que estas rotinas de software não fazem parte das fontes do compilador, mas sim serão encontradas na biblioteca correspondente, como o clib, ou glib para o compilador GNU. Ver http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions

Se se quiser um maior controlo, deve avaliar cuidadosamente o que precisa. Alguns dos métodos típicos são a interpolação de tabelas de consulta, a chamada de montagem (que é muitas vezes lenta), ou outros esquemas de aproximação como Newton-Raphson para raízes quadradas.

Se você quer uma implementação em software, não em hardware, o local para procurar uma resposta definitiva a esta pergunta é o Capítulo 5 de receitas numéricas. Minha cópia está em uma caixa, então eu não posso dar detalhes, mas a versão curta (Se eu me lembro deste direito) é que você toma tan(theta/2) como sua operação primitiva e computa os outros de lá. A computação é feita com uma aproximação de séries, mas é algo que converge Muito mais rapidamente que um Taylor série.

Desculpa não me lembrar de mais sem pôr a mão no livro.

Como muitas pessoas salientaram, depende da implementação. Mas tanto quanto eu entendo sua pergunta, você estava interessado em um realsoftware implementação de funções matemáticas, mas apenas não conseguiu encontrar um. Se for esse o caso, então aqui está você:

• Obter o código-fonte da glibc de http://ftp.gnu.org/gnu/glibc/
• Veja o ficheiro dosincos.c localizado na pasta glibc root \sysdeps\ieeee754\dbl-64
• da mesma forma você pode encontre implementações do resto da biblioteca de matemática, basta procurar o arquivo com o nome apropriado

Você também pode ter uma olhada nos arquivos com a extensão .tbl, seu conteúdo não é nada mais do que tabelas enormes de pré-compostos 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 usam eles apenas fazem uma pesquisa rápida, que é Muito mais rápido. BTW, eles usam Série de alfaiate para calcular seno e cosseno.

Espero que isto ajude.

Vou tentar responder pelo caso de {[[0]} num programa C, compilado com o compilador C do GCC num processador x86 actual (digamos um duo Intel Core 2).

Na linguagem C, A Biblioteca C padrão inclui funções matemáticas comuns, não incluídas na própria língua (ex. pow, sin e cos para potência, seno e cosseno, respectivamente). Cujos cabeçalhos estão incluídos em matemática.h.

Agora num sistema GNU / Linux, estas funções de bibliotecas são fornecidas por glibc (GNU libc ou GNU C Library). Mas o compilador GCC quer que você ligue para a biblioteca de matemática (libm.so) Usando a bandeira do compilador -lm para permitir a utilização destas funções matemáticas. Não sei porque não faz parte da biblioteca C. estas seriam uma versão de software das funções de ponto flutuante, ou "soft-float".

De lado: a razão para ter as funções matemáticas separadas é histórica, e foi meramente destinada a reduzir o tamanho do executável programas em muito sistemas Unix antigos, possivelmente antes de bibliotecas compartilhadas estarem disponíveis, tanto quanto sei.

Agora o compilador pode otimizar a função biblioteca C padrão {[[0]} (fornecida por libm.so) para ser substituída por uma chamada para uma instrução nativa para a função sin (built-in do seu CPU/FPU), que existe como uma instrução FPU (FSIN para x86/x87) em processadores mais recentes como a série Core 2 (isto é correto, tanto quanto o i486DX). Isso dependeria da otimização bandeiras passadas para o compilador gcc. Se o compilador foi dito para escrever código que seria executado em qualquer i386 ou processador mais recente, ele não faria tal otimização. A bandeira -mcpu=486 informaria o compilador de que era seguro fazer tal otimização.

Agora, se o programa executado a versão do software do pecado() função, seria fazê-lo com base em uma CORDIC (Coordenadas de Rotação Computador DIgital) ou BKM algoritmo, ou mais, provavelmente, de uma tabela ou cálculo de séries de potência que é comumente usado agora para calcular tais funções transcendentais. [Src: http://en.wikipedia.org/wiki/Cordic#Application]

Qualquer recente (desde 2, 9 x aprox.) versão do gcc também oferece uma versão embutida do sin, __builtin_sin() que será usado para substituir a chamada padrão para a versão da biblioteca C, como uma otimização.

Tenho a certeza que isso é claro como lama, mas espero que te dê mais informações do que esperavas, e montes de saltos. off points to learn more yourself.

Não há nada como atingir a fonte e ver como alguém realmente fez isso em uma biblioteca de uso comum; vamos olhar para uma implementação da biblioteca C em particular. Escolhi a uLibC.

Aqui está a função do pecado:

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

Que parece tratar de alguns casos especiais,e depois realiza alguma redução de argumentos para mapear a entrada para o intervalo [-pi/4, pi/4], (dividindo o argumento em duas partes, um grande parte e cauda) antes de chamar

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

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

Sempre que tal função é 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, gráficos informáticos)
• a avaliação de uma série que converge para o valor desejado --- provavelmente não uma série de taylor, mais provavelmente algo baseado numa quadratura chique como Clenshaw-Curtis.

Se não houver suporte de hardware, então o compilador provavelmente usa o último método, emitindo apenas código de montagem (sem símbolos de depuração), ao invés de usar uma biblioteca c --- tornando difícil para você rastrear o código real para baixo em seu depurador.

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

Calcular seno / cosseno / tangente é, na verdade, muito fácil de fazer através do código usando a série Taylor. Escrever um leva uns 5 segundos.

Todo o processo pode ser resumido com esta equação aqui: http://upload.wikimedia.org/math/5/4/6/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<b; i++)
        c *= a;
    return c;
}

double _fact(double x) {
    double ret = 1;
    for (int i=1; i<=x; i++) 
        ret *= i;
    return ret;
}

double _sin(double x) {
    double y = x;
    double s = -1;
    for (int i=3; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _cos(double x) {
    double y = 1;
    double s = -1;
    for (int i=2; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _tan(double x) {
     return (_sin(x)/_cos(x));  
}

Não uses a série Taylor. Os polinômios de Chebyshev são mais rápidos e precisos, como apontado por algumas pessoas acima. Aqui está uma implementação (originalmente da ROM ZX Spectrum): https://albertveli.wordpress.com/2015/01/10/zx-sine/

Se queres pecado, então asm Volátil ("fsin": "=t"(vsin): "0" (xrads))); se você quer cos então asm Volátil ("fcos": "=t"(vcos): "0" (xrads)); se você quer sqrt então asm Volátil ("fsqrt": "=t"(vsqrt): "0" (Valor)); então, por que usar código impreciso quando as instruções da máquina vai fazer.