calcul d'une norme en entier

Version imprimable

Je dispose de
1- X entier codé sur n bits
2 - Y entier codé sur n bits

je sais de plus que la norme du vecteur (X,Y) peut être codée sur n bits
(par contre X^2 tout comme Y^2 peut dépasser n bits!)

quelqu'un aurrait-il/elle un algorithme éfficace pour calculer la valeur entière de
sqrt(X^2 + Y^2) Sans passer par des floatting point

pour info n=40 dans mon cas, le nombre de bits effectivement utiliés va de 20 a 38 ou 39.

Actuellement j'ai créé des buffers de 80 bits pour stocker X^2 et Y^2
un buffer de 81 bits pour la somme et je cacule le sqrt avec la méthode
"de la potence " adaptée à la base 2.

Cela me semble un peu lourd!!

Merci de vos réponses

Voici mon code qui fonctionne mais me semble non optimal

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
 
#define Nbits   40
unsigned long long L1, L2;
unsigned short  W1[Nbits], W2[Nbits],
                W1_sqr[2*Nbits], W2_sqr[2*Nbits],
                W_sm[2*(Nbits+1)]; // en principe 2*Nbits+1 suffisent mais un nombre paire est + facile a programmer dans 'pendule'
void Bit40( unsigned long long L, unsigned short * W)
   {
   short i;
   for (i=0; i < Nbits; i++)
      {
      W[i] = L & 1;
      L = L>> 1;
      }
   }
void LSQR( unsigned short * W,unsigned short * W2)
   {
   short i;
   memset(W2,0,2*Nbits*sizeof(short));
   for (i=0; i< Nbits; i++) if ( W[i] )
      {
      short retenue =0;
      short k;
      for (k=0; k < Nbits; k++)
         {
         short u= W2[i+k] + W[k] + retenue;
         if ( u <= 1)
            {
            W2[i+k] = u;
            retenue = 0;
            }
         else
            {
            W2[i+k] = u - 2;
            retenue = 1;
            }
         }
      }
   }
void Add2( unsigned short * W1_sqr, unsigned short * W2_sqr,unsigned short * W_sm)
   {
   short k,   retenue = 0;
   memset(W_sm,0,2*(Nbits+1)*sizeof(short));
   for (k=0; k < 2*Nbits; k++)
      {
      short m = retenue+W1_sqr[k] + W2_sqr[k] ;
      if ( m <=1)
         {
         W_sm[k] = m;
         retenue = 0;
         }
      else
         {
         W_sm[k] = m-2;
         retenue = 1;
         }
      }
   W_sm[2*Nbits]=retenue; //W_sm[2*Nbits+1] est toujours nul
   }
void potence (unsigned short *W_sm, unsigned long long * square, unsigned long long * reste)
   {
   *square=0;
   * reste =0;
   short i = Nbits;
   while (1)
      {
      long long L;
      short j=2*i+1;
      short k=2*i;
      * reste = (* reste) * 4  + (W_sm[j] << 1) +  W_sm[k];
      L = 4* (*square) + 1;
      if ( L <= (* reste) )
         {
         (* reste) -= L;
         (*square) = 2*(*square)  +1;
         }
      else
         (* square) =  2 * ( * square);
      i--;
      if ( i < 0) break;
      }
   }
void  Norme(unsigned long long L1, unsigned long long L2, unsigned long long * L, unsigned long long * R)
   {
   Bit40(L1,W1);
   Bit40(L2,W2);
   LSQR(W1,W1_sqr);
   LSQR(W2,W2_sqr);
   Add2(W1_sqr,W2_sqr,W_sm);
   potence (W_sm,  L,R);
   }

Note: J'utilise les unsigned short car au final le SW doit aller sur un DSP de type C5510 qui ne connait pas d'entité de 8 bits.

27/05/2008, 22h56
Zavonen

racine(u) = u ^1/2
racine (u) = racine(1 + (u-1))^1/2
Je suggère donc d'utiliser le développement limite de (1+v)^s où v=(u-1)
Malheureusement ce développement limité n'est valable que pour v <1 (me semble-t-il de mémoire).
Alors si par exemple x >y
Commencer par mettre x en facteur et écrire:
(x^2+y^2) ^1/2 = x(1+(y/x)^2)^1/2
et utiliser le truc précédent avec v=y/x
28/05/2008, 13h35
Jean-Marc.Bourguet

