(*********************************************************************** Mathematica-Compatible Notebook This notebook can be used on any computer system with Mathematica 3.0, MathReader 3.0, or any compatible application. The data for the notebook starts with the line of stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. ***********************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 53237, 1596]*) (*NotebookOutlinePosition[ 53902, 1620]*) (* CellTagsIndexPosition[ 53858, 1616]*) (*WindowFrame->Normal*) Notebook[{ Cell[BoxData[ FrameBox[\(PROGRAMMING\ \ TIPS\)]], "Text", TextAlignment->Center, FontFamily->"Times", FontSize->18, FontWeight->"Bold", FrameBoxOptions->{BoxFrame->0.5, BoxMargins->{{1, 1}, {2, 2}}}], Cell[CellGroupData[{ Cell["What Is a Prime Number?", "Title"], Cell[TextData[{ "The notion of primality is, in many respects, elusive. In this column we \ describe several ways in which ", StyleBox["Mathematica", FontSlant->"Italic"], " can help us answer the age-old question, \[OpenCurlyDoubleQuote]Is n \ prime?\[CloseCurlyDoubleQuote]" }], "Subtitle"], Cell["by Stan Wagon", "Subtitle"], Cell[TextData[{ "Although it is one of the most elementary concepts of mathematics, the \ notion of primality still holds many secrets. In this article we will delve \ into several algorithms that determine primality. It is important that such \ an algorithm be fast, and for that one must give up on the notion of absolute \ correctness. Still, there is an algorithm that is very fast and is perfect \ for all integers up to ", Cell[BoxData[ \(TraditionalForm\`10\^16\)]], ", and possibly much farther. This algorithm is used by the built-in ", StyleBox["PrimeQ", "Input"], " command and we will discuss it in detail here. Along the way we will \ learn about efficient approaches to linear recurrences, as well as ways to \ deal with the situation where absolute certainty regarding a number's \ primality is required. " }], "Text"], Cell[CellGroupData[{ Cell[TextData[{ CounterBox["Section"], ". Probable Primes" }], "Section"], Cell[TextData[{ "Prime numbers have many beautiful properties. One of the most famous is \ that, if p is an odd prime, then ", Cell[BoxData[ \(TraditionalForm\`2\^\(p - 1\) \[Congruent] \ 1\ \(\((mod\ p)\).\)\)]], "This is known as Fermat's little theorem and is in fact more general: if \ gcd(a, p) = 1, then ", Cell[BoxData[ \(TraditionalForm\`a\^\(p - 1\) \[Congruent] \ 1\ \(\((mod\ p)\).\)\)]], " The proof is simple and can be found in any number theory textbook. \ Because modular powers can be computed very quickly (one uses the base-2 \ digits of p \[Dash] 1 to break the problem down into steps; by reducing mod p \ at each step the numbers stay small), this is a fine initial test for \ primality." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[{ \(\(p = Prime[10000]; \)\), \(PowerMod[2, p - 1, p]\)}], "Input"], Cell[BoxData[ \(1\)], "Output"] }, Closed]], Cell[TextData[{ "A composite, odd n for which ", Cell[BoxData[ \(TraditionalForm\`\(2\^n\[Dash]1 \[Congruent] \ 1\ \((mod\ n)\)\ \)\)]], "is called a ", StyleBox["2-pseudoprime", FontSlant->"Italic"], ". The following shows that there are none among the first 100 odd integers \ beyond one million." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Select[Range[10\^6 + 1, 10\^6 + 200, 2], PowerMod[2, # - 1, #] == 1\ &] == Select[Range[10\^6 + 1, 10\^6 + 200, 2], PrimeQ]\)], "Input"], Cell[BoxData[ \(True\)], "Output"] }, Closed]], Cell["\<\ Of course, 2-pseudoprimes do exist. There are 23 of them under \ 10000.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Select[Range[3, 10000, 2], \n\t \(! PrimeQ[#]\)\ && \ PowerMod[2, \ # - 1, #] == 1\ &]\)], "Input"], Cell[BoxData[ \({341, 561, 645, 1105, 1387, 1729, 1905, 2047, 2465, 2701, 2821, 3277, 4033, 4369, 4371, 4681, 5461, 6601, 7957, 8321, 8481, 8911}\)], "Output"] }, Closed]], Cell[TextData[{ "Similarly one can define ", StyleBox["b-pseudoprimes", FontSlant->"Italic"], " to be odd composite n for which ", Cell[BoxData[ \(TraditionalForm\`b\^n \[Congruent] 1\ \((mod\ n)\)\)]], "." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Select[Range[3, 10000, 2], \n\t \(! PrimeQ[#]\)\ && \ PowerMod[13, \ # - 1, #] == 1\ &]\)], "Input"], Cell[BoxData[ \({21, 85, 105, 231, 357, 427, 561, 1099, 1785, 1891, 2465, 3605, 5149, 5185, 5565, 6601, 7107, 8841, 8911, 9577, 9637}\)], "Output"] }, Closed]], Cell[TextData[{ "Pseudoprime tests are very fast, and so they are perfect for an initial \ investigation into the primality question for a large integer n. But the \ failure rate is too great and so more sophisticated methods must be used for \ something as important as ", StyleBox["PrimeQ", "Input"], ". Note that 561 (= 3\[CenterDot]11\[CenterDot]17) is both a 2-psp and a \ 13-psp. In fact, 561 is a b-pseudoprime for every b relatively prime to 561. \ Such numbers are called Carmichael numbers (it is now known that there are \ infinitely many of them) show that primality cannot be verified by just \ trying lots of b-pseudoprime tests.\n\nThis sort of test has nothing to do \ with probability theory. If a number passes several pseudoprime tests then \ one might say that it is probably prime, but the procedure is totally \ deterministic. Sometimes randomness is introduced into the choice of b, but \ then different runs might yield different results. To ensure repeatability, \ the built-in prime tester makes no use of random numbers." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell[TextData[{ CounterBox["Section"], StyleBox[". Strong Pseudoprimes", "Section"] }], "Section"], Cell[TextData[{ " A very pretty and easy enhancement of the pseudoprime test comes by \ considering square roots. This leads to the notion of ", StyleBox["strong pseudoprime", FontSlant->"Italic"], ". Suppose ", Cell[BoxData[ \(TraditionalForm\`b\^\(n - 1\) \[Congruent] 1\ \((mod\ n)\)\)]], " where n is odd. Write ", Cell[BoxData[ \(TraditionalForm\`n\ - \ 1\)]], " as ", Cell[BoxData[ \(TraditionalForm\`\(2\^s\) m\)]], " and observe that ", Cell[BoxData[ \(TraditionalForm\`b\^\(n - 1\)\)]], " is the last term of the sequence ", Cell[BoxData[ \(TraditionalForm\`b\^m\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\^\(2 m\)\)]], ", . . . , ", Cell[BoxData[ FormBox[ SuperscriptBox["b", RowBox[{\(2\^s\), AdjustmentBox["m", BoxMargins->{{-0.25, 0.25}, {0, 0}}]}]], TraditionalForm]]], ", where each term of the sequence is obtained by squaring its predecessor. \ But if n is indeed prime then not only is the last term equal to 1 (modulo \ n), but the last term that is not 1 must be ", Cell[BoxData[ \(TraditionalForm\`\(-1\)\)]], ". This is because the only square roots of +1 modulo a prime p are \ \[PlusMinus]1 (proof: ", Cell[BoxData[ \(TraditionalForm \`x\^2\ - \ 1\ = \ \((x\ - \ 1)\) \((x\ + \ 1)\)\)]], ", and if the right side is 0, then p divides it, in which case p divides \ one of the two factors so x \[Congruent] \[PlusMinus]1). So if p is prime the \ b-sequence just given either consists entirely of 1s or looks like { . . . , \ \[Dash]1, +1, . . . , +1}. Call n a ", StyleBox["b-strong pseudoprime", FontSlant->"Italic"], " (b-spsp) if its b-sequence looks like this but n is composite. Here is \ code to find such numbers. First we isolate an odd-part routine that returns \ ", Cell[BoxData[ \(TraditionalForm\`{m, \ s}\)]], " as just defined." }], "Text"], Cell[BoxData[ \(OddPart[n_] := Module[{m = n, s = 0}, While[EvenQ[m], m /= 2; \(s++\)]; {m, s}]\)], "Input"], Cell[CellGroupData[{ Cell[BoxData[ \(OddPart\ /@\ {100, 101, 102}\)], "Input"], Cell[BoxData[ \({{25, 2}, {101, 0}, {51, 1}}\)], "Output"] }, Closed]], Cell[TextData[{ "Now we can look at some b-sequences. In the code that follows ", StyleBox["NestList", "Input"], " is used to iterate the mod-n squaring function the number of times \ specified by ", StyleBox["oddPart", FontFamily->"Courier", FontWeight->"Bold"], ". For legibility we replace ", Cell[BoxData[ \(TraditionalForm\`n\ - \ 1\)]], " by ", Cell[BoxData[ \(TraditionalForm\`\(-1\)\)]], "." }], "Text"], Cell[BoxData[ \(spspSequence[b_, n_] := Module[{m, s}, \n\t\t{m, s} = OddPart[n - 1]; \n\t\t NestList[Mod[#\^2, n]&, PowerMod[b, m, n], s] /. n - 1 \[Rule] \(-1\)] \)], "Input"], Cell["\<\ When we look at the 1000th or 10000th prime we see the proper \ behavior.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(spspSequence[2, 7919]\)], "Input"], Cell[BoxData[ \({1, 1}\)], "Output"] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ \(spspSequence[2, 104729]\)], "Input"], Cell[BoxData[ \({36639, \(-1\), 1, 1}\)], "Output"] }, Closed]], Cell["\<\ The smallest 2-pseudoprime, 341, is quickly unmasked by this test, \ since the 2-sequence has a 1 preceded by a 32.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(spspSequence[2, 341]\)], "Input"], Cell[BoxData[ \({32, 1, 1}\)], "Output"] }, Closed]], Cell[TextData[ "Here is complete code to recognize b-strong pseudoprimes. We treat n = 1 as \ a special case. We need only check that \[Dash]1 is a member of the \ b-sequence with its last entry dropped, for its presence anywhere guarantees \ that the sequence has the form typical of primes."], "Text"], Cell[BoxData[{ \(StrongPseudoprime[b_, n_] := Module[{bSeq, m, s}, \n\t\t{m, s} = OddPart[n - 1]; \n\t\t bSeq = NestList[PowerMod[#, 2, n]\ &, PowerMod[b, m, n], s]; \n\t\t Union[bSeq] == {1} \[Or] MemberQ[Drop[bSeq, \(-1\)], n - 1]] /; n > 1\n\), \(\(StrongPseudoprime[_, 1] = False; \)\)}], "Input"], Cell[TextData[{ "We can now select all the 2-spsps below 10,000, where the term ", StyleBox["2-spsp", FontSlant->"Italic"], " is used for an odd integer that passes the 2-strong pseudoprime test but \ is not in fact prime." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Select[Range[10000], \[InvisibleSpace]\(! \((PrimeQ[#])\)\) \[And] StrongPseudoPrime[2, #]\ &]\)], "Input"], Cell[BoxData[ \({2047, 3277, 4033, 4681, 8321}\)], "Output"] }, Closed]], Cell[TextData[ "We see that this test too is fooled by some rather small integers. But in \ this situation the combination of several strong pseudoprime tests leads to a \ very strong test. For example, the first integer that is a b-strong \ pseudoprime for b = 2, 3, 5, and 7 is 3,215,031,751 [PSW]\n\nDaniel \ Bleichenbacher [Bleichenbacher, 1996] has made a detailed investigation of \ integers that are b-spsps for several b. His largest example is \ 68528663395046912244223605902738356719751082784386681071, which is a b-spsp \ for all b \[LessEqual] 100. Here is the proof, where the psp is given in \ factored form."], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Table[ StrongPseudoPrime[i, 18215745452589259639*4337082250616490391*867416450123298079], {i, 101}]\)], "Input"], Cell[BoxData[ \({True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, False}\)], "Output"] }, Closed]], Cell[TextData[{ "In such a case we say that 101 is the first witness to the compositeness \ of n. Unlike the case for pseudoprimes, all composites have a first witness \ that is relatively small. Well, this is under the assumption of the extended \ Riemann hypothesis, but the evidence for the truth of this statement is good. \ More precisely, under ERH every composite integer n is not a b-strong \ pseudoprime for some b less than ", Cell[BoxData[ FormBox[ RowBox[{"2", RowBox[{"(", RowBox[{"log", " ", FormBox[\(\(n)\)\^2\), "TraditionalForm"]}]}]}], TraditionalForm]]], ". In short, ERH \[Implies] \[OpenCurlyDoubleQuote]primes are in \ \[ScriptCapitalP]\[CloseCurlyDoubleQuote], where \[ScriptCapitalP] is the \ class of problems solvable in polynomial time. See [Bach and Shallit, 1996].\n\ \n", StyleBox["Mathematica", FontSlant->"Italic"], "'s ", StyleBox["PrimeQ", "Input"], " function starts off (after some checking for small divisors) by using a \ 2-spsp test. If the input fails the test, it is proved composite. If it \ passes, then a 3-spsp test is used, and then one other test, to be described \ in section 4." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell["3. Perrin's Pretty Primality Test", "Section"], Cell[TextData[{ "As a warmup to the more important Lucas test, we discuss a simple and very \ pretty test. Consider the recurrence defined by the code that follows; it was \ in fact studied by Lucas, but it is now known as the Perrin recurrence, after \ R. Perrin who examined it in 1899. The code that follows serves as the \ mathematical definition, but also works as a recursive definition in ", StyleBox["Mathematica", FontSlant->"Italic"], ". In such cases it is essential to cache the values as they are computed, \ and this is done by the ", Cell[BoxData[ FormBox[ StyleBox[\(a\_n = a\_\(n - 2\) + a\_\(n - 3\)\), FontFamily->"Courier", FontWeight->"Bold"], TraditionalForm]]], " phrase. While subscripted variables do not work in all situations, they \ can be used here, as opposed to the more familiar ", StyleBox["a[n_]", "Input"], " setup." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[{ \(a\_0 = 3; a\_1 = 0; a\_2 = 2; \n a\_n_ := \(a\_n = a\_\(n - 2\) + a\_\(n - 3\)\)\), \({a\_0, \ a\_1, \ a\_2, \ a\_3, \ a\_4, \ a\_5}\)}], "Input"], Cell[BoxData[ \({3, 0, 2, 3, 2, 5}\)], "Output"] }, Closed]], Cell[TextData[{ "Note that to clear the definition one would use ", StyleBox["Clear[Subscript]", "Input"], ". In order to allow ", StyleBox["a", "Input"], " to work directly as a function too, it is convenient to add the following \ two lines." }], "Text"], Cell[BoxData[{ \(a[n_] := a\_n\), \(\(Attributes[a]\ = \ Listable; \)\)}], "Input"], Cell["Here are the first 20 values.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(a[Range[0, 20]]\)], "Input"], Cell[BoxData[ \({3, 0, 2, 3, 2, 5, 5, 7, 10, 12, 17, 22, 29, 39, 51, 68, 90, 119, 158, 209, 277}\)], "Output"] }, Closed]], Cell[TextData[{ "Start from the left and count 0, 1, 2, 3, . . . . You will notice that at \ the primes, and only there, the a-value is divisible by the index (with the \ exception of ", Cell[BoxData[ \(TraditionalForm\`\(a\_1)\)\)]], ". We can check it." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(pr\ = \ Prime[Range[20]]; \na[pr]/pr\)], "Input"], Cell[BoxData[ \({1, 1, 1, 1, 2, 3, 7, 11, 28, 120, 197, 892, 2479, 4148, 11687, 56010, 271913, 461529, 2270882, 6599404}\)], "Output"] }, Closed]], Cell[TextData[{ "For a proof, we will need some elementary facts about linear recurrences. \ If p(x) is the characteristic polynomial of the recurrence \[LongDash] p(x) = \ ", Cell[BoxData[ \(TraditionalForm\`x\^3\)]], "\[Dash] x \[Dash] 1 in this case \[LongDash] then ", Cell[BoxData[ \(TraditionalForm\`a\_n\)]], "is equal to ", Cell[BoxData[ \(TraditionalForm \`\(c\_1\) \[Alpha]\^n\ + \ \(c\_2\) \[Beta]\^n\ + \ \(c\_3\) \[Gamma]\^n\)]], ", where \[Alpha], \[Beta], and \[Gamma] are the three roots of the \ polynomial and the coefficients ", Cell[BoxData[ \(TraditionalForm\`c\_i\)]], " are determined by the initial conditions. The roots of a cubic are easy \ to find and we can then solve for the coefficients. We suppress the large \ symbolic output of these two steps." }], "Text"], Cell[BoxData[ \({\[Alpha], \[Beta], \[Gamma]} = x /. Solve[x\^3 - x - 1 == 0, x]; \n soln = Solve[ Thread[\n\t\t\t\ \ {3, 0, 2} == c\_1\ \[Alpha]\^{0, 1, 2} + c\_2\ \[Beta]\^{0, 1, 2} + c\_3\ \[Gamma]\^{0, 1, 2}], {c\_1, c\_2, c\_3}]; \)], "Input"], Cell["\<\ Checking the numerical values shows that that the coefficients are \ surely {1, 1, 1}.\ \>", "Text"], Cell[CellGroupData[{ Cell["soln //N //Chop", "Input"], Cell[BoxData[ \({{c\_1.` \[Rule] 1.00000000000000022`, c\_2.` \[Rule] 0.999999999999999644`, c\_3.` \[Rule] 1.00000000000000022`}}\)], "Output"] }, Closed]], Cell[TextData[{ "Well, it is easy to see that setting each ", Cell[BoxData[ \(TraditionalForm\`c\_i\)]], "= 1 does indeed work. Setting n = 0 yields ", Cell[BoxData[ \(TraditionalForm\`c\_1\)]], "+ ", Cell[BoxData[ \(TraditionalForm\`c\_2\)]], " + ", Cell[BoxData[ \(TraditionalForm\`c\_3\)]], " = 3, which is satisfied. Setting n = 1 yields ", Cell[BoxData[ \(TraditionalForm\`0\ = \ \(a\_1 = \ c\_1\)\)]], "\[Alpha] + ", Cell[BoxData[ \(TraditionalForm\`c\_2\)]], " \[Beta] + ", Cell[BoxData[ \(TraditionalForm\`c\_3\)]], "\[Gamma] = \[Alpha] + \[Beta] + \[Gamma], which holds because 0 is the \ coefficient of ", Cell[BoxData[ \(TraditionalForm\`x\^2\)]], " in the polynomial. And it is easy to see, by expanding ", Cell[BoxData[ FormBox[ RowBox[{"(", RowBox[{"\[Alpha]", " ", "+", " ", "\[Beta]", " ", "+", " ", FormBox[\(\(\[Gamma])\)\^2\), "TraditionalForm"]}]}], TraditionalForm]]], ", that ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^2\)]], "+ ", Cell[BoxData[ \(TraditionalForm\`\[Beta]\^2\)]], "+ ", Cell[BoxData[ FormBox[ RowBox[{ FormBox["", "TraditionalForm"], \(\[Gamma]\^2\)}], TraditionalForm]]], " = 2. So we have that ", Cell[BoxData[ \(TraditionalForm\`\(\(a\_n\ = \)\ \)\)]], Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^n\ + \ \[Beta]\^n\ + \ \[Gamma]\^n\)]], ". In fact, for n \[GreaterEqual] 10, it is true that ", Cell[BoxData[ \(TraditionalForm\`a\_n\)]], " is simply the nearest integer to ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^n\)]], ", where \[Alpha] is the real root, which can be simplified to the \ following form." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\[Alpha] = \(\((\@27 - \@23)\)\^\(1/3\) + \((\@27 + \@23)\)\^\(1/3\)\)\/\(2\^\(1/3\)\ \@3\); \n{a[Range[9, 20]]\ , Round[\[Alpha]\^Range[9, 20]]}\)], "Input"], Cell[BoxData[ \({{12, 17, 22, 29, 39, 51, 68, 90, 119, 158, 209, 277}, {13, 17, 22, 29, 39, 51, 68, 90, 119, 158, 209, 277}}\)], "Output"] }, Closed]], Cell[TextData[{ "Now, if p is prime, then ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^p\ + \ \[Beta]\^p\ + \ \[Gamma]\^p\)]], " and ", Cell[BoxData[ FormBox[ RowBox[{"(", RowBox[{"\[Alpha]", " ", "+", " ", "\[Beta]", " ", "+", " ", FormBox[\(\(\[Gamma])\)\^p\), "TraditionalForm"]}]}], TraditionalForm]]], " are in the subring of algebraic integers of the number field determined \ by the roots of p(x). Therefore the binomial theorem and elementary facts \ about divisibility of ", Cell[BoxData[ FormBox[ RowBox[{"(", "\[NegativeThinSpace]", GridBox[{ { StyleBox["p", FontSize->10]}, { StyleBox["i", FontSize->10]} }, RowSpacings->0.3], "\[NegativeThinSpace]", ")"}], TraditionalForm]], "Input", FontWeight->"Plain"], " by p tell us that ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^p\ + \ \[Beta]\^p\ + \ \[Gamma]\^p\)]], " \[Congruent] ", Cell[BoxData[ FormBox[ RowBox[{"(", RowBox[{"\[Alpha]", " ", "+", " ", "\[Beta]", " ", "+", " ", FormBox[\(\(\[Gamma])\)\^p\ \((mod\ p)\)\), "TraditionalForm"]}]}], TraditionalForm]]], " in this ring. But ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\ + \ \[Beta]\ + \ \[Gamma]\)]], " = 0, so ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^p\ + \ \[Beta]\^p\ + \ \[Gamma]\^p\)]], " has the form p\[Rho] where \[Rho] is in the aforementioned ring. But ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^p\ + \ \[Beta]\^p\ + \ \[Gamma]\^p\)]], " is an integer \[LongDash] it equals ", Cell[BoxData[ \(TraditionalForm\`a\_p\)]], "\[LongDash] and \[Rho] is an integer linear combination of elements of an \ integral basis (see [Bach and Shallit, 1996, \[Section]6.7]) so \[Rho] must \ in fact be an integer too. Therefore ", Cell[BoxData[ \(TraditionalForm\`a\_p\)]], "is an integer divisible by p.\n\nWhenever one has a simple formula that \ holds for primes, it is natural to ask if it holds for any composites. Call n \ a ", StyleBox["Perrin pseudoprime", FontSlant->"Italic"], " if n is composite and n divides ", Cell[BoxData[ \(TraditionalForm\`a\_n\)]], ". They exist, but they are not easy to find. Indeed, none was known until \ 1982, when Adams and Shanks found the first one. Here are the results of a \ quick search up to 100." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Select[Range[2, 100], \[InvisibleSpace]\(! PrimeQ[#]\) && Mod[a\_#, #] == 0\ &]\)], "Input"], Cell[BoxData[ \({}\)], "Output"] }, Closed]], Cell[TextData[{ "The recursive definition makes even a modest search up to 1000 difficult \ because of recursion limits. One could program the definition in a ", StyleBox["Do", "Input"], "-loop. But it is much more efficient to take an approach based on \ matrices. First we digress to discuss how to raise matrices to very large \ powers, and that requires a digression on the raising of integers to very \ large powers modulo n. Recall that modular powers for integers are very fast \ when ", StyleBox["PowerMod", "Input"], " is used." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(PowerMod[1007, 34000000000, 100001]\)], "Input"], Cell[BoxData[ \(23354\)], "Output"] }, Closed]], Cell[TextData[{ "The algorithm underlying this is simple. As an example, suppose the \ exponent is 26. Write 26 in base 2 and view the goal, ", Cell[BoxData[ \(TraditionalForm\`b\^n\)]], ", as being ", Cell[BoxData[ \(TraditionalForm\`b\^11010\)]], ". Now start with 1 and scan the base-2 digits from left to right. When a 0 \ is read, square the current value modulo m. When a 1 is read, square and then \ multiply by b, again reducing modulo m. The sequence will look like the \ following, where the exponent is in base 2: ", Cell[BoxData[ \(TraditionalForm\`b\^0\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\^1\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\^10\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\^11\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\^110\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\^1100\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\^1101\)]], ", ", Cell[BoxData[ \(TraditionalForm\`b\^11010\)]], ".\n\nThe fact that a modular reduction takes place after each digit is \ scanned means that the numbers never get large. If we wished to implement \ this ourselves, here is a fast way to do it. ", StyleBox["Fold", "Input"], " is used to scan the base-2 digits and do the right thing." }], "Text"], Cell[BoxData[ \(PowerModManual[b_, n_, m_] := Fold[Mod[#1\^2*If[#2 == 1, b, 1], m]\ &, \n\t\t1, IntegerDigits[n, 2]] \)], "Input"], Cell[CellGroupData[{ Cell[BoxData[ \({PowerModManual[123, 456, 789], PowerMod[123, 456, 789]}\)], "Input"], Cell[BoxData[ \({699, 699}\)], "Output"] }, Closed]], Cell[TextData[{ "Now it is simple to generalize the idea to modular powers of integer \ matrices. We start with ", StyleBox["mat", "Input"], ", as opposed to the identity matrix." }], "Text"], Cell[BoxData[ \(MatrixPowerMod[mat_, \ n_, \ m_]\ := \ Fold[\n\ \ Mod[If[#2\ == \ 1, #1.#1.mat, #1.#1], m]\ &, \n\ \ mat, Rest[IntegerDigits[n, 2]]]\)], "Input"], Cell[TextData[{ "Here is an example that uses the matrix ", Cell[BoxData[ FormBox[ StyleBox[ RowBox[{"(", "\[NegativeThinSpace]", GridBox[{ { StyleBox["0", FontSize->10], StyleBox["1", FontSize->10]}, { StyleBox["1", FontSize->10], StyleBox["1", FontSize->10]} }], "\[NegativeThinSpace]", ")"}], GridBoxOptions->{RowSpacings->0.5}], TraditionalForm]]], ", which defines the Fibonacci recurrence, in that ", Cell[BoxData[ FormBox[ StyleBox[ RowBox[{"(", "\[NegativeThinSpace]", GridBox[{ { StyleBox["0", FontSize->10], StyleBox["1", FontSize->10]}, { StyleBox["1", FontSize->10], StyleBox["1", FontSize->10]} }], "\[NegativeThinSpace]", ")"}], GridBoxOptions->{RowSpacings->0.5}], TraditionalForm]]], Cell[BoxData[ FormBox[ StyleBox[ RowBox[{"(", "\[NegativeThinSpace]", GridBox[{ { StyleBox["a", FontSize->10]}, { StyleBox["b", FontSize->10]} }], "\[NegativeThinSpace]", ")"}], GridBoxOptions->{RowSpacings->0.5}], TraditionalForm]]], " = ", Cell[BoxData[ FormBox[ StyleBox[ RowBox[{"(", "\[NegativeThinSpace]", GridBox[{ { StyleBox["b", FontSize->10]}, { StyleBox[\(a + b\), FontSize->10]} }], "\[NegativeThinSpace]", ")"}], GridBoxOptions->{RowSpacings->0.5}], TraditionalForm]]], ". Using \[Infinity] as the modulus means that we get the unreduced \ numbers. Taking the 100th power should yield the 100th Fibonacci number." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(MatrixPowerMod[{{0, 1}, {1, 1}}, 100, \[Infinity]]\)], "Input"], Cell[BoxData[ \({{218922995834555169026, 354224848179261915075}, { 354224848179261915075, 573147844013817084101}}\)], "Output"] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ \(Fibonacci[100]\)], "Input"], Cell[BoxData[ \(354224848179261915075\)], "Output"] }, Closed]], Cell[TextData[{ "Now we can formulate the Perrin recurrence in terms of multiplication by ", Cell[BoxData[ TagBox[ StyleBox[ RowBox[{"(", GridBox[{ {"0", "1", "0"}, {"0", "0", "1"}, {"1", "1", "0"} }], ")"}], FontFamily->"Times", FontSize->9, GridBoxOptions->{RowSpacings->0.5}], (MatrixForm[ #]&)]]], ". Here, for example, is how to get ", Cell[BoxData[ \(TraditionalForm\`a\_7919\)]], " (mod 7919). Because 7919 is prime, we expect 0." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{"Mod", "[", RowBox[{ RowBox[{ RowBox[{"First", "[", RowBox[{"MatrixPowerMod", "[", RowBox[{ StyleBox[ RowBox[{"(", GridBox[{ {"0", "1", "0"}, {"0", "0", "1"}, {"1", "1", "0"} }], ")"}], FontFamily->"Times", FontSize->9, GridBoxOptions->{RowSpacings->0.5}], ",", "7919", ",", "7919"}], "]"}], "]"}], ".", \({3, 0, 2}\)}], ",", "7919"}], "]"}]], "Input"], Cell[BoxData[ \(0\)], "Output"] }, Closed]], Cell[TextData[{ "Of course, we can formulate this method for any linear recurrence. The ", StyleBox["LinearRecurrence", "Input"], " function that follows accepts the coefficients of the recurrence ({1, 1, \ 0} for Perrin's rule), the initial condition ({3, 0, 2} for Perrin), the \ desired index n, and a modulus m (which defaults to \[Infinity], giving the \ unreduced numbers). It forms the matrix, takes the mod-m power efficiently, \ and picks out the correct element of the matrix product with the initial \ vector. We also arrange the function to handle a list in its third argument." }], "Text"], Cell[TextData[ "LinearRecurrence[coeffs_, init_, n_Integer, m_:\[Infinity]] := \n \ Mod[First[MatrixPowerMod[\n Append[RotateRight /@ \n \ Drop[IdentityMatrix[Length[coeffs]], -1],\n coeffs], n, m]] . init, m] /; \ n \[GreaterEqual] Length[coeffs]\n\nLinearRecurrence[coeffs_, init_, \ n_Integer, m_:\[Infinity]] := Mod[\ninit\[LeftDoubleBracket]n+1\ \[RightDoubleBracket],m] /; n < Length[coeffs]\t\t\n\n\ LinearRecurrence[coeffs_, init_, n_List, m_:\[Infinity]] := \n \ LinearRecurrence[coeffs, init, #, m]& /@ n "], "Input", InitializationCell->True], Cell["As a check we look at the millionth prime.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(LinearRecurrence[{1, 1, 0}, {3, 0, 2}, Prime[10\^6], Prime[10\^6]]\)], "Input"], Cell[BoxData[ \(0\)], "Output"] }, Closed]], Cell["And now we can search for Perrin pseudorpimes.", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Select[Range[2, 500], \[InvisibleSpace]\(! PrimeQ[#]\) && LinearRecurrence[{1, 1, 0}, {3, 0, 2}, #, #] == 0\ &]\)], "Input"], Cell[BoxData[ \({}\)], "Output"] }, Closed]], Cell["\<\ The first one is quite large. Since I know where it is, I can focus \ the search very well!\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Select[Range[271200, 271500], \[InvisibleSpace]\(! PrimeQ[#]\) && \n\t\t LinearRecurrence[{1, 1, 0}, {3, 0, 2}, #, #] == 0\ &]\)], "Input"], Cell[BoxData[ \({271441}\)], "Output"] }, Closed]], Cell[TextData[{ "It is striking, especially when compared to other pseudoprime tests, how \ large the first failure of the Perrin test is. We can see that 271441, which \ is ", Cell[BoxData[ \(TraditionalForm\`521\^2\)]], ", is not a 2-pseudoprime:" }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(PowerMod[2, 271440, 271441]\)], "Input"], Cell[BoxData[ \(63042\)], "Output"] }, Closed]], Cell["\<\ So we might be tempted into thinking that Perrin + 2-psp has the \ potential of being a good primality test. But in fact there are many Perrin \ pseudoprimes that are even 2-strong pseudoprimes. The first such is \ 27,664,033.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[{ \(n = 3037*9109; \nLinearRecurrence[{1, 1, 0}, {3, 0, 2}, n, n]\), \(StrongPseudoprime[2, n]\)}], "Input"], Cell[BoxData[ \(0\)], "Output"] }, Closed]], Cell[BoxData[ \(True\)], "Input", FontWeight->"Plain"], Cell[TextData[{ "One can go farther with Perrin's idea. For example, one can run the \ recurrence backwards as well as forwards and discover that ", Cell[BoxData[ \(TraditionalForm\`a\_\(-p\) \[Congruent] 1\ \((mod\ p)\)\)]], "when p is prime. This and other variations, as well as references to \ papers with more detailed investigations, can be found in chapter 18 of \ [Wagon, 1998]. Such ideas lead to a Perrin test that is even stronger, in the \ sense that the first failure is quite large. But let us leave Perrin and move \ on to Lucas." }], "Text"] }, Closed]], Cell[CellGroupData[{ Cell[TextData[{ CounterBox["Section"], ". The Lucas Primality Test" }], "Section"], Cell[TextData[{ StyleBox["Mathematica", FontSlant->"Italic"], "'s ", StyleBox["PrimeQ", "Input"], " uses the 2- and 3-strong pseudoprime tests together with a Lucas test \ that we will now describe. No composite is known that passes these three \ tests, and all cases up to ", Cell[BoxData[ \(TraditionalForm\`10\^16\)]], " have been checked.\n\nA Lucas sequence is the sequence of integers \ obtained from the Fibonacci-lik recurrence ", Cell[BoxData[ \(TraditionalForm\`V\_0\ = \ 2\)]], ", ", Cell[BoxData[ \(TraditionalForm\`V\_1\ = \ P\)]], ", ", Cell[BoxData[ \(TraditionalForm\`V\_n\ = \ P\ V\_\(n - 1\) - \ Q\ V\_\(n - 2\)\)]], ", where \[CapitalDelta], defined to be ", Cell[BoxData[ \(TraditionalForm\`P\^2 - 4 Q\)]], ", is congruent to 0 or 1 (mod 4) and is not a perfect square. Since the \ special case Q = 1 will be our primary interest, we make that the default. As \ with the Perrin sequence, the V-sequence is related to the roots of a \ polynomial: if \[Alpha] and \[Beta] are the roots of the quadratic ", Cell[BoxData[ \(TraditionalForm\`\(x\^2 - P\[ThinSpace]x + Q\ \)\)]], ", then ", Cell[BoxData[ \(TraditionalForm\`V\_n\)]], "= ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^n + \[Beta]\^n\)]], "." }], "Text"], Cell[BoxData[ \(\(Lucas[P_, Q_: 1]\)[0]\ := \ 2; \n\(Lucas[P_, Q_: 1]\)[1]\ := \ P; \n\(Lucas[P_, Q_: 1]\)[n_]\ := \ \(\(Lucas[P, Q]\)[n]\ = \n\t\t\ P\ \ \(Lucas[P, Q]\)[n - 1]\ - \ Q\ \(Lucas[P, Q]\)[n - 2]\)\)], "Input"], Cell[TextData[{ "Because of recursion limits, the preceding definition will not work well \ except for quite small values of n. We will use the matrix method for serious \ work. We next summarize the facts that lead to an efficient prime test using \ the V-sequence.\n\nT", StyleBox["HEOREM ", FontSize->10], StyleBox["1", FontSize->11], ". Suppose p is prime, P is relatively prime to p, \[CapitalDelta] is not \ a square modulo p (where \[CapitalDelta] = ", Cell[BoxData[ \(TraditionalForm\`\(P\^2 - 4\ \)\)]], "), and ", Cell[BoxData[ \(TraditionalForm\`{V\_n}\)]], "\nis the Lucas sequence that starts with (2, P).\n\n(a). ", Cell[BoxData[ \(TraditionalForm\`V\_\(p + 1\) \[Congruent] 2\ \((mod\ p)\)\)]], "\n\n(b). ", Cell[BoxData[ \(TraditionalForm\`V\_\(2 i\)\ = \ V\_i\%2\ - \ 2\)]], "\n\n(c). ", Cell[BoxData[ \(TraditionalForm\`V\_\(2 i + 1\) = \(V\_i\) V\_\(i + 1\) - P\)]], "\n\n(d). ", Cell[BoxData[ FormBox[ RowBox[{ RowBox[{ FormBox[\(\ \(2 V\_\(i + 1\)\ = \)\ \), "TraditionalForm"], "P", " ", \(V\_i\)}], " ", "+", " ", \(\[CapitalDelta]\ \@\(\((V\_i\%2\ - \ 4)\)/\[CapitalDelta]\)\)}], TraditionalForm]]], "\n\nWe will not prove these identities; the reader can find proofs in \ [Bressoud, 1989, chap. 12]. These identities are, in essence, equivalent to \ the following result, a quadratic field analog of Fermat's little theorem \ that puts the primality test we are about to describe on a basis very similar \ to that of ordinary pseudoprimes (see [Bleichenbacher, 1996]).\n\nT", StyleBox["HEOREM ", FontSize->10], "2. Let p, P, and \[CapitalDelta] be as in Theorem 1. Let \[Alpha] be the \ larger root of ", Cell[BoxData[ \(TraditionalForm\`x\^2 - P\[ThinSpace]x + 1\ = \ 0\)]], "; \[Alpha] = ", Cell[BoxData[ FormBox[ RowBox[{ FormBox[ RowBox[{\(\[InvisiblePrefixScriptBase]\^1\), AdjustmentBox[\(/\_2\), BoxMargins->{{-0.3, 0}, {0, 0}}]}], "TraditionalForm"], \((P + \@\[CapitalDelta])\)}], TraditionalForm]]], ". Then ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\^n\)]], "\[Congruent] ", Cell[BoxData[ \(TraditionalForm\`\[Alpha]\&_\ \ \((mod\ p)\)\)]], ".\n\nOur first task is to define a routine to get the least P that works. \ All the work that follows will be restricted to odd inputs. Further, if P = \ \[PlusMinus]1 or \[PlusMinus]2, then ", Cell[BoxData[ \(TraditionalForm\`V\_n\)]], " is a periodic sequence, so we start our search for P at 3." }], "Text"], Cell[BoxData[ \(LucasParameter[n_?OddQ] := Module[{P = 3}, While[JacobiSymbol[P\^2 - 4, n] == 1 || GCD[n, P] \[NotEqual] 1, \(P++\)]; P]\)], "Input", GroupPageBreakWithin->True], Cell[TextData[{ "We can use this routine to illustrate Theorem 2. In the code that follows \ we use ", StyleBox["PowerMod", "Input"], " to rationalize the 2 in the denominator. But because ", Cell[BoxData[ FormBox[ SqrtBox[ StyleBox["k", FontFamily->"Courier", FontWeight->"Bold"]], TraditionalForm]]], " has the internal form ", StyleBox["Power[k,\[ThinSpace]Rational[1,\[ThinSpace]2]]", "Input"], ", we must hide the square root from the rule that transforms rationals. \ Then we add the square root back at the end." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(n = 61; P = LucasParameter[n]; \[CapitalDelta] = P\^2 - 4; \n \[Alpha] = x /. Last[ Solve[x\^2 - P\ x + 1 == 0, x] /. P\^2 - 4 \[Rule] \[CapitalDelta]] \)], "Input"], Cell[BoxData[ \(1\/2\ \((5 + \@21)\)\)], "Output"] }, Closed]], Cell[CellGroupData[{ Cell[BoxData[ \(\(\(ExpandAll[{\[Alpha], \[Alpha]\^n /. \@21 \[RuleDelayed] \(-\@21\)}] /. \@21 \[RuleDelayed] hide\) /. Rational[m_, 2] \[RuleDelayed] PowerMod[2, \ \(-1\), \ n]\) /. hide \[Rule] \@21\)], "Input"], Cell[BoxData[ \({31 + 31\ \@21, 31 + 31\ \@21}\)], "Output"] }, Closed]], Cell[TextData[{ "The reader can check some small odd composite values of n to see that the \ equation of Theorem 2 fails (to do this, add \"", StyleBox[ "/.\[ThinSpace]m_Integer\[ThinSpace]\[RuleDelayed]\[ThinSpace]Mod[m,\ \[ThinSpace]n]", "Input"], "\" to the end of the code). The Lucas test that follows is essentially a \ computationally efficient was of checking the congruence of Theorem 2.\n\nThe \ value of P that ", StyleBox["LucasParameter", "Input"], " gets will typically be quite small." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(Union[LucasParameter\ /@\ Range[3, 1000, 2]]\)], "Input"], Cell[BoxData[ \({3, 4, 5, 6, 9, 11, 15, 17, 21, 27, 29}\)], "Output"] }, Closed]], Cell["\<\ Now we can verify the Lucas condition (Theorem 1(a)) for primality \ in the case of 7.\ \>", "Text"], Cell[BoxData[ \(n\ = \ 7; \n ls = \(\(Lucas[LucasParameter[n]]\)\ [#]&\)\ /@\ Range[0, n + 1]; \n DisplayForm[GridBox[{Range[0, \ n + 1], ls, Mod[ls, n]}]]\)], "Input"], Cell[CellGroupData[{ Cell["Code for GridBox below", "Subsubsection"], Cell[CellGroupData[{ Cell[BoxData[ \(n\ = \ 7; \t ls = \(\(Lucas[LucasParameter[n]]\)\ [#]&\)\ /@\ Range[0, n + 1]; \n DisplayForm[ StyleBox[GridBox[\n \t\t\t{Prepend[Range[0, \ n + 1], "\"], \n\t\t\t\t Prepend[ls, "\"], \n\t\ \ \ \ \t Prepend[\tMod[ls, n], \ "\"]}, \n\tGridFrame -> 1, \ RowLines \[Rule] 1, ColumnAlignments -> {Center, Right}, \ \ ColumnLines\ -> \ {1, 0}], \ FontSize -> 9, \ Background -> GrayLevel[0.85]]]\)], "Input"], Cell[BoxData[ TagBox[ StyleBox[GridBox[{ {"n", "0", "1", "2", "3", "4", "5", "6", "7", "8"}, {"V", "2", "3", "7", "18", "47", "123", "322", "843", "2207"}, {\(V\ \((mod\ 7)\)\), "2", "3", "0", "4", "5", "4", "0", "3", "2"} }, ColumnAlignments->{Center, Right}, GridFrame->True, RowLines->True, ColumnLines->{True, False}], FontSize->9, Background->GrayLevel[0.850004]], DisplayForm]], "Output"] }, Closed]], Cell["\<\ And here is the last term of the V-sequence for the first 25 odd \ primes.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\(Mod[\(Lucas[LucasParameter[#]]\)\ [# + 1], #]\ &\)\ /@\ Prime[Range[2, 26]]\)], "Input"], Cell[BoxData[ \({2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2}\)], "Output"] }, Closed]], Cell[TextData[{ "But much more is true, and some ideas of the theory of strong pseudoprimes \ can be used [Baillie and Wagstaff, 1980]. Theorem 1(b) means that we can get \ ", Cell[BoxData[ \(TraditionalForm\`V\_\(n + 1\)\)]], " by first computing ", Cell[BoxData[ \(TraditionalForm\`V\_s\)]], " where s is the odd part of n+1, and then repeatedly squaring and \ subtracting 2 mod-n. Moreover, if n is really a prime and ", Cell[BoxData[ \(TraditionalForm\`V\^2 - \ 2 \[Congruent] \ 2\ \((mod\ n)\)\)]], ", then ", Cell[BoxData[ \(TraditionalForm\`V\^2\ \[Congruent] \ 4\)]], ", and this means that V must be one of \[PlusMinus]2 (mod n); this is \ analogous to the situation in the strong pseudoprime test, where the fact \ that the modular square root of 1 is one of \[PlusMinus]1 is used to great \ advantage. Further, the V-values that occur right after those just mentioned \ are important too: still assuming p is prime, if ", Cell[BoxData[ \(TraditionalForm\`V\_\(2 i\) \[Congruent] \ \(\[PlusMinus]2\)\)]], ", then ", Cell[BoxData[ \(TraditionalForm \`\(V\_\(2 i + 1\)\ \[Congruent] \ \(\[PlusMinus]P\)\ \)\)]], "; this follows from Theorem 1(d). And, finally, these next V-values are \ easily computed from Theorem 1(c) \[LongDash] ", Cell[BoxData[ \(TraditionalForm\`V\_\(2 i + 1\) = \(V\_i\) V\_\(i + 1\) - P\)]], " \[LongDash] so we can easily check this condition on the next V-value.\n\n\ Now we can form a strong Lucas pseudoprime test as follows:\n\n 1. Given \ n, let P be the Lucas parameter as defined earlier.\n 2. Write n + 1 as ", Cell[BoxData[ \(TraditionalForm\`s\ 2\^r\)]], ", where s is odd.\n 3. Get ", Cell[BoxData[ \(TraditionalForm\`V\_s\)]], "and ", Cell[BoxData[ \(TraditionalForm\`V\_\(s + 1\)\)]], " modulo n by the matrix method.\n 4. Form the two V-sequences indexed \ by {s, 2s, 4s, 8s, . . . , n+1} and {s+1, 2s+1, 4s+1, 8s+1, . . . , n+2} \ modulo n.\n 5. Check that ", Cell[BoxData[ \(TraditionalForm\`v\_\(n + 1\) = 2\)]], ", and that either the first sequence consists of all 2s or the first 2 in \ the sequence is preceded by a ", Cell[BoxData[ \(TraditionalForm\`\(-2\)\)]], ".\n 6. Check that the first 2 or ", Cell[BoxData[ \(TraditionalForm\`\(-2\)\)]], " in the first sequence is matched by a P or ", Cell[BoxData[ \(TraditionalForm\`\(-P\)\)]], " in the next sequence.\n\nIt would be more efficient to not compute the \ whole V-sequence, but to check the condition term-by-term. But, in part \ because we wish to look at the whole sequences here, we simply follow that \ approach. Exercise: Program the Lucas test more efficiently, checking each \ term of the sequence rather than precomputing the entire sequence. Here is \ code that implements these ideas. The ", StyleBox["lspspSequence", "Input"], " function uses ", StyleBox["LinearRecurrence", "Input"], " to get started and then ", StyleBox["NestList", "Input"], " to continue, and returns two lists of V-values, the first indexed by {s, \ 2s, 4s, 8s, . . . , n+1} and the second indexed by {s+1, 2s+1, 4s+1, 8s+1, . \ . . , n+2}. " }], "Text"], Cell[BoxData[ \(lspspSequence[P_, n_?OddQ] := Module[{oddData = OddPart[n + 1]}, \n\t\t Transpose[ NestList[ Mod[{#\[LeftDoubleBracket]1\[RightDoubleBracket]\^2 - 2, #\[LeftDoubleBracket]1\[RightDoubleBracket]\ # \[LeftDoubleBracket]2\[RightDoubleBracket] - P}, n]&, \n\t\t\t\t\ LinearRecurrence[{\(-1\), P}, {2, P}, oddData\[LeftDoubleBracket]1\[RightDoubleBracket] + {0, 1}, n], oddData\[LeftDoubleBracket]2\[RightDoubleBracket]] /. { n - 1 \[Rule] \(-1\), n - 2 \[Rule] \(-2\)}]]\)], "Input"], Cell[TextData[{ "Here is an example. The parameter P is 5 and a ", Cell[BoxData[ \(TraditionalForm\`\(-2\)\)]], " shows up with the corresponding value in the sequence of next terms being \ the expected ", Cell[BoxData[ \(TraditionalForm\`7919\ - \ 5\)]], ", or 7914." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(P\ = \ LucasParameter[7919]; \n{P, lspspSequence[LucasParameter[7919], 7919]}\)], "Input"], Cell[BoxData[ \({5, {{3883, 7830, 0, \(-2\), 2}, {4721, 7072, 4107, 7914, 5}}}\)], "Output"] }, Closed]], Cell["\<\ Now we can write some code to check the strong Lucas pseudoprime \ conditions in general. We will use the term strong Lucas pseudoprime do \ describe a composite integer that passes the test as just described, but the \ reader should be aware that the terminology in the literature is not \ standard. Others have used the term QF-psp and strong QF-psp to emphasize \ that this idea is a pseudoprime construction in a certain quadrtic number \ field.\ \>", "Text"], Cell[BoxData[ \(StrongLucasPSP[n_] := Module[{P = LucasParameter[n], lSeq, pos}, \n\t\t lSeq = lspspSequence[P, n]; \n\t\t lSeq\[LeftDoubleBracket]1, \(-1\)\[RightDoubleBracket] == 2 \[And] \n \t\t\((lSeq\[LeftDoubleBracket]{1, 2}, 1\[RightDoubleBracket] == { 2, P} || \n\t\t\t\t\t\ \((pos = Position[lSeq\[LeftDoubleBracket]1\[RightDoubleBracket], \(-2\)]; \n\t\t\t\t\t\t pos \[NotEqual] {} \[And] lSeq\[LeftDoubleBracket]2, pos\[LeftDoubleBracket]1, 1 \[RightDoubleBracket]\[RightDoubleBracket] == n - P) \))\)]\)], "Input"], Cell["\<\ We first check that a dozen random primes under 10 million pass the \ test.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(12 == Length[Select[Table[Prime[Random[Integer, {2, 10\^7}]], {12}], StrongLucasPSP]]\)], "Input"], Cell[BoxData[ \(True\)], "Output"] }, Closed]], Cell[TextData[{ "And now we search for strong Lucas pseudoprimes according to this \ definition. We avoid a ", StyleBox["Select[Range[]]", "Input"], " construction since that requires forming a huge range; instead we use a ", StyleBox["Do", "Input"], "-loop. We add a ", StyleBox["Print", "Input"], " statement to monitor progress since this is a time-consuming computation, \ and such a trace eases the task of stopping and restarting if desired." }], "Text"], Cell[BoxData[ \(LucasPSPs = {}; \n Do[If[Mod[i + 1, 2000] == 0, Print[{i, LucasPSPs}]]; \n\t\ \ \ \ \ If[GCD[i, 6] == 1 \[And] \[InvisibleSpace]\(! PrimeQ[i]\) \[And] StrongLucasPSP[i], \n\t\t\ \ \ \ \ \ AppendTo[LucasPSPs, i]], {i, 5, 100000, 2}]; \)], "Input"], Cell[CellGroupData[{ Cell["LucasPSPs", "Input"], Cell[BoxData[ \({989, 3239, 5777, 10877, 27971, 29681, 30739, 31631, 39059, 72389, 73919, 75077}\)], "Output"] }, Closed]], Cell["\<\ So there are 12 strong Lucas psps under 100,000. But they all would \ have their compositeness detected by a 2-spsp test! Indeed, even the simple \ 2-psp test suffices.\ \>", "Text"], Cell[CellGroupData[{ Cell[BoxData[ \(\((1 == PowerMod[2, # - 1, #]\ &)\)\ /@\ LucasPSPs\)], "Input"], Cell[BoxData[ \({False, False, False, False, False, False, False, False, False, False, False, False}\)], "Output"] }, Closed]], Cell[TextData[{ "Now, what the built-in ", StyleBox["PrimeQ", "Input"], " function does on input n is, in essence, the following:\n\n 1. If n \ has a small divisor, return ", StyleBox["False", "Input"], ".\n 2. If n is proved to be composite by the 2-strong pseudoprime \ test, return ", StyleBox["False", "Input"], ".\n 3. If n is proved to be composite by the 3-strong pseudoprime \ test, return ", StyleBox["False", "Input"], ".\n 4. If n is proved to be composite by the strong Lucas pseudoprime \ test, return ", StyleBox["False", "Input"], ".\n 5. Return ", StyleBox["True.", "Input"], "\n\nSo, if ", StyleBox["PrimeQ", "Input"], " declares that n is prime, is n definitely prime? Good question! The \ common belief is that such a combination of simple tests will ", StyleBox["not", FontSlant->"Italic"], " truly characterize primality (note that if such a finite combination of \ tests did characterize primality, that would mean that the primes are in the \ complexity class ", StyleBox["\[ScriptCapitalP]", FontWeight->"Bold"], "). On the other hand, Daniel Bleichenbacher [Bleichenbacher, 1996] has \ computed a list of all composite integers under ", Cell[BoxData[ \(TraditionalForm\`10\^16\)]], " that are both 2- and 3-strong pseudoprimes; there are 52,593 of them and \ he checked that ", StyleBox["PrimeQ", "Input"], " is correct (i.e., ", StyleBox["False", "Input"], ") for all of them. In other words, under ", Cell[BoxData[ \(TraditionalForm\`10\^16\)]], " there is no counterexample to the assertion that primality is determined \ by the conjunction of the 2-spsp test, the 3-spsp test, and the strong Lucas \ pseudoprime test. ", StyleBox["Mathematica", FontSlant->"Italic"], "'s", StyleBox[" ", FontSlant->"Italic"], StyleBox["PrimeQ", "Input"], " is actually a little different than the algorithm described here in that \ a slightly different method is used for the choice of Lucas parameter. But I \ checked that all the 2- and 3-spsps under ", Cell[BoxData[ \(TraditionalForm\`10\^16\)]], " fail the Lucas test as presented here too. So, for either method of \ choosing the Lucas parameter, we can be sure that the prime-query algorithm \ based on 2-spsp, 3-spsp, and the Lucas test tells the truth up to ", Cell[BoxData[ \(TraditionalForm\`10\^16\)]], "." }], "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[TextData[{ CounterBox["Section"], ". Certifiable Primality" }], "Section"], Cell[CellGroupData[{ Cell["References", "Subsection"], Cell[TextData[{ "[Bach and Shallit, 1996] Eric Bach and Jeffrey Shallit, ", StyleBox["Algorithmic Number Theory, Volume 1, Efficient Algorithms", FontSlant->"Italic"], ", MIT Press, Cambridge.\n\n[Baillie and Wagstaff, 1980] R. Baillie and S. \ S. Wagstaff, Lucas pseudoprimes, ", StyleBox["Math. Comp.", FontSlant->"Italic"], " 35, 1391\[Dash]1417.\n\n[Bleichenbacher, 1996] D. Bleichenbacher, \ Efficiency and Security of Cryptosystems based on Number Theory, Ph. D \ dissertation, Swiss Federal Institute of Technology, ETH Diss. No. 11404, Z\ \[UDoubleDot]rich.\n\n[Bressoud, 1989] D. M. Bressoud, ", StyleBox["Factorization and Primality Testing", FontSlant->"Italic"], ", Springer, New York.\n\n[Schlafly and Wagon, 1994] A. Schlafly and S. \ Wagon, Carmichael's conjecture is valid below ", Cell[BoxData[ \(TraditionalForm\`10\^\(10, 000, 000\)\)]], ", ", StyleBox["Mathematics of Computation", FontSlant->"Italic"], " 63, 415\[Dash]419.\n\n[Wagon, 1998] Stan Wagon, ", StyleBox["Mathematica in Action", FontSlant->"Italic"], ", Springer/TELOS, New York.\n\n\n" }], "Text"] }, Open ]] }, Open ]] }, Open ]] }, FrontEndVersion->"X 3.0", ScreenRectangle->{{0, 1600}, {0, 1280}}, AutoGeneratedPackage->None, WindowSize->{645, 656}, WindowMargins->{{Automatic, 245}, {Automatic, 60}} ] (*********************************************************************** Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. ***********************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[1709, 49, 215, 7, 67, "Text"], Cell[CellGroupData[{ Cell[1949, 60, 40, 0, 104, "Title"], Cell[1992, 62, 302, 7, 151, "Subtitle"], Cell[2297, 71, 33, 0, 64, "Subtitle"], Cell[2333, 73, 849, 15, 122, "Text"], Cell[CellGroupData[{ Cell[3207, 92, 78, 3, 52, "Section"], Cell[3288, 97, 750, 16, 104, "Text"], Cell[CellGroupData[{ Cell[4063, 117, 89, 2, 43, "Input"], Cell[4155, 121, 35, 1, 27, "Output"] }, Closed]], Cell[4205, 125, 338, 10, 50, "Text"], Cell[CellGroupData[{ Cell[4568, 139, 165, 2, 51, "Input"], Cell[4736, 143, 38, 1, 27, "Output"] }, Closed]], Cell[4789, 147, 95, 3, 32, "Text"], Cell[CellGroupData[{ Cell[4909, 154, 126, 2, 43, "Input"], Cell[5038, 158, 170, 2, 43, "Output"] }, Closed]], Cell[5223, 163, 238, 8, 32, "Text"], Cell[CellGroupData[{ Cell[5486, 175, 127, 2, 43, "Input"], Cell[5616, 179, 158, 2, 43, "Output"] }, Closed]], Cell[5789, 184, 1066, 16, 194, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[6892, 205, 102, 3, 32, "Section"], Cell[6997, 210, 1955, 52, 176, "Text"], Cell[8955, 264, 122, 3, 27, "Input"], Cell[CellGroupData[{ Cell[9102, 271, 61, 1, 27, "Input"], Cell[9166, 274, 62, 1, 27, "Output"] }, Closed]], Cell[9243, 278, 450, 15, 50, "Text"], Cell[9696, 295, 204, 4, 66, "Input"], Cell[9903, 301, 97, 3, 32, "Text"], Cell[CellGroupData[{ Cell[10025, 308, 54, 1, 27, "Input"], Cell[10082, 311, 40, 1, 27, "Output"] }, Closed]], Cell[CellGroupData[{ Cell[10159, 317, 56, 1, 27, "Input"], Cell[10218, 320, 55, 1, 27, "Output"] }, Closed]], Cell[10288, 324, 139, 3, 32, "Text"], Cell[CellGroupData[{ Cell[10452, 331, 53, 1, 27, "Input"], Cell[10508, 334, 44, 1, 27, "Output"] }, Closed]], Cell[10567, 338, 303, 4, 68, "Text"], Cell[10873, 344, 346, 6, 107, "Input"], Cell[11222, 352, 245, 6, 50, "Text"], Cell[CellGroupData[{ Cell[11492, 362, 145, 3, 27, "Input"], Cell[11640, 367, 64, 1, 27, "Output"] }, Closed]], Cell[11719, 371, 633, 9, 140, "Text"], Cell[CellGroupData[{ Cell[12377, 384, 156, 4, 59, "Input"], Cell[12536, 390, 697, 9, 139, "Output"] }, Closed]], Cell[13248, 402, 1227, 26, 176, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[14512, 433, 52, 0, 32, "Section"], Cell[14567, 435, 915, 19, 70, "Text"], Cell[CellGroupData[{ Cell[15507, 458, 176, 3, 70, "Input"], Cell[15686, 463, 52, 1, 70, "Output"] }, Closed]], Cell[15753, 467, 266, 7, 70, "Text"], Cell[16022, 476, 93, 2, 70, "Input"], Cell[16118, 480, 45, 0, 70, "Text"], Cell[CellGroupData[{ Cell[16188, 484, 48, 1, 70, "Input"], Cell[16239, 487, 121, 2, 70, "Output"] }, Closed]], Cell[16375, 492, 280, 7, 70, "Text"], Cell[CellGroupData[{ Cell[16680, 503, 70, 1, 70, "Input"], Cell[16753, 506, 145, 2, 70, "Output"] }, Closed]], Cell[16913, 511, 850, 21, 70, "Text"], Cell[17766, 534, 292, 5, 70, "Input"], Cell[18061, 541, 110, 3, 70, "Text"], Cell[CellGroupData[{ Cell[18196, 548, 32, 0, 70, "Input"], Cell[18231, 550, 166, 3, 70, "Output"] }, Closed]], Cell[18412, 556, 1808, 58, 70, "Text"], Cell[CellGroupData[{ Cell[20245, 618, 201, 4, 70, "Input"], Cell[20449, 624, 151, 2, 70, "Output"] }, Closed]], Cell[20615, 629, 2553, 65, 70, "Text"], Cell[CellGroupData[{ Cell[23193, 698, 119, 2, 70, "Input"], Cell[23315, 702, 36, 1, 70, "Output"] }, Closed]], Cell[23366, 706, 559, 11, 70, "Text"], Cell[CellGroupData[{ Cell[23950, 721, 68, 1, 70, "Input"], Cell[24021, 724, 39, 1, 70, "Output"] }, Closed]], Cell[24075, 728, 1331, 40, 70, "Text"], Cell[25409, 770, 149, 3, 70, "Input"], Cell[CellGroupData[{ Cell[25583, 777, 89, 1, 70, "Input"], Cell[25675, 780, 44, 1, 70, "Output"] }, Closed]], Cell[25734, 784, 197, 5, 70, "Text"], Cell[25934, 791, 183, 3, 70, "Input"], Cell[26120, 796, 2143, 62, 70, "Text"], Cell[CellGroupData[{ Cell[28288, 862, 83, 1, 70, "Input"], Cell[28374, 865, 139, 2, 70, "Output"] }, Closed]], Cell[CellGroupData[{ Cell[28550, 872, 47, 1, 70, "Input"], Cell[28600, 875, 55, 1, 70, "Output"] }, Closed]], Cell[28670, 879, 594, 19, 70, "Text"], Cell[CellGroupData[{ Cell[29289, 902, 639, 17, 70, "Input"], Cell[29931, 921, 35, 1, 70, "Output"] }, Closed]], Cell[29981, 925, 608, 9, 70, "Text"], Cell[30592, 936, 572, 9, 70, "Input", InitializationCell->True], Cell[31167, 947, 58, 0, 70, "Text"], Cell[CellGroupData[{ Cell[31250, 951, 102, 2, 70, "Input"], Cell[31355, 955, 35, 1, 70, "Output"] }, Closed]], Cell[31405, 959, 62, 0, 70, "Text"], Cell[CellGroupData[{ Cell[31492, 963, 162, 3, 70, "Input"], Cell[31657, 968, 36, 1, 70, "Output"] }, Closed]], Cell[31708, 972, 115, 3, 70, "Text"], Cell[CellGroupData[{ Cell[31848, 979, 176, 3, 70, "Input"], Cell[32027, 984, 42, 1, 70, "Output"] }, Closed]], Cell[32084, 988, 273, 7, 70, "Text"], Cell[CellGroupData[{ Cell[32382, 999, 60, 1, 70, "Input"], Cell[32445, 1002, 39, 1, 70, "Output"] }, Closed]], Cell[32499, 1006, 250, 5, 70, "Text"], Cell[CellGroupData[{ Cell[32774, 1015, 130, 2, 70, "Input"], Cell[32907, 1019, 35, 1, 70, "Output"] }, Closed]], Cell[32957, 1023, 60, 2, 70, "Input"], Cell[33020, 1027, 569, 10, 70, "Text"] }, Closed]], Cell[CellGroupData[{ Cell[33626, 1042, 86, 3, 32, "Section"], Cell[33715, 1047, 1328, 36, 158, "Text"], Cell[35046, 1085, 258, 5, 75, "Input"], Cell[35307, 1092, 2683, 67, 450, "Text"], Cell[37993, 1161, 210, 5, 47, "Input"], Cell[38206, 1168, 593, 15, 69, "Text"], Cell[CellGroupData[{ Cell[38824, 1187, 216, 5, 49, "Input"], Cell[39043, 1194, 54, 1, 43, "Output"] }, Closed]], Cell[CellGroupData[{ Cell[39134, 1200, 247, 4, 70, "Input"], Cell[39384, 1206, 64, 1, 33, "Output"] }, Closed]], Cell[39463, 1210, 524, 11, 104, "Text"], Cell[CellGroupData[{ Cell[40012, 1225, 77, 1, 27, "Input"], Cell[40092, 1228, 73, 1, 27, "Output"] }, Closed]], Cell[40180, 1232, 110, 3, 32, "Text"], Cell[40293, 1237, 181, 3, 59, "Input"], Cell[CellGroupData[{ Cell[40499, 1244, 47, 0, 42, "Subsubsection"], Cell[CellGroupData[{ Cell[40571, 1248, 524, 10, 139, "Input"], Cell[41098, 1260, 542, 15, 54, "Output"] }, Closed]], Cell[41655, 1278, 98, 3, 32, "Text"], Cell[CellGroupData[{ Cell[41778, 1285, 117, 2, 27, "Input"], Cell[41898, 1289, 116, 2, 27, "Output"] }, Closed]], Cell[42029, 1294, 3242, 72, 482, "Text"], Cell[45274, 1368, 669, 12, 87, "Input"], Cell[45946, 1382, 303, 9, 50, "Text"], Cell[CellGroupData[{ Cell[46274, 1395, 116, 2, 43, "Input"], Cell[46393, 1399, 99, 2, 27, "Output"] }, Closed]], Cell[46507, 1404, 473, 8, 86, "Text"], Cell[46983, 1414, 734, 14, 107, "Input"], Cell[47720, 1430, 99, 3, 32, "Text"], Cell[CellGroupData[{ Cell[47844, 1437, 136, 3, 47, "Input"], Cell[47983, 1442, 38, 1, 27, "Output"] }, Closed]], Cell[48036, 1446, 477, 11, 86, "Text"], Cell[48516, 1459, 297, 5, 75, "Input"], Cell[CellGroupData[{ Cell[48838, 1468, 26, 0, 27, "Input"], Cell[48867, 1470, 121, 2, 27, "Output"] }, Closed]], Cell[49003, 1475, 192, 4, 50, "Text"], Cell[CellGroupData[{ Cell[49220, 1483, 83, 1, 27, "Input"], Cell[49306, 1486, 125, 2, 43, "Output"] }, Closed]], Cell[49446, 1491, 2424, 60, 356, "Text"] }, Closed]] }, Closed]], Cell[CellGroupData[{ Cell[51919, 1557, 83, 3, 32, "Section"], Cell[CellGroupData[{ Cell[52027, 1564, 32, 0, 45, "Subsection"], Cell[52062, 1566, 1135, 25, 320, "Text"] }, Open ]] }, Open ]] }, Open ]] } ] *) (*********************************************************************** End of Mathematica Notebook file. ***********************************************************************)