ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cal/cnt.c
(Generate patch)

Comparing ray/src/cal/cnt.c (file contents):
Revision 1.3 by greg, Mon Apr 11 18:08:19 2022 UTC vs.
Revision 1.4 by greg, Wed Apr 20 19:13:12 2022 UTC

# Line 13 | Line 13 | static const char      RCSid[] = "$Id$";
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 < char  buf[256];
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(
24 < char  *b,
25 < int  i
26 < )
50 > tack(char *b, long i)
51   {
52          char  *cp;
53          char  *res;
# Line 34 | Line 58 | int  i
58                  *cp++ = '0';
59          else
60                  do {
61 <                        *cp++ = i%10 + '0';
62 <                        i /= 10;
61 >                        *cp++ = i%10L + '0';
62 >                        i /= 10L;
63                  } while (i);
64          res = cp--;
65   #define c i
# Line 51 | Line 75 | int  i
75  
76   /* Loop over dimensions, spitting out buffer after each increment */
77   static void
78 < loop(
55 < int  *n,
56 < char  *b
57 < )
78 > loop(long *n, char *b)
79   {
80 <        int  i;
80 >        long  i;
81  
82          if (n[0] == 0) {
83                  *b = '\0';
# Line 68 | Line 89 | char  *b
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 bit count table */
132 +                int     b;
133 +                buf[i] = 0;
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 +        b = 0;                                  /* browse the leaves */
194 +        while (ski >= 8-buf[*brp]) {            /* buf contains bit counts */
195 +                pos += 8;
196 +                ski -= 8-buf[*brp++];
197 +        }
198 +        while (ski >= 0) {                      /* target zero leaf bit */
199 +                pos++;
200 +                ski -= !(*brp & 1<<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, random() % 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(
74 < int *n
75 < )
224 > shuffle(long *n)
225   {
226 <        int  sub[MAXDIM];
227 <        int  ndim;
228 <        int  alen;
80 <        int  *myshuf;
81 <        int  i, j;
226 >        long  alen;
227 >        uint16  *myshuf;
228 >        int  i;
229  
230 <        alen = 1;               /* allocate shuffle index array */
231 <        for (j = 0; n[j]; j++)
232 <                if ((alen *= n[j]) < 0)
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 <        myshuf = (int *)malloc(alen*sizeof(int));
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;
95 –                                /* get unique starting point */
96 –        srandom((long)time(0));
252                                  /* perform Fisher-Yates shuffle */
253          for (i = 0; i < alen-1; i++) {
254                  int     ix = random()%(alen-i) + i;
# Line 102 | Line 257 | int *n
257                  myshuf[ix] = ndx;
258          }
259                                  /* put randomly indexed output */
260 <        for (i = alen; i--; ) {
261 <                int     aval = myshuf[i];
262 <                for (j = 0; n[j+1]; j++) {
263 <                        printf("\t%d", aval % n[j]);
109 <                        aval /= n[j];
110 <                }
111 <                printf("\t%d\n", aval);
112 <        }
113 <        free(myshuf);
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(
119 < int  argc,
120 < char  *argv[]
121 < )
268 > main(int argc, char *argv[])
269   {
270          char  *prog = argv[0];
271          int  doshuffle = 0;
272 <        int  n[MAXDIM];
272 >        long  n[MAXDIM];
273          int  a;
274  
275          argv++; argc--;
# Line 133 | Line 280 | char  *argv[]
280                  argv++; argc--;
281          }
282          for (a = 0; a < argc; a++)
283 <                if ((n[a] = atoi(argv[a])) <= 1)
283 >                if ((n[a] = atol(argv[a])) <= 1)
284                          goto userr;
285          n[a] = 0;
286          if (!a)

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines