Skip to content

Instantly share code, notes, and snippets.

@yutopio
Created May 24, 2013 14:43
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save yutopio/5644022 to your computer and use it in GitHub Desktop.
Save yutopio/5644022 to your computer and use it in GitHub Desktop.
Enumerates some pseudoprimes.
using System;
using System.IO;
unsafe class Program
{
const int size = 1000000000;
static bool[] sieve = null;
static void Main()
{
var start = DateTime.Now;
CreateSieve();
var end = DateTime.Now;
Console.WriteLine("Sieve creation: {0}", end - start);
for (ulong i = 0; i < 100; i++)
{
var j = 0;
for (ulong k = i * 10000000 + 1; j < 5000000; j++, k += 2)
// if (!IsPrime(k) && MillarRabin(k))
if (!sieve[k] && MillarRabin(k))
Console.WriteLine(k);
}
}
static void CreateSieve()
{
sieve = new bool[size];
fixed (bool* pt = sieve)
{
ulong* init = (ulong*)pt;
init[0] = 0x0100010001010000;
for (int i = 1; i < size / 8; i++)
init[i] = 0x0100010001000100;
int upper = (int)Math.Ceiling(Math.Sqrt(size));
if (upper % 2 == 1) upper++;
for (int i = 3; i < upper; i++)
if (pt[i])
for (int j = i * i; j < size; j += i)
pt[j] = false;
}
}
/* static bool IsPrime(ulong n)
{
if (n < 2) return false;
else if (n == 2) return true;
else if ((n & 1) == 0) return false;
var j = (ulong)Math.Sqrt(n);
for (ulong i = 3; i <= j; i += 2)
if (n % i == 0) return false;
return true;
}
*/
static bool MillarRabin(ulong n)
{
if (n < 2) return false;
else if (n == 2) return true;
else if ((n & 1) == 0) return false;
var oddPower = n - 1;
while ((oddPower & 1) == 0) oddPower >>= 1;
foreach (var @base in new ulong[] { 2, 3, 5, 13 })
{
var power = oddPower;
var temp = PowMod(@base, power, n);
if (temp == 1) continue;
while (power != n - 1 && temp != n - 1)
{
temp = (temp * temp) % n;
power <<= 1;
}
if (power == n - 1) return false;
}
return true;
}
static ulong PowMod(ulong @base, ulong power, ulong modulo)
{
ulong ret = 1;
@base %= modulo;
while (power != 0)
{
if ((power & 1) == 1) ret = (ret * @base) % modulo;
@base = (@base * @base) % modulo;
power >>= 1;
}
return ret;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment