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