%
% \iffalse
%
%<*pkg>
%% Copyright (c) 2010 by Lars Hellstrom. All rights reserved.
%% See the file "license.terms" for information on usage and redistribution
%% of this file, and for a DISCLAIMER OF ALL WARRANTIES.
%</pkg>
%<*driver>
\documentclass{tclldoc}
\usepackage{amsmath,amsfonts}
\usepackage{url}
\newcommand{\Tcl}{\Tcllogo}
\begin{document}
\DocInput{numtheory.dtx}
\end{document}
%</driver>
% \fi
%
% \title{Number theory package}
% \author{Lars Hellstr\"om}
% \date{30 May 2010}
% \maketitle
%
% \begin{abstract}
% This package provides a command to test whether an integer is a
% prime, but may in time come to house also other number-theoretic
% operations.
% \end{abstract}
%
% \tableofcontents
%
% \section*{Preliminaries}
%
% \begin{tcl}
%<*pkg>
package require Tcl 8.5
% \end{tcl}
% \Tcl~8.4 is seriously broken with respect to arithmetic overflow,
% so we require 8.5. There are (as yet) no explicit 8.5-isms in the
% code, however.
% \begin{tcl}
package provide math::numtheory 1.0
namespace eval ::math::numtheory {
namespace export isprime
}
%</pkg>
% \end{tcl}
% \setnamespace{math::numtheory}
%
% \Tcl lib has its own test file boilerplate.
% \begin{tcl}
%<*test>
source [file join\
[file dirname [file dirname [file join [pwd] [info script]]]]\
devtools testutilities.tcl]
testsNeedTcl 8.5
testsNeedTcltest 2
testing {useLocal numtheory.tcl math::numtheory}
%</test>
% \end{tcl}
%
% And the same is true for the manpage.
% \begin{tcl}
%<*man>
[manpage_begin math::numtheory n 1.0]
[copyright "2010 Lars Hellstr\u00F6m\
<Lars dot Hellstrom at residenset dot net>"]
[moddesc {Tcl Math Library}]
[titledesc {Number Theory}]
[category Mathematics]
[require Tcl [opt 8.5]]
[require math::numtheory [opt 1.0]]
[description]
[para]
This package is for collecting various number-theoretic operations,
though at the moment it only provides that of testing whether an
integer is a prime.
[list_begin definitions]
%</man>
% \end{tcl}
%
%
% \section{Primes}
%
% The first (and so far only) operation provided is |isprime|, which
% tests if an integer is a prime.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::isprime] [arg N] [
opt "[arg option] [arg value] ..."
]]
The [cmd isprime] command tests whether the integer [arg N] is a
prime, returning a boolean true value for prime [arg N] and a
boolean false value for non-prime [arg N]. The formal definition of
'prime' used is the conventional, that the number being tested is
greater than 1 and only has trivial divisors.
[para]
To be precise, the return value is one of [const 0] (if [arg N] is
definitely not a prime), [const 1] (if [arg N] is definitely a
prime), and [const on] (if [arg N] is probably prime); the latter
two are both boolean true values. The case that an integer may be
classified as "probably prime" arises because the Miller-Rabin
algorithm used in the test implementation is basically probabilistic,
and may if we are unlucky fail to detect that a number is in fact
composite. Options may be used to select the risk of such
"false positives" in the test. [const 1] is returned for "small"
[arg N] (which currently means [arg N] < 118670087467), where it is
known that no false positives are possible.
[para]
The only option currently defined is:
[list_begin options]
[opt_def -randommr [arg repetitions]]
which controls how many times the Miller-Rabin test should be
repeated with randomly chosen bases. Each repetition reduces the
probability of a false positive by a factor at least 4. The
default for [arg repetitions] is 4.
[list_end]
Unknown options are silently ignored.
%</man>
% \end{tcl}
%
%
% \subsection{Trial division}
%
% As most books on primes will tell you, practical primality
% testing algorithms typically start with trial division by a list
% of small known primes to weed out the low hanging fruit. This is
% also an opportunity to handle special cases that might arise for
% very low numbers (e.g.\ $2$ is a prime despite being even).
%
% \begin{proc}{prime_trialdivision}
% This procedure is meant to be called as
% \begin{quote}
% |prime_trialdivision| \word{$n$}
% \end{quote}
% from \emph{within} a procedure that returns |1| if $n$ is a prime
% and |0| if it is not. It does not return anything particular, but
% \emph{it causes its caller to return provided} it is able to
% decide what its result should be. This means one can slap it in
% as the first line of a primality checker procedure, and then on
% lines two and forth worry only about the nontrivial cases.
% \begin{tcl}
%<*pkg>
proc ::math::numtheory::prime_trialdivision {n} {
if {$n<2} then {return -code return 0}
% \end{tcl}
% Integers less than $2$ aren't primes.\footnote{
% Well, at least as one usually defines the term for integers.
% When considering the concept of prime in more general rings,
% one may have to settle with accepting all associates of primes
% as primes as well.
% } This saves us many worries by excluding negative numbers from
% further considerations.
% \begin{tcl}
if {$n<4} then {return -code return 1}
% \end{tcl}
% Everything else below \(2^2 = 4\) (i.e., $2$ and $3$) are primes.
% \begin{tcl}
if {$n%2 == 0} then {return -code return 0}
% \end{tcl}
% Remaining even numbers are then composite.
% \begin{tcl}
if {$n<9} then {return -code return 1}
% \end{tcl}
% Now everything left below \(3^2 = 9\) (i.e., $5$ and $7$) are
% primes. Having decided those, we can now do trial division with
% $3$, $5$, and $7$ in one go.
% \begin{tcl}
if {$n%3 == 0} then {return -code return 0}
if {$n%5 == 0} then {return -code return 0}
if {$n%7 == 0} then {return -code return 0}
% \end{tcl}
% Any numbers less that \(11^2 = 121\) not yet eliminated are
% primes; above that we know nothing.
% \begin{tcl}
if {$n<121} then {return -code return 1}
}
%</pkg>
% \end{tcl}
% This procedure could be extended with more primes, pushing the
% limit of what can be decided further up, but the returns are
% diminishing, so we might be better off with a different method
% for testing primality. No analysis of where the cut-off point
% lies have been conducted (i.e., $7$ as last prime for trial
% division was picked arbitrarily), but note that the optimum
% probably depends on what distribution the input values will have.
%
% \begin{tcl}
%<*test>
test prime_trialdivision-1 "Trial division of 1" -body {
::math::numtheory::prime_trialdivision 1
} -returnCodes 2 -result 0
test prime_trialdivision-2 "Trial division of 2" -body {
::math::numtheory::prime_trialdivision 2
} -returnCodes 2 -result 1
test prime_trialdivision-3 "Trial division of 6" -body {
::math::numtheory::prime_trialdivision 6
} -returnCodes 2 -result 0
test prime_trialdivision-4 "Trial division of 7" -body {
::math::numtheory::prime_trialdivision 7
} -returnCodes 2 -result 1
test prime_trialdivision-5 "Trial division of 101" -body {
::math::numtheory::prime_trialdivision 101
} -returnCodes 2 -result 1
test prime_trialdivision-6 "Trial division of 105" -body {
::math::numtheory::prime_trialdivision 105
} -returnCodes 2 -result 0
% \end{tcl}
% Note that extending the number of primes for trial division is
% likely to change the results in the following two tests ($121$
% is composite, $127$ is prime).
% \begin{tcl}
test prime_trialdivision-7 "Trial division of 121" -body {
::math::numtheory::prime_trialdivision 121
} -returnCodes 0 -result ""
test prime_trialdivision-8 "Trial division of 127" -body {
::math::numtheory::prime_trialdivision 127
} -returnCodes 0 -result ""
%</test>
% \end{tcl}
% \end{proc}
%
%
% \subsection{Pseudoprimality tests}
%
% After trial division, the next thing tried is usually to test the
% claim of Fermat's little theorem: if $n$ is a prime, then \(a^{n-1}
% \equiv 1 \pmod{n}\) for all integers $a$ that are not multiples of
% $n$, in particular those \(0 < a < n\); one picks such an $a$ (more
% or less at random) and computes $a^{n-1} \bmod n$. Numbers that
% pass are said to be \emph{(Fermat) pseudoprimes (to base $a$)}.
% Most composite numbers quickly fail this test.
% (One particular class that fails are the powers of primes, since
% the group of invertible elements in $\mathbb{Z}_n$ for \(n = p^{k+1}\)
% is cyclic\footnote{
% The easiest way to see that it is cyclic is probably to exhibit
% an element of order $(p -\nobreak 1) p^k$. A good start is to
% pick a primitive root $a$ of $\mathbb{Z}_p$ and compute its order
% modulo $p^{k+1}$; this has to be a number on the form $(p
% -\nobreak 1) p^r$. If \(r=k\) then $a$ is a primitive root and we're
% done, otherwise $(p +\nobreak 1) a$ will be a primitive root
% because $p+1$ can be shown to have order $p^k$ modulo $n$ and the
% least common multiple of $(p -\nobreak 1) p^r$ and $p^k$ is
% $(p -\nobreak 1) p^k$. To exhibit the order of $p+1$, one may
% use induction on $k$ to show that \( (1 +\nobreak p)^N \equiv 1
% \pmod{p^{k+1}}\) implies \(p^k \mid N\); in \((1 +\nobreak p)^N =
% \sum_{i=0}^N \binom{N}{i} p^i\), the induction hypothesis implies
% all terms with \(i>1\) vanish modulo $p^{k+1}$, leaving just
% \(1+Np \equiv 1 \pmod{p^{k+1}}\).
% } of order $(p -\nobreak 1) p^k$ rather than order $p^{k+1}-1$.
% Therefore it is only to bases $a$ of order dividing $p-1$ (i.e., a
% total of $p-1$ out of $p^{k+1}-1$) that prime powers are
% pseudoprimes. The chances of picking one of these are generally
% rather slim.)
%
% Unfortunately, there are also numbers (the so-called \emph{Carmichael
% numbers}) which are pseudoprimes to every base $a$ they are coprime
% with. While the above trial division by $2$, $3$, $5$, and $7$ would
% already have eliminated all Carmichael numbers below \(29341 = 13
% \cdot 37 \cdot 61\), their existence means that there is a
% population of nonprimes which a Fermat pseudoprimality test is very
% likely to mistake for primes; this would usually not be acceptable.
%
% \begin{proc}{Miller--Rabin}
% The Miller--Rabin test is a slight variation on the Fermat test,
% where the computation of $a^{n-1} \bmod n$ is structured so that
% additional consequences of $n$ being a prime can be tested.
% Rabin~\cite{Rabin}
% proved that any composite $n$ will for this test be revealed as
% such by at least $3/4$ of all bases $a$, thus making it a valid
% probabilistic test. (Miller~\cite{Miller} had first designed it as
% a deterministic polynomial algorithm, but in that case the proof
% that it works relies on the generalised Riemann hypothesis.)
%
% Given natural numbers $s$ and $d$ such that \(n-1 = 2^s d\), the
% computation of $a^{n-1}$ is organised as $(a^d)^{2^s}$, where the
% $s$ part is conveniently performed by squaring $s$ times. This is
% of little consequence when $n$ is not a pseudoprime since one
% will simply arrive at some \(a^{n-1} \not\equiv 1 \pmod{n}\), but
% in the case that $n$ is a pseudoprime these repeated squarings will
% exhibit some $x$ such that \(x^2 \equiv 1 \pmod{n}\), and this
% makes it possible to test another property $n$ must have if it is
% prime, namely that such an \(x \equiv \pm 1 \pmod{n}\).
%
% That implication is of course well known for real (and complex)
% numbers, but even though what we're dealing with here is rather
% residue classes modulo an integer, the proof that it holds is
% basically the same. If $n$ is a prime, then the residue class
% ring $\mathbb{Z}_n$ is a field, and hence the ring
% $\mathbb{Z}_n[x]$ of polynomials over that field is a Unique
% Factorisation Domain. As it happens, \(x^2 \equiv 1 \pmod{n}\) is
% a polynomial equation, and $x^2-1$ has the known factorisation
% \((x -\nobreak 1) (x +\nobreak 1)\). Since factorisations are
% unique, and any zero $a$ of $x^2-1$ would give rise to a factor
% $x-a$, it follows that \(x^2 \equiv 1 \pmod{n}\) implies \(x
% \equiv 1 \pmod{n}\) or \(x \equiv -1 \pmod{n}\), as claimed.
% But this assumes $n$ is a prime.
%
% If instead \(n = pq\) where \(p,q > 2\) are coprime, then there
% will be additional solutions to \(x^2 \equiv 1 \pmod{n}\).
% For example, if \(x \equiv 1 \pmod{p}\) and \(x \equiv -1
% \pmod{q}\) (and such $x$ exist by the Chinese Remainder Theorem),
% then \(x^2 \equiv 1 \pmod{p}\) and \(x^2 \equiv 1 \pmod{q}\),
% from which follows \(x^2 \equiv 1 \pmod{pq}\), but \(x \not\equiv
% 1 \pmod{n}\) since \(x-1 \equiv -2 \not\equiv 0 \pmod{q}\), and
% \(x \not\equiv -1 \pmod{n}\) since \(x+1 \equiv 2 \not\equiv 0
% \pmod{p}\). The same argument applies when \(x \equiv -1 \pmod{p}\)
% and \(x \equiv 1 \pmod{q}\), and in general, if $n$ has $k$
% distinct odd prime factors then one may construct $2^k$ distinct
% solutions \(0<x<n\) to \(x^2 \equiv 1 \pmod{n}\). It is thus not
% too hard to imagine that a ``random'' $a^d$ squaring to $1$
% modulo $n$ will be one of the nonstandard square roots of~$1$
% when $n$ is not a prime, even if the above is not a proof that
% at least $3/4$ of all $a$ are witnesses to the compositeness
% of~$n$.
%
% Getting down to the implementation, the actual procedure has the
% call syntax
% \begin{quote}
% |Miller--Rabin| \word{n} \word{s} \word{d} \word{a}
% \end{quote}
% where all arguments should be integers such that \(n-1 = d2^s\),
% \(d,s \geq 1\), and \(0 < a < n\). The procedure computes
% $(a^d)^{2^s} \mod n$, and if in the course of doing this the
% Miller--Rabin test detects that $n$ is composite then this procedure
% will return |1|, otherwise it returns |0|.
%
% The first part of the procedure merely computes \(x = a^d \bmod n\),
% using exponentiation by squaring. $x$, $a$, and $d$ are modified in
% the loop, but $xa^d \bmod n$ would be an invariant quantity.
% Correctness presumes the initial \(d \geq 1\).
% \begin{tcl}
%<*pkg>
proc ::math::numtheory::Miller--Rabin {n s d a} {
set x 1
while {$d>1} {
if {$d & 1} then {set x [expr {$x*$a % $n}]}
set a [expr {$a*$a % $n}]
set d [expr {$d >> 1}]
}
set x [expr {$x*$a % $n}]
% \end{tcl}
% The second part will $s-1$ times square $x$, while checking each
% value for being \(\equiv \pm 1 \pmod{n}\). For most part, $-1$
% means everything is OK (any subsequent square would only
% yield~$1$) whereas $1$ arrived at without a previous $-1$ signals
% that $n$ cannot be prime. The only exception to the latter is
% that $1$ before the first squaring (already \(a^d \equiv 1
% \pmod{n}\)) is OK as well.
% \begin{tcl}
if {$x == 1} then {return 0}
for {} {$s>1} {incr s -1} {
if {$x == $n-1} then {return 0}
set x [expr {$x*$x % $n}]
if {$x == 1} then {return 1}
}
% \end{tcl}
% There is no need to square $x$ the $s$th time, because if at this
% point \(x \not\equiv -1 \pmod{n}\) then $n$ cannot be a prime; if
% \(x^2 \not\equiv 1 \pmod{n}\) it would fail to be a pseudoprime
% and if \(x^2 \equiv 1 \pmod{n}\) then $x$ would be a nonstandard
% square root of $1 \pmod{n}$, but it is not necessary to find out
% which of these cases is at hand.
% \begin{tcl}
return [expr {$x != $n-1}]
}
%</pkg>
% \end{tcl}
%
% As for testing, the minimal allowed value of $n$ is $3$, which
% is a prime.
% \begin{tcl}
%<*test>
test Miller--Rabin-1.1 "Miller--Rabin 3" -body {
list [::math::numtheory::Miller--Rabin 3 1 1 1]\
[::math::numtheory::Miller--Rabin 3 1 1 2]
} -result {0 0}
% \end{tcl}
% To exercise the first part of the procedure, one may consider the
% case \(s=1\) and \(d = 2^2+2^0 = 5\), i.e., \(n=11\). Here, \(2^5
% \equiv -1 \pmod{11}\) whereas \(4^5 \equiv 1^5 \equiv 1
% \pmod{11}\). A bug on the lines of not using the right factors in
% the computation of $a^d$ would most likely end up with something
% different here.
% \begin{tcl}
test Miller--Rabin-1.2 "Miller--Rabin 11" -body {
list [::math::numtheory::Miller--Rabin 11 1 5 1]\
[::math::numtheory::Miller--Rabin 11 1 5 2]\
[::math::numtheory::Miller--Rabin 11 1 5 4]
} -result {0 0 0}
% \end{tcl}
% $27$ will on the other hand be exposed as composite by most bases,
% but $1$ and $-1$ do not spot it. It is known from the argument
% about prime powers above that at least one of $2$ and \(8 = (3
% +\nobreak 1) \cdot 2\) is a primitive root of $1$ in
% $\mathbb{Z}_{27}$; it turns out to be $2$.
% \begin{tcl}
test Miller--Rabin-1.3 "Miller--Rabin 27" -body {
list [::math::numtheory::Miller--Rabin 27 1 13 1]\
[::math::numtheory::Miller--Rabin 27 1 13 2]\
[::math::numtheory::Miller--Rabin 27 1 13 3]\
[::math::numtheory::Miller--Rabin 27 1 13 4]\
[::math::numtheory::Miller--Rabin 27 1 13 8]\
[::math::numtheory::Miller--Rabin 27 1 13 26]
} -result {0 1 1 1 1 0}
% \end{tcl}
% Taking \(n = 65 = 1 + 2^6 = 5 \cdot 13\) instead focuses on the
% second part of the procedure. By carefully choosing the base, it
% is possible to force the result to come from:
% \begin{tcl}
test Miller--Rabin-1.4 "Miller--Rabin 65" -body {
% \end{tcl}
% The first |return|
% \begin{tcl}
list [::math::numtheory::Miller--Rabin 65 6 1 1]\
% \end{tcl}
% the second |return|, first iteration
% \begin{tcl}
[::math::numtheory::Miller--Rabin 65 6 1 64]\
% \end{tcl}
% the third |return|, first iteration---\(14 \equiv 1 \pmod{13}\)
% but \(14 \equiv -1 \pmod{5}\)
% \begin{tcl}
[::math::numtheory::Miller--Rabin 65 6 1 14]\
% \end{tcl}
% the second |return|, second iteration
% \begin{tcl}
[::math::numtheory::Miller--Rabin 65 6 1 8]\
% \end{tcl}
% the third |return|, second iteration---\(27 \equiv 1 \pmod{13}\)
% but \(27^2 \equiv 2^2 \equiv -1 \pmod{5}\)
% \begin{tcl}
[::math::numtheory::Miller--Rabin 65 6 1 27]\
% \end{tcl}
% the final |return|
% \begin{tcl}
[::math::numtheory::Miller--Rabin 65 6 1 2]
} -result {0 0 1 0 1 1}
% \end{tcl}
% There does however not seem to be any \(n=65\) choice of $a$ which
% would get a |0| out of the final |return|.
%
% An $n$ which allows fully exercising the second part of the
% procedure is \(17 \cdot 257 = 4369\), for which \(s=4\)
% and \(d=273\). In order to have \(x^{2^{s-1}} \equiv -1
% \pmod{n}\), it is necessary to have \(x^8 \equiv -1\) modulo both
% $17$ and $257$, which is possible since the invertible elements
% of $\mathbb{Z}_{17}$ form a cyclic group of order $16$ and the
% invertible elements of $\mathbb{Z}_{257}$ form a cyclic group of
% order $256$. Modulo $17$, an element of order $16$ is $3$,
% whereas modulo $257$, an element of order $16$ is $2$.
%
% There is an extra complication in that what the caller can
% specify is not the $x$ to be repeatedly squared, but the $a$
% which satisfies \(x \equiv a^d \pmod{n}\). Since \(d=273\) is
% odd, raising something to that power is an invertible operation
% modulo both $17$ and $257$, but it is necessary to figure out
% what the inverse is. Since \(273 \equiv 1 \pmod{16}\), it turns
% out that \(a^d \equiv a \pmod{17}\), and \(x=3\) becomes \(a=3\).
% From \(273 \equiv 17 \pmod{256}\), it instead follows that \(x
% \equiv a^d \pmod{257}\) is equivalent to \(a \equiv x^e
% \pmod{257}\), where \(17e \equiv 1 \pmod{256}\). This has the
% solution \(e = 241\), so the $a$ which makes \(x=2\) is \(a
% = 2^{241} \bmod 257\). However, since \(x=2\) was picked on
% account of having order $16$, hence \(2^{16} \equiv 1
% \pmod{257}\), and \(241 \equiv 1 \pmod{16}\), it again turns out
% that \(x=2\) becomes \(a=2\).
%
% For \(a = 2\), one may observe that \(a^{2^1} \equiv 4
% \pmod{257}\), \(a^{2^2} \equiv 16 \pmod{257}\), \(a^{2^3} \equiv
% -1 \pmod{257}\), and \(a^{2^4} \equiv 1 \pmod{257}\). For
% \(a=3\), one may observe that \(a^{2^1} \equiv 9 \pmod{17}\),
% \(a^{2^2} \equiv 13 \pmod{17}\), \(a^{2^3} \equiv -1 \pmod{17}\),
% and \(a^{2^4} \equiv 1 \pmod{17}\). For solving simultaneous
% equivalences, it is furthermore useful to observe that \(2057
% \equiv 1 \pmod{257}\) and \(2057 \equiv 0 \pmod{17}\) whereas
% \(2313 \equiv 1 \pmod{17}\) and \(2313 \equiv 0 \pmod{257}\).
% \begin{tcl}
test Miller--Rabin-1.5 "Miller--Rabin 17*257" -body {
% \end{tcl}
% In order to end up at the first |return|, it is necessary to take
% \(a \equiv 1 \pmod{17}\) and \(a \equiv 1 \pmod{257}\); the
% solution \(a=1\) is pretty obvious.
% \begin{tcl}
list [::math::numtheory::Miller--Rabin 4369 4 273 1]\
% \end{tcl}
% In order to end up at the second |return| of the first iteration,
% it is necessary to take \(a \equiv -1 \pmod{17}\) and
% \(a \equiv -1 \pmod{257}\); the solution \(a \equiv -1 \pmod{n}\)
% is again pretty obvious.
% \begin{tcl}
[::math::numtheory::Miller--Rabin 4369 4 273 4368]\
% \end{tcl}
% Hitting the third |return| at the first iteration can be achieved
% with \(a \equiv -1 \pmod{17}\) and \(a \equiv 1 \pmod{257}\);
% now a solution is \(a \equiv 2057 - 2313 \equiv 4113 \pmod{n}\).
% \begin{tcl}
[::math::numtheory::Miller--Rabin 4369 4 273 4113]\
% \end{tcl}
% Hitting the second |return| at the second iteration happens if
% \(a^2 \equiv -1\) modulo both prime factors, i.e., for \(a \equiv
% 16 \pmod{257}\) and \(a \equiv 13 \pmod{17}\). This has the
% solution \(a \equiv 16 \cdot 2057 + 13 \cdot 2313 \equiv 1815
% \pmod{n}\).
% \begin{tcl}
[::math::numtheory::Miller--Rabin 4369 4 273 1815]\
% \end{tcl}
% To hit the third |return| at the second iteration, one may keep
% \(a \equiv 16 \pmod{257}\) but take \(a \equiv 1 \pmod{17}\). This
% has the solution \(a \equiv 16 \cdot 2057 + 1 \cdot 2313 \equiv 273
% \pmod{n}\).
% \begin{tcl}
[::math::numtheory::Miller--Rabin 4369 4 273 273]\
% \end{tcl}
% Hitting the second |return| at the third and final iteration happens
% if \(a^4 \equiv -1\) modulo both prime factors, i.e., for \(a \equiv
% 4 \pmod{257}\) and \(a \equiv 9 \pmod{17}\). This has the
% solution \(a \equiv 4 \cdot 2057 + 9 \cdot 2313 \equiv 2831
% \pmod{n}\).
% \begin{tcl}
[::math::numtheory::Miller--Rabin 4369 4 273 2831]\
% \end{tcl}
% And as before, to hit the third |return| at the third and final
% iteration one may keep the above \(a \equiv 9 \pmod{17}\) but
% change the other to \(a \equiv 1 \pmod{257}\). This has the
% solution \(a \equiv 1 \cdot 2057 + 9 \cdot 2313 \equiv 1029
% \pmod{n}\).
% \begin{tcl}
[::math::numtheory::Miller--Rabin 4369 4 273 1029]\
% \end{tcl}
% To get a |0| out of the fourth |return|, one takes \(a \equiv
% 2 \pmod{257}\) and \(a \equiv 3 \pmod{17}\); this means \(a \equiv
% 2 \cdot 2057 + 3 \cdot 2313 \equiv 2315 \pmod{n}\).
% \begin{tcl}
[::math::numtheory::Miller--Rabin 4369 4 273 2315]\
% \end{tcl}
% Finally, to get a |1| out of the fourth |return|, one may take
% \(a \equiv 1 \pmod{257}\) and \(a \equiv 3 \pmod{17}\); this means
% \(a \equiv 1 \cdot 2057 + 3 \cdot 2313 \equiv 258 \pmod{n}\).
% \begin{tcl}
[::math::numtheory::Miller--Rabin 4369 4 273 258]
} -result {0 0 1 0 1 0 1 0 1}
% \end{tcl}
% It would have been desirable from a testing point of view to also
% find a value of $a$ that would make \(a^{n-1} \equiv -1
% \pmod{n}\), since such an $a$ would catch an implementation error
% of running the squaring loop one step too far, but that does not
% seem possible; picking \(n=pq\) such that both $p-1$ and $q-1$
% are divisible by some power of $2$ implies that $n-1$ is
% divisible by the same power of $2$.
% \end{proc}
%
% A different kind of test is to verify some exceptional numbers and
% boundaries that the |isprime| procedure relies on. First, $1373653$
% appears prime when \(a=2\) or \(a=3\), but \(a=5\) is a witness to
% its compositeness.
% \begin{tcl}
test Miller--Rabin-2.1 "Miller--Rabin 1373653" -body {
list\
[::math::numtheory::Miller--Rabin 1373653 2 343413 2]\
[::math::numtheory::Miller--Rabin 1373653 2 343413 3]\
[::math::numtheory::Miller--Rabin 1373653 2 343413 5]
} -result {0 0 1}
% \end{tcl}
% $25326001$ is looking like a prime also to \(a=5\), but \(a=7\)
% exposes it.
% \begin{tcl}
test Miller--Rabin-2.2 "Miller--Rabin 25326001" -body {
list\
[::math::numtheory::Miller--Rabin 25326001 4 1582875 2]\
[::math::numtheory::Miller--Rabin 25326001 4 1582875 3]\
[::math::numtheory::Miller--Rabin 25326001 4 1582875 5]\
[::math::numtheory::Miller--Rabin 25326001 4 1582875 7]
} -result {0 0 0 1}
% \end{tcl}
% $3215031751$ is a tricky composite that isn't exposed even by
% \(a=7\), but \(a=11\) will see through it.
% \begin{tcl}
test Miller--Rabin-2.3 "Miller--Rabin 3215031751" -body {
list\
[::math::numtheory::Miller--Rabin 3215031751 1 1607515875 2]\
[::math::numtheory::Miller--Rabin 3215031751 1 1607515875 3]\
[::math::numtheory::Miller--Rabin 3215031751 1 1607515875 5]\
[::math::numtheory::Miller--Rabin 3215031751 1 1607515875 7]\
[::math::numtheory::Miller--Rabin 3215031751 1 1607515875 11]
} -result {0 0 0 0 1}
% \end{tcl}
% Otherwise the lowest composite that these four will fail for is
% $118670087467$.
% \begin{tcl}
test Miller--Rabin-2.4 "Miller--Rabin 118670087467" -body {
list\
[::math::numtheory::Miller--Rabin 118670087467 1 59335043733 2]\
[::math::numtheory::Miller--Rabin 118670087467 1 59335043733 3]\
[::math::numtheory::Miller--Rabin 118670087467 1 59335043733 5]\
[::math::numtheory::Miller--Rabin 118670087467 1 59335043733 7]\
[::math::numtheory::Miller--Rabin 118670087467 1 59335043733 11]
} -result {0 0 0 0 1}
%</test>
% \end{tcl}
%
%
% \subsection{Putting it all together}
%
% \begin{proc}{isprime}
% The user level command for testing primality of an integer $n$ is
% |isprime|. It has the call syntax
% \begin{quote}
% |math::numtheory::isprime| \word{n}
% \begin{regblock}[\regstar]\word{option}
% \word{value}\end{regblock}
% \end{quote}
% where the options may be used to influence the exact algorithm
% being used. The call returns
% \begin{description}
% \item[0] if $n$ is found to be composite,
% \item[1] if $n$ is found to be prime, and
% \item[on] if $n$ is probably prime.
% \end{description}
% The reason there might be \emph{some} uncertainty is that the
% primality test used is basically a probabilistic test for
% compositeness---it may fail to find a witness for the
% compositeness of a composite number $n$, even if the probability
% of doing so is fairly low---and to be honest with the user, the
% outcomes of ``definitely prime'' and ``probably prime'' return
% different results. Since |on| is true when used as a boolean, you
% usually need not worry about this fine detail. Also, for \(n <
% 10^{11}\) (actually a little more) the primality test is
% deterministic, so you only encounter the ``probably prime''
% result for fairly high $n$.
%
% At present, the only option that is implemented is |-randommr|,
% which controls how many rounds (by default 4) of the |Miller--Rabin|
% test with random bases are run before returing |on|. Other options
% are silently ignored.
%
% \begin{tcl}
%<*pkg>
proc ::math::numtheory::isprime {n args} {
prime_trialdivision $n
% \end{tcl}
% Implementation-wise, |isprime| begins with |prime_trialdivision|,
% but relies on the Miller--Rabin test after that. To that end, it
% must compute $s$ and $d$ such that \(n = d 2^s + 1\); while this
% is fairly quick, it's nice not having to do it more than once,
% which is why this step wasn't made part of the |Miller--Rabin|
% procedure.
% \begin{tcl}
set d [expr {$n-1}]; set s 0
while {($d&1) == 0} {
incr s
set d [expr {$d>>1}]
}
% \end{tcl}
% The deterministic sequence of Miller--Rabin tests combines
% information from \cite{PSW80,Jaeschke}, but most of these
% numbers may also be found on Wikipedia~\cite{Wikipedia}.
% \begin{tcl}
if {[Miller--Rabin $n $s $d 2]} then {return 0}
if {$n < 2047} then {return 1}
if {[Miller--Rabin $n $s $d 3]} then {return 0}
if {$n < 1373653} then {return 1}
if {[Miller--Rabin $n $s $d 5]} then {return 0}
if {$n < 25326001} then {return 1}
if {[Miller--Rabin $n $s $d 7] || $n==3215031751} then {return 0}
if {$n < 118670087467} then {return 1}
% \end{tcl}
% \(3215031751 = 151 \cdot 751 \cdot 28351\) is a Carmichael
% number~\cite[p.\,1022]{PSW80}.
%
% Having exhausted this list of limits below which |Miller--Rabin|
% for \(a=2,3,5,7\) detects all composite numbers, we now have to
% resort to picking bases at random and hoping we find one which
% would reveal a composite $n$. In the future, one might want to
% add the possibility of using a deterministic test (such as the
% AKR~\cite{CL84} or AKS~\cite{AKS04} test) here instead.
%
% \begin{tcl}
array set O {-randommr 4}
array set O $args
for {set i $O(-randommr)} {$i >= 1} {incr i -1} {
if {[Miller--Rabin $n $s $d [expr {(
% \end{tcl}
%
% The probabilistic sequence of Miller--Rabin tests employs
% \Tcl's built-in pseudorandom number generator |rand()| for
% choosing bases, as this does not seem to be an application that
% requires high quality randomness. It may however be observed
% that since by now \(n > 10^{11}\), the space of possible bases $a$
% is always several times larger than the state space of |rand()|,
% so there may be a point in tweaking the PRNG to avoid some less
% useful values of $a$.
%
% It is a trivial observation that the intermediate $x$ values
% computed by |Miller--Rabin| for \(a=a_1a_2\) are simply the
% products of the corresponding values computed for \(a=a_1\) and
% \(a=a_2\) respectively---hence chances are that if no
% compositeness was detected for \(a=a_1\) or \(a=a_2\) then it
% won't be detected for \(a=a_1a_2\) either. There is a slight
% chance that something interesting could happen if \(a_1^{d2^k}
% \equiv -1 \equiv a_2^{d2^k} \pmod{n}\) for some \(k>0\), since in
% that case \((a_1a_2)^{d2^k} \equiv 1 \pmod{n}\) whereas no direct
% conclusion can be reached about $(a_1a_2)^{d2^{k-1}}$, but this
% seems a rather special case (and cannot even occur if \(n
% \equiv 3 \pmod{4}\) since in that case \(s=1\)), so it seems
% natural to prefer $a$ that are primes. Generating only prime $a$
% would be much work, but avoiding numbers divisible by $2$ or $3$
% is feasible.
%
% First turn |rand()| back into the integer it internally is, and
% adjust it to be from $0$ and up.
% \begin{tcl}
(round(rand()*0x100000000)-1)
% \end{tcl}
% Then multiply by $3$ and set the last bit. This has the effect
% that the range of the PRNG is now $\{1,3,7,9,13,15,\dotsb,
% 6n +\nobreak 1, 6n +\nobreak 3, \dotsb \}$.
% \begin{tcl}
*3 | 1)
% \end{tcl}
% Finally add $10$ so that we get $11$, $13$, $17$, $19$, \dots
% \begin{tcl}
+ 10
}]]} then {return 0}
}
% \end{tcl}
% That ends the |for| loop for |Miller--Rabin| with random bases.
% At this point, since the number in question passed the requested
% number of Miller--Rabin rounds, it is proclaimed to be ``probably
% prime''.
% \begin{tcl}
return on
}
%</pkg>
% \end{tcl}
%
% Tests of |isprime| would mostly be asking ``is $n$ a prime'' for
% various interesting $n$. Several values of $n$ should be the same
% as the previous tests:
% \begin{tcl}
%<*test>
test isprime-1.1 "1 is not prime" -body {
::math::numtheory::isprime 1
} -result 0
test isprime-1.2 "0 is not prime" -body {
::math::numtheory::isprime 0
} -result 0
test isprime-1.3 "-2 is not prime" -body {
::math::numtheory::isprime -2
} -result 0
test isprime-1.4 "2 is prime" -body {
::math::numtheory::isprime 2
} -result 1
test isprime-1.5 "6 is not prime" -body {
::math::numtheory::isprime 6
} -result 0
test isprime-1.6 "7 is prime" -body {
::math::numtheory::isprime 7
} -result 1
test isprime-1.7 "101 is prime" -body {
::math::numtheory::isprime 101
} -result 1
test isprime-1.8 "105 is not prime" -body {
::math::numtheory::isprime 105
} -result 0
test isprime-1.9 "121 is not prime" -body {
::math::numtheory::isprime 121
} -result 0
test isprime-1.10 "127 is prime" -body {
::math::numtheory::isprime 127
} -result 1
test isprime-1.11 "4369 is not prime" -body {
::math::numtheory::isprime 4369
} -result 0
test isprime-1.12 "1373653 is not prime" -body {
::math::numtheory::isprime 1373653
} -result 0
test isprime-1.13 "25326001 is not prime" -body {
::math::numtheory::isprime 25326001
} -result 0
test isprime-1.14 "3215031751 is not prime" -body {
::math::numtheory::isprime 3215031751
} -result 0
% \end{tcl}
% To get consistent results for large non-primes, it is necessary
% to reduce the number of random rounds and\slash or reset the PRNG.
% \begin{tcl}
test isprime-1.15 "118670087467 may appear prime, but isn't" -body {
expr srand(1)
list\
[::math::numtheory::isprime 118670087467 -randommr 0]\
[::math::numtheory::isprime 118670087467 -randommr 1]
} -result {on 0}
% \end{tcl}
% However, a few new can be added. On~\cite[p.\,925]{Jaeschke} we
% can read that \(p=22 \mkern1mu 754 \mkern1mu 930 \mkern1mu 352
% \mkern1mu 733\) is a prime, and $p (3p -\nobreak 2)\) is a
% composite number that looks prime to |Miller--Rabin| for all
% \(a \in \{2,3,5,7,11,13,17,19,23,29\}\).
% \begin{tcl}
test isprime-1.16 "Jaeschke psi_10" -body {
expr srand(1)
set p 22754930352733
set n [expr {$p * (3*$p-2)}]
list\
[::math::numtheory::isprime $p -randommr 25]\
[::math::numtheory::isprime $n -randommr 0]\
[::math::numtheory::isprime $n -randommr 1]
} -result {on on 0}
% \end{tcl}
% On the same page it is stated that \(p=137 \mkern1mu 716 \mkern1mu
% 125 \mkern1mu 329 \mkern1mu 053\) is a prime such that
% $p (3p -\nobreak 2)\) is a composite number that looks prime to
% |Miller--Rabin| for all \(a \in
% \{2,3,5,7,11,13,17,19,23,29,31\}\).
% \begin{tcl}
test isprime-1.17 "Jaeschke psi_11" -body {
expr srand(1)
set p 137716125329053
set n [expr {$p * (3*$p-2)}]
list\
[::math::numtheory::isprime $p -randommr 25]\
[::math::numtheory::isprime $n -randommr 0]\
[::math::numtheory::isprime $n -randommr 1]\
[::math::numtheory::isprime $n -randommr 2]
} -result {on on on 0}
% \end{tcl}
% RFC~2409~\cite{RFC2409} lists a number of primes (and primitive
% generators of their respective multiplicative groups). The
% smallest of these is defined as \(p = 2^{768} - 2^{704} - 1 +
% 2^{64} \cdot \left( [2^{638} \pi] + 149686 \right)\) (where the
% brackets probably denote rounding to the nearest integer), but
% since high precision (roughly $200$ decimal digits would be
% required) values of \(\pi = 3.14159\dots\) are a bit awkward to
% get hold of, we might as well use the stated hexadecimal digit
% expansion for~$p$. It might also be a good idea to verify that
% this is given with most significant digit first.
% \begin{tcl}
test isprime-1.18 "OAKLEY group 1 prime" -body {
set digits [join {
FFFFFFFF FFFFFFFF C90FDAA2 2168C234 C4C6628B 80DC1CD1
29024E08 8A67CC74 020BBEA6 3B139B22 514A0879 8E3404DD
EF9519B3 CD3A431B 302B0A6D F25F1437 4FE1356D 6D51C245
E485B576 625E7EC6 F44C42E9 A63A3620 FFFFFFFF FFFFFFFF
} ""]
expr srand(1)
list\
[::math::numtheory::isprime 0x$digits]\
[::math::numtheory::isprime 0x[string reverse $digits]]
} -result {on 0}
% \end{tcl}
%
% A quite different thing to test is that the tweaked PRNG really
% produces only \(a \equiv 1,5 \pmod{6}\).
% \begin{tcl}
test isprime-2.0 "PRNG tweak" -setup {
namespace eval ::math::numtheory {
rename Miller--Rabin _orig_Miller--Rabin
proc Miller--Rabin {n s d a} {
expr {$a>7 && $a%6!=1 && $a%6!=5}
}
}
} -body {
::math::numtheory::isprime 118670087467 -randommr 500
} -result on -cleanup {
namespace eval ::math::numtheory {
rename Miller--Rabin ""
rename _orig_Miller--Rabin Miller--Rabin
}
}
%</test>
% \end{tcl}
% \end{proc}
%
%
% \section*{Closings}
%
% \begin{tcl}
%<*man>
[list_end]
[keywords {number theory} prime]
[manpage_end]
%</man>
% \end{tcl}
%
% \begin{tcl}
%<test>testsuiteCleanup
% \end{tcl}
%
%
% \begin{thebibliography}{9}
%
% \bibitem{AKS04}
% Manindra Agrawal, Neeraj Kayal, and Nitin Saxena:
% PRIMES is in P,
% \textit{Annals of Mathematics} \textbf{160} (2004), no. 2,
% 781--793.
%
% \bibitem{CL84}
% Henri Cohen and Hendrik W. Lenstra, Jr.:
% Primality testing and Jacobi sums,
% \textit{Mathematics of Computation} \textbf{42} (165) (1984),
% 297--330.
% \texttt{doi:10.2307/2007581}
%
% \bibitem{RFC2409}
% Dan Harkins and Dave Carrel.
% \textit{The Internet Key Exchange (IKE)},
% \textbf{RFC 2409} (1998).
%
% \bibitem{Jaeschke}
% Gerhard Jaeschke: On strong pseudoprimes to several bases,
% \textit{Mathematics of Computation} \textbf{61} (204), 1993,
% 915--926.
% \texttt{doi:\,10.2307/2153262}
%
% \bibitem{Miller}
% Gary L. Miller:
% Riemann's Hypothesis and Tests for Primality,
% \textit{Journal of Computer and System Sciences} \textbf{13} (3)
% (1976), 300--317. \texttt{doi:10.1145/800116.803773}
%
% \bibitem{PSW80}
% C.~Pomerance, J.~L.~Selfridge, and S.~S.~Wagstaff~Jr.:
% The pseudoprimes to $25 \cdot 10^9$,
% \textit{Mathematics of Computation} \textbf{35} (151), 1980,
% 1003--1026.
% \texttt{doi: 10.2307/2006210}
%
% \bibitem{Rabin}
% Michael O. Rabin:
% Probabilistic algorithm for testing primality,
% \textit{Journal of Number Theory} \textbf{12} (1) (1980),
% 128--138. \texttt{doi:10.1016/0022-314X(80)90084-0}
%
% \bibitem{Wikipedia}
% Wikipedia contributors:
% Miller--Rabin primality test,
% \textit{Wikipedia, The Free Encyclopedia}, 2010.
% Online, accessed 10 September 2010,
% \url{http://en.wikipedia.org/w/index.php?title=Miller%E2%80%93Rabin_primality_test&oldid=383901104}
%
% \end{thebibliography}
%
\endinput
