Diviser efficacement la valeur non signée par une puissance de deux, en arrondissant au maximum

Je veux implémenter une division entière non signée par une puissance arbitraire de deux, en arrondissant efficacement. Donc ce que je veux, mathématiquement, c’est ceiling(p/q) 0 . En C, l’implémentation de strawman, qui ne tire pas parti du domaine restreint de q pourrait ressembler à la fonction 1 suivante :

 /** q must be a power of 2, although this version works for any q */ uint64_t divide(uint64_t p, uint64_t q) { uint64_t res = p / q; return p % q == 0 ? res : res + 1; } 

… bien sûr, je ne souhaite pas réellement utiliser la division ou le mod au niveau de la machine, car cela prend beaucoup de cycles, même sur du matériel moderne. Je cherche une réduction de force qui utilise des changements et / ou une autre opération (s) bon marché – en tirant parti du fait que q est une puissance de 2.

Vous pouvez supposer que nous avons une fonction efficace lg(unsigned int x) , qui renvoie le journal de base 2 de x , si x est une puissance de deux.

Un comportement indéfini convient si q est égal à zéro.

Veuillez noter que la solution simple: (p+q-1) >> lg(q) ne fonctionne pas en général – essayez-la avec p == 2^64-100 and q == 256 2 par exemple.

Détails de la plateforme

Je suis intéressé par les solutions en C, susceptibles de bien fonctionner sur diverses plates-formes, mais pour des raisons de concrétisme, l’atsortingbution de la prime et parce que toute discussion définitive sur la performance doit inclure une architecture cible, je vais être précis sur comment je vais les tester:

  • Skylake CPU
  • gcc 5.4.0 avec les indicateurs de compilation -O3 -march=haswell

Utiliser les commandes intégrées à gcc (telles que les instructions bitscan / lead zero) convient très bien, et en particulier, j’ai implémenté la fonction lg() je disais disponible comme suit:

 inline uint64_t lg(uint64_t x) { return 63U - (uint64_t)__builtin_clzl(x); } inline uint32_t lg32(uint32_t x) { return 31U - (uint32_t)__builtin_clz(x); } 

J’ai vérifié que ceux-ci comstacknt jusqu’à une seule instruction bsr , au moins avec -march=haswell , malgré l’implication apparente d’une soustraction. Vous êtes bien entendu libre d’ignorer ces éléments et d’utiliser les autres éléments intégrés que vous souhaitez dans votre solution.

Référence

J’ai écrit un repère pour les réponses existantes et je partagerai et mettrai à jour les résultats au fur et à mesure des changements.

Écrire une bonne référence pour une petite opération potentiellement en ligne est assez difficile. Lorsque du code est inséré dans un site d’appel, une grande partie du travail de la fonction peut disparaître, en particulier lorsqu’elle est dans une boucle 3 .

Vous pouvez simplement éviter tout le problème en ligne en vous assurant que votre code n’est pas en ligne: déclarez-le dans une autre unité de compilation. J’ai essayé avec le binary de bench , mais les résultats sont plutôt inutiles. Presque toutes les implémentations liées à 4 ou 5 cycles par appel, mais même une méthode factice qui ne fait rien d’autre que return 0 prend le même temps. Donc, vous ne mesurez généralement que la charge d’ call + ret . De plus, vous n’utiliserez presque jamais vraiment les fonctions de ce type – à moins que vous ne fassiez des bêtises, elles seront disponibles pour l’inline et cela changera tout .

Ainsi, les deux tests sur lesquels je me concentrerai le plus souvent feront appel à la méthode testée de manière répétée dans une boucle, permettant l’optimisation en ligne, l’optimisation entre fonctions, le levage de la boucle et même la vectorisation.

Il existe deux types de référence: la latence et le débit. La principale différence est que, dans le repère de latence, chaque appel à divide dépend de l’appel précédent. Par conséquent, les appels ne peuvent généralement pas se chevaucher 4 :

 uint32_t bench_divide_latency(uint32_t p, uint32_t q) { uint32_t total = p; for (unsigned i=0; i < ITERS; i++) { total += divide_algo(total, q); q = rotl1(q); } return total; } 

Notez que le total dépend donc de la sortie de chaque appel de division et qu’il s’agit également d’une entrée dans l’appel de divide .

La variante de débit, en revanche, n’alimente pas la sortie d’une division dans la suivante. Cela permet de superposer le travail d’un appel à un autre (à la fois par le compilateur, mais surtout le processeur) et même de permettre la vectorisation:

 uint32_t bench_divide_throughput(uint32_t p, uint32_t q) { uint32_t total = p; for (unsigned i=0; i < ITERS; i++) { total += fname(i, q); q = rotl1(q); } return total; } 

Notez que nous introduisons ici le compteur de boucles comme dividende – ceci est variable , mais cela ne dépend pas de l’appel de division précédent.

En outre, chaque référence a trois types de comportement pour le diviseur, q :

  1. Diviseur constant à la compilation. Par exemple, un appel à divide(p, 8) . Ceci est courant dans la pratique et le code peut être beaucoup plus simple lorsque le diviseur est connu au moment de la compilation.
  2. Diviseur invariant. Ici, le diviseur n’est pas connu au moment de la compilation, mais rest constant pour toute la boucle de benchmarking. Cela permet un sous-ensemble des optimisations que fait la constante de compilation.
  3. Diviseur variable. Le diviseur change à chaque itération de la boucle. Les fonctions de référence ci-dessus illustrent cette variante en utilisant une instruction “rotation à gauche 1” pour modifier le diviseur.

En combinant tout, vous obtenez un total de 6 points de repère distincts.

Résultats

Global

Afin de choisir un meilleur algorithme global, j’ai examiné chacun des 12 sous-ensembles des algorithmes proposés: (latency, throughput) x (constant a, invariant q, variable q) x (32-bit, 64-bit) et un score de 2, 1 ou 0 par sous-test comme suit:

  • Le ou les meilleurs algorithmes (dans la limite de 5% de tolérance) reçoivent un score de 2.
  • Les algorithmes “assez proches” (pas plus de 50% plus lent que le meilleur) reçoivent un score de 1.
  • Les algorithmes restants ont un score de zéro.

