src/cal/cnt.c

#ifndef lint
static const char       RCSid[] = "$Id: cnt.c,v 1.5 2022/04/20 20:56:01 greg Exp $";
#endif
/*
 *  cnt.c - simple counting program.
 *
 *      2/1/88
 *
 *  Added -s (shuffle) option April 2022
 */

#include  <stdio.h>
#include  <time.h>
#include  "random.h"

#ifndef uint16
#define uint16  unsigned short          /* 16-bit unsigned integer */
#endif
#undef uby8
#define uby8  unsigned char             /* 8-bit unsigned integer */

#define MAXDIM  50

#define NLEVELS 9                       /* number of tree levels */
#define BRORDER 6                       /* branches/level */

/* Tree branch structure for quick occupancy search */
/* with 9 levels & 6 branches per level, we can store 1.94 Gbits in 259 MBytes (4.5% overhead) */
const struct {
        long    capacity;               /* slots/branch this level */
        long    skip_bytes;             /* bytes until next branch */
        int     cntr_siz;               /* occupancy counter size */
} tree_br[NLEVELS] = {
        {248L, 32L, 1},
        {248L*6, 32L*6+2, 2},
        {248L*6*6, (32L*6+2)*6+2, 2},
        {248L*6*6*6, ((32L*6+2)*6+2)*6+2, 2},
        {248L*6*6*6*6, (((32L*6+2)*6+2)*6+2)*6+3, 3},
        {248L*6*6*6*6*6, ((((32L*6+2)*6+2)*6+2)*6+3)*6+3, 3},
        {248L*6*6*6*6*6*6, (((((32L*6+2)*6+2)*6+2)*6+3)*6+3)*6+3, 3},
        {248L*6*6*6*6*6*6*6, ((((((32L*6+2)*6+2)*6+2)*6+3)*6+3)*6+3)*6+4, 4},
        {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},
};

char  buf[256];                         /* buffer for ordered array output */


/* Encode integer in string and return pointer to end */
static char *
tack(char *b, long i)
{
        char  *cp;
        char  *res;

        *b++ = '\t';
        cp = b;
        if (i == 0)
                *cp++ = '0';
        else
                do {
                        *cp++ = i%10L + '0';
                        i /= 10L;
                } while (i);
        res = cp--;
#define c i
        while (cp > b) {                /* reverse string */
                c = *cp;
                *cp-- = *b;
                *b++ = c;
        }
#undef c
        return(res);
}


/* Loop over dimensions, spitting out buffer after each increment */
static void
loop(long *n, char *b)
{
        long  i;

        if (n[0] == 0) {
                *b = '\0';
                puts(buf);
                return;
        }
        for (i = 0; i < n[0]; i++)
                loop(n+1, tack(b, i));
}


/* Print out shuffled value */
static void
print_shuf(long *n, long aval)
{
        int     i;

        for (i = 0; n[i+1]; i++) {
                printf("\t%ld", aval % n[i]);
                aval /= n[i];
        }
        printf("\t%ld\n", aval);
}


/* Allocate and prepare occupancy tree */
static uby8 *
tree_alloc(long alen)
{
        uby8  *troot;
        double  bytes_per_bit;
        int  i;
        int  ht = 0;
                                        /* how tall does our tree need to be? */
        while (tree_br[ht].capacity*BRORDER < alen)
                if (++ht >= NLEVELS) {
                        fputs("Array too large to shuffle\n", stderr);
                        exit(1);
                }
        bytes_per_bit = 1.;             /* figure out tree size (with overhead) */
        for (i = ht; i >= 0; i--)
                bytes_per_bit += (double)tree_br[i].cntr_siz;
        bytes_per_bit += (double)tree_br[ht].skip_bytes;
        bytes_per_bit /= (double)tree_br[ht].capacity;
        troot = (uby8 *)calloc((long)(alen*bytes_per_bit)+2, 1);
        if (troot == NULL) {
                fputs("Not enough memory for shuffle\n", stderr);
                exit(1);
        }
        *troot = ht;                    /* first byte is tree height */
        for (i = 256; i--; ) {          /* assign 0-bit count table */
                int     b;
                buf[i] = 8;
                for (b = i; b; b >>= 1)
                        buf[i] -= b&1;
        }
        return(troot);
}


