#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;

    // Find Y bounding range
    double min_y = pts[1];
    double max_y = pts[1];
    for (int i = 1; i < M; ++i) {
        double y = pts[i * stride + 1];
        if (y < min_y) min_y = y;
        if (y > max_y) max_y = y;
    }
    double y_range = max_y - min_y;
    double corner_bonus = (y_range > 0.0) ? (0.1 * y_range < 2.0 ? 0.1 * y_range : 2.0) : 2.0;

    int best_idx = 0;
    double best_score = -1e30;
    double best_x = 1e30;

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

        double xi = pts[i * stride + 0];
        double yi = pts[i * stride + 1];

        double vin_x = xi - pts[prev * stride + 0];
        double vin_y = yi - pts[prev * stride + 1];
        double vout_x = pts[next * stride + 0] - xi;
        double vout_y = pts[next * stride + 1] - yi;

        double l_in = sqrt(vin_x * vin_x + vin_y * vin_y);
        double l_out = sqrt(vout_x * vout_x + vout_y * vout_y);

        double sharpness = 0.0;
        if (l_in > 1e-7 && l_out > 1e-7) {
            double dot = vin_x * vout_x + vin_y * vout_y;
            double cos_ang = dot / (l_in * l_out);
            if (cos_ang > 1.0) cos_ang = 1.0;
            if (cos_ang < -1.0) cos_ang = -1.0;
            sharpness = 1.0 - cos_ang; // 0 for straight line, 1 for 90 deg corner, 2 for 180 deg turnaround
        }

        double score = yi + sharpness * corner_bonus;

        // Highest score wins (back of the object). Tie-break: smallest X for deterministic consistency across layers
        if (score > best_score + 1e-6) {
            best_score = score;
            best_x = xi;
            best_idx = i;
        } else if (fabs(score - best_score) <= 1e-6) {
            if (xi < best_x) {
                best_x = xi;
                best_idx = i;
            }
        }
    }

    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;
}