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 // 9 10 // 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. 29 // 30 31 using System; 32 33 namespace Mono.Math.Prime { 34 35 #if INSIDE_CORLIB 36 internal 37 #else 38 public 39 #endif PrimalityTest(BigInteger bi, ConfidenceFactor confidence)40 delegate bool PrimalityTest (BigInteger bi, ConfidenceFactor confidence); 41 42 #if INSIDE_CORLIB 43 internal 44 #else 45 public 46 #endif 47 sealed class PrimalityTests { 48 PrimalityTests()49 private PrimalityTests () 50 { 51 } 52 53 #region SPP Test 54 GetSPPRounds(BigInteger bi, ConfidenceFactor confidence)55 private static int GetSPPRounds (BigInteger bi, ConfidenceFactor confidence) 56 { 57 int bc = bi.BitCount(); 58 59 int Rounds; 60 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; 74 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"); 92 } 93 } 94 Test(BigInteger n, ConfidenceFactor confidence)95 public static bool Test (BigInteger n, ConfidenceFactor confidence) 96 { 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); 102 } 103 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> RabinMillerTest(BigInteger n, ConfidenceFactor confidence)126 public static bool RabinMillerTest (BigInteger n, ConfidenceFactor confidence) 127 { 128 int bits = n.BitCount (); 129 int t = GetSPPRounds (bits, confidence); 130 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; 135 136 BigInteger.ModulusRing mr = new BigInteger.ModulusRing (n); 137 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); 144 145 // still here ? start at round 1 (round 0 was a == 2) 146 for (int round = 0; round < t; round++) { 147 148 if ((round > 0) || (y == null)) { 149 BigInteger a = null; 150 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)); 156 157 y = mr.Pow (a, r); 158 } 159 160 if (y == 1) 161 continue; 162 163 for (int j = 0; ((j < s) && (y != n_minus_1)); j++) { 164 165 y = mr.Pow (y, 2); 166 if (y == 1) 167 return false; 168 } 169 170 if (y != n_minus_1) 171 return false; 172 } 173 return true; 174 } 175 SmallPrimeSppTest(BigInteger bi, ConfidenceFactor confidence)176 public static bool SmallPrimeSppTest (BigInteger bi, ConfidenceFactor confidence) 177 { 178 int Rounds = GetSPPRounds (bi, confidence); 179 180 // calculate values of s and t 181 BigInteger p_sub1 = bi - 1; 182 int s = p_sub1.LowestSetBit (); 183 184 BigInteger t = p_sub1 >> s; 185 186 187 BigInteger.ModulusRing mr = new BigInteger.ModulusRing (bi); 188 189 for (int round = 0; round < Rounds; round++) { 190 191 BigInteger b = mr.Pow (BigInteger.smallPrimes [round], t); 192 193 if (b == 1) continue; // a^t mod p = 1 194 195 bool result = false; 196 for (int j = 0; j < s; j++) { 197 198 if (b == p_sub1) { // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1 199 result = true; 200 break; 201 } 202 203 b = (b * b) % bi; 204 } 205 206 if (result == false) 207 return false; 208 } 209 return true; 210 } 211 212 #endregion 213 214 // TODO: Implement the Lucus test 215 // TODO: Implement other new primality tests 216 // TODO: Implement primality proving 217 } 218 } 219