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