-
-
Save Hermann-SW/83deac0d8977f33b56b548641c26c590 to your computer and use it in GitHub Desktop.
/* 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; | |
} |
Now allows to compute values for all n<2^31 for CPUs with 64bit unsigned long type.
Functionality verified by comparing against table of first 100k entries of "Period in residues modulo n in iteration of x^2 + 1 starting at 0."
https://oeis.org/A248218
using this command:
for((i=1; i<=100000; ++i)); do echo $i
./analyze $i 1 | grep -v : | head -1 | cut -f1 -d\ ; done > 100k
Switch to "unsigned long" was needed for computing a(357911) = 54670
(357911 = 71^3).
Output for n=357911 from comments section reports period for starting at 0, followed by all other cycle periods. Number of cycles is c=19. m is maximal number of function evaluations until cycle is reached, l is maximal period length, lm is maximal sum over all starting values (average values are reported as well):
$ ./analyze 357911 1 | grep -v :
54670 4970 4970 4970 4970 4970 4970 70 770 4970 70 70 70 70 70 11 70 5 2
c 19
m 7 1.8
l 54670 41496.4
lm 54677 41498.2
lm1 54676
lm2 54675
lm3 4970
lm4 54670
$
Identifying values on period 2 cycle:
$ ./analyze 357911 1 | grep : | grep " 2$"
117047: 19 0 2
240863: 19 0 2
$
(117047^2+1)%357911 = 240863
(240863^2+1)%357911 = 117047
Identifying values on last cycle of length 70 (cycle with number 17):
$ ./analyze 357911 1 | grep ": 17 0" | wc --lines
70
$ ./analyze 357911 1 | grep ": 16 0" | wc --lines
11
$
Identifying values not on cycle with period 11, but entering it after some function applications:
$ ./analyze 357911 1 | grep ": 16 [1-9]" | wc --lines
41
$
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
$
Fix lint reported issues (revision 4):
https://www.raspberrypi.org/forums/viewtopic.php?f=33&t=317340&p=1899285#p1899265
Fix 36(!) clang-tidy reported issues:
https://www.raspberrypi.org/forums/viewtopic.php?f=33&t=317340&p=1901054#p1901054
https://www.raspberrypi.org/forums/viewtopic.php?f=33&t=312344