Skip to content

Instantly share code, notes, and snippets.

@Hermann-SW
Last active Jun 22, 2021
Embed
What would you like to do?
Determines period of "x->(x^2+a) mod n" for start value s (arbitrary precision, gmplib)
/* 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;
}
@Hermann-SW

This comment has been minimized.

Copy link
Owner Author

@Hermann-SW Hermann-SW commented Jun 22, 2021

$ ./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

@Hermann-SW

This comment has been minimized.

Copy link
Owner Author

@Hermann-SW Hermann-SW commented Jun 22, 2021

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
$ 
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment