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

Если показатель п равен одному из этих чисел, то наименьший делитель числа

3*2" +1, больший единицы, больше наибольшего простого модуля, использованного для составления таблицы. Вот примеры.

FactorInteger [3*2‘Л36+1] {{206158430209, 1} } FactorInteger [3*2л92+1]

{{1132314641089,1},{13119392730926401,1}} FactorInteger [3*2л96+1] {{392840481939253,1},{605040718740253,1}} Factorlnteger [3*2л108+1]

{{805213, 1}, {1781311693519, 1},{678750386188027, 1}}

Таким образом, нетривиальные делители чисел 3*2я+1, показатели которых вы­

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

J+ -4J;

P rin t[ " Y e s * " ,Y e s ,N o = " ,N o ] )

Вот результаты счета.

prog4//Timing

Yes14451 ; No60549 {1893.61 Second,Null)

Как видите, в пределах первых 300 000 придется проверить только 14451 число, кратное 4. Правда, отсев остальных чисел при этом у меня занял примерно 1893 се­ кунды, т.е. чуть больше получаса.

Хотя, как и в случае чисел 5*2я + 1, время на построение таблицы окупается, пове­ дение программы может показаться несколько неожиданным из-за временных затрат на построение таблицы и возможных временных затрат на вычисление функции PrimeQ в случае показателей, кратных 4. Действительно, сначала очень быстро печа­ таются первые значения л, а затем после п = 534 (на моем компьютере) происходит весьма существенная задержка. Она связана как с отсутствием искомых значений п вплоть до л = 2 208, так и с тем, что для части чисел (чуть более 11 %) приходится вы­ полнять полную проверку, которая теперь уже требует весьма ощутимых временных за­ трат, причем зачастую весьма непредсказуемых для показателей л, кратных 4. Затем вы­ дача замедляется, но поведение программы выглядит вполне стабильно, пока не будет

j

t-Sort[DeleteCases[Array[f3,nnprime,3],Null]]; While[j<-nn.

If[fx[j,t]. If[MM3nPrime01[j],Print[j]]];

Эта программа существенно быстрее первоначальной, но если хотите по ней найти л —20 909, понадобится время и ангельское терпение.

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

такое основание а, чтобы для установления простоты числа 3-2я+ 1 достаточно было

Но все дело в том, что если показатель л делится на 4, то выбрать основание а не­ сколько сложнее. Если показатель л делится на 4, в качестве основания а можно взять

любой квадратичный невычет г по модулю 3*2"+1. Поэтому если п делится на 4,то сначала нужно найти (желательно небольшой) квадратичный невычет г по модулю

З'Г +1, а затем проверить, выполняется ли сравнение г2-"4 =-1 (mod *-2"+1).

Но как найти квадратичный вычет? Для этого можно использовать квадратичный за­

кон взаимности.

Е ст Р>\ и 0>1 — положительные нечетные взаимно простые числа, то

Прежде всего заметим, что если Р = 3 2*+1, то •£—i = 3-2*~* — число четное

2

прпиЖ

Я - 1 (? -1

Поэтому (-1) 2 2 = 1 и потому (*)- . Как видим, для поиска квадратичного

невычета г можно использовать функцию JacobiSymbol. В качестве г можно брать

небольшие нечетные числа и вычислять для них символ Якоби. (Для квадратичного невычета символ Якоби равен -1.) Как только мы найдем квадратичный невычет г,

можем проверить, выполняется ли сравнение r ~ 'k =-1 (mod к - 2я+1). Теперь легко написать функцию, вычисляющую квадратичный невычет.

ММЗпг [пп_]

Module[{г=5,

JS=JacobiSymbol[5,nn]},

While[JS!=-l && r<nn, r+=2;

If[GCD[r,nn]~l,

JS=JacobiSymbol[r,nn]]];r]

Чтобы запрограммировать критерий простоты, воспользуемся функцией PowerMod.

MM3nPrime [n_] : =

Module!{t = 3*2Лп+1),

If[Mod[n,4]!=0,

PowerMod[5, (t—1)/2,t]==t-l,

PowerMod[ММЗпг[t],(t-1)/2,t]==t-l]]

Теперь, наконец, можем написать окончательную версию программы.

M3nv01 [ п_] : =

Block[{j=l, nnprime=25000/primen=Prime[nnprime],nn=Min[n,primen]}, While [MMkn[3,j]<= primen &&• j<=n, {If[MM3nPrime[j], Print[j]],

j + + ) l ;

i f [j <n; t=Sort[DeleteCases[Array[f3,nnprime,3],Null]]; While[j<=nn,

If[fx[j,t], If[MM3nPrime[j],Print[j]]];

j++]]]

Запуская программу с параметром nnprime = 10000 и 25000, убеждаемся, что

при этих (и всех промежуточных) значениях генерация таблицы начинается после вы­ вода показателя 12, но до вывода показателя 18. Затем программа довольно быстро добирается до показателя 3912, после вывода которого “уходит в себя” аж до вывода показателя 20 909. (Заметьте, что, в отличие от предыдущей версии программы, пе­ рерыв между показателями 534 и 2 208 не слишком заметен — эту “достоприме­

чательность” вы запросто можете пропустить, если выйдете выпить чашечку кофе!)

превосходящих 300 000, число 3-2Я-М вляется составным.

Резюме

Мы рассмотрели определение частного (функция Quotient) и остатка (функция Mod), а также возведение в степень в модулярной арифметике (функция PowerMod) и

основу основ модулярной арифметики — китайскую теорему об остатках (функция ChineseRemainder). Кроме того, мы затронули вопрос об извлечении квадратных

корней (функции SqrtMod, SqrtModList и JacobiSymbol). Но мы не ограничились корнями второй степени, а научились находить первообразные корни по модулю п (функция PrimitiveRoot), а также показатели (функция MultiplicativeOrder).

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

Задачи

Задача 7.1. Найти длину периодов дробей 1/19, 1/41, 1/(13*37), 1/(17*29) 1/(7*23*31), 1/(11*13*17), 1/(2*11*13), 1/(4*53*73).

Задача 7.2. Вычислить квадратные корни из 7 по модулю 17.

Задача 7.3. Решить квадратное уравнение 2дг2+5х+4 = 0 в поле вычетов по модулю 11.

Задача 7.4. Найти последние 50 цифр чисел 123402, З3’ , 44* 53’ (в случае, если1 число записывается менее чем 50 цифрами, вычислить число).

Это утверждение называется теоремой Эйлера и очень часто малой теоремой Фер­ ма, который знал ее частный случай для простых модулей т. Эта теорема позволяет упростить вычисление остатков степеней чисел, взаимно простых с модулем. Рассмот­ рим пример.

Пример 8.1. “Бытие Бога”. Во второй половине XIX века некоторых авторов при­

влекали числа 22, З3’ , 44* Вот что пишет об этих числах известный немецкий попу­ ляризатор математики Вальтер Литцман в середине XX столетия:

В одной, изданной в 1874 году книге под названием “Бьггие Бога” (автор книги — Крениг)

рассматривается ряд чисел: 22 , З3 , 44 О последнем из этих чисел автор говорит: “Представьте себе отрезок такой длины, что световому лучу понадобился бы квинтилли­ он лет, чтобы пройти этот путь. Затем представьте себе шар с диаметром, равным этому отрезку, наполненный типографской краской. Всей этой краски не хватило бы, чтобы четко напечатать это число даже самыми мелкими цифрами, какие только существуют” X. Маурер исследовал эти числа с точки зрения теории чисел. Обозначив для простоты

число х г через 3Jt и введя затем общее обозначение яа = х( х) (где х и п — целые по­

ложительные числа), он показал, как можно найти последние цифры этих числовых ве­

Что вы чувствуете по мере чтения этого отрывка? Я чувствую интерес и нарастаю­ щее недоверие. Я, например, в захвате от сравнения с шаром. Но я сомневаюсь, за­

канчивается ли 109 — на 9 392 745 289. И даже если последние п цифр числа ”9 по­

вторяются во всех последующих числах я+,9, я+29, и так далее, то я полагаю, что это как-то связано с девяткой и с тем, что эти числа записываются в десятичной сис­ теме счисления. Едва ли такими свойствами обладает и любой другой .ряд 1х = х, 2д ,

V и так далее, где х — целое. И действительно, * 2 = 2 заканчивается на 2, а

2а* = 4 — на 4. Уже последняя цифра не повторяется во всех последующих числах

лое, то тут явное вранье. Вальтер Литцман, правда, говорит не о таких свойствах, а об аналогичных свойствах, но в чем отличие, не указывает явно. Такими свойствами не обладают, например, числа а = 3, 4, 8 и многие другие. Поэтому в чем именно состо­ ит аналогия, для меня пока не совсем понятно. И мне, естественно, хочется прояс­

нить этот вопрос и еще сильнее хочется вычислить последние п цифр числа "9.

Давайте начнем с вычисления последних п цифр числа я9 . Как это сделать с помо­ щью системы Mathematica?

Вот определение функции ях .

GodF[n_,x_J :=Module[{t = 1), Do[t = x^t, {i, n}]; t]

Функция яа возрастает очень быстро, значительно быстрее факториалов, поэтому, если нужно вычислить, скажем, последние к знаков какого-нибудь значения этой функции, необходимо использовать PowerMod.

GodFLastk[n_,x_, kj := If [n==l,Mod [x, 10лк],PowerMod [x,GodF[n-1, x],IC'k]]

Пусть к = 50. Тогда мы в мгновение ока сможем вычислить последние 50 знаков чисел ‘9, 29 и 39 . Но вот вычислить последние даже 10 знаков числа 49 так легко не удастся. Этот пример показывает, что в некоторых случаях применение функции powerMod может продвинуть вычисления всего лишь на один шаг!

Но функция Эйлера (и малая теорема Ферма) позволяет быстро вычислить по­ следние 50 знаков числа 49. Программа, конечно, будет другая.

PowerMod[9,PowerMod[9,9Л9,EulerPhi[10Л50]],10А50]

30363975097419408973491530163140828233401045865289

Можно ли быстро вычислить последние 50 знаков числа 59 по такой программе?

PowerMod[9,PowerMod[9,9Л9Л9,EulerPhi[10л50]],10л50]

Едва ли.

В чем же причина такого медленного продвижения, почему мы продвигаемся только на один шаг? Это происходит потому, что мы на самом деле устраняем быст­

рый рост функции пх =хс 'х) только на последнем шаге. Почему это происходит? Потому, что мы все время связаны с определением функции GodF. Именно она ис­ пользуется внутри функции GodFLastk. Даже фактически отказавшись от явного применения функции GodFLastk, мы все равно стремимся записать выражение для

пх= хс 'х) с некоторыми упрощениями для последнего или предпоследнего возведения

в степень. Конечно, это так естественно использовать выражение для ях = хг '*)! Однако нужно иметь в виду, что использование подобных “прямолинейных” опре­ делений часто весьма неэффективно. Если уж мы решили повышать эффективность, то делать это нужно на всю “глубину” выражения. Правда, это неизбежно приведет к появлению новых функций. Нам, например, придется отказаться от функции GodFLastk. Как ни удивительно, но вместо нее нам понадобится только одна функ­ ция, притом довольно простая. Давайте найдем ее определение. Обозначим через <7(л,

х, к) остаток от деления пх на к. (Таким образом, G(n, х, К Г ) — это последние т

| Предположим теперь, что х и к взаимно просты. Тогда по малой теореме Ферма это выражение можно записать так: PowerMod [х, G(/i, х, <р(&)), к ]. Так что в этом случае С(л+1, х, к) = PowerMod[х, <7(и, х, <р(&)), к]. Этр и есть нужное нам рекуррентное соотношение. С учетом этого соотношения определение нужной нам функции выглядит так.

GodFk[n_,x_, k j :=

If[n==l,Mod[х , k],

If[GCD[x,k]==1,PowerMod[x,GodFk[n-1,x, EulerPhi[k]],k],

PowerMod[x,GodF[n-l,x],k]]

];

Теперь легко написать программу.

Do[Print ["n=", n, ":”,GodFk [n, x=9,10Ak] ],{n, k=50} ]

Вот полученные результаты (последние п цифр я выделил полужирным и кое-где Добавил для наглядности пробелы).

ДшшЙчч;TKTOKjpb ш ш т ййшшпярвшл®мдчгасжмма f»/»» члсш 3. Bor так эаяпияст*

«й mpwpwwa.

Ш)Ц^ariwittt[",i№f*="’„un„ "^"'„©адакКв^х-З,, l®*kj I» 4m»fc=5®B 1

Вот что получится (кое-где для наглядности добавлены пробелы, а повторяющиеся цифры выделены полужирным).

Теперь понятно, что имел в виду Вальтер Литцман под аналогичными свойствами. Начиная с некоторого п = п0 последняя цифра числа не изменяется, причем начиная

е п = л0 + 1 повторяются уже две цифры, с я = л1)+ 2 - три цифры, с п — п()+ к — к

цифр и т. д. (Через н()(х) удобно обозначать наименьшее п0 с таким свойством.) Для

некоторых л*, например для х = 9, п0(х) = 1, а для некоторых х л5(х)>1. Если //„(*)>!,

то нельзя утверждать, что последние п цифр числа пх повторяются во всех последую­

щих числах "+,.v, л+2д:, и т.д. Если /;0(х)>1, то повторение последних п цифр нач­

нется не с числа пх , а несколько позднее! Да уж, разбираться в популярных книжках иногда приходится с системой Mathematical

Пример 8.2. График функции Эйлера.

Давайте теперь построим график функции Эйлера. Сначала используем функции Table и EulerPhi для построения таблицы tl (точнее, списка) значений функции Эйлера.

tlTable[EulerPhi[k],{k,1,п=10л3}];

Теперь можем использовать функцию ListPlot для построения графика.

ListPlot[tl,PlotRange->All]

100C

800

600

400

200^

200 400 600 800 1000

А вот график для n = 100000.

20000 40000 60000 80000 100000