Общие соображения
Заметим, что:
-
Величины
С = lg (X * p1 * p2 * … * pk) и П = cnt(X) * lg p1 * lg p2 * … * lg pk
— это сумма и произведение величин (a1+1) lg p1, (a2+1) lg p2,…, (ak+1) lg pk.
Не будь задача дискретной, в соответствии с теоремами о средних оптимум достигался бы при равенстве этих величин, (а с учётом потенцирования – и величин вида piai+1). -
Наш случай – дискретное приближение к оптимуму, так что эти величины должны хотя бы стремиться к равенству. А поскольку основания растут, то показатели не должны возрастать (a1>=a2>=…>=ak).
-
Сумма и произведение величин идут “в одном флаконе”, и алгоритмы линейного динамического программирования здесь неприменимы.
-
Так как показатели степеней канонического разложения не возрастают, то требуемое оптимальное число X можно представить в виде произведения X=Pk…P2P1, где Pk=p1p2…pk— это произведения последовательных простых чисел (праймориалы), которые могут повторяться и которые можно протабулировать.
И поскольку праймориал при возрастании k растёт быстрее факториала, то их произведение для числа X подбирается очень быстро.
_Известно, что
Every highly composite number is a product of primorials_
Алгоритм и анализ результатов
Алгоритм по п.4 рекурсивный, параметром является частное от деления числа на уже готовую часть цепочки.
Для поиска текущего оптимального праймориала (звена цепочки) ведётся перебор по убыванию в границах от “жадного” звена для текущего частного до “жадного” звена для квадратного корня из этого частного.
Перебор обрывается, когда целевая функция (количество множителей) начинает падать.
Такая схема обеспечивает решение задачи.
Дополнительный контроль невозрастания звеньев в цепочке избавляет от повторения пройденного.
Для случая N=1018 оптимальное решение нашлось за 118 шагов и совпало с предыдущим ответом:
X = P12P4P2P2P1P1P1P1 = 28 * 34 * 52 * 72 * 11 * 13 … 37 = 897612484786617600,
cnt(X) = (8+1) (4+1) (2+1)2 (1+1)8 = 103680.
При этом порядок величин 29=512, 35=243, 53=125, 73=343, 112=121, …, 372=1369 примерно одинаков, что подтверждает теорию оптимума.
Для случая N=1036 (cnt(X)=127401984 множителей) понадобилось всего 168 шагов. Т.е. сложность алгоритма возрастает не быстрее логарифма. Хотя разрядность данных меняется.
Исходный код на PHP использует библиотеку BCMath. границу возможностей не тестировал (скорее всего, первым подведёт int для cnt(X)):
function prods(&$prod){
$p = array( "2", "3", "5", "7", "11", "13", "17", "19", "23", "29",
"31", "37", "41", "43", "47", "53", "59", "61", "67", "71",
"73", "79", "83", "89", "97",);
$str="1";
foreach($p as $val){
$str = bcmul($str,$val);
array_push($prod,$str);
}
return 1;
}
function greedy_prod($n){
global $p_prod;
foreach($p_prod as $key=>$val){
if(bccomp($val,$n)>0) break;
$i_m=$key+1;
}
if (!isset($i_m)) $i_m=0;
return $i_m;
}
function show_chain($chain){
$s = "[";
foreach($chain as $len) $s = $s . sprintf("%d.",$len);
$s = $s."]";
return $s;
}
function multi_prod($chain){
global $p_prod;
$m_prod = "1";
foreach($chain as $val) $m_prod= bcmul($m_prod, $p_prod[$val-1]);
return $m_prod;
}
function multi_count($chain){
$cnt = array();
for ($i=0; $i<$chain[0]; $i++) array_push($cnt, 0);
foreach($chain as $key=>$val){
for ($i=0; $i<$val; $i++) $cnt[$i]++;
}
$count=1;
foreach($cnt as $val) $count *= $val+1;
return $count;
}
function maxi_count($n, &$main_chain, &$latch_chain){
global $maxi_counter, $p_prod;
$maxi_counter++;
$m_prod = multi_prod($main_chain);
$nn = bcdiv($n,$m_prod);
$long_chain = greedy_prod($nn);
if(count($main_chain)>0){
$last = array_pop($main_chain); array_push($main_chain,$last);
$long_chain = min($long_chain, $last);
}
$short_chain = greedy_prod(bcsqrt($nn));
$short_chain = max($short_chain,1);
if ($long_chain==0) {
$max_count = multi_count($main_chain);
if($max_count> multi_count($latch_chain)) $latch_chain=$main_chain;
}
else {
$max_count = 0;
for ($cur_chain=$long_chain; $cur_chain>$short_chain-1; $cur_chain--){
array_push($main_chain,$cur_chain);
$counter = maxi_count($n, $main_chain, $latch_chain);
array_pop($main_chain);
if($counter>$max_count) $max_count = $counter;
else break;
}
}
$s1 = show_chain($main_chain); $s2=show_chain($latch_chain);
print("<tr><td align="right"> $m_prod</td><td align="right">$max_count</td><td>$s2</td><td>$s1</td></tr>");
return $max_count;
}
//$N="1000000000000000000";
$N="1000000000000000000000000000000000000";
$p_prod = array();
if (prods($p_prod)>0){
print("<br>p_prod:<br>");
foreach($p_prod as $key=>$val) printf(" %d=>%s", $key, $val);
}
$maxi_counter = 0; $chain=array(); $optimal_chain=array();
print("<br><br>maxi_count:<br>N=$N");
print('<table border="2">');
print('<tr><th>current number</th><th>multi_prod</th><th>optimal chain</th><th>current chain</th></tr>');
$max_count=maxi_count($N, $chain, $optimal_chain);
print('</table>');
printf("<br>Количество последовательностей = %d<br>Цепочка произведений = %s,<br>Kоличество множителей = %d<br>Минимальное число = %s", $maxi_counter, show_chain($optimal_chain), $max_count, multi_prod($optimal_chain));
Результаты для N = 1036:
<pre>
p_prod:
0=>2 1=>6 2=>30 3=>210 4=>2310 5=>30030 6=>510510 7=>9699690 8=>223092870
9=>6469693230 10=>200560490130 11=>7420738134810 12=>304250263527210
13=>13082761331670030 14=>614889782588491410 15=>32589158477190044730
16=>1922760350154212639070 17=>117288381359406970983270
18=>7858321551080267055879090 19=>557940830126698960967415390
20=>40729680599249024150621323470 21=>3217644767340672907899084554130
22=>267064515689275851355624017992790 23=>23768741896345550770650537601358310
24=>2305567963945518424753102147331756070
maxi_count:
N=1000000000000000000000000000000000000
current number multi_prod optimal chain current chain
713062256890366523119516128040749300 56623104 [24.3.] [24.3.]
855674708268439827743419353648899160 67108864 [24.2.2.] [24.2.2.]
570449805512293218495612902432599440 62914560 [24.2.2.] [24.2.1.1.]
285224902756146609247806451216299720 62914560 [24.2.2.] [24.2.1.]
142612451378073304623903225608149860 67108864 [24.2.2.] [24.2.]
23768741896345550770650537601358310 67108864 [24.2.2.] [24.]
616919031242227216631491481563344900 63700992 [24.2.2.] [23.5.]
673002579536975145416172525341830800 94371840 [23.4.2.1.] [23.4.2.1.]
336501289768487572708086262670915400 94371840 [23.4.2.1.] [23.4.2.]
897336772715966860554896700455774400 99090432 [23.4.1.1.1.1.] [23.4.1.1.1.1.]
448668386357983430277448350227887200 99090432 [23.4.1.1.1.1.] [23.4.1.1.1.]
224334193178991715138724175113943600 99090432 [23.4.1.1.1.1.] [23.4.1.1.]
112167096589495857569362087556971800 99090432 [23.4.1.1.1.1.] [23.4.1.]
56083548294747928784681043778485900 99090432 [23.4.1.1.1.1.] [23.4.]
961432256481393064880246464774044000 100663296 [23.3.3.1.1.] [23.3.3.1.1.]
480716128240696532440123232387022000 100663296 [23.3.3.1.1.] [23.3.3.1.]
240358064120348266220061616193511000 100663296 [23.3.3.1.1.] [23.3.3.]
576859353888835838928147878864426400 94371840 [23.3.3.1.1.] [23.3.2.2.1.]
288429676944417919464073939432213200 94371840 [23.3.3.1.1.] [23.3.2.2.]
769145805185114451904197171819235200 100663296 [23.3.3.1.1.] [23.3.2.1.1.1.1.]
384572902592557225952098585909617600 100663296 [23.3.3.1.1.] [23.3.2.1.1.1.]
192286451296278612976049292954808800 100663296 [23.3.3.1.1.] [23.3.2.1.1.]
96143225648139306488024646477404400 100663296 [23.3.3.1.1.] [23.3.2.1.]
48071612824069653244012323238702200 100663296 [23.3.3.1.1.] [23.3.2.]
8011935470678275540668720539783700 100663296 [23.3.3.1.1.] [23.3.]
267064515689275851355624017992790 100663296 [23.3.3.1.1.] [23.]
579755234179442444545257054963143400 84934656 [23.3.3.1.1.] [22.6.2.]
773006978905923259393676073284191200 95551488 [23.3.3.1.1.] [22.6.1.1.1.]
386503489452961629696838036642095600 95551488 [23.3.3.1.1.] [22.6.1.1.]
193251744726480814848419018321047800 95551488 [23.3.3.1.1.] [22.6.1.]
96625872363240407424209509160523900 95551488 [23.3.3.1.1.] [22.6.]
891931129506834530069626238404836000 113246208 [22.5.3.1.1.] [22.5.3.1.1.]
445965564753417265034813119202418000 113246208 [22.5.3.1.1.] [22.5.3.1.]
222982782376708632517406559601209000 113246208 [22.5.3.1.1.] [22.5.3.]
535158677704100718041775743042901600 106168320 [22.5.3.1.1.] [22.5.2.2.1.]
267579338852050359020887871521450800 106168320 [22.5.3.1.1.] [22.5.2.2.]
713544903605467624055700990723868800 113246208 [22.5.3.1.1.] [22.5.2.1.1.1.1.]
356772451802733812027850495361934400 113246208 [22.5.3.1.1.] [22.5.2.1.1.1.]
178386225901366906013925247680967200 113246208 [22.5.3.1.1.] [22.5.2.1.1.]
89193112950683453006962623840483600 113246208 [22.5.3.1.1.] [22.5.2.1.]
44596556475341726503481311920241800 113246208 [22.5.3.1.1.] [22.5.2.]
7432759412556954417246885320040300 113246208 [22.5.3.1.1.] [22.5.]
851388805438342051430097773022798000 104857600 [22.5.3.1.1.] [22.4.4.2.]
567592536958894700953398515348532000 100663296 [22.5.3.1.1.] [22.4.4.1.1.]
283796268479447350476699257674266000 100663296 [22.5.3.1.1.] [22.4.4.1.]
141898134239723675238349628837133000 104857600 [22.5.3.1.1.] [22.4.4.]
608134861027387179592926980730570000 98304000 [22.5.3.1.1.] [22.4.3.3.]
729761833232864615511512376876684000 113246208 [22.5.3.1.1.] [22.4.3.2.2.]
973015777643819487348683169168912000 125829120 [22.4.3.2.1.1.1.] [22.4.3.2.1.1.1.]
486507888821909743674341584584456000 125829120 [22.4.3.2.1.1.1.] [22.4.3.2.1.1.]
243253944410954871837170792292228000 125829120 [22.4.3.2.1.1.1.] [22.4.3.2.1.]
121626972205477435918585396146114000 125829120 [22.4.3.2.1.1.1.] [22.4.3.2.]
20271162034246239319764232691019000 125829120 [22.4.3.2.1.1.1.] [22.4.3.]
675705401141541310658807756367300 125829120 [22.4.3.2.1.1.1.] [22.4.]
3217644767340672907899084554130 125829120 [22.4.3.2.1.1.1.] [22.]
790130551223459534127080290097448600 71663616 [22.4.3.2.1.1.1.] [21.8.1.]
395065275611729767063540145048724300 71663616 [22.4.3.2.1.1.1.] [21.8.]
623787277281678579574010755340091000 84934656 [22.4.3.2.1.1.1.] [21.7.3.]
748544732738014295488812906408109200 99532800 [22.4.3.2.1.1.1.] [21.7.2.2.]
998059643650685727318417208544145600 111476736 [22.4.3.2.1.1.1.] [21.7.2.1.1.1.]
499029821825342863659208604272072800 111476736 [22.4.3.2.1.1.1.] [21.7.2.1.1.]
249514910912671431829604302136036400 111476736 [22.4.3.2.1.1.1.] [21.7.2.1.]
124757455456335715914802151068018200 111476736 [22.4.3.2.1.1.1.] [21.7.2.]
20792909242722619319133691844669700 111476736 [22.4.3.2.1.1.1.] [21.7.]
513707169526088242002126504397722000 94371840 [22.4.3.2.1.1.1.] [21.6.4.1.]
256853584763044121001063252198861000 94371840 [22.4.3.2.1.1.1.] [21.6.4.]
880640862044722700575074007538952000 123863040 [22.4.3.2.1.1.1.] [21.6.3.2.1.1.]
440320431022361350287537003769476000 123863040 [22.4.3.2.1.1.1.] [21.6.3.2.1.]
220160215511180675143768501884738000 123863040 [22.4.3.2.1.1.1.] [21.6.3.2.]
587093908029815133716716005025968000 113246208 [22.4.3.2.1.1.1.] [21.6.3.1.1.1.1.]
293546954014907566858358002512984000 113246208 [22.4.3.2.1.1.1.] [21.6.3.1.1.1.]
146773477007453783429179001256492000 113246208 [22.4.3.2.1.1.1.] [21.6.3.1.1.]
73386738503726891714589500628246000 113246208 [22.4.3.2.1.1.1.] [21.6.3.1.]
36693369251863445857294750314123000 123863040 [22.4.3.2.1.1.1.] [21.6.3.]
528384517226833620345044404523371200 111476736 [22.4.3.2.1.1.1.] [21.6.2.2.2.1.]
264192258613416810172522202261685600 111476736 [22.4.3.2.1.1.1.] [21.6.2.2.2.]
704512689635778160460059206031161600 119439360 [22.4.3.2.1.1.1.] [21.6.2.2.1.1.1.1.]
352256344817889080230029603015580800 119439360 [22.4.3.2.1.1.1.] [21.6.2.2.1.1.1.]
176128172408944540115014801507790400 119439360 [22.4.3.2.1.1.1.] [21.6.2.2.1.1.]
88064086204472270057507400753895200 119439360 [22.4.3.2.1.1.1.] [21.6.2.2.1.]
44032043102236135028753700376947600 119439360 [22.4.3.2.1.1.1.] [21.6.2.2.]
7338673850372689171458950062824600 119439360 [22.4.3.2.1.1.1.] [21.6.2.]
1223112308395448195243158343804100 123863040 [22.4.3.2.1.1.1.] [21.6.]
869350594582610871080521776673068000 100663296 [22.4.3.2.1.1.1.] [21.5.5.1.1.]
434675297291305435540260888336534000 100663296 [22.4.3.2.1.1.1.] [21.5.5.1.]
217337648645652717770130444168267000 100663296 [22.4.3.2.1.1.1.] [21.5.5.]
592739041760871048463992120458910000 98304000 [22.4.3.2.1.1.1.] [21.5.4.3.]
711286850113045258156790544550692000 113246208 [22.4.3.2.1.1.1.] [21.5.4.2.2.]
948382466817393677542387392734256000 125829120 [22.4.3.2.1.1.1.] [21.5.4.2.1.1.1.]
474191233408696838771193696367128000 125829120 [22.4.3.2.1.1.1.] [21.5.4.2.1.1.]
237095616704348419385596848183564000 125829120 [22.4.3.2.1.1.1.] [21.5.4.2.1.]
118547808352174209692798424091782000 125829120 [22.4.3.2.1.1.1.] [21.5.4.2.]
19757968058695701615466404015297000 125829120 [22.4.3.2.1.1.1.] [21.5.4.]
508062035795032327254850388964780000 106168320 [22.4.3.2.1.1.1.] [21.5.3.3.2.]
677416047726709769673133851953040000 117964800 [22.4.3.2.1.1.1.] [21.5.3.3.1.1.1.]
338708023863354884836566925976520000 117964800 [22.4.3.2.1.1.1.] [21.5.3.3.1.1.]
169354011931677442418283462988260000 117964800 [22.4.3.2.1.1.1.] [21.5.3.3.1.]
84677005965838721209141731494130000 117964800 [22.4.3.2.1.1.1.] [21.5.3.3.]
609674442954038792705820466757736000 115605504 [22.4.3.2.1.1.1.] [21.5.3.2.2.2.]
812899257272051723607760622343648000 127401984 [21.5.3.2.2.1.1.1.] [21.5.3.2.2.1.1.1.]
406449628636025861803880311171824000 127401984 [21.5.3.2.2.1.1.1.] [21.5.3.2.2.1.1.]
203224814318012930901940155585912000 127401984 [21.5.3.2.2.1.1.1.] [21.5.3.2.2.1.]
101612407159006465450970077792956000 127401984 [21.5.3.2.2.1.1.1.] [21.5.3.2.2.]
16935401193167744241828346298826000 127401984 [21.5.3.2.2.1.1.1.] [21.5.3.2.]
2822566865527957373638057716471000 127401984 [21.5.3.2.2.1.1.1.] [21.5.3.]
94085562184265245787935257215700 127401984 [21.5.3.2.2.1.1.1.] [21.5.]
40729680599249024150621323470 127401984 [21.5.3.2.2.1.1.1.] [21.]
746835726498886408969432055023615800 71663616 [21.5.3.2.2.1.1.1.] [20.9.2.]
995780968665181878625909406698154400 80621568 [21.5.3.2.2.1.1.1.] [20.9.1.1.1.]
497890484332590939312954703349077200 80621568 [21.5.3.2.2.1.1.1.] [20.9.1.1.]
248945242166295469656477351674538600 80621568 [21.5.3.2.2.1.1.1.] [20.9.1.]
124472621083147734828238675837269300 80621568 [21.5.3.2.2.1.1.1.] [20.9.]
974133556302895316047085289161238000 99532800 [21.5.3.2.2.1.1.1.] [20.8.3.2.]
649422370868596877364723526107492000 95551488 [21.5.3.2.2.1.1.1.] [20.8.3.1.1.]
324711185434298438682361763053746000 95551488 [21.5.3.2.2.1.1.1.] [20.8.3.1.]
162355592717149219341180881526873000 99532800 [21.5.3.2.2.1.1.1.] [20.8.3.]
779306845042316252837668231328990400 104509440 [21.5.3.2.2.1.1.1.] [20.8.2.2.1.1.]
389653422521158126418834115664495200 104509440 [21.5.3.2.2.1.1.1.] [20.8.2.2.1.]
194826711260579063209417057832247600 104509440 [21.5.3.2.2.1.1.1.] [20.8.2.2.]
519537896694877501891778820885993600 95551488 [21.5.3.2.2.1.1.1.] [20.8.2.1.1.1.1.]
259768948347438750945889410442996800 95551488 [21.5.3.2.2.1.1.1.] [20.8.2.1.1.1.]
129884474173719375472944705221498400 95551488 [21.5.3.2.2.1.1.1.] [20.8.2.1.1.]
64942237086859687736472352610749200 95551488 [21.5.3.2.2.1.1.1.] [20.8.2.1.]
32471118543429843868236176305374600 104509440 [21.5.3.2.2.1.1.1.] [20.8.2.]
5411853090571640644706029384229100 104509440 [21.5.3.2.2.1.1.1.] [20.8.]
657967402064236309961627783029959000 75497472 [21.5.3.2.2.1.1.1.] [20.7.5.]
717782620433712338139957581487228000 106168320 [21.5.3.2.2.1.1.1.] [20.7.4.2.1.]
358891310216856169069978790743614000 106168320 [21.5.3.2.2.1.1.1.] [20.7.4.2.]
957043493911616450853276775316304000 113246208 [21.5.3.2.2.1.1.1.] [20.7.4.1.1.1.1.]
478521746955808225426638387658152000 113246208 [21.5.3.2.2.1.1.1.] [20.7.4.1.1.1.]
239260873477904112713319193829076000 113246208 [21.5.3.2.2.1.1.1.] [20.7.4.1.1.]
119630436738952056356659596914538000 113246208 [21.5.3.2.2.1.1.1.] [20.7.4.1.]
59815218369476028178329798457269000 113246208 [21.5.3.2.2.1.1.1.] [20.7.4.]
512701871738365955814255415348020000 99532800 [21.5.3.2.2.1.1.1.] [20.7.3.3.1.]
256350935869182977907127707674010000 99532800 [21.5.3.2.2.1.1.1.] [20.7.3.3.]
615242246086039146977106498417624000 111476736 [21.5.3.2.2.1.1.1.] [20.7.3.2.2.1.]
307621123043019573488553249208812000 111476736 [21.5.3.2.2.1.1.1.] [20.7.3.2.2.]
820322994781385529302808664556832000 119439360 [21.5.3.2.2.1.1.1.] [20.7.3.2.1.1.1.1.]
410161497390692764651404332278416000 119439360 [21.5.3.2.2.1.1.1.] [20.7.3.2.1.1.1.]
205080748695346382325702166139208000 119439360 [21.5.3.2.2.1.1.1.] [20.7.3.2.1.1.]
102540374347673191162851083069604000 119439360 [21.5.3.2.2.1.1.1.] [20.7.3.2.1.]
51270187173836595581425541534802000 119439360 [21.5.3.2.2.1.1.1.] [20.7.3.2.]
8545031195639432596904256922467000 119439360 [21.5.3.2.2.1.1.1.] [20.7.3.]
284834373187981086563475230748900 119439360 [21.5.3.2.2.1.1.1.] [20.7.]
503151542755004237029480069375851000 67108864 [21.5.3.2.2.1.1.1.] [20.6.6.]
928895155855392437592886281924648000 110100480 [21.5.3.2.2.1.1.1.] [20.6.5.2.1.1.]
464447577927696218796443140962324000 110100480 [21.5.3.2.2.1.1.1.] [20.6.5.2.1.]
232223788963848109398221570481162000 110100480 [21.5.3.2.2.1.1.1.] [20.6.5.2.]
619263437236928291728590854616432000 100663296 [21.5.3.2.2.1.1.1.] [20.6.5.1.1.1.1.]
309631718618464145864295427308216000 100663296 [21.5.3.2.2.1.1.1.] [20.6.5.1.1.1.]
154815859309232072932147713654108000 100663296 [21.5.3.2.2.1.1.1.] [20.6.5.1.1.]
77407929654616036466073856827054000 100663296 [21.5.3.2.2.1.1.1.] [20.6.5.1.]
38703964827308018233036928413527000 110100480 [21.5.3.2.2.1.1.1.] [20.6.5.]
738893873975880348085250451530970000 92160000 [21.5.3.2.2.1.1.1.] [20.6.4.4.]
633337606265040298358786101312260000 106168320 [21.5.3.2.2.1.1.1.] [20.6.4.3.2.]
844450141686720397811714801749680000 117964800 [21.5.3.2.2.1.1.1.] [20.6.4.3.1.1.1.]
422225070843360198905857400874840000 117964800 [21.5.3.2.2.1.1.1.] [20.6.4.3.1.1.]
211112535421680099452928700437420000 117964800 [21.5.3.2.2.1.1.1.] [20.6.4.3.1.]
105556267710840049726464350218710000 117964800 [21.5.3.2.2.1.1.1.] [20.6.4.3.]
760005127518048358030543321574712000 115605504 [21.5.3.2.2.1.1.1.] [20.6.4.2.2.2.]
506670085012032238687028881049808000 113246208 [21.5.3.2.2.1.1.1.] [20.6.4.2.2.1.1.]
253335042506016119343514440524904000 113246208 [21.5.3.2.2.1.1.1.] [20.6.4.2.2.1.]
126667521253008059671757220262452000 115605504 [21.5.3.2.2.1.1.1.] [20.6.4.2.2.]
21111253542168009945292870043742000 115605504 [21.5.3.2.2.1.1.1.] [20.6.4.2.]
3518542257028001657548811673957000 117964800 [21.5.3.2.2.1.1.1.] [20.6.4.]
16754963128704769797851484161700 117964800 [21.5.3.2.2.1.1.1.] [20.6.]
557940830126698960967415390 119439360 [21.5.3.2.2.1.1.1.] [20.]
1 127401984 [21.5.3.2.2.1.1.1.] []
Количество последовательностей = 168
Цепочка произведений = [21.5.3.2.2.1.1.1.],
Kоличество множителей = 127401984
Минимальное число = 812899257272051723607760622343648000
</pre>Расширенная постановка задачи
Итак, мы нашли оптимальное решение задачи, которому соответствует минимальное среди n<N. В приведённом примере отношение N/n=1.23. Но если бы потребовалось найти все возможные числа до N с полученным количеством делителей, то возникла бы комбинаторная задача на тему возможных обменов простых множителей.
Оценка на количество делителей числа и сверхсоставные числа
Предыстория и мотивация
Иногда бывает так, что асимптотика решения задачи зависит от количества делителей числа во входе. Вы, возможно, слышали такое утверждение: «количество делителей числа $n$ — это примерно кубический корень из $n$». Действительно ли у чисел может быть так много делителей?
На самом деле, это, конечно, не так, но оценка в кубический корень дает нужный порядок величин для грубых оценок на числах, с которыми мы имеем дело в реальной жизни (отличие не больше, чем в $4$ раза при $n le 10^{15}$, и не больше, чем в $10$ раз при $n le 10^{18}$). Давайте разберемся, сколько все таки на самом деле делителей у чисел, и как этим пользоваться, а также придумаем новую более точную оценку.
Субполиномиальность и сверхсоставные числа
Будем обозначать количество делителей числа $n$ как $d(n)$. На самом деле $d(n)$ — это субполиномиальная величина, то есть для достаточно больших $n$ она меньше $n^{varepsilon}$ для сколь угодно малого положительного $varepsilon$. С доказательством можно, к примеру, ознакомиться, в этой статье. То есть для ооочень больших $n$ можно оценивать $d(n)$ не только как кубический корень, но и как корень четвертой, пятой… любой степени! Однако нас не очень сильно интересуют эти теоретические оценки, потому что по настоящему они достигаются только для невероятно огромных чисел, с которыми мы не работаем. Нас интересуют более практичные оценки.
Максимальное количество делителей неразрывно связано с таким понятием как «сверхсоставные числа» (highly composite numbers). Это такие числа, у которых больше делителей, чем у всех меньших чисел. Первые сверхсоставные числа — это $1, 2, 4, 6, 12, 24ldots$ Тогда, к примеру, если мы хотим понять, какое максимальное количество делителей есть у чисел, не превосходящих $1000$, можно посмотреть в список сверхсоставных чисел, найти там самое больше число, не превосходящее $1000$ (это будет $840$), и посмотреть, сколько у него делителей ($32$). На удивление сверхсоставные числа встречаются не очень часто: существует всего $156$ сверхсоставных чисел, не превосходящих $10^{18}$. К примеру, сверхсоставные числа можно найти как последовательность на OEIS: oeis.org/A002182.
Быстрая генерация больших сверхсоставных чисел
Для генерации сверхсоставных чисел есть этот замечательный скрипт. По ссылке есть скрипт, способный очень быстро вычислить все сверхсоставные числа до MAXN (при MAXN$=10^{100}$ программа работает меньше секунды), а также ниже представлен список всех всех сверхсоставных чисел до $10^{18}$ с их количествами делителей и разложениями.
Давайте разберемся, как же этому коду удается так быстро находить такие большие сверхсоставные числа. Для начала придумаем парочку очевидных алгоритмов поиска сверхсоставных чисел. Давайте идти по числам по возрастанию, находить их количества делителей, и если все меньшие числа имели меньше делителей, то новое число будет сверхсоставным. Количество делителей числа можно находить аналогично проверке на простоту за $O(sqrt n)$ (все делители разбиваются на пары вида ($x$, $frac{n}{x}$), и один из этих делителей будет не больше корня из $n$). В частности, из этого следует самая простая оценка на количество делителей числа: $d(n) le 2 cdot sqrt n$. Тогда мы найдем все сверхсоставные числа до n за $O(n cdot sqrt n)$, что очень долго.
Более умный способ основан на полезной идее о том, что вместо того, чтобы для каждого числа перебирать делители, мы можем действовать наоборот: для каждого числа будем перебирать, какие числа на него делятся. Эта идея похожа на идею решета Эратосфена и, как известно, работает за $O(n log n)$, потому что числа, делящиеся на $k$ — это $k$, $2k$, $3k$, $ldots$ Их $frac{n}{k}$ штук. Всего мы получим $frac{n}{1} + frac{n}{2} + ldots + frac{n}{n} = O(n log n)$ (сумма гармонического ряда). В частности, из этого следует, что в среднем у чисел от $1$ до $n$ как раз таки очень мало делителей: $O(log n)$. Но это все еще долго, такой алгоритм будет работать только для $n$ порядка $10^8$.
Для того, чтобы придумать эффективный алгоритм, давайте подробнее изучим нашу задачу. Давайте поймем, чему же равно количество делителей числа $n$. Пусть $n = p_1^{a_1} cdot p_2^{a_2} cdot ldots cdot p_k^{a_k}$. Тогда у $n$ есть ровно $(a_1 + 1) cdot (a_2 + 1) cdot ldots cdot (a_k + 1)$ делителей, потому что любой делитель $n$ должен состоять только из тех простых, которые есть в $n$, и при этом степени должны быть не больше, чем степени в $n$. Поэтому у нас есть $a_1 + 1$ вариантов степени $p_1$ (от $0$ до $a_1$), $a_2 + 1$ вариантов степени $p_2$ и так далее. Мы хотим, чтобы это произведение было как можно больше.
Несложно заметить, глядя на список по ссылке выше, что степени вхождения простых $2, 3, 5, 7 ldots$ в сверхсоставные числа не возрастают. Действительно, если есть, к примеру, число $40 = 2^3 cdot 3^0 cdot 5^1$, то можно поменять местами степени тройки и пятерки и получить число $24 = 2^3 cdot 3^1 cdot 5^0$, у которого будет ровно столько же делителей, но само число будет меньше.
Собственно, ровно на этой идее и основан алгоритм выше. Он начинает с числа $1$ и постепенно генерирует все возможные числа, не превосходящие MAXN, у которых степени вхождения простых не возрастают. Так как сумма степеней простых не больше логарифма MAXN, таких невозрастающих наборов степеней будет не очень много. После чего, зная разложение, легко посчитать для каждого из них количество делителей. Однако это еще не ответ. Это лишь кандидаты на то, чтобы быть сверхсоставными. Остается отсортировать их и пройтись по порядку, выбирая все числа, у которых больше делителей, чем у всех предыдущих.
Таблица для повседневного использования
Список сверхсоставных чисел — это, конечно, хорошо, но он не очень удобен для использования на практике. На практике у нас есть некоторое ограничение на $n$, и мы хотим узнать, какое максимальное количество делителей может быть у числа с такими ограничениями. На этот случай есть такая замечательная табличка:
| $le N$ | $n$ | Факторизация | $d(n)$ |
|---|---|---|---|
| 20 | 12 | $2^2 cdot 3^1$ | 6 |
| 50 | 48 | $2^4 cdot 3^1$ | 10 |
| 100 | 60 | $2^2 cdot 3^1 cdot 5^1$ | 12 |
| $10^3$ | 840 | $2^3 cdot 3^1 cdot 5^1 cdot 7^1$ | 32 |
| $10^4$ | 7560 | $2^3 cdot 3^3 cdot 5^1 cdot 7^1$ | 64 |
| $10^5$ | 83160 | $2^3 cdot 3^3 cdot 5^1 cdot 7^1 cdot 11^1$ | 128 |
| $10^6$ | 720720 | $2^4 cdot 3^2 cdot 5^1 cdot 7^1 cdot 11^1 cdot 13^1$ | 240 |
| $10^7$ | 8648640 | $2^6 cdot 3^3 cdot 5^1 cdot 7^1 cdot 11^1 cdot 13^1$ | 448 |
| $10^8$ | 73513440 | $2^5 cdot 3^3 cdot 5^1 cdot 7^1 cdot 11^1 cdot 13^1 cdot 17^1$ | 768 |
| $10^9$ | 735134400 | $2^6 cdot 3^3 cdot 5^2 cdot 7^1 cdot 11^1 cdot 13^1 cdot 17^1$ | 1344 |
| $10^{11}$ | 97772875200 | $2^6 cdot 3^3 cdot 5^2 cdot 7^2 cdot 11^1 cdot 13^1 cdot 17^1 cdot 19^1$ | 4032 |
| $10^{12}$ | 963761198400 | $2^6 cdot 3^4 cdot 5^2 cdot 7^1 cdot 11^1 cdot 13^1 cdot 17^1 cdot 19^1 cdot 23^1$ | 6720 |
| $10^{15}$ | 866421317361600 | $2^6 cdot 3^4 cdot 5^2 cdot 7^1 cdot 11^1 cdot 13^1 cdot 17^1 cdot 19^1 cdot 23^1 cdot 29^1 cdot 31^1$ | 26880 |
| $10^{18}$ | 897612484786617600 | $2^8 cdot 3^4 cdot 5^2 cdot 7^2 cdot 11^1 cdot 13^1 cdot 17^1 cdot 19^1 cdot 23^1 cdot 29^1 cdot 31^1 cdot 37^1$ | 103680 |
Исходный код таблицы в $LaTeX$ доступен по ссылке
По первому столбцу можно понять ограничения, по второму столбцу можно увидеть наименьшее число с наибольшим количеством делителей, потом его разложение на простые, и последний столбец, наконец, говорит нам то самое количество делителей. Если у вас нет этой таблицы, ее можно, к примеру, получить по двум последовательностям OEIS: oeis.org/A066150, oeis.org/A066151.
Численный анализ оценки в кубический корень
Для точности рекомендуется пользоваться этой таблицей и не полагаться на оценку в кубический корень, однако если этой таблички нет под рукой, то эта оценка даст вам примерное представление о максимальном количестве делителей.
Чтобы лучше понимать эту оценку, давайте проанализируем графики:
На данном графике на оси $X$ отложены десятичные логарифмы $n$ (то есть количество цифр в числе $n$), красный график — это кубический корень, а синий график — это максимальное количество делителей среди чисел, меньших $n$ (обратите внимание, что все числа по оси $Y$ умножаются на миллион). Из-за экспоненциального роста не видно, что происходит для чисел, меньших $10^{11}$, но и так видно, что для больших $n$ кубический корень уходит сильно выше, но все равно отношение не больше $10$ для обозримых чисел.
Давайте ограничим наш график на числа, меньшие $10^{12}$, чтобы посмотреть внимательнее на эту часть графика:
Видно, что для маленьких чисел количество делителей все таки больше кубического корня, а где-то в районе между $10^{10}$ и $10^{11}$ эти графики меняются местами и уже никогда не встречаются вновь. Как я говорил раньше, сильный всплеск отношения между этими графиками происходит уже на очень больших числах, и для $n le 10^{15}$ отношение между этими величинами не превосходит $4$.
Но если мы смотрим на такие графики, то мы видим только то, что происходит для больших чисел, потому что все маленькие величины становятся просто точкой. Для того, чтобы так не происходило, прологарифмируем также и ось $Y$:
Теперь происходящее стало еще понятнее.
Более точная практическая оценка
Как мы уже говорили раньше, оценка в кубический корень очень проста и дает правильный порядок величин, когда нет доступа к таблице. Кроме того, кубический корень достаточно легко вычислить даже без компьютера. Однако если у вас нет таблицы, но при этом есть компьютер, я рекомендую вам другую очень точную оценку на максимальное количество делителей: $exp((ln ln n)^{1.8})$.
Посмотрите, как удивительно точно следуют друг за другом зеленый и синий графики! До $10^{12}$ они идут настолько ровно рядом, что пропадает проблема с тем, что кубический корень достаточно сильно меньше $d(n)$ для малых $n$.
Наша новая оценка отличается от реального значения не больше, чем в $1.5$ раза для $n le 10^{14}$, и не больше, чем в $2.5$ раза для $n le 10^{18}$.
Еще одним удивительным фактом является то, что все три графика пересекаются примерно в одном и том же месте между $10^{11}$ и $10^{12}$.
$begingroup$
While playing around with programming, a problem suddenly came across my mind.
Given $n$ let’s say $2^{64}$ (64-bit unsigned integer). Find a positive integer $x$, $1 leq x leq n$, such that $x$ has maximum divisors.
My attempt was:
Let $n = p_1^{a_1} times p_2^{a_2} times … times p_k^{a_k}.$
The number of divisors is given by:
$$(a_1 + 1) times (a_2 + 1) times … times (a_k + 1)$$
We need to find the maximum of these products. Since $(a_1 + 1)_{max} implies {a_1}_{max}$, the product becomes: $a_1 times a_2 times …. times a_k$
Let $P = a_1 times a_2 times …. times a_k$.
The maximum of $P$ means $a_i = a_j$ where $1 leq i, j leq k$.
So I guess my question is how can we find a number that is less than $n$ which has this property?
Thanks,
asked Mar 11, 2011 at 17:17
$endgroup$
2
$begingroup$
Sequence A002182 gives the highly composite numbers, where the number of divisors sets a record.
You want the largest under $2^{64}$. Unfortunately the table doesn’t go that high, but the references do. Somewhat in the spirit of Mark Eicenlaub’s answer, if you think of $$n = p_1^{a_1} times p_2^{a_2} times … times p_k^{a_k} text{ and } d(n)=(a_1 + 1) times (a_2 + 1) times … times (a_k + 1)$$ and think of adding one to the exponent of $p_i$, $log (n)$ increases by $log (p_i)$ and $log (d(n))$ increases by $log(a_i+2)-log(a_i+1)$.
So the figure of merit for an increase is $$frac{log(a_i+2)-log(a_i+1)}{log (p_i)}$$ and you can just look through the primes to find which you should add in. This will miss some, where taking out one and adding another gets you a step up.
answered Mar 11, 2011 at 19:01
Ross MillikanRoss Millikan
369k27 gold badges252 silver badges444 bronze badges
$endgroup$
7
$begingroup$
I don’t know the exact answer to the question, but here is a way to get an approximate solution.
We are maximizing $a_1a_2a_3ldots$ subject to a constraint $p_1^{a_1}p_2^{a_2}p_3^{a_3}ldots leq n$. Let’s treat the $a_i$ as continuous variables that we can differentiate. In that case we might as well set the constraint $=n$ rather than $leq n$.
Take logarithms of both products. Now we are maximizing $ln a_1 + ln a_2 + ln a_3 ldots$ subject to $a_1 ln p_1 + a_2 ln p_2 + a_3 ln p_3 ldots = ln n$.
Introduce a Langrange multiplier and set the gradients of the two sums to be proportional to each other.
$$ frac{1}{a_1} hat{a_1} + frac{1}{a_2} hat{a_2} + frac{1}{a_3} hat{a_3 }+ ldots = lambda (ln p_1 hat{a_1} + ln p_2 hat{a_2} + ln p_3 hat{a_3} + ldots) $$
Dotting both sides with $hat{a_i}$ gives
$$frac{1}{a_i} = lambda ln p_i$$
or
$$a_i = frac{1}{lambda ln p_i}$$
$lambda$ is then chosen to satisfy the original constraint.
This should specify the maximum $a_i$ treating $a_i$ as continuous. To get the $a_i$ to be integers you would round some up and some down. Since you’re writing code, you could simply brute-force search this much smaller space of possible integer solutions.
answered Mar 11, 2011 at 18:33
Mark EichenlaubMark Eichenlaub
6,5433 gold badges32 silver badges40 bronze badges
$endgroup$
2
$begingroup$
This answer might be useful if you are loking for an actual algorithm to solve the problem (and not a thoretical approach).
If $x$ is an integer in the range $[1,n]$ with maximum number of divisors write
$$ x = q_1^{a_1} q_2^{a_2} dots q_k^{a_k} $$
with the $q_i$ distinct primes and
$$ a_1 ge a_2 ge cdots ge a_k. quad(*)$$
Then we see that the integer
$$ x’ = 2^{a_1} 3^{a_2} dots p_k^{a_k} quad (**)$$
has the same number of divisors and is $le x$. So we can limit our search to the integers of the form $(**)$ which verify $(*)$. This makes the search really easy, and fast in small ranges: start with $(a_1,dots,a_k) = (t,0,dots,0)$ and $k$ such that $2cdot3cdotdotscdot p_{k+1} > n$, and $2^t le n < 2^{t+1}$ and go through all the $k$-tuples verifying $(*)$ and $(**)$ in decreasing lexicographic order, keeping their product maximal and $<n$.
You can find a C++ program that performs this task in this thread at TopCoder, though you need to adapt it to arrive up to $2^{64}$, since it is limited to $nle 2^{64}/47$.
EDIT: I have written a PARI-GP version which can arrive much further, for $n=2^{64}$ the probably non-unique integer with maximum number of divisors is:
$$ 18401055938125660800 = 2^7cdot 3^4cdot 5^2cdot 7^2cdot 11 cdot 13 cdot 17 cdot 19 cdot 23 cdot 29 cdot 31 cdot 37 cdot 41 $$ with 184320 divisors.
answered Mar 11, 2011 at 19:00
$endgroup$
2
You must log in to answer this question.
Not the answer you’re looking for? Browse other questions tagged
.
Not the answer you’re looking for? Browse other questions tagged
.
Сайт переезжает. Большинство статей уже перенесено на новую версию.
Скоро добавим автоматические переходы, но пока обновленную версию этой статьи можно найти там.
Теория чисел
- Простые числа
- Разложение на простые множители
- Решето Эратосфена
- Линейное решето Эратосфена*
- НОД и НОК
- Алгоритм Евклида
- Расширенный алгоритм Евклида*
- Операции по модулю
- Быстрое возведение в степень
- Деление по простому модулю*
Простые числа
Простым называется натуральное число, которое делится только на единицу и на себя. Единица при этом простым числом не считается. Составным числом называют непростое число, которое еще и не единица.
Примеры простых чисел: (2), (3), (5), (179), (10^9+7), (10^9+9).
Примеры составных чисел: (4), (15), (2^{30}).
Еще одно определение простого числа: (N) — простое, если у (N) ровно два делителя. Эти делители при этом равны (1) и (N).
Проверка на простоту за линию
С точки зрения программирования интересно научиться проверять, является ли число (N) простым. Это очень легко сделать за (O(N)) – нужно просто проверить, делится ли оно хотя бы на одно из чисел (2, 3, 4, ldots, N-1) . (N > 1) является простым только в случае, если оно не делится на на одно из этих чисел.
def is_prime(n):
if n == 1:
return False
for i in range(2, n): # начинаем с 2, так как на 1 все делится; n не включается
if n % i == 0:
return False
return True
for i in range(1, 10):
print(i, is_prime(i))(1, False)
(2, True)
(3, True)
(4, False)
(5, True)
(6, False)
(7, True)
(8, False)
(9, False)Проверка на простоту за корень
Алгоритм можно ускорить с (O(N)) до (O(sqrt{N})).
Пусть (N = a times b), причем (a leq b). Тогда заметим, что (a leq sqrt N leq b).
Почему? Потому что если (a leq b < sqrt{N}), то (ab leq b^2 < N), но (ab = N). А если (sqrt{N} < a leq b), то (N < a^2 leq ab), но (ab = N).
Иными словами, если число (N) равно произведению двух других, то одно из них не больше корня из (N), а другое не меньше корня из (N).
Из этого следует, что если число (N) не делится ни на одно из чисел (2, 3, 4, ldots, lfloorsqrt{N}rfloor), то оно не делится и ни на одно из чисел (lceilsqrt{N}rceil + 1, ldots, N-2, N-1), так как если есть делитель больше корня (не равный (N)), то есть делитель и меньше корня (не равный 1). Поэтому в цикле for достаточно проверять числа не до (N), а до корня.
def is_prime(n):
if n == 1:
return False
# Удобно вместо for i in range(2, n ** 0.5) писать так:
i = 2
while i * i <= n:
if n % i == 0:
return False
i += 1
return True
for i in [1, 2, 3, 10, 11, 12, 10**9+6, 10**9+7]:
print(i, is_prime(i))(1, False)
(2, True)
(3, True)
(10, False)
(11, True)
(12, False)
(1000000006, False)
(1000000007, True)Разложение на простые множители
Любое натуральное число можно разложить на произведение простых, и с такой записью очень легко работать при решении задач. Разложение на простые множители еще называют факторизацией.
[11 = 11 = 11^1] [100 = 2 times 2 times 5 times 5 = 2^2 times 5^2] [126 = 2 times 3 times 3 times 7 = 2^1 times 3^2 times 7^1]
Рассмотрим, например, такую задачу:
Условие: Нужно разбить (N) людей на группы равного размера. Нам интересно, какие размеры это могут быть.
Решение: По сути нас просят найти число делителей (N). Нужно посмотреть на разложение числа (N) на простые множители, в общем виде оно выглядит так:
[N= p_1^{a_1} times p_2^{a_2} times ldots times p_k^{a_k}]
Теперь подумаем над этим выражением с точки зрения комбинаторики. Чтобы «сгенерировать» какой-нибудь делитель, нужно подставить в степень (i)-го простого число от 0 до (a_i) (то есть (a_i+1) различное значение), и так для каждого. То есть делитель (N) выглядит ровно так: [M= p_1^{b_1} times p_2^{b_2} times ldots times p_k^{b_k}, 0 leq b_i leq a_i] Значит, ответом будет произведение ((a_1+1) times (a_2+1) times ldots times (a_k + 1)).
Алгоритм разложения на простые множители
Применяя алгоритм проверки числа на простоту, мы умеем легко находить минимальный простой делитель числа N. Ясно, что как только мы нашли простой делитель числа (N), мы можем число (N) на него поделить и продолжить искать новый минимальный простой делитель.
Будем перебирать простой делитель от (2) до корня из (N) (как и раньше), но в случае, если (N) делится на этот делитель, будем просто на него делить. Причем, возможно, нам понадобится делить несколько раз ((N) может делиться на большую степень этого простого делителя). Так мы будем набирать простые делители и остановимся в тот момент, когда (N) стало либо (1), либо простым (и мы остановились, так как дошли до корня из него). Во втором случае надо еще само (N) добавить в ответ.
Напишем алгоритм факторизации:
def factorize(n):
factors = []
i = 2
while i * i <= n: # перебираем простой делитель
while n % i == 0: # пока N на него делится
n //= i # делим N на этот делитель
factors.append(i)
i += 1
# возможно, в конце N стало большим простым числом,
# у которого мы дошли до корня и поняли, что оно простое
# его тоже нужно добавить в разложение
if n > 1:
factors.append(n)
return factors
for i in [1, 2, 3, 10, 11, 12, 10**9+6, 10**9+7]:
print(i, '=', ' x '.join(str(x) for x in factorize(i)))1 =
2 = 2
3 = 3
10 = 2 x 5
11 = 11
12 = 2 x 2 x 3
1000000006 = 2 x 500000003
1000000007 = 1000000007Задание
За сколько работает этот алгоритм?
.
.
.
.
Решение
За те же самые (O(sqrt{N})). Итераций цикла while с перебором делителя будет не больше, чем (sqrt{N}). Причем ровно (sqrt{N}) операций будет только в том случае, если (N) – простое.
А итераций деления (N) на делители будет столько, сколько всего простых чисел в факторизации числа (N). Понятно, что это не больше, чем (O(log{N})).
Задание
Докажите, что число (N) имеет не больше, чем (O(log{N})) простых множителей в факторизации.
Разные свойства простых чисел*
Вообще, про простые числа известно много свойств, но почти все из них очень трудно доказать. Вот еще некоторые из них:
- Простых чисел, меньших (N), примерно (frac{N}{ln N}).
- N-ое простое число равно примерно (Nln N).
- Простые числа распределены более-менее равномерно. Например, если вам нужно найти какое-то простое число в промежутке, то можно их просто перебрать и проверить — через несколько сотен какое-нибудь найдется.
- Для любого (N ge 2) на интервале ((N, 2N)) всегда найдется простое число (Постулат Бертрана)
- Впрочем, существуют сколь угодно длинные отрезки, на которых простых чисел нет. Самый простой способ такой построить – это начать с (N! + 2).
- Есть алгоритмы, проверяющие число на простоту намного быстрее, чем за корень.
- Максимальное число делителей равно примерно (O(sqrt[3]{n})). Это не математический результат, а чисто эмпирический — не пишите его в асимптотиках.
- Максимальное число делителей у числа на отрезке ([1, 10^5]) — 128
- Максимальное число делителей у числа на отрекзке ([1, 10^9]) — 1344
- Максимальное число делителей у числа на отрезке ([1, 10^{18}]) — 103680
- Наука умеет факторизовать числа за (O(sqrt[4]{n})), но об этом как-нибудь в другой раз.
- Любое число больше трёх можно представить в виде суммы двух простых (гипотеза Гольдбаха), но это не доказано.
Решето Эратосфена
Часто нужно не проверять на простоту одно число, а найти все простые числа до (N). В этом случае наивный алгоритм будет работать за (O(Nsqrt N)), так как нужно проверить на простоту каждое число от 1 до (N).
Но древний грек Эратосфен предложил делать так:
Запишем ряд чисел от 1 до (N) и будем вычеркивать числа: * делящиеся на 2, кроме самого числа 2 * затем деляющиеся на 3, кроме самого числа 3 * затем на 5, затем на 7, и так далее и все остальные простые до n. Таким образом, все незачеркнутые числа будут простыми — «решето» оставит только их.
Красивая визуализация
Задание
Найдите этим способом на бумажке все простые числа до 50, потом проверьте с программой:
N = 50
prime = [1] * (N + 1)
prime[0], prime[1] = 0, 0
for i in range(2, N + 1): # можно и до sqrt(N)
if prime[i]:
for j in range(2 * i, N + 1, i): # идем с шагом i, можно начиная с i * i
prime[j] = 0
for i in range(1, N + 1):
if prime[i]:
print(i)2
3
5
7
11
13
17
19
23
29
31
37
41
43
47У этого алгоритма можно сразу заметить несколько ускорений.
Во-первых, число (i) имеет смысл перебирать только до корня из (N), потому что при зачеркивании составных чисел, делящихся на простое (i > sqrt N), мы ничего не зачеркнем. Почему? Пусть существует составное (M leq N), которое делится на %i%, и мы его не зачеркнули. Но тогда (i > sqrt N geq sqrt M), а значит по ранее нами доказанному утверждению (M) должно делиться и на простое число, которое меньше корня. Но это значит, что мы его уже вычеркнули.
Во-вторых, по этой же самое причине (j) имеет смысл перебирать только начиная с (i^2). Зачем вычеркивать (2i), (3i), (4i), …, ((i-1)i), если они все уже вычеркнуты, так как мы уже вычеркивали всё, что делится на (2), (3), (4), …, ((i-1)).
Асимптотика
Такой код будет работать за (O(N log log N)) по причинам, которые мы пока не хотим объяснять формально.
Гармонический ряд
Научимся оценивать асимптотику величины (1 + frac{1}{2} + ldots + frac{1}{N}), которая нередко встречается в задачах, где фигурирует делимость.
Возьмем (N) равное (2^i – 1) и запишем нашу сумму следующим образом: [left(frac{1}{1}right) + left(frac{1}{2} + frac{1}{3}right) + left(frac{1}{4} + ldots + frac{1}{7}right) + ldots + left(frac{1}{2^{i – 1}} + ldots + frac{1}{2^i – 1}right)]
Каждое из этих слагаемых имеет вид [frac{1}{2^j} + ldots + frac{1}{2^{j + 1} – 1} le frac{1}{2^j} + ldots + frac{1}{2^j} = 2^j frac{1}{2^j} = 1]
Таким образом, наша сумма не превосходит (1 + 1 + ldots + 1 = i le 2log_2(2^i – 1)). Тем самым, взяв любое (N) и дополнив до степени двойки, мы получили асимптотику (O(log N)).
Оценку снизу можно получить аналогичным образом, оценив каждое такое слагаемое снизу значением (frac{1}{2}).
Попытка объяснения асимптотики** (для старших классов)
Мы знаем, что гармонический ряд (1 + frac{1}{2} + frac{1}{3} + ldots + frac{1}{N}) это примерно (log N), а значит [N + frac{N}{2} + frac{N}{3} + ldots + frac{N}{N} sim N log N]
А что такое асимптотика решета Эратосфена? Мы как раз ровно (frac{N}{p}) раз зачеркиваем числа делящиеся на простое число (p). Если бы все числа были простыми, то мы бы как раз получили (N log N) из формули выше. Но у нас будут не все слагаемые оттуда, только с простым (p), поэтому посмотрим чуть более точно.
Известно, что простых чисел до (N) примерно (frac{N}{log N}), а значит допустим, что k-ое простое число примерно равно (k ln k). Тогда
[sum_{substack{2 leq p leq N \ text{p is prime}}} frac{N}{p} sim frac{1}{2} + sum_{k = 2}^{frac{N}{ln N}} frac{N}{k ln k} sim int_{2}^{frac{N}{ln N}} frac{N}{k ln k} dk =N(lnlnfrac{N}{ln N} – lnln 2) sim N(lnln N – lnlnln N) sim N lnln N]
Но вообще-то решето можно сделать и линейным.
Задание
Решите 5 первых задач из этого контеста:
https://informatics.msk.ru/mod/statements/view.php?id=34271
Линейное решето Эратосфена*
Наша цель — для каждого числа до (N) посчитать его минимальный простой делитель. Будем хранить его в массиве min_d. Параллельно будем хранить и список всех найденных простых чисел primes – это ровно те числа (x), у которых (min_d[x] = x).
Основное утверждение такое:
Пусть у числа (M) минимальный делитель равен (a). Тогда, если (M) составное, мы хотим вычеркнуть его ровно один раз при обработке числа (frac{M}{a}).
Мы также перебираем число (i) от (2) до (N). Если (min_d[i]) равно 0 (то есть мы не нашли ни один делитель у этого числа еще), значит оно простое – добавим в primes и сделаем (min_d[i] = i).
Далее мы хотим вычеркнуть все числа (i times k) такие, что (k) – это минимальный простой делитель этого числа. Из этого следует, что необходимо и достаточно перебрать (k) в массиве primes, и только до тех пор, пока (k < min_d[i]). Ну и перестать перебирать, если (i times k > N).
Алгоритм пометит все числа по одному разу, поэтому он корректен и работает за (O(N)).
N = 30
primes = []
min_d = [0] * (N + 1)
for i in range(2, N + 1):
if min_d[i] == 0:
min_d[i] = i
primes.append(i)
for p in primes:
if p > min_d[i] or i * p > N:
break
min_d[i * p] = p
print(i, min_d)
print(min_d)
print(primes)2 [0, 0, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
3 [0, 0, 2, 3, 2, 0, 2, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
4 [0, 0, 2, 3, 2, 0, 2, 0, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
5 [0, 0, 2, 3, 2, 5, 2, 0, 2, 3, 2, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0]
6 [0, 0, 2, 3, 2, 5, 2, 0, 2, 3, 2, 0, 2, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0]
7 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 0, 2, 0, 2, 3, 0, 0, 0, 0, 0, 3, 0, 0, 0, 5, 0, 0, 0, 0, 0]
8 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 0, 0, 0, 3, 0, 0, 0, 5, 0, 0, 0, 0, 0]
9 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 2, 0, 0, 3, 0, 0, 0, 5, 0, 3, 0, 0, 0]
10 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 2, 0, 2, 3, 0, 0, 0, 5, 0, 3, 0, 0, 0]
11 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 0, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 0, 5, 0, 3, 0, 0, 0]
12 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 0, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 2, 5, 0, 3, 0, 0, 0]
13 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 2, 5, 2, 3, 0, 0, 0]
14 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 0]
15 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 2]
16 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 0, 2, 0, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 2]
17 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 0, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 2]
18 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 0, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 2]
19 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 2]
20 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 2]
21 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 2]
22 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 0, 2, 5, 2, 3, 2, 0, 2]
23 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 0, 2]
24 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 0, 2]
25 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 0, 2]
26 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 0, 2]
27 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 0, 2]
28 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 0, 2]
29 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 29, 2]
30 [0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 29, 2]
[0, 0, 2, 3, 2, 5, 2, 7, 2, 3, 2, 11, 2, 13, 2, 3, 2, 17, 2, 19, 2, 3, 2, 23, 2, 5, 2, 3, 2, 29, 2]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]Этот алгоритм работает асимптотически быстрее, чем обычное решето. Но на практике, если писать обычное решето Эратсфена с оптимизациями, то оно оказывается быстрее линейнего. Также линейное решето занимает гораздо больше памяти – ведь в обычном решете можно хранить просто (N) бит, а здесь нам нужно (N) чисел и еще массив primes.
Зато один из «побочных эффектов» алгоритма — он неявно вычисляет факторизацию всех чисел от (1) до (N). Ведь зная минимальный простой делитель любого числа от (1) до (N) можно легко поделить на это число, посмотреть на новый минимальный простой делитель и так далее.
НОД и НОК
Введем два определения.
Наибольший общий делитель (НОД) чисел (a_1, a_2, ldots, a_n) — это максимальное такое число (x), что все (a_i) делятся на (x).
Наименьшее общее кратное (НОК) чисел (a_1, a_2, ldots, a_n) — это минимальное такое число (x), что (x) делится на все (a_i).
Например, * НОД(18, 30) = 6 * НОД(60, 180, 315) = 15 * НОД(1, N) = 1 * НОК(12, 30) = 6 * НОК(1, 2, 3, 4) = 12 * НОК(1, (N)) = (N)
Зачем они нужны? Например, они часто возникают в задачах.
Условие: Есть (N) шестеренок, каждая (i)-ая зацеплена с ((i-1))-ой. (i)-ая шестеренка имеет (a_i) зубчиков. Сколько раз нужно повернуть полносьтю первую шестеренку, чтобы все остальные шестеренки тоже вернулись на изначальное место?
Решение: Когда одна шестеренка крутится на 1 зубчик, все остальные тоже крутятся на один зубчик. Нужно найти минимальное такое число зубчиков (x), что при повороте на него все шестеренки вернутся в изначальное положение, то есть (x) делится на все (a_i), то есть это НОК((a_1, a_2, ldots, a_N)). Ответом будет (frac{x}{a_1}).
Еще пример задачи на применение НОД и НОК:
Условие: Город — это прямоугольник (n) на (m), разделенный на квадраты единичного размера. Вертолет летит из нижнего левого угла в верхний правый по прямой. Вертолет будит людей в квартале, когда он пролетает строго над его внутренностью (границы не считаются). Сколько кварталов разбудит вертолёт?
Решение: Вертолет пересечет по вертикали ((m-1)) границу. С этим ничего не поделать — каждое считается как новое посещение какого-то квартала. По горизонтали то же самое — ((n-1)) переход в новую ячейку будет сделан.
Однако еще есть случай, когда он пересекает одновременно обе границы (то есть пролетает над каким-нибудь углом) — ровно тот случай, когда нового посещения квартала не происходит. Сколько таких будет? Ровно столько, сколько есть целых решений уравнения (frac{n}{m} = frac{x}{y}). Мы как бы составили уравнение движения вертолёта и ищем, в скольки целых точках оно выполняется.
Пусть (t = НОД(n, m)), тогда (n = at, m = bt).
Тогда (frac{n}{m} = frac{a}{b} = frac{x}{y}). Любая дробь с натуральными числителем и знаменателем имеет ровно одно представление в виде несократимой дроби, так что (x) должно делиться на (a), а (y) должно делиться на (b). А значит, как ответ подходят ((a, b), (2a, 2b), (3a, 3b), cdots, ((t-1)a, (t-1)b)). Таких ответов ровно (t = НОД(n, m))
Значит, итоговый ответ: ((n-1) + (m-1) – (t-1)).
Кстати, когда (НОД(a, b) = 1), говорят, что (a) и (b) взаимно просты.
Алгоритм Евклида
Осталось придумать, как искать НОД и НОК. Понятно, что их можно искать перебором, но мы хотим хороший быстрый способ.
Давайте для начала научимся искать (НОД(a, b)).
Мы можем воспользоваться следующим равенством: [НОД(a, b) = НОД(a, b – a), b > a]
Оно доказывается очень просто: надо заметить, что множества общих делителей у пар ((a, b)) и ((a, b – a)) совпадают. Почему? Потому что если (a) и (b) делятся на (x), то и (b-a) делится на (x). И наоборот, если (a) и (b-a) делятся на (x), то и (b) делится на (x). Раз множства общих делитей совпадают, то и максимальный делитель совпадает.
Из этого равенства сразу следует следующее равенство: [НОД(a, b) = НОД(a, b operatorname{%} a), b > a]
(так как (НОД(a, b) = НОД(a, b – a) = НОД(a, b – 2a) = НОД(a, b – 3a) = ldots = НОД(a, b operatorname{%} a)))
Это равенство дает идею следующего рекурсивного алгоритма:
[НОД(a, b) = НОД(b operatorname{%} a, a) = НОД(a operatorname{%} , (b operatorname{%} a), b operatorname{%} a) = ldots]
Например: [НОД(93, 36) = ] [= НОД(36, 93spaceoperatorname{%}36) = НОД(36, 21) = ] [= НОД(21, 15) = ] [= НОД(15, 6) = ] [= НОД(6, 3) = ] [= НОД(3, 0) = 3]
Задание:
Примените алгоритм Евклида и найдите НОД чисел: * 1 и 500000 * 10, 20 * 18, 60 * 55, 34 * 100, 250
По-английски наибольший общий делитель — greatest common divisor. Поэтому вместо НОД будем в коде писать gcd.
def gcd(a, b):
if b == 0:
return a
return gcd(b, a % b)
print(gcd(1, 500000))
print(gcd(10, 20))
print(gcd(18, 60))
print(gcd(55, 34))
print(gcd(100, 250))
print(gcd(2465473782, 12542367456))1
10
6
1
50
6Вообще, в C++ такая функция уже есть в компиляторе g++ — называется __gcd. Если у вас не Visual Studio, то, скорее всего, у вас g++. Вообще, там много всего интересного.
А за сколько оно вообще работает?
Задание
Докажите, что алгоритм Евклида для чисел (N), (M) работает за (O(log(N+M))).
Кстати, интересный факт: самыми плохими входными данными для алгоритма Евклида являются числа Фибоначчи. Именно там и достигается логарифм.
Как выразить НОК через НОД
(НОК(a, b) = frac{ab}{НОД(a, b)})
По этой формуле можно легко найти НОК двух чисел через их произведение и НОД. Почему она верна?
Посмотрим на разложения на простые множители чисел a, b, НОК(a, b), НОД(a, b).
[ a = p_1^{a_1}times p_2^{a_2}timesldotstimes p_n^{a_n} ] [ b = p_1^{b_1}times p_2^{b_2}timesldotstimes p_n^{b_n} ] [ ab = p_1^{a_1+b_1}times p_2^{a_2+b_2}timesldotstimes p_n^{a_n+b_n} ]
Из определений НОД и НОК следует, что их факторизации выглядят так: [ НОД(a, b) = p_1^{min(a_1, b_1)}times p_2^{min(a_2, b_2)}timesldotstimes p_n^{min(a_n, b_n)} ] [ НОК(a, b) = p_1^{max(a_1, b_1)}times p_2^{max(a_2, b_2)}timesldotstimes p_n^{max(a_n, b_n)} ]
Тогда посчитаем (НОД(a, b) times НОК(a, b)): [ НОД(a, b)НОК(a, b) = p_1^{min(a_1, b_1)+max(a_1, b_1)}times p_2^{min(a_2, b_2)+max(a_2, b_2)}timesldotstimes p_n^{min(a_n, b_n)+max(a_n, b_n)} =] [ = p_1^{a_1+b_1}times p_2^{a_2+b_2}timesldotstimes p_n^{a_n+b_n} = ab]
Формула доказана.
Как посчитать НОД/НОК от более чем 2 чисел
Для того, чтобы искать НОД или НОК у более чем двух чисел, достаточно считать их по цепочке:
(НОД(a, b, c, d, ldots) = НОД(НОД(a, b), c, d, ldots))
(НОК(a, b, c, d, ldots) = НОК(НОК(a, b), c, d, ldots))
Почему это верно?
Ну просто множество общих делителей (a) и (b) совпадает с множеством делителей (НОД(a, b)). Из этого следует, что и множество общих делителей (a), (b) и еще каких-то чисел совпадает с множеством общих делителей (НОД(a, b)) и этих же чисел. И раз совпадают множества общих делителей, то и наибольший из них совпадает.
С НОК то же самое, только фразу “множество общих делителей” надо заменить на “множество общих кратных”.
Задание
Решите задачи F, G, H, I из этого контеста:
https://informatics.msk.ru/mod/statements/view.php?id=34271
Расширенный алгоритм Евклида*
Очень важным для математики свойством наибольшего общего делителя является следующий факт:
Для любых целых (a, b) найдутся такие целые (x, y), что (ax + by = d), где (d = gcd(a, b)).
Из этого следует, что существует решение в целых числах, например, у таких уравнений: * (8x + 6y = 2) * (4x – 5y = 1) * (116x + 44y = 4) * (3x + 11y = -1)
Мы сейчас не только докажем, что решения у таких уравнений существуют, но и приведем быстрый алгоритм нахождения этих решений. Здесь нам вновь пригодится алгоритм Евклида.
Рассмотрим один шаг алгоритма Евклида, преобразующий пару ((a, b)) в пару ((b, a operatorname{%} b)). Обозначим (r = a operatorname{%} b), то есть запишем деление с остатком в виде (a = bq + r).
Предположим, что у нас есть решение данного уравнения для чисел (b) и (r) (их наибольший общий делитель, как известно, тоже равен (d)): [bx_0 + ry_0 = d]
Теперь сделаем в этом выражении замену (r = a – bq):
[bx_0 + ry_0 = bx_0 + (a – bq)y_0 = ay_0 + b(x_0 – qy_0)]
Tаким образом, можно взять (x = y_0), а (y = (x_0 – qy_0) = (x_0 – (a operatorname{/} b)y_0)) (здесь (/) обозначает целочисленное деление).
В конце алгоритма Евклида мы всегда получаем пару ((d, 0)). Для нее решение требуемого уравнения легко подбирается — (d * 1 + 0 * 0 = d). Теперь, используя вышесказанное, мы можем идти обратно, при вычислении заменяя пару ((x, y)) (решение для чисел (b) и (a operatorname{%} b)) на пару ((y, x – (a / b)y)) (решение для чисел (a) и (b)).
Это удобно реализовывать рекурсивно:
def extended_gcd(a, b):
if b == 0:
return a, 1, 0
d, x, y = extended_gcd(b, a % b)
return d, y, x - (a // b) * y
a, b = 3, 5
res = extended_gcd(a, b)
print("{3} * {1} + {4} * {2} = {0}".format(res[0], res[1], res[2], a, b))3 * 2 + 5 * -1 = 1Но также полезно и посмотреть, как будет работать расширенный алгоритм Евклида и на каком-нибудь конкретном примере. Пусть мы, например, хотим найти целочисленное решение такого уравнения: [116x + 44y = 4] [(2times44+28)x + 44y = 4] [44(2x+y) + 28x = 4] [44x_0 + 28y_0 = 4] Следовательно, [x = y_0, y = x_0 – 2y_0] Будем повторять такой шаг несколько раз, получим такие уравнения: [116x + 44y = 4] [44x_0 + 28y_0 = 4, x = y_0, y = x_0 – 2y_0] [28x_1 + 16y_1 = 4, x_0 = y_1, y_0 = x_1 – y_1] [16x_2 + 12y_2 = 4, x_1 = y_2, y_1 = x_2 – y_2] [12x_3 + 4y_3 = 4, x_2 = y_3, y_2 = x_3 – y_3] [4x_4 + 0y_4 = 4, x_3 = y_4, y_3 = x_4 – 3 y_4] А теперь свернем обратно: [x_4 = 1, y_4 = 0] [x_3 = 0, y_3 =1] [x_2 = 1, y_2 =-1] [x_1 = -1, y_1 =2] [x_0 = 2, y_0 =-3] [x = -3, y =8]
Действительно, (116times(-3) + 44times8 = 4)
Задание
Решите задачу J из этого контеста:
https://informatics.msk.ru/mod/statements/view.php?id=34273
Операции по модулю
Выражение (a equiv b pmod m) означает, что остатки от деления (a) на (m) и (b) на (m) равны. Это выражение читается как «(a) сравнимо (b) по модулю (m)».
Еще это можно опрделить так: (a) сравнимо c (b) по модулю (m), если ((a – b)) делится на (m).
Все целые числа можно разделить на классы эквивалентности — два числа лежат в одном классе, если они сравнимы по модулю (m). Говорят, что мы работаем в «кольце остатков по модулю (m)», и в нем ровно (m) элементов: (0, 1, 2, cdots, m-1).
Сложение, вычитение и умножение по модулю определяются довольно интуитивно — нужно выполнить соответствующую операцию и взять остаток от деления.
С делением намного сложнее — поделить и взять по модулю не работает. Об этом подробнее поговорим чуть дальше.
a = 30
b = 50
mod = 71
print('{} + {} = {} (mod {})'.format(a, b, (a + b) % mod, mod))
print('{} - {} = {} (mod {})'.format(a, b, (a - b) % mod, mod)) # на C++ это может не работать, так как модуль от отрицательного числа берется странно
print('{} - {} = {} (mod {})'.format(a, b, (a - b + mod) % mod, mod)) # на C++ надо писать так, чтобы брать модулю от гарантированно неотрицательного числа
print('{} * {} = {} (mod {})'.format(a, b, (a * b) % mod, mod))
# print((a / b) % mod) # а как писать это, пока неясно30 + 50 = 9 (mod 71)
30 - 50 = 51 (mod 71)
30 - 50 = 51 (mod 71)
30 * 50 = 9 (mod 71)Задание
Посчитайте: * (2 + 3 pmod 5) * (2 * 3 pmod 5) * (2 ^ 3 pmod 5) * (2 – 4 pmod 5) * (5 + 5 pmod 6) * (2 * 3 pmod 6) * (3 * 3 pmod 6)
Для умножения (в C++) нужно ещё учитывать следующий факт: при переполнении типа всё ломается (разве что если вы используете в качестве модуля степень двойки).
intвмещает до (2^{31} – 1 approx 2 cdot 10^9).long longвмещает до (2^{63} – 1 approx 8 cdot 10^{18}).long long longв плюсах нет, при попытке заиспользовать выдает ошибкуlong long long is too long.- Под некоторыми компиляторами и архитектурами доступен
int128, но не везде и не все функции его поддерживают (например, его нельзя вывести обычными методами).
Зачем нужно считать ответ по модулю
Очень часто в задаче нужно научиться считать число, которое в худшем случае гораздо больше, чем (10^{18}). Тогда, чтобы не заставлять вас писать длинную арифметику, автор задачи часто просит найти ответ по модулю большого числа, обычно (10^9 + 7)
Кстати, вместо того, чтобы писать (1000000007) удобно просто написать (1e9 + 7). (1e9) означает (1 times 10^9)
int mod = 1e9 + 7; # В C++
cout << mod;1000000007N = 1e9 + 7 # В питоне такое число становится float
print(N)
print(int(N))1000000007.0
1000000007Быстрое возведение в степень
Задача: > Даны натуральные числа (a, b, c < 10^9). Найдите (a^b) (mod (c)).
Мы хотим научиться возводить число в большую степень быстро, не просто умножая (a) на себя (b) раз. Требование на модуль здесь дано только для того, чтобы иметь возможность проверить правильность алгоритма для чисел, которые не влезают в int и long long.
Сам алгоритм довольно простой и рекурсивный, постарайтесь его придумать, решая вот такие примеры (прямо решать необязательно, но можно придумать, как посчитать значение этих чисел очень быстро):
- (3^2)
- (3^4)
- (3^8)
- (3^{16})
- (3^{32})
- (3^{33})
- (3^{66})
- (3^{132})
- (3^{133})
- (3^{266})
- (3^{532})
- (3^{533})
- (3^{1066})
Да, здесь специально приведена такая последовательность, в которой каждое следующее число легко считается через предыдущее: его либо нужно умножить на (a=3), либо возвести в квадрат. Так и получается рекурсивный алгоритм:
- (a^0 = 1)
- (a^{2k}=(a^{k})^2)
- (a^{2k+1}=a^{2k}times a)
Нужно только после каждой операции делать mod: * (a^0 pmod c = 1) * (a^{2k} pmod c = (a^{k} pmod c)^2 pmod c) * (a^{2k+1} pmod c = ((a^{2k}pmod c) times a) pmod c)
Этот алгоритм называется быстрое возведение в степень. Он имеет много применений: * в криптографии очень часто надо возводить число в большую степень по модулю * используется для деления по простому модулю (см. далее) * можно быстро перемножать не только числа, но еще и матрицы (используется для динамики, например)
Асимптотика этого алгоритма, очевидно, (O(log c)) – за каждые две итерации число уменьшается хотя бы в 2 раза.
Задание
Решите задачу K из этого контеста:
https://informatics.msk.ru/mod/statements/view.php?id=34271
Задание
Решите как можно больше задач из практического контеста:
https://informatics.msk.ru/mod/statements/view.php?id=34273
Деление по модулю*
Давайте все-таки научимся не только умножать, но и делить по простому модулю. Вот только что это значит?
(a / b) = (a times b^{-1}), где (b^{-1}) – это обратный элемент к (b).
Определение: (b^{-1}) – это такое число, что (bb^{-1} = 1)
Утверждение: в кольце остатков по простому модулю (p) у каждого остатка (кроме 0) существует ровно один обратный элемент.
Например, обратный к (2) по модулю (5) это (3) ((2 times 3 = 1 pmod 5)))
Задание
Найдите обратный элемент к: * числу (3) по модулю (5) * числу (3) по модулю (7) * числу (1) по модулю (7) * числу (2) по модулю (3) * числу (9) по модулю (31)
Давайте докажем это утверждение: надо заметить, что если каждый ненулевой остаток (1, 2, ldots, (p-1)) умножить на ненулевой остаток (a), то получатся числа (a, 2a, ldots, (p-1)a) – и они все разные! Они разные, потому что если (xa = ya), то ((x-y)a = 0), а значит ((x – y) a) делится на (p), (a) – ненулевой остаток, а значит (x = y), и это не разные числа. И из того, что все числа получились разными, это все ненулевые, и их столько же, следует, что это ровно тот же набор чисел, просто в другом порядке!
Из этого следует, что среди этих чисел есть (1), причем ровно один раз. А значит существует ровно один обратный элемент (a^{-1}). Доказательство закончено.
Это здорово, но этот обратный элемент еще хочется быстро находить. Быстрее, чем за (O(p)).
Есть несколько способов это сделать.
Через малую теорему Ферма
Малая теорема Ферма: > (a^{p-1} = 1 pmod p), если (p) – простое, (a neq 0 pmod p)).
Доказательство: В предыдущем пункте мы выяснили, что множества чисел (1, 2, ldots, (p-1)) и (a, 2a, ldots, (p-1)a) совпадают. Из этого следует, что их произведения тоже совпадают по модулю: ((p-1)! = a^{p-1} (p-1)! pmod p).
((p-1)!neq 0 pmod p) а значит на него можно поделить (это мы кстати только в предыдущем пункте доказали, поделить на число – значит умножить на обратный к нему, который существует).
А значит, (a^{p – 1} = 1 pmod p).
Как это применить Осталось заметить, что из малой теоремы Ферма сразу следует, что (a^{p-2}) – это обратный элемент к (a), а значит мы свели задачу к возведению (a) в степень (p-2), что благодаря быстрому возведению в степень мы умеем делать за (O(log p)).
Обобщение У малой теоремы Ферма есть обобщение для составных (p):
Теорема Эйлера: > (a^{varphi(p)} = 1 pmod p), (a) – взаимно просто с (p), а (varphi(p)) – это функция Эйлера (количество чисел, меньших (p) и взаимно простых с (p)).
Доказывается теорема очень похоже, только вместо ненулевых остатков (1, 2, ldots, p-1) нужно брать остатки, взаимно простые с (p). Их как раз не (p-1), а (varphi(p)).
Для нахождения обратного по этой теореме достаточно посчитать функцию Эйлера (varphi(p)) и найти (a^{-1} = a^{varphi(p) – 1}).
Но с этим возникают большие проблемы: посчитать функцию Эйлера сложно. Более того, на предполагаемой невозможности быстро ее посчитать построены некоторые криптографические алгоритм типа RSA. Поэтому быстро делить по составному модулю этим способом не получится.
Через расширенный алгоритм Евклида
Этим способом легко получится делить по любому модулю! Рекомендую.
Пусть мы хотим найти (a^{-1} pmod p), (a) и (p) взаимно простые (а иначе обратного и не будет существовать).
Давайте найдем корни уравнения
[ax + py = 1]
Они есть и находятся расширенным алгоритмом Евклида за (O(log p)), так как (НОД(a, p) = 1), ведь они взаимно простые.
Тогда если взять остаток по модулю (p):
[ax = 1 pmod p]
А значит, найденный (x) и будет обратным элементом к (a).
То есть надо просто найти (x) из решения того уравнения по модулю (p). Можно брать по модулю прямо походу решения уравнения, чтобы случайно не переполниться.
В данной статье мы поговорим о том, как найти все делители числа. Начнем с доказательства теоремы, с помощью которой можно задать вид всех делителей определенного числа. Далее возьмем примеры нахождения всех нужных делителей и покажем, как именно определить, сколько делителей имеет конкретное число. В последнем пункте подробно рассмотрим примеры задач на нахождение общих делителей нескольких чисел.
Как найти все делители числа
Чтобы понять материал, изложенный в данном пункте, нужно хорошо знать, что вообще из себя представляют кратные числа и делители. Здесь мы поговорим только о поиске делителей натуральных чисел, т.е. целых положительных. Этим можно ограничиться, поскольку свойство делимости гласит, что делители целого отрицательного числа аналогичны делителям целого положительного, которое будет противоположным по отношению к этому числу. Также сразу уточним, что у нуля есть бесконечно большое число делителей, и находить их смысла не имеет, поскольку в итоге все равно получится 0.
Если речь идет о простом числе, то его можно разделить только на единицу и на само себя. Значит, у любого простого числа a есть всего 4 делителя, два из которых больше 0 и два меньше: 1, -1, a, -a. Возьмем простое число 7: у него есть делители 7, -7, 1 и -1, и все. Еще один пример: 367 – тоже простое число, которое можно разделить лишь на 1, -1, 367 и -367.
Сложнее определить все делители составного числа. Сформулируем теорему, которая лежит в основе данного действия.
Допустим, у нас есть выражение, означающее каноническое разложение числа на простые множители, вида a=p1s1·p2s2·…·pnsn. Тогда натуральными делителями числа a будут следующие числа: d=p1t2·p2t2·…·pntn, где t1=0, 1, …, s1, t2=0, 1, …, s2, …, tn=0, 1, …, sn.
Перейдем к доказательству этой теоремы. Зная основное определение делимости, мы можем утверждать, что a можно разделить на d, если есть такое число q, что делает верным равенство a=d·q, т.е. q=p1(s1−t1)·p2(s2-t2)·…·pn(sn-tn).
Любое число, делящее a, будет иметь именно такой вид, поскольку, согласно свойствам делимости, других простых множителей, кроме p1, p2, …, pn, оно иметь не может, а их показатели в данном случае не превысят s1, s2, …, sn.
Учитывая доказательство этой теоремы, мы можем сформировать схему нахождения всех положительных делителей данного числа.
Для этого нужно выполнить следующие действия:
- Выполнить каноническое разложение на простые множители и получить выражение вида a=p1s1·p2s2·…·pnsn.
- Найти все значения d=p1t2·p2t2·…·pntn, где числа t1, t2, …, tn будут принимать независимо друг от друга каждое из значений t1=0, 1, …, s1, t2=0, 1, …, s2, …, tn=0, 1, …, sn.
Самым трудным в таком расчете является именно перебор всех комбинаций указанных значений. Разберем подробно решения нескольких задач, чтобы наглядно показать применение данной схемы на практике.
Условие: найти все делители 8.
Решение
Разложим восьмерку на простые множители и получим 8=2·2·2. Переведем разложение в каноническую форму и получим 8=23. Следовательно, a=8, p1=2, s1=3.
Поскольку все делители восьмерки будут значениями p1t1=2t1, то t1 может принять значения нуля, единицы, двойки, тройки. 3 будет последним значением, ведь s1=3. Таким образом, если t1=0, то 2t1=20=1, если 1, то 2t1=21=2, если 2, то 2t1=22=4, а если 3, то 2t1=23=8.
Для нахождения делителей удобно все полученные значения оформлять в виде таблицы:
| t1 | 2t1 |
| 0 | 20=1 |
| 1 | 21=2 |
| 2 | 22=4 |
| 3 | 23=8 |
Значит, положительными делителями восьмерки будут числа 1, 2, 4 и 8, а отрицательными −1, −2, −4 и −8.
Ответ: делителями данного числа будут ±1, ±2, ±4, ±8.
Возьмем пример чуть сложнее: в нем при разложении числа получится не один, а два множителя.
Условие: найдите все делители числа 567, являющиеся натуральными числами.
Решение
Начнем с разложения данного числа на простые множители.
56718963217133337
Приведем разложение к каноническому виду и получим 567=34·7. Затем перейдем к вычислению всех натуральных множителей. Для этого будем присваивать t1 и t2 значения 0, 1, 2, 3, 4 и 0, 1, вычисляя при этом значения 3t1·7t2. Результаты будем вносить в таблицу:
| t1 | t2 | 3t1·7t2 |
| 0 | 0 | 30·70=1 |
| 0 | 1 | 30·71=7 |
| 1 | 0 | 31·70=3 |
| 1 | 1 | 31·71=21 |
| 2 | 0 | 32·70=9 |
| 2 | 1 | 32·71=63 |
| 3 | 0 | 33·70=27 |
| 3 | 1 | 33·71=189 |
| 4 | 0 | 34·70=81 |
| 4 | 1 | 34·71=567 |
Ответ: натуральными делителями 567 будут числа 27, 63, 81, 189, 1, 3, 7, 9, 21 и 567.
Продолжим усложнять наши примеры – возьмем четырехзначное число.
Условие: найти все делители 3 900, которые будут больше 0.
Решение
Проводим разложение данного числа на простые множители. В каноническом виде оно будет выглядеть как 3 900=22·3·52·13. Теперь приступаем к нахождению положительных делителей, подставляя в выражение 2t1·3t2·5t3·13t4 значения t1, равные 0, 1 и 2, t2=0,1, t3=0, 1, 2, t4=0, 1. Результаты представляем в табличном виде:
| t1 | t2 | t3 | t4 | 2t1·3t2·5t3·13t4 |
| 0 | 0 | 0 | 0 | 20·30·50·130=1 |
| 0 | 0 | 0 | 1 | 20·30·50·131=13 |
| 0 | 0 | 1 | 0 | 20·30·51·130=5 |
| 0 | 0 | 1 | 1 | 20·30·51·131=65 |
| 0 | 0 | 2 | 0 | 20·30·52·130=25 |
| 0 | 0 | 2 | 1 | 20·30·52·131=325 |
| 0 | 1 | 0 | 0 | 20·31·50·130=3 |
| 0 | 1 | 0 | 1 | 20·31·50·131=39 |
| 0 | 1 | 1 | 0 | 20·31·51·130=15 |
| 0 | 1 | 1 | 1 | 20·31·51·131=195 |
| 0 | 1 | 2 | 0 | 20·31·52·130=75 |
| 0 | 1 | 2 | 1 | 20·31·52·131=975 |
| t1 | t2 | t3 | t4 | 2t1·3t2·5t3·13t4 |
| 1 | 0 | 0 | 0 | 21·30·50·130=2 |
| 1 | 0 | 0 | 1 | 21·30·50·131=26 |
| 1 | 0 | 1 | 0 | 21·30·51·130=10 |
| 1 | 0 | 1 | 1 | 21·30·51·131=130 |
| 1 | 0 | 2 | 0 | 21·30·52·130=50 |
| 1 | 0 | 2 | 1 | 21·30·52·131=650 |
| 1 | 1 | 0 | 0 | 21·31·50·130=6 |
| 1 | 1 | 0 | 1 | 21·31·50·131=78 |
| 1 | 1 | 1 | 0 | 21·31·51·130=30 |
| 1 | 1 | 1 | 1 | 21·31·51·131=390 |
| 1 | 1 | 2 | 0 | 21·31·52·130=150 |
| 1 | 1 | 2 | 1 | 21·31·52·131=1950 |
| t1 | t2 | t3 | t4 | 2t1·3t2·5t3·13t4 |
| 2 | 0 | 0 | 0 | 22·30·50·130=4 |
| 2 | 0 | 0 | 1 | 22·30·50·131=52 |
| 2 | 0 | 1 | 0 | 22·30·51·130=20 |
| 2 | 0 | 1 | 1 | 22·30·51·131=260 |
| 2 | 0 | 2 | 0 | 22·30·52·130=100 |
| 2 | 1 | 0 | 1 | 22·30·52·131=1300 |
| 2 | 1 | 0 | 0 | 22·31·50·130=12 |
| 2 | 1 | 0 | 1 | 22·31·50·131=156 |
| 2 | 1 | 1 | 0 | 22·31·51·130=60 |
| 2 | 1 | 1 | 1 | 22·31·51·131=780 |
| 2 | 1 | 2 | 0 | 22·31·52·130=300 |
| 2 | 1 | 2 | 1 | 22·31·52·131=3900 |
Ответ: делителями числа 3 900 будут:195, 260, 300, 325, 390, 650, 780, 975, 75, 78, 100, 130, 150, 156, 13,15, 20, 25, 26, 30, 39, 50,52, 60, 65, 1, 2, 3, 4, 5, 6, 10, 12, 1 300, 1 950, 3 900
Как определить количество делителей конкретного числа
Чтобы узнать, сколько положительных делителей у конкретного числа a, каноническое разложение которого выглядит как a=p1s1·p2s2·…·pnsn, нужно найти значение выражения (s1+1) ·(s2+1) ·…·(sn+1). О количестве наборов переменных t1, t2, …, tn мы можем судить по величине записанного выражения.
Покажем на примере, как это вычисляется. Определим, сколько будет натуральных делителей у числа 3 900, которое мы использовали в предыдущей задаче. Каноническое разложение мы уже записывали: 3 900=22·3·52·13. Значит, s1=2, s2=1, s3=2, s4=1. Теперь подставим значения s1, s2, s3 и s4 в выражение (s1+1) ·(s2+1) ·(s3+1) ·(s4+1) и вычислим его значение. Имеем (2+1)·(1+1)·(2+1)·(1+1)=3·2·3·2=36. Значит, это число имеет всего 36 делителей, являющихся натуральными числами. Пересчитаем то количество, что у нас получилось в предыдущей задаче, и убедимся в правильности решения. Если учесть и отрицательные делители, которых столько же, сколько и положительных, то получится, что у данного числа всего будет 72 делителя.
Условие: определите, сколько делителей имеет 84.
Решение
Раскладываем число на множители.
844221712237
Записываем каноническое разложение: 84=22·3·7. Определяем, сколько у нас получится положительных делителей: (2+1)·(1+1)·(1+1) =12. Для учета отрицательных нужно умножить это число на 2:2·12=24.
Ответ: всего у 84 будет 24 делителя – 12 положительных и 12 отрицательных.
Как вычислить общие делители нескольких чисел
Зная свойства наибольшего общего делителя, можно утверждать, что количество делителей некоторого набора целых чисел будет совпадать с количеством делителей НОД тех же чисел. Это будет справедливо не только для двух чисел, но и для большего их количества. Следовательно, чтобы вычислить все общие делители нескольких чисел, надо определить их наибольший общий множитель и найти все его делители.
Разберем пару таких задач.
Условие: сколько будет натуральных общих делителей у чисел 140 и 50? Вычислите их все.
Решение
Начнем с вычисления НОД (140, 50).
Для этого нам потребуется алгоритм Евклида:
140=50·2+40, 50=40·1+10, 40=10·4, значит, НОД (50, 140)=10.
Далее выясним, сколько положительных делителей есть у десяти. Разложим его на простые множители и получим 20·50=1, 20·51=5, 21·50=2 и 21·51=10. Значит, все натуральные общие делители исходного числа – это 1, 2, 5 и 10, а всего их четыре.
Ответ: данные числа имеют четыре натуральных делителя, равные 10, 5, 2 и 1.
Условие: выясните, сколько общих положительных делителей есть у чисел 585, 315, 90 и 45.
Решение
Вычислим их наибольший общий делитель, разложив число на простые множители. Поскольку 90=2·3·3·5, 45=3·3·5, 315=3·3·5·7 и 585=3·3·5·13, то таким делителем будет 5: НОД (90, 45, 315, 585) =3·3·5=32·5.
Чтобы узнать количество этих чисел, нужно выяснить, сколько положительных делителей имеет НОД.
Считаем:
НОД (90, 45, 315, 585) =32·5:(2+1)·(1+1) =6.
Ответ: у данных чисел шесть общих делителей.
