Rabin-Miller Algorithm

First, I will attach the matrix67 explanation:

Certificate -----------------------------------------------------------------------------------------------------------------------------------------------------

The work of Miller and Rabin has taken a revolutionary step in the Fermat test and established the legendary Miller-Rabin algorithm. The new test is based on the following theorem: If P is a prime number, X is a positive integer smaller than P, and x ^ 2 mod p = 1, then either x = 1 or X = p-1. This is obvious, because x ^ 2 mod p = 1 is equivalent to P which can divide x ^ 2-1, that is, P can divide (x + 1) (x-1 ). Since P is a prime number, it is only possible that the X-1 can be divisible by P (at this time x = 1) or x + 1 can be divisible by P (at this time x = PM ).
The following example shows how the theorem is applied to the Fermat test. As mentioned above, 341 can pass a 2-based Fermat test, because 2 ^ 340 mod 341 = 1. If 341 is a prime number, 2 ^ 170 mod 341 may only be 1 or 340. If 2 ^ 170 mod 341 is equal to 1, we can continue to view the result of dividing 2 ^ 85 by 341. We found that 2 ^ 85 mod 341 = 32, this result removed the prime crown on the head of 341, and the real face behind the mask appeared, the attempt to pretend to be a prime number to interact with our plain mm was revealed.
This is Miller-Rabin's method for the testing of the kernel. Continuously extract factor 2 from the exponent N-1 and express n-1 as D * 2 ^ r (where D is an odd number ). What we need to calculate is the remainder of D * 2 ^ r of a divided by N. So a ^ (D * 2 ^ (r-1) is either equal to 1 or n-1. If a ^ (D * 2 ^ (r-1) is equal to 1, the theorem continues to apply to a ^ (D * 2 ^ (R-2, the A ^ d mod n = 1 or N-1 is obtained after an I satisfies a ^ (D * 2 ^ I) mod n = n-1 or 2 in the final index. In this way, the Fermat theorem is enhanced in the following form:
Extract factor 2 as much as possible and express n-1 as D * 2 ^ r. If n is a prime number, or a ^ d mod n = 1, or there is an I so that a ^ (D * 2 ^ I) mod n = N-1 (0 <= I <R) (note that I can be equal to 0, this will unify the situation of a ^ d mod n = n-1 to the end)
Miller-Rabin is also an uncertain algorithm. We call the sum of Miller-Rabin tests based on A as the base-strong Pseudo Prime (strong Pseudo Prime ). The first base-2 strong Pseudo Prime Number is 2047. The first strong Pseudo Prime Number base on 2 and 3 is as large as 1 373 653.
The code of the Miller-Rabin algorithm is also very simple: Calculate the values of D and R (which can be accelerated by bitwise operations), and then calculate the values of a ^ d mod N in two parts, finally, it is squared to R times. The code of the program is simpler than imagined. I write a copy and put it below. Although I have already switched to C, I believe there are still many people who do not understand C. Let me write Pascal again. The isprime function returns whether the specified base number A and N can pass the test. If the function returns false, n is not a prime number. If the function returns true, n is very likely a prime number.Note that the data range of this Code is limited to longint. You may need to change them to int64 or high-precision computing.
function pow( a, d, n:longint ):longint;
begin
   if d=0 then exit(1)
   else if d=1 then exit(a)
   else if d and 1=0 then exit( pow( a*a mod n, d div 2, n) mod n)
   else exit( (pow( a*a mod n, d div 2, n) * a) mod n);
end;

function IsPrime( a,n:longint ):boolean;
var
   d,t:longint;
begin
   if n=2 then exit(true);
   if (n=1) or (n and 1=0) then exit(false);
   d:=n-1;
   while d and 1=0 do d:=d shr 1;
   t:=pow( a, d, n );
   while ( d<>n-1 ) and ( t<>1 ) and ( t<>n-1 ) do
   begin
      t:=(t * t)mod n;
      d:=d shl 1;
   end;
   exit( (t=n-1) or (d and 1=1) );
end;
Miller-Rabin is currently the most widely used algorithm for the determination of large numbers. Generally, the base number is still randomly selected, but when the number of tests is not large, there are some tips to select the base number. For example, if the number to be tested is less than 4 759 123 141, you only need to test three bases 2, 7, and 61. Of course, the more you test, the larger the correct range. If you use the first seven prime numbers (2, 3, 5, 7, 11, 13, and 17) each time, all the numbers that do not exceed 341 550 071 728 320 are correct. If the base number is 2, 3, 7, 61, and 24251, the unique strong Pseudo Prime Number in 10 ^ 16 is 46 856 248 255 981. These conclusions make the Miller-Rabin algorithm very useful in Oi. It is generally believed that the accuracy rate of Miller-Rabin testing is acceptable. The loss rate of the random selection of K bases for testing is approximately 4 ^ (-K ).

The Miller-Rabin algorithm is an RP algorithm. RP is a kind of time complexity, mainly for the determination problem. One algorithm is the RP algorithm, which indicates that it can be completed in polynomial time and can accurately judge the case where the answer is no, but at the same time, it may also make the right sentence wrong (the error probability cannot exceed 1/2 ). The RP algorithm is based on randomization, so running this algorithm multiple times can reduce the error rate. There are other probabilistic testing algorithms, such as the solovay-strassen algorithm. In addition, some native testing algorithms need to know some auxiliary information (such as n-1 prime factor) in advance, or the number to be tested must meet certain conditions (for example, the number to be tested must be in the form of 2 ^ N-1 ). The AKs algorithm has made a sensation in the world over the past few years. It is the first Polynomial-based, deterministic, and non-conditional prime judgment algorithm. At that time, a paper was published with the title primes is in P, and the whole world went crazy. There were several mm in our class that came to the initial tide. The algorithm is mainly based on the following fact: N is a prime number when and only when (X-a) ^ n squared (x ^ N-A) (mod N ). Note that X is an unknown number in the polynomial, and each side of the equation is a polynomial. For example, when a = 1, the proposition is equivalent to the following conclusion: WHEN n is a prime number, the n + 1 line in the Yang Hui triangle can be divisible by N except the 1 at both ends.

Certificate ------------------------------------------------------------------------------------------------------------------------------------------------

The quick power here is ugly...

There is no problem above. We will mainly study the sentence exit (t = N-1) or (D and 1 = 1;

Or is followed by a ^ d mod p = 1 or P-1

In the above section, it means that a ^ d is in the process of repeating the square step, for example, mod p = p-1. why P-1? Why is it no longer or T = 1?

In fact, if a ^ d is in the process of repeating the square for a mod p = 1, then x ^ 2 mod p = 1?

Then x = 1 or P-1, And because X is not equal to 1, we first tested whether a ^ d mod p = 1 or P-1,

If it is not, it will first appear in the Process of repeating the square, instead of 1, and if there is 1, it must have appeared before. Otherwise, it will continue to push forward, A ^ d mod p = 1, which is in conflict with the WHILE LOOP

So there is no problem with such an algorithm.

My code:

 1 var p,n,m:int64; t:longint; 2 procedure init; 3  begin 4    readln(p); 5  end; 6 procedure mul(x,y:int64;var z:int64); 7  var tmp:int64; 8  begin 9    tmp:=0;10    while y>0 do11     begin12      if odd(y) then tmp:=(tmp+x) mod p;13      y:=y>>1;14      x:=(x+x) mod p;15     end;16    z:=tmp;17  end;18 19 function check(x:int64):boolean;20   var cs,y,z:int64;21       j:longint;22   begin23     cs:=n;y:=1;24     while cs>0 do25       begin26         if odd(cs) then mul(y,x,y);27         cs:=cs>>1;28         mul(x,x,x);29       end;30     z:=n;31     while (z<>p-1) and (y<>1) and (y<>p-1) do32       begin33         mul(y,y,y);34         z:=z<<1;35       end;36   exit((odd(z)) or (y=p-1));37   end;38 39 function isprime(p:int64):boolean;40  var i:longint;41  begin42   if (not(odd(p))) or (p=1) then exit(false);43   n:=p-1;m:=0;44   while n and 1=0 do n:=n>>1;45   for i:=1 to 10 do46    if not(check(random(p-2)+2)) then exit(false);47   exit(true);48  end;49 procedure main;50  begin51   if isprime(p) then writeln(‘Yes‘) else writeln(‘No‘);52  end;53 begin54   assign(input,‘input.txt‘);assign(output,‘output.txt‘);55   reset(input);rewrite(output);56   readln(t);57   while t>0 do58    begin59     dec(t);60     init;61     main;62    end;63   close(input);close(output);64 end.65       

View code