Determines period of "x->(x^2+a) mod n" for start value s (arbitrary precision, gmplib)

period.c

/*  gcc -O3 -Wall -Wextra -pedantic -fomit-frame-pointer -m64 -o period period.c -lgmp

$ ./period 
period [n [a [s]]]
  Determines period of x->(x^2+a) mod n for start value s.
$
*/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <inttypes.h>
#include <time.h>
#include <sys/time.h>

#include <gmp.h>

void
determine_period_start_at(mpz_t n, mpz_t a, mpz_t s)
{
  unsigned long c = 0;
  mpz_t x, z, t;

  mpz_inits (t, NULL);
  mpz_init_set (x, s);
  mpz_init_set (z, s);

  do
  {
    ++c;

    mpz_mul (t, x, x);
    mpz_add (t, t, a);
    mpz_mod (x, t, n);

    mpz_mul (t, z, z);
    mpz_add (t, t, a);
    mpz_mod (z, t, n);

    mpz_mul (t, z, z);
    mpz_add (t, t, a);
    mpz_mod (z, t, n);
  }
  while (mpz_cmp (x, z) != 0);

  gmp_printf ("%ld [%Zd]", c, x);
}

int
main (int argc, char *argv[])
{
  mpz_t n,a,s;
  mpz_inits (n, a, s, NULL);

  if (argc > 1)
  {
    struct timespec t0, t1;

    mpz_set_str (n, argv[1], 0);
    mpz_set_str (a, argc>2 ? argv[2] : "1", 0);
    mpz_set_str (s, argc>3 ? argv[3] : "0", 0);

    gmp_printf ("n=%Zd a=%Zd s=%Zd: ", n, a, s);
    clock_gettime(CLOCK_REALTIME, &t0);
    determine_period_start_at (n, a, s);
    clock_gettime(CLOCK_REALTIME, &t1);

    printf(" %ldns\n", (t1.tv_sec-t0.tv_sec) * 1000000000 + t1.tv_nsec - t0.tv_nsec);
  }
  else
  {
    printf("period [n [a [s]]]\n");
    printf("  Determines period of x->(x^2+a) mod n for start value s.\n");
  }

  return 0;
}

$ ./period 
period [n [a [s]]]
  Determines period of x->(x^2+a) mod n for start value s.
$ 
$ ./period $((71*71*71*71))
n=25411681 a=1 s=0: 3881570 [2443684] 483587929ns
$ 
$ ./period $((71*71*71*71*71))
n=1804229351 a=1 s=0: 275591470 [1171381010] 37944932965ns
$ 

Default values a=1 and s=0 for sequence A248218:
https://oeis.org/A248218

Took less than an hour to compute a(128,100,283,921) = a(71^6) = 19,566,994,370 on Intel i7:

$ ./period $((71*71*71*71*71*71))
n=128100283921 a=1 s=0: 19566994370 [71536325699] 3528148502656ns
$