Добавил:
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

книги / Mathematica 5. ╨б╨░╨╝╨╛╤Г╤З╨╕╤В╨╡╨╗╤М

.pdf
Скачиваний:
1
Добавлен:
19.11.2023
Размер:
33.75 Mб
Скачать

Функция A rra y [ f , 1 0 0 0 0 0 , 4] генерирует следую щ ий список из 100000 элементов

{f [ 4 ] ,

f [ 5 ] ,

...}. ф ункция D e le t e C a s e s [A rra y [ f , 1 0 0 0 0 0 , 4 ] , N u ll]

удалйет из

него элементы ,

равные N u ll. Затем функция S o r t

сортирует этот список, а функция

Union рассматривает его как множество и тем самым исключает

из него

повторяю ­

щиеся элементы .

 

 

 

Из

массива, который первоначально содержал

сто элементов,

получается список

всего из 15 элементов.

t l= U n io n [ S o r t [ D e le t e C a s e s [ A r r a y [ f ,1 0 0 , 4 ] , N u l l ] ] ]

{ 3 , 5 , 7 , 1 1 , 2 3 , 2 9 , 3 7 ,4 3 , 7 3 ,8 3 , 1 3 1 ,1 7 9 , 1 9 1 ,2 3 9 , 2 5 1 }

Понятно, что это очень мало. Н о даже если бы мы исходили из тысячи простых чисел, у нас бы получился список из 80 элементов.

t l= U n io n [ S o r t [ D e le t e C a s e s [ A r r a y [ f , 1 0 0 0 , 4 ] , N u l l ] ]]

{ 3 ,5 , 7 ,1 1 , 2 3 ,2 9 , 3 7 ,4 3 , 4 7 ,5 3 , 7 3 ,7 9 , 8 3 ,1 1 3 , 1 3 1 ,1 7 9 , 1 9 1 ,1 9 7 , 2 3 3 ,2 3 9 , 2 5 1 ,3 5 9 ,3 9 7 ,4 1 9 ,4 3 1 ,4 4 3 ,4 6 1 ,4 8 7 ,4 9 1 ,5 4 7 ,5 5 7 ,5 7 1 ,5 7 7 ,6 0 1 ,6 5 9 ,6 8 3 ,7 1 9 ,7 4 3 ,7 6 1 ,857, 9 1 1 ,9 4 1 , 1 0 1 3 , 1 0 1 9 , 1 0 3 1 ,1 1 0 3 , 1 2 2 3 , 1 3 2 1 ,1 4 3 9 , 1 4 5 1 , 14 9 9 , 1 5 1 1 ,1 5 5 9 , 15 8 3 ,1 8 1 1 ,1 9 3 1 ,2 0 0 3 ,2 0 3 9 ,2 0 6 3 ,2 3 3 9 ,2 3 5 1 ,2 3 9 9 ,2 4 5 9 ,2 5 4 3 ,2 6 9 9 ,2 8 1 9 ,2 9 0 3 ,2 9 3 9 ,2 9 6 3 ,3 0 2 3 ,3 2 9 9 ,3 3 5 9 ,3 4 9 1 ,3 5 3 9 ,3 6 2 3 ,3 7 7 9 ,3 8 0 3 ,3 8 5 1 ,3 8 6 3 ,3 9 1 1 }

Если ж е взять сто тысяч простых чисел,

получится список из 4536 элементов. К о­

нечно, это тож е не

очень много. Кроме того,

не

все элементы этого списка годятся

для отсеивания. Н е

годятся, например, р =

3 ,

5 ,

7. Ведь если простое число р попа­

ло в нашу таблицу, это означает лиш ь, что

М р делится на некоторое простое число,

возможно на себя, — и тогда он о является

простым числом!

П оэтом у для отсеивания

пригодны лиш ь те

элементы наш ей таблицы, для

которых

М р больш е самого боль­

шого простого модуля q , использованного

для построения наш ей таблицы. Правда,

число М ерсенна М

растет гораздо быстрее, чем

л -е простое число р п . П оэтом у для

большинства элем ентов из полученного списка выполняется неравенство М

> q , где

q — самый больш ой простой модуль, использованный для построения наш ей

табли­

цы. Однако проблема заключается

в том, что в наш ей

таблице слиш ком мало чисел.

Даже если для построения таблицы мы используем миллион простых

чисел, таблица

будет содержать всего лиш ь 37330

