An Elementary Method For Fast Modular Exponentiation With Factored Modulus

1. Introduction

An important problem at the intersection of cryptography and number theory is the Modular Exponentiation Problem. This is the problem of computing $a^{n}\pmod{m}$ for positive integers $a,n,m$ . It happens that solving for $n$ in the equation $a^{n}\equiv b\pmod{m}$ , called the Discrete Logarithm Problem (DLP), is computationally hard. This makes modular exponentiation useful in many cryptosystems, most notably the RSA public key cryptosystem and the Diffie-Hellman key exchange. The standard algorithm to solve the Modular Exponentiation Problem is the Repeated Squaring Algorithm, which runs in $O(\log n)$ steps (in this paper, a step is a modular multiplication). To our knowledge, there are no existing significant asymptotic improvements to this algorithm.

In this paper, we present an efficient algorithm to solve the Modular Exponentiation Problem when the factorization of the modulus is known. Such a situation is potentially practically useful, as modular exponentiation is generally performed for encryption and thus the user has the freedom of choosing $m$ and knowings its factorization. In the 1984 paper [1], Bach proves the existence of a Hensel lifting type algorithm that can lift solutions to the discrete logarithm from modulo $p$ to modulo $p^{e}$ in polynomial time. One can intuit that in a similar way, there exists a speedup to modular exponentiation when the factorization of the modulus has some prime power in it; this is what our algorithm provides. The algorithm is quite elementary: it hinges on the binomial theorem and the clever recursive computation of inverses and binomial coefficients. It also depends on a set of parameters: If $m$ can be factored as $\prod_{i=1}^{k}p_{i}^{e_{i}}$ , we choose for each $i$ an integer parameter $t_{i}\in[1,e_{i}]$ . We then compute $a^{n}\pmod{m}$ in $O(\max(e_{i}/t_{i})+\sum_{i=1}^{k}t_{i}\log p_{i})$ steps (again, steps are modular multiplications). We additionally provide an $O(1)$ memory algorithm for memory-sensitive scenarios, where we define $O(1)$ memory to mean the memory of an integer modulo $m$ , or $\log m$ bits. For general $m$ we make asymptotic improvements as measured by specific metrics. For a certain infinite family of $m$ we achieve a complexity of $O(\sqrt{\log m})$ . Furthermore, our algorithm profits from the development of other fast algorithms. For example, if one were to discover a faster modular exponentiation algorithm, we could use this as the modExp function in our algorithm. This, in turn, makes our algorithm faster than the default one.

In section 2 of this paper, we lay out some preliminary results that our algorithm and analysis require. In section 3, we present the algorithm with pseudocode. In section 4, we mathematically analyze our algorithm. Because its complexity depends on number-theoretic properties of $m$ , we require results from analytic number theory to estimate the “average” complexity. We discuss families of $m$ for which we make significant improvements. Then, in section 5, we present an version of our algorithm for matrices, linear recurrences, and $\mathbb{Z}/p^{k}\mathbb{Z}$ ring extensions. Finally, in section 6, we test our algorithm for prime powers against Python’s built-in $\mathrm{pow}$ function to show our algorithm does practically. A Python implementation of the algorithm and programmatical analysis for it can be found at [2].

2. Preliminaries

The complexity of our algorithm depends on the number-theoretic properties of the modulus, so we will need multiple preliminary definitions to continue with this analysis. In general, we use the variable $p$ for a prime. We denote the canonical decomposition $m$ as $m=\prod p_{i}^{e_{i}}$ . By convention, we let $a\pmod{b}$ denote the least residue of $a$ modulo $b$ . We let $\nu_{p}(n)$ denote the $p$ -adic valuation of $n$ , and we use this notation interchangeably with $\nu(n,p)$ . We let $\zeta(\bullet)$ be the Riemann Zeta Function and let $\varphi(\bullet)$ denote Euler’s Totient Function. Additionally, we let $\vartheta(\bullet)$ denote Chebyshev’s Function. We must define the radical of an integer $n$ because this is intimately related to the optimal complexity that our algorithm may achieve:

Definition 2.1.

If $n=\prod_{i=1}^{k}p_{i}^{e_{i}}$ , define the radical of $n$ as $\operatorname{rad}n=\prod_{i=1}^{k}p_{i}$ .

A useful notation that will assist our analysis is the following:

Definition 2.2.

Let $n=\prod_{i=1}^{k}p_{i}^{e_{i}}$ , and define the multiset $S=\{x_{1},x_{2},\cdots,x_{k}\}$ . We define $H_{S}(n):=\max\left(\frac{e_{i}}{x_{i}}\right)$ , and set $H(n):=H_{\{1,1,\ldots,1\}}(n)$ .

Our algorithm hinges on the following theorem, due to Euler:

Theorem 2.3.

(Euler) Let $a,n$ be relatively prime positive integers. Then $a^{\varphi(n)}\equiv 1\pmod{n}$ .

Another result due to Euler which we need to use in our analysis is the Euler Product formula:

Theorem 2.4.

(Euler) Let $a(n):\mathbb{N}\to\mathbb{C}$ be a multiplicative function, i.e., $a(mn)=a(m)a(n)$ is true when $\gcd(m,n)=1$ . Then

\sum_{n\geq 1}\frac{a(n)}{n^{s}}=\prod_{p\leavevmode\nobreak\ \mathrm{prime}}\sum_{k=0}^{\infty}\frac{a(p^{k})}{p^{ks}}.

In particular, we have the following Euler Product for $\zeta$ :

\zeta(s)=\prod_{p\leavevmode\nobreak\ \text{prime}}\frac{1}{1-p^{-s}}.

To acquire precise asymptotics for specific sums, we need the following asymptotics:

Theorem 2.5.

(Stirling’s Approximation) $\log(x!)=x\log x-x+O(\log x)$

Theorem 2.6.

(Chebyshev) $\vartheta(x)=O(x)$

To approximate average order summations, we will need Abel’s summation formula, referred to as partial summation:

Theorem 2.7.

(Abel) Let $(a_{n})^{\infty}_{n=0}$ be a sequence of complex numbers. Define $A(t)=\sum_{0\leq n\leq t}a_{n}$ . Fix real numbers $x<y$ and let $\phi$ be a continuously differentiable function on $[x,y]$ . Then

\sum_{x<n\leq y}a_{n}\phi(n)=A(y)\phi(y)-A(x)\phi(x)-\int_{x}^{y}A(u)\phi^{\prime}(u)\mathrm{d}u.

