EXTRACTION D'UNE RACINE CARRÉE

BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019 - 4 août 2005 à 12:06
cs_Charolais Messages postés 1 Date d'inscription vendredi 20 octobre 2006 Statut Membre Dernière intervention 20 octobre 2006 - 20 oct. 2006 à 21:56
Cette discussion concerne un article du site. Pour la consulter dans son contexte d'origine, cliquez sur le lien ci-dessous.

https://codes-sources.commentcamarche.net/source/33073-extraction-d-une-racine-carree

cs_Charolais Messages postés 1 Date d'inscription vendredi 20 octobre 2006 Statut Membre Dernière intervention 20 octobre 2006
20 oct. 2006 à 21:56
En réponse à notre ami COSMOBOB, je me permet d'émettre 2 critiques sur l'algo :
- D'une part, le résultat d'une racine carré d'un 32 bits sera codé sur 16 bits (si ton "int" est bien un 32 bits, cale cela dépend des compilos)
- D'autre part, le test "milieucarre < milieu" , bien que corrigeant le résultat de calcul pour certaines valeurs ne corrige pas tout ! (quand on tend vers 2^32 ça tend aussi vers des résultats systématique faux !)
Il suffit de borner la variable "droite" avec la racine de la valeur max du type d'entée, qui se trouve être aussi la valeur max de sortie (on calcul une racine carré je rappel) donc avec 2^16 = 65536.
Qu'est-ce que vous en pensez-vous ?
Il y a peut-être encore plus optimum, mais au moins, je pense que cet algo fonctionne parfaitement cette fois-ci
Enfin, pour un sujet noté "Niveau de la source : Débutant", je trouve qu'il y a beaucoup de commentaires. C'est pas si "Débutant" que cela de faire un algo de racine carré, vu le nombre d'algos qui existent, qui soit rapide et fiable...


typedef unsigned int u16;
typedef unsigned long u32;

u16 mysqrt(u32 a)
{
u16 gauche 0, milieu 0;
u16 droite;
u32 milieucarre = 0;

if (a < 65536)
{
droite = (u16)value;
}
else
{
droite = 65535;
}

if (a <= 1)
return ((u16)a);
while ((droite - gauche) > 1)
{
milieu = (u16)(((u32)droite + gauche) / 2);
milieucarre = ((u32)milieu * milieu);
if (milieucarre > a)
{
droite = milieu;
}
else
{
gauche = milieu;
}
}
return (gauche);
}
cs_Kirua Messages postés 3006 Date d'inscription dimanche 14 avril 2002 Statut Membre Dernière intervention 31 décembre 2008
3 sept. 2005 à 03:24
Bon, ça me chagrinait cosmobob que tu dises que mon algo était lent, parce que je voyais pas trop pq, alors j'ai benchmarké avec le high frequency timer de windows, et le benchmark donne ton algo gagnant (4% plus rapide que le mien) sur 10 000 000 de naturels compris entre 0 et MAX_RAND (=2^16-1 chez moi). Je pense pas que tu aies pu voir une telle différence à l'oeil nu, mais en attendant t'as qd même gagné, bravo ^^.
cs_Kirua Messages postés 3006 Date d'inscription dimanche 14 avril 2002 Statut Membre Dernière intervention 31 décembre 2008
3 sept. 2005 à 03:10
arnaud, je l'ai trouvée en griffonant sur un bout de papier.
quant à "c'est une trouvaille ça", ben je trouve que oui, dans la mesure où je ne risque vraiment pas d'overflow, que je ne bosse qu'avec des entiers, que c'est élégant, que c'est du sûr, que c'est indirect... la recherche dichotomique je l'ai faite aussi, t'en fais pas, et c'est très bien si tu veux une approximation avec virgule flottante, mais c'était pas la question, comme j'ai précisé: fallait le naturel égal ou inférieur.
au passage, ma recherche dichotomique était environ 25 fois plus lente que la version de la librairie standard, qui, je suppose, fait appel au FPU d'Intel, donc bon ^^.
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
6 août 2005 à 09:35
Pour FSQRT tu vas demander à Intel, c'est en dur dans les transistors de la fpu, je ne suis pas chimiste.
Arnaud16022 Messages postés 1329 Date d'inscription vendredi 15 août 2003 Statut Membre Dernière intervention 16 juin 2010 2
6 août 2005 à 01:50
oh et pour Kirua : t'es pas bete toi ^^ tu la tires d'ou cette méthode?
Arnaud16022 Messages postés 1329 Date d'inscription vendredi 15 août 2003 Statut Membre Dernière intervention 16 juin 2010 2
6 août 2005 à 01:47
pour Vecchio : j'avais fait ca sur TI83, mais c'était en basic...lol
ca marchait tres bien, mais j'ai pas pu comparer question perfs avec la fonction "normale" (c'etait franchement lent...3-4 secondes)
JCDjcd -> mdr le gars qui se la pete :D tu peux preciser ?
Brunews ou tout autre: et le coprocesesseur, lui, il fait comment pour le faire, son FSQRT ?

++
ad
cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
5 août 2005 à 15:17
En dix itérations de l'algorithme de Newton (méthode des tangentes) on obtient une bonne approximation.
cosmobob Messages postés 700 Date d'inscription mardi 30 décembre 2003 Statut Membre Dernière intervention 27 janvier 2009 4
5 août 2005 à 12:28
et recorrigée sans avoir besoin des int 64:

int mysqrt(int a)
{
int gauche 0, droite a, milieu = 0, milieucarre = 0;
if (a <= 1)
return a;
while (droite - gauche > 1)
{
milieu = (droite + gauche) / 2;
milieucarre = milieu*milieu;
if (milieucarre < milieu || milieucarre > a)
{
droite = milieu;
}
else
{
gauche = milieu;
}
}
return gauche;
}
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
5 août 2005 à 12:24
Va bon ainsi, résultats exacts.
cosmobob Messages postés 700 Date d'inscription mardi 30 décembre 2003 Statut Membre Dernière intervention 27 janvier 2009 4
5 août 2005 à 12:18
version corrigée:

int mysqrt(int a)
{
int gauche = 0;
int droite = a;
if (a <= 1)
return a;
while (droite - gauche > 1)
{
int milieu = (droite + gauche) / 2;
__int64 milieucarre = (__int64)milieu*milieu;
if (milieucarre > a)
{
droite = milieu;
}
else
{
gauche = milieu;
}
}
return gauche;
}
cosmobob Messages postés 700 Date d'inscription mardi 30 décembre 2003 Statut Membre Dernière intervention 27 janvier 2009 4
5 août 2005 à 12:12
non, mais ya un overflow intermédiaire ;) (le milieu²)
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
5 août 2005 à 12:07
v = 99986984;
ultoa(RacineCosmo(v), buf, 10);
MessageBox(0, buf, "racine cosmo", 0);

j'obtiens 7639941 au lieu de 9999.
Aurais-je mal copier-coller ?
cosmobob Messages postés 700 Date d'inscription mardi 30 décembre 2003 Statut Membre Dernière intervention 27 janvier 2009 4
5 août 2005 à 11:50
a part pour a = 1 (ca rentre pas ds la boucle), oui
(t'as continuellement gauche² <= a et droite² > a, donc quand droite-gauche<=1, t'es assuré que gauche est le plus grand entier tel que gauche² <= a)
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
5 août 2005 à 11:22
tu es certain d'avoir des résultats corrects ?
cosmobob Messages postés 700 Date d'inscription mardi 30 décembre 2003 Statut Membre Dernière intervention 27 janvier 2009 4
5 août 2005 à 10:46
extraction d'une racine carré d'un entier en un temps en log(n) :

int mysqrt(int a)
{
int gauche = 0;
int droite = a;

while (droite - gauche > 1)
{
int milieu = (droite + gauche) / 2;
if (milieu*milieu > a)
{
droite = milieu;
}
else
{
gauche = milieu;
}
}
return gauche;
}
cosmobob Messages postés 700 Date d'inscription mardi 30 décembre 2003 Statut Membre Dernière intervention 27 janvier 2009 4
5 août 2005 à 10:17
salut,
tu parles d'une curiosité, découvrir que (n+1)²-n² = 2n+1 lol
par contre c'est pas tres tres rapide comme methode !!

sinon je signale que x->x²-a est une fonction croissante, une recherche dichotomique marchera donc tres bien !!!!
cs_Kirua Messages postés 3006 Date d'inscription dimanche 14 avril 2002 Statut Membre Dernière intervention 31 décembre 2008
5 août 2005 à 01:56
pour trouver la partie entière d'une racine carrée sans utiliser sqrt, j'avais écrit cet algo il y a qq mois, en découvrant une curiosité par hasard, je te tape le code avec le commentaire:

int Racine(unsigned int n)
{
//on évacue d'abord un cas gênant
if(n == 0) return 0;
//on peut énumérer tous les carrés
//naturels très simplement en ajoutant
//le prochain impair dans une liste:
//0
//0 + 1 = 1
//1 + 3 = 4
//4 + 5 = 9
//9 + 7 = 16
//...
//On va parcourir la liste des carrés naturels, et lorsqu'on a dépassé
//le nombre demandé, on sait quelle est sa racine carrée naturelle.
//Illustration: tous les nombres de 9 à 15 compris doivent retourner 3
//comme réponse, donc on peut passer tout de suite de 9 à 16. Si n est
//>= à 9 mais < 16, alors la réponse est 3 (et on retourne Ret).

unsigned int Ret = 0;
unsigned int Carre = 1;
unsigned int ListeImpairs = 1;

for(; Carre <= n; Ret++)
{
ListeImpairs += 2;
Carre += ListeImpairs;
}

return Ret;
}

j'avais écrit ça pour un exercice sur le site de prologin.
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 22:33
Si si, on est pas futé mais enfin je pense qu'on avait tout de même compris. Pour autant c'est les vacances pour pas mal de monde et donc on a un peu de temps pour pousser les choses un peu plus loin.
sqrt() c'est "plus con tu meurs" comme fonction (encore que on peut dire macro sur compilo correct car code mis inline):
FLD qword ptr[param]
FSQRT
tu vois, ils ne se sont pas arrachés un neurone sur ce coup.
Skaaar Messages postés 18 Date d'inscription jeudi 1 juillet 2004 Statut Membre Dernière intervention 21 avril 2006
4 août 2005 à 22:10
Je pense que beaucoup ici n'ont pas compris le but de la source...
Il s'agissait d'expliquer comment on peut calculer une racine carrée à la main, pas d'annalyser les perfs d'une méthode qui est évidement la plus efficace d'un point de vue algorythmique (sqrt()). Ceux qui ont concu ce language à mon avis ont passé quelques heures à concevoir cette fonction pour quelle soit la plus efficace possible.
A bon entendeur,
Ciao ++
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 20:36
oui c'est pas con comme idée, tu peux commencer par 1/4 valeur et suivant inf ou sup tu incr ou décrémentes sur le même principe. La dichotomie nouvelle forme...
vecchio56 Messages postés 6535 Date d'inscription lundi 16 décembre 2002 Statut Membre Dernière intervention 22 août 2010 14
4 août 2005 à 20:24
l'utilisation de sqrt n'est pas immédiate, tu a regardé la mise en oeuvre?
Je peux faire une boucle mais avec une incrémentation variable (je monte de plus en plus, dès que je dépasse je descend un peu puis de plus en plus, et je finis par converger assez vite, mais ca j'arrive pas trop à le faire)
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 20:13
ben là que commence le problème du choix:
- soit sqrt() va bon.
- soit tu fais par boucle jusqu'à trouver, portable mais ne sert plus à rien car rédhibitoire question perf.
vecchio56 Messages postés 6535 Date d'inscription lundi 16 décembre 2002 Statut Membre Dernière intervention 22 août 2010 14
4 août 2005 à 19:44
C'étati une classe protable au debut, si je commence a mettre de l'asm dedans c'est pas bon
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 19:40
Va voir les RSQRTxx dans manuels Intel (P3 mini).
1er calcul: x0 = RSQRTSS(a)
résultat corrigé: x1 = 0.5 * x0 * (3 - (a * x0)) * x0)
vecchio56 Messages postés 6535 Date d'inscription lundi 16 décembre 2002 Statut Membre Dernière intervention 22 août 2010 14
4 août 2005 à 19:17
Justement je cherche a calculer rapidement la partie entiere de la racine carrée d'un entier, mais j'y arrive pas
Pour cette source: http://www.cppfrance.com/code.aspx?ID=19105
Si quelqu'un a une piste
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 16:31
OUPS suis allé trop vite:
enlever tmps += bnTicksResult(); de chaque boucle do.
et placer tmps = bnTicksResult(); sous chaque boucle.
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 15:57
y a même quasi rien à gagner en asm, les multiplications et divisions en boucle plombent les perfs.
En tout cas c'est sympa comme algo, tu peux nous en sortir d'autres.
cs_Tripo Messages postés 5 Date d'inscription mardi 11 mai 2004 Statut Membre Dernière intervention 4 août 2005
4 août 2005 à 15:16
wep, suffit de rajouter un test pour vérifier que le chiffre en entré est supérieur ou égal à 0.
Pour la faire en asm, oui on y gagnerais en vitesse mais comme je le répète, ce n'est pas le but...
skone007 Messages postés 166 Date d'inscription mercredi 24 avril 2002 Statut Membre Dernière intervention 23 juin 2009
4 août 2005 à 15:02
Ouaih en effet c'est plus rapide mais si on fait racine de -1 ca donne quoi faudrait gerer ce cas et c'est parfais et même le faire en ASM ca serai encore plus rapide ...
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 14:07
arf, j'ai oublié à mettre en haut:
unsigned __int64 tmps, tmps1;
pour les temps de calcul.
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 14:05
Alors après tests, sqrt() est énormément + rapide.
J'ai enlevé les String(eries) pour essayer d'accélérer mais rien n'y fait.
bnPrecis.dll dans mes sources pour bnTicksStart() et bnTicksResult().

