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 // 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
  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
  49                 private PrimalityTests ()
  50                 {
  51                 }
  52
  53                 #region SPP Test
  54
  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
  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>
 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
 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 }