чисел. Это значит,

что с помощ ью

такой таблицы

мы сможем отбраковать м енее 4% чисел. Соответственно, и бы стродействие програм­ мы мы повысим м енее, чем на 4 %. Так что в данном случае услож нение программы

не оправдано, потому что для больш инства чисел (более 96%) будет применяться кри­ терий Люка—Л емера.

Но оказывается, эф ф ективны е критерии простоты сущ ествуют не только для чисел

Мерсенна! Эффективны е критерии сущ ествуют, например, и для чисел вида к - 2я + 1 .

Простые числа вида к •2" +1

Для чисел вида

к Т +1

сущ ествует весьма

эф ф ективны й

критерий

простоты .

Он состоит в проверке сравнения a r 'k = - 1 (m od

к - Т + 1) для

некоторого

а. Однако

выбор этого а зависит от к.

Оказывается, нуж но различать случаи, когда к

не делится

на 3 и делится на 3.

Случай, когда к не делится на 3, совсем прост.

 

Модулярная арифметика...

201

Случай, когда кне делится на 3

Если к не делится на 3, то, как заметил Прот (Proth) в 1878 году, числа Л-2" +1 яв­

ляются простыми при fc£2" + 1 тогда и только тогда, когда З2 * s - l (mod &-2"+l). Впрочем, эта традиционная формулировка критерия Прота несколько некорректна.

В ней, например, (молчаливо!) исключается случай п = 1, к =

1. Для п =

1, к = 1 име­

ем к - 1<3 = 2+1 = 2" + 1, хотя 3г"'*=3н0 (mod

*-2"+1),

причем

к - 2Я+ 1 = 3 -

простое число. Дело в том, что если число к ■2" +1

простое, то в традиционной фор-,

мулировке критерия Прота число 3 должно быть квадратичным невычетом по модулю

Jk• 2" +1. Но при п = 1, к ~ 1 число

к-2" -И = 1

Впрочем, п = 1, к =

1 — единствен­

ное решение уравнения *-2" + 1 =

3 в целых

числах, отличных от

0. (Уравнение

к •2" +1 = 3 в целых числах имеет всего два решения: указанное только что решение п — 1, к — 1 и и = 0, к = 2.)

Разыскивая простые числа вида к ■2" +1, конечно, нельзя ограничить перебор по п только простыми числами, но все равно, многие значения отсеиваются тривиальным образом.

Пример 7.7. Простые числа вида 5-2"+1. Рассмотрим, например, случай к = 5.1

Тогда числа вида 5-2*+1 делятся на 3 при четных л. Так что перебор можем вестиj

только по нечетным л. ] Используя функцию PowerMod, очень легко запрограммировать критерий простоты. |

M k n P rim etk _ ,n _ ]:«M odule( { t = k*2An + l),

I f lk < » 2 An+l,Pow erM od[ 3 , ( t - 1 ) / 2 ,t ] = = t - l ,P r i m e Q [ t ] ]]

А вот и программа для перебора по нечетным л.

M k n tk _ ,n l_ ]: « B lo c k [(n = l) .

W hile (n < « n l, {If[M k n P rim e[k ,n ], P r in t [ n ] ] , n + = 2 )]];

Теперь можем выполнить перебор.

Mkn(5,30Q00Q]

Так можно получить (если хватит терпения) следующие значения л: l , 3 ,

7, 13,

1 5 , 2 5 ,

3 9 , 55, 75,

85, 127, 1947, 3313, 4687, 5 9 4 7 , 1 3165,

23473,

2 6 6 0 7 ,

125413, 209787,

240937.

 

Впрочем, при дневном счете терпение обычно улетучивается после печати числа 13165, если не раньше. Конечно, метод перебора можно усовершенствовать, заметив,

что если остаток от деления и на 3 равен 2, то число 5-2’ +1 делится на 7. Так что числа, даюшие в остатке 2 при делении на 3, можно при переборе пропустить. Можно также пропустить числа, дающие в остатке 1 при делении на 10, поскольку тогда чис­ ло 5 2*+1 делится на 11. Можно, конечно, пойти и дальше. Можно, например, взять несколько небольших нечетных простых чисел д, для которых существует решение

сравнения 5-2“+1*0 (mod f). Пусть г — показатель 2 по модулю

q:

2's i (mod q).

(Такое s можно найти с помощью функции M u ltip lic a tiv e O r d e r [2 ,

q].) Тогда при

