Dans le calcul numérique, il est souvent nécessaire de mettre les nombres à l'échelle pour être dans une plage sûre.
Par exemple, calculer la distance euclidienne: sqrt (a ^ 2 + b ^ 2) . Ici, si la magnitude de a ou b est trop petite / grande, alors un sous-dépassement / dépassement peut se produire.
Une approche courante pour résoudre ce problème consiste à divisez les nombres par le plus grand nombre de magnitude. Cependant, cette solution est:
J'ai donc pensé qu'au lieu de diviser par le plus grand nombre de magnitude, multiplions-le par un nombre réciproque proche de la puissance de 2. Cela semble une meilleure solution, car:
Donc, j'aimerais créer une petite fonction utilitaire, qui a une logique comme celle-ci (par ^ , je veux dire l'exponentiation):
void getScaler(double value, double &scaler, double &scalerReciprocal) {
int e = <exponent of value>;
if (e<-1022) { scaler=2^-1022; scalerReciprocal = 2^1022; }
} else if (e>1022) { scaler=2^1022; scalerReciprocal = 2^-1022; }
} else { scaler=2^e; scalerReciprocal = 2^(2046-e); }
}
Cette fonction doit renvoyer un scaler et un scalerReciprocal normalisés, les deux étant des nombres de puissance de 2, où scaler est proche de la valeur , et scalerReciprocal est l'inverse de scaler.
Les exposants maximum autorisés pour scaler / scaleReciprocal sont -1022..1022 (je ne veux pas travailler avec un scaler sous-normal, car les nombres sous-normaux peuvent être lents ).
Quel serait un moyen rapide de le faire? Cela peut-il être fait avec des opérations pures en virgule flottante? Ou devrais-je extraire l'exposant de la valeur et utiliser de simples if s pour faire la logique? Y a-t-il une sorte d'astuce pour faire la comparaison avec (-) 1022 rapidement (car la plage est symétrique)?
Remarque: scaler n'a pas besoin d'être le plus proche possible puissance de 2. Si une logique en a besoin, scaler peut être à une petite puissance de 2 de la valeur la plus proche.
3 Réponses :
Vous pouvez utiliser
int get_exp(double *d) {
long long *l = (long long *) d;
return ((*l & (0x7ffLL << 52) )>> 52)-1023 ;
}
La valeur renvoyée est la partie fractionnaire de x et exp est l'exposant (moins le décalage).
Alternativement, le code suivant obtient la partie exposante d'un double.
double frexp (double x, int* exp);
frexp fait plus que nécessaire ici. Et en fait moins, car je dois bloquer exp, puis je dois reconvertir pour obtenir un double . Je ne pense pas que frexp vraiment utilisable dans mon cas, car j'ai besoin de vitesse. Je préfère extraire l'exposant à la main, si c'est la voie à suivre (c'est juste un memcpy , un décalage et un masque).
@geza: frexp n'est qu'un décalage et un masque, standardisés et documentés pour le rendre portable. Si vous souhaitez ajuster l'exposant, associez-le avec ldexp (notez que ldexp ajoute à l'exposant plutôt que de le remplacer)
@BenVoigt: malheureusement non. Consultez le code source. Il gère les nombres nan / inf et sous-normaux. Et il fait une certaine logique, qui pourrait faire partie de ma logique if (exp <-1022) .. . Donc, avec frexp , j'aurais du code redondant. Je ne dis pas que c'est très lent. Mais quand même, il est préférable (pour moi) d'extraire l'exposant manuellement, si c'est la voie à suivre.
La fonction s = get_scale (z) calcule la "puissance de fermeture de 2". Puisque les bits de fraction de s
sont nuls, l'inverse de s est juste une soustraction d'entiers (peu coûteuse): voir la fonction inv_of_scale .
Sur x86, get_scale et inv_of_scale compilez en un assemblage assez efficace avec clang.
Le compilateur clang traduit les opérateurs ternaires en minsd et maxsd ,
voir également le commentaire de Peter Cordes .
Avec gcc, il est légèrement plus efficace de
traduire ces fonctions en intrinsèques x86
code ( get_scale_x86 et inv_of_scale_x86 ), voir Godbolt .
Notez que C autorise explicitement le poinçonnage de type
via une union, alors que C ++ (c ++ 11) n'a pas une telle autorisation
Bien que gcc 8.2 et clang 7.0 ne se plaignent pas de l'union, vous pouvez améliorer
la portabilité C ++ en en utilisant l'astuce memcpy au lieu de la
astuce syndicale. Une telle modification du code devrait être triviale.
Le code doit gérer correctement les sous-normaux.
/* Test how well get_scale vectorizes: */
void get_scale_vec(double * __restrict__ t, double * __restrict__ x){
int n = 1024;
int i;
for (i = 0; i < n; i++){
x[i] = get_scale(t[i]);
}
}
La sortie semble correcte:
Portable code:
x pow_of_2 inverse pow2*inv x*inverse
4.940656e-324 2.225074e-308 4.494233e+307 1.000000e+00 2.220446e-16
1.099790e-320 2.225074e-308 4.494233e+307 1.000000e+00 4.942713e-13
1.100000e-300 7.466109e-301 1.339386e+300 1.000000e+00 1.473324e+00
1.100000e-05 7.629395e-06 1.310720e+05 1.000000e+00 1.441792e+00
7.000000e-01 5.000000e-01 2.000000e+00 1.000000e+00 1.400000e+00
1.700000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.700000e+00
1.231000e+02 6.400000e+01 1.562500e-02 1.000000e+00 1.923437e+00
1.100000e+300 6.696929e+299 1.493222e-300 1.000000e+00 1.642544e+00
1.790000e+308 4.494233e+307 2.225074e-308 1.000000e+00 3.982882e+00
-1.099790e-320 2.225074e-308 4.494233e+307 1.000000e+00 -4.942713e-13
-7.000000e-01 5.000000e-01 2.000000e+00 1.000000e+00 -1.400000e+00
-1.700000e+00 1.000000e+00 1.000000e+00 1.000000e+00 -1.700000e+00
-1.231000e+02 6.400000e+01 1.562500e-02 1.000000e+00 -1.923437e+00
-1.100000e+307 5.617791e+306 1.780059e-307 1.000000e+00 -1.958065e+00
x86 specific SSE code:
x pow_of_2 inverse pow2*inv x*inverse
4.940656e-324 2.225074e-308 4.494233e+307 1.000000e+00 2.220446e-16
1.099790e-320 2.225074e-308 4.494233e+307 1.000000e+00 4.942713e-13
1.100000e-300 7.466109e-301 1.339386e+300 1.000000e+00 1.473324e+00
1.100000e-05 7.629395e-06 1.310720e+05 1.000000e+00 1.441792e+00
7.000000e-01 5.000000e-01 2.000000e+00 1.000000e+00 1.400000e+00
1.700000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.700000e+00
1.231000e+02 6.400000e+01 1.562500e-02 1.000000e+00 1.923437e+00
1.100000e+300 6.696929e+299 1.493222e-300 1.000000e+00 1.642544e+00
1.790000e+308 4.494233e+307 2.225074e-308 1.000000e+00 3.982882e+00
-1.099790e-320 2.225074e-308 4.494233e+307 1.000000e+00 -4.942713e-13
-7.000000e-01 5.000000e-01 2.000000e+00 1.000000e+00 -1.400000e+00
-1.700000e+00 1.000000e+00 1.000000e+00 1.000000e+00 -1.700000e+00
-1.231000e+02 6.400000e+01 1.562500e-02 1.000000e+00 -1.923437e+00
-1.100000e+307 5.617791e+306 1.780059e-307 1.000000e+00 -1.958065e+00
Vectorisation
La fonction get_scale doit vectoriser avec des compilateurs prenant en charge la vectorisation automatique. Le morceau suivant de
code vectorise très bien avec clang (pas besoin d'écrire du code intrinsèque SSE / AVX).
XXX
Malheureusement, gcc ne trouve pas les instructions vmaxpd et vminpd .
Merci d'avoir répondu! Sur la base de votre solution, j'en ai trouvé une (peut-être) plus rapide.
Re: union type-punning: GNU C ++ le prend explicitement en charge en tant qu'extension d'ISO C ++. Voir gcc.gnu.org/onlinedocs/gcc/Optimize-Options .html # Type-punnin g et gcc.gnu.org/onlinedocs/gcc/... . Je pense que MSVC le prend également en charge, mais IDK s'il est documenté. Pourtant, pas d'inconvénient à utiliser memcpy pour de bons compilateurs, tant qu'il s'agit de la largeur totale du type.
@PeterCordes: pour autant que je sache, MSVC n'utilise pas l'optimisation basée sur des règles d'aliasing strictes. Par exemple, ce code recharge * a deux fois inutilement.
@geza: l'alias stricte est lié mais distinct du poinçonnage de type union * (int *) & my_float est toujours UB dans GNU C, bien que cela puisse souvent fonctionner dans la pratique pour des cas simples comme celui-ci où il n'y a pas également d'opérations FP mélangées.
@PeterCordes: sûr :) Je voulais dire qu'avec MSVC, toutes sortes de punitions de type fonctionnent, pas seulement basées sur les unions, comme avec gcc / clang. Je n'ai jamais vu cela documenté nulle part, c'est juste basé sur mon expérience.
@geza: Je vois. Oui, je pense qu'il est vrai que MSVC prend également intentionnellement en charge le casting de pointeurs pour le poinçonnage de type (car beaucoup de code non portable hérité l'utilise). En général, nous devons examiner attentivement une optimisation ou une situation de travail manquée et conclure qu'elle est officiellement prise en charge, mais je pense que dans ce cas, nous pouvons. (Basé sur un usage répandu de l'idiome, et que si MSVC allait casser quelque chose, il le briserait probablement.)
@PeterCordes: Re: punition de type d'union: Merci pour les liens. J'ai fait quelques expériences avec union contre memcpy . Avec clang, MSVC et gcc, je n'ai trouvé aucune différence entre l'assembly généré pour ces deux cas. Avec icc, il y avait une petite différence
Sur la base de la réponse de wim, voici une autre solution, qui peut être plus rapide, car elle contient une instruction de moins. La sortie est un peu différente, mais remplit toujours les conditions.
L'idée est d'utiliser des opérations sur les bits pour corriger les cas de bordure: mettez un 01 au lsb de l'exposant, quel que soit son valeur. Donc, exposant:
00 à 01 ), ou réduit de moitié (lorsque 10-> 01) ou 1/4 (lorsque 11-> 01) Donc, cette routine modifiée fonctionne (et je pense que c'est plutôt cool que le problème puisse être résolu avec seulement 2 instructions fast asm ):
#include<stdio.h>
#include<stdint.h>
#include<immintrin.h>
/* gcc -Wall -m64 -O3 -march=sandybridge dbl_scale.c */
union dbl_int64{
double d;
uint64_t i;
};
double get_scale(double t){
union dbl_int64 x;
uint64_t and_i;
uint64_t or_i;
/* 0xFEDCBA9876543210 */
and_i = 0x7FD0000000000000ull;
or_i = 0x0010000000000000ull;
x.d = t;
x.i = (x.i & and_i)|or_i; /* Set fraction bits to zero, take absolute value */
return x.d;
}
double get_scale_x86(double t){
__m128d x = _mm_set_sd(t);
__m128d x_and = _mm_castsi128_pd(_mm_set1_epi64x(0x7FD0000000000000ull));
__m128d x_or = _mm_castsi128_pd(_mm_set1_epi64x(0x0010000000000000ull));
x = _mm_and_pd(x, x_and);
x = _mm_or_pd(x, x_or);
return _mm_cvtsd_f64(x);
}
/* Compute the inverse 1/t of a double t with all zero fraction bits */
/* and exponent between the limits of function get_scale */
/* A single integer subtraction is much less expensive than a */
/* floating point division. */
double inv_of_scale(double t){
union dbl_int64 x;
/* 0xFEDCBA9876543210 */
uint64_t inv_mask = 0x7FE0000000000000ull;
x.d = t;
x.i = inv_mask - x.i;
return x.d;
}
double inv_of_scale_x86(double t){
__m128i inv_mask = _mm_set1_epi64x(0x7FE0000000000000ull);
__m128d x = _mm_set_sd(t);
__m128i x_i = _mm_sub_epi64(inv_mask, _mm_castpd_si128(x));
return _mm_cvtsd_f64(_mm_castsi128_pd(x_i));
}
int main(){
int n = 14;
int i;
/* Several example values, 4.94e-324 is the smallest subnormal */
double y[14] = { 4.94e-324, 1.1e-320, 1.1e-300, 1.1e-5, 0.7, 1.7, 123.1, 1.1e300,
1.79e308, -1.1e-320, -0.7, -1.7, -123.1, -1.1e307};
double z, s, u;
printf("Portable code:\n");
printf(" x pow_of_2 inverse pow2*inv x*inverse \n");
for (i = 0; i < n; i++){
z = y[i];
s = get_scale(z);
u = inv_of_scale(s);
printf("%14e %14e %14e %14e %14e\n", z, s, u, s*u, z*u);
}
printf("\nx86 specific SSE code:\n");
printf(" x pow_of_2 inverse pow2*inv x*inverse \n");
for (i = 0; i < n; i++){
z = y[i];
s = get_scale_x86(z);
u = inv_of_scale_x86(s);
printf("%14e %14e %14e %14e %14e\n", z, s, u, s*u, z*u);
}
return 0;
}
Quelles 3 instructions comptez-vous? Le vôtre le réduit à 2, juste ANDPS / ORPS au lieu de ANDPS / MINPD / MAXPD. Ou si vous comptez en fait la mise à l'échelle de 2 valeurs en fonction de la magnitude maximale, alors vous aurez besoin de AND + AND + (sélectionnez l'exposant le plus élevé avec MAXPD ou AVX512 VPMAXUQ) + OU + PSUBQ, puis appliquez-le aux deux entrées avec 2x MULPD ou VFMADD ... si vous ou le compilateur pouvez contracter la première étape après la normalisation dans un FMA.
Et BTW, évitez les intrinsèques pour le scalaire; ils sont nuls car il n'y a aucun moyen de dire au compilateur que vous voulez un vecteur avec un élément supérieur non défini, c'est-à-dire pas d'équivalent scalaire-> 128 de __m256 _mm256_castps128_ps256 (__m128 a) . Malheureusement, à part clang, la plupart des compilateurs gaspillent en fait une instruction s'étendant à zéro pour _mm_set_sd (t) . Comment fusionner un scalaire dans un vecteur sans que le compilateur ne gaspille une instruction de remise à zéro des éléments supérieurs? Limitation de conception dans les éléments intrinsèques d'Intel? . Utilisez simplement la version union type-pun, je pense que tous les principaux compilateurs x86 su
Oh, je viens de regarder votre lien Godbolt. C'est étrange que les compilateurs scalaires n'utilisent pas ANDPS / ORPS , même dans une boucle, et extraient réellement les regs GP avec movq . Il n'y a pas de registres XMM préservés des appels dans x86-64 System V, ils ne pouvaient donc pas hisser les constantes, mais les utiliser à partir de la mémoire serait toujours une victoire. Espérons que les compilateurs vectoriseront automatiquement la version C pure.
@PeterCordes: oui, c'est en fait juste deux :) J'étais content d'avoir pu supprimer une instruction de la version de wim, je n'ai pas pris la peine de vérifier qu'en fait la première instruction n'est pas vraiment nécessaire. Merci pour le lien, je vais le vérifier!
Vous voudrez peut-être ajouter une version de votre fonction qui prend 2 entrées et renvoie un facteur d'échelle à utiliser pour les deux, car comme je l'ai dit, vous pouvez combiner cette recherche de magnitude maximale avec la remise à zéro du significande.
Belle solution. Je pense que votre fonction get_scale est un peu plus bancale que la mienne, néanmoins, elle peut bien fonctionner pour votre application.
@wim: rappelez-vous que le cas d'utilisation de l'OP est juste de mettre à l'échelle les valeurs à une magnitude où les mettre au carré ne débordera pas à l'infini ou ne sous-dépassera pas en dénormal / zéro, pour des choses comme sqrt (a ^ 2 + b ^ 2) < / code> et similaire. Laisser le plus grand dans la plage [1..2) au lieu de [0.5 .. 1) est tout à fait correct. Utiliser OR pour appliquer un minimum est une très bonne idée pour ce cas d'utilisation.
@PeterCordes: Avec l'idée ET / OU, les valeurs normales sont mises à l'échelle dans la plage [0,5 ... 8,0), ce qui est parfait pour calculer de manière fiable une hypoténuse, en effet.
Je ne pense pas que cela répondra à votre question, mais si vous recherchez des performances accrues, le moyen le plus rapide de multiplier et de diviser par 2 est de changer de bit. Dans le traitement du signal, vous soustrayez le décalage du signal et, comme vous l'avez dit, mettez-le à l'échelle de sorte que la plage soit [0,1] et vous concevez le filtre pour avoir une magnitude maximale de 1. De plus, je ne comprends pas: vous ne mettez à l'échelle que si l'exposant est inférieur à -1022 ou supérieur à 1022: pourquoi? quelle gamme devrait-il être?
@FrancescoBoi: Non, la mise à l'échelle se produirait toujours. Mais
scalerne devrait avoir que des exposants entre-1022..1022(ce qui représente presque toute la plage. Seules les valeurs limites problématiques sont éliminées).Êtes-vous intéressé par le C pur portable qui compilera efficacement pour x86, ou êtes-vous également intéressé par le C avec des intrinsèques pour les instructions SIMD comme AVX512
_mm512_getexp_pd(extraire l'exposant sous forme dedouble) et_mm512_scalef_pdqui faitdst [63: 0]: = tmp_src1 [63: 0] * POW (2, FLOOR (tmp_src2 [63: 0]) ))(c'est-à-dire ajoute la partie entière d'un double au champ d'exposant d'un autre.)@PeterCordes: Je m'intéresse surtout au C pur, ou peut-être aux extensions, qui sont largement disponibles. AVX n'est toujours pas si répandu. Mais merci de l'avoir mentionné, je ne connaissais pas
_mm512_getexp_pd/_mm512_scalef_pd. Ils semblent similaires à l'ancienFXTRACT/FSCALE.En C pur, si vous finissez par serrer quelque chose, écrivez-le pour qu'il puisse compiler efficacement en
minsd/maxsd(ou vectoriser automatiquement avec leurs versions SIMD). Ce sont SSE2 et donc la ligne de base pour x86-64, et je pense que la plupart des jeux d'instructions SIMD contiennent des instructions min / max. Voir Quelle est l'instruction qui donne le minimum et le maximum de FP sans branche sur x86? . par exemple.a. Vous pouvez simplement serrer séparément le détartreur et son réciproque. Bien sûr, si vous faites des choses avec des exposants entiers, vous devrez peut-être serrer dans la plage d'entiers avant de remplir un double.@PeterCordes: oui, merci pour l'info. J'ai une pince, puis je dois créer une double valeur basée sur la pince. Je m'attends à ce qu'il puisse y avoir un petit truc ici. Si la pince ne modifie pas la valeur, j'ai besoin d'une simple opération
et(effacez simplement le significand). Si la pince se modifie, un peu plus se produit. Il existe peut-être un moyen extrêmement intelligent de compacter cesifen quelque chose d'intelligent.Pouvez-vous supposer que C ++
doubleest représenté par IEEE 754 64 bits en virgule flottante binaire?@PatriciaShanahan: Il existe une macro standard prédéfinie pour vous dire si c'est le cas ou non.
@PatriciaShanahan: oui.
@PeterCordes En ce qui concerne l'opérateur ternaire et
minsd,minpd,maxsd,maxpd, clang fait mieux travail que gcc. Lien GodboltIl me vient à l'esprit qu'une autre solution est d'utiliser une version extrême du split Veltkamp-Dekker:
double get_scale (double x) {double d = x * (0x1p52 + 1); double t = d - x; retour d - t; }. Ce ne sont que trois instructions arithmétiques ordinaires, mais elles ne gèrent pas les valeurs proches de l'extrémité supérieure de la plage finie.@EricPostpischil: c'est une belle astuce pour obtenir une valeur de puissance de 2 proche. Sur les plates-formes où la manipulation de nombres à virgule flottante en tant qu'entiers est lente, cette astuce peut faire partie de la solution. J'ai besoin d'exposants serrés (cela peut être facilement fait avec une comparaison), et de la réciproque également (je ne vois pas de moyen rapide d'avoir cela sans manipuler la valeur comme un entier).