dist.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4
   5
   6 static int loops;
   7
   8 struct rec {
   9         int d;
  10         int x;
  11 };
  12
  13
  14 static int comp(struct rec *r1, struct rec *r2)
  15 {
  16         return r1->x - r2->x;
  17 }
  18
  19 static int comp2(struct rec *r1, struct rec *r2)
  20 {
  21         return r1->d - r2->d;
  22 }
  23
  24 void shuffle(struct rec *d, int n)
  25 {
  26         int i;
  27
  28         for (i=0;i<n;i++)
  29                 d[i].x = random();
  30
  31         qsort(d, n, sizeof(*d), comp);
  32 }
  33
  34 int count(struct rec *d, int n, int x)
  35 {
  36         int i, ret=0;
  37
  38         for (i=0;i<n;i++)
  39                 if (d[i].d == x) ret++;
  40         return ret;
  41 }
  42
  43
  44 static double avg(struct rec *d, int n, int loops)
  45 {
  46         int sum=0;
  47         int i;
  48
  49         for (i=0;i<loops;i++) {
  50                 shuffle(d, 2*n);
  51                 sum += abs(n/2 - count(d, n, 1));
  52         }
  53         return sum/(1.0*i);
  54 }
  55
  56 static double countr(struct rec *d, int n, int loops, int r)
  57 {
  58         int sum=0;
  59         int i;
  60
  61         for (i=0;i<loops;i++) {
  62                 shuffle(d, 2*n);
  63                 if (count(d, n, 1) == r) sum++;
  64         }
  65         return sum/(1.*loops);
  66 }
  67
  68
  69 static double sqr(double x)
  70 {
  71         return x*x;
  72 }
  73
  74 static double fact(int n)
  75 {
  76         int i;
  77         double ret=1;
  78         for (i=2;i<=n;i++)
  79                 ret *= i;
  80         return ret;
  81 }
  82
  83 static double binomial(int n, int k)
  84 {
  85         return fact(n) / (fact(k) * fact(n-k));
  86 }
  87
  88 static double P0(int n)
  89 {
  90         return sqr(fact(n/2)) / fact(n);
  91 }
  92
  93 static double P1(int n)
  94 {
  95         return P0(n) * (n/2) * (n/2);
  96 }
  97
  98 static double P2(int n)
  99 {
 100         return P0(n) * (n/2) * (n/2 - 1)/2 * binomial(n/2, 2);
 101 }
 102
 103 static double Pn(int n, int r)
 104 {
 105         return sqr(binomial(n/2, r)) * sqr(fact(n/2)) / fact(n);
 106 }
 107
 108 static double En(int n)
 109 {
 110         return (n/2) * (sqr(fact(n/2))/fact(n)) * binomial(n-1, n/2-1);
 111 }
 112
 113
 114 /*
 115   > p;
 116
 117                                            2     2
 118                              binomial(n, i)  (n!)
 119                              ---------------------
 120                                      (2 n)!
 121
 122
 123 > 4*sum(p*(n/2 - i),i=0..n/2);
 124
 125                         2                       2            2
 126                     (n!)  binomial(n, 1 + 1/2 n)  (1 + 1/2 n)
 127                   2 ------------------------------------------
 128                                      (2 n)! n
 129
 130
 131 This is the expected number of elements in the wrong cell after the
 132 hypercube phase of a 4 way internal sort.
 133
 134  */
 135
 136
 137 /* take n random elements. Sample k of those elements. How far is the median
 138    of the n elements from the median of the k elements? Answer in terms
 139    of the k elements (so the answer is at most k/2) */
 140 static int median_deviation(int n, int k)
 141 {
 142         struct rec *d;
 143         int i, j;
 144         double s;
 145         double sum;
 146
 147         srandom(getpid() ^ time(NULL));
 148
 149         d = (struct rec *)malloc(sizeof(d[0])*n);
 150
 151         for (i=0;i<n;i++) {
 152                 d[i].d = i;
 153         }
 154
 155         sum = 0;
 156         for (j=0;j<loops;j++) {
 157                 shuffle(d, n);
 158
 159                 qsort(d, k, sizeof(*d), comp2);
 160
 161                 sum += abs(n/2 - d[k/2].d);
 162         }
 163         sum /= j;
 164
 165         free(d);
 166
 167         return sum;
 168 }
 169
 170
 171 int main(int argc, char *argv[])
 172 {
 173         struct rec *d;
 174         int i;
 175         int N = atoi(argv[1]);
 176         int k = atoi(argv[2]);
 177
 178         loops = atoi(argv[3]);
 179
 180         printf("median_deviation(N, k) = %g\n", median_deviation(N, k));
 181 }