Faster algorithm to find how many numbers are not divisible by a given set of numbers

Question

I am trying to solve an online judge problem: http://opc.iarcs.org.in/index.php/problems/LEAFEAT

The problem in short:

If we are given an integer L and a set of N integers s1,s2,s3..sN, we have to find how many numbers there are from 0 to L-1 which are not divisible by any of the 'si's.

For example, if we are given, L = 20 and S = {3,2,5} then there are 6 numbers from 0 to 19 which are not divisible by 3,2 or 5.

L <= 1000000000 and N <= 20.

I used the Inclusion-Exclusion principle to solve this problem:

/*Let 'T' be the number of integers that are divisible by any of the 'si's in the 
given range*/

for i in range 1 to N
  for all subsets A of length i
    if i is odd then:
      T += 1 + (L-1)/lcm(all the elements of A)
    else
      T -= 1 + (L-1)/lcm(all the elements of A)
return T

Here is my code to solve this problem

#include <stdio.h>

int N;
long long int L;
int C[30];

typedef struct{int i, key;}subset_e;
subset_e A[30];
int k;

int gcd(a,b){
int t;
    while(b != 0){
            t = a%b;
            a = b;
            b = t;
    }

    return a;
}

long long int lcm(int a, int b){
    return (a*b)/gcd(a,b);
}

long long int getlcm(int n){
  if(n == 1){
    return A[0].key;
  }
  int  i;
  long long int rlcm = lcm(A[0].key,A[1].key);
  for(i = 2;i < n; i++){
    rlcm = lcm(rlcm,A[i].key);
  }
  return rlcm;
}

int next_subset(int n){
  if(k == n-1 && A[k].i == N-1){
    if(k == 0){
      return 0;
    }
    k--;
  }
  while(k < n-1 && A[k].i == A[k+1].i-1){
    if(k <= 0){
      return 0;
    }
    k--;
  }
  A[k].key = C[A[k].i+1];
  A[k].i++;
  return 1;
}

int main(){
  int i,j,add;
  long long int sum = 0,g,temp;
  scanf("%lld%d",&L,&N);
  for(i = 0;i < N; i++){
    scanf("%d",&C[i]);
  }
  for(i = 1; i <= N; i++){
    add = i%2;
    for(j = 0;j < i; j++){
      A[j].key = C[j];
      A[j].i = j;
    }
    temp = getlcm(i);
    g = 1 + (L-1)/temp;
    if(add){
      sum += g;
    } else {
      sum -= g;
    }
    k = i-1;
    while(next_subset(i)){
      temp = getlcm(i);
      g = 1 + (L-1)/temp;
      if(add){
        sum += g;
      } else {
        sum -= g;
      }
    }
  }
  printf("%lld",L-sum);
  return 0;
}

The next_subset(n) generates the next subset of size n in the array A, if there is no subset it returns 0 otherwise it returns 1. It is based on the algorithm described by the accepted answer in this stackoverflow question.

The lcm(a,b) function returns the lcm of a and b. The get_lcm(n) function returns the lcm of all the elements in A. It uses the property : LCM(a,b,c) = LCM(LCM(a,b),c)

When I submit the problem on the judge it gives my a 'Time Limit Exceeded'. If we solve this using brute force we get only 50% of the marks.

As there can be upto 2^20 subsets my algorithm might be slow, hence I need a better algorithm to solve this problem.

EDIT:

After editing my code and changing the function to the Euclidean algorithm, I am getting a wrong answer, but my code runs within the time limit. It gives me a correct answer to the example test but not to any other test cases; here is a link to ideone where I ran my code, the first output is correct but the second is not.

Is my approach to this problem correct? If it is then I have made a mistake in my code, and I'll find it; otherwise can anyone please explain what is wrong?

Answer 1

You could also try changing your lcm function to use the Euclidean algorithm.

int gcd(int a, int b) {
    int t;

    while (b != 0) {
        t = b;
        b = a % t;
        a = t;
    }

    return a;
}

int lcm(int a, int b) {
    return (a * b) / gcd(a, b);
}

At least with Python, the speed differences between the two are pretty large:

>>> %timeit lcm1(103, 2013)
100000 loops, best of 3: 9.21 us per loop
>>> %timeit lcm2(103, 2013)
1000000 loops, best of 3: 1.02 us per loop

Answer 2

Typically, the lowest common multiple of a subset of k of the s_i will exceed L for k much smaller than 20. So you need to stop early.

Probably, just inserting

if (temp >= L) {
    break;
}

after

while(next_subset(i)){
  temp = getlcm(i);

will be sufficient.

Also, shortcut if there are any 1s among the s_i, all numbers are divisible by 1.

I think the following will be faster:

unsigned gcd(unsigned a, unsigned b) {
    unsigned r;
    while(b) {
        r = a%b;
        a = b;
        b = r;
    }
    return a;
}

unsigned recur(unsigned *arr, unsigned len, unsigned idx, unsigned cumul, unsigned bound) {
    if (idx >= len || bound == 0) {
        return bound;
    }
    unsigned i, g, s = arr[idx], result;
    g = s/gcd(cumul,s);
    result = bound/g;
    for(i = idx+1; i < len; ++i) {
        result -= recur(arr, len, i, cumul*g, bound/g);
    }
    return result;
}

unsigned inex(unsigned *arr, unsigned len, unsigned bound) {
    unsigned i, result = bound, t;
    for(i = 0; i < len; ++i) {
        result -= recur(arr, len, i, 1, bound);
    }
    return result;
}

call it with

unsigned S[N] = {...};
inex(S, N, L-1);

You need not add the 1 for the 0 anywhere, since 0 is divisible by all numbers, compute the count of numbers 1 <= k < L which are not divisible by any s_i.

Answer 3

Create an array of flags with L entries. Then mark each touched leaf:

for(each size in list of sizes) {
    length = 0;
    while(length < L) {
        array[length] = TOUCHED;
        length += size;
    }
}

Then find the untouched leaves:

for(length = 0; length < L; length++) {
    if(array[length] != TOUCHED) { /* Untouched leaf! */ }
}

Note that there is no multiplication and no division involved; but you will need up to about 1 GiB of RAM. If RAM is a problem the you can use an array of bits (max. 120 MiB).

This is only a beginning though, as there are repeating patterns that can be copied instead of generated. The first pattern is from 0 to S1*S2, the next is from 0 to S1*S2*S3, the next is from 0 to S1*S2*S3*S4, etc.

Basically, you can set all values touched by S1 and then S2 from 0 to S1*S2; then copy the pattern from 0 to S1*S2 until you get to S1*S2*S3 and set all the S3's between S3 and S1*S2*S3; then copy that pattern until you get to S1*S2*S3*S4 and set all the S4's between S4 and S1*S2*S3*S4 and so on.

Next; if S1*S2*...Sn is smaller than L, you know the pattern will repeat and can generate the results for lengths from S1*S2*...Sn to L from the pattern. In this case the size of the array only needs to be S1*S2*...Sn and doesn't need to be L.

Finally, if S1*S2*...Sn is larger than L; then you could generate the pattern for S1*S2*...(Sn-1) and use that pattern to create the results from S1*S2*...(Sn-1) to S1*S2*...Sn. In this case if S1*S2*...(Sn-1) is smaller than L then the array doesn't need to be as large as L.