/* 13× faster version of hadamardsin.js (see also hadamardsin.py). */
#include <stddef.h>
#include <stdalign.h>
#include <setjmp.h>
#include <stdio.h>
#include <string.h>

#include "kmregion.h"

/* Some basic array arithmetic functions */

typedef double sc;              /* scalar */

sc *concat(kmregion *km, sc *xs, size_t nx, sc *ys, size_t ny)
{
  sc *result = km_new(km, (nx + ny) * sizeof(*result));
  memcpy(result, xs, nx * sizeof(*xs));
  memcpy(result + nx, ys, ny * sizeof(*ys));
  return result;
}

sc *add(kmregion *km, sc *xs, sc *ys, size_t n)
{
  sc *result = km_new(km, n * sizeof(*result));
  for (size_t i = 0; i < n; i++) result[i] = xs[i] + ys[i];
  return result;
}

sc *sub(kmregion *km, sc *xs, sc *ys, size_t n)
{
  sc *result = km_new(km, n * sizeof(*result));
  for (size_t i = 0; i < n; i++) result[i] = xs[i] - ys[i];
  return result;
}

sc *zeros(kmregion *km, size_t n)
{
  sc *result = km_new(km, n * sizeof(*result));
  memset(result, 0, n * sizeof(*result)); /* Assuming IEEE-754! */
  return result;
}

sc sum(sc *xs, size_t n)
{
  sc result = 0;
  for (size_t i = 0; i < n; i++) result += xs[i];
  return result;
}

/* Fast Hadamard transform. */

sc *hada(kmregion *km, sc *xs, size_t len)
{
  if (len < 2) return xs;
  size_t half = len / 2;
  sc *left = hada(km, xs, half);
  sc *right = hada(km, xs + half, half);
  return concat(km,
                add(km, left, right, half), half,
                sub(km, left, right, half), half);
}

void show(sc *xs, size_t len)
{
  printf("[");
  for (size_t i = 0; i < len; i++) printf("%s%+.1f", i ? " " : "", xs[i]);
  printf("] ");
}

int main()
{
  kmregion *km = km_start(km_libc_disc, 0);

  for (int i = 0; i < 16; i++) {
    sc *a = zeros(km, 16);
    a[i] = 1;
    show(a, 16); show(hada(km, a, 16), 16); printf("\n");
  }

  enum { n = 1024 * 256 };
  sc *b = zeros(km, n);
  b[0] = 1;
  b[7] = 3;
  printf("%f\n", (float)sum(hada(km, b, n), n));
  km_end(km);
  return 0;
}