DESIGN OF HIGH-PERFORMANCE AND SCALABLE MONTGOMERY MODULAR

MULTIPLIER CIRCUITS

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of U.S. application Serial No. 63/014,390 filed April

23, 2020, the disclosure of which is hereby incorporated in its entirety by reference herein.

TECHNICAL FIELD

[0002] The present invention relates to a cryptographic systems, and more particularly, to a

Montgomery modular multiplier including both fixed-precision and scalable designs.

BACKGROUND

[0003] High-speed and area- efficient hardware implementation of modular multiplication

(MM) has been a subject of great interest. The design of a modular multiplication significantly affects the performance of many cryptosystems like Rivest-Shamir-Adleman (RSA) [23], elliptic-curve cryptography (ECC) [11], somewhat homomorphic encryption (SHE) [17], and fully homomorphic encryption (FHE) [6, 22]. Based on the hardness of the underlying mathematical problem, a cryptosystem usually deals with large integers from hundreds of bits to thousands of bits. More concretely, the performance of a variable bitwidth N of the modulus is an important concern. Therefore, efficient hardware realization of an MM algorithm (dealing with large integers) is a key challenge.

[0004] a, Basic Modular Multiplication

[0005] To better explain algorithms, some notation is provided. An N- bit integer X will be represented as in radix- 2 and in radix- 2 ^m where n _m = [N/m] and X _i = X[im + m — 1: im ] denotes bits of X belonging to word i. X is padded with (n _m * m — N) zeros to the left of the most significant bit (MSB). We refer to n _m as the word count and m as the word width.

[0006] Given a multiplicand X, a multiplier Y, and a modulus M, the modular multiplication

(MM) targets to compute Z = XYmodM . Realization 1 presents the basic N -bit modular multiplication, where we decompose the multiplier Y into n _m = [N/m] words. In the algorithm, we scan multiplier Y from the most significant word to the least significant word. In particular, we left shift the intermediate result Z ^(i+r) from the (i+l)-st iteration by m bits and add it to the product of X and Y _i . Next, the quotient is computed as t he floor of divided by M, and Z ⁽ⁱ⁾ is calculated by subtracting q ⁽ⁱ⁾ M from T ⁽ⁱ⁾ .

[0007] It is easy to show that Z ⁽ⁱ⁾ lies in the range [0, M) :

[0008] We can unfold the expression of Z® until Z ⁽ⁿ m ⁾ as shown below.

[0009] Although this realization is simple and immediately correct (i.e., it does not need a correction step as found in more advanced modular multiplication algorithms), it relies on doing large- bitwidth divisions to compute q ⁽ⁱ⁾ which is quite expensive. This realization is also referred as interleaved modular multiplication (IMM), especially when m is 1 or 2 and multiplication can be easily optimized.

[0010] b. Montgomery Modular Multiplication

[0011] Montgomery modular multiplication [14] is widely used because it replaces time- consuming divisions of the basic modular multiplication with multiplications and shifts.

[0012] Let M be an N-bit odd positive integer. Instead of directly computing Z = XYmodM, the operation can be performed efficiently in Montgomery domain with Montgomery modular multiplication. Assume R is a power of 2 and R > M. We can transform X (in the ordinary domain) into X (in the Montgomery domain) as below: where R ¹ is the modular multiplicative inverse of R with R ¹ ∈ (0, M) and RR ¹ modM = 1. The Montgomery form modM is represented by X, Y, and R ^-1 as follows: [0013] Notice that the transformation between X and can also be fitted into operation

XYR ^-1 modM with different X and Y, as shown below:

[0014]

[0015] Therefore, the design of Z = XYR ^-1 modM determines the performance of

Montgomery MM.

[0016] Depending on whether the bitwidth N is fixed or variable, the design of the

Montgomery MM falls into two categories: fixed-precision [3, 7, 12] and scalable [10, 25, 26] Montgomery MM. Given a multiplicand X , a multiplier Y , and a modulus M , the modular multiplication computes Z = XYmodM . A fixed-precision Montgomery MM with radix 2 ^m (FPR2tm) assumes the bitwidth N of the modulus is fixed and computes intermediate results related to m bits of the multiplier iteratively. When N is known, some optimization techniques, like Karatsuba multiplication and Toom-Cook multiplication [7], are possible to simplify the three required multiplications in each iteration. When N is variable, a scalable Montgomery MM with radix 2 ^m (MWR2tm) focuses on designing some ingenious modules, which process w bits of the multiplicand and the modulus, and m bits of the multiplier in each iteration.

[0017] When m = 1, the multiplication with 1 bit of the multiplier (radix 2 design) is simply a parallel collection of 2-input AND logic. Many prior art references focused on improving the radix 2 design, either for fixed-precision (FPR2tl) or scalable (MWR2tl) versions. The disadvantage of radix 2 design is its high cycle latency because it has to compute the intermediate results one multiplier bit a time. When m > 1, the cycle latency can be reduced to approximately 1/m. However, the design becomes quite costly in terms of area and delay of the resulting circuit, which hinders the development of a high-radix Montgomery MM.

[0018] c. Fixed-precision Montgomery Modular Multiplication [0019] Realization 2 shows the fixed-precision radix-2 ^m Montgomery modular multiplication

(FPR2tm), where we scan a group of m-bits of multiplier Y _i in the i-th iteration. Let d = [N/m] denote the number of iterations. Bazout’s identity implies that rr ^_1 — MM' = 1 when r and M are positive and relatively prime with r ^-1 ∈ (0, M) and M' ∈ (0, r). Both r ^-1 and M' can be computed through the extended Euclidean algorithm.

[0020] In line 3 of realization 2, the product X Y _i is added to the intermediate result obtained in the (i — l)-st iteration to produce the current intermediate result Z ⁽ⁱ⁾ . We compute <7 ⁽ⁱ⁾ in line 4 so that Z ⁽ⁱ⁾ + q ⁽ⁱ⁾ M becomes a multiple of r, and therefore, the computation of (Z ⁽ⁱ⁾ + <7 ⁽ⁱ⁾ M) /r in line 5 is reduced to a simple right shift (Z ⁽ⁱ⁾ + <7 ⁽ⁱ⁾ M) by m bits. To achieve this goal, we set q ⁽ⁱ⁾ = Z ⁽ⁱ⁾ M' so that where º _r denotes the congruence relationship modulo r. Notice that is just one of many possible choices for Adding/subtracting any multiples of r to/ from will still meet the requirement In other words, is unique up to a modulo on r. Consequently, valid choices of form a set as shown below:

[0021] Line 5 of realization 2 can be rewritten as in (8) below, which establishes a relation between Using this relation, we can recursively unfold

[0022] Consider the last iteration

Defining we have:

[0023] The Bezout’s identity implies that RR ^_1 — MM ^-1 = 1 since /? is even and M is odd.

Therefore, . [0024] The range 'modr can be estimated in equation (11). Since

M is a N- bit odd integer, we must have Consequently, the range of Z^ ^d~ ^ is [0,2 M), and only one correction in line 8-10 is needed to adjust Z back to the range [0, M).

[0025] In realization 2, we have to compute first and + q^M next. The critical path in each iteration thus encompasses three multiplications and two additions. When m = 1 , the multiplication with 1 bit of multiplier XY[i] can be simplified as a parallel collection of 2-input AND logic, as shown in realization 3. The quotient q^ is simply the least significant bit (LSB) of Z ⁽ⁱ⁾ in line 4. The critical path reduces to only two additions, which can be even optimized if we use carry- save compression. [0026] The disadvantage of FPR2tl design is its high cycle latency because it has to compute the intermediate results about one multiplier bit a time. When m > 1, the cycle latency can be reduced to approximately 1/m since we scan m bits of multiplier in one iteration. However, two challenges hinders the development of a high-radix Montgomery MM:

[0027] · Computationally expensive operations. In each iteration, the multiplications involved with XY _i and q^M are costly but necessary. The critical path grows up very fast when m increases. Although the cycle latency can be reduced to approximately 1/m, the time latency (product of cycle latency and clock period) does not reduce too much as we expect.

[0028] · Computation dependency. We have to compute quotient first and update intermediate result Z ⁽ⁱ⁾ next. This data dependency unavoidably cause a long critical path, and the situation becomes worse when the m increases.

[0029] d. Scalable Montgomery Modular Multiplication

[0030] References [8, 27] proved that the correction step is not immediately necessary.

Following equation (11) when r = 2 (m = 1), when inputs X and Y are in [0, M) we need d = N iterations to ensure Z in the range [0,2 M). When inputs X and Y are in [0,2 M), we just need to set R = 2 ⁿ⁺² and use d = N + 2 iterations in FPR2tl. When we need to do multiple modular multiplications, like when doing modular exponentiation, this setup allows us to start the next modular multiplication without doing the correction step. Therefore, we present various algorithms for scalable Montgomery MM to perform Z = XYR ^~1 modM where inputs X and Y and output Z are in [0,2M).

[0031] Based on the FPR2tl algorithm, reference [26] presented a radix- 2 scalable

Montgomery modular multiplication. Originally, the authors described their algorithms as a multiple word radix- 2 Montgomery multiplication (MWR2MM). To use a consistent naming convention throughout our article, we refer to the scalable Montgomery MM with radix- 2 as MWR2tl. Realization 4 shows the MWR2tl algorithm. Each outer loop, which goes through N + 2 iterations, scans 1 bit of the multiplier Y . The i-th outer loop iteration is equivalent to computing the result of the i-th iteration of the FPR2tl. MWR2tl uses an inner loop, which iterates e times, to scan a word of w > 2 bits of the multiplicand X in each iteration, e is equal to \(N + 2) /w] + 1 in binary form (and \(N + 1 )/w] + 1 in carry-save form) when X, Y E [0,2 M).

[0032] We compute based the least significant bit (LSB) in the 0-th inner loop iteration and use it in the remaining inner loop iterations. The w-bit additions in lines 5 and 9 result in a (w + 1) -bit integer. We assign the most significant bit to C _a or C _b and the remaining w least significant bits to Z _j . Finally, because we had to right shift Z by 1 bit in line 5 of the FPR2tl algorithm, the new z†} _t in line 10 of the MWR2tl realization is calculated as the concatenation of zj^ [0] and

Z _j \ ^L [w — 1: 1], Note that realization 4 also work for the case where w = 1. In that case, the only change is to modify line 10 of the 4 realization to read as z†} _t = z†^ [0] . [0033] Once we finish the 1-st inner loop iteration of the i-th outer loop iteration (i.e., i,j =

1), we have computed Z ₀ , which is required for the 0-th inner loop iteration of the (i + l)-st outer loop iteration (i.e., i + 1 ,j = 0). Therefore, if we use N + 2 processing elements (PEs) to do the outer loop iterations (one PE per outer loop iteration), we can start the next PE two cycles after the current PE. Exploiting this observation, we can design multiple-PE hardware that will have an issue latency (L) of 2. The clock period of a PE is lower bounded by two w-bit additions in binary form (lines 5 and

9).

[0034] Reference [26] presented a design of a radix 2 scalable modular multiplication

(MWR2tl) with the issue latency of L = 2. Later references such as [9, 10, 25] improved MWR2tl by reducing the issue latency to one (L = 1) by using elaborate methods. In spite of this low issue latency, the problem remains that, with a radix 2 design, we can only issue one bit of the multiplier every L cycles. As a result, we need at least LN cycles to finish the N- bit modular multiplication. This deficiency cannot be overcome no matter how many PEs are instantiated.

[0035] The scalable architecture of Montgomery MM provide us the capability when the bitwidth N of modulus M is variable. A PE only needs e cycles to complete the assigned outer loop iteration, and we can reuse it every e cycles, etc. As shown in realization 4, the value of e is determined by bitwidth N of modulus and word width w of multiplicand X. Ideally, we prefer at least p = [e /L] PEs so that the 1-st PE becomes free when the p-th PE generates Z ₀ . We can thus reuse the 1-st PE to compute the (p + l)-st outer loop iteration.

[0036] In contrast, we can only instantiate p PEs for a presumed bitwidth N of modulus M.

When the actual bitwidth is larger than the presumed, we do not have enough PEs. The 1-st PE is still busy when the p-th PE generates Z ₀ , and we need a memory management unit to buffer intermediate results until 1-st PE becomes free. A memory management unit could be designed as an internal first- in first-out (FIFO) queue or an external RAM module, whose memory size can be very large. Consequently, the actual bitwidth can be larger than the presumed, as long as we can buffer all intermediate results. The non-decompos ability challenge is thus solved. [0037] The problem in the fixed-precision Montgomery MM also exists in the scalable

Montgomery MM. Many references focus on improving radix-2 or radix-4 scalable Montgomery MM by reducing issue latency or critical path. High-radix design with fewer iteration is not well studied, since multiplication is an expensive operation and elongates the critical path.

SUMMARY

[0038] Disclosed herein are three designs and implementations of highly efficient

Montgomery modular multiplier circuits as summarized below.

[0039] A first design is that of a scalable high-radix (i.e., 2 ^m ) Montgomery Modular (MM)