Je me batirais des fonctions de calcul en 96 ou 128 bits (a partir de ce qui est disponible, 32 ou 64 bits -- as-tu verifier que tu disposais bien de long long sur ton DSP?) et pour le calcul de la racine carre, j'utiliserais une approximation (si tu utilises 2 long long pour faire des entiers de 128 bits, prendre la racine carree de la partie haute et la decaler de 32 bits -- si c'est nul, prendre la racine carree de la partie basse) puis quelques iterations de N-R:
x_n+1 = (x_n + a/x_n)/2 devrait te donner toute la precision dont tu as besoin.

Le plus dur, c'est d'ecrire et de tester la division faite a partir de mots plus petit. J'ai donne l'algo de base ici par le passe (http://www.developpez.net/forums/sho...40&postcount=7), il y a des algos plus complexes mais je doute de leur utilite dans le cas 2 ou 4 chiffres (de 64 ou 32 bits).
28/05/2008, 18h54
j.p.mignot

Citation:

as-tu verifier que tu disposais bien de long long sur ton DSP?

Le C5510 supporte un type 'long long' qui manipule des entitées entières signées ou non de 40 bits.

Citation:

x_n+1 = (x_n + a/x_n)/2

Je voulais éviter la méthode de newton qui implique implicitement le float.
Par ailleurs même si avec des arrondis liés à une division entière la série converge vers la bonne valeur, compte tenu que a peut occuper jusqu'à 81 bits,
cette méthode continue à maintenir les problèmes de manipulation des nombres assez grands ( hors du codage max normalement soutenu par mon DSP)
.
L'extraction de sqrt(n) {n sur 82 bits} par la méthode de la potence en base 2 comme proposées implique
41 comparaisons et éventuellement soustractions d'entiers
41 shift left et incrément de 1
pour garentir la partie entière de sqrt(n)

Je ne suis pas sur que cela soit nécessairement plus lourd que 5 ou 6 ittétations de newtown.
28/05/2008, 22h23
Jean-Marc.Bourguet

Citation:

Envoyé par j.p.mignot

Le C5510 supporte un type 'long long' qui manipule des entitées entières signées ou non de 40 bits.

J'avais pensé à l'absence, pas à la présence d'un type non standard (la norme C demande que long long fasse au moins 64 bits) :-(

Citation:

Je voulais éviter la méthode de newton qui implique implicitement le float.
Par ailleurs même si avec des arrondis liés à une division entière la série converge vers la bonne valeur, compte tenu que a peut occuper jusqu'à 81 bits,

En entier (avec troncature, pas arrondi), si j'ai bonne mémoire NR converge vers l'entier inferieur ou égal à la solution.

Citation:

cette méthode continue à maintenir les problèmes de manipulation des nombres assez grands ( hors du codage max normalement soutenu par mon DSP)

Ton problème exige que tu manipules de tels nombres. La question est donc comment tu les représentes. Choisir la base 2 plutôt qu'une autre (2^16 et 2^20 sont ceux qui viennent à l'esprit) a l'avantage de la simplicité, mais certainement pas celui des perf.

Citation:

L'extraction de sqrt(n) {n sur 82 bits} par la méthode de la potence en base 2 comme proposées implique
41 comparaisons et éventuellement soustractions d'entiers
41 shift left et incrément de 1
pour garentir la partie entière de sqrt(n)

Je ne suis pas sur que cela soit nécessairement plus lourd que 5 ou 6 ittétations de newtown.

