model.c

   1 /*
   2  * Copyright (c) 2017 Mohamed Aslan <maslan@sce.carleton.ca>
   3  *
   4  * Permission to use, copy, modify, and distribute this software for any
   5  * purpose with or without fee is hereby granted, provided that the above
   6  * copyright notice and this permission notice appear in all copies.
   7  *
   8  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
   9  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  10  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  11  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  12  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  13  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  14  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  15  */
  16
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <math.h>
  20 #include <stdint.h>
  21 #include "utils.h"
  22 #include "model.h"
  23
  24 static uint64_t
  25 nCr_r(uint32_t n, uint32_t r)
  26 {
  27         if (n == r)
  28                 return 1;
  29         if (r == 1)
  30                 return n;
  31         return nCr_r(n - 1, r) + nCr_r(n - 1, r - 1);
  32 }
  33
  34 static uint64_t
  35 nCr(uint32_t n, uint32_t r)
  36 {
  37         int              i, j;
  38         uint64_t        *C[101];
  39
  40         for (i = 0 ; i < 101 ; i++)
  41                 C[i] = (uint64_t *)xcalloc(101, sizeof(uint64_t));
  42         for (i = 1 ; i <= n ; i++) {
  43                 for (j = 1 ; j <= r ; j++) {
  44                         if (i == j)
  45                                 C[i][j] = 1;
  46                         else if (j == 1)
  47                                 C[i][j] = i;
  48                         else
  49                                 C[i][j] = C[i - 1][j] + C[i - 1][j - 1];
  50                 }
  51         }
  52         for (i = 0 ; i < 101 ; i++)
  53                 free(C[i]);
  54         return C[n][r];
  55 }
  56
  57 double
  58 rw2phi(int r, int w, int n)
  59 {
  60         if (r + w <= n)
  61                 return 1.0 - ((double)nCr(n - w, r) / (double)nCr(n, r));
  62         else
  63                 return 1.0;
  64 }
  65
  66 void
  67 phi2rw(double phi, int n, int *r, int *w)
  68 {
  69         int             i, i_min = 1, j, j_min = 1;
  70         double          d, d_min = 1.1;
  71
  72         for (i = 1 ; i < n ; i++) {
  73                 for (j = i ; j <= n - i ; j++) {
  74                         d = fabs(rw2phi(i, j, n) - phi);
  75                         if (d < d_min) {
  76                                 d_min = d;
  77                                 i_min = i;
  78                                 j_min = j;
  79                         }
  80                 }
  81         }
  82         *r = i_min;
  83         *w = j_min;
  84 }