Le plus rapide de la mise en œuvre du sinus, du cosinus et de la racine carrée en C++ (n'a pas besoin d'être exacte)

Je suis googler la question pour la dernière heure, mais il y a seulement des points de la Série de Taylor ou un exemple de code qui est soit trop lent ou ne compile pas du tout. Eh bien, la plupart des réponses que j'ai trouvé sur Google "Google, c'est déjà demandé", mais malheureusement ce n'est pas...

Je suis de profilage de mon jeu sur le bas de gamme Pentium 4 et découvert que ~85% du temps d'exécution est perdu sur le calcul de sinus, cosinus et de la racine carrée (à partir de la bibliothèque C++ standard dans Visual Studio), et cela semble être très CPU dépendant (sur mon I7 les mêmes fonctions a obtenu que 5% de temps d'exécution, et le jeu est waaaaaaaaaay plus rapide). Je ne peut pas optimiser ces trois fonctions, ni de calculer le sinus et le cosinus en un seul passage (il y interdépendantes), mais je n'ai pas besoin de trop de résultats précis pour ma simulation, afin que je puisse vivre avec plus rapide rapprochement.

Donc, la question: Quels sont le moyen le plus rapide pour calculer le sinus, le cosinus et la racine carrée de float en C++?

MODIFIER
Table de recherche sont de plus en plus douloureux comme résultant de Cache Miss, c'est bien plus coûteux moderne de la CPU de la Série de Taylor. Les Processeurs sont juste tellement vite ces jours-ci, et le cache n'est pas.

J'ai fait une erreur, je l'ai bien que j'ai besoin de calculer plusieurs factorielles pour les Séries de Taylor, et je vois maintenant, ils peuvent être mis en œuvre sous forme de constantes.

De sorte que la mise à jour question: est-il rapide d'optimisation pour la racine carrée ainsi?

EDIT2

Je suis en utilisant la racine carrée de calculer la distance, pas de normalisation - ne peut pas utiliser rapide inverse de la racine carrée de l'algorithme (comme indiqué dans le commentaire: http://en.wikipedia.org/wiki/Fast_inverse_square_root

EDIT3

Je ne peux pas fonctionner sur des carrés des distances, j'ai besoin de distance exacte pour les calculs

Il a été résolu il y a plusieurs années - utilisation précalculées tableau pour obtenir le sinus/cosinus numéros.
stackoverflow.com/questions/3688649/...
Pour l'inverse de la racine carrée (ce qui est fréquent, puisqu'il est impliqué dans le vecteur de normalisation), il est bien connu de formule (en.wikipedia.org/wiki/Fast_inverse_square_root), mais honnêtement, il est un peu dépassé, et probablement 1.0/sqrt(x) va permettre à certaines d'optimisation du compilateur.
jetez un oeil à ce pour le sinus et le cosinus: devmaster.net/forums/topic/4648-fast-and-accurate-sinecosine
Ce n'est pas le moyen le plus rapide depuis récemment, la CPU sont waaaay plus rapide maintenant, alors il y a des années et le cache n'est pas beaucoup plus rapide. Je veux dire que le Cache est beaucoup plus douloureux, puis en calculant les sinus à l'aide de Tylor Série
Quel est exactement le problème avec les Séries de Taylor? Il ressemble exactement à ce dont vous avez besoin. Il vous permet de calculer sin, cos, tan jusqu'à ce que la précision que vous voulez. Si vous avez de la difficulté à la mettre en œuvre, puis de publier que d'une question distincte. Sinon, les gens ont suggéré une table de recherche qui peut être très bon, mais les recherches peuvent être coûteux aussi. Heureusement, la mise en œuvre est si rapide que vous pouvez faire les tests.
Vous avez vraiment besoin d'expliquer davantage sur ce que vous êtes en utilisant les résultats. Il existe une variété de différentes approches à différents compromis, et en fonction de ce que vous faites, vous pouvez également être en mesure d'éviter entièrement les opérations, l'utilisation SIMD, ou amortir leurs coûts à travers un grand nombre de calculs. Il est impossible de vous conseiller sur microoptimisation sans connaître les détails du code spécifique à microoptimised.
Pour obtenir un numéro de look-up table) est plus lent que le calcul de la série de Taylor? Et tout le profilage des résultats pour le prouver?
Que diriez-dectecting le CPU et l'utilisation d'un natif de l'instruction sur un processeur moderne avec une table de recherche ou d'autres code optimisé sur des machines plus anciennes.
Je suis assez sûr que je l'ai vu la comparaison des articles, mais ne peuvent pas les trouver dès maintenant
Si vous utilisez le calcul de la distance à des fins de comparaison avec une autre distance, un peu de travail avec le carré de la distance et vous pouvez le faire sans sqrt entièrement.
Comme je l'ai poser dessus, pouvez-vous expliquer pourquoi vous êtes le calcul de la distance? par exemple, si vous êtes à faire des calculs afin de la comparer avec quelque chose, ne faites pas de la place de la racine carrée de la valeur à comparer avec la place. De même, si vous expliquez ce que vous faites avec sin/cos il y a peut-être des moyens pour éviter ces opérations.
J'ai besoin de distance exacte, que les au carré de la distance de truc ne peut pas être utilisé dans mon cas
L'approche habituelle pour la racine carrée est d'obtenir une estimation pour la racine carrée (ou à l'inverse de la racine carrée) à l'aide d'un natif du PROCESSEUR de l'instruction, d'une table de recherche, Carmack la fonction etc, généralement à 1 partie en 64, puis l'utiliser dans de Newton-Raphson jusqu'à ce que vous avez suffisamment de précision (IIRC 4 tours pour plein précision). C'est effectivement ce que la bibliothèque C++ aurez à faire, donc si vous avez besoin d'un résultat exact que les chances sont que l'application que vous utilisez actuellement est déjà optimale et la seule façon d'améliorer serait d'utiliser SIMD (ESS, etc), ou de modification de votre algorithme pour éviter le fonctionnement / de réduire sa fréquence.
Quand vous posez ce genre de question, vous devez spécifier des conditions beaucoup plus de précision. Avez-vous des informations sur la distribution des numéros pour lesquels vous aurez à calculer sin/cos/sqrt (dire qu'ils sont tous proches de 0)? Avez-vous des contraintes spécifiques sur la précision (par exemple, le péché(0) absolument être 0)? etc. Toute information supplémentaire donne un moyen d'améliorer la solution.
Des tables pour les fonctions trigonométriques sont rarement une bonne idée sur n'importe quel processeur de ces 15 dernières années. Qu'est-ce que votre spécification minimale? Si votre spécification minimale est un Pentium 4 pensez à activer le SSE2 de génération de code et l'utilisation de l'ESS sqrt intrinsèques. Le moteur d'exécution lien de l'ESS des versions optimisées de fonctions trigonométriques ainsi ce qui peut être assez d'un facteur que vous n'avez pas besoin de chercher plus loin.
Si la précision est d'aucun intérêt que ce soit, vous pouvez utiliser 0, rapprochant sine, que sa sortie sera dans l'intervalle [-1, 1]. Huhu -- désolé, blague -- juste la vitesse et la précision ont tendance à être liés ensemble, et c'est un peu ambigu quant à quel degré d'approximation est acceptable. Lut ne sont généralement pas si bénéfique ici. J'ai trouvé ces bits tripoter nombre magique des solutions telles que l'id est rapide rsqrt de toujours offrir quelques petits avantages, même si je suis dans le C++03 ère (legacy/plate-forme) et un peu en retard sur la dernière optimiseurs/libs standard. Il y a sin/cos variantes de ces derniers aussi bien.

InformationsquelleAutor PiotrK | 2013-09-06

c++math optimization trigonometry x86

16

Le moyen le plus rapide est à pré-calculer les valeurs d'une utilisation d'un tableau comme dans cet exemple:

Créer sine table de recherche en C++

MAIS si vous insistez sur le calcul au moment de l'exécution, vous pouvez utiliser le développement en série de Taylor de l'expansion de sinus ou cosinus...

Pour en savoir plus sur les séries de Taylor... http://en.wikipedia.org/wiki/Taylor_series

L'un des éléments-clés pour obtenir que cela fonctionne bien est pré-calculer les factorielles et les tronquer à un bon nombre de termes. Le factorielles grandir dans le dénominateur très rapidement, de sorte que vous n'avez pas besoin de transporter plus de quelques termes.

Aussi...Ne multipliez pas vos x^n depuis le début à chaque fois...par exemple, multiplier x^3 par x un autre deux fois, alors que par deux pour calculer les exposants.
- Et comme pour la racine carrée?
- Il l'a dit au début de son post qu'il n'aime pas les séries de Taylor, mais n'a pas expliqué pourquoi.
- Mon erreur, je pensais que je dois calculer la Factorielle plusieurs fois, mais j'ai raté que je peux utiliser précalculées constante
- +1 pour precomputing factorielles, qui devrait enregistrer un tas de cycles de cpu par lui-même.
- Il y a un lien intéressant à la question des commentaires qui montre une méthode encore mieux que la série de Taylor.
- Ce n'est pas la meilleure façon de calculer sin/cos en termes d'efficacité. Il y a de vieux réponses dans stackoverflow qui avait déjà discuté de cela dans les moindres détails. Aussi - GSL bibliothèque scientifique GNU, qui est la norme numérique de la bibliothèque utilisée partout, aussi ne les utilise pas. Sachant que la meilleure procédure numérique vous permet de balance de précision/vitesse de façon plus précise.
- La série de Taylor avec 3 ou 4 itérations est sensiblement plus rapide qu'un lookuptable aujourd'hui, bien que ce n'est pas la manière la plus rapide possible (ce serait une approximation parabolique, je pense que le code est sur Flipcode 10 ans ou si il ya...).
- rpc.sh/7rksf mais cela fonctionne plus lentement que cmath sine
InformationsquelleAutor KeatsKelleher
39

Tout d'abord, le développement en série de Taylor n'est PAS le meilleur/le plus rapide à mettre en œuvre sinus/cosinus. Il est également pas la façon dont les professionnels des bibliothèques de mise en œuvre de ces fonctions trigonométriques, et de connaître le meilleur numérique de la mise en œuvre permet d'affiner la précision pour obtenir la vitesse de manière plus efficace. En outre, ce problème a déjà été largement discuté dans StackOverflow. Voici juste un exemple.

Deuxième, la grande différence que vous voyez entre les anciens/nouveaux PC est dû au fait que la modernité de l'architecture Intel a explicite de l'assemblée du code pour calculer elementar fonctions trigonométriques. Il est assez difficile de les battre sur la vitesse d'exécution.

Enfin, parlons de ce code sur votre ancien PC. Vérifier gsl gnu scientific library (ou recettes ou numérique) de la mise en œuvre, et vous verrez qu'ils sont fondamentalement utilisation de Tchebychev Approximation de la Formule.

Approximation de Chebyshev converge plus rapidement, de sorte que vous devez évaluer les moins de termes. Je ne vais pas écrire de mise en œuvre en détail ici, car il y a déjà de très belles réponses postées sur StackOverflow. Vérifier celui-ci par exemple . Juste régler le nombre de termes de cette série de modifier l'équilibre entre la précision et de vitesse.

Par la route: de la règle à 0 pour ce genre de problème: si vous voulez des détails de mise en œuvre de certaines de fonction ou de la méthode numérique, vous devriez jeter un oeil sur la GSL code avant toute action - GSL est LA NORME numérique de la bibliothèque.

EDIT: vous pouvez améliorer le temps d'exécution en incluant agressif à virgule flottante options d'optimisation de gcc/icc. Cela permettra de diminuer la précision, mais il semble que c'est exactement ce que vous voulez.

EDIT2: Vous pouvez essayer de faire un gros péché de la grille et utilisez gsl routine (gsl_interp_cspline_periodic pour spline avec des conditions) pour spline de la table (la spline permettra de réduire les erreurs par rapport à une interpolation linéaire => vous avez besoin de moins de points sur votre table => moins de cache miss)!

InformationsquelleAutor Vivian Miranda
25

Ici est la garantie la plus rapide possible de la fonction sinus en C++:
```
double FastSin(double x)
{
    return 0;
}
```
Oh, tu voulais une précision meilleure que |1.0|? Bien lecture sur.

Ingénieurs dans les années 1970, de formidables découvertes dans ce domaine, mais les nouveaux programmeurs sont tout simplement pas au courant que ces méthodes existent, parce qu'ils ne sont pas instruits dans le cadre de la norme en sciences informatiques programmes d'études.

Vous devez commencer par comprendre que il n'est pas "parfait" de la mise en œuvre de ces fonctions pour toutes les applications. Par conséquent, superficielle réponses à des questions comme "qui est le plus rapide" sont garantis d'être mauvais.

La plupart des gens qui posent cette question ne comprennent pas l'importance de la compromis entre performances et précision. En particulier, vous allez avoir à faire des choix quant à l'exactitude des calculs avant de faire quoi que ce soit d'autre. Combien d'erreur pouvez-vous tolérer dans le résultat? 10^-4? 10^-16?

À moins que vous pouvez quantifier l'erreur dans une méthode, ne l'utilisez pas.

Aucune personne n'utilise les séries de Taylor seul à rapprocher des êtres transcendantaux dans le logiciel. À l'exception de certains cas bien spécifiques, la série de Taylor en général l'approche de la cible lentement à travers le commun des gammes d'entrée.

Les algorithmes que vos grands-parents utilisée pour calculer les êtres transcendantaux efficacement, sont désignés collectivement comme CORDIC et ont été assez simple pour être mis en œuvre dans le matériel. Voici une bien documenté CORDIC mise en œuvre en C. CORDIC implémentations, en général, nécessite une très petite table de recherche, mais la plupart des implémentations ne nécessitent même pas de matériel multiplicateur être disponible. La plupart des CORDIC implémentations de vous permettre de traiter des performances de précision, y compris celui que j'ai lié.

Il y a eu beaucoup d'améliorations incrémentales à l'origine CORDIC algorithmes au fil des ans. Par exemple, l'année dernière, des chercheurs au Japon ont publié un l'article sur une amélioration de la CORDIC avec de meilleurs angles de rotation, ce qui réduit les opérations nécessaires.

Si vous avez du matériel multiplicateurs assis autour de (et vous en avez certainement), ou si vous ne pouvez pas se permettre un tableau de recherche comme CORDIC l'exige, vous pouvez toujours utiliser un Polynôme de tchebychev de faire la même chose. Polynômes de tchebychev besoin multiplie, mais c'est rarement un problème sur le matériel moderne. Nous aimons les polynômes de Tchebychev parce que ils ont hautement prévisible erreurs maximales pour un rapprochement. Le maximum de le dernier terme d'un polynôme de Tchebychev, au sein de votre gamme d'entrée, les limites de l'erreur dans le résultat. Et cette erreur est plus petit que le nombre de termes est importante, plus. Voici un exemple d'un polynôme de Tchebychev de donner une condition sine rapprochement à travers une vaste gamme, en ignorant la symétrie naturelle de la fonction sinus et seulement à la résolution des le rapprochement problème en mettant plus de coefficients à elle.

Nous aussi comme les polynômes de Tchebychev, car l'erreur dans l'approximation est également répartie entre les différents résultats. Si vous êtes à l'écriture de plugins audio ou de faire un traitement du signal numérique, polynômes de Tchebychev vous donner un bon marché et prévisible de tramage effet "gratuitement".

Si vous voulez trouver vos propres coefficients du polynôme de Tchebychev à travers une gamme spécifique, de nombreuses bibliothèques de mathématiques appeler le processus de trouver les coefficients "De tchebychev ajustement" ou quelque chose comme ça.

Racines carrées, alors comme aujourd'hui, sont généralement calculées avec une variante de la Algorithme de Newton-Raphson, généralement avec un nombre fixe d'itérations. Habituellement, quand quelqu'un développe un "inédite" algorithme pour faire des racines carrées, il s'agit simplement de Newton-Raphson dans le déguisement.

De Newton-Raphson, CORDIC, et les polynômes de Tchebychev de vous laisser sur le compromis vitesse-précision, de sorte que la réponse peut être tout aussi imprécis que vous le souhaitez.

Enfin, lorsque vous avez terminé tous votre fantaisie d'analyse comparative et de la micro-optimisation, assurez-vous que votre "rapide" est effectivement plus rapide que la version de bibliothèque. Voici un exemple typique de la bibliothèque de la mise en œuvre de la fnis() délimitée sur le domaine de -pi/4, pi/4. Et il n'est tout simplement pas ce maudit lent.

Il y a des gens qui ont consacré leur vie à la résolution de ces problèmes de façon efficace, et ils ont produit quelques résultats fascinants. Lorsque vous êtes prêt à rejoindre la vieille école, ramasser une copie de Numérique Recettes.

tl:dr; go google "sine rapprochement" ou "cosinus rapprochement" ou "racine carrée rapprochement" ou "théorie de l'approximation."
- Pour float/double, la plupart des plates-formes efficaces de matériel sqrt. Sur x86, matériel sqrt est plus rapide que tout ce que vous pourriez faire cuire vous-même, à l'exception de l'utilisation du matériel rapide approximative réciproque sqrt instruction. Je suppose que sans un matériel FPU, ou si elle est très lente sqrt mais rapidement se multiplier, NR pourrait être une victoire.
- Le matériel x86 lui-même est en train de faire un Newton-Raphson itératif rapprochement.
- Il n'a pas d'importance la façon dont le matériel est câblé; tout ce qui importe est la façon dont il est rapide par rapport à un FP multiplier (ou fused multiply-add). Le sqrt instruction est une boîte noire qui crache correctement arrondies sqrt résultats très rapide (par exemple, sur Skylake avec 12 cycle de latence, un par 3 cycle de traitement). Vous ne pouvez pas battre ce avec un Newton-Raphson itération de départ avec rsqrtps (approx réciproque sqrt). À l'aide de tout rsqrtps (ce qui donne 12 bits de précision) est plus rapide, ou si vous avez besoin de la racine au lieu de l'inverse, x * approx_rsqrt(x) est un peu plus rapide que sqrt(x).
- À moins que vous goulot d'étranglement sur uop débit plutôt que sqrt temps de latence, dans le cas où l'utilisation de la plaine sqrtps est plus rapide encore que rsqrtps + fmaddps, parce qu'il sqrtps décode à un seul uop (la table de la recherche + de Newton-Raphson qui se passe à l'intérieur de la cloison de l'unité, n'est pas conduit par le microcode qui serait en concurrence avec d'autres instructions pour l'exécution des ressources).
- Vous êtes un x86 expert et votre avis est correct de ce point de vue. Mon code a exécuter sur toutes les plateformes imaginables, alors j'ai tendance à penser de manière algorithmique avant de tirer l'assembleur. "Le meilleur optimiseur est entre vos oreilles" -Abrash
- C'est certainement un luxe pour vectoriser manuellement pour une seule ISA où tout le MATÉRIEL qui vous est efficace sqrt. Je ne sais pas ce que HW sqrt est comme sur les BRAS, haut de gamme ou PowerPC (utilisé dans quelques grappes HPC). Il est certainement important d'avoir de bons algorithmes, mais lorsque la constante de facteurs de question, il est important de comprendre les coûts relatifs et de savoir quand un rapprochement est en fait plus lentement. sqrt est spécial parce que c'est l'une des fonctions que la norme IEEE nécessite de produire un correctement, le résultat arrondi, avec + - * /, de sorte qu'il est généralement mis en œuvre dans le matériel.
- Si le BRAS en question a NEON capacités, alors il peut faire sqrt. Du pouce et de l'Thumb2 ont pas intégré dans sqrt. Aucune idée sur PPC.
- Ok, mais la question est de savoir comment rapide que le matériel sqrt est. P5 (pentium 586) a x87 fsqrt, mais il en faut 70 cycles d'horloge, vs fadd / fmul de 3 cycle de latence. (P5 est une commande de la CPU, mais il peut se chevauchent fsqrt avec entier instructions pour 69 de ces cycles d'horloge, mais seulement par 2 de 70 cycles avec d'autres FP travail). Un rapide approximative sqrt serait certainement utile d'examiner sur P5.
- Pour "assurez-vous que votre "rapide" est effectivement plus rapide que la version de bibliothèque", +1
- Si vous déclarez FastSine comme constexpr, il sera encore plus rapide.
- davidhigh, vous allez référence pour nous.
InformationsquelleAutor johnwbyrd
19

Pour la racine carrée, il y a une approche appelée décalage de bits.

Un nombre à virgule défini par la norme IEEE-754 est à l'aide de certains bits représentent décrire fois de multiples basé sur 2. Certains bits sont pour représenter la valeur de base.
```
float squareRoot(float x)
{
  unsigned int i = *(unsigned int*) &x;

  //adjust bias
  i  += 127 << 23;
  //approximation of square root
  i >>= 1;

  return *(float*) &i;
}
```
C'est une constante de temps de calcul de la carrés de la racine
- Quel est le but de prendre l'adresse de i puis moulage float* avant d'être déréférencé il de nouveau?
- en.wikipedia.org/wiki/...
- type de beaucoup les jeux de mots. Sachant que la norme ieee représentation de float, cela peut fonctionner.
- Exprimés par la référence est plus lisible pour ce genre de hack.
InformationsquelleAutor BigTailWolf
10

Basé sur l'idée de http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648 et certains manuels de réécriture pour améliorer la performance dans un micro de référence j'ai fini avec la suite de cosinus de mise en œuvre qui est utilisé dans le HPC de simulation physique qui est un goulot d'étranglement par la répétition de cos appels sur un grand nombre de l'espace. C'est assez précis et beaucoup plus vite qu'une table de recherche, et plus particulièrement la division est nécessaire.
```
template<typename T>
inline T cos(T x) noexcept
{
    constexpr T tp = 1./(2.*M_PI);
    x *= tp;
    x -= T(.25) + std::floor(x + T(.25));
    x *= T(16.) * (std::abs(x) - T(.5));
    #if EXTRA_PRECISION
    x += T(.225) * x * (std::abs(x) - T(1.));
    #endif
    return x;
}
```
Le compilateur Intel, au moins, est également assez intelligent dans la vectorisation de cette fonction dans une boucle.

Si EXTRA_PRECISION est défini, l'erreur maximale est d'environ 0.00109 de la plage -à π π, en supposant T est double comme il est généralement défini dans la plupart des C++ implémentations. Sinon, l'erreur maximale est d'environ de 0,056 pour la même gamme.
- il y a une division dans la première ligne
- Oui, mais c'est une constante de compilation de la division qui est infiniment bon marché au moment de l'exécution 😛
- J'aimerais voir une référence à la bibliothèque standard de cosinus. stackoverflow.com/questions/824118/why-is-floor-so-slow
- Je n'ai pas accès à une licence de le compilateur Intel, où la différence était la plus grande. La raison d'être comme je l'ai dit qu'il a réussi à vectoriser la fonction ci-dessus et en l'entourant de code, alors qu'il ne l'ai pas fait en tant que bien pour le std::cos. Le lien pour "plancher est lent" montre également comment -ffast-math aide à soulager la question quelque peu. La CPI fait par défaut.
- mesurer par vous-même: godbolt.org/g/YhbQTk
InformationsquelleAutor milianw
5

Pour les architectures x86, le matériel FP racine carrée instructions sont rapides (sqrtss est sqrt Scalaire Simple précision). Simple précision est plus rapide que la double précision, donc certainement utiliser float au lieu de double pour le code où vous pouvez vous permettre d'utiliser moins de précision.

Pour les 32 bits de code, vous avez généralement besoin de compilateur options pour l'amener à faire FP mathématiques avec des instructions SSE, plutôt que de x87. (par exemple,-mfpmath=sse)

Pour obtenir C sqrt() ou sqrtf() fonctions inline comme juste sqrtsd ou sqrtss, vous avez besoin de compiler avec -fno-math-errno. Avoir les fonctions mathématiques ensemble errno sur NaN est généralement considérée comme une erreur de conception, mais la norme l'exige. Sans cette option, gcc inlines, mais alors il ne le compare+les branches pour voir si le résultat est NaN, et si oui, appelle une fonction de la bibliothèque alors il peut errno. Si votre programme ne vérifie pas errno après les fonctions mathématiques, il n'y a pas de danger à l'aide de -fno-math-errno.

Vous n'avez pas besoin de "dangereux" pièces de -ffast-math pour obtenir sqrt et quelques autres fonctions inline mieux ou pas du tout, mais -ffast-math peut faire une grande différence (en permettant par exemple le compilateur à l'auto-vectorisation dans les cas où cela change le résultat, parce que FP mathématiques n'est pas associatif.

par exemple avec gcc6.3 compilation float foo(float a){ return sqrtf(a); }
```
foo:    # with -O3 -fno-math-errno.
    sqrtss  xmm0, xmm0
    ret
```
```
foo:   # with just -O3
    pxor    xmm2, xmm2   # clang just checks for NaN, instead of comparing against zero.
    sqrtss  xmm1, xmm0
    ucomiss xmm2, xmm0
    ja      .L8          # take the slow path if 0.0 > a
    movaps  xmm0, xmm1
    ret

.L8:                     # errno-setting path
    sub     rsp, 24
    movss   DWORD PTR [rsp+12], xmm1   # store the sqrtss result because the x86-64 SysV ABI has no call-preserved xmm regs.
    call    sqrtf                      # call sqrtf just to set errno
    movss   xmm1, DWORD PTR [rsp+12]
    add     rsp, 24
    movaps  xmm0, xmm1    # extra mov because gcc reloaded into the wrong register.
    ret
```
gcc du code pour l'NaN cas semble bien trop compliqué; il n'a même pas utiliser le sqrtf valeur de retour! De toute façon, c'est le genre de dégâts que vous obtenez en fait sans -fno-math-errno, pour chaque sqrtf() site d'appel dans votre programme. Surtout que c'est juste du code de ballonnements, et aucun des .L8 bloc sera exécuté lors de la prise de la racine carrée d'un nombre >= 0.0, mais il y a encore plusieurs instructions supplémentaires dans la voie rapide.

Si vous savez que votre entrée à sqrt est non-zéro, vous pouvez utiliser le rapide mais très approximative réciproque sqrt instruction, rsqrtps (ou rsqrtss pour la version scalaire). Un Newton-Raphson itération le met à peu près la même précision que le matériel de simple précision sqrt instruction, mais pas tout à fait.

sqrt(x) = x * 1/sqrt(x), pour x!=0, de sorte que vous pouvez calculer une racine carrée avec rsqrt et de se multiplier. Ce sont à la fois rapides, même sur P4 (était-ce toujours d'actualité en 2013)?

Sur P4, il peut être utile d'utiliser rsqrt + itération de Newton pour remplacer une seule sqrt, même si vous n'avez pas besoin de diviser quelque chose en elle.

Voir aussi une réponse que j'ai écrit récemment au sujet de la manipulation des zéros lors du calcul de sqrt(x) x*rsqrt(x), avec une Itération de Newton. J'ai inclus une discussion d'erreur d'arrondi si vous voulez convertir la valeur FP pour un entier, et des liens vers d'autres questions pertinentes.

P4:
- sqrtss: 23c temps de latence, pas de pipeline
- sqrtsd: 38c temps de latence, pas de pipeline
- fsqrt (x87): 43 ° c temps de latence, pas de pipeline
- rsqrtss /mulss: 4c + 6c de latence. Peut-être l'un par 3c débit, car ils n'ont apparemment pas besoin de la même unité d'exécution (mmx vs fp).
- SIMD paniers versions sont un peu plus lentes
Skylake:
- sqrtss/sqrtps: 12c de latence, un par 3c débit
- sqrtsd/sqrtpd: 15-16c de latence, un par 4-6c débit
- fsqrt (x87): 14-21cc de latence, un par 4-7c débit
- rsqrtss /mulss: 4c + 4c temps de latence. Un par 1c débit.
- SIMD 128b vecteur versions sont à la même vitesse. 256b vecteur versions sont un peu plus élevé de temps de latence, près de la moitié du débit. Le rsqrtss version a plein rendement pour 256b vecteurs.
Avec une Itération de Newton, la rsqrt version n'est pas beaucoup, si en plus vite.

Numéros de Agner de Fog, d'expérimentation, de. Voir son microarch guides de comprendre ce qui rend le code d'une course rapide ou lente. Également voir les liens à la x86 la balise wiki.

IDK la meilleure façon de calculer sin/cos. J'ai lu que le matériel fsin /fcos (et seulement légèrement plus lente fsincos qui fait les deux à la fois) ne sont pas le moyen le plus rapide, mais IDK ce qui est.

InformationsquelleAutor Peter Cordes

QT est rapide implémentations de sinus (qFastSin) et cosinus (qFastCos) qui utilise table avec interpolation. Je l'utilise dans mon code et ils sont plus rapides que std:sin/cos et suffisamment précis pour que ce dont j'ai besoin (erreur ~= 0.01% je dirais):

https://code.woboq.org/qt5/qtbase/src/corelib/kernel/qmath.h.html#_Z8qFastSind

#define QT_SINE_TABLE_SIZE 256


inline qreal qFastSin(qreal x)
{
   int si = int(x * (0.5 * QT_SINE_TABLE_SIZE / M_PI)); //Would be more accurate with qRound, but slower.
   qreal d = x - si * (2.0 * M_PI / QT_SINE_TABLE_SIZE);
   int ci = si + QT_SINE_TABLE_SIZE / 4;
   si &= QT_SINE_TABLE_SIZE - 1;
   ci &= QT_SINE_TABLE_SIZE - 1;
   return qt_sine_table[si] + (qt_sine_table[ci] - 0.5 * qt_sine_table[si] * d) * d;
}

inline qreal qFastCos(qreal x)
{
   int ci = int(x * (0.5 * QT_SINE_TABLE_SIZE / M_PI)); //Would be more accurate with qRound, but slower.
   qreal d = x - ci * (2.0 * M_PI / QT_SINE_TABLE_SIZE);
   int si = ci + QT_SINE_TABLE_SIZE / 4;
   si &= QT_SINE_TABLE_SIZE - 1;
   ci &= QT_SINE_TABLE_SIZE - 1;
   return qt_sine_table[si] - (qt_sine_table[ci] + 0.5 * qt_sine_table[si] * d) * d;
}

La LUT et la licence peut être trouvé ici: https://code.woboq.org/qt5/qtbase/src/corelib/kernel/qmath.cpp.html#qt_sine_table

InformationsquelleAutor Adriel Jr

J'utilise le code suivant pour calculer les fonctions trigonométriques en quadruple précision. La constante N détermine le nombre de bits de précision requis (par exemple N=26 donnera simple précision). Selon le niveau de précision désiré, le précalculées de stockage peut être petit et tient dans le cache. Il nécessite seulement l'addition et de la multiplication des opérations et est également très facile à vectoriser.

L'algorithme de pré-calcule sin et cos valeurs de 0.5^i, i=1,...,N. on peut Alors combiner ces valeurs précalculées, pour calculer sin et cos pour n'importe quel angle jusqu'à une résolution de 0,5^N

template <class QuadReal_t>
QuadReal_t sin(const QuadReal_t a){
const int N=128;
static std::vector<QuadReal_t> theta;
static std::vector<QuadReal_t> sinval;
static std::vector<QuadReal_t> cosval;
if(theta.size()==0){
#pragma omp critical (QUAD_SIN)
if(theta.size()==0){
theta.resize(N);
sinval.resize(N);
cosval.resize(N);
QuadReal_t t=1.0;
for(int i=0;i<N;i++){
theta[i]=t;
t=t*0.5;
}
sinval[N-1]=theta[N-1];
cosval[N-1]=1.0-sinval[N-1]*sinval[N-1]/2;
for(int i=N-2;i>=0;i--){
sinval[i]=2.0*sinval[i+1]*cosval[i+1];
cosval[i]=sqrt(1.0-sinval[i]*sinval[i]);
}
}
}
QuadReal_t t=(a<0.0?-a:a);
QuadReal_t sval=0.0;
QuadReal_t cval=1.0;
for(int i=0;i<N;i++){
while(theta[i]<=t){
QuadReal_t sval_=sval*cosval[i]+cval*sinval[i];
QuadReal_t cval_=cval*cosval[i]-sval*sinval[i];
sval=sval_;
cval=cval_;
t=t-theta[i];
}
}
return (a<0.0?-sval:sval);
}

InformationsquelleAutor Dhairya

C'est une implémentation de la Série de Taylor de la méthode précédemment donné en akellehe réponse.

unsigned int Math::SIN_LOOP = 15;
unsigned int Math::COS_LOOP = 15;
//sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...
template <class T>
T Math::sin(T x)
{
T Sum       = 0;
T Power     = x;
T Sign      = 1;
const T x2  = x * x;
T Fact      = 1.0;
for (unsigned int i=1; i<SIN_LOOP; i+=2)
{
Sum     += Sign * Power / Fact;
Power   *= x2;
Fact    *= (i + 1) * (i + 2);
Sign    *= -1.0;
}
return Sum;
}
//cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...
template <class T>
T Math::cos(T x)
{
T Sum       = x;
T Power     = x;
T Sign      = 1.0;
const T x2  = x * x;
T Fact      = 1.0;
for (unsigned int i=3; i<COS_LOOP; i+=2)
{
Power   *= x2;
Fact    *= i * (i - 1);
Sign    *= -1.0;
Sum     += Sign * Power / Fact;
}
return Sum;
}

InformationsquelleAutor hkBattousai

1

Sur 100000000 test, milianw réponse est 2 fois plus lent que std::cos mise en œuvre.
Cependant, vous pouvez gérer l'exécuter plus rapidement en faisant les étapes suivantes:

->utiliser float

->ne pas utiliser de l'étage, mais static_cast

->n'utilisez pas d'abs mais ternaire conditionnelle

->utiliser des #define constante pour la division

->utilisation de la macro pour éviter l'appel de fonction
```
//1 /(2 * PI)
#define FPII 0.159154943091895
//PI /2
#define PI2 1.570796326794896619
#define _cos(x)         x *= FPII;\
x -= .25f + static_cast<int>(x + .25f) - 1;\
x *= 16.f * ((x >= 0 ? x : -x) - .5f);
#define _sin(x)         x -= PI2; _cos(x);
```
Sur 100000000 appel à std::cos et _ des _cos(x), std::cos exécuter sur ~14s vs ~3s pour _ des _cos(x) (un peu plus pour _sin(x))
- tous vos commentaires me font me demander si vous avez réellement compilé avec le compilateur les optimisations activées. Plus particulièrement, cette fonction doit obtenir inline, donc à l'aide d'une macro ou pas ne devrait faire aucune différence que ce soit.
- "L'appel d'une fonction en ligne peut ou ne peut pas générer un appel de fonction, qui, en général, subit une très petite quantité de surcharge. L'exacte des situations dans lesquelles une fonction en ligne est insérée varier selon le compilateur; plus de faire un effort de bonne foi pour inline petites fonctions (au moins lorsque l'optimisation est activée), mais il n'est pas nécessaire qu'ils le fassent (C99, §6.7.4):" (stackoverflow.com/questions/5226803/...)
- À droite, c'est le compilateur dépend, quel compilateur utilisez-vous? L'assemblée, ce n'est pas tout appel de fonction sur clang ou gcc: godbolt.org/g/UjAKBh j'irais jusqu'à affirmer que vous pouvez signaler un bug de votre compilateur, si ce n'est pas l'in-lining cette fonction. Similaires le compilateur doit faire l'constante devision pour vous, pas besoin d'obfusquer le code de définir la constante...
- Le Code a été compilé avec visual studio 2015. J'ai fait le test de 1,5 ans, et je ne me souviens pas si c'était avec l'optimisation activé ou pas (je suis d'accord il faut que j'écrive ce moment de répondre). Cependant, toutes les conditions dont j'ai été conservés parce qu'ils ont moins de temps d'exécution. En plus, il me semble que les personnes handicapées de l'optimisation est plus restrictive de la condition de permis d'optimisation, donc, cela signifie que le code précédent fonctionne dans plus de cas (si l'optimisation est désactivé).
- En ce qui concerne microsoft comportement avec fonction inline, miscrosoft dit: "le compilateur ne pas incorporer une fonction si son adresse est prise ou si elle est trop grande pour la ligne." (msdn.microsoft.com/en-us/library/cx3b23a3.aspx) qui est vraiment pas clair. Nous pouvons aussi avoir un oeil pour le mot-clé __forceinline
InformationsquelleAutor Hugo Zevetel

Partage mon code, c'est un 6ème polynomiale de degré, rien de spécial, mais disposés de façon à éviter pows. Sur Core i7 c'est 2,3 fois plus lente que celle de la mise en œuvre, bien qu'un peu plus vite pour [0..2*PI] gamme. Pour un vieux processeur, cela pourrait être une alternative à la norme sin/cos.

/*
On [-1000..+1000] range with 0.001 step average error is: +/- 0.000011, max error: +/- 0.000060
On [-100..+100] range with 0.001 step average error is:   +/- 0.000009, max error: +/- 0.000034
On [-10..+10] range with 0.001 step average error is:     +/- 0.000009, max error: +/- 0.000030
Error distribution ensures there's no discontinuity.
*/
const double PI          = 3.141592653589793;
const double HALF_PI     = 1.570796326794897;
const double DOUBLE_PI   = 6.283185307179586;
const double SIN_CURVE_A = 0.0415896;
const double SIN_CURVE_B = 0.00129810625032;
double cos1(double x) {
if (x < 0) {
int q = -x / DOUBLE_PI;
q += 1;
double y = q * DOUBLE_PI;
x = -(x - y);
}
if (x >= DOUBLE_PI) {
int q = x / DOUBLE_PI;
double y = q * DOUBLE_PI;
x = x - y;
}
int s = 1;
if (x >= PI) {
s = -1;
x -= PI;
}
if (x > HALF_PI) {
x = PI - x;
s = -s;
}
double z = x * x;
double r = z * (z * (SIN_CURVE_A - SIN_CURVE_B * z) - 0.5) + 1.0;
if (r > 1.0) r = r - 2.0;
if (s > 0) return r;
else return -r;
}
double sin1(double x) {
return cos1(x - HALF_PI);
}

InformationsquelleAutor Josh

Alors permettez-moi de reformuler ce que, cette idée vient de rapprocher les cosinus & sinus fonctions sur un intervalle [-pi/4,+pi/4], délimitée à l'erreur à l'aide de l'algorithme de Remez. Ensuite, à l'aide de la gamme réduite flotteur reste et une LUT pour les sorties cos & sinus de l'quotient entier, le rapprochement peut être déplacé à n'importe quel angulaire argument.

Son tout à fait unique et j'ai pensé qu'il pourrait être élargi de manière à faire un algorithme plus efficace en termes de bornée erreur.

void sincos_fast(float x, float *pS, float *pC){
float cosOff4LUT[] = { 0x1.000000p+00,  0x1.6A09E6p-01,  0x0.000000p+00, -0x1.6A09E6p-01, -0x1.000000p+00, -0x1.6A09E6p-01,  0x0.000000p+00,  0x1.6A09E6p-01 };
int     m, ms, mc;
float   xI, xR, xR2;
float   c, s, cy, sy;
//Cody & Waite's range reduction Algorithm, [-pi/4, pi/4]
xI  = floorf(x * 0x1.45F306p+00 + 0.5);
xR  = (x - xI * 0x1.920000p-01) - xI*0x1.FB5444p-13;
m   = (int) xI;
xR2 = xR*xR;
//Find cosine & sine index for angle offsets indices
mc = (  m  ) & 0x7;     //two's complement permits upper modulus for negative numbers =P
ms = (m + 6) & 0x7;     //two's complement permits upper modulus for negative numbers =P, note phase correction for sine.
//Find cosine & sine
cy = cosOff4LUT[mc];     //Load angle offset neighborhood cosine value 
sy = cosOff4LUT[ms];     //Load angle offset neighborhood sine value 
c = 0xf.ff79fp-4 + xR2 * (-0x7.e58e9p-4);               //TOL = 1.2786e-4
//c = 0xf.ffffdp-4 + xR2 * (-0x7.ffebep-4 + xR2 * 0xa.956a9p-8);  //TOL = 1.7882e-7
s = xR * (0xf.ffbf7p-4 + x2 * (-0x2.a41d0cp-4));   //TOL = 4.835251e-6
//s = xR * (0xf.fffffp-4 + xR2 * (-0x2.aaa65cp-4 + xR2 * 0x2.1ea25p-8));  //TOL = 1.1841e-8
*pC = c*cy - s*sy;     
*pS = c*sy + s*cy;
}
float sqrt_fast(float x){
union {float f; int i; } X, Y;
float ScOff;
uint8_t e;
X.f = x;
e = (X.i >> 23);           //f.SFPbits.e;
if(x <= 0) return(0.0f);
ScOff = ((e & 1) != 0) ? 1.0f : 0x1.6a09e6p0;  //NOTE: If exp=EVEN, b/c (exp-127) a (EVEN - ODD) := ODD; but a (ODD - ODD) := EVEN!!
e = ((e + 127) >> 1);                            //NOTE: If exp=ODD,  b/c (exp-127) then flr((exp-127)/2)
X.i = (X.i & ((1uL << 23) - 1)) | (0x7F << 23);  //Mask mantissa, force exponent to zero.
Y.i = (((uint32_t) e) << 23);
//Error grows with square root of the exponent. Unfortunately no work around like inverse square root... :(
//Y.f *= ScOff * (0x9.5f61ap-4 + X.f*(0x6.a09e68p-4));        //Error = +-1.78e-2 * 2^(flr(log2(x)/2))
//Y.f *= ScOff * (0x7.2181d8p-4 + X.f*(0xa.05406p-4 + X.f*(-0x1.23a14cp-4)));      //Error = +-7.64e-5 * 2^(flr(log2(x)/2))
//Y.f *= ScOff * (0x5.f10e7p-4 + X.f*(0xc.8f2p-4 +X.f*(-0x2.e41a4cp-4 + X.f*(0x6.441e6p-8))));     //Error =  8.21e-5 * 2^(flr(log2(x)/2))
//Y.f *= ScOff * (0x5.32eb88p-4 + X.f*(0xe.abbf5p-4 + X.f*(-0x5.18ee2p-4 + X.f*(0x1.655efp-4 + X.f*(-0x2.b11518p-8)))));   //Error = +-9.92e-6 * 2^(flr(log2(x)/2))
//Y.f *= ScOff * (0x4.adde5p-4 + X.f*(0x1.08448cp0 + X.f*(-0x7.ae1248p-4 + X.f*(0x3.2cf7a8p-4 + X.f*(-0xc.5c1e2p-8 + X.f*(0x1.4b6dp-8))))));   //Error = +-1.38e-6 * 2^(flr(log2(x)/2))
//Y.f *= ScOff * (0x4.4a17fp-4 + X.f*(0x1.22d44p0 + X.f*(-0xa.972e8p-4 + X.f*(0x5.dd53fp-4 + X.f*(-0x2.273c08p-4 + X.f*(0x7.466cb8p-8 + X.f*(-0xa.ac00ep-12)))))));    //Error = +-2.9e-7 * 2^(flr(log2(x)/2))
Y.f *= ScOff * (0x3.fbb3e8p-4 + X.f*(0x1.3b2a3cp0 + X.f*(-0xd.cbb39p-4 + X.f*(0x9.9444ep-4 + X.f*(-0x4.b5ea38p-4 + X.f*(0x1.802f9ep-4 + X.f*(-0x4.6f0adp-8 + X.f*(0x5.c24a28p-12 ))))))));   //Error = +-2.7e-6 * 2^(flr(log2(x)/2))
return(Y.f);
}

Cela n'explique pas comment vous avez calculé les constantes (ce qui est l'algorithme utilisé), il y a une erreur sur x2 qui doit lire xR2, et c'est l'erreur calculée est beaucoup plus grand que les autres réponses. Aussi, j'ai comparé et deux fois plus lent que milianw réponses. Les commentaires n'ont aucun sens (la dernière ligne de commentaire a déclaré que l'erreur est inférieure à la décommenté code, pourquoi?)

InformationsquelleAutor nimig18

Vous devez vous connecter pour publier un commentaire.