Par conséquent, le score total maximum est de 24, mais aucun algorithme ne l’a atteint. Voici le total des résultats:

 ╔═══════════════════════╦═══════╗ ║ Algorithm ║ Score ║ ╠═══════════════════════╬═══════╣ ║ divide_user23_variant ║ 20 ║ ║ divide_chux ║ 20 ║ ║ divide_user23 ║ 15 ║ ║ divide_peter ║ 14 ║ ║ divide_chrisdodd ║ 12 ║ ║ stoke32 ║ 11 ║ ║ divide_chris ║ 0 ║ ║ divide_weather ║ 0 ║ ╚═══════════════════════╩═══════╝ 

Ainsi, pour les besoins de ce code de test spécifique, avec ce compilateur spécifique et sur cette plate-forme , user2357112 “variant” (avec ... + (p & mask) != 0 ) donne les meilleurs résultats, à égalité avec la suggestion de chux (qui fait code identique).

Voici tous les sous-scores qui se résument à ce qui précède:

 ╔══════════════════════════╦═══════╦════╦════╦════╦════╦════╦════╗ ║ ║ Total ║ LC ║ LI ║ LV ║ TC ║ TI ║ TV ║ ╠══════════════════════════╬═══════╬════╬════╬════╬════╬════╬════╣ ║ divide_peter ║ 6 ║ 1 ║ 1 ║ 1 ║ 1 ║ 1 ║ 1 ║ ║ stoke32 ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║ ║ divide_chux ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║ ║ divide_user23 ║ 8 ║ 1 ║ 1 ║ 2 ║ 2 ║ 1 ║ 1 ║ ║ divide_user23_variant ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║ ║ divide_chrisdodd ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║ ║ divide_chris ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ ║ divide_weather ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ ║ 64-bit Algorithm ║ ║ ║ ║ ║ ║ ║ ║ ║ divide_peter_64 ║ 8 ║ 1 ║ 1 ║ 1 ║ 2 ║ 2 ║ 1 ║ ║ div_stoke_64 ║ 5 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 1 ║ ║ divide_chux_64 ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║ ║ divide_user23_64 ║ 7 ║ 1 ║ 1 ║ 2 ║ 1 ║ 1 ║ 1 ║ ║ divide_user23_variant_64 ║ 10 ║ 2 ║ 2 ║ 2 ║ 1 ║ 2 ║ 1 ║ ║ divide_chrisdodd_64 ║ 6 ║ 1 ║ 1 ║ 2 ║ 0 ║ 0 ║ 2 ║ ║ divide_chris_64 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ ║ divide_weather_64 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ 0 ║ ╚══════════════════════════╩═══════╩════╩════╩════╩════╩════╩════╝ 

Ici, chaque test est nommé comme XY, avec X dans {Latency, Throughput} et Y dans {Constant Q, Invariant Q, Variable Q}. Ainsi, par exemple, LC est “Test de latence à q constant”.

Une parsing

Au plus haut niveau, les solutions peuvent être grossièrement divisées en deux catégories: rapide (les 6 premiers) et lent (les 2 derniers). La différence est plus grande: tous les algorithmes rapides étaient les plus rapides sur au moins deux sous-tests et, en général, quand ils ne finissaient pas premiers, ils entraient dans la catégorie “assez proche” (leurs seules exceptions étant les vectorisations ayant échoué dans le cas de stoke et chrisdodd ). Les algorithmes lents ont cependant marqué 0 (même pas proche) à chaque test. Vous pouvez donc éliminer la plupart des algorithmes lents.

Vectorisation automatique

Parmi les algorithmes rapides , la capacité de vectorisation automatique est un facteur de différenciation important.

Aucun des algorithmes n’a été capable de s’auto-vectoriser dans les tests de latence, ce qui est logique puisque les tests de latence sont conçus pour intégrer leurs résultats directement à l’itération suivante. Vous ne pouvez donc vraiment calculer les résultats que de manière sérielle.

Pour les tests de débit, cependant, de nombreux algorithmes ont été capables de s’auto-vectoriser pour le cas constant Q et le cas invariant Q. Dans ces deux tests, le diviseur q est invariant de la boucle (et dans le premier cas, il s’agit d’une constante de compilation). Le dividende est le compteur de boucle, il est donc variable, mais prévisible (et en particulier un vecteur de dividendes peut être calculé de manière sortingviale en ajoutant 8 au vecteur de saisie précédent: [0, 1, 2, ..., 7] + [8, 8, ..., 8] == [8, 9, 10, ..., 15] ).

Dans ce scénario, gcc a pu vectoriser peter , stoke , chux , user23 et user23_variant . Il n’a pas été en mesure de vectoriser chrisdodd pour une raison quelconque, probablement parce qu’il incluait une twig (mais les conditions conditionnelles n’empêchent pas la vectorisation de manière ssortingcte, de nombreuses autres solutions comportant des éléments conditionnels mais toujours vectorisés). L’impact a été énorme: les algorithmes vectorisés ont montré une amélioration de 8 fois le débit par rapport aux variantes qui ne l’étaient pas mais qui étaient par ailleurs rapides.

La vectorisation n’est pas gratuite! Voici les tailles de fonction pour la variante “constante” de chaque fonction, avec le Vec? colonne indiquant si une fonction vectorisée ou non:

 Size Vec? Name 045 N bench_c_div_stoke_64 049 N bench_c_divide_chrisdodd_64 059 N bench_c_stoke32_64 212 Y bench_c_divide_chux_64 227 Y bench_c_divide_peter_64 220 Y bench_c_divide_user23_64 212 Y bench_c_divide_user23_variant_64 

La tendance est claire – les fonctions vectorisées prennent environ 4 fois la taille des fonctions non vectorisées. C’est à la fois parce que les boucles principales sont plus grandes (les instructions vectorielles ont tendance à être plus grandes et plus nombreuses), et parce que la configuration de la boucle, et en particulier le code post-boucle, est beaucoup plus volumineux: par exemple, la version vectorisée requirejs une réduction de additionnez toutes les sums partielles dans un vecteur. Le nombre de boucles étant fixe et un multiple de 8, aucun code de fin n’est généré – mais s’il était variable, le code généré serait encore plus volumineux.

