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

8 6 6 5 7 8 8 9 9 5 3 5 0 2 7 2 2 5 3 9 5 4 3 4 5 1 6 5 9 8 3 9 1 7 3 3 6 4 2 7 5 0 7 1 5 4 3 3 1 7 4 9 3 8 6 3 7 7 9 5 7 6 5 0 2 2 3 3 0 7

1 6 8 9 5 8 6 3 7 1 9 7 1 9 2 1 1 0 5 7 8 7 3 7 8 5 7 3 3 6 9 4 3 2 1 2 4 5 7 7 1 5 5 2 1 2 7 5 5 1 3 9 9 8 3 1 7 7 8 5 4 7 6 7 1 6 7 8 5 9

1 2 9 9 6 4 5 0 6 7 2 9 6 2 7 4 8 3 7 3 6 5 3 0 2 2 1 5 2 3 4 3 2 0 5 0 7 4 7 8 3 4 0 9 2 7 9 0 5 6 5 3 7 1 2 7 3 8 3 2 6 4 0 5 3 5 9 0 9 7

6 9 9 6 3 5 1 3 4 3 5 9 7 7 5 3 7 9 9 2 8 3 6 8 0 7 5 2 8 1 7 5 4 8 3 8 2 7 2 4 4 7 8 1 4 4 5 3 6 9 4 0 9 7 9 9 7 2 3 0 4 7 1 8 4 1 7 6 2 5

8 9 4 4 7 9 5 1 5 4 0 1 8 0 7 2 6 2 4 2 8 3 6 5 9 7 6 1 4 2 9 1 8 8 3 4 8 9 6 7 9 1 8 8 1 5 3 7 7 2 8 5 4 7 6 7 8 1 0 7 4 9 6 6 1 6 1 2 6 6

1 8 5 4 7 6 2 6 6 6 8 5 3 2 3 5 5 2 9 0 0 5 5 7 1 8 8 8 4 9 1 6 7 9 8 8 5 5 4 7 0 0 6 8 4 7 3 5 8 2 6 8 5 0 8 9 7 3 9 1 8 7 0 0 8 5 1 0 7 5

4 0 2 8 1 8 8 5 3 9 2 5 3 4 9 0 5 2 9 1 2 2 8 8 2 0 3 9 7 1 9 7 2 4 0 3 2 2 3 5 7 8 7 0 0 6 0 7 3 2 8 3 8 7 7 3 5 8 2 8 2 6 1 7 0 0 4 3 1 5

0 6 0 2 2 5 0 4 0 6 6 0 1 9 6 1 6 5 6 9 9 4 3 9 7 5 4 3 6 1 0 2 6 8 5 5 2 6 6 3 7 4 0 3 6 6 8 2 9 0 6 1 9 0 1 7 4 9 2 3 4 9 4 3 2 4 1 7 8 7 9 9 3 5 9 6 8 1 4 2 2 6 2 7 1 7 7 2 8 9 } '

Вычисление, как видите, заняло более четырех минут (260,64 секунды). Но гораздо быстрее результат можно получить вот так.

P ow erM od[ 9 , 9 Л9 , 1 0 Л1 0 0 0 ] //T im in g

{ 0 . 0 1 5 S eco n d , 1 6 4 3 7 8 9 4 7 5 7 4 7 7 8 4 4 7 3 1 1 4 2 7 8 0 7 6 7 0 3 7 2 6 7 0 0 ...2 6 2 7 1 7 7 2 8 9 }

(Здесь середина пропущена.) Это вычисление заняло лишь 0,015 секунды. Это более чем в тысячу (точнее, в 17376) раз быстрее!

Пример 7.4. Графики функции PowerMod.

Давайте теперь построим несколько графиков функции PowerMod. Поскольку это функция двух аргументов, построим изображения поверхности z = PowerMod [х, у), Для этого используем функцию

P lo t3 D [ P o w e r M o d [I n t e g e r P a r t [x ]

, I n t e g e r P a r t [ у ] , 1 0 0 0 ] , {x , 1 , 1 0 } , {у , 1 , 1 0 } ] ;

10

А вот вид издалека.

Plot3D(PowerMod[IntegerPart[x]

, I n t e g e r P a r t [ y ] , 1 0 0 0 ] , ( x , 1 , 2 5 } , {y , 1 , 2 5 } ] ;

Китайская теорема об остатках —

функция ChineseRemainder