We compare our algorithm to the standard repeated squaring method for modular exponentiation (there is no general asymptotic improvement to this method, as far as we know). We may compute $a^{n}\pmod{m}$ in $O(\log n)$ steps by this method. By Euler’s theorem, we can reduce this to $O(\log(\varphi(m)))$ which is considered $O(\log m)$ .

For a ring $R$ , we let $GL_{n}(R)$ denote the general linear group over $R$ . We will need the well-known fact that $|GL_{n}(\mathbb{F}_{p})|=\prod_{i=0}^{n-1}(p^{n}-p^{i})$ . We will also need the following theorem due to Lagrange:

Theorem 2.8.

(Lagrange) Let $G$ be a group of order $n$ . Then for every $a\in G$ , $a^{n}$ is the identity.

To compute the order of general linear groups over $\mathbb{Z}/T\mathbb{Z}$ , we must first compute the order of $GL_{n}(\mathbb{Z}/p^{k}\mathbb{Z})$ . We do this with the following lemma:

Lemma 2.9.

For a prime $p$ and a $k\geq 0$ we have

|GL_{n}(\mathbb{Z}/p^{k}\mathbb{Z})|=p^{(k-1)n^{2}}|GL_{n}(\mathbb{F}_{p})|=p^{(k-1)n^{2}}\prod_{i=0}^{n-1}(p^{n}-p^{i}).

Proof.

We proceed by induction on $k$ , with $k=1$ trivial. Notice the natural ring homomorphism $\varphi:\mathbb{Z}/p^{k}\mathbb{Z}\to\mathbb{Z}/p^{k-1}\mathbb{Z}$ given by $a+p^{k}\mathbb{Z}\mapsto a+p^{k-1}\mathbb{Z}$ . This induces the surjection $\varphi_{n}:GL_{n}(\mathbb{Z}/p^{k}\mathbb{Z})\to GL_{n}(\mathbb{Z}/p^{k-1}\mathbb{Z})$ . Hence $|GL_{n}(\mathbb{Z}/p^{k}\mathbb{Z})|=|\ker(\varphi_{n})||GL_{n}(\mathbb{Z}/p^{k-1}\mathbb{Z})|$ . Since there are $p$ choices for each entry of a matrix in the kernel of $\varphi_{n}$ , the size of the kernel is $p^{n^{2}}$ , and we have the desired result. ∎

This allows us to prove the following important lemma:

Lemma 2.10.

If $T=p_{1}^{t_{1}}p_{2}^{t_{2}}\cdots p_{k}^{t_{k}}$ is the canonical factorization of some $T\in\mathbb{N}$ ,

|GL_{n}(\mathbb{Z}/T\mathbb{Z})|=\prod_{j=1}^{k}p_{j}^{(t_{j}-1)n^{2}}\prod_{i=0}^{n-1}(p_{j}^{n}-p_{j}^{i}).

Proof.

By the Chinese Remainder Theorem, we have the ring isomorphism

\mathbb{Z}/T\mathbb{Z}\to\prod_{j=1}^{k}\mathbb{Z}/p_{j}^{t_{j}}\mathbb{Z},

and thus there is a corresponding isomorphism

GL_{n}(\mathbb{Z}/T\mathbb{Z})\to\prod_{j=1}^{k}GL_{n}(\mathbb{Z}/p_{j}^{t_{j}}\mathbb{Z}).

This implies the desired result by Lemma 2.9. ∎

Remark 2.11.

It follows from Lemma 2.10 that we may compute $|GL_{n}(\mathbb{Z}/T\mathbb{Z})|$ in $O(kn)$ steps, which is insignificant.

We use standard convention (from [3]) for the asymptotic notations $O(\bullet),o(\bullet),\ll,\gg,\sim$ , and $\Omega(\bullet)$ . We use $<O(\bullet)$ and $>O(\bullet)$ rather unconventionally: $f(x)<O(g(x))$ provided there is a function $h(x)=O(g(x))$ such that $f(x)<h(x)$ for sufficiently large $x$ . We use a similar definition for $>O(\bullet)$ .

3. The Algorithm

Our main theorem is the following:

Theorem 3.1.

Let $m$ be a positive integer with known factorization $p_{1}^{e_{1}}p_{2}^{e_{2}}\cdots p_{k}^{e_{k}}$ and let $a,n$ be positive integers such that $\gcd(a,m)=1$ . For each $1\leq i\leq k$ , choose an integer parameter $t_{i}\in[1,e_{i}]$ . With $\mathcal{T}=\{t_{1},t_{2},\cdots,t_{k}\}$ , we may compute $a^{n}\pmod{m}$ in

O\left(H_{\mathcal{T}}(m)+\sum_{i=1}^{k}t_{i}\log p_{i}\right)

steps.

First, we need a definition that will simplify our calculation of modular inverses.

Definition 3.2.

For some $a$ and $m=\prod p_{i}^{e_{i}}$ , define $\{u_{a},v_{a}\}$ to be the inverse pair of $a$ modulo $m$ , where $v_{a}=\prod p_{i}^{\nu(a,p_{i})}$ and $u_{a}\equiv(a/v_{a})^{-1}\pmod{m}$ . Define the inverse pair of $0$ as $\{0,0\}$ .

Example 1.

Let us compute the inverse pair of $12$ modulo $20$ as an example. We have $v_{12}=2^{\nu(12,2)}\cdot 5^{\nu(12,5)}=4$ . Then, $u_{12}\equiv(12/4)^{-1}\equiv 7\pmod{20}$ . So the inverse pair of $12$ modulo $20$ is $\{4,7\}$ . Observe also that $4\cdot 7^{-1}\equiv 12\pmod{20}$ .

Notice that $u_{a}$ always exists because $a/v_{a}$ is clearly invertible modulo $m$ . Intuitively, inverse pairs are inverses equipped with an extra parameter that allows for computation of fractions with denominator not necessarily relatively prime to $m$ .

Consider the following lemma that allows for recursive computation of modular inverses.

Lemma 3.3.

If all of $1,2,\dots,a$ are invertible modulo $m$ ,

a^{-1}\equiv-(m\pmod{a})^{-1}\Big{\lfloor}\frac{m}{a}\Big{\rfloor}\pmod{m}.

Proof.

Note $m=a\lfloor\frac{m}{a}\rfloor+(m\pmod{a})$ . Working in $\mathbb{Z}/m\mathbb{Z}$ ,