De plus, malgré l’importante amélioration de l’exécution, la vectorisation de gcc est en réalité médiocre. Voici un extrait de la version vectorisée de la routine de Peter:

  on entry: ymm4 == all zeros on entry: ymm5 == 0x00000001 0x00000001 0x00000001 ... 4007a4: c5 ed 76 c4 vpcmpeqd ymm0,ymm2,ymm4 4007ad: c5 fd df c5 vpandn ymm0,ymm0,ymm5 4007b1: c5 dd fa c0 vpsubd ymm0,ymm4,ymm0 4007b5: c5 f5 db c0 vpand ymm0,ymm1,ymm0 

Ce bloc fonctionne indépendamment sur 8 éléments DWORD provenant de ymm2 . Si nous considérons que x est un élément DWORD unique de ymm2 et que y la valeur entrante de ymm1 ces instructions de foud correspondent à:

  x == 0 x != 0 x = x ? 0 : -1; // -1 0 x = x & 1; // 1 0 x = 0 - x; // -1 0 x = y1 & x; // y1 0 

Ainsi, les trois premières instructions pourraient simplement être remplacées par les premières, car les états sont identiques dans les deux cas. Il s’agit donc de deux cycles ajoutés à cette chaîne de dépendance (qui ne fait pas partie d’une boucle) et de deux cycles supplémentaires. Il est évident que les phases d’optimisation de gcc ont en quelque sorte une faible interaction avec le code de vectorisation, car de telles optimisations sortingviales manquent rarement dans le code scalaire. L’examen des autres versions vectorisées montre de la même manière que beaucoup de performances ont chuté sur le sol.

Branches vs sans twigs

Presque toutes les solutions compilées en code sans twig, même si le code C avait des twigs conditionnelles ou explicites. Les portions conditionnelles étaient suffisamment petites pour que le compilateur décide généralement d’utiliser le déplacement conditionnel ou une variante. Une exception est chrisdodd qui a compilé avec une twig (en vérifiant si p == 0 ) dans tous les tests de débit, mais aucun des tests de latence. Voici un exemple typique du test de débit constant q :

 0000000000400e60 : 400e60: 89 f8 mov eax,edi 400e62: ba 01 00 00 00 mov edx,0x1 400e67: eb 0a jmp 400e73  400e69: 0f 1f 80 00 00 00 00 nop DWORD PTR [rax+0x0] 400e70: 83 c2 01 add edx,0x1 400e73: 83 fa 01 cmp edx,0x1 400e76: 74 f8 je 400e70  400e78: 8d 4a fe lea ecx,[rdx-0x2] 400e7b: c1 e9 03 shr ecx,0x3 400e7e: 8d 44 08 01 lea eax,[rax+rcx*1+0x1] 400e82: 81 fa 00 ca 9a 3b cmp edx,0x3b9aca00 400e88: 75 e6 jne 400e70  400e8a: c3 ret 400e8b: 0f 1f 44 00 00 nop DWORD PTR [rax+rax*1+0x0] 

La twig en 400e76 ignore le cas où p == 0 . En fait, le compilateur aurait pu simplement éplucher la première itération (en calculant explicitement son résultat), puis éviter complètement le saut, car après cela, il peut prouver que p != 0 . Dans ces tests, la twig est parfaitement prévisible, ce qui pourrait donner un avantage au code compilé à l’aide d’une twig (car le code de comparaison et de twig est essentiellement hors ligne et presque libre), et explique en grande partie pourquoi chrisdodd gagne le débit, variable q cas.

Résultats détaillés du test

Vous trouverez ici des résultats de tests détaillés et des détails sur les tests eux-mêmes.

Latence

Les résultats ci-dessous testent chaque algorithme sur 1e9 itérations. Les cycles sont calculés simplement en multipliant le temps / appel par la fréquence d’horloge. Vous pouvez généralement supposer que quelque chose comme 4.01 est identique à 4.00 , mais les écarts plus importants tels que 5.11 semblent être réels et reproductibles.

Les résultats de divide_plusq_32 utilisent (p + q - 1) >> lg(q) mais ne sont affichés qu’à titre de référence, car cette fonction échoue pour les grands p + q . Les résultats pour dummy sont une fonction très simple: return p + q et vous permet d’estimer le temps système de référence 5 (l’addition elle-même devrait prendre un cycle au maximum).

 ============================== Bench: Comstack-time constant Q ============================== Function ns/call cycles divide_peter_32 2.19 5.67 divide_peter_64 2.18 5.64 stoke32_32 1.93 5.00 stoke32_64 1.97 5.09 stoke_mul_32 2.75 7.13 stoke_mul_64 2.34 6.06 div_stoke_32 1.94 5.03 div_stoke_64 1.94 5.03 divide_chux_32 1.55 4.01 divide_chux_64 1.55 4.01 divide_user23_32 1.97 5.11 divide_user23_64 1.93 5.00 divide_user23_variant_32 1.55 4.01 divide_user23_variant_64 1.55 4.01 divide_chrisdodd_32 1.95 5.04 divide_chrisdodd_64 1.93 5.00 divide_chris_32 4.63 11.99 divide_chris_64 4.52 11.72 divide_weather_32 2.72 7.04 divide_weather_64 2.78 7.20 divide_plusq_32 1.16 3.00 divide_plusq_64 1.16 3.00 divide_dummy_32 1.16 3.00 divide_dummy_64 1.16 3.00 ============================== Bench: Invariant Q ============================== Function ns/call cycles divide_peter_32 2.19 5.67 divide_peter_64 2.18 5.65 stoke32_32 1.93 5.00 stoke32_64 1.93 5.00 stoke_mul_32 2.73 7.08 stoke_mul_64 2.34 6.06 div_stoke_32 1.93 5.00 div_stoke_64 1.93 5.00 divide_chux_32 1.55 4.02 divide_chux_64 1.55 4.02 divide_user23_32 1.95 5.05 divide_user23_64 2.00 5.17 divide_user23_variant_32 1.55 4.02 divide_user23_variant_64 1.55 4.02 divide_chrisdodd_32 1.95 5.04 divide_chrisdodd_64 1.93 4.99 divide_chris_32 4.60 11.91 divide_chris_64 4.58 11.85 divide_weather_32 12.54 32.49 divide_weather_64 17.51 45.35 divide_plusq_32 1.16 3.00 divide_plusq_64 1.16 3.00 divide_dummy_32 0.39 1.00 divide_dummy_64 0.39 1.00 ============================== Bench: Variable Q ============================== Function ns/call cycles divide_peter_32 2.31 5.98 divide_peter_64 2.26 5.86 stoke32_32 2.06 5.33 stoke32_64 1.99 5.16 stoke_mul_32 2.73 7.06 stoke_mul_64 2.32 6.00 div_stoke_32 2.00 5.19 div_stoke_64 2.00 5.19 divide_chux_32 2.04 5.28 divide_chux_64 2.05 5.30 divide_user23_32 2.05 5.30 divide_user23_64 2.06 5.33 divide_user23_variant_32 2.04 5.29 divide_user23_variant_64 2.05 5.30 divide_chrisdodd_32 2.04 5.30 divide_chrisdodd_64 2.05 5.31 divide_chris_32 4.65 12.04 divide_chris_64 4.64 12.01 divide_weather_32 12.46 32.28 divide_weather_64 19.46 50.40 divide_plusq_32 1.93 5.00 divide_plusq_64 1.99 5.16 divide_dummy_32 0.40 1.05 divide_dummy_64 0.40 1.04 

