Aritmética Computacional

Francisco Rodríguez Henríquez Aritmética Computacional

Francisco Rodríguez HenríquezCINVESTAV

e-mail: [email protected]

Aritmética Computacional


Fairy Tale: Chinese Emperor used to count his army by giving a series of tasks.

1. All troops should form groups of 3. Report back the number of soldiers that were not able to do this.

2. Now form groups of 5. Report back.3. Now form groups of 7. Report back.4. Etc.At the end, if product of all group numbers is sufficiently

large, can ingeniously figure out how many troops.

Chinese Remainder Theorem




mod 3:

N mod 3 = 1



mod 5:

N mod 5 = 2



mod 7:

N mod 7 = 2



Secret inversion formula (for N < 105 = 3·5·7):N a (mod 3)N b (mod 5)N c (mod 7)

Implies that N = (-35a + 21b + 15c) mod 105.So in our case a = 1, b = 2, c = 2 gives:N = (-35·1 + 21·2 + 15·2) mod 105

= (-35 + 42 + 30) mod 105= 37 mod 105 = 37



CRT: Example

1. Find three numbers l,m,n with following properties– l 1(mod 3), l 0(mod 5), l 0(mod 7)– m0(mod 3), m 1(mod 5), m 0(mod 7)– n 0(mod 3), n 0(mod 5), n 1(mod 7)

2. Then y = al+bm +cn [secret formula] satisfies– y al+bm +cn (mod 3)

a·1+0 + 0 (mod 3) a (mod 3) – Similarly, y b (mod 5) – Similarly, y c (mod 7)

3. This will imply x y (mod 3·5·7)


Find three numbers l,m,n: Standard trick.EG, to find l :a) Multiply together all modulii different from 3.

Result: 5·7 = 35b) Find an inverse of this number mod 3: In this

case it’s easy. 35 2(mod 3) so find an inverse of 2 [2 or anything congruent to 2(mod 3)]. Practice shows that should choose inverse of smallest magnitude: –1.

c) l is the product of (a) and (b): l = -35l is 0 mod 5 and 7 since it’s divisible by 5·7. But (c)

guarantees that it’s 1 modulo 3!

CRT: Example


Similarly, m = 21 and n = 15. So our solution to all three congruences is:

x = -35a + 21b + 15cIf we want to guarantee a solution between 0 and

104, just computex mod 105 .

The same tricks can be generalized to prove:

CRT: Example


THM (CRT): Let m1, m2, … , mn be pairwise relatively prime positive integers. Then there is a unique solution x in [0,m1·m2···mn-1] to the system of congruences:

x a1 (mod m1 )

x a2 (mod m2 )

x an (mod mn )



CRT: Conversion Algorithm

Step 1. Compute using multi-precision arithmetic.

Step 2. Compute the multiplicative inverses of modulo mi for 1 ≤ i ≤ n, i.e., compute the constants ci

such that,

Step 3. Compute u by performing the sum (in multiprecision arithmetic):

inii m

MmmmmmM 121

imM

.1for ,mod1 nimcmM

iii

MucmMuc

mMuc

mMu nn

n

mod222

111


CRT: Conversion Algorithm

Theorem. Given the moduli m1, m2,…, mn and

the remainders u1, u2,…, un the number u can

be computed in O(n2).


CRT: Mixed-Radix Conversion Algorithm

Step 1. Compute constants cij for 1 ≤ i < j ≤ n such that,

Step 2. Compute

Step 3. Compute

jiij mmc mod1

nnnnnnnn mcvcvcvuv

mcvcvuvmcvuv

muv

mod

,mod,mod

,mod

,112211

223213133

212122

111

121213121 nn mmmvmmvmvvu



Computation of u using the above formula also requires O(n2) arithmetic operations. We now define Vij for

0 ≤ i < j ≤ n such that Voi = ui for 1 ≤ i ≤ n. These Vij are the temporary values of vj resulting from the operations in Step 2 of the mixed-radix conversion algorithm. This way, we build a triangular table of values with diagonal entries Vi = Vi-1,j for 0 ≤ i ≤ n. The entries of this table are named multiplied differences.



An Example: For n = 4, it can be given as follows,

Where [mi] stands for modulo mi.

4342324344241214244140104144404

3231213233130103133303

2120102122202

1101

mcVVVmcVVVmcVVVmuVmcVVVmcVVVmuV

mcVVVmuVmuV


Finite fields: Arithmetic operations

FP finite field operations : Addition, subtraction,

multiplication, Squaring, inversion, exponentiation and Primality Testing


Arithmetic Operations in GFp

Operation Bit Complexity

