Mandelbrot

mandelbrot.c

// Author: Will Killian
//
// Modified/Generalized from:
//   https://code.it4i.cz/jansik/mandelbrot/-/blob/master/mandelbrot-real-fma-omp.c

#include <stdio.h>
#include <stdlib.h>

#define SIZE 512

static const double input[SIZE] __attribute__((aligned(64))) =
{ -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
  -2.00000, -1.92742, -1.85484,   -1.78226, -1.70968, -1.63710, -1.56452,
  -1.49194, -1.41935, -1.34677,   -1.27419, -1.20161, -1.12903, -1.05645,
  -0.98387, -0.91129, -0.83871,   -0.76613, -0.69355, -0.62097, -0.54839,
  -0.47581, -0.40323, -0.33065,   -0.25806, -0.18548, -0.11290, -0.04032,
   0.03226,  0.10484,  0.17742,    0.25000,
};

static double output[SIZE] __attribute__((aligned(64)));

void mandelbrot (int iters, int size) {
  for (int i = 0; i < size; ++i) {
    double a = 0.0;
    double x = input[i];
    for (int iter = 0; iter < iters; ++iter) {
      #ifdef USE_FMA
      #pragma clang fp contract(fast)
      #endif
      a = (a * a) + x;
    }
    output[i] = a;
  }
}

int main (int argc, char** argv) {
  int iters = 2000000;
  int size = 32;
  
  if (argc != 1) {
    if (argc != 3) {
      puts("Usage: ./mandelbrot <iters> <size>");
      exit(-1);
    }
    int i = atoi(argv[1]);
    int s = atoi(argv[2]);
    if (i < 0) {
      puts("Using default iteration count");
    } else {
      iters = i;
    }
    if (s < 0 || s > SIZE) {
      puts("Using default size");
    } else {
      size = s;
    }
  }
  
  mandelbrot (iters, size);
  
  printf ("Outputs:\n");
  for (int i = 0; i < size; i += 4) {
    printf ("%d: %f %f %f %f\n", i, output[i+0], output[i+1], output[i+2], output[i+3]);
  }
  return 0;
}