Comment c calcule-t-il sin() et les autres fonctions mathématiques?

OK kiddies, c'est l'heure des pros.... C'est l'une de mes plus grandes plaintes inexpérimenté ingénieurs en logiciel. Ils viennent dans le calcul des fonctions transcendantales à partir de zéro (en utilisant la série de Taylor) comme si personne n'avait jamais fait ces calculs auparavant dans leurs vies. Pas vrai. Il s'agit d'un problème bien défini qui a été abordé des milliers de fois par des ingénieurs logiciels et matériels très intelligents et dont la solution est bien définie. Fondamentalement, la plupart des fonctions transcendantes utilisation Polynômes de tchebychev pour les calculer. Les polynômes utilisés dépendent des circonstances. D'abord, la bible sur cette question est un livre appelé "Computer Approximations" par Hart et Cheney. Dans ce livre, vous pouvez décider si vous avez un adder matériel, un multiplicateur, un diviseur, etc, et décider quelles opérations sont les plus rapides. par exemple, si vous aviez un diviseur vraiment rapide, la façon la plus rapide de calculer sine pourrait être P1(x)/P2(x) où P1, P2 sont des polynômes de Chebyshev. Sans le diviseur rapide, il pourrait être juste P (x), où P a beaucoup plus de termes que P1 ou P2....donc ce serait plus lent. Donc, la première étape est de déterminer votre matériel et ce qu'il peut faire. Ensuite, vous choisissez la combinaison appropriée de polynômes de Chebyshev (est généralement de la forme cos(ax) = aP(x) pour cosinus par exemple, encore une fois où P est un polynôme de Chebyshev). Ensuite, vous décidez quelle précision décimale vous voulez. par exemple, si vous voulez une précision de 7 chiffres, vous cherchez dans le tableau approprié dans le livre que j'ai mentionné, et il vous donnera (pour précision = 7,33) un nombre N = 4 et un nombre polynomial 3502. N est l'ordre du polynôme (donc c'est p4.x ^ 4 + p3.x^3 + p2.x^2 + p1.x + p0), parce que N = 4. Ensuite, vous regardez la valeur réelle des valeurs p4,p3,p2,p1,p0 au dos du livre sous 3502 (elles seront en virgule flottante). Ensuite, vous mettez en œuvre votre algorithme dans le logiciel sous la forme: (((p4.x + p3).x + p2).x + p1).x + p0 ....et c'est comme ça qu'on calcule cosine à 7 décimales sur ce matériel.

Note que la plupart des implémentations matérielles d'opérations transcendantales dans un FPU impliquent généralement un certain microcode et des opérations comme celle-ci (cela dépend du matériel). Les polynômes de Chebyshev sont utilisés pour la plupart des transcendentals mais pas tous. par exemple, la racine carrée est plus rapide pour utiliser une double itération de la méthode de Newton raphson en utilisant d'abord une table de recherche. Encore une fois, ce livre "Computer Approximations" vous le dira.

si vous prévoyez d'implémenter ces fonctions, je recommande à quiconque que recevoir une copie de ce livre. C'est vraiment la bible pour ces types d'algorithmes. Notez qu'il existe des groupes de moyens alternatifs pour calculer ces valeurs comme les cordics, etc, mais ceux-ci ont tendance à être meilleurs pour des algorithmes spécifiques où vous n'avez besoin que d'une faible précision. Pour garantir la précision à chaque fois, les polynômes de chebyshev sont la solution. Comme je l'ai dit, problème bien défini. A été résolu depuis 50 ans maintenant.....et c'est comment il est fait.

