pollard rho algorithm in c++ never finds factors

Question

I am trying to implement the pollard rho algorithm to factor large numbers. I am using the gmp package. The algorithm is taken from "Introduction to algorithms". The symptoms is that the while loop never breaks. My thoughts is that I am missing something in the implementation.

My code is as follows:

---

#include <iostream>
#include <gmp.h>
#include <stdio.h>
#include <math.h>
#include <gmpxx.h>

using namespace std;
//int numberToFactor;
int compare1;
int compare2;
int factor[100] ;
int k;
int i;



int main (){
    mpz_t numberToFactor;
    mpz_init(numberToFactor);
    mpz_t D;
    mpz_init(D);
    mpz_t x;
    mpz_init(x);
    mpz_t y;
    mpz_init(y);
    mpz_t random;
    mpz_init(random);
    mpz_t largest;
    mpz_init(largest);
    mpz_t square;
    mpz_init(square);
    mpz_t yMinusx;
    mpz_init(yMinusx);

    unsigned long int exp = 2;
    unsigned long int one = 1;


    gmp_randstate_t state;
    gmp_randinit_mt(state);


    while(stdin){
            cout<<"enter number";
        cin >> numberToFactor;

        //Is it a prime number?
        int c = mpz_probab_prime_p(numberToFactor, 5);

//Start of Pollard Rho algorithm
        i = 1;
        //Generate a random number random.
        mpz_urandomm(x, state, numberToFactor);

        //y = x;
        mpz_set(y, x);
        cout<<"y="<<y<<"\n";
        k = 2;



        while(true){
            i++;
            //square = x².
            mpz_pow_ui(square, x, exp);
            cout<<"x="<<x<<"\n";
            cout<<"x²="<<square<<"\n";

            //square = x²-1
            mpz_sub_ui(square, square, one);
            cout<<"x²-1="<<square<<"\n";

            //x= (x²-1)mod numberToFactor
            mpz_mod(x, square, numberToFactor); 
            cout<<"x = x²-1 mod numberToFactor="<<x<<"\n";

            //y-x
            mpz_sub(yMinusx, y, x);
            //abs(y-x)
            mpz_abs(yMinusx, yMinusx);
            cout<<"y-x="<<yMinusx<<"\n";

            //D=gcd(abs(y-x),numberToFactor).
            mpz_gcd(D, yMinusx, numberToFactor);
            cout<<"gcd(abs(y-x), numberToFactor)="<<D<<"\n";

            compare1 = mpz_cmp_ui(D, one);
            compare2 = mpz_cmp(D, numberToFactor);
            //is D!=1 AND D!=numberToFactor?
            if(compare1 != 0 && compare2 != 0){
                cout<<"D"<<D<< "\n";
                //is D a prime?
                if(mpz_probab_prime_p(D, 5) == 2){
                    //save D?
                    cout<<"Found one prime divisor!"<<"\n";
                    cout<<"prime = "<<D<<"\n";
                }
                //is the remainder a prime?
                //numberToFactor = numberTofactor/D
                mpz_div(numberToFactor , numberToFactor, D);
                if(mpz_probab_prime_p(numberToFactor, 5) == 2){

                    cout<<"The remainder is a prime!"<<"\n";
                    cout<<"prime = "<<numberToFactor<<"\n";
                }
                compare1 = mpz_cmp_ui(numberToFactor, one);

                if( 0 == compare1 ){
                    break;
                }
            }
            if(i == k){
                //y = x;
                mpz_set(y, x);
                k = 2*k;
            }
        }
        cout<<"the task is completed!"<<"\n";




    }

}

end of code

and this is a bit of what I get in the terminal before I kill the execution.

---

.
.
.
x=4
x²=16
x²-1=15
x = x²-1 mod numberToFactor=5
y-x=0
gcd(abs(y-x), numberToFactor)=10
x=5
x²=25
x²-1=24
x = x²-1 mod numberToFactor=4
y-x=1
gcd(abs(y-x), numberToFactor)=1
x=4
x²=16
x²-1=15
x = x²-1 mod numberToFactor=5
y-x=0
gcd(abs(y-x), numberToFactor)=10
x=5
x²=25
x²-1=24
x = x²-1 mod numberToFactor=4
y-x=1
gcd(abs(y-x), numberToFactor)=1
x=4
x²=16
x²-1=15
x = x²-1 mod numberToFactor=5
y-x=0
gcd(abs(y-x), numberToFactor)=10
x=5
x²=25
x²-1=24
x = x²-1 mod numberToFactor=4
y-x=1
gcd(abs(y-x), numberToFactor)=1
x=4
x²=16
x²-1=15
x = x²-1 mod numberToFactor=5
y-x=0
gcd(abs(y-x), numberToFactor)=10
x=5
x²=25
x²-1=24
x = x²-1 mod numberToFactor=4
y-x=1
gcd(abs(y-x), numberToFactor)=1
x=4
x²=16
x²-1=15
x = x²-1 mod numberToFactor=5
y-x=0
gcd(abs(y-x), numberToFactor)=10
x=5
x²=25
x²-1=24
x = x²-1 mod numberToFactor=4
y-x=1
gcd(abs(y-x), numberToFactor)=1
x=4
x²=16
x²-1=15
x = x²-1 mod numberToFactor=5
y-x=0
gcd(abs(y-x), numberToFactor)=10
x=5
x²=25
x²-1=24
x = x²-1 mod numberToFactor=4
y-x=1
gcd(abs(y-x), numberToFactor)=1
.
.
.

---

I am sorry of how the printouts look, I was in a bit of a hurry. I would be grateful for any help.

//Helen

Answer 1

using System;
using System.Numerics;

namespace Simple_Pollard_Rho_Floyd
{   class Program
    {   static void Main()
        {   /* Variables ****************************************************************/
            BigInteger f, n, x, y;

            /* Assign Data **************************************************************/
            n = 339850094323758691; // 764567807 * 444499613 (~71 ms)
            x = y = 2;
            Console.WriteLine("\n\t Simple Pollard Rho Floyd Factoring\n\n {0}", n);

            /* Algorithm ****************************************************************/
            do
            {   x = ((x * x) + 1) % n;
                y = ((y * y) + 1) % n;      // Floyd
                y = ((y * y) + 1) % n;
                f = BigInteger.GreatestCommonDivisor(BigInteger.Abs(y - x), n);
            } while (f == 1);

            /* Output Data **************************************************************/
            Console.WriteLine("\n factor = {0}\n", f);
            Console.Read();
        } // end of Main
    } // end Program
} // end namespace

Answer 2

I'm not sure if this is the only Problem, but your program is not computing the value of y.

You use x = f(x) = x²-1, so you have also to compute y = f(f(y)) = y⁴ - 2y².