Хитрый китаец Сунь Цюуоколо 2000 лет назад (конечно, это могло быть и 200 лет до нашей эры, и 200 лет после начала нашей эры) открыл правило “гай-йен” (“боль­ шое обобщение”), которое является частным случаем целочисленного аналога интер­ поляционной формулы Лагранжа для полиномов. (Правда, примерно в то же самое вре­ мя, в 100 г. н. э., греческий математик Никомах знал точно тот же частный случай.) Полностью формула была сформулирована и доказана впервые, по-видимому, в 1734 году Л. Эйлером. Сформулируем китайскую (или греко-китайскую, как назвал ее Акритас, автор одного из широко известных учебников по компьютерной алгебре?) теорему об остатках и укажем соответствующее правило (формулу).

Пусть m,, /гг,, тГ — попарно взаимно простые модули (попарно взаимно простые положи­ тельные целые числа), т = т, пи ... тг — их произведение, а о, и,, и2, иг — произ­

вольные целые числа. Тогда существует ровно одно целое число и9 такое, что айи<а+т9 удов­

летворяющее всем сравнениям и = и{(то6т().Более того,

и=X w c^ M<(modtn),

где с, — произведение всех модулей, кроме т.%: с{ = т/mi9 а г/.=crl(modmf) (т.е. dc.= l(modm()).

Эта теорема имеет многочисленные применения в математике, информатике (машинная арифметика) и криптологии (секретные ключи). Для нахождения числа м, существование которого гарантировал хитрый китаец Сунь Цю, в пакете теории чисел существует “китайская” функция C h in e seR em a in d e r . C h in e seR em a in d e r [список^], список_2] — это такое наименьшее целое неотрицательное число и, что Mod [и, список_2] = список_ 7. Чтобы ею воспользоваться, нужно сначала загрузить пакет тео­ рии чисел.

1Встречается также написание Сунь Цу.

< < N u m b e rT h e o ry 'N u m b er T h eo ry F u n c tio n s'

После загрузки пакета можно вызвать эту функцию.

ChineseRemainder[{0,1/2}, {4,9,121}]

244

Полученный ответ означает, что 244s0(mod4), 244sl(mod9) и 244 s2(modl21). С помощью функции Mod это легко проверить.

Mod[244,{4,9,121}]

{ 0 , 1, 2 }

Пример 7.5. “Задача о корзинке с яйцами”. Есть одна очень старая задача, в которой используется китайская теорема об остатках. Вот современная формулировка задачи. Крестьянка (наверняка отличница местной школы, а может даже победительница олимпиад по поднятию тяжестей) несла по проселочной дороге на базар яйца. Однако по дороге очень быстро ехал экипаж с сеньором, и экипаж нечаянно зацепил корзин­ ку, в которой были яйца. Корзинка опрокинулась, и яйца разбились. Сеньор пожелал возместить убытки и спросил, сколько яиц было в корзинке. Крестьянка ответила, что не помнит, однако вспомнила, что когда вынимала яйца парами, то оставалось одно яйцо, когда тройками, то — два яйца, когда семерками — то 6 яиц и, вообще, когда она вынимала по р, ( р{ — /-е простое число) яиц, то в корзинке оставалось ровно /

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

C h in e s e R e m a in d e r [Range[100],Prime[Range[100]] ]

7088952796061978701604564987060183578666885084349794336112315062885608

9433674441002051192538320421552012038430556484751827147448338964951322

2800615000316327826942553040750737082239796728321817952808494505090787

967742393

Ну и корзинка, которая вмещала столько яиц!

Корни в системе остаточных классов

Задача о корзинке с яйцами представляет собой задачу о решении системы срав­ нений вида M3M((modm() с попарно взаимно простыми модулями щ . Конечно, она

допускает дальнейшие обобщения. Если не считать линейных сравнений и их систем, то наиболее простой из таких задач окажется задача об извлечении квадратного корня

в системе остаточных классов, т.е. задача о решении сравнения хг =(1 (modл ).

/

Квадратный корень по модулю —

функции SqrtMod и SqrtModList

Конечно, квадратных корней в системе остаточных классов — решений сравнения i дг =d (modл) — может быть более одного. Поэтому в системе Mathematica предуемот-!

рено две функции для нахождения таких решений. Функция SqrtM od находит наи­ меньший вычет, а функция SqrtM odL ist — список всех вычетов.

Наименьший квадратный корень по модулю — функция SqrtMod

Выражение SqrtMod[d, п] представляет собой наименьший неотрицательный квадратный корень из d по модулю я, т.е. такой наименьший неотрицательный вычет т,

что т2 =d(modn). Но это, конечно, в том случае, если такое т существует. Для этого d, понятно, должно быть полным квадратом по модулю я, т.е. символ Якоби jacobiSymbol [d, п] должен быть равным 1. Это условие также достаточно, если я является простым числом. Для простых я система Mathematica использует алгоритм Шенкса (Shanks). Для примера найдем самое маленькое неотрицательное целоечисло, такое, что его квадрат равен 3 по модулю 11.

SqrtMod [3,11]

5

Этот результат легко проверить.

Mod[Range [1 1 ]А2 ,1 1 ] 11,4/9,5,3/3/5, 9, 4,1,0}

Как видите, только квадраты чисел 5 и 6 по модулю 11 равны 3.

Если квадратный корень из d по модулю я не существует, SqrtMod [d, п] останет­ ся невычисленным.

SqrtMod[3, 5]

PowerMod::root The equation xA2

3 (mod 5) has no integer solutions. More...

SqrtMod[3, 5]

Убедимся, например, что квадратный корень из 3 по модулю 5 не существует.

Mod[{0,1,2,3,4}А2,5]

10,1,4,4,1} Ну и, конечно, нельзя не пройти мимо вычисления квадратного корня из 2.

SqrtMod[2,10А64+57] //Timing (0.016

Second,876504467496681643735926111996546100401033611976777074909122865

}

Функция SqrtMod может вычислять квадратные корни не только по простому мо­ дулю, но и по составному.

SqrtMod[3,11А512] //Timing {0. Second,

2718318953660034486015737562192386420778710413038788143739003923841504

0704483070407436551431775335493107616260852428027933235604314335079592

1736056000587594734102740265815589321182911694313197466445321220556586

2095453218619638905081162811397143598915650966649758986637069077078836

4312646366381249473813239422043085278619070206444242341251162079379755

9465277881216754388848763002962044844839330650726351213084503486673574

1383763679625506479581066834653760431240067274371648124587492530060166

7779537400688115698515767318620929348306180}

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

используется тогда, когда нужно найти наименьшее т , удовлетворяющее хотя бы од­ ному из сравнений кт (mod/i) . 3 является наименьшим натуральным числом,

удовлетворяющим хотя бы одному из сравнений 2w=5(mod/z) и 2m=8(modn).

Table[PowerMod[2,m,11],{m,10}]

{2, 4/ 8 , 5 ,1 0 , 9, 7 , 3 , 6 ,1 }

Вот как это можно вычислить.

M u ltip lic a tiv e O r d e r [ 2 , 1 1 , { 5 , 8 } ]

3

Пример 7.6. Длина периода систематической дроби по основанию Ь. Вот функция, которая вычисляет длину периода систематической дроби, представляющей рацио­ нальное число г по основанию Ь.

DigitCycleLength [r_Rational, b__Integer?Positive] := MultiplicativeOrder[b,FixedPoint[Quotient[#,GCD[#,b]]&,Denominator[r]]

]

Вот пример. Длина периода дроби 1/49 в десятичной системе равна 42.

DigitCycleLength[ 1 /4 9 ,1 0 ]

42

А вот и подтверждение.

N [1/49,151]

0.020408163265306122448979591836734693877551

0 2 0 4 0 8 1 6 3 2 6 5 3 0 6 1 2 2 4 4 8 9 7 9 5 9 1 8 3 6 7 3 4 6 9 3 8 7 7 5 5 1

020408163265306122448979591836734693877551

0 2 0 4 0 8 1 6 3 2 6 5 3 0 6 1 2 2 4 4 8 9 7 9 5 9

Что такое первообразный корень по модулю п

Решая сравнение ктs i (modп) относительно к, можно заметить, что при некоторых

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

не существует. Первообразные корни существуют только по модулям 2, 4, р' и 2 /?'

(здесь р — произвольное нечетное простое число). Для вычисления первообразных кор­

ней по модулю п в системе Mathematica предусмотрена функция PrimitiveRoot [п].

P rim itiveR oot [ 5 ]

2

Если первообразного корня не существует, функция остается невычисленной.

,P rim itiveR oot [1 1 * 1 3 ] P rim itiv eR o o t [1 4 3 ]

Следует помнить, что PrimitiveRoot использует Factorlnteger как подпро­ грамму, так что она может не справиться с вычислениями в случае очень большого параметра.

Критерии простоты чисел специального вида

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

П р о с т ы е ч и с л а М е р с е н н а , т е с т Л к ж а - Л е м е р а

Наибольшее простое число поймано решетом Эратосфена.

Заголовок статьи из газеты Suddeutsche Zeitung от 17 ноября 1978 года, рассказывающей о нахождении калифорнийскими школьниками, Никкелом и Ноллом,

6533-значного простого числа Мерсенна 221701 —1

Конечно, решето Эратосфена годится для поиска простых чисел Мерсенна не бо­ лее чем консервная банка для плавания через Атлантический океан. Гораздо более пригодна для этой цели функция PrimeQ. С ее помощью мы нашли все те 1индексы простых чисел Мерсенна, которые не превышают 5000. Дальнейшее продвижение оказалось практически невозможным. Однако есть другой метод проверки простоты чисел этого вида. Его опубликовал Э. Люка (Lukas2) в 1878 году, а в 1930 году усовер­ шенствовал Д. X. Лемер. Он основан на следующей теореме.

Если р>2 — простое число, то число Мерсенна М р = 2Р-1 является простым тогда и только тогда, когда L p_2 = 0, где последовательность L. определяется так:

1+ = 4 , Ll41=(L *-2)m od(2'-l).

Совершенно элементарное доказательство этой теоремы имеется во многих книгах; в частности, его можно найти во втором томе трехтомника Кнута. Чтобы сократить доказательство, можно также заметить, что критерий Люка основан на малой теореме

Ферма для чисел из полей Q(>/5 ) или Q(>/3). Тест Люка позволил проверить вручную все простые числа до Л/127. Когда для проверки простоты стали применяться ЭВМ,

началась настоящая гонка. Очередные новые простые числа Мерсенна рассматрива­ лись как свидетельство быстрого роста возможностей ЭВМ. Чтобы убедиться, что ? М тх является составным, Д. Дж. Уилер потратил около 100 часов машинного времени

1LLIAC I — самой быстрой ЭВМ выпуска пятидесятых годов. Гурвицу на IBM 7090 для этой же цели потребовалось лишь 5 часов в 1962 году. В 1964 году Гиллис сделал это на ILLIAC II за 49 минут. А в 1971 году Такерман потратил на это только 3 минуты на IBM 360/91. (Но то были суперЭВМ, а я на своем весьма слабеньком ПК на это потратил 4,125 секунды.) Проверка простоты Л/|1213 на ILLIAC II заняла 2 часа 15 ми­

нут, и всего лишь 8 минут на IBM 360/91. (Но это столько нужно суперЭВМ, а моему старенькому ПК нужно всего лишь 9,672 секунды.) Найденное очередное простое число Мерсенна становилось визиткой фирмы. Например, доказательству простоты числа

2 Встречается также написание Лукаш и Лукас.

д/11213 Иллинойский университет посвятил выпуск фирменного конверта. IBM в долгу

не осталась: она выпустила конверт, посвященный доказательству простоты М]9921

(это 6002-значное число, найденное Б. Такерманом 4 марта 1971 года, довольно долго было рекордсменом). Проверка простоты этого числа и у меня заняла изрядное время — 42,094 секунды.

С помощью функций Mod и PowerMod несложно запрограммировать тест Люка- Лемера. Сначала дадим нужные определения.

Мп[п_] :в2 Ап-1 ;

LLn [n_,mn__] : =Module [ {t = 4}, Do[t = Mod [PowerMod [t, 2,mn] -2, mn], {i/

n}]; t]; MPrime[p_]:=LLn[p-2,Mn[p]]==0;

Поскольку единственным простым числом Мерсенна Мр с четным индексом р

является М2, можем ограничиться поиском только нечетных индексов р. (Конечно,

программа при этом пропустит число М2 = 3, но разве мы можем о нем забыть?) Вот

нужная нам программа поиска нечетных индексов р простых чисел Мерсенна Мр.

Block[{p=3}, While [р<45000/ {If[MPrime[р],Print[р]], p=NextPrime [р] }]]

Даже на весьма слабеньком компьютере всего за несколько часов вы можете найти все нечетные индексы р первых 24 простых чисел Мерсенна, известных до 1980 года:

3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 42 53, 442 3, 9689, 9941, 1 12 1 3, 19937. Здесь, конечно, только 23 числа, поскольку индекс первого простого числа Мерсенна М2 = 3 четный: р = 2.

Надеюсь, вы о нем не забыли, так что не забудьте и отпечатать специальные кон­ верты!

Заметьте, что данную программу очень легко распараллелить и организовать поиск на нескольких ЭВМ. Кроме того, данную программу легко модифицировать так, что­ бы поиск продолжался после уже проверенных чисел. Впрочем, возможно, программу придется модифицировать так, чтобы она время от времени подавала “признаки жиз­ ни”, даже если не успевает найти простое число Мерсенна. Это можно сделать, на­ пример, так.

Block[{р=19937, i=l},

While [р<45000, {If[MPrime[р],Print[р]], p=NextPrime[р], 'If[Mod[i++,500]==0, -Print["***p=",p]]}]]

Так что если не поленитесь потратить ночку-другую на вычисления, то найдете и другие простые числа р, при которых р-е число Мерсенна Мр = 2р-1 является про­

стым. В частности, вы сможете найти 6533-значное простое число Мерсенна А/21701

(найдено калифорнийскими школьниками Никкелом и Ноллом в 1978 году), 6987значное простое число Мерсенна М23209 (найдено Ноллом), 13395-значное простое

число Мерсенна МШ91 (это 27-е простое число Мерсенна найдено Д. Словинским в

1979 году), 25962-значное простое число Мерсенна МШАЪ (это число Мерсенна най­

дено Дэвидом Словинским в январе 1983 года). Если хотите, проверьте, не пропустил ли я чего-нибудь и сможете ли вы продвинуться дальше. У меня проверка простоты Мгт заняла 53,547 секунды, М23209 — 64,25 секунды, МШ91 — 350,562 секунды (ого,

почти 6 минут!), М86243 — 1845,34 секунды (это чуть больше получаса). Думаю, что

дальнейшее продвижение по данной программе весьма затруднительно. В данном случае

либо нужен компьютер со значительно более высоким быстродействием, либо нужно изменить программу так, чтобы большинство составных чисел отсеивалось еще до применения критерия Люка. Казалось бы, перед применением критерия Люка можно

проверить, что число Мр не удовлетворяет сравнению аМр =а(тодМр) при некото­

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

программы. Время проверки сравнения аМр =а(то6Мр) может оказаться, например,

примерно равным времени проверки по критерию Люка-Лемера. Тогда, понятно, та­ кая дополнительная проверка только замедляет выполнение программы. Так что

XX столетия Дэвид Словинский. Вместе с Полем Гейджем он получил следующие ин­ дексы простых чисел Мерсенна: 756 839, 859 433, 1 257 787. Но степени 1 398 269 и 2 976 221 были найдены Жоэлем Арменгодом (Joel Armengaud) и Гордоном Спенсом (Gordon Spence) на персональных компьютерах по программе, разработанной Георгом Вольтманом (George Woltman) в 1996 году. Предварительная проверка Л/2976221 на ком­

пьютере с процессором Pentium с тактовой частотой 100 МГц заняла 15 суток в августе 1998 года! 27 января 1998 года Роланд Кларксон (Roland Clarkson) обнаружил простоту Мтхугп а 1 июня 1999 года Найян Хьяратвала (Nayan Hayratwala) — простоту Мб972593 .

Впрочем, люди стремятся к рекордам, и после 1983 года числа проверялись не подряд, так что могли что-то и пропустить! Правда, добровольцы тщательно потом все прос­ мотрели до миллиона!

Как из этого видно, критерий Люка-Лемера применяется только после некоторого предварительного отбора, которому уделяется огромное внимание. Давайте попробуем найти хотя бы какой-нибудь очень простой (а, значит, и довольно грубый) метод от­

оно делится на простое число q, если р является показателем 2 по простому модулю q. Так что такие индексы р чисел Мерсенна можно отсеять сразу, если М р ф q. Поэтому

для предварительного отсева нам понадобится таблица показателей 2 по простым мо- \ дулям q. Причем, чтобы уменьшить таблицу, из нее можно вычеркнуть все составные |

(показатель 2 по найденному простому модулю q, если он существует). Если оказыва­ ется, что показатель 2 по простому модулю q существует, проверяется, является ли он простым числом. Если оказывается, что найденный показатель 2 по простому модулю q является простым числохг, он возвращается в качестве значения функции. (В против­ ном случае функция возвращает Null.) С помощью этой функции легко строится нужная нам таблица простых показателей 2 по простым модулям q.

[[Sort I DeleteCases [Array [f, 100000,4] ,N ull] ] ] ;