2010-06-21 Atsushi Enomoto <atsushi@ximian.com>
[mcs.git] / class / Mono.Security / Mono.Math.Prime / PrimalityTests.cs
blobcd3953e9955aa5c68d01637761d46b69d576e12e
1 //
2 // Mono.Math.Prime.PrimalityTests.cs - Test for primality
3 //
4 // Authors:
5 // Ben Maurer
6 //
7 // Copyright (c) 2003 Ben Maurer. All rights reserved
8 //
11 // Permission is hereby granted, free of charge, to any person obtaining
12 // a copy of this software and associated documentation files (the
13 // "Software"), to deal in the Software without restriction, including
14 // without limitation the rights to use, copy, modify, merge, publish,
15 // distribute, sublicense, and/or sell copies of the Software, and to
16 // permit persons to whom the Software is furnished to do so, subject to
17 // the following conditions:
18 //
19 // The above copyright notice and this permission notice shall be
20 // included in all copies or substantial portions of the Software.
21 //
22 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
23 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
24 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
25 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
26 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
27 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
28 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
31 using System;
33 namespace Mono.Math.Prime {
35 #if INSIDE_CORLIB
36 internal
37 #else
38 public
39 #endif
40 delegate bool PrimalityTest (BigInteger bi, ConfidenceFactor confidence);
42 #if INSIDE_CORLIB
43 internal
44 #else
45 public
46 #endif
47 sealed class PrimalityTests {
49 private PrimalityTests ()
53 #region SPP Test
55 private static int GetSPPRounds (BigInteger bi, ConfidenceFactor confidence)
57 int bc = bi.BitCount();
59 int Rounds;
61 // Data from HAC, 4.49
62 if (bc <= 100 ) Rounds = 27;
63 else if (bc <= 150 ) Rounds = 18;
64 else if (bc <= 200 ) Rounds = 15;
65 else if (bc <= 250 ) Rounds = 12;
66 else if (bc <= 300 ) Rounds = 9;
67 else if (bc <= 350 ) Rounds = 8;
68 else if (bc <= 400 ) Rounds = 7;
69 else if (bc <= 500 ) Rounds = 6;
70 else if (bc <= 600 ) Rounds = 5;
71 else if (bc <= 800 ) Rounds = 4;
72 else if (bc <= 1250) Rounds = 3;
73 else Rounds = 2;
75 switch (confidence) {
76 case ConfidenceFactor.ExtraLow:
77 Rounds >>= 2;
78 return Rounds != 0 ? Rounds : 1;
79 case ConfidenceFactor.Low:
80 Rounds >>= 1;
81 return Rounds != 0 ? Rounds : 1;
82 case ConfidenceFactor.Medium:
83 return Rounds;
84 case ConfidenceFactor.High:
85 return Rounds << 1;
86 case ConfidenceFactor.ExtraHigh:
87 return Rounds << 2;
88 case ConfidenceFactor.Provable:
89 throw new Exception ("The Rabin-Miller test can not be executed in a way such that its results are provable");
90 default:
91 throw new ArgumentOutOfRangeException ("confidence");
95 public static bool Test (BigInteger n, ConfidenceFactor confidence)
97 // Rabin-Miller fails with smaller primes (at least with our BigInteger code)
98 if (n.BitCount () < 33)
99 return SmallPrimeSppTest (n, confidence);
100 else
101 return RabinMillerTest (n, confidence);
104 /// <summary>
105 /// Probabilistic prime test based on Rabin-Miller's test
106 /// </summary>
107 /// <param name="n" type="BigInteger.BigInteger">
108 /// <para>
109 /// The number to test.
110 /// </para>
111 /// </param>
112 /// <param name="confidence" type="int">
113 /// <para>
114 /// The number of chosen bases. The test has at least a
115 /// 1/4^confidence chance of falsely returning True.
116 /// </para>
117 /// </param>
118 /// <returns>
119 /// <para>
120 /// True if "this" is a strong pseudoprime to randomly chosen bases.
121 /// </para>
122 /// <para>
123 /// False if "this" is definitely NOT prime.
124 /// </para>
125 /// </returns>
126 public static bool RabinMillerTest (BigInteger n, ConfidenceFactor confidence)
128 int bits = n.BitCount ();
129 int t = GetSPPRounds (bits, confidence);
131 // n - 1 == 2^s * r, r is odd
132 BigInteger n_minus_1 = n - 1;
133 int s = n_minus_1.LowestSetBit ();
134 BigInteger r = n_minus_1 >> s;
136 BigInteger.ModulusRing mr = new BigInteger.ModulusRing (n);
138 // Applying optimization from HAC section 4.50 (base == 2)
139 // not a really random base but an interesting (and speedy) one
140 BigInteger y = null;
141 // FIXME - optimization disable for small primes due to bug #81857
142 if (n.BitCount () > 100)
143 y = mr.Pow (2, r);
145 // still here ? start at round 1 (round 0 was a == 2)
146 for (int round = 0; round < t; round++) {
148 if ((round > 0) || (y == null)) {
149 BigInteger a = null;
151 // check for 2 <= a <= n - 2
152 // ...but we already did a == 2 previously as an optimization
153 do {
154 a = BigInteger.GenerateRandom (bits);
155 } while ((a <= 2) && (a >= n_minus_1));
157 y = mr.Pow (a, r);
160 if (y == 1)
161 continue;
163 for (int j = 0; ((j < s) && (y != n_minus_1)); j++) {
165 y = mr.Pow (y, 2);
166 if (y == 1)
167 return false;
170 if (y != n_minus_1)
171 return false;
173 return true;
176 public static bool SmallPrimeSppTest (BigInteger bi, ConfidenceFactor confidence)
178 int Rounds = GetSPPRounds (bi, confidence);
180 // calculate values of s and t
181 BigInteger p_sub1 = bi - 1;
182 int s = p_sub1.LowestSetBit ();
184 BigInteger t = p_sub1 >> s;
187 BigInteger.ModulusRing mr = new BigInteger.ModulusRing (bi);
189 for (int round = 0; round < Rounds; round++) {
191 BigInteger b = mr.Pow (BigInteger.smallPrimes [round], t);
193 if (b == 1) continue; // a^t mod p = 1
195 bool result = false;
196 for (int j = 0; j < s; j++) {
198 if (b == p_sub1) { // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
199 result = true;
200 break;
203 b = (b * b) % bi;
206 if (result == false)
207 return false;
209 return true;
212 #endregion
214 // TODO: Implement the Lucus test
215 // TODO: Implement other new primality tests
216 // TODO: Implement primality proving