/* Get number of slots available at this branch location */
static long
get_avail(const uby8 *ctrp, int lvl)
{
        long    cnt = 0;
        int     n = tree_br[lvl].cntr_siz;

        while (--n > 0) {               /* LSB first */
                cnt |= ctrp[n];
                cnt <<= 8;
        }
        cnt |= ctrp[0];

        return(tree_br[lvl].capacity - cnt);
}


/* Increment branch occupancy counter */
static void
incr_counter(uby8 *ctrp, int n)
{
        n = tree_br[n].cntr_siz;

        while (n-- > 0)                 /* LSB first */
                if (++(*ctrp++))
                        break;
}


/* Skip to and allocate a leaf from tree */
static long
eat_nth_leaf(uby8 *brp, long ski)
{
        int     lvl = *brp++;                   /* tree height in first byte */
        long    pos = 0;
        int     b;

        while (lvl >= 0) {                      /* descend to leaves */
                long  navail;
                b = 0;                          /* select each branch */
                while (ski >= (navail = get_avail(brp, lvl))) {
                        if (++b >= BRORDER) {
                                fputs("Shuffle tree error!\n", stderr);
                                exit(1);
                        }
                        pos += tree_br[lvl].capacity;
                        ski -= navail;
                        brp += tree_br[lvl].skip_bytes;
                }
                incr_counter(brp, lvl);         /* we intend to eat one */
                brp += tree_br[lvl--].cntr_siz; /* drop a level */
        }
        while (ski >= buf[*brp]) {              /* browse the leaves */
                pos += 8;
                ski -= buf[*brp++];             /* buf contains 0-bit counts */
        }
        b = 0;                                  /* find target bit in byte */
        while ((ski -= !(*brp & 1<<b)) >= 0) {
                pos++;
                b++;
        }
        *brp |= 1<<b;                           /* eat it */
        return(pos);                            /* & return leaf's slot# */
}


/* Shuffle all possible output strings and spit out randomly (tree version) */
static void
big_shuffle(long *n, long alen)
{
        uby8  *tree_root;
                                        /* size and allocate holder tree */
        tree_root = tree_alloc(alen);

        while (alen > 0)                /* allocate and print random array entries */
                print_shuf(n, eat_nth_leaf(tree_root, irandom(alen--)));

        free(tree_root);                /* all done */
}


/* Shuffle all possible output strings and spit out randomly */
static void
shuffle(long *n)
{
        long  alen;
        uint16  *myshuf;
        int  i;

        alen = 1;               /* compute shuffle size */
        for (i = 0; n[i]; i++) {
                if (alen*n[i] <= alen) {
                        fputs("Array too large to count!\n", stderr);
                        exit(1);
                }
                alen *= n[i];
        }
                                /* get unique starting point */
        srandom((long)time(0));

        if (alen > 1L<<16) {    /* use large shuffle method? */
                big_shuffle(n, alen);
                return;
        }
        myshuf = (uint16 *)malloc(alen*sizeof(uint16));
        if (myshuf == NULL) {
                fputs("Insufficient memory for shuffle\n", stderr);
                exit(1);
        }
        for (i = alen; i--; )   /* initialize in any order */
                myshuf[i] = i;
                                /* perform Fisher-Yates shuffle */
        for (i = 0; i < alen-1; i++) {
                int     ix = irandom(alen-i) + i;
                int     ndx = myshuf[i];
                myshuf[i] = myshuf[ix];
                myshuf[ix] = ndx;
        }
                                /* put randomly indexed output */
        for (i = alen; i--; )
                print_shuf(n, (long)myshuf[i]);

        free(myshuf);           /* all done */
}


