Initial commit.

Project YACS (DES Simulator for an MM1 queue), linear congruential method, random walk for a Markov Chain and the Vector-Matrix Product code.

linear-congruential.c

@@ -0,0 +1,44 @@
+/*
+Code from Rosetta Code: https://rosettacode.org/wiki/Linear_congruential_generator#C
+*/
+#include <stdio.h>
+
+/* always assuming int is at least 32 bits */
+int rand();
+int rseed = 0;
+
+inline void srand(int x)
+{
+	rseed = x;
+}
+
+#ifndef MS_RAND
+#define RAND_MAX ((1U << 31) - 1)
+
+inline int rand()
+{
+	return rseed = (rseed * 1103515245 + 12345) & RAND_MAX;
+}
+
+#else /* MS rand */
+
+#define RAND_MAX_32 ((1U << 31) - 1)
+#define RAND_MAX ((1U << 15) - 1)
+
+inline int rand()
+{
+	return (rseed = (rseed * 214013 + 2531011) & RAND_MAX_32) >> 16;
+}
+
+#endif/* MS_RAND */
+
+int main()
+{
+	int i;
+	printf("rand max is %d\n", RAND_MAX);
+
+	for (i = 0; i < 100; i++)
+		printf("%d\n", rand());
+
+	return 0;
+}

random-walk.c

@@ -0,0 +1,58 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License version 2 as
+ * published by the Free Software Foundation;
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+ * Following code simulates a Discrete Time Markov Chain (DTMC). 
+ * Input is a DTMC (stochastic matrix).
+ * Output is a uniformised probability vector.
+ * Author: Ricardo M. Czekster ([email protected])
+ * Date: 25/11/2019
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <time.h>
+
+#define ORDER 3         // model size
+#define MAXRUNS 1000000 // number of runs
+
+int main(int argc, char *argv[]) {
+   float D[ORDER][ORDER] = {
+      { 0.470588235, 0.176470588, 0.352941176 },
+      { 0.470588235, 0.0        , 0.529411765 },
+      { 0.176470588, 0.117647059, 0.705882353 }
+   };
+   int i;
+   float r;             // variable for saving random number
+   int runs = 0;
+   int state = 0; // starts at 0 state (could be any -- selected at random)
+   int visits[ORDER];
+   float acc;
+   srand(time(NULL));            // set initial random seed -- based on now
+   for (i = 0; i < ORDER; i++)   // states initialisation
+      visits[i] = 0;
+   while (runs++ < MAXRUNS) {
+      r = (float) rand()/RAND_MAX;
+      acc = 0.0;
+      for (i = 0; i < ORDER; i++) {
+         acc += D[state][i];
+         if (r < acc) {
+            visits[i]++;
+            state = i;
+            break;
+         }
+      }
+   }
+   for (i = 0; i < ORDER; i++)
+      printf("%d=%f\n", i, (float)visits[i]/runs);
+}

vector-matrix-product.c

@@ -0,0 +1,76 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License version 2 as
+ * published by the Free Software Foundation;
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+ * Following code does a Vector-Matrix Product (VMP) for a Discrete Time Markov Chain (DTMC). 
+ * Input is a DTMC (stochastic matrix).
+ * Output is a uniformised probability vector.
+ * Author: Ricardo M. Czekster ([email protected])
+ * Date: 25/11/2019
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <time.h>
+#include <string.h>     // memcpy
+#include <math.h>       // pow, sqrt
+
+#define ORDER 3         // model size
+#define MAXRUNS 1000000 // number of runs
+#define RESIDUE 1e-10   // residual difference between two iterations
+
+/** Checks for convergence, i.e., all positions from two vectors must be greater than the RESIDUE constant. */
+int converge(float *r1, float *r2) {
+   int i;
+   for (i = 0; i < ORDER; i++)
+      if (sqrt(pow(r1[i] - r2[i], 2)) > RESIDUE)
+         return 0;
+   return 1;
+}
+
+/** Performs a vector-matrix product. */
+void multiply(float *v, float m[ORDER][ORDER]) {
+   int i,j;
+   float aux[ORDER];
+   for (i = 0; i < ORDER; i++) {
+      aux[i] = 0;
+      for (j = 0; j < ORDER; j++)
+         aux[i] += m[j][i] * v[j];
+   }
+   memcpy(v, aux, sizeof(float) * ORDER);
+}
+
+int main(int argc, char *argv[]) {
+   /*Change to *your* matrix here or build a text file reader containing the input DTMC matrix.*/
+   float D[ORDER][ORDER] = {
+      { 0.470588235, 0.176470588, 0.352941176 },
+      { 0.470588235, 0.0        , 0.529411765 },
+      { 0.176470588, 0.117647059, 0.705882353 }
+   };
+   int i;
+   int runs = 0;
+   float pvec[ORDER];
+   float old[ORDER];
+   pvec[0] = 1.0;  // initialise first position
+   for (i = 1; i < ORDER; i++)
+      pvec[i] = 0.0;
+   do {
+      memcpy(old, pvec, sizeof(float) * ORDER); // copy to old array
+      multiply(pvec, D);
+      if (++runs > MAXRUNS)
+         break;
+   } while (!converge(old, pvec));
+   printf("Number of iterations: %d\n", runs-1);
+   for (i = 0; i < ORDER; i++)
+      printf("%d=%f\n", i, pvec[i]);
+}