Number theory algorithms. Most divisors on the segment

Question

I'm looking for an efficient algorithm to solve the following problem. Let d(n) denote number of positive divisors of n where n is positive integer. We're given some 1 <= a <= b <= 10^18 and the task is to find maximum value of d on segment [a..b] and (may be we need more complex algorithm for this part) to find number which maximizes value of d.

Some time ago i found the following code in free access: http://ideone.com/qvxPj

unsigned long long n, res;
int p, primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 51, 53, 59, 61, 67, 71};

unsigned long long mul(unsigned long long a, unsigned long long b){
    unsigned long long res = 0;

    while (b){
        if (b & 1LL) res = (res + a);
        if (res >= n) return 0;
        a = (a << 1LL);
        b >>= 1LL;
    }

    return res;
}

void backtrack(int i, int lim, unsigned long long val, unsigned long long r){
    if (r > res) res = r;
    if (i == p) return;

    int d;
    unsigned long long x = val;

    for (d = 1; d <= lim; d++){
        x = mul(x, primes[i]);
        if (x == 0) return;
        backtrack(i + 1, d, x, r * (d + 1));
    }
}

int main(){    
    p = sizeof(primes) / sizeof(int);

    while (scanf("%llu", &n) != EOF){
        res = 0;
        backtrack(0, 100, 1, 1);
        printf("Maximum number of divisors of any number less than %llu = %llu\n", n, res);
    }
    return 0;
}

I'll be very glad if someone explains me how it works because (as for me) this program runs extremely fast.

Thanks in advance for any help.

Answer 1

It iterates over all numbers like this:

num = P1^D1 * P2^D2 * P3^D3 * ... * Ps^Ds
constraints:
  Pi <= 71
  1 <= Di <= 100
  sequence (Pi) is a sorted list of first s primes
  sequence (Di) is nonincreasing
  num <= n

Let's check the first constraint. Suppose that the minimal optimal number has prime factor q > 71. If any prime p <= 71 is not used in this number, then we can replace q with p in same power. Obviously, the number of divisors will stay same, but the number would decrease -> contradiction. Then there is no unused prime less than 71. But product of all primes up to 71 is already so huge, that the number we consider must be larger than 64-bit n. That's not possible.

Now let's explain the second and the third constraints. Suppose that our minimal optimal number has a prime q in its factorization, but doesn't have some prime p, where p < q. Then we can replace q with p in same order, the number would have same number of divisors, but it would become less -> contradiction. It means the all the primes in factorization of the sought-for optimal (minimal) number must be exactly the first s primes. There can be no holes in the set of primes used. BTW, Di <= 100 is obvious, because even 2^100 does not fit 64-bit integer already.

Now we have to explain the fourth constraint. Suppose that D[i] < D[i+1] for some i. Then we can replace P[i]^D[i] * P[i+1]^D[i+1] with P[i]^D[i+1] * P[i+1]^D[i], and the number would become smaller. For instance, replace 5^2 * 7^3 with 5^3 * 7^2: number of divisors is same, but the result is smaller. Obviously, if we search the minimal optimal number, we can safely assume this condition too.

Now let's consider the code. mul is a small function that calculates product of a and b. It is calculated by a funny binary procedure. The main reason for this procedure is: if the product is greater than n, then the function returns 0. This procedure is only a guard against overflow that may happen otherwise.

Finally, we got to backtrack. This is a usual recursive search. val is the current number, r is its number of divisors, i shows the index of the prime we are going to add now, lim limits power of each prime to 100. At the very start you see update of current optimal answer (stored in res), and hard stopping condition (all primes used).

Then there is a loop that checks every power for the current prime number. It starts with power 1, since zero powers are forbidden. It maintains current number in x and multiplies it by Pi each iteration to increment power. If x becomes greater than n, it stops immediately. Finally, it calls itself in order to search for the next prime.