The Montessori Algorithm for determining prime numbers of the ferma Theorem

Http://hi.baidu.com/xiaohanhoho/blog/item/f941e9948ed9f842d0135e72.html

Conventions:
X % Y is x modulo y, that is, the remainder of X divided by Y. When x <Y, X % Y = x, all the modulo operation objects are integers.
X ^ y indicates the Power Y of X.
The multiplication operator has a higher priority than multiplication, division, and modulo. addition and subtraction have the lowest priority.
When we see x ^ y/z, we calculate the multiplication side and then the division.
A/B is called dividing a by B. It is also called dividing a by B.
If a % B = 0, A can be divisible by B, or B can be divisible by.
A * B represents a multiplied by B or a multiplied by B, B multiplied by a, B multiplied by ...... It's the same as TMD!

Review Elementary School Mathematics
Public factor: two different natural numbers A and B. If there is a natural number C, which can be divisible by A or B, then C is the public factor of A and B.
Public multiple: two different natural numbers A and B. If a natural number C can be divisible by A or B, C is the public multiple of A and B.
Ry number: two different natural numbers. They have only one public factor of 1, and are called ry.

Ferma is a French mathematician and translated as "ferma". His profile is as follows. I don't know. I was shocked.
Http://www.cmr.com.cn/BasicStudy/LearnColumn/Maths/shuxuejiashi/j12.htm

Ferma's theorem:
If n is an arbitrary positive integer, P is a prime number, and N cannot be divisible by P (obviously N and P are mutually qualitative), then:
N ^ P % P = N (that is, the power P of n divided by the remainder of P is N)

However, I checked a lot of information and found the formula as follows:
(N ^ (P-1) % P = 1

After analysis, the two formulas are actually the same and can be deformed from each other. The original formula is:
(N ^ P-N) % P = 0 (that is, the p power of N minus N can be divisible by P, because The FeiMa theorem knows that the P power of N is divided by the remainder of P is N)

Put N to A, N ^ p becomes your n * (N ^ (P-1), then (N ^ P-N) % P = 0: (N * (N ^ (P-1)-1) % P = 0
Note the above formula, meaning: N * (N ^ (P-1)-1) can be divisible by P

And because N * (N ^ (P-1)-1) must be able to divide N (this does not bother !)
Therefore, N * (N ^ (P-1)-1) is a public multiple of N and P, primary knowledge ^_^

Because the premise is that N and P are mutually qualitative, and the minimum public multiple of the number is their product, so there must be a positive integer m to make the equation true:
N * (N ^ (P-1)-1) = M * n * P
About n on both sides, which is simplified:
N ^ (P-1)-1 = m * P
Because M is an integer, obviously:
(N ^ (P-1)-1) % P = 0
That is:
N ^ (P-1) % P = 1
========================================================== ====
Product modulus decomposition Formula

First, there is a theorem. If X % z = 0, that is, X can be divisible by Z, then:
(X + y) % z = Y % Z
You don't need to prove this...

If there are three positive integers x, y, and z, then there must be: (x * Y) % z = (X % Z) * (Y % z) % Z

It took a long time to prove it. We need to discuss the situation in detail:

1. When X and Y are larger than Z, there must be integers A and B to make the following equations true:
X = z * I + a (1)
Y = z * j + B (2)
Needless to say, this is the nature of the modulus operation!
(1) and (2) are substituted into (x * Y) modz: (z * I + a) (z * j + B) % Z
Multiply and then extract the first three items, and then convert them into: (z * I * j + I * A + I * B) + a * B) % Z (3)
Because Z * (z * I * j + I * A + I * B) is an integer multiple of Z ...... Dizzy, again.
(3) can be reduced to: (A * B) % Z
And because a = x % Z and B = Y % Z are substituted into the formula above, it becomes the original formula.

2. When X is larger than Z and Y is smaller than Z, the conversion is the same:
X = z * I +
(X * Y) % Z:
(Z * I * Y + A * Y) % Z
According to the theorem, the conversion is as follows: (A * Y) % Z
Because a = x % Z and Y = Y % Z are substituted into the above formula, the original formula is obtained.
Similarly, when X is smaller than Z and Y is larger than Z, the original formula is also true.

3. When X is smaller than Z and Y is also less than Z, x = x % Z, y = Y % Z, so the original formula is true.
========================================================== ==================
Fast Multiplication Algorithm

For example, to calculate 2 ^ 13, the traditional method requires 12 multiplications.

/* Calculate n ^ p */
Unsigned power (unsigned N, unsigned P)
{
For (INT I = 0; I <p; I ++) N * = N;
Return N;
}

Damn multiplication, it's time to optimize it! Save the result of 2*2 to see if it is 4*4*4*4*4*4*4*2.
Save the 4*4 Results: 16*16*16*2
A total of five operations are performed, which are 2*2, 4*4, and 16*16*16*2.

In this analysis, our algorithm needs to calculate a multiplication that is less than half of the total.
To illustrate this algorithm, let's take another example 2 ^ * 2*2*2*2*2*2*2*2.
Separated by two: (2*2) * (2*2) * (2*2) * 2
If we use 2*2 for calculation, then the index can be divided by 2. If there is no excess of one, we will multiply it later.
Separate them by two, and divide the index by 2: (2*2) * (2*2) * (2*2) * 2
In fact, the 2*2 in the last bracket is the rest of this time, so we will multiply it again later.
Now that the index is 1, you can calculate the final result: 16*4*2 = 128

The optimized algorithm is as follows:
Unsigned power (unsigned N, unsigned P)
{
Unsigned main = N; // use main to save the result
Unsigned odd = 1; // odd is used to calculate the "remaining" product.
While (p> 1)
{// Always calculate until the index is less than or equal to 1
If (P % 2 )! = 0)
{// If the exponent P is an odd number, it indicates that there will be an extra number after calculation. Here, we will multiply it into the result.
Odd * = Main; // multiply the "remaining"
}
Main * = Main; // subject multiplication party
P/= 2; // divide the index by 2
}
Return main * odd; // Finally, multiply the subject and the "remaining" as the result.
}

Is it perfect? No, not enough! See it? Main is unnecessary, and we can have faster code to judge the odd number. We need to know that division or modulo operation is very inefficient, so we can use an even number to optimize the code, that is, the limit bit in the even binary notation must be 0!

Perfect version:
Unsigned power (unsigned N, unsigned P)
{// Calculate the power of N to the power of P
Unsigned odd = 1; // oddk is used to calculate the product of "remaining"
While (p> 1)
{// Always calculate until the index is less than or equal to 1
If (P & 1 )! = 0)
{// Determine whether P is an odd number, and the even number must be 0
Odd * = N; // If the odd is an odd number, multiply the "remaining"
}
N * = N; // subject Multiplication
P/= 2; // divide the index by 2
}
Return N * odd; // Finally, multiply the subject and the "remaining" as the result.
}
========================================================== ======================
Montgomery "Fast Power modulo Algorithm

