himada-dispatch 0.1.1

Adaptive SIMD dispatch for Himada — auto-selects fastest kernel at runtime
#ifndef HIMADA_H
#define HIMADA_H

#include <stddef.h>
#include <stdint.h>

#ifdef __cplusplus
extern "C" {
#endif

// Opaque dispatch handle
typedef struct HwdnaDispatch HwdnaDispatch;

// ---------------------------------------------------------------------------
// Dispatch lifecycle
// ---------------------------------------------------------------------------

/// Create a dot-product dispatch. `name` may be NULL (default name used).
HwdnaDispatch* himada_dispatch_create(const char* name);

/// Destroy a dispatch handle. Safe to call with NULL.
void himada_dispatch_destroy(HwdnaDispatch* handle);

/// Write the selected kernel name into `buf` (null-terminated, up to `max_len` bytes).
void himada_dispatch_selected_name(const HwdnaDispatch* handle, char* buf, size_t max_len);

// ---------------------------------------------------------------------------
// Handle-based operations (reuse dispatch across calls)
// ---------------------------------------------------------------------------

/// Compute dot product: result = a · b. Both slices must have `n` elements.
double himada_dispatch_dot(HwdnaDispatch* handle, const double* a, const double* b, size_t n);

// ---------------------------------------------------------------------------
// One-shot operations  (create + destroy dispatch internally)
// ---------------------------------------------------------------------------

/// Dot product (f64): result = a · b
double himada_dot(const double* a, const double* b, size_t n);

/// Dot product (f32): result = a · b
float himada_dot_f32(const float* a, const float* b, size_t n);

/// Reduce sum (f64): result = sum(a[0..n])
double himada_reduce_sum(const double* a, size_t n);

/// Reduce sum (f32): result = sum(a[0..n])
float himada_reduce_sum_f32(const float* a, size_t n);

/// Reduce max (f64): result = max(a[0..n])
double himada_reduce_max(const double* a, size_t n);

/// Reduce max (f32): result = max(a[0..n])
float himada_reduce_max_f32(const float* a, size_t n);

/// Absolute max (f64): result = max(|a[0]|, ..., |a[n-1]|)
double himada_abs_max(const double* a, size_t n);

/// Absolute max (f32): result = max(|a[0]|, ..., |a[n-1]|)
float himada_abs_max_f32(const float* a, size_t n);

/// Argmax (f64): return index of maximum element
size_t himada_argmax(const double* a, size_t n);

/// Argmax (f32): return index of maximum element
size_t himada_argmax_f32(const float* a, size_t n);

/// Euclidean distance (f64): sqrt(sum((a[i] - b[i])^2))
double himada_euclidean_distance(const double* a, const double* b, size_t n);

/// Euclidean distance (f32): sqrt(sum((a[i] - b[i])^2))
float himada_euclidean_distance_f32(const float* a, const float* b, size_t n);

/// Cosine similarity (f64): dot(a,b) / (|a| * |b|)
double himada_cosine_similarity(const double* a, const double* b, size_t n);

/// Cosine similarity (f32): dot(a,b) / (|a| * |b|)
float himada_cosine_similarity_f32(const float* a, const float* b, size_t n);

/// Softmax (f64): output[i] = exp(input[i] - max) / sum(exp(input - max))
void himada_softmax(const double* input, double* output, size_t n);

/// Softmax (f32): output[i] = exp(input[i] - max) / sum(exp(input - max))
void himada_softmax_f32(const float* input, float* output, size_t n);

/// Hadamard product (f64): c[i] = a[i] * b[i]
void himada_hadamard_product(const double* a, const double* b, double* c, size_t n);

/// Hadamard product (f32): c[i] = a[i] * b[i]
void himada_hadamard_product_f32(const float* a, const float* b, float* c, size_t n);

/// Add (f64): c[i] = a[i] + b[i]
void himada_add(const double* a, const double* b, double* c, size_t n);

/// Add (f32): c[i] = a[i] + b[i]
void himada_add_f32(const float* a, const float* b, float* c, size_t n);

/// Subtract (f64): c[i] = a[i] - b[i]
void himada_sub(const double* a, const double* b, double* c, size_t n);

/// Subtract (f32): c[i] = a[i] - b[i]
void himada_sub_f32(const float* a, const float* b, float* c, size_t n);

/// Multiply (f64): c[i] = a[i] * b[i]
void himada_mul(const double* a, const double* b, double* c, size_t n);

/// Multiply (f32): c[i] = a[i] * b[i]
void himada_mul_f32(const float* a, const float* b, float* c, size_t n);

/// Negate (f64): c[i] = -a[i]
void himada_negate(const double* a, double* c, size_t n);

/// Negate (f32): c[i] = -a[i]
void himada_negate_f32(const float* a, float* c, size_t n);

/// Clamp (f64): c[i] = a[i].clamp(lo, hi)
void himada_clamp(const double* a, double lo, double hi, double* c, size_t n);

/// Clamp (f32): c[i] = a[i].clamp(lo, hi)
void himada_clamp_f32(const float* a, float lo, float hi, float* c, size_t n);

/// Matrix multiply (f64): C += A * B (flat N×N, row-major). C must be zeroed or pre-initialized.
void himada_matmul(const double* a, const double* b, double* c, size_t n);

/// Matrix multiply (f32): C += A * B (flat N×N, row-major). C must be zeroed or pre-initialized.
void himada_matmul_f32(const float* a, const float* b, float* c, size_t n);

/// Fused matmul + bias + ReLU (f64): C = ReLU(A * W + bias)
///   A ∈ ℝ^{M×K}, W ∈ ℝ^{K×N}, bias ∈ ℝ^{N}, C ∈ ℝ^{M×N}
void himada_matmul_bias_relu(const double* a, const double* w, const double* bias,
                             double* c, size_t m, size_t k, size_t n);

/// Memchr: find first occurrence of `byte` in `data[0..n]`.
/// Returns index or -1 if not found.
int64_t himada_memchr(uint8_t byte, const uint8_t* data, size_t n);

/// Fused dot-clamp-reduce: sum(clamp(a[i] * b[i], 0, 1))
double himada_fused_dot_clamp_reduce(const double* a, const double* b, size_t n);

#ifdef __cplusplus
}
#endif

#endif // HIMADA_H