À la question des courbes de Bézier, de la vitesse Arduino et d'un site intéressant, ou comment j'ai passé le week-end

«N'importe qui peut résoudre le paradoxe gris avec des dauphins, et vous essayez de le faire sans dauphins. "

En fait, j'avais prévu de passer le week-end d'une manière légèrement différente, d'aller à Copter Huck (pas que j'étais un fan de copters, juste pour voir ce que les jeunes inventaient, pour passer du temps comme ça), mais la sœur aînée était catégoriquement contre. Bien sûr, j'ai insisté (c'est-à-dire que j'ai ri deux fois et dit: "Eh bien, peut-être ... ce sera amusant, de toute façon"), mais elle était implacable, et quand ma femme a pris son parti, il n'y avait aucune chance de voyage. Bon, d'accord, "je n'en voulais pas vraiment", mais je me suis assis un peu sur un casse-tête drôle du domaine de la programmation, que j'ai moi-même inventé, dont je fais rapport.

(Remarque nécessaire - le week-end précédent était destiné, c'est toujours comme ça - écrire un programme nécessite quelques heures, rédiger un rapport à ce sujet et cinq jours de voyage dans les transports publics ne sont pas terminés.)

Dans un article récent, l'auteur a abordé le problème de l'accélération (entre autres) du calcul des courbes de Bézier (KB) sur MK avec des paramètres relativement faibles. Eh bien, en fait, ces paramètres sont au niveau de l'ordinateur central moyen des années 70, mais à l'heure actuelle, ils sont considérés comme clairement insuffisants. À la suite de certaines actions, l'auteur a réussi à accélérer quelque peu les calculs, à mon avis, n'est clairement pas suffisant, j'ai donc décidé d'écrire comment cela devrait être fait en première approximation. Je connais parfaitement la recette universelle pour résoudre les problèmes de vitesse - prendre MK avec une fréquence plus élevée ou passer à une autre famille, mais je viens de l'époque où nous avons appris à nous débrouiller avec ce que nous avons, simplement parce qu'il n'y avait rien d'autre, du mot du tout. À l'heure actuelle, l'approche est dépassée, mais il me semblait qu'elle ne serait pas sans intérêt pour les lecteurs modernes de Habr.

Nous formulons le problème - nous voulons calculer le plus rapidement possible les coordonnées des points de la courbe de Bézier définis par les points extrêmes A et B et le foyer imaginaire C. La formule de calcul du point P sur la courbe est donnée par

$$

$$(1) P = T * T * B + 2 * T * (1-T) * C + (1-T) * (1-T) * A$$

où T varie de 0 à 1 inclus. (Sur Wiki, ils écrivent que cette formule était secrète à un moment donné, c'était étrange comme ça, mais tout est possible). Il est clair que nous ne le prendrons pas sous une forme complexe, au lieu de cela, nous rechercherons séparément les coordonnées X et Y. Nous allons estimer la complexité du calcul en utilisant cette formule, simplement en comptant le nombre de signes d'opérations arithmétiques dans cette expression - 7 multiplications et 5 additions (=> 7 * 5 + ) Il est possible qu'un bon compilateur (et maintenant tous les compilateurs sont bons et optimisent parfaitement si vous ne les interdisez pas explicitement) réduira les coûts à 7 * 3 +, bien qu'il serait préférable de l'aider en calculant (1-T) à l'avance. De manière générale, un bon compilateur peut généralement faire des merveilles si toutes les valeurs de la formule sont représentées par des constantes, mais nous supposons que toutes les valeurs sont statiquement indéfinies.

Première partie, mathématiques

Nous commençons le processus d'optimisation, pour lequel nous développons les crochets et regroupons les termes en T (peut-être qu'un jour le compilateur pourra le faire pour nous, mais jusqu'à présent cette partie du travail est affectée à l'intelligence naturelle), obtenant

$$

$$(2) P = T * T * (B + A-2 * C) + T * 2 * (C-A) + A$$

=> 5 * 5 +, ce qui est clairement meilleur que la valeur initiale de 7 * 5 +, mais une amélioration relative de 7 * 3 + devrait encore être envisagée.

Si nous prenons le temps d'exécution de l'opération d'addition comme un, le temps de multiplication ne sera exactement pas inférieur à un, en règle générale, plus, mais combien cela dépend de la mise en œuvre de MK (j'ai écrit d'abord sur l'architecture, mais ce n'est pas entièrement vrai). Lorsqu'il n'y a pas de multiplicateur matériel sur le cristal, le temps d'exécution de la multiplication sera dix (30+) fois supérieur à un, et lorsqu'il est présent, il sera plusieurs fois (1-6). Par conséquent, nous pouvons croire avec confiance que le remplacement de la multiplication par l'addition donne presque toujours un gain (et souvent significatif) de temps d'exécution. Eh bien, nous remarquerons immédiatement que la transition des nombres à virgule fixe à un point flottant (nous laissons de côté la preuve de ce fait) conduit à une augmentation du temps d'exécution de plus de 20 fois pour l'addition (l'alignement est très influent ici), mais seulement à une légère augmentation pour la multiplication . Par conséquent, pour les nombres à virgule flottante, les temps d'addition et de multiplication diffèrent peu, en particulier en termes relatifs (nous pouvons nous attendre à un maximum de 2 fois), mais ils diffèrent toujours et ne sont pas en faveur de la multiplication, donc il y a un gain ici.

