1 | #ifndef __HALTON_H
|
---|
2 | #define __HALTON_H
|
---|
3 |
|
---|
4 | #include <iostream>
|
---|
5 | using namespace std;
|
---|
6 |
|
---|
7 | class Halton2 {
|
---|
8 | static float _invBases[2];
|
---|
9 | float _prev[2];
|
---|
10 |
|
---|
11 | float halton(float baseRec, float prev) const {
|
---|
12 | float r = 1 - prev - 1e-10f;
|
---|
13 | if (baseRec < r)
|
---|
14 | return prev + baseRec;
|
---|
15 | float h = baseRec;
|
---|
16 | float hh;
|
---|
17 | do {
|
---|
18 | hh = h;
|
---|
19 | h *= baseRec;
|
---|
20 | } while (h >= r);
|
---|
21 | return prev + hh + h - 1;
|
---|
22 | }
|
---|
23 |
|
---|
24 | public:
|
---|
25 |
|
---|
26 | void Reset() {
|
---|
27 | _prev[0] =_prev[1] = 0;
|
---|
28 | }
|
---|
29 |
|
---|
30 | Halton2() {
|
---|
31 | _invBases[0] = 1.0f/2;
|
---|
32 | _invBases[1] = 1.0f/3;
|
---|
33 | Reset();
|
---|
34 | }
|
---|
35 |
|
---|
36 | void
|
---|
37 | GetNext(float &a, float &b) {
|
---|
38 | a = halton(_invBases[0], _prev[0]);
|
---|
39 | b = halton(_invBases[1], _prev[1]);
|
---|
40 | _prev[0] = a;
|
---|
41 | _prev[1] = b;
|
---|
42 | }
|
---|
43 | };
|
---|
44 |
|
---|
45 |
|
---|
46 | /**
|
---|
47 | * Assert whether the argument is a prime number.
|
---|
48 | * @param number the number to be checked
|
---|
49 | */
|
---|
50 | inline bool IsPrime(const int number) {
|
---|
51 | bool isIt = true;
|
---|
52 | for(int i = 2; i < number; i++) {
|
---|
53 | if(number % i == 0) {
|
---|
54 | isIt = false;
|
---|
55 | break;
|
---|
56 | }
|
---|
57 | }
|
---|
58 | if(number == 2) {
|
---|
59 | isIt = false;
|
---|
60 | }
|
---|
61 | return isIt;
|
---|
62 | }
|
---|
63 |
|
---|
64 | /**
|
---|
65 | * Find the nth prime number.
|
---|
66 | * @param index the ordinal position in the sequence
|
---|
67 | */
|
---|
68 | inline int FindPrime(const int index) {
|
---|
69 | if(index < 1) {
|
---|
70 | cerr<<"FindPrime: The argument must be non-negative."<<endl;
|
---|
71 | return -1;
|
---|
72 | }
|
---|
73 | int prime = 1;
|
---|
74 | int found = 1;
|
---|
75 | while(found != index) {
|
---|
76 | prime += 2;
|
---|
77 | if(IsPrime(prime) == true) {
|
---|
78 | found++;
|
---|
79 | }
|
---|
80 | }
|
---|
81 | return prime;
|
---|
82 | }
|
---|
83 |
|
---|
84 | struct HaltonSequence {
|
---|
85 | public:
|
---|
86 | int index;
|
---|
87 |
|
---|
88 | HaltonSequence():index(1) {}
|
---|
89 |
|
---|
90 | void Reset() {
|
---|
91 | index = 1;
|
---|
92 | }
|
---|
93 |
|
---|
94 | void GenerateNext() {
|
---|
95 | index++;
|
---|
96 | }
|
---|
97 |
|
---|
98 | /**
|
---|
99 | * Returns the nth number in the sequence, taken from a specified dimension.
|
---|
100 | * @param index the ordinal position in the sequence
|
---|
101 | * @param dimension the dimension
|
---|
102 | */
|
---|
103 |
|
---|
104 | double GetNumber(const int dimension) {
|
---|
105 | int base = FindPrime(dimension);
|
---|
106 | if(base == 1) {
|
---|
107 | base++; //The first dimension uses base 2.
|
---|
108 | }
|
---|
109 | double remainder;
|
---|
110 | double output = 0.0;
|
---|
111 | double fraction = 1.0 / (double)base;
|
---|
112 | int N1 = 0;
|
---|
113 | int copyOfIndex = index;
|
---|
114 | if((base >= 2) && (index >= 1)) {
|
---|
115 | // changed by matt
|
---|
116 | //if(base >= 2 & index >= 1) {
|
---|
117 | while(copyOfIndex > 0) {
|
---|
118 | N1 = (copyOfIndex / base);
|
---|
119 | remainder = copyOfIndex % base;
|
---|
120 | output += fraction * remainder;
|
---|
121 | copyOfIndex = (int)(copyOfIndex / base);
|
---|
122 | fraction /= (double)base;
|
---|
123 | }
|
---|
124 | return output;
|
---|
125 | }
|
---|
126 | else {
|
---|
127 | cerr<<"Error generating Halton sequence."<<endl;
|
---|
128 | exit(1);
|
---|
129 | }
|
---|
130 | }
|
---|
131 | };
|
---|
132 |
|
---|
133 | extern Halton2 halton2;
|
---|
134 |
|
---|
135 | #endif
|
---|