struct point { float32 x; float32 y; }

float32 distance(const point &p1, const point &p2) {
    float32 x1 = p1.x; float32 x2 = p2.x;
    float32 y1 = p1.y; float32 y2 = p2.y;
    float32 dx = (x2 - x1); float32 dy = (y2 - y1);
    return sqrt(dx * dx + dy * dy);
}

point add(const point &p1, const point &p2) {
    return new point(p1.x + p2.x, p1.y + p2.y);
}

point divide(const point &p, float32 div) {
    if(div == 0.0f) return new point(0.0f, 0.0f); // edge case...
    return new point(p.x / div, p.y / div);
}

type temp_center = pair<point, float32>;
pair<uint32, float32> mk_acc(uint32 idx, float32 dist) { return new pair<uint32, float32>(idx, dist); }
temp_center mk_temp_center(point p, float32 counter) { return new temp_center(p, counter); }

@Export
void kmeans(const vec<point> &points, vec<point> &centers) {
    vec<temp_center> temp{centers.size()};
    project p in points into temp {
        pair<uint32, float32> acc;
        @NoVectorize
        accumulate (c, idx) in centers with mk_acc(0u, INF) into acc {
            float32 dist = distance(p, c);
            if(dist < acc.second) continue with mk_acc(idx, dist);
            else continue;
        }

        uint32 closest_idx = acc.first;
        continue[closest_idx] with (curr) -> mk_temp_center(add(p, curr.first), curr.second + 1.0f);
    }

    project c in temp into centers {
        continue with divide(c.first, c.second);
    }
}