int
main(int argc, char *argv[])
{
        char  *prog = argv[0];
        int  doshuffle = 0;
        long  n[MAXDIM];
        int  a;

        argv++; argc--;
        if (argc <= 0)
                goto userr;
        if (argv[0][0] == '-' && argv[0][1] == 's') {
                doshuffle = 1;
                argv++; argc--;
        }
        for (a = 0; a < argc; a++)
                if ((n[a] = atol(argv[a])) <= 1)
                        goto userr;
        n[a] = 0;
        if (!a)
                goto userr;

        if (doshuffle)
                shuffle(n);
        else
                loop(n, buf);

        return(0);
userr:
        fputs("Usage: ", stderr);
        fputs(prog, stderr);
        fputs(" [-s] N0 [N1 ..]\n", stderr);
        return(1);
}
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
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	1.6	static const char RCSid[] = "$Id: cnt.c,v 1.5 2022/04/20 20:56:01 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* cnt.c - simple counting program.
6			*
7			* 2/1/88
8	greg	1.3	*
9			* Added -s (shuffle) option April 2022
10	greg	1.1	*/
11
12			#include <stdio.h>
13	greg	1.3	#include <time.h>
14			#include "random.h"
15	greg	1.1
16	greg	1.4	#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	greg	1.3	#define MAXDIM 50
23	greg	1.1
24	greg	1.4	#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			{248L6, 32L6+2, 2},
36			{248L66, (32L6+2)6+2, 2},
37			{248L666, ((32L6+2)6+2)6+2, 2},
38			{248L6666, (((32L6+2)6+2)6+2)6+3, 3},
39			{248L66666, ((((32L6+2)6+2)6+2)6+3)6+3, 3},
40			{248L666666, (((((32L6+2)6+2)6+2)6+3)6+3)6+3, 3},
41			{248L6666666, ((((((32L6+2)6+2)6+2)6+3)6+3)6+3)6+4, 4},
42			{248L66666666, (((((((32L6+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	greg	1.1
47
48	greg	1.3	/* Encode integer in string and return pointer to end */
49			static char *
50	greg	1.4	tack(char *b, long i)
51	greg	1.1	{
52	greg	1.3	char *cp;
53	greg	1.1	char *res;
54
55			*b++ = '\t';
56			cp = b;
57			if (i == 0)
58			*cp++ = '0';
59			else
60			do {
61	greg	1.4	*cp++ = i%10L + '0';
62			i /= 10L;
63	greg	1.1	} 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	greg	1.3	/* Loop over dimensions, spitting out buffer after each increment */
77	schorsch	1.2	static void
78	greg	1.4	loop(long n, char b)
79	greg	1.1	{
80	greg	1.4	long i;
81	greg	1.1
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	greg	1.3
91
92	greg	1.4	/* 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)(alenbytes_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	greg	1.5	for (i = 256; i--; ) { /* assign 0-bit count table */
132	greg	1.4	int b;
133	greg	1.5	buf[i] = 8;
134	greg	1.4	for (b = i; b; b >>= 1)
135	greg	1.5	buf[i] -= b&1;
136	greg	1.4	}
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	greg	1.5	while (ski >= buf[brp]) { / browse the leaves */
194	greg	1.4	pos += 8;
195	greg	1.5	ski -= buf[brp++]; / buf contains 0-bit counts */
196	greg	1.4	}
197	greg	1.5	b = 0; /* find target bit in byte */
198			while ((ski -= !(*brp & 1<<b)) >= 0) {
199	greg	1.4	pos++;
200	greg	1.5	b++;
201	greg	1.4	}
202	greg	1.5	brp \|= 1<<b; / eat it */
203			return(pos); /* & return leaf's slot# */
204	greg	1.4	}
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	greg	1.6	print_shuf(n, eat_nth_leaf(tree_root, irandom(alen--)));
217	greg	1.4
218			free(tree_root); /* all done */
219			}
220
221
222	greg	1.3	/* Shuffle all possible output strings and spit out randomly */
223			static void
224	greg	1.4	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	greg	1.3	exit(1);
235	greg	1.4	}
236			alen *= n[i];
237			}
238			/* get unique starting point */
239			srandom((long)time(0));
240	greg	1.3
241	greg	1.4	if (alen > 1L<<16) { /* use large shuffle method? */
242			big_shuffle(n, alen);
243			return;
244			}
245			myshuf = (uint16 )malloc(alensizeof(uint16));
246	greg	1.3	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	greg	1.6	int ix = irandom(alen-i) + i;
255	greg	1.3	int ndx = myshuf[i];
256			myshuf[i] = myshuf[ix];
257			myshuf[ix] = ndx;
258			}
259			/* put randomly indexed output */
260	greg	1.4	for (i = alen; i--; )
261			print_shuf(n, (long)myshuf[i]);
262
263			free(myshuf); /* all done */
264	greg	1.3	}
265
266
267			int
268	greg	1.4	main(int argc, char *argv[])
269	greg	1.3	{
270			char *prog = argv[0];
271			int doshuffle = 0;
272	greg	1.4	long n[MAXDIM];
273	greg	1.3	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	greg	1.4	if ((n[a] = atol(argv[a])) <= 1)
284	greg	1.3	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			}