I’ve been reading some about the Erdos heuristic to construct Carmichael numbers, first described by Erdos in 1956. It is based on Korselt’s criterion, that N is a Carmichael number if N is composite, squarefree and such that if p|N, then (p-1)|(N-1). The best reference I’ve found to understand the basics of Erdos’s heuristic is a paper by Alford, Grantham, Hayman, and Shallue: “Constructing Carmichael numbers through improved subset-product algorithms”, at arXiv: https://arxiv.org/abs/1203.6664 .
The basic Pari/gp code is copied below with comments. However, I needed a more sophisticated method to find a Carmichael number for this value of Lambda:
? lambda=367567200 \\ Start with a colossally abundant number Lambda
%502 = 367567200
? factor(lambda) \\ Lambda has as prime factors the primes from 2 to 17
%503 =
[ 2 5]
[ 3 3]
[ 5 2]
[ 7 1]
[11 1]
[13 1]
[17 1]
? count=0;forprime(p=18,lambda+1,if((lambda%(p-1))==0,count=count+1;print(p," ",count))) \\ count p s.t. (p-1)|Lambda but p coprime to Lambda
[snip 344 lines] \\ There are a total of 344 such p
? vec=vector(count); \\ define the vector vec to contain the primes above; count=344.
? pr=1;count=0;forprime(p=18,lambda+1,if((lambda%(p-1))==0,count=count+1;pr=pr*p;vec[count]=p;print(p," ",count))); \\ populate the vector
[snip 344 lines]
? rec=0;for(iter=1,10000000, num=2^count+1+random(2^count);d=digits(num,2);pr=1;for(j=1,count,if(d[j+1],pr=pr*vec[j]));if((pr%lambda)==1,if(omega(pr)>rec,rec=omega(pr);print(pr," ",omega(pr))))) \\ search for a subset of the 344 primes with product 1 modulo Lambda
\\ Check that if pr is printed and has 2 or more prime factors, pr is a Carmichael number by Korselt's criterion:
\\ The primes p dividing pr are such that p-1 divides Lambda.
\\ Lambda divides pr-1 , so (p-1) divides (pr-1). Clearly pr is squarefree. If pr has 2 or more prime factors,
\\ then pr is a Carmichael number by Korselt's criterion.
---------------------- more advanced method, using a distance function -----------
omeg2(a)={indic=vector(7);for(j=1,7,if(!((a%Q[j])==1), indic[j]=1));m=0;for(j=1,7,if(indic[j],m=j));if(indic[1]&&indic[2]&&indic[3]&&indic[4]&&indic[5]&&indic[6]&&indic[7],m=0);return(m)}
This is from Definition 1.2 of the paper by Alford,Grantham, Hayman, and Shallue.
The idea is to find residues 'a' at a distance zero from 1 (mod Lambda), as measured by omeg2. One should think of the set S of primes among the 344 as varying over time, and "converging" on a solution:
? rec=10000000;for(iter=0,3700000, r=random(2^26);num=bitxor(num1d,r);d=digits(num,2);pr=1;for(j=1,count,if(d[j+1],pr=(pr*vec[j])));if(1,if(omeg2(pr)<1,if((
pr%lambda)<rec, rec=pr%lambda; print(pr%lambda," ",omeg2(pr)," ",r)))))
7490741 0 30039465
7138591 0 33364677
3345961 0 5701645
908291 0 57347721
470941 0 20734282
278821 0 3503003
207121 0 57979429
11281 0 39237696
8291 0 51175419
2231 0 20270848
1 0 27951752 \\ pr mod lambda is 1: pr could be a solution
\\ define: num1f=bitxor(num1d,27951752);
? rec=10000000;for(iter=0,0, r=0;num=bitxor(num1f,r);d=digits(num,2);pr=1;for(j=1,count,if(d[j+1],pr=(pr*vec[j])));if(1,if(omeg2(pr)<1,if((pr%lambda)<rec, r
ec=pr%lambda; print(pr%lambda," ",omeg2(pr)," ",pr)))))
1 0 1596801821246215215354150928470348085731488866618488240016890068110834294005836195843083760201819285043624222454304260383649471511260259062651806532465739265474163777909654773275087933709057884038920286248499087515898888381843855058557965381537860953388042519444732256055085357630776490339701159229635081392002550005267013162644954895574924103089191627344677142896313524681088763356747091773627633054616953589777173693008935824601362690678107598301015071365434109165379742620324834341210643035560978212679381806564772575121781922420769841930910895626114231655314947560769898207931472066357277272282989955029155024650476290253417897455317722337252200848545531234515826808372134772372224231068409475847704921561861994790639992927388759816654823564328959786592554476280164092163519201
time = 1 ms.
? pr
%611 = 1596801821246215215354150928470348085731488866618488240016890068110834294005836195843083760201819285043624222454304260383649471511260259062651806532465739265474163777909654773275087933709057884038920286248499087515898888381843855058557965381537860953388042519444732256055085357630776490339701159229635081392002550005267013162644954895574924103089191627344677142896313524681088763356747091773627633054616953589777173693008935824601362690678107598301015071365434109165379742620324834341210643035560978212679381806564772575121781922420769841930910895626114231655314947560769898207931472066357277272282989955029155024650476290253417897455317722337252200848545531234515826808372134772372224231068409475847704921561861994790639992927388759816654823564328959786592554476280164092163519201
? omega(pr)
time = 28 ms.
%612 = 185 \\ pr has 185 distinct prime factors
? pr%2
%613 = 1
? Mod(2,pr)^(pr-1) \\ Checking whether pr is a Fermat base 2 pseudoprime
time = 12 ms.
%614 = Mod(1, 1596801821246215215354150928470348085731488866618488240016890068110834294005836195843083760201819285043624222454304260383649471511260259062651806532465739265474163777909654773275087933709057884038920286248499087515898888381843855058557965381537860953388042519444732256055085357630776490339701159229635081392002550005267013162644954895574924103089191627344677142896313524681088763356747091773627633054616953589777173693008935824601362690678107598301015071365434109165379742620324834341210643035560978212679381806564772575121781922420769841930910895626114231655314947560769898207931472066357277272282989955029155024650476290253417897455317722337252200848545531234515826808372134772372224231068409475847704921561861994790639992927388759816654823564328959786592554476280164092163519201) \\ pr is a base 2 pseudoprime
? forprime(p=2,97, if(gcd(p,pr)==1, rem= lift( Mod(p,pr)^(pr-1) ); print(p," ",rem)))
2 1
3 1
5 1
7 1
11 1
13 1
17 1
37 1
41 1
47 1
53 1
59 1
73 1
83 1
97 1 \\ This means that pr is a base p pseudoprime for all primes p from 2 to 97
\\ that are coprime to pr.
Finally,
? for(iter=1,100,p=randomprime(pr); if(gcd(p,pr)==1, rem= lift( Mod(p,pr)^(pr-1) ); print(iter," ",rem)))
1 1
2 1
3 1
4 1 \\ This shows that pr is a base p pseudoprime for 100 random prime bases from
5 1 \\ 2 to pr.
6 1
7 1 \\ Conclusion: pr is probably a Carmichael number.
8 1
9 1
10 1
11 1
12 1
13 1
14 1
15 1
16 1
17 1
18 1
19 1
20 1
21 1
22 1
23 1
24 1
25 1
26 1
27 1
28 1
29 1
30 1
31 1
32 1
33 1
34 1
35 1
36 1
37 1
38 1
39 1
40 1
41 1
42 1
43 1
44 1
45 1
46 1
47 1
48 1
49 1
50 1
51 1
52 1
53 1
54 1
55 1
56 1
57 1
58 1
59 1
60 1
61 1
62 1
63 1
64 1
65 1
66 1
67 1
68 1
69 1
70 1
71 1
72 1
73 1
74 1
75 1
76 1
77 1
78 1
79 1
80 1
81 1
82 1
83 1
84 1
85 1
86 1
87 1
88 1
89 1
90 1
91 1
92 1
93 1
94 1
95 1
96 1
97 1
98 1
99 1
100 1
Addendum:
I chose a larger Lambda = 6,983,776,800
? log(pr)/log(10)
%724 = 2593.2872840518554785045839089406295506 \\ pr has 2594 digits
? for(n=1,100,p=randomprime(pr\10^2000);if(gcd(p,pr)==1,print(n," ",lift(Mod(p,pr)^(pr-1)))))
1 1
2 1
3 1 \\ pr is probably a Carmichael number with 549 prime factors
4 1
5 1
6 1
7 1
8 1
9 1
10 1
11 1
12 1
13 1
14 1
15 1
16 1
17 1
18 1
19 1
20 1
21 1
22 1
23 1
24 1
25 1
26 1
27 1
28 1
29 1
30 1
31 1
32 1
33 1
34 1
35 1
36 1
37 1
38 1
39 1
40 1
41 1
42 1
43 1
44 1
45 1
46 1
47 1
48 1
49 1
50 1
51 1
52 1
53 1
54 1
55 1
56 1
57 1
58 1
59 1
60 1
61 1
62 1
63 1
64 1
65 1
66 1
67 1
68 1
69 1
70 1
71 1
72 1
73 1
74 1
75 1
76 1
77 1
78 1
79 1
80 1
81 1
82 1
83 1
84 1
85 1
86 1
87 1
88 1
89 1
90 1
91 1
92 1
93 1
94 1
95 1
96 1
97 1
98 1
99 1
100 1
time = 42,009 ms.
? pr
%726 = 19376889019188860475758567487505704842683046126738922749256746129804275722743763143853003709024313490696819220866847652680920553692408893587159828887475685813220844239111193690447018821260621240658858329597752726220092392592598684285279158744217265809174012038328200551207562219853890768486270905661752410310013656382541283834397034052957540806955793482516594798786321762556886824684994671230076738649018398082038220013573467705837232670243490324041026613520364361599279976542247025410249979472280206477804978113169919349104178223757715984364157895303227690276593573651241576725784041398440144036329653672987923509502258980779900247857017440421168751190204448181089280625418993450500767573000108004093186753497088765550899711596127958661224385684398801660482504131544892569062561118261630967660101587099306578168113770744072219450898573970822424045606101804167038092643802215922523890401784850508163374914836286201279094340217553351164417456880173386882245474383822715206219252342731033879668819098316801433449270690976578483703337088282016358865448945883473403353063953403358458191582516954409747364784253325101066490191185714878476620273339177086215698314906050003282657763452304520498207265110601674056069141376453053738200114597257098045688497833969816190328804640040615736355883582486555542953218219848142716893501080625712765766603723917482673920577227144168017168327092815722603691426987661027607200742230019690928973677865225904604763572369100766603764573153553328705084508018877416781629612260580241920869948881001864604977748266502699401575264804931030975548191214509753606662540642250057106327325292100424757632178315790830730455300930088175013715507733803399882778982202692704041020777336271613965419964220418821183628500540385584490199948663914349773219445482546068233259968577545241919161251479573484874451102602938710809845297009908059024538206031571726385403028849433335144155400024758850195164054475578553780827074807748045598857223249855316059837495442616158527122286452971800936178356782833296586359454051929164489991307822124179533529674207889530939396071525186343021157751603403493291574439091928237658734043006532447410849487728268818988920410215780155595100654865196472812873679202349263619365098620018181137461556099765072354175857485098412284536010977943151053691292114658140755542980632118327148051752280372132477159715755851126061959750303344344184121215093680671026938782944126072309332111019393638999380869841697102828831467708939163648410917769814796381909184753614154630429283909819471943508031365754291322843319000498905465532045166819743815547981904193253433601