Un peu de maths ? calcul d'une racine carree - Page 1

#define EPSILON 0.000000000000000001 double m_abs(double x) {   if(x < 0)     return (-x);   else     return x; } double racine_carre(double v) {   double r, s;   s = v;   r = (v + 1) / 2;   while(EPSILON < m_abs(s - r))     {       s = r;       r = 0.5*(r + v/r);     }   return r; }

voila pour l'instant ma fonction de calcul d'une racine.

bon le gros du pb c'est que je voudrais obtenir la meme precision que la fonction sqrt de la lib maths c ...
quelqu'un aurait-il une idee pour ameliorer la precision car la malgre toutes mes tentatives je reste encore trop loin du resultat donne par sqrt

Tiens, on venait d'en parler là:
topics/30469-pedrom-080/3#61

hum ca parle pas trop de la racine

Tu peux encore faire des tests multiplicatifs une fois que ta racine est bien approchée, par exemple, si aprx^2 - sqrt^2 = delta, alors sqrt # aprx - 0.5*delta/aprx (tu peux itérer racine+=(valeur-racine^2)/racine plusieurs fois). Je ne garantis pas que ça soit vraiment plus précis, mais c'est à essayer (et la division n'a plus besoin d'être aussi précise avec cette méthode).

Sinon fais attention, ta routine va être très lente si x n'est pas de l'ordre de 1... (remarque, si c'est pour TI, de toute façon elle va être très lente )

Bah je sais pas trop si tu cherches la vitesse, mais si c'est pas super important, il suffit d'itérer ta boucle jusqu'à ce que ça soit constant.

Je crois que là t'obtiens un truc très proche de sqrt (enfin égal tout le temps, ou presque). Maintenant c'est vrai que c'est lent. (Je crois que c'est 3 fois plus lent que sqrt )

Il existe une méthode d'extraction des racines carrées à la main, en base 10 : à chaque itération, on a un chiffre en plus derrière la virgule (si ce n'est pas un carré exact).
Je me suis amusé l'année dernière à l'adapter en base 2, où les calculs sont bien plus simples, directement écrit en Asm x86.
Si ça t'intéresse et si je retrouve ma vieille feuille de brouillon, je t'envoie mon code, et si je ne retrouve pas la feuille, je pense pouvoir le refaire en moins d'une heure (le plus dur étant de retrouver la méthode originale en base 10 ...).

@++

L'algorithme en base 10, il est là.

pas con le test d'egalite seulement qd je le fait ca ne donne qd meme pas la valeur de sqrt ...
y'a pas moyen de connaitre l'algo exact utilise pour cette fct ?

C'est quoi, exactement, le pb? La fonction que t'as mis, c'est celle de la libc ou la tienne? ("voila pour l'instant _ma_ fonction de calcul d'une racine")

c la mienne

or on me demande d'atteindre la meme precision que sqrt, chose que je ne parvient pas a avoir

Bah tu dois pas en être très loin, non? C'est à cause des erreurs d'arrondi au cours de la division? Ou alors c'est à cause de ton EPSILON qui est trop grand?

[pre] #include <stdio.h> #include <stdlib.h> #include <math.h> float my_sqrt(float a) { float n = a; float old = -1; while(old != n) { old = n; n = 0.5 * (n + (a/n)); } return n; } int main() { for(n= 1000; n<10000000; n++) { if((float)my_sqrt((float)n) != (float)sqrt((float)n)) printf("%d - ", n); } return 0; } [pre] Chez moi j'obtiens aucun affichage, donc ça veut dire qu'il n'y a pas de différences. Enfin en enlevant le cast en float dans le test d'égalité, j'obtiens pas pareil par contre.

ben le pb c que le epsilon meme en le modifiant la precision ne change pas ...

je crois que g compris, il arrete la comparaison a la limite des floats g l'impression

v bidouiller avec une union

Personnellement, je ne vois absolument pas où un arrondi pourrait intervenir...
Tu peux pas essayer de débugger en lui donnant en pâture un truc du style "v = 1+16*EPSILON" et voir où ça foire?

1.4142135623730949234300169337075203657150268554687500
1.4142135623730951454746218587388284504413604736328125
0.0000000000000002220446049250313080847263336181640625


voila ce que ca me donne pour racine de 2

la premiere ligne c mon resultat
la seconde c sqrt
la derniere c la difference entre les 2

et bon la on voit tres bien ou ca s'arrete, et on arrive pas a passer cette limite ...

Oui, mais ça doit vouloir dire que quelque part, il y a un arrondi qui se fait : avec un débuggeur pas à pas, on devrait pouvoir trouver où... (mais ça doit pas être très visible avec racine(2), donc c'est pour ça qu'avec racine(1+16*EPSILON) on devrait pouvoir voir ça tout de suite puisque en arrondissant vers un float, on passe à 1.0 tout rond)

Fin: ; le resultat est dans bl ; c'est une valeur entiere arrondie par defaut Fin_Programme: mov ah, 0x4C int 0x21 section .data section .bss

Bon, voilà ce que ça donne en Asm x86 (je viens de le refaire, pas envie de rentrer chez moi pour rechercher une feuille que je ne retrouverai sans doute jamais) :	org	0x0100


section	.text

Debut:
	; on suppose que la valeur entiere dont on veut la racine est dans ax
	; cette valeur est strictement inferieure a 65536
	; donc le resultat est strictement inférieur a 256, il tient dans bl

	xor	bx, bx			; test de nullité
	test	ax, ax
	jz	Fin

	mov	cl, 7			; initialisation
	mov	dx, ax

Taille:
	rcl	dx, 1			; taille du resultat
	jc	Taille_OK
	rcl	dx, 1
	jc	Taille_OK
	dec	cl
	jmp	Taille

Taille_OK:
	inc	bx			; preparation premiere iteration
	mov	dx, ax
	ror	dx, cl
	ror	dx, cl
	dec	dx
	rcl	ax, 3
	rcl	dx, 1
	rcl	ax, 1
	rcl	dx, 1

Iteration:
	rol	bx, 2			; recherche de valeur de bit
	inc	bx
	cmp	bx, dx
	ja	Suite

	sub	dx, bx			; bit a 1
	inc	bx

Suite:
	ror	bx, 1			; preparation iteration suivante
	rcl	ax, 1
	rcl	dx, 1
	rcl	ax, 1
	rcl	dx, 1
	dec	cl
	jnz	Iteration

N'ayant pas NAsm ici, je ne peux pas l'assembler pour le tester sous Debug.com et vérifier que je ne me sois pas gourré quelque part ...
Bien sûr, on peux largement améliorer ça (continuer après la virgule, prendre la racine d'une valeur non entière), mais ce n'était qu'un premier jet.

@++

Je crois que Vark parle de racine carrée flottante, pas entière...

c cool mais je dois coder en C

mais qu'est-ce qui pourrait causer l'arrondis, passke la je nage en plein mystere

Moi aussi, et je pense que la meilleure façon de voir pourquoi est d'appliquer ce que j'ai dit au ./17. Ou alors tu peux regarder le code que génère ton compilo, si tu sais comment marche la programmation en ASM avec le FPU.

euh oue mais je c pas utiliser de debugger

bon g teste ma fonction marche parfaitement avec des float
visiblent la couille vient des double ...

Bon, j'ai compilé ton prog et voilà ce que m'affiche VS : 1.41421356237309490000
Donc en fait c'est tout simplement que les "double" ne vont pas plus loin que ça En fait l'erreur doit être simplement de DBL_EPSILON ou un truc du style, pas plus. C'est juste ton printf qui extrapole des chiffres au-delà de la précision maximale.

printf("%d \n", (int)(sqrt(marmotte) == racine_carre(marmotte)));

cette ligne me renvoit 0 avec des doubles ... dc ma valeur n'est pas la bonne ...

Donc tu veux _la_ valeur la plus proche, et tu ne peux pas tolérer une erreur de DBL_EPSILON? Fallait le dire...
Je pense que pour ce genre de choses, racine+=0.6*(valeur-racine^2)/racine (attention, 0.6, pas 0.5 ) devrait pas mal arranger ton histoire, et si tu as un doute tu peux encore ré-itérer. Mieux, tu peux voir si la valeur après cette itération est bien la même qu'avant, et sinon tu sauras que tu n'es pas sûr de ta valeur.

la g eu droit a une boucle infinie :/

argl, je comprends le problème. En l'occurrence, on a 2 - vark_sqrt(2)^2 == ansi_sqrt(2)^2 - 2, donc c pas vraiment facile de voir qui a raison Il faudrait essayer avec un deuxième float pour multiplier la précision par 2, je te tiens au courant de ce que je trouve.

Il suffit de taper ça dans la TI-89:
14142135623730949234300169337075203657150268554687500²=199999999999999964539536283296691068389630300640621776844330428901663765373086789622902870178222656250000
14142135623730951454746218587388284504413604736328125²=200000000000000027343234630647692806884916507957232351955020204815893780647684252471663057804107666015625
2*10^104-ans(2)=35460463716703308931610369699359378223155669571098336234626913210377097129821777343750000
2*10^104-ans(2)=-27343234630647692806884916507957232351955020204815893780647684252471663057804107666015625
Donc c'est la deuxième réponse la plus précise.