Factorisation (fixed code at http://e-maxx.ru/algo/factorization)

factorization.cpp

/*
Compiled against:

$ gcc --version
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 8.0.0 (clang-800.0.42.1)
Target: x86_64-apple-darwin16.7.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin
*/

#include <algorithm>
#include <cmath>
#include <iostream>
#include <map>
#include <vector>

using std::vector;
using std::map;

using std::cout;
using std::abs;

//! Возвращает true, если n четное
template<class T>
bool even(const T &n) {
    // return n % 2 == 0;
    return (n & 1) == 0;
}

//! Делит число на 2
template<class T>
void bisect(T &n) {
    // n /= 2;
    n >>= 1;
}

//! Умножает число на 2
template<class T>
void redouble(T &n) {
    // n *= 2;
    n <<= 1;
}

//! Возвращает true, если n - точный квадрат простого числа
template<class T>
bool perfect_square(const T &n) {
    T sq = (T) ceil(sqrt((double) n));
    return sq * sq == n;
}

//! Вычисляет корень из числа, округляя его вниз
template<class T>
T sq_root(const T &n) {
    return (T) floor(sqrt((double) n));
}

//! Возвращает количество бит в числе
template<class T>
unsigned bits_in_number(T n) {
    if (n == 0)
        return 1;
    unsigned result = 0;
    while (n) {
        bisect(n);
        ++result;
    }
    return result;
}

//! Возвращает значение k-го бита числа (биты нумеруются с нуля)
template<class T>
bool test_bit(const T &n, unsigned k) {
    return (n & (T(1) << k)) != 0;
}

//! Умножает a *= b (mod n)
template<class T>
void mulmod(T &a, T b, const T &n) {
    a *= b;
    a %= n;
}

template<>
void mulmod(int &a, int b, const int &n) {
    a = int((((long long) a) * b) % n);
}

template<>
void mulmod(unsigned &a, unsigned b, const unsigned &n) {
    a = unsigned((((unsigned long long) a) * b) % n);
}

template<>
void mulmod(unsigned long long &a, unsigned long long b, const unsigned long long &n) {
    if (a >= n)
        a %= n;
    if (b >= n)
        b %= n;
    unsigned long long res = 0;
    while (b)
        if (!even(b)) {
            res += a;
            while (res >= n)
                res -= n;
            --b;
        } else {
            redouble(a);
            while (a >= n)
                a -= n;
            bisect(b);
        }
    a = res;
}

template<>
void mulmod(long long &a, long long b, const long long &n) {
    bool neg = false;
    if (a < 0) {
        neg = !neg;
        a = -a;
    }
    if (b < 0) {
        neg = !neg;
        b = -b;
    }
    unsigned long long aa = a;
    mulmod<unsigned long long>(aa, (unsigned long long) b, (unsigned long long) n);
    a = (long long) aa * (neg ? -1 : 1);
}


//! Вычисляет a^k (mod n)
template<class T, class T2>
T powmod(T a, T2 k, const T &n) {
    T res = 1;
    while (k)
        if (!even(k)) {
            mulmod(res, a, n);
            --k;
        } else {
            mulmod(a, a, n);
            bisect(k);
        }
    return res;
}

//! Переводит число n в форму q*2^p
template<class T>
void transform_num(T n, T &p, T &q) {
    T p_res = 0;
    while (even(n)) {
        ++p_res;
        bisect(n);
    }
    p = p_res;
    q = n;
}

//! Алгоритм Евклида
template<class T, class T2>
T gcd(const T &a, const T2 &b) {
    return (a == 0) ? b : gcd(b % a, a);
}

//! Вычисляет jacobi(a,b) - символ Якоби
template<class T>
T jacobi(T a, T b) {

#pragma warning (push)
#pragma warning (disable: 4146)

    if (a == 0)
        return 0;
    if (a == 1)
        return 1;

    if (a < 0) {
        if ((b & 2) == 0) {
            return jacobi(-a, b);
        } else {
            return -jacobi(-a, b);
        }
    }

    T e, a1;
    transform_num(a, e, a1);

    T s;
    if (even(e) || (b & 7) == 1 || (b & 7) == 7)
        s = 1;
    else
        s = -1;
    if ((b & 3) == 3 && (a1 & 3) == 3)
        s = -s;
    if (a1 == 1)
        return s;
    return s * jacobi(b % a1, a1);

#pragma warning (pop)

}

//! Вычисляет pi(b) первых простых чисел. Возвращает вектор с простыми и в pi - pi(b)
template<class T, class T2>
const std::vector<T> &get_primes(const T &b, T2 &pi) {

    static std::vector<T> primes;
    static T counted_b;

    if (counted_b >= b)
        pi = T2(std::upper_bound(primes.begin(), primes.end(), b) - primes.begin());
    else {

        if (counted_b == 0) {
            primes.push_back(2);
            counted_b = 2;
        }

        T first = counted_b == 2 ? 3 : primes.back() + 2;
        for (T cur = first; cur <= b; ++ ++cur) {
            bool cur_is_prime = true;
            for (auto iter = primes.begin(), end = primes.end(); iter != end; ++iter) {
                const T &div = *iter;
                if (div * div > cur)
                    break;
                if (cur % div == 0) {
                    cur_is_prime = false;
                    break;
                }
            }
            if (cur_is_prime)
                primes.push_back(cur);
        }

        counted_b = b;
        pi = (T2) primes.size();

    }

    return primes;

}

//! Тривиальная проверка n на простоту, перебираются все делители до m.
//! Результат: 1 - если n точно простое, p - его найденный делитель, 0 - если неизвестно
template<class T, class T2>
T2 prime_div_trivial(const T &n, T2 m) {

    if (n == 2 || n == 3)
        return 1;
    if (n < 2)
        return 0;
    if (even(n))
        return 2;

    T2 pi;
    const vector<T2> &primes = get_primes(m, pi);

    for (auto iter = primes.begin(), end = primes.end(); iter != end; ++iter) {
        const T2 &div = *iter;
        if (div * div > n)
            break;
        else if (n % div == 0)
            return div;
    }

    if (n < m * m)
        return 1;
    return 0;

}

template<class T, class T2>
bool miller_rabin(T n, T2 b) {
    // сначала проверяем тривиальные случаи
    if (n == 2)
        return true;
    if (n < 2 || even(n))
        return false;

    // проверяем, что n и b взаимно просты (иначе это приведет к ошибке)
    // если они не взаимно просты, то либо n не просто, либо нужно увеличить b
    if (b < 2)
        b = 2;
    for (T g; (g = gcd(n, b)) != 1; ++b)
        if (n > g)
            return false;

    // разлагаем n-1 = q*2^p
    T n_1 = n;
    --n_1;
    T p, q;
    transform_num(n_1, p, q);

    // вычисляем b^q mod n, если оно равно 1 или n-1, то n простое (или псевдопростое)
    T rem = powmod(T(b), q, n);
    if (rem == 1 || rem == n_1)
        return true;

    // теперь вычисляем b^2q, b^4q, ... , b^((n-1)/2)
    // если какое-либо из них равно n-1, то n простое (или псевдопростое)
    for (T i = 1; i < p; i++) {
        mulmod(rem, rem, n);
        if (rem == n_1)
            return true;
    }

    return false;
}

template<class T, class T2>
bool lucas_selfridge(const T &n, T2 unused) {
    // сначала проверяем тривиальные случаи
    if (n == 2)
        return true;
    if (n < 2 || even(n))
        return false;

    // проверяем, что n не является точным квадратом, иначе алгоритм даст ошибку
    if (perfect_square(n))
        return false;

    // алгоритм Селфриджа: находим первое число d такое, что:
    // jacobi(d,n)=-1 и оно принадлежит ряду { 5,-7,9,-11,13,... }
    T2 dd;
    for (T2 d_abs = 5, d_sign = 1;; d_sign = -d_sign, ++ ++d_abs) {
        dd = d_abs * d_sign;
        T g = gcd(n, d_abs);
        if (1 < g && g < n)
            // нашли делитель - d_abs
            return false;
        if (jacobi(T(dd), n) == -1)
            break;
    }

    // параметры Селфриджа
    T2
            p = 1,
            q = (p * p - dd) / 4;

    // разлагаем n+1 = d*2^s
    T n_1 = n;
    ++n_1;
    T s, d;
    transform_num(n_1, s, d);

    // алгоритм Лукаса
    T
            u = 1,
            v = p,
            u2m = 1,
            v2m = p,
            qm = q,
            qm2 = q * 2,
            qkd = q;
    for (unsigned bit = 1, bits = bits_in_number(d); bit < bits; bit++) {
        mulmod(u2m, v2m, n);
        mulmod(v2m, v2m, n);
        while (v2m < qm2)
            v2m += n;
        v2m -= qm2;
        mulmod(qm, qm, n);
        qm2 = qm;
        redouble(qm2);
        if (test_bit(d, bit)) {
            T t1, t2;
            t1 = u2m;
            mulmod(t1, v, n);
            t2 = v2m;
            mulmod(t2, u, n);

            T t3, t4;
            t3 = v2m;
            mulmod(t3, v, n);
            t4 = u2m;
            mulmod(t4, u, n);
            mulmod(t4, (T) dd, n);

            u = t1 + t2;
            if (!even(u))
                u += n;
            bisect(u);
            u %= n;

            v = t3 + t4;
            if (!even(v))
                v += n;
            bisect(v);
            v %= n;
            mulmod(qkd, qm, n);
        }
    }

    // точно простое (или псевдо-простое)
    if (u == 0 || v == 0)
        return true;

    // довычисляем оставшиеся члены
    T qkd2 = qkd;
    redouble(qkd2);
    for (T2 r = 1; r < s; ++r) {
        mulmod(v, v, n);
        v -= qkd2;
        if (v < 0) v += n;
        if (v < 0) v += n;
        if (v >= n) v -= n;
        if (v >= n) v -= n;
        if (v == 0)
            return true;
        if (r < s - 1) {
            mulmod(qkd, qkd, n);
            qkd2 = qkd;
            redouble(qkd2);
        }
    }

    return false;
}

template<class T>
bool baillie_pomerance_selfridge_wagstaff(T n) {
    // сначала проверяем на тривиальные делители - например, до 29
    int div = prime_div_trivial(n, 29);
    if (div == 1)
        return true;
    if (div > 1)
        return false;

    // тест Миллера-Рабина по основанию 2
    if (!miller_rabin(n, 2))
        return false;

    // сильный тест Лукаса-Селфриджа
    return lucas_selfridge(n, 0);
}

template<class T>
bool isprime(T n) {
    return baillie_pomerance_selfridge_wagstaff(n);
}

template<class T>
T pollard_p_1(T n) {
    // параметры алгоритма, существенно влияют на производительность и качество поиска
    const T b = 13;
    const T q[] = {2, 3, 5, 7, 11, 13};

    // несколько попыток алгоритма
    T a = 5 % n;
    for (int j = 0; j < 10; j++) {

        // ищем такое a, которое взаимно просто с n
        while (gcd(a, n) != 1) {
            mulmod(a, a, n);
            a += 3;
            a %= n;
        }

        // вычисляем a^M
        for (size_t i = 0; i < sizeof q / sizeof q[0]; i++) {
            T qq = q[i];
            T e = (T) floor(log((double) b) / log((double) qq));
            T aa = powmod(a, powmod(qq, e, n), n);
            if (aa == 0)
                continue;

            // проверяем, не найден ли ответ
            T g = gcd(aa - 1, n);
            if (1 < g && g < n)
                return g;
        }

    }

    // если ничего не нашли
    return 1;
}

template<class T>
T pollard_rho(T n, unsigned iterations_count = 100000) {
    T
            b0 = rand() % n,
            b1 = b0,
            g;
    mulmod(b1, b1, n);
    if (++b1 == n)
        b1 = 0;
    g = gcd(abs(b1 - b0), n);
    for (unsigned count = 0; count < iterations_count && (g == 1 || g == n); count++) {
        mulmod(b0, b0, n);
        if (++b0 == n)
            b0 = 0;
        mulmod(b1, b1, n);
        ++b1;
        mulmod(b1, b1, n);
        if (++b1 == n)
            b1 = 0;
        g = gcd(abs(b1 - b0), n);
    }
    return g;
}

template<class T>
T pollard_bent(T n, unsigned iterations_count = 19) {
    T
            b0 = rand() % n,
            b1 = (b0 * b0 + 2) % n,
            a = b1;
    for (unsigned iteration = 0, series_len = 1;
         iteration < iterations_count; iteration++, series_len *= 2) {
        T g = gcd(b1 - b0, n);
        for (unsigned len = 0; len < series_len && (g == 1 && g == n); len++) {
            b1 = (b1 * b1 + 2) % n;
            g = gcd(abs(b1 - b0), n);
        }
        b0 = a;
        a = b1;
        if (g != 1 && g != n)
            return g;
    }
    return 1;
}

template<class T>
T pollard_monte_carlo(T n, unsigned m = 100) {
    T b = rand() % (m - 2) + 2;

    static std::vector<T> primes;
    static T m_max;
    if (primes.empty())
        primes.push_back(3);
    if (m_max < m) {
        m_max = m;
        for (T prime = 5; prime <= m; ++ ++prime) {
            bool is_prime = true;
            for (auto iter = primes.begin(), end = primes.end(); iter != end; ++iter) {
                T div = *iter;
                if (div * div > prime)
                    break;
                if (prime % div == 0) {
                    is_prime = false;
                    break;
                }
            }
            if (is_prime)
                primes.push_back(prime);
        }
    }

    T g = 1;
    for (size_t i = 0; i < primes.size() && g == 1; i++) {
        T cur = primes[i];
        while (cur <= n)
            cur *= primes[i];
        cur /= primes[i];
        b = powmod(b, cur, n);
        g = gcd(abs(b - 1), n);
        if (g == n)
            g = 1;
    }

    return g;
}

template<class T, class T2>
T ferma(const T &n, T2 unused) {
    T2
            x = sq_root(n),
            y = 0,
            r = x * x - y * y - n;
    for (;;)
        if (r == 0)
            return x != y ? x - y : x + y;
        else if (r > 0) {
            r -= y + y + 1;
            ++y;
        } else {
            r += x + x + 1;
            ++x;
        }
}

template<class T, class T2>
void factorize(const T &n, std::map<T, unsigned> &result, T2 unused) {
    if (n == 1);
    else
        // проверяем, не простое ли число
    if (isprime(n))
        ++result[n];
    else
        // если число достаточно маленькое, то его разлагаем простым перебором
    if (n < 1000 * 1000) {
        T div = prime_div_trivial(n, 1000);
        ++result[div];
        factorize(n / div, result, unused);
    } else {
        // число большое, запускаем на нем алгоритмы факторизации
        T div;
        // сначала идут быстрые алгоритмы Полларда
        div = pollard_monte_carlo(n);
        if (div == 1)
            div = pollard_rho(n);
        if (div == 1)
            div = pollard_p_1(n);
        if (div == 1)
            div = pollard_bent(n);
        // придётся запускать 100%-ый алгоритм Ферма
        if (div == 1)
            div = ferma(n, unused);
        // рекурсивно обрабатываем найденные множители
        factorize(div, result, unused);
        factorize(n / div, result, unused);
    }
}

int main() {
    typedef long long base;

    base n = 1024 * 3 * 3 * 3 * 7 * 7;

    map<base, unsigned> m;
    factorize(n, m, (long long) 0);

    for (auto &i : m) {
        cout << i.first << "^" << i.second << " * ";
    }
    cout << '1';

    return 0;
}