mcs/class/corlib/Mono.Math.Prime/PrimalityTests.cs

   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 // Copyright (C) 2004 Novell, Inc (http://www.novell.com)
  12 //
  13 // Permission is hereby granted, free of charge, to any person obtaining
  14 // a copy of this software and associated documentation files (the
  15 // "Software"), to deal in the Software without restriction, including
  16 // without limitation the rights to use, copy, modify, merge, publish,
  17 // distribute, sublicense, and/or sell copies of the Software, and to
  18 // permit persons to whom the Software is furnished to do so, subject to
  19 // the following conditions:
  20 //
  21 // The above copyright notice and this permission notice shall be
  22 // included in all copies or substantial portions of the Software.
  23 //
  24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  25 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  26 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  27 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  28 // LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  29 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  30 // WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  31 //
  32
  33 using System;
  34
  35 namespace Mono.Math.Prime {
  36
  37 #if INSIDE_CORLIB
  38         internal
  39 #else
  40         public
  41 #endif
  42         delegate bool PrimalityTest (BigInteger bi, ConfidenceFactor confidence);
  43
  44 #if INSIDE_CORLIB
  45         internal
  46 #else
  47         public
  48 #endif
  49         sealed class PrimalityTests {
  50
  51                 private PrimalityTests ()
  52                 {
  53                 }
  54
  55                 #region SPP Test
  56
  57                 private static int GetSPPRounds (BigInteger bi, ConfidenceFactor confidence)
  58                 {
  59                         int bc = bi.BitCount();
  60
  61                         int Rounds;
  62
  63                         // Data from HAC, 4.49
  64                         if      (bc <= 100 ) Rounds = 27;
  65                         else if (bc <= 150 ) Rounds = 18;
  66                         else if (bc <= 200 ) Rounds = 15;
  67                         else if (bc <= 250 ) Rounds = 12;
  68                         else if (bc <= 300 ) Rounds =  9;
  69                         else if (bc <= 350 ) Rounds =  8;
  70                         else if (bc <= 400 ) Rounds =  7;
  71                         else if (bc <= 500 ) Rounds =  6;
  72                         else if (bc <= 600 ) Rounds =  5;
  73                         else if (bc <= 800 ) Rounds =  4;
  74                         else if (bc <= 1250) Rounds =  3;
  75                         else                 Rounds =  2;
  76
  77                         switch (confidence) {
  78                                 case ConfidenceFactor.ExtraLow:
  79                                         Rounds >>= 2;
  80                                         return Rounds != 0 ? Rounds : 1;
  81                                 case ConfidenceFactor.Low:
  82                                         Rounds >>= 1;
  83                                         return Rounds != 0 ? Rounds : 1;
  84                                 case ConfidenceFactor.Medium:
  85                                         return Rounds;
  86                                 case ConfidenceFactor.High:
  87                                         return Rounds << 1;
  88                                 case ConfidenceFactor.ExtraHigh:
  89                                         return Rounds << 2;
  90                                 case ConfidenceFactor.Provable:
  91                                         throw new Exception ("The Rabin-Miller test can not be executed in a way such that its results are provable");
  92                                 default:
  93                                         throw new ArgumentOutOfRangeException ("confidence");
  94                         }
  95                 }
  96
  97                 /// <summary>
  98                 ///     Probabilistic prime test based on Rabin-Miller's test
  99                 /// </summary>
 100                 /// <param name="bi" type="BigInteger.BigInteger">
 101                 ///     <para>
 102                 ///         The number to test.
 103                 ///     </para>
 104                 /// </param>
 105                 /// <param name="confidence" type="int">
 106                 ///     <para>
 107                 ///     The number of chosen bases. The test has at least a
 108                 ///     1/4^confidence chance of falsely returning True.
 109                 ///     </para>
 110                 /// </param>
 111                 /// <returns>
 112                 ///     <para>
 113                 ///             True if "this" is a strong pseudoprime to randomly chosen bases.
 114                 ///     </para>
 115                 ///     <para>
 116                 ///             False if "this" is definitely NOT prime.
 117                 ///     </para>
 118                 /// </returns>
 119                 public static bool RabinMillerTest (BigInteger bi, ConfidenceFactor confidence)
 120                 {
 121                         int Rounds = GetSPPRounds (bi, confidence);
 122
 123                         // calculate values of s and t
 124                         BigInteger p_sub1 = bi - 1;
 125                         int s = p_sub1.LowestSetBit ();
 126
 127                         BigInteger t = p_sub1 >> s;
 128
 129                         int bits = bi.BitCount ();
 130                         BigInteger a = null;
 131                         BigInteger.ModulusRing mr = new BigInteger.ModulusRing (bi);
 132
 133                         // Applying optimization from HAC section 4.50 (base == 2)
 134                         // not a really random base but an interesting (and speedy) one
 135                         BigInteger b = mr.Pow (2, t);
 136                         if (b != 1) {
 137                                 bool result = false;
 138                                 for (int j=0; j < s; j++) {
 139                                         if (b == p_sub1) {         // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
 140                                                 result = true;
 141                                                 break;
 142                                         }
 143
 144                                         b = (b * b) % bi;
 145                                 }
 146                                 if (!result)
 147                                         return false;
 148                         }
 149
 150                         // still here ? start at round 1 (round 0 was a == 2)
 151                         for (int round = 1; round < Rounds; round++) {
 152                                 while (true) {                     // generate a < n
 153                                         a = BigInteger.GenerateRandom (bits);
 154
 155                                         // make sure "a" is not 0 (and not 2 as we have already tested that)
 156                                         if (a > 2 && a < bi)
 157                                                 break;
 158                                 }
 159
 160                                 if (a.GCD (bi) != 1)
 161                                         return false;
 162
 163                                 b = mr.Pow (a, t);
 164
 165                                 if (b == 1)
 166                                         continue;              // a^t mod p = 1
 167
 168                                 bool result = false;
 169                                 for (int j = 0; j < s; j++) {
 170
 171                                         if (b == p_sub1) {         // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
 172                                                 result = true;
 173                                                 break;
 174                                         }
 175
 176                                         b = (b * b) % bi;
 177                                 }
 178
 179                                 if (!result)
 180                                         return false;
 181                         }
 182                         return true;
 183                 }
 184
 185                 public static bool SmallPrimeSppTest (BigInteger bi, ConfidenceFactor confidence)
 186                 {
 187                         int Rounds = GetSPPRounds (bi, confidence);
 188
 189                         // calculate values of s and t
 190                         BigInteger p_sub1 = bi - 1;
 191                         int s = p_sub1.LowestSetBit ();
 192
 193                         BigInteger t = p_sub1 >> s;
 194
 195
 196                         BigInteger.ModulusRing mr = new BigInteger.ModulusRing (bi);
 197
 198                         for (int round = 0; round < Rounds; round++) {
 199
 200                                 BigInteger b = mr.Pow (BigInteger.smallPrimes [round], t);
 201
 202                                 if (b == 1) continue;              // a^t mod p = 1
 203
 204                                 bool result = false;
 205                                 for (int j = 0; j < s; j++) {
 206
 207                                         if (b == p_sub1) {         // a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
 208                                                 result = true;
 209                                                 break;
 210                                         }
 211
 212                                         b = (b * b) % bi;
 213                                 }
 214
 215                                 if (result == false)
 216                                         return false;
 217                         }
 218                         return true;
 219                 }
 220
 221                 #endregion
 222
 223                 // TODO: Implement the Lucus test
 224                 // TODO: Implement other new primality tests
 225                 // TODO: Implement primality proving
 226         }
 227 }