We will use this operation later: (x ^ y) % Z

The problem is that when X and Y are large, how can a 32-bit integer variable calculate the result effectively?
Considering the final optimization code above and the product modulo decomposition formula mentioned above, I think you may snap your head and say, "oh, I understand !".

The following explanation is intended for students who have not yet made such an action. X ^ y can be seen as multiplication of Y x, that is, there is a product modulo decomposition formula, so we can separate the process of multiplication of Y x and then modulo.For example: (17 ^ 25) % 29 can be divided into :( (17*17) % 29 * (17*17) % 29 *......
If we use the above Code to optimize this process, we will get the famous "Montgomery" Fast power mode algorithm:
Unsigned Montgomery (unsigned N, unsigned P, unsigned m)
{// Quickly calculate the value of (N ^ e) % m, which is very similar to the power algorithm.
Unsigned r = n % m; // The r here can not be saved
Unsigned k = 1;
While (p> 1)
{
If (P & 1 )! = 0)
{
K = (k * r) % m; // directly modulo
}
R = (R * r) % m; // same as above
P/= 2;
}
Return (R * k) % m; // same as above
}

The above code can also be optimized. The following is the speed version of Montgomery:

Unsigned Montgomery (unsigned N, unsigned P, unsigned m)
{// Quickly calculate the value of (N ^ e) % m
Unsignedk = 1;
N % = m;
While (P! = 1)
{
If (0! = (P & 1) k = (k * n) % m;
N = (N * n) % m;
P> = 1;
}
Return (N * k) % m;
}
========================================================== ==================
How can I determine whether a number is a prime number?

Dummies:
Bool isprime (unsigned N)
{
If (n <2)
{// The number smaller than 2 is neither a combination nor a prime number.
Throw 0;
}
For (unsigned I = 2; I <n; ++ I)
{// Returns the prime number if all vertices are less than it.
If (N % I = 0)
{// If the sum is exceeded, the sum is used.
Return false;
}
}
Return true;
}

Dividing a number by a number that is larger than half of a number must be different. Do you still need to judge it ??

The following are primary school students' practices:
Bool isprime (unsigned N)
{
If (n <2)
{// The number smaller than 2 is neither a combination nor a prime number.
Throw 0;
}
For (unsigned I = 2; I <n/2 + 1; ++ I)
{// Returns the prime number if they are not divided by any decimal number.
If (0 = n % I)
{// The sum is exceeded.
Return false;
}
}
Return true; // there is no division, prime number
}

A sum can be obtained by multiplying two or more prime numbers. If a number cannot be divisible by all the prime numbers smaller than half of it, it is impossible to divide all the total numbers smaller than half of it. It is useful to create a prime number table.

The following is the practice of middle school students:
Bool isprime2 (unsigned N)
{
If (n <2)
{// The number smaller than 2 is neither a combination nor a prime number.
Throw 0;
}
Static unsigned aprimelist [] = {// prime number table
1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41,
43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,113,
193,241,257,337,353,401,433,449,577,593,641,
673,769,881,929,977,100 9, 1153,120 1, 1217,124 9,
1297,1361, 1409,148 9, 1553,160 1, 1697,177 7, 1873,
1889,201 7, 2081,211 3, 2129,216 1, 2273,241 7, 2593,
2609,265 7, 2689,275 3, 2801,283 3, 2897,304 1, 3089,
3121,313 7, 3169,321 7, 3313,332 9, 3361,345 7, 3617,
3697,376 1, 3793,388 9, 4001,404 9, 4129,417 7, 4241,
4273,428 9, 4337,448 1, 4513,456 1, 4657,467 3, 4721,
4801,481 7, 4993,500 9, 5153,523 3, 5281,529 7, 5393,
5441,552 1, 5569,585 7, 5953,611 3, 6257,633 7, 6353,
6449,648 1, 6529,657 7, 6673,668 9, 6737,683 3, 6961,
6977,705 7, 7121,729 7, 7393,745 7, 7489,753 7, 7649,
7681,779 3, 7841,787 3, 7937,801 7, 8081,816 1, 8209,
8273,835 3, 8369,851 3, 8609,864 1, 8689,873 7, 8753,
8849,892 9, 9041,913 7, 9281,937 7, 9473,952 1, 9601,
9649,969 7, 9857
};

Const int nlistnum = sizeof (inclumelist)/sizeof (unsigned); // calculate the number of elements in the prime list
For (unsigned I = 2; I <nlistnum; ++ I)
{
If (n/2 + 1 <aprimelist [I])
{
Return true;
}
If (0 = n % else melist [I])
{
Return false;
}
}
/* Because the number of elements in the prime number table is limited, you can only use the dumb method for the numbers that cannot be determined by the prime number table */
For (unsigned I = systmelist [nListNum-1]; I <n/2 + 1; I ++)
{
If (0 = n % I)
{// The sum is exceeded.
Return false;
}
}
Return true;
}
It's still too bad. Now we need to judge the large prime number. The prime number table is used for the top p! Of course, we can use a dynamic prime number table for optimization, which is the practice of college students. However, dynamic EPO table strategies are complex and inefficient. Therefore, let's jump to the practice of experts directly:
According to the ferma theorem mentioned above, for the prime numbers N and P of the two mutually qualitative, there must be: N ^ (P-1) % P = 1
Let's use this property to determine the prime number. Of course, you will worry that when P is very large, the multiplication will be very troublesome. Don't worry! Didn't we have a fast idempotent algorithm above? Make good use of the joy that the big mathematician of Montgomery brings to us!

The algorithm IDEA is as follows:
For N, obtain any prime number from the prime number table and perform a trojan test on it. If n fails to be tested after many prime numbers are obtained, n is considered as the prime number. Of course, the more tests, the more accurate, but generally 50 tests are enough. In addition, we use the "primary school student" algorithm to construct an array containing 500 prime numbers. The Division test of Q will greatly improve the pass rate. The method is as follows:

Bool isprime3 (unsigned N)
{
If (n <2)
{// The number smaller than 2 is neither a combination nor a prime number.
Throw 0;
}
Static unsigned parameter melist [] = {
2, 3, 5, 7, 11, 17, 19, 23, 29, 31, 41,
43, 47, 53, 59, 67, 71, 73, 79, 83, 89, 97
};
Const int nlistnum = sizeof (inclumelist)/sizeof (unsigned );
For (INT I = 0; I <nlistnum; ++ I)
{// Judge the current prime number based on the number in the prime number table
If (1! = Montgomery (mongomelist [I], n-1, n) // Montgomery Algorithm
{
Return false;
}
}
Return true;
}

OK. This is the practice of experts.
Wait, what? It seems a bit strange. Let's take a look at this number of 29341. It is equivalent to 13*37*61. It is obviously a combination, but it has passed the test !! Oh, sorry, I forgot to add the numbers 13, 37, and 61 to the prime number table. I did it on purpose. I just want to explain that the ferma test is not completely reliable.
Now we have discovered an important point. The ferma theorem is a necessary condition for prime numbers, not a sufficient condition. This is not a prime number, but there are still a lot of numbers that can pass the ferma test. In mathematics, they are referred to as the carmic number. Now mathematicians have found all the Karl mcks within 10 ^ 16, the maximum value is 9585921133193329. We must find more effective testing methods. Mathematicians studied and extended the ferma's small theorem, and summed up a variety of fast and effective Prime Number testing methods. Currently, the fastest algorithm is rapinmiller's testing algorithm, the following describes the rapinmiller test.
========================================================== ======================================
Rapinmiller Test

The rapinmiller test is an uncertain algorithm. It can only be used to determine from the probability sense that a number may be a prime number, but it cannot be ensured. AlgorithmThe process is as follows:
1. Select t random number a, and a <n is true.
2. Locate R and m so that n = 2 * r * m + 1 is true.
Fast way to get R and M: N is represented by binary number B, so c = B-1. Because N is an odd number (the prime number is an odd number), the forward bit of C is 0. statistics are collected from the 0 of the forward bit of C until the first 1 is encountered. At this time, the number of 0 is R, and m is the value of B shifted right to R..
3. If a ^ m % n = 1, pass the test of a for N and then perform the test of the next.
4. If a ^ m % n! = 1, then let I iterate from 0 to R and perform the following test
5. If a ^ (2 ^ I) * m) % N = N-1 then pass a test on N; otherwise the test of the next I
6. If I = R and has not passed the test, the test of N for this a fails, indicating that N is the sum.
7. Perform the next A-to-n test until a specified number

Verification shows that when T is a prime number and A is a random number with an average distribution, the test efficiency is 1/(4 ^ t ). If T> 8, the probability of test errors will be less than 10 ^ (-5), which is sufficient for general applications. If the required prime number is extremely large or a higher degree of assurance is required, you can increase the T value appropriately. The following code is used:

Bool rabbinmillertest (unsigned N)
{
If (n <2)
{// The number smaller than 2 is neither a combination nor a prime number.
Throw 0;
}

Const unsigned nprimelistsize = sizeof (g_aprimelist)/sizeof (unsigned); // calculates the number of elements in the prime table.
For (INT I = 0; I <nprimelistsize; ++ I)
{// Judge the current prime number based on the number in the prime number table
If (n/2 + 1 <= g_w.melist [I])
{// If it is already smaller than the number of the current prime number table, it must be a prime number
Return true;
}
If (0 = n % g_1_melist [I])
{// If the remainder is 0, it indicates that it must not be a prime number.
Return false;
}
}
// Locate R and m so that n = 2 ^ r * m + 1;
Int r = 0, M = n-1; // (n-1) must be the sum
While (0 = (M & 1 ))
{
M> = 1; // shifts one digit to the right.
R ++; // count the number of shifts to the right
}
Const unsigned ntestcnt = 8; // indicates the number of tests performed.
For (unsigned I = 0; I <ntestcnt; ++ I)
{// Use a random number for testing,
Int A = g_1_melist [rand () % nprimelistsize];
If (1! = Montgomery (a, m, n ))
{
Int J = 0;
Int e = 1;
For (; j <r; ++ J)
{
If (n-1 = Montgomery (a, m * E, n ))
{
Break;
}
E <= 1;
}
If (j = r)
{
Return false;
}
}
}
Return true;
}