Counting pairs of relatively prime numbers Problem from Hacker Earth:


  Inverted GCD:
  
  Given an array a of \$N\$ numbers , you have to find the number of pair of indices \$i\$ and \$j\$ that satisfy the following relation:
  
  
  \$i &lt; j\$ 
  \$a_i &gt; a_j\$
  gcd\$(a_i , a_j) = 1\$
  
  
  Input:
  
  The first line of the input contains a single integer N - denoting the size of the array.
  The next line contains N space separated integers denoting the array a.
  
  Output:
  Output a single integer - the answer to the problem stated above.
  
  Constraints:
  \$1 \le N \le 10^5\$
  \$1 \le a_i \le 10^5\$


For this problem I've submitted the following solution in Python:

from fractions import gcd
n = int(input())
count = 0;
a = [int(x) for x in input().split()]
for i in range(n):
    for j in range(n):
        if(a[i]&gt;a[j] and i&lt;j and gcd(a[i] , a[j]) == 1):
            count +=1;
print(count)


I'm not sure what is the complexity of gcd() function in Python, but it clear due to two for loop that the time complexity is minimum \$O(n^2)\$.I can't understand how to reduce it. Because here we have to check these conditions between every pair of integers. I've seen such question here but there is no condition for gcd(), so can't relate these two questions. I'm newbie to Python.
tl;dr TLE beyond input size 40000 to 50000 (10 seconds for input size 100000)

Conditions (1) and (2) are transitive but independent and one-sided. This means they can only prune half of all candidates on average each, and so their combination only reduces the average candidate set to one fourth (i.e. half of a half). Such a small constant factor is not worth it to invest in (as in 'devising a special data structure'), unless overall performance is already close to a passing grade. Asymptotic performance doesn't change at all, of course. Hence it is best to take the doubling in performance that comes from index pruning (condition (2), which is essentially free) and test condition (1) as pairs are processed.

Condition (3) induces a squarish relation all on its own already, and it is not transitive. This means that there is no efficiency gain from precomputing this relation, and memoisation is likely to be completely ineffective in the face of random inputs.

The evaluation of the GCD for a pair of numbers takes time logarithmic in the smaller of the two inputs. This time can be cut by a considerable constant factor with a bit of precomputation, by computing integer-sized 'bit vectors' representing the presence of a handful of the smallest primes in the (radical) factorisation of a number. 

This way a simple binary AND can prove the absence of shared small factors and thus avoid the computation of the actual GCD in about 89% of all cases (assuming factorisation by the first 32 primes, i.e. up to 131). In the remaining 11%, only the factorisation 'remainders' need to be GCDed, i.e. what remains after dividing out the small factors/radicals represented by the bit vectors. This saves some iterations of Euclid's algorithm.

This can give you an overall speed boost of an order of magnitude or more. On top of the factor 2 speed gain from index pruning mentioned by coderodde and some general efficiency tuning this might just be enough to beat the time limit...

Note: one might be tempted to extend the 'bit vectors' from 32 bits to 64 bits to utilise the bigger machine words on 64-bit machines, but this is not worth it. It doubles the up-front factorisation effort but it improves the rejection rate only from 88.7% to 90.3%, i.e. by a completely negligible 1.8%. 

Here's a small table with the pass/reject rates for bit masks involving the given number of the first primes:

 8 primes: pass 17,102%, rej. 82,898%                 # up to 4 extra radicals
12 primes: pass 14,872%, rej. 85,128% // boost 2,690% # up to 3 extra radicals
16 primes: pass 13,609%, rej. 86,391% // boost 1.484%
32 primes: pass 11,300%, rej. 88,700% // boost 2,673%
64 primes: pass  9,683%, rej. 90,317% // boost 1,823% # up to 1 extra radical


The 'up to X extra radicals' bit says how many different prime factors can remain in the factorisations of numbers not exceeding 100,000 after the the given number of the first small primes has been factored out. 

P.S.: the bounds on the number of numbers is the same as the bound on the size of the numbers, which means that the factorisations can be pre-computed wholesale without any divisions (by borrowing the logic behind the Sieve of Eratosthenes). The partial factoring for the bit vectors takes less than a millisecond, and factoring the input remainders takes less than 10 milliseconds for 100,000 inputs. Hence the stuff can always be computed wholesale and it isn't even worth it to have different code paths for small and big input batches.

--- UPDATE ---

For getting a better handle on the problem I needed a closer fix on some numbers, and so I coded things up in C++. It turned out that computing the gcd of the remainders (after factoring out the small primes used in the 'bit vectors') is orders of magnitude too slow, even though it occurs only for 0.35% of all candidates. This necessitated replacing the partial factorisation remainders with lists of the (up to two) remaining radicals.

'Radical' Optimisation

For numbers less than 107113 = 43 * 47 * 53, factoring out the 13 smallest primes leaves at most two distinct prime factors ('radicals'). Moreover, only factors up to half the limit need be considered since the only numbers with bigger prime factors are the primes, which are coprime to all numbers but themselves and self-pairings get excluded via criteria (1) and (2) both. This means that 2 x 16 bits = 32 bits are sufficient for storing the 'list' of remaining radicals for a number. 

This size reduction is beneficial for L1 cache utilisation and it would allow to stuff both radicals into a single 32-bit integer; however, it collides with an important optimisation of the intersection operation. This optimisation is the reduction of branchings to the absolute minimum of two to three comparisons, by always treating the two-element lists as full and by replacing all but one of the remaining non-comparison branches with arithmetic (poor man's CMOV). This makes things easier for the CPU's branch prediction logic and reduces the amount of incurred branch mis-prediction penalties, but it requires the use of sentinel values for 'missing' radicals. With the only safe sentinel for a number a[i] being a[i] itself this means that the slots for the radicals must be able to hold full-sized numbers. 

Elimination of one particularly nasty branch required adding a third radical slot (which always holds a sentinel value) but it resulted in a speedup of 10%. Without this explanation it might seem odd that the code uses three slots even though there can be at most two radicals. First, the data:

uint16_t sieve[1 + LIMIT];  // bit masks of partial factorisations of 1 .. LIMIT
uint16_t v[1 + LIMIT];      // v[i] = sieve[a[i]]
uint32_t f[1 + LIMIT][3];   // remaining factors (radicals) for each input (3rd is sentinel)
uint8_t skip[1 + LIMIT];    // for each a[j]: how many successors are not less


The code for generating these tables is elementary and of little interest here (especially as it takes less than 10 milliseconds even in C#). Hence I'm skipping it in the interest of brevity. Here is the meat of the code, the function that does the counting as specified by the task:

uint64_t count_matches (uint32_t const a[], uint32_t count)
{
    uint64_t result = 0;

    prepare_tables(a, count);

    for (uint32_t i = 0; i &lt; count; ++i)
    {
        uint32_t a_i = a[i];
        unsigned v_i = v[i];

        for (uint32_t j = i + 1; j &lt; count; ++j)
        {
            if ((v_i &amp; v[j]) == 0)
            {
                if (a_i &gt; a[j])
                {
                    if (uint32_t d = f[i][0] - f[j][0])
                    {
                        uint32_t di = (d &gt;&gt; 31) &amp; 1;
                        uint32_t ci = di;
                        uint32_t cj = di ^ 1;

                        if ((d = f[i][ci] - f[j][cj]) != 0)
                        {
                            di = (d &gt;&gt; 31) &amp; 1;
                            ci += di;
                            cj += di ^ 1;

                            if ((d = f[i][ci] - f[j][cj]) != 0)
                                ++result;
                        }
                    }
                }
                else j += skip[j];
            }
        }
    }

    return result;
}


On my laptop this clocks a bit less than 19 seconds for random inputs of size 100,000, compared to 22 seconds for a 'branchy' version of the intersection test (which has no use for the third slot):

                    int32_t d = int32_t(f[i][0] - f[j][0]);

                    if (d &gt; 0)
                    {
                        d = int32_t(f[i][0] - f[j][1]);

                        result += d &gt; 0 || d &amp;&amp; f[i][1] != f[j][1];
                    }
                    else if (d != 0)
                    {
                        d = int32_t(f[i][1] - f[j][0]);

                        result += d &lt; 0 || d &amp;&amp; f[i][1] != f[j][1];
                    }


Of course, in Python the optimisation is bound to get killed by the interpreter overhead, because the code contains a lot of explicit little operations instead of pushing processing into the engine. But in any case it should be beneficial to utilise the fact the 'radical factor set intersection' involves only lists of at most two elements here instead of general sets.

Skip Counts

The skip count thing hasn't been mentioned yet - it's a small optimisation that costs virtually nothing and improves speed by about 10%. It utilises the fact that a[i] &lt;= a[j] &amp;&amp; a[j] &lt;= a[j + k] implies a[i] &lt;= a[j + k], and storing/precomputing this is cheap because it involves only a[j] and its context in the input sequence instead of depending on pairings of a[i] and a[j]. 

Half of all skip counts are 0, of course, and half of the rest is 1. The rest peters off quickly because the probability of a skip count having value k is 1 / ((k+1)*(k+2)), which means that skip counts of four bits or so already realise most of the performance gain that can be gotten this way. 

Four bits would probably be optimum if one wanted to stuff several different values - like input number, small factor bit vector and skip count - into a single integer-sized variable, especially in languages where array indexings carry a heavy cost (Python, Java, C#). In my sample code I used a separate array for simplicity, and hence allocated one whole byte for each skip count.

Numbers and Summary

The reason I sat down and got coding was that I wanted to get a fix on some numbers, in order to get a better handle on the problem and a better feeling for the lay of the land. Here are those numbers. 

For a random input vector of 100,000 numbers - the limit given in the task description - and enumerating only the 4,999,950,000 pairs of numbers with i &lt; j, it took the following amounts of time to count the following things:


11.0 s for 2512386875 inversions
21.5 s for 3029211382 coprime pairs
17.3 s for 1527936260 candidates (a[i] &gt; a[j] and no common factors among the first 16 primes)
15.7 s for 1527936260 candidates (dito, with byte-sized skip counts)
85.7 s for 1522549277 matches (with gcd applied to remainders)
18.6 s for 1522549277 matches (with list intersection for remainder radicals)


The most important message here is that the number of inversions (i.e. criterion (2) in the task description), the number of coprime pairs (criterion (3)) and the number of results to be counted are all the same order of magnitude as the number of pairs with i &lt; j. The percentages are 50.2%, 60.6% and 30.5% respectively. 

Also, just counting the number of inversions takes more than five times as long as the two seconds admissible as per the task description, even though my laptop is likely quite a bit faster than the hardware used by the robot judge.

Even if I found some fiendishly clever algorithm or data structure which can prune everything that doesn't fulfil the task conditions then the time would still be 30.5% * 19 s which is about 6 seconds. In other words, on a maximum input file the code would still take three times as long as the time limit.

For comparison, counting the matches for a vector of 1000 random inputs - as per coderodde's 'lightened' specification - takes less than 2 milliseconds with this code. I'd hate to think how long coderodde's solution would take on a maximum input file as per the actual task description; extrapolation from the given results puts it at more than half an hour.

In the first two paragraphs of this post I argumented why the intersection of the task criteria must leave a number of pairings that is still square in the number of inputs; now we know that the actual number is about 30% of 'j-pruned' pairs or 15% of all input pairs. Moreover, the C++ timings suggest that counting each result pair individually is simply not an option, since the optimisation potential that's still left is far less than an order of magnitude (which is why I didn't bother to create an optimised version).

Any admissible solution must count whole swathes of pairs in one go instead of processing each pair individually, but at the moment I cannot offer any specific advice as to how one might accomplish that. One seedling of an idea might lie in the fact that any prime must a priori be coprime to any input number other than itself, so that only criteria (1) and (2) need be applied for counting it. Similar considerations apply to sets of input numbers that have disjoint sets of prime factors.

I took a peek at the submission table at HackerEarth, and the only successful submissions I could see were in C++ with timings between 3.3 seconds and 12 seconds, with codes sizes ranging from 10KB to several hundred KB. This suggests that real solutions might be quite involved, compared to the simple code that has been discussed in this topic so far. I did not look at any of the code, however, because I intend to have a stab at the problem myself first...

UPDATE #2

250,000,000 factor set intersections per second

For an input vector of 10^5 values the current algorithm has to filter 0.5 * (10^5)^2 pairs of numbers - that's about 5 billion (American-style). Performance can be roughly doubled by looking only at inversions and by removing all branches but one from the set intersection code.

The standard algorithm for inversion counting by merge sorting divides the input into non-inverted and inverted runs and tallies the latter:

long count_inversions_m1 (List&lt;int&gt; values)
{
    long inversions = 0;
    int n = values.Count;
    var a = values.ToArray();
    var b = new int[n];

    for (int w = 1; w &lt; n; w &lt;&lt;= 1)
    {
        for (int beg = 0; beg &lt; n; beg += w + w)
        {
            int mid = Math.Min(beg + w, n);
            int end = Math.Min(mid + w, n);

            for (int i = beg, j = mid, k = beg; k &lt; end; ++k)
            {
                if (i &lt; mid &amp;&amp; (j &gt;= end || a[i] &lt;= a[j]))
                {
                    b[k] = a[i++];
                }
                else // a[k] &gt; a[j] for all k in [i .. mid - 1]
                {
                    inversions += mid - i;
                    b[k] = a[j++];
                }
            }
        }

        var t = a;
        a = b;
        b = t;
    }

    return inversions;
}


In the a[i] &gt; a[j] branch, the run from a[i] to a[mid - 1] is inverted with respect to a[j]; all inversions in the input pass that branch eventually. Replacing the line

 inversions += mid - i;


with a small loop that performs factor set intersection on the values involved gives an algorithm that is a tad faster than the earlier one, because it never has to look at any non-inverted pairs. Simply dropping the intersection code from the earlier version in there would work, but it can be made faster by using branch-less comparison. 

That is, instead of branching based on the result of comparing factors, a value is computed that is 0 if the two lists of non-bitmapped prime factors (at most two elements each) have some value in common, and 1 if they do not. I chose to abuse the fact that, given a pair of unequal values a and b, either a - b or b - a must be negative but no difference is negative if the values are equal. Hence the high bit of (a - b) | (b - a) is 1 if the values are different, and 0 if they are equal (as long as the difference does not exceed INT32_MAX).

A further speedup comes from ‘compressing’ factors down to 16 bits. This is possible because the factors are only ever compared for equality, and any unique mapping of values will do just as well for that. There are only 9592 primes up to 100000, and so there are many viable schemes for assigning each prime a unique 16-bit proxy value. I chose p &lt; UINT_MAX ? p : (ushort)(p + 3) for simplicity. This ‘compression’ from 32 bits to 16 bits leads to better cache utilisation and reduces register pressure for 32-bit architectures (since a whole pair fits into a single machine word).

Here's the final set intersection code, ready for replacing the inversion branch in the above procedure (C++ version shown, but the C# is similar):

uint32_t a_j = b[k] = a[j++];                // x 815,024
unsigned f_j = factor_bits[a_j];
unsigned bj0 = big_factors[a_j][0];
unsigned bj1 = big_factors[a_j][1];

for (uint32_t l = i; l &lt; mid; ++l)
{
    uint32_t a_i = a[l];                     // x 2,512,386,875

    if ((f_j &amp; factor_bits[a_i]) == 0)   
    {
        unsigned bi0 = big_factors[a_i][0];  // x 1,527,936,260
        unsigned di0 = ((bi0 - bj0) | (bj0 - bi0)) &amp; ((bi0 - bj1) | (bj1 - bi0));
        unsigned bi1 = big_factors[a_i][1];
        unsigned di1 = ((bi1 - bj0) | (bj0 - bi1)) &amp; ((bi1 - bj1) | (bj1 - bi1));

        result += (di0 &amp; di1) &gt;&gt; 31;         // 1,522,549,277 x 1 vs 5,386,983 x 0
    }
}


The 16-bit factor bitmaps are kept in a separate array to improve cache efficiency, because the bitmap AND must occur for every inversion pair whereas comparison of the up to two extra factors occurs less often. I've annotated some code lines with execution counts for an input of 100,000 random values, to put things a bit into perspective.

The resulting code comes closer to the goal, but not close enough. The C++ version still needs 10 seconds for a maximum-size (non-degenerate) input, and for the C# edition it's 12.5 seconds. I submitted the C# but, as expected, it failed all (non-degenerate) test cases where the input size exceeded 40,000 significantly. The performance target seems almost within reach, but it seems hopeless to try squeezing another five-fold performance increase out of the current approach.

I've chosen to complete the write-up despite the failure because some of the techniques may be useful for other challenges. Thanks to Kaushal28 for bringing this neat problem to our attention - I had a lot of fun trying to solve it and I learnt a lot along the way.