randcounts.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, optarg); 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", errstr, optarg); 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", errstr, optarg); 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, optarg); 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, weights[i]); 82 if (weights[i] > 1.0) 83 errx(1, "weight %d is greater than 1 (%g)", i, weights[i]); 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 seed for each ms */ 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] / n); 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 }