Maintenant, cela étant dit, il existe des techniques par lesquelles les polynômes de Chebyshev peuvent être utilisés pour obtenir un résultat de précision unique avec un polynôme de faible degré (comme l'exemple pour cosinus ci-dessus). Ensuite, il existe d'autres techniques pour interpoler entre les valeurs pour augmenter la précision sans avoir à aller à un polynôme beaucoup plus grand, comme "Gal's Accurate Tables Method". Cette dernière technique est celle à laquelle se réfère le post faisant référence à la littérature ACM. Mais en fin de compte, les polynômes de Chebyshev sont ce qui est utilisé pour obtenir 90% du chemin.

de Profiter de.

fonctions comme sine et cosine sont mis en œuvre dans le microcode à l'intérieur des microprocesseurs. Les puces Intel, par exemple, ont des instructions d'assemblage pour celles-ci. Un compilateur C générera du code qui appelle ces instructions d'assemblage. (Par contre, un compilateur Java ne le fera pas. Java évalue les fonctions trig dans le logiciel plutôt que le matériel, et donc il tourne beaucoup plus lentement.)

puces ne pas utiliser la série Taylor pour calculer les fonctions de trig, au moins pas entièrement. Tout d'abord ils utilisent CORDIC , mais ils peuvent également utiliser une série de Taylor courte pour polir le résultat de CORDIC ou pour des cas spéciaux tels que le calcul sinus avec une haute précision relative pour de très petits angles. Pour plus d'explications, voir cette StackOverflow answer .

Oui, il y a aussi des algorithmes logiciels pour calculer sin . Fondamentalement, le calcul de ce genre de choses avec un ordinateur numérique est habituellement fait en utilisant méthodes numériques comme approximant le série Taylor représentant la fonction.

méthodes numériques peuvent approximer des fonctions à un montant arbitraire de précision et puisque la quantité de précision que vous avez dans un nombre flottant est fini, ils conviennent ces tâches assez bien.

C'est une question complexe. Les processeurs Intel de la famille x86 ont une implémentation matérielle de la fonction sin() , mais ils font partie du FPU X87 et ne sont plus utilisés en mode 64 bits (où les registres SSE2 sont utilisés à la place). Dans ce mode, une implémentation logicielle est utilisée.

il y a plusieurs implémentations de ce genre là-bas. L'un est en fdlibm et est utilisé en Java. Pour autant que je sache, l'implémentation glibc contient des parties de fdlibm, et d'autres pièces fournies par IBM.

implémentations logicielles de fonctions transcendantales telles que sin() utilisent généralement des approximations par polynômes, souvent obtenues à partir de la série Taylor.

utiliser série taylor et essayer de trouver la relation entre les termes de la série afin que vous ne calculez pas les choses encore et encore

voici un exemple pour 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;}

en utilisant ceci nous pouvons obtenir le nouveau terme de la somme en utilisant celui déjà utilisé (nous évitons le factoriel et x^2p)

explication http://img514.imageshack.us/img514/1959/82098830.jpg

les polynômes de Chebyshev, comme mentionné dans une autre réponse, sont les polynômes où la plus grande différence entre la fonction et le polynôme est aussi petite que possible. C'est un excellent début.

Dans certains cas, l'erreur maximale est pas ce qui vous intéresse, mais le maximum de l'erreur relative. Par exemple, pour la fonction sine, l'erreur près de x = 0 devrait être beaucoup plus petite que pour les valeurs plus grandes; vous voulez un petit relatif erreur. Ainsi, on calcule le polynôme de Chebyshev pour sin x / x, et on multiplie ce polynôme par X.

ensuite, vous devez trouver comment évaluer le polynôme. Vous voulez l'évaluer de manière à ce que les valeurs intermédiaires soient petites et donc que les erreurs d'arrondissement soient petites. Sinon, les erreurs d'arrondi peuvent devenir beaucoup plus grande que les erreurs dans le polynôme. Et avec des fonctions comme la fonction sinus, si vous êtes négligent, alors il peut être possible que l' résultat que vous calculez pour sin x est plus grand que le résultat pour sin y même si x < Y. Il faut donc choisir avec soin l'ordre de calcul et le calcul des limites supérieures de l'erreur d'arrondissement.

par exemple, sin x = x - x^3/6 + x^5 / 120 - x^7 / 5040... Si vous calculez naïvement sin x = x * (1 - x^2/6 + x^4/120 - x^6/5040...), alors que la fonction entre parenthèses est en baisse, et il se produira que si y est le prochain plus grand nombre à x, alors parfois sin y sera plus petit que sin X. Calculez plutôt sin x = x-x^3 * (1/6-x^2 / 120 + x^4/5040...) où cela ne se fait pas.

pour calculer les polynômes de Chebyshev, vous devez généralement arrondir les coefficients à double précision, par exemple. Mais alors qu'un polynôme de Chebyshev est optimal, le polynôme de Chebyshev avec des coefficients arrondis à la double précision n'est pas le polynôme optimal avec des coefficients de double précision!