(m\pmod{a})^{-1}\cdot\left(m-\Big{\lfloor}\frac{m}{a}\Big{\rfloor}\right)=\frac{m-\lfloor\frac{m}{a}\rfloor}{m-a\lfloor\frac{m}{a}\rfloor}=\frac{-\lfloor\frac{m}{a}\rfloor}{-a\lfloor\frac{m}{a}\rfloor}=\frac{1}{a},

as desired. ∎

We now extend this result from a recursive computation of modular inverses to a recursive computation of inverse pairs.

Lemma 3.4.

We may linearly compute $\{u_{n},v_{n}\}$ given the inverse pairs $\bigcup_{1\leq i<n}\{u_{i},v_{i}\}$ .

Proof.

We claim that the following formula holds:

u_{n}=\begin{cases}u_{m\%n}\cdot\frac{-\lfloor\frac{m}{n}\rfloor}{v_{m\%n}}\pmod{m}&v_{n}=1\\ u_{n/v_{n}}\pmod{m}&v_{n}>1.\end{cases}

Note that $\gcd(n/v_{n},m)=1$ by definition, so that

u_{n/v_{n}}\equiv\left(\frac{n}{v_{n}}\right)^{-1}\equiv u_{n}\pmod{m}.

Suppose $v_{n}=1$ , so that $m$ and $n$ are relatively prime. Decompose $m=qn+r$ via the division algorithm. We wish to show that

u_{n}\equiv u_{r}\cdot\frac{q}{v_{r}}\equiv\frac{v_{r}}{r}\cdot\frac{q}{v_{r}}\equiv\frac{q}{r}\pmod{m}.

Note that

u_{n}\equiv n^{-1},

so we wish to show that

qn+r\equiv 0\pmod{m},

which is obvious. There are no inversion issues as $\gcd(r,m)=1$ by the Euclidean Algorithm, and $\gcd(v_{r},m)=1$ as $v_{r}\mid r$ . ∎

We are now equipped to prove Theorem 3.1. The algorithm is as follows.

Proof.

First, define $T=p_{1}^{t_{1}}p_{2}^{t_{2}}\cdots p_{i}^{t_{i}}$ and $\Phi=\varphi(T)$ . Decompose $n=M\Phi+r$ with the division algorithm. Then,

a^{n}\equiv a^{M\Phi+r}\equiv a^{r}\cdot(a^{\Phi})^{M}\pmod{m}.

Now, $a^{r}$ can be computed with the standard Repeated Squaring Algorithm. The complexity of this is $O(\log r)$ . By the division algorithm, $r<\Phi<T$ , so this is $O(\log T)$ . For computation of the second term, first notice that by Euler’s theorem, $a^{\Phi}\equiv 1\pmod{T}$ so there exists some integer $s$ with $a^{\Phi}=Ts+1$ . Now, we expand using the binomial theorem:

(a^{\Phi})^{M}=(Ts+1)^{M},

=1\cdot\binom{M}{0}+Ts\cdot\binom{M}{1}+\cdots+(Ts)^{i}\binom{M}{i}+\cdots+(Ts)^{M}\binom{M}{M}.

Let $\ell=1+\text{max}(\lfloor e_{i}/t_{i}\rfloor)$ . Only the first $\ell$ terms need to be computed, as all the later terms are $0$ modulo $m$ . This is because for a prime $p_{i}$ ,

\nu_{p_{i}}(T^{\ell})=\nu_{p_{i}}(\left(p_{i}^{t_{i}}\right)^{\ell})=\nu_{p_{i}}(p_{i}^{t_{i}(\lfloor e_{i}/t_{i}\rfloor+1)})=t_{i}(\lfloor e_{i}/t_{i}\rfloor+1)\geq t_{i}(e_{i}/t_{i})=e_{i},

so $p_{i}^{e_{i}}\mid T^{\ell}$ for all $i$ , and thus $m\mid T^{\ell}$ and $T^{\ell}\equiv 0\pmod{m}$ .

Additionally, we can linearly compute the $\binom{M}{i}$ terms with the identity

\binom{M}{i+1}=\frac{M!}{(i+1)!(M-i-1)!},

=\frac{M-i}{i+1}\frac{M!}{i!(M-i)!}=\frac{M-i}{i+1}\binom{M}{i}.

Notice that in our binomial recursion, the one inverse we need to calculate $\binom{M}{i}$ is the inverse of $i$ . Using Lemma 3.4, because we are performing a linear computation, we can compute the inverse pairs of $0$ through $\ell-1$ modulo $m$ in $O(\ell)$ time and $O(\ell)$ space. This is the only time we use more than $O(1)$ memory in this algorithm. From here, we can use our identity to compute $\binom{M}{0}$ to $\binom{M}{\ell-1}$ in $O(\ell)$ time and space. This is the final aspect of our calculation, and we may now recover $a^{n}\pmod{m}$ . It is now clear that the number of steps is $O(\ell+\log T)$ and the space complexity is $O(\ell)$ . Because $\ell=1+H_{\mathcal{T}}(m)$ and

\log T=\log\left(\prod_{i=1}^{k}p_{i}^{t_{i}}\right)=\sum_{i=1}^{k}\log(p_{i}^{t_{i}})=\sum_{i=1}^{k}t_{i}\log p_{i},

this concludes the proof of our main theorem. ∎

Example 2.

We present the details of the algorithm in the special case of $k=1$ , $p_{1}=p$ , $e_{1}=e$ , and $t_{1}=1$ (so $T=p$ ). Our goal is to compute $a^{n}\pmod{p^{e}}$ . We let $n=(p-1)m+r$ , so that

a^{n}\equiv(a^{p-1})^{m}a^{r}.

Notice that $a^{p-1}\equiv 1\pmod{p}$ , so we let $a^{p-1}=1+sp$ . Compute $a^{r}$ via the standard modExp algorithm. Now compute

(a^{p-1})^{m}=(1+sp)^{m}\equiv 1+\binom{m}{1}(sp)^{1}+\binom{m}{2}(sp)^{2}+\cdots+\binom{m}{e-1}(sp)^{e-1}\pmod{p^{e}},

by the Binomial Theorem. We can recursively compute the $\binom{m}{i}$ coefficients modulo $p^{e}$ by the method outlined above. Thus we may compute each $\binom{m}{i}(sp)^{i}$ term for $0\leq i\leq e-1$ . Finally, compute $a^{n}\equiv(a^{p-1})^{m}a^{r}\pmod{p^{e}}$ .

Example 3.

As a particular case, let us compute $7^{123}\pmod{11^{3}}$ . Decompose $123=12\cdot 10+3$ . Compute $7^{3}\equiv 343\pmod{11^{3}}$ . Write $s=\frac{7^{10}-1}{11}$ . Notice $s$ is an integer by Euler. Then

