J'utilise une bibliothèque de points fixe de Anthony Williams décrite dans l'article du Dr Dobb " optimisant les applications à forte intensité de maths avec arithmétique à point fixe " pour calculer la distance entre deux points géographiques à l'aide du Méthode de ligne RHumb . Cela fonctionne assez bien lorsque la distance entre les points est significative (supérieure à quelques kilomètres), mais est très pauvre à des distances plus petites. Le pire des cas étant lorsque les deux points sont égaux ou presque égaux, le résultat est de 194 mètres, tandis que j'ai besoin de précision d'au moins 1 mètre sur des distances> = 1 mètre. par comparaison avec un La mise en oeuvre de point flottant à double précision, j'ai localisé le problème à la fonction Corriger le résultat pour Le Notez que Utilisation de la La méthode de la ligne Rhumb est une recommandation corporelle de normes pour cette application afin de ne peut donc pas être modifiée et, dans tous les cas, une fonction de racine carrée plus précise est susceptible d'être requise ailleurs dans l'application ou dans les applications futures. conclusion
À la lumière de la solution et de l'analyse de Justin Peel, et de la comparaison avec l'algorithme dans "l'art négligé de l'arithmétique de point fixe" , j'ai adapté ce dernier comme suit: tandis que cela donne une précision beaucoup plus grande, l'amélioration J'avais besoin de ne pas être atteint. Le format Q36.28 seul à propos de la précision dont j'ai besoin, mais il n'est pas possible d'effectuer une SQRT () sans perte de quelques bits de précision. Cependant, une pensée latérale offre une meilleure solution. Mon application teste la distance calculée contre une limite de distance. La solution plutôt évidente à l'arrière est de tester le carré de la distance contre le carré de la limite! corrigé :: sqrt () , qui fonctionne mal à de petites valeurs: fixe :: sqrt (0) est trivial en le traitant comme un cas particulier, mais cela ne résoudra pas le problème des petites distances non nulles, où l'erreur commence à 194 mètres. et converge vers zéro avec une distance croissante. J'ai probablement besoin d'au moins un ordre d'amélioration de la magintude de la précision envers zéro. corrigé :: sqrt () algorithim est brièvement expliqué à la page 4 de l'article lié ci-dessus, mais Je me lance de le suivre, laissez-la seule déterminer s'il est possible de l'améliorer. Le code de la fonction est reproduit ci-dessous: m_nval est la valeur de représentation du point fixe interne, c'est un int.4_t < / Code> Et la représentation utilise format Q36.28 ( fixe_resolution_shift < / code> = 28). La représentation elle-même a suffisamment de précision pour au moins 8 décimales et une fraction de l'arc équatorial est bonne pour des distances d'environ 0,14 mètre. La limitation n'est donc pas la représentation à point fixe. fixe :: sqrt () algorithme pour de petites valeurs non nulles tout en conservant sa convergence délimitée et déterministe? < p>
4 Réponses :
Il y a plusieurs années, j'ai travaillé sur un programme de démonstration pour un petit ordinateur que notre tenue avait construite. L'ordinateur disposait d'une instruction carrée intégrée et nous avons construit un programme simple pour démontrer l'ordinateur effectuant 16 bits ajouter / soustraire / multiplier / division / racine carrée sur un TTY. Hélas, il s'est avéré qu'il y avait un bogue grave dans l'instruction racine carrée, mais nous avions promis de démo la fonction. Nous avons donc créé une gamme de carrés des valeurs 1-255, puis a utilisé une recherche simple pour correspondre à la valeur saisie de l'une des valeurs de matrice. L'index était la racine carrée.
Malheureusement, j'ai besoin de meilleure précision sur une plage plus large qui serait réalisable avec une recherche.
étant donné que IE P> sqrt (AB) = sqrt (a) sqrt (a) sqrt (b) code>, alors ne pouvez pas simplement piéger le cas où votre numéro est petit et la déplace par un nombre donné de bits , calculez la racine et décalage de la moitié de la moitié du nombre de bits pour obtenir le résultat? sqrt(n) = sqrt(n.2^k)/sqrt(2^k)
= sqrt(n.2^k).2^(-k/2)
Solution très intelligente et efficace.
Je ne sais pas comment vous obtenez les chiffres de Voici ce que je fais: p> corrigé :: sqrt () code> montré dans la table. 0.00000000: sqrtfix0()=0.00003052 / err=3.051758e-05 sqrt()=0.00000000 sqrtfix2()=0.00000000 / err=0.000000e+00
0.00001000: sqrtfix0()=0.00311279 / err=4.948469e-05 sqrt()=0.00316228 sqrtfix2()=0.00316207 / err=2.102766e-07
0.00002000: sqrtfix0()=0.00445557 / err=1.656955e-05 sqrt()=0.00447214 sqrtfix2()=0.00447184 / err=2.974807e-07
0.00003000: sqrtfix0()=0.00543213 / err=4.509667e-05 sqrt()=0.00547723 sqrtfix2()=0.00547720 / err=2.438148e-08
0.00004000: sqrtfix0()=0.00628662 / err=3.793423e-05 sqrt()=0.00632456 sqrtfix2()=0.00632443 / err=1.262553e-07
0.00005000: sqrtfix0()=0.00701904 / err=5.202484e-05 sqrt()=0.00707107 sqrtfix2()=0.00707086 / err=2.060551e-07
0.00006000: sqrtfix0()=0.00772095 / err=2.501943e-05 sqrt()=0.00774597 sqrtfix2()=0.00774593 / err=3.390476e-08
0.00007000: sqrtfix0()=0.00836182 / err=4.783859e-06 sqrt()=0.00836660 sqrtfix2()=0.00836649 / err=1.086198e-07
0.00008000: sqrtfix0()=0.00894165 / err=2.621519e-06 sqrt()=0.00894427 sqrtfix2()=0.00894409 / err=1.777289e-07
2.00000000: sqrtfix0()=1.41419983 / err=1.373327e-05 sqrt()=1.41421356 sqrtfix2()=1.41421356 / err=1.851493e-09
34359700000.00000000: sqrtfix0()=185363.69654846 / err=5.097361e-06 sqrt()=185363.69655356 sqrtfix2()=185363.69655356 / err=1
.164153e-09
La mise en œuvre initiale a évidemment quelques problèmes. Je suis devenu frustré par la tentative de les réparer avec la manière dont le code est actuellement fait et fini par y aller avec une approche différente. Je pourrais probablement réparer l'original maintenant, mais j'aime bien mon chemin de toute façon.
Je traite le numéro d'entrée comme étant dans Q64 pour démarrer, ce qui est identique à partir de 28, puis sur le dos de 14 ensuite (le SQRT la refuse ). Cependant, si vous faites cela, la précision est limitée à 1/2 ^ 14 = 6.1035e-5 car les 14 derniers bits seront 0. Pour remédier à cela, je transfère ensuite Au fait, une grosse erreur dans le code d'origine est que le terme où m = 0 n'est jamais considéré comme si le bit ne peut jamais être réglé. Quoi qu'il en soit, voici le code. Profitez! avec la sortie du programme étant A et < Code> reste correctement et pour continuer à remplir les chiffres, je fais la boucle à nouveau. Le code peut être rendu plus efficace et plus propre, mais je laisserai cela à quelqu'un d'autre. La précision indiquée ci-dessous est à peu près aussi bonne que possible avec le Q36.28. Si vous comparez le point fixe sqrt avec le point flottant SQRT du numéro d'entrée après avoir été tronqué par point fixe (le convertir en point fixe et en arrière), les erreurs sont autour de 2E-9 (je n'ai pas fait cela dans le code ci-dessous, mais il nécessite une ligne de changement). C'est juste en ligne avec la meilleure précision du Q36.28 qui est 1/2 ^ 28 = 3.7529e-9.
C'est exactement ce que j'ai demandé, et plus ou moins un remplaçant d'abandon du corps de code existant. Malheureusement, mon estimation de la précision requise était incorrecte et même avec cette vaste amélioration, elle est insuffisante. Cela améliore la précision dans laquelle SQRT () est utilisé ailleurs, alors je suis susceptible de le retenir. Je vais le regarder plus loin comme une diligence raisonnable, mais si vous dites que c'est la limitation de la performance, dans ce cas, je devrai utiliser le point STD: SQRT () et flottant.
J'ai comparé les résultats de cela avec l'algorithme dans " L'art négligé de l'arithmétique de point fixe " , et cela donne des résultats identiques et est probablement la version plus efficace / nettoyeuse que vous faisiez référence. Vous avez au moins m'a fait comprendre les limites de ce qui peut être atteint avec le Q36.28. Merci.
Clifford - L'URL à l'article a été remplie (demandant une connexion UBM Techweb, vous étiez probablement connecté et n'a pas fait l'expérience de cette expérience). J'ai essayé de trouver l'article ailleurs mais non - Google Cache semblait être la meilleure alternative. Merci de référencer l'article en premier lieu.
@Dan, je suis allé au lien d'origine et je n'ai pas eu de problème. Je n'ai jamais utilisé UBM et je ne pense pas que je suis connecté à DDJ, donc je ne sais pas pourquoi il y a un problème.
J'utilise le code pendant quelques années; J'ai téléchargé la bibliothèque probablement avant l'exigence de connexion. Obtenez-le de Anthony's propre site web
@MarkRansom - OK, je devais regarder dans ceci ... J'avais des biscuits de Dddobbs / UBM / DDJ, dès que je les ai effacés, cela m'a laissé passer sans me forcer à me connecter. Nice, donnez à vos utilisateurs enregistrés plus d'un obstacle.
Comment avez-vous obtenu les chiffres de
fixes :: sqrt () montré dans la table? Quel compilateur + OS avez-vous utilisé? Je ne reçois pas les mêmes chiffres pour quoi que ce soit, mais une racine carrée de 0. Ni avec GCC (DJGPP / DOS), ni avec Open WatCom (Windows). Mes résultats sont tous différents d'environ 10 ^ -5 à 10 ^ -6 pourx de la table, jamais 10 ^ -7 ou 10 ^ -9. Avez-vous utilisé plus de 28 morceaux fractionnaires lorsque vous remplissez la table? Comment avez-vous converti / depuisfixe et quelle est la taille de votre type de point flottant (BTW, était-ce double ou long double)?@ALEX: Le test a été effectué à l'aide de VC ++, mais des résultats identiques ont été obtenus à l'aide de Bras Realview. La conversion a été réalisée avec la fonction de membre fixe :: to_double (). Je vais ajouter le code de test à la question dès que je peux - il n'est pas disponible pour moi à présent. Tous les types de test sont fixes ou doubles. L'erreur a été calculée par
fabs (std :: sqrt (x) - sqrt (fixe (x)))@ALEX: Correction ... La conversion était avec
corrigé :: as_double () .as_double () exécutex / (double) 1 << fixe_resolution_shift