1 | #include "common.h"
|
---|
2 | #include "Halton.h"
|
---|
3 |
|
---|
4 |
|
---|
5 |
|
---|
6 | void Halton::TestHalton(int n, int dim)
|
---|
7 | {
|
---|
8 | Halton h(dim);
|
---|
9 |
|
---|
10 | for (int i = 1; i <= n; ++ i)
|
---|
11 | {
|
---|
12 | std::cout << "halton " << i << " of dim " << dim << " = " << h.GetNext() << std::endl;
|
---|
13 | }
|
---|
14 | }
|
---|
15 |
|
---|
16 |
|
---|
17 | void Halton::TestPrime()
|
---|
18 | {
|
---|
19 | int dim = 2;
|
---|
20 |
|
---|
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)
|
---|
29 | {
|
---|
30 | mBase = FindPrime(dim);
|
---|
31 | }
|
---|
32 |
|
---|
33 |
|
---|
34 | float Halton::GetNext()
|
---|
35 | {
|
---|
36 | ++ mIndex;
|
---|
37 |
|
---|
38 | float result = .0f;
|
---|
39 | float fraction = 1.0f / (float)mBase;
|
---|
40 | int idx = mIndex;
|
---|
41 |
|
---|
42 | while (idx > 0)
|
---|
43 | {
|
---|
44 | int digit = idx % mBase;
|
---|
45 | result += fraction * (float)digit;
|
---|
46 |
|
---|
47 | idx = (idx - digit) / mBase;
|
---|
48 | fraction /= (float)mBase;
|
---|
49 | }
|
---|
50 |
|
---|
51 | //mIndex %= 100000000;
|
---|
52 | return result;
|
---|
53 | }
|
---|
54 |
|
---|
55 |
|
---|
56 | bool Halton::IsPrime(int n)
|
---|
57 | {
|
---|
58 | bool isPrime = true;
|
---|
59 |
|
---|
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;
|
---|
70 | }
|
---|
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 | }
|
---|