Montgomery Reduction for Gaussian Integers ^†

Abstract

Modular arithmetic over integers is required for many cryptography systems. Montgomery reduction is an efficient algorithm for the modulo reduction after a multiplication. Typically, Montgomery reduction is used for rings of ordinary integers. In contrast, we investigate the modular reduction over rings of Gaussian integers. Gaussian integers are complex numbers where the real and imaginary parts are integers. Rings over Gaussian integers are isomorphic to ordinary integer rings. In this work, we show that Montgomery reduction can be applied to Gaussian integer rings. Two algorithms for the precision reduction are presented. We demonstrate that the proposed Montgomery reduction enables an efficient Gaussian integer arithmetic that is suitable for elliptic curve cryptography. In particular, we consider the elliptic curve point multiplication according to the randomized initial point method which is protected against side-channel attacks. The implementation of this protected point multiplication is significantly faster than comparable algorithms over ordinary prime fields.

1. Introduction

Montgomery reduction is an efficient method that performs modulo reduction after an integer multiplication. The reduction algorithm uses multiplication and simple bit-wise operations instead of an expensive division [1]. Public-key cryptography based on the Rivest-Shamir-Adleman (RSA) system benefits from the efficient Montgomery modulo reduction [2,3]. Similarly, elliptic curve cryptography (ECC) systems over prime fields apply Montgomery reduction to support arbitrary prime curves [4,5,6,7,8,9,10,11,12,13,14].

In this work, we consider the modulo reduction for Gaussian integers. Gaussian integers are complex numbers such that the real and imaginary parts are integers. Arithmetic over Gaussian integers for the RSA system was proposed in [15,16,17,18]. In [17,18], it was shown for the RSA system that a significant complexity reduction can be achieved with Gaussian integers. Due to the isomorphism between Gaussian integer rings and ordinary integer rings, this arithmetic is suitable for many cryptography systems. The Rabin cryptography system and elliptic curve cryptography over Gaussian integers were considered in [19,20,21,22]. Moreover, coding applications over Gaussian integers use Gaussian integer arithmetic [23,24,25,26,27,28].

To reduce the complexity of the modular reduction, we investigate Montgomery reduction algorithms for finite Gaussian integer rings. However, it is not trivial to generalize Montgomery reduction to Gaussian integers. The final step of the reduction algorithm is based on the total order of integers which does not exist for complex numbers. In [29], a Montgomery reduction algorithm was proposed that utilizes the absolute value to measure the size of a Gaussian integer. This study is an extension of our paper [29]. In this work, we present a new approach that aims on reducing the complexity of the reduction presented in [29]. In this algorithm, the absolute value is replaced by the Manhattan weight. The calculation of the Manhattan weight requires only a single addition, whereas the calculation of the absolute value requires addition and two squaring operations. On the other hand, the reduction using the Manhattan weight may not always obtain a unique solution. There are at most two possible solutions that are congruent. For many applications uniqueness is not an issue for the intermediate results in calculations, which require many subsequent reduction steps. For the final result, the correct representative can be calculated using absolute values.

The proposed concept of Montgomery reduction was applied in [22] for the efficient calculation of the elliptic curve point multiplication. In particular, an area-efficient coprocessor was designed in [22]. This processor uses the proposed Montgomery modular arithmetic over Gaussian integers. It is shown in [21,22] that a key representation with a Gaussian integer expansion is beneficial for the calculation of the point multiplication. This key representation reduces the computational complexity and the memory requirements of a secure hardware implementation, which is protected against side-channel attacks [30,31,32,33,34,35,36,37]. In this work, we provide the theoretical justification for the reduction algorithms applied in [22]. Furthermore, we present new results for protected implementations of the point multiplication. In particular, we consider the randomized initial point method proposed in [38]. This method is protected against several side-channel attacks, i.e., differential power analysis (DPA), zero-value point attacks (ZPA), and refined power analysis (RPA) [32,38]. We present synthesis results for a field-programmable gate array (FPGA) based on the architecture proposed in [22]. The protected implementation of the point multiplication is significantly faster with Gaussian integers. Moreover, it requires less memory and fewer flip-flops than the PM algorithms over ordinary prime fields reported in [39,40].

This publication is organized as follows. We review the Montgomery multiplication of ordinary integers and discuss some properties of Gaussian integers in Section 2. In Section 3, Montgomery reduction based on the absolute value is investigated. We present a low-complexity reduction algorithm based on the Manhattan weight in Section 4. To demonstrate that the arithmetic over Gaussian integers is useful, we consider the complexity for the elliptic curve point multiplication according to the randomized initial point method in Section 5. In Section 6, we discuss the results and conclude our work.

2. Preliminaries and Problem Statement

In this section, we briefly discuss Montgomery reduction of ordinary integers. Moreover, we review some properties of Gaussian integer rings.

2.1. Montgomery Reduction of Ordinary Integers

Montgomery reduction is considered in several ECC and RSA cryptography systems to decrease the complexity of the arithmetic in prime fields [9,11,14] as well as in rings [2,3]. The Montgomery multiplication uses the arithmetic over a ring ${\mathcal{R}}_n$ that is isomorphic to the integer ring ${\mathbb{Z}}_n=\left\{xmodn:x=0,\dots ,n-1,x\in \mathbb{Z}\right\}$ [1]. Montgomery reduction algorithm reduces the complexity of the modulo reduction. The expensive calculation $modn$ is replaced by $modR$, where $R>n$ is a power of two. Hence, the modulo reduction in the ring ${\mathcal{R}}_n$ is a simple bitwise AND operation with $R-1$.

For the Montgomery arithmetic, an element $x\in {\mathbb{Z}}_n$ is mapped to the Montgomery domain as

$$X=xRmodn.$$

(1)

Addition in the Montgomery domain is isomorphic to ordinary integer modular arithmetic, because

$$(X+Y)modn=(xR+yR)modn=(x+y)Rmodn.$$

(2)

The multiplication in the Montgomery domain needs Montgomery reduction function $\mu \left(Z\right)=Z{R}^{-1}modn$ described in Algorithm 1. This function defines the inverse mapping, since $\mu \left(X\right)=xR{R}^{-1}modn=xmodn$. Using the reduction function, we obtain the product

$$\mu \left(XY\right)=XY{R}^{-1}modn=xyRmodn$$

(3)

in the Montgomery domain.

Algorithm 1 Montgomery reduction according to [1].

input: Z, with $0\le Z<nR$, $n^{\prime }=-n^{-1}modR$, and $R=2^l\ge n$

output:$M=\mu \left(Z\right)=Z{R}^{-1}modn$

