IdentifiantMot de passe
Loading...
Mot de passe oublié ?Je m'inscris ! (gratuit)
Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Algorithmes et structures de données Discussion :

calcul d'une norme en entier


Sujet :

Algorithmes et structures de données

Vue hybride

Message précédent Message précédent   Message suivant Message suivant
  1. #1
    Membre émérite
    Inscrit en
    Juin 2005
    Messages
    644
    Détails du profil
    Informations professionnelles :
    Secteur : Industrie

    Informations forums :
    Inscription : Juin 2005
    Messages : 644
    Par défaut calcul d'une norme en entier
    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 C : Sélectionner tout - Visualiser dans une fenêtre à part
    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.

  2. #2
    Rédacteur
    Avatar de Zavonen
    Profil pro
    Inscrit en
    Novembre 2006
    Messages
    1 772
    Détails du profil
    Informations personnelles :
    Âge : 77
    Localisation : France

    Informations forums :
    Inscription : Novembre 2006
    Messages : 1 772
    Par défaut
    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
    Ce qu'on trouve est plus important que ce qu'on cherche.
    Maths de base pour les nuls (et les autres...)

  3. #3
    Expert confirmé

    Inscrit en
    Novembre 2005
    Messages
    5 145
    Détails du profil
    Informations forums :
    Inscription : Novembre 2005
    Messages : 5 145
    Par défaut
    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).

  4. #4
    Membre émérite
    Inscrit en
    Juin 2005
    Messages
    644
    Détails du profil
    Informations professionnelles :
    Secteur : Industrie

    Informations forums :
    Inscription : Juin 2005
    Messages : 644
    Par défaut
    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.


    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.

  5. #5
    Expert confirmé

    Inscrit en
    Novembre 2005
    Messages
    5 145
    Détails du profil
    Informations forums :
    Inscription : Novembre 2005
    Messages : 5 145
    Par défaut
    Citation Envoyé par j.p.mignot Voir le message
    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) :-(

    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.

    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.

    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?

  6. #6
    Membre émérite
    Inscrit en
    Juin 2005
    Messages
    644
    Détails du profil
    Informations professionnelles :
    Secteur : Industrie

    Informations forums :
    Inscription : Juin 2005
    Messages : 644
    Par défaut
    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

Discussions similaires

  1. [Débutant] Calcul d'un vecteur à partir d'une norme
    Par Globoxx dans le forum MATLAB
    Réponses: 6
    Dernier message: 31/10/2013, 17h43
  2. recuperer les minimum d'une séquence d'entiers?
    Par novice12 dans le forum Algorithmes et structures de données
    Réponses: 8
    Dernier message: 25/01/2005, 03h44
  3. Recuperer un champ calculé dans une variable....
    Par vijeo dans le forum MS SQL Server
    Réponses: 4
    Dernier message: 21/12/2004, 14h57
  4. calcul dans une requête
    Par blaz dans le forum Langage SQL
    Réponses: 8
    Dernier message: 22/12/2003, 10h31
  5. Réponses: 3
    Dernier message: 28/09/2003, 10h46

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo