static inline int mandel_scalar1(double cr, double ci, int count) {
  double zr = cr, zi = ci;
  int i, cmp;

  for (i = 0; i < count; ++i) {
    double rr = zr*zr;
    double ii = zi*zi;

    cmp = rr + ii <= 4.;
    if (!cmp) break;

    double nzr = rr - ii;
    double nzi = 2. * zr * zi;

    zr = nzr + cr;
    zi = nzi + ci;
  }
  return cmp;
}

void mandelbrot_scalar1(double x0, double y0, double x1, double y1, int width, int height, int maxIter, unsigned char* output) {
  double dx = (x1 - x0) / width;
  double dy = (y1 - y0) / width;

#pragma omp parallel for default(none) firstprivate(x0, y0, dx, dy, width, height, maxIter) shared(output) schedule(guided)
  for (int j = 0; j < height; ++j) {
    unsigned char* row = output + j*(width/8);
    for (int I = 0; I < width/8; ++I) {
      unsigned char p = 0;
      for (int i = 0; i < 8; ++i) {
        double x = x0 + (I*8+i)*dx;
        double y = y0 + j*dy;

        int b = mandel_scalar1(x, y, maxIter);
        p = (p << 1) | b;
      }
      row[I] = p;
    }
  }
}



static inline int iter_scalar2(double* pzr, double* pzi, double* pn2, const double* pcr, double ci) {
  for (int k = 0; k < 8; ++k) {
    double rr = pzr[k] * pzr[k];
    double ii = pzi[k] * pzi[k];
    double ri = pzr[k] * pzi[k];

    pn2[k] = rr + ii;

    pzr[k] = rr - ii + pcr[k];
    pzi[k] = ri + ri + ci;
  }
}

static inline unsigned char mand8_scalar2(int to_prune, const double* pcr, double ci, int maxIter) {
  double pzr[8], pzi[8], pn2[8];
  for (int k = 0; k < 8; ++k) {
    pzr[k] = pcr[k];
    pzi[k] = ci;
  }
  int i = 0;

  if (to_prune) {
    for (i = 0; i < maxIter-3; i += 4) {
      for (int j = 0; j < 4; ++j) {
        iter_scalar2(pzr, pzi, pn2, pcr, ci);
      }
      int do_break = 1;
      for (int k = 0; k < 8; ++k) {
        do_break &= !(pn2[k] < 4.);
      }
      if (do_break) return 0;
    }
    for (; i < maxIter; i++) {
      iter_scalar2(pzr, pzi, pn2, pcr, ci);
    }
  } else {
    for (i = 0; i < maxIter-5; i += 6) {
      for (int j = 0; j < 6; ++j) {
        iter_scalar2(pzr, pzi, pn2, pcr, ci);
      }
    }
    for (; i < maxIter; i++) {
      iter_scalar2(pzr, pzi, pn2, pcr, ci);
    }
  }
  unsigned char byte = 0;
  for (int k = 0; k < 8; ++k) {
    byte = (byte << 1) | (pn2[k] <= 4.);
  }
  return byte;
}


void mandelbrot_scalar2(double x0, double y0, double x1, double y1, int width, int height, int maxIter, unsigned char* output) {
  double dx = (x1 - x0) / width;
  double dy = (y1 - y0) / height;

  double X[width];
  for (int i = 0; i < width; ++i) {
    X[i] = x0 + i * dx;
  }

#pragma omp parallel for default(none) firstprivate(y0, dy, width, height, maxIter) shared(output, X) schedule(guided)
  for (int j = 0; j < height; ++j) {
    unsigned char* row = output + j*(width/8);
    double y = y0 + j*dy;
    int to_prune = 0;
    for (int I = 0; I < width/8; ++I) {
      double* px = &X[I*8];
      unsigned char byte = mand8_scalar2(to_prune, px, y, maxIter);
      to_prune = !byte;
      row[I] = byte;
    }
  }
}