Skip to content

Instantly share code, notes, and snippets.

@Hermann-SW
Last active August 12, 2021 15:20
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Hermann-SW/83deac0d8977f33b56b548641c26c590 to your computer and use it in GitHub Desktop.
Save Hermann-SW/83deac0d8977f33b56b548641c26c590 to your computer and use it in GitHub Desktop.
Analyze Pollard Rho cycles built by "x->(x^2+a) mod n"
/* gcc analyze.c -Wall -pedantic -Wextra -O3 -o analyze */
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>
typedef unsigned long unsig;
#define fmt "%lu"
unsig n,a,*A,*L,*M,l;
unsig dfs(unsig x)
{
unsig ret;
if (A[x]==0)
{
A[x]=l;
ret = dfs((x*x+a)%n);
}
else
{
if (A[x]==l)
{
unsig y=(x*x+a)%n, c=1;
M[y]=0;
while (y!=x)
{
y=(y*y+a)%n;
++c;
M[y]=0;
}
L[l]=c;
ret = l;
}
else
{
ret = A[x];
}
}
return ret;
}
unsig correct(unsig x, unsig c)
{
if (A[x]!=c)
{
A[x]=c;
M[x]=correct((x*x+a)%n, c)+1;
}
return M[x];
}
unsig correct2(unsig x)
{
if (M[x]==n)
{
M[x]=correct2((x*x+a)%n)+1;
}
return M[x];
}
int main(int argc, char *argv[])
{
unsig x;
assert(argc==3);
assert(strtol(argv[1], NULL, 0) >= 0);
assert(strtol(argv[2], NULL, 0) >= 0);
n=(unsig) strtol(argv[1], NULL, 0);
a=(unsig) strtol(argv[2], NULL, 0);
assert(n>2);
A=malloc(sizeof(unsig)*n); assert(A);
M=malloc(sizeof(unsig)*n); assert(M);
L=malloc(sizeof(unsig)*1000); assert(L);
l=1;
memset(A, 0, sizeof(unsig)*n);
for(x=0; x<n; ++x) { M[x]=n; }
memset(L, 0, 1000);
for(x=0; x<n; ++x)
{
if (A[x]==0)
{
unsig c = dfs(x);
if (c==l)
{
++l;
correct2(x);
}
else
{
correct(x, c);
}
}
}
for(x=1; x<l; ++x) { printf(fmt" ", L[x]); }
printf("\n");
for(x=0; x<n; ++x)
{
printf(fmt": "fmt" "fmt" "fmt"\n", x, A[x], M[x], L[A[x]]+M[x]);
}
{
unsig mm=0;
unsig sm=0;
unsig ml=0;
unsig sl=0;
unsig mlm=0;
unsig slm=0;
for(x=0; x<n; ++x)
{
sm+=M[x];
if (M[x]>mm) { mm=M[x]; }
sl+=L[A[x]];
if (L[A[x]]>ml) { ml=L[A[x]]; }
slm+=(L[A[x]]+M[x]);
if (L[A[x]]+M[x]>mlm) { mlm=L[A[x]]+M[x]; }
}
printf("c "fmt"\n", l-1);
printf("m "fmt" %.1f\n", mm, (1.0*sm)/n);
printf("l "fmt" %.1f\n", ml, (1.0*sl)/n);
printf("lm "fmt" %.1f\n", mlm, (1.0*slm)/n);
printf("lm1 "fmt"\n", L[A[1]]+M[1]);
printf("lm2 "fmt"\n", L[A[2]]+M[2]);
printf("lm3 "fmt"\n", L[A[3]]+M[3]);
printf("lm4 "fmt"\n", L[A[4]]+M[4]);
}
return 0;
}
@Hermann-SW
Copy link
Author

Hermann-SW commented Jun 20, 2021

Previous comment 357911 = 71^3.
If having 26.9g free ram, computing for n = 71^5 = 1804229351 takes some time, but works:

$ ./analyze 1804229351 1 | grep -v :
275591470 25053770 25053770 25053770 25053770 25053770 25053770 352870 3881570 25053770 352870 352870 352870 352870 352870 54670 352870 4970 4970 4970 4970 4970 4970 4970 770 70 11 70 70 70 70 70 70 5 2 
c 35
m 7 1.8
l 275591470 209183296.4
lm 275591477 209183298.2
lm1 275591476
lm2 275591475
lm3 25053770
lm4 275591470
$

@Hermann-SW
Copy link
Author

@Hermann-SW
Copy link
Author

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment