#include <stdint.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#ifdef _WIN32
#define EXPORT __declspec(dllexport)
#else
#define EXPORT __attribute__((visibility("default")))
#endif

EXPORT int find_sharpest_seam_c(const double* pts, int M, int stride, double scale) {
    if (M < 3) return 0;

    int64_t x_buf[1024];
    int64_t y_buf[1024];
    int64_t* x = (M <= 1024) ? x_buf : (int64_t*)malloc(M * sizeof(int64_t));
    int64_t* y = (M <= 1024) ? y_buf : (int64_t*)malloc(M * sizeof(int64_t));

    for (int i = 0; i < M; ++i) {
        x[i] = (int64_t)llround(pts[i * stride + 0] * scale);
        y[i] = (int64_t)llround(pts[i * stride + 1] * scale);
    }

    int64_t signed_area = (x[M - 1] * y[0] - x[0] * y[M - 1]);
    for (int i = 0; i < M - 1; ++i) {
        signed_area += (x[i] * y[i + 1] - x[i + 1] * y[i]);
    }

    if (signed_area == 0) {
        if (M > 1024) { free(x); free(y); }
        return 0;
    }

    int best_idx = 0;
    double min_cos_concave = 2.0;
    int best_concave_idx = -1;
    double min_cos_overall = 2.0;

    for (int i = 0; i < M; ++i) {
        int prev = (i == 0) ? M - 1 : i - 1;
        int next = (i == M - 1) ? 0 : i + 1;

        int64_t vin_x = x[i] - x[prev];
        int64_t vin_y = y[i] - y[prev];
        int64_t vout_x = x[next] - x[i];
        int64_t vout_y = y[next] - y[i];

        int64_t l2_in = vin_x * vin_x + vin_y * vin_y;
        int64_t l2_out = vout_x * vout_x + vout_y * vout_y;

        if (l2_in == 0 || l2_out == 0) continue;

        int64_t cross = vin_x * vout_y - vin_y * vout_x;
        int is_concave = (signed_area > 0) ? (cross < 0) : (cross > 0);

        double dot = (double)(vin_x * vout_x + vin_y * vout_y);
        double denom = sqrt((double)l2_in * (double)l2_out);
        double cos_ang = dot / denom;

        if (is_concave) {
            if (cos_ang < min_cos_concave) {
                min_cos_concave = cos_ang;
                best_concave_idx = i;
            }
        }
        if (cos_ang < min_cos_overall) {
            min_cos_overall = cos_ang;
            best_idx = i;
        }
    }

    if (M > 1024) { free(x); free(y); }

    if (best_concave_idx >= 0) return best_concave_idx;
    return best_idx;
}

EXPORT void align_path_seam_sharpest_c(const double* in_pts, int M, int stride, double* out_pts, double scale) {
    int best_idx = find_sharpest_seam_c(in_pts, M, stride, scale);
    if (best_idx == 0) {
        memcpy(out_pts, in_pts, (M + 1) * stride * sizeof(double));
        return;
    }
    int split = M - best_idx;
    memcpy(out_pts, in_pts + (best_idx * stride), split * stride * sizeof(double));
    memcpy(out_pts + (split * stride), in_pts, (best_idx + 1) * stride * sizeof(double));
}

// Fast 64-bit hash for 3 float32 values (12 bytes)
static inline uint64_t hash_v3(const float* v) {
    const uint32_t* u = (const uint32_t*)v;
    uint64_t h = (uint64_t)u[0] * 0x9e3779b97f4a7c15ULL;
    h ^= (uint64_t)u[1] * 0xbf58476d1ce4e5b9ULL;
    h = (h ^ (h >> 30)) * 0x94d049bb133111ebULL;
    h ^= (uint64_t)u[2] * 0x4bfb2f35a09a447fULL;
    h = (h ^ (h >> 27)) * 0x94d049bb133111ebULL;
    return h ^ (h >> 31);
}

static inline int v3_eq(const float* a, const float* b) {
    const uint32_t* ua = (const uint32_t*)a;
    const uint32_t* ub = (const uint32_t*)b;
    return (ua[0] == ub[0]) && (ua[1] == ub[1]) && (ua[2] == ub[2]);
}

EXPORT int weld_vertices_c(
    const float* in_verts,
    int num_in_verts,
    float* out_unique_verts,
    uint32_t* out_inv_indices
) {
    if (num_in_verts <= 0) return 0;

    uint32_t table_size = 1024;
    while (table_size < (uint32_t)num_in_verts * 2) {
        table_size <<= 1;
    }
    uint32_t mask = table_size - 1;

    uint32_t* table = (uint32_t*)calloc(table_size, sizeof(uint32_t));
    if (!table) return -1;

    int unique_count = 0;

    for (int i = 0; i < num_in_verts; ++i) {
        const float* v = in_verts + (i * 3);
        uint64_t h = hash_v3(v);
        uint32_t slot = (uint32_t)(h & mask);

        while (1) {
            uint32_t entry = table[slot];
            if (entry == 0) {
                int uid = unique_count++;
                out_unique_verts[uid * 3 + 0] = v[0];
                out_unique_verts[uid * 3 + 1] = v[1];
                out_unique_verts[uid * 3 + 2] = v[2];
                table[slot] = (uint32_t)(uid + 1);
                out_inv_indices[i] = (uint32_t)uid;
                break;
            }
            uint32_t uid = entry - 1;
            const float* existing_v = out_unique_verts + (uid * 3);
            if (v3_eq(v, existing_v)) {
                out_inv_indices[i] = (uint32_t)uid;
                break;
            }
            slot = (slot + 1) & mask;
        }
    }

    free(table);
    return unique_count;
}

EXPORT double polygon_signed_area_2d_c(const double* pts, int N, int stride) {
    if (N < 3) return 0.0;
    double area = pts[(N - 1) * stride + 0] * pts[1] - pts[0] * pts[(N - 1) * stride + 1];
    for (int i = 0; i < N - 1; ++i) {
        area += pts[i * stride + 0] * pts[(i + 1) * stride + 1] - pts[(i + 1) * stride + 0] * pts[i * stride + 1];
    }
    return area * 0.5;
}