void racine(__int64 nbr, int pres, char *psz)
{
char *c;
int precision = 0;
__int64 res, tmp, courant;
// compteurs
int cpt, tmpVar, compteur=0;
precision = pres;
courant = (__int64) premiers(nbr);
nbr = apresPremiers(nbr);
cpt = carreInferieur((int) courant);
res = cpt;
// stockage pour l'affichage
//var.getIntValue(cpt);
//renvoie = renvoie + var;
ultoa(cpt, psz, 10);
c = psz;
while(*c) c++;
courant = courant - (cpt*cpt);
while(nbr != 0) {
cpt = premiers(nbr);
courant = (courant*100) + cpt;
nbr = apresPremiers(nbr);
tmp = res*2;
tmpVar = (int) ((courant/10)/tmp);
while(courant < (((tmp*10)+tmpVar) * tmpVar) ) tmpVar--;
res = (res*10) + tmpVar;
//var.getIntValue(tmpVar);
//renvoie = renvoie + var;
ultoa(tmpVar, c, 10);
while(*c) c++;
courant -= ((tmp*10)+tmpVar) * tmpVar;
}
if(precision != 0) {
*c++ ','; *c 0;
while(compteur != precision) {
courant *= 100;
tmp = res*2;
tmpVar = (int) ((courant/10)/tmp);
while(courant < (((tmp*10)+tmpVar) * tmpVar) ) tmpVar--;
res = (res*10) + tmpVar;
//var.getIntValue(tmpVar);
//renvoie = renvoie + var;
ultoa(tmpVar, c, 10);
while(*c) c++;
courant -= ((tmp*10)+tmpVar) * tmpVar;
compteur++;
}
}
}

int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE x, PSTR y, int z)
{
char buf[40];
int n, i;
tmps = 0;
tmps1 = 0;
i = 1000;
bnTicksStart();
do {
racine(99986984, 0, buf);
tmps += bnTicksResult();
} while(--i);
MessageBox(0, buf, szappname, 0);

i = 1000;
bnTicksStart();
do {
n = (int) sqrt(99986984);
itoa(n, buf,10);
tmps1 += bnTicksResult();
} while(--i);
MessageBox(0, buf, szappname, 0);

_ui64toa(tmps, buf, 10);
MessageBox(0, buf, szappname, 0);
_ui64toa(tmps1, buf, 10);
MessageBox(0, buf, szappname, 0);

return 0;
}
asmanur Messages postés 230 Date d'inscription mercredi 11 février 2004 Statut Membre Dernière intervention 4 août 2005
4 août 2005 à 12:58
Je me suis ptet planté dans le test lol je fais rarement ce genre de truc. J'ai fait calculer 10 racines chacun et apparrament sqrt a un peu de mal à suivre.
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 12:55
ah je n'aurais pas cru, je vais tester sur 32 et 64 bits.
asmanur Messages postés 230 Date d'inscription mercredi 11 février 2004 Statut Membre Dernière intervention 4 août 2005
4 août 2005 à 12:41
en effet c'est plus rapide que sqrt
cs_Tripo Messages postés 5 Date d'inscription mardi 11 mai 2004 Statut Membre Dernière intervention 4 août 2005
4 août 2005 à 12:40
Je ne sais pas s'il est plus rapide...
Le but n'était pas d'améliorer mais de montrer comment extraire une racine "à la main".

Pour ce qui est du String retourner, c'est juste une question d'affichage, le programme pet être modifier pour retourner autre chose, je souhaitais juste montrer le traitement, pas faire un programme tip top.
asmanur Messages postés 230 Date d'inscription mercredi 11 février 2004 Statut Membre Dernière intervention 4 août 2005
4 août 2005 à 12:35
ton programme retourne une string ce qui est totalement inutile le programmeur en a peu faire de voir une virgule a la place d'un point car il utilise le nombre.
Quant à savoir si c'est plus rapide que sqrt je teste
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 12:13
Tu es certain que c'est plus performant que sqrt() ???
BruNews Messages postés 21040 Date d'inscription jeudi 23 janvier 2003 Statut Modérateur Dernière intervention 21 août 2019
4 août 2005 à 12:06
SVP, éviter les accents dans noms de fichiers et dossiers, ça empêche la visualisation des fichiers dans la liste.
Rejoignez-nous