всех

(mod s) 5-2*+ls0 (mod g). Если g не является числом вида

5-2*+1, все чис­

ла мте (mod $) можно опустить при переборе. Если же д является числом м |и 5-2"+1,

то нужно оставить только то число яя/ (mod s), при котором 5-2* +1 = д. Эта идея вполне реализуема, но не слишком ли она усложняет перебор?

Глава7

Давайте попробуем. Вот функция, которая по заданному номеру п простого числа ^ находит / — решение сравнения 5-2'+1 = 0 (mod д )и 5 - показатель 2 по модулю q, если такое решение существует.

f [п_] : =

Block[{q= Prime[n], s=MultiplicativeOrder[2,q],

t=MultiplicativeOrder[2,q,{-PowerMod[5,-1,q]}]}, If[NumberQ[t],{s,t}]]

Теперь с помощью этой функции можно сгенерировать таблицу остатков и соот­ ветствующих им модулей.

t=Sort[DeleteCases[Array[f,100,4]/Null]]

{{3,2}, {10,1},{11,5},{12,9}, {18,11}, {20,3}, {23,14}, {28,20}, {36,31}, {50 ,17}, {51,22}, {52,31},{58,23},{60,8},{66,18},{82,14},{83, 56},{100,26},{ 106,6}, {130,19}, {131,111}, {138,121},{148,118},{162, 66},{172,47},{178,1 29},{180,114},{183,105},{191,124},{196,9},{200,52},{204,42},{210,183}, {224,64 },{226, 102},{231,139},{243, 117},{251,41},{260, 156}, {268, 194}, {2 92,265}, {316,231}, {346,242},{348,208},{372,17},{378,45},{388,110},{418 ,355}, {420, 352},{442,422},{4 60, 182},{4 66, 192},{490, 359}, {508,164}, {522 ,212},{540,314},{546,94},{556,305},{562,166}}

Чтобы в таблице сначала шли меньшие модули s (они помогают отбраковать больше кандидатов), таблицу пришлось отсортировать, предварительно вычеркнув из нее эле­ менты N U L L , которые соответствуют тем простым модулям q, для которых сравнение

5-2'+ 1 = 0 (mod q) неразрешимо.

Каждая пара в этой таблице описывает арифметическую прогрессию из показате­ лей, для которых выполняется некоторое тождественное сравнение. Пусть, например, в таблице имеется пара (s, /). Она была получена для некоторого простого числа q.

Ее наличие в таблице означает, что 5• 2Ь+' + 1= 0 (mod q) для любых натуральных к.

Теперь давайте не будем спешить с написанием программы, а предварительно оце­ ним количество чисел, которые может помочь отбраковать такая таблица. Для этого сначала составим программу, которая подсчитывает количество чисел, отбракованных с помощью таблицы. Для этой программы нам понадобится функция, отбраковываю­ щая заданное число по заданной таблице. Вот код этой функции.