Addition a + b mod n O(lg a + lg b) = O (lg n)

Subtraction a – b mod n O(lg a + lg b) = O (lg n)

Multiplication a*b mod n O(lg a lg b) = O ((lg n)2)

Inversion a-1 mod n O ((lg n)3)

Exponentiation ak mod n O ((lg n)3)


Modular Addition and Subtraction


Modular Addition

Input: A modulus p, and integers a, b in [0, p-1]Output: c = (a + b) mod p.1. C0 = Add(a0, b0);

2. For i from 1 to t-1do: Ci = Add_with_carry(ai, bi);

3. If the carry bit is set, then subtract p from c = (ct-1,…, c2,c1,c0). (why??)

4. If c ≥ p then c -= p; (why??)5. Return(c);


Modular Subtraction

Input: A modulus p, and integers a, b in [0, p-1]Output: c = (a - b) mod p.1. C0 = Subtract(a0, b0);2. For i from 1 to t-1do: Ci = Subtract_with_borrow(ai, bi);

3. If the carry bit is set, then add p to c = (ct-1,…, c2,c1,c0). (why??)

4. Return(c);


Modular Multiplication


Modular Multiplication

Computation of c = ab mod n can be performed by using:• Classical: Normal integer multiplication followed by

reduction• Blakley’s method: The multiplication steps are

interleaved with reduction steps.• Montgomery’s method: Uses predominantly modulo

2j arithmetic.


Modular Multiplication: Classical Method


Integer Multiplication

We perform the operations radix W = 2w: wordsize of the computer:

We define (Carry, Sum) pairs. Our notation is:

1

0021

1

0021

s

j

iiss

s

j

iiss

WbWbbbb

WaWaaaa

:jiij abt

12,,1,0for :1,,1,0for :,

sittsibaba

i

ii



01234567

30313233

20212223

10111213

00010203

0123

0123

tttttttttttt

tttttttt

ttttbbbbaaaa



1. for i = 0 to s-1 do:2. C:= 0

3. for j = 0 to s-1 do:4. (C, S) := ti+j + ajbi + C;5. ti+j := S;

6. end7. ti+j+1:= C;

8. end



i j Step (C, S) Partial t0 0 t0 + a0b0 + C

0 + 87 + 0

(0, *)(5, 6)

000000000006

1 t1 + a1b0 + C

0 + 37 + 3 (3,3) 000036

2 t2 + a2b0 + C

0 + 37 + 3 (2, 4) 000436

002436




3 + 85 + 0

(0, *)(4, 3) 002436

1 t2 + a1b1 + C

4 + 45 + 4 (2, 8) 002836

2 t3 + a2b1 + C

2 + 35 + 2 (1, 9) 009836

019836




8 + 88 + 0

(0, *)(7, 2) 019236

1 t3 + a1b2 + C

9 + 48 + 7 (4, 8) 018236

2 t4 + a2b2 + C

1 + 38 + 4 (2, 9) 098236

298236



This algorithm requires s2 = (k/w)2 inner product steps: (C, S) := ti+j+ajbi+C;

In other words, O(k2) bit operations.The variables ti+j, aj, bi, C and S each hold a single-

word, or a w-bit number.Notice that from the main operation in the loop we

obtain a double-word, or a 2w-bit number since:

1212121212 2 WWWWW


A straightforward modification of the multiplication algorithm gives the following algorithm for squaring. There are roughly ½ fewer multiplication operations.

Integer Squaring


Integer Squaring [Guajardo and Paar]

Input: An integer a [0, p-1], a = (at-1 at-2 … a1 a0)Output: c = a2.1. for i from 0 to 2t-1 do: ci = 0;2. for i from 0 to t-1 do

3. (uv) = c2i + ai2;

4. C2i=v; C1= u; C2 = 0;5. for j from i+1 to t-1 do

6. (uv) = ci+j + ai aj + C1; C1 = u;7. (uv) = v + ai aj + C2; ci+j = v ; C2 = u;

8. (uv) = C1+C2, C2 = u;9. (uv) = ci+t + v; ci+t= v;10.ci+t+1 = C2 + u;

11. return (c);


Integer Squaring [Classical]

Input: An integer a [0, p-1], a = (at-1 at-2 … a1 a0)Output: c = a2.1. r0 = r1 = r2 = 0;2. for k from 0 to 2(t-1) do

3. For each elmt. of {(i, j)| i+j = k, 0 ≤ i ≤ j < t} do4. (uv) = ai aj;

5. If (i < j) then (uv) << 1; r2 = AddC(r2, 0);

6. r0 = Add(r0, v); r1 = AddC(r1, u); r2 = AddC(r2, 0);

