[2839] | 1 | #include "common.h"
|
---|
| 2 | #include "Halton.h"
|
---|
| 3 |
|
---|
| 4 |
|
---|
| 5 |
|
---|
[3227] | 6 | void Halton::TestHalton(int n, int dim)
|
---|
| 7 | {
|
---|
| 8 | Halton h(dim);
|
---|
[2839] | 9 |
|
---|
[3227] | 10 | for (int i = 1; i <= n; ++ i)
|
---|
| 11 | {
|
---|
| 12 | std::cout << "halton " << i << " of dim " << dim << " = " << h.GetNext() << std::endl;
|
---|
| 13 | }
|
---|
| 14 | }
|
---|
[2839] | 15 |
|
---|
| 16 |
|
---|
[3227] | 17 | void Halton::TestPrime()
|
---|
| 18 | {
|
---|
| 19 | int dim = 2;
|
---|
[2839] | 20 |
|
---|
[3227] | 21 | for (int i = 1; i < 10; ++ i)
|
---|
| 22 | {
|
---|
| 23 | std::cout << "prime for dim " << i << " " << Halton::FindPrime(i) << std::endl;
|
---|
| 24 | }
|
---|
| 25 | }
|
---|
| 26 |
|
---|
| 27 |
|
---|
| 28 | Halton::Halton(int dim): mIndex(0)
|
---|
[2839] | 29 | {
|
---|
[3227] | 30 | mBase = FindPrime(dim);
|
---|
| 31 | }
|
---|
[2839] | 32 |
|
---|
| 33 |
|
---|
[3227] | 34 | float Halton::GetNext()
|
---|
| 35 | {
|
---|
| 36 | ++ mIndex;
|
---|
[2839] | 37 |
|
---|
[3227] | 38 | float result = .0f;
|
---|
| 39 | float fraction = 1.0f / (float)mBase;
|
---|
| 40 | int idx = mIndex;
|
---|
| 41 |
|
---|
| 42 | while (idx > 0)
|
---|
[2839] | 43 | {
|
---|
[3227] | 44 | int digit = idx % mBase;
|
---|
| 45 | result += fraction * (float)digit;
|
---|
| 46 |
|
---|
| 47 | idx = (idx - digit) / mBase;
|
---|
| 48 | fraction /= (float)mBase;
|
---|
[2839] | 49 | }
|
---|
| 50 |
|
---|
[3227] | 51 | mIndex %= 100000000;
|
---|
| 52 | return result;
|
---|
[2839] | 53 | }
|
---|
| 54 |
|
---|
| 55 |
|
---|
[3227] | 56 | bool Halton::IsPrime(int n)
|
---|
[2839] | 57 | {
|
---|
[3227] | 58 | bool isPrime = true;
|
---|
[2839] | 59 |
|
---|
[3227] | 60 | for (int i = 2; i < n; ++ i)
|
---|
| 61 | {
|
---|
| 62 | if (n % i == 0)
|
---|
| 63 | {
|
---|
| 64 | isPrime = false;
|
---|
| 65 | break;
|
---|
| 66 | }
|
---|
| 67 | }
|
---|
| 68 |
|
---|
| 69 | return isPrime;
|
---|
[2839] | 70 | }
|
---|
[3227] | 71 |
|
---|
| 72 |
|
---|
| 73 | int Halton::FindPrime(int idx)
|
---|
| 74 | {
|
---|
| 75 | if (idx < 1)
|
---|
| 76 | {
|
---|
| 77 | std::cerr << "error: cannot find " << idx << "th prime" << std::endl;
|
---|
| 78 | return 0;
|
---|
| 79 | }
|
---|
| 80 |
|
---|
| 81 | // only even prime number
|
---|
| 82 | if (idx == 1) return 2;
|
---|
| 83 |
|
---|
| 84 | int number = 3;
|
---|
| 85 | int primeFound = 1;
|
---|
| 86 |
|
---|
| 87 | while (1)
|
---|
| 88 | {
|
---|
| 89 | if (IsPrime(number)) ++ primeFound;
|
---|
| 90 | if (primeFound == idx) break;
|
---|
| 91 |
|
---|
| 92 | number += 2;
|
---|
| 93 | }
|
---|
| 94 |
|
---|
| 95 | return number;
|
---|
| 96 | }
|
---|
| 97 |
|
---|
| 98 |
|
---|
| 99 | void HaltonSequence::TestHalton(int n, int dim)
|
---|
| 100 | {
|
---|
| 101 | HaltonSequence h(dim);
|
---|
| 102 | float *haltons = new float[dim];
|
---|
| 103 |
|
---|
| 104 | for (int i = 1; i <= n; ++ i)
|
---|
| 105 | {
|
---|
| 106 | h.GetNext(haltons);
|
---|
| 107 | std::cout << "halton " << i << " of dim " << dim << " = ";
|
---|
| 108 |
|
---|
| 109 | for (int j = 0; j < dim; ++ j)
|
---|
| 110 | std::cout << haltons[j] << " ";
|
---|
| 111 |
|
---|
| 112 | std::cout << std::endl;
|
---|
| 113 | }
|
---|
| 114 | }
|
---|
| 115 |
|
---|
| 116 |
|
---|
| 117 | HaltonSequence::HaltonSequence(int dim)
|
---|
| 118 | {
|
---|
| 119 | for (int i = 0; i < dim; ++ i)
|
---|
| 120 | {
|
---|
| 121 | mHaltons.push_back(Halton(i + 1));
|
---|
| 122 | }
|
---|
| 123 | }
|
---|
| 124 |
|
---|
| 125 |
|
---|
| 126 | void HaltonSequence::GetNext(float *numbers)
|
---|
| 127 | {
|
---|
| 128 | for (size_t i = 0; i < mHaltons.size(); ++ i)
|
---|
| 129 | {
|
---|
| 130 | numbers[i] = mHaltons[i].GetNext();
|
---|
| 131 | }
|
---|
| 132 | }
|
---|