1: $t=Z{n}^{\prime }modR$           // bitwise AND with $R-1$ 2: $q=(Z+tn)\mathrm{div}R$              // shift right by l 3: if ($q\ge n$) then 4:  $M=q-n$ 5: else 6:  $M=q$ 7: end if

Typically, all variables are mapped at the beginning of the calculation into the Montgomery domain using this reduction function as

$$X=\mu \left(x{R}^2\right)=x{R}^2{R}^{-1}modn=xRmodn,$$

(4)

since only one precomputed value $R^{2}modn$ is required.

We illustrate the Montgomery reduction of ordinary integers with an example. The binary representation of positive integers is denoted by brackets ${(·)}_2$ with the subscript 2.

Example 1.

We consider the product of two numbers $x=14$ and $y=7$ from ${\mathbb{Z}}_{29}$. The two integers x and y are mapped to the integers $X=xRmodn=13$ and $Y=yRmodn=21$ with $R=32=2^5>n=29$ in the Montgomery domain. Hence, all elements of the Montgomery domain can be represented with five binary digits. The product $xymodn=98mod29=11$ corresponds to $xyRmodn=4$ in the Montgomery domain.

With Algorithm 1, we can reduce the product $Z=XY=273$ as follows. With $n^{\prime }=11$, we first calculate $t=Z{n}^{\prime }modR=3003mod32=27$. Note that the modulo reduction of the binary representation is a simple bitwise AND operation with $R-1=31$ or equivalently the truncation of all bits except the five least significant bits. The integer 3003 has the binary representation ${\left(101110111011\right)}_2$ and the five least significant bits ${\left(11011\right)}_2$ represent the value $t=27$.

Next, we calculate the quotient $q=(Z+tn)\mathrm{div}R=1056\mathrm{div}32=33$. For $R=32$, the division by R without remainder in the binary representation is a simple right-shift by five bits or equivalently the truncation of the five least significant bits. The integer 1056 has the binary representation ${\left(10000100000\right)}_2$. The division without remainder results in ${\left(100001\right)}_2$ which is the binary representation of 33. The final result in the Montgomery domain is $M=\mu \left(Z\right)=q-29=4$ because $q=33>n=29$. To obtain the result in ${\mathbb{Z}}_{29}$, we apply Montgomery reduction for the inverse mapping $M^{\prime }=\mu \left(M\right)=\mu \left(4\right)=11$.

2.2. Rings of Gaussian Integers

The modulo arithmetic over Gaussian integers is based the modulo function

$$xmod\pi =x-\left[\frac{x{\pi }^{*}}{\pi {\pi }^{*}}\right]·\pi ,$$

(5)

where x and $\pi$ are Gaussian integers. ${\pi }^{*}$ is the conjugate of the Gaussian integer $\pi$. The brackets $\left[·\right]$ denote rounding to the closest Gaussian integer [23], i.e., for a complex number $x=c+d\mathrm{i}$, we have $\left[x\right]=\left[c\right]+\left[d\right]\mathrm{i}$. The set

$${\mathcal{G}}_p=\left\{xmod\pi :x=0,\dots ,p-1,x\in \mathbb{Z}\right\}$$

(6)

is a finite ring. For primes p of the form $p\equiv 1mod4$, ${\mathcal{G}}_p$ is a finite field which is isomorphic to the prime field $GF\left(p\right)$ over ordinary integers [23]. Such fields are suitable for elliptic curve cryptography [21,22]. A prime of the form $p\equiv 1mod4$ is the sum of two perfect squares, i.e., $p=\pi {\pi }^{*}={\left|a\right|}^2+{\left|b\right|}^2$ with the Gaussian integer $\pi =a+b\mathrm{i}$.

Moreover, for integers $n=cd$ such that c and d are primes of the form $c\equiv d\equiv 1mod4$ and $c\ne d$, the set ${\mathcal{G}}_n$ is a ring that is isomorphic to the ring ${\mathbb{Z}}_n$ over ordinary integers [26]. Such Gaussian integers are suitable for the RSA public-key system [17,18].

Consider Montgomery reduction in Algorithm 1. Lines 3 to 7 of this algorithm uniquely determine the smallest integer that is congruent to $Z{R}^{-1}modn$. However, complex numbers cannot be totally ordered [24]. Hence, it is not possible to directly apply this reduction algorithm to Gaussian integers. This reduction problem is not regarded in [17,18].

In this work, we consider two reduction strategies using different norms. One algorithm uses the absolute value $\left|x\right|=\sqrt{{\left|c\right|}^2+{\left|d\right|}^2}$ of the Gaussian integer $x=c+d\mathrm{i}$. The second approach is based on the Manhattan weight $∥x∥=\left|c\right|+\left|d\right|$. Calculating the Manhattan weight is less complex than calculating the absolute value, because only addition is required, whereas squaring is necessary to determine $\left|x\right|$. Note that

$$\left|x\right|\le ∥x∥$$

(7)

holds for any x which follows from squaring both sides of the inequality.

2.3. Finding Primes of the Form $p=a^2+b^2$

An algorithm for primes of the form $p\equiv 1mod4$ to find $a,b$ such that $p=a^2+b^2$ is proposed in [23]. This algorithm consists of two steps.

• Find x such that $x^2\equiv -1modp$, which can be done using the Tonelli-Shanks algorithm [41].
• Use the Euclidean algorithm to determine the greatest common divisor of p and x. The first two remainders less than $\sqrt{p}$ are a and b.

On the other hand, there exist many primes of the form $p\equiv 1mod4$ such that $p=a^2+{(a-1)}^2$ or $p=a^2+1$. Exploiting this observation we can search for a such that the sum $p=a^2+1$ or $p=2a^2-2a+1$ is prime. This can significantly reduce the search complexity.

Table 1. Examples for primes of the form $p=a^2+b^2$ with $b=1$ or $b=a-1$.

bits	p	a	b
156	8 $\times {10}^{46}+74\times {10}^{24}+17113$	$2\times {10}^{23}+93$	$2\times {10}^{23}+92$
169	$4\times {10}^{50}+216\times {10}^{25}+2917$	$2\times {10}^{25}+54$	1
170	$8\times {10}^{50}+6\times {10}^{26}+113$	$2\times {10}^{25}+8$	$2\times {10}^{25}+7$
190	$8\times {10}^{56}+3\times {10}^{30}+2813$	$2\times {10}^{28}+38$	$2\times {10}^{28}+37$
256	$8\times {10}^{76}+2516\times {10}^{38}+197821$	$2\times {10}^{38}+315$	$2\times {10}^{38}+314$
382	$9\times {10}^{114}+384\times {10}^{57}+4097$	$3\times {10}^{57}+64$	1