8. ck = r0; r0 = r1; r1 = r2; r2 = 0;

9. c2t-1 = r0;11. return (c);


Reduction

Given t, the computation of R which satisfiest = Qn + R

With R < n. Here t is a 2k-bit number and n is a k-bit number.

The number t and n are positive, so are the results Q and R.

Since we are not interested in the quotient, steps of the division algorithm can be simplified.


Reduction

Two algorithms of interest:

• Restoring Division

• Non-restoring division


Restoring Division

1. R0 := t;2. n := 2kn;

3. for i = 1 to k do:4. Ri := Ri-1-n;5. if Ri<0 then Ri := Ri-1;

6. n := n/2;

6. end7. Return Rk;


Restoring Division: An example

• We give an example of the restoring division algorithm for computing 3019 mod 53, where,

3019 = (101111001011)2

53 = (110101)2

The result is:51 = (110011)2



R0 101111 001011 tn 110101 Subtract

-000110 Negative RemainderR1 101111 001011 Restoren/2 11010 100000 Subtract

+10100 100000 Positive RemainderR2 10100 101011 Not restoren/2 1101 010000 Subtract

+0111 010000 Positive rem.R3 0111 011011 Not restoren/2 110 101000 Subtract



+000 110000 Positive remainderR4 000 110011 Not Restoren/2 11 0101n/2 1 10101n/2 0 110101 Subtract

- 000010 Negative RemainderR5 110011 RestoreR 110011 Final Remainder


Non restoring Division Algorithm

• The non-restoring division algorithm allows a negative remainder.

• Suppose Ri:=Ri-1-n< 0, then the restoring algorithm assigns Ri:=Ri-

1 and performs a subtraction with the shifted n, obtaining Ri+1:= Ri-

n/2 = Ri-1-n/2;

• However, if Ri = Ri-1 – n < 0, then the non-restoring algorithm lets

Ri remain negative and adds the shifted n in the following cycle. Thus it obtains,

Ri+1:= Ri+n/2 = (Ri-1-n)+n/2 = Ri-1-n/2;

i.e., the same value (!!)


Non-Restoring Division Algorithm

1. R0 := t;2. n := 2kn;

3. for i = 1 to k do:4. if Ri-1<0 then Ri := Ri-1-n;

5. else Ri := Ri-1+n;

6. n := n/2;

6. end7. Return Rk;


Non-Restoring Division Algorithm

• Since the remainder is allowed to stay negative, we use 2’s complement coding to represent such numbers.

• Also, note that the nonrestoring division algorithm may require a final restoration cycle in which a negative remainder is corrected by adding the last value of n back to it.

• Example Computation of 51 = 3019 mod 53.



R0 101111 001011 tn 110101 Subtract

1111010 Negative Remaindern/2 011010 100000 addR2 010100 100000 Positive remaindern/2 1101 010000 SubtractR3 0111 010000 Positive remaindern/2 110 101000 SubtractR4 0000 110 Positive remaindern/2 011 010100n/2 01 101010



n/2 110101 subtractR5 1 111110 Negative Remainder

n 0 110101 Add (restore)R 110011 Final Remainder


Barrett Reduction

Barrett reduction computes r = x mod m given x and m. The algorithm requires the precomputation of the quantity,

It is advantageous if many reductions are performed with a single modulus. Typically, the radix b is chosen to be a power of two closed to the word-size of the processor.

Barrett reduction is based on the following fact:

Given

pb k2

121 /1//

as, written becan ,0 and

kkk bpbbxQ

pxpRRQpx


Barrett Reduction

Input: positive integers x = (x2k-1 … x1x0), p = (pk-1 … p1p0)

Output: x mod p.1. 2. 3. if r < 0 then 4. While r ≥ p do: r= r-p; 5. Return(r);

;//ˆ 11 kk bbxq ;modˆmod 11 kk bpqbxr

1 kbrr

pbbxpkpb

kk

b

22 ,0,1log,,3


Barrett Reduction

Example: Let b = 4, k = 3, x = (313221)b, and p = (233)b (i.e.,

x = 3561, and p = 47). Then = |46/p| = 87 = (1113)b,

|x/bk-1| = |(313221)b/42| = (3132)b,

|x/bk-1| = (3132)b (1113)b = (10231302)b

Hence q = (1023)b,

r1 = (3221)b (why??)

r2 = (1023)b (233)b mod b4 =(3011)b, and r = r1 – r2 = (210)b

Thus x mod p = (210)b = 36


Barrett Reduction : Computational efficiency

• All divisions performed in the algorithm are

