dkh.numeric.prime source code

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