Par exemple pour sin (x), où vous avez besoin de coefficients pour x, x^3, x^5, x^7 etc. vous faites ce qui suit: Calculez la meilleure approximation de sin x avec un polynôme (ax + bx^3 + cx^5 + dx^7) avec une précision supérieure à la double, puis arrondissez a à la double précision, donnant A. la différence entre a et a serait assez grande. Calculer maintenant la meilleure approximation de (sin x - Ax) avec un polynôme (b x^3 + cx^5 + dx^7). Vous obtenez des coefficients différents, parce qu'ils s'adaptent à la différence entre a et A. Tour b à double précision B. puis approximatif (sin x-Ax-BX^3) avec un polynôme cx^5 + dx^7 et ainsi de suite. Vous obtiendrez un polynôme qui est presque aussi bon que le polynôme Chebyshev original, mais beaucoup mieux que Chebyshev arrondi à double précision.

ensuite, vous devez tenir compte des erreurs d'arrondissement dans le choix du polynôme. Vous avez trouvé un polynôme avec un minimum d'erreur dans le polynôme en ignorant l'erreur d'arrondi, mais vous voulez optimiser le polynôme plus erreur d'arrondi. Une fois que vous avez le polynôme de Tchebychev, vous pouvez calculer les limites de l'erreur d'arrondi. Dire que f (x) est une fonction P (x) est le polynôme, et E (x) est l'erreur d'arrondi. Vous ne voulez pas optimiser | f (x) - P (x) |, vous souhaitez optimiser | f (x) - P (x) +/- E (x) |. Vous obtiendrez un polynôme légèrement différent qui essaie de maintenir les erreurs polynomiales en bas là où l'erreur d'arrondi est grande, et détend les erreurs polynomiales un peu là où l'erreur d'arrondi est petite.

Tout cela, vous obtiendrez facilement les erreurs d'arrondi au plus 0.55 fois le dernier morceau,+, -,*, / ont des erreurs d'arrondi au plus de 0,50 fois le dernier bit.

pour sin spécifiquement, en utilisant L'expansion de Taylor vous donnerait:

sin (x) := x - x^3/3! + x^5/5! - x^7j / 7! + ... (1)

vous continueriez à ajouter des termes jusqu'à ce que la différence entre eux soit inférieure à un niveau de tolérance accepté ou juste pour un nombre fini d'étapes (plus rapide, mais moins précis). Un exemple serait quelque chose comme:

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;
}

Note: (1) fonctionne en raison de la aproximation sin(x)=x pour les petits angles. Pour des angles plus grands, vous devez calculer de plus en plus de termes pour obtenir des résultats acceptables. Vous pouvez utiliser un argument while et continuer pour une certaine précision:

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;

}

la mise en œuvre effective des fonctions de bibliothèque dépend du compilateur et/ou du fournisseur de bibliothèque spécifique. Que ce soit en matériel ou en logiciel, Qu'il s'agisse d'une extension Taylor ou non, etc. vont varier.

je me rends compte que ça ne sert à rien.

concernant la fonction trigonométrique telle que sin() , cos() , tan() il n'a pas été fait mention, après 5 ans, d'un aspect important des fonctions de trig de haute qualité: réduction de gamme .

Une première étape dans l'une de ces fonctions est de réduire l'angle en radians, à une distance de 2*π intervalle. Mais π est irrationnel donc des réductions simples comme x = remainder(x, 2*M_PI) introduire l'erreur comme M_PI , ou machine pi, est un rapprochement de π. Alors, comment faire x = remainder(x, 2*π) ?

premières bibliothèques ont utilisé la précision étendue ou la programmation artisanale pour donner des résultats de qualité, mais encore sur une gamme limitée de double . Quand une grande valeur a été demandée comme sin(pow(2,30)) , les résultats étaient dénués de sens ou 0.0 et peut-être avec un indicateur d'erreur fixé à quelque chose comme TLOSS perte totale de précision ou PLOSS perte partielle de précision.

bon intervalle la réduction de grandes valeurs à un intervalle comme-π à π est un problème difficile qui rivalise avec les défis de la fonction trig de base, comme sin() , lui-même.

un bon rapport est réduction D'Argument pour les arguments énormes: bon au dernier bit (1992). Il couvre bien la question: discute le besoin et comment les choses étaient sur diverses plates-formes (SPARC, PC, HP, 30+ autres) et fournit un algorithme de solution donne la qualité résultats pour tous "1519350920 double de -DBL_MAX à DBL_MAX .

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