カーマイケル数の列挙

カーマイケル数を列挙してみる。カーマイケル数といえば、フェルマーテストの偽陽性を引き起こすクセモノ。あらかじめカーマイケル数を列挙しておけば、フェルマーテストと併用することで、正確な素数判定ができる。

いろいろな関数を使うことになった。

• 最大公約数
• 最小公倍数
• 素因数分解
• 繰り返し二乗法

など基本的な数学系処理の練習にいいかも。。

さて、問題のカーマイケル数の部分だが、カーマイケル関数を使って、

となる最小のmを求める。このmをλ(n)と表記する。ここで、aとnはcoprimeである。カーマイケル数は、

となるnなので、(n-1)がλ(n)で割り切れるとき、nはカーマイケル数の候補となる。さらに、ここから素数を除く必要がある。素数の場合は、φ(n) = n-1となるので、それを除けばいい。逆に、カーマイケル数は、合成数なので、totient functionはn-1より小さくなるため、必ずφ(n) != n-1となり、この除去で偽陰性となることはない。

map<int, int>factorize(int x) {
  map<int, int> ret;

  for (int i = 2; i*i <= x; i++) {
      while (x % i == 0) {
          ret[i]++;
          x /= i;
      }
  }

  if (x != 1)
      ret[x] = 1;

  return ret;
}

int mpow(int x, int p) {
  int ret = 1;

  while (p) {
      if (p & 1)
          ret = ret * x;
      x = x * x;
      p >>= 1;
  }

  return ret;
}

int gcd(int a, int b) {
  if (b)
      return gcd(b, a % b);
  return a;
}

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

int carmichaelFunction(int x) {
  map<int, int> fact = factorize(x);

  int ret = 1;
  for (map<int, int>::iterator itr = fact.begin(); itr != fact.end(); itr++) {
      int p = itr->first;
      int k = itr->second;

      int lambda = (p > 2 || k < 3) ? mpow(p, k-1) * (p-1) : mpow(p, k-2);
      ret = lcm(ret, lambda);
  }

  return ret;
}

int eulerFunction(int x) {
  map<int, int> fact = factorize(x);

  double ret = x;
  for (map<int, int>::iterator itr = fact.begin(); itr != fact.end(); itr++)
      ret *= (1 - 1./itr->first);

  return ret + 0.5;
}

int main() {

  for (int i = 2; i < 100000; i++) {
      int lambda = carmichaelFunction(i);
      int phi = eulerFunction(i);

      if ((i-1) % lambda == 0 && phi != i-1) {
          printf("%d\n", i);
      }
  }

  return 0;
}

カーマイケル数は、561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361...と続く。