CluE  1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
primegenerator.cpp
Go to the documentation of this file.
1 #include <set>
2 
3 #include "primegenerator.h"
4 
5 namespace CluE
6 {
7 
8 unsigned long long PrimeGenerator::getPrime(unsigned long long lowerBound, unsigned long long upperBound, double errorProbability)
9 {
10  if(lowerBound < 2)
11  lowerBound = 2;
12  if(lowerBound > upperBound)
13  throw "TODO";
14  if(errorProbability <= 0 || errorProbability > 0.25)
15  throw "TODO";
16 
17  std::set<unsigned long long> checkedNumbers;
19  std::uniform_int_distribution<unsigned long long> dis(lowerBound, upperBound);
20 
21  // Sample from the given interval and check if the number is prime
22  while(checkedNumbers.size() < upperBound - lowerBound + 1)
23  {
24  unsigned long long candidate = 0;
25  while(candidate == 0 || checkedNumbers.count(candidate) > 0)
26  candidate = dis(gen);
27  checkedNumbers.insert(candidate);
28 
29  if(isMillerRabinPrime(candidate, errorProbability))
30  return candidate;
31  }
32 
33  throw "TODO";
34 }
35 
36 bool PrimeGenerator::isMillerRabinPrime(unsigned long long candidate, double errorProbability)
37 {
38  if(candidate == 2)
39  return true;
40  else if((candidate & 1) == 0) // mod 2
41  return false;
42 
43  if(errorProbability <= 0 || errorProbability > 0.25)
44  throw "TODO";
45 
46  // Write candidate = 2^s * r + 1
47  unsigned long long r = candidate-1;
48  int s = 0;
49  while((r & 1) == 0) // mod 2
50  {
51  s++;
52  r /= 2;
53  }
54 
55  // Random number generator
56  std::set<unsigned long long> usedNumbers;
58  std::uniform_int_distribution<unsigned long long> dis(2, candidate-1);
59 
60  double currentErrorProbability = 0.5;
61  for(unsigned long i = 1; errorProbability < currentErrorProbability; i++)
62  {
63  // Sample an unused number
64  unsigned long long randomBase = 0;
65  while(randomBase == 0 || usedNumbers.count(randomBase) > 0)
66  randomBase = dis(gen);
67  usedNumbers.insert(randomBase);
68 
69  // Calculate y=randomBase^r mod candidate
70  unsigned long long y = potMod(randomBase, r, candidate);
71 
72  // 1-Condition derived from Femat's little theorem
73  bool onlyOnes = (y == 1 || y == candidate-1);
74 
75  // Calculate y^2 mod candidate, (y^2 mod candidate)^2 mod candidate, ...
76  for(int j = 1; j <= s; j++)
77  {
78  y = (y*y) % candidate;
79  if(onlyOnes && y != 1)
80  return false;
81  else if(!onlyOnes)
82  {
83  if(y == 1)
84  return false;
85  else
86  onlyOnes = (y == candidate-1);
87  }
88  }
89  if(y != 1)
90  return false;
91 
92  if(i >= candidate/4)
93  return true;
94  currentErrorProbability *= 0.5;
95  }
96 
97  return true;
98 }
99 
100 unsigned long long PrimeGenerator::potMod(unsigned long long x, unsigned long long k, unsigned long long m)
101 {
102  unsigned long long result = 1;
103  unsigned long long base = x;
104  while(k > 0)
105  {
106  if((k & 1) == 1)
107  result = (result * base) % m;
108  k >>= 1;
109  base = (base * base) % m;
110  }
111 
112  return result;
113 }
114 
115 }
Encapsulates an STL random generator.
bool isMillerRabinPrime(unsigned long long candidate, double errorProbability)
Miller-Rabin primality test.
unsigned long long potMod(unsigned long long x, unsigned long long b, unsigned long long m)
Calculates "x^b mod m".
static RandomGenerator getRandomGenerator()
Definition: randomness.h:23
unsigned long long getPrime(unsigned long long lowerBound, unsigned long long upperBound, double errorProbability)
Returns a number [lowerBound, upperBound] which is prime with probablity 1-errorProbability.