Problème de précision ?

Bonjour,

Je dois chercher dans un triangles équilatéral de 273 de côté, les points à l'intérieur du triangle pour lesquels la distance aux trois sommets de celui-ci sont des valeurs entières.

Voici les résultats :

http://img225.imageshack.us/content_...ng&via=mupload

Le problème est le suivant. Une fois trouvé les valeurs des coordonnées donnant le résultat attendu, je reinjecte ces coordonnées dans les équations de cercles afin de vérifier que j'ai bien calculé ce qu'il fallait. Cependant, comme vous pourrez le remarquer sur l'image ci dessus. Pour un des résultats, je ne récupère pas une valeur entière.

Je voulais savoir si cela venait d'un problème de précision ou si cela pourrait venir de ma fonction qui détermine si une valeur est entière ( unsigned int isInteger( float x ) ) ?

Quelques explications supplémentaires pour s'imprégner du code :

Les résultats sont organisés par groupe de 3 lignes
ligne 1 : (x,y) les coordonnées du point qui donne des distances entières aux trois sommets.
lignes 2 : Les distances en question aux 3 sommets
lignes 3 : Vérification des valeurs trouvées. J...'ai réinjécté le couple (x,y) dans les equations de cercles afin de voir que mes calculs étaient justes.

En effet je retrouve bien les valeurs déterminées par l'algorithme mais pas complétement pour une des valeurs.

J'ai pris 3 points A(0,0) B(273/2, 273*V3/2) C(0,273). Et j'ai déduis les équations de cercles de centre A, B, C.

Concrètement ce que je fais :

Je fais varier les rayons des cercles de centre A et B par pas de 1. Donc double for imbriquées.
Ensuite je regarde s'il y a intersection entre les cercles par un tas de calculs ennuyeux.
Je renvois la solution qui appartient au triangle équilatéral.
Du point d'intersection trouvé, je calcul la distance au point C.
Si cette valeur est entière je retourne le résultat sinon je continue.
Pour obtenir les autre résultats, il suffit de faire des symétries... Mais cela donnerait de toutes façons des combinaisons différentes des valeurs déjà trouvées.

Voici le code source :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

#include<stdio.h> #include<stdlib.h> #include<math.h>   float checkRa( float x, float y ) { return sqrt( x*x + y*y ); }   float checkRb( float x, float y ) { return sqrt( (x-273/2)*(x-273/2) + (y-273*sqrt(3)/2)*(y-273*sqrt(3)/2) ); }   float checkRc ( float x, float y ) { return sqrt( (x-273)*(x-273) + y*y ); }   float DELTA( float A, float B, float C ) { return sqrt( B*B-4*A*C ); }   float x( float A, float B, float D ) { return (-B+D)/(2*A); }   float y( float N, float x, float xa, float ya, float xb, float yb ) { return N-x*(xb-xa)/(yb-ya); }   float N( float ra, float rb, float xa, float xb, float ya, float yb ) { return (ra*ra-rb*rb-xa*xa+xb*xb-ya*ya+yb*yb)/(2*(yb-ya)); }   float A( float xa, float ya, float xb, float yb ) { return (xb-xa)*(xb-xa)/((yb-ya)*(yb-ya)) + 1; }   float B( float N, float xa, float ya, float xb, float yb ) { return 2*yb*(xb-xa)/(yb-xa) - 2*N*(xb-xa)/(yb-xa) - 2*xb; }   float C( float xb, float yb, float N, float rb ) { return xb*xb+yb*yb+N*N-rb*rb-2*yb*N; }   float distanceAB( float xa, float ya, float xb, float yb ) { return sqrt( (xb-xa)*(xb-xa)+(yb-ya)*(yb-ya) ); }   float isInteger( float x ) { if ( x-(int)x == 0 ) { return 1; }   else return 0; }   int main (void) { float xa = 0; float ya = 0; float xb = 273/2; float yb = 273*sqrt(3)/2; float xc = 273; float yc = 0;   size_t i,j; for ( i=0 ; i<=273 ; i++ ) { for ( j=0 ; j<=273 ; j++ ) { float _N = N( i, j, xa, xb, ya, yb); float _A = A( xa, ya, xb, yb ); float _B = B( _N, xa, ya, xb, yb ); float _C = C( xb, yb, _N, j );   if ( DELTA( _A, _B, _C ) > 0 ) { float _x = x( _A, _B, DELTA( _A, _B, _C ) ); float _y = y( _N, _x, xa, ya, xb, yb );   if ( isInteger( distanceAB( _x, _y, xc, yc ) ) == 1 ) { printf("\n x=%.20f, y=%.20f\n %d - %d - %f\n %.20f - %.20f - %.20f\n ", _x, _y, i, j, distanceAB( _x, _y, xc, yc ), checkRa(_x, _y), checkRb(_x, _y), checkRc(_x, _y) ); } } } }   return 0; }

---

Merci d'avance pour vos lumières.

Qia,

En utilisant des réels (float), tous les calculs sont des approximations (entre 2 réels quelconques différents, il y a une infinité de réels ...).

Tu peux améliorer un peu en utilisant des double, mais ce problème est irréductible. :(

Ok, du coup je sors du cadre du topic mais sait-on jamais... Auriez vous une autre idée pour résoudre se problème puisque celui-ce n'est pas viable...

Bonsoir,

Ce que vous pouvez faire c'est faire une structure nombre avec un champ int mantisse et un autre champ int exposant, si vos nombres n'ont pas trop de décimales, ca pourrait peut-être convenir

Selon droggo il n'y a pas de moyen...
Je parlais de la méthode pour déterminer les valeurs.

Citation:

---

Envoyé par darkwall_37

Selon droggo il n'y a pas de moyen...
Je parlais de la méthode pour déterminer les valeurs.

---

Ce problème est récurrent et on le retrouve dans d'autres domaines. Par exemple, dans les transactions bancaires, il est inconcevable d'avoir une erreur due à l'arrondi.
Il existe une librairie qui calcule les nombres décimaux en précision absolue. Elle travaille en fait comme dans le type decimal(n, m). Chaque chiffre est écrit sur une position et les calculs se font comme au primaire, chiffre par chiffre, colonne par colonne.

Toi, ton problème est que tu as des fractions et des racines. Ben tu peux te lancer à créer une librairie pour gérer les fractions. Une structure "s_frac" contenant un numerateur, un denominateur. Puis tu codes les opérations "addition", "soustraction", "multiplication" et "division" etc.

Et tu fais de même avec les racines...

Si je comprend bien le problème :
Soit un triangle équilatéral ABC de coté L entier. Trouver tous les points M situés sur ou à l'intérieur du triangle pour lesquels la distance de M aux sommets A, B et C soit un entier.

Il s'agit d'un problème de valeurs entières et il faut impérativement ne pas faire d'approximations liés aux calcul sur flottants (et notamment utiliser la racine carrée) qui peuvent engendrer de fausses solutions et des solutions approchées.
On doit donc formuler le problème en terme d'entiers.

En conservant les coordonnées que tu utilises pour les sommets A, B, C et en appelant da, db, dc la distance d'un point M (x,y) à ces sommets, on a :

Citation:

---

da^2 = x^2 + y^2
db^2 = da^2 + L^2 - Lx -LyV3
dc^2 = da^2 + L^2 -2Lx

---

On cherche x et y pour avoir da, db, dc entiers.
Si da est entier, da^2 est entier , L^2 est entier et da^2 + L^2 est entier.
Pour avoir aussi db entier, il faut donc que Lx+LyV3 soit entier : Lx+LyV3 = p
Pour avoir aussi dc entier, il faut donc que 2Lx soit entier : 2Lx = m
x et y doivent donc vérifier

Citation:

---

x = m/(2L)
y =(p-m/2)/(LV3)= V3*(2p-m)/(6L)

---

Alors, on doit avoir :

Citation:

---

3(Lda)^2 = m^2 + p^2 -pm
db^2 = da^2 + L^2 - p
dc^2 = da^2 + L^2 - m

---

On a maintenant à résoudre un problème entre entiers et il faut déterminer p, m, da, db, dc qui satisfont ce système.

On peut se limiter, compte tenu des symétries du triangle à rechercher les points sur 1/6 du triangle , le triangle AOH où O est le centre du triangle et H le milieu de AC. Alors on a 0 <= m <= L^2 et m/2 <= p <= m.
Pour ce triangle AOH, on a 0 <= da <= L/V3 et da <= dc <= db
On peut envisager la démarche suivante :

Code:

---

1 2 3 4 5 6

Pour da = 0...~L/V3 pour da <= dc <=L soit m = da^2 + L^2 - dc^2 pour dc <= db <= L soit p = da^2 + L^2 - db^2 vérifiant m/2 <= p <= m Si on a un tel couple (p,m), vérifier que 3(Lda)^2 = m^2 + p^2 -pm Si c'est le cas, on a une solution exacte.

---

Ce qui correspond au code suivant :

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

#include <stdio.h> #define L 273 int main(void) { int da,db,dc ; long m,p; long S[L+1]; int damax; int i; double const V3 = 1.7320508076; double v1; for(da=0; da<=L ; da++) S[da] = da*da; damax = L*V3 + 0.001; for(da=0; da<=damax ; da++) { v1 = 3.0*S[L]*S[da]; for(dc=da; dc<L; dc++) { m = S[da]+S[L]-S[dc]; for(db=dc; db<L ; db++) { p = S[da]+S[L]-S[db]; if(p>=m/2 && p<=m && v1 == (double)m*m +(double)p*(p-m)) { printf("m = %d p = %d\nda = %3d db = %3d dc = %3d\n", m,p,da, db, dc); printf("x = %d/%d [=%f]\ny = (%d/%d)sqrt(3) [=%f]\n\n", m,2*L, 0.5*m/L, (2*p-m),6*L, V3*(2*p-m)/6/L); } } } } return 0; }

---

donnant ces résultats très rapidement:

Citation:

---

m = 35490 p = 17745
da = 65 db = 247 dc = 208
x = 35490/546 [=65.000000]
y = (0/1638)sqrt(3) [=0.000000]

m = 49713 p = 40674
da = 97 db = 208 dc = 185
x = 49713/546 [=91.049451]
y = (31635/1638)sqrt(3) [=33.451421]

m = 65520 p = 32760
da = 120 db = 237 dc = 153
x = 65520/546 [=120.000000]
y = (0/1638)sqrt(3) [=0.000000]

---

Note :
Les long (32 bits min) sont largement suffisants pour les données m et p en L^2 mais insuffisants pour les expressions de v1 et m*m +p*(p-m) de l'ordre de L^4. Si on en dispose, on peut utiliser des long long (64 bits min). Sinon, le code présenté utilise des doubles qui ont une mantisse couramment supérieure à 50 bits ce qui assure une représentation exacte de ces entiers.

Wahou, je pensais le topic abandonné ! Merci notament à Diogene pour ton implication dans ce problème. Je n'ai pas pensé à cette approche, jolie et vraiment merci beaucoup !