% 
% \iffalse
% 
%<*pkg_common>
%% 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_common>
%<*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.1.1
namespace eval ::math::numtheory {
   namespace export isprime
}
%</pkg>
% \end{tcl}
% Additional procedures are placed into a separate file primes.tcl,
% sourced from the primary.
% \begin{tcl}
%<*pkg_primes>
# primes.tcl --
#     Provide additional procedures for the number theory package
#
namespace eval ::math::numtheory {
    variable primes {2 3 5 7 11 13 17}
    variable nextPrimeCandidate 19
    variable nextPrimeIncrement  1 ;# Examine numbers 6n+1 and 6n+5

    namespace export firstNprimes primesLowerThan primeFactors uniquePrimeFactors factors \
                     totient moebius legendre jacobi gcd lcm \
                     numberPrimesGauss numberPrimesLegendre numberPrimesLegendreModified
}

%</pkg_primes>
% \end{tcl}
% \setnamespace{math::numtheory}
% 
% \Tcl lib has its own test file boilerplate.
% We require tcltest 2.1 to allow the definition of a custom matcher
% comparing lists of integers
% \begin{tcl}
%<*test_primes>
# -*- tcl -*-
# primes.test --
#    Additional test cases for the ::math::numtheory package
#
# Note:
#    The tests assume tcltest 2.1, in order to compare
#    list of integer results

# -------------------------------------------------------------------------

%</test_primes>
%<*test_common>
source [file join\
  [file dirname [file dirname [file join [pwd] [info script]]]]\
  devtools testutilities.tcl]

testsNeedTcl     8.5
testsNeedTcltest 2.1

%</test_common>
%<*test_primes>
support {
    useLocal math.tcl math
}

%</test_primes>
%<*test_common>
testing {
    useLocal numtheory.tcl math::numtheory
}

%</test_common>
% \end{tcl}
% and a bit more for the additional tests. This is where tcltest 2.1
% is required
% \begin{tcl}
%<*test_primes>
proc matchLists { expected actual } {
   set match 1
   foreach a $actual e $expected {
      if { $a != $e } {
         set match 0
         break
      }
   }
   return $match
}
customMatch equalLists matchLists

%</test_primes>
% \end{tcl}
% 
% And the same is true for the manpage.
% \begin{tcl}
%<*man>
[comment {
    __Attention__ This document is a generated file.
    It is not the true source.
    The true source is

        numtheory.dtx

    To make changes edit the true source, and then use
    
        sak.tcl docstrip/regen modules/math

    to update all generated files.
}]
[vset VERSION 1.1.1]
[manpage_begin math::numtheory n [vset VERSION]]
[keywords {number theory}]
[keywords prime]
[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 [vset VERSION]]]

[description]
[para]
This package is for collecting various number-theoretic operations, with
a slight bias to prime numbers.

[list_begin definitions]
%</man>
% \end{tcl}
% 
% 
% \section{Primes}
% 
% The first 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}
% Then we have |firstNprimes|, which returns a list containing
% the first |n| primes.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::firstNprimes] [arg N]]
Return the first N primes

[list_begin arguments]
[arg_def integer N in]
Number of primes to return
[list_end]

[call [cmd math::numtheory::primesLowerThan] [arg N]]
Return the prime numbers lower/equal to N

[list_begin arguments]
[arg_def integer N in]
Maximum number to consider
[list_end]

[call [cmd math::numtheory::primeFactors] [arg N]]
Return a list of the prime numbers in the number N

[list_begin arguments]
[arg_def integer N in]
Number to be factorised
[list_end]

%</man>
% \end{tcl}
% Similarly |primesLowerThan| returns a list of the prime numbers
% which are less than |n|, or equal to it.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::primesLowerThan] [arg N]]
Return the prime numbers lower/equal to N

[list_begin arguments]
[arg_def integer N in]
Maximum number to consider
[list_end]

%</man>
% \end{tcl}
% Then |primeFactors| returns the list of the prime numbers in |n|.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::primeFactors] [arg N]]
Return a list of the prime numbers in the number N

[list_begin arguments]
[arg_def integer N in]
Number to be factorised
[list_end]

%</man>
% \end{tcl}
% And |uniquePrimeFactors| does the same, with duplicates removed.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::uniquePrimeFactors] [arg N]]
Return a list of the [emph unique] prime numbers in the number N

[list_begin arguments]
[arg_def integer N in]
Number to be factorised
[list_end]

%</man>
% \end{tcl}
% |factors| returns all factors of |n|, not just just primes.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::factors] [arg N]]
Return a list of all [emph unique] factors in the number N, including 1 and N itself
[list_begin arguments]
[arg_def integer N in]
Number to be factorised
[list_end]

%</man>
% \end{tcl}
% |totient| computes the Euler totient for |n|
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::totient] [arg N]]
Evaluate the Euler totient function for the number N (number of numbers
relatively prime to N)

[list_begin arguments]
[arg_def integer N in]
Number in question
[list_end]

%</man>
% \end{tcl}
% |moebius| computes the Moebious function on |n|.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::moebius] [arg N]]
Evaluate the Moebius function for the number N

[list_begin arguments]
[arg_def integer N in]
Number in question
[list_end]

%</man>
% \end{tcl}
% |legendre| computes the Legendre symbol (|a/p|).
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::legendre] [arg a] [arg p]]
Evaluate the Legendre symbol (a/p)

[list_begin arguments]
[arg_def integer a in]
Upper number in the symbol
[arg_def integer p in]
Lower number in the symbol (must be non-zero)
[list_end]

%</man>
% \end{tcl}
% |jacobi| compute the Jacobi symbol (|a/p|).
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::jacobi] [arg a] [arg b]]
Evaluate the Jacobi symbol (a/b)

[list_begin arguments]
[arg_def integer a in]
Upper number in the symbol
[arg_def integer b in]
Lower number in the symbol (must be odd)
[list_end]

%</man>
% \end{tcl}
% |gcd| computes the greatest common divisor of |n| and |m|.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::gcd] [arg m] [arg n]]
Return the greatest common divisor of [term m] and [term n]

[list_begin arguments]
[arg_def integer m in]
First number
[arg_def integer n in]
Second number
[list_end]

%</man>
% \end{tcl}
% |lcm| computes the least common multiple of |n| and |m|.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::lcm] [arg m] [arg n]]
Return the lowest common multiple of [term m] and [term n]

[list_begin arguments]
[arg_def integer m in]
First number
[arg_def integer n in]
Second number
[list_end]

%</man>
% \end{tcl}
% |numberPrimesGauss| estimates the number of primes below |n|
% using a formula by Gauss.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::numberPrimesGauss] [arg N]]
Estimate the number of primes according the formula by Gauss.

[list_begin arguments]
[arg_def integer N in]
Number in question
[list_end]

%</man>
% \end{tcl}
% |numberPrimesLegendre| estimates the number of primes below |n|
% using a formula by Legendre.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::numberPrimesLegendre] [arg N]]
Estimate the number of primes according the formula by Legendre.

[list_begin arguments]
[arg_def integer N in]
Number in question
[list_end]

%</man>
% \end{tcl}
% |numberPrimesLegendreModified| estimates the number of primes below
% |n| using Legendre's modified formula.
% \begin{tcl}
%<*man>
[call [cmd math::numtheory::numberPrimesLegendreModified] [arg N]]
Estimate the number of primes according the modified formula by Legendre.

[list_begin arguments]
[arg_def integer N in]
Number in question
[list_end]

%</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
}

# Add the additional procedures
#
source [file join [file dirname [info script]] primes.tcl]
%</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 {Add-ons}
%
% A number of additional functions around factoring numbers
%
% \begin{tcl}
%<*pkg_primes>
# ComputeNextPrime --
#     Determine the next prime
#
# Arguments:
#     None
#
# Result:
#     None
#
# Side effects:
#     One prime added to the list of primes
#
# Note:
#     Using a true sieve of Erathostenes might be faster, but
#     this does work. Even computing the first ten thousand
#     does not seem to be slow.
#
proc ::math::numtheory::ComputeNextPrime {} {
    variable primes
    variable nextPrimeCandidate
    variable nextPrimeIncrement

    while {1} {
        #
        # Test the current candidate
        #
        set sqrtCandidate [expr {sqrt($nextPrimeCandidate)}]

        set isprime 1
        foreach p $primes {
            if { $p > $sqrtCandidate } {
                break
            }
            if { $nextPrimeCandidate % $p == 0 } {
                set isprime 0
                break
            }
        }

        if { $isprime } {
            lappend primes $nextPrimeCandidate
        }

        #
        # In any case get the next candidate
        #
        if { $nextPrimeIncrement == 1 } {
            set nextPrimeIncrement 5
            set nextPrimeCandidate [expr {$nextPrimeCandidate + 4}]
        } else {
            set nextPrimeIncrement 1
            set nextPrimeCandidate [expr {$nextPrimeCandidate + 2}]
        }

        if { $isprime } {
            break
        }
    }
}

# firstNprimes --
#     Return the first N primes
#
# Arguments:
#     number           Number of primes to return
#
# Result:
#     List of the first $number primes
#
proc ::math::numtheory::firstNprimes {number} {
    variable primes

    while { [llength $primes] < $number } {
        ComputeNextPrime
    }

    return [lrange $primes 0 [expr {$number-1}]]
}

# primesLowerThan --
#     Return the primes lower than some threshold
#
# Arguments:
#     threshold        Threshold for the primes
#
# Result:
#     List of primes lower/equal to the threshold
#
proc ::math::numtheory::primesLowerThan {threshold} {
    variable primes

    while { [lindex $primes end] < $threshold } {
        ComputeNextPrime
    }

    set n 0
    foreach p $primes {
        if { $p > $threshold } {
            break
        } else {
            incr n
        }
    }
    return [lrange $primes 0 [expr {$n-1}]]
}

# primeFactors --
#     Determine the prime factors of a number
#
# Arguments:
#     number           Number to factorise
#
# Result:
#     List of prime factors
#
proc ::math::numtheory::primeFactors {number} {
    variable primes

    #
    # Make sure we have enough primes
    #
    primesLowerThan [expr {sqrt($number)}]

    set factors {}

    set idx 0

    while { $number > 1 } {
        set p [lindex $primes $idx]
        if {$p == {}} {
            lappend factors $number
            break
        }
        if { $number % $p == 0 } {
            lappend factors $p
            set number [expr {$number/$p}]
        } else {
            incr idx
        }
    }

    return $factors
}

# uniquePrimeFactors --
#     Determine the unique prime factors of a number
#
# Arguments:
#     number           Number to factorise
#
# Result:
#     List of unique prime factors
#
proc ::math::numtheory::uniquePrimeFactors {number} {
    return [lsort -unique -integer [primeFactors $number]]
}

# totient --
#     Evaluate the Euler totient function for a number
#
# Arguments:
#     number           Number in question
#
# Result:
#     Totient of the given number (number of numbers
#     relatively prime to the number)
#
proc ::math::numtheory::totient {number} {
    set factors [uniquePrimeFactors $number]

    set totient 1

    foreach f $factors {
        set totient [expr {$totient * ($f-1)}]
    }

    return $totient
}

# factors --
#     Return all (unique) factors of a number
#
# Arguments:
#     number           Number in question
#
# Result:
#     List of factors including 1 and the number itself
#
# Note:
#     The algorithm for constructing the power set was taken from
#     wiki.tcl.tk/2877 (algorithm subsets2b).
#
proc ::math::numtheory::factors {number} {
    set factors [primeFactors $number]

    #
    # Iterate over the power set of this list
    #
    set result [list 1 $number]
    for {set n 1} {$n < [llength $factors]} {incr n} {
        set subsets [list [list]]
        foreach f $factors {
            foreach subset $subsets {
                lappend subset $f
                if {[llength $subset] == $n} {
                    lappend result [Product $subset]
                } else {
                    lappend subsets $subset
                }
            }
        }
    }
    return [lsort -unique -integer $result]
}

# Product --
#     Auxiliary function: return the product of a list of numbers
#
# Arguments:
#     list           List of numbers
#
# Result:
#     The product of all the numbers
#
proc ::math::numtheory::Product {list} {
    set product 1
    foreach e $list {
        set product [expr {$product * $e}]
    }
    return $product
}

# moebius --
#     Return the value of the Moebius function for "number"
#
# Arguments:
#     number         Number in question
#
# Result:
#     The product of all the numbers
#
proc ::math::numtheory::moebius {number} {
    if { $number < 1 } {
        return -code error "The number must be positive"
    }
    if { $number == 1 } {
        return 1
    }

    set primefactors [primeFactors $number]
    if { [llength $primefactors] != [llength [lsort -unique -integer $primefactors]] } {
        return 0
    } else {
        return [expr {(-1)**([llength $primefactors]%2)}]
    }
}

# legendre --
#     Return the value of the Legendre symbol (a/p)
#
# Arguments:
#     a              Upper number in the symbol
#     p              Lower number in the symbol
#
# Result:
#     The Legendre symbol
#
proc ::math::numtheory::legendre {a p} {
    if { $p == 0 } {
        return -code error "The number p must be non-zero"
    }

    if { $a % $p == 0 } {
        return 0
    }

    #
    # Just take the brute force route
    # (Negative values of a present a small problem, but only a small one)
    #
    while { $a < 0 } {
        set a [expr {$p + $a}]
    }

    set legendre -1
    for {set n 1} {$n < $p} {incr n} {
        if { $n**2 % $p == $a } {
            set legendre 1
            break
        }
    }

    return $legendre
}

