Need Explanation of this code

Question

The code outputs number of pairs of (i,j) which satisfy the condition

(2^j-1) % (2^i-1) == 0

where the

1<=i<j<=n

n is the number user input under which number of (i,j) pairs are to be found. Code works perfectly, but what is the logic behind this code is hard to understand.

P.S: t is a variable which will let user to input numbers more than one at a time.

    #include<stdio.h>
    #include<math.h>
    int main()
    {   
        int t;
        long n,sum,ans; 
        scanf("%d",&t);
        while(t--)
        {
            scanf("%ld",&n);
            int nrt=(int)sqrt(n);
            sum=0;
            for(int i=1;i<=nrt;i++)
            {
                 sum+=n/i;
            }
            ans=2*sum-nrt*nrt-n;
            printf("%ld\n",ans);
        }
    return 0;
    }

Answer 1

It's a very interesting problem. If you tried some small inputs, you will get a rough idea on the code.

I used a very simple code to generate all valid pairs when n = 10, and here is what I got:

1 2
1 3
1 4
1 5
1 6
1 7
1 8
1 9
1 10
2 4
2 6
2 8
2 10
3 6
3 9
4 8
5 10

Surprise? We could see a very obvious pattern here: when i, j satisfies j = k * i, where k is an integer and j = k * i < n, i, j is a valid pair. It's totally nothing to do with the original equation, only depends on n.

Actually this is not surprise, since (2^(nk) - 1) = ((2^k)^n - 1) = (a^n - 1), where a = 2^k, therefore we could apply factoring rule, which gives (a^n - 1) = (a - 1)(a^(n - 1) + a^(n - 2) + .. + 1), and thus could be factorized by (a - 1), i.e. (2^(nk) - 1) % (2^k - 1) == 0.

Now the problem becomes how to count this number efficiently. By the condition, we have j > i. Previous we know j = k * i. Therefore, k must be in the range of [2, n / i]. For each i, we have totally (n / i) - 2 + 1 = (n / i) - 1 valid choices for k. Therefore, the total valid pairs will be sigma((n / i) - 1, 1 <= i <= n).

As for how to transform the equation into the code you gave, please refer to @MOehm's answer.

Answer 2

Let's take a brute-force approach to the problem and print the results*:

 #############################   2^^1 -1 == 1
  -#-#-#-#-#-#-#-#-#-#-#-#-#-#   2^^2 -1  == 3
   --#--#--#--#--#--#--#--#--#   2^^3 -1  == 7
    ---#---#---#---#---#---#--   2^^4 -1  == 15
     ----#----#----#----#----#   2^^5 -1  == 31
      -----#-----#-----#-----#   2^^6 -1  == 63
       ------#------#------#--   2^^7 -1  == 127
        -------#-------#------   2^^8 -1  == 255
         --------#--------#---   2^^9 -1  == 511
          ---------#---------#   2^^10 -1 == 1023
           ----------#--------   2^^11 -1 == 2047
            -----------#------   2^^12 -1 == 4095
             ------------#----   2^^13 -1 == 8191
              -------------#--   2^^14 -1 == 16383
               --------------#   2^^15 -1 == 32767
                --------------   2^^16 -1 == 65535
                 -------------   2^^17 -1 == 131071
                           ...   ...

A hash mark denotes the cases where your condition is met. A nice pattern emerges: Every one of your numbers is divisible by 1, every other is divisible by 3, every third is divisible by 7 and so on. Every ith number is divisible by 2^^i - 1.**

With this insight, we can code your function as:

int f(int n)
{
    int sum = 0;
    int i;

    for (i = 1; i <= n; i++) sum += (n - i) / i;

    return sum;
}

We can substitute (n - i) / i with n / i - 1 and move the common subtrahend -1 into the return value:

int g(int n)
{
    int sum = 0;
    int i;

    for (i = 1; i <= n; i++) sum += n / i;

    return sum - n;
}

Now let's look at the sum ∑(1, n: n / i). For example:

∑(i = 1, 9: 9 / i) = 9 + 4 + 3 + 2 + 1 + 1 + 1 + 1 + 1

We can get at the same sum by looking at it from right to left and by counting how often each summand occurs:

∑(i = 1, 9: 9 / i) = 5*1 + 1*2 + 1*3 + 1*4 + 1*9

We can get at this representation easily:

∑(i = 1, n: n / i) = ∑(1, n: i * (n / i - n / (i + 1))

This really is just another way to write the sum; you can see that by grouping the summands differently, so that they share the same denominator:

∑(i = 1, N: i * (n / i - n / (i + 1))
    = n + ∑(i = 1, n: ((i + 1) - i)  * n / (i + 1))
    = n + ∑(i = 1, n: n / (i + 1)) - (N + 1) * (n / (N + 1))
    = n + ∑(i = 2, n + 1: n / i) - c
    = ∑(i = 1, n: n / i) - c

The additional term c = (N + 1) * (n / (N + 1)) is a correction term because only half of the term for i = n + 1 is used. When summing over the whole range, n / (n + 1) is zero and vanishes. It doesn't vanish when summing only part of the array, as we'll see later.

If we split the sum into a head and tail at s = sqrt(n), we get:

∑(i = 1, n: n / i) = ∑(i = 1, s: n / i) + ∑(s + 1, n: n / i)

Let's represent the head in the original fashion and the tail in the "count the summands" fashion, for example:

∑(i = 1, 9: 9 / i) = (9 + 4 + 3)   +   (5*1 + 1*2)

For any n:

∑(i = 1, n: n / i) 
    = ∑(i = 1, s: n / i) + ∑(1, s - 1: i * (n / i - n / (i + 1))
    = ∑(i = 1, s: n / i) + ∑(1, s: n / i) - s * (n / s)

All divisions are integer divisions (and that's why there have to be parentheses sometimes) and n / s == s, so:

∑(1, n: n / i) = ∑(i = 1, s: n / i) + ∑(i = 1, s: n / i) - s * (n / s)
               = 2 * ∑(i = 1, s: n / i) - s²

And that yields your original function:

int h(int n)
{
    int nrt = sqrt(n);
    int sum = 0;
    int i;

    for(i = 1; i <= nrt; i++) sum += n/i;

    return 2 * sum - nrt * nrt - n;
}

where ∑(1, n: n / i) in g above has been replaced with 2 * ∑(i = 1, s: n / i) - s².

*) I've stolen D's power operator ^^ here, in order not to confuse the old C buffs who take ^ at face value, i.e. as xor.

**) I don't know, why the pattern shows. There is probably a good explanation, but for now, I trust my pattern-matching skills. Blindly. Edit @nevets's answer has the explanation for this pattern.

Answer 3

The variable i works from 1 to nrt, nrt being the square root of n converted explicitily into an int value. Each time the loop works, sum is added by the result of (n/i). The code then prints the ans (long type) with is calculated as (twice the sum-nrt square-n).