Реализации алгоритмов/P-1 метод Полларда

P-1 метод Полларда (читается как п-1 метод Полларда) — один из методов факторизации целых чисел.

С/С++

Ниже приведен заготовка функции на С++, реализующая первую и вторую стадии алгоритма. Некоторые пояснения:

Класс Т — класс, с которым производятся вычисления, это может быть как стандартный числовой класс, например long int, так и пользовательский класс, для которого необходимо перегрузить функции и операции '+', '-', '+=', '-=', взятия остатка по модулю n '%', '%='.
Массив q — массив простых чисел меньше либо равных b (B в предыдущих обозначениях).
Массив p — массив простых чисел между границами $~B1$ и $~B2=B1^{2}$ .
Необходимо определить так же функции mulmod(T a,T b,T n)— перемножение чисел a,b по модулю n &nbsp $~mulmod(a,b,n)=ab\mod n$ , gsd(T a,T b) — поиск НОД чисел a и b, powmod(T a,T b, T n) — возведение числа a в степень b по модулю n $~powmod(a,b,n)=a^{b}\mod n$ .
Первая стадия выполняется для нескольких а.

template <class T>
T pollard_p_1_1 (T n)
{
    // параметры алгоритма, существенно влияют на производительность и качество поиска
    const T b1 = 13;
    const T q[] = { 2, 3, 5, 7, 11, 13 };
    const T b1=13;
    const T b2=169; //b2=b1^2
    T p[]= [13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
     101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167];
    T *D;
    size_t i;
    // несколько попыток алгоритма
    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 (i = 0; i < sizeof q / sizeof q[0]; i++)
        {
            T qq = q[i];
            T e = (T) floor (log ((double)b1) / 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;
        }

    }
    T apowM=aa;//apowM=a^M посчитанное в первой стадии
    for (i = 0; i< (sizeof q / sizeof p[0])-1; i++)
    {
        p[i]= p[i+1]-p[i]; //переопределяем массив p
    }
    for ( i = 0; i< (sizeof q / sizeof p[0])-1; i++)
    {
        aa = mulmod(aa, powmod(aa, p[i], n), n);
        g = gsd (aa-1, n);
        if (g !=1 ) 
            return g;
    }
    
    // если ничего не нашли
     return 1;

}