Added carmichael function

Author: fredbr

math/binary_gcd.cpp

@@ -0,0 +1,29 @@
+#include <bits/stdc++.h>
+
+// important!
+#pragma GCC target("bmi")
+
+template<typename T>
+T gcd(T a, T b) {
+    auto ctz = [] (T const& x) { return __builtin_ctzll(x); };
+
+    int shift = ctz(a|b);
+    b >>= ctz(b);
+
+    while (a) {
+        a >>= ctz(a);
+        if (a < b) std::swap(a, b);
+        a -= b;
+    }
+
+    return b << shift;
+}
+
+int main() {
+    using ll = uint64_t;
+    
+    ll a, b;
+    std::cin >> a >> b;
+
+    std::cout << gcd<ll>(a, b) << "\n";
+}

math/carmichael.cpp

@@ -0,0 +1,281 @@
+#include <bits/stdc++.h>
+
+#include <bits/stdc++.h>
+
+using namespace std;
+
+using u128 = __uint128_t;
+using s128 = __int128_t;
+using u64 = uint64_t;
+
+inline u128 mul(u128 x, u128 y, u128 c)
+{
+    u128 ans = 0;
+    static const u64 shift = 25;
+    while (y) {
+        ans += x*(y&((1<<shift)-1));
+        x = (x<<shift)%c;
+        y >>= shift;
+    }
+    return ans % c;
+}
+
+u128 HI(u128 x) { return x>>64; };    
+u128 LO(u128 x) { return u64(x); };
+    
+namespace Prime {
+
+    int const num_tries = 20;
+
+    u128 power(u128 a, u128 b, u128 c) {
+        u128 ans = 1;
+        while (b) {
+
+            if (b&1) ans = mul(ans, a, c);
+            a = mul(a, a, c);
+            b >>= 1;
+        }
+        return ans;
+    }
+
+    bool witness(u128 base, u128 s, u128 d, u128 n) {
+        u128 x = power(base, d, n);
+        u128 y;
+        while (s--) {
+
+            y = mul(x, x, n);
+            if (y == 1 and x != 1 and x != n-1) return false;
+            x = y;
+        }
+        if (x != 1) return false;
+        return true;
+    }
+
+    mt19937_64 rd;
+
+    bool is_prime(u128 x) {
+        if (x == 2 or x == 3 or x == 1) return true;
+        if (x % 2 == 0 or x % 3 == 0) return false;
+
+        u128 s = 1, d = x>>1;
+        while (!(d&1)) s++, d>>= 1;
+
+        uniform_int_distribution<u64> rng(2, HI(x)? u64(-1ull): (u64(x)-1));
+
+        for (int i = 0; i < num_tries; i++) {
+            u64 p = rng(rd);
+            if (!witness(p, s, d, x)) return false;
+        }
+        return true;
+    }
+} //end namespace Prime
+
+namespace Factor {
+
+    int const num_tries = 10;
+
+    u128 n, niv_n;
+    void setn(u128 n_) {
+        // calculates multiplicative inverse of n mod 2^128
+        // to use in modular reduction (how?)
+        n = n_;
+        niv_n = 1;
+        u128 x = n;
+        for (int i = 1; i <= 126; i++)
+        {
+            niv_n*=x;
+            x*=x; 
+        } 
+    }
+  
+    struct u256 {
+        u128 hi, lo;
+        
+        static u256 mul128(u128 x,u128 y) {
+            u128 t1 = LO(x)*LO(y);
+            u128 t2 = HI(x)*LO(y) + HI(y)*LO(x) + HI(t1);
+            return { HI(x)*HI(y) + HI(t2) , (t2<<64) + LO(t1) };
+        }
+    };
+
+    u128 mmul(u128 x,u128 y) {   
+        u256 m = u256::mul128(x,y);
+        // some kind of montgomery reduction?
+        u128 ans = m.hi - u256::mul128(m.lo*niv_n, n).hi;
+        if (s128(ans) < 0) ans += n;
+        return ans;
+    }
+
+    inline u128 f(u128 x, u128 inc) {
+        return mmul(x, x) + inc;
+    }
+
+    inline u128 gcd(u128 a, u128 b) {
+        auto ctz = [] (u128 x) {
+            if (!HI(x)) return __builtin_ctzll(x);
+            return 64 + __builtin_ctzll(x>>64);
+        };
+
+        int shift = ctz(a|b);
+        b >>= ctz(b);
+        while (a) {
+            a >>= ctz(a);
+            if (a < b) swap(a, b);
+            a -= b;
+        }
+        return b << shift;
+    }
+
+    inline u128 rho(u128 seed, u128 n, u64 inc) {
+        static const int step=1<<9;
+        
+        setn(n);
+
+        auto sub = [&] (u128 x, u128 y) { return x>y ? x-y:y-x; };
+
+        u128 y = f(seed, inc);
+
+        for (int l=1;; l<<=1)
+        {
+            u128 x = y;
+            for (int i = 1; i <= l; i++) y = f(y, inc);
+            for (int k = 0; k < l; k += step)
+            {
+                int d = min(step, l-k);
+                u128 g = seed, y0 = y;
+
+                for (int i = 1; i <= d; i++) {
+                    y = f(y, inc);
+                    g = mmul(g, sub(x,y));
+                }
+                
+                g = gcd(g,n);
+                
+                if (g==1) continue;
+                if (g!=n) return g;
+                
+                u128 y = y0;
+                
+                while (gcd(sub(x,y), n) == 1) y = f(y, inc);
+                
+                return gcd(sub(x,y), n) % n;
+            }
+        }
+        return 0;
+    }
+
+    mt19937_64 rd;
+
+    u128 rho(u128 x) {  
+        if (x <= 3) return 0;
+        if (x%2 == 0) return 2;
+        if (x%3 == 0) return 3;
+
+        uniform_int_distribution<u64> rng(2, HI(x)? u64(-1ull): (u64(x)-1));
+
+        for (u128 i = 2; i < num_tries; i++) {
+            u128 ans = rho(rng(rd), x, i);
+            if (ans != 0 and ans != x) return ans;
+        }
+        return 0;
+    }
+
+    void factor(u128 x, vector<u128>& primes) {
+
+        if (Prime::is_prime(x)) primes.push_back(x);
+        else {
+            u128 fac = rho(x);
+            factor(fac, primes);
+            factor(x/fac, primes);
+        }
+    }
+
+    vector<u128> factor(u128 x) {
+        vector<u128> v;
+        factor(x, v);
+        sort(v.begin(), v.end());
+        return v;
+    }
+
+    vector<pair<u128,u128>> factor_pairs(u128 x) {
+        auto v = factor(x);
+        vector<pair<u128, u128>> res;
+        
+        for (int i = 0; i < (int)v.size(); i++) {
+            if (i == 0 or res.back().first != v[i])
+                res.emplace_back(v[i], 1);
+            else
+                res.back().second++;
+        }
+        return res;
+    }
+} // end namespace Factor
+
+template <typename T> 
+void read(T &x)
+{
+    #define gc (c=getchar())
+
+    char c;
+    while(gc < '0');
+    x = c - '0';
+    while(gc >= '0') x = x*10 + c - '0';
+
+    #undef gc
+}
+
+template <typename T> 
+void write(T x, char c = '\n')
+{
+    static char st[40];
+    int top=0;
+    do {st[++top] = '0' + x%10;} while (x /=10);
+    do {putchar(st[top]);} while (--top);
+    if (c != 0) putchar(c);
+}
+
+template<typename T>
+T gcd(T a, T b) {
+    if (b == 0) return a;
+    return gcd(b, a%b);
+}
+
+template<typename T>
+T lcm(T a, T b) {
+    return a/gcd(a, b)*b;
+}
+
+template<typename T>
+T power(T a, T b) {
+    T ans = 1;
+    while (b) {
+        if (b&1) ans *= a;
+        a *= a;
+        b >>= 1;
+    }
+    return ans;
+}
+
+template<typename T>
+T carmichael(pair<T, T> const& x) {
+    T ans = power(x.first, x.second-1) * (x.first-1);
+    return ans;
+}
+template<typename T>
+T carmichael(T x) {
+    auto factors = Factor::factor_pairs(x);
+    T ans = 1;
+
+    for (auto const& i : factors) {
+        ans = lcm(ans, carmichael(i));
+    }
+
+    return ans;
+}
+
+int main() {
+    u128 x;
+    read(x);
+
+    write(carmichael(x));
+}

math/determinant.cpp

@@ -47,6 +47,10 @@ T determinant(std::vector<std::vector<T>> A) {
             std::swap(A[best], A[i]);
         }
 
+        if (A[i][i] == T()) {
+            return T();
+        }
+
         ans *= T(A[i][i]);
         divide_row<T>(A[i], A[i][i]);
 

math/pollard_rho.cpp

@@ -120,8 +120,6 @@ int main()
 	cin >> x;
 
 	vector<ll> fac;
-	for (int i = 0; i < 100; i++)
-		Factor::factor(fac, x), fac.clear();
 
 	Factor::factor(fac, x);
 	for (ll i : fac) cout << i << " ";