# jacobi --
#     Return the value of the Jacobi symbol (a/b)
#
# Arguments:
#     a              Upper number in the symbol
#     b              Lower number in the symbol
#
# Result:
#     The Jacobi symbol
#
# Note:
#     Implementation adopted from the Wiki - http://wiki.tcl.tk/36990
#     encoded by rmelton 9/25/12
#     Further references:
#     http://en.wikipedia.org/wiki/Jacobi_symbol
#     http://2000clicks.com/mathhelp/NumberTh27JacobiSymbolAlgorithm.aspx
#
proc ::math::numtheory::jacobi {a b} {
    if { $b<=0 || ($b&1)==0 } {
        return 0;
    }

    set j 1
    if {$a<0} {
        set a [expr {0-$a}]
        set j [expr {0-$j}]
    }
    while {$a != 0} {
        while {($a&1) == 0} {
            ##/* Process factors of 2: Jacobi(2,b)=-1 if b=3,5 (mod 8) */
            set a [expr {$a>>1}]
            if {(($b & 7)==3) || (($b & 7)==5)} {
                set j [expr {0-$j}]
            }
        }
        ##/* Quadratic reciprocity: Jacobi(a,b)=-Jacobi(b,a) if a=3,b=3 (mod 4) */
        lassign [list $a $b] b a
        if {(($a & 3)==3) && (($b & 3)==3)} {
            set j [expr {0-$j}]
        }
        set a [expr {$a % $b}]
    }
    if {$b==1} {
        return $j
    } else {
        return 0
    }
}

# gcd --
#     Return the greatest common divisor of two numbers n and m
#
# Arguments:
#     n              First number
#     m              Second number
#
# Result:
#     The greatest common divisor
#
proc ::math::numtheory::gcd {n m} {
    #
    # Apply Euclid's good old algorithm
    #
    if { $n > $m } {
        set t $n
        set n $m
        set m $t
    }

    while { $n > 0 } {
        set r [expr {$m % $n}]
        set m $n
        set n $r
    }

    return $m
}

# lcm --
#     Return the lowest common multiple of two numbers n and m
#
# Arguments:
#     n              First number
#     m              Second number
#
# Result:
#     The lowest common multiple
#
proc ::math::numtheory::lcm {n m} {
    set gcd [gcd $n $m]
    return [expr {$n*$m/$gcd}]
}

# numberPrimesGauss --
#     Return the approximate number of primes lower than the given value based on the formula by Gauss
#
# Arguments:
#     limit            The limit for the largest prime to be included in the estimate
#
# Returns:
#     Approximate number of primes
#
proc ::math::numtheory::numberPrimesGauss {limit} {
    if { $limit <= 1 } {
        return -code error "The limit must be larger than 1"
    }
    expr {$limit / log($limit)}
}

# numberPrimesLegendre --
#     Return the approximate number of primes lower than the given value based on the formula by Legendre
#
# Arguments:
#     limit            The limit for the largest prime to be included in the estimate
#
# Returns:
#     Approximate number of primes
#
proc ::math::numtheory::numberPrimesLegendre {limit} {
    if { $limit <= 1 } {
        return -code error "The limit must be larger than 1"
    }
    expr {$limit / (log($limit) - 1.0)}
}

# numberPrimesLegendreModified --
#     Return the approximate number of primes lower than the given value based on the
#     modified formula by Legendre
#
# Arguments:
#     limit            The limit for the largest prime to be included in the estimate
#
# Returns:
#     Approximate number of primes
#
proc ::math::numtheory::numberPrimesLegendreModified {limit} {
    if { $limit <= 1 } {
        return -code error "The limit must be larger than 1"
    }
    expr {$limit / (log($limit) - 1.08366)}
}
%</pkg_primes>
%\end{tcl}
% \begin{tcl}
%<*test_primes>
test first-few-primes-1 "First 10 primes" -match equalLists -body {
    ::math::numtheory::firstNprimes 10
} -result {2 3 5 7 11 13 17 19 23 29}

test first-few-primes-2 "First 12 primes" -match equalLists -body {
    ::math::numtheory::firstNprimes 12
} -result {2 3 5 7 11 13 17 19 23 29 31 37}

test first-few-primes-3 "First 20 primes" -match equalLists -body {
    ::math::numtheory::firstNprimes 20
} -result {2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71}

test primes-lower-than-1 "Primes lower/equal 101" -match equalLists -body {
    ::math::numtheory::primesLowerThan 101
} -result {2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101}

test primes-lower-than-2 "Primes lower/equal 2" -match equalLists -body {
    ::math::numtheory::primesLowerThan 2
} -result {2}

test primes-lower-than-3 "Primes lower/equal 4" -match equalLists -body {
    ::math::numtheory::primesLowerThan 4
} -result {2 3}

test primes-lower-than-4 "Primes lower/equal 102" -match equalLists -body {
    ::math::numtheory::primesLowerThan 102
} -result {2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101}

test prime-factors-1 "Prime factors 100" -match equalLists -body {
    ::math::numtheory::primeFactors 100
} -result {2 2 5 5}

test prime-factors-2 "Unique prime factors 100" -match equalLists -body {
    ::math::numtheory::uniquePrimeFactors 100
} -result {2 5}

test prime-factors-3 "Prime factors 2900" -match equalLists -body {
    ::math::numtheory::primeFactors 2900
} -result {2 2 5 5 29}

test prime-factors-4 "Unique prime factors 2900" -match equalLists -body {
    ::math::numtheory::uniquePrimeFactors 2900
} -result {2 5 29}

test prime-factors-5 "Prime factors 964" -match equalLists -body {
    ::math::numtheory::primeFactors 964
} -result {2 2 241}

test prime-factors-6 "Prime factors 960" -match equalLists -body {
    ::math::numtheory::primeFactors 960
} -result {2 2 2 2 2 2 3 5}

test prime-factors-7 "Prime factors 914" -match equalLists -body {
    ::math::numtheory::primeFactors 914
} -result {2 457}

test totient-1 "Totient 15" -body {
    ::math::numtheory::totient 15
} -result 8

test totient-2 "Totient 30" -body {
    ::math::numtheory::totient 30
} -result 8

test totient-3 "Totient 35" -body {
    ::math::numtheory::totient 35
} -result 24

test totient-4 "Totient 105" -body {
    ::math::numtheory::totient 105
} -result 48

test factors-1 "All factors 30" -match equalLists -body {
    ::math::numtheory::factors 30
} -result {1 2 3 5 6 10 15 30}

test factors-1 "All factors 128" -match equalLists -body {
    ::math::numtheory::factors 128
} -result {1 2 4 8 16 32 64 128}

test factors-1 "All factors 250" -match equalLists -body {
    ::math::numtheory::factors 250
} -result {1 2 5 10 25 50 125 250}

test moebius-1 "Moebius for first 19 numbers" -match equalLists -body {
    set result {}
    for {set n 1} {$n < 20} {incr n} {
        lappend result [::math::numtheory::moebius $n]
    }
    set result
} -result {1 -1 -1 0 -1 1 -1 0 0 1 -1 0 -1 1 1 0 -1 0 -1}

test legendre-1 "Legendre symbol (-1/3)" -body {
    ::math::numtheory::legendre -1 3
} -result -1
test legendre-2 "Legendre symbol (-3/7)" -body {
    ::math::numtheory::legendre -3 7
} -result 1

test jacobi-1 "Jacobi symbol (6/7)" -body {
    ::math::numtheory::jacobi 6 7
} -result -1

test jacobi-2 "Jacobi symbol (6/9)" -body {
    ::math::numtheory::jacobi 6 9
} -result 0

test jacobi-3 "Jacobi symbol (3/11)" -body {
    ::math::numtheory::jacobi 3 11
} -result 1

test gcd-1 "Greatest common divisor 2 and 3" -body {
    ::math::numtheory::gcd 2 3
} -result 1

test gcd-2 "Greatest common divisor 20 and 12" -body {
    ::math::numtheory::gcd 20 12
} -result 4

test gcd-3 "Greatest common divisor 600 and 125" -body {
    ::math::numtheory::gcd 600 125
} -result 25

test lcm-1 "Lowest common multiple 3 and 4" -body {
    ::math::numtheory::lcm 3 4
} -result 12

test lcm-2 "Lowest common multiple 12 and 20" -body {
    ::math::numtheory::lcm 12 20
} -result 60

test number-primes "Exercise prime estimators" -match equalLists -body {
    set estimate1 [::math::numtheory::numberPrimesGauss 1000]
    set estimate2 [::math::numtheory::numberPrimesLegendre 1000]
    set estimate3 [::math::numtheory::numberPrimesLegendreModified 1000]
    set result [list [expr {int($estimate1)}] [expr {int($estimate2)}] [expr {int($estimate3)}]]
} -result {144 169 171}
%</test_primes>
%\end{tcl}
% 
% \section*{Closings}
% 
% \begin{tcl}
%<*man>
[list_end]

[vset CATEGORY {math :: numtheory}]
[include ../common-text/feedback.inc]
[manpage_end]
%</man>
% \end{tcl}
% 
% \begin{tcl}
%<test_common>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