Big integer power modulo

I. Problem description:

Calculate (a ^ power) % m, where power is a non-negative big integer, A, M is an integer greater than 1.

Ii. Problem Analysis:

Obviously, power is a big integer, so the overflow of Power computing must be considered. How can we avoid overflow? This can be achieved by reducing the power and successive modulo operations. A natural idea is to divide power into the sum of two integers, power = N1 + N2, then a ^ power = (a ^ N1) * (a ^ N2 ). usually use the bipartite, n1 = n2 or | n1-n2 | = 1. This involves the calculation of (a * B) % m.

It is not hard to prove that (A * B) % m = (a % m) * (B % m) % m.

Proof: If a = PM + R1, B = QM + R2 is set, R1 = A % m, R2 = B % m,

Then (A * B) % m = (R1 * R2) % m = (a % m) * (B % m) % m. The certificate is complete.

In particular, when a = B, (a ^ 2) % m = (a % m) ^ 2) % m; this is a simple but key conclusion.

Iii. Algorithm Design:

Algorithm 1: grouping policy:

(1) If power is an odd number: Make power = 2 k + 1, K is a non-negative integer, then

A ^ (2 k + 1) = (a ^ K) ^ 2 * A; A ^ (2 k + 1) % m = (a ^ K % m) ^ 2% M * A) % m
(2) If power is an even number: If power is equal to 2 K and K is a non-negative integer

A ^ (2 k) = (a ^ K) ^ 2; a ^ (2 k) % m = (a ^ K % m) ^ 2% M
(3) If power = 1: returns a % m; If power = 0: returns 1% M.

Based on the above (1)/(2)/(3), the corresponding recursive solving program can be designed. The time complexity is T (n) = T (n/2) + C = O (logn), and the space complexity is O (1). The disadvantage is that there is a certain amount of recursive call overhead.

Algorithm 2: Binary decomposition of Integers

Separate the big integer power in binary format:

Power = A [n] * 2 ^ N + A [N-1] * 2 ^ (N-1) +... + A [1] * 2 + A [0]

Where: A [I] is 0 or 1 (I =,..., n), then

A ^ power = a ^ (A [n] * 2 ^ N) * A ^ (A [N-1] * 2 ^ (N-1 ))*... * A ^ (A [1] * 2) * A ^ A [0]
Obviously:

(1) If a [I] = 0, a [I] * 2 ^ I = 0, a ^ (A [I] * 2 ^ I) = 1, this product item can be deleted;

(2) If a [I] = 1, a [I] * 2 ^ I = 2 ^ I, a ^ (2 ^ I) % m = (a ^ (2 ^ (I-1) % m) ^ 2% M.

Make temp = a ^ (2 ^ (I-1) % m, then a ^ (2 ^ I) % m = (temp * temp) % m.

For example, if power = 26 = 16 + 8 + 2 = (11010) 2, a ^ 26 = a ^ (2 ^ 4 + 2 ^ 3 + 2 ^ 1 );

Computing a ^ power is actually the product of computing power in binary representation where all bits correspond to 1.

(A ^ 26) % m = (a ^ 16% m) * (a ^ 8% m) * (a ^ 2% M) % m

However, a ^ 8% M = (a ^ 4% m) × (a ^ 4% m) % m can be solved successively by dynamic programming. For simplicity, (a ^ I) % m is called "modulo product item ".

Algorithm Description:

Make temp = a ^ (2 ^ I) % m, and I is the position corresponding to the current binary bit of power,

Temp indicates(Modulo)Product

Step 1: Initialize temp as a % m, result = 1;

Step 2: Perform binary decomposition on power. If power> = 1, perform the modulo 2 operation. Otherwise, go to step 3.

[1] if the remainder is 1, the binary value at this position is 1. In the product, add the temp item: Result = (result * temp) % m;

The product item corresponding to the next binary is temp = (temp * temp) % m

[2] if the remainder is 0, the binary value at this position is 0. The temp item is not required in the product, and the result remains unchanged,

The product item corresponding to the next binary is temp = (temp * temp) % m

Step 3: The product items corresponding to all binary bits are calculated and the algorithm ends.

For example, the calculation process of result = (3 ^ 26) % 5 is as follows: 26 = (11010) 2;

(1) initialization: temp = 3 ^ 1% 5 = 3;, result = 1;

(2) If the binary value of detection 26 is 0, no product item is added. Result = 1, temp = (3 ^ 2) % 5 = (temp * temp) % 5 = 4

(3) If the secondary low binary value of detection 26 is 1, add the product item, result = (result * temp) % 3 = 4, temp = (3 ^ 4) % 5 = (4*4) % 5 = 1;

(4) If the secondary low binary value of detection 26 is 0, no product item is added. Result = 4, temp = (3 ^ 8) % 5 = (1*1) % 5 = 1;

(5) If the secondary low binary value of detection 26 is 1, add the product item, result = (result * temp) % 5 = 4, temp = (3 ^ 16) % 5 = 1;

(6) If the secondary low binary value of detection 26 is 1, add the product item, result = (result * temp) % 5 = 4, temp = (3 ^ 32) % 5 = 1.

As all bits are detected, 3 ^ 26% 5 = 4. We can see from the above,

3 ^ 26% 5 = (3 ^ 16) % 5) * (3 ^ 8) % 5) * (3 ^ 2) % 5) % 5. 3 ^ 16% 5, 3 ^ 8% 5, and 3 ^ 2% 5 are obtained through dynamic programming.

 

The complete procedure is as follows:

The program is compiled and run in ubuntu10.10 gcc4.4.5.
$ Gcc-g-wall bintmode. C runtime. C-o bintmode # compile the connection
$ Bintmode or GDB bintmode # Run

 

/** Bintmode. c: This program calculates (a ^ power) % Mod. power is a big integer * basic formula: (A * B) mod m = (a mod m) * (B mod m) mod m */# include <stdio. h> # include <assert. h> int modmrec (INT, Int, INT); int modmdyn (INT, Int, INT); void testmodmrec (INT); void testmodmdyn (INT ); void testvalid (INT (* Fun) (int A, int power, int mode); void testinvalid (INT (* Fun) (int A, int power, int mode )); int main () {printf ("******** solve with recursive technology :******** \ N "); testvalid (modmrec);/* testinvalid (modmrec); */defruntime (testmodmrec ); printf ("********* solve with binary decomposition technology: ********** \ n"); testvalid (modmdyn ); /* testinvalid (modmdyn); */defruntime (testmodmdyn); Return 0;}/** modmrec: recursive solution (a ^ power) % mod, power is a big integer */INT modmrec (int A, int power, int mod) {assert (A> = 1 & Power> = 0 & mod> = 1 ); if (power = 0) {return 1% MOD;} If (power = 1) {return % MOD;} If (Power % 2! = 0) {int temp = modmrec (A, power/2, MoD); Return (temp * A) % MOD;} else {int temp = modmrec (, power/2, MoD); Return (temp * temp) % mod ;}/ ** modmdyn: use binary decomposition technology to solve (a ^ power) % mod, power is a big integer */INT modmdyn (int A, int power, int mod) {assert (A> = 1 & Power> = 0 & mod> = 1 ); int result = 1; int temp = A % MOD; while (power> = 1) {If (Power % 2! = 0) {result = (result * temp) % MOD;} temp = (temp * temp) % MOD; power/= 2;} return result ;} void testvalid (INT (* Fun) (int A, int power, int mode) {int base [5] = {2, 3, 5, 7, 11}; int I, K; for (I = 0; I <5; I ++) {for (k = 0; k <20; k ++) {printf ("% d ^ % d mod % d = % d. \ n ", base [I], K, 10, (* Fun) (base [I], K, 10) ;}} void testinvalid (INT (* Fun) (int A, int power, int mode) {printf ("0 ^ 3 mod 5"); (* Fun) (0, 3, 5 ); printf ("3 ^-1 mod 5"); (* Fun) (3,-1, 5); printf ("3 ^ 5 mod 0"); (* Fun) (3, 5, 0);} void testmodmrec (int n) {modmrec (2, N, 10);} void testmodmdyn (int n) {modmdyn (2, n, 10 );}

# Include <stdio. h> # include <time. h> # define max_scale 10 long defscale [max_scale] = {1, 10,100,100 0, 10000,100 000, 1000000,100 00000, 100000000,100 0000000};/* runtime: measure the running time of the function that fun points to * Fun: pointer to the test function; * scale: array of problem scales */void Runtime (void (* Fun) (Long Scale ), long * scale, int Len); void defruntime (void (* Fun) (Long Scale); void Runtime (void (* Fun) (Long Scale), long * scale, int Len) {int I; clock_t start, end; for (I = 0; I <Len; I ++) {start = clock (); (* Fun) (scale [I]); End = clock (); printf ("scale: % 12ld \ t run time: % 8.4f (MS ). \ n ", scale [I], (double) (end-Start) * 1000/clocks_per_sec) ;}} void defruntime (void (* Fun) (Long Scale )) {Runtime (fun, defscale, max_scale );}

 

Iv. Summary:Scientific computing algorithms
1. it is often necessary to read relevant materials, understand the corresponding mathematical nature and conclusions, and simplify the problem. For example, in this example, (A * B) % m = (a % m) * (B % m) % m modulo properties to avoid multiplication overflow.
2. The sub-governance method and binary decomposition method are preferred. The division method halved the number to be solved, while the binary decomposition sought a solution from a special perspective of the number.