illustrates some examples of such primes with sufficient bit lengths for ECC applications.

3. Montgomery Reduction for Gaussian Integers

In this section, we consider some important properties of Gaussian integers and derive Montgomery reduction for Gaussian integers using the absolute value. First, we show that the set ${\mathcal{G}}_n$ is symmetric and that the absolute value of any element of this ring is bounded.

Lemma 1.

For any $x\in {\mathcal{G}}_n$ we have the following symmetry

$$-x,\mathrm{i}x,-\mathrm{i}x\in {\mathcal{G}}_n.$$

(8)

The absolute value $\left|x\right|$ is bounded by

$$\left|x\right|<\frac{\left|\pi \right|}{\sqrt{2}}.$$

(9)

Moreover, for any $x^{\prime }\notin {\mathcal{G}}_n$ with $x=x^{\prime }mod\pi$ we have

$$\left|x\right|<\left|x^{\prime }\right|.$$

(10)

Proof. 

Consider the quotient $\frac{x}{\pi }=\frac{x{\pi }^{*}}{n}=c+d\mathrm{i}$, where c and d are the real part and the imaginary part of $\frac{x}{\pi }$. For $x\in {\mathcal{G}}_n$ we have $x=xmod\pi$. Hence, the rounded quotient satisfies

$$\left[\frac{x{\pi }^{*}}{\pi {\pi }^{*}}\right]=0.$$

This implies $\left[c\right]=\left[d\right]=0$ and $\left|c\right|<1/2$, $\left|d\right|<1/2$. Hence, we can bound the absolute value of the quotient $\frac{x}{\pi }$ as

$$\left|\frac{x}{\pi }\right|<\frac{1}{\sqrt{2}}.$$

Multiplying both sides by $\left|\pi \right|$ results in (9). Note that $x,-x,ix$, or $-ix$ have the same absolute value. Hence, (8) holds.

Finally, consider a Gaussian integer $x^{\prime }\notin {\mathcal{G}}_n$ which is congruent to x, i.e., $x=x^{\prime }mod\pi$. We use the notation $\frac{x^{\prime }}{\pi }=c^{\prime }+d^{\prime }\mathrm{i}$. Note that $x=x^{\prime }mod\pi$ implies $c=c^{\prime }-\left[c^{\prime }\right]$ and $d=d^{\prime }-\left[d^{\prime }\right]$. Moreover, at least one rounded value $\left[c^{\prime }\right]$ or $\left[d^{\prime }\right]$ is nonzero because $x^{\prime }\notin {\mathcal{G}}_n$, where $\left[c^{\prime }\right]\ne 0$ results in $\left|c\right|<1/2<\left|c^{\prime }\right|$. Similarly, $\left[d^{\prime }\right]\ne 0$ implies $\left|d\right|<1/2<\left|d^{\prime }\right|$. Hence, we have (10). □

The next example demonstrates that $\left|x\right|<\frac{\left|\pi \right|}{\sqrt{2}}$ is not a sufficient condition for $x\in {\mathcal{G}}_n$.

Example 2.

We consider the finite field ${\mathcal{G}}_{29}$ for the Gaussian integer $\pi =5+2\mathrm{i}$. All elements of this field are depicted in Figure 1. Consider $q^{\prime }=3$ with $q=3mod\pi =-2-2\mathrm{i}$, i.e., $q^{\prime }\notin {\mathcal{G}}_{29}$ and $q\in {\mathcal{G}}_{29}$. Both Gaussian integers satisfy $\left|q^{\prime }\right|<\frac{\left|\pi \right|}{\sqrt{2}}$ and $\left|q\right|<\frac{\left|\pi \right|}{\sqrt{2}}\approx 3.8079$. Hence, it is not possible to determine the representative based on the bound (9). The representative $q=-2-2\mathrm{i}$ follows from (10), as $\left|q\right|\approx 2.8284<\left|q^{\prime }\right|=3$.

Next, we consider now Montgomery reduction. The reduction for Gaussian integers can be performed according to Algorithm 2, where

$$\begin{array}{ccc}\hfill {\alpha }^{\prime }& =& \underset{\alpha \in \{0,\pm 1,\pm \mathrm{i}\}}{argmin}\left|q-\alpha \pi \right|,\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\alpha }^{″}& =& \underset{\alpha \in \{\pm 1,\pm \mathrm{i},\pm 1\pm \mathrm{i}\}}{argmin}\left|q-\alpha \pi \right|.\hfill \end{array}$$

(12)

Steps 1 and 2 in Algorithm 2 are applied separately for the real- and imaginary part. We demonstrate that Algorithm 2 always obtains a precision reduction resulting in the correct representative.

Algorithm 2 Montgomery reduction for Gaussian integers

input:$Z=XY$, ${\pi }^{\prime }=-{\pi }^{-1}modR$, $R=2^l>\frac{\left|\pi \right|}{\sqrt{2}}$

output:$M=\mu \left(Z\right)=Z{R}^{-1}mod\pi$

1: $t=Z{\pi }^{\prime }modR$           // bitwise AND of $\mathrm{Re}$, $\mathrm{Im}$ with $R-1$ 2: $q=(Z+t\pi )\mathrm{div}R$                // shift $\mathrm{Re}$, $\mathrm{Im}$ right by l 3: if ($\left|q\right|<\frac{\sqrt{2}-1}{\sqrt{2}}\left|\pi \right|$) then 4:  $M=q$ 5: else if ($\left|q\right|<\frac{\left|\pi \right|}{\sqrt{2}}$) then 6:  determine ${\alpha }^{\prime }$ according to (11) 7:  $M=q-{\alpha }^{\prime }\pi$ 8: else 9:  determine ${\alpha }^{″}$ according to (12) 10:  $M=q-{\alpha }^{″}\pi$ 11: end if

Proposition 1.

For M and Z according to Algorithm 2, we have

$$M=Z{R}^{-1}mod\pi .$$

(13)

Proof. 

We consider the sum $Z+t\pi$ in step 2 of Algorithm 2, where

$$\begin{array}{ccc}\hfill Z+t\pi & \equiv & Z+Z{\pi }^{\prime }\pi modR\hfill \\ & \equiv & Z+Z(-{\pi }^{-1})\pi modR\hfill \\ & \equiv & Z-ZmodR\equiv 0,\hfill \end{array}$$

implies that R divides $Z+t\pi$. Considering the corresponding quotient $q=(Z+t\pi )\mathrm{div}R$, we have