Débit

Voici les résultats pour les tests de débit. Notez que beaucoup des algorithmes ici ont été auto-vectorisés, donc la performance est relativement très bonne pour ceux-ci: une fraction de cycle dans de nombreux cas. Un résultat est que, contrairement à la plupart des résultats de latence, les fonctions 64 bits sont considérablement plus lentes, car la vectorisation est plus efficace avec des tailles d’éléments plus petites (bien que l’écart soit plus important que prévu).

 ============================== Bench: Comstack-time constant Q ============================== Function ns/call cycles stoke32_32 0.39 1.00 divide_chux_32 0.15 0.39 divide_chux_64 0.53 1.37 divide_user23_32 0.14 0.36 divide_user23_64 0.53 1.37 divide_user23_variant_32 0.15 0.39 divide_user23_variant_64 0.53 1.37 divide_chrisdodd_32 1.16 3.00 divide_chrisdodd_64 1.16 3.00 divide_chris_32 4.34 11.23 divide_chris_64 4.34 11.24 divide_weather_32 1.35 3.50 divide_weather_64 1.35 3.50 divide_plusq_32 0.10 0.26 divide_plusq_64 0.39 1.00 divide_dummy_32 0.08 0.20 divide_dummy_64 0.39 1.00 ============================== Bench: Invariant Q ============================== Function ns/call cycles stoke32_32 0.48 1.25 divide_chux_32 0.15 0.39 divide_chux_64 0.48 1.25 divide_user23_32 0.17 0.43 divide_user23_64 0.58 1.50 divide_user23_variant_32 0.15 0.38 divide_user23_variant_64 0.48 1.25 divide_chrisdodd_32 1.16 3.00 divide_chrisdodd_64 1.16 3.00 divide_chris_32 4.35 11.26 divide_chris_64 4.36 11.28 divide_weather_32 5.79 14.99 divide_weather_64 17.00 44.02 divide_plusq_32 0.12 0.31 divide_plusq_64 0.48 1.25 divide_dummy_32 0.09 0.23 divide_dummy_64 0.09 0.23 ============================== Bench: Variable Q ============================== Function ns/call cycles stoke32_32 1.16 3.00 divide_chux_32 1.36 3.51 divide_chux_64 1.35 3.50 divide_user23_32 1.54 4.00 divide_user23_64 1.54 4.00 divide_user23_variant_32 1.36 3.51 divide_user23_variant_64 1.55 4.01 divide_chrisdodd_32 1.16 3.00 divide_chrisdodd_64 1.16 3.00 divide_chris_32 4.02 10.41 divide_chris_64 3.84 9.95 divide_weather_32 5.40 13.98 divide_weather_64 19.04 49.30 divide_plusq_32 1.03 2.66 divide_plusq_64 1.03 2.68 divide_dummy_32 0.63 1.63 divide_dummy_64 0.66 1.71 

a Au moins en spécifiant unsigned, nous évitons toute la boîte de vers liés au comportement de décalage à droite des entiers signés en C et C ++.

0 Bien entendu, cette notation ne fonctionne pas réellement en C où / tronque le résultat et le ceiling ne fait rien. Alors considérez cette pseudo-notation plutôt que le simple C.

1 Je suis également intéressé par les solutions où tous les types sont uint32_t plutôt que uint64_t .

2 En général, tout p et qp + q >= 2^64 pose un problème en raison du débordement.

3 Cela dit, la fonction devrait être dans une boucle, car la performance d’une fonction microscopique qui prend une demi-douzaine de cycles n’a d’importance que si elle est appelée dans une boucle assez étroite.

4 Ceci est un peu une simplification – seul le dividende p dépend du résultat de l’itération précédente, de sorte que certains travaux liés au traitement de q peuvent toujours se chevaucher.

5 Utilisez toutefois ces estimations avec prudence – les frais généraux ne sont pas simplement additifs. Si le temps système est égal à 4 cycles et qu’une fonction f prend 5, il est probablement inexact de dire que le coût du travail réel dans f est 5 - 4 == 1 , en raison du chevauchement de l’exécution.

Cette réponse concerne ce qui est idéal dans asm; ce que nous aimerions convaincre le compilateur d’émettre pour nous. (Je ne suggère pas réellement d’utiliser inline asm, sauf comme sharepoint comparaison lors de l’parsing comparative des résultats du compilateur. https://gcc.gnu.org/wiki/DontUseInlineAsm ).

J’ai réussi à obtenir une très bonne sortie asm de C pur pour ceil_div_andmask , voir ci-dessous. (C’est pire qu’un CMOV sur Broadwell / Skylake, mais probablement bon sur Haswell. Néanmoins, la version user23 / chux est encore meilleure dans les deux cas.) C’est un des rares cas où j’ai eu le compilateur à émettre. le asm je voulais.


Il semble que l’idée générale de Chris Dodd concernant le return ((p-1) >> lg(q)) + 1 avec la gestion des cas spéciaux pour d = 0 soit l’une des meilleures options. C’est-à-dire que sa mise en œuvre optimale dans asm est difficile à battre avec une mise en œuvre optimale de toute autre chose. Les chux (p >> lg(q)) + (bool)(p & (q-1)) aussi des avantages (comme une latence plus faible à partir de p-> résultat), et plus de CSE lorsque le même q est utilisé pour plusieurs divisions. Voir ci-dessous pour une parsing de latence / débit de ce que gcc en fait.

Si le même e = lg(q) est réutilisé pour plusieurs dividendes, ou si le même dividende est réutilisé pour plusieurs diviseurs, différentes implémentations peuvent donner une plus grande place à l’expression. Ils peuvent également vectoriser efficacement avec AVX2 .

Les twigs ne coûtent pas cher et sont très efficaces si elles prédisent très bien . Par conséquent, il est préférable de d==0 twig sur d==0 si elle n’est presque jamais prise. Si d==0 n’est pas rare, l’asym sans twigs se comportera mieux en moyenne. Dans l’idéal, nous pouvons écrire quelque chose en C qui permettra à gcc de faire le bon choix lors de l’optimisation guidée par le profil et de le comstackr correctement dans les deux cas.

Étant donné que les meilleures implémentations asm sans twig n’apportent pas beaucoup de temps de latence par rapport à une implémentation en twig, le mode sans twig est probablement la solution à moins que la twig ne fonctionne de la même manière, peut-être 99% du temps. Cela pourrait être probable pour les twigs sur p==0 , mais probablement moins probable pour les twigs sur p & (q-1) .


Il est difficile de guider gcc5.4 dans l’émission de valeurs optimales. Ceci est mon travail en cours sur Godbolt ).

Je pense que les séquences optimales pour Skylake pour cet algorithme sont les suivantes. (Présenté en tant que fonctions autonomes pour l’ABI AMD64 SysV, mais parle de débit / latence en supposant que le compilateur émettra quelque chose de similaire dans une boucle, sans RET attaché).


La twig continue de calculer d-1 pour détecter d == 0 , au lieu d’un test et d’une twig séparés. Réduit bien le nombre d’uop, en particulier sur la famille SnB où JC peut fusionner avec SUB.

 ceil_div_pjc_branch: xor eax,eax ; can take this uop off the fast path by adding a separate xor-and-return block, but in reality we want to inline something like this. sub rdi, 1 jc .d_was_zero ; fuses with the sub on SnB-family tzcnt rax, rsi ; tzcnt rsi,rsi also avoids any false-dep problems, but this illustrates that the q input can be read-only. shrx rax, rdi, rax inc rax .d_was_zero: ret 
  • Uops du domaine fondu: 5 (sans compter ret) et l’un d’eux est un xor-zéro (pas d’unité d’exécution)
  • Latence HSW / SKL avec prédiction de twig réussie:
    • (d == 0): Aucune dépendance des données sur d ou q ne rompt la chaîne dep . (contrôlez la dépendance de d pour détecter les erreurs de prédiction et retirer la twig).
    • (d! = 0): q-> résultat: tzcnt + shrx + inc = 5c
    • (d! = 0): d-> résultat: sub + shrx + inc = 3c
  • Débit: probablement juste goulot d’étranglement sur le débit uop

J’ai essayé mais je n’ai pas réussi à obtenir que gcc se twig sur CF à partir de la soustraction, mais il souhaite toujours faire une comparaison séparée. Je sais que gcc peut être amené à créer des twigs sur CF après avoir soustrait deux variables, mais peut-être que cela échoue si l’une est une constante de compilation. (IIRC, cela comstack typiquement à un test CF avec vars non signés: foo -= bar; if(foo>bar) carry_detected = 1; )


Sans twig avec ADC / SBB pour gérer le cas d==0 . La gestion du zéro ajoute une seule instruction au chemin critique (par rapport à une version sans traitement spécial pour d == 0), mais convertit également une autre d’un INC en une sbb rax, -1 pour que CF annule le -= -1 L’utilisation d’un CMOV est moins chère avant Broadwell, mais nécessite des instructions supplémentaires pour le configurer.

 ceil_div_pjc_asm_adc: tzcnt rsi, rsi sub rdi, 1 adc rdi, 0 ; d? d-1 : d. Sets CF=CF shrx rax, rdi, rsi sbb rax, -1 ; result++ if d was non-zero ret 
  • Uops de domaine fondu: 5 (sans compter ret) sur SKL. 7 sur HSW
  • Latence SKL:
    • q-> résultat: tzcnt + shrx + sbb = 5c
    • d-> résultat: sub + adc + shrx (dep sur q commence ici) + sbb = 4c
  • Débit: TZCNT s’exécute sur p1. Les CFF, ADC et SHRX ne fonctionnent que sur p06. Donc, je pense que nous avons un goulot d’étranglement sur 3 uops pour p06 par itération, ce qui donne au mieux une itération par 1.5c .

Si q et d sont prêts en même temps, notez que cette version peut exécuter SUB / ADC en parallèle avec la latence 3c de TZCNT. Si les deux proviennent de la même ligne de cache cache-miss, c’est certainement possible. Dans tous les cas, introduire le dep sur q le plus tard possible dans la chaîne de dépendance d-> résultat est un avantage.

Obtenir cela depuis C semble peu probable avec gcc5.4. Il existe un élément insortingnsèque pour l’add-with-carry, mais gcc en fait un gâchis total. Il n’utilise pas d’opérandes immédiats pour ADC ou SBB, et stocke le report dans un registre entier entre chaque opération. gcc7, clang3.9 et icc17 en font un code terrible.

 #include  // comstacks to completely horrible code, putting the flags into integer regs between ops. T ceil_div_adc(T d, T q) { T e = lg(q); unsigned long long dm1; // unsigned __int64 unsigned char CF = _addcarry_u64(0, d, -1, &dm1); CF = _addcarry_u64(CF, 0, dm1, &dm1); T shifted = dm1 >> e; _subborrow_u64(CF, shifted, -1, &dm1); return dm1; } # gcc5.4 -O3 -march=haswell mov rax, -1 tzcnt rsi, rsi add rdi, rax setc cl xor edx, edx add cl, -1 adc rdi, rdx setc dl shrx rdi, rdi, rsi add dl, -1 sbb rax, rdi ret 

