[dev-c++] calcul de la racine carrée d'un réel

Soyez le premier à donner votre avis sur cette source.

Snippet vu 8 526 fois - Téléchargée 19 fois

Contenu du snippet

Cette source permet de calculer une racine carrée par la méthode de Newton avec une approximation par division de l'exposant binaire faisant partie de la représentation en virgule flottante (ma réelle contribution).

Source / Exemple :


#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>

#define u_int unsigned int

#define iter 10

// Ce type de structure est peu employé mais parfait pour cette routine :
// Elle permet de déclarer deux variable de type/taille différente sur le même espace de
// mémoire.

union flint {
	float f;
	u_int i;
};

float rac(float x)
{
	float y1;
 	u_int b;
	flint ya;
	int i;

	if (x<=0) return 0;

//	on duplique le nombre avec un type adéquat pour les opérations binaires
	ya.f=x;

//	en tant que float un nombre est sous la forme m * 2 ^ e avec m < 1
//	la racine carrée = (m * 2 ^ e) ^ 1/2 et comme m < 1
//	m * 2 ^ (e/2) en est une approximation
//	
//	le type float est grossièrement codé ainsi (sur 32bits) :
//	0 01010101 11100011100000001110000
//	^ exposant      mantisse
//	'-signe

//	première opération :
//	déplace l'exposant (opérateur >>):
//		0exposant11100011100000001110000
//	=>	000000000000000000000000exposant
	b = ya.i>>23;

//	l'exposant est augmenté de 127 en mémoire pour permettre la représentation de nombre
//	dont l'exposant serait négatif.
//	il faut donc soustraire 127 avant de diviser par 2
//	b=(b-127)/2+127=(b+127)/2
//	=>	000000000000000000000000demiexpo
	b=(b+127)/2;

//	finalement on "réassemble" le nombre avec le nouvel exposant :

//	d'abord on annule les bit de l'ancien exposant
//		00101010111100011100000001110000 & 10000000011111111111111111111111
//		               x                            0x807fffff
//	=>	00000000011100011100000001110000

//	<< repositionne les bit au bon endroit
//	000000000000000000000000demiexpo => 0demiexpo00000000000000000000000

//	le | ou binaire permet "d'ajouter" notre exposant
//		00000000011100011100000001110000 | 0demiexpo00000000000000000000000
//	=>	0demiexpo11100011100000001110000
	ya.i = (ya.i&0x807fffff)|(b<<23);

//	méthode de Newton
	ya.f=(ya.f+x/ya.f)*0.5;
	ya.f=(ya.f+x/ya.f)*0.5;
	ya.f=(ya.f+x/ya.f)*0.5;
	return ya.f;
}

int main(int argc, char *argv[])
{
	float a = 54444.000549, b;
	clock_t t;
	int i, n=10000000;

	printf("%f\n", sqrtf(a));
	printf("%f\n", rac(a));

	t = clock();
	for(i=0;i<n;++i) sqrtf(a);
	printf("%i\n", clock()-t);

	t = clock();
	for(i=0;i<n;++i) rac(a);
	printf("%i\n", clock()-t);

	getchar();
}

A voir également

Ajouter un commentaire

Commentaires

ctx_man
Messages postés
285
Date d'inscription
mardi 28 décembre 2004
Statut
Membre
Dernière intervention
20 janvier 2013

En l'occurence sur ce genre d'opération, tous les compilos que j'ai testé détectent s'il va y avoir perte de donnée et te préviennent. Mais je suis d'accord sur le fait qu'il ne faut pas présumer de ce que va coder le compilo, je l'ai fais ainsi parce que j'ai vérifié ce que ca donnais avec le mien, mais comme je le disais également, pour en être véritablement certain, il faut caster, exactement comme tu le montres.
BruNews
Messages postés
21042
Date d'inscription
jeudi 23 janvier 2003
Statut
Modérateur
Dernière intervention
21 août 2019
18
ya.c[3] = (ya.c[3] + 127) / 2;
On ne doit surtout pas présupposer sur ce qui est codé par le compilo, quel qu'il soit.
Soit on vérifie le listing ASM et si calculs faits sur DWORD alors ok mais si faits sur AL que nenni.
Dans tous les cas je déconseille fort cette manière, aucune assurance que ira encore quand on change de compilo ou simplement de version de compilo.