fx[n_, : =

Block[{1= Length[t],i=l,c=(i<=l)}, While[c,

{s,r}=t [[i]]; answer=(Mod[n,s]!=r); i++; c=(i<=l)&&answer];

answer]

Эта функция принимает значение True, если число п не было отсеяно с помощью таблицы t, т.е. если оно подлежит дальнейшему испытанию. Теперь можем составить программу подсчета отбракованных чисел.

R<rog:=

Block[ {nn=200000,Yes=0,No=0,j=l},

While[j<=nn,

If[fx[j,t ] ,Yes++,No++];

j+=2 ];

Print["Yes=",Yes,No=",No]]

Модулярная арифметика...

203

Конечно, количество отбракованных чисел и время отбраковки зависят от размера таблицы /. Сначала давайте возьмем совсем небольшую таблицу.

t=Sort[DeleteCases[Array[f,10,4]/Null]] {{3,2},{10,1},{11,5},{12,9},{18,11},{20, 3}, {28,20}, {36,31}}

Теперь запускаем программу.

prog//Timing

Yes= 27274 ; No= 72726 {5.188 Second,Null}

Как видите, за 5 секунд удалось отбраковать более 72 % чисел! Давайте попробуем увеличить таблицу.

t=Sort[DeleteCases[Array[f,25,4],Null]]

{{3,2},{10,1},{11,5},{12,9},{18,11},{20,3}, {23,14}, {28,20},{36,31},{51 ,22},{52,31},{58,23},{60,8},{66,18},{82,14},{100,26},{106,6}}

Опять считаем процент отбракованных чисел и замеряем время:

prog//Timing

Yes= 23152 ; No= 76848 {6.985 Second,Null}

Время отбраковки возросло менее, чем на 2 секунды, а отбраковали мы почти 77 % чисел! Значит, таблицу стоит увеличить.

t=Sort[DeleteCases[Array[f,100,4],Null]];

Опять измеряем время:

prog//Timing

Yes= 1,8329 ; No= 81671 {14.047 Second,Null}

Как видим, таблицу увеличивать стоило! Время отбраковки, конечно, возросло, притом более чем вдвое, но 14 секунд по сравнению с часами счета роли не играют. Так что увеличение таблицы окупает себя! Увеличиваем таблицу еще.

t=Sort[DeleteCases[Array[f,1000,4],Null]];

Теперь снова измеряем время и процент отбраковки:

prog//Timing

Yes= 12827 ; No= 87173 {74.172 Second,Null}

Конечно, время возросло: теперь отбраковка занимает более одной минуты. Но про­ верять по-настоящему придется менее 13 % чисел! Потратив дополнительно на пред­ варительную отбраковку чуть более минуты, мы почти в 8 раз уменьшим количество чисел, которые надо будет проверить на простоту. Оказывается, увеличение таблицы оправдало себя и на этот раз. Поэтому еще раз увеличиваем таблицу.

t=Sort[DeleteCases[Array[f, 10000,4], Null]]///Timing {89.969 Second,'Null}

На этот раз сама генерация таблицы занимает более минуты. Но и это время оку­ пается.

prog//Timing

Yes= 9821 ; No= 90179 {527.781 Second,Null}

204

Гпава 7

Количество чисел, которые подлежат проверке на простоту, мы уменьшили более чем в 10 раз! Может быть, стоит увеличить таблицу еще? Ну, на этот раз, во всяком случае, уж точно не в 10 раз. Ведь ждать генерации таблицы несколько часов — заня­ тие весьма неприятное! Поэтому лучше сначала попытаться увеличить таблицу, ска­ жем, в 2,5 раза.

t=Sort[DeleteCases[Array[f,25000,4],Null]];//Timing {638.234 Second,Null}

На генерацию таблицы понадобилось более десяти с половиной минут. Заметно, что время генерации таблицы растет нелинейно с ростом таблицы. Наверное, пора ос­ тановиться. Проверяем:

prog//Timing

Yes= 8985 ; No= 91015

{1201.44 Second,Null}

На отбраковку потрачено более двадцати минут. Это на 9,39063 минут больше, чем в предыдущем случае. Несмотря на все это, дополнительно отбраковать удалось менее одного процента чисел. Действительно, пора остановиться. (Надо признать, конечно, что это довольно грубый, притом субъективный критерий останова.)

Теперь можем приступить к написанию программы. Пусть ее заголовок (имя и список параметров) будет М5п[п_]. (Параметр п — верхний конец диапазона, в кото­ ром происходит перебор.) Прежде всего нужно решить следующий вопрос: как испо­

льзовать таблицу так, чтобы не пропустить ни одного простого числа вида 5*2я+1 даже в том случае, если оно использовалось при построении таблицы. Конечно, мож­ но заблаговременно составить список п для таких чисел. Но это значит, что список п для таких чисел нужно составить по другой программе, решающей фактически ту же задачу, решением которой мы сейчас как раз и заняты. При решении этой вспомога­ тельной задачи, правда, можно ограничиться перебором на меньших начальных отрез­ ках натурального ряда. Конечно, этот способ решения вполне корректен. Но такой способ ведет к программе, которая выглядит несколько искусственно. Фактически по­ лучится не программа, а подглядывание в раздел ответов в случае небольших зна­ чений п. Ничего плохого в этом нет. Просто этот способ реализовать нельзя, если мы не знаем нужную нам часть ответа для решаемой задачи. Поэтому поступим иначе. Мы вообще сделаем вид, что не знаем даже начальных значений п. Сначала в прог­

рамме найдем все те начальные значения п, при которых числа вида 5-2л+1 не пре­ восходят тех простых чисел q, которые использовались для построения таблицы. Хотя

это и

будет чистый перебор, перебирать придется недолго, поскольку числа вида

5 2л +1

с ростом п растут гораздо быстрее, чем л-е простое число рп. Например, если

мы будем строить таблицу, используя первые 25 000 простых чисел, то перебирать без дополнительной отбраковки придется лишь первый десяток нечетных значений п. За­ тем можно будет построить таблицу и запустить процесс предварительной отбраковки. Поскольку проверка малых значений выполняется весьма быстро, о времени ее вы­ полнения можно не беспокоиться. Может возникнуть и противоположный вопрос: не слишком ли рано мы запускаем отбраковку? Но мы видели, что даже отбраковка 100 000 нечетных чисел с помощью огромной таблицы, построенной с использованием первых 25 000 простых чисел, занимает всего лишь чуть более 20 минут. Так что на предварительную отбраковку одного числа требуется совсем немного времени. Ну а выигрыш отбраковка дает уже при весьма малых значениях п. Поэтому не стоит бес­ покоиться из-за того, что на некотором, весьма малом отрезке отбраковка может оказаться менее эффективной, чем прямая проверка. Для таблицы, построенной с использованием первых 25 000 простых чисел, перерасход времени, например, не превысит 12 секунд. Так что беспокойство по поводу того, что временные затраты

Модулярная арифметика...

205

на отбраковку малых значений п могут быть больше, чем на прямую проверку, не­ обоснованно. Но перед запуском отбраковки нужно сгенерировать таблицу t. Даже если таблица строится с использованием первых 10 000 простых чисел, генерация таб­ лицы занимает более минуты. Хотя это время (и даже большее) окупается, поведение программы из-за временных затрат на построение таблицы может показаться не­ сколько неожиданным. Действительно, сначала очень быстро печатаются первые зна­ чения л, затем происходит весьма заметная (но необъяснимая, если не знать причину) задержка из-за генерации таблицы, после чего опять довольно быстро печатается сле­ дующая серия небольших значений п и лишь затем (на моем компьютере после п = 5947) опять происходит весьма существенная задержка. Она связана как с отсутст­ вием искомых значений п вплоть до п = 13165, так и с тем, что для части чисел (менее 10 %) приходится выполнять полную проверку, которая теперь уже требует весьма ощутимых временных затрат. В приведенной ниже программе учтено все, о чем говорилось выше.

М5п[п_]:=

Block[{j=l, nnprime=25000,primen=Prime[nnprime],nn=Min[n,primen] },

While [5*2Aj+l<= primen && j<=n, {If[MknPrime[5 ,j ] , P r i n t [ j ] ] , j+=2}]; I f [j<n,

t=Sort[DeleteCases[Array[f,nnprime,4],Null]];

While[j<=nn,

I f [ f x [j , t], I f [MknPrime[5,j ], Print[ j ]]];

j+=2]]]

Эта программа существенно быстрее первоначальной, и даже за короткую летнюю ночь она сможет на весьма слабеньком ПК (по нынешним представлениям, конечно) найти п = 26607.

Ну и, наконец, мораль этой истории. Хотя поначалу казалось, что предварительная отбраковка сильно усложнит программу, оказалось, что в системе Mathematica есть все необходимые функции для того, чтобы компактно записать построение необходи­ мой таблицы и тем самым ускорить выполнение программы более чем в 10 раз!

Случай Ас = 3

 

 

Критерий простоты чисел

вида 3-2" + 1

несколько сложнее. Впрочем, вся слож­

ность заключается в выборе

основания

я, для которого проверяется сравнение

^ * = - 1 (mod *-2" + 1). Если п не делится на 4, то в качестве а достаточно взять 5.

Иными словами, если п не делится на 4, то 3-2” + 1 будет простым тогда и только то­

гда, когда 5^* =-1 (mod к-2п+1).

Вот как можно реализовать этот способ. Сначала определим вспомогательную j

функцию:

' j

MMkn[k_,n_] :=

 

k*2An+l

 

Теперь запишем критерий простоты.

MM3nPrime01[п_]:=

Module!{t = 3*2An+l},

If[Mod[nr4]!=0,PowerMod[5,(t-1)/ 2 , t]==t-l,PrimeQ[t]]]

После этого можем написать программу поиска тех л, при которых число 3• 2я+1 является простым.

ММЗп[п_]:=

Block[{i=l}, While [i<=n, {If[MM3nPrime01[i],Print[i]], i++}]3;

206

Глава 7

Конечно, данная программа работает существенно быстрее, чем написанная нами ранее для той же цели. Но все же давайте попробуем применить предварительный отсев.

Конечно, метод перебора можно усовершенствовать, если найти арифметическую прогрессию таких показателей л, для которых число 3-2"+1 делится на некоторое простое число q. Тогда показатели п, являющиеся членами этой прогрессии, при пе­ реборе можно пропустить, если 3-2п+1 ф q. Иными словами, можно пропустить те

показатели п из этой прогрессии, для которых число 3-2"+1>0. Можно, конечно, пойти и дальше. Можно, например, взять несколько небольших нечетных простых чисел q, для которых существует решение сравнения 3-2'+1 = 0 (mod q). Пусть s —

показатель 2 по модулю q: 2'=1 (mod q). (Такое s можно найти с помощью функции M ultiplicati-veOrder [2, q].) Тогда при всех rmt (mod s) число 3- 2я +1 = 0 (mod q).

Если q не является числом вида 3*2я+ 1 , все числа rmt (mod s) можно опустить при

переборе. Если же q является числом вида 3*2" + 1, то нужно оставить только то число

л=/(mod s), при котором 3-2" + 1 = q. Иными словами, мы сможем отбраковать все те показатели п, которые принадлежат хотя бы одной из построенных нами указанным выше способом арифметических прогрессий. Как мы видели на примере чисел вида 5-2"+ 1, эта идея значительно повышает быстродействие программы и практически не усложняет перебор. Однако в данном случае таблицы получатся, конечно, совсем другие, и потому оптимальные параметры (главный из них — размер таблицы) при­ дется подбирать заново.

Сначала определим функцию, которая по заданному номеру п простого числа q находит / — решение сравнения 3- 2' +1 = 0 (mod q) и s — показатель 2 по модулю q, если такое решение существует.

f3[nj :=

Block[{q= Prime[n], s=MultiplicativeOrder[2,q],

t=MultiplicativeOrder[ 2 , q,{-PowerMod[3,-1,q]}]}, If[NumberQ[t],{s,t}]]

Теперь с помощью этой функции можно сгенерировать таблицу остатков и соот­ ветствующих им модулей.

t=Sort[DeleteCases[Array[f3,100,3],Null]]

{{3,1},{4,3},{10,7},{12,2},{18,14},{28,9},{36,28},{39,29},{48,5},{51,4 2}, {52,9}, {58,37}, {60,24},{66, 60},{73,10},{82,51},{96,6},{99,14},{100, 81}, {102, 60},{106, 89},{130, 123},{136/72},{138,28},{148,135},{156,22},{ 162,142},{166,98},{172,59},{178,159},{180,34},{183,48},{196,113},{204, 91}, {210,62},{22 6,67},{231,197},{243, 146},{268,25},{2 92,281},{316,225} ,{346, 21 }, {348, 148}, {372,320},{378,52},{388, 311},{418,109},{420,226}, { 442,345},{460, 119},{4 66,249},{490, 415},{508,245}, {522,486},{540, 166}, { 546,390}, {556, 399}}

Чтобы в таблице сначала шли меньшие модули s (они помогают отбраковать больше кандидатов), ,таблицу пришлось отсортировать, предварительно вычеркнув из нее эле­ менты N U L L , которые соответствуют тем простым модулям q, для которых сравнение

3-241 = 0 (mod q) неразрешимо.

Каждая пара в этой таблице описывает арифметическую прогрессию из показате­ лей, для которых выполняется некоторое тождественное сравнение. Пусть, например, в таблице имеется пара (s, /). Она была получена для некоторого простого числа q.

Ее наличие в таблице означает, что 3*2*,+| +1 s 0 (mod q) для любых натуральных к.

Модулярная арифметика...

207

Теперь давайте предварительно оценим количество чисел, которые может помочь от­ браковать такая таблица. Для этого сначала составим программу, которая подсчитывает количество чисел, отбракованных с помощью таблицы. Для этой программы нам пона­ добится функция, отбраковывающая заданное число по заданной таблице. Вот код этой функции.

fx[n_,tj :=

Block[{1= Length[t],i=l,c=(i<=l)}, While[c,

{s,r}=t[[i]]; answer^ (Mod [n,js] !=r); i++; c=(i<=l)&&answer];

answer]

Эта функция принимает значение True, если число п не было отсеяно с помощью таблицы t, т.е. если оно подлежит дальнейшему испытанию. Теперь можем составить программу подсчета отбракованных чисел.

,prog:= Block[{nn=100000,Yes=0,No=0, j=l},

While[j<=nn,

If[fx[j,t],Yes++,No++];

j++];

Print["Yes=",Yes,No=",No]]

Заметьте, что программа отличается от аналогичной программы для чисел вида 5 2* + 1. Отличие связано с тем, что на этот раз перебор не ограничивается только не­

четными показателями. Конечно, как и для чисел вида 5-2' + 1, количество отбрако­ ванных чисел и время отбраковки зависят от размера таблицы /. Сначала давайте сге­ нерируем совсем небольшую таблицу.

t=Sort[DeleteCases[Array[f3,10,3],Null]]

{{3,1},(4,3),{10,7),{12,2),(18,14),{28,9],{36,28)}

Теперь запускаем программу.

prog//Timing

Yes= 33650 ; No= 66350 {4.937 Second,Null}

Как видите, менее чем за 5 секунд удалось отбраковать более 66 % чисел! Давайте попробуем увеличить таблицу.

t=Sort[DeleteCases(Array[f3,25,3],Null]]

{{3,1},{4,3},{10,7},{12,2},{18,14},{28,9}, {36,28}, {39,29},{48,5},{51,4 2},{52,9},{58,37},{60,24},{66,60}, {82,51}, {100,81}}

Опять считаем процент отбракованных чисел и замеряем время.

prog//Timing

Yes= 25645 ; No= 74355 {7.188 Second,Null}

Время отбраковки возросло чуть более чем на 2 секунды, а отбраковали мы более 74 % чисел! Значит, таблицу стоит увеличить.

t*Sort[DeleteCases[Array[f3,100,3], Null]];

208

Гпава 7

Опять измеряем время.

p rog//T im in g

Yes= 21055 ; N0= 78945 {14.765 Second,Null}

Как видим, таблицу увеличивать стоило! Время отбраковки, конечно, возросло, притом более чем вдвое, но 15 секунд по сравнению с часами счета роли не играют. Так что увеличение таблицы окупает себя! Увеличиваем таблицу еще.

t=Sort[DeleteCases[Array[f3,1000,3],Null]];

Теперь снова измеряем время и процент отбраковки.

prog//Timing

Yes= 15388 ; N0= 84612 (84.64 Second,Null]

Конечно, время возросло: теперь отбраковка занимает более одной минуты. Но про­ верять по-настоящему придется чуть более 15 % чисел! Потратив дополнительно на предварительную отбраковку немногим более одной минуты, мы почти в 8 раз умень­ шим количество чисел, которые надо будет проверить на простоту. Так что увеличе­ ние таблицы оправдало себя и на этот раз. Поэтому еще раз увеличиваем таблицу.

t=Sort[DeleteCases[Array[f3, 10000,3],Null]];//Timing (90.703 Second,Null}

На этот раз сама генерация таблицы занимает более полутора минут. Но и это время окупается.

prog//Timing

Yes= 12123 ; N0= 87877

(630.734 Second,Null}

 

Количество чисел, которые подлежат проверке на простоту, мы уменьшили более

 

чем в 8 раз! Может быть, стоит увеличить таблицу еще? Ну, на этот раз, во всяком

 

случае, уж точно не в 10 раз. Ведь ждать генерации таблицы несколько часов — заня­

 

тие весьма неприятное! Поэтому лучше сначала попытаться увеличить таблицу, ска­

 

жем, в 2,5 раза.

 

t=Sort[DeleteCases[Array[f3,25000,3],Null]];//Timing

 

(629.781 Second,Null}

 

На генерацию таблицы понадобилось почти десять с половиной минут. Заметно,

I

что время генерации таблицы растет нелинейно с ростом таблицы. Наверное, нужно

t

проверить, не пора ли остановиться.

| prog//Timing

Yes= 11129 ; N0= 88871 (1405.2 Second,Null}

На отбраковку потрачено почти двадцать четыре минуты. Это существенно боль­ ше, чем в предыдущем случае. Несмотря на это, дополнительно отбраковать удалось чуть более одного процента чисел. Кажется, пора остановиться. (Конечно, это до­ вольно грубый и притом субъективный критерий останова.)

Теперь можем приступить к написанию программы. Пусть ее заголовок (имя и список параметров) будет М3п[п_]. (Параметр п — верхний конец диапазона, в кото­ ром происходит перебор.) Прежде всего нужно решить ^следующий вопрос: как ис­

пользовать таблицу так, чтобы не пропустить ни одного простого числа вида 3-2" + 1 даже в том случае, если оно использовалось при построении таблицы. Конечно, можно

Модулярная арифметика...

209

заблаговременно составить список л для таких чисел. Но это Значит, что список л для таких чисел нужно составить по другой программе, решающей фактически ту же зада­ чу, решением которой мы сейчас как раз и заняты. При решении этой вспомога­ тельной задачи, правда, можно ограничиться перебором на меньших начальных отрез­ ках натурального ряда. Конечно, этот способ решения вполне корректен. Но он ведет к программе, которая выглядит несколько искусственно. Фактически получится не программа, а подглядывание в раздел ответов в случае небольших значений п. Ничего плохого в этом нет. Просто этот способ реализовать нельзя, если мы не знаем нужную нам часть ответа для решаемой задачи. Поэтому мы поступим иначе. Мы вообще сде­ лаем вид, что не знаем даже начальных значений п. Сначала найдем в программе все

те начальные значения п, при которых числа вида 3-2" +1 не превосходят тех простых чисел q, которые использовались для построения таблицы. Хотя это и будет чистый перебор, перебирать придется недолго, поскольку числа вида 3• 2" +1 с ростом л рас­ тут гораздо быстрее, чем л-е простое число рп. Например, если мы будем строить

таблицу, используя первые 25 000 простых чисел, то перебирать без дополнительной отбраковки придется лишь три-четыре десятка значений л. Затем можно будет по­ строить таблицу и запустить процесс предварительной отбраковки. Поскольку провер­ ка малых значений выполняется весьма быстро, о времени ее выполнения можно не беспокоиться. Конечно, вполне законным является и вопрос: не слишком ли рано мы запускаем отбраковку? Но мы видели, что даже отбраковка 100 000 чисел с помощью] огромной таблицы, построенной с использованием первых 25 000 простых чисел, за- i нимает не более 24 минут. Так что на предварительную отбраковку одного числа тре­ буется совсем немного времени. Ну а выигрыш отбраковка дает уже при весьма малых значениях л. Поэтому не стоит беспокоиться из-за того, что на некотором, весьма ма­ лом отрезке отбраковка может оказаться менее эффективной, чем прямая проверка. Для таблицы, построенной с использованием первых 25 000 простых чисел, перерас­ ход времени, например, не превысит десятка-двух секунд. Так что беспокойство по поводу того, что временные затраты на отбраковку малых значений л могут быть больше, чем на прямую проверку, необоснованно. Но перед запуском отбраковки нужно сгенерировать таблицу t. Даже если таблица строится с использованием пер­ вых 10 000 простых чисел, генерация таблицы занимает немного более полторы минуты.

Кроме того, в данной программе есть еще один источник непредсказуемых вре­ менных затрат. Это показатели л, кратные 4. Давайте сначала посмотрим, что это бу­ дут за показатели. Вот нужная нам программа.

prog401:=

Block[{nn=8000,Yes=0,No=0,j=4},

While[j<=nn,

If[fx[j,t],Print[j,", "];Yes++,No++]; j+=4];

Print["Yes=",Yes, N o = " , N o ]]

В пределах первых восьми тысяч будет найдено 401 число, кратное 4, такое, что его не удастся отсеять с помощью таблицы, для генерации которой использовались 25 000 простых чисел. (Отсеяно будет 1599 чисел.) Далее приведен список тех чисел, которые кратны 4 и выдерживают отсев.

36, 92, 96, 108, 128, 132, 164, 188, 200,

216,

240, 252, 276, 296,

308,

312, 332,

336,

368,

372,

396,

404,

408,

416,

440,

452,

488,

492,

512,

516, 540,

560,

584,

600,

620, 636, 660, 668, 696, 704,

708,

732,

768,

800, 828,

836,

884,

888,

912, 920, 948, 956, 992, 996,

1016,

 

1020,

1028, 1052, 1056,

1088,

1092,

1124,

1128,

1136,

1140,

1176,

 

1196,

1212, 1232, 1280,

1304,

1332,

1376,

1388,

1412,

1424,

1428,

 

1476,

1484, 1496, 1548,

1568,

1596,

1608,

1664,

1668,

1680,

1716,

 

210

 

 

 

 

 

 

 

 

 

 

 

 

Гпава 7