dkh.numeric.prime source code

1 module dkh.numeric.prime;
2 
3 /// xの約数一覧を返す
4 T[] divisorList(T)(T x) {
5     import std.algorithm : sort;
6     T[] res;
7     for (T i = 1; i*i <= x; i++) {
8         if (x%i == 0) {
9             res ~= i;
10             if (i*i != x) res ~= x/i;
11         }
12     }
13     sort(res);
14     return res;
15 }
16 
17 ///
18 unittest {
19     import std.range, std.algorithm;
20     assert(equal(divisorList(1), [1]));
21     assert(equal(divisorList(2), [1, 2]));
22     assert(equal(divisorList(4), [1, 2, 4]));
23     assert(equal(divisorList(24), [1, 2, 3, 4, 6, 8, 12, 24]));
24 }
25 
26 /// xの素因数一覧を返す
27 T[] factorList(T)(T x) {
28     T[] res;
29     for (T i = 2; i*i <= x; i++) {
30         while (x % i == 0) {
31             res ~= i;
32             x /= i;
33         }
34     }
35     if (x > 1) res ~= x;
36     // res is sorted!
37     return res;
38 }
39 
40 ///
41 unittest {
42     import std.range, std.algorithm;
43     assert(equal(factorList(1), new int[0]));
44     assert(equal(factorList(2), [2]));
45     assert(equal(factorList(4), [2, 2]));
46     assert(equal(factorList(24), [2, 2, 2, 3]));
47 }
48 
49 import dkh.numeric.primitive;
50 
51 /// Millar-Rabin Test
52 bool isPrime(ulong n) {
53     import dkh.int128;
54     if (n <= 1) return false;
55     if (n == 2) return true;
56     if (n % 2 == 0) return false;
57     ulong d = n-1;
58     while (d % 2 == 0) d /= 2;
59     ulong[] alist = [2,3,5,7,11,13,17,19,23,29,31,37];
60     foreach (a; alist) {
61         if (n <= a) break;
62         ulong y = ulongPowMod(a, d, n);
63         ulong t = d;
64         while (t != n-1 && y != 1 && y != n-1) {
65             y = mul128(y, y).mod128(n);
66             t <<= 1;
67         }
68         if (y != n-1 && t % 2 == 0) {
69             return false;
70         }
71     }
72     return true;
73 }
74 
75 ///
76 unittest {
77     assert(!isPrime(0));
78     assert(!isPrime(1));
79     assert(isPrime(2));
80     assert(isPrime(10^^9 + 7));
81 }