Determines for "x->(x^2+1)%n" number of steps needed by Pollard-Rho factorization algorithm, with s on cycle of start value 0

search_cycle.c

//  gcc -Wall -Wextra -pedantic search_cycle.c -o search_cycle
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

typedef unsigned long unsig;

unsig gcd(unsig a, unsig b)
{
  if (a<b)
  {
    if (a==0) return b; else b%=a;
  }
  while (b!=0)
  {
    a%=b;
    if (a==0) return b;
    b%=a;
  }
  return a;
}

unsig dif(unsig a, unsig b)
{
  return (a<b) ? b-a : a-b;
}

unsig f(unsig x, unsig a, unsig n)
{
  return (x*x + a) %n;
}

int main(int argc, char *argv[])
{
  unsig n=atoi(argv[1]);
  unsig s=atoi(argv[2]);
  unsig x,y,c;
  assert(argc == 3);

  x=s;
  do
  {
    for(y=f(x,1,n), c=1; gcd(dif(x, y), n)==1; y=f(y,1,n), ++c)  {}
    printf("%lu: %lu [%lu]\n",x,c, gcd(dif(x,y), n));
    x = f(x,1,n);
  }
  while (x!=s);
}

Works for n<2^32 on machines with "sizeof(unsigned long)==8" (because of square step in "f()").
Allows to be used for n=71^5 (for https://oeis.org/A248218):

$ echo "l(71^5)/l(2)" | bc -ql
30.74873559752341030734
$ 

Find values on start value 0 cycle using analyze.c:
https://gist.github.com/Hermann-SW/83deac0d8977f33b56b548641c26c590

$ ./analyze $((71*71*71)) 1  | grep ": 1 0 " | head -3
4: 1 0 54670
6: 1 0 54670
16: 1 0 54670
$ 

Take any of the values on cycle identified, here 6:

$ ./search_cycle $((71*71*71)) 6 | head
6: 11 [71]
37: 11 [71]
1370: 11 [71]
87346: 11 [71]
92841: 11 [71]
238580: 11 [71]
40516: 11 [71]
166411: 11 [71]
331030: 11 [71]
323764: 11 [71]
$ 

Period of cycle:

$ ./search_cycle $((71*71*71)) 6 | wc --lines
54670
$

For all values on cycle, Pollard-Rho factorization algorithm stops after 11 steps, with divisor 71 found:

$ ./search_cycle $((71*71*71)) 6 | cut -f2 -d: | sort -u
 11 [71]
$