title: C++ - Brent-Pollard rho algorithm infinite loop
body: <p>I have this following function. I got it from two sites and tried to adapt it into my own but it didn't work very well.</p>
  
  <p>When I test <code>unsigned long int max - 2</code> or to but it as a number <code>4294967293</code> it puts the following code into an infinite loop where <code>ys</code> will keep coming back to the same value.</p>
  
  <p>I am very new to factorisation algorithms and step by step help into why I am getting the infinite loop would be great.</p>
  
  <p>The following code is just my "rho" function. I have another function called <code>gdc</code> which is the same as every other <code>gdc</code> recursive function out there.</p>
  
  <pre><code>unsigned long int rho(unsigned long int n) 
  {
      if (n == 1) return n;
      if (n % 2 == 0) return 2;
  
      unsigned long int y = rand() % n;
      unsigned long int x;
      unsigned long long int ys = y;
      unsigned long int c;
      do
          c = rand() % n;
      while (c == 0 || c == n - 2);
      unsigned long int m = 1000;
      unsigned long int d = 1;
      unsigned long int q = 1;
      unsigned long int r = 1;
      while (d == 1)
      {
          x = y;
          for (int i = 0; i &lt; r; i++)
          {
              y = y * y % n;
              y += c;
              if (y &lt; c)
                  y += (std::numeric_limits&lt;unsigned long&gt;::max() - n) + 1;
              y %= n;
          }
          int j = 0;
          while (j &lt; r &amp;&amp; d == 1)
          {
              ys = y;
              for (int i = 0; i &lt; m &amp;&amp; i &lt; (r-j); i++)
              {
                  y = y * y % n;
                  y += c;
                  if (y &lt; c)
                      y += (std::numeric_limits&lt;unsigned long&gt;::max() - n) + 1;
                  y %= n;
                  q *= ((x&gt;y) ? x - y : y - x) % n;
              }
              d = gcd(q, n);
              j += m;
          }
          r *= 2;
      }
      if (d == n)
      {
          do
          {
              ys = ys * ys % n;
              std::cout &lt;&lt; ys &lt;&lt; std::endl;
              ys += c;
              if (ys &lt; c)
                  ys += (std::numeric_limits&lt;unsigned long&gt;::max() - n) + 1;
              ys %= n;
              d = gcd( ((x&gt;ys) ? x - ys : ys - x) , n);
          } while (d == 1);
      }
      return d;
  }
  </code></pre>
  
  <p>The examples I adapted my code from :</p>
  
  <ul>
  <li><a href="http://xn--2-umb.com/09/12/brent-pollard-rho-factorisation" rel="nofollow">Brent-Pollard ρ factorisation</a></li>
  <li><a href="https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/" rel="nofollow">Pollard Rho Brent Integer Factorization</a></li>
  </ul>
  
  <h2>EDIT</h2>
  
  <p>I have done as Amd has suggested and tidied up my code and moved repetitive lines into helper functions. However I am still getting an infinite loop for the <code>d==n</code> part near the bottom. For some reason, <code>f(ys)</code> ends up essentially returning the same thing that it had returned earlier so it keeps cycling through a series of values.</p>
  
  <pre><code>uint64_t rho(uint64_t n)
  {
      if (n == 1) return n;
      if (n % 2 == 0) return 2;
  
      uint64_t y = rand() % n;
      uint64_t x;
      unsigned long long int ys = y;
      uint64_t c;
      do c = rand() % n; while (c == 0 || c == n - 2);
      uint64_t m = 1000;
      uint64_t d = 1;
      uint64_t q = 1;
      uint64_t r = 1;
      do 
      {
          x = y;
          for (int i = 0; i &lt;= r; i++)
              y = f(y, c, n);
  
          int j = 0;
          do 
          {
              ys = y;
              for (int i = 0; i &lt;= min(m, r-j); i++)
              {
                  y = f(y, c, n);
                  q *= (abs(x,y) % n);
              }
              d = gcd(q, n);
              j += m;
          } while (j &lt; r &amp;&amp; d == 1);
          r *= 2;
      } while (d == 1);
      if (d == n)
      {
          do
          {
              ys = f(ys, c, n);
              d = gcd(abs(x, ys), n);
          } while (d == 1);
      }
      return d;
  }
  </code></pre>