En revenant au paragraphe précédent, nous constatons que pour PT, la note 5 * 5 + ne devrait pas avoir un avantage significatif sur 7 * 3 +, mais nous avons encore des réserves. Faites attention au fait que nous devons calculer l'ensemble des points sur la courbe de Bézier lorsque le paramètre T change, et tous les autres paramètres de la courbe sont fixes (mais pas constants, mais désolé), alors le reste de la formule peut être calculé à l'avance et obtenir

$$

$$(3) P = T * T * A1 + T * B1 + A$$

=> 3 * 2 +, où $$$A1 = A + B-2 * C$et $$$B1 = 2 * (C-A)$déjà bon, mais si vous vous souvenez du plan de Horner et écrivez

$$

$$(4) P = T * (T * A1 + B1) + A$$

=> 2 * 2 +, alors par rapport à la décision «sur le front» nous devons gagner plus de 2 fois, presque 3, et ces optimisations sont tout à fait évidentes.

Vérifions la théorie avec la pratique (bien que cela soit complètement redondant, nous sommes confiants dans nos estimations, mais soudain j'ai sous-estimé le compilateur), pour lequel nous devons mesurer le temps réel d'exécution de différentes options sur du matériel réel. Eh bien, il se trouve que chez moi, j'ai beaucoup de toutes sortes de cartes de débogage pour MK de diverses sociétés (y compris des raretés comme les débogages de Luminary Micro ou Intel Edisson, essayez d'en acheter une maintenant), mais il n'y a pas une seule carte Arduino («Eh bien nous n'avons pas d'ananas »). Cela semblerait être une impasse, mais il existe des options - un site très intéressant tinkercad.com vient à notre aide, sur lequel vous pouvez construire votre circuit sur une planche à pain en utilisant le module Arduino, écrire un croquis et l'exécuter immédiatement. Dans le même temps, vous pouvez définir des points d'arrêt, exécuter le programme étape par étape et même (une chose sans précédent pour un vrai Arduino) afficher les valeurs des variables aux points d'arrêt.

Nous nous tournons vers ce site et commençons à mesurer. Pour commencer, nous vérifions nos hypothèses sur le temps d'exécution des opérations et, après avoir éliminé les circonstances environnantes, nous obtenons les données suivantes pour les nombres entiers:

8 + 8 => 8 - 1 temps, 16 + 16 => 16 - 2,
8 * 8 => 16 - 2, 16 * 16 => 16 - 14 (la seule chose qui s'est avérée inattendue, j'ai pensé obtenir 4 * 2 + 4 * 2 = 16, il y a des optimisations intéressantes),
8/8 => 8 - 230, 16/16 => 16 - 230.

Faites attention aux deux derniers chiffres, il est clair d'après eux que l'opération de division est interdite si nous voulons vraiment compter rapidement. Maintenant (enfin) nous mesurons le temps nécessaire pour effectuer des opérations sur le nombre de PT avec une mantisse de 24 bits
a + b - 126 (et dépend fortement des opérandes), a * b - 140, a / b - 482.
Les données obtenues correspondent bien à nos hypothèses théoriques, il est clair qu'il y a une implémentation matérielle à bord de ce MK: pour la multiplication, pour la division, pas pour les opérations, PT.

Maintenant, nous commençons à mesurer le temps de calcul complet. Nous fixons les valeurs A = 140, B = 120, C = 70 et construisons 170 points uniformément répartis sur le bureau d'études. Pourquoi précisément ces valeurs - elles ont été données dans le message spécifié lors de l'évaluation des performances. Voici les algorithmes et le temps d'exécution du test correspondant.

Formule (1) => 20 ms ou 1 900 cycles d'horloge par échantillon
Formule (1) => 18 ms ou 1660 cycles d'horloge par échantillon (considérer séparément 1-T)
Formule (2) => 16 ms ou 1540 cycles d'horloge par échantillon
Formule (3) => 10 ms ou 923 cycles d'horloge par échantillon
Formule (4) => 8 ms ou 762 mesures par comptage

On peut voir que la réduction du temps d'exécution qui en résulte (de 20 ms à 8 ms) correspond bien à celle attendue et nous avons pu accélérer les calculs de plus de 2 fois. A noter qu'en plus de considérations tout à fait évidentes et mathématiques, ne dépassant pas le cours du lycée, nous n'en avions pas besoin.

Et maintenant, parlons de quoi faire si le résultat n'est pas suffisant, et nous avons déjà tout éliminé des formules de calcul. Ils m'ont écrit ici (dans les commentaires d'un autre article) qu'en général tout problème peut être réduit à l'informatique avec FT et, malgré la controverse évidente de l'hypothèse (essayez de le faire pour la solution numérique des équations de Navier-Stokes), dans ce cas particulier cette recommandation est applicable Bien que, comme toujours, il existe des nuances.

Deuxième partie, Informatique

Une fois les modifications de l'algorithme épuisées, seules les structures de données restent et nous entrons dans le sol des nombres à virgule fixe. Ici, nous trouverons de nombreux pièges auxquels nous n'avons pas pensé pour la plage et la précision du PT (en général, pour le PT, il faut penser à ces problèmes, mais ici tout est plus simple, beaucoup a été fait pour nous). Il est nécessaire de mener une petite étude du problème pour déterminer la représentation nécessaire de FT (sélectionnée dans le poste 9.7 susmentionné, à en juger par les résultats, elle est clairement insuffisante), mais je propose de prendre un chemin légèrement différent. Soit dit en passant, si nous ne faisons pas 170 pas sur l'intervalle, mais 128 (je ne vois aucune raison de nous interdire cette étape), cette idée nous conviendrait parfaitement. Si nous prenons en compte le fait que les constantes définissant la KB sont données par des entiers, et que le seul paramètre T peut être représenté par une fraction de la forme et / et nous utiliserons le résultat pour le rendu à l'écran, c'est-à-dire qu'il se traduira en coordonnées entières, alors nous pouvons il suffit de tout faire en nombres entiers, qui traitent beaucoup plus rapidement.

Nous utilisons uniquement la dernière formule et la réécrivons dans la nouvelle notation

$$

$$(5) P = u / U * (u / U * A1 + B1) + A$$

(=> 2 * 2 + 2 /), où A1 et B1 sont calculés de la même manière que pour PT. De toute évidence, tous les nombres sont des nombres entiers et les opérations correspondantes doivent être effectuées beaucoup plus rapidement. Afin de ne pas perdre de précision lors de l'opération de division entière (2/3 = 1! = 1.5) et de faire la division au tout dernier moment, on transforme légèrement la formule en la forme

$$

$$(6) P = ((et * A1 + B1 * Et) / Et * Et + A * Et) / Et$$

(=> 4 * 2 + 2 /). Tous les nombres FT, donc nous implémentons cet algorithme et obtenons ... vous voici, grand-mère et le jour de Yuryev ... 1869 cycles, mais c'est bien pire que pour FT, nous sommes partis de cela, une sorte de poubelle, car les entiers sont beaucoup plus rapides.

Nous commençons le débriefing et il s'avère que le simple changement de type de variables ne suffit pas. Premièrement, nous devons utiliser des nombres non pas 8 ou même 16, mais 32 bits, sinon un débordement se produira et des nombres longs, bien que plus rapides que PT, mais pas autant que pour compenser les défauts de l'algorithme. Deuxièmement, ces défauts sont dans nous avons de nouveau calculé des constantes sur chaque mesure - nous les supprimons par un calcul préliminaire B2 = B1 * I, A2 = A * I * I. Ensuite, nous obtenons

$$

$$(7) P = ((et * A1 + B2) * et + A2) / Et / Et$$

(=> 2 * 2 + 2 /) avec un résultat de 1684 est meilleur que le précédent, mais nous ne nous en sommes pas encore éloignés.

On exclut le calcul d'une autre constante And2 = Et * Et on obtient

$$

$$(8) P = ((et * A1 + B2) * et + A2) / II$$

(=> 2 * 2 + 1 /), avec un temps d'exécution de 956 cycles - mais c'est intéressant, l'exclusion d'une opération a conduit à une augmentation significative de la productivité.

C'est ce qui nous ralentit - la division, car c'est une opération qui prend beaucoup de temps, mais nous avons une astuce intéressante pour y faire face. Pour calculer l'expression 1 / Et nous pouvons effectuer des transformations élémentaires 1 / = 1 / * ( / ) = 1 * ( / ) / . Si nous choisissons le degré de deux comme H, alors la division par H peut être remplacée par des décalages, et si l'exposant est un multiple de 8, alors même des décalages ne seront pas nécessaires. Et la valeur de N / A devra être calculée honnêtement, mais une seule fois, après quoi seule la multiplication restera dans le cycle de calcul.

Faites attention au fait que nous avons fait une conversion pas tout à fait correcte et remplacé le N / A par sa valeur arrondie K pour passer aux opérations exclusivement avec des entiers. L'inexactitude consiste dans la perte de précision et des recherches supplémentaires doivent être effectuées pour prouver l'applicabilité de cette approche à notre cas. Nous écrivons H / I sous la forme (K * I + d) / I = K + (d / I), où q est inférieur à I. Ensuite, l'erreur absolue en allant de H / I à K sera d / I, et l'erreur relative sera d / I I / (K + d / I)> = d / I / (K + 1) ~ d / I / K, à condition que K >> 1 (ce n'est pas un décalage). Il s'ensuit que la valeur de H doit être choisie aussi grande que possible, car l'erreur de calcul absolue est égale à A * d / I / K> = A * 1 / N / I. Si nous voulons que l'erreur ne soit pas supérieure à l'unité, nous devons résister à la condition A / K <= 1, puis K> = A, nous convertissons K * I> = A * I, ce qui signifie H> = A * I, alors nous ne le faisons pas perdre en précision. Dans notre cas, A <= 256 et I <= 256, nous obtenons H> = 2 ** 16, ce qui est tout à fait acceptable. Évidemment, dans les formules ci-dessus, les modules des nombres originaux doivent être utilisés.

Nous notons pour l'avenir que si nous arrondissons non pas vers le bas, mais vers l'entier le plus proche, alors les exigences sont quelque peu réduites et H devrait être suffisamment moitié moins, bien qu'il y ait des nuances.

Dans tous les cas, nous pouvons fournir la précision requise et obtenir l'algorithme suivant: H = 2 ** 16; K = [N / A] (I <256); 0 <= et <= AND;

$$

$$(9) P = ((((et * A1 + B2) * et + A2) * K) >> 16) * K) >> 16$$

(=> 4 * 2 + 2 >> 16) où toutes les opérations sont effectuées sur des entiers longs. Nous implémentons cet algorithme et obtenons 583 cycles d'horloge ... mais cela est déjà proche de l'idéal, mais pas encore idéal.

Viennent ensuite les petits paramètres pour un MK spécifique - travailler avec des variables globales est plus rapide. qu'avec les locaux, mais encore plus rapide avec les registres locaux, ce qui entraîne une réduction du temps à 506 cycles d'horloge.

De plus, nous notons que la dernière multiplication avant le décalage peut être effectuée avec des nombres de 16 bits, ce qui donnera 504 - une bagatelle, mais agréable.

Au total, nous avons accéléré les calculs par rapport à l'implémentation «front» en 1900/504 - plus de 3 fois, et nous n'avons pas du tout perdu le mot. C'est le résultat que j'appelle l'optimisation du temps, et non 20% reçu dans le message d'origine.

Est-il possible d'atteindre des indicateurs encore meilleurs - c'est possible, mais c'est le sujet du prochain billet.