Tu oublies des choses dans tes comptes: j'en suis déjà à 5 shifts par itération rien que pour potence, il faut compter aussi les opérations faites pour la transformation (déjà 80 shift) et celles pour les calculs avant la racine carrée (il y une boucle où tu passes déjà potentiellement 3200 fois -- 2*40*40). J'ai comme l'impression qu'utiliser des chiffres de 16 ou 20 bits et NR serait quand même moins coûteux au total. Naturellement, une évaluation précises en cycles ou des mesures -- sur le DSP, pas sur un PC -- sont nécessaires pour confirmer ou infirmer, mais la situation me semble moins claire que tu ne le penses. Si on démarre NR par une estimation dans une table, il ne faut pas beaucoup (on double le nombre de chiffres par itération, si tu parts avec 5 bits de précision, il en faut donc 3) -- mais la méthode la plus rapide est peut-être un mélange -- travailler par chiffres de 16 bits mais un calcul de racine carrée par potence par exemple.

Il n'y a pas de bibliothèques donnant une racine carrée entière pour ton DSP?
29/05/2008, 19h49
j.p.mignot

Citation:

Il n'y a pas de bibliothèques donnant une racine carrée entière pour ton DSP?

of course! mais ssi n est un entier reconu codable sur 40 bit max.
Ici je suis obligé d'envisager des entiers pouvant aller jusqu'à 81 bits.
( somme de 2 carrés d'entiers 40 bits max donc somme de 2 entiers 80 bits d'où 81)

Là je dois moi même définir toutes les opérations (+, *, /, sqr, sqrt, ...)

a> * de 2 <40 bits> donne au max de 40*40 sommes de bits
( à faire sur X et sur Y ) et 2 résultats codés 80 bits

b> + de 2 <80 bits> donne 80 sommes de bits et un résultat sur 81 bits

c> Sqrt(somme)

'à prioris' l'opérarion (a) à faire sur X et sur Y me semble plus lourde que sqrt

Citation:

Envoyé par j.p.mignot

Ici je suis obligé d'envisager des entiers pouvant aller jusqu'à 81 bits.
( somme de 2 carrés d'entiers 40 bits max donc somme de 2 entiers 80 bits d'où 81)

Là je dois moi même définir toutes les opérations (+, *, /, sqr, sqrt, ...)

a> * de 2 <40 bits> donne au max de 40*40 sommes de bits
( à faire sur X et sur Y ) et 2 résultats codés 80 bits

b> + de 2 <80 bits> donne 80 sommes de bits et un résultat sur 81 bits

c> Sqrt(somme)

'à prioris' l'opérarion (a) à faire sur X et sur Y me semble plus lourde que sqrt

Si tu es sûr que la racine carrée tient sur 40 bits, tu sais que la somme tient sur 80 bits.

Pour le carré (codé dans le browser):
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 typedef unsigned long long u40; #define MASK20 (1ULL<<20 - 1) struct u80 { u40 l; u40 h; }; void square(u40 i, u80* r) { u40 il = u40 & MASK20; u40 ih = u40 >> 20; u40 acc = il*il; u40 tmp = acc & MASK20; acc >>= 20; acc += 2* ((il*ih) & MASK20); r->l = tmp + (acc & MASK20) << 20; acc >>= 20; r->h = acc + ih*ih; }
L'addition est du même genre. La division, c'est le plus compliqué et j'ai donné l'algo ci-dessus. Tu as peut-être intérêt à travailler avec 6 chiffres sur 16 bits et profiter de l'accumulateur sur 40 bits pour simplifier le traitement des reports plutôt que de gérer les reports en utilisant quatre chiffres sur 20 bits ou deux chiffres sur 40 que tu redécouperais tout le temps en deux.

Version 6 chiffres sur 16 bits.
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 typedef unsigned short u16; typedef unsigned long long u40; #define NDIGITS 6 #define MASK16 ((1 << 16) - 1) typedef struct u96 { u16 d[NDIGITS]; } u96[1]; void square(u96 r, u40 l) { u40 l = a & MASK16; u40 m = (a >> 16) & MASK16; u40 h = a >> 32; u40 acc = l*l; r->d[0] = acc & MASK16; acc >>= 16; acc += 2 * l * m; r->d[1] = acc & MASK16; acc >>= 16; acc += 2 * h * l + 2 * m * m; r->d[2] = acc & MASK16; acc >>= 16; acc += 2 * h * m; r->d[3] = acc & MASK16; acc >>= 16; acc += m * m; r->d[4] = acc & MASK16; acc >>= 16; r->d[5] = acc; }