Iterating over numbers with many divisors

Question

I would like to iterate over the numbers with more than $D$ divisors in a large range $[x, x+N]$. Current values I'm working with are $D=626$ and $x\approx N\approx10^{11}.$

At the moment I'm using a factorization sieve (PARI/GP's forfactored) to iterate through the whole interval and skipping anything with at most $D$ divisors. (The factorization sieve computes the number of divisors 'for free'.)

If $D$ were really big, I might be able to restrict to multiples of 2, or 6, or 12. But even odd numbers this size can have over a thousand divisors, so that's not workable. But it seems like much better should be possible, since numbers with this many divisors are quite rare. There should be some solution which loops over exponents for the first few (100? 10000?) primes, but I'm not sure what stopping criteria to use. Thoughts?

Answer 1

I have written a script for this in C++, which you can find at the bottom of this answer. It seems it runs for your use case in $\approx 75$ seconds. It should be faster for you than forfactored, please let me know if anyone finds any more optimisation that I can do to the code.

I took your idea of looping over exponents, but in a reverse way (largest prime to lowest prime); and my stopping criteria uses 'Highly composite numbers', i.e., numbers $n$ where $d(n)$, the number of divisors of $n$ , increases to a record. See https://oeis.org/A002182 and https://oeis.org/A002183 for more explanation.

My script finds numbers in $[1,N]$ with divisors more than $D$. The high level idea is:

• Run linear time sieve till $\left\lceil\sqrt N\right\rceil$, to find the prime numbers in $\left[2,\left\lceil\sqrt N\right\rceil\right]$. Let them be $p_1,p_2 \dots p_k$

• Call recursive function solve(k, 1, 1)

• Details of solve(ind, val, divisors) are:

  • If ind = 0 (terminating step) : OUTPUT(val) only if divisors $> D$

  • Else find the largest value of A002182 $\le \frac{N}{\text{val}}\ = $ A002182[i]. Also find the corresponding A002183[i]

  • If A002183[i] $\times$ divisors $>D$, run a loop $t$ for exponents $p_\text{ind}^0, p_\text{ind}^1,p_\text{ind}^2,\dots p_\text{ind}^s \le \frac{N}{\text{val}}$ and call solve(k-1, val * (p_ind)^t, divisors*(t+1)). This is the stopping criteria currently, and weeds out any possibilities which can't make it to more than $D$ divisors.

The runtime complexity of this is somewhere between $\Omega(\sqrt N)$ and $O(N)$, and definitely depends highly on $D$ (where on high $D$, its closer to $\Theta(\sqrt N)$), but I have no idea currently how it depends on that. I'm running some time tests to get some idea, will edit the post if I get it.

Here are the current results :

Rows represent $N$ with values :1000,10000,100000, ... 1000000000,10000000000

Columns represent D with values : 100,200,300,400, ... 1800,1900,2000

0.000014,0.000013,0.000012,0.000012,0.000016,0.000010,0.000013,0.000012,0.000011,0.000013,0.000012,0.000013,0.000012,0.000012,0.000016,0.000012,0.000011,0.000010,0.000012,0.000010,
0.000161,0.000020,0.000019,0.000018,0.000017,0.000019,0.000021,0.000020,0.000021,0.000022,0.000019,0.000020,0.000019,0.000018,0.000022,0.000019,0.000020,0.000025,0.000021,0.000020,
0.003108,0.000296,0.000038,0.000039,0.000034,0.000043,0.000033,0.000034,0.000034,0.000034,0.000038,0.000035,0.000039,0.000036,0.000032,0.000037,0.000036,0.000035,0.000038,0.000037,
0.059675,0.007385,0.001908,0.000523,0.000089,0.000092,0.000091,0.000089,0.000085,0.000086,0.000089,0.000090,0.000092,0.000082,0.000089,0.000086,0.000090,0.000086,0.000094,0.000091,
0.895572,0.144097,0.045533,0.014107,0.007312,0.002999,0.001725,0.000264,0.000252,0.000251,0.000271,0.000253,0.000263,0.000257,0.000261,0.000259,0.000256,0.000258,0.000257,0.000262,
11.268128,2.709907,0.881835,0.306602,0.179348,0.079676,0.050591,0.022058,0.015891,0.009958,0.006493,0.003759,0.003304,0.000740,0.000725,0.000734,0.000733,0.000755,0.000753,0.000750,
^C,41.459918,17.274902,5.737790,3.582669,1.678419,1.142178,0.549022,0.430282,0.279197,0.202219,0.113233,0.085952,0.068420,0.051472,0.028519,0.028648,0.018177,0.017932,0.011946,

