dkh.numeric.primitive source code

1 module dkh.numeric.primitive;
2 
3 import std.traits, std.bigint;
4 
5 /// lcm
6 T lcm(T)(in T a, in T b) {
7     import std.numeric : gcd;
8     return a / gcd(a,b) * b;
9 }
10 
11 ///
12 unittest {
13     assert(lcm(2, 4) == 4);
14     assert(lcm(3, 5) == 15);
15     assert(lcm(1, 1) == 1);
16     assert(lcm(0, 100) == 0);
17 }
18 
19 /// 高速累乗
20 Unqual!T pow(T, U)(T x, U n)
21 if (!isFloatingPoint!T && (isIntegral!U || is(U == BigInt))) {
22     return pow(x, n, T(1));
23 }
24 
25 /// ditto
26 Unqual!T pow(T, U, V)(T x, U n, V e)
27 if ((isIntegral!U || is(U == BigInt)) && is(Unqual!T == Unqual!V)) {
28     Unqual!T b = x, v = e;
29     Unqual!U m = n;
30     while (m) {
31         if (m & 1) v *= b;
32         b *= b;
33         m /= 2;
34     }
35     return v;
36 }
37 
38 unittest {
39     assert(pow(3, 5) == 243);
40     assert(pow(3, 5, 2) == 486);
41 }
42 
43 ///
44 T powMod(T, U, V)(T x, U n, V md)
45 if (isIntegral!U || is(U == BigInt)) {
46     T r = T(1);
47     Unqual!U m = n;
48     while (m) {
49         if (m & 1) r = (r*x)%md;
50         x = (x*x)%md;
51         m >>= 1;
52     }
53     return r % md;
54 }
55 
56 unittest {
57     immutable int B = 3;
58     assert(powMod(5, B, 100) == 25); //125 % 100
59 }
60 
61 //todo: consider binary extgcd
62 /// a*T[0]+b*T[1]=T[2], T[2]=gcd
63 T[3] extGcd(T)(in T a, in T b) 
64 if (!isIntegral!T || isSigned!T) //unsignedはNG
65 {
66     if (b==0) {
67         return [T(1), T(0), a];
68     } else {
69         auto e = extGcd(b, a%b);
70         return [e[1], e[0]-a/b*e[1], e[2]];
71     }
72 }
73 
74 ///
75 unittest {
76     import std.numeric : gcd;
77     foreach (i; 0..100) {
78         foreach (j; 0..100) {
79             auto e = extGcd(i, j);
80             assert(e[2] == gcd(i, j));
81             assert(e[0] * i + e[1] * j == e[2]);
82         }
83     }
84 }
85 
86 /// calc inverse, (x * invMod(x)) % md == 1
87 T invMod(T)(T x, T md) {
88     auto r = extGcd!T(x, md);
89     assert(r[2] == 1);
90     auto z = r[0];
91     return (z % md + md) % md;
92 }
93 
94 unittest {
95     import std.numeric : gcd;
96     foreach (i; 1..100) {
97         foreach (j; 1..i) {
98             if (gcd(i, j) != 1) continue;
99             auto k = invMod(j, i);
100             assert(j * k % i == 1);
101         }
102     }
103 }