LXR iproute2/netem/stats.c

   1/* SPDX-License-Identifier: NIST-PD */
   2/*
   3 * Experimental data  distribution table generator
   4 * Taken from the uncopyrighted NISTnet code (public domain).
   5 *
   6 * Rread in a series of "random" data values, either
   7 * experimentally or generated from some probability distribution.
   8 * From this, report statistics.
   9 */
  10
  11#include <stdio.h>
  12#include <stdlib.h>
  13#include <math.h>
  14#include <malloc.h>
  15#include <sys/types.h>
  16#include <sys/stat.h>
  17
  18void
  19stats(FILE *fp)
  20{
  21        struct stat info;
  22        double *x;
  23        int limit;
  24        int n=0, i;
  25        double mu=0.0, sigma=0.0, sumsquare=0.0, sum=0.0, top=0.0, rho=0.0;
  26        double sigma2=0.0;
  27
  28        fstat(fileno(fp), &info);
  29        if (info.st_size > 0) {
  30                limit = 2*info.st_size/sizeof(double);  /* @@ approximate */
  31        } else {
  32                limit = 10000;
  33        }
  34        x = (double *)malloc(limit*sizeof(double));
  35
  36        for (i=0; i<limit; ++i){
  37                fscanf(fp, "%lf", &x[i]);
  38                if (feof(fp))
  39                        break;
  40                sumsquare += x[i]*x[i];
  41                sum += x[i];
  42                ++n;
  43        }
  44        mu = sum/(double)n;
  45        sigma = sqrt((sumsquare - (double)n*mu*mu)/(double)(n-1));
  46
  47        for (i=1; i < n; ++i){
  48                top += ((double)x[i]-mu)*((double)x[i-1]-mu);
  49                sigma2 += ((double)x[i-1] - mu)*((double)x[i-1] - mu);
  50
  51        }
  52        rho = top/sigma2;
  53
  54        printf("mu =    %12.6f\n", mu);
  55        printf("sigma = %12.6f\n", sigma);
  56        printf("rho =   %12.6f\n", rho);
  57        /*printf("sigma2 = %10.4f\n", sqrt(sigma2/(double)(n-1)));*/
  58        /*printf("correlation rho = %10.6f\n", top/((double)(n-1)*sigma*sigma));*/
  59}
  60
  61
  62int
  63main(int argc, char **argv)
  64{
  65        FILE *fp;
  66
  67        if (argc > 1) {
  68                fp = fopen(argv[1], "r");
  69                if (!fp) {
  70                        perror(argv[1]);
  71                        exit(1);
  72                }
  73        } else {
  74                fp = stdin;
  75        }
  76        stats(fp);
  77        return 0;
  78}
  79