Math/NumberTheory/factorize_brent.h

• View this file on GitHub
• Last update: 2022-11-27 00:04:55+08:00
• Link: https://judge.yosupo.jp/problem/factorize
• Link: https://www.spoj.com/problems/FACT1/

Depends on

• Misc/int128.h

Code

// This seems to only work fine for random test cases
// - AC (20 digits) https://www.spoj.com/problems/FACT1/
// - TLE https://judge.yosupo.jp/problem/factorize

// brent {{{
#include "../../Misc/int128.h"

i128 brute(i128 n) {
    for (int i=2; i*i <= n; ++i) {
        if (n % i == 0) return i;
    }
    return n;
}

// returns random int in range [1, n-1]
i128 rand(i128 n) {
    return rng() % (n-1) + 1;
}

// return one divisor of n
i128 brent(i128 n) {
    if (n % 2 == 0) return 2;
    if (n < 200) return brute(n);

    i128 y = rand(n), c = rand(n), m = rand(n);
    auto sq_c = [c, n] (auto y) {
        return ((y*y) % n + c) % n;
    };

    i128 g = 1, r = 1, q = 1, x, ys;

    while (g == 1) {
        x = y;
        for (int i = 0; i < r; ++i) y = sq_c(y);

        int64_t k = 0;
        while (k < r && g == 1) {
            ys = y;
            for (int i = 0; i < min(m, r-k); ++i) {
                y = sq_c(y);
                q = q * my_abs(x - y) % n;
            }
            g = gcd(q, n);
            k += m;
        }
        r *= 2;
    }

    if (g == n) {
        while (true) {
            ys = sq_c(ys);
            g = gcd(my_abs(x - ys), n);
            if (g > 1) break;
        }
    }
    return g;
}

void factorize(i128 n, vector<i128>& fs) {
    if (n == 1) return;

    auto divisor = brent(n);
    if (divisor == n) {
        fs.push_back(divisor);
        return;
    }

    factorize(divisor, fs);
    factorize(n/divisor, fs);
}

vector<i128> factorize(i128 a) {
    if (a == 1) return {};

    vector<i128> res;
    factorize(a, res);
    sort(res.begin(), res.end());
    return res;
}
// }}}

#line 1 "Math/NumberTheory/factorize_brent.h"
// This seems to only work fine for random test cases
// - AC (20 digits) https://www.spoj.com/problems/FACT1/
// - TLE https://judge.yosupo.jp/problem/factorize

// brent {{{
#line 1 "Misc/int128.h"
// i128 helper functions {{{
using i128 = __int128_t;
i128 str2i128(std::string str) {
    i128 ret = 0;
    bool minus = false;
    for (auto c : str) {
        if (c == '-')
            minus = true;
        else
            ret = ret * 10 + c - '0';
    }
    return minus ? -ret : ret;
}
std::istream &operator>>(std::istream &is, i128 &x) {
    std::string s;
    return is >> s, x = str2i128(s), is;
}
std::ostream &operator<<(std::ostream &os, const i128 &x) {
    i128 tmp = x;
    if (tmp == 0) return os << 0;
    std::vector<int> ds;
    if (tmp < 0) {
        os << '-';
        while (tmp) {
            int d = tmp % 10;
            if (d > 0) d -= 10;
            ds.emplace_back(-d), tmp = (tmp - d) / 10;
        }
    } else {
        while (tmp) ds.emplace_back(tmp % 10), tmp /= 10;
    }
    std::reverse(ds.begin(), ds.end());
    for (auto i : ds) os << i;
    return os;
}
i128 my_abs(i128 n) {
    if (n < 0) return -n;
    return n;
}
i128 gcd(i128 a, i128 b) {
    if (b == 0) return a;
    return gcd(b, a % b);
}
// Count trailing zeroes
int ctz128(i128 n) {
    if (!n) return 128;
 
    if (!static_cast<uint64_t>(n)) {
        return __builtin_ctzll(static_cast<uint64_t>(n >> 64)) + 64;
    } else {
        return __builtin_ctzll(static_cast<uint64_t>(n));
    }
}
// }}}

#line 7 "Math/NumberTheory/factorize_brent.h"

i128 brute(i128 n) {
    for (int i=2; i*i <= n; ++i) {
        if (n % i == 0) return i;
    }
    return n;
}

// returns random int in range [1, n-1]
i128 rand(i128 n) {
    return rng() % (n-1) + 1;
}

// return one divisor of n
i128 brent(i128 n) {
    if (n % 2 == 0) return 2;
    if (n < 200) return brute(n);

    i128 y = rand(n), c = rand(n), m = rand(n);
    auto sq_c = [c, n] (auto y) {
        return ((y*y) % n + c) % n;
    };

    i128 g = 1, r = 1, q = 1, x, ys;

    while (g == 1) {
        x = y;
        for (int i = 0; i < r; ++i) y = sq_c(y);

        int64_t k = 0;
        while (k < r && g == 1) {
            ys = y;
            for (int i = 0; i < min(m, r-k); ++i) {
                y = sq_c(y);
                q = q * my_abs(x - y) % n;
            }
            g = gcd(q, n);
            k += m;
        }
        r *= 2;
    }

    if (g == n) {
        while (true) {
            ys = sq_c(ys);
            g = gcd(my_abs(x - ys), n);
            if (g > 1) break;
        }
    }
    return g;
}

void factorize(i128 n, vector<i128>& fs) {
    if (n == 1) return;

    auto divisor = brent(n);
    if (divisor == n) {
        fs.push_back(divisor);
        return;
    }

    factorize(divisor, fs);
    factorize(n/divisor, fs);
}

vector<i128> factorize(i128 a) {
    if (a == 1) return {};

    vector<i128> res;
    factorize(a, res);
    sort(res.begin(), res.end());
    return res;
}
// }}}

Back to top page