Multiplication circuit replacing the integer multiplications in each iteration of the Montgomery MM algorithm (related to the product of m bits of the multiplier and the multiplicand) with carry-save compressions completely eliminating costly multiplications. The disclosed Montgomery MM decomposes the multiplicand itself using a radix of 2 ^W with w > 2m, thereby achieving a scalable design, which can deliver an issue latency of one cycle.

[0040] A second design is that of a high-performance iterative Montgomery modular multiplication circuit of fixed-precision, wherein the computations of the quotient and intermediate result in each iteration are done in parallel. This parallelism breaks the data dependency in each iteration, and thus, greatly reduces the overall circuit latency for performing a modular multiplication. This realization also replaces required multiplications and additions in each iteration with novel encoding and compressions.

[0041] The third design is that of a high-performance, yet scalable, Montgomery Modular

Multiplication circuit wherein the potential overflow of intermediate results is simply corrected by a 2-bit addition, resulting in a homogeneous decomposition of intermediate results, which in turn allows us the utilize a custom processing element to decomposes the multiplicand itself using a radix of 2 ^W with w > 3m and deliver an issue latency of one clock cycle. The design of the processing element can be pipelined, wherein the optimal number of processing elements can be reduced or alternatively the circuit throughput can be increased. [0042] The foregoing summary is illustrative only and is not intended to be in any way limiting. In addition to the illustrative aspects, embodiments, and features described above, further aspects, embodiments, and features will become apparent by reference to the drawings and the following detailed description.

BRIEF DESCRIPTION OF THE DRAWINGS

[0043] For a further understanding of the nature, objects, and advantages of the present disclosure, reference should be had to the following detailed description, read in conjunction with the following drawings, wherein like reference numerals denote like elements and wherein:

[0044] FIGURES la, lb, and lc. Embodiments of Montgomery modular multiplier circuits.

[0045] FIGURE Id, le, and If. Decoders that is used in the Montgomery modular multiplier circuits.

[0046] FIGURE 2. Example of z-th outer loop iteration when N=8, w=8, and m=2.

[0047] FIGURE 3. Data Dependency Graph of the MWR2tm Algorithm.

[0048] FIGURES 4a and 4b. Schedule for N = 6, w = 4, and m = 2 with (a) p = 3 and (b) p = 2.

[0049] FIGURE 5. Internal Design of the Processing Element (PE) for the MWR2tm realization with L = 1.

[0050] FIGURE 6. Proposed Radix 2 ^m Scalable Architecture of Montgomery Modular

Multiplication.

[0051] FIGURE 7. Special Compression for 5-to 2 in PE.

[0052] FIGURE 8. PEc Diagram for C Task.

[0053] FIGURE 9. Special Compression for 8-to-2 in PEc. [0054] FIGURES 10a and 10b. Time Latency of the Proposed Montgomery MM Design

Presuming N = 1,024 with Actual Values of (a) N < 1,024 and (b) N > 1,024.

[0055] FIGURE 11. Throughput of the Proposed Montgomery MM Design Presuming N =

1,024 and with Actual Values of N ³ 1,024.

[0056] FIGURE 12. Example of Encode when m = 8.

[0057] FIGURE 13. Model to group k = 2 bits.

[0058] FIGURE 14. Example of EncodeSN when m = 8.

[0059] FIGURE 15. Computing q ⁽ⁱ⁾ .

[0060] FIGURE 16. 8-bit signed multiplication with radix-4 modified Booth Coding.

[0061] FIGURE 17. Data model to accurately represent Z by Z _s , Z _c , and Z _carry .

[0062] FIGURE 18. Compute Z ⁽ⁱ⁾ .

[0063] FIGURE 19. Block-level diagram for the hardware implementation of the proposed

Montgomery modular multiplication.

[0064] FIGURE 20. Data Dependency Graph.

[0065] FIGURES 21a and 21b. Schedule for N = 8, w = 6, and m = 2 with (a) p = 3 and (b) p = 2

[0066] FIGURE 22. Internal Design of the Processing Element (PE) for the MWR2tm realization with L = 1.

[0067] FIGURE 23. Proposed Radix 2 ^m Scalable Architecture of Montgomery Modular

Multiplication [0068] FIGURE 24. Overflow when computing by truncation when m = 4.

[0069] FIGURE 25. Data model for overflow correction.

[0070] FIGURE 26. PEc Diagram for C task.

[0071] FIGURE 27. Table 1: Area and Performance of the Proposed Montgomery Modular

Multiplication Design for w and m, when the Modulus M Has N = 1024 Bits

[0072] FIGURE 28. Table 2: Comparisons of Different Montgomery Modular Multiplication with N = 1024 and N = 2048.

[0073] FIGURE 29. Table 3: Area and Performance of Different Montgomery Modular

Multiplication Designs on Virtex-7 FPGA.

DETAILED DESCRIPTION

[0074] Reference will now be made in detail to presently preferred embodiments and methods of the present invention, which constitute the best modes of practicing the invention presently known to the inventors. The Figures are not necessarily to scale. However, it is to be understood that the disclosed embodiments are merely exemplary of the invention that may be embodied in various and alternative forms. Therefore, specific details disclosed herein are not to be interpreted as limiting, but merely as a representative basis for any aspect of the invention and/or as a representative basis for teaching one skilled in the art to variously employ the present invention.

[0075] It is also to be understood that this invention is not limited to the specific embodiments and methods described below, as specific components and/or conditions may, of course, vary. Furthermore, the terminology used herein is used only for the purpose of describing particular embodiments of the present invention and is not intended to be limiting in any way. [0076] It must also be noted that, as used in the specification and the appended claims, the singular form "a," "an," and "the" comprise plural referents unless the context clearly indicates otherwise. For example, reference to a component in the singular is intended to comprise a plurality of components.

[0077] The term “comprising” is synonymous with “including,” “having,” “containing,” or

“characterized by.” These terms are inclusive and open-ended and do not exclude additional, unrecited elements or method steps.

[0078] The phrase “consisting of’ excludes any element, step, or ingredient not specified in the claim. When this phrase appears in a clause of the body of a claim, rather than immediately following the preamble, it limits only the element set forth in that clause; other elements are not excluded from the claim as a whole.

[0079] The phrase “consisting essentially of’ limits the scope of a claim to the specified materials or steps, plus those that do not materially affect the basic and novel characteristic(s) of the claimed subject matter.

[0080] With respect to the terms “comprising,” “consisting of,” and “consisting essentially of,” where one of these three terms is used herein, the presently disclosed and claimed subject matter can include the use of either of the other two terms.

[0081] It should also be appreciated that integer ranges explicitly include all intervening integers. For example, the integer range 1-10 explicitly includes 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10. Similarly, the range 1 to 100 includes 1, 2, 3, 4. . . . 97, 98, 99, 100. Similarly, when any range is called for, intervening numbers that are increments of the difference between the upper limit and the lower limit divided by 10 can be taken as alternative upper or lower limits. For example, if the range is 1.1. to 2.1 the following numbers 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, and 2.0 can be selected as lower or upper limits. [0082] When referring to a numerical quantity, in a refinement, the term “less than” includes a lower non-included limit that is 5 percent of the number indicated after “less than.” A lower non- includes limit means that the numerical quantity being described is greater than the value indicated as a lower non-included limited. For example, “less than 20” includes a lower non-included limit of 1 in a refinement. Therefore, this refinement of “less than 20” includes a range between 1 and 20. In another refinement, the term “less than” includes a lower non-included limit that is, in increasing order of preference, 20 percent, 10 percent, 5 percent, 1 percent, or 0 percent of the number indicated after “less than.”

[0083] The term “electronic component” refers is any physical entity in an electronic device or system used to affect electron states, electron flow, or the electric fields associated with the electrons. Examples of electronic components include, but are not limited to, capacitors, inductors, resistors, thyristors, diodes, transistors, etc. Electronic components can be passive or active.

[0084] The term “electronic device” or “system” refers to a physical entity formed from one or more electronic components to perform a predetermined function on an electrical signal.

[0085] It should be appreciated that in any figures for electronic devices, a series of electronic components connected by lines (e.g., wires) indicates that such electronic components are in electrical communication with each other. Moreover, when lines directed connect one electronic component to another, these electronic components can be connected to each other as defined above.

[0086] The processes, methods, or algorithms disclosed herein can be deliverable to/implemented by a processing device, controller, or computer, which can include any existing programmable electronic control unit or dedicated electronic control unit. Similarly, the processes, methods, or algorithms can be stored as data and instructions executable by a controller or computer in many forms including, but not limited to, information permanently stored on non-writable storage media such as ROM devices and information alterably stored on writeable storage media such as floppy disks, magnetic tapes, CDs, RAM devices, and other magnetic and optical media. The processes, methods, or algorithms can also be implemented in an executable software object. Alternatively, the processes, methods, or algorithms can be embodied in whole or in part using suitable hardware components, such as Application Specific Integrated Circuits (ASICs), Field-Programmable Gate Arrays (FPGAs), state machines, controllers or other hardware components or devices, or a combination of hardware, software and firmware components.

[0087] Throughout this application, where publications are referenced, the disclosures of these publications in their entireties are hereby incorporated by reference into this application to more fully describe the state of the art to which this invention pertains.

[0088] Abbreviation:

[0089] “FPGAs” means Field-Programmable Gate Arrays.

[0090] “PE” means processing element.

[0091] The terms “processing element” and “processing component” refers to a computer processor (e.g., CPU), programable logic array, Field-Programmable Gate Array, or any assembly of logic circuits designed to implement the predetermined steps set forth below.

[0092] Referring to Figures la, lb, and lc, embodiments of a Montgomery modular multiplier circuit are provided. The Montgomery modular multiplier 10 includes at least one digital processing component 12, which produces a modular multiplication result by (i) decomposing the multiplier into a group of words, each word containing a number of bits, and (ii) calculating one or more intermediate results based on the one or more intermediate results from a previous iteration and a product of the multiplicand and each word of the multiplier in an outer loop iteration. Montgomery modular multiplier circuit 10 can also include includes a plurality of processing elements (PEs) 14 ¹ and postprocessing processing element 16. The processing elements (PEs) 14 ¹ include addition block 18 and compressor bock 22 and post-processing processing element 16 includes addition block 18, where i is an integer label for the processing elements. Typically, i runs from 1 to the total number of processing elements. [0093] In a variation, the one or more intermediate results in each outer loop iteration is decomposed into a number of parts where some of the said parts are moved to a next iteration for processing. In a refinement, carry-save compressors 24 within compressor block 22 are used to replace multiplication operations to generate one or more intermediate results. In a further refinement, at least one digital processing component 12 is configured such that (i) the multiplicand is first left-shifted by some number of bits and is then decomposed into a second group of words, each word containing a second number of bits, and (ii) the one or more intermediate results in each outer loop iteration are calculated based on the one or more intermediate results from two previous outer loop iterations and products of each word of the multiplicand and one word of the multiplier in an inner iteration loop.

[0094] In a variation, the at least one digital processing component 12 is configured such that modular multiplication result is produced by (i) decomposing the multiplier into a group of words, each word containing a number of bits, and (ii) independently calculating two sets of one or more intermediate results based on the two sets of one or more intermediate results from a previous iteration and a product of the multiplicand and each word of the multiplier in an outer loop iteration. In a refinement, encoding circuit 20 are used to speed up the calculation of any expression that may be modeled as the product of a sum of two integers with another integer.

[0095] In another variation, at least one digital processing component 12 is configured such that (i) the one or more intermediate results in each outer loop iteration are partitioned into two sets which can be calculated independently, (ii) first encoding circuits are used to speed up the calculation of any expression that may be modeled as the product of a sum of two integers with another integer, and (iii) second encoding circuits are used to decompose each of the intermediate results in the set of one or more intermediate results into two corresponding intermediate results.

[0096] In another variation, the at least one digital processing component 12 is configured such that (i) the multiplicand is first left shifted by some number of bits and is then decomposed into a third group of words, each word containing a third number of bits, (ii) the one or more intermediate results in each outer loop iteration are calculated based on the two sets of one or more intermediate results from two previous outer loop iterations and products of each word of the multiplicand and one word of the multiplier in an inner iteration loop, and (iii) second encoding circuits are used to decompose each of the intermediate results in the set of one or more intermediate results into two corresponding intermediate results.

[0097] Referring to Figure la, Montgomery modular multiplier circuit calculates an output value Z congruent to XT/? ^_1 modM, wherein X is an (N + l)-bit wide multiplicand and Y is an (N + l)-bit wide multiplier with X and Y being integer numbers in the range [0,2 M), M is an N-bit wide modulus, R is an integer that is a power of 2 such that R > M, and R ^-1 is a modular multiplicative inverse of R with respect to M, the Montgomery modular multiplier circuit using a first input r ^-1 e [0, M) and a second input M' E [0, r)such that rr ^_1 — MM' = 1 with r = 2 ^m . The Montgomery modular multiplier circuit 10 includes at least one digital processing component configured 12 to: a) decompose multiplier Y into n _m = \(N + m + 2 )/m] words of m bits each, pad (n _m * m — N — 1) zero’s to the left of Y such that produce output value Z by iteratively calculating and combining XY _t r.

