計算機プログラムの構造と解釈（SICP）学習メモ
蒋　いつ峰 2008/08/03   |  プログラミング基礎
第一章　練習1.21〜1.28 2008/08/03  |  PK：SICP 練習 1.21〜1.28

SICP第一章の練習1.21〜1.28、素数問題とフェルマ検査。


問題

練習1.21：smallest-divisor 方法で、199、1999、19999 の最少因子を求める。
練習1.22：上記 smallest-divisor で素数を求める時の時間を計る。
練習1.23：改良されたsmallest-divisor。
練習1.24：フェルマ検査。
練習1.25：expmod に、fast-expt を直接使う場合。
練習1.26：square を使わない時のコスト。
練習1.27：Carmichael 数。
練習1.28：Miller-Rabin検査（騙されないフェルマ検査）。


解答

各練習の解答ソースを個別で書きます（単独で実行できない）。

＊これからのソースは Tango ライブラリを利用します。コンパイルするには、Tango をインストールしなければ成りません。Tango については こちら を参照します。

練習1.21の解答

 

  1. /** 
  2.  * 練習1.21。 
  3.  * smallest-divisorで最少因子を求める。 
  4.  */  
  5. float _sqrt;  
  6.  
  7. int smallestDivisor(int n) {  
  8.     _sqrt = sqrt(cast(float) n);  
  9.     return findDivisor(n, 2);  
  10. }  
  11.   
  12. int findDivisor(int n, int testDivisor) {  
  13.     if (testDivisor > _sqrt) {  
  14.         return n;  
  15.     }  
  16.     if (n % testDivisor == 0) {  
  17.         return testDivisor;  
  18.     }  
  19.     testDivisor++;  
  20.     return findDivisor(n, testDivisor);  
  21. }  
  22.       
  23. void main() {  
  24.     int ret;  
  25.   
  26.     // for 1.21  
  27.     Stdout("---- 1.21 ----").newline;  
  28.     ret = smallestDivisor(199);  
  29.     Stdout.formatln("smallestDivisor(199) = {}", ret);  
  30.   
  31.     ret = smallestDivisor(1999);  
  32.     Stdout.formatln("smallestDivisor(1999) = {}", ret);  
  33.   
  34.     ret = smallestDivisor(19999);  
  35.     Stdout.formatln("smallestDivisor(19999) = {}", ret);  
  36. }  

実行結果。

---- 1.21 ----
smallestDivisor(199) = 199
smallestDivisor(1999) = 1999
smallestDivisor(19999) = 7


練習1.22の解答

 
  1. /** 
  2.  * 練習1.22。 
  3.  * smallestDivisor の時間コストを計る。 
  4.  */  
  5.   
  6. // 素数かどうか  
  7. bool prime(int n) {  
  8.     int sd = smallestDivisor(n);  
  9.     return sd == n;  
  10. }  
  11.   
  12. // 時間を計る  
  13. StopWatch _sw;  
  14. bool timedPrimeTest(int n) {  
  15.     _sw.start();  
  16.     if (prime(n)) {  
  17.         _sw.stop();  
  18.         ulong elapsed = _sw.microsec();  
  19.         Stdout.formatln("{} ***, elapsed: {}", n, elapsed);  
  20.         return true;  
  21.     }  
  22.     _sw.stop();  
  23.     Stdout.(n).newline;  
  24.     return false;  
  25. }  
  26.   
  27. void searchForPrimes(int start) {  
  28.     int s;  
  29.     if (start % 2 == 0) {  
  30.         s = start + 1;  
  31.     } else {  
  32.         s = start;  
  33.     }  
  34.     int max = s * 10;  
  35.     for (int i = s; i += 2; i < max) {  
  36.         if (timedPrimeTest(i)) {  
  37.             return;  
  38.         }  
  39.     }  
  40. }  
  41.   
  42. void main() {  
  43.     // for 1.22  
  44.     Stdout.newline;  
  45.     Stdout("---- 1.22 ----").newline;  
  46.     searchForPrimes(1000);  
  47.     searchForPrimes(10000);  
  48.     searchForPrimes(100000);  
  49.     searchForPrimes(1000000);  
  50. }  

実行結果。

---- 1.22 ----
1003
1005
1007
1009 ***, elapsed: 3
10003
10005
10007 ***, elapsed: 7
100003 ***, elapsed: 19
1000003 ***, elapsed: 53

smallest-divisor で素数を求める時の時間コストは O(平方根(n)) です。10 の平方根は 3.16 です。上記結果、O(10,000)/O(1,000) = 7/3 = 2.33、O(1,000,000)/O(100,000) = 53/19 = 2.79。n が大きい時、予想結果に近いことが分かります。


練習1.23の解答

 
  1. /** 
  2.  * 練習1.23。 
  3.  * 改良されたsmallest-divisor 
  4.  */  
  5. int smallestDivisor2(int n) {  
  6.     _sqrt = sqrt(cast(float) n);  
  7.     return findDivisor2(n, 2);  
  8. }  
  9.   
  10. int findDivisor2(int n, int testDivisor) {  
  11.     if (testDivisor > _sqrt) {  
  12.         return n;  
  13.     }  
  14.     if (n % testDivisor == 0) {  
  15.         return testDivisor;  
  16.     }  
  17.     if (testDivisor == 2) {  
  18.         testDivisor = 3;  
  19.     } else {  
  20.         testDivisor = testDivisor + 2;  
  21.     }  
  22.     return findDivisor2(n, testDivisor);  
  23. }  
  24.   
  25. bool prime2(int n) {  
  26.     int sd = smallestDivisor2(n);  
  27.     return sd == n;  
  28. }  
  29.   
  30. bool timedPrimeTest2(int n) {  
  31.     _sw.start();  
  32.     if (prime2(n)) {  
  33.         _sw.stop();  
  34.         ulong elapsed = _sw.microsec();  
  35.         Stdout.formatln("{} ***, elapsed: {}", n, elapsed);  
  36.         return true;  
  37.     }  
  38.     _sw.stop();  
  39.     Stdout(n).newline;  
  40.     return false;  
  41. }  
  42.   
  43. void searchForPrimes2(int start) {  
  44.     int s;  
  45.     if (start % 2 == 0) {  
  46.         s = start + 1;  
  47.     } else {  
  48.         s = start;  
  49.     }  
  50.     int max = s * 10;  
  51.     int count = 0;  
  52.     for (int i = s; i += 2; i < max) {  
  53.         if (timedPrimeTest2(i)) {  
  54.             // 三つまで探す  
  55.             if (++count == 3) {  
  56.                 return;  
  57.             }  
  58.         }  
  59.     }  
  60. }  
  61.   
  62. void main() {  
  63.     // for 1.23  
  64.     Stdout.newline;  
  65.     Stdout("---- 1.23 ----").newline;  
  66.     searchForPrimes2(1000);  
  67.     searchForPrimes2(10000);  
  68.     searchForPrimes2(100000);  
  69.     searchForPrimes2(1000000);  
  70. }  

実行結果。

---- 1.23 ----
1003
1005
1007
1009 ***, elapsed: 1
1011
1013 ***, elapsed: 2
1015
1017
1019 ***, elapsed: 2
10003
10005
10007 ***, elapsed: 3
10009 ***, elapsed: 3
10011
10013
10015
10017
10019
10021
10023
10025
10027
10029
10031
10033
10035
10037 ***, elapsed: 3
100003 ***, elapsed: 7
100005
100007
100009
100011
100013
100015
100017
100019 ***, elapsed: 7
100021
100023
100025
100027
100029
100031
100033
100035
100037
100039
100041
100043 ***, elapsed: 8
1000003 ***, elapsed: 30
1000005
1000007
1000009
1000011
1000013
1000015
1000017
1000019
1000021
1000023
1000025
1000027
1000029
1000031
1000033 ***, elapsed: 21
1000035
1000037 ***, elapsed: 22

改良される前（練習1.22）より、半分に近い早くなりました。他の判定を入れて、時間がかかるので、半分には行かないですね。


練習1.24の解答

 
  1. /** 
  2.  * 練習1.24 フェルマ検査。 
  3.  * Scheme 言語では計算できるらしいですが。 
  4.  * n が非常に大きい時（100000付近）、base の exp 乗はあまりにも大きくて、計算出来なくなる。 
  5.  */  
  6. int expmod(int base, int exp, int m) {  
  7.     if (exp == 0) {  
  8.         return 1;  
  9.     }  
  10.     if (exp % 2 == 0) {  
  11.         return square(expmod(base, exp / 2, m)) % m;  
  12.   
  13.     } else {  
  14.         return (base * expmod(base, exp - 1, m)) % m;  
  15.     }  
  16. }  
  17.   
  18. int square(int n) {  
  19.     return n * n;  
  20. }  
  21.   
  22. bool fermatTest(int n) {  
  23.     int a = Kiss.shared.toInt(1, n - 1) + 1;  
  24.     int ret = expmod(a, n, n);  
  25.     return ret == a;  
  26. }  
  27.   
  28. bool fastPrime(int n, int times) {  
  29.     if (times == 0) {  
  30.         return true;  
  31.     }  
  32.     if (fermatTest(n)) {  
  33.         return fastPrime(n, times - 1);  
  34.     } else {  
  35.         return false;  
  36.     }  
  37. }  
  38.   
  39. void searchForPrimes3(int start) {  
  40.     int s;  
  41.     if (start % 2 == 0) {  
  42.         s = start + 1;  
  43.     } else {  
  44.         s = start;  
  45.     }  
  46.     int max = s * 10;  
  47.     int count = 0;  
  48.     for (int i = s; i += 2; i < max) {  
  49.         _sw.start();  
  50.         if (fastPrime(i, 3)) {  
  51.             _sw.stop();  
  52.             ulong elapsed = _sw.microsec();  
  53.             Stdout.formatln("{} ***, elapsed: {}", i, elapsed);  
  54.             // 三つまで探す  
  55.             if (++count == 3) {  
  56.                 return;  
  57.             }  
  58.         }  
  59.         _sw.stop();  
  60.         Stdout(i).newline;  
  61.     }  
  62. }  
  63.   
  64. void main() {  
  65.     // for 1.24  
  66.     Stdout.newline;  
  67.     Stdout("---- 1.24 ----").newline;  
  68.     searchForPrimes3(10);  
  69.     searchForPrimes3(100);  
  70.     searchForPrimes3(1000);  
  71.     searchForPrimes3(10000);  
  72. }  

実行結果。

---- 1.24 ----
13 ***, elapsed: 5
13
15
17 ***, elapsed: 4
17
19 ***, elapsed: 2
103 ***, elapsed: 3
103
105
107 ***, elapsed: 3
107
109 ***, elapsed: 3
1003
1005
1007
1009 ***, elapsed: 4
1009
1011
1013 ***, elapsed: 4
1013
1015
1017
1019 ***, elapsed: 4
10003
10005
10007 ***, elapsed: 5
10007
10009 ***, elapsed: 5
10009
10011
10013
10015
10017
10019
10021
10023
10025
10027
10029
10031
10033
10035
10037 ***, elapsed: 4

コストの増加は大体 O(log(n)) と成っていることが分かります。


練習1.25の解答

fast-expt(base, exp) が大きすぎな数になるため、上手く行かないでしょう。


練習1.26の解答

expmod(base, exp/2, m) * expmod(base, exp/2, m) の所は、余計な計算をしているため。

an = an/2 * an/2 = an/4 * an/4 * an/4 * an/4 * = ... = a * a * ... * a (n 回) となってしまい、コストは O(n) です。


練習1.27の解答

 
  1. /** 
  2.  * 練習1.27。
  3.  * Carmichael 数。
  4.  */  
  5. void carmichael(int n) {  
  6.     bool fermat = true;  
  7.     int ret;  
  8.     for (int i = 1; i < n; i++) {  
  9.         ret = expmod(i, n, n);  
  10.         if (ret != i) {  
  11.             fermat = false;  
  12.             break;  
  13.         }  
  14.     }  
  15.     if (fermat) {  
  16.         Stdout.formatln("{} ***", n);  
  17.     } else {  
  18.         Stdout(n).newline;  
  19.     }  
  20. }  
  21.   
  22. void main() {  
  23.     // for 1.27  
  24.     Stdout.newline;  
  25.     Stdout("---- 1.27 ----").newline;  
  26.     carmichael(561);  
  27.     carmichael(1105);  
  28.     carmichael(1729);  
  29.     carmichael(2465);  
  30.     carmichael(2821);  
  31.     carmichael(6601);  
  32.     // 6789 は Carmichael でない  
  33.     carmichael(6789);  
  34. }  

実行結果。

---- 1.27 ----
561 ***
1105 ***
1729 ***
2465 ***
2821 ***
6601 ***
6789

完全に騙されましたね。6789は Carmichael 数でないため、フェルマ検査は騙されていない。

練習1.28の解答は割愛します。

閲覧  |  コメント  |  目次

 
ヘルプ  |  ご利用規約  |  相互リンク  |  問合せ
リンクはご自由に、問合せはお気軽に
©2007 Uprush