$$\begin{array}{ccc}\hfill qmod\pi & \equiv & (Z+t\pi )\mathrm{div}Rmod\pi \hfill \\ & \equiv & (Z+t\pi )R^{-1}mod\pi \hfill \\ & \equiv & Z{R}^{-1}+t\pi R^{-1}mod\pi \hfill \\ & \equiv & Z{R}^{-1}mod\pi .\hfill \end{array}$$

This shows that q is congruent to $Z{R}^{-1}mod\pi$ or $Z{R}^{-1}=q-\alpha \pi$. The absolute value of $Z{R}^{-1}mod\pi$ is bounded, i.e.,

$$\left|Z{R}^{-1}\right|=\left|XY\right|R^{-1}\le \frac{{\left|\pi \right|}^2}{2R}<\frac{R\left|\pi \right|}{\sqrt{2}R}=\frac{\left|\pi \right|}{\sqrt{2}}.$$

As shown in Example 2, the quotient q may not be the correct representative $Z{R}^{-1}mod\pi$ even for $\left|q\right|<\frac{\left|\pi \right|}{\sqrt{2}}$. The congruent values $q-\alpha \pi$ with $\alpha \in \{\pm 1,\pm \mathrm{i}\}$ are also possible candidates. However, for any $x\in {\mathcal{G}}_n$ and $\alpha \in \{\pm 1,\pm \mathrm{i}\}$, the lower bound

$$\left|x-\alpha \pi \right|\ge \left|\left|\pi \right|-\left|x\right|\right|>\frac{\sqrt{2}-1}{\sqrt{2}}\left|\pi \right|,$$

follows from the triangle inequality and Lemma 1. Consequently, the quotient q is the unique solution for $\left|q\right|<\frac{\sqrt{2}-1}{\sqrt{2}}\left|\pi \right|$. Furthermore, for any $x\in {\mathcal{G}}_n$ the condition $\left|\alpha \right|>\sqrt{2}$ implies

$$\left|x-\alpha \pi \right|\ge \left|\left|\alpha \pi \right|-\left|x\right|\right|>\frac{\left|\pi \right|}{\sqrt{2}}.$$

Hence, only the candidates according to (11) are possible for $\frac{\sqrt{2}-1}{\sqrt{2}}\left|\pi \right|\le \left|q\right|<\frac{\left|\pi \right|}{\sqrt{2}}$. Using Lemma 1, it follows that the correct candidate can be found by minimizing the absolute value among all candidates according to (11).

If $\left|q\right|\ge \frac{\left|\pi \right|}{\sqrt{2}}$, the result M is calculated as $M=q-{\alpha }^{″}\pi$ with ${\alpha }^{″}$ according to (12). To demonstrate that (12) is correct, we derive the bound $\left|\alpha \right|\le \sqrt{2}$. Consider again the quotient $q=(Z+t\pi )\mathrm{div}R$. We aim to compensate the offset $t\pi R^{-1}$ by $\alpha \pi$, where the absolute values are bounded by

$$\left|\alpha \pi \right|=\left|t\pi R^{-1}\right|=\left|t\right|\left|\pi \right|R^{-1}\le \sqrt{2}\left|\pi \right|.$$

This follows from $\left|t\right|\le \sqrt{R^2+R^2}=\sqrt{2}R$ due to the reduction $modR$ and implies the bound $\left|\alpha \right|\le \sqrt{2}$. □

Example 3.

As an example for Montgomery reduction, we consider the product of two numbers $x=6$ and $y=7$ from ${\mathbb{Z}}_{29}$ or the corresponding field ${\mathcal{G}}_{29}$ of Gaussian integers with $\pi =5+2\mathrm{i}$. The two integers x and y are mapped to the Gaussian integers $X=xRmod\pi =2-\mathrm{i}$ and $Y=yRmod\pi =-2$ with $R=8$ in the Montgomery domain. The product $xymodp=42mod29=13$ corresponds to $xyRmod\pi =-\mathrm{i}$ in the Montgomery domain.

With Algorithm 2, we can reduce the product $Z=XY=-4+2\mathrm{i}$ as follows. With ${\pi }^{\prime }=-1+2\mathrm{i}$, we first calculate $t=Z{\pi }^{\prime }modR=-10\mathrm{i}mod8=-2\mathrm{i}$. Note that the modulo reduction is a simple bitwise AND operation with $R-1=7$.

Next, we calculate the quotient $q=(Z+t\pi )\mathrm{div}R=-8\mathrm{i}\mathrm{div}8=-\mathrm{i}$. For $R=8=2^3$, the division by R without remainder is a simple right-shift by three bits. The quotient $M=\mu \left(Z\right)=q=-\mathrm{i}$ is the final result without further reduction steps because $\left|q\right|=1<\frac{\sqrt{2}-1}{\sqrt{2}}\left|\pi \right|=1.5773$.

The product $xymodp=13$ corresponds to $xymod\pi =1+\mathrm{i}$ in the field ${\mathcal{G}}_{29}$. To obtain this result, we apply Montgomery reduction for the inverse mapping $M^{\prime }=\mu \left(M\right)=M{R}^{-1}$. With $M=-\mathrm{i}$, we have $t=Z{\pi }^{\prime }modR=2+\mathrm{i}mod8=2+\mathrm{i}$ and the result $M^{\prime }=q=(Z+t\pi )\mathrm{div}R=8+8\mathrm{i}\mathrm{div}8=1+\mathrm{i}$. Again, no reduction is required because the condition $\left|q\right|=\sqrt{2}<\frac{\sqrt{2}-1}{\sqrt{2}}\left|\pi \right|=1.5773$ is satisfied.

We illustrate the final reduction procedure of Algorithm 2 in the following example. There are up to eight possible values $\alpha$ according to (11) and (12). However, not all potential values have to be considered for the reduction. The number of valid candidates $\alpha$ can be reduced based on the signs of the real and imaginary part of q.

Example 4.

We consider again the finite field ${\mathcal{G}}_{29}$ of Gaussian integers with $\pi =5+2\mathrm{i}$. Figure 2 depicts all possible quotients q in the Montgomery domain after the first two steps of Algorithm 2. For $q=3$ we have $\frac{\sqrt{2}-1}{\sqrt{2}}\left|\pi \right|\approx 1.5773<\left|q\right|=3<\frac{\left|\pi \right|}{\sqrt{2}}\approx 3.8079$. Thus, ${\alpha }^{\prime }$ should be determined according to (11). However, there are only two possible values for α, i.e., $\alpha \in \{0,1\}$, as q is in the first quadrant. Calculating $\left|q-0\right|=3$ and $\left|q-\pi \right|\approx 2.3852$ leads to ${\alpha }^{\prime }=1$. Consequently, we obtain the representative $M=q-{\alpha }^{\prime }\pi =-2-2\mathrm{i}$.

Algorithm 2 achieves the desired precision reduction, where the final reduction step always results in the correct representative $M=qmod\pi$. This result satisfies the bound $\left|M\right|\le \frac{\left|\pi \right|}{\sqrt{2}}$. The number of possible candidates $\alpha$ is restricted depending on the absolute value of q, as illustrated in (11) and (12).

Nonetheless, the reduction according to Algorithm 2 can be rather complex due to the squaring to determine the absolute values. This complexity may not be required in all applications of the Montgomery multiplication. For applications in cryptography, it is adequate to find the unique representative $Z{R}^{-1}mod\pi$ only in the last calculation step, whereas for intermediates results a reduction that achieves $M\equiv Z{R}^{-1}mod\pi$ is sufficient. This motivates a low complexity Montgomery reduction based on the Manhattan weight, which is considered in the next section.

4. Precision Reduction for Gaussian Integers Using the Manhattan Weight

In this section, we present a low complexity precision reduction for Gaussian integers based on the Manhattan weight, i.e., the absolute value in the reduction algorithm is replaced by the Manhattan weight. The calculation of the absolute value requires squaring, whereas the Manhattan weight requires only addition. On the other hand, this algorithm may not obtain a unique solution. However, there are at most two possible solutions that are congruent. In ECC or RSA systems that require many reduction steps, this algorithm can be used for the intermediate results.

Without loss of generality we consider $\pi =a+b\mathrm{i}$ with $a>b\ge 1$. The precision reduction is described in Algorithm 3, where

$$\widehat{\alpha }=\underset{\alpha \in \{\pm 1,\pm \mathrm{i},\pm 1\pm \mathrm{i}\}}{argmin}∥q-\alpha \pi ∥,$$

(14)

and $N=a-1$. We demonstrated that this algorithm always obtains $M\equiv Z{R}^{-1}mod\pi$, where $∥M∥\le N$.

Algorithm 3 Precision reduction for Gaussian integers

input:$Z=XY$, ${\pi }^{\prime }=-{\pi }^{-1}modR$, $R=2^l>N$

output:$M\equiv \mu \left(Z\right)=Z{R}^{-1}mod\pi$

1: $t=Z{\pi }^{\prime }modR$           // bitwise AND of $\mathrm{Re}$, $\mathrm{Im}$ with $R-1$ 2: $q=(Z+t\pi )\mathrm{div}R$              // shift $\mathrm{Re}$ and $\mathrm{Im}$ right by l 3: if ($∥q∥\le N$) then 4:  $M=q$ 5: else 6:  determine $\widehat{\alpha }$ according to  (14) 7:  $M=q-\widehat{\alpha }\pi$ 8: end if

First, we note that using the symmetry according to Lemma 1, we can restrict the following derivations to the elements of the first quadrant. Next, we consider some important properties of the Manhattan weight.

Lemma 2.

For $x\in {\mathcal{G}}_p$ the Manhattan weight is upper bounded by

$$∥x∥\le a-1=N.$$

(15)

Proof. 

Let c and d be the real part and the imaginary part of $\frac{x}{\pi }$. For $x\in {\mathcal{G}}_p$ we have $x=xmod\pi$ and consequently

$$\left[\frac{x{\pi }^{*}}{\pi {\pi }^{*}}\right]=0.$$

(16)

This implies $\left[c\right]=\left[d\right]=0$ and $\left|c\right|<1/2$, $\left|d\right|<1/2$. Hence, we consider $(1/2+1/2\mathrm{i})\pi$ to upper bound the real and imaginary parts of x as

$$\begin{array}{ccc}\hfill \left|\mathrm{Re}\left\{\mathrm{x}\right\}\right|& <& \left|\frac{a-b}{2}\right|,\hfill \\ \hfill \left|\mathrm{Im}\left\{\mathrm{x}\right\}\right|& <& \left|\frac{a+b}{2}\right|.\hfill \end{array}$$

Note that either a is odd and b is even or vice versa because p is an odd prime. Furthermore, the real- and imaginary parts of x are integers. Consequently, we have

$$\begin{array}{ccc}\hfill \left|\mathrm{Re}\left\{\mathrm{x}\right\}\right|& \le & \left|\frac{a-b}{2}-\frac{1}{2}\right|,\hfill \\ \hfill \left|\mathrm{Im}\left\{\mathrm{x}\right\}\right|& \le & \left|\frac{a+b}{2}-\frac{1}{2}\right|.\hfill \end{array}$$

We can bound the Manhattan weight of x as

$$∥x∥\le ∥\frac{a-b-1}{2}+\frac{a+b-1}{2}\mathrm{i}∥.$$

With $a>b\ge 1$, the Manhattan weight of x is upper bounded by

$$∥x∥\le \frac{a-b-1}{2}+\frac{a+b-1}{2}=a-1.$$

(17)

Hence (15) holds. □

Next, we consider bounds on the Manhattan weight for the sum and product of two elements $x,y\in {\mathcal{G}}_p$. Note that we consider arithmetic without modulo reduction.

Lemma 3.

For $x,y\in {\mathcal{G}}_p$ we have the upper bounds

$$\begin{array}{ccc}\hfill ∥x+y∥& \le & 2N,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill ∥xy∥& \le & N^2,\hfill \end{array}$$

(19)

for arithmetic without modulo reduction.

Proof. 

The bound on the sum follows from the triangle inequality and (15), i.e.,

$$∥x+y∥\le ∥x∥+∥y∥\le 2N.$$

(20)

Without loss of generality we consider two elements $x=c+d\mathrm{i}$ and $y=e+f\mathrm{i}$ from the first quadrant for the product in (19). This implies $∥x∥=c+d\le N$ or $d\le N-c$. Similarly, we have $f\le N-e$. For the product $xy$ we have

$$xy=(ce-df)+(ed+cf)\mathrm{i}.$$

(21)

First, consider the absolute value of the imaginary part $\mathrm{Im}\left\{\mathrm{xy}\right\}$

$$\left|\mathrm{Im}\left\{\mathrm{xy}\right\}\right|=ed+cf\le e(N-c)+c(N-e)=eN+cN-2ce.$$

To determine the maximum value, we consider the bivariate function $g(e,c)=eN+cN-2ce$ and its partial derivatives

$$\begin{array}{ccc}\hfill \frac{{\partial }^{}g(e,c)}{\partial e}& =& N-2c=0,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \frac{{\partial }^{}g(e,c)}{\partial c}& =& N-2e=0.\hfill \end{array}$$

(23)

