[355] | 1 | #ifndef __HALTON_H
|
---|
| 2 | #define __HALTON_H
|
---|
| 3 |
|
---|
[492] | 4 | #include <iostream>
|
---|
| 5 | using namespace std;
|
---|
| 6 |
|
---|
[355] | 7 | class Halton2 {
|
---|
[492] | 8 | static float _invBases[2];
|
---|
| 9 | float _prev[2];
|
---|
| 10 |
|
---|
| 11 | float halton(float baseRec, float prev) const {
|
---|
[497] | 12 | float r = 1 - prev - 1e-10f;
|
---|
[492] | 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() {
|
---|
[497] | 31 | _invBases[0] = 1.0f/2;
|
---|
| 32 | _invBases[1] = 1.0f/3;
|
---|
[492] | 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 | };
|
---|
[355] | 44 |
|
---|
[492] | 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;
|
---|
[355] | 56 | }
|
---|
[492] | 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 {
|
---|
[355] | 85 | public:
|
---|
[492] | 86 | int index;
|
---|
[355] | 87 |
|
---|
[492] | 88 | HaltonSequence():index(1) {}
|
---|
[355] | 89 |
|
---|
[492] | 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;
|
---|
[497] | 114 | if((base >= 2) && (index >= 1)) {
|
---|
| 115 | // changed by matt
|
---|
| 116 | //if(base >= 2 & index >= 1) {
|
---|
[492] | 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;
|
---|
[355] | 125 | }
|
---|
[492] | 126 | else {
|
---|
| 127 | cerr<<"Error generating Halton sequence."<<endl;
|
---|
| 128 | exit(1);
|
---|
[355] | 129 | }
|
---|
[492] | 130 | }
|
---|
[355] | 131 | };
|
---|
| 132 |
|
---|
| 133 | extern Halton2 halton2;
|
---|
| 134 |
|
---|
| 135 | #endif
|
---|