Point rapide fixe pow, log, exp et sqrt
https://stackoverflow.com/questions/4657468
-
09-10-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
J'ai une classe de point fixe (10,22) et j'ai besoin d'un prisonnier de guerre, un sqrt, une exp et une fonction de journal.
Hélas, je ne sais pas où commencer même à ce sujet. Quelqu'un peut-il me fournir des liens vers des articles utiles ou, mieux encore, me fournir un code?
Je suppose qu'une fois que j'ai une fonction exp alors il devient relativement facile à mettre en œuvre pow et sqrt comme ils deviennent.
pow( x, y ) => exp( y * log( x ) )
sqrt( x ) => pow( x, 0.5 )
Il est juste les exp et log fonctions que je trouve difficile (comme si je me souviens que quelques-uns de mes règles de journaux, je ne me souviens pas beaucoup d'autre à leur sujet).
On peut supposer que, il y aurait également une méthode plus rapide pour sqrt et pow donc des pointeurs sur ce front serait apprécié même si sa juste dire utiliser les méthodes i schéma ci-dessus.
S'il vous plaît noter:. Cette plate-forme doit être croix et dans le plus pur code C / C ++, donc je ne peux pas utiliser tout assembleur Optimisations
La solution
Une solution très simple consiste à utiliser une approximation déterminée par des tables décent. Vous ne avez pas réellement besoin de beaucoup de données si vous réduisez vos entrées correctement. exp(a)==exp(a/2)*exp(a/2), ce qui signifie que vous vraiment besoin de exp(x) calculate pour 1 < x < 2. Au cours de cette gamme, une approximation runga-Kutta donnerait des résultats raisonnables avec ~ 16 entrées IIRC.
De même, sqrt(a) == 2 * sqrt(a/4) == sqrt(4*a) / 2 qui signifie que vous devez seulement les entrées table pour 1 < a < 4. Log (a) est un peu plus difficile: log(a) == 1 + log(a/e). Ceci est une itération assez lente, mais log (1024) est seulement 6,9 de sorte que vous n'aurez pas beaucoup d'itérations.
Vous userait un algorithme similaire « -premier entier » pour pow: pow(x,y)==pow(x, floor(y)) * pow(x, frac(y)). Cela fonctionne parce que pow(double, int) est trivial (diviser pour régner).
[modifier] Pour la composante intégrale de log(a), il peut être utile de stocker une 1, e, e^2, e^3, e^4, e^5, e^6, e^7 de table afin que vous pouvez réduire log(a) == n + log(a/e^n) par une simple recherche binaire hardcoded d'un dans ce tableau. L'amélioration de 7 à 3 étapes est pas si grand, mais il signifie que vous suffit de diviser une fois par e^n au lieu des temps de n par e.
[modifier 2]
Et pour ce terme dernier log(a/e^n), vous pouvez utiliser log(a/e^n) = log((a/e^n)^8)/8 - chaque itération produit 3 plus de bits par consultation de table . Cela permet de conserver votre taille code et petite table. Ceci est du code généralement pour les systèmes embarqués, et ils n'ont pas de grandes caches.
[modifier 3]
C'est Stil pas intelligent de mon côté. log(a) = log(2) + log(a/2). Vous pouvez simplement stocker la valeur de point fixe log2=0.30102999566, compter le nombre de zéros non significatifs, a de changement de vitesse dans la plage utilisée pour votre table de consultation, et multiplier ce changement de vitesse (nombre entier) par le log2 constante à virgule fixe. Peut-être aussi bas que 3 instructions.
Utilisation e pour l'étape de réduction vous donne juste une « belle » constante de log(e)=1.0 mais qui est fausse optimisation. ,30102999566 est tout aussi bien une constante de 1,0; les deux sont des 32 constantes de bits en virgule fixe 10,22. En utilisant 2 la constante de réduction de la plage vous permet d'utiliser un décalage de bits pour une division.
Vous obtenez toujours le truc de modifier 2, log(a/2^n) = log((a/2^n)^8)/8. Au fond, cela vous obtient un résultat (a + b/8 + c/64 + d/512) * 0.30102999566 - avec b, c, d dans l'intervalle [0,7]. a.bcd est vraiment un nombre octal. Pas une surprise puisque nous avons utilisé 8 comme la puissance. (L'astuce fonctionne aussi bien avec une puissance de 2, 4 ou 16).
[modifier 4]
avait encore une extrémité ouverte. pow(x, frac(y) est juste pow(sqrt(x), 2 * frac(y)) et nous avons un 1/sqrt(x) décent. Cela nous donne l'approche beaucoup plus efficace. Dis binaire frac(y)=0.101, à savoir 1/2 plus 1/8. Ensuite, cela signifie x^0.101 est (x^1/2 * x^1/8). Mais x^1/2 est juste sqrt(x) et x^1/8 est (sqrt(sqrt(sqrt(x))). Enregistrement d'une opération plus, Newton-Raphson NR(x) nous donne 1/sqrt(x) donc nous calculons 1.0/(NR(x)*NR((NR(NR(x))). Nous ne inverser le résultat final, ne pas utiliser la fonction sqrt directement.
Autres conseils
Voici un exemple C Clay mise en œuvre de S. Turner de log base point fixe 2 algorithme [ 1 ]. L'algorithme ne nécessite aucune sorte de table de consultation. Cela peut être utile sur les systèmes où les contraintes de mémoire sont serrés et le processeur ne dispose pas d'un FPU, comme est le cas avec de nombreux micro-contrôleurs. Connexion de base e et une base log 10 sont alors également pris en charge en utilisant la propriété des logarithmes que, pour toute base n
logₘ(x)
logₙ(x) = ───────
logₘ(n)
où, pour cet algorithme, m est égal à 2.
Une belle caractéristique de cette mise en œuvre est qu'il supporte une précision variable: la précision peut être déterminée lors de l'exécution, au détriment de la gamme. La façon dont je l'ai implémenté, le processeur (ou compilateur) doit être capable de faire des mathématiques 64 bits pour la tenue des résultats intermédiaires. Il peut être facilement adaptée pour ne pas nécessiter le support 64 bits, mais la gamme sera réduite.
Lors de l'utilisation de ces fonctions, x devrait être une valeur à virgule fixe mis à l'échelle en fonction de la
precision spécifié. Par exemple, si precision est 16, alors x devrait être mis à l'échelle par 2 ^ 16 (65536). Le résultat est une valeur en virgule fixe avec le même facteur d'échelle comme entrée. Une valeur de retour de INT32_MIN représente l'infini négatif. Une valeur de retour de INT32_MAX indique une erreur et errno sera réglé sur EINVAL, indiquant que la précision d'entrée est invalide.
#include <errno.h>
#include <stddef.h>
#include "log2fix.h"
#define INV_LOG2_E_Q1DOT31 UINT64_C(0x58b90bfc) // Inverse log base 2 of e
#define INV_LOG2_10_Q1DOT31 UINT64_C(0x268826a1) // Inverse log base 2 of 10
int32_t log2fix (uint32_t x, size_t precision)
{
int32_t b = 1U << (precision - 1);
int32_t y = 0;
if (precision < 1 || precision > 31) {
errno = EINVAL;
return INT32_MAX; // indicates an error
}
if (x == 0) {
return INT32_MIN; // represents negative infinity
}
while (x < 1U << precision) {
x <<= 1;
y -= 1U << precision;
}
while (x >= 2U << precision) {
x >>= 1;
y += 1U << precision;
}
uint64_t z = x;
for (size_t i = 0; i < precision; i++) {
z = z * z >> precision;
if (z >= 2U << precision) {
z >>= 1;
y += b;
}
b >>= 1;
}
return y;
}
int32_t logfix (uint32_t x, size_t precision)
{
uint64_t t;
t = log2fix(x, precision) * INV_LOG2_E_Q1DOT31;
return t >> 31;
}
int32_t log10fix (uint32_t x, size_t precision)
{
uint64_t t;
t = log2fix(x, precision) * INV_LOG2_10_Q1DOT31;
return t >> 31;
}
Le code de cette mise en œuvre vit également Github , ainsi que d'un programme échantillon / test qui illustre comment utiliser cette fonction pour calculer et logarithmes d'affichage de numéros lus depuis l'entrée standard.
[ 1 ] CS Turner, "Un rapide binaire Logarithme algorithme" , IEEE signal Processing Mag. , p. 124140, septembre 2010.
Un bon point de départ est livre de Jack Crenshaw, < em> "Math Toolkit pour la programmation en temps réel" . Il a une bonne discussion des algorithmes et implémentations pour différentes fonctions transcendantes.
Vérifier mon point fixe sqrt mise en œuvre en utilisant uniquement les opérations entier. Il était amusant d'inventer. Assez vieux maintenant.
Sinon vérifier CORDIC ensemble des algorithmes. C'est la façon de mettre en œuvre toutes les fonctions que vous avez énumérés et les fonctions trigonométriques.
