Optimize me! (C, performance) -- followup to bit-twiddling question

Question

Thanks to some very helpful stackOverflow users at Bit twiddling: which bit is set?, I have constructed my function (posted at the end of the question).

Any suggestions -- even small suggestions -- would be appreciated. Hopefully it will make my code better, but at the least it should teach me something. :)

Overview

This function will be called at least 10¹³ times, and possibly as often as 10¹⁵. That is, this code will run for months in all likelihood, so any performance tips would be helpful.

This function accounts for 72-77% of the program's time, based on profiling and about a dozen runs in different configurations (optimizing certain parameters not relevant here).

At the moment the function runs in an average of 50 clocks. I'm not sure how much this can be improved, but I'd be thrilled to see it run in 30.

Key Observation

If at some point in the calculation you can tell that the value that will be returned will be small (exact value negotiable -- say, below a million) you can abort early. I'm only interested in large values.

This is how I hope to save the most time, rather than by further micro-optimizations (though these are of course welcome as well!).

Performance Information

• smallprimes is a bit array (64 bits); on average about 8 bits will be set, but it could be as few as 0 or as many as 12.
• q will usually be nonzero. (Notice that the function exits early if q and smallprimes are zero.)
• r and s will often be 0. If q is zero, r and s will be too; if r is zero, s will be too.
• As the comment at the end says, nu is usually 1 by the end, so I have an efficient special case for it.
• The calculations below the special case may appear to risk overflow, but through appropriate modeling I have proved that, for my input, this will not occur -- so don't worry about that case.
• Functions not defined here (ugcd, minuu, star, etc.) have already been optimized; none take long to run. pr is a small array (all in L1). Also, all functions called here are pure functions.
• But if you really care... ugcd is the gcd, minuu is the minimum, vals is the number of trailing binary 0s, __builtin_ffs is the location of the leftmost binary 1, star is (n-1) >> vals(n-1), pr is an array of the primes from 2 to 313.
• The calculations are currently being done on a Phenom II 920 x4, though optimizations for i7 or Woodcrest are still of interest (if I get compute time on other nodes).
• I would be happy to answer any questions you have about the function or its constituents.

What it actually does

Added in response to a request. You don't need to read this part.

The input is an odd number n with 1 < n < 4282250400097. The other inputs provide the factorization of the number in this particular sense:

smallprimes&1 is set if the number is divisible by 3, smallprimes&2 is set if the number is divisible by 5, smallprimes&4 is set if the number is divisible by 7, smallprimes&8 is set if the number is divisible by 11, etc. up to the most significant bit which represents 313. A number divisible by the square of a prime is not represented differently from a number divisible by just that number. (In fact, multiples of squares can be discarded; in the preprocessing stage in another function multiples of squares of primes <= lim have smallprimes and q set to 0 so they will be dropped, where the optimal value of lim is determined by experimentation.)

q, r, and s represent larger factors of the number. Any remaining factor (which may be greater than the square root of the number, or if s is nonzero may even be less) can be found by dividing factors out from n.

Once all the factors are recovered in this way, the number of bases, 1 <= b < n, to which n is a strong pseudoprime are counted using a mathematical formula best explained by the code.

Improvements so far

• Pushed the early exit test up. This clearly saves work so I made the change.
• The appropriate functions are already inline, so __attribute__ ((inline)) does nothing. Oddly, marking the main function bases and some of the helpers with __attribute ((hot)) hurt performance by almost 2% and I can't figure out why (but it's reproducible with over 20 tests). So I didn't make that change. Likewise, __attribute__ ((const)), at best, did not help. I was more than slightly surprised by this.

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

You seem to be wasting much time doing divisions by the factors. It is much faster to replace a division with a multiplication by the reciprocal of divisor (division: ~15-80(!) cycles, depending on the divisor, multiplication: ~4 cycles), IF of course you can precompute the reciprocals.

While this seems unlikely to be possible with q, r, s - due to the range of those vars, it is very easy to do with p, which always comes from the small, static pr[] array. Precompute the reciprocals of those primes and store them in another array. Then, instead of dividing by p, multiply by the reciprocal taken from the second array. (Or make a single array of structs.)

Now, obtaining exact division result by this method requires some trickery to compensate for rounding errors. You will find the gory details of this technique in this document, on page 138.

EDIT:

After consulting Hacker's Delight (an excellent book, BTW) on the subject, it seems that you can make it even faster by exploiting the fact that all divisions in your code are exact (i.e. remainder is zero).

It seems that for every divisor d which is odd and base B = 2^word_size, there exists a unique multiplicative inverse d⃰ which satisfies the conditions: d⃰ < B and d·d⃰ ≡ 1 (mod B). For every x which is an exact multiple of d, this implies x/d ≡ x·d⃰ (mod B). Which means you can simply replace a division with a multiplication, no added corrections, checks, rounding problems, whatever. (The proofs of these theorems can be found in the book.) Note that this multiplicative inverse need not be equal to the reciprocal as defined by the previous method!

How to check whether a given x is an exact multiple of d - i.e. x mod d = 0 ? Easy! x mod d = 0 iff x·d⃰ mod B ≤ ⌊(B-1)/d⌋. Note that this upper limit can be precomputed.

So, in 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
}

Assuming the pr[] array becomes an array of struct prime:

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

the while(smallprimes) loop in your code becomes:

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;
}

And for the mulinv() function:

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;
   }
}

Note you can replace ulong with any other unsigned type - just use the same type consistently.

The proofs, whys and hows are all available in the book. A heartily recommended read :-).

Answer 2

If your compiler supports GCC function attributes, you can mark your pure functions with this attribute:

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

This attribute indicates to the compiler that the result of the function depends only on its argument(s). This information can be used by the optimiser.

Is there a reason why you've opencoded vals() instead of using __builtin_ctz() ?

Answer 3

It is still somewhat unclear, what you are searching for. Quite frequently number theoretic problems allow huge speedups by deriving mathematical properties that the solutions must satisfiy.

If you are indeed searching for the integers that maximize the number of non-witnesses for the MR test (i.e. oeis.org/classic/A141768 that you mention) then it might be possible to use that the number of non-witnesses cannot be larger than phi(n)/4 and that the integers for which have this many non-witnesses are either are the product of two primes of the form

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

or they are Carmichael numbers with 3 prime factors. I'd think above some limit all integers in the sequence have this form and that it is possible to verify this by proving an upper bound for the witnesses of all other integers. E.g. integers with 4 or more factors always have at most phi(n)/8 non-witnesses. Similar results can be derived from you formula for the number of bases for other integers.

As for micro-optimizations: Whenever you know that an integer is divisible by some quotient, then it is possible to replace the division by a multiplication with the inverse of the quotient modulo 2^64. And the tests n % q == 0 can be replaced by a test

n * inverse_q < max_q,

where inverse_q = q^(-1) mod 2^64 and max_q = 2^64 / q. Obviously inverse_q and max_q need to be precomputed, to be efficient, but since you are using a sieve, I assume this should not be an obstacle.