CMOV pour corriger le résultat entier : pire latence de q-> résultat, car elle est utilisée plus tôt dans la chaîne d-> result.

 ceil_div_pjc_asm_cmov: tzcnt rsi, rsi sub rdi, 1 shrx rax, rdi, rsi lea rax, [rax+1] ; inc preserving flags cmovc rax, zeroed_register ret 
  • Uops de domaine fondu: 5 (sans compter ret) sur SKL. 6 sur HSW
  • Latence SKL:
    • q-> résultat: tzcnt + shrx + lea + cmov = 6c (pire que ADC / SBB de 1c)
    • d-> résultat: sub + shrx (dep sur q commence ici) + lea + cmov = 4c
  • Débit: TZCNT s’exécute sur p1. LEA est p15. CMOV et SHRX sont p06. SUB est p0156. En théorie, uniquement goulot d’étranglement sur le débit uop du domaine fondu, donc une itération par 1.25c . Avec de nombreuses opérations indépendantes, les conflits de ressources liés au vol de p1 ou de p06 par SUB ou LEA ne devraient pas être un problème de débit, car à 1 iter par 1.25c, aucun port n’est saturé avec des uops pouvant être exécutés uniquement sur ce port.

CMOV pour obtenir un opérande pour SUB : J’espérais pouvoir trouver un moyen de créer un opérande pour une instruction ultérieure qui produirait un zéro si nécessaire, sans dépendance d’entrée sur q, e ou le résultat SHRX. Cela aiderait si d est prêt avant q ou en même temps.

Cela n’atteint pas cet objective et nécessite un mov rdx,-1 supplémentaire de 7 octets mov rdx,-1 dans la boucle.

 ceil_div_pjc_asm_cmov: tzcnt rsi, rsi mov rdx, -1 sub rdi, 1 shrx rax, rdi, rsi cmovnc rdx, rax sub rax, rdx ; res += d ? 1 : -res ret 

Version à faible temps de latence pour les CPU pré-BDW avec une CMOV coûteuse , en utilisant SETCC pour créer un masque pour AND.

 ceil_div_pjc_asm_setcc: xor edx, edx ; needed every iteration tzcnt rsi, rsi sub rdi, 1 setc dl ; d!=0 ? 0 : 1 dec rdx ; d!=0 ? -1 : 0 // AND-mask shrx rax, rdi, rsi inc rax and rax, rdx ; zero the bogus result if d was initially 0 ret 

Encore 4c latence à partir de d-> résultat (et 6 à partir de q-> résultat), car le SETC / DEC se produit en parallèle avec le SHRX / INC. Nombre total d’uop: 8. La plupart de ces insns peuvent fonctionner sur n’importe quel port, il devrait donc être 1 iter par 2 horloges.

Bien sûr, pour les versions antérieures à HSW, vous devez également remplacer SHRX.