DWORD v = (DWORD) ya.c[3];
ya.c[3] = (BYTE) ((v + 127) / 2);
C'est pas épuisant à taper et on est assuré du bon résultat.
ctx_man
Messages postés
285
Date d'inscription
mardi 28 décembre 2004
Statut
Membre
Dernière intervention
20 janvier 2013

Personnellement je ne pense pas que gcc soit moins bon en optimisation que vc++.
Il y a des tonnes de flags qui peuvent jouer et tellement de choses peuvent influer sur le résultat du bench présent dans ta fonction "main".
Concernant le dépassement de l'opération, normalement ca passe correctement. ya.[3] + 127 fait plus que 255, oui, mais il n'est pas stocké, donc il se retrouve soit en registre, soit en pile, une fois divisé par 2, le nombre ne peux plus faire plus de 8 bits, et donc tiens sur un octet, et là il est stocké dans ya.c[3]. Pour en être absolument certain, tu peux caster.
Après, l'utilisation de ce tableau est surtout une facilité d'écriture du programme (plus simple qu'un masque), question perf ca doit jouer a peine, pas facile de vérifier. J'ai repris le code tel qu'il est actuellement et j'ai ajouté une seconde fonction en utilisant le tableau de char, chez moi, 6 fois sur 10, la fonction utilisant le char est plus rapide de quelques millisecondes, mais c'est pas non plus loins devant. Faudrait refaire le teste encore et encore pour savoir si en moyenne une fonction est plus rapide que l'autre. Le mieux étant de décompiler pour voir l'assembleur résultant et compter le nombre de cycles consommés. Une comparaison par appel de l'horloge ne fait que donner un ordre d'idée.
Si tu veux continuer a optimiser encore plus, il faut passer en assembleur surtout que chaque compilateur va donner un assembleur différent et donc des performances différentes (avec un compilateur une fonction sera plus rapide que l'autre, et ce sera l'inverse sur un autre compilateur). Mais bon, de toute façon tu n'atteindra pas le score de la fonction standard qui se contente d'utiliser l'instruction assembleur du coprocesseur mathématique ou de SSE2 et qui donc te sort le résultat en une seule instruction.
cs_Jhep
Messages postés
68
Date d'inscription
dimanche 31 mars 2002
Statut
Membre
Dernière intervention
18 janvier 2010

bon je viens de tester la vitesse et contrairement à mon attente l'utilisation des char augmente le temps d'exécution. mais bon gcc n'est pas forcement aussi bien que vc++ en terme d'optimisation.
mais en y repensant ca ne devrait même pas marcher car l'opération ya.c[3] = (ya.c[3] + 127) / 2 dépasse le plus grand nombre sur 1 octet : 255...
ctx_man
Messages postés
285
Date d'inscription
mardi 28 décembre 2004
Statut
Membre
Dernière intervention
20 janvier 2013

Oui oui, on s'en tamponne du signe, je l'ai conservé uniquement pour la pédagogie. C'était un bit de perdu lors du déplacement binaire du coup je préférai montrer comment le conserver et le restaurer. Mais je suis d'accord, vu que c'est une racine carrée le signe est forcément positif, donc 0, or ce sont des 0 qui sont insérés dans les décalages. (Encore que, je crois me souvenir que quand tu décales, le dernier bit sorti est posé dans un flag de retenu et est ré-inséré si tu re-décale dans l'autre sens, mais je suis pas catégorique la dessus, ca demande teste)

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.