[0098] Still referring to Figure la, the at least one digital processing component 12 includes a plurality of processing elements (PEs) 14 ¹ such that the i-th PE produces first intermediate results based on one or more intermediate results from the (i-2)-nd PE, one or more second intermediate results from the (i-l)-st PE, and the value of XY _i V. The at least one digital processing component 12 also includes a post-processing processing element 16 which takes as input the intermediate results of the last two PEs and produces the output value Z. A weighted sum of first intermediate results produced by the i-th PE is congruent to X r ^~l .

[0099] In a variation, the at least one digital processing component 12 is further configured to: a) left shift multiplicand X by m bits to produce X' = Xr, pad (e * w — N — m — 1) zero’s to the left of X' , and decompose X' into e = [(iV + 2 m + 2 )/w] words of w bits each such that b) produce output value Z by iteratively calculating and combining Xj Y _t . [0100] In a refinment thereof, the at least one digital processing component 12 comprises: a) a plurality of processing elements (PEs) such that the i-th PE produces first intermediate results based on one or more intermediate results from the (i-2)-nd PE, one or more intermediate results from the (i-l)-st PE, and the value of XY _L r in e cycles, where in each cycle the i- th PE receives w bits of each intermediate result from the (i-2)-nd PE and the (i-l)-st PE, and produces w bits of the first intermediate results; and b) a post-processing PE which takes intermediate results of the last two PEs and produces the output value Z in e cycles.

[0101] In a refinement, the data initiation latency is 1 cycle and wherein the Montgomery modular multiplier circuit is scalable.

[0102] In a variation, two or more first intermediate results and two of the first intermediate results are decomposed into e = [(iV + 2m + 2 )/w] words of w bits, where m > 2 and w > 2m. In a refinment, at least one of the first intermediate results has zero’s at bit positions iw to iw + w- m- 1 for 0 < i < e — 1. In another refinment, at least one of the first intermediate results has zero’s at bit positions iw + w- m to iw + w- 1 for 0 < i £ e — 1.

[0103] In another variation, integer multiplication operations are unrolled into partial products, addition operations that are required to sum partial products and other results are replaced with k-to-2 compressors where k indicates numbers of integers that need to be compressed successively. In a refinement, the at least one digital processing component is further configured to pad (n _m * m — (N + 1)) zero’s to the left of multiplier Y and decompose the padded Y into n _m = [(iV + 2m + 2 )/m] words of m bits each, such that Y = E _j 2 ^im with Z _j = Y[im + m — 1: m] for 0 < i £ n _m — 1.

The at least one digital processing component includes: a) a plurality of processing elements (PEs), the i-th PE producing first intermediate results based on intermediate results from the (i-2)-nd PE and the (i-l)-st PE and such that a weighted sum of intermediate results is congruent to Y r ^-i ; and b) a post-processing PE taking intermediate results of the last and next to the last PEs and producing the output value Z. In a refinment thereof, there are two or more first intermediate results and at least two of the first intermediate results are decomposed into e = \(N + 2m + 2 )/w] words of w bits, where m > 2 and w > 2m. At least one of the first intermediate results can have zero’s at bit positions iw to iw +w-m-l for 0 < i £ e — 1. Moreover, the at least one of the first intermediate results can have zero’s at bit positions iw+w-m to iw+w-1 for 0 < i £ e — 1. In a refinment, one of the intermediate results is a 1-bit value. In yet another refinement, a Booth coding technique is used to speed up the compressors.

[0104] In a refinement, multiplicand X is left shifted by m bits to produce X' = Xr, padded with (e * w — m — (iV + 1)) zero’s on the left, and the padded X' is decomposed into e words of w bits each with such that The PEs compute their first intermediate results by repeatedly calculating the product of Xj Y _t for different i and j indices, thereby, enabling design of PEs that are independent of value of N, and producing a scalable Montgomery modular multiplier circuit. In a further refinement, the data initiation latency is 1 cycle. In a further refinement, the Montgomery modular multiplier circuit is scalable. Characteristically, the post-processing processing elements can repeatedly receive intermediate results of the last two PEs and produce the partial output value Z.

[0105J With reference to Figure lb, a Montgomery modular multiplier circuit is provided. The

Montgomery modular multiplier circuit 10 calculates an output value Z congruent to XT/? ^_1 modM, wherein X is an N- bit multiplicand and Y is an N- bit multiplier with X and Y being integer numbers in the range [0, M), M is a N- bit wide modulus, R is an integer that is a power of 2 such that R > M, and R ^-1 is a modular multiplicative inverse of R with respect to M. Characteristically, the circuit uses a first additional input r ^-1 ∈ [0, M) and a second additional input M" E [0, r ² ) such that r ² r ^-1 — MM" = 1 with r = 2 ^m . Montgomery modular multiplier circuit 10 includes at least one digital processing component 12 configured to: a) decompose multiplier Y into n _m = \N /m\ words of m bits each, pad (n _m * m — N) zero’s to the left of Y such that b) produce output value Z by iteratively calculating and combining values of XY{r ² .

[0106] Still referring to Figure lb, the Montgomery modular multiplier circuit 10 includes a digital processing element 14 ¹ configured to produce two intermediate results independently of one another in the i-th iteration, wherein the two intermediate results are calculated based on the value of XY _t r ² and the corresponding two intermediate results produced in the (i-l)-st iteration. The Montgomery modular multiplier circuit also includes a post-processing processing element 16 which takes the intermediate results of the last iteration as input and produces the output value Z. In a refinement, one intermediate result produced in the z-th iteration is congruent to is the multiplicative inverse of r ^i_1 modM.

[0107] Referring to Figure lc, in a the Montgomery modular multiplier circuit, the digital processing element(s) 14 ¹ produces two groups of intermediate results independently of one another in the z-th iteration. Moreover, the two groups of intermediate results produced in z-th iteration are calculated based on the value of XY{r ² and the corresponding two groups of intermediate results produced in the (i-l)-st iteration wherein each group contains three intermediate results and F _j = F _j — F _j [m — 1 ]r + Y _i -ilm — 1] , Encoder circuts 20 and 28 accelerate the calculation of intermediate results in the i-th iteration. Finally, a post-processing processing element which takes the intermediate results of the last iteration as input and produces the output value Z.

[0108] In a variation, a weighted sum of first intermediate results produced in the i-th iteration is congruent to 2f is the multiplicative inverse of r ^i-1 modM.

[0109] In another variation of the Montgomery modular multiplier circuit of figure lc, X is an

(iV + l)-bit multiplicand and Y is an (iV + l)-bit multiplier with X and Y being integer numbers in the range [0,2 M) wherein the number of iterations d is at least \(N + 4) /m] and R is r ^d~2 ; and postprocessing processing element 16 takes the intermediate results of the last iteration as input and produces the output value Z. [0110] In another variation, the at least one digital processing component 12 includes a plurality of processing elements 14 ¹ such that the i-th PE produces first intermediate results based on one or more intermediate results from the (i-2)-nd PE, one or more second intermediate results from the (i-l)-st PE, and the value of XY _i r ² . The at least one digital processing component 12 also includes a post-processing element which takes as input the intermediate results of the last two PEs and produces the output value Z. In a refinement, there are two or more first intermediate results and two of the first intermediate results are decomposed into e = \(N + 2m + 3 )/w] words where w > 3m.

[0111] In a variation of Figure lc, the Montgomery modular multiplier circuit is configured such that at least one of the first intermediate results has zero’s at bit positions iw to iw + w — m — 1 for 0 < i < e — 1 and at least one of the first intermediate results has a tail bit at bit positions iw for 0 < i < e — 1.

[0112] In another variation of Figure lc, the Montgomery modular multiplier circuit is configured such that at least one of the first intermediate results has zero’s at bit positions iw + w — m to iw + w — 1 for 0 < i £ e — 1 and the weight of bit positions iw for 0 < i £ e — 1 is — 2 ^lw .

[0113] In another variation of Figure lc, the Montgomery modular multiplier circuit is configured such that multiplicand X is left shifted by 2m bits to produce X' = Xr ² , padded with (e * w — N — m — 1)) zero’s on the left, and the padded X' is decomposed into e words of w bits each with such that such that the PEs compute their first group of intermediate results by repeatedly calculating the product of XjY _i for different i and j indices, thereby, enabling design of PEs that are independent of value of N, and producing a scalable Montgomery modular multiplier circuit. In a refinement, the data initiation latency is 1 cycle and the Montgomery modular multiplier circuit is scalable. In another refinement, the post-processing processing element 16 repeatedly receives intermediate results of the last two PEs and producing the partial output value Z. In a refinement, the processing elements are designed into a pipeline with L stages where the data initiation latency is L cycle. [0114] With reference to Figure Id, an encoder circuit is provided. The endocoder 30 calculates (t + l)-bit outputs d and b from t-bit inputs a and b, wherein t ³ 3. The encoder circuit ensuring that v _& + vg = a + b where and vg = . The encoder includes at least one digital processing component 32 configured to: a) calculate a 3-bit intermediate result d as d[ 2: 0] = a[t — 1: t — 2] + b[t — 1: t — 2] + (a[t — 3]ANDd[t — 3]); b) if t ³ 4, calculate d[t — 4: 0] = a[t — 4: 0] and b[t — 4: 0] = b[t — 4: 0]; c) calculate d[t — 3] = a[t — 3]XORh[t — 3]; d[t — 2] = d[0]; d[t — 1] = rf [1] ; d[t] = d[ 2]; and d) calculate b[t — 1: t — 3] = 0; b[t\ = d[l],

[0115] In a refinement, v^ _[t —1 _: 0 _] + ^v B _[t -1 _: 0 _] congruent to a + hmod2 ^f with v^ _[t —1 _:0] ⁼

[0116] In a variation, encoder 30 calculates 3-bit outputs d and b from 2-bit inputs a and b, wherein the encoder circuit ensures that + vg = a + b with v _& = 4d[2] — 2d[l] + d[0] and vg = 4b [2] — 2 b [1] + b [0] . The encoder includes at least one digital processing component 32 configured to calculate a 3-bit intermediate result d as d[2: 0] = a + b; calculate d = d; and calculate b\l \ 0] = 0; b[ 2] = d[l], In a variation, i7 _a[i: o _] + ^v B _[i-. o _] congruent to a + dmod4 with i7 _a[i: o _] ^{= —} 2d[l] +

[0117] With reference to Figure le, an encoder circuit is provided. The encoder 34 calculates

1-bit outputs c _out and e and a 2-bit output d from 2-bit inputs a and b, and 1-bit input c _in . The encoder circuit ensures that 4 c _out + 4e — 2d[l] + d[0] = a + b + c _in . The encoder 34 includes at least one digital processing component 36 configured to alculate a 2-bit value as 2c [1] + d[0] = a[0] + d[0] + c _in ; calculate c _out = a[l]ANDd[l] ; calculate d[l] = a[l]XORd[l]XORc[l] ; and calculate e = a[l]XORd[l]ORc[l]. [0118] With reference to Figure If, encoder circuit 38 calculates (t + l)-bit outputs a and b from t-bit inputs a and b and 1-bit input c, where t ³ 2. The encoder circuit ensuring that + ^v b = b[2i + l]2 ^2l+1 . The encoder includes one digital processing component 40 configured to: a) encode a\ 1: 0], b[ 1: 0], and c to generate outputs c) encode ^ to generate outputs d) calculate e) calculate b[2i + 2] = for 0 < i < t/2 — 1. In a refinement, if parameter t is odd, inputs a and b are padded by one 0 on the left. In another refinement, i½ _[ί-ΐ:q] + ^v b _[t -i _·. o _] is congruent to

[0119] In a refinment of the variation set ofrth above,, s equal to

[0120] In a refinment of the variation set ofrth above, is equal to

[0121] Additional details of the invention are set forth below.

[0122] 1. Design of a Scalable Montgomery Modular Multiplier Circuit

[0123] 1.1 Fixed-precision Optimization

[0124] Depending on the value of m, the high-radix Montgomery MM falls into one of two categories: (i) using explicit and expensive multiplication units when m is large, and (ii) replacing the explicit multiplications with simplified logic expressions when m is small. When m is large, e.g., m = 16, prior art references have mainly focused on optimizing the MM time latency (product of cycle latency and clock period) and hardware cost (chip area) when they have two multiplication units [13, 24] or more than two multiplication units [16]. Considering area and delay, the bitwidth of the multiplication units cannot be made very large and are commonly chosen as 16 bits to maintain a good area-latency trade-off. When m is small, e.g., m = 2 , the multiplication unit can be explicitly simplified, for example by using modified Booth coding to reduce partial products. It would be desirable to develop a parametric high-radix (fixed-precision or scalable) Montgomery MM to support arbitrary m.

[0125] 1,1,1 Operation Rescheduling

[0126] The value of Z is updated twice in each iteration of realization 2, which in turn necessitates separate hardware resources and increases the data dependence. Fortunately, the FPR2tm realization can be rescheduled to address this shortcoming. More precisely, the accumulation of XY _i+1 into Z ^(t+1) in the (i+l)-st iteration can be moved to the i-th iteration. Consequently, we can compute q^M and XY _L r in parallel and accumulate them into Z^K

[0127] Algorithm 5 is a modified version of the FPR2tm realization incorporating the rescheduling technique, where d is the number of iterations.

