GopherProxy

	trandcounts.c - numtools - perform numerical operations on vectors and matrices…
	git clone git://src.adamsgaard.dk/numtools
	Log
	Files
	Refs
	README
	LICENSE
	---
	trandcounts.c (2893B)
	---
	1 #include <stdio.h>
	2 #include <stdlib.h>
	3 #include <err.h>
	4 #include <limits.h>
	5 #include <string.h>
	6 #include <math.h>
	7 #include <sys/time.h>
	8 #include <unistd.h>
	9
	10 #include "util.h"
	11
	12 static void
	13 usage(void)
	14 {
	15 errx(1, "usage: randcounts [-d delimstr] [-n] [-N num] [-p prec]"
	16 "[-r repeats] [-R] [-s seed] weight1 [weight2 ...]");
	17 }
	18
	19 int
	20 main(int argc, char *argv[])
	21 {
	22 int i, ch, nbins, prec = 17, finalnl = 1, s = 0;
	23 long j, seed, *counts = NULL, n = 1, r, repeats = 1, ratio = 0;
	24 double val = 0.0, weightsum = 0.0, *weights = NULL;
	25 char *delimstr = "\t";
	26 const char *errstr;
	27 struct timeval t1;
	28
	29 if (pledge("stdio", NULL) == -1)
	30 err(2, "pledge");
	31
	32 while ((ch = getopt(argc, argv, "d:nN:r:Rp:s:")) != -1) {
	33 switch (ch) {
	34 case 'd':
	35 delimstr = optarg;
	36 break;
	37 case 'n':
	38 finalnl = 0;
	39 break;
	40 case 'N':
	41 n = strtonum(optarg, 0, LONG_MAX, &errstr);
	42 if (errstr != NULL)
	43 errx(1, "bad num value, %s: %s", errstr,…
	44 break;
	45 case 'r':
	46 repeats = strtonum(optarg, 0, LONG_MAX, &errstr);
	47 if (errstr != NULL)
	48 errx(1, "bad repeats value, %s: %s", err…
	49 break;
	50 case 'R':
	51 ratio = 1;
	52 break;
	53 case 'p':
	54 prec = strtonum(optarg, 0, INT_MAX, &errstr);
	55 if (errstr != NULL)
	56 errx(1, "bad precision value, %s: %s", e…
	57 break;
	58 case 's':
	59 seed = strtonum(optarg, 0, INT_MAX, &errstr);
	60 if (errstr != NULL)
	61 errx(1, "bad seed value, %s: %s", errstr…
	62 break;
	63 default:
	64 usage();
	65 }
	66 }
	67 argc -= optind;
	68 argv += optind;
	69 if (argc < 1)
	70 usage();
	71
	72 nbins = argc;
	73 if (!(weights = calloc(nbins, sizeof(double))) \|\|
	74 !(counts = calloc(nbins, sizeof(long))))
	75 err(1, "calloc");
	76
	77 for (i = 0; i < nbins; i++) {
	78 if (!sscanf(argv[i], "%lf", &weights[i]))
	79 errx(1, "bad weight value: %s", argv[i]);
	80 if (weights[i] <= 0.0)
	81 errx(1, "weight %d is not positive (%g)", i, wei…
	82 if (weights[i] > 1.0)
	83 errx(1, "weight %d is greater than 1 (%g)", i, w…
	84 }
	85
	86 for (i = 0; i < nbins; i++)
	87 weightsum += weights[i];
	88 if (fabs(weightsum - 1.0) > 1e-3)
	89 errx(1, "weights do not sum to 1 (%g)", weightsum);
	90
	91 if (s)
	92 #ifdef __OpenBSD__
	93 srand48_deterministic(seed);
	94 #else
	95 srand48(seed);
	96 else {
	97 gettimeofday(&t1, NULL);
	98 srand48(t1.tv_sec * t1.tv_usec); /* Caveat: identical se…
	99 }
	100 #endif
	101
	102 for (r = 0; r < repeats; r++) {
	103 for (i = 0; i < nbins; i++)
	104 counts[i] = 0;
	105 for (j = 0; j < n; j++) {
	106 val = drand48();
	107 weightsum = 0.0;
	108 for (i = 0; i < nbins; i++) {
	109 weightsum += weights[i];
	110 if (val <= weightsum) {
	111 counts[i]++;
	112 break;
	113 }
	114 }
	115 }
	116 for (i = 0; i < nbins; i++) {
	117 if (ratio)
	118 printf("%.*g", prec, (double)counts[i] /…
	119 else
	120 printf("%ld", counts[i]);
	121 if (i < nbins - 1)
	122 fputs(delimstr, stdout);
	123 else
	124 if (finalnl)
	125 putchar('\n');
	126 }
	127 }
	128
	129 free(counts);
	130 free(weights);
	131
	132 return 0;
	133 }