Nous pouvons obtenir que gcc5.4 émette quelque chose de presque aussi bon: (utilise toujours un sub rdi, 1 distinct au lieu de définir un masque basé sur sub rdi, 1 , mais sinon les mêmes instructions que ci-dessus). Voir sur Godbolt .

 T ceil_div_andmask(T p, T q) { T mask = -(T)(p!=0); // TEST+SETCC+NEG T e = lg(q); T nonzero_result = ((p-1) >> e) + 1; return nonzero_result & mask; } 

Lorsque le compilateur sait que p est non nul, il en profite pour créer du code agréable:

 // http://stackoverflow.com/questions/40447195/can-i-hint-the-optimizer-by-giving-the-range-of-an-integer #if defined(__GNUC__) && !defined(__INTEL_COMPILER) #define assume(x) do{if(!(x)) __builtin_unreachable();}while(0) #else #define assume(x) (void)(x) // still evaluate it once, for side effects in case anyone is insane enough to put any inside an assume() #endif T ceil_div_andmask_nonzerop(T p, T q) { assume(p!=0); return ceil_div_andmask(p, q); } # gcc5.4 -O3 -march=haswell xor eax, eax # gcc7 does tzcnt in-place instead of wasting an insn on this sub rdi, 1 tzcnt rax, rsi shrx rax, rdi, rax add rax, 1 ret 

Chux / user23_variant

seule une latence de 3c à partir de p-> résultat et une constante q peuvent beaucoup CSE.

 T divide_A_chux(T p, T q) { bool round_up = p & (q-1); // comstacks differently from user23_variant with clang: AND instead of return (p >> lg(q)) + round_up; } xor eax, eax # in-place tzcnt would save this xor edx, edx # target for setcc tzcnt rax, rsi sub rsi, 1 test rsi, rdi shrx rdi, rdi, rax setne dl lea rax, [rdx+rdi] ret 

Faire le SETCC avant TZCNT permettrait un TZCNT sur place, en sauvegardant le xor eax,eax . Je n’ai pas regardé comment cela s’inscrit dans une boucle.

  • Uops de domaine fondu: 8 (sans compter ret) sur HSW / SKL
  • Latence HSW / SKL:
    • q-> résultat: (tzcnt + shrx (p) | sub + test (p) + setne) + lea (ou add) = 5c
    • d-> résultat: test (dep sur q commence ici) + setne + lea = 3c . (la chaîne de fils shrx-> est plus courte et donc pas le chemin critique)
  • Débit: Probablement juste goulot d’étranglement sur le front-end, à un iter par 2c . Saving the xor eax,eax should speed this up to one per 1.75c (but of course any loop overhead will be part of the bottleneck, because frontend bottlenecks are like that).
 uint64_t exponent = lg(q); uint64_t mask = q - 1; // v divide return (p >> exponent) + (((p & mask) + mask) >> exponent) // ^ round up 

The separate computation of the “round up” part avoids the overflow issues of (p+q-1) >> lg(q) . Depending on how smart your comstackr is, it might be possible to express the “round up” part as ((p & mask) != 0) without branching.

The efficient way of dividing by a power of 2 for an unsigned integer in C is a right shift — shifting right one divides by two (rounding down), so shifting right by n divides by 2 n (rounding down).

Now you want to round up rather than down, which you can do by first adding 2 n -1, or equivalently subtracting one before the shift and adding one after (except for 0). This works out to something like:

 unsigned ceil_div(unsigned d, unsigned e) { /* compute ceil(d/2**e) */ return d ? ((d-1) >> e) + 1 : 0; } 

The conditional can be removed by using the boolean value of d for addition and subtraction of one:

 unsigned ceil_div(unsigned d, unsigned e) { /* compute ceil(d/2**e) */ return ((d - !!d) >> e) + !!d; } 

Due to its size, and the speed requirement, the function should be made static inline. It probably won’t make a different for the optimizer, but the parameters should be const. If it must be shared among many files, define it in a header:

 static inline unsigned ceil_div(const unsigned d, const unsigned e){... 

Efficiently dividing unsigned value by a power of two, rounding up

[Re-write] given OP’s clarification concerning power-of-2 .

The round-up or ceiling part is easy when overflow is not a concern. Simple add q-1 , then shift.

Otherwise as the possibility of rounding depends on all the bits of p smaller than q , detection of those bits is needed first before they are shifted out.

 uint64_t divide_A(uint64_t p, uint64_t q) { bool round_up = p & (q-1); return (p >> lg64(q)) + round_up; } 

This assumes code has an efficient lg64(uint64_t x) function, which returns the base-2 log of x , if x is a power-of-two.`

My old answer didn’t work if p was one more than a power of two (whoops). So my new solution, using the __builtin_ctzll() and __builtin_ffsll() functions 0 available in gcc (which as a bonus, provides the fast logarithm you mentioned!):

 uint64_t divide(uint64_t p,uint64_t q) { int lp=__builtin_ffsll(p); int lq=__builtin_ctzll(q); return (p>>lq)+(lp<(lq+1)&&lp); } 

Note that this is assuming that a long long is 64 bits. It has to be tweaked a little otherwise.

The idea here is that if we need an overflow if and only if p has fewer trailing zeroes than q . Note that for a power of two, the number of trailing zeroes is equal to the logarithm, so we can use this builtin for the log as well.

The &&lp part is just for the corner case where p is zero: otherwise it will output 1 here.

0 Can't use __builtin_ctzll() for both because it is undefined if p==0 .

If the dividend/divisor can be guaranteed not to exceed 63 (or 31) bits, you can use the following version mentioned in the question. Note how p+q could overflow if they use all 64 bit. This would be fine if the SHR instruction shifted in the carry flag, but AFAIK it doesn’t.

 uint64_t divide(uint64_t p, uint64_t q) { return (p + q - 1) >> lg(q); } 

If those constraints cannot be guaranteed, you can just do the floor method and then add 1 if it would round up. This can be determined by checking if any bits in the dividend are within the range of the divisor. Note: p&-p extracts the lowest set bit on 2s complement machines or the BLSI instruction

 uint64_t divide(uint64_t p, uint64_t q) { return (p >> lg(q)) + ( (p & -p ) < q ); } 

Which clang comstacks to:

  bsrq %rax, %rsi shrxq %rax, %rdi, %rax blsiq %rdi, %rcx cmpq %rsi, %rcx adcq $0, %rax retq 

That's a bit wordy and uses some newer instructions, so maybe there is a way to use the carry flag in the original version. Lets see: The RCR instruction does but seems like it would be worse ... perhaps the SHRD instruction... It would be something like this (unable to test at the moment)

  xor edx, edx ;edx = 0 (will store the carry flag) bsr rcx, rsi ;rcx = lg(q) ... could be moved anywhere before shrd lea rax, [rsi-1] ;rax = q-1 (adding p could carry) add rax, rdi ;rax += p (handle carry) setc dl ;rdx = carry flag ... or xor rdx and setc shrd rax, rdx, cl ;rax = rdx:rax >> cl ret 

1 more instruction, but should be compatible with older processors (if it works ... I'm always getting a source/destination swapped - feel free to edit)

Addenda:

I've implemented the lg() function I said was available as follows:

 inline uint64_t lg(uint64_t x) { return 63U - (uint64_t)__builtin_clzl(x); } inline uint32_t lg32(uint32_t x) { return 31U - (uint32_t)__builtin_clz(x); } 

The fast log functions don't fully optimize to bsr on clang and ICC but you can do this:

 #if defined(__x86_64__) && (defined(__clang__) || defined(__INTEL_COMPILER)) static inline uint64_t lg(uint64_t x){ inline uint64_t ret; //other comstackrs may want bsrq here __asm__("bsr %0, %1":"=r"(ret):"r"(x)); return ret; } #endif #if defined(__i386__) && (defined(__clang__) || defined(__INTEL_COMPILER)) static inline uint32_t lg32(uint32_t x){ inline uint32_t ret; __asm__("bsr %0, %1":"=r"(ret):"r"(x)); return ret; } #endif 

There has already been a lot of human brainpower applied to this problem, with several variants of great answers in C along with Peter Cordes’s answer which covers the best you could hope for in asm, with notes on trying to map it back to C .

So while the humans are having their kick at the can, I thought see what some brute computing power has to say! To that end, I used Stanford’s STOKE superoptimizer to try to find good solutions to the 32-bit and 64-bit versions of this problem.

Usually, a superoptimizer is usually something like a brute force search through all possible instruction sequences until you find the best one by some mesortingc. Of course, with something like 1,000 instructions that will quickly spiral out of control for more than a few instructions 1 . STOKE, on the hand, takes a guided randomized approach: it randomly makes mutations to an existing candidate program, evaluating at each step a cost function that takes both performance and correctness into effect. That’s the one-liner anyway – there are plenty of papers if that stoked your curiosity.

So within a few minutes STOKE found some pretty interesting solutions. It found almost all the high-level ideas in the existing solutions, plus a few unique ones. For example, for the 32-bit function, STOKE found this version:

 neg rsi dec rdi pext rax, rsi, rdi inc eax ret 

It doesn’t use any leading/trailing-zero count or shift instructions at all. Pretty much, it uses neg esi to turn the divisor into a mask with 1s in the high bits, and then pext effectively does the shift using that mask. Outside of that sortingck it’s using the same sortingck that user QuestionC used: decrement p, shift, increment p – but it happens to work even for zero dividend because it uses 64-bit registers to distinguish the zero case from the MSB-set large p case.

I added the C version of this algorithm to the benchmark, and added it to the results. It’s competitive with the other good algorithms, tying for first in the “Variable Q” cases. It does vectorize, but not as well as the other 32-bit algorithms, because it needs 64-bit math and so the vectors can process only half as many elements at once.

Even better, in the 32-bit case it came up with a variety of solutions which use the fact that some of the intuitive solutions that fail for edge cases happen to “just work” if you use 64-bit ops for part of it. Par exemple:

 tzcntl ebp, esi dec esi add rdi, rsi sarx rax, rdi, rbp ret 

That’s the equivalent of the return (p + q - 1) >> lg(q) suggestion I mentioned in the question. That doesn’t work in general since for large p + q it overflows, but for 32-bit p and q this solution works great by doing the important parts in 64-bit. Convert that back into C with some casts and it actually figures out that using lea will do the addition in one instruction 1 :

 stoke_32(unsigned int, unsigned int): tzcnt edx, esi mov edi, edi ; goes away when inlining mov esi, esi ; goes away when inlining lea rax, [rsi-1+rdi] shrx rax, rax, rdx ret 

So it’s a 3-instruction solution when inlined into something that already has the values zero-extended into rdi and rsi . The stand-alone function definition needs the mov instructions to zero-extend because that’s how the SysV x64 ABI works .


For the 64-bit function it didn’t come up with anything that blows away the existing solutions but it did come up with some neat stuff, like:

  tzcnt r13, rsi tzcnt rcx, rdi shrx rax, rdi, r13 cmp r13b, cl adc rax, 0 ret 

That guy counts the trailing zeros of both arguments, and then adds 1 to the result if q has fewer trailing zeros than p , since that’s when you need to round up. Intelligent!

In general, it understood the idea that you needed to shl by the tzcnt really quickly (just like most humans) and then came up with a ton of other solutions to the problem of adjusting the result to account for rounding. It even managed to use blsi and bzhi in several solutions. Here’s a 5-instruction solution it came up with:

 tzcnt r13, rsi shrx rax, rdi, r13 imul rsi, rax cmp rsi, rdi adc rax, 0 ret 

It’s a basically a “multiply and verify” approach – take the truncated res = p \ q , multiply it back and if it’s different than p add one: return res * q == p ? ret : ret + 1 . Cool. Not really better than Peter’s solutions though. STOKE seems to have some flaws in its latency calculation – it thinks the above has a latency of 5 – but it’s more like 8 or 9 depending on the architecture. So it sometimes narrows in solutions that are based on its flawed latency calculation.


1 Interestingly enough though this brute force approach reaches its feasibility around 5-6 instructions: if you assume you can sortingm the instruction count to say 300 by eliminating SIMD and x87 instructions, then you would need ~28 days to try all 300 ^ 5 5 instruction sequences at 1,000,000 candidates/second. You could perhaps reduce that by a factor of 1,000 with various optimizations, meaning less than an hour for 5-instruction sequences and maybe a week for 6-instruction. As it happens, most of the best solutions for this problem fall into that 5-6 instruction window…

2 This will be a slow lea , however, so the sequence found by STOKE was still optimal for what I optimized for, which was latency.

You can do it like this, by comparing dividing n / d with dividing by (n-1) / d .

 #include  int main(void) { unsigned n; unsigned d; unsigned q1, q2; double actual; for(n = 1; n < 6; n++) { for(d = 1; d < 6; d++) { actual = (double)n / d; q1 = n / d; if(n) { q2 = (n - 1) / d; if(q1 == q2) { q1++; } } printf("%u / %u = %u (%f)\n", n, d, q1, actual); } } return 0; } 

Sortie du programme:

 1 / 1 = 1 (1.000000) 1 / 2 = 1 (0.500000) 1 / 3 = 1 (0.333333) 1 / 4 = 1 (0.250000) 1 / 5 = 1 (0.200000) 2 / 1 = 2 (2.000000) 2 / 2 = 1 (1.000000) 2 / 3 = 1 (0.666667) 2 / 4 = 1 (0.500000) 2 / 5 = 1 (0.400000) 3 / 1 = 3 (3.000000) 3 / 2 = 2 (1.500000) 3 / 3 = 1 (1.000000) 3 / 4 = 1 (0.750000) 3 / 5 = 1 (0.600000) 4 / 1 = 4 (4.000000) 4 / 2 = 2 (2.000000) 4 / 3 = 2 (1.333333) 4 / 4 = 1 (1.000000) 4 / 5 = 1 (0.800000) 5 / 1 = 5 (5.000000) 5 / 2 = 3 (2.500000) 5 / 3 = 2 (1.666667) 5 / 4 = 2 (1.250000) 5 / 5 = 1 (1.000000) 

Mettre à jour

I posted an early answer to the original question, which works, but did not consider the efficiency of the algorithm, or that the divisor is always a power of 2. Performing two divisions was needlessly expensive.

I am using MSVC 32-bit comstackr on a 64-bit system, so there is no chance that I can provide a best solution for the required target. But it is an interesting question so I have dabbled around to find that the best solution will discover the bit of the 2**n divisor. Using the library function log2 worked but was so slow. Doing my own shift was much better, but still my best C solution is

 unsigned roundup(unsigned p, unsigned q) { return p / q + ((p & (q-1)) != 0); } 

My inline 32-bit assembler solution is faster, but of course that will not answer the question. I steal some cycles by assuming that eax is returned as the function value.

 unsigned roundup(unsigned p, unsigned q) { __asm { mov eax,p mov edx,q bsr ecx,edx ; cl = bit number of q dec edx ; q-1 and edx,eax ; p & (q-1) shr eax,cl ; divide p by q, a power of 2 sub edx,1 ; generate a carry when (p & (q-1)) == 0 cmc adc eax,0 ; add 1 to result when (p & (q-1)) != 0 } } ; eax returned as function value 

This seems efficient and works for signed if your comstackr is using arithmetic right shifts (usually true).

 #include  int main (void) { for (int i = -20; i <= 20; ++i) { printf ("%5d %5d\n", i, ((i - 1) >> 1) + 1); } return 0; } 

Use >> 2 to divide by 4, >> 3 to divide by 8, &ct. Efficient lg does the work there.

You can even divide by 1! >> 0