-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgenerate_abundant_carmichael_numbers.pl
293 lines (217 loc) · 17.1 KB
/
generate_abundant_carmichael_numbers.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
#!/usr/bin/perl
# Erdos construction method for Carmichael numbers:
# 1. Choose an even integer L with many prime factors.
# 2. Let P be the set of primes d+1, where d|L and d+1 does not divide L.
# 3. Find a subset S of P such that prod(S) == 1 (mod L). Then prod(S) is a Carmichael number.
# Alternatively:
# 3. Find a subset S of P such that prod(S) == prod(P) (mod L). Then prod(P) / prod(S) is a Carmichael number.
use 5.020;
use warnings;
use ntheory qw(:all);
use List::Util qw(shuffle);
use experimental qw(signatures);
# Modular product of a list of integers
sub vecprodmod ($arr, $mod) {
#~ if ($mod > ~0) {
#~ my $prod = Math::GMPz->new(1);
#~ foreach my $k(@$arr) {
#~ $prod = ($prod * $k) % $mod;
#~ }
#~ return $prod;
#~ }
if ($mod < ~0) {
my $prod = 1;
foreach my $k(@$arr) {
$prod = mulmod($prod, $k, $mod);
}
return $prod;
}
my $prod = 1;
foreach my $k (@$arr) {
$prod = Math::Prime::Util::GMP::mulmod($prod, $k, $mod);
}
Math::GMPz->new($prod);
}
# Primes p such that p-1 divides L and p does not divide L
sub lambda_primes ($L) {
#grep { ($L % $_) != 0 } grep { $_ > 2 and is_prime($_) } map { $_ + 1 } divisors($L);
grep { ($_ > 2) and (($L % $_) != 0) and is_prime($_) } map { ($_ >= ~0) ? (Math::GMPz->new($_)+1) : ($_ + 1) } divisors($L);
}
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89);
#my @prefix = (3,5,17,23, 29);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 617, 2003, 2549);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 353, 617);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 419, 449, 617);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003);
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 617, 1409, 2003, 2549, 3137, 9857, 10193, 16073, 68993, 202049, 275969, 1500929, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2549, 3137, 4019, 4289, 10193, 16073, 21977, 38459, 50513, 52529, 76649, 93809, 97553]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2549, 8009, 9857, 23297, 50513, 68993, 275969, 375233, 1500929, 3232769, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2297, 2549, 3329, 8009, 10193, 16073, 23297, 50177, 93809, 202049, 275969, 656657, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897, 1500929, 2475089, 8044037, 18386369]
foreach my $n(
#qw(327976184723917106922696038500316850197097366397296955613721508063986316871814680735 685142249888262836361512024427161900061736398403953340277064230345667415945220868055415 1447705574013899373231874907614593094830449009827553408005436718720395249892251694201091895 3154550445776286734272255423692198353635548392414238876043846610091741249515216441664179239205 7498366409610233567365151142116355486591698528768645808356223392188068950097669481835754051590285 19443264100119335640177836911507709776732274285097098581067687255943662787603256966400110255773609005 6752147266492648336087178462400067734188627847011899231029509812418215 8622492059311111925183326896484886496558877760634195318024684030458060555 11183372200926512166962774984740897786036864455542551327478015187504104539835 15757371431105455643250549953499924980525942017859454820416523399193283296627515 22895460689396227049643049082435390996704193751949787854065208499027840629999779295 4065021504672018226118572207226948264489296099895 2126006246943465532260013264379693942327901860245085 1226705604486379612114027653547083404723199373361414045 837839927864197275073880887372657965425945172005845792735 8400431334513514794031775534168508142396705 3334971239801865373230614887064897732531491885 1444042546834207706608856246099100718186135986205 11304017840824088923079122979486080913570205802842055 8692789719593724381847845571224796222535488262385540295 8075601649502569950736648535667835690735468595756166934055 9020447042494370634972836414340972466551518421459638465339435 13972672468823780113572923605814166350688302034840979982810784815 21895177758646863437968771290310798671528569288595815633064499805105 45739026337813297721916763225459258424823181243876658857471740092864345 108721665604982208684996146186916657275804701816694818104210326200738548065 207093460403925930428246384004848306979999376211770028014102314959732822821308183615 14116844787629511288659115050445171405 6338463309645650568607942657649881960845 )
#qw( 269616807542440095880838902934479765646815 39623553877269279114765582013213715065835439560755)
#qw(80457526688110964090941908021307082343637031298345387811989511961940980765)
#qw(156400339552952614158261462316373572492517490659424507963838717745)
#qw(30300520893268092108675921598650019435254761657479255573181729185)
#qw(23727894199896704861923196240133139730035052198495893166156405)
#qw(5560414251522943074360489383944833466511209118626077678018703909135127547157897173895)
#qw(14807987984119898011183681250338047627675234576733585022035)
#qw(1469693990779095715245182961386962460712186859726612101336192430649765)
#qw(10094456962219244277998564908104312836457170747058383045)
#qw(48976162706891785856850267556951183909509591565)
#qw(232256510885289043333840440634089959783515)
#qw(76901951875930663768872450745768943476147523769184184245)
#qw(93426807299183714704340654128147124952637393355)
qw(100567069213330155763552910794560952586261995)
#qw(2291345950693302094766820752641350945890207963412715)
#qw(327976184723917106922696038500316850197097366397296955613721508063986316871814680735 685142249888262836361512024427161900061736398403953340277064230345667415945220868055415 1447705574013899373231874907614593094830449009827553408005436718720395249892251694201091895 3154550445776286734272255423692198353635548392414238876043846610091741249515216441664179239205 7498366409610233567365151142116355486591698528768645808356223392188068950097669481835754051590285 19443264100119335640177836911507709776732274285097098581067687255943662787603256966400110255773609005 6752147266492648336087178462400067734188627847011899231029509812418215 8622492059311111925183326896484886496558877760634195318024684030458060555 11183372200926512166962774984740897786036864455542551327478015187504104539835 15757371431105455643250549953499924980525942017859454820416523399193283296627515 22895460689396227049643049082435390996704193751949787854065208499027840629999779295 80175683972122511383731338517508443634187270777707283918484112701538301509380205 125635296784315975338307007456935731174771453308667313900264604603310518465198781235 224007734166435384028201394295716408684617501249353820684171790007702654423449426942005 467952156673683517234912712683751577742165960109900131409234869326090845090585852881848445 988782907051493271917370561900767083769196673712218977667713278886029955676407907139345764285 2154557954465203839507950454381771475533079552018925152337947234692659273418892829656634420377015 5121384257763789526510398230065470797342130095148985087107300576864451092916708256093820017236164655 13279749380381506242241462610559765777508143336721318330869230395809521683933024508051275304693374950415 4065021504672018226118572207226948264489296099895 2126006246943465532260013264379693942327901860245085 1226705604486379612114027653547083404723199373361414045 837839927864197275073880887372657965425945172005845792735 19560340012826425693128516437578862635 8400431334513514794031775534168508142396705 3334971239801865373230614887064897732531491885 1444042546834207706608856246099100718186135986205 11304017840824088923079122979486080913570205802842055 8692789719593724381847845571224796222535488262385540295 8075601649502569950736648535667835690735468595756166934055 9020447042494370634972836414340972466551518421459638465339435 13972672468823780113572923605814166350688302034840979982810784815 21895177758646863437968771290310798671528569288595815633064499805105 45739026337813297721916763225459258424823181243876658857471740092864345 108721665604982208684996146186916657275804701816694818104210326200738548065 17277031840122951876810012573270045985 99436691737018840272248469415851811397274754265692619259832323874449185818279122553922438597548102645740113148322901071530942568237363765081022282605920356808140142566166169841228255378791958207751846269805134897319182266879091065192147432933981465795852869698389629543929 895029662324906581290508473212082154386870063145499265957750747193917121550330382107855869816530471914306758448054432544850014056704511249494281565735889131630069423238061694740895526664506415827974368274516019210769959584178698677794519043838767173628471680155205055524904929 8306770296037457980957209139881334474864541056053378687353884684706744805108616276343010327767219309836681025156393188448752980460274568906556427211594787030658674317072450588890251382973284045299430111955783174295155994900762502428610931245867598138445845663520458120326642646049 78058720471863992647054894287464900060302092303733599525064454382189280933605667148795268050028559854535291593394626791852931757385200124014910746507356213727099562557529818183801692245799950173678744762048494488851580884082465235321656920917417819706975611700101744956709460944922453 749441775250366193404374040053950505478960388208146289040143826523399286243548010295583368548324203163393334588181811828579997802655306390667158077217127007993882900114843784382680047251925321617489628460427595587464028068075748724323228097728128487006672847932676853329367534532200471253 207093460403925930428246384004848306979999376211770028014102314959732822821308183615 14116844787629511288659115050445171405 6338463309645650568607942657649881960845 83562900542610185478464992069940404422726948660826302385114906951342076455803610529344640789 50221303226108721472557460234034183058058896145156607733454059077756587949937969928136129114189 32191855367935690463909332010015911340215752429045385557144051868841972875910238723935258762195149 73447150516193077029324563940188527242280274750433560293457441337357999076671648487217809)
#qw(1568442344968558041435703706985298315 622671610952517542449974371673163431055 269616807542440095880838902934479765646815)
#qw(80457526688110964090941908021307082343637031298345387811989511961940980765)
#qw(183946261763273755985808210039437380929387193716936037921324457691523453096923545)
#qw(80457526688110964090941908021307082343637031298345387811989511961940980765)
#qw(68671670497867034860945549416314930790009427315 39623553877269279114765582013213715065835439560755 28806323668774765916434578123606370852862364560668885 22152062901287794989738190577053299185851158347154372565 20579266435296361545466779046082514943655726104506412112885 23727894199896704861923196240133139730035052198495893166156405 30300520893268092108675921598650019435254761657479255573181729185 39299775598568715464952670313449075207525425869750594478416702752945 55373383818383320090118312471649746967403325050478587620089134178899505 80457526688110964090941908021307082343637031298345387811989511961940980765 117387531437953896608684243803087033139366428664285920817692697952471890936135 183946261763273755985808210039437380929387193716936037921324457691523453096923545)
#qw(100567069213330155763552910794560952586261995)
#qw(39623553877269279114765582013213715065835439560755)
#qw(23727894199896704861923196240133139730035052198495893166156405 30300520893268092108675921598650019435254761657479255573181729185 39299775598568715464952670313449075207525425869750594478416702752945 55373383818383320090118312471649746967403325050478587620089134178899505 80457526688110964090941908021307082343637031298345387811989511961940980765)
) {
my @prefix = factor($n);
#my @prefix = (3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019);
#my @prefix = (5, 7, 13, 17, 19, 23, 37, 59, 67, 73, 89, 97, 109, 163, 193, 199, 233, 257, 349, 353, 397, 433, 487, 523, 577, 727, 769, 929, 1153, 1277, 1297, 1409, 1453, 1459, 1567, 1783, 2089, 2113, 2179, 2377, 2593);
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2297, 2549, 3137, 9857, 10193, 68993, 88397, 93809, 5850209, 8044037, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 2003, 2549, 3137, 4019, 4289, 10193, 16073, 21977, 38459, 50513, 52529, 76649, 93809, 97553]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2549, 8009, 9857, 23297, 50513, 68993, 275969, 375233, 1500929, 3232769, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 197, 353, 449, 617, 1409, 2003, 2297, 2549, 3329, 8009, 10193, 16073, 23297, 50177, 93809, 202049, 275969, 656657, 18386369]
#~ [3, 5, 17, 23, 29, 53, 83, 89, 113, 257, 353, 449, 617, 1409, 2003, 2297, 2549, 3137, 3329, 4019, 7547, 9857, 10193, 16073, 17837, 23297, 68993, 88397, 93809, 202049, 896897, 1500929, 2475089, 8044037, 18386369]
#my @prefix = (3, 5, 17, 23, 29);
my $prefix_prod = Math::GMPz->new(vecprod(@prefix));
my sub method_1 ($L) { # smallest numbers first
(vecall { ($L % ($_-1)) == 0 } @prefix) or return;
my @P = lambda_primes($L);
#@P = grep{$_ < 1e5 } @P;
#~ @P = grep {
#~ (($prefix_prod % $_) == 0)
#~ ? 1
#~ : Math::Prime::Util::GMP::gcd(Math::Prime::Util::GMP::totient(Math::Prime::Util::GMP::mulint($prefix_prod, $_)), Math::Prime::Util::GMP::mulint($prefix_prod, $_)) eq '1'
#~ } @P;
#vecprodmod(@P, 3*5*17*23) == 0 or return;
#vecprodmod(\@P, 3*5*17*23*29) == 0 or return;
if (@prefix) {
vecprodmod(\@P, $prefix_prod) == 0 or return;
}
#return if (vecprod(@P) < ~0);
#@P = grep { $_ > $prefix[-1] } @P;
@P = grep { gcd($prefix_prod, $_) == 1 } @P;
say "# Testing: $L -- ", scalar(@P);
my $n = scalar(@P);
my @orig = @P;
my $max = 1e6;
my $max_k = scalar(@prefix)>>1;
my $L_rem = invmod($prefix_prod, $L);
foreach my $k (1 .. @P) {
#next if (binomial($n, $k) > 1e6);
next if ($k > $max_k);
@P = @orig;
my $count = 0;
forcomb {
if (vecprodmod([@P[@_]], $L) == $L_rem) {
say vecprod(@P[@_], $prefix_prod);
}
lastfor if (++$count > $max);
} $n, $k;
next if (binomial($n, $k) < $max);
@P = reverse(@P);
$count = 0;
forcomb {
if (vecprodmod([@P[@_]], $L) == $L_rem) {
say vecprod(@P[@_], $prefix_prod);
}
lastfor if (++$count > $max);
} $n, $k;
@P = shuffle(@P);
$count = 0;
forcomb {
if (vecprodmod([@P[@_]], $L) == $L_rem) {
say vecprod(@P[@_], $prefix_prod);
}
lastfor if (++$count > $max);
} $n, $k;
}
my $B = vecprodmod([@P, $prefix_prod], $L);
my $T = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P));
foreach my $k (1 .. @P) {
#next if (binomial($n, $k) > 1e6);
last if ($k > $max_k);
@P = @orig;
my $count = 0;
forcomb {
if (vecprodmod([@P[@_]], $L) == $B) {
my $S = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
say vecprod($prefix_prod, $T / $S) if ($T != $S);
}
lastfor if (++$count > $max);
} $n, $k;
next if (binomial($n, $k) < $max);
@P = reverse(@P);
$count = 0;
forcomb {
if (vecprodmod([@P[@_]], $L) == $B) {
my $S = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
say vecprod($prefix_prod, $T / $S) if ($T != $S);
}
lastfor if (++$count > $max);
} $n, $k;
@P = shuffle(@P);
$count = 0;
forcomb {
if (vecprodmod([@P[@_]], $L) == $B) {
my $S = Math::GMPz->new(Math::Prime::Util::GMP::vecprod(@P[@_]));
say vecprod($prefix_prod, $T / $S) if ($T != $S);
}
lastfor if (++$count > $max);
} $n, $k;
}
}
use Math::GMPz;
#my %seen;
#my $count = 0;
#method_1(205058304)
method_1(carmichael_lambda($prefix_prod));
}
__END__
while (<>) {
chomp(my $n = $_);
#$n > 1e5 or next;
#$n *= 656656;
#$n *= 78848;
#$n *= 1232;
#$n = lcm($n, 29586292736);
#$n = lcm($n, 1232);
#$n = lcm($n, 5789168, 147090944);
#$n = lcm($n, 78848);
#$n = lcm($n, 656656);
#$n = lcm($n, 9193184);
#$n = lcm($n, 10506496);
#$n = lcm($n, 36772736);
#~ next if ($n < 707981814540);
#~ next if ($n > 44351949725003712);
#next if ($n < 1e6);
#next if ($n < 1e8);
#next if ($n < 7813080); # for 2^64
#next if ($n < ~0);
#next if ($n < 17125441200);
#next if (length($n) > 45);
# if (++$count % 1000 == 0) {
#say "Testing: $n";
# $count = 0 ;
# }
#say "Testing: $n";
if ($n > ~0) {
$n = Math::GMPz->new($n);
}
next if $seen{$n}++;
method_1($n);
}