[0128] The value of d influences the upper bound of output Z . Indeed, we must find the minimum value of d to guarantee that the output Z lies in [0,2 M) as explained next. Based on the equation in line 4 of the algorithm 5.1.1, we may explicitly compute Z ⁽ⁱ⁾ in terms of X, M, r, Y _t , ···, Y _i , by recursively unfolding the Z ⁽ⁱ⁾ until Z* ^-1) = 0 in (12). Note that q^ is 0.

Pn+P

[0129] In each iteration, we can scan m bits of the multiplier Y . When i ³ m , the multiplier

Y is completely scanned, and we can find a tight upper bound of Z ⁽ⁱ⁾ as follows:

[0130] Therefore, the upper bound of Z ^(d ^ in the last iteration is

(14)

Following the equation 14, 4 Mr so that we can guarantee that In other words, if d ³ the output Z will always lie in [0,2 M) as required.

[0131] 1,1,2 Intermediate Result Decomposition

[0132] During the computation of q^ ^l \ we only need the m least significant bits of Zl ^i-1 l because of the operation Zl ^i-1 lmodr. When w > 2m, we can decompose Z into ZM and ZR as shown in (15) below. Since w — m ³ m, we only need ZR (more precisely ZR ₀ ) to compute q.

[0133] This decomposition seems of little value for the FPR2tm algorithm. However, it is an important observation in order to support the MWR2tm algorithm. As shown in reference [25], the MWR2tl realization achieves L = 1 when w > 2 and m = 1. This article generalizes the aforesaid result and proves that the MWR2tm realization can also achieve L = 1 when w > 2m for m > 1. We describe how the word decomposition is incorporated into the FPR2tm algorithm. In the 0-th iteration of realization 5, we can decompose Z* ^-1 ! into ZM* ^-1 ! and ZR^ ^_1 \ Since the w — m least significant bits of ZMi ^-1 ! are 0, we have and therefore we only need to use ZR^ ¹ ^

[0136] We observe that ZM* ^-1 ! /r is not necessary to the required computations in the 0-th iteration and can be moved to the next iteration. We define /r and decompose may thus be rewritten as:

[0137]

[0138] As before is not necessary to compute q^\ i.e., [0139]

[0140] Again, is not necessary to the required computations in the 1-st iteration and can be moved to the second iteration. Since = 0, we can define = 0 and Z R _jiew ⁼ ZR^ ^~r> = 0 . Now then, for an arbitrary iteration i , we define ^ ( ^ _W decompose the new \

[0141] Realization 6 is a modified version of the FPR2tm realization incorporating rescheduling and decomposition techniques. Note that the decomposition is quite easy to do in hardware.

[0142] 1.1.3 Elimination of Multiplication Operations

[0143] Even after rescheduling the original realization to remove some data dependencies, the performance of the FPR2tm realization continues to be limited by multiplications. Implementation of a multiplication operation involves generating partial products and summing them together. The timing critical paths of this operation are the carry chains for the addition of partial products. We can use a ternary tree arrangement of CSAs to speed up the compression. Meanwhile, if we use the carry- save form to pass partial results from one iteration to next and only deploy the addition at the end of the algorithm, the performance will improve further.

[0144] Instead of using Z ⁽ⁱ⁾ , we can use the carry-save form ZS ⁽ⁱ⁾ and ZC ⁽ⁱ⁾ in line 4 of realization 6 and eliminate the final addition as detailed below.

[0145]

[0146] ^{1 J}

[0148] where ZS ⁽ⁱ⁾ » m = |ZS ⁽ⁱ⁾ /rJ, meaning that we right shift by m bits and thus ignore the least m significant bits of ZS ⁽ⁱ⁾ . As a result, Z ⁽ⁱ⁾ /r can be rewritten as:

[0149]

[0150] Based on the modular definition, we have

[0151]

[0152] Since Z ⁽ⁱ⁾ is a multiple of r, (ZS ⁽ⁱ⁾ rnodr + ZC ⁽ⁱ⁾ modr) can only be either 0 or r.

When ZC ⁽ⁱ⁾ rnodr = 0, we have to choose ZS ⁽ⁱ⁾ rnodr = 0 since ZS ⁽ⁱ⁾ rnodr cannot be r. When ZC ⁽ⁱ⁾ modr ¹ 0, we have to choose ZS ⁽ⁱ⁾ rnodr = r — ZC ⁽ⁱ⁾ modr ¹ 0 since ZC ⁽ⁱ⁾ modr cannot be negative. We may define the indicator c ⁽ⁱ⁾ to represent the (ZS ⁽ⁱ⁾ rnodr + ZC ⁽ⁱ⁾ rnodr) /r as follows (Notice that it is also correct to use ZS ⁽ⁱ⁾ rnodr to define c ⁽ⁱ⁾ ). An indicator l(x) = 1 when the condition x is true and 0 otherwise.

[0153] In summary, Z^/r in line 5 of realization 6 may be replaced by the summation of

ZS^ » m, ZC^ » m, and . Moreover, ZS^ » m and ZC^ » m can be decomposed into (ZSM^ ^l \ZSR^) and (ZCM ⁽ⁱ \ ZCR ⁽ⁱ ^). respectively, and moved forward to the next iteration.

[0154] Realization 7 is a modified version of the FPR2tm realization utilizing operation rescheduling, word decomposition, and CSA compression. Each iteration now consists of three compressions and one addition instead of three multiplications as was the case in realization 2. The first compression in line 5 compresses the 5-operand addition into 2 (temporary) integers (TS^ ^l \ Next, the required computation is performed into two steps: a compression and an addition. In particular, a second compression in line 6 takes O(logm) units of delay to compress the 2m-operand addition into 2 integers Notice that products each produce m partial product terms that need to be added together. The m least significant bits of qS^ and qC^ are summed up in line 7, where the sum is generated by a logarithmic carry propagate adder (CPA), like the Kogge-Stone adder (KSA), requiring 0(logm ) levels of logic. Finally, the third compression in line 8 compresses the (2m + 2)-operand addition into 2 integers Notice that multiplication with r = 2 ^m is simply done by left shifting the operand and that the products q^M and XY{r, each produce m partial product terms which must be summed together with TS^ and TC^K

[0155J 1.1.4 Booth coding

[0156] Because of the modular property to compute is not always identical to use qS ⁽ⁱ⁾ modr and qC^modr unless the sum of them is always less than r. Although it is possible to replace qZ ⁾ with the resulting (qS ^d> modr)M and (qC^modr)M would cause double area compared with q^M. Instead, we choose to complete the modular addition Consequently, we can use radix-4 modified Booth coding to reduce the number of partial products for each multiplication performed in line 6 of realization 5.1.3 from 2m to m + 2. Subsequently, we save m — 2 CSAs and also reduce the delay. Modified Booth coding is not compatible with the carry-save form because it can achieve the correct result only after the final addition and ignore potential overflows. This is why we do not use Booth coding to improve the third compression in line 8 of realization 7.

[0157] 1.2 Radix-2 ^m Scalable Montgomery Design

[0158] To extend the fixed-precision realization 7 to a scalable one (see realization 8 below), we need to iteratively compute w bits of intermediate results in realization 7 (see line 8). Moreover, we must find the maximum values of the intermediate results in order to derive e in terms of m. Equation (12) provides a way to compute the intermediate result Z ⁽ⁱ⁾ as shown next.

[0159]

[0160]

[0161] Equation (24) shows that the maximum possible value of intermediate result Z ⁽ⁱ⁾ is always less than 2 ^N+m+2 so that it can be presented in N + m + 2 bits. Note that lines 10 and 11 of the realization 7 performs a division by r using a right shift by m bits. The maximum possible value of the intermediate result is the value just before the division, which can be presented in N + 2m + 2 bits. Therefore, minimum value of e can be set as . Notice that that w > 2m, utilizing the triangle inequality 1, we can derive a somewhat pessimistic lower bound on e which is e = + 1 as was done for realization 4.

[0162] Partial results ZSR^ ¹ ^ and ZCR ^(l ^ in the i-th iteration of realization 7 may be explicitly represented as in equation (25), where the most m significant bits of ZSR _{j (i} —1) and ZCR _{j (i} -i) are 0. [0163]

[0164] Similarly, partial results ZSM^ ^l~2 ^ /r and ZCM^ ^l~2) /r may explicitly be represented as in equation (26), where again the least w — m significant bits of ZSM ²⁾ and ZCM _j ^{1 2)} are 0. Note that ZSM^ ^l~2 ^ and ZCM^ ^l~2 ^ are multiples of r = 2 ^m because w — m ³ m.

[0165]

[0166]

[0167]

[0168]

[0169]

[0170]

[0171] Finally can also be represented explicitly as follows:

[0172]

[0173] The expression + XY _t r in line 8 of realization 7 can thus be written as: [° ¹⁷⁴ 1

[0175] In the context of the MWR2tm realization where N is variable, when Y _i is provided in the i-th iteration, we need an inner loop to scan w bits of X' and M in successive iterations. Consider any iteration i of the outer loop. In the 0-th inner loop iteration (j = 0), let us redefine intermediate results ( TS,TC ) as the compression of ZSR _Q + ZCR _Q + ZSM ₀ '/r + ZCM ₀ '/r + c · 1 (j == 0), where

ZSR _Q and ZCR _Q denote ZSR^ and ZCR^ in (current) iteration i whereas ZSM _Q and ZCM _Q ' denote ZSM _Q ¹ ' ⁾ and ZCM^ ¹ ' ⁾ in (previous) iteration i — 1. We can guarantee that only TS and TC need w — m + 1 bits as shown in figure 7. From equation (28), TS and TC contain necessary information to allow calculation of q in the 0-th inner loop iteration.

[0176] Next, we compress TS + TC + X _Q YI + qM ₀ into intermediate results (OS, OC) .

Assume the bitwidths of OS and OC are w + x where x > 0. The w least significant bits of OS and OC are actually ZS^ and ZC^ (parts of ZS^ and ZC^ in line 8 of realization 7), respectively. Consequently, we can compute c^ ^l \ ZSR^ \ and ZCR^ \ The most x significant bits of OS and OC will affect the computation of ZSR^ and ZCR^ and so on. These are named as FBS and FBC and used in the next inner loop iteration. Similarly, when we get to the 1-st inner loop iteration (j = 1), we can update ( TS, TC ) as the compression of ZSR ₁ + ZCR ₁ + ZSM[/r + ZCM[/r + c · 1 (j == 0). Furthermore, we update (OS, OC) as the compression of TS + TC + X[Y _i + qM _t + FBS + FBC. As before, the w least significant bits of OS and OC are actually ZS^ and ZC^ such that we can partition and assign them to ZSM ₀ , ZCM ₀ , ZSR _i, and ZCR ₁ . FBS and FBC are still updated as the x most significant bits of OS and OC. The procedure can proceed until all e inner loop iterations are completed. [0177] One important question is how to set the bitwidths of FBS and FBC to guarantee that the bitwidths of OS and OC are large enough to avoid overflows. In another word, assume FBS and FBC each have x bits, the sum of OS and OC must be precisely representable with w + x bits. Equation (31) captures the constraints.

[0178]

[0179]

[0180]

[0181] Based on realization 7, we describe a radix- 2 ^m scalable Montgomery modular multiplication (MWR2tm) with an issue latency of L = 1. Because the final result Z is in the range [0,2 M), Z can also be represented in e words.

[0182] The last thing is to collect all carry-save form results together and compute the final

Z = XYR ^~1 modM. Once we have computed , we are able to compress their sum into two (w + 1) -bit integers (PS, PC) . The w least significant bits PS[w — 1: 0] and PC[w — 1: 0] cooperate with c (working as the carry-in bit) to compute the last word Z ₀ . The carry-out bit during the addition automatically works as carry-in in the next cycle. The two most significant bits are named as DO 1 and DO 2 to be added with ZSM[/r, and ZCM[/r. We thus compute Z _x and update DO 1 and DO 2. The procedure continues until we collect all e words of Z. We can guarantee it suffices to use two 1-bit DO 1 and DO 2 in (34). Detail of the compression is in the figure 9.

[0183] Figure 2 shows an example of the i-th outer loop iteration when N = 8, w = 8, m =

= 2. The modulus M has 8 bits so that M ₁ are 0. Because X E [0,2 M), X' = Xr has iV + l + m = 11 bits and the least m = 2 significant bits are 0.

[0184] [0185] 1,3, Scheduling and Performance Estimation

[0186] Figure 2 shows the data dependency graph of realization 8. Node A denotes tasks of 0- th inner loop iteration for the MWR2tm algorithm. Node B denotes tasks of inner loop iteration when j ¹ 0. Node C denotes the task in lines 25-31 of the MWR2tm algorithm. We can start the A task of PE i as soon as we finish the A task of PE i — 1. The issue latency is thus L = 1.

[0187] Each PE needs e cycles to execute one A task and (e — 1) B tasks. Ideally, we need PEs. In the context of a scalable design, the bitwidth N of modulus M is variable such that we cannot know how many PEs are required in advance. As we mentioned earlier, after e cycles, a PE becomes free and may thus be reused. Assume p PEs are employed. Let k denote the number of rounds we use all p PEs, resulting in kp outer loop iterations in algorithm 5.2. From (14), it is easy to see that k = guarantees that output Z will lie in [0,2 M) as required.

[0188] Figure 4 shows an example when N = 6, w = 4, m = 2, and e = 3. When we have p = 3 PEs (p > e) in figure 4(a), we need k = 2 rounds. PE 1 finishes the computation about Y ₀ and becomes free when we need to issue Y ₃ . Therefore, we need kp cycles to finish all A tasks. When we only have 2 PEs (p < e) in figure 4(b), we need k = 3 rounds. When we want to issue Y ₂ , PE 1 is still busy and we have to wait. Therefore, we need (k — l)e + p cycles to finish all A tasks. When all A tasks are done, we can start the C tasks and collect all carry-save form results, which needs e + 1 cycles.

[0189] Based on the magnitude of p and e, the number of cycles required to perform Z =

XYR ^~1 modM (where inputs X and Y and output Z are in [0,2 M)) by using the MWR2tm realization is as follows: where

A simpler (upper bound) equation for the cycle count is: cycle_cnt < /cmax(e,p) + min(e, p) + 1 (37) where

[0190] Note that we use radix-2 ^w for decomposing the multiplicand into e words of w bits each, and radix r = 2 ^m for decomposing the multiplier.

[0191] To start a required iteration of realization 7 on a new PE, we must collect enough information to compute the least m significant bits of TS^ and TC^ in line 5 so as to compute q^k When w > 2m, after the first cycle, we finish the computation about the w — m least significant bits because we need to right shift by m bits. After the second cycle, we finish the computation about the 2 w — m ³ m least significant bits. The procedure goes on until all e inner loop iterations are done. So the issue latency is L = 1.

[0192] Unfortunately, when m £ w < 2m so that w — m < m , the MWR2tm realization cannot support an issue latency of L = 1. Instead, we have 2 w — m ³ m and thus can only start computations on a new PE with an issue latency of 2 cycles. Generally speaking, the issue latency of realization 8 is L = [(2m) /w]. One may modify realization 8 and the corresponding data dependency graph to support various combination of w and m, some resulting in an issue latency of 1 cycle, others to an issue latency of 2 cycles or more.

[0193] Other optimization methods (such as increasing the cycle latency or increasing the circuit area as in [9, 10]) may be used to achieve L = 1 even when m < w < 2m. In view of the above discussion, setting w > 2m provides the highest performance realization (i.e., realization 7) in terms of circuit area (hardware usage), cycle latency, and issue latency (L = 1). [0194] 1.4. Hardware Implementation

[0195] As we can observe, tasks A and B have many similarities and thus can be implemented in one PE by using a control signal CT as depicted in figure 5. In the A task ( CT = 1), the 5-to-2 compression handles variable c, the (2m+4)-to-2 compression masks inputs FBS and FBC to be 0 and forces the least m bits of outputs OS and OC to be 0, and finally q is calculated. In the B task (CT = 0), the 5-to-2 compression masks c and thus forces it to be 0, the (2m+4)-to-2 compression accepts inputs FBS and FBC and generates OS and OC, and finally the stored value of q is used. Notice that ZSMi, ZSRi, etc. are represented by m bits, w m bits, etc. instead of w bits, since the remaining bits are known to be zeros.

[0196] Figure 7 shows the 5-to-2 compression of c · 1 (j == 0), ZSM-/r, ZCM[/r, ZSR^, and

ZCRi by using two (w — m)-bit CSAs. Recall that the w — 2m least significant bits of ZSM[/r are 0.

[0197] Unlike A and B tasks, the C task needs only one special processing element (named

PEc). Figure 8 shows the implementation of the PEc. When CT = 1, we force feedback signals DO 1 and DO 2 to 0 and set carry signal as c. When CT = 0, we accept feedback signals DO 1 and DO 2 and set carry signal as the carry-out bit from the previous cycle.

[0198] Although there are 8 integers, we can compress them by using two w-bit CSAs as in figure 9. DZSRi stands for the delayed ZSR _L (it is delayed by one clock cycle). The main delay of the C task is thus a w-bit addition, where we can use a KSA with a gate delay of log(w) to speed up the computation.

[0199] We show the block-level architecture of the proposed radix-2 ^m scalable Montgomery modular multiplication in figure 6. The required ZSM[ and ZCM[ in the i-th outer loop iteration in realization 8 are actually ZSM _t and ZCM _t which are generated in the (i — 2)-nd outer loop iteration. Similarly, ZSM[ and ZCM[ required in the PEc are ZSM _t and ZCM _t generated by the next to the last PE. To compute Z _j , we need ZSM[, ZCM[, ZSM _t , ZCMi, ZSR^ ZCRi, and c. However, they are not generated in the same cycle. More precisely, ZSM _t and ZCM _t are generated in the (i + l)-st inner loop iteration whereas ZSM[ , ZCM[, ZSR^ ZCR^ and c are generated in the i-th inner loop iteration and need to be delayed by one clock cycle as shown in figure 8. When p < e as in figure 4(b), a first- in first-out queue with a depth of e — p slots is necessary to buffer the outputs.

[0200] 1.5. Hardware Emulation Results

[0201] 1.5.1 Analysis of Area and Time

[0202] The MWR2tm realization was coded in the Verilog hardware description language and implemented using the Xilinx Vivado Virtex 7 xc7v585tffgl 157-3. Given a chip area constraint and the word size m of multiplier, there are many possible kernel configurations with the number p of PEs and the word size w of multiplicand for a designer- specified bitwidth N of the modulus M.

[0203] The number of cycles is a function of e and p as shown in equation (35). When e < p, the cycle count is kp + e + 1. A PE becomes free before we issue another word of the multiplier, meaning that we have instantiated too many PEs. The reduction of p can improve resource utilization until p = e. Note also that the change of p does not affect cycle count since kp is roughly a constant. When e > p, the cycle count is ke + p + 1. A PE is still busy when we want to issue another word, so we have to wait for PEs to become available. Increasing p causes k to be reduced and thus lowers the cycle count until p becomes equal to e. Therefore, for a designer- specified bitwidth N of the modulus M, the optimum value of p is e.

[0204] Table 1 presents our implementations for different w and m combinations, (see, Figure

N +2771+2

27) As we discussed, the optimum p is e = , which is almost independent of m when w >

W

2m, and thus, decreases with the increase of w. Meanwhile, increasing w means that a PE is more complex, and hence, it will require more LUTs to do the required compressions and addition. Consequently, different implementations of the Montgomery MM with a fixed m have nearly the same area (in terms of the number of LUTs they need). The clock period of the kernel in figure 6. is bounded by the critical paths in both the PEs and the PEc. By increasing w, more computations are implemented as purely combinational logic, and thus the area decreases. However, when w becomes too large, the latency of the KSA addition in the PEc dominates the clock period such that the Montgomery MM time latency increases whereas the cycle latency is somewhat reduced. The placement and routing inside a PE also become more complex and require more area. Based on the product of area and time, the configurations (32,2), (64,4), and (128,8) for (w, m) achieve the best trade-off. When we run multivariate linear regression, we can estimate the area (LUT + FF) to scale as follows:

LUT_count = 4.34 wmp + 6107 (39)

[0205] In the context of a scalable Montgomery MM design, N is a designer- specified variable, and therefore, we cannot optimize the architecture for a given N a priori. Instead, we presume an N value under constraints on the time latency and area. Choices of w, m, and p determine the latency and area for a presumed N as shown in table 1.

[0206] Figure 10 shows the time latency of three scalable Montgomery modular multiplications with a presumed value of N = 1,024. When the actual (designer-specified) N is less than or equal to 1,024, the difference of time latency is small as shown in figure 10(a). However, when the actual N is larger than 1,024 where e > p, a PE is still busy when we want to issue another word of the multiplier such that we have to wait. Figure 10(b) shows the curve of time latency when N > 1,024. When m increases, the slope decreases. Consequently, for actual N = 16,384, the time latency of the configuration w = 128, m = 8, and p = 9 is only 170.64 ps whereas the configuration w = 32 , m = 2 , and p = 33 results in 332.2 p s time latency. Figure 11 shows the corresponding throughput of the same configurations for the presumed N = 1,024.

[0207] 1,5,2 Performance Comparison

[0208] Table 2 reports the performance of different Montgomery modular multiplications in terms of time and area, (see, figure 28). Notice that different radix (different m values) result in different products of area and time. Ideally, it would be fair to compare our design with other scalable designs under the same m value. However, to the best of the authors’ knowledge, there is no high- radix scalable design, which does not use DSP’s. Strictly speaking, we do not know the complexity of a DSP (used as a multiplier) in terms of an equivalent number of FUTs. Although we cannot compare area with DSP-based design [1, 15], the worst latency of our design when m=2 is still smaller than their latency. Meanwhile, many references [13, 16, 20] (regardless of whether they are scalable or not) focus on optimizations for a given radix (i.e., m is fixed). When m is small, the required multiplications can be explicitly optimized. When m is large, the multiplication cost becomes non- negligible, and many designs make use of the well-optimized DSP modules on FPGA. It is difficult or inefficient to extend these designs to support a parameterized m value.

[0209] References [4, 5, 21] provide parameterized but non-scalable designs. Recall that a scalable design optimized for a presumed value of N , for example N = 1,024, can still support Montgomery modular multiplication for any N > 1024, e.g., N = 2048 and N = 4096, whereas a non-scalable design cannot support cases when N > 1024 . Evidently, there are time and area overheads to convert a non-scalable design to a scalable one (if this conversion is possible at all), for example to perform decomposition of the intermediate results and employ complex connections among PEs. When m = 2 and m = 4, the area and time product of our proposed scalable design is 29.95% and 17.44% larger than the reference [21]. However, our design’s overhead is constant since each PE has a fixed number of inputs and is only connected with at most two neighboring PEs. Recall that the critical delay of a PE is that of performing three compressions and one m-bit addition. Therefore, the clock period of a PE scales as O(logm). As a result, when m = 8, our scalable design has a smaller product of area and time than the non-scalable designs in [4, 21].

[0210] 4, Inventive Design of a High-Performance, Fixed-Precision Montgomery

Modular Multiplier Circuit

[0211] Feedforward Montgomery MM [19, 28] is a variant of Montgomery MM, which is able to compute in parallel. However, feedforward Montgomery MM can only guarantee that

Z ^(d_1) in the last iteration falls in the range We cannot easily adjust Z back to the range

[0, M) with a constant number of corrections. Note that the value of d can change for different variants of the Montgomery MM. In contrast, this thesis proposes an efficient fixed-precision Montgomery

MM. For the sake of clarity, we first present a basic version of proposed Montgomery MM, in which we parallelize the computation while guaranteeing that Z^ ^d~r) lies in the range [0,2 M). The basic version shows the correctness of proposed Montgomery MM. Next, we present an advanced version of the proposed Montgomery MM, in which multiplication and addition operations are replaced with encoding and compression operations. This advanced version significantly reduces the critical path delay while guaranteeing that Z^ ^d~r) lies in the range (— M, 2 M).

[0212] 2,1 Basic version,

[0213] Realization 9 shows the basic version of proposed Montgomery MM. In this algorithm, q9 ⁾ and Z ⁽ⁱ⁾ are independent of each other and can thus be computed in parallel. Based on BAOzout’s identity, we can find r ^-1 E (0, M) and M” E (0, r ² ) such that r ² r ^_1 — MM” = 1.

[0214] Let us prove the correctness of realization 9 by mathematical induction. At the beginning Assume and We would like to derive a q^∈ [0, r) such that Z ⁽ⁱ⁾ + q^M º _r 0. Note that this condition must be met in every iteration of the realization to ensure that the computational result is correct. To do this we first define: as in line 3 of realization 9. Next we show that

[0215] Furthermore, equation (40) implies the following:

[0216] Because both q4 ^{l 1)} and are valid choices in [0, r), due to uniqueness of up to a modulo on r, we conclude = g ⁽ⁱ⁾ modr. Therefore, g^ — is a multiple of r i.e., gii) _ _q a- 1) = _r o (42)

[0217] In line 4 of realization 9, we have Because is 2 m bits and is m bits, q^ is also m bits. With this choice for q^ ^l \ and based on line 5 of realization 9, we can now prove that + q^M º _r 0. From equation (40), we know And the proof is complete.

[0218] Consequently, the three multiplications in realization 9 may be done in parallel. The critical path for computing Z ⁽ⁱ⁾ is thus reduced to one multiplication and two additions. We can also unfold the expression of Z ⁽ⁱ⁾ until as shown below:

Because Z^ ¹⁾ = 0, we have q ^{( 1)} = q^ = 0.

[0219] When i > [N/m] — 1 , we have å _k ^l ₌₀ Y _k ^{rk =} Y · Because q^ E [0, r) , we can estimate the upper bound of Z ⁽ⁱ⁾ as follows:

[0220] When (i - l)m > JV , we have M < 2 ^N £ 2 ^(i_1)m = r ^i_1 such that Z ⁽ⁱ⁾ < 2 M .

Therefore, when we set d = [N/m] + 2, we have Z^ ^d_1 ^ E [0,2 M). Compared with realization 2, realization 9 only needs two additional iterations.

[0221] 4,2 Advanced version

[0222] Although realization 9 breaks the data dependency between q^ and Z ⁽ⁱ⁾ and thus reduces the computation latency, multiplications and additions still hinder the development of a very high-performance Montgomery modular multiplication. Realization 10 shows an advanced version of the proposed Montgomery MM, which replaces all multiplications and additions in each iteration with encoding and compression operations. The result is an realization which can be implemented in hardware with much lower hardware area and computation latency than any and all prior art designs.

[0223] In this algorithm, Z^ ^r) is represented by three values: an (N + 2m + l)-bit signed number zf \ an (iV + 2m + l)-bit signed number Z^ ¹ \ and a 1-bit unsigned value Z^ _ar ]. _y . Note that we do not explicitly perform the addition.

[0224] The operation (Z^ ^modr ² ) is simply replaced with ¹ ' ⁾ as follows:

₇ a-1)

Notice that the difference b is represented by an m-bit unsigned number q _{s (} i—1) , an m-bit unsigned number q _{c (} i—1) , and a 1-bit unsigned value q (i-i) . carry·

[0225] Consider computing the expression Z = (A + B) * C, which shows up in each iteration of the Montgomery MM. Normally, we would need to first perform the full-size addition D = A + B followed by the multiplication Z = D * C (during which we can utilize the modified Booth coding to cut the number of generated partial product in half). This serial sequence of add and multiply operations increases the latency per iteration of algorithms 3.3 and 6.1. Notice that the hardware resource cost to compute the expression Z = A * C + B * C directly is costly so that is not an option. Instead, we propose a novel encoding technique (called Encode), which removes the necessity of performing the full-bitwidth addition of A and B when computing (A + B) * C, it does so by grouping and encoding every k bits of A and B and using these encoded k- bit groups to generate partial products of D * C. Although the Encode function is good for multiplication, the results of Encode are presented in a special format, which does not correspond to either signed number or unsigned number representations. Therefore, we propose an auxiliary version of Encode (call EncodeSN for Encode in Signed Numbers), which results are standard signed numbers. It will be shown that the sum of any values encoded by the EncodeSN function is equal to the sum of these values when they are encoded by the Encode function so that EncodeSN is used to perform addition and subtraction whereas Encode is used to do D * C.

[0226] We also point out that the Encode function when k = 2 also results in removal of half of the partial products of D * C, similar to using a radix-4 modified Booth coding when doing D * C. In this thesis, we discuss and show the implementation of this encoding approach for the k = 2 case. Extension to higher values of k is straight-forward.

[0227] 2,2,1 EncodeSN

[0228] Using a similar procedure as the one we did for realization 9, we determine so that Z ^(l~1] + qh- ^!j /vf = _r o. Recall that q^ ^l~r) is only unique up to a modulo on r as proved in (7). This observation motivates our encoding techniques, name Encode and EncodeSN. EncodeSN calculates (t + l)-bit outputs d and b from t-bit inputs a and b. The encoding process consists of the following steps: a) calculate a 3-bit intermediate result d as d[ 2: 0] = a[t — 1: t — 2] + b[t — 1: t — 2] + (a[t — 3]ANDh[t — 3]); b) calculate d[t — 4: 0] = a[t — 4: 0] and b[t — 4: 0] = b[t — 4: 0]; c) calculate d[t — 3] = a[t — 3]XORh[t — 3] ; d[t — 2] = d[0] ; d[t — 1] = d[l] ; d[t] = d[ 2]; d) calculate b[t — 1: t — 3] = 0; b[t] = d[l].

The values of outputs d and b are defined as ensuring that v _& + vg = a + b.

[0229] We first perform EncodeSN(¾ (i—1) , (i— l) ). Although the outputs are m + 1 bits, we truncate outputs and assign the least m bits and q^ qu _arry ' ^s simply the qu _arry itself. Because q _s [m]2 ^m and q _c [m]2 ^m are multiples of r = 2 ^m , the difference between =

We use the notation = EncodeSN(qd ^i_1) ) to represent the said encoding process.

[0230] Figure 12 shows the case when m = 8, where the superscript (i — 1) is omitted for brevity. The signs + and — stand for positive and negative weights respectively. As an example, = 418 and = —94. The difference = 512 = 2 * 256 = 2 * 2 ⁸ = 2 r.

[0231] 2,2,2 Encode

[0232] To explain the relationship between EncodeSN and Encode, we must carefully analyze the various steps of a 2-bit encoding as shown in Figure 13. Considering 2-bit unsigned numbers a = 2a _! + a ₀ and b = 2 b _t + b ₀ , and 1-bit unsigned number c _in . The goal is to compute c _out , d, and e such that

CL + b + C _jjj — 4c _0Uf — 2di + dg + 4e (51)

Note that d _j = d[i] for i = 0,1.

[0233] We define 2 * + d ₀ = a ₀ + b ₀ + c _in for step (i) of figure 13. Table 3 presents the truth table for different and c _x so that a ₁ + b ₁ + c ₁ = 2 d ₂ + d _t + 2 c _out in step (ii). The overflow due to a _t = b _t is captured in c _out = a-,ANϋZ) _! . TABLE 3

Truth table for « _÷ , X , c _¾

[0234] In step (iii), we split d ₁ as 2 d ₁ — cO Notice that d ₂ and cannot be 1 at the same time, therefore, we only need e = d ₂ 0Rd ₁ to represent cl ₂ + cl ₁ in step (iv). In conclusion, realization 11 shows the 2-bit encoding, which takes three inputs a, b, and c _in and generates outputs e, d, and c _{out .} The critical path delay is equal to that for a 2-bit adder. To avoid confusions, XOR and AND have the same precedence, which are higher than OR.

[0235] Taking a second look a = EncodeSN(q^ O, we actually do 2-bit encoding for the leading 2 bits o with c _in = qf [m — SJAND and ignore the outputs c _out and e. Note that m — 3] and are adjusted accordingly. Instead of only encoding the leading 2 bits of and we can encode every 2 bits. [0236] Realization 12 calculates (t + l)-bit outputs a and b for two t-bit inputs a and b and a 1-bit input c. Without loss of generality, we assume the bitwidth t for inputs a and b is even. Otherwise, we can pad a 0 to the left of a and b. The output c _out in a 2-bit encoding is used as the input c _in for the next 2-bit encoding. For example, the c _out from Encode2bit(a[t — 3: t — 4], b[t — 3: t — 4], c _in ) serves as the c _in for Encode2bit(a[t — 1: t — 2], b[t — 1: t — 2], c _in ).

[0237] Recall that the sum of inputs is equal to the sum of outputs for a 2-bit encoding as shown in (51). The goal of the Encode procedure is to represent the sum a + b + c in a different format. The values of a and b are defined as ensuring that Vg _, + VE = CL + b + c.

[0238] If we also apply (a, b ) = EncodeSN(a, b ), we must have [0239] We can perform qu _arry )· Similar to the computation of we can truncate outputs and assign the least m bits to ¾ (i—1) and ¾ (i— l) . Because of equation (53), we have:

We use the notation = Encode(g^ ^i_1 ^) to represent the said encoding process.

[0240] Although realization 12 is presented by a for loop, all 2-bit encodings can be preformed in parallel. Figure 14 shows the example when m = 8. In step (ii), we group 2 bits and compute c _out for each group. In step (iii), we perform 2-bit encoding and generate outputs e and d for each group. c _out is already computed and used in step (ii). In step (iv), we combine results into two numbers. Output e when encoding the leading 2 bits is ignored. In step (v), we fill empty slots with 0. Notice that we get all information in step (iii). Steps (iv) and (v) are only shown to explain the computation to the reader. The critical path delay to perform Encode is roughly equal to that of a 2-bit adder.

[0241] An important observation is that we can group 2-bit parts of again. As shown in step (v) of figure 14, some bits are constant 0. Their value can be treated as (—2) * 4 ³ + 2 * 4 ² + l * 4 ¹ + (—2) * 4° = —94. The value range for a group is [—2,2]. Like radix-4 modified Booth coding, we can use a group to generate a partial product directly. When we compute we only have m/2 partial products. There is no need to do an m-bit addition for q^ ¹ ' ⁾ + q^ ¹ ' ⁾ + qu _arry ^an >' more. Since the range of a group is [—2,2], the range of m-bit — 1 ),^(r — 1)] as shown below: [0242] 2.2.3 Computation of q®

[0243] Recall the computation of = (Z^ ¹⁾ modr ² )M"modr ² in realization 9, resulting

[0244] If we add/subtract any multiples of r ² to/from (Z^ ^i_1 ^modr ² ), the above equation will still hold. In other words, it is safe to replace (Z^ ^i_1 ^modr ² ) with z[ ^l ^ or z[ ^l ^ = Encod e(Z^ ¹ ' ⁾ ), because the difference between any two of and z[ ^{l 1} ' ⁾ is a multiple of r ² . The benefit is that we will only have m partial products if we compute Z^ . We thus have:

[0245] Equation (57) implies:

[0246] Considering equation (50) and the uniqueness of up to a modulo on r, we have

[0247] Following equation (43) which was related to realization 9, the obtained q^ for realization 10 should meet the following requirement:

[0248] Figure 15 shows compression of — q ^(i_1) for the least significant 2m bits. In step (i), we use CSA to compress m + 3 numbers into two unsigned numbers Notice that each partial product has a tail bit. The aforesaid m + 3 numbers consist of three parts: m numbers

~(i — l-! (i — l-! (I — l-! from Z _L , 2 numbers from q _s ^J and q _c , and 1 number constructed by collecting the tail bits and qu _arry We then use a balanced ternary tree arrangement of CSAs to perform the compression

(Rp’ Rp) ⁼ Compress(Z^ ^M" — q6 ^-1) ) with O(logm) number of full adder cell delays. Here, we do the fixed -bitwidth compression and collect the least significant 2m bits of the two numbers after compression because any overflow is a multiple of r ² .

[0249] Because of equation (59), the sum of the least significant m bits can only be 0 or r.

When qp ⁽ [m — 1] = 1 or qp [m — 1] = 1, the following equality holds: qp[m — 1: 0] + qp [m — 1: 0] = r . On the other hand, when qP[m — 1] = qP[m — 1] = 0 , we have qP[m — 1: 0] + qp [m — 1: 0] = 0. Therefore, we may define qu _arry ^as follows:

[0250] In step (ii), the computation of is reduced to a right shift by m bits:

[0251] Note that the resulting q6 ⁾ = (qp, qp , qp _rry ) is valid since it meets equation (60).

Consequently, the computation of q6 ⁾ _re lies on encoding and compression operations, rather then the more expensive multiplication and addition operations.

[0252] 2,2,4 Error Analysis [0253] Before doing the error analysis for Z^l , let us briefly discuss the optimization to compute AY _j based on the modified Booth coding. Instead of computing AY _j , we scan Y _j and Y^_ ₁ [m — 1] and encode them as Y _j , that is,

Y _j i[m - 1] + im] (63) 1 + im])4 ^fc

Note that U_ _c = 0.

[0254] Following the radix-4 modified Booth coding, we can further divide Y _j into m/2 groups if we group every 2 bits. In case m is an odd number, we can sign-extend Y _j by 1 bit and generate [m/2] groups. Alternatively, we can directly apply radix-2 ^m modified Booth coding to Y _j . Because the range of — 2Y[2k + 1 + im ] + Y[2k + im ] + Y[2k — 1 + im ] is [—2,2], we can use it to generate a partial product. Consequently, the product of XY^ has m/2 partial products, whereas XY^ has m/2 + 1 partial products. Moreover, the range of Y _j falls in [— r/2, r/2] as shown next:

[0255] We can further compute the range for å^= ₀ ¾ ^as follows:

[0256] When ( i + l)m — 1 > N, we have Y[(t + l)m — 1] = 0 so that [0257] Since , we can unfold the expression until

Z* ^-1) = 0 as shown below:

Note that because Z^ ¹⁾ = 0, we have

[0258] When i > [(iV + 1 )/m] — 1, the following equality holds:

[0259] We can derive upper and lower bounds for

[0260] When meeting the condition M/r ^i_1 < 4/3, we have — M < 2 M. Equation (70) given below simplifies this condition. Because log y ¹ < N + 2, this condition is met when N < (i — l)m. Therefore, when we set d = [N/m] + 2, we have Z^ ^d~ ^ E (— M, 2M) . We only need one correction, either addition or subtraction, to adjust Z back to [0, M) . [0261] We can also compute the maximum size of any intermediate result as follows.

Recall that

[0262] In other words, independent of the iteration count i, any can be represented by an

(iV + 2m + l)-bit signed number.

[0263] 2,2,5 Computation of Z ⁽ⁱ⁾

[0264] The last thing to do is to efficiently compute At this time, it is not safe to add or subtract a multiple of r. Our proposed encoding and modified Booth coding save half of the partial products and reduce the hardware area. After compression, the partial products are compressed into two signed numbers. We can do an addition (while ignoring any overflows) to get the correct result. Figure 16 shows an example to compute 8-bit signed multiplication for C = AB with radix-4 modified Booth coding. For A = 108 and B = 125, the correct product C is 13,500. Once we get all 4 partial products, we compress them into two 16-bit signed numbers CS and CC. Although CS + CC = —52,036 ¹ 13,500, we can perform the 16-bit addition while ignoring the overflow bit, thereby obtaining the correct result.

[0265] The said addition of the two signed numbers seems necessary. We show next that in fact this costly add operation can be avoided and the computation of Z^ can be accomplished by using only encoding and compression operations.

[0266] As we discussed in equation (71), Z^ is an (iV + 2m + 1) -bit signed number. We want to accurately represent Z^ by a triplet: an (iV + 2m + 2)-bit signed number Z^K an (iV + 2m + 2) -bit signed number z£\ and a 1-bit unsigned value Z^ _rry . Figure 15 shows our general data model, where the leading 2 bits of Z _s are same and the leading 2 bits of Z _c are 0. Because the resulting Z is an (N + 2m + l)-bit signed number, Z[N + 2m + 1] must be the sign extension.

[0267] We define the carry from the addition for the least significant (N + 2m) bits as C _in .

Table 4 shows all possible combinations of 5 _! and C _in . Notice that it is impossible to have = 0 and C _in = 1 because the resulting Z becomes larger than 2 ^N+2m , which would violate the fact that Z is a (iV + 2m + l)-bit signed number. For the remaining three cases, can never have any overflows, therefore, the sum Z ⁽ⁱ⁾ can be correctly presented.

TABLE 4

Truth table for the proposed data model

[0268] We again prove the result by using mathematical induction. The base case is ^r) =

Z _V ¹ ' ⁾ + z£ ¹ - ⁾ + Z^H _y = 0. Assume Z^ ^i_1 ^ = (zf Z^ Z^ _ar l _y ) fits in the data model of figure 15, if we show that also fits in the data model, then we will have proven that the sum Z ⁽ⁱ⁾ can be correctly presented.

[0269] Figure 18 shows the computation of Z ⁽ⁱ⁾ . Both q^ ^L~r> M and XY{r ² generate m/2 partial products. Z ^< ' ^1~1) is represented by Z^ ¹ \ z£ ^{l 1} \ and Z^ _a ^l _r ^ _y . Because Z ⁽ⁱ⁾ is an (iV + 2m + 1) -bit signed number, is an ( [N + 3m + 1) -bit signed number, which is also a multiple of r. To sum all partial products together, the first thing to do is to sign extend the numbers so that all of them have same bitwidth. Here, we deliberately extend them into N + 3m + 2 bits in step (i) of Figure 18. For sake of brevity, all signs are labelled as S although different numbers could have different sign values. [0270] The sign bits in a number can be converted into l’s followed by a S in step (ii). Any overflows can be ignored because Z^V only has (iV + 3m + 1) bits. S is the Boolean inversion of S. We further merge all constant l’s into a number called prefix in step (iii) after ignoring any overflows. Because the prefix is a constant number which is only dependent on m, the real implementation can skip steps (i) and (ii), and start at step (iii). When m > 4, we can prove that the maximum value for dashed polygon in step (iii) is less than 2 ^N+3m+1 as shown below:

[0271] Since the dashed polygon is less than 2 ^N+3m+1 , the compression of the m + 3 numbers will result in two (iV + 3m + 1) -bit unsigned numbers and Z ⁽ⁱ⁾ in step (iv) regardless of the compression strategy. Note that Z _S [N + 3m] and Z _C [N + 3m] cannot be 1 at the same time; otherwise Z ⁽ⁱ⁾ + Z ⁽ⁱ⁾ will be larger than 2 ^N+3m+1 . As a result, the sum of Z _S [N + 3m], Z _C [N + 3m], and the leading 2 bits of the prefix will be a 2-bit number with the same value, which is Z _S [N +

3m]XNORZ _c [iV + 3m] . Since Z _S [N + 3m] is always the same as Z _S [N + 3m + 1] , we can hide Z _S [N + 3m + 1] and Z _C [N + 3m + 1] in a hardware implementation.

[0272] Because Z^ ^i_1 ^ + XY{r ² + q^ ^l~ ^M º _r 0, in step (v), the operation (Z^ ^i_1 ^ + XY _i r ² + can be simplified as shown below:

Because fits the data model in Figure 17, we thus correctly present without doing an addition. [0273] When m < 4, we can modify the data model of Figure 17 to have N + 2m + 3 bits, sign extend all numbers into (iV + 3m + 3) bits in step (i), and prove the dashed polygon is less than 2 ^N+3m +2 j _{n slC} p ( jjjj The remaining proofs are similar to what we did above and are thus omitted here for brevity.

[0274] 2,3 Hardware Implementation

[0275] Figure 19 shows the block-level diagram for the hardware implementation of the proposed Montgomery modular multiplier. With the help of Encode, EncodeSN, and compression, the critical paths to update (jqs _> qc _> S _carry ) ^an d ( Z _s ,Z _c ,Z _carry ) correspond to those for doing a 2-bit addition and compression of m + 3 numbers.

[0276] The critical path of the proposed Montgomery MM thus lies in the addition module for the final summation and correction step. Considering the trade-off between time and area, we use a high performance logarithmic Brent-Kung adder (BKA) [2] and specify a 3-cycle timing constraint. Note that we need to perform one final summation and do at most one correction. We only need to use BKA twice in 6 consecutive cycles, which justifies the 3-cycle timing constraint. The control signal CT is used to control the data flow of BKA. When CT = 0, we compute Z _s + Z _c + Z _carry and store it in Z. When CT = 1, we check the sign of Z. If Z < 0, we compute Z + M. If Z > 0, we compute Z — M. If the computation result is negative, we do not need any correction and simply choose the original Z. Otherwise, the result, which is positive, is stored into Z. In conclusion, the total clock cycle count to do the Montgomery MM is only

CycleCount = u + 6 (74) where u = [N/m] + 2 is the number of iterations.

[0277] 2,4 Complexity Analysis

[0278] Table 5 shows complexity analysis for different algorithms in terms of area and time.

The iteration count, which remains the same for all algorithms, is equal to 0(N /m). Assume we use a textbook multiplier to perform an N X m multiplication with area and time complexities of O(iVm) . For algorithms 3.3 and 6.1, the critical path in each iteration is bounded by the multiplications such that the area and time complexity per iteration of the realization are O(Nrri). Meanwhile, the clock period (latency per iteration) for realization 10 is reduced a lot due to replacement of additions and multiplications with compression and encoding operations. More precisely, since encoding only takes a constant delay and compression has O(logm) time complexity, the clock period is O(logm) time. Although encoding technique and modified Booth coding reduce the number of partial products by a factor of two, the area complexity remains O(Nrri). As a result, the area-time product complexity is reduced from 0(N ³ m) in algorithms 3.3 and 6.1 to 0(N ² \ogm) in realization 10.

TABLE 5

Complexity Analysis of Montgomery MM Algorithms realization 2 OiN/m) OiNrn) O(Nrn) 0{N' ^:> m ) realization 9 0(N/rn) OTTVm) OiNrnj OiN^ ^' ni} realization Ϊ0 0(N/rn) O ( log mj 0(7Vm) 0(W ² iogrn)

[0279] 2,5 Hardware Emulation Results

[0280] The proposed fixed-precision Montgomery MM design was coded in the Verilog hardware description language and implemented using the default flow of the Xilinx Vivado Virtex-7 xc7v585tffgl 157-3 device, including synthesis and implementation. Instead of optimizing the realization for a specific value of m, our design maintains the flexibility to implement different m values, reflecting the trade-off between latency and area.

[0281] References [21] and [4] present two Montgomery MM designs that are parametric on m. When m doubles, the area roughly doubles whereas the cycle count is cut in half. However, the clock period also increases a lot. The period when m = 8 is more than double that of m = 2.

Consequently, the the area-time product (ATP) increases a lot when m increases. [0282] Because of the modified Booth coding and the proposed Encode/EncodeSN, we eliminate half of the partial products. Meanwhile, the critical path of our design only depends on m so that the clock period does not change by much. The increase in the clock period is mainly due to the larger fan-out count and more challenging placement and routing problems when N changes from 1,024 to 2,048. Considering the configuration N = 1,024 and m = 8, compared with [21] and [4], our design reduces the computation latency by 49.53% and 46.09%, respectively, and the area by 31.97% and 36.33%, respectively. As a result, we save more than 65% in terms of the ATP metric. The advantages of our design in other configurations are also evident. When increasing m, the clock period gradually increases so that the savings in computation latency decrease accordingly. Based on the experimental results for our proposed Montgomery MM, we achieve the minimum value of the ATP metric for m = 8. Finally, reference [18] provide a Montgomery MM design specifically optimized for the configuration N = 2,048 and m = 4. Compared with its synthesis results, our flexible design (which is parameterizable in m) achieves 29.77% reduction in the ATP. (see, Table 6, Figure 29).

[0283] 3 Inventive Design of a High-Performance, Scalable Montgomery Modular

Multiplier Circuit

[0284] It is desirable to design a high-performance, yet scalable, architecture for doing

Montgomery modular multiplication. In the realization 8, both the quotient and intermediate result for each iteration are unsigned numbers so that it is not necessary to worry about the overflow problem mentioned in figure 16. Instead, both quotient and intermediate result of the realization 10 are signed numbers. This property gives rise to two problems:

• Once we apply intermediate result decomposition in section 5.1.2, the resulting ZM _e-t is a signed number whereas the lower significant words are unsigned numbers. We need a special task to complete the last iteration of the inner loop.

• Following the data model in figure 17, we can prevent any overflows of applying compression on the signed numbers by representing results in w + m + 2 bits. In the last iteration of the inner loop, we need to fit the results in w + m bits but it would not be safe to simply truncate the leading 2 bits.

[0285] 3,1 Variant for a Different Input Range

[0286] As we discussed, the scalable design targets to compute Z = XYR ^~1 modM when inputs X and Y and output Z all lie in the range [0,2 M) . This can be done easily by adjusting the number of iterations.

[0287] Intermediate result Z^ satisfies the following relation:

When (i + 1) > \(N + 2)/m], we have åj( ₌₀ Y _k ^ ^k = Y so that

[0288] Based on equation (76), we can calculate the range of Z^ as follows:

Therefore, when 4 M /r ^d~2 < 1/3 or 12M < r ^d~2 , we have∈ (— M, M) . We can set d =

\(N + 4)/m] + 2 to ensure 12M < 16 M < 2 ^N+A < 2 ^m(d_2) = r ^d_2 . Realization 13 shows the details. [0289] Since we prove Z^ ^d ^∈ (— M, M) , the final result Z = Z^ ^d ^ 4 M e (0,2 M) .

Equation (78) generalizes the constraint: md ³ N A 2 m A 4 (78)

When the number of iteration d meets the above equation, the final result Z lies in (0,2 M).

[0290] 3.2 Homogeneous Decomposition of the Intermediate Result [0291] Taking a closer look at realization 13, it can be seen that we only need the 2m least significant bits of zf ¹ ' ⁾ and Z^ ¹ ' ⁾ to compute q^ ^l k We may decompose both Z^ ¹ ' ⁾ and z£ ^{l 1} ' ⁾ into e words where each word has w bits (see the discussion of e below). In this case, the most significant word is a signed number whereas the remaining (lower significant) words are unsigned numbers. Instead, we extend the decomposition for Z (Z _s or Z _c in any iteration) into e words as shown below:

Notice that Z[i] = Z[N + 2m] for i ³ N + 2m to perform sign extension and Z[— 1] = 0.

[0292] In this case, Z is homogeneously decomposed into words Z£) for 0 < i < e — 1. Each

ZE _t is a w-bit signed number and a 1-bit unsigned tail number such that: Z = åfrd (ZE _j )2 ^wl . We can further decompose Z£) into two parts ZME _t and ZRE _t as shown in (80). The least w — m bits of ZMEi are 0 and the most m bits of ZRE^ are 0. When w ³ 3m and w — m ³ 2m, we only need ZRE (more precisely ZRE ₀ ) to compute q.

[0293] Realization 14 shows the Montgomery modular multiplication after the aforesaid homogeneous decomposition. We simply place Z^ _rry in the following bit of ZCRE^ , say

[0294] 3,3 Scalable Architecture [0295] Similar to realization 8, we use an inner loop in realization 15 to do the operation

Compress^XY _i r ² + qMBM + ZREMB + ZMEMB / _r ) by scanning w bits of X' = Xr ² , M , ZREMB _{^ an} d ZMEMB i _n each iteration.

[0296] 3.3.1 Inner Loop Computations [0297] In the first iteration of the inner loop (J = 0), we need to compute the sum of ZSRE ₀ ,

ZCRE ₀ , ZSME ₀ /r, ZCME ₀ /r, X^Y _i , and qM ₀ . Notice that they are signed numbers. We can use two (w + x)-bit signed numbers OS and OC (x > 0) to accurately represent the sum. Based on the least w bits, we can compute ZSRE ₀ and ZCRE ₀ as in (81). Notice that OS + OC is a multiple of r, we can embed Z _carry = OS[m — l]OROC[m — 1] to the tail bit of ZCRE ₀ . The leading x bits of OS and OC are stored in FBS and FBC for future reference.

[0298] In the second iteration (j = 1), we need to compute the sum of ZSRE ₁ , ZCRE ₁ , FBS, and FBC . Again, we represent them into two (w + x)-bit signed number OS and OC. Based on the least w bits and the homogeneous decomposition in (80), we can compute ZSRE^ZCRE _t and ZSME ₀ /ZCME ₀ . The leading x bits of OS and OC are again stored in FBS and FBC for future reference.

[0299] The next question is what bitwidth w + x for OS and OC is sufficient to accurately represent intermediate results in each iteration of the inner loop. Following the data model in figure 17 and optimization in figure 18, we can prove that setting x = m + 3 will suffice. Since OS\w + m + 2] and OS[w + m + 1] are the same, we may hide OS[w + m + 2] .

[0300] 3,3,2 Scheduling and Computation of e [0301] Figure 20 shows the data dependency graph, where we can issue a new word of multiplier V) every cycle. In the A task, we compute quotient q and intermediate result ZSRE ₀ and ZCRE _Q . In each of the B tasks, we compute ZSRE _j and ZCRE _j , ZSME _j-1 , and ZCME _j-1 .

[0302] Figure 20(a) shows the schedule when we have 3 PEs to finish the computation in 3 rounds. In the A task of the second round of PE, we generate ZSRE ₀ and ZCRE ₀ and also generate ZSME _e and ZCME _e-1 for the first round. Since we generate w bits every cycle, we actually try to fit intermediate result into ew + m bits. Based on the number of PE instances, we can compute the number of cycles ass shown below:

/7rp + e + 1 if e < p cycle_cnt (83) \ke + p + 1 otherwise

[0303] 3,3,3 Overflow Correction

[0304] Ideally, we can generate w bits of the intermediate result until the last iteration j = e — 1, in which we compute ZSRE _g ^, ZCRE _g ^, ZSME _e-2 , and ZCME _e-2 . If we use the whole m + 2 bits of FBS and FBC to represent ZSME _e-1 and ZCME _e-1 , we can accurately represent the data without overflow. However, it is dangerous to directly take m bits to fit ZSME _g ^ and ZCME _e-1 . Figure 24 shows an example. If we finish the addition, the sum is a negative number, but if we directly use the 4 least significant bits, ZSME _g ^ and ZCME _g ^ become positive and thus cause overflow. We need a method to correct the potential overflow.

[0305] Fortunately, we just need a 2-bit addition to correct the overflow. Figure 25 shows our data model. Assume we have two (N + l)-bit signed number Z _s and Z _c and we know the sum Z is an N- bit signed number, we just need one addition on the leading 2 bits. Since Z is a N- bit signed number, Z[N] must be the sign extension and Z[N ] = Z[N — 1],

[0306] Once we do the addition and ignore any overflow, we can deduce the possible carry

C _in for the least N — l bits as shown in table 7. When (c ₁₍ c ₀ ) = (0,0), C _in has to be 0 and Z[N ] = Z[N — 1] = 0. C _in cannot be 1, otherwise Z[N ] ¹ Z[N — 1] . Similarly, C _in is 1 when (c _lt c ₀ ) = (1,0). The case when (c _± , c ₀ ) = (0,1) is impossible since Z[N\ is always different from Z[N — 1] no matter what C _in is. When (c^ c ₀ ) = (1,1), C _in can be 0 or 1. Except for the case (c^ c ₀ ) = (0,1), the remaining cases will not cause overflow and thus accurately represent the final sum Z.

TABLE 7

Truth table for the overflow correction

[0307] We can use the data model in figure 25 to compute ZSME _g ^ and ZCME _e-1 . redEquation (84) shows that the intermediate result Z^ (after division by r) is an (N + 2m + 2)-bit signed number.

[0308] We thus require ew + m ³ (N + 3m + 2) + 1.

[0309] Assume we instantiate p PEs and go through them k rounds, we equivalently finish kp iterations of the outer loop. Based on equation (78), we have Since e is large enough, we just need to add FBS[m — 1: m — 2\ and FBC[m — 1: m — 2] and then assign the summation results to ZSME _e-1 /ZCME _e-1 .

[0310] 3,3,4 Hardware Implementation

[0311] Figure 22 shows the implementation of processing elements to do tasks A and B. In task A ( CT = 1), we perform the following operations:

• compute Z _L and q and compute q _s , q _c , and q _{carry '■>}

• compute q and compress qM _Q + X _Q Y^ + ZSME _g /r + ZCME _g /r + ZSRE _Q + ZCRE _Q ;

• generate results ZSRE ₀ and ZCRE ₀ ;

• generate results ZSRE _g ^ and ZCRE _e-1 for the previous iteration based on the value stored in FBS and FBC; and

• update FBS and FBC.

[0312] In task B (CT = 0), we perform the following operations:

• compress qM _j + XjY^ + ZSMEj /r + ZCMEj /r + ZSRE _j + ZCRE _j + FBS + FBC based on the stored q;

• generate results ZSRE _j , ZCRE _j , ZSME _j-1 , and ZCME _j-1 ; and

• update FBS and FBC.

[0313] Unlike the A and B tasks, the C task needs only one special processing element (named

PEc). Figure 26 shows the implementation of the PEc. When CT = 1, we force feedback signals DO 1 and DO 2 to 0 and set the carry signal as 0. When CT = 0, we accept feedback signals DO 1 and DO 2 and set the carry signal as the carry-out bit from the previous cycle.

[0314] We also show the block-level architecture of the proposed high-performance, radix-2 ^m scalable Montgomery modular multiplication in figure 23. The required ZSMEj and ZCMEj in the i- th outer loop iteration in realization 8 are actually ZSME _j and ZCME _j which are generated in the ( i — 2)-nd outer loop iteration. Similarly, ZSMEj and ZCMEj required in the PEc are ZSME _j and ZCME _j generated by the next to the last PE. To compute Z _j , we need ZSMEj , ZCMEj , ZSME _j , ZCME _j , ZSRE _j , ZCRE _j . However, they are not generated in the same cycle. More precisely, ZSME _j and ZCME _j are generated in the (J + l)-st inner loop iteration whereas ZSMEj , ZCMEj , ZSRE _j , and ZCRE _j are generated in the j- th inner loop iteration and need to be delayed by one clock cycle as shown in figure 26. When p < e as in figure 21(b), a first- in first-out queue with a depth of e — p slots is necessary to buffer the outputs.

[0315J 3.3.5 Pipelining

[0316] The design in section 7.3.4 finishes each iteration of the inner loop in one cycle. Instead, we can pipeline the processing element in figure 22 into 2 stages or more. Once the connections among different PEs are adjusted accordingly, we can achieve an issue latency (L) of 2 or more. Equation (87) generalizes the number of cycles. if e £ Lp

(87) otherwise

[0317] Although the number of cycles increases when L is more than 1, the required number of PEs (which is optimal for a specified bitwidth N ) is reduced to e/L which results in area savings. If we balance the delay in each stage of the pipeline, the impact on the full circuit latency will be minimized.

[0318] Additional details of some aspects of the invention are set forth in B. Zhang, Z. Cheng and ML Pedram, "High-Radix Design of a Sealable Montgomery Modular Multiplier with Low Latency," in IEEE Transactions on Computers , dot: 10.1109/TC.2G21.3052999; the entire disclosure of which is hereby incorporated by reference.

[0319] While exemplary embodiments are described above, it is not intended that these embodiments describe all possible forms of the invention. Rather, the words used in the specification are words of description rather than limitation, and it is understood that various changes may be made without departing from the spirit and scope of the invention. Additionally, the features of various implementing embodiments may be combined to form further embodiments of the invention. [0320] References:

[0321] [1] Dorian Amiet, Andreas Curiger, and Paul Zbinden. Flexible fpga-based architectures for curve point multiplication over gf (p). In 2016 Euromicro Conference on Digital System Design (DSD), pages 107-114. IEEE, 2016.

[0322] [2] Richard P Brent and Hsiang T Kung. A regular layout for parallel adders. IEEE

Computer Architecture Letters, 31(03):260-264, 1982.

[0323] [3] Viktor Bunimov and Manfred Schimmler. Area and time efficient modular multiplication of large integers. In Proceedings IEEE International Conference on Application- Specific Systems, Architectures, and Processors. ASAP 2003, pages 400-409. IEEE, 2003.

[0324] [4] Serdar Siier Erdem, Tugrul Yamk, and Anil Celebi. A general digit-serial architecture for montgomery modular multiplication. IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 25(5): 1658-1668, 2017.

[0325] [5] Sahar Fatemi, Maryam Zare, Amir Farzad Khavari, and Mohammad Maymandi-

Nejad. Efficient implementation of digit-serial montgomery modular multiplier architecture. IET Circuits, Devices & Systems, 13(7):942-949, 2019.

[0326] [6] Craig Gentry. Fully homomorphic encryption using ideal lattices. In Proceedings of the forty-first annual ACM symposium on Theory of computing, pages 169-178, 2009.

[0327] [7] Zhen Gu and Shuguo Li. A division-free toom-cook multiplication-based montgomery modular multiplication. IEEE Transactions on Circuits and Systems II: Express Briefs, 66(8): 1401-1405, 2018.

[0328] [8] Gael Hachez and Jean-Jacques Quisquater. Montgomery exponentiation with no final subtractions: Improved results. In International Workshop on Cryptographic Hardware and Embedded Systems, pages 293-301. Springer, 2000. [0329] [9] David Harris, Ram Krishnamurthy, Mark Anders, Sanu Mathew, and Steven Hsu.

An improved unified scalable radix-2 montgomery multiplier. In 17th IEEE Symposium on Computer Arithmetic (ARITH’05), pages 172-178. IEEE, 2005.

[0330] [10] Miaoqing Huang, Kris Gaj, and Tarek El-Ghazawi. New hardware architectures for montgomery modular multiplication algorithm. IEEE Transactions on computers , 60(7):923- 936, 2010.

[0331] [11] Neal Koblitz. Elliptic curve cryptosystems. Mathematics of computation,

48(177):203-209, 1987.

[0332] [12] Ciaran Mclvor, Maire McLoone, and John V McCanny. Fast montgomery modular multiplication and rsa cryptographic processor architectures. In The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003, volume 1, pages 379-384. IEEE, 2003.

[0333] [13] NIEL Michael. Montgomery multiplication method, February 17 2015. US

Patent 8,959,134.

[0334] [14] Peter L Montgomery. Modular multiplication without trial division.

Mathematics of computation, 44(170):519-521, 1985.

[0335] [15] Miguel Morales-Sandoval and Arturo Diaz-Perez. Scalable gf (p) montgomery multiplier based on a digit-digit computation approach. IET Computers & Digital Techniques, 10(3): 102-109, 2016.

[0336] [16] Michael Andrew Moshier and Jeff Furlong. Scalable, faster method and apparatus for montgomery multiplication, September 282010. US Patent 7,805,479.

[0337] [17] Michael Naehrig, Kristin Lauter, and Vinod Vaikuntanathan. Can homomorphic encryption be practical? In Proceedings of the 3rd ACM workshop on Cloud computing security workshop, pages 113-124, 2011. [0338] [18] Bien-Cuong Nguyen and Cong-Kha Pham. An efficient hardware implementation of radix- 16 montgomery multiplication. In 2019 IEEE 8th Global Conference on Consumer Electronics (GCCE), pages 1121-1122. IEEE, 2019.

[0339] [19] Holger Orup. Simplifying quotient determination in high-radix modular multiplication. In Proceedings of the 12th Symposium on Computer Arithmetic, pages 193-199. IEEE, 1995.

[0340] [20] Erdem Ozcan and Serdar S Erdem. A fast digit based montgomery multiplier designed for fpgas with dsp resources. Microprocessors and Microsystems, 62:12-19, 2018.

[0341] [21] Jeng-Shyang Pan, Pengfei Song, and Chun-Sheng Yang. Efficient digit-serial modular multiplication algorithm on fpga. IET Circuits, Devices & Systems, 12(5):662-668, 2018.

[0342] [22] Ronald L Rivest, Len Adleman, Michael L Dertouzos, et al. On data banks and privacy homomorphisms. Foundations of secure computation, 4(11): 169-180, 1978.

[0343] [23] Ronald L Rivest, Adi Shamir, and Leonard Adleman. A method for obtaining digital signatures and public-key cryptosystems. Communications of the ACM, 21(2): 120-126, 1978.

[0344] [24] Sudhir K Satpathy, Raghavan Kumar, Sanu K Mathew, and Vikram B Suresh.

Parallel computation techniques for accelerated cryptographic capabilities, December 3 2019. US Patent 10,498,532.

[0345] [25] Ming-Der Shieh and Wen-Ching Lin. Word-based montgomery modular multiplication algorithm for low-latency scalable architectures. IEEE transactions on computers, 59(8): 1145—1151, 2010.

[0346] [26] Alexandre F Tenca and £etin K KoV. A scalable architecture for modular multiplication based on montgomery’s algorithm. IEEE Transactions on computers, 52(9): 1215— 1221, 2003. [0347] [27] Colin D Walter. Montgomery’s multiplication technique: How to make it smaller and faster. In International Workshop on Cryptographic Hardware and Embedded Systems, pages 80-93. Springer, 1999.

[0348] [28] Tao Wu, ShuGuo Li, and LiTian Liu. Fast rsa decryption through high-radix scalable montgomery modular multipliers. Science China Information Sciences, 58(6): 1-16, 2015.