This results in a maximum for $c=e=N/2$ and the bound $\left|\mathrm{Im}\left\{\mathrm{xy}\right\}\right|=ed+cf\le N^2/2$. Due to symmetry this bound also holds for the absolute value of the real part. Hence, we have

$$∥xy∥\le \frac{N^2}{2}+\frac{N^2}{2}=N^2.$$

(24)

□

The following proposition demonstrates that Algorithm 3 results in the desired reduction.

Proposition 2.

For M and Z according to Algorithm 3, we have

$$\begin{array}{ccc}\hfill M& \equiv & Z{R}^{-1}mod\pi ,\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill ∥M∥& \le & N.\hfill \end{array}$$

(26)

Proof. 

The first two steps are identical in Algorithm 3 and Algorithm 2. Hence, the first statement (25) follows from the same arguments as in the proof of Proposition 1.

For $∥q∥\le N$, we immediately have (26). For $∥q∥>N$, we consider again the corresponding quotient $q=(Z+t\pi )\mathrm{div}R$. The Manhattan weight of $Z{R}^{-1}mod\pi$ is bounded by

$$∥\frac{Z}{R}∥=\frac{∥Z∥}{R}=\frac{∥XY∥}{R}\le \frac{N^2}{R}\le \frac{NR}{R}=N,$$

(27)

where we have used (24) and the assumption $R\ge N$.

Similar to the proof of Proposition 1, we aim to compensate $t\pi R^{-1}$ with $\alpha \pi$. From Proposition 1 follows that $∥\alpha ∥\le \sqrt{2}$. Hence, the minimization in (14) over $\alpha \in \{\pm 1,\pm \mathrm{i},\pm 1\pm \mathrm{i}\}$ is sufficient. From (27) follows that there is at least one solution with $∥M∥=∥q-\widehat{\alpha }\pi ∥\le N$. Consequently, the minimization in (14) will find a solution satisfying (26). □

Finally, we can now describe the final reduction step of Algorithm 3. Note that there are eight possible values for $\alpha$ according to (14), but not all potential values have to be considered. The reduction procedure is demonstrated in the following example.

Example 5.

Again, we consider the finite field ${\mathcal{G}}_{29}$ of Gaussian integers with $\pi =5+2\mathrm{i}$, where all possible values of q after the first two steps of Algorithm 3 are depicted in Figure 2.

For instance, consider the product of $x=19$ and $y=7$ from ${\mathbb{Z}}_{29}$. The two integers x and y are mapped to the Gaussian integers $X=xRmod\pi =2-2\mathrm{i}$ and $Y=yRmod\pi =-2$ with $R=8$ in the Montgomery domain. The product $xymodp=133mod29=17$ corresponds to $xyRmod\pi =3-\mathrm{i}$ in the Montgomery domain.

With Algorithm 2, we can reduce the product $Z=XY=-4+4\mathrm{i}$ as follows. We calculate $t=Z{\pi }^{\prime }modR=4+4\mathrm{i}$ and the quotient $q=(Z+t\pi )\mathrm{div}R=1+4\mathrm{i}$. A valid solution of the reduction is a Gaussian integer M with $∥M∥\le N$. We have $∥q∥=5>N=4$. Hence, further reduction is required. As q is in the first quadrant, we have three possible values for α, i.e., $\alpha \in \{1,\mathrm{i},1+\mathrm{i}\}$. Calculating $q-\alpha \pi$ results in $-4+2\mathrm{i},3-\mathrm{i},-2-3\mathrm{i}$, where only the solution $M=3-\mathrm{i}$ satisfies the condition $∥M∥\le N=4$. Hence, we choose $M=3-\mathrm{i}$ as the final result.

The final reduction step results in a Gaussian integer x with $∥x∥\le N$. Hence, Algorithm 3 achieves the desired precision reduction using the Manhattan weight. However, the result may not always be an element of the ring. For instance, for some $\pi$ the value $x=a-1$ is a ring element. Alternatively, the point $x^{\prime }=x-\pi =-1-b\mathrm{i}$ can be a ring element. Both values are congruent and can satisfy $∥x∥=a-1\le N$, $∥x^{\prime }∥=b+1\le N$ depending on $\pi$. The ring element is the value with the smallest absolute value. There exist at most two congruent values with Manhattan weight less or equal N. Consequently, the correct representation can be selected by minimizing the absolute value of these two congruent values. Depending on the application this step might be required once to obtain the final result $x=qmod\pi$.

5. Elliptic Curve Point Multiplication

In this section, we consider the elliptic curve point multiplication to demonstrate that the arithmetic over Gaussian integers enables an efficient calculation. Calculations of the elliptic curve point multiplication are prone to side-channel attacks. For example, an attacker can estimate the secret key based on the required power consumption of the different arithmetic operations over time. There are many known attacks on the point multiplication like timing attacks, simple power analysis, differential power analysis, refined power analysis, and zero-value point attacks [32,38].

The concept for the point multiplication over Gaussian integers was introduced in [22], where an area-efficient coprocessor design was proposed with an arithmetic unit that enables Montgomery modular arithmetic over Gaussian integers. It is shown in [22] that a key representation with a Gaussian integer expansion is beneficial to reduce the computational complexity and the memory requirements of a secure hardware implementation considering timing attacks and simple power analysis. This architecture supports different point multiplication algorithms. In contrast to [22], we consider the randomized initial point method which is additionally protected against differential power analysis, refined power analysis, and zero-value point attacks [32,38]. We refer to [22] for the architecture and the implementation details.

In the following, we consider the elliptic curve

$$y^2=x^3+\alpha x+\beta ,$$

(28)

which is recommended for prime fields $GF\left(p\right)$ [32,42]. The parameters $\alpha$ and $\beta$ are constant coefficients. The pair x and y defines the coordinates of a point $P(x,y)$ on the curve. The one-way function of elliptic-curve cryptography is the point multiplication, i.e., $kP$, where P is a point on the elliptic curve and k is an integer. The point multiplication is typically implemented based on the binary expansion $k={\sum }_{i=0}^{r-1}{k}_i{2}^i$ of the integer k with binary digits $k_i\in \{0,1\}$, where r is the length in bits and $k_{r-1},k_{r-2},\dots ,k_0$ is the binary representation of the key (the integer k). This method requires two different point operations, the point addition (ADD) and the point doubling operation (DBL) [32]. Let $P(x_1,y_1)$ and $Q(x_2,y_2)$ be two distinct points on an elliptic curve (28), then the sum $S(x_3,y_3)=P(x_1,y_1)+Q(x_2,y_2)$ is calculated as

$$x_3={\left({\displaystyle \frac{y_2-y_1}{x_2-x_1}}\right)}^2-x_1-x_2,y_3=\left({\displaystyle \frac{y_2-y_1}{x_2-x_1}}\right)(x_1-x_3)-y_1.$$

(29)

Similarly, the point doubling operation $S(x_3,y_3)=2P(x_1,y_1)$ is defined as

$$x_3={\left({\displaystyle \frac{3x_1^2+\alpha }{2y_1}}\right)}^2-2x_1,y_3=\left({\displaystyle \frac{3x_1^2+\alpha }{2y_1}}\right)(x_1-x_3)-y_1.$$

(30)

The calculations in (29) and (30) are typically performed using arithmetic over prime fields. However, we consider arithmetic over Gaussian integer fields, where $x_i,y_i,\alpha ,\beta \in {\mathcal{G}}_p$.

As an alternative to the binary expansion key representation, $\tau$-adic expansions of the integer k with a non-binary basis $\tau$ were proposed to speed up the point multiplication (PM) and to improve the resistance against attacks [38,43,44,45,46,47]. The $\tau$-adic expansion results in the point multiplication

$$kP=\sum _{i=0}^{l-1}{\kappa }_i{\tau }^{i}P=\tau (\dots \tau (\tau {\kappa }_{l-1}P+{\kappa }_{l-2}P)+\dots )+{\kappa }_{0}P,$$

(31)

where we consider digits ${\kappa }_i$ from a Gaussian integer field.

To demonstrate that a $\tau$-adic expansion with arithmetic over Gaussian integers can reduce the computational complexity, we consider an example based on the curve (28) with $\beta =0$. For the products ${\kappa }_{i}P(x,y)$ we use the endomorphism $\mathrm{i}P(x,y)=P(-x,\mathrm{i}y)$ for prime fields [32,48]. Note that $P(-x,\mathrm{i}y)$ is a point on the curve (28) if $P(x,y)$ is a valid point. The negation of any point $-P(x,y)$ is $P(x,-y)$ according to [32]. Hence, for the products ${\kappa }_{i}P(x,y)$ with ${\kappa }_i\in \{0,\pm 1,\pm \mathrm{i}\}$ we have

$$\begin{array}{ccc}\hfill P& =& P(x,y),\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill -P& =& P(x,-y),\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill \mathrm{i}P& =& P(-x,\mathrm{i}y),\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill -\mathrm{i}P& =& P(-x,-\mathrm{i}y).\hfill \end{array}$$

(35)

Using this endomorphism, we can calculate other products $\tau P$, e.g., $\tau P=(2+\mathrm{i})P=2P+\mathrm{i}P$ requires a DBL, the mapping $\mathrm{i}P$, and an ADD operation. Several other endomorphisms for different curves over prime fields can be fund in [32,48].

For the point multiplication, we consider the randomized initial point method according to Algorithm 4. This method introduces a random point R into the calculation of the point multiplication. Consequently, all intermediate results in the calculation of the point multiplication become randomized. The algorithm starts with a precomputation phase, where the points $S_{l-1},\dots ,S_1,S_0$ are computed in steps 1 to 4 and stored in a memory. The point addition in step 5 results in the point $Q=S_{l-1}+R$. This point is multiplied with $\tau$ and then added to the corresponding precomputed point $S_j$ in the loop in steps 6 to 9 of this algorithm. Due to the second product of the precomputations $-(\tau -1)·R$, we obtain the point $Q+R$ after each iteration of this loop. Correspondingly, the final result is the point $Q+R=k·P+R$. Hence, we obtain the correct result of the point multiplication $k·P$ by subtracting R from Q in step 10. The point addition in step 8 is computed in each iteration since the $\tau$-adic expansion of the key according to [22] excludes all zero elements. This prevents SPA and timing attacks. Furthermore, all stored points $S_{l-1},\dots ,S_1,S_0$ are randomized due to the subtracting the random point $-(\tau -1)·R$. Thus, the resistance against DPA, RPA, and ZPA is increased.

We can use the reduction based on the Manhattan weight in Algorithm 3 for all interim results in the point multiplication. The aim of the reduction of interim results is finding a Gaussian integer $\tilde{x}$ that has a Manhattan weight satisfying $∥\tilde{x}∥\le N$ and is congruent to the actual representative x. The reduction algorithm can be stopped once a value $\tilde{x}=q-\alpha \pi$ with $∥\tilde{x}∥\le N$ is found. Hence, not all offsets in (14) have to be considered in every reduction.

Table 2. Primes of the form $p=a^2+b^2$ suitable for elliptic curve cryptography (ECC) applications and the percentage of occurrences of offset reduction steps.

${log}_2\left(\mathit{p}\right)$	a	b	No Reduction	1 Reduction	2 Reductions	3 Reductions
188	$2^{94}-150$	1	$91.7\%$	$4.4\%$	$2.4\%$	$1.5\%$
189	$2^{94}-94$	$a-1$	$91.9\%$	$5.0\%$	$1.5\%$	$1.6\%$
198	$2^{99}-34$	1	$91.6\%$	$4.1\%$	$2.5\%$	$1.6\%$
199	$2^{99}-37$	$a-1$	$91.6\%$	$5.3\%$	$1.5\%$	$1.7\%$
208	$2^{104}-120$	1	$91.5\%$	$4.0\%$	$2.7\%$	$1.8\%$
209	$2^{104}-10$	$a-1$	$91.6\%$	$4.2\%$	$2.7\%$	$1.6\%$

presents numerical results on the number of required reduction steps for different field sizes. The four columns on the right in Table 2 provide the percentage for the number of required offset reduction steps after an arithmetic operation, where no reduction is required if the result q already satisfies $∥q∥\le N$. These results illustrate that sequential processing of the offset reduction is suitable for the point multiplication because $91\%$ of all operations require no reduction since $∥q∥\le N$ is satisfied. About $4-5\%$ of all operations require a single reduction step. Two or three calculation steps were required in approximately $4\%$ of all cases.

Algorithm 4 Point multiplication according to the randomized initial point method from [38].

input:P, $k={\left({\kappa }_{l-1},\dots ,{\kappa }_1,{\kappa }_0\right)}_{\tau }$

output:$k·P$

1: $R=randompoint$                    // randomized initial point 2: for ($j=l-1$ down to 0) do 3:  $S_j={\kappa }_j·P-(\tau -1)·R$         // store all precomputed points in memory 4: end for 5: $Q=S_{l-1}+\tau ·R$;              // ADD ${\kappa }_{l-1}·P-(\tau -1)·R$ and $\tau ·R$ 6: for ($j=l-2$ down to 0) do 7:  $Q=\tau ·Q$               // use DBL and ADD operations since $\tau \in \mathbb{Z}\left[\mathrm{i}\right]$ 8:  $Q=Q+S_j$     // ADD Q and ${\kappa }_j·P-(\tau -1)·R$ that is read from the memory 9: end for 10: return$Q-R$

