Add Barnsley Fern in C (#832)

* Add Barnsley Fern in C

* Add comments and fix typos

Co-authored-by: James Schloss <jrs.schloss@gmail.com>

Author: mahdisarikhani

CONTRIBUTORS.md

@@ -56,4 +56,5 @@ This file lists everyone, who contributed to this repo and wanted to show up her
 - Ishaan Verma
 - Delphi1024
 - ntindle
-- Ridham177
+- Mahdi Sarikhani
+- Ridham177

contents/barnsley/barnsley.md

@@ -125,6 +125,8 @@ The biggest differences between the two code implementations is that the Barnsle
 {% method %}
 {% sample lang="jl" %}
 [import, lang:"julia"](code/julia/barnsley.jl)
+{% sample lang="c" %}
+[import, lang:"c"](code/c/barnsley.c)
 {% sample lang="java" %}
 [import, lang:"java"](code/java/Barnsley.java)
 {% endmethod %}

contents/barnsley/code/c/barnsley.c

@@ -0,0 +1,108 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+struct matrix {
+    double xx, xy, xz,
+           yx, yy, yz,
+           zx, zy, zz;
+};
+
+struct point2d {
+    double x, y;
+};
+
+struct point3d {
+    double x, y, z;
+};
+
+struct point3d matmul(struct matrix mat, struct point3d point)
+{
+    struct point3d out = {
+        mat.xx * point.x + mat.xy * point.y + mat.xz * point.z,
+        mat.yx * point.x + mat.yy * point.y + mat.yz * point.z,
+        mat.zx * point.x + mat.zy * point.y + mat.zz * point.z
+    };
+    return out;
+}
+
+// This function reads in the Hutchinson operator and corresponding
+// probabilities and returns a randomly selected transform
+// This works by choosing a random number and then iterating through all
+// probabilities until it finds an appropriate bin
+struct matrix select_array(struct matrix *hutchinson_op, double *probabilities,
+                           size_t num_op)
+{
+    // random number to be binned
+    double rnd = (double)rand() / RAND_MAX;
+
+    // This checks to see if a random number is in a bin, if not, that
+    // probability is subtracted from the random number and we check the next
+    // bin in the list
+    for (size_t i = 0; i < num_op; ++i) {
+        if (rnd < probabilities[i]) {
+            return hutchinson_op[i];
+        }
+        rnd -= probabilities[i];
+    }
+}
+
+// This is a general function to simulate a chaos game
+//  - output_points: pointer to an initialized output array
+//  - num: the number of iterations
+//  - initial_point: the starting point of the chaos game
+//  - hutchinson_op: the set of functions to iterate through
+//  - probabilities: the set of probabilities corresponding to the likelihood
+//      of choosing their corresponding function in hutchingson_op
+//  - nop: the number of functions in hutchinson_op
+void chaos_game(struct point2d *output_points, size_t num,
+                struct point2d initial_point, struct matrix *hutchinson_op,
+                double *probabilities, size_t nop)
+{
+    // extending point to 3D for affine transform
+    struct point3d point = {initial_point.x, initial_point.y, 1.0};
+
+    for (size_t i = 0; i < num; ++i) {
+        point = matmul(select_array(hutchinson_op, probabilities, nop), point);
+        output_points[i].x = point.x;
+        output_points[i].y = point.y;
+    }
+}
+
+int main()
+{
+    struct matrix barnsley_hutchinson[4] = {
+        {
+            0.0, 0.0, 0.0,
+            0.0, 0.16, 0.0,
+            0.0, 0.0, 1.0
+        },
+        {
+            0.85, 0.04, 0.0,
+            -0.04, 0.85, 1.60,
+            0.0, 0.0, 1.0
+        },
+        {
+            0.2, -0.26, 0.0,
+            0.23, 0.22, 1.60,
+            0.0, 0.0, 1.0
+        },
+        {
+            -0.15, 0.28, 0.0,
+            0.26, 0.24, 0.44,
+            0.0, 0.0, 1.0
+        }
+    };
+
+    double barnsley_probabilities[4] = {0.01, 0.85, 0.07, 0.07};
+    struct point2d output_points[10000];
+    struct point2d initial_point = {0.0, 0.0};
+    chaos_game(output_points, 10000, initial_point, barnsley_hutchinson,
+               barnsley_probabilities, 4);
+    FILE *f = fopen("barnsley.dat", "w");
+    for (size_t i = 0; i < 10000; ++i) {
+        fprintf(f, "%f\t%f\n", output_points[i].x, output_points[i].y);
+    }
+    fclose(f);
+
+    return 0;
+}