find all Carmichael numbers in a given interval [a,b) - C

I am working in a math software with different features one of them to be to find all Carmichael numbers in a given interval [a,b)

This is my code, but I don't know if I have done it correctly or not cause I can't test it since the smallest Carmichael number is 560 which is too big for my pc to process.

#include <stdio.h>

int main() {

  unsigned int begin, end;

  printf("Write an int (begin):\n");
  scanf("%d", &begin);

  printf("Write an int (end):\n");
  scanf("%d", &end);

  int i;

  for( int i=begin; i<end; i++ ) {

    long unsigned int a_nr = i-1;

    int a[a_nr];

    for( int j=0; j<a_nr; j++ ) {
      a[j] = j;
    }

    unsigned long c_nr[a_nr];

    for( int k=0; k<a_nr; k++ ) {
      unsigned long current_c_nr;
      int mod;
      for( int l=0; l<i; l++ ) {
        current_c_nr= current_c_nr * a[k];
      }
      mod = current_c_nr%i;
      if( mod==a[k] && mod!=a[k] ) {
        c_nr[k] = i;
      }

    }

  }

  return 0;
}

If it is not correct, where is the mistake?

Thank you

P.S Overflow should be prevented.

Solution

When you say "This is my code, but I don't know if I have done it correctly or not cause I can't test it since the smallest Carmichael number is 560 which is too big for my pc to process" then the conclusion is -- you haven't done it correctly. You should be able to process 561 (560 must be a typo) in a small fraction of a second. Even if your algorithm is in principle correct, if it can't handle the smallest Carmichael number then it is useless.

n is Carmichael if and only if it is composite and, for all a with 1 < a < n which are relatively prime to n, the congruence a^(n-1) = 1 (mod n) holds. To use this definition directly, you need:

1) An efficient way to test if a and n are relatively prime

2) An efficient way to compute a^(n-1) (mod n)

For the first -- use the Euclidean algorithm for greatest common divisors. It is most efficiently computed in a loop, but can also be defined via the simple recurrence gcd(a,b) = gcd(b,a%b) with basis gcd(a,0) = a. In C this is just:

unsigned int gcd(unsigned int a, unsigned int b){
    return b == 0? a : gcd(b, a%b);
}

For the second point -- almost the worst possible thing you can do when computing a^k (mod n) is to first compute a^k via repeated multiplication and to then mod the result by n. Instead -- use exponentiation by squaring, taking the remainder (mod n) at intermediate stages. It is a divide-and-conquer algorithm based on the observation that e.g. a^10 = (a^5)^2 and a^11 = (a^5)^2 * a. A simple C implementation is:

unsigned int modexp(unsigned int a, unsigned int p, unsigned int n){
    unsigned long long b;
    switch(p){
        case 0:
            return 1;
        case 1:
            return a%n;
        default:
            b = modexp(a,p/2,n);
            b = (b*b) % n;
            if(p%2 == 1) b = (b*a) % n;
            return b;
        }
}

Note the use of unsigned long long to guard against overflow in the calculation of b*b.

To test if n is Carmichael, you might as well first test if n is even and return 0 in that case. Otherwise, step through numbers, a, in the range 2 to n-1. First check if gcd(a,n) == 1 Note that if n is composite then you must have at least one a before you reach the square root of n with gcd(a,n) > 1). Keep a Boolean flag which keeps track of whether or not such an a has been encountered and if you exceed the square root without finding such an a, return 0. For those a with gcd(a,n) == 1, compute the modular exponentiation a^(n-1) (mod n). If this is ever different from 1, return 0. If your loop finishes checking all a below n without returning 0, then the number is Carmichael, so return 1. An implementation is:

int is_carmichael(unsigned int n){
    int a,s;
    int factor_found = 0;
    if (n%2 == 0) return 0;
    //else:
    s = sqrt(n);
    a = 2;
    while(a < n){
        if(a > s && !factor_found){
            return 0;
        }
        if(gcd(a,n) > 1){
            factor_found = 1;
        }
        else{
            if(modexp(a,n-1,n) != 1){
                return 0;
            }
        }
        a++;
    }
    return 1; //anything that survives to here is a carmichael
}

A simple driver program:

int main(void){
    unsigned int n;
    for(n = 2; n < 100000; n ++){
    if(is_carmichael(n)) printf("%u\n",n);
    }
    return 0; 
}

output:

C:\Programs>gcc carmichael.c

C:\Programs>a
561
1105
1729
2465
2821
6601
8911
10585
15841
29341
41041
46657
52633
62745
63973
75361

This only takes about 2 seconds to run and matches the initial part of this list.

This is probably a somewhat practical method for checking if numbers up to a million or so are Carmichael numbers. For larger numbers, you should probably get yourself a good factoring algorithm and use Korseldt's criterion as described in the Wikipedia entry on Carmichael numbers.