On the other hand, the different execution times of the reduction steps may cause concerns about side-channel vulnerabilities. However, the probability for additional reduction steps is relativity low. Moreover, the randomized initial point method also randomizes the occurrence of these additional steps.

Next, we demonstrate that determining the point multiplication using Gaussian integers reduces the computational complexity. The number of required point operations per iteration of the point multiplication is similar to the results in [22]. However, some additional computations are required to subtract the random point. As in [22], we estimate the computational complexity for the point multiplication in terms of multiplication equivalent operations M.

Table 3. Complexity results for Algorithm 4 with different $\tau$-adic expansions in comparison with ordinary integer expansions from [38] for a binary key of length $r=163$.

Reference	${\left|\mathit{\tau }\right|}^2$ or $2^{\mathit{w}}$	l	M per PM incl. Precomputations
proposed	5	$0.430r$	2372
[38]	4	$0.506r$	3054
proposed	17	$0.245r$	1866
[38]	16	$0.2515r$	2763

illustrates two examples for the complexity with different bases $\tau$. The parameter r denotes the maximum key length in bits and the value l is the number of iterations per point multiplication. A larger value of ${\left|\tau \right|}^2$ reduces the number of iterations and speeds up the computation compared with a conventional binary PM, but requires more storage for the precomputed points.

The results are compared with the point multiplication according to randomized initial point but using ordinary integer expansions as proposed in [38] for comparable base values w. For instance, we can compare the complex base $\tau =4+\mathrm{i}$ with ${\left|\tau \right|}^2=17$ digits and the base $w=4$ with $2^w=16$ digits from [38]. For these values, the number M of multiplications is reduced by $32.5\%$. This reduction results mainly from the simpler calculation of $\tau ·Q$ for Gaussian integers due to the used endomorphism. In [21], a similar PM with $\tau$-adic expansions over Gaussian integers was proposed. This algorithm additionally reduces the number of precomputed and stored points. However, this method considers only timing and simple power analysis attacks. It is more vulnerable to statistical attacks than the randomized initial point method.

To illustrate the efficiency of Gaussian integer arithmetic, we consider latencies for the point multiplication according to Algorithm 4 using both norms in

Table 4. Field-programmable gate array (FPGA) synthesis results for the point multiplication according to Algorithm 4 using the absolute value (abs.) and the Manhattan weight (Man).

$\mathit{\tau }$	r	LUT	FF	RAM [kbit]	Slices	DSP	${\mathit{f}}_{\mathbf{clk}}$ [$\mathbf{MHz}$]	Protected PM Latency [$\mathbf{ms}$
								abs.	Man
$2+\mathrm{i}$	189	1540	521	17.3	426	4	227	7.59	6.60
$4+\mathrm{i}$	189	1540	521	19.1	426	4	227	5.68	4.93
$4+\mathrm{i}$	253	1891	678	20.4	532	4	212	7.92	6.87

. The table provides results for a Xilinx Virtex-7 field-programmable gate array (FPGA) based on the architecture from [22]. The hardware requirements are represented by the number of look-up tables (LUT), flip-flops (FF), slices, and digital signal processor (DSP) units, as well as by the RAM size. The maximum clock frequency is denoted by $f_{clk}$. These point multiplications are significantly faster than the results from [39,40]. For instance, the design in [40] considers key lengths up to $r=256$ and the unprotected PM. This design is optimized for the use of DSP units which reduces the number of LUT. The number 1990 of LUT is similar to the proposed design with key length $r=253$ using DSP units. It employs much more flip-flops (1786) and memory (234 kbit). An unprotected PM of length $r=256$ requires 23.5ms whereas the proposed protected PM requires only 6.87ms for $r=253$.

Note that the latency values for the reduction algorithms depend on the hardware architecture. The calculation of the Manhattan weight requires one addition which is performed in a single clock cycle with the architecture from [22]. For the absolute value, two additional multiplications are required which need four clock cycles each. Hence, the latency for calculating the Manhattan weight is only 11% of the time for the absolute value. The results in Table 4 illustrate the impact of this complexity reduction on the total latency of a point multiplication. Applying the Manhattan weight reduces the latency of a PM by 13% compared with the reduction algorithm from [29].

6. Conclusions

Gaussian integers are suitable for RSA [17,18] and ECC applications [22]. Furthermore, Gaussian integers have applications in coding theory [23,26]. The Montgomery multiplication for Gaussian integers was previously proposed in [17,18]. However, no reduction algorithm was advised.

The generalization of Montgomery reduction to Gaussian integers is not trivial, because complex numbers cannot be totally ordered. In this work, we have presented two algorithms for Montgomery reduction for Gaussian integers which use different norms. The first algorithm utilizes the absolute value to measure the size of a Gaussian integer. This algorithm can uniquely determine the correct representative. The second algorithm is based on the Manhattan weight of Gaussian integers, which reduces the computational complexity for the modulo reduction. However, two congruent solutions may occur. It is not possible to determine the correct candidate based on the Manhattan weight. The correct representative is the candidate with the smaller absolute value. Hence, the first algorithm is required to determine the final result, e.g., for the inverse mapping from the Montgomery domain to the original field representation.

It is shown in [21,22] that a key representation with a Gaussian integer expansion is beneficial to reduce the computational complexity and the memory requirements of elliptic curve point multiplication algorithms that are protected against side-channel attacks. However, only timing and simple power analysis attacks were considered. This PM is vulnerable to statistical attacks [38]. As an alternative, we have shown that Gaussian integer expansions can be used with the randomized initial point method. This method is protected against statistical attacks like DPA, RPA, and ZPA. The FPGA implementation of this protected PM with Gaussian integer expansions is significantly faster than the PM algorithms over ordinary prime fields reported in [39,40].

However, Gaussian integer fields can only be constructed for primes of the form $pmod4=1$, hence a generalization of this work to Eisenstein integers could be beneficial. Eisenstein integers are complex numbers of the form $x=a+b\omega$, where $\omega =\frac{1}{2}·\left(1+\sqrt{-3}\right)$, and Eisenstein integer fields can be constructed for primes of the form $pmod6=1$. An elliptic curve point multiplication using Eisenstein integers was considered in [49] showing similar properties as the point multiplication over Gaussian integers. Up to now no efficient modulo operation for Eisenstein integers is known. We believe that a generalization of the Montgomery modular multiplication to Eisenstein integers would be a promising direction for further research.