p-1 Pollard Factorization algorithm

pollard_base.cpp

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <vector>

using namespace std;

typedef unsigned long long ull;

ull gcd (ull a, ull b) {
    return b ? gcd (b, a % b) : a;
}

vector<int> generateAllPrimesBelow(int b) {
    int n = 7, h = 5, s = 25;

    int firstPrimes[] = {2, 3, 5, 7, 11, 13, 17, 17};
    vector<int> primes (firstPrimes, firstPrimes + sizeof(firstPrimes) / sizeof(int));
    primes.resize(100);

    step3:
    {
        primes[n] += 2;
        int k = 1;

        if (primes[n] > b) {
            primes.resize(n);
            return primes;
        }

        if (primes[n] > s) {
            s = s + h;
            h = h + 1;
            s = s + h;
        }

        while (primes[k] <= h) {
            if (primes[n] % primes[k] == 0) {
                goto step3;
            }
            k += 1;
        }

        n += 1;
        primes[n] = primes[n-1];
        goto step3;
    }

    return primes;
}

ull m_pow(unsigned int x, int a) {
    ull result = x;
    for (int i = 1; i < a; i++) {
        result *= x;
    }
    return result;
}

/**
 * Tries to find a prime factor
 * @return p where p | m and p is prime
 */
long long factorizeNumber(ull m, ull B = 10^6, int numberOfIterations = 10) {
    srand(time(NULL));

    vector<int> primes = generateAllPrimesBelow(B);

    ull k = 1;
    for (vector<int>::iterator it = primes.begin() ; it != primes.end(); ++it) {
        int alpha = (*it <= 31) ? rand() % 10 : 1;
        k *=  *it ^ alpha;
    }

    ull p;
    int a;
    int step = -1;
    while(1) {
        step += 1;
        a = primes[step];

        if (gcd(a, m) > 1) { return a; }

        p = gcd(m, (m_pow(a, k) - 1) % m);
        if ((m > p) & (p > 1)) {
            return p;
        }
        if (step > numberOfIterations) {
            return -1;
        }
    }
}

int main() {
    int m, B, c;

    cout << "Целое составное число m: ";
    cin >> m;
    cout << "Граница B: ";
    cin >> B;

    long long factor = factorizeNumber(m, B);
    if (factor == -1) {
        cout << "Failed to find factor!" << endl;
    } else {
        cout << factor << endl;
    }

    return 0;
}