#include <stddef.h>
#include <float.h>
#include <assert.h>

#define SIMD_SIZE 8

inline double sqr(double x)
{
    return x*x;
}

size_t foo(double x[], double y[], size_t n, double tourx, double toury)
{
    double closedist[SIMD_SIZE];
    size_t closeidx[SIMD_SIZE];
    double dist[SIMD_SIZE];
    size_t i, j, simd_size;
    assert(n>0);
    for (i=0; i<SIMD_SIZE; i++) {
        closedist[i] = DBL_MAX;
        closeidx[i] = i;
    }
    for (j = SIMD_SIZE; j < n-SIMD_SIZE; j += SIMD_SIZE) {
      for (i=0; i<SIMD_SIZE; i++)
        dist[i] = sqr(x[i+j]-tourx)+sqr(y[i+j]-toury);
      for (i=0; i<SIMD_SIZE; i++) {
        if (dist[i]<closedist[i]) {
          closedist[i]=dist[i];
        }
      }
      for (i=0; i<SIMD_SIZE; i++) {
        if (dist[i]<closedist[i]) {
          closeidx[i] = j+i;
        }
      }
    }
    for (i=0; i<n-j; i++)
      dist[i] = sqr(x[i+j]-tourx)+sqr(y[i+j]-toury);
    for (i=0; i<n-j; i++) {
      if (dist[i]<closedist[i]) {
        closedist[i]=dist[i];
      }
    }
    for (i=0; i<n-j; i++) {
      if (dist[i]<closedist[i]) {
        closeidx[i] = j+i;
      }
    }
    simd_size=SIMD_SIZE/2;
    for (i=0; i<simd_size; i++)
      if (closedist[i]>closedist[i+simd_size]) {
        closeidx[i]=closeidx[i+simd_size];
        closedist[i]=closedist[i+simd_size];
      }
    simd_size=SIMD_SIZE/4;
    for (i=0; i<simd_size; i++)
      if (closedist[i]>closedist[i+simd_size]) {
        closeidx[i]=closeidx[i+simd_size];
        closedist[i]=closedist[i+simd_size];
      }
    simd_size=SIMD_SIZE/8;
    double cd = closedist[0];
    size_t ci = closeidx[0];
    if (cd>closedist[1]) {
        ci=closeidx[1];
        cd=closedist[1];
      }
    return ci;
}