Finally the code:

#include<bits/stdc++.h>
using namespace std;

int sqrtint(long long N, long long l, long long r) {
    if (l == r) return l;
    else if(((__int128)1) * ((l+r)/2) * ((l+r)/2) >= N) return sqrtint(N,  l, (l+r)/2);
    else return sqrtint(N, (l+r)/2+1, r);
}

int main() {
    clock_t start, end;

    long long N, D;
    cin>>N>>D;
    start = clock();

    int rt_N = sqrtint(N, 1, N);
    static vector<int> d(rt_N+1), prm;
    for(int i=0; i<=rt_N; i++) {
        d[i] = i;
    }

    for(int i=2; i<=rt_N; i++) {
        if(d[i] == i) {
            prm.push_back(i);
        }
        for(int j=0; j < prm.size() && i*prm[j] <= rt_N && prm[j] <= d[i]; j++) {
            d[i*prm[j]] = prm[j];
        }
    }

    long long A002182[] = {
        1,2,4,6,12,24,36,48,60,120,180,240,360,720,840,
        1260,1680,2520,5040,7560,10080,15120,20160,25200,
        27720,45360,50400,55440,83160,110880,166320,
        221760,277200,332640,498960,554400,665280,720720,
        1081080,1441440,2162160,2882880,3603600,4324320,
        6486480,7207200,8648640,10810800,14414400,17297280,
        21621600,32432400,36756720,43243200,61261200,
        73513440,110270160,122522400,147026880,183783600,
        245044800,294053760,367567200,551350800,698377680,
        735134400,1102701600,1396755360,2095133040,2205403200,
        2327925600,2793510720,3491888400,4655851200,5587021440,
        6983776800,10475665200,13967553600,20951330400,
        27935107200,41902660800,48886437600,64250746560,
        73329656400,80313433200,97772875200,128501493120,
        146659312800,160626866400,240940299600,293318625600,
        321253732800,481880599200,642507465600,963761198400
    };
    int A002183[] = {
        1,2,3,4,6,8,9,10,12,16,18,20,24,30,32,36,40,48,
        60,64,72,80,84,90,96,100,108,120,128,144,160,168,
        180,192,200,216,224,240,256,288,320,336,360,384,
        400,432,448,480,504,512,576,600,640,672,720,768,
        800,864,896,960,1008,1024,1152,1200,1280,1344,
        1440,1536,1600,1680,1728,1792,1920,2016,2048,2304,
        2400,2688,2880,3072,3360,3456,3584,3600,3840,4032,
        4096,4320,4608,4800,5040,5376,5760,6144,6720
    };

    struct call_stack_element{
        int index;
        long long val;
        int divisors;
    };

    static vector<call_stack_element> call_stack;
    static vector<long long> ans;
    call_stack.push_back((call_stack_element){((int)prm.size())-1, 1LL, 1});

    while(!call_stack.empty()) {
        auto tmp = call_stack.back();
        call_stack.pop_back();

        if(tmp.index == -1 && tmp.divisors > D) {
            ans.push_back(tmp.val);
        }
        else if(tmp.index >= 0) {
            int pntr = upper_bound(
                A002182, A002182 + sizeof(A002182)/sizeof(long long), N/tmp.val
            ) - A002182;
            if(A002183[pntr-1] * tmp.divisors > D || N/tmp.val >= 963761198400)
            {
                long long prime_power = 1LL;
                for(int exp=0; tmp.val*prime_power <= N; exp++) {
                    call_stack.push_back(
                        (call_stack_element){tmp.index-1, tmp.val*prime_power, tmp.divisors*(exp+1)}
                    );
                    prime_power *= prm[tmp.index];
                }
            }
        }
    }

    sort(ans.begin(), ans.end());
    for(auto it : ans) {
        cout<<it<<',';
    }
    end = clock();
    // cout<<fixed<<setprecision(6)<<((long double)(end-start))/CLOCKS_PER_SEC;

    return 0;
}