カーマイケル数の列挙

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

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

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

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

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

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

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

view plainprint?

1. map<int, int>factorize(int x) {  
2.   map<int, int> ret;  
3.   
4.   for (int i = 2; i*i <= x; i++) {  
5.       while (x % i == 0) {  
6.           ret[i]++;  
7.           x /= i;  
8.       }  
9.   }  
10.   
11.   if (x != 1)  
12.       ret[x] = 1;  
13.   
14.   return ret;  
15. }  
16.   
17. int mpow(int x, int p) {  
18.   int ret = 1;  
19.   
20.   while (p) {  
21.       if (p & 1)  
22.           ret = ret * x;  
23.       x = x * x;  
24.       p >>= 1;  
25.   }  
26.   
27.   return ret;  
28. }  
29.   
30. int gcd(int a, int b) {  
31.   if (b)  
32.       return gcd(b, a % b);  
33.   return a;  
34. }  
35.   
36. int lcm(int a, int b) {  
37.   return (long long int) a * b / gcd(a, b);  
38. }  
39.   
40. int carmichaelFunction(int x) {  
41.   map<int, int> fact = factorize(x);  
42.   
43.   int ret = 1;  
44.   for (map<int, int>::iterator itr = fact.begin(); itr != fact.end(); itr++) {  
45.       int p = itr->first;  
46.       int k = itr->second;  
47.   
48.       int lambda = (p > 2 || k < 3) ? mpow(p, k-1) * (p-1) : mpow(p, k-2);  
49.       ret = lcm(ret, lambda);  
50.   }  
51.   
52.   return ret;  
53. }  
54.   
55. int eulerFunction(int x) {  
56.   map<int, int> fact = factorize(x);  
57.   
58.   double ret = x;  
59.   for (map<int, int>::iterator itr = fact.begin(); itr != fact.end(); itr++)  
60.       ret *= (1 - 1./itr->first);  
61.   
62.   return ret + 0.5;  
63. }  
64.   
65. int main() {  
66.   
67.   for (int i = 2; i < 100000; i++) {  
68.       int lambda = carmichaelFunction(i);  
69.       int phi = eulerFunction(i);  
70.   
71.       if ((i-1) % lambda == 0 && phi != i-1) {  
72.           printf("%d\n", i);  
73.       }  
74.   }  
75.   
76.   return 0;  
77. }  

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...と続く。