[C/ASM] Optimiser au cycle près : la référence
Posté le 20/06/2021 22:24
J'adore programmer des add-ins, je pense que ça ne vous aura pas échappé. Le sujet est très large, mais s'il y a une problématique spécifique que je considère ma «spécialité», c'est d'exploiter toutes les astuces logicielles et matérielles pour améliorer les performances.
Ce topic est une référence de techniques d'optimisation d'add-ins, de tous les styles et niveaux. Ce post contient un index de tous les sujets qui me paraissent utiles (sauf ceux que j'ai oubliés), et j'écrirai des bouts au gré des opportunités et des découvertes.
Le sujet étant large, je ne m'attarderai pas sur les notions générales de C/assembleur dans un premier temps. Cela dit, seules les catégories «Bas niveau» risquent d'être un peu rudes ; en cas de doute, demandez dans les commentaires
Tutoriels qui vous aideront à comprendre les concepts :
La référence de l'optimisation au cycle près
Identifier les besoins d'optimisation
Il est très tentant, mais inutile, d'optimiser du code qui ne limite pas les performances du programme. Identifier les parties critiques et les liens faibles qui monopolisent les ressources est toujours la première étape.
Optimisations algorithmiques
On ne peut pas sous-estimer l'importance des optimisations algorithmiques. En un sens, l'algorithmique est la science de calculer ce qu'on veut avec la bonne
méthode, indépendamment de l'implémentation. Avec l'exception du dessin sur Graph 90+E et des programmes très intensifs en calcul, l'immense majorité des améliorations du
code sont des fourmis au pied de la montagne de l'algorithmique.
Optimisations de calcul et d'implémentation
Vous allez voir que nos chères calculatrices ne sont pas si bien équipées pour calculer que ça.
Optimisations sur la mémoire
En général, plus vous progressez en C, plus vous regardez la mémoire de près. Je vais jusqu'à dire que la mémoire est
le concept fondamental et omniprésent qui guide l'implémentation de tout programme en C (avec bien sûr l'algorithmique qui guide la conception). Sur la calculatrice, c'est encore plus vrai, et de toutes les optimisations de code très peu concurrenceront une gestion supérieure de la mémoire.
Optimisations sur le code assembleur
C'est dans le boucles critiques, exécutées plusieurs millions de fois par seconde, où chaque cycle compte, que ce topic trouve son nom. Vous avez optimisé les algos, les méthodes de calcul, conçu le programme avec la meilleure distribution possible de la mémoire, et maintenant vous voulez écrire en assembleur la version la plus rapide possible. Voilà comment.
Citer : Posté le 21/06/2021 22:17 | #
Calcul en point fixe et en point flottant (Programmation générale — ★☆☆)
Format des types en point flottant
En plus des types entiers (int et compagnie), le langage C spécifie deux types de nombres «à virgule» : float et double. Le premier est nommé comme ça parce que c'est un nombre à «virgule flottante», et le second parce que c'est un float avec double précision.
Quoi qu'on fasse, pour un ordinateur un nombre c'est une suite de chiffres, du style 37492. Si le nombre était entier ça s'arrêterait là, mais pour un nombre décimal on veut pouvoir indiquer la position de la virgule. Le principe du nombre en point flottant c'est de stocker, en plus de la suite de chiffres, la position de la virgule (sous la forme d'une puissance de 2, ou pour l'exemple ci-dessous, de 10).
Par exemple (conceptuellement), le nombre en point flottant «valeur=37492, puissance=10²» va représenter le nombre décimal 37492×10² = 3749200. Le nombre «valeur=37492, puissance=10⁻³» représenterait 37492×10⁻³ = 37.492 qui a 3 chiffres après la virgule. Généralement on peut mettre des puissances 10 assez fines, imaginez-vous 10⁻³⁸ pour une float ou 10⁻³⁰⁰ pour un double, et pareil dans les puissances positives. La plage est immense, on peut vraiment tout représenter.
À cela, on ajoute un bit pour spécifier si le nombre est positif ou négatif (le complément à 2 ne marchant pas vraiment), et donc on obtient un nombre qui a 3 champs :
Tout ça est compacté ensemble dans 32 bits (pour un float) ou 64 bits (pour un double) selon le format fixé par la norme IEEE 754. (Il existait d'autres formats avant, et c'était un enfer, mais aujourd'hui c'est fini et il n'existe plus que l'IEEE 754, ce qui est vraiment vraiment bien.)
La complexité de calculer en point flottant
Imaginons que vous voulez additionner deux nombres. Pour des entiers, c'est vraiment facile : il suffit de les jeter dans un circuit d'addition, qui fait la même chose qu'un humain qui pose l'addition au papier. Les entiers négatifs sont encodés d'une façon particulière (complément à 2) conçue pour qu'ajouter un nombre négatif ne demande aucune action particulière, le circuit est identique. En assembleur, ça prend 1 seule instruction qui occupe l'ALU (unité de calcul arithmétique) pendant 1 seul cycle d'horloge, on ne peut pas faire plus simple.
Pour des nombres en point flottant, c'est beaucoup plus compliqué.
Vous ne serez donc pas surpris quand je dirai qu'électroniquement, c'est beaucoup plus compliqué, et ça ne prend pas qu'1 cycle. C'est assez compliqué pour qu'un coprocesseur complet y soit dédié d'ailleurs ; on appelle ça un FPU (pour Floating-Point Unit ou unité de calcul en point flottant).
Grâce à des techniques électroniques (pipelining), un processeur moderne avec un FPU à qui on donne un gros paquet d'additions à faire va réussir à en compléter une par cycle en moyenne (il triche, il en exécute plusieurs en même temps). Donc on arrive à s'en sortir, mais c'est plus difficile de réaliser cette moyenne qu'avec des entiers.
Point flottant logiciel sur la calculatrice
Sur la calculatrice, vous serez peut-être surpris de savoir qu'il n'y a pas de FPU. Le processeur est un SuperH 4AL-DSP, qui est variante du SuperH 4A où le FPU est remplacé par un DSP (un autre type de coprocesseur plutôt orienté vers le traitement du signal). Pour une calculatrice, ça peut paraître surprenant étant donné la grande quantité de calcul, mais c'est surtout parce que CASIO préfère utiliser des nombres en point flottant en base 10, pour éviter les problèmes d'approximation de la base 2 (c'est un sujet à part entière, je ne détaille pas). Notez juste qu'en Basic CASIO toutes les variables sont des nombres en point flottant en base 10, sauf en mode BASE où c'est des entiers en binaire.
Je pense que vous voyez venir la tragédie. Toute la complexité des nombres en point flottant ne peut pas être cachée dans le circuit électronique d'un FPU s'il n'y a pas de FPU. Et donc tout doit être fait en logiciel, et on parle de point flottant logiciel (ou soft-fp pour les anglophones).
Voilà le code assembleur pour ajouter deux entiers (stockés dans les registres r4 et r5).
Et voilà le code assembleur pour ajouter deux float (pris dans _addsub_fs.o de libgcc ; vous pouvez ignorer les nombres en début de ligne).
0: 61 42 mov.l @r4,r1
2: e3 01 mov #1,r3
4: 31 36 cmp/hi r3,r1
6: 89 01 bt c <__fpadd_parts+0xc>
8: a0 85 bra 116 <__fpadd_parts+0x116>
a: 60 43 mov r4,r0
c: 62 52 mov.l @r5,r2
e: 32 36 cmp/hi r3,r2
10: 8b 52 bf b8 <__fpadd_parts+0xb8>
12: 60 13 mov r1,r0
14: 88 04 cmp/eq #4,r0
16: 8f 02 bf.s 1e <__fpadd_parts+0x1e>
18: 60 23 mov r2,r0
1a: a0 80 bra 11e <__fpadd_parts+0x11e>
1c: 88 04 cmp/eq #4,r0
1e: 88 04 cmp/eq #4,r0
20: 8d 4a bt.s b8 <__fpadd_parts+0xb8>
22: 88 02 cmp/eq #2,r0
24: 8d 43 bt.s ae <__fpadd_parts+0xae>
26: 60 13 mov r1,r0
28: 88 02 cmp/eq #2,r0
2a: 89 45 bt b8 <__fpadd_parts+0xb8>
2c: 2f 86 mov.l r8,@-r15
2e: 2f 96 mov.l r9,@-r15
30: 52 42 mov.l @(8,r4),r2
32: 57 52 mov.l @(8,r5),r7
34: 61 23 mov r2,r1
36: 50 43 mov.l @(12,r4),r0
38: 31 78 sub r7,r1
3a: 41 11 cmp/pz r1
3c: 8f 4d bf.s da <__fpadd_parts+0xda>
3e: 53 53 mov.l @(12,r5),r3
40: e8 1f mov #31,r8
42: 31 87 cmp/gt r8,r1
44: 8d 4f bt.s e6 <__fpadd_parts+0xe6>
46: 21 18 tst r1,r1
48: 8d 0a bt.s 60 <__fpadd_parts+0x60>
4a: 67 1b neg r1,r7
4c: 69 33 mov r3,r9
4e: 49 7d shld r7,r9
50: e7 ff mov #-1,r7
52: 68 73 mov r7,r8
54: 48 1d shld r1,r8
56: 61 87 not r8,r1
58: 21 38 tst r3,r1
5a: 67 7a negc r7,r7
5c: 63 93 mov r9,r3
5e: 23 7b or r7,r3
60: 51 41 mov.l @(4,r4),r1
62: 57 51 mov.l @(4,r5),r7
64: 31 70 cmp/eq r7,r1
66: 8d 43 bt.s f0 <__fpadd_parts+0xf0>
68: 21 18 tst r1,r1
6a: 61 03 mov r0,r1
6c: 8d 02 bt.s 74 <__fpadd_parts+0x74>
6e: 31 38 sub r3,r1
70: 61 33 mov r3,r1
72: 31 08 sub r0,r1
74: 41 11 cmp/pz r1
76: 8d 02 bt.s 7e <__fpadd_parts+0x7e>
78: e3 00 mov #0,r3
7a: 61 1b neg r1,r1
7c: e3 01 mov #1,r3
7e: d7 32 mov.l 148 <__fpadd_parts+0x148>,r7 ! 3ffffffe
80: 16 31 mov.l r3,@(4,r6)
82: 63 13 mov r1,r3
84: 73 ff add #-1,r3
86: 16 22 mov.l r2,@(8,r6)
88: 33 76 cmp/hi r7,r3
8a: 8d 0a bt.s a2 <__fpadd_parts+0xa2>
8c: 16 13 mov.l r1,@(12,r6)
8e: 72 ff add #-1,r2
90: 31 1c add r1,r1
92: 63 13 mov r1,r3
94: 73 ff add #-1,r3
96: 33 76 cmp/hi r7,r3
98: 65 23 mov r2,r5
9a: 8f f9 bf.s 90 <__fpadd_parts+0x90>
9c: 72 ff add #-1,r2
9e: 16 13 mov.l r1,@(12,r6)
a0: 16 52 mov.l r5,@(8,r6)
a2: e1 03 mov #3,r1
a4: 26 12 mov.l r1,@r6
a6: 60 63 mov r6,r0
a8: 69 f6 mov.l @r15+,r9
aa: 00 0b rts
ac: 68 f6 mov.l @r15+,r8
ae: 88 02 cmp/eq #2,r0
b0: 8d 04 bt.s bc <__fpadd_parts+0xbc>
b2: 60 43 mov r4,r0
b4: 00 0b rts
b6: 00 09 nop
b8: 00 0b rts
ba: 60 53 mov r5,r0
bc: 26 12 mov.l r1,@r6
be: 60 63 mov r6,r0
c0: 51 41 mov.l @(4,r4),r1
c2: 16 11 mov.l r1,@(4,r6)
c4: 61 43 mov r4,r1
c6: 71 08 add #8,r1
c8: 62 16 mov.l @r1+,r2
ca: 16 22 mov.l r2,@(8,r6)
cc: 61 12 mov.l @r1,r1
ce: 52 51 mov.l @(4,r5),r2
d0: 16 13 mov.l r1,@(12,r6)
d2: 51 41 mov.l @(4,r4),r1
d4: 21 29 and r2,r1
d6: 00 0b rts
d8: 16 11 mov.l r1,@(4,r6)
da: 61 73 mov r7,r1
dc: 31 28 sub r2,r1
de: e8 1f mov #31,r8
e0: 31 87 cmp/gt r8,r1
e2: 8f 25 bf.s 130 <__fpadd_parts+0x130>
e4: 69 03 mov r0,r9
e6: 32 77 cmp/gt r7,r2
e8: 89 17 bt 11a <__fpadd_parts+0x11a>
ea: 62 73 mov r7,r2
ec: af b8 bra 60 <__fpadd_parts+0x60>
ee: e0 00 mov #0,r0
f0: 30 3c add r3,r0
f2: 16 11 mov.l r1,@(4,r6)
f4: e1 03 mov #3,r1
f6: 16 22 mov.l r2,@(8,r6)
f8: 40 11 cmp/pz r0
fa: 16 03 mov.l r0,@(12,r6)
fc: 8d d3 bt.s a6 <__fpadd_parts+0xa6>
fe: 26 12 mov.l r1,@r6
100: 61 03 mov r0,r1
102: 41 01 shlr r1
104: c9 01 and #1,r0
106: 21 0b or r0,r1
108: 69 f6 mov.l @r15+,r9
10a: 72 01 add #1,r2
10c: 60 63 mov r6,r0
10e: 68 f6 mov.l @r15+,r8
110: 16 13 mov.l r1,@(12,r6)
112: 00 0b rts
114: 16 22 mov.l r2,@(8,r6)
116: 00 0b rts
118: 00 09 nop
11a: af a1 bra 60 <__fpadd_parts+0x60>
11c: e3 00 mov #0,r3
11e: 8f c9 bf.s b4 <__fpadd_parts+0xb4>
120: 60 43 mov r4,r0
122: 52 41 mov.l @(4,r4),r2
124: 51 51 mov.l @(4,r5),r1
126: 32 10 cmp/eq r1,r2
128: 89 c4 bt b4 <__fpadd_parts+0xb4>
12a: d0 08 mov.l 14c <__fpadd_parts+0x14c>,r0 ! 0 <__fpadd_parts>
12c: 00 0b rts
12e: 00 09 nop
130: 62 1b neg r1,r2
132: 49 2d shld r2,r9
134: e2 ff mov #-1,r2
136: 68 23 mov r2,r8
138: 48 1d shld r1,r8
13a: 61 87 not r8,r1
13c: 21 08 tst r0,r1
13e: 62 2a negc r2,r2
140: 60 93 mov r9,r0
142: 20 2b or r2,r0
144: af 8c bra 60 <__fpadd_parts+0x60>
146: 62 73 mov r7,r2
148: 3f ff addv r15,r15
14a: ff fe .word 0xfffe
14c: 00 00 .word 0x0000
...
00000150 <___addsf3>:
150: 2f 86 mov.l r8,@-r15
152: 4f 22 sts.l pr,@-r15
154: d8 0e mov.l 190 <___addsf3+0x40>,r8 ! 0 <__fpadd_parts>
156: 7f c8 add #-56,r15
158: 1f 51 mov.l r5,@(4,r15)
15a: 65 f3 mov r15,r5
15c: 2f 42 mov.l r4,@r15
15e: 75 08 add #8,r5
160: 48 0b jsr @r8
162: 64 f3 mov r15,r4
164: 65 f3 mov r15,r5
166: 64 f3 mov r15,r4
168: 75 18 add #24,r5
16a: 48 0b jsr @r8
16c: 74 04 add #4,r4
16e: d0 09 mov.l 194 <___addsf3+0x44>,r0 ! 0 <__fpadd_parts>
170: 66 f3 mov r15,r6
172: 65 f3 mov r15,r5
174: 64 f3 mov r15,r4
176: 76 28 add #40,r6
178: 75 18 add #24,r5
17a: 40 0b jsr @r0
17c: 74 08 add #8,r4
17e: 64 03 mov r0,r4
180: d0 05 mov.l 198 <___addsf3+0x48>,r0 ! 0 <__fpadd_parts>
182: 40 0b jsr @r0
184: 00 09 nop
186: 7f 38 add #56,r15
188: 4f 26 lds.l @r15+,pr
18a: 00 0b rts
18c: 68 f6 mov.l @r15+,r8
18e: 00 09 nop
Vous voyez la complexité ? Non ? Voilà le code pour ajouter deux double (pris dans _addsub_df.o de libgcc).
0: 61 42 mov.l @r4,r1
2: e0 01 mov #1,r0
4: 31 06 cmp/hi r0,r1
6: 89 01 bt c <__fpadd_parts+0xc>
8: a0 d8 bra 1bc <__fpadd_parts+0x1bc>
a: 60 43 mov r4,r0
c: 67 52 mov.l @r5,r7
e: 37 06 cmp/hi r0,r7
10: 8b 6d bf ee <__fpadd_parts+0xee>
12: 60 13 mov r1,r0
14: 88 04 cmp/eq #4,r0
16: 8f 02 bf.s 1e <__fpadd_parts+0x1e>
18: 60 73 mov r7,r0
1a: a0 d7 bra 1cc <__fpadd_parts+0x1cc>
1c: 88 04 cmp/eq #4,r0
1e: 88 04 cmp/eq #4,r0
20: 8d 65 bt.s ee <__fpadd_parts+0xee>
22: 88 02 cmp/eq #2,r0
24: 8d 5e bt.s e4 <__fpadd_parts+0xe4>
26: 60 13 mov r1,r0
28: 88 02 cmp/eq #2,r0
2a: 89 60 bt ee <__fpadd_parts+0xee>
2c: 2f 86 mov.l r8,@-r15
2e: 2f 96 mov.l r9,@-r15
30: 2f a6 mov.l r10,@-r15
32: 2f b6 mov.l r11,@-r15
34: 2f c6 mov.l r12,@-r15
36: 2f d6 mov.l r13,@-r15
38: 2f e6 mov.l r14,@-r15
3a: 51 42 mov.l @(8,r4),r1
3c: 5b 52 mov.l @(8,r5),r11
3e: 6a 13 mov r1,r10
40: 59 43 mov.l @(12,r4),r9
42: 3a b8 sub r11,r10
44: 4a 11 cmp/pz r10
46: 58 44 mov.l @(16,r4),r8
48: 50 53 mov.l @(12,r5),r0
4a: 8f 63 bf.s 114 <__fpadd_parts+0x114>
4c: 57 54 mov.l @(16,r5),r7
4e: ec 3f mov #63,r12
50: 3a c7 cmp/gt r12,r10
52: 8d 65 bt.s 120 <__fpadd_parts+0x120>
54: 2a a8 tst r10,r10
56: 8d 68 bt.s 12a <__fpadd_parts+0x12a>
58: 6e a3 mov r10,r14
5a: 7e e0 add #-32,r14
5c: 4e 11 cmp/pz r14
5e: 89 01 bt 64 <__fpadd_parts+0x64>
60: a0 e2 bra 228 <__fpadd_parts+0x228>
62: 6b ab neg r10,r11
64: 6c eb neg r14,r12
66: 6b 03 mov r0,r11
68: 4b cd shld r12,r11
6a: ec 00 mov #0,r12
6c: 4e 11 cmp/pz r14
6e: 89 01 bt 74 <__fpadd_parts+0x74>
70: a0 d1 bra 216 <__fpadd_parts+0x216>
72: 6d a3 mov r10,r13
74: ed ff mov #-1,r13
76: 4d ed shld r14,r13
78: ee 00 mov #0,r14
7a: 6d d7 not r13,r13
7c: 6e e7 not r14,r14
7e: 2d 09 and r0,r13
80: 2e 79 and r7,r14
82: 54 41 mov.l @(4,r4),r4
84: 2d eb or r14,r13
86: 55 51 mov.l @(4,r5),r5
88: 2d d8 tst r13,r13
8a: e7 ff mov #-1,r7
8c: 67 7a negc r7,r7
8e: 34 50 cmp/eq r5,r4
90: 60 c3 mov r12,r0
92: 8f 4f bf.s 134 <__fpadd_parts+0x134>
94: 27 bb or r11,r7
96: 00 08 clrt
98: 63 83 mov r8,r3
9a: 33 7e addc r7,r3
9c: 62 93 mov r9,r2
9e: 32 0e addc r0,r2
a0: 16 41 mov.l r4,@(4,r6)
a2: 16 12 mov.l r1,@(8,r6)
a4: 16 23 mov.l r2,@(12,r6)
a6: 16 34 mov.l r3,@(16,r6)
a8: e1 03 mov #3,r1
aa: 26 12 mov.l r1,@r6
ac: d1 73 mov.l 27c <__fpadd_parts+0x27c>,r1 ! 1fffffff
ae: 32 12 cmp/hs r1,r2
b0: 8b 01 bf b6 <__fpadd_parts+0xb6>
b2: 32 16 cmp/hi r1,r2
b4: 89 01 bt ba <__fpadd_parts+0xba>
b6: a0 92 bra 1de <__fpadd_parts+0x1de>
b8: 60 63 mov r6,r0
ba: 6e f6 mov.l @r15+,r14
bc: 64 23 mov r2,r4
be: 6d f6 mov.l @r15+,r13
c0: 44 01 shlr r4
c2: 6c f6 mov.l @r15+,r12
c4: 65 33 mov r3,r5
c6: 51 62 mov.l @(8,r6),r1
c8: 60 33 mov r3,r0
ca: 6b f6 mov.l @r15+,r11
cc: 45 25 rotcr r5
ce: 6a f6 mov.l @r15+,r10
d0: c9 01 and #1,r0
d2: 20 5b or r5,r0
d4: 69 f6 mov.l @r15+,r9
d6: 71 01 add #1,r1
d8: 16 04 mov.l r0,@(16,r6)
da: 60 63 mov r6,r0
dc: 68 f6 mov.l @r15+,r8
de: 16 43 mov.l r4,@(12,r6)
e0: 00 0b rts
e2: 16 12 mov.l r1,@(8,r6)
e4: 88 02 cmp/eq #2,r0
e6: 8d 04 bt.s f2 <__fpadd_parts+0xf2>
e8: 60 43 mov r4,r0
ea: 00 0b rts
ec: 00 09 nop
ee: 00 0b rts
f0: 60 53 mov r5,r0
f2: 26 12 mov.l r1,@r6
f4: 60 63 mov r6,r0
f6: 51 41 mov.l @(4,r4),r1
f8: 16 11 mov.l r1,@(4,r6)
fa: 61 43 mov r4,r1
fc: 71 08 add #8,r1
fe: 62 16 mov.l @r1+,r2
100: 16 22 mov.l r2,@(8,r6)
102: 62 16 mov.l @r1+,r2
104: 16 23 mov.l r2,@(12,r6)
106: 61 12 mov.l @r1,r1
108: 52 51 mov.l @(4,r5),r2
10a: 16 14 mov.l r1,@(16,r6)
10c: 51 41 mov.l @(4,r4),r1
10e: 21 29 and r2,r1
110: 00 0b rts
112: 16 11 mov.l r1,@(4,r6)
114: 6a b3 mov r11,r10
116: 3a 18 sub r1,r10
118: ec 3f mov #63,r12
11a: 3a c7 cmp/gt r12,r10
11c: 8f 61 bf.s 1e2 <__fpadd_parts+0x1e2>
11e: 6e a3 mov r10,r14
120: 31 b7 cmp/gt r11,r1
122: 89 50 bt 1c6 <__fpadd_parts+0x1c6>
124: 61 b3 mov r11,r1
126: e9 00 mov #0,r9
128: e8 00 mov #0,r8
12a: 54 41 mov.l @(4,r4),r4
12c: 55 51 mov.l @(4,r5),r5
12e: 34 50 cmp/eq r5,r4
130: 8d b2 bt.s 98 <__fpadd_parts+0x98>
132: 00 08 clrt
134: 24 48 tst r4,r4
136: 8d 3c bt.s 1b2 <__fpadd_parts+0x1b2>
138: 65 83 mov r8,r5
13a: 65 73 mov r7,r5
13c: 35 8a subc r8,r5
13e: 67 03 mov r0,r7
140: 37 9a subc r9,r7
142: 47 11 cmp/pz r7
144: 8f 3c bf.s 1c0 <__fpadd_parts+0x1c0>
146: 63 5a negc r5,r3
148: 62 73 mov r7,r2
14a: 63 53 mov r5,r3
14c: e7 00 mov #0,r7
14e: d4 4c mov.l 280 <__fpadd_parts+0x280>,r4 ! fffffff
150: e5 ff mov #-1,r5
152: 00 08 clrt
154: 16 71 mov.l r7,@(4,r6)
156: 35 3e addc r3,r5
158: e7 ff mov #-1,r7
15a: 37 2e addc r2,r7
15c: 16 12 mov.l r1,@(8,r6)
15e: 37 46 cmp/hi r4,r7
160: 16 23 mov.l r2,@(12,r6)
162: 8d a1 bt.s a8 <__fpadd_parts+0xa8>
164: 16 34 mov.l r3,@(16,r6)
166: 37 42 cmp/hs r4,r7
168: 8b 01 bf 16e <__fpadd_parts+0x16e>
16a: a0 81 bra 270 <__fpadd_parts+0x270>
16c: e7 fe mov #-2,r7
16e: d0 44 mov.l 280 <__fpadd_parts+0x280>,r0 ! fffffff
170: 71 ff add #-1,r1
172: e4 ff mov #-1,r4
174: e5 ff mov #-1,r5
176: ea fe mov #-2,r10
178: 43 00 shll r3
17a: 00 09 nop
17c: 42 24 rotcl r2
17e: 00 08 clrt
180: 69 33 mov r3,r9
182: 39 5e addc r5,r9
184: 67 23 mov r2,r7
186: 37 4e addc r4,r7
188: 37 02 cmp/hs r0,r7
18a: 68 13 mov r1,r8
18c: 8f f4 bf.s 178 <__fpadd_parts+0x178>
18e: 71 ff add #-1,r1
190: 37 06 cmp/hi r0,r7
192: 8f 68 bf.s 266 <__fpadd_parts+0x266>
194: 39 a6 cmp/hi r10,r9
196: e1 03 mov #3,r1
198: 16 23 mov.l r2,@(12,r6)
19a: 60 63 mov r6,r0
19c: 16 34 mov.l r3,@(16,r6)
19e: 16 82 mov.l r8,@(8,r6)
1a0: 26 12 mov.l r1,@r6
1a2: 6e f6 mov.l @r15+,r14
1a4: 6d f6 mov.l @r15+,r13
1a6: 6c f6 mov.l @r15+,r12
1a8: 6b f6 mov.l @r15+,r11
1aa: 6a f6 mov.l @r15+,r10
1ac: 69 f6 mov.l @r15+,r9
1ae: 00 0b rts
1b0: 68 f6 mov.l @r15+,r8
1b2: 00 08 clrt
1b4: 35 7a subc r7,r5
1b6: 67 93 mov r9,r7
1b8: af c3 bra 142 <__fpadd_parts+0x142>
1ba: 37 0a subc r0,r7
1bc: 00 0b rts
1be: 00 09 nop
1c0: 62 7a negc r7,r2
1c2: af c4 bra 14e <__fpadd_parts+0x14e>
1c4: e7 01 mov #1,r7
1c6: e0 00 mov #0,r0
1c8: af af bra 12a <__fpadd_parts+0x12a>
1ca: e7 00 mov #0,r7
1cc: 8f 8d bf.s ea <__fpadd_parts+0xea>
1ce: 60 43 mov r4,r0
1d0: 52 41 mov.l @(4,r4),r2
1d2: 51 51 mov.l @(4,r5),r1
1d4: 32 10 cmp/eq r1,r2
1d6: 89 88 bt ea <__fpadd_parts+0xea>
1d8: d0 2a mov.l 284 <__fpadd_parts+0x284>,r0 ! 0 <__fpadd_parts>
1da: 00 0b rts
1dc: 00 09 nop
1de: af e0 bra 1a2 <__fpadd_parts+0x1a2>
1e0: 00 09 nop
1e2: 7e e0 add #-32,r14
1e4: 4e 11 cmp/pz r14
1e6: 8f 2a bf.s 23e <__fpadd_parts+0x23e>
1e8: 61 ab neg r10,r1
1ea: 6c eb neg r14,r12
1ec: 61 93 mov r9,r1
1ee: 41 cd shld r12,r1
1f0: ec 00 mov #0,r12
1f2: 4e 11 cmp/pz r14
1f4: 8f 2e bf.s 254 <__fpadd_parts+0x254>
1f6: 6d a3 mov r10,r13
1f8: ed ff mov #-1,r13
1fa: 4d ed shld r14,r13
1fc: ee 00 mov #0,r14
1fe: 6d d7 not r13,r13
200: 6e e7 not r14,r14
202: 2d 99 and r9,r13
204: 2e 89 and r8,r14
206: 2d eb or r14,r13
208: 2d d8 tst r13,r13
20a: e8 ff mov #-1,r8
20c: 68 8a negc r8,r8
20e: 28 1b or r1,r8
210: 69 c3 mov r12,r9
212: af 8a bra 12a <__fpadd_parts+0x12a>
214: 61 b3 mov r11,r1
216: de 1c mov.l 288 <__fpadd_parts+0x288>,r14 ! 7fffffff
218: 7d e1 add #-31,r13
21a: 4e dd shld r13,r14
21c: ed ff mov #-1,r13
21e: 4d ad shld r10,r13
220: 2d eb or r14,r13
222: ee ff mov #-1,r14
224: af 29 bra 7a <__fpadd_parts+0x7a>
226: 4e ad shld r10,r14
228: 6c 03 mov r0,r12
22a: 7b 1f add #31,r11
22c: 3c cc add r12,r12
22e: 4c bd shld r11,r12
230: 6d ab neg r10,r13
232: 6b 73 mov r7,r11
234: 4b dd shld r13,r11
236: 2b cb or r12,r11
238: 6c 03 mov r0,r12
23a: af 17 bra 6c <__fpadd_parts+0x6c>
23c: 4c dd shld r13,r12
23e: 6c 93 mov r9,r12
240: 71 1f add #31,r1
242: 3c cc add r12,r12
244: 4c 1d shld r1,r12
246: 6d ab neg r10,r13
248: 61 83 mov r8,r1
24a: 41 dd shld r13,r1
24c: 21 cb or r12,r1
24e: 6c 93 mov r9,r12
250: af cf bra 1f2 <__fpadd_parts+0x1f2>
252: 4c dd shld r13,r12
254: de 0c mov.l 288 <__fpadd_parts+0x288>,r14 ! 7fffffff
256: 7d e1 add #-31,r13
258: 4e dd shld r13,r14
25a: ed ff mov #-1,r13
25c: 4d ad shld r10,r13
25e: 2d eb or r14,r13
260: ee ff mov #-1,r14
262: af cc bra 1fe <__fpadd_parts+0x1fe>
264: 4e ad shld r10,r14
266: 89 01 bt 26c <__fpadd_parts+0x26c>
268: af 88 bra 17c <__fpadd_parts+0x17c>
26a: 43 00 shll r3
26c: af 94 bra 198 <__fpadd_parts+0x198>
26e: e1 03 mov #3,r1
270: 35 76 cmp/hi r7,r5
272: 8b 01 bf 278 <__fpadd_parts+0x278>
274: af 19 bra aa <__fpadd_parts+0xaa>
276: e1 03 mov #3,r1
278: af 79 bra 16e <__fpadd_parts+0x16e>
27a: 00 09 nop
27c: 1f ff mov.l r15,@(60,r15)
27e: ff ff .word 0xffff
280: 0f ff mac.l @r15+,@r15+
282: ff ff .word 0xffff
284: 00 00 .word 0x0000
286: 00 00 .word 0x0000
288: 7f ff add #-1,r15
28a: ff ff .word 0xffff
0000028c <___adddf3>:
28c: 2f 86 mov.l r8,@-r15
28e: 4f 22 sts.l pr,@-r15
290: d8 0f mov.l 2d0 <___adddf3+0x44>,r8 ! 0 <__fpadd_parts>
292: 7f b4 add #-76,r15
294: 1f 51 mov.l r5,@(4,r15)
296: 65 f3 mov r15,r5
298: 1f 73 mov.l r7,@(12,r15)
29a: 75 10 add #16,r5
29c: 1f 62 mov.l r6,@(8,r15)
29e: 2f 42 mov.l r4,@r15
2a0: 48 0b jsr @r8
2a2: 64 f3 mov r15,r4
2a4: 65 f3 mov r15,r5
2a6: 64 f3 mov r15,r4
2a8: 75 24 add #36,r5
2aa: 48 0b jsr @r8
2ac: 74 08 add #8,r4
2ae: d0 09 mov.l 2d4 <___adddf3+0x48>,r0 ! 0 <__fpadd_parts>
2b0: 66 f3 mov r15,r6
2b2: 65 f3 mov r15,r5
2b4: 64 f3 mov r15,r4
2b6: 76 38 add #56,r6
2b8: 75 24 add #36,r5
2ba: 40 0b jsr @r0
2bc: 74 10 add #16,r4
2be: 64 03 mov r0,r4
2c0: d0 05 mov.l 2d8 <___adddf3+0x4c>,r0 ! 0 <__fpadd_parts>
2c2: 40 0b jsr @r0
2c4: 00 09 nop
2c6: 7f 4c add #76,r15
2c8: 4f 26 lds.l @r15+,pr
2ca: 00 0b rts
2cc: 68 f6 mov.l @r15+,r8
2ce: 00 09 nop
Ajouter deux entiers, ici, prend environ 8 ns (si on ignore le parallélisme bas-niveau qui est un sujet beaucoup plus compliqué). Ajouter deux float quelconques m'a pris autour de 2-3 µs (~300 fois plus), et ajouter deux double autour de 3-4 µs (~400 fois plus).
Et là ce n'est qu'une addition, ça ne s'améliore pas quand il s'agit de multiplier (ordre d'idée en double : dans les 5 µs), et vous pouvez oublier la division si vous tenez aux performances (ordre d'idée en double : dans les 10-11 µs). Diviser est lent aussi pour les entiers (et donne lieu à d'autres optimisations), mais ça n'a rien à voir avec ce niveau-là.
Soit dit en passant, en plus de prendre du temps ce code prend aussi de la place, et la lib mathématique en point flottant logiciel (OpenLibm pour nous) prend aussi beaucoup de place. Chaque fonction occupe vite 4 kio de code sinon plus, donc les add-ins grossissent vite (et puis ça ruine un peu le cache).
Les nombres en point fixe
Les nombres en point fixe sont une autre option pour représenter des nombres décimaux. Comme son nom l'indique, la virgule est fixée à une certaine position, ce qui élimine l'exposant. Par exemple, pour un point fixe «de puissance 10⁻³», on stockera la valeur 2874 pour représenter le nombre décimal 2.874. La valeur de l'exposant, -3, n'est pas représentée dans la mémoire mais est implicite.
L'avantage c'est qu'éliminer l'exposant signifie qu'on perd un champ sur les 3, et du coup on peut récupérer le complément à 2 et éliminer le bit de signe. Ultimement, le format complexe à trois champs redevient un entier.
Les nombres en point fixe sont vraiment aussi simples qu'écrire 1383 centimes dans un entier pour signifier 13.83 €. Il n'y a pas d'arnaque, on décide simplement qu'on utilise la valeur entière X pour représenter le nombre décimal X/100 (ou X/1000, ou autre puissance de 10), et ça marche exactement pareil avec les nombres négatifs.
Notez quand même que les puissances de 10 ça ne marche pas très bien sur un processeur qui compte en binaire, on préfère utiliser des puissances de 2. Du coup, dans la suite de cette partie, je vais utiliser du binaire.
Les entiers du processeur SuperH font 32 bits, et on a le choix sur la position de la virgule. Ici je vais prendre peut-ête le format le plus évident, où je mets la virgule au milieu, ce qui me laisse 16 chiffres pour la partie entière et 16 pour la partie décimale. On dit que j'utilise le «format de point fixe 16:16».
Comme il y a 16 bits de partie décimale, l'entier x représentera la valeur décimale x/2¹⁶, c'est-à-dire x/65536. Les entiers sur 32 bits vont de -2147483648 à 2147483647, et donc le format en point fixe va aller de -2147483648/2¹⁶ = -32768, à 2147483647/2¹⁶ = 32767.99998.
Comme on va le voir, les opérations sont pas tout à fait les mêmes que les entiers, mais ça se passe assez bien.
Opérations en point fixe
Pour faire les opérations en point fixe sans se planter, je vous conseille de poser les étapes suivantes :
Par exemple, pour additionner. Si j'ai deux nombre décimaux x et y, alors :
La multiplication montre un exemple où il ne suffit pas d'utiliser l'opération équivalente sur les entiers :
On remarque aussi que la multiplication de 2¹⁶x et 2¹⁶y est un nombre entier de 64 bits dont on récupère les 32 bits du milieu, donc il faut bien faire la multiplication 64 bits. On écrit donc, en réalité, z = ((int64_t)x * (int64_t)y) >> 16.
(Le compilateur est assez intelligent pour faire une multiplication 32×32→64 bits au lieu de 64×64→64 bits, et utilise même l'instruction SuperH xtrct pour sortir les 32 bits du milieu.)
Implémentation en pratique
Je ne voudrais pas trop m'étendre sur le code sinon il y en a pour des pages, mais voici les grandes idées. Je vous conseille de donner un nom différent au type des nombres en point fixe, parce que si vous utilisez int vous allez tout le temps confondre les vrais entiers avec les nombres en point fixe, et croyez-moi vous n'avez pas envie de debugger ça.
Par exemple, on peut appeler ce type fixed.
Notez que le compilateur ne vous empêchera pas d'additioner un int et un fixed. Le nommage ne sert qu'à vous pour bien vous souvenir de ce que représente chaque variable. En C++, si vous enrobez l'entier dans une classe, vous pourrez empêcher les interactions de ce genre, et vous ne le regretterez pas ; je vous le conseille. Vous pourriez enrober l'entier dans une structure en C et obtenir le même effet, mais ça rendrait le code beaucoup plus verbeux, donc c'est pas forcément intéressant.
Pour obtenir un fixed à partir d'un entier ou d'un nombre décimal, multipliez-le simplement par 2¹⁶ = 65536.
fixed y = 1.625 * 65536; // y = 1.625 en point fixe 16:16
Vous pourriez être tentés de faire << 16 à la place de la multiplication, mais ça donne des warnings (en C) ou des erreurs (en C++) pour les nombres négatifs. C'est parce que le C/C++ n'impose pas que le complément à 2 soit utilisé pour les entiers, et si le complément à 2 n'est pas utilisé, un décalage de bits d'un nombre négatif vers la gauche serait illégal/insensé. Ne vous embêtez pas, écrivez *65536, le compilateur optimisera ça en un décalage de bits absolument toutes les fois.
Notez que le 1.625 * 65536 est calculé par le compilateur, vous n'allez pas vous retrouver avec une multiplication de double pour initialiser votre nombre en point fixe (ce serait fort !).
Pour additioner et multiplier, ça donne les opérations que j'ai mentionnées plus haut ; ça marche aussi pour soustraire, diviser, le négatif...
x + y;
# Soustraction
x - y;
# Opposé
-x;
# Multiplication
((int64_t)x * (int64_t)y) >> 16;
# Division (y ≠ 0)
((int64_t)x * 65536) / y;
# Modulo (y ≠ 0)
x % y;
# Incrémenter
x + (1 * 65536) // et pas x+1
Je vous conseille d'enrober ça dans des fonctions pour ne pas avoir à le réécrire à chaque fois. Dans ce cas, il est de bon goût de les avoir en static inline dans un header, pour que le compilateur les insère sur place au lieu de faire un vrai appel de fonction (qui prendrait aussi longtemps que l'opération en elle-même).
Pendant que j'y suis, je vous invite très fortement à avoir une fonction fix() pour générer les constantes, histoire d'éviter d'écrire x+1 par erreur (ça m'est arrivé, c'est aussi chiant à debugger que ça en a l'air).
{
return n * 65536;
}
static inline fixed mul_fixed(fixed x, fixed y)
{
return ((int64_t)x * (int64_t)y) >> 16;
}
Vous ne pouvez pas écrire x*y à moins de faire du C++, il faut se contenter la notation mul_fixed(x,y) et du boost énorme de performances qui l'accompagne par rapport aux nombres en point flottant.
Avantages, inconvénients, conclusion
Voilà pour les grandes idées du calcul en point fixe. Les avantages sont assez clairs, vous allez avoir un code beaucoup plus compact et beacoup plus rapide, ce qui devient vite très important pour les programmes qui calculent beaucoup ou les jeux qui ont pas mal de physique.
Il y a plusieurs inconvénients qu'il convient de mesurer pour éviter les déconvenues.
Et dans un moment de publicité qui fera rager Yatis, j'observerai qu'on peut implémenter en C++ une classe de point fixe qui élimine les deux derniers points (ie. pas de risque de mélanger avec des entiers, et syntaxe naturelle avec les opérateurs normaux) sans compromis sur l'optimisation du code résultant. En fait, le C++ donne même des outils pour exiger que certains calculs soient faits à la compilation, ce qui est encore plus avantageux que la version C.
J'ai joué avec ça dans une bibliothèque numérique pour mon moteur de jeu en C++ (un truc plus gros avec différentes tailles de point fixes, et de l'arithmétique de vecteurs/matrices/etc), ça existera sûrement de façon réutilisable un jour. Honnêtement c'est assez agréable à utiliser.
Citer : Posté le 22/06/2021 19:17 | #
Les forces et les limites de l'optimisation dans le compilateur (Programmation générale — ★★☆)
Le compilateur optimise déjà le code C/C++ qu'on lui donne à compiler, alors pourquoi est-ce qu'on s'embêterait à le faire nous-même ? La réponse est assez simple : parce qu'on peut faire mieux. Le compilateur n'est pas à l'aise avec toutes les optimisations du code.
Voici quelques idées générales pour vous aider à déterminer si le compilateur fera une optimisation :
n*256 est égal à n<<8 pour tout entier. Il n'y a besoin d'aucun contexte, aucune hypothèse sur la valeur de n, aucune astuce pour que cette transformation soit valide. Le compilateur la fera systématiquement. Par contre, dérouler une boucle pour faire le travail de 2 tours en un seul n'est possible que si le nombre de tours est pair. Selon la façon dont la fonction est écrite, le compilateur peut ne pas être capable de déterminer si c'est le cas (et aussi ne pas avoir envie de dérouler même s'il l'est).
Le compilateur n'aura pas de mal à voir que dans y = condition ? -1 : x + 1, la valeur de x n'est nécessaire que dans une branche, et saura réordonner le code pour ne pas calculer x si la condition est vraie (en supposant que x ne soit pas utilisé ailleurs). Il n'aura pas non plus de difficultés à voir quelles valeurs calculées dans des boucles sont constantes et peuvent être calculées une seule fois avant le début de la boucle. Mais il ne réalisera probablement pas que plusieurs parties indépendantes de plusieurs fonctions différentes sont exécutées sous les mêmes conditions, même si les conditions sont calculées de façon identique. L'optimisation sur plusieurs fonctions (inter-procédurale) est assez exotique déjà (et sur plusieurs fichiers encore plus ; ça s'appelle de la LTO, mais gint ne l'active pas).
En particulier quand il s'agit de faire du calcul sur les entiers et les booléens, vous pouvez faire confiance au compilateur pour utiliser les bonnes astuces. Par exemple, je parie que vous n'avez pas vu venir la technique ci-dessous (que GCC utilise). Par contre, pas de garanties sur l'optimisation du pipeline ou du parallélisme.
and r4, r5
shll r5
# bt/bf/etc
Ce sera le sujet d'à peu près toute la fin de cette partie. Les astuces du jeu d'instruction SuperH, GCC les connaît. Les astuces d'ordonnancement des instructions spécifiques au SuperH 4AL-DSP, probablement pas. GCC n'ira pas utiliser les registres spéciaux de copie du CPU pour remplacer memcpy() vers la XRAM/YRAM par exemple. GCC ne sait pas non plus à quelle vitesse la RAM répond sur la calculatrice, puisque ce n'est pas lié au processeur mais au MPU. De façon générale plus c'est spécifique moins c'est le job du compilateur d'optimiser.
Quitte à caricaturer, moins il y a besoin de connaître le programme pour optimiser, plus le compilateur est susceptible de s'en charger. Si quelqu'un qui ne connaît pas votre programme ni la plateforme spécifiquement lit le code et repère une optimisation, c'est généralement dans le domaine de ce que le compilateur peut faire. S'il faut connaître en profondeur le fonctionnement du code ou le fonctionnement du matériel, alors il n'y a que vous qui avez «l'autorité» pour effectuer les transformations.
Dans les optimisations que vous pouvez faire mais que le compilateur ne fait pas, on trouve le plus souvent :
En gros tout ce qui est dans l'index de ce topic se rattache à un sujet d'optimisation que le compilateur ne fera pas pour vous.
Citer : Posté le 23/06/2021 09:58 | #
Hey ! Super, bravo !
J’ai juste une petite question :
Par exemple, pour un point fixe «de puissance 10⁻³», on stockera la valeur 2874 pour représenter le nombre décimal 2.874. La valeur de l'exposant, -3, n'est pas représentée dans la mémoire mais est implicite.
Pourquoi dans 2874 la valeur -3 est implicite ?
Le nombre pourrait tout autant être 287,4, et à ce moment là ce ne sera pas -3 mais le nombre stocke dans la mémoire sera quand meme 2874…
Citer : Posté le 23/06/2021 10:01 | #
Je devine que Lephé utilise la notation scientifique, avec par défaut un chiffre avant la virgule. Je suis peut être complètement à côté de la plaque.
Citer : Posté le 23/06/2021 10:05 | #
Mais, si tu veux 287,4 tu fais comment vu que c’est un chiffre avant la virgule ?
Citer : Posté le 23/06/2021 10:06 | #
Très bonne question ! C'est justement parce que je choisis d'utiliser -3 comme exposant. Le principe du type en point fixe, c'est de faire comme le point flottant sauf que tous les nombres ont le même exposant. Le programmeur a le choix de cet exposant, mais doit ensuite s'y tenir.
Si j'avais choisi -4, mon entier aurait dû contenir 28740 pour représenter 2.874. De la même façon, si j'avais choisi -1 et utilisé la valeur 287.4, j'aurais aussi eu l'entier 2874 dans la mémoire.
Comment je fais la différence ? Dans la conception du type ou du programme. En C, comme le compilateur ne fait pas la différence entre un int et un nombre en point fixe, et ne connaît pas la valeur de l'exposant, c'est dans la conception du programme. En C++, on pourrait dire au compilateur quelle est la valeur de l'exposant pour qu'il vérifie qu'on ne se trompe pas.
Ajouté le 23/06/2021 à 10:08 :
Tu ne « peux pas », tu n'as que des nombres avec exactement 3 chiffres après la virgule.
Ici, pas de panique, tu peux écrire 287400 ce qui représente 287.400 = 287.4. Mais ça ne marcherais pas si tu voulais 4 chiffres après la virgule par exemple.
Citer : Posté le 23/06/2021 10:25 | #
Pour faire simple, et c'est ce que Lephe explique, ça revient exactement à compter en centimes.
Cette baguette coûte 110 centimes, ce croissant 150 centimes. Au total j'ai payé 260 centimes, mais vu que je sais que ce sont des centimes, je sais que ça représente 2,60 €
Et si je veux plus de précision, par exemple sur le prix de l'essence, je peux compter en millièmes d'euros. Le litre de gasoil coûte 1345 millièmes d'euros, si j'en prends 10 litres ça donne 13450 millièmes d'euro, et vu que je sais que je compte en millièmes, ça fait 13,45 €
Mais dans tous les cas il faut se fixer une convention et s'y tenir, parce que si je ne te dis pas dans quel unité est le nombre 4235, tu ne saura pas si c'est 4235 €, 42,35 € ou 4,235 €.
Citer : Posté le 23/06/2021 10:26 | #
C'est plus clair comme ça, merci Darks. <3
Citer : Posté le 23/06/2021 16:35 | #
Super tuto, bien construit et fourni
Je me demandais juste : l'objectif de tout ceci (la virgule fixe) est d'éviter toute opération sur des float/double qui soit gourmande non ?
Alors je ne comprends pas trop cette ligne-ci :
Si j'ai bien compris, il serait plus aisé (ça prendrait moins de temps) d'exécuter ce code si on calculait directement le résultat (pour éviter de faire une opération avec une virgule flottante), mais un problème qui surviendrait serait de se souvenir ce qu'on y a fait...
Encore jamais fait d'add-in, mais ce genre de tuto me dire super envie de commencer ! (mais bon, vu tout ce que je dois finir en ce moment...)
Citer : Posté le 23/06/2021 16:40 | #
Exactement, c'est beaucoup plus simple, et donc rapide, de faire des opérations en point fixe. Lorsque ni la plage élevée ni la précision élevée des float n'est nécessaire, c'est juste un énorme boost en performances. Comme application, prenez Windmill, qui utilise quasiment exclusivement du point fixe (et à l'époque toute la conversion en point fixe avait été un gros facteur de progrès).
Comme je l'ai mentionné ensuite, le compilateur précalcule le résultat à ta place. C'est le genre de choses qu'il peut faire automatiquement et sans poser de questions. Tu peux être sûr que dès -O1 c'est identique à fixed y = 106496.
Lance-toi un jour, tu ne le regretteras pas !
Citer : Posté le 23/06/2021 18:14 | #
Mais …
fixed x = 1 * 65536; // x= 1.0 en point fixe 16:16
Attend, 1 x 65536 = 65536 non ?
Et pourquoi 2 puissance 16 ? Il y a un chiffre donc 1 bit ?!
Citer : Posté le 23/06/2021 18:24 | #
Dans cet exemple je compte non pas en centièmes, ni en millièmes, mais en 65536èmes. Donc oui, en effet, ça fait une valeur entière égale à 65536, qui représente autant de 65536èmes, pour un total de 1.
De façon générale, tu peux écrire
et sauf cas problématiques (overflows et subtilités de précision), ça te donne <valeur> comptée en 65536èmes.
Pourquoi 2¹⁶ ? Parce que j'ai choisi de compter en 65536èmes. J'aurais pu prendre un autre diviseur, mais pour que le code soit performant il faut une puissance de 2, et pour diverses autres raisons 2¹⁶ ça marche bien.
Citer : Posté le 24/06/2021 09:48 | #
Merci beaucoup pour le travail derrière l'article. La lecture est fluide et le contenu aussi précis qu'intéressant.
Ce à quoi je ne m'attendais pas : "Notez juste qu'en Basic CASIO toutes les variables sont des nombres en point flottant en base 10".
Est-il possible d'avoir un microcontrolleur disposant d'un FPU qui ne sache traiter que les float et non les double ? L'ordre de grandeur des temps présentés dans l'article entre l'addition de deux entiers et celui d'un flottant pour une calculatrice qui n'a pas de FPU me fait penser à celui déjà rencontré sur certains cibles entre les float et les double émulés.
Dans le cas d'une target avec FPU, les nombres en point fixe font-ils concurrence aux nombre en virgule flottante ? La réponse est peut-être hardware specific...
Hâte de lire les prochains articles !
La Planète Casio est accueillante : n'hésite pas à t'inscrire pour laisser un message ou partager tes créations !
Citer : Posté le 24/06/2021 10:28 | #
Merci beaucoup ! Oui en Basic tout est en point flottant, ce qui explique la notation scientifique et pourquoi certaines fonctions comme MOD() ne marchent plus quand les valeurs sont trop grandes (exposant > 0).
C'est possible oui, quoique je ne connais pas d'exemple. Peut-être sur une machine un peu ancienne ?
Pour les machines avec FPU, non, on utilise souvent du point flottant sans se poser de question. Le pipelining fait que ça ne pose pas de problèmes de performances particuliers pour les calculs mondains, et les gros algos numériques sont de toute façon déjà optimisés. Quand le problème devient complexe les questions portent plus sur les algorithmes, la gestion de la mémoire/des caches/etc que sur le format numérique.
D'ailleurs les flottants ont déjà des problèmes bien connus de précision et de stabilité numérique, qui ne feraient qu'empirer fabuleusement avec la plage limitée du point fixe, donc ce serait beaucoup d'effort d'ingénierie en plus.
Par contre, dans les shaders qui tournent sur les GPU (pas que pour les questions d'effets de lumière ; aujourd'hui des shaders de différents types sont utilisés dans tout ce qui touche au GPU, y compris le rendu graphique et le calcul), on utilise beaucoup des float au lieu de double, probablement pour une combinaison de place occupée en mémoire et de vitesse de calcul.
Le format est plus modeste que des articles bien présentés et très accessibles, et je partirai sans doute dans des trucs plus compliqués sur certains sujets, donc n'hésitez pas à demander si d'autres choses ne sont pas claires
Citer : Posté le 24/06/2021 15:21 | #
Notons que dans le cas des GPU, les performances sont souvent exprimées sous la forme de trois valeurs en flops : half, single, double.
Et effectivement plus on réduit la taille, plus les performances sont élevées
Exemple ici avec la RX 5700XT : https://www.techpowerup.com/gpu-specs/radeon-rx-5700-xt.c3339
Citer : Posté le 24/06/2021 19:39 | #
Ah oui quand même ça fait un ratio de x16 entre le float et le double ! Je ne m'y attendais pas. Le ratio de x2 entre le half et le float correspond plus à ce que j'imaginais.
Citer : Posté le 27/06/2021 00:08 | #
Précalcul des divisions (Maths — ★★☆)
Introduction
Exploration mathématique et preuve
Implémentation
Cette technique est longue mais c'est très progressif. Ne soyez pas intimidé·e ! Si c'est trop dur regardez juste la partie «Implémentation», vous aurez le code sans les maths.
Ajout plus tard : Sur injonction de Dark Storm j'ai fait plus de tests, et on dirait que le compilateur fait parfois cette optimisation tout seul ! Regardez bien le code assembleur avant de vous lancer, ce serait dommage d'optimiser pour rien.
Motivation pour éviter les divisions
Les opérations arithmétiques sur les entiers ne coûtent pas très cher ; les additions, soustractions, négations, décalages de bits, et autres carry et overflow sont des opérations basiques qui occupent l'ALU pendant 1 seul cycle. Les multiplications prennent plus longtemps ; 1 cycle pour les multiplications 16 bits × 16 bits → 32 bits (très rares en C, on s'amuse rarement à ça) et 3 cycles pour les multiplications 32 bits × 32 bits → 64 bits. Le multiplieur est séparé de l'ALU donc on peut quand même travailler pendant les 2 cycles de délai tant qu'on n'a pas besoin du résultat de la multiplication, donc en général ça ne pose aucun problème de performances.
La division, par contre, est une autre affaire. Diviser est de loin l'opération la plus compliquée, et la méthode classique ressemble à ça.
r=0
Pour i=31 à 0:
Si q × 2^i ≤ p:
r ← r + 2^i
p ← p - q × 2^i
C'est exactement comme on le fait à l'école primaire ; typiquement pour diviser 84934 par 23 on commence à regarder combien de fois 23000 passe dans 84934 (→ 3) ; on soustrait 3 fois 23000, il reste 15394 et maintenant on cherche combien de fois 2300 passe dedans (→ 6) ; et ainsi de suite jusqu'à trouver le résultat 3692.
En binaire, c'est un tout petit peu plus simple parce que chaque chiffre ne peut être que 0 ou 1, donc on ne demande pas « combien de fois » chaque multiple du diviseur passe, mais simplement si le multiple passe. C'est matérialisé par la condition « q × 2^i ≤ p ».
Dans le processeur SuperH, la méthode est similaire. On calcule un bit à la fois, et donc ça prend autant d'opérations qu'il y a de bits dans le résultat. Spécifiquement, la séquence d'instruction est la suivante :
Quand le résultat fait 8 ou 16 bits on peut récupérer T lors du div1 suivant (par un mécanisme bien foutu), mais GCC divise toujours les int de 32 bits et pour ceux-là il faut récupérer T à la main avec rotcl. La fonction _sdivsi3 de GCC qui calcule la division entre deux entiers signés ressemble donc à ça :
2: 60 53 mov r5,r0
4: 20 08 tst r0,r0
6: 89 48 bt 9a <div0>
8: e2 00 mov #0,r2
a: 21 27 div0s r2,r1
c: 33 3a subc r3,r3
e: 31 2a subc r2,r1
10: 23 07 div0s r0,r3
12: 41 24 rotcl r1
14: 33 04 div1 r0,r3
16: 41 24 rotcl r1
18: 33 04 div1 r0,r3
1a: 41 24 rotcl r1
1c: 33 04 div1 r0,r3
1e: 41 24 rotcl r1
20: 33 04 div1 r0,r3
22: 41 24 rotcl r1
24: 33 04 div1 r0,r3
26: 41 24 rotcl r1
28: 33 04 div1 r0,r3
2a: 41 24 rotcl r1
2c: 33 04 div1 r0,r3
2e: 41 24 rotcl r1
30: 33 04 div1 r0,r3
32: 41 24 rotcl r1
34: 33 04 div1 r0,r3
36: 41 24 rotcl r1
38: 33 04 div1 r0,r3
3a: 41 24 rotcl r1
3c: 33 04 div1 r0,r3
3e: 41 24 rotcl r1
40: 33 04 div1 r0,r3
42: 41 24 rotcl r1
44: 33 04 div1 r0,r3
46: 41 24 rotcl r1
48: 33 04 div1 r0,r3
4a: 41 24 rotcl r1
4c: 33 04 div1 r0,r3
4e: 41 24 rotcl r1
50: 33 04 div1 r0,r3
52: 41 24 rotcl r1
54: 33 04 div1 r0,r3
56: 41 24 rotcl r1
58: 33 04 div1 r0,r3
5a: 41 24 rotcl r1
5c: 33 04 div1 r0,r3
5e: 41 24 rotcl r1
60: 33 04 div1 r0,r3
62: 41 24 rotcl r1
64: 33 04 div1 r0,r3
66: 41 24 rotcl r1
68: 33 04 div1 r0,r3
6a: 41 24 rotcl r1
6c: 33 04 div1 r0,r3
6e: 41 24 rotcl r1
70: 33 04 div1 r0,r3
72: 41 24 rotcl r1
74: 33 04 div1 r0,r3
76: 41 24 rotcl r1
78: 33 04 div1 r0,r3
7a: 41 24 rotcl r1
7c: 33 04 div1 r0,r3
7e: 41 24 rotcl r1
80: 33 04 div1 r0,r3
82: 41 24 rotcl r1
84: 33 04 div1 r0,r3
86: 41 24 rotcl r1
88: 33 04 div1 r0,r3
8a: 41 24 rotcl r1
8c: 33 04 div1 r0,r3
8e: 41 24 rotcl r1
90: 33 04 div1 r0,r3
92: 41 24 rotcl r1
94: 31 2e addc r2,r1
96: 00 0b rts
98: 60 13 mov r1,r0
0000009a <div0>:
9a: 00 0b rts
9c: e0 00 mov #0,r0
9e: 00 09 nop
Le coût de cette procédure comprend le test pour la division par 0, l'appel de la fonction (pas négligeable, surtout que les conventions d'appel obligent à sauvegarder des registres), et bien sûr 32 coups de div1 et rotcl, le tout comptant probablement dans les ~80 cycles. C'est pas la mort, mais dès qu'on met ça dans une boucle ça commence à piquer pas mal.
Le principe de cette technique d'optimisation est de remplacer les divisions par des multiplication 32×32 → 64 bits, qui prend en tout autour de ~10 cycles, dans les cas où on connaît le diviseur à l'avance.
Multiplication en point fixe
Imaginons que je veux diviser un entier n par une constance c=10. Je peux diviser par 10... ou je peux multiplier par 1/10. Bien sûr 1/10 ce n'est pas un nombre entier, mais 2³²/10 c'est presque un nombre entier, si je peux multiplier n par ça j'obtiendrais à peu près 2³² × n/10, et ensuite c'est une affaire de diviser par 2³² (décalage de bits) pour récupérer le résultat.
Il faut faire attention quand même, parce que 2³²/10 n'est pas exactement un entier, et donc on n'aura qu'une approximation. L'erreur d'approximation se fait multiplier par n, et il faut qu'elle soit assez petite pour que le résultat reste correct.
Dans cette technique, je vais dérouler les calculs pour vous montrer comment on approche une question mathématique de ce genre, et comment on peut trouver la solution par le raisonnement plutôt que par le test.
Je présuppose que les signes sont traités séparément (donc les nombres sont positifs). Si c est une puissance de 2, on peut faire un décalage de bits, donc on ne s'embête pas avec cette méthode. Pour l'approximation, je vous demande d'admettre qu'il faut arrondir vers le haut sinon la méthode ne marche pas du tout.
Il y a pas mal d'entiers et de nombres décimaux dans les calculs, donc je vais utiliser les notations courtes pour les parties entières : ⌊x⌋ représente la partie entière inférieure de x, ⌈x⌉ sa partie entière supérieure, et {x} sa partie décimale (on a donc x = ⌊x⌋ + {x}, et aussi ⌊x⌋ ≤ x ≤ ⌈x⌉).
Et voilà les variables du problème :
Ce qu'on veut faire, c'est profiter qu'on connaisse la valeur de c pour calculer ⌊n/c⌋ en remplaçant la division par une multiplication par x. La valeur de x qu'on a choisi c'est ⌈2^m / c⌉ = 2^m/c + ϵ, et donc la multiplication va donner
À la fin on va re-diviser par 2^m (avec un décalage de bits) et on perdra au passage la partie décimale de ce résultat. Autrement dit, on obtiendra ⌊nx/2^m⌋, et toute la technique consiste à trouver la bonne valeur de x (plus spécifiquement la bonne valeur de m) qui fait que ⌊nx/2^m⌋ = ⌊n/c⌋.
Pourquoi le format 0:32 n'est pas suffisant
L'exemple introductif prend x = ⌈2³²/c⌉, ce qui revient, si vous avez lu la technique sur le calcul en point fixe, à représenter 1/c en point fixe 0:32. Il est facile de voir que ce n'est pas assez précis ; l'erreur causée par la partie entière est dans l'intervalle [0:2⁻³²[, donc quand on multiplie par un n qui peut être n'importe quelle valeur entre 0 et 2³²-1, l'erreur se fait multiplier par 2³²-1 (on regarde le pire cas) et on obtient une erreur qui peut être à peu près n'importe quel nombre dans [0:1[. C'est trop, puisque ça peut parfaitement nous faire passer d'un entier à l'autre (par exemple une erreur de 0.3 transforme un résultat correct de 88/10 = 8.8 tronqué à 8, en un résultat incorrect de 88/10 + 0.3 = 9.1 tronqué à 9).
La raison pour laquelle ça ne marche pas, c'est que quand on calcule ⌈2³²/c⌉ on perd trop de précision. Par exemple si je prend 1 million comme constante, j'obtiens 4295, un nombre de 13 bits, alors que j'ai 32 bits à ma disposition. Vous pourriez être tenté·e de penser que quand on divise par 1 million le résultat est forcément petit, donc il y a peu de bits à calculer, donc c'est normal si on a besoin de peu de bits de précision. C'est raisonnable, mais malheureusement ce n'est pas vrai, et on va voir pourquoi.
Une approximation correcte avec m = 32+k+1
Lorsqu'on calcule le produit nx, il y a deux facteurs qui contribuent à la partie décimale du résultat, mis en valeur dans l'équation présentée dans l'introduction :
À la fin on veut diviser par 2^m (par décalage de bits, ce qui arrondit vers le bas) et obtenir ⌊n/c⌋, ce qui signifie que ces deux termes doivent disparaître (le troisième terme ⌊2^m × n/c⌋ nous donnant précisément ⌊n/c⌋ après le décalage de bits). Pour ça, il faut que {2^m × n/c} + nϵ < 2^m.
La difficulté c'est que dans la division, {n/c} peut être très proche de 1 (du genre 0.99 si jamais n/c = 5.99) et il faut que l'erreur soit malgré tout assez petite pour éviter qu'on déborde sur le multiple de 2^m supérieur.
L'astuce pour s'en sortir est de remarquer que {n/c} ne peut pas être n'importe quelle valeur : c'est soit 1/c, soit 2/c, etc, jusqu'à (c-1)/c. Par exemple si c=10, {n/c} ne peut être que 0.1, 0.2, ... ou 0.9 (et donc pas 0.99 par exemple). Si on arrive à garantir que nϵ < 2^m/c alors on réussira à éviter les erreurs, puisque même dans le pire cas on aura {2^m × n/c} + nϵ ≤ 2^m(c-1)/c + nϵ < 2^m, et donc tout disparaîtra sans laisser de traces lors du décalage de bits final.
On veut donc que le résultat de la multiplication soit précis à nϵ < 2^m/c, ce qui montre bien que plus c est grand (par exemple 1 million) plus la précision nécessaire est grande. Ainsi, contrairement à ce que notre intuition nous dit, même quand le diviseur est grand on a besoin d'absolument tous les bits disponibles dans x pour maintenir une précision élevée.
Maintenant qu'on connaît la précision nécessaire, on peut déterminer la bonne valeur de m. Comme n peut aller jusqu'à 2³²-1, pour avoir nϵ < 2^m/c il faut que ϵ < 2^(m-32)/c. On sait de plus que c est coincé entre 2^k et 2^(k+1), et qu'il faudra donc une précision de 1/2^(k+1) pour atteindre 1/c (comme on est en binaire on n'a pas d'intermédiaire entre k bits et k+1 bits ; on est obligés de prendre k+1). Par conséquent, il faut que ϵ < 2^(m-32-k-1).
On se rappelle que 0 ≤ ϵ < 1 compte tenu de la façon dont on a formé notre calcul, et donc pour que ça marche à tous les coups on est obligés de prendre au moins m = 32+k+1. Et ça, contrairement à notre idée initiale de m = 32, ça va passer.
On peut revenir à l'exemple de c = 1 million pour comprendre l'importance de m = 32+k+1. Avec m=32, on avait x = 4295, un nombre pitoyable de 13 bits qui laissait 19 bits sur 32 inutilisés. Avec m=32+k+1, on magnifie suffisament 1/c pour exploiter toute la place disponible, et on obtient x = 4503599628, un nombre qui ne gâche aucun des 32 bits.
Le 33ème bit de l'approximation
Si vous avez le regard affûté, vous aurez remarqué que x = 4503599628 ça ne tient pas dans 32 bits. En fait, m = 32+k+1 ça nous oblige à calculer 33 bits de 1/c, et donc dans l'absolu il nous manque 1 bit pour représenter x.
Cependant, le 33ème bit vaut toujours 1. C'est le principe en fait ; on a choisi m = 32+k+1 spécifiquement pour amener le premier bit non nul de 1/c à la 33ème position et stocker toutes les décimales binaires importantes dans les 32 bits disponibles.
Du coup on peut facilement stocker non pas x mais x - 2³² et ensuite faire une addition pendant le calcul pour rajouter le 2³². Le code est le suivant :
n_sur_c = ((uint64_t)n_sur_c + n) >> (k + 1);
Vous pouvez voir qu'on décale n*x de 32+k+1 bits, comme promis ; on s'arrête en chemin après 32 pour rajouter n, ce qui revient à ajouter 2³²×n au résultat (ce qui compense le 2³² qui n'était pas présent dans x lors de la multiplication).
Vérifions que ça marche
Maintenant qu'on a bien choisi notre valeur, on peut (re)vérifier que tout marche bien. Puisque m = 32+k+1, notre approximation de 1/c est
x ∈ [0..2³²[ comme promis. Et comme on a arrondi vers le haut, on a bien le ϵ qu'on veut.
Maintenant, on peut prendre n'importe quelle valeur de n entre 0 et 2³²-1, et calculer le produit avec l'addition spéciale.
De ce résultat, on veut éliminer q durant le décalage de 32+k+1 bits vers la droite, et garder uniquement ⌊n/c⌋. On veut donc que q < 2^(32+k+1), et comme {n/c} peut atteindre (c-1)/c, on veut avoir une garantie forte que
Et tout se goupille puisque ϵ < 1, c < 2^(k+1), et n < 2³². Autrement dit, une fois le décalage de 32+k+1 bits vers la droite effectué, on récupère exactement ⌊n/c⌋.
Et voilà !
Le cas particulier où m = 32+k suffit
Le 33ème bit est un peu frustrant quand même, il nous oblige à faire une addition en plus. Séparer le décalage de bits en deux est inévitable (on verra pourquoi quand on comparera le code assembleur), mais l'addition c'est toujours un cycle en plus dont on aimerait bien se passer. Heureusement, il se trouve qu'on peut parfois l'éviter.
Il n'y a pas de miracles en fait. On a choisi m = 32+k+1 pour représenter un maximum de bits de 1/c, à tel point que le bit de poids fort (la k+1ème décimale) est arrivé à la 33ème position pour laisser la place au 32+k+1ème. Vous allez me dire, ok, mais si le 32+k+1ème bit c'est 0, ça ne sert à rien. Et oui, c'est exact. Si le 32+k+1ème bit est 0, on peut prendre m = 32+k et éviter d'avoir un 33ème bit implicite (parfois ça marche aussi par pure chance).
Dans ce cas, le code ressemble à ça :
Il y a une autre situation très courante où c'est suffisant. Si n est un entier signé, alors la plage de valeurs possibles pour n n'est pas [0:2³²[, mais [2⁻³¹:2³¹[. Une fois les signes gérés, on se retrouve avec [0:2³¹[. Comme n est 2 fois plus petit qu'avant, l'erreur nϵ a 2 fois plus de marge, ce qui est précisément ce dont on a besoin pour réduire la précision de 1/c d'un bit. Et donc là aussi, on peut prendre m = 32 + k.
La preuve est encore plus simple ; on prend
Puis on considère n'importe quel n ∈ [0..2³¹[, ce qui nous donne le produit
et une fois encore q < 2^(32+k) puisque {n/c} ≤ (c-1)/c et nϵ < 2^(32+k)/c grâce à la borne plus stricte sur n.
Code de précalcul et code de division
D'abord voici la fonction qui calcule k (le logarithme en base 2 de c) ; elle divise simplement c par 2 jusqu'à atteindre 0.
{
int log = -1;
while(c) {
log++;
c >>= 1;
}
return log;
}
Cas général — Dans le cas général (m = 32 + k + 1), on calcule x de la façon suivante :
int k = u32log2(c);
uint32_t x = ((1ull << (32+k+1)) + (c-1)) / c;
Le fait d'ajouter (c-1) avant de diviser par c produit un arrondi vers le haut au lieu d'un arrondi vers le bas (c'est une technique classique).
Ensuite la multiplication se fait comme ceci :
n_sur_c = ((uint64_t)n_over_c + n) >> (k + 1);
Par exemple, avec c = 1729 :
n_sur_1729 = ((uint64_t)n * 0x2f3b5f81) >> 32;
n_sur_1729 = ((uint64_t)n_sur_10 + n) >> 11;
Cas restreint — Dans le cas restreint (m = 32 + k), on calcule x comme ça :
int k = u32log2(c);
uint32_t x = ((1ull << (32+k)) + (c-1)) / c;
La multiplication se fait en un seul coup :
Par exemple avec n=10, qui a la chance de marcher avec tous les n (même non signés) parce que son x normal est pair :
n_sur_10 = ((uint64_t)n * 0xcccccccd) >> 35;
Avec ça, vous pouvez tout précalculer ; mais continuez à lire, il y a des options plus rapides.
Programme pour générer les inverses et expérimenter avec les formules
Le programme ci-dessous calcule les valeurs de k et x pour les deux méthodes (m = 32+k+1, et m = 32+k), et vérifie toutes les divisions possibles sur 32 bits pour déterminer si m = 32+k suffit. Les tests sont parallélisés sur 4 threads et prennent moins de 20 secondes sur mon modeste Core i3-8100.
Programme : fixed-point-inverse.c
Compilez le programme avec les optis et la bibliothèque de threads.
Le programme rapporte des erreurs même pour les méthodes qui échouent complètement, donc n'hésitez pas à jouer avec et tester des variations. Par exemple, m = 32+k+1 peut être plus rapide en assembleur si k+1 vaut 2, 8 ou 16 parce que le jeu d'instructions peut faire ces décalages de bits sans spécifier la valeur dans un registre.
Voilà ce que donne le programme sur un diviseur où m = 32+k suffit sur toutes les entrées :
Inverse of 10 (k=3) with 32+k+1 places: 0x9999999a
success!
Inverse of 10 (k=3) with 32+k places: 0xcccccccd
success!
Et voilà un exemple sur un diviseur où m = 32+k ne marche pas sur tout le monde (mais quand même sur tous les nombres signés, puisqu'on a prouvé que c'était toujours le cas) :
Inverse of 1729 (k=10) with 32+k+1 places: 0x2f3b5f81
success!
Inverse of 1729 (k=10) with 32+k places: 0x979dafc1
sample error: bffff98f -> 001c6d90 (instead of 001c6d8f)
sample error: fffff9aa -> 0025e76b (instead of 0025e76a)
errors: 956331, but none in signed inputs
N'hésitez pas à noter ces résultats dans des fonctions dans un header, histoire pas avoir à copier les constantes partout :
static inline uint32_t div_10(uint32_t n) {
return ((uint64_t)n * 0xcccccccd) >> 35;
}
static inline uint32_t div_1729(uint32_t n) {
uint32_t r = ((uint64_t)n * 0x2f3b5f81) >> 32;
return ((uint64_t)r + n) >> 11;
}
Vous voilà paré·e pour la version C de cette technique.
Comment le compilateur C++ peut tout calculer à votre place
Le compilateur C++ est très fort. Déterminer la valeur de x est un truc très facile qu'il peut parfaitement faire à votre place. Déterminer si m = 32+k suffit est plus difficile, mais on peut se contenter de détecter le cas où x est pair et laisser l'addition dans les quelques autres même si elle n'est pas nécessaire. L'avantage de tout ça c'est qu'il n'y a rien à précalculer pour nous, le compilateur fait tout.
On a besoin de trois outils du C++ pour faire ça :
Voilà comment ça se passe. On commence par définir les fonctions utiles pour calculer ; le u32log2(), la fonction fixdiv_k1() qui calcule x dans le cas où m = 32+k+1, et la fonction fixdiv_k() qui calcule x dans le cas où m = 32+k.
{
int log = -1;
while(c) {
log++;
c >>= 1;
}
return log;
}
constexpr uint32_t fixdiv_k1(uint32_t c)
{
int k = u32log2(c);
return ((1ull << (32+k+1)) + (c-1)) / c;
}
constexpr uint32_t fixdiv_k(uint32_t c)
{
int k = u32log2(c);
return ((1ull << (32+k)) + (c-1)) / c;
}
J'ai mis constexpr sur les trois, ce qui indique au compilateur que les fonctions peuvent être évaluées (exécutées) à la compilation. Bien sûr, le compilateur ne me croit pas sur parole ; il y a des règles qu'on doit respecter, et évidemment pour que la fonction puisse être évaluée à la compilation il faut que tous les paramètres soient connus à la compilation.
Ensuite, vient l'outil vraiment magique, le template. La façon dont ça marche est que je crée une fonction fixdiv mais je dis qu'il y a un paramètre c qui est à déterminer à la compilation (le paramètre du template). Chaque fois qu'on appelera fixdiv, le compilateur générera une copie de la fonction où c sera hardcodé (naturellement il faut que c soit une constante).
constexpr uint32_t fixdiv(uint32_t n)
{
int const k = u32log2(c);
uint32_t const x = fixdiv_k1(c);
if constexpr ((x & 1) == 0) {
uint32_t const x = fixdiv_k(c);
return ((uint64_t)n * x) >> (32 + k);
}
else {
uint32_t r = ((uint64_t)n * x) >> 32;
return ((uint64_t)r + n) >> (k + 1);
}
}
Le corps de la fonction implémente la technique complète de précalcul. Le compilateur commence par calculer x avec fixdiv_k1 (c'est-à-dire m = 32+k+1). Si x est pair, on décide d'utiliser m = 32+k donc on recalcule x et on utilise la méthode directe. Sinon, on reste sur la méthode avec l'addition supplémentaire.
Le if constexpr indique que la condition peut (doit) être calculée à la compilation (ce qui possible puisque tout ne dépend que de la valeur de c), et donc le compilateur ne génère qu'une des deux branches. Grâce aux optimisations, tous les appels à fixdiv_k1 et fixdiv_k sont remplacés par les résultats, si x est pair le calcul initial de fixdiv_k1 est entièrement éliminé, et on obtient vraiment la même chose qu'en C, mais sans avoir à travailler.
Cerises sur le gâteau, comme la fonction générée est aussi constexpr, elle est inlinée par défaut, ce qui veut dire qu'on a exactement le même résultat que si on avait précalculé la constante à la main et copié le code de la méthode appropriée à l'endroit de la division ; et si n est une constante le compilateur optimise aussi le calcul comme avec la division normale !
Programme : fixed-point-inverse-template.cpp
Comparaison du code assembleur
On peut obtenir le code assembleur de la division en définissant une fonction quelconque qui l'utilise (ou en prenant la version C, les résultats sont les mêmes) :
{
return fixdiv<10>(n);
}
Et on obtient le code avec objdump -d :
% sh-elf-objdump -d fpit.o
(...)
00000000 <__Z6div_10m>:
0: d1 02 mov.l c <__Z6div_10m+0xc>,r1 ! cccccccd
2: 34 15 dmulu.l r1,r4
4: 00 0a sts mach,r0
6: 40 09 shlr2 r0
8: 00 0b rts
a: 40 01 shlr r0
(...)
La fonction commence par charger la valeur x = 0xcccccccd (dans r1), multiplie n (r4) par x, récupère les 32 bits du haut (mach) puis décale de 3 bits en deux temps (il n'y a aucune instruction qui décale de 3 bits d'un coup, l'autre option aurait été de charger la valeur 3 dans un registre et d'utiliser shld).
Temps total, sauf erreur de ma part : 9 cycles (3 pour le mov.l avec la dépendance RAW, 3 pour dmul.u avec la dépendance sur mach, ensuite 1 par calcul ; peut-être que sts et shlr2 parallélisent par forwarding et descendent le total à 8).
C'est le bon moment pour expliquer où est passé le décalage de bits par 32. En fait lorsqu'on multiplie 32 bits × 32 bits → 64 bits, le résultat est séparé en deux registres : mach (32 bits du haut) et macl (32 bits du bas). Le SuperH n'a pas de registres de 64 bits, donc on ne peut manipuler qu'une moitié à la fois. Le compilateur transforme donc le décalage de 32 bits du résultat vers la droite en le fait de lire mach et d'ignorer macl. Si on voulait un mélange des bits des deux ce serait assez casse-pieds, heureusement le cas ne se présente pas avec cette méthode.
Et voilà, avec un peu d'effort de preuve, du précalcul si on fait du C ou un template constexpr si on fait du C++, on a réduit d'un facteur 8-10 le temps nécessaire pour diviser par des constantes. Diviser par 10 en particulier est super utile pour les fonctions d'affichage de nombres en décimal.
Citer : Posté le 30/06/2021 11:28 | #
Excellente explication, merci
Y'a un truc que je ne comprends pas, c'est pourquoi le compilo ne le fait déjà pas tout seul ?
Il sait très bien remplacer *4 par <<2 par exemple. Donc remplacer /10 par la dizaine d'instructions dont tu parles, je ne vois pas où est la difficulté supplémentaire.
uint32_t result = n / 10;
Y'a pas besoin de contexte, aucune astuce liée au matériel, un poil de calcul préliminaire pour trouver les bonnes constantes, mais je suppose que tu n'es pas le premier à avoir optimisé des divisions entières. Le template ne me semble pas si spécifique que ça, donc je suis surpris qu'il ne soit pas déjà intégré dans le compilateur.
Citer : Posté le 30/06/2021 11:53 | #
Alors j'ai pas inventé la méthode non. Enfin je l'ai réinventée (si j'avais su...), mais ensuite j'ai retrouvé des traces assez facilement jusque dans des vieux papiers de recherche.
Il faut voir d'abord que la méthode n'est pas entièrement directe. D'une part là je n'ai présenté que la version non signée, sachant que la version signée nécessite test et inversion de signe en plus. Et d'autre part j'ai éludé la division par 0, qui a une sémantique différente selon les architectures (exception en point flottant sur x86 par exemple). Et du coup tu peux probablement pas faire l'opti directement sur l'IR parce que tu sais pas exactement quel comportement tu veux obtenir, il faudrait la faire dans le back-end (ce qui ouvre tout un paquet d'explications parce que tu te doutes que le backend SuperH est pas le plus perfectionné).
Ensuite, ce n'est pas forcément plus rapide sur toutes les architectures. Ici la division prend 80 cycles, mais sur un x86 par exemple tu peux payer du genre 3-6 cycles pour la multiplication et 12-40 cycles pour la division donc c'est pas évident qu'avec le branchement signé et le surcoût de charger une grosse constante (au lieu par exemple de 10) ce soit plus rapide.
Enfin, les compilateurs l'implémentent parfois, et même si je n'ai pas d'exemple hyper concret sous la main, le manuel d'optimisation x86 d'Intel en parle justement pour les divisions 128-bits (page 360) :
Modern compilers can transform expressions of integer division in high-level language code with a constant divisor into assembly sequences that use IMUL/MUL to replace IDIV/DIV instructions. Typically, compilers will replace a divisor value that is within the range of 32-bits if the divisor value is known at compile time. If the divisor value is not known at compile time or the divisor is greater than those repre-sented by 32-bits, DIV or IDIV will be generated.
The latency of a DIV instruction with a 128-bit dividend is quite long. For dividend values greater than 64-bits, the latency can range from 70-90 cycles.
The basic technique that a compiler employs to transform integer division into 128-bit integer multiplica-tion is based on the congruence principle of modular arithmetic. It can be easily extended to handle larger divisor values to take advantage of fast 128-bit IMUL/MUL operation.
La section 20.2.6.2 référence ce paragraphe dans le contexte de toutes les divisions, et l'exemple 14-21 montre le cas exact de la division par 10 pour convertir les donnés numériques en texte (le "%d" essentiellement).
Après ça reste des trucs assez tordus, tu peux sans doute trouver facilement un compilo qui le fait mais tu ne peux pas raisonnablement espérer qu'un compilateur tiré au hasard le fasse.
Sur SuperH spécifiquement, GCC a une option de stratégie de division, mais ça ne concerne que les cas génériques :
Set the division strategy to be used for integer division
operations. STRATEGY can be one of:
'call-div1'
Calls a library function that uses the single-step division
instruction 'div1' to perform the operation. Division by zero
calculates an unspecified result and does not trap. This is
the default except for SH4, SH2A and SHcompact.
'call-fp'
Calls a library function that performs the operation in double
precision floating point. Division by zero causes a
floating-point exception. This is the default for SHcompact
with FPU. Specifying this for targets that do not have a
double precision FPU defaults to 'call-div1'.
'call-table'
Calls a library function that uses a lookup table for small
divisors and the 'div1' instruction with case distinction for
larger divisors. Division by zero calculates an unspecified
result and does not trap. This is the default for SH4.
Specifying this for targets that do not have dynamic shift
instructions defaults to 'call-div1'.
Il utilise call-table de base, c'est pour ça que diviser par 0 produit un TLB miss (l'accès à la table donne n'importe quoi).
Donc pour répondre à ta question, le compilo ne le fait pas tout seul probablement pour une combinaison de (1) c'est pas toujours rentable selon les architectures, (2) la méthode s'alourdit si on veut exactement la même sémantique que la division, sans hypothèses, et (3) implémentation incomplètes, back-ends pas très développées, etc.
La morale c'est que si tu veux que ton code fasse exactement ce que tu lui dis, même que c'est quelque chose d'exotique, il vaut mieux être le plus explicite possible. Si l'opti est tellement dépendante du back-end que le compilateur n'en implémente rien tu as aussi à gagner à le faire en assembleur (et les parties assembleur de ce topic montreront qu'il y a pas mal à gagner, sur SH4 en tous cas).
Citer : Posté le 30/06/2021 17:10 | #
Ceci étant dit, typiquement sur la division par 10 l'opti de passer par des templates est grosso-modo équivalente à l'opti de la call-table ? Ou du moins c'est pas si horrible de faire un / 10 dans un bout de code que ce que tu laissais imaginer au début.