(***********************************************************************
Mathematica-Compatible Notebook
This notebook can be used on any computer system with Mathematica 4.0,
MathReader 4.0, or any compatible application. The data for the notebook
starts with the line containing 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[ 28029, 812]*)
(*NotebookOutlinePosition[ 28917, 843]*)
(* CellTagsIndexPosition[ 28821, 837]*)
(*WindowFrame->Normal*)
Notebook[{
Cell[TextData[{
"MATHEMATICA IN ACTION ",
StyleBox["for issue 8-4\n of ",
FontColor->RGBColor[1, 0, 0]],
StyleBox["Mathematica in Education and Research",
FontSlant->"Italic",
FontColor->RGBColor[1, 0, 0]]
}], "Text",
TextAlignment->Center,
FontSize->16],
Cell["Certified Primes", "Title",
Evaluatable->False,
AspectRatioFixed->True],
Cell["\<\
An elegant recursive method for getting large integers that are \
certifiably prime outperforms traditional ways of generating mere \
\[OpenCurlyDoubleQuote]probable primes\[CloseCurlyDoubleQuote].\
\>", \
"Subtitle",
Evaluatable->False,
TextAlignment->Left,
AspectRatioFixed->True,
FontWeight->"Plain",
FontSlant->"Plain",
FontTracking->"Plain",
FontVariations->{"Underline"->False,
"Outline"->False,
"Shadow"->False}],
Cell["by Stan Wagon", "Subtitle",
CellMargins->{{Inherited, Inherited}, {18, Inherited}},
Evaluatable->False,
TextAlignment->Left,
AspectRatioFixed->True,
FontFamily->"Helvetica"],
Cell[TextData[{
"Suppose you want a random 100-digit prime. The obvious way to proceed is \
to use ",
StyleBox["PrimeQ", "Input",
FontWeight->"Plain"],
", as in the following code, which generates a random d-digit prime by \
choosing the first prime after a random d-digit integer."
}], "Text"],
Cell[BoxData[{
\(\(PrimeAfter[n_] :=
2 /; \(-2\) \[LessEqual] n \[LessEqual] 1;\)\), "\n",
\(\(PrimeAfter[n_] :=
Module[{j = n + 1 + Mod[n, 2]},
While[\[InvisibleSpace]\[Not] PrimeQ[j], j += 2];
j];\)\n\t\t\), "\n",
\(\(RandomPrime[d_] :=
Module[{n}, \n\ \ While[
n = PrimeAfter[Random[Integer, 10\^\(d - {1, 0}\)]];
Log[10. , n] > d]; n];\)\)}], "Input"],
Cell[CellGroupData[{
Cell[BoxData[
\(RandomPrime[5]\)], "Input"],
Cell[BoxData[
\(11549\)], "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[BoxData[
\(rp = RandomPrime[50]\)], "Input"],
Cell[BoxData[
\(32061110146573155597903708933477946846369567284001\)], "Output"]
}, Open ]],
Cell[TextData[{
"The problem with this approach is that it is conceivable (though no \
examples are known; see [Wagon 1999]) that ",
StyleBox["PrimeQ", "Input",
FontWeight->"Plain"],
" would return ",
StyleBox["True", "Input",
FontWeight->"Plain"],
" on a composite integer n. Thus it is of interest to find a method that \
produces integers that are certifiably prime. The primary motivation is \
esthetic, since it is annoying to have to say \[OpenCurlyDoubleQuote]very \
probably prime\[CloseCurlyDoubleQuote] when one wishes to say \
\[OpenCurlyDoubleQuote]prime\[CloseCurlyDoubleQuote]; in practical \
applications to cryptography one tends to not worry about the chance of \
getting a nonprime because this probability can be made very small with a few \
more pseudoprime tests. One could use the routines in the standard package ",
StyleBox["PrimeQ`", "Input",
FontWeight->"Plain"],
" to validate the output of ",
StyleBox["RandomPrime", "Input",
FontWeight->"Plain"],
". The ",
StyleBox["ProvablePrimeQ", "Input",
FontWeight->"Plain"],
" function in that package uses elliptic curve methods and works reasonably \
well, but it is slow and its theoretical underpinnings are quite advanced."
}], "Text"],
Cell[BoxData[
\(Needs["\"]\)], "Input"],
Cell[CellGroupData[{
Cell[BoxData[
\(ProvablePrimeQ[rp] // Timing\)], "Input"],
Cell[BoxData[
\({22.683333333333337`\ Second, True}\)], "Output"]
}, Open ]],
Cell["\<\
The method that we will describe in this article (due to U. Maurer \
[Maurer 1995]) produces certified primes very quickly and is both simple and \
elegant.\
\>", "Text"],
Cell[CellGroupData[{
Cell[TextData[{
CounterBox["Section"],
". The Theory"
}], "Section"],
Cell[TextData[{
"Maurer\[CloseCurlyQuote]s work relies on the following result, discovered \
by H. C. Pocklingon in 1914. The main idea is similar to the idea of Pratt \
certificates discussed in an earlier column of mine [Wagon 1997], but with \
the important new point that the witness, b, need only work for a partial \
factorization of n. We give a simplified version of Pocklington\
\[CloseCurlyQuote]s theorem where the partial factorization consists of a \
single prime q; this is adequate to the task at hand. For the general \
version, which handles the case that the factored part of ",
Cell[BoxData[
\(TraditionalForm\`n - 1\)]],
" has several primes (as opposed to just one), see [Bressoud and Wagon, \
2000]."
}], "Text"],
Cell[TextData[{
"P",
StyleBox["OCKLINGTON\[CloseCurlyQuote]S", "SmallText"],
" T",
StyleBox["HEOREM", "SmallText"],
". Suppose ",
Cell[BoxData[
\(TraditionalForm\`n - 1\)]],
" has the prime factor q, where ",
Cell[BoxData[
\(TraditionalForm\`q > \@n\)]],
" and the following hold:"
}], "Text",
CellTags->"Thm: Pratt Certificate Shorter"],
Cell[TextData[{
"\t1. ",
Cell[BoxData[
\(TraditionalForm\`2\^\(n - 1\) \[Congruent] 1\ \((mod\ n)\)\)]],
"\n\n\t2. ",
Cell[BoxData[
FormBox[
RowBox[{
RowBox[{"gcd",
RowBox[{"(",
RowBox[{
RowBox[{
StyleBox[\(2\^\(\((n - 1)\)/q\)\),
FontFamily->"Times"], "-", "1"}], ",", " ", "n"}],
")"}]}], "=", "1"}], TextForm]]],
"."
}], "Text"],
Cell["Then n is prime.", "Text"],
Cell[TextData[{
"P",
StyleBox["ROOF", "SmallText"],
". Suppose n is not prime and p is its smallest prime divisor; then ",
Cell[BoxData[
\(TraditionalForm\`p < \@n\)]],
" and p is odd (because, by (1), n is odd). Let d be the mod-p order of 2. \
Then, by standard results, d divides ",
Cell[BoxData[
\(TraditionalForm\`p - 1\)]],
" (recall that Fermat\[CloseCurlyQuote]s little theorem tells us that ",
Cell[BoxData[
\(TraditionalForm\`2\^\(p - 1\) \[Congruent] 1\ \((mod\ p)\)\)]],
"). Further, by (1), d divides ",
Cell[BoxData[
\(TraditionalForm\`n - 1\)]],
". But d does not divide ",
Cell[BoxData[
\(TraditionalForm\`\((n - 1)\)/q\)]],
", because doing so would violate (2) (for then ",
Cell[BoxData[
FormBox[
RowBox[{
StyleBox[\(2\^\(\((n - 1)\)/q\)\),
FontFamily->"Times"], "\[Congruent]", \(1\ \((mod\ p)\)\)}],
TextForm]]],
" and so p would divide both n and ",
Cell[BoxData[
\(TraditionalForm\`2\^\(\((n - 1)\)/q\) - 1\)]],
"). Therefore q divides d. But d divides ",
Cell[BoxData[
\(TraditionalForm\`p - 1\)]],
", so q divides ",
Cell[BoxData[
\(TraditionalForm\`p - 1\)]],
", which means that ",
Cell[BoxData[
\(TraditionalForm\`p - 1 \[GreaterEqual] q \[GreaterEqual] \@n\)]],
", a contradiction. \[Square]\n\nOf course, any integer b can serve as a \
witness in place of 2; but we may as well restrict to the choice of 2 \
because, if n is truly prime, the chance that 2 is not a witness is very \
small (about ",
Cell[BoxData[
\(TraditionalForm\`1/\@n\)]],
")."
}], "Text"],
Cell[TextData[{
"And now here is Maurer\[CloseCurlyQuote]s method of generating certified \
primes. Suppose we want p, a d-digit prime that is certifiably prime. Select \
a prime q having ",
Cell[BoxData[
\(TraditionalForm\`\[LeftCeiling]d\/2\[RightCeiling] + 1\)]],
" many digits (q should be a ",
StyleBox["certified",
FontSlant->"Italic"],
" prime: this can be done recursively). Then choose a random integer ",
Cell[BoxData[
\(TraditionalForm\`r < q\)]],
" (choosing repeatedly) until ",
Cell[BoxData[
\(TraditionalForm\`p = r\ q + 1\)]],
" can be declared probably prime by ",
StyleBox["PrimeQ", "Input",
FontWeight->"Plain"],
" (or your favorite pseudoprime test). Because q is prime and greater than \
",
Cell[BoxData[
\(TraditionalForm\`\@p\)]],
", we have enough of a factorization of ",
Cell[BoxData[
\(TraditionalForm\`p - 1\)]],
" to apply Pocklington\[CloseCurlyQuote]s Theorem and check that 2 is a \
witness to p\[CloseCurlyQuote]s primality. Below ",
Cell[BoxData[
\(TraditionalForm\`10\^16\)]],
" we can simply use ",
StyleBox["PrimeQ", "Input",
FontWeight->"Plain"],
", which is known to tell the truth in that domain."
}], "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell[TextData[{
CounterBox["Section"],
". The Implementation"
}], "Section"],
Cell[TextData[{
"Because ",
StyleBox["PrimeQ", "Input",
FontWeight->"Plain"],
" is known to be 100% correct below ",
Cell[BoxData[
\(TraditionalForm\`10\^16\)]],
", we can consider primes below that bound to be self-certifying. However, \
experiments show that it is faster to lower the crossover point to ",
Cell[BoxData[
\(TraditionalForm\`10\^7\)]],
". The ",
StyleBox["QuickCompositeQ", "Input",
FontWeight->"Plain"],
" function performs a fast check for compositeness with a 2-pseudoprime \
test, preceded by a check for small divisors. When looking for the random \
multiplier r, we take some care (by looking at the mod-3 class of q) to put \
it in the right mod 6 class so that ",
Cell[BoxData[
\(TraditionalForm\`1 +
q\ r \[Congruent] \(\[PlusMinus]1\)\ \ \((mod\ 6)\)\)]],
", and so has a chance of being prime. We do this by choosing r and then \
subtracting ",
Cell[BoxData[
\(TraditionalForm\`Mod[r, 6]\)]],
" to get a multiple of 6; then we add the small correction (0 or 4 or 2) \
based on ",
Cell[BoxData[
\(TraditionalForm\`Mod[q, 3]\)]],
". But if this fails to yield a prime value of ",
Cell[BoxData[
\(TraditionalForm\`1 + r\ q\)]],
" there is really no need to pick another large random number r; it is \
sufficient to add 6 to r\[CloseCurlyQuote]s current value and try again (and \
again and again). The use of the decimal point in the bounds in the first ",
StyleBox["Random", "Input",
FontWeight->"Plain"],
" speeds things up by avoiding rational arithmetic; and the ",
StyleBox["If", "Input",
FontWeight->"Plain"],
" statement that follows avoids the mod-3 reduction of ",
StyleBox["q", "Input",
FontWeight->"Plain"],
" whenever the random choice from ",
Cell[BoxData[
\(TraditionalForm\`{0, 1}\)]],
" is 0. ",
StyleBox["Mod[q,3,2]", "Input",
FontWeight->"Plain"],
" uses the new modular offset feature in version 4; in this case the \
modulus is 3, so the offset of 2 gives a residue in ",
Cell[BoxData[
\(TraditionalForm\`{2, 3, 4}\)]],
"."
}], "Text"],
Cell[BoxData[
\(primeProd = Times@@Prime[Range[3, 150]]; \n\n
QuickCompositeQ[n_] :=
1 < GCD[n, primeProd] < n || PowerMod[2, n - 1, n] \[NotEqual] 1; \n\n
Options[CertifiedPrime] = {ShowCertificate \[Rule] False}; \n\n
CertifiedPrime[d_, opts___] :=
\(If[\(ShowCertificate /. \[InvisibleSpace]{opts}\) /.
\[InvisibleSpace]Options[CertifiedPrime], Flatten, First]\)[
PrimeAndCertificate[d]]; \n\n
PrimeAndCertificate[d_] :=
Module[\[IndentingNewLine]{p, r, q,
q1 = PrimeAndCertificate[Ceiling[d\/2] + 1]},
\[IndentingNewLine]q = q1\[LeftDoubleBracket]1\[RightDoubleBracket];
\[IndentingNewLine]r =
Random[Integer, {10. \^\(d - 1\), 10. \^d}\/q];
\[IndentingNewLine]r = r - Mod[r, 6] + adjust; \n\t\t
While[r = r + 6; p = 1 + q\ r; \[IndentingNewLine]\ \
QuickCompositeQ[p] || GCD[p, PowerMod[2, r, p] - 1] \[NotEqual] 1];
\[IndentingNewLine]{p, q1}]; \[IndentingNewLine]\n\n
PrimeAndCertificate[d_ /; d < 8] := {q\_0 = RandomPrime[d];
adjust = 1 - PowerMod[q\_0, \(-1\), 6]; q\_0}\)], "Input"],
Cell[CellGroupData[{
Cell[BoxData[
\(CertifiedPrime[80]\)], "Input"],
Cell[BoxData[
\(716015818822597323106370182923784258251457719053500230359529255810823709\
66160659\)], "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[BoxData[
\(cert = CertifiedPrime[80, ShowCertificate \[Rule] True]\)], "Input"],
Cell[BoxData[
\({12922585523681472466982853331520255260838539501843702770149527805065244\
858734559, 14656315658677753155342729356277390667109, 2435069850302794312841,
757342063721, 9204449}\)], "Output"]
}, Open ]],
Cell["\<\
And we can easily check that all the proper conditions hold \
throughout the certificate.\
\>", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(Table[{p, q} =
cert\[LeftDoubleBracket]{i,
i + 1}\[RightDoubleBracket]; \n\t\ \ \ \ \ \ \ \ \ \ \ {q > \@p,
PowerMod[2, p - 1, p],
GCD[p, PowerMod[2, \((p - 1)\)/q, \ p] - 1]}, \ \n\t{i, \
Length[cert] - 1}]\)], "Input"],
Cell[BoxData[
\({{True, 1, 1}, {True, 1, 1}, {True, 1, 1}, {True, 1, 1}}\)], "Output"]
}, Open ]],
Cell["We can easily get quite large primes.", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(CertifiedPrime[300]\)], "Input"],
Cell[BoxData[
\(717594999266282283078146644939214989156028944821874092350708699524052503\
703682086625283362917371193676459708295513249755943105366317479770956943728794\
872152425301125657665231181489266649994289975210281535941619488026900700832600\
149496476194967307644263896877277578316553770365520324099902579353733281\)], \
"Output"]
}, Open ]]
}, Open ]],
Cell[CellGroupData[{
Cell[TextData[{
CounterBox["Section"],
". Timing Surprises"
}], "Section"],
Cell[TextData[{
"Now for a surprise: Maurer\[CloseCurlyQuote]s method is ",
StyleBox["faster",
FontSlant->"Italic"],
" than standard methods of generating probable primes, while at the same \
time having the great advantage that there is no doubt about the status of \
the output!\n",
"\n",
"We will take several steps to make a fair timing comparison. The \
randomness in the process makes a true comparison difficult, but the results \
in this section will provide solid evidence that Maurer\[CloseCurlyQuote]s \
algorithm is just as fast as, and often faster than, standard methods for \
getting a random probable prime. We will optimize the ",
StyleBox["CertifiedPrime", "Input",
FontWeight->"Plain"],
" code in several ways:\n\n\[Bullet] There is no need to keep track of the \
certificate (the smaller primes q). For security reasons this is sound, since \
knowledge of the factors of ",
Cell[BoxData[
\(TraditionalForm\`p - 1\)]],
" can help one factor a number of the form ",
Cell[BoxData[
\(TraditionalForm\`p\ p\_1\)]],
". The prime returned will be a certified prime, but neither the user nor \
anyone else will be able to check its certificate since it is too difficult \
to find the large prime factor of ",
Cell[BoxData[
\(TraditionalForm\`p - 1\)]],
". \n\n\[Bullet] We replace the recursive method with an iterative method: \
first we generate the list of digit-sizes that will show up; then we reverse \
this list and work upward, using ",
StyleBox["Fold", "Input",
FontWeight->"Plain"],
" with ",
StyleBox["#1", "Input",
FontWeight->"Plain"],
" representing q, the certified prime in hand, and ",
StyleBox["#2", "Input",
FontWeight->"Plain"],
" representing d, the number of digits.\n\n\[Bullet] We will choose r to be \
in the right congruence class modulo 30, so that ",
Cell[BoxData[
\(TraditionalForm\`1 + q\ r\)]],
" is congruent to one of ",
Cell[BoxData[
\(TraditionalForm\`{1, 7, 11, 13, 17, 19, 23, 29}\)]],
" modulo 30. In fact, if we note the mod-30 class of the first prime \
generated ",
Cell[BoxData[
\(TraditionalForm\`\((q\_0)\)\)]],
", then we can arrange that each subsequent q has the same mod-30 class. \
This is handled with the ",
StyleBox["adjust", "Input",
FontWeight->"Plain"],
" variable, which is used to adjust the random r upward after it is reduced \
downward to be 0 (mod 30).\n\n\[Bullet] ",
StyleBox["Quotient", "Input",
FontWeight->"Plain"],
" is faster than ",
StyleBox["Ceiling", "Input",
FontWeight->"Plain"],
".\n\n\[Bullet] The gcd part of the quick test for compositeness is done in \
two pieces, the first using ",
Cell[BoxData[
\(TraditionalForm\`7\[CenterDot]11\[CenterDot]13\[CenterDot]17\
\[CenterDot]19\)]],
", thus allowing it to more often work with machine integers.\n\nTo be \
scrupulously fair we should avoid ",
StyleBox["PrimeQ", "Input",
FontWeight->"Plain"],
" in the generation of the smallest random primes. But this is done only \
once for each call, and the timing issue on this point is insignificant."
}], "Text"],
Cell[BoxData[{
\(\(p7To150 =
Times @@ Prime[Range[7, 150]];\)\[IndentingNewLine]\), "\n",
\(\(CertifiedPrime[d_]\ := \
Module[{q0, adjust, r,
p, \n\ \ \ \ \ \ \ \ sizes =
Rest[Reverse[
FixedPointList[If[# > 10, Quotient[# + 1, 2] + 1, #] &,
d]]]}, \[IndentingNewLine]\[IndentingNewLine] (*\
Get\ q0, \
the\ smallest\ of\ the\ certified\ primes\ \
*) \[IndentingNewLine]q0 =
Random[Integer, \ {10. \^\(First[sizes] - 1\), \
10. \^First[sizes]}]; \[IndentingNewLine]q0 += \
1 - Mod[q0, 2]; \[IndentingNewLine]While[\(! PrimeQ[q0]\),
q0 += 2]; \[IndentingNewLine]\[IndentingNewLine]If[d > 10,
adjust =
1 + PowerMod[\(-q0\), \(-1\),
30030]]; \[IndentingNewLine]Fold[\((\ \ r = \((\(\((# +
Mod[\(-#\), 30030] + adjust)\) &\)\ [
Random[Integer, {10. \^\(#2 - 1\), 10. \^#2}\/\(\(\ \
\)\(#1\)\)]])\); \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ While[
p = 1 + #1\ *r; \[IndentingNewLine]\((1 < GCD[p, p7To150] <
p)\)\ \[Or] \[IndentingNewLine]\ \((PowerMod[2,
p - 1, p] \[NotEqual]
1)\) \[Or] \ \[IndentingNewLine]\((GCD[p,
PowerMod[2, r, p]\ - \ 1] \[NotEqual] 1)\),
r += 30030]; \[IndentingNewLine]p)\) &\ , q0,
Rest[sizes]]];\)\)}], "Input"],
Cell[CellGroupData[{
Cell[BoxData[
\(CertifiedPrime[100] // Timing\)], "Input"],
Cell[BoxData[
\({0.7000000000000028`\ Second,
958259328023165216549412278282311999891698008437309946015842129539007853\
5356058560119325336825576359}\)], "Output"]
}, Open ]],
Cell["\<\
These optimizations yield a nicely fast routine. Here is the \
average time needed for a 100-digit prime, based on 1000 trials.\
\>", "Text"],
Cell[CellGroupData[{
Cell[BoxData[{
\(\(trials = 1000;\)\), "\n",
\(Plus @@ Table[\(Timing[CertifiedPrime[100]]\)\[LeftDoubleBracket]1\
\[RightDoubleBracket], {trials}]\/trials\)}], "Input"],
Cell[BoxData[
\(1.1244666666666667`\ Second\)], "Output"]
}, Open ]],
Cell[TextData[{
"If we wish to compare this with standard methods of generating random \
probable primes (sometimes known as \[OpenCurlyDoubleQuote]industrial grade \
primes\[CloseCurlyDoubleQuote]), we have several choices (see [Wagon 1997, \
1999] for the definitions).\n\n1. Use the built-in ",
StyleBox["PrimeQ", "Input",
FontWeight->"Plain"],
", which calls n prime if it passes the 2- and 3-strong pseudoprime tests \
and a Lucas pseudoprime test.\n\n2. Use a method that uses only a single \
2-strong psp test.\n\n3. Use 2- and 3- strong psp tests.\n\n4. Use 2-, 3-, \
and 5-strong psp tests.\n\nMethod 2 is surely the fastest and will be \
adequately reliable for large integers, where the probability of running into \
a 2-strong psp is very small. Methods 3 and 4 add some extra security, while \
method 1 is probably the best since not a single counterexample is currently \
known. "
}], "Text"],
Cell[TextData[{
"Coding the strong pseudoprime test gives us an opportunity to show off two \
new and useful version 4 functions: ",
StyleBox["NestWhile", "Input",
FontWeight->"Plain"],
" and ",
StyleBox["IntegerExponent", "Input",
FontWeight->"Plain"],
". To perform a b-strong pseudoprime test on n one first forms ",
Cell[BoxData[
\(TraditionalForm\`s = b\^m\ \((mod\ n)\)\)]],
" where m is the odd part of ",
Cell[BoxData[
\(TraditionalForm\`n - 1\)]],
". Then one squares repeatedly (mod n) until the exponent reaches ",
Cell[BoxData[
\(TraditionalForm\`n - 1\)]],
". If ",
Cell[BoxData[
\(TraditionalForm\`n - 1\)]],
" shows up before the penultimate step, n passes the test; if it does not, \
n fails and is definitely composite. Thus we wish to iterate the modular \
squaring operation, exiting the ",
StyleBox["Nest", "Input",
FontWeight->"Plain"],
" sequence as soon as a condition is met. There was no simple way to do \
this in version 3. But now we can use:"
}], "Text"],
Cell[BoxData[
StyleBox[
\(n - 1 \[Equal]
NestWhile[PowerMod[#, 2, n]&, s, # \[NotEqual] n - 1\ &, 1, r - 1]\),
"Input"]], "Input"],
Cell[TextData[{
"The third argument tells us under what condition to continue the \
iteration; the fourth argument gives the maximum number of iterations. If the \
result of this is ",
Cell[BoxData[
\(TraditionalForm\`n - 1\)]],
", the number passes; otherwise it fails."
}], "Text"],
Cell[TextData[{
"So here is code that returns ",
StyleBox["True", "Input",
FontWeight->"Plain"],
" if n passes the b-strong pseudoprime test, followed by the ",
StyleBox["PrimeAfterSPSP", "Input",
FontWeight->"Plain"],
" and ",
StyleBox["RandomPrimeSPSP", "Input",
FontWeight->"Plain"],
" functions that make use of the ",
StyleBox["StrongPseudoprimeTest", "Input",
FontWeight->"Plain"],
" code. We split the gcd step into two as before."
}], "Text"],
Cell[BoxData[{
\(\(Attributes[StrongPseudoprimeTest] =
Listable;\)\n\), "\[IndentingNewLine]",
\(\(StrongPseudoprimeTest[b_Integer, n_]\ := \
Module[\[IndentingNewLine]{r = IntegerExponent[n - 1, 2], s}, \n\t
s = PowerMod[b, \(n - 1\)\/2\^r, n]; \n\t
s \[Equal] 1 \[Or]
s \[Equal] n - 1 \[Or] \n\t\ \ \ \ \ n - 1 \[Equal]
NestWhile[PowerMod[#, 2, n] &, s, # \[NotEqual] n - 1\ &, 1,
r - 1]];\)\)}], "Input"],
Cell["\<\
Here are some examples; 2047 is the first composite that passes the \
2-strong pseudoprime test; 1373653 is the first composite that passes both \
the 2- and 3-strong psp tests.\
\>", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(StrongPseudoprimeTest[
2, {101, 541, \ 7919, \ 2047, \ 1373653}]\)], "Input"],
Cell[BoxData[
\({True, True, True, True, True}\)], "Output"]
}, Open ]],
Cell["But 2047 is not a 3-strong pseudoprime.", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(StrongPseudoprimeTest[
3, {101, 541, \ 7919, \ 2047, \ 1373653}]\)], "Input"],
Cell[BoxData[
\({True, True, True, False, True}\)], "Output"]
}, Open ]],
Cell[TextData[{
"Now we modify the ",
StyleBox["PrimeAfter", "Input",
FontWeight->"Plain"],
" and ",
StyleBox["RandomPrime", "Input",
FontWeight->"Plain"],
" code so that the 2-strong pseudoprime test is used."
}], "Text"],
Cell[BoxData[{
\(\(PrimeAfterSPSP2[n_] :=
Module[{j = n + 1 + Mod[n, 2]}, \[IndentingNewLine]While[
1 < GCD[j, 15015] < j \[Or]
1 < GCD[j, p7To150] <
j\ \[Or] \ \[IndentingNewLine]\[InvisibleSpace]\[Not] \
StrongPseudoprimeTest[2, j]\ \ \[InvisibleSpace], \[IndentingNewLine]j +=
2]; \[IndentingNewLine]j];\)\[IndentingNewLine]\), "\
\[IndentingNewLine]",
\(\(RandomPrimeSPSP2[d_]\ := \
Module[{n},
While[n =
PrimeAfterSPSP2[
Random[Integer,
10\^\(d - {1, 0}\)]]; \ \[IndentingNewLine]\ \ Log[10. ,
n] > d]; \ n];\)\)}], "Input"],
Cell[CellGroupData[{
Cell[BoxData[
\(Plus @@ Table[\(Timing[RandomPrimeSPSP2[100]]\)\[LeftDoubleBracket]1\
\[RightDoubleBracket], {trials}]\/trials\)], "Input"],
Cell[BoxData[
\(1.0269500000000011`\ Second\)], "Output"]
}, Open ]],
Cell[TextData[{
"The timing experiment shows that the generation of certified primes takes \
essentially the same amount of time as the generation of probably primes. And \
when more sophisticated probable prime generators are used, such as ",
StyleBox["PrimeQ", "Input",
FontWeight->"Plain"],
", then the certified-prime algorithm takes less time!"
}], "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell["References", "Subsection"],
Cell[TextData[{
"[Bressoud and Wagon 2000] D. Bressoud and S. Wagon, ",
StyleBox["A Course in Computational Number Theory,",
FontSlant->"Italic"],
" Springer-Verlag, New York.\n[Maurer 1995] U. Maurer, Fast generation of \
prime numbers and secure public-key cryptographic parameters, ",
StyleBox["Journal of Cryptology",
FontSlant->"Italic"],
" 9, 123\[Dash]155.\n[Wagon 1997] S. Wagon, What is a prime number?, ",
StyleBox["Mathematica in Education and Research",
FontSlant->"Italic"],
" 6:4, 54\[Dash]61.\n[Wagon 1999] S. Wagon, ",
StyleBox["Mathematica in Action",
FontSlant->"Italic"],
", Springer/TELOS, New York.\n"
}], "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell["About the Author", "Subsection"],
Cell[TextData[{
"Stan Wagon uses ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" extensively in his research, teaching, and exposition, and is especially \
appreciative of the new ways ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" allows us to look at mathematical objects. He is author or coauthor of \
several books, including the just-released second edition of ",
StyleBox["Mathematica in Action",
FontSlant->"Italic"],
", ",
StyleBox["Animating Calculus",
FontSlant->"Italic"],
", ",
StyleBox["VisualDSolve",
FontSlant->"Italic"],
", and ",
StyleBox["Which Way Did the Bicycle Go?",
FontSlant->"Italic"],
". He teaches a high-altitude course, ",
StyleBox["Rocky Mountain Mathematica",
FontSlant->"Italic"],
", every summer in the mountains of Colorado.\n\nStan Wagon\nDepartment of \
Mathematics and Computer Science\nMacalester College\nSt. Paul, MN 55105\n\
wagon@macalester.edu"
}], "Text"]
}, Open ]]
},
FrontEndVersion->"4.0 for Microsoft Windows",
ScreenRectangle->{{0, 1024}, {0, 695}},
AutoGeneratedPackage->None,
CellGrouping->Manual,
WindowSize->{839, 673},
WindowMargins->{{7, Automatic}, {Automatic, 0}}
]
(***********************************************************************
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->{
"Thm: Pratt Certificate Shorter"->{
Cell[6302, 194, 370, 13, 33, "Text",
CellTags->"Thm: Pratt Certificate Shorter"]}
}
*)
(*CellTagsIndex
CellTagsIndex->{
{"Thm: Pratt Certificate Shorter", 28703, 830}
}
*)
(*NotebookFileOutline
Notebook[{
Cell[1717, 49, 276, 9, 60, "Text"],
Cell[1996, 60, 81, 2, 105, "Title",
Evaluatable->False],
Cell[2080, 64, 447, 14, 93, "Subtitle",
Evaluatable->False],
Cell[2530, 80, 190, 5, 62, "Subtitle",
Evaluatable->False],
Cell[2723, 87, 304, 7, 52, "Text"],
Cell[3030, 96, 443, 10, 112, "Input"],
Cell[CellGroupData[{
Cell[3498, 110, 47, 1, 30, "Input"],
Cell[3548, 113, 39, 1, 29, "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[3624, 119, 53, 1, 30, "Input"],
Cell[3680, 122, 84, 1, 29, "Output"]
}, Open ]],
Cell[3779, 126, 1251, 26, 128, "Text"],
Cell[5033, 154, 66, 1, 30, "Input"],
Cell[CellGroupData[{
Cell[5124, 159, 61, 1, 30, "Input"],
Cell[5188, 162, 69, 1, 29, "Output"]
}, Open ]],
Cell[5272, 166, 180, 4, 33, "Text"],
Cell[CellGroupData[{
Cell[5477, 174, 73, 3, 53, "Section"],
Cell[5553, 179, 746, 13, 90, "Text"],
Cell[6302, 194, 370, 13, 33, "Text",
CellTags->"Thm: Pratt Certificate Shorter"],
Cell[6675, 209, 470, 16, 71, "Text"],
Cell[7148, 227, 32, 0, 33, "Text"],
Cell[7183, 229, 1647, 45, 147, "Text"],
Cell[8833, 276, 1235, 33, 112, "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell[10105, 314, 84, 4, 53, "Section"],
Cell[10192, 320, 2121, 54, 166, "Text"],
Cell[12316, 376, 1145, 20, 446, "Input"],
Cell[CellGroupData[{
Cell[13486, 400, 51, 1, 30, "Input"],
Cell[13540, 403, 116, 2, 29, "Output"]
}, Open ]],
Cell[CellGroupData[{
Cell[13693, 410, 88, 1, 30, "Input"],
Cell[13784, 413, 215, 3, 48, "Output"]
}, Open ]],
Cell[14014, 419, 113, 3, 33, "Text"],
Cell[CellGroupData[{
Cell[14152, 426, 289, 6, 75, "Input"],
Cell[14444, 434, 90, 1, 29, "Output"]
}, Open ]],
Cell[14549, 438, 53, 0, 33, "Text"],
Cell[CellGroupData[{
Cell[14627, 442, 52, 1, 30, "Input"],
Cell[14682, 445, 342, 5, 67, "Output"]
}, Open ]]
}, Open ]],
Cell[CellGroupData[{
Cell[15073, 456, 82, 4, 53, "Section"],
Cell[15158, 462, 3140, 72, 489, "Text"],
Cell[18301, 536, 1566, 28, 384, "Input"],
Cell[CellGroupData[{
Cell[19892, 568, 62, 1, 30, "Input"],
Cell[19957, 571, 174, 3, 48, "Output"]
}, Open ]],
Cell[20146, 577, 151, 3, 33, "Text"],
Cell[CellGroupData[{
Cell[20322, 584, 178, 3, 63, "Input"],
Cell[20503, 589, 61, 1, 29, "Output"]
}, Open ]],
Cell[20579, 593, 916, 15, 261, "Text"],
Cell[21498, 610, 1049, 27, 90, "Text"],
Cell[22550, 639, 155, 4, 30, "Input"],
Cell[22708, 645, 295, 7, 52, "Text"],
Cell[23006, 654, 484, 14, 52, "Text"],
Cell[23493, 670, 501, 9, 166, "Input"],
Cell[23997, 681, 201, 4, 52, "Text"],
Cell[CellGroupData[{
Cell[24223, 689, 103, 2, 30, "Input"],
Cell[24329, 693, 64, 1, 29, "Output"]
}, Open ]],
Cell[24408, 697, 55, 0, 33, "Text"],
Cell[CellGroupData[{
Cell[24488, 701, 103, 2, 30, "Input"],
Cell[24594, 705, 65, 1, 29, "Output"]
}, Open ]],
Cell[24674, 709, 240, 8, 33, "Text"],
Cell[24917, 719, 695, 15, 192, "Input"],
Cell[CellGroupData[{
Cell[25637, 738, 143, 2, 42, "Input"],
Cell[25783, 742, 61, 1, 29, "Output"]
}, Open ]],
Cell[25859, 746, 371, 7, 52, "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell[26267, 758, 32, 0, 47, "Subsection"],
Cell[26302, 760, 669, 15, 109, "Text"]
}, Open ]],
Cell[CellGroupData[{
Cell[27008, 780, 38, 0, 47, "Subsection"],
Cell[27049, 782, 964, 27, 185, "Text"]
}, Open ]]
}
]
*)
(***********************************************************************
End of Mathematica Notebook file.
***********************************************************************)