Optimiser moi! (C, les performances) -- suivi de peu-tourner à la question

Question

Grâce à quelques très utile stackOverflow utilisateurs à Peu se tourner: qui bit est défini?, J'ai construit ma fonction (posté à la fin de la question).

Toutes les suggestions -- même les petites suggestions seraient appréciées. J'espère qu'il va faire de mon mieux, mais au moins il devrait m'apprendre quelque chose. :)

Vue d'ensemble

Cette fonction sera appelée à au moins 10^{à 13} fois, et éventuellement aussi souvent que 10^{à 15}. C'est, ce code sera exécuté par mois en toute probabilité, de sorte que toute la performance conseils seraient utiles.

Cette fonction comptes pour 72-77% du programme du temps, basé sur le profilage et une douzaine de courses dans différentes configurations (optimisation de certains paramètres qui ne sont pas pertinentes ici).

Au moment où la fonction s'exécute dans une moyenne de 50 horloges. Je ne sais pas combien cela peut être amélioré, mais je serais ravis de les voir courir dans 30.

La Clé De L'Observation

Si, à un certain point dans le calcul, vous pouvez dire que la valeur qui sera retournée sera petite (valeur exacte négociable -- dire, en dessous d'un million de) , vous pouvez abandonner tôt. Je ne suis intéressé dans de grandes valeurs.

C'est de cette façon j'ai l'espoir de sauver le plus de temps, plutôt que par d'autres micro-optimisations (si ce sont bien sûr les bienvenus aussi!).

Information Sur Le Rendement

• smallprimes est un tableau de bits (64 bits); la moyenne est d'environ 8 bits sera réglé, mais il pourrait être aussi peu que 0 ou 12.
• q sera généralement différente de zéro. (Notez que la fonction se termine tôt si q et smallprimes sont à zéro.)
• r et s sera souvent 0. Si q est égal à zéro, r et s sera trop; si r est égal à zéro, s sera trop.
• Comme le commentaire dit à la fin, nu est généralement de 1 à la fin, j'ai donc efficace d'un cas particulier pour elle.
• Les calculs ci-dessous le cas particulier peut apparaître pour des risques de débordement, mais par le biais de modélisation appropriées j'ai prouvé que, pour mon entrée, cela ne se fera pas -- ne vous inquiétez pas au sujet de cette affaire.
• Fonctions qui ne sont pas définis ici (ugcd, minuu, étoiles, etc.) ont déjà été optimisés; aucun temps pour s'exécuter. pr est un petit tableau (tous en L1). Aussi, toutes les fonctions appelées ici sont des fonctions pures.
• Mais si vous vous souciez vraiment... ugcd est le pgcd, minuu est le minimum, vals-les-bains est le nombre de fuite binaire 0, __builtin_ffs est l'emplacement le plus à gauche binaire 1, la star est de (n-1) >> vals(n-1), pr est un tableau de nombres premiers de 2 à 313.
• Les calculs sont en cours sur un Phenom II x4 920, même si des optimisations pour le i7 ou Woodcrest sont encore de l'intérêt (si je reçois de temps de calcul sur les autres nœuds).
• Je serais heureux de répondre à toutes les questions que vous avez au sujet de la fonction ou de ses mandants.

Ce qu'il fait réellement

Ajouté en réponse à une demande. Vous n'avez pas besoin de lire cette partie.

L'entrée est un nombre impair n, avec 1 < n < 4282250400097. Les autres entrées de fournir la factorisation du nombre dans ce sens:

smallprimes&1 est défini si le nombre est divisible par 3, smallprimes&2 est défini si le nombre est divisible par 5, smallprimes&4 est défini si le nombre est divisible par 7, smallprimes&8 est défini si le nombre est divisible par 11, etc. jusqu'à le bit le plus significatif qui représente 313. Un nombre divisible par le carré d'un nombre premier n'est pas représenté de manière différente d'un nombre divisible par juste ce nombre. (En fait, les multiples de places peut être mis au rebut; dans l'étape de prétraitement dans une autre fonction des multiples de carrés de nombres premiers <= lim ont smallprimes et q définie sur 0, de sorte qu'ils seront supprimées, où la valeur optimale de la mfr est déterminé par l'expérience.)

q, r et s représentent les plus grands facteurs de la le nombre. Tout en restant le facteur (qui peut être supérieure à la racine carrée du nombre, ou si s est non nulle peut-être même moins) peuvent être obtenus en divisant les facteurs de n.

Une fois que tous les facteurs sont récupérés de cette façon, le nombre de bases, 1 <= b < n, dans laquelle n est un fort pseudoprime sont comptées à l'aide d'une formule mathématique mieux expliqué par le code.

Des améliorations jusqu'à présent

• Poussé le début de la sortie de test. De toute évidence, cela permet d'économiser du travail, donc j'ai fait le changement.
• Les fonctions appropriées sont déjà en ligne, alors __attribute__ ((inline)) ne fait rien. Curieusement, le marquage de la fonction principale bases et de certaines aides avec __attribute ((hot)) nuire à la performance de près de 2% et je ne peux pas comprendre pourquoi (mais c'est reproductible avec plus de 20 tests). Donc je n'ai pas d'effectuer le changement. De même, __attribute__ ((const)), au mieux, n'a pas aidé. J'ai été plus que légèrement surpris par cette.

Code

ulong bases(ulong smallprimes, ulong n, ulong q, ulong r, ulong s)
{
    if (!smallprimes & !q)
        return 0;

    ulong f = __builtin_popcountll(smallprimes) + (q > 1) + (r > 1) + (s > 1);
    ulong nu = 0xFFFF;  // "Infinity" for the purpose of minimum
    ulong nn = star(n);
    ulong prod = 1;

    while (smallprimes) {
        ulong bit = smallprimes & (-smallprimes);
        ulong p = pr[__builtin_ffsll(bit)];
        nu = minuu(nu, vals(p - 1));
        prod *= ugcd(nn, star(p));
        n /= p;
        while (n % p == 0)
            n /= p;
        smallprimes ^= bit;
    }
    if (q) {
        nu = minuu(nu, vals(q - 1));
        prod *= ugcd(nn, star(q));
        n /= q;
        while (n % q == 0)
            n /= q;
    } else {
        goto BASES_END;
    }
    if (r) {
        nu = minuu(nu, vals(r - 1));
        prod *= ugcd(nn, star(r));
        n /= r;
        while (n % r == 0)
            n /= r;
    } else {
        goto BASES_END;
    }
    if (s) {
        nu = minuu(nu, vals(s - 1));
        prod *= ugcd(nn, star(s));
        n /= s;
        while (n % s == 0)
            n /= s;
    }

    BASES_END:
    if (n > 1) {
        nu = minuu(nu, vals(n - 1));
        prod *= ugcd(nn, star(n));
        f++;
    }

    // This happens ~88% of the time in my tests, so special-case it.
    if (nu == 1)
        return prod << 1;

    ulong tmp = f * nu;
    long fac = 1 << tmp;
    fac = (fac - 1) / ((1 << f) - 1) + 1;
    return fac * prod;
}

Answer 1

Vous semblez perdre beaucoup de temps à faire des divisions par les facteurs. Il est beaucoup plus rapide pour remplacer une division avec une multiplication par l'inverse du diviseur (division: ~15-80(!) cycles, selon le diviseur, la multiplication: ~4 cycles), SI bien sûr, vous pouvez précalculer les inverses.

Bien que cela semble peu probable, possible, avec q, r, s - en raison de la répartition de ces revendeurs à valeur ajoutée, il est très facile à faire avec p, qui vient toujours de la petite, statique pr[] array. Précalculer les inverses de ces nombres premiers et de les stocker dans un autre tableau. Alors, au lieu de diviser par p, multiplier par l'inverse prises à partir de la deuxième matrice. (Ou faire un seul tableau de structures.)

Maintenant, l'obtention de division exacte résultat obtenu par cette méthode nécessite une ruse pour compenser les erreurs d'arrondi. Vous trouverez les détails croustillants de cette technique dans le présent document, à la page 138.

EDIT:

Après consultation du Hacker Plaisir (un excellent livre, d'ailleurs) sur le sujet, il semble que vous pouvez le rendre encore plus rapide en exploitant le fait que toutes les divisions de votre code sont exacts (c'est à dire le reste est égal à zéro).

Il semble que, pour chaque diviseur d , ce qui est surprenant et la base B = 2^word_size, il existe un unique inverse multiplicatif d⃰ qui satisfait les conditions: d⃰ < B et d·d⃰ ≡ 1 (mod B). Pour tous les x qui est un multiple exact de d, ce qui implique x/d ≡ x·d⃰ (mod B). Ce qui signifie que vous pouvez simplement remplacer une division avec une multiplication, pas d'ajout de corrections, les chèques, les problèmes d'arrondi, peu importe. (Les preuves de ces théorèmes peuvent être trouvés dans le livre.) Notez que cette multiplicatif inverse n'a pas besoin d' être égale à la réciproque tel que défini par la méthode précédente!

Comment vérifier si une donnée x est un multiple exact de d - d. x mod d = 0 ? Facile! x mod d = 0 mfi x·d⃰ mod B ≤ ⌊(B-1)/d⌋. Notez que cette limite supérieure peut être précalculées.

Ainsi, dans le code:

unsigned x, d;
unsigned inv_d = mulinv(d);          //precompute this!
unsigned limit = (unsigned)-1 / d;   //precompute this!

unsigned q = x*inv_d;
if(q <= limit)
{
   //x % d == 0
   //q == x/d
} else {
   //x % d != 0
   //q is garbage
}

En supposant que l' pr[] tableau devient un tableau d' struct prime:

struct prime {
   ulong p;
   ulong inv_p;  //equal to mulinv(p)
   ulong limit;  //equal to (ulong)-1 / p
}

l' while(smallprimes) boucle dans votre code devient:

while (smallprimes) {
    ulong bit = smallprimes & (-smallprimes);
    int bit_ix = __builtin_ffsll(bit);
    ulong p = pr[bit_ix].p;
    ulong inv_p = pr[bit_ix].inv_p;
    ulong limit = pr[bit_ix].limit;
    nu = minuu(nu, vals(p - 1));
    prod *= ugcd(nn, star(p));
    n *= inv_p;
    for(;;) {
        ulong q = n * inv_p;
        if (q > limit)
            break;
        n = q;
    }
    smallprimes ^= bit;
}

Et pour l' mulinv() fonction de:

ulong mulinv(ulong d) //d needs to be odd
{
   ulong x = d;
   for(;;)
   {
      ulong tmp = d * x;
      if(tmp == 1)
         return x;
      x *= 2 - tmp;
   }
}

Remarque: vous pouvez remplacer ulong avec n'importe quel autre type non signé - il suffit d'utiliser le même type de manière cohérente.

Les preuves, pourquois et comments sont tous disponibles dans le livre. Un vivement recommandé de lire :-).

Answer 2

Si votre compilateur GCC supporte les attributs des fonctions, vous pouvez marquer vos fonctions pures avec cet attribut:

ulong star(ulong n) __attribute__ ((const));

Cet attribut indique au compilateur que le résultat de la fonction ne dépend que de son argument(s). Ces informations peuvent être utilisées par l'optimiseur.

Est-il une raison pour laquelle vous avez opencoded vals() au lieu d'utiliser __builtin_ctz() ?

Answer 3

Il est encore un peu flou, ce que vous êtes recherche pour. Très souvent, des problèmes théoriques liés aux nombres permettent de grandes accélérations en dérivant des propriétés mathématiques que les solutions doivent satisfiy.

Si vous êtes à la recherche pour les entiers que de maximiser le nombre de non-témoins pour le M. de test (c'est à dire oeis.org/classic/A141768 que vous mentionnez), alors qu'il pourrait être possible d'utiliser que le nombre de non-témoins ne peuvent pas être plus grand que phi(n)/4 et que les entiers pour lesquels un grand nombre de non-témoins sont sont le produit de deux nombres premiers de la forme

(k+1)*(2k+1)

ou ils sont les nombres de Carmichael avec 3 facteurs premiers. Je pensais ci-dessus, la limite de tous les nombres entiers dans l'ordre ce formulaire, et qu'il est possible de vérifier cela en faisant preuve d'une limite supérieure pour les témoins de tous les autres entiers. E. g. entiers avec 4 ou plus des facteurs de toujours avoir dans la plupart des phi(n)/8 non-témoins. Des résultats similaires peuvent être dérivées à partir de, vous de la formule pour le nombre de bases pour les autres entiers.

Comme pour les micro-optimisations: Chaque fois que vous savez qu'un entier est divisible par certains quotient, alors il est possible de remplacer la division par une multiplication par l'inverse du quotient modulo 2^64. Et les tests n % q == 0 peut être remplacé par un essai

n * inverse_q < max_q,

où inverse_q = q^(-1) mod 2^64 et max_q = 2^64 / q. Évidemment inverse_q et max_q doivent être précalculées, pour être efficace, mais depuis que vous êtes à l'aide d'un tamis, je suppose que ce ne doit pas être un obstacle.

Answer 4

Petite optimisation, mais:

ulong f;
ulong nn;
ulong nu = 0xFFFF;  // "Infinity" for the purpose of minimum
ulong prod = 1;

if (!smallprimes & !q)
    return 0;

// no need to do this operations before because of the previous return
f = __builtin_popcountll(smallprimes) + (q > 1) + (r > 1) + (s > 1);
nn = star(n);

BTW: vous devez éditer votre post pour ajouter star() et d'autres fonctions que vous utilisez définition

Answer 5

Essayez de remplacer ce modèle (r et q trop):

n /= p; 
while (n % p == 0) 
  n /= p; 

Avec ceci:

ulong m;
  ...
m = n / p; 
do { 
  n = m; 
  m = n / p; 
} while ( m * p == n); 

Dans mes tests limités, j'ai eu une petite accélération (10%) de l'élimination du modulo.

Aussi, si p, q ou r sont constants, le compilateur remplace les divisions par des multiplications. Si il y a peu de choix pour p, q ou r, ou si certains d'entre eux sont de plus en plus fréquentes, vous pourriez obtenir quelque chose par la spécialisation de la fonction de ces valeurs.