simple right-shifts of the base b representation.

• Since the k+1 MSBs of x/bk-1| are not needed

to determine q (why??), only a partial multiple-

precision multiplication is necessary.


Reduction

The arithmetic in Barrett reduction can be reduced by choosing b to be a power of 2. For primes p of special form, there exist very fast modular reduction techniques [For example, see “Software Implementation of the NIST Elliptic Curves Over Prime Fields”, Brown, Hankerson, López and Menezes].


Modular Multiplication: Blakley’s Method


Blakley’s Method

Let ai and bi represent the bits of the k-bit numbers a and b, respectively. The product t (2k-bit number) can be written as,

This formulation yields the shift-add multiplication algorithm. Blakley’s algorithm uses this formulation and furthermore reduces the partial product modulo n at each step.

ik

ii

k

i

ii bababat 22

1

0

1

0


Blakley’s Method

1. R := 0;

2. For i = 0 to k-1do

3. R := 2R + ak-1-ib;

4. R := R mod n;

5. End

6. Return R;


Blakley’s Method

Assuming that 0 ≤ a, b, R ≤ n-1, the new R will be in the range 0 ≤ R ≤ 3n – 3

SinceAt most two subtraction will be needed to bring the

new R to the range [0, n - 1]. Thus we can use While (R ≥ n) R -= n;Blakley’s algorithm computes the remainder R in k

steps, where at each step one left shift, one addition, and at most two subtractions are performed; the operands involved in these computations are of length k bits.

331122: nnnbarR j


Modular Multiplication: Montgomery’s Method


Montgomery’s Method

This method replaces division by n operations with division by r = 2k. Assuming n is a k-bit integer, i.e., 2k-1 < n < 2k

We assign r = 2k. Now, we perform the mapping of the integers a [0, n-1] to the integers [0, n-1] using the one-to-one mapping

We call the n-residue of a.

nraa mod: a

a



We now define the Montgomery product of two n-residues as

Also we need n’ such that rr-1-nn’ = 1; r-1 and n’ are computed by using the extended Euclid’s algorithm.

nrbabaoMon mod),(Pr 1

; else )( then if

/: .3mod: 2.;: .1,Pr

ureturnnureturnnurnmtu

rntmbat

baoMon



This routine requires only modulo r arithmetic, which is efficiently accomplished on a computer if r = 2j.

Theorem 1. If c = ab mod n thenProof:

);,(Pr baoMonc

),(Prmod

mod

modmod

1

1

baoMonnrba

nrrbra

nrbanrcc



Theorem 2.

Proof:

)1,(Pr coMonc

)1,(Prmod1

mod

mod

1

1

1

coMonnrc

nrc

nrrcc



MonPro procedure can be utilized to compute c: =ab mod n as follows:ModMul(a, b, n) /* n is odd (why???) */1. Compute n’ using EEA.2. 3. 4. 5. 6. Return c;

nraa modnrbb mod

baoMonc ,Pr: 1,Pr: coMonc



Since preprocessing operations such as,

• Computation of n’ and,

• Conversion from ordinary to n-residue

• Conversion from n-residue to ordinary

Are time consuming, it is not a good idea to use Montgomery’s method for a single modular multiplication. However, it is very suitable for modular exponentiation.



MonPro procedure can be utilized to compute c: =Me mod n as follows:ModExp(M, e, n) /* n is odd (why???) */1. Compute n’ using EEA.2. 3. 4. for i=k-1 down to 0 do 5. 6. If ei = 1 then

7. C := 8. Return C;

nrMM modnrC mod1

ccoMonc ,Pr: CMoMonc ,Pr:

This function uses the Binary method that will be discuss in detail later. AnyOther exponential algorithmwill work as well.

);1,(Pr CoMon



Example: Computation of 710 mod 13r = 2k=16. Since 16*9-13*11 = 1, we have r-1 = 9, n´= 11.

M = 7, thus C = 1, thusHence,

813mod167mod: nrCM

813mod167mod: nrCM

8M and 3 C

ei Step 5 Step 6

1 MonPro(3, 3) = 3 MonPro(8, 3) = 80 MonPro(8, 8) = 41 MonPro(4, 4) = 1 MonPro(8, 1) = 70 MonPro(7, 7) = 12



Step 7: C = MonPro(12, 1) = 4Computation of MonPro(3, 3):t := 3*3 = 9;m := 9*11 mod 16 = 3u := (9+3*13)/16=48/16=3

Computation of MonPro(8, 1):t := 8*1 = 8;m :=8*11 mod 16 = 8u :=(8+8*13)/16=112/16=7

Documents

Aritmética Computacional