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

Commenter

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

S'inscrire maintenant

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.

Du même auteur (cs_Jhep)