7^{120}\equiv(7^{10})^{12}\equiv(1+11s)^{12}=\binom{12}{0}+\binom{12}{1}\cdot 11s+\binom{12}{2}\cdot(11s)^{2}\equiv 23\pmod{11^{3}}.

Therefore, $7^{123}\equiv 343\cdot 23\equiv 1234\pmod{11^{3}}$ .

Notice that this algorithm is $O(H_{\mathcal{T}}(m))$ space, which is not always ideal. For more memory-sensitive scenarios, we may compute inverses directly rather than recursively. This gives rise to the following theorem:

Theorem 3.5.

Let $m$ be a positive integer with known factorization $p_{1}^{e_{1}}p_{2}^{e_{2}}\cdots p_{k}^{e_{k}}$ and let $a,n$ be positive integers such that $\gcd(a,m)=1$ . For each $1\leq i\leq k$ , choose an integer parameter $t_{i}\in[1,e_{i}]$ . With $\mathcal{T}=\{t_{1},t_{2},\cdots,t_{k}\}$ , we may compute $a^{n}\pmod{m}$ in

O\left(H_{\mathcal{T}}(m)\log(H_{\mathcal{T}}(m))+\sum_{i=1}^{k}t_{i}\log p_{i}\right)

steps and $O(1)$ memory.

Proof.

Identical to Theorem 3.1, except instead of recursively computing inverse pairs, compute the modular inverses $u_{n}$ with standard inversion algorithms: $u_{n}\equiv(n/v_{n})^{-1}\pmod{m}$ . If, for example, we use the Extended Euclidean Algorithm then calculating the inverses of $1$ to $\ell$ has complexity $O(\sum_{n=1}^{\ell}\log(n/v_{n}))=O(\sum_{n=1}^{\ell}\log n)=O(\ell\log\ell)$ We do not give pseudocode for this algorithm, instead see [2]. ∎

3.1. Pseudocode

Define $\textsc{modExp}(a,b,c)$ to be the standard Repeated Squaring Algorithm for computing $a^{b}\pmod{c}$ . We use the integer division notation $a//b=\lfloor\frac{a}{b}\rfloor$ . The pseudocode builds up to our algorithm, $\textsc{FastModExp}(a$ , $n$ , $P$ , $E$ , $T)$ . Here $a$ and $n$ are positive integers, as usual. $P$ and $E$ are equal sized arrays, $P$ is an array of primes and $E$ is an array of positive integers. $T$ also has the same length as $P$ and $E$ , and is an array of non-negative integer parameters such that $T[i]\leq E[i]$ for all $i$ . We store the inverse pairs in a 2-dimensional array $L$ . This algorithm computes $a^{n}$ modulo $m=\prod_{i=1}^{\mathrm{len}(P)}P[i]^{E[i]}$ .

Algorithm 1 Computes

\nu(n,p)

function

\nu

(

n

,

p

)

if

n\%p=0

then

return

1+\nu(n//p,p)

end if

return

0

end function

Algorithm 2 Computes the inverse pair of

i

given inverse pairs of

0

through

i-1

function nextInversePair(

i

,

m

,

L

,

P

)

v\leftarrow 1

j\leftarrow 0

while

j<\text{len}(P)

do

v\leftarrow v\cdot P[j]^{\nu(i,P[j])}

j\leftarrow j+1

end while

if

v=1

then

u\leftarrow(L[m\%i][0]\cdot((m-m//i)//L[m\%i][1]))\%m

end if

if

v\neq 1

then

u\leftarrow L[i//v][0]

end if

return

[u,v]

end function

Algorithm 3 Computes the inverse pairs of

0

through

i

function generateInversePairs(

i

,

m

,

P

)

L\leftarrow[[0,0],[1,1]]

j\leftarrow 2

while

j<i+1

do

L\leftarrow L+\textsc{nextInversePair}(j\text{, }m\text{, }L\text{, }P)

j\leftarrow j+1

end while

return

L

end function

Algorithm 4 Computes

a^{n}

modulo

m=\prod P[i]^{E[i]}

function FastModExp(

a

,

n

,

P

,

E

,

T

)

t\leftarrow 1

\phi\leftarrow 1

m\leftarrow 1

i\leftarrow 0

while

i<\text{len}(P)

do

\text{temp}\leftarrow P[i]^{T[i]-1}

\phi=\phi\cdot\text{temp}\cdot(P[i]-1)

t=t\cdot\text{temp}\cdot P[i]

m=m\cdot P[i]^{E[i]}

i\leftarrow i+1

end while

r\leftarrow n\%\phi

q\leftarrow(n-r)//\phi

c\leftarrow\textsc{modExp}(a,\phi,m)-1

\text{sum}\leftarrow 0

\text{choose}=1

\text{cExp}\leftarrow 1

\ell\leftarrow 0

i\leftarrow 0

while

i<\text{len}(P)

do

if

\ell<E[i]//T[i]

then

\ell\leftarrow E[i]//T[i]

end if

i\leftarrow i+1

end while

\text{inverses}\leftarrow\textsc{generateInversePairs}(\ell\text{, }m\text{, }P)

i\leftarrow 0

while

i<\textsc{min}(\ell\text{, }q+1)

do

\text{sum}\leftarrow(\text{sum}+(\text{choose}\cdot\text{cExp}))\%m

\text{cExp}\leftarrow(\text{cExp}\cdot c)\%m

\text{choose}\leftarrow(((\text{choose}\cdot(q-i))\%m)//\text{inverses}[i+1][1]\cdot\text{inverses}[i+1][0])\%m

i\leftarrow i+1

end while

\text{ar}\leftarrow\textsc{modExp}(a,r,m)

return

(\text{sum}\cdot\text{ar})\%m

end function

4. Mathematical Analysis

It’s important to optimally choose the parameters $t_{i}$ . For general $m$ , the optimal choice is of the form $t_{i}=O(1)$ for all $i$ . This is due to the following theorem:

Theorem 4.1.

[4] We have that

\lim_{n\to\infty}\frac{1}{n}\sum_{k\leq n}H(k)=C,

where $C=1+\sum_{k\geq 2}\left(1-\frac{1}{\zeta(k)}\right)\approx 1.705$ is Niven’s constant.

In other words, for general $m$ we expect $H(m)$ to be $O(1)$ . The choice of $t_{i}=O(1)\forall i$ is then clear. Let us proceed with some analysis of this choice.

By some metrics, we make asymptotic improvement to the standard $O(\log m)$ repeated squaring methods with $t_{i}=1$ (i.e., $T=\operatorname{rad}m$ ). But by other metrics, we may not. First, let us sweep the $H(m)$ term under the rug (we can do so due to Theorem 4.1). If we just compute $\sum_{m\leq x}\frac{\log m}{\log\operatorname{rad}m}$ or $\sum_{m\leq x}(\log m-\log\operatorname{rad}m)$ , we get $\Omega(x)$ . However, the sum $\sum_{m\leq x}\frac{m}{\operatorname{rad}m}$ is not $O(x)$ , and is not even $O(x(\log x)^{A})$ for any $A$ .

Theorem 4.2.

For any $x\in\mathbb{N}$ , we have that

\sum_{m\leq x}\frac{\log m}{\log\operatorname{rad}m}=\Omega(x)=\sum_{m\leq x}(\log m-\log\operatorname{rad}m)

Proof.

By the simple estimate that $\frac{a}{b}\leq a-b+1$ for $a\geq b\geq 1$ , we have that

\sum_{m\leq x}\frac{\log m}{\log\operatorname{rad}m}\leq O(x)+\sum_{m\leq x}(\log m-\log\operatorname{rad}m).

Note that, by Stirling’s approximation,

\sum_{m\leq x}\log m=\log(x!)=x\log x+O(x),

and

\sum_{m\leq x}\log\operatorname{rad}m=\sum_{m\leq x}\log\left(\prod_{p\mid m}p\right)=\sum_{m\leq x}\sum_{p\mid m}\log p.

Swapping the order of the summation, we may rewrite this sum as

\sum_{p\leq x}\log p\left\lfloor\frac{x}{p}\right\rfloor=\sum_{p\leq x}\frac{x}{p}\log p-O\left(\sum_{p\leq x}\log p\right)=x\sum_{p\leq x}\frac{\log p}{p}-O(\vartheta(x)).

Thus, applying Merten’s first theorem and Chebyshev’s asymptotics for $\vartheta$ , we have that

\sum_{m\leq x}\log\operatorname{rad}m=x\log x-O(x).

Therefore,

\sum_{m\leq x}\frac{\log m}{\log\operatorname{rad}m}=O(x).

On the other hand, $\log m\geq\log\operatorname{rad}m$ , so that

\sum_{m\leq x}\frac{\log m}{\log\operatorname{rad}m}\geq x.

Therefore,

\sum_{m\leq x}\frac{\log m}{\log\operatorname{rad}m}=\Omega(x),

as desired. ∎

Theorem 4.3.

For any $x\in\mathbb{N}$ and $A>0$ , we have that

\sum_{m\leq x}\frac{m}{\operatorname{rad}m}\neq O(x(\log x)^{A})

Proof.

Notice that $f(m):=\frac{m}{\mathrm{rad}(m)}$ is multiplicative, and $f(p^{k})=p^{k-1}$ for prime $p$ and $k\in\mathbb{N}$ . Therefore, the Dirichlet Series of $f$ is

F(s)=\sum_{n=1}^{\infty}\frac{f(n)}{n^{s}}=\prod_{p}\left(1+\sum_{k\geq 1}\frac{p^{k-1}}{p^{ks}}\right).

This expression simplifies to

F(s)=\prod_{p}\left(1+\frac{1}{p^{s}-p}\right).

Thus

\frac{F(s)}{(\zeta(s))^{A+1}}=\prod_{p}\left(\left(1+\frac{1}{p^{s}-p}\right)(1-p^{-s})^{A+1}\right)

by the Euler product for $\zeta$ . Sending $s\to 1^{+}$ on the real axis, the product tends to infinity. Because $\lim_{s\to 1^{+}}\zeta(s)^{A+1}(1-s)^{A+1}=1$ , $F(s)\neq O\left(\frac{1}{(s-1)^{A+1}}\right)$ . Let $S(x)=\sum_{n\leq x}f(n)$ and suppose for the sake of contradiction that $S(x)=O(x(\log x)^{A})$ . For $s>1$ we have that, by partial summation,

F(s)=\int_{1^{-}}^{\infty}\frac{\mathrm{d}S(x)}{x^{s}}=\left[\frac{S(x)}{x^{s}}\right]_{1^{-}}^{\infty}+\int_{1}^{\infty}\frac{S(x)}{x^{s+1}}\mathrm{d}x.

Note that because $S(x)=O(x(\log x)^{A})$ the first term on the RHS vanishes. Thus $F(s)$ is on the order of

\int_{1}^{\infty}\frac{(\log x)^{A}}{x^{s}}\mathrm{d}x.

To arrive at the desired contradiction, we must show that this integral is $O\left(\frac{1}{(s-1)^{A+1}}\right)$ . In order to do this, we induct on $A$ (we may assume $A$ is an integer). For $A=0$ , the result is obvious. For higher $A$ , we may integrate by parts. Set $u=(\log x)^{A}$ and $\mathrm{d}v=x^{-s}\mathrm{d}x$ . Then

\int\frac{(\log x)^{A}}{x^{s}}\mathrm{d}x=\frac{(\log x)^{A}x^{1-s}}{1-s}-\frac{A}{1-s}\int\frac{(\log x)^{A-1}}{x^{s}}\mathrm{d}x,

which implies the result by the induction hypothesis, as both terms are $O(x(\log x)^{A})$ . ∎

In fact, there is a more precise estimate than Theorem 4.3.

Theorem 4.4.

[5] For any $x\in\mathbb{N}$ we have that

\sum_{m\leq x}\frac{m}{\operatorname{rad}m}=x\exp\left((1+o(1))\sqrt{\frac{8\log x}{\log\log x}}\right)

By the metric given by Theorem 4.2, we make an $O(1)$ improvement “on average” to the repeated squaring method. However, by Theorem 4.3, we “expect” $m\gg(\log m)^{A}(\operatorname{rad}m)$ (very loosely), so that we “expect” $\log m\gg\log\operatorname{rad}m+A\log\log m$ . By this metric, we do make asymptotic improvement to the repeated squaring method.

Nonetheless, it is particularly fruitful to work with smooth $m$ rather than general $m$ . We have the following corollary of Theorem 3.1:

Corollary 4.5.

Let $m$ be a positive integer with known factorization $p_{1}^{e_{1}}p_{2}^{e_{2}}\cdots p_{k}^{e_{k}}$ such that all primes $p_{i}$ have the same bit length as an integer $P$ , and all $e_{i}$ are such that $e_{i}\sim k\log P$ . Let $a,n$ be positive integers such that $\gcd(a,m)=1$ . We may then compute $a^{n}\pmod{m}$ in $O(\sqrt{\log m})$ steps.

Proof.

The idea is to use Theorem 3.1 and set $H_{\mathcal{T}}(n)\approx\mathsf{c}$ for some $\mathsf{c}$ , and choose $t_{i}$ such that $\frac{e_{i}}{t_{i}}\approx\mathsf{c}$ for all $i$ . In other words, take $t_{i}=\lfloor\frac{e_{i}}{\mathsf{c}}\rfloor$ , for example. Then $H_{\mathcal{T}}(n)=\mathsf{c}+O(1)$ , and

\sum_{i=1}^{k}t_{i}\log p_{i}=\frac{1}{\mathsf{c}}\sum_{i=1}^{k}(e_{i}\log p_{i})+O\left(\sum_{i=1}^{k}\log p_{i}\right)=\frac{\log m}{\mathsf{c}}+O(\log\operatorname{rad}m).

Notice that

\sqrt{\log m}=\sqrt{\log P\sum_{i=1}^{k}e_{i}}\sim k\log P,

so we may choose $\mathsf{c}$ on the order of $\sqrt{\log m}$ (as $\mathsf{c}<e_{i}$ is necessary and sufficient). With this choice, we compute $a^{n}\pmod{m}$ in

O(\sqrt{\log m})+O(\log\operatorname{rad}m)

steps. Note that

\log\operatorname{rad}m\sim k\log P\sim\sqrt{\log m}

so that our modular exponentiation can be completed in $O(\sqrt{\log m})$ steps. ∎

Remark 4.6.

This shows that in general it is best to choose all $t_{i}$ as $O(1)$ . When faced with the problem of optimizing these parameters, we can hence choose a relatively strong upper bound. On the other hand, it seems very difficult to compute the exactly optimal multiset $\mathcal{T}$ . Note that this is important because a smart choice of $t_{i}$ can make for a big performance boost (see Figure 3). If one could better understand the constant factors at play, then, modifying the above calculations, they could theoretically find the optimal choice of $\mathcal{T}$ .

In other words, when $m$ has large exponents relative to its prime factors, our algorithm makes large improvements (as $\log\operatorname{rad}m\ll\log m$ in the case where the exponents are on the order of $\log P$ ). In particular, our algorithm does well for prime powers. We have the following cororollary:

Cororollary 4.7.

Let $m=p^{\ell}$ where $\ell=O(\log p)$ . Let $a,n$ be positive integers such that $\gcd(a,m)=1$ . We may then compute $a^{n}\pmod{m}$ in $O(\log p)$ steps.

Proof.

This trivially follows from Corollary 4.5 with $k=1$ . ∎

This is quite impressive as, if we pick $\ell$ on the order of $\log p$ , we may modular exponentiate modulo $p^{\ell}$ in the same number of steps (asymptotically) as modulo $p$ . Furthermore, because each operation (modular multiplication or addition) can be taken to be $O((\log m)^{2})$ , we have the following result:

Corollary 4.8.

Let the modulus $m$ be a positive integer with known factorization $p_{1}^{e_{1}}p_{2}^{e_{2}}\cdots p_{k}^{e_{k}}$ such that all primes $p_{i}$ have the same bit length as an integer $P$ , and all $e_{i}$ are such that $e_{i}\sim k\log P$ . If the standard modular exponentiation algorithm takes $T$ time, our algorithm for the same values takes $cT^{5/6}$ time, for some constant $c$ and given unit of time.

Proof.

For convenience, denote $O((\log m)^{e})$ simply as $f(e)$ . The standard algorithm takes $f(1)$ steps and hence $f(3)$ time. Our algorithm takes $f(1/2)$ steps and hence $f(5/2)$ time. Then, if $T_{1}$ is the time it takes our algorithm to run and $T_{2}$ is the time for the standard algorithm to run, we have for some constants $c_{1},c_{2}$ that

(\log m)^{5/2}\sim c_{1}T_{1},

(\log m)^{3}\sim c_{2}T_{2}.

Hence $T_{1}\sim c_{2}^{5/6}c_{1}^{-1}T_{2}^{5/6}$ , so we obtain the desired result by taking the appropriate $c\sim c_{2}^{5/6}c_{1}^{-1}$ . ∎

Remark 4.9.

One may think that for such smooth $m$ , a CRT-type approach may also be fast. However, it is well-known (see e.g. [6]) that CRT is $O(\ell^{2})$ , where the moduli have bit length $\ell$ . Hence this is a non-issue.

5. General Modular Exponentiation

We may extend the ideas of our method to prove a more general theorem about matrix modular exponentiation, which translates over to linear recurrences. In this section, a step is a matrix multiplication.

Theorem 5.1.

Let $m=p_{1}^{e_{1}}p_{2}^{e_{2}}\cdots p_{k}^{e_{k}}$ be the canonical factorization of some $m\in\mathbb{N}$ . Choose parameters $t_{i}\in[1,e_{i}]$ for each $1\leq i\leq k$ . Let $d\in\mathbb{N}$ . Let $A\in GL_{d}(\mathbb{Z}/m\mathbb{Z})$ and $n\in\mathbb{N}$ . Then we may compute $A^{n}$ in $O\left(\max(e_{i}/t_{i})+\sum_{j=1}^{k}(t_{j}-1)d^{2}\log p_{j}+\sum_{j=1}^{k}\sum_{i=0}^{d-1}\log(p_{j}^{d}-p_{j}^{i})\right)$ steps.

Proof.

Let $T=\prod_{j=1}^{k}p_{j}^{t_{j}}$ . The algorithm is identical to that of Theorem 3.1, except we decompose the exponent modulo $|GL_{d}(\mathbb{Z}/T\mathbb{Z})|$ instead of $\varphi(T)$ . This works because of Lagrange’s theorem 2.8. The complexity is hence $O\left(\max(e_{i}/t_{i})+\log|GL_{d}(\mathbb{Z}/T\mathbb{Z})|\right)$ . By Lemma 2.10, this is the desired complexity. ∎

Remark 5.2.

A near-identical result holds for the special linear group instead. This is less practically useful but still important to note.

This asymptotic may be slightly difficult to work with, but a trivial corollary is the following:

Corollary 5.3.

Let $m=p_{1}^{e_{1}}p_{2}^{e_{2}}\cdots p_{k}^{e_{k}}$ be the canonical factorization of some $m\in\mathbb{N}$ . Choose parameters $t_{i}\in[1,e_{i}]$ for each $1\leq i\leq k$ . Let $d\in\mathbb{N}$ . Let $A\in GL_{d}(\mathbb{Z}/m\mathbb{Z})$ and $n\in\mathbb{N}$ . Then we may compute $A^{n}$ in $O\left(\max(e_{i}/t_{i})+d^{2}\sum_{i=1}^{k}t_{i}\log p_{i}\right)$ steps.

Proof.

This follows from Theorem 5.1 and the obvious estimate $\log(p_{j}^{d}-p_{j}^{i})<\log(p_{j}^{d})=d\log p_{j}$ . Indeed, the double sum is upper-bounded by $d^{2}\sum_{j=1}^{k}\log p_{j}$ and so the result follows. ∎

Remark 5.4.

With $d=1$ , we get Theorem 3.1 for prime powers. Also, this shows we can modular exponentiate matrices of size $O(1)$ in the same complexity as natural numbers.

Another corollary is that we may quickly compute the residue of large elements of a linear recurrent sequence modulo prime powers:

Corollary 5.5.

Let $m=p_{1}^{e_{1}}p_{2}^{e_{2}}\cdots p_{k}^{e_{k}}$ be the canonical factorization of some $m\in\mathbb{N}$ . Choose parameters $t_{i}\in[1,e_{i}]$ for each $1\leq i\leq k$ . Let $(u_{n})^{\infty}_{n=0}$ be a sequence of elements of $\mathbb{Z}/m\mathbb{Z}$ related with a linear recurrence relation of degree $d$ :

u_{n+d}=c_{d-1}u_{n+d-1}+\cdots+c_{0}u_{n}.

Suppose that $\gcd(c_{0},m)=1$ . Given sufficient initial terms, we may compute any element $u_{N}$ in $O\left(\max(e_{i}/t_{i})+\sum_{j=1}^{k}(t_{j}-1)d^{2}\log p_{j}+\sum_{j=1}^{k}\sum_{i=0}^{d-1}\log(p_{j}^{d}-p_{j}^{i})\right)$ steps.

Proof.

Linear recurrence relations of order $d$ can be represented as matrix powers. In particular, the respective matrices are in $GL_{d}(\mathbb{Z}/m\mathbb{Z})$ , so the result follows immediately by Theorem 5.1. One must make sure that the matrix produced has determinant invertible in $\mathbb{Z}/m\mathbb{Z}$ (and hence in $\mathbb{Z}/T\mathbb{Z}$ ). This determinant has magnitude $c_{0}$ , so we require $\gcd(c_{0},m)=1$ . This condition is met, so we are done. ∎

Remark 5.6.

This shows that we may compute the $n$ th Fibonacci number modulo $m$ in the same number of steps as we perform a modular exponentiation modulo $m$ . For example, modulo $p^{k}$ we may do it in $O(k+\log p)$ steps.

Fiduccia [7] and, more recently, Bostan and Mori [8] provide the state-of-the-art results for this problem in the case of a sequences over a general ring. The amount of steps taken is on the order of $M(d)\log n$ , where $M(d)=O(d\log d\log\log d)$ is the number of operations to multiply two polynomials in the ring. We achieve a stronger bound for the case where the ring is $\mathbb{Z}/m\mathbb{Z}$ . Our bound does not depend on $n$ because of the reduction we make modulo $|GL_{d}(\mathbb{Z}/T\mathbb{Z})|$ . In order to reduce $n$ in the same manner modulo $m$ for the bounds given by Fiduccia and by Bostan and Mori, the best possible reduction is by $|GL_{d}(\mathbb{Z}/m\mathbb{Z})|$ . By Lemma 2.10, $M(d)\log|GL_{d}(\mathbb{Z}/m\mathbb{Z})|\approx M(d)d^{2}\log m\gg O(d^{3}\log d\log m)$ . If $\omega\in[2,3]$ is the exponent for matrix multiplication, we take under $O(d^{\omega}\max(e_{i}/t_{i})+d^{2+\omega}\log T)$ steps of the same complexity as the steps taken by Fiduccia. Therefore Fiduccia’s algorithm is better as a function of $d$ , but we are better as a function of $m$ . It is often the case that $m\gg d$ , so we make significant improvement. Indeed, the key optimization that Fiduccia makes (involving the characteristic polynomial) only affects the complexity as a function of $d$ .

We can also apply our algorithm to ring extensions:

Theorem 5.7.

Let $p$ be a prime and $k$ and $n$ be natural numbers. Consider a finite ring extension $R=\mathbb{Z}/p^{k}\mathbb{Z}[\alpha_{1},\alpha_{2},\cdots,\alpha_{n}]$ . Consider the corresponding field extension $F=\mathbb{F}_{p}[\alpha_{1},\alpha_{2},\cdots,\alpha_{k}]$ . We can exponentiate in $R$ in $O(k+\log|F|)$ steps, where each step is an operation in $R$ .

Proof.

The algorithm is exactly that of Theorem 3.1 for prime powers and $t=1$ , except that we decompose the exponent modulo $|F|$ . This works due to Lagrange’s Theorem 2.8. ∎

A nice example of this theorem is that we provide fast modular exponentiation for Gaussian Integers modulo prime powers! This case is also related to Theorem 5.1 due to the bijection $a+bi\leftrightarrow\begin{pmatrix}a&-b\\ b&a\end{pmatrix}$ .

6. Programmatical Analysis

In this section, we test only integer modular exponentiation of our algorithm, as we prove that the generalizations have similar complexities.

Let us begin by testing our algorithm for $m=p^{k}$ . We will test against Python’s built-in $\mathrm{pow}$ function. We will iterate over primes $p$ , and choose $k$ randomly in a small interval around $\log p$ . In particular, we choose $k$ uniformly at random in the interval $[\log p-\sqrt{\log p},\log p+\sqrt{\log p}]$ . We choose $a\in[p^{k}/2,p^{k}]$ randomly (such that $\gcd(a,p)=1$ ) as well. We choose $n$ randomly in the interval $[p^{k}/2,p^{k}]$ . We choose $t_{i}=1$ for simplicity. We then compare the runtime for computation of $a^{n}\pmod{p^{k}}$ via our method ( $\mathcal{R}_{1}$ ) and Python’s built-in $\mathrm{pow}$ function ( $\mathcal{R}_{2}$ ) over $1000$ iterations.

In Figure 1 we iterate over $0\leq n\leq 3.5\cdot 10^{4}$ , plotting $p_{n+0.7\cdot 10^{5}}$ versus the respective ratio $\frac{\mathcal{R}_{2}}{\mathcal{R}_{1}}$ .

Refer to caption — Figure 1. Ratio for small primes

In Figure 2 we iterate over $0\leq n\leq 3.5\cdot 10^{4}$ , plotting $p_{n+3.5\cdot 10^{5}}$ versus the respective ratio $\frac{\mathcal{R}_{2}}{\mathcal{R}_{1}}$ .

As seen in the figures, there is a high variance in the ratio over small primes, whereas it steadies out for larger primes. The large jump in the second figure at about $n=16000$ is likely because Python has different optimizations for smaller calls of the $\mathrm{pow}$ function. We still do not have an explanation for the random jumps in data. Nonetheless, a basic implementation of our algorithm makes significant improvements to the highly optimized $\mathrm{pow}$ function.

The below figure shows how the problem of optimizing $\mathcal{T}$ is important.

For a specific choice of computing $a^{n}\pmod{p^{e}}$ with $p=101,e=200,a=13,$ and $n=\lfloor\frac{p^{e}}{3}\rfloor$ (chosen arbitrarily), we vary the choice of $\mathcal{T}$ over the set $\{\{1\},\{2\},\cdots,\{50\}\}$ . We plot this against the runtime for computation to create the blue curve. As seen in the figure, the optimal value of $\mathcal{T}$ yields a time approximately $5$ times faster than $\mathcal{T}=\{1\}$ . The orange curve demonstrates how the graph follows a curve of the form $at+b/t$ , as indicated by Corollary 4.5. The $R^{2}$ value is $96.5\%$ .

Recall that we achieve a complexity of $O(\sqrt{\log m})$ for an infinite family of $m$ , by Corollary 4.5. We aim to show this empirically. Because the complexity of the repeated squaring algorithm is $O(\log m)$ , if we graph the ratio $r=\mathcal{R}_{2}/\mathcal{R}_{1}$ of the runtime of python’s built-in modular exponentiation function to ours, we expect to see $r\propto\sqrt{\log m}$ . With $y=r$ and $x=\log m$ , we anticipate a graph of the form $y=c\sqrt{x}$ .

The setup is as follows. Let $P(n)$ be the first prime greater than $10^{n}$ . For $10\leq n<50,$ we let $p=P(n)$ and $e=n$ . Pick $t_{i}=1$ for simplicity. Then, $m=p^{e}$ is in the desired family. For each $n$ , we compute the ratio $r$ and plot it against $\log_{10}m\approx n^{2}$ .

The desired $y=c\sqrt{x}$ curve appears. Interestingly, the ratio for even $n$ is on average $24.8$ percent higher than for odd $n$ . The curves fit quite well: the $R^{2}$ values for even and odd $n$ are $94.9$ and $85.2$ percent respectively.

One might notice three values of odd $n$ appear to fit on the curve for even $n$ . Those values are $21$ , $39$ , and $45$ . We do not know why this phenomenon occurs or why the ratio is higher in general for even $n$ .

Remark 6.1.

In Python, we have achieved computation time more than $200$ times faster than the built-in $\mathrm{pow}$ for specific large values of $m$ .

See [2] for the Python code used to create these graphs.

7. Conclusion

We presented a fast algorithm for modular exponentiation when the modulus is known. We also presented a variant of this algorithm which uses less memory. We analyzed this algorithm in the general case, then shifted our focus to the specific case where the modulus has large prime exponents. We showed particular interest in the case where the modulus is a prime power, and we analyzed this case programmatically, testing it against Python’s built-in $\mathrm{pow}$ function. We also presented a stronger version of our algorithm for matrix modular exponentiation, which applies to the computation of large terms in linear recurrent sequences modulo some $m$ .

This algorithm has potential practical use in cryptography. Fast modular exponentiation is vital in the fast encryption of classic algorithms such as RSA and the Diffie-Hellman Key Exchange. It is even used in quantum algorithms: modular exponentiation is the bottleneck of Shor’s algorithm. If one could construct a cryptosystem in which it is useful to have a known modulus with large prime exponents, our algorithm would be applicable to its encryption process. For example, a variant of Takagi’s cryptosystem [9] with larger exponents has such properties. Additionally, work has been done on using matrix exponentiation and linear recurrences for error-correcting codes. For example, Matsui’s 2012 paper [10] uses linear sequences for a decoding algorithm. It is quite possible that our algorithm is potentially useful for such an algorithm.

There are a couple of things that we wish to do with this work going forward. We’d like to find a framework for programmatically testing general moduli (not only prime powers). Additionally, we hope to to make further progress on the front of optimizing $\mathcal{T}$ in practice. Furthermore, we want to come up with explanations for some of the phenomena that we see in the figures in section 5. Finally, we aim to implement further optimizations to our algorithm such as Montgomery Reduction.

8. Acknowledgments

We thank Dr. Simon Rubinstein-Salzedo for the useful discussion we had during the creation of this paper. We also thank Eva Goedhart and Nandita Sahajpal for awarding us the Lehmer Prize for this work at the 2023 West Coast Number Theory conference. We finally thank Dr. John Gorman from Jesuit High School Portland for inspiring this paper.

An Elementary Method For Fast Modular Exponentiation With Factored Modulus

Abstract.

1. Introduction

2. Preliminaries

Definition 2.1.

Definition 2.2.

Theorem 2.3.

Theorem 2.4.

Theorem 2.5.

Theorem 2.6.

Theorem 2.7.

Theorem 2.8.

Lemma 2.9.

Proof.

Lemma 2.10.

Proof.

Remark 2.11.

3. The Algorithm

Theorem 3.1.

Definition 3.2.

Example 1.

Lemma 3.3.

Proof.

Lemma 3.4.

Proof.

Proof.

Example 2.

Example 3.

Theorem 3.5.

Proof.

3.1. Pseudocode

4. Mathematical Analysis

Theorem 4.1.

Theorem 4.2.

Proof.

Theorem 4.3.

Proof.

Theorem 4.4.

Corollary 4.5.

Proof.

Remark 4.6.

Cororollary 4.7.

Proof.

Corollary 4.8.

Proof.

Remark 4.9.

5. General Modular Exponentiation

Theorem 5.1.

Proof.

Remark 5.2.

Corollary 5.3.

Proof.

Remark 5.4.

Corollary 5.5.

Proof.

Remark 5.6.

Theorem 5.7.

Proof.

6. Programmatical Analysis

Remark 6.1.

7. Conclusion

8. Acknowledgments

References