Multiplication rapide / division par 2 pour les flotteurs et les doubles (C / C++)
Dans le logiciel que j'écris, je fais des millions de multiplication ou de division par 2 (ou puissances de 2) de mes valeurs. Je voudrais vraiment que ces valeurs soient int afin que je puisse accéder aux opérateurs bitshift
int a = 1;
int b = a<<24
Cependant, je ne peux pas, et je dois rester avec des doubles.
Ma question Est: comme il existe une représentation standard des doubles (signe, exposant, mantisse), existe-t-il un moyen de jouer avec l'exposant pour obtenir des multiplications/divisions rapides par une puissance de 2?
Je peux même supposer que le nombre de bits va être fixé (le logiciel fonctionnera sur des machines qui auront toujours des doubles de 64 bits)
PS: et oui, l'algorithme ne fait que ces opérations. C'est le goulot d'étranglement (c'est déjà multithread).
Edit: ou suis-je complètement trompé et les compilateurs intelligents optimisent déjà les choses pour moi?
Résultats temporaires (avec Qt pour mesurer le temps, exagéré, mais je m'en fous):
#include <QtCore/QCoreApplication>
#include <QtCore/QElapsedTimer>
#include <QtCore/QDebug>
#include <iostream>
#include <math.h>
using namespace std;
int main(int argc, char *argv[])
{
QCoreApplication a(argc, argv);
while(true)
{
QElapsedTimer timer;
timer.start();
int n=100000000;
volatile double d=12.4;
volatile double D;
for(unsigned int i=0; i<n; ++i)
{
//D = d*32; // 200 ms
//D = d*(1<<5); // 200 ms
D = ldexp (d,5); // 6000 ms
}
qDebug() << "The operation took" << timer.elapsed() << "milliseconds";
}
return a.exec();
}
Exécute suggèrent que D = d*(1<<5); et D = d*32; exécuter dans le même temps (200 ms) alors que D = ldexp (d,5); est beaucoup plus lent (6000 ms). Je sais que c'est un micro benchmark, et que soudainement, ma RAM a explosé parce que Chrome a soudainement demandé de calculer Pi dans mon dos chaque fois que je cours ldexp(), donc ce benchmark ne vaut rien. Mais je vais le garder quand même.
De l'autre, j'ai du mal à faire reinterpret_cast<uint64_t *> parce qu'il y a une violation const (semble le mot-clé volatile interfère)
8 réponses
Vous pouvez supposer en toute sécurité le formatage IEEE 754, dont les détails peuvent devenir assez gnarley (esp. quand vous entrez dans les sous-normaux). Dans les cas courants, cependant, cela devrait fonctionner:
const int DOUBLE_EXP_SHIFT = 52;
const unsigned long long DOUBLE_MANT_MASK = (1ull << DOUBLE_EXP_SHIFT) - 1ull;
const unsigned long long DOUBLE_EXP_MASK = ((1ull << 63) - 1) & ~DOUBLE_MANT_MASK;
void unsafe_shl(double* d, int shift) {
unsigned long long* i = (unsigned long long*)d;
if ((*i & DOUBLE_EXP_MASK) && ((*i & DOUBLE_EXP_MASK) != DOUBLE_EXP_MASK)) {
*i += (unsigned long long)shift << DOUBLE_EXP_SHIFT;
} else if (*i) {
*d *= (1 << shift);
}
}
EDIT: après avoir fait un peu de timing, cette méthode est étrangement plus lente que la méthode double sur mon compilateur et ma machine, même dépouillée au code exécuté minimum:
double ds[0x1000];
for (int i = 0; i != 0x1000; i++)
ds[i] = 1.2;
clock_t t = clock();
for (int j = 0; j != 1000000; j++)
for (int i = 0; i != 0x1000; i++)
#if DOUBLE_SHIFT
ds[i] *= 1 << 4;
#else
((unsigned int*)&ds[i])[1] += 4 << 20;
#endif
clock_t e = clock();
printf("%g\n", (float)(e - t) / CLOCKS_PER_SEC);
Dans le DOUBLE_SHIFT se termine en 1.6 secondes, avec une boucle interne de
movupd xmm0,xmmword ptr [ecx]
lea ecx,[ecx+10h]
mulpd xmm0,xmm1
movupd xmmword ptr [ecx-10h],xmm0
Versus 2,4 secondes sinon, avec une boucle interne de:
add dword ptr [ecx],400000h
lea ecx, [ecx+8]
Vraiment inattendu!
EDIT 2: mystère résolu! L'un des changements pour VC11 est maintenant qu'il vectorise toujours les boucles à virgule flottante, forçant efficacement /arch: SSE2, bien que VC10, même avec /arch:SSE2 est encore pire avec 3,0 secondes avec une boucle interne de:
movsd xmm1,mmword ptr [esp+eax*8+38h]
mulsd xmm1,xmm0
movsd mmword ptr [esp+eax*8+38h],xmm1
inc eax
VC10 Sans / arch: SSE2 (même avec /arch:SSE) est de 5,3 secondes... avec 1 / 100ème des itérations!!, boucle intérieure:
fld qword ptr [esp+eax*8+38h]
inc eax
fmul st,st(1)
fstp qword ptr [esp+eax*8+30h]
Je savais que la pile FP x87 était aweful, mais 500 fois pire est un peu ridicule. Vous ne verrez probablement pas ces types d'accélérations convertir, c'est-à-dire les ops matriciels en SSE ou int hacks, car c'est le pire des cas de chargement dans la pile FP, en faisant un op, et en stockant, mais c'est un bon exemple pour pourquoi x87 n'est pas le chemin à PARCOURIR pour quoi que ce soit concerner.
C'est l'une de ces choses très spécifiques à l'application. Il peut aider dans certains cas et pas dans d'autres. (Dans la grande majorité des cas, une multiplication simple est toujours la meilleure.)
"L'intuitif" façon de faire c'est juste pour extraire les bits d'un entier de 64 bits et ajouter la valeur de décalage directement dans l'exposant. (cela fonctionnera tant que vous ne frappez pas NAN ou INF)
Donc quelque chose comme ceci:
union{
uint64 i;
double f;
};
f = 123.;
i += 0x0010000000000000ull;
// Check for zero. And if it matters, denormals as well.
Notez que ce code N'est pas compatible C dans tout façon, et est montré juste pour illustrer l'idée. Toute tentative de mise en œuvre doit être effectuée directement dans assembly ou SSE intrinsics.
Cependant, dans la plupart cas la surcharge de déplacer les données de L'unité FP vers l'unité entière (et retour) coûtera beaucoup plus cher que de simplement faire une multiplication pure et simple. C'est particulièrement le cas pour l'ère pré-SSE où la valeur doit être stockée à partir du FPU x87 dans la mémoire, puis relue dans l'entier inscrire.
À L'ère SSE, L'entier SSE et FP SSE utilisent les mêmes registres ISA (bien qu'ils aient toujours des fichiers de registre séparés). Selon le Agner Fog , Il y a une pénalité de cycle 1 à 2 pour déplacer des données entre les unités D'exécution SSE entières et SSE FP. Donc, le coût est beaucoup mieux que l'ère x87, mais il est toujours là.
Dans l'ensemble, cela dépendra de ce que vous avez d'autre sur votre pipeline. Mais dans la plupart des cas, la multiplication sera toujours plus rapide. J'ai couru dans cette exactement le même problème avant donc je parle de l'expérience de première main.
Maintenant, avec les instructions AVX 256 bits qui ne prennent en charge que les instructions FP, il y a encore moins d'incitation à jouer des tours comme celui-ci.
Que diriez-vous de ldexp ?
Tout compilateur à moitié décent générera du code optimal sur votre plate-forme.
Mais comme le souligne @ Clinton, le simple fait de l'écrire de la manière "évidente" devrait faire aussi bien. Multiplier et diviser par des puissances de deux est un jeu d'enfant pour un compilateur moderne.
Grignoter directement la représentation en virgule flottante, en plus d'être non portable, ne sera presque certainement pas plus rapide (et pourrait bien être plus lent).
Et bien sûr, vous ne devriez pas perdez même du temps à penser à cette question à moins que votre outil de profilage ne vous le dise. Mais le genre de personnes qui écoutent ce conseil n'en aura jamais besoin, et ceux qui en ont besoin n'écouteront jamais.
[mise à jour]
OK, donc je viens d'essayer ldexp avec G++ 4.5.2. L'en-tête cmath l'inline comme un appel à __builtin_ldexp, qui à son tour...
...émet un appel à la fonction libm ldexp. J'aurais pensé que ce builtin serait trivial à optimiser, mais je suppose que les développeurs GCC ne a obtenu autour d'elle.
Donc, multiplier par 1 << p est probablement votre meilleur pari, comme vous l'avez découvert.
Le moyen Le plus rapide de le faire est probablement:
x *= (1 << p);
Ce genre de chose peut simplement être fait en appelant une instruction machine pour ajouter p à l'exposant. Dire au compilateur d'extraire à la place les bits avec un masque et de faire quelque chose manuellement rendra probablement les choses plus lentes, pas plus rapides.
Rappelez-vous, C / c++ n'est pas un langage d'assemblage. L'utilisation d'un opérateur bitshift ne compile pas nécessairement en une opération d'assemblage bitshift, pas nécessairement en utilisant la multiplication compiler à la multiplication. Il y a toutes sortes de choses étranges et merveilleuses qui se passent comme quels registres sont utilisés et quelles instructions peuvent être exécutées simultanément que je ne suis pas assez intelligent pour comprendre. Mais votre compilateur, avec de nombreuses années de connaissances et d'expérience et beaucoup de puissance de calcul, est beaucoup mieux à faire ces jugements.
P. S. gardez à l'esprit que si vos doubles sont dans un tableau ou une autre structure de données plate, votre compilateur peut être vraiment intelligent et utilisez SSE à plusieurs 2, ou même 4 doubles en même temps. Cependant, faire beaucoup de décalage de bits va probablement confondre votre compilateur et empêcher cette optimisation.
Quelles autres opérations cet algorithme nécessite-t-il? Vous pourriez être capable de casser vos flotteurs en paires int (signe/mantisse et magnitude), de faire votre traitement et de les reconstituer à la fin.
Multipliant par 2 peut être remplacé par un ajout: x *= 2 est équivalent à x += x.
La Division par 2 peut être remplacée par la multiplication par 0,5. La Multiplication est généralement beaucoup plus rapide que la division.
Bien qu'il y ait peu/aucun avantage pratique à traiter les puissances de deux spécialement pour le flotteur de types doubles, il y a un cas pour cela pour les typesdouble-double . Double-double multiplication et division est compliquée en général, mais est trivial pour multiplier et diviser par une puissance de deux.
Par exemple pour
typedef struct {double hi; double lo;} doubledouble;
doubledouble x;
x.hi*=2, x.lo*=2; //multiply x by 2
x.hi/=2, x.lo/=2; //divide x by 2
, En fait j'ai surchargé << et >> pour doubledouble, de sorte que c'est l'analogue des nombres entiers.
//x is a doubledouble type
x << 2 // multiply x by four;
x >> 3 // divide x by eight.
Selon ce que vous multipliez, si vous avez des données suffisamment récurrentes, une table de recherche peut fournir de meilleures performances, au détriment de la mémoire.