|
Hi,
Miguel suggested that I post this code to help with development on the JIT. It is a replacement for the current Mono.Math.BigInteger. The code will only work if things such as unsafe code and long emulation are correctly implemented, and its speed is directly related to the speed of long emulation.
This code finds the larger than a predefined number. It should output the number 4512701454269052790334614871457068366211465637830644056043775988879311685387672935588935849974585007504790954267419831655667405131168814386004100494057429, and then halt. Just press enter, and the program will quit.
The code takes ~ 3.5 seconds to execute using the new JIT on Miguel’s P4 1.8ghz, and 1.4 seconds on the old JIT with the same computer. It takes .46 seconds to execute on the MS JIT on a 1 ghz AMD.
The attached BigInteger code is not appropriate for production use, as it still has bugs that need to be worked out.
Sincerely,
|
using System;
namespace BigIntegerNewJitTest
{
/// <summary>
/// Summary description for Class1.
/// </summary>
class Class1
{
/// <summary>
/// The main entry point for the application.
/// </summary>
[STAThread]
static void Main(string[] args)
{
Console.WriteLine(BigInteger.BigInteger.NextHightestPrime(new
BigInteger.BigInteger(Start)));
Console.Read();
}
public static uint[] Start {
get {
return new uint[] {
0x5629a00f,
0x3b672419, 0x85891da8, 0x2d63ff0c, 0xd8b91375, 0xfb57f659,
0x07e2de17,
0x7de561be, 0xc29d6912, 0x198cb7fd, 0x251d33ad, 0x6bc20a0a,
0xdd1e4060,
0x809d5fb2, 0x20fcd816, 0xafc2ddb2
};
}
}
}
}
using System;
using System.Diagnostics;
using System.Security.Cryptography;
namespace BigInteger {
public class BigInteger {
#region Data Storage
/// <summary>
/// The Length of this BigInteger
/// </summary>
uint length = 1;
/// <summary>
/// The data for this BigInteger
/// </summary>
uint[] data;
// /// <summary>
// /// The Sign of this BigInteger
// /// </summary>
// Sign sign;
#endregion
#region Constants
/// <summary>
/// Default length of a BigInteger in bytes
/// </summary>
const uint DEFAULT_LEN = 20;
/// <summary>
/// Table of primes below 2000.
/// </summary>
/// <remarks>
/// <para>
/// This table was generated using Mathematica 4.1 using
the following function:
/// </para>
/// <para>
/// <code>
/// PrimeTable[x_] := Prime[Range[1, PrimePi[x]]]
/// PrimeTable[6000]
/// </code>
/// </para>
/// </remarks>
public static readonly uint[] smallPrimes = {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, 101, 103, 107, 109, 113, 127, 131, 137,
139, 149, 151,
157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223,
227, 229, 233,
239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307,
311, 313, 317,
331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397,
401, 409, 419,
421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487,
491, 499, 503,
509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593,
599, 601, 607,
613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677,
683, 691, 701,
709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787,
797, 809, 811,
821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883,
887, 907, 911,
919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997,
1009, 1013, 1019,
1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087,
1091, 1093, 1097,
1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,
1187, 1193, 1201,
1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
1283, 1289, 1291,
1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,
1381, 1399, 1409,
1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471,
1481, 1483, 1487,
1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559,
1567, 1571, 1579,
1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637,
1657, 1663, 1667,
1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747,
1753, 1759, 1777,
1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867,
1871, 1873, 1877,
1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973,
1979, 1987, 1993,
1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063,
2069, 2081, 2083,
2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143,
2153, 2161, 2179,
2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269,
2273, 2281, 2287,
2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
2371, 2377, 2381,
2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447,
2459, 2467, 2473,
2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579,
2591, 2593, 2609,
2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683,
2687, 2689, 2693,
2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753,
2767, 2777, 2789,
2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857,
2861, 2879, 2887,
2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969,
2971, 2999,
3001,
3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083,
3089, 3109, 3119,
3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209,
3217, 3221, 3229,
3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319,
3323, 3329, 3331,
3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
3433, 3449, 3457,
3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529,
3533, 3539, 3541,
3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617,
3623, 3631, 3637,
3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719,
3727, 3733, 3739,
3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833,
3847, 3851, 3853,
3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929,
3931, 3943, 3947,
3967, 3989,
//
4001, 4003, 4007, 4013, 4019, 4021,
4027, 4049, 4051, 4057, 4073,
//
4079, 4091, 4093, 4099, 4111, 4127,
4129, 4133, 4139, 4153, 4157, 4159, 4177,
//
4201, 4211, 4217, 4219, 4229, 4231,
4241, 4243, 4253, 4259, 4261, 4271, 4273,
//
4283, 4289, 4297, 4327, 4337, 4339,
4349, 4357, 4363, 4373, 4391, 4397, 4409,
//
4421, 4423, 4441, 4447, 4451, 4457,
4463, 4481, 4483, 4493, 4507, 4513, 4517,
//
4519, 4523, 4547, 4549, 4561, 4567,
4583, 4591, 4597, 4603, 4621, 4637, 4639,
//
4643, 4649, 4651, 4657, 4663, 4673,
4679, 4691, 4703, 4721, 4723, 4729, 4733,
//
4751, 4759, 4783, 4787, 4789, 4793,
4799, 4801, 4813, 4817, 4831, 4861, 4871,
//
4877, 4889, 4903, 4909, 4919, 4931,
4933, 4937, 4943, 4951, 4957, 4967, 4969,
//
4973, 4987, 4993, 4999,
//
5003, 5009, 5011, 5021, 5023, 5039,
5051, 5059, 5077,
//
5081, 5087, 5099, 5101, 5107, 5113,
5119, 5147, 5153, 5167, 5171, 5179, 5189,
//
5197, 5209, 5227, 5231, 5233, 5237,
5261, 5273, 5279, 5281, 5297, 5303, 5309,
//
5323, 5333, 5347, 5351, 5381, 5387,
5393, 5399, 5407, 5413, 5417, 5419, 5431,
//
5437, 5441, 5443, 5449, 5471, 5477,
5479, 5483, 5501, 5503, 5507, 5519, 5521,
//
5527, 5531, 5557, 5563, 5569, 5573,
5581, 5591, 5623, 5639, 5641, 5647, 5651,
//
5653, 5657, 5659, 5669, 5683, 5689,
5693, 5701, 5711, 5717, 5737, 5741, 5743,
//
5749, 5779, 5783, 5791, 5801, 5807,
5813, 5821, 5827, 5839, 5843, 5849, 5851,
//
5857, 5861, 5867, 5869, 5879, 5881,
5897, 5903, 5923, 5927, 5939, 5953, 5981,
//
5987
};
public enum Sign : int {
Negative = -1,
Zero = 0,
Positive = 1
};
#endregion
#region Constructors
public BigInteger () {
data = new uint[DEFAULT_LEN];
}
public BigInteger ( Sign sign, uint len ) {
this.data = new uint[len];
this.length = len;
}
public BigInteger( BigInteger bi ) {
this.data = (uint[])bi.data.Clone();
this.length = bi.length;
}
public BigInteger( BigInteger bi, uint len ) {
this.data = new uint[len];
for( uint i = 0; i < bi.length; i++ )
this.data[i] = bi.data[i];
this.length = bi.length;
}
#endregion
#region Conversions
public BigInteger ( uint ui ) {
data = new uint[2];
data[0] = ui;
}
public BigInteger ( int i ) {
data = new uint[1];
data[0] = (uint)(i > 0 ? i : -i);
// NEGPROBLEM
//sign = (Sign)Math.Sign(i);
}
public BigInteger ( ulong ul ) {
data = new uint[2];
length = 2;
data[0] = (uint)(ul & 0xFFFFFFFF);
data[1] = (uint)(ul >> 32);
this.Normalize();
}
public BigInteger(byte[] inData, Sign sign) {
length = (uint)inData.Length >> 2;
int leftOver = inData.Length & 0x3;
if(leftOver != 0) // length not multiples of 4
length++;
data = new uint[length];
for(int i = inData.Length - 1, j = 0; i >= 3; i -= 4, j++) {
data[j] = (uint)(
(inData[i-3] << 24) +
(inData[i-2] << 16) +
(inData[i-1] << 8) +
(inData[i ] << 0)
);
}
if(leftOver == 1)
data[length-1] = (uint)inData[0];
else if(leftOver == 2)
data[length-1] = (uint)((inData[0] << 8) + inData[1]);
else if(leftOver == 3)
data[length-1] = (uint)((inData[0] << 16) + (inData[1]
<< 8) + inData[2]);
this.Normalize();
}
public BigInteger(uint[] inData, Sign sign) {
length = (uint)inData.Length;
data = new uint[length];
for(int i = (int)length - 1, j = 0; i >= 0; i--, j++)
data[j] = inData[i];
this.Normalize();
}
public BigInteger(uint[] inData) {
length = (uint)inData.Length;
data = new uint[length];
for(int i = (int)length - 1, j = 0; i >= 0; i--, j++)
data[j] = inData[i];
this.Normalize();
}
public static implicit operator BigInteger (uint value) {
return (new BigInteger (value));
}
public static implicit operator BigInteger (int value) {
return (new BigInteger (value));
}
public static implicit operator BigInteger (ulong value) {
return (new BigInteger (value));
}
#endregion
#region Operators
public static BigInteger operator + (BigInteger bi1, BigInteger bi2) {
if (bi1 == 0) {
return new BigInteger(bi2);
}
else if (bi2 == 0) {
return new BigInteger(bi1);
}
else {
return Kernel.AddSameSign( bi1, bi2 );
}
}
public static BigInteger operator - (BigInteger bi1, BigInteger bi2) {
if (bi1 == 0) {
// NEGPROBLEM: what if bi2 != 0
return new BigInteger( bi2 );
}
else if (bi2 == 0) {
return new BigInteger(bi1);
}
else {
switch (Kernel.CompareSameSign(bi1, bi2)) {
case Sign.Zero:
return 0;
case Sign.Positive:
return Kernel.Subtract( bi1, bi2 );
case Sign.Negative:
// NEGPROBLEM
BigInteger ret = Kernel.Subtract( bi2,
bi1 );
//ret.sign = (Sign)(-1 *
(int)ret.sign);
return ret;
default:
throw new Exception();
}
}
}
public static int operator % (BigInteger bi, int i) {
if ( i > 0 )
return (int)Kernel.DwordMod(bi, (uint)i);
else
return -(int)Kernel.DwordMod(bi, (uint)-i);
}
public static uint operator % (BigInteger bi, uint ui) {
return Kernel.DwordMod(bi, (uint)ui);
}
public static BigInteger operator % ( BigInteger bi1, BigInteger bi2) {
return Kernel.multiByteDivide(bi1, bi2)[1];
}
public static BigInteger operator / (BigInteger bi, int i) {
BigInteger ret;
if ( i > 0 ) {
ret = Kernel.DwordDiv( bi, (uint)i );
}
else {
ret = Kernel.DwordDiv( bi, (uint)-i );
// NEGPROBLEM
//ret.sign = (Sign)(-(uint)ret.sign);
}
return ret;
}
public static BigInteger operator / ( BigInteger bi1, BigInteger bi2) {
return Kernel.multiByteDivide(bi1, bi2)[0];
}
public static BigInteger operator * (BigInteger bi1, BigInteger bi2) {
//
// Validate pointers
//
if ( bi1.data.Length < bi1.length ) throw new
IndexOutOfRangeException( "bi1 out of range" );
if ( bi2.data.Length < bi2.length ) throw new
IndexOutOfRangeException( "bi2 out of range" );
BigInteger ret = new BigInteger( Sign.Positive, bi1.length +
bi2.length );
if (bi1 == 0 || bi2 == 0) return ret;
Kernel.Multiply(bi1.data, 0, bi1.length, bi2.data, 0,
bi2.length, ret.data, 0);
ret.Normalize();
return ret;
}
public static BigInteger operator * (BigInteger bi, int i) {
if ( i == 0 ) return 0;
if ( i == 1 ) return new BigInteger (bi);
uint ii = i > 0 ? (uint)i : (uint)-i;
Sign s = (Sign)(i / ii);
return Kernel.MultiplyByDword( bi, ii, s );
}
// public static BigInteger operator - ( BigInteger bi ) {
// BigInteger ret = new BigInteger(bi);
// ret.sign = (Sign)((int)bi.sign * -1);
// return ret;
// }
public static BigInteger operator <<(BigInteger bi1, int shiftVal) {
return Kernel.LeftShift(bi1, shiftVal);
}
public static BigInteger operator >>(BigInteger bi1, int shiftVal) {
return Kernel.RightShift(bi1, shiftVal);
}
#endregion
#region Random
private static RandomNumberGenerator rng;
private static RandomNumberGenerator Rng {
get {
if (rng == null)
rng = new RNGCryptoServiceProvider();
return rng;
}
}
/// <summary>
/// Generates a new, random BigInteger of the specified length.
/// </summary>
/// <param name="bits">The number of bits for the new number.</param>
/// <param name="rng">A random number generator to use to obtain the
bits.</param>
/// <returns>A random number of the specified length.</returns>
public static BigInteger genRandom ( int bits, RandomNumberGenerator
rng ) {
int dwords = bits >> 5;
int remBits = bits & 0x1F;
if (remBits != 0)
dwords++;
BigInteger ret = new BigInteger( Sign.Positive, (uint)dwords );
byte[] random = new byte [dwords << 2];
rng.GetBytes( random );
Buffer.BlockCopy(random, 0, ret.data, 0, (int)dwords << 2);
if (remBits != 0) {
uint mask = (uint)(0x01 << (remBits-1));
ret.data[dwords-1] |= mask;
mask = (uint)(0xFFFFFFFF >> (32 - remBits));
ret.data[dwords-1] &= mask;
}
else
ret.data[dwords-1] |= 0x80000000;
ret.Normalize();
return ret;
}
/// <summary>
/// Generates a new, random BigInteger of the specified length using
the default RNG crypto service provider.
/// </summary>
/// <param name="bits">The number of bits for the new number.</param>
/// <returns>A random number of the specified length.</returns>
public static BigInteger genRandom ( int bits ) {
return genRandom( bits, Rng);
}
/// <summary>
/// Randomizes the bits in "this" from the specified RNG.
/// </summary>
/// <param name="rng">A RNG.</param>
public void randomize ( RandomNumberGenerator rng ) {
int bits = (int)length;
int dwords = bits >> 5;
int remBits = bits & 0x1F;
if (remBits != 0)
dwords++;
byte[] random = new byte [dwords << 2];
rng.GetBytes( random );
Buffer.BlockCopy(random, 0, data, 0, (int)dwords << 2);
if (remBits != 0) {
uint mask = (uint)(0x01 << (remBits-1));
data[dwords-1] |= mask;
mask = (uint)(0xFFFFFFFF >> (32 - remBits));
data[dwords-1] &= mask;
}
else
data[dwords-1] |= 0x80000000;
Normalize();
}
/// <summary>
/// Randomizes the bits in "this" from the default RNG.
/// </summary>
public void randomize () {
randomize( Rng );
}
#endregion
#region Bitwise
public int bitCount() {
this.Normalize();
uint value = data[length - 1];
uint mask = 0x80000000;
uint bits = 32;
while(bits > 0 && (value & mask) == 0) {
bits--;
mask >>= 1;
}
bits += ((length - 1) << 5);
return (int)bits;
}
/// <summary>
/// Tests if the specified bit is 1.
/// </summary>
/// <param name="bitNum">The bit to test. The least significant bit is
0.</param>
/// <returns>True if bitNum is set to 1, else false.</returns>
public bool testBit (uint bitNum) {
uint bytePos = bitNum >> 5; // divide by 32
byte bitPos = (byte)(bitNum & 0x1F); // get the lowest 5
bits
uint mask = (uint)1 << bitPos;
return ((this.data[bytePos] | mask) == this.data[bytePos]);
}
#endregion
#region Compare
public static bool operator == ( BigInteger bi1, uint ui ) {
if( bi1.length != 1 ) bi1.Normalize();
return bi1.length == 1 && bi1.data[0] == ui;
}
public static bool operator != ( BigInteger bi1, uint ui ) {
if( bi1.length != 1 ) bi1.Normalize();
return !(bi1.length == 1 && bi1.data[0] == ui);
}
public static bool operator == (BigInteger bi1, BigInteger bi2) {
return Kernel.Compare(bi1, bi2) == 0;
}
public static bool operator !=( BigInteger bi1, BigInteger bi2) {
return Kernel.Compare(bi1, bi2) != 0;
}
public static bool operator > (BigInteger bi1, BigInteger bi2) {
return Kernel.Compare(bi1, bi2) > 0;
}
public static bool operator < (BigInteger bi1, BigInteger bi2) {
return Kernel.Compare(bi1, bi2) < 0;
}
public static bool operator >= (BigInteger bi1, BigInteger bi2) {
return Kernel.Compare(bi1, bi2) >= 0;
}
public static bool operator <= (BigInteger bi1, BigInteger bi2) {
return Kernel.Compare(bi1, bi2) <= 0;
}
public Sign Compare (BigInteger bi) {
return Kernel.Compare( this, bi );
}
#endregion
#region Formatting
public string ToString (uint radix) {
if(radix < 2 || radix > 36)
throw (new ArgumentException("Radix must be >= 2 and
<= 36"));
return ToString( radix, "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"
);
}
public string ToString( uint radix, string charSet ) {
if ( charSet.Length < radix )
throw new ArgumentException("charSet length less than
radix", "charSet");
string result = "";
BigInteger a = new BigInteger(this);
uint rem = 0;
if(a == 0)
result = "0";
else {
while(a.length > 1 || (a.length == 1 && a.data[0] !=
0)) {
rem = Kernel.SingleByteDivideInPlace( a, radix
);
result = charSet[(int)rem] + result;
}
}
return result;
}
#endregion
#region Misc
/// <summary>
/// Normalizes this by setting the length to the actual number of
/// uints used in data and by setting the sign to Sign.Zero if the
/// value of this is 0.
/// </summary>
private void Normalize() {
// Normalize length
while(length > 0 && data[length-1] == 0) length--;
// Check for zero
if (length == 0) {
length++;
}
}
#endregion
#region Object Impl
public override int GetHashCode () {
uint val = 0;
for (uint i = 0; i < this.length; i++) {
val ^= this.data[i];
}
unchecked {
return (int)val;
}
}
public override string ToString() {
return ToString( 10 );
}
public override bool Equals (object o) {
if ( o == null ) return false;
if ( o is int) return this == (int)o;
return Kernel.Compare(this, (BigInteger)o) == 0;
}
#endregion
#region Number Theory
public BigInteger gcd( BigInteger bi ) {
return Kernel.gcd( this, bi );
}
public BigInteger modInverse( BigInteger mod ) {
return Kernel.modInverse(this, mod);
}
#endregion
#region Prime Testing
/// <summary>
/// Returns the number of Strong Pseudo Prime tests that must be
performed
/// to make the chance of failing less than 2^-80
/// </summary>
/// <param name="bc">Number of bits in the number</param>
/// <returns>The number of rounds of the SPP test that should be
done</returns>
private static int SppItrsForBitCount( uint bc ) {
// Data from HAC, 4.49
if ( bc <= 100 ) return 27;
else if ( bc <= 150 ) return 18;
else if ( bc <= 200 ) return 15;
else if ( bc <= 250 ) return 12;
else if ( bc <= 300 ) return 9;
else if ( bc <= 350 ) return 8;
else if ( bc <= 400 ) return 7;
else if ( bc <= 500 ) return 6;
else if ( bc <= 600 ) return 5;
else if ( bc <= 800 ) return 4;
else if ( bc <= 1250 ) return 3;
else return 2;
}
public bool isProbablePrime() {
BigInteger thisVal = this;
// test for divisibility by primes < 2000
for(int p = 0; p < smallPrimes.Length; p++) {
if(thisVal % smallPrimes[p] == 0)
return false;
}
return
Kernel.SmallPrimeSppTest(thisVal,
SppItrsForBitCount(thisVal.length * 32));
}
#endregion
#region Prime Number Generation
/// <summary>
/// Generates the smallest prime >= bi
/// </summary>
/// <param name="bi">A BigInteger</param>
/// <returns>The smallest prime >= bi. More mathematically, if bi is
prime: bi, else Prime[PrimePi[bi] + 1].</returns>
public static BigInteger NextHightestPrime( BigInteger bi ) {
// We copy the value to a new variable, because we will be
// modifying the bits. The length is increased incase there is
not
// an prime >= bi with the same length as bi (very unlikely,
but
// not worth the risk)
BigInteger curVal = new BigInteger( bi, bi.length + 1 );
// All primes > 2 are odd.
curVal.data[0] |= 0x1;
const uint primeProd1 = 3u* 5u * 7u * 11u * 13u * 17u * 19u *
23u * 29u;
uint pMod1 = curVal % primeProd1;
while(true) {
if ( pMod1 % 3 == 0 ) goto biNotPrime;
if ( pMod1 % 5 == 0 ) goto biNotPrime;
if ( pMod1 % 7 == 0 ) goto biNotPrime;
if ( pMod1 % 11 == 0 ) goto biNotPrime;
if ( pMod1 % 13 == 0 ) goto biNotPrime;
if ( pMod1 % 17 == 0 ) goto biNotPrime;
if ( pMod1 % 19 == 0 ) goto biNotPrime;
if ( pMod1 % 23 == 0 ) goto biNotPrime;
if ( pMod1 % 29 == 0 ) goto biNotPrime;
for(int p = 9; p < smallPrimes.Length; p++) {
if(curVal % smallPrimes[p] == 0)
goto biNotPrime;
}
if ( Kernel.SmallPrimeSppTest(curVal,
SppItrsForBitCount(curVal.length * 32)) ) return curVal;
biNotPrime:
pMod1 += 2;
if (pMod1 >= primeProd1) pMod1 -= primeProd1;
curVal.Incr2();
}
}
/// <summary>
/// Increments this by two
/// </summary>
private void Incr2() {
unchecked {
int i = 0;
data[0] += 2;
// If there was no carry, nothing to do
if (data[0] < 2) {
// Account for the first carry
data[++i]++;
// Keep adding until no carry
while(data[i++] == 0x0)
data[i]++;
// See if we increased the data length
if( length == (uint)i )
length++;
}
}
}
#endregion
public sealed class ModulusRing {
BigInteger mod, constant;
public ModulusRing( BigInteger mod ) {
this.mod = mod;
// calculate constant = b^(2k) / m
uint i = mod.length << 1;
constant = new BigInteger( Sign.Positive, i + 1 );
constant.data[i] = 0x00000001;
constant = constant / mod;
}
public BigInteger BarrettReduction(BigInteger x) {
BigInteger n = mod;
uint k = n.length,
kPlusOne = k+1,
kMinusOne = k-1;
BigInteger q3;
//
// Validate pointers
//
if ( x.data.Length < x.length ) throw new
IndexOutOfRangeException( "x out of range" );
if( x.length < kPlusOne )
return x;
// q1 = x / b^(k-1)
// q2 = q1 * constant
// q3 = q2 / b^(k+1), Needs to be accessed with an
offset of kPlusOne
// TODO: We should the method in HAC p 604 to do this
(14.45)
q3 = new BigInteger(Sign.Positive, x.length -
kMinusOne + constant.length);
Kernel.Multiply(x.data, kMinusOne, x.length -
kMinusOne, constant.data, 0, constant.length, q3.data, 0);
// r1 = x mod b^(k+1)
// i.e. keep the lowest (k+1) words
uint lengthToCopy = (x.length > kPlusOne) ? kPlusOne :
x.length;
x.length = lengthToCopy;
x.Normalize();
// r2 = (q3 * n) mod b^(k+1)
// partial multiplication of q3 and n
BigInteger r2 = new BigInteger(Sign.Positive,
kPlusOne);
Kernel.MultiplyMod2p32pmod(q3.data, (int)kPlusOne,
(int)q3.length - (int)kPlusOne, n.data, 0, (int)n.length, r2.data, 0, (int)kPlusOne);
r2.Normalize();
if (r2 < x ) {
Kernel.MinusEq(x, r2);
}
else {
BigInteger val = new BigInteger(Sign.Positive,
kPlusOne + 1);
val.data[kPlusOne] = 0x00000001;
BigInteger aa = x - r2;
BigInteger bb = aa + val;
Kernel.MinusEq( val, r2 );
Kernel.PlusEq( x, val );
}
while(x >= n)
Kernel.MinusEq(x, n);
return x;
}
public BigInteger Multiply( BigInteger a, BigInteger b ) {
// XXX: is this right?
if (a.length >= mod.length << 1)
a %= mod;
if (b.length >= mod.length << 1)
b %= mod;
if (a.length >= mod.length)
a = BarrettReduction( a );
if (b.length >= mod.length)
b = BarrettReduction( b );
return BarrettReduction( a * b);
}
public BigInteger Pow(BigInteger b, BigInteger exp) {
BigInteger resultNum = new BigInteger( (BigInteger)1,
mod.length << 1 );
BigInteger tempNum;
tempNum = b % mod; // ensures (tempNum * tempNum) <
b^(2k)
tempNum = new BigInteger( tempNum, mod.length << 1 );
int totalBits = exp.bitCount();
int count = 0;
uint[] wkspace = new uint[mod.length << 1];
// perform squaring and multiply exponentiation
for(int pos = 0; pos < exp.length; pos++) {
uint mask = 0x01;
for(int index = 0; index < 32; index++) {
if((exp.data[pos] & mask) != 0) {
Array.Clear(wkspace, 0,
wkspace.Length);
Kernel.Multiply(resultNum.data, 0, resultNum.length, tempNum.data, 0, tempNum.length,
wkspace, 0);
resultNum.length +=
tempNum.length;
uint[] t = wkspace;
wkspace = resultNum.data;
resultNum.data = t;
resultNum =
BarrettReduction(resultNum);
}
mask <<= 1;
Kernel.SquarePositive(tempNum, ref
wkspace);
tempNum = BarrettReduction(tempNum);
if(tempNum == 1) {
return resultNum;
}
count++;
if(count == totalBits)
break;
}
}
return resultNum;
}
#region Two Pow
/// <summary>
/// Calculates 2^exp (mod mod)
/// </summary>
/// <param name="exp">The exponent</param>
/// <returns>2^exp (mod mod)</returns>
///
public BigInteger TwoPow(BigInteger exp) {
if((mod.data[0] & 1) == 1) return OddModTwoPow(exp);
// TODO: make tests for this
else return EvenModTwoPow(exp);
}
private unsafe BigInteger EvenModTwoPow(BigInteger exp) {
exp.Normalize();
uint[] wkspace = new uint[mod.length << 1 + 1];
BigInteger resultNum = new BigInteger( 2, mod.length
<< 1 +1);
uint value = exp.data[exp.length - 1];
uint mask = 0x80000000;
// Find the first bit of the exponent
while((value & mask) == 0)
mask >>= 1;
//
// We know that the first itr will make the val 2,
// so eat one bit of the exponent
//
mask >>= 1;
uint wPos = exp.length - 1;
do {
value = exp.data[wPos];
do {
Kernel.SquarePositive(resultNum, ref
wkspace);
if( !(resultNum.length < mod.length) )
resultNum =
BarrettReduction(resultNum);
if( (value & mask) != 0 ) {
//
// resultNum = (resultNum * 2)
% mod
//
fixed ( uint* u =
resultNum.data) {
//
// Double
//
uint* uu = u;
uint* uuE = u +
resultNum.length;
uint x, carry = 0;
while (uu < uuE) {
x = *uu;
*uu = (x << 1)
| carry;
carry = x >>
(32 - 1);
uu++;
}
// subtraction inlined
because we know it is square
if( carry != 0 ||
resultNum >= mod ) {
uu = u;
uint c = 0;
uint[] s =
mod.data;
uint i = 0;
do {
uint a
= s[i];
if (
((a += c) < c) | ( (*(uu++) -= a) > ~a ) )
c = 1;
else
c = 0;
i++;
} while ( uu <
uuE );
}
}
}
} while ( (mask >>= 1) > 0 );
mask = 0x80000000;
} while( wPos-- > 0 );
return resultNum;
}
private unsafe BigInteger OddModTwoPow(BigInteger exp) {
uint[] wkspace = new uint[mod.length << 1 + 1];
BigInteger resultNum =
Montgomery.ToMont((BigInteger)2, this.mod);
resultNum = new BigInteger( resultNum, mod.length << 1
+1);
uint mPrime = Montgomery.Inverse ( mod.data[0] );
//
// TODO: eat small bits, the ones we can do with no
modular reduction
//
uint pos = (uint)exp.bitCount() - 2;
do {
Kernel.SquarePositive(resultNum, ref wkspace);
resultNum = Montgomery.Reduce( resultNum, mod,
mPrime);
if( exp.testBit(pos) ) {
//
// resultNum = (resultNum * 2) % mod
//
fixed ( uint* u = resultNum.data) {
//
// Double
//
uint* uu = u;
uint* uuE = u +
resultNum.length;
uint x, carry = 0;
while (uu < uuE) {
x = *uu;
*uu = (x << 1) | carry;
carry = x >> (32 - 1);
uu++;
}
// subtraction inlined because
we know it is square
if( carry != 0 || resultNum
>= mod ) {
fixed( uint* s =
mod.data ) {
uu = u;
uint c = 0;
uint* ss = s;
do {
uint a
= *ss++;
if (
((a += c) < c) | ( (*(uu++) -= a) > ~a ) )
c = 1;
else
c = 0;
} while ( uu <
uuE );
}
}
}
}
} while ( pos-- > 0 );
resultNum = Montgomery.Reduce( resultNum, mod, mPrime);
return resultNum;
}
#endregion
#region Pow Small Base
public BigInteger Pow(uint b, BigInteger exp) {
if ( b != 2 ) {
if((mod.data[0] & 1) == 1) return OddPow(b,
exp);
// TODO: make tests for this
else return EvenPow(b, exp);
}
else {
if((mod.data[0] & 1) == 1) return
OddModTwoPow(exp);
// TODO: make tests for this
else return EvenModTwoPow(exp);
}
}
private unsafe BigInteger OddPow(uint b, BigInteger exp) {
exp.Normalize();
uint[] wkspace = new uint[mod.length << 1 + 1];
BigInteger resultNum =
Montgomery.ToMont((BigInteger)b, this.mod);
resultNum = new BigInteger( resultNum, mod.length << 1
+1);
uint mPrime = Montgomery.Inverse ( mod.data[0] );
uint value = exp.data[exp.length - 1];
uint mask = 0x80000000;
while((value & mask) == 0)
mask >>= 1;
//
// We know that the first itr will make the val b
//
mask >>= 1;
uint wPos = exp.length - 1;
do {
value = exp.data[wPos];
do {
//
// r = r ^ 2 % m
//
Kernel.SquarePositive(resultNum, ref
wkspace);
resultNum = Montgomery.Reduce(
resultNum, mod, mPrime);
if( (value & mask) != 0 ) {
//
// r = r * b % m
//
// TODO: Is Unsafe really
speeding things up?
fixed ( uint* u =
resultNum.data) {
uint i = 0;
ulong mc = 0;
do {
mc +=
(ulong)u[i] * (ulong)b;
u[i] =
(uint)(mc & 0xFFFFFFFF);
mc >>= 32;
} while (++i <
resultNum.length);
if(resultNum.length <
mod.length){
if (mc != 0) {
u[i] =
(uint)mc;
resultNum.length++;
while(resultNum >= mod)
Kernel.MinusEq( resultNum, mod);
}
}
else if( mc != 0 ) {
//
// First, we
estimate the quotient by dividing
// the first
part of each of the numbers. Then
// we correct
this, if necessary, with a subtraction.
//
uint cc =
(uint)mc;
// We would
rather have this estimate overshoot,
// so we add
one to the divisor
uint
divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u[i -1]) /
(mod.data[mod.length-1] + 1));
uint t;
i = 0;
mc = 0;
do {
mc +=
(ulong)mod.data[i] * (ulong)divEstimate;
t =
u[i];
u[i]
-= (uint)mc;
mc >>=
32;
if(
u[i] > t) mc++;
i++;
} while ( i <
resultNum.length );
cc -= (uint)mc;
if( cc != 0 ) {
uint
sc = 0, j = 0;
uint[]
s = mod.data;
do {
uint a = s[j];
if ( ((a += sc) < sc) | ( (u[j] -= a) > ~a ) ) sc = 1;
else sc = 0;
j++;
}
while ( j < resultNum.length );
cc -=
sc;
}
while(resultNum >= mod)
Kernel.MinusEq( resultNum, mod);
}
else {
while(resultNum >= mod)
Kernel.MinusEq( resultNum, mod);
}
}
}
} while ( (mask >>= 1) > 0 );
mask = 0x80000000;
} while( wPos-- > 0 );
resultNum = Montgomery.Reduce( resultNum, mod, mPrime);
return resultNum;
}
private unsafe BigInteger EvenPow(uint b, BigInteger exp) {
exp.Normalize();
uint[] wkspace = new uint[mod.length << 1 + 1];
BigInteger resultNum = new BigInteger( (BigInteger)b,
mod.length << 1 + 1);
uint value = exp.data[exp.length - 1];
uint mask = 0x80000000;
while((value & mask) == 0)
mask >>= 1;
//
// We know that the first itr will make the val b
//
mask >>= 1;
uint wPos = exp.length - 1;
do {
value = exp.data[wPos];
do {
//
// r = r ^ 2 % m
//
Kernel.SquarePositive(resultNum, ref
wkspace);
if( !(resultNum.length < mod.length) )
resultNum =
BarrettReduction(resultNum);
if( (value & mask) != 0 ) {
//
// r = r * b % m
//
// TODO: Is Unsafe really
speeding things up?
fixed ( uint* u =
resultNum.data) {
uint i = 0;
ulong mc = 0;
do {
mc +=
(ulong)u[i] * (ulong)b;
u[i] =
(uint)(mc & 0xFFFFFFFF);
mc >>= 32;
} while (++i <
resultNum.length);
if(resultNum.length <
mod.length){
if (mc != 0) {
u[i] =
(uint)mc;
resultNum.length++;
while(resultNum >= mod)
Kernel.MinusEq( resultNum, mod);
}
}
else if( mc != 0 ) {
//
// First, we
estimate the quotient by dividing
// the first
part of each of the numbers. Then
// we correct
this, if necessary, with a subtraction.
//
uint cc =
(uint)mc;
// We would
rather have this estimate overshoot,
// so we add
one to the divisor
uint
divEstimate = (uint) ((((ulong)cc << 32) | (ulong) u[i -1]) /
(mod.data[mod.length-1] + 1));
uint t;
i = 0;
mc = 0;
do {
mc +=
(ulong)mod.data[i] * (ulong)divEstimate;
t =
u[i];
u[i]
-= (uint)mc;
mc >>=
32;
if(
u[i] > t) mc++;
i++;
} while ( i <
resultNum.length );
cc -= (uint)mc;
if( cc != 0 ) {
uint
sc = 0, j = 0;
uint[]
s = mod.data;
do {
uint a = s[j];
if ( ((a += sc) < sc) | ( (u[j] -= a) > ~a ) ) sc = 1;
else sc = 0;
j++;
}
while ( j < resultNum.length );
cc -=
sc;
}
while(resultNum >= mod)
Kernel.MinusEq( resultNum, mod);
}
else {
while(resultNum >= mod)
Kernel.MinusEq( resultNum, mod);
}
}
}
} while ( (mask >>= 1) > 0 );
mask = 0x80000000;
} while( wPos-- > 0 );
return resultNum;
}
#endregion
}
public sealed class Montgomery {
public static uint Inverse( uint n ) {
uint y = n, z;
Debug.Assert( (n & 1) == 1, "n Must be odd");
while ((z = n * y) != 1)
y *= 2 - z;
return (uint)-y;
}
public static BigInteger ToMont( BigInteger n, BigInteger m ) {
n.Normalize(); m.Normalize();
n <<= (int)m.length * 32;
n %= m;
return n;
}
public static unsafe BigInteger Reduce( BigInteger n,
BigInteger m, uint mPrime ) {
BigInteger A = n;
fixed( uint* a = A.data, mm = m.data ) {
for( uint i = 0; i < m.length; i++ ) {
// The mod here is taken care of by
the CPU,
// since the multiply will overflow.
uint u_i = a[0] * mPrime /* % 2^32 */;
//
// A += u_i * m;
// A >>= 32
//
// mP = Position in mod
// aSP = the source of bits from a
// aDP = destination for bits
uint* mP = mm, aSP = a, aDP = a;
ulong c = (ulong)u_i * (ulong)*(mP++)
+ *(aSP++);
c >>= 32;
uint j = 1;
// Multiply and add
for( ; j < m.length; j++ ) {
c += (ulong)u_i *
(ulong)*(mP++) + *(aSP++);
*(aDP++) = (uint)c;
c >>= 32;
}
// Account for carry
// TODO: use a better loop here, we
dont need the ulong stuff
for( ; j < A.length; j++ ) {
c += *(aSP++);
*(aDP++) = (uint)c;
c >>= 32;
if( c == 0) {j++; break;}
}
// Copy the rest
for( ; j < A.length; j++ ) {
*(aDP++) = *(aSP++);
}
*(aDP++) = (uint)c;
}
while( A.length > 1 && a[A.length-1] == 0 )
A.length--;
}
if ( A >= m ) Kernel.MinusEq( A, m );
return A;
}
public static BigInteger Reduce( BigInteger n, BigInteger m ) {
return Reduce(n, m, Inverse(m.data[0]));
}
}
/// <summary>
/// Low level functions for the BigInteger
/// </summary>
private sealed class Kernel {
#region Addition/Subtraction
/// <summary>
/// Adds two numbers with the same sign.
/// </summary>
/// <param name="bi1">A BigInteger</param>
/// <param name="bi2">A BigInteger</param>
/// <returns>bi1 + bi2</returns>
public static BigInteger AddSameSign ( BigInteger bi1,
BigInteger bi2 ) {
uint[] x, y;
uint yMax, xMax, i = 0;
// x should be bigger
if (bi1.length < bi2.length){
x = bi2.data;
xMax = bi2.length;
y = bi1.data;
yMax = bi1.length;
}
else {
x = bi1.data;
xMax = bi1.length;
y = bi2.data;
yMax = bi2.length;
}
BigInteger result = new BigInteger( Sign.Positive,
xMax + 1 );
uint[] r = result.data;
ulong sum = 0;
// Add common parts of both numbers
do {
sum = ((ulong)x[i]) +
((ulong)y[i]) + sum;
r[i] = (uint)sum;
sum >>= 32;
} while ( ++i < yMax );
// Copy remainder of longer number while carry
propagation is required
bool carry = (sum != 0);
if (carry){
if ( i < xMax ) {
do {
carry = ((r[i] = x[i] + 1) ==
0);
} while ( ++i < xMax && carry);
}
if (carry) {
r[i] = 1;
result.length = ++i;
return result;
}
}
// Copy the rest
if( i < xMax ) {
do {
r[i] = x[i];
} while ( ++i < xMax );
}
result.Normalize();
return result;
}
public static BigInteger Subtract ( BigInteger big, BigInteger
small ) {
BigInteger result = new BigInteger ( Sign.Positive,
big.length );
uint [] r = result.data, b = big.data, s = small.data;
uint i = 0, c = 0;
do {
uint x = s[i];
if ( ((x += c) < c) | ( (r[i] = b[i] - x) > ~x
) )
c = 1;
else
c = 0;
} while ( ++i < small.length );
if ( i == big.length ) goto fixup;
if ( c == 1 ) {
do {
r[i] = b[i] - 1;
} while (b[i++] == 0 && i < big.length);
if ( i == big.length ) goto fixup;
}
do {
r[i] = b[i];
} while (++i < big.length);
fixup:
result.Normalize();
return result;
}
public static void MinusEq ( BigInteger big, BigInteger small
) {
uint[] b = big.data, s = small.data;
uint i = 0, c = 0;
do {
uint x = s[i];
if ( ((x += c) < c) | ( (b[i] -= x) > ~x ) )
c = 1;
else
c = 0;
} while ( ++i < small.length );
if ( i == big.length ) goto fixup;
if ( c == 1 ) {
do {
b[i]--;
} while (b[i++] == 0 && i < big.length);
if ( i == big.length ) goto fixup;
}
fixup:
// Normalize length
while(big.length > 0 && big.data[big.length-1]
== 0) big.length--;
// Check for zero
if (big.length == 0) {
big.length++;
}
}
public static void PlusEq ( BigInteger bi1, BigInteger bi2 ) {
uint[] x, y;
uint yMax, xMax, i = 0;
bool flag = false;
// x should be bigger
if (bi1.length < bi2.length){
flag = true;
x = bi2.data;
xMax = bi2.length;
y = bi1.data;
yMax = bi1.length;
}
else {
x = bi1.data;
xMax = bi1.length;
y = bi2.data;
yMax = bi2.length;
}
uint[] r = bi1.data;
ulong sum = 0;
// Add common parts of both numbers
do {
sum = ((ulong)x[i]) +
((ulong)y[i]) + sum;
r[i] = (uint)sum;
sum >>= 32;
} while ( ++i < yMax );
// Copy remainder of longer number while carry
propagation is required
bool carry = (sum != 0);
if (carry){
if ( i < xMax ) {
do {
carry = ((r[i] = x[i] + 1) ==
0);
} while ( ++i < xMax && carry);
}
if (carry) {
r[i] = 1;
bi1.length = ++i;
return;
}
}
// Copy the rest
if( flag && i < xMax - 1 ) {
do {
r[i] = x[i];
} while ( ++i < xMax);
}
bi1.length = xMax + 1;
bi1.Normalize();
}
#endregion
#region Compare
/// <summary>
/// Compares two BigInteger
/// </summary>
/// <param name="bi1">A BigInteger</param>
/// <param name="bi2">A BigInteger</param>
/// <returns>The sign of bi1 - bi2</returns>
public static Sign Compare( BigInteger bi1, BigInteger bi2 ) {
//
// Step 2. Compare the lengths
//
uint l1 = bi1.length, l2 = bi2.length;
while(l1 > 0 && bi1.data[l1-1] == 0) l1--;
while(l1 > 0 && bi2.data[l2-1] == 0) l2--;
if (l1 == 0 && l2 == 0) return Sign.Zero;
// bi1 len < bi2 len
if (l1 < l2) return Sign.Negative;
// bi1 len > bi2 len
else if (l1 > l2) return Sign.Positive;
//
// Step 3. Compare the bits
//
// bi1 and bi2 have same sign, length
uint pos = l1 - 1;
do {
if (bi1.data[pos] < bi2.data[pos])
return Sign.Negative;
else if (bi1.data[pos] > bi2.data[pos])
return Sign.Positive;
} while (pos-- > 0);
// same number
return Sign.Zero;
}
public static Sign CompareSameSign( BigInteger bi1, BigInteger
bi2 ) {
uint l1 = bi1.length, l2 = bi2.length;
while(l1 > 0 && bi1.data[l1-1] == 0) l1--;
while(l2 > 0 && bi2.data[l2-1] == 0) l2--;
if (l1 == 0 && l2 == 0) return Sign.Zero;
// bi1 len < bi2 len
if (l1 < l2) return Sign.Negative;
// bi1 len > bi2 len
else if (l1 > l2) return Sign.Positive;
// bi1 and bi2 have same sign, length
else {
uint pos = bi1.length - 1;
do {
if (bi1.data[pos] < bi2.data[pos])
return Sign.Negative;
else if (bi1.data[pos] > bi2.data[pos])
return Sign.Positive;
} while (pos-- > 0);
// same number
return Sign.Zero;
}
}
#endregion
#region Division
#region Dword
/// <summary>
/// Performs n / d and n % d in one operation.
/// </summary>
/// <param name="n">A BigInteger, upon exit this will hold n /
d</param>
/// <param name="d">The divisor</param>
/// <returns>n % d</returns>
public static uint SingleByteDivideInPlace ( BigInteger n,
uint d) {
ulong r = 0;
uint i = n.length;
while (i-- > 0) {
r <<= 32;
r |= n.data[i];
n.data[i] = (uint)(r / d);
r %= d;
}
n.Normalize();
return (uint)r;
}
public static uint DwordMod ( BigInteger n, uint d ) {
ulong r = 0;
uint i = n.length;
while (i-- > 0) {
r <<= 32;
r |= n.data[i];
r %= d;
}
return (uint)r;
}
public static BigInteger DwordDiv ( BigInteger n, uint d ) {
BigInteger ret = new BigInteger ( Sign.Positive,
n.length );
ulong r = 0;
uint i = n.length;
while (i-- > 0) {
r <<= 32;
r |= n.data[i];
ret.data[i] = (uint)(r / d);
r %= d;
}
ret.Normalize();
return ret;
}
public static BigInteger[] DwordDivMod( BigInteger n, uint d,
Sign dSign) {
BigInteger ret = new BigInteger ( Sign.Positive ,
n.length );
ulong r = 0;
uint i = n.length;
while (i-- > 0) {
r <<= 32;
r |= n.data[i];
ret.data[i] = (uint)(r / d);
r %= d;
}
ret.Normalize();
BigInteger rem = new BigInteger( (uint)r);
rem.Normalize();
return new BigInteger[] {ret, rem};
}
#endregion
#region BigNum
public static BigInteger[] multiByteDivide(BigInteger bi1,
BigInteger bi2) {
if( CompareSameSign(bi1, bi2) == Sign.Negative )
return new BigInteger[2] { 0, new
BigInteger(bi1) };
bi1.Normalize(); bi2.Normalize();
if( bi2.length == 1 )
return DwordDivMod(bi1, bi2.data[0],
Sign.Positive);
uint remainderLen = bi1.length + 1;
int divisorLen = (int)bi2.length + 1;
uint mask = 0x80000000;
uint val = bi2.data[bi2.length - 1];
int shift = 0;
int resultPos = (int)bi1.length - (int)bi2.length;
while(mask != 0 && (val & mask) == 0) {
shift++; mask >>= 1;
}
BigInteger quot = new BigInteger ( Sign.Positive,
bi1.length - bi2.length + 1 );
BigInteger rem = (bi1 << shift);
uint [] remainder = rem.data;
bi2 = bi2 << shift;
int j = (int)(remainderLen - bi2.length);
int pos = (int)remainderLen - 1;
uint firstDivisorByte = bi2.data[bi2.length-1];
ulong secondDivisorByte = bi2.data[bi2.length-2];
while(j > 0) {
ulong dividend = ((ulong)remainder[pos] << 32)
+ (ulong)remainder[pos-1];
ulong q_hat = dividend /
(ulong)firstDivisorByte;
ulong r_hat = dividend %
(ulong)firstDivisorByte;
do {
if(q_hat == 0x100000000 ||
(q_hat * secondDivisorByte) >
((r_hat << 32) + remainder[pos-2])) {
q_hat--;
r_hat +=
(ulong)firstDivisorByte;
if(r_hat < 0x100000000)
continue;
}
} while (false);
//
// At this point, q_hat is either exact, or
one too large
// (more likely to be exact) so, we attempt to
multiply the
// divisor by q_hat, if we get a borrow, we
just subtract
// one from q_hat and add the divisor back.
//
uint t;
uint dPos = 0;
int nPos = pos - divisorLen + 1;
ulong mc = 0;
uint uint_q_hat = (uint)q_hat;
do {
mc += (ulong)bi2.data[dPos] *
(ulong)uint_q_hat;
t = remainder[nPos];
remainder[nPos] -= (uint)mc;
mc >>= 32;
if( remainder[nPos] > t) mc++;
dPos++; nPos++;
} while ( dPos < divisorLen );
nPos = pos - divisorLen + 1;
dPos = 0;
// Overestimate
if( mc != 0 ) {
q_hat--;
ulong sum = 0;
do {
sum =
((ulong)remainder[nPos]) + ((ulong)bi2.data[dPos]) + sum;
remainder[nPos] = (uint)sum;
sum >>= 32;
dPos++; nPos++;
} while ( dPos < divisorLen );
}
quot.data[resultPos--] = (uint)q_hat;
pos--;
j--;
}
quot.Normalize();
rem.Normalize();
BigInteger[] ret = new BigInteger[2] { quot, rem };
if( shift != 0 )
ret[1] >>= shift;
return ret;
}
#endregion
#endregion
#region Shift
public static BigInteger LeftShift (BigInteger bi, int n) {
if( n == 0) return new BigInteger( bi, bi.length + 1);
int w = n >> 5;
n &= ((1 << 5) - 1);
BigInteger ret = new BigInteger( Sign.Positive,
bi.length + 1 + (uint)w);
uint i = 0, l = bi.length;
if (n != 0) {
uint x, carry = 0;
while (i < l) {
x = bi.data[i];
ret.data[i + w] = (x << n) | carry;
carry = x >> (32 - n);
i++;
}
ret.data[i + w] = carry;
}
else {
while (i < l) {
ret.data[i + w] = bi.data[i];
i++;
}
}
ret.Normalize();
return ret;
}
public static BigInteger RightShift( BigInteger bi, int n) {
if( n == 0) return new BigInteger( bi );
int w = n >> 5;
int s = n & ((1 << 5) - 1);
BigInteger ret = new BigInteger( Sign.Positive,
bi.length - (uint)w + 1);
uint l = (uint)ret.data.Length - 1;
if (s != 0) {
uint x, carry = 0;
while (l-- > 0) {
x = bi.data[l + w];
ret.data[l] = (x >> n) | carry;
carry = x << (32 - n);
}
}
else {
while (l-- > 0)
ret.data[l] = bi.data[l + w];
}
ret.Normalize();
return ret;
}
#endregion
#region Multiply
public static BigInteger MultiplyByDword ( BigInteger n, uint
f, Sign s ) {
BigInteger ret = new BigInteger( Sign.Positive,
n.length + 1 );
uint i = 0;
ulong c = 0;
do {
c += (ulong)n.data[i] * (ulong)f;
ret.data[i] = (uint)(c & 0xFFFFFFFF);
c >>= 32;
} while (++i < n.length);
ret.data[i] = (uint)c;
ret.Normalize();
return ret;
}
/// <summary>
/// Multiplies the data in x[xOffset:xOffset+xLen] by
/// y[yOffset:yOffset+yLen] and puts it into
/// d[dOffset:dOffset+xLen+yLen].
/// </summary>
/// <remarks>
/// This code is unsafe! It is the caller's responsibility to
make
/// sure that it is safe to access x[xOffset:xOffset+xLen],
/// y[yOffset:yOffset+yLen], and d[dOffset:dOffset+xLen+yLen].
/// </remarks>
public static unsafe void Multiply ( uint[] x, uint xOffset,
uint xLen, uint[] y, uint yOffset, uint yLen, uint[] d, uint dOffset ) {
fixed (uint* xx = x, yy = y, dd = d) {
uint* xP = xx + xOffset,
xE = xP + xLen,
yB = yy + yOffset,
yE = yB + yLen,
dB = dd + dOffset;
for(; xP < xE; xP++, dB++) {
if(*xP == 0) continue;
ulong mcarry = 0;
uint* dP = dB;
for(uint* yP = yB; yP < yE; yP++,
dP++) {
mcarry += ((ulong)*xP *
(ulong)*yP) + (ulong)*dP;
*dP = (uint)(mcarry &
0xFFFFFFFF);
mcarry >>= 32;
}
if(mcarry != 0)
*dP = (uint)mcarry;
}
}
}
/// <summary>
/// Multiplies the data in x[xOffset:xOffset+xLen] by
/// y[yOffset:yOffset+yLen] and puts the low mod words into
/// d[dOffset:dOffset+mod].
/// </summary>
/// <remarks>
/// This code is unsafe! It is the caller's responsibility to
make
/// sure that it is safe to access x[xOffset:xOffset+xLen],
/// y[yOffset:yOffset+yLen], and d[dOffset:dOffset+mod].
/// </remarks>
public static unsafe void MultiplyMod2p32pmod ( uint[] x, int
xOffset, int xLen, uint[] y, int yOffest, int yLen, uint[] d, int dOffset, int mod ) {
fixed (uint* xx = x, yy = y, dd = d) {
uint* xP = xx + xOffset,
xE = xP + xLen,
yB = yy + yOffest,
yE = yB + yLen,
dB = dd + dOffset,
dE = dB + mod;
for(; xP < xE; xP++, dB++) {
if(*xP == 0) continue;
ulong mcarry = 0;
uint* dP = dB;
for(uint* yP = yB; yP < yE && dP < dE;
yP++, dP++) {
mcarry += ((ulong)*xP *
(ulong)*yP) + (ulong)*dP;
*dP = (uint)mcarry;
mcarry >>= 32;
}
if(mcarry != 0 && dP < dE)
*dP = (uint)mcarry;
}
}
}
public static unsafe void SquarePositive(BigInteger bi, ref
uint[] wkSpace) {
Debug.Assert(wkSpace.Length >= bi.length << 1);
uint [] t = wkSpace;
wkSpace = bi.data;
uint [] d = bi.data;
uint dl = bi.length;
bi.data = t;
fixed (uint* dd = d, tt = t) {
uint* ttE = tt + t.Length;
// Clear the dest
for(uint* ttt = tt; ttt < ttE; ttt++)
*ttt = 0;
uint* dP = dd, tP = tt;
for(uint i = 0; i < dl; i++, dP++) {
if (*dP == 0)
continue;
ulong mcarry = 0;
uint bi1val = *dP;
uint* dP2 = dP + 1, tP2 = tP + 2*i + 1;
for(uint j = i + 1; j < dl; j++,
tP2++, dP2++) {
// k = i + j
mcarry += ((ulong)bi1val *
(ulong)*dP2) + *tP2;
*tP2 = (uint)mcarry;
mcarry >>= 32;
}
if(mcarry != 0)
*tP2 = (uint)mcarry;
}
// Double t. Inlined for speed.
tP = tt;
uint x, carry = 0;
while (tP < ttE) {
x = *tP;
*tP = (x << 1) | carry;
carry = x >> (32 - 1);
tP++;
}
if ( carry != 0 ) *tP = carry;
// Add in the diagnals
dP = dd;
tP = tt;
for (uint* dE = dP + dl; (dP < dE); dP++,
tP++ ) {
ulong val = (ulong)*dP * (ulong)*dP +
*tP;
*tP = (uint)val;
val >>= 32;
*(++tP) += (uint)val;
if (*tP < (uint)val) {
uint* tP3 = tP;
// Account for the first carry
(*++tP3)++;
// Keep adding until no carry
while((*tP3++) == 0x0)
(*tP3)++;
}
}
bi.length <<= 1;
// Normalize length
while( tt[bi.length-1] == 0 && bi.length > 1 )
bi.length--;
}
}
public static bool Double( uint[] u, int l ) {
uint x, carry = 0;
uint i = 0;
while (i < l) {
x = u[i];
u[i] = (x << 1) | carry;
carry = x >> (32 - 1);
i++;
}
if ( carry != 0 ) u[l] = carry;
return carry != 0;
}
#endregion
#region Number Theory
public static BigInteger gcd(BigInteger a, BigInteger b) {
BigInteger x = a;
BigInteger y = b;
BigInteger g = y;
while(x.length > 1) {
g = x;
x = y % x;
y = g;
}
if ( x == 0 ) return g;
// TODO: should we have something here if we can
convert to long?
//
// Now we can just do it with single precision. I am
using the extended gcd method,
// as it should be faster.
//
uint yy = x.data[0];
uint xx = y % yy;
int t = 0;
while ( ((xx | yy) & 1) == 0) {
xx >>= 1; yy >>= 1; t++;
}
while ( xx != 0 ){
while ( (xx & 1) == 0 ) xx >>= 1;
while ( (yy & 1) == 0 ) yy >>= 1;
if ( xx >= yy )
xx = (xx - yy) >> 1;
else
yy = (yy - xx) >> 1;
}
return yy << t;
}
public static BigInteger modInverse(BigInteger bi, BigInteger
modulus) {
BigInteger[] p = { 0, 1 };
BigInteger[] q = new BigInteger[2]; // quotients
BigInteger[] r = { 0, 0 }; // remainders
int step = 0;
BigInteger a = modulus;
BigInteger b = bi;
ModulusRing mr = new ModulusRing( modulus );
while(b.length > 1 || (b.length == 1 && b.data[0] !=
0)) {
if(step > 1) {
p[0] += modulus;
BigInteger pval = (p[0] - (p[1] *
q[0]) % modulus) % modulus;
p[0] = p[1];
p[1] = pval;
}
BigInteger[] divret = multiByteDivide(a, b);
q[0] = q[1];
r[0] = r[1];
q[1] = divret[0];
r[1] = divret[1];
a = b;
b = divret[1];
step++;
}
if(r[0].length > 1 || (r[0].length == 1 &&
r[0].data[0] != 1))
throw (new ArithmeticException("No inverse!"));
// TODO: clean this up
p[0] += modulus;
BigInteger t = (p[1] * q[0]) % modulus;
return (p[0] - t) % modulus;
//BigInteger result = ((p[0] - (p[1] * q[0])) %
modulus);
// NEGPROBLEM
//if(result.isNeg)
// result += modulus; // get the least positive
modulus
}
#endregion
#region Prime Testing
/// <summary>
/// Probabilistic prime test based on Rabin-Miller's test
/// </summary>
/// <param name="bi" type="BigInteger.BigInteger">
/// <para>
/// The number to test.
/// </para>
/// </param>
/// <param name="confidence" type="int">
/// <para>
/// The number of chosen bases. The test
has a
/// 1/4^confidence chance of falsely
returning True.
/// </para>
/// </param>
/// <returns>
/// <para>
/// True if "this" is a strong pseudoprime
to randomly chosen
/// bases.
/// </para>
/// <para>
/// False if "this" is definitely NOT
prime.
/// </para>
/// </returns>
public static bool RabinMillerTest(BigInteger bi, int
confidence) {
BigInteger thisVal;
thisVal = bi;
// calculate values of s and t
BigInteger p_sub1 = thisVal - 1;
int s = 0;
for(int index = 0; index < p_sub1.length; index++) {
uint mask = 0x01;
for(int i = 0; i < 32; i++) {
if((p_sub1.data[index] & mask) != 0) {
goto end_loop;
}
mask <<= 1;
s++;
}
}
end_loop:
BigInteger t = p_sub1 >> s;
int bits = thisVal.bitCount();
BigInteger a = null;
RandomNumberGenerator rng =
RandomNumberGenerator.Create();
ModulusRing mr = new ModulusRing( thisVal );
for(int round = 0; round < confidence; round++) {
bool done = false;
while(!done) { // generate a < n
a = genRandom( bits, rng );
uint byteLen = a.length;
// make sure "a" is not 0
if (a > 1 && a < thisVal)
done = true;
}
BigInteger gcdTest = a.gcd(thisVal);
if (gcdTest != 1) return false;
BigInteger b = mr.Pow(a, t);
if (b == 1) continue; // a^t mod p = 1
bool result = false;
for(int j = 0; j < s; j++) {
if(b == p_sub1) { //
a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
result = true;
break;
}
b = (b * b) % thisVal;
}
if(result == false)
return false;
}
return true;
}
public static bool SmallPrimeSppTest( BigInteger bi, int
confidence) {
BigInteger thisVal;
thisVal = bi;
// calculate values of s and t
BigInteger p_sub1 = thisVal - 1;
int s = 0;
for(int index = 0; index < p_sub1.length; index++) {
uint mask = 0x01;
for(int i = 0; i < 32; i++) {
if((p_sub1.data[index] & mask) != 0) {
goto end_loop;
}
mask <<= 1;
s++;
}
}
end_loop:
BigInteger t = p_sub1 >> s;
ModulusRing mr = new ModulusRing( thisVal );
for(int round = 0; round < confidence; round++) {
uint a = smallPrimes[round];
BigInteger b = mr.Pow(a, t);
if (b == 1) continue; // a^t mod p = 1
bool result = false;
for(int j = 0; j < s; j++) {
if(b == p_sub1) { //
a^((2^j)*t) mod p = p-1 for some 0 <= j <= s-1
result = true;
break;
}
b = (b * b) % thisVal;
}
if(result == false)
return false;
}
return true;
}
#endregion
}
}
}
