ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cal/cnt.c
Revision: 1.6
Committed: Thu Apr 21 02:52:40 2022 UTC (2 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.5: +3 -3 lines
Log Message:
perf(cnt): improved randomness of cnt -s option with new irandom(modulus) call

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: cnt.c,v 1.5 2022/04/20 20:56:01 greg Exp $";
3 #endif
4 /*
5 * cnt.c - simple counting program.
6 *
7 * 2/1/88
8 *
9 * Added -s (shuffle) option April 2022
10 */
11
12 #include <stdio.h>
13 #include <time.h>
14 #include "random.h"
15
16 #ifndef uint16
17 #define uint16 unsigned short /* 16-bit unsigned integer */
18 #endif
19 #undef uby8
20 #define uby8 unsigned char /* 8-bit unsigned integer */
21
22 #define MAXDIM 50
23
24 #define NLEVELS 9 /* number of tree levels */
25 #define BRORDER 6 /* branches/level */
26
27 /* Tree branch structure for quick occupancy search */
28 /* with 9 levels & 6 branches per level, we can store 1.94 Gbits in 259 MBytes (4.5% overhead) */
29 const struct {
30 long capacity; /* slots/branch this level */
31 long skip_bytes; /* bytes until next branch */
32 int cntr_siz; /* occupancy counter size */
33 } tree_br[NLEVELS] = {
34 {248L, 32L, 1},
35 {248L*6, 32L*6+2, 2},
36 {248L*6*6, (32L*6+2)*6+2, 2},
37 {248L*6*6*6, ((32L*6+2)*6+2)*6+2, 2},
38 {248L*6*6*6*6, (((32L*6+2)*6+2)*6+2)*6+3, 3},
39 {248L*6*6*6*6*6, ((((32L*6+2)*6+2)*6+2)*6+3)*6+3, 3},
40 {248L*6*6*6*6*6*6, (((((32L*6+2)*6+2)*6+2)*6+3)*6+3)*6+3, 3},
41 {248L*6*6*6*6*6*6*6, ((((((32L*6+2)*6+2)*6+2)*6+3)*6+3)*6+3)*6+4, 4},
42 {248L*6*6*6*6*6*6*6*6, (((((((32L*6+2)*6+2)*6+2)*6+3)*6+3)*6+3)*6+4)*6+4, 4},
43 };
44
45 char buf[256]; /* buffer for ordered array output */
46
47
48 /* Encode integer in string and return pointer to end */
49 static char *
50 tack(char *b, long i)
51 {
52 char *cp;
53 char *res;
54
55 *b++ = '\t';
56 cp = b;
57 if (i == 0)
58 *cp++ = '0';
59 else
60 do {
61 *cp++ = i%10L + '0';
62 i /= 10L;
63 } while (i);
64 res = cp--;
65 #define c i
66 while (cp > b) { /* reverse string */
67 c = *cp;
68 *cp-- = *b;
69 *b++ = c;
70 }
71 #undef c
72 return(res);
73 }
74
75
76 /* Loop over dimensions, spitting out buffer after each increment */
77 static void
78 loop(long *n, char *b)
79 {
80 long i;
81
82 if (n[0] == 0) {
83 *b = '\0';
84 puts(buf);
85 return;
86 }
87 for (i = 0; i < n[0]; i++)
88 loop(n+1, tack(b, i));
89 }
90
91
92 /* Print out shuffled value */
93 static void
94 print_shuf(long *n, long aval)
95 {
96 int i;
97
98 for (i = 0; n[i+1]; i++) {
99 printf("\t%ld", aval % n[i]);
100 aval /= n[i];
101 }
102 printf("\t%ld\n", aval);
103 }
104
105
106 /* Allocate and prepare occupancy tree */
107 static uby8 *
108 tree_alloc(long alen)
109 {
110 uby8 *troot;
111 double bytes_per_bit;
112 int i;
113 int ht = 0;
114 /* how tall does our tree need to be? */
115 while (tree_br[ht].capacity*BRORDER < alen)
116 if (++ht >= NLEVELS) {
117 fputs("Array too large to shuffle\n", stderr);
118 exit(1);
119 }
120 bytes_per_bit = 1.; /* figure out tree size (with overhead) */
121 for (i = ht; i >= 0; i--)
122 bytes_per_bit += (double)tree_br[i].cntr_siz;
123 bytes_per_bit += (double)tree_br[ht].skip_bytes;
124 bytes_per_bit /= (double)tree_br[ht].capacity;
125 troot = (uby8 *)calloc((long)(alen*bytes_per_bit)+2, 1);
126 if (troot == NULL) {
127 fputs("Not enough memory for shuffle\n", stderr);
128 exit(1);
129 }
130 *troot = ht; /* first byte is tree height */
131 for (i = 256; i--; ) { /* assign 0-bit count table */
132 int b;
133 buf[i] = 8;
134 for (b = i; b; b >>= 1)
135 buf[i] -= b&1;
136 }
137 return(troot);
138 }
139
140
141 /* Get number of slots available at this branch location */
142 static long
143 get_avail(const uby8 *ctrp, int lvl)
144 {
145 long cnt = 0;
146 int n = tree_br[lvl].cntr_siz;
147
148 while (--n > 0) { /* LSB first */
149 cnt |= ctrp[n];
150 cnt <<= 8;
151 }
152 cnt |= ctrp[0];
153
154 return(tree_br[lvl].capacity - cnt);
155 }
156
157
158 /* Increment branch occupancy counter */
159 static void
160 incr_counter(uby8 *ctrp, int n)
161 {
162 n = tree_br[n].cntr_siz;
163
164 while (n-- > 0) /* LSB first */
165 if (++(*ctrp++))
166 break;
167 }
168
169
170 /* Skip to and allocate a leaf from tree */
171 static long
172 eat_nth_leaf(uby8 *brp, long ski)
173 {
174 int lvl = *brp++; /* tree height in first byte */
175 long pos = 0;
176 int b;
177
178 while (lvl >= 0) { /* descend to leaves */
179 long navail;
180 b = 0; /* select each branch */
181 while (ski >= (navail = get_avail(brp, lvl))) {
182 if (++b >= BRORDER) {
183 fputs("Shuffle tree error!\n", stderr);
184 exit(1);
185 }
186 pos += tree_br[lvl].capacity;
187 ski -= navail;
188 brp += tree_br[lvl].skip_bytes;
189 }
190 incr_counter(brp, lvl); /* we intend to eat one */
191 brp += tree_br[lvl--].cntr_siz; /* drop a level */
192 }
193 while (ski >= buf[*brp]) { /* browse the leaves */
194 pos += 8;
195 ski -= buf[*brp++]; /* buf contains 0-bit counts */
196 }
197 b = 0; /* find target bit in byte */
198 while ((ski -= !(*brp & 1<<b)) >= 0) {
199 pos++;
200 b++;
201 }
202 *brp |= 1<<b; /* eat it */
203 return(pos); /* & return leaf's slot# */
204 }
205
206
207 /* Shuffle all possible output strings and spit out randomly (tree version) */
208 static void
209 big_shuffle(long *n, long alen)
210 {
211 uby8 *tree_root;
212 /* size and allocate holder tree */
213 tree_root = tree_alloc(alen);
214
215 while (alen > 0) /* allocate and print random array entries */
216 print_shuf(n, eat_nth_leaf(tree_root, irandom(alen--)));
217
218 free(tree_root); /* all done */
219 }
220
221
222 /* Shuffle all possible output strings and spit out randomly */
223 static void
224 shuffle(long *n)
225 {
226 long alen;
227 uint16 *myshuf;
228 int i;
229
230 alen = 1; /* compute shuffle size */
231 for (i = 0; n[i]; i++) {
232 if (alen*n[i] <= alen) {
233 fputs("Array too large to count!\n", stderr);
234 exit(1);
235 }
236 alen *= n[i];
237 }
238 /* get unique starting point */
239 srandom((long)time(0));
240
241 if (alen > 1L<<16) { /* use large shuffle method? */
242 big_shuffle(n, alen);
243 return;
244 }
245 myshuf = (uint16 *)malloc(alen*sizeof(uint16));
246 if (myshuf == NULL) {
247 fputs("Insufficient memory for shuffle\n", stderr);
248 exit(1);
249 }
250 for (i = alen; i--; ) /* initialize in any order */
251 myshuf[i] = i;
252 /* perform Fisher-Yates shuffle */
253 for (i = 0; i < alen-1; i++) {
254 int ix = irandom(alen-i) + i;
255 int ndx = myshuf[i];
256 myshuf[i] = myshuf[ix];
257 myshuf[ix] = ndx;
258 }
259 /* put randomly indexed output */
260 for (i = alen; i--; )
261 print_shuf(n, (long)myshuf[i]);
262
263 free(myshuf); /* all done */
264 }
265
266
267 int
268 main(int argc, char *argv[])
269 {
270 char *prog = argv[0];
271 int doshuffle = 0;
272 long n[MAXDIM];
273 int a;
274
275 argv++; argc--;
276 if (argc <= 0)
277 goto userr;
278 if (argv[0][0] == '-' && argv[0][1] == 's') {
279 doshuffle = 1;
280 argv++; argc--;
281 }
282 for (a = 0; a < argc; a++)
283 if ((n[a] = atol(argv[a])) <= 1)
284 goto userr;
285 n[a] = 0;
286 if (!a)
287 goto userr;
288
289 if (doshuffle)
290 shuffle(n);
291 else
292 loop(n, buf);
293
294 return(0);
295 userr:
296 fputs("Usage: ", stderr);
297 fputs(prog, stderr);
298 fputs(" [-s] N0 [N1 ..]\n", stderr);
299 return(1);
300 }