// simd module  — TezzNative 2026 high-performance vector ops
//
// Width tiers:
//   v4f_*  — 4-wide f32 (SSE baseline, backward-compatible)
//   v8f_*  — 8-wide f32 (AVX2 tier, 2x throughput)
//   v8i_*  — 8-wide i32
//
// AI kernels: relu, relu6, leaky, gelu_approx, sigmoid_approx, tanh_approx, fmadd
// Bulk helpers: map_relu, map_relu6, map_sigmoid, map_gelu, map_tanh
// Reductions: v8f_reduce_sum, v8f_reduce_max, v8f_reduce_min
// BLAS-style: dot_product, axpby, scale_inplace

typedef [float;4] Vec4f
typedef [float;8] Vec8f
typedef [int;4]   Vec4i
typedef [int;8]   Vec8i

// ─────────────────────────────────────────────
// 4-wide ops (backward compatible API)
// ─────────────────────────────────────────────

fn v4f_add(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = a[0] + b[0]
    out[1] = a[1] + b[1]
    out[2] = a[2] + b[2]
    out[3] = a[3] + b[3]
  ret 0

fn v4f_sub(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = a[0] - b[0]
    out[1] = a[1] - b[1]
    out[2] = a[2] - b[2]
    out[3] = a[3] - b[3]
  ret 0

fn v4f_mul(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = a[0] * b[0]
    out[1] = a[1] * b[1]
    out[2] = a[2] * b[2]
    out[3] = a[3] * b[3]
  ret 0

fn v4f_div(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = a[0] / b[0]
    out[1] = a[1] / b[1]
    out[2] = a[2] / b[2]
    out[3] = a[3] / b[3]
  ret 0

fn v4f_dot(a:*float, b:*float) -> float:
  unsafe:
    ret a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + a[3]*b[3]

fn v4i_add(out:*int, a:*int, b:*int) -> int:
  unsafe:
    out[0] = a[0] + b[0]
    out[1] = a[1] + b[1]
    out[2] = a[2] + b[2]
    out[3] = a[3] + b[3]
  ret 0

fn v4i_mul(out:*int, a:*int, b:*int) -> int:
  unsafe:
    out[0] = a[0] * b[0]
    out[1] = a[1] * b[1]
    out[2] = a[2] * b[2]
    out[3] = a[3] * b[3]
  ret 0

fn v4f_add_strided(out:*float, a:*float, b:*float, stride:int) -> int:
  unsafe:
    out[0]        = a[0]        + b[0]
    out[stride]   = a[stride]   + b[stride]
    out[stride*2] = a[stride*2] + b[stride*2]
    out[stride*3] = a[stride*3] + b[stride*3]
  ret 0

fn v4f_sub_strided(out:*float, a:*float, b:*float, stride:int) -> int:
  unsafe:
    out[0]        = a[0]        - b[0]
    out[stride]   = a[stride]   - b[stride]
    out[stride*2] = a[stride*2] - b[stride*2]
    out[stride*3] = a[stride*3] - b[stride*3]
  ret 0

fn v4f_mul_strided(out:*float, a:*float, b:*float, stride:int) -> int:
  unsafe:
    out[0]        = a[0]        * b[0]
    out[stride]   = a[stride]   * b[stride]
    out[stride*2] = a[stride*2] * b[stride*2]
    out[stride*3] = a[stride*3] * b[stride*3]
  ret 0

fn v4f_div_strided(out:*float, a:*float, b:*float, stride:int) -> int:
  unsafe:
    out[0]        = a[0]        / b[0]
    out[stride]   = a[stride]   / b[stride]
    out[stride*2] = a[stride*2] / b[stride*2]
    out[stride*3] = a[stride*3] / b[stride*3]
  ret 0

fn v4i_add_strided(out:*int, a:*int, b:*int, stride:int) -> int:
  unsafe:
    out[0]        = a[0]        + b[0]
    out[stride]   = a[stride]   + b[stride]
    out[stride*2] = a[stride*2] + b[stride*2]
    out[stride*3] = a[stride*3] + b[stride*3]
  ret 0

fn v4i_mul_strided(out:*int, a:*int, b:*int, stride:int) -> int:
  unsafe:
    out[0]        = a[0]        * b[0]
    out[stride]   = a[stride]   * b[stride]
    out[stride*2] = a[stride*2] * b[stride*2]
    out[stride*3] = a[stride*3] * b[stride*3]
  ret 0

fn v4f_load(out:*float, p:*float) -> int:
  unsafe:
    out[0] = p[0]
    out[1] = p[1]
    out[2] = p[2]
    out[3] = p[3]
  ret 0

fn v4f_store(p:*float, v:*float) -> int:
  unsafe:
    p[0] = v[0]
    p[1] = v[1]
    p[2] = v[2]
    p[3] = v[3]
  ret 0

fn v4i_load(out:*int, p:*int) -> int:
  unsafe:
    out[0] = p[0]
    out[1] = p[1]
    out[2] = p[2]
    out[3] = p[3]
  ret 0

fn v4i_store(p:*int, v:*int) -> int:
  unsafe:
    p[0] = v[0]
    p[1] = v[1]
    p[2] = v[2]
    p[3] = v[3]
  ret 0

fn v4f_hsum(p:*float) -> float:
  unsafe:
    ret p[0] + p[1] + p[2] + p[3]

// ─────────────────────────────────────────────
// 8-wide f32 ops (AVX2 tier — 2x throughput)
// ─────────────────────────────────────────────

fn v8f_add(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = a[0] + b[0]
    out[1] = a[1] + b[1]
    out[2] = a[2] + b[2]
    out[3] = a[3] + b[3]
    out[4] = a[4] + b[4]
    out[5] = a[5] + b[5]
    out[6] = a[6] + b[6]
    out[7] = a[7] + b[7]
  ret 0

fn v8f_sub(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = a[0] - b[0]
    out[1] = a[1] - b[1]
    out[2] = a[2] - b[2]
    out[3] = a[3] - b[3]
    out[4] = a[4] - b[4]
    out[5] = a[5] - b[5]
    out[6] = a[6] - b[6]
    out[7] = a[7] - b[7]
  ret 0

fn v8f_mul(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = a[0] * b[0]
    out[1] = a[1] * b[1]
    out[2] = a[2] * b[2]
    out[3] = a[3] * b[3]
    out[4] = a[4] * b[4]
    out[5] = a[5] * b[5]
    out[6] = a[6] * b[6]
    out[7] = a[7] * b[7]
  ret 0

fn v8f_div(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = a[0] / b[0]
    out[1] = a[1] / b[1]
    out[2] = a[2] / b[2]
    out[3] = a[3] / b[3]
    out[4] = a[4] / b[4]
    out[5] = a[5] / b[5]
    out[6] = a[6] / b[6]
    out[7] = a[7] / b[7]
  ret 0

// Fused multiply-add: out[i] += a[i] * b[i]
fn v8f_fmadd(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = out[0] + a[0] * b[0]
    out[1] = out[1] + a[1] * b[1]
    out[2] = out[2] + a[2] * b[2]
    out[3] = out[3] + a[3] * b[3]
    out[4] = out[4] + a[4] * b[4]
    out[5] = out[5] + a[5] * b[5]
    out[6] = out[6] + a[6] * b[6]
    out[7] = out[7] + a[7] * b[7]
  ret 0

// Fused multiply-sub: out[i] -= a[i] * b[i]
fn v8f_fmsub(out:*float, a:*float, b:*float) -> int:
  unsafe:
    out[0] = out[0] - a[0] * b[0]
    out[1] = out[1] - a[1] * b[1]
    out[2] = out[2] - a[2] * b[2]
    out[3] = out[3] - a[3] * b[3]
    out[4] = out[4] - a[4] * b[4]
    out[5] = out[5] - a[5] * b[5]
    out[6] = out[6] - a[6] * b[6]
    out[7] = out[7] - a[7] * b[7]
  ret 0

fn v8f_load(out:*float, p:*float) -> int:
  unsafe:
    out[0] = p[0]
    out[1] = p[1]
    out[2] = p[2]
    out[3] = p[3]
    out[4] = p[4]
    out[5] = p[5]
    out[6] = p[6]
    out[7] = p[7]
  ret 0

fn v8f_store(p:*float, v:*float) -> int:
  unsafe:
    p[0] = v[0]
    p[1] = v[1]
    p[2] = v[2]
    p[3] = v[3]
    p[4] = v[4]
    p[5] = v[5]
    p[6] = v[6]
    p[7] = v[7]
  ret 0

fn v8f_zero(p:*float) -> int:
  unsafe:
    p[0] = 0.0
    p[1] = 0.0
    p[2] = 0.0
    p[3] = 0.0
    p[4] = 0.0
    p[5] = 0.0
    p[6] = 0.0
    p[7] = 0.0
  ret 0

fn v8f_fill(p:*float, v:float) -> int:
  unsafe:
    p[0] = v
    p[1] = v
    p[2] = v
    p[3] = v
    p[4] = v
    p[5] = v
    p[6] = v
    p[7] = v
  ret 0

fn v8f_scale(out:*float, a:*float, s:float) -> int:
  unsafe:
    out[0] = a[0] * s
    out[1] = a[1] * s
    out[2] = a[2] * s
    out[3] = a[3] * s
    out[4] = a[4] * s
    out[5] = a[5] * s
    out[6] = a[6] * s
    out[7] = a[7] * s
  ret 0

fn v8f_hsum(p:*float) -> float:
  unsafe:
    ret p[0] + p[1] + p[2] + p[3] + p[4] + p[5] + p[6] + p[7]

// ─────────────────────────────────────────────
// 8-wide int ops
// ─────────────────────────────────────────────

fn v8i_add(out:*int, a:*int, b:*int) -> int:
  unsafe:
    out[0] = a[0] + b[0]
    out[1] = a[1] + b[1]
    out[2] = a[2] + b[2]
    out[3] = a[3] + b[3]
    out[4] = a[4] + b[4]
    out[5] = a[5] + b[5]
    out[6] = a[6] + b[6]
    out[7] = a[7] + b[7]
  ret 0

fn v8i_mul(out:*int, a:*int, b:*int) -> int:
  unsafe:
    out[0] = a[0] * b[0]
    out[1] = a[1] * b[1]
    out[2] = a[2] * b[2]
    out[3] = a[3] * b[3]
    out[4] = a[4] * b[4]
    out[5] = a[5] * b[5]
    out[6] = a[6] * b[6]
    out[7] = a[7] * b[7]
  ret 0

fn v8i_zero(p:*int) -> int:
  unsafe:
    p[0] = 0
    p[1] = 0
    p[2] = 0
    p[3] = 0
    p[4] = 0
    p[5] = 0
    p[6] = 0
    p[7] = 0
  ret 0

fn v8i_fill(p:*int, v:int) -> int:
  unsafe:
    p[0] = v
    p[1] = v
    p[2] = v
    p[3] = v
    p[4] = v
    p[5] = v
    p[6] = v
    p[7] = v
  ret 0

// ─────────────────────────────────────────────
// High-throughput kernels for numeric compute paths
// ─────────────────────────────────────────────

// ReLU: out[i] = max(0, a[i])  unrolled x8
fn v8f_relu(out:*float, a:*float) -> int:
  unsafe:
    v0:float = a[0]
    if v0 < 0.0:
      v0 = 0.0
    out[0] = v0
    v1:float = a[1]
    if v1 < 0.0:
      v1 = 0.0
    out[1] = v1
    v2:float = a[2]
    if v2 < 0.0:
      v2 = 0.0
    out[2] = v2
    v3:float = a[3]
    if v3 < 0.0:
      v3 = 0.0
    out[3] = v3
    v4:float = a[4]
    if v4 < 0.0:
      v4 = 0.0
    out[4] = v4
    v5:float = a[5]
    if v5 < 0.0:
      v5 = 0.0
    out[5] = v5
    v6:float = a[6]
    if v6 < 0.0:
      v6 = 0.0
    out[6] = v6
    v7:float = a[7]
    if v7 < 0.0:
      v7 = 0.0
    out[7] = v7
  ret 0

// ReLU6: clamp(a[i], 0, 6)
fn v8f_relu6(out:*float, a:*float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      v:float = a[i]
      if v < 0.0:
        v = 0.0
      if v > 6.0:
        v = 6.0
      out[i] = v
    i = i + 1
  ret 0

// Leaky ReLU: x >= 0 ? x : alpha*x
fn v8f_leaky(out:*float, a:*float, alpha:float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      v:float = a[i]
      if v < 0.0:
        out[i] = alpha * v
      else:
        out[i] = v
    i = i + 1
  ret 0

// Elementwise max / min
fn v8f_max(out:*float, a:*float, b:*float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      v:float = a[i]
      w:float = b[i]
      if v > w:
        out[i] = v
      else:
        out[i] = w
    i = i + 1
  ret 0

fn v8f_min(out:*float, a:*float, b:*float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      v:float = a[i]
      w:float = b[i]
      if v < w:
        out[i] = v
      else:
        out[i] = w
    i = i + 1
  ret 0

// Clamp all 8 to [lo, hi]
fn v8f_clamp(out:*float, a:*float, lo:float, hi:float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      v:float = a[i]
      if v < lo:
        v = lo
      if v > hi:
        v = hi
      out[i] = v
    i = i + 1
  ret 0

// AXPBY: y[i] = alpha * x[i] + beta * y[i]  (BLAS-style, 8 elements)
fn v8f_axpby(y:*float, x:*float, alpha:float, beta:float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      y[i] = alpha * x[i] + beta * y[i]
    i = i + 1
  ret 0

// Fast sigmoid: 0.5 * x / (1 + |x|) + 0.5  (Pade approx, no exp needed)
// Max error vs true sigmoid: ~1.5% — good enough for training
fn v8f_sigmoid_approx(out:*float, a:*float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      x:float = a[i]
      ax:float = x
      if ax < 0.0:
        ax = 0.0 - ax
      out[i] = 0.5 * (x / (1.0 + ax)) + 0.5
    i = i + 1
  ret 0

// Fast GELU: x * sigmoid(1.702 * x), <0.5% error
fn v8f_gelu_approx(out:*float, a:*float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      x:float = a[i]
      s:float = 1.702 * x
      as2:float = s
      if as2 < 0.0:
        as2 = 0.0 - as2
      sig:float = 0.5 * (s / (1.0 + as2)) + 0.5
      out[i] = x * sig
    i = i + 1
  ret 0

// Fast tanh: 2*sigmoid(2x) - 1
fn v8f_tanh_approx(out:*float, a:*float) -> int:
  i:int = 0
  while i < 8:
    unsafe:
      x:float = a[i]
      s:float = 2.0 * x
      as3:float = s
      if as3 < 0.0:
        as3 = 0.0 - as3
      sig:float = 0.5 * (s / (1.0 + as3)) + 0.5
      out[i] = 2.0 * sig - 1.0
    i = i + 1
  ret 0

// ─────────────────────────────────────────────
// Bulk reduction ops (n elements, 8-wide unrolled)
// ─────────────────────────────────────────────

// Sum n floats with 8 parallel accumulators + scalar tail
fn v8f_reduce_sum(p:*float, n:int) -> float:
  acc0:float = 0.0
  acc1:float = 0.0
  acc2:float = 0.0
  acc3:float = 0.0
  acc4:float = 0.0
  acc5:float = 0.0
  acc6:float = 0.0
  acc7:float = 0.0
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      acc0 = acc0 + p[i + 0]
      acc1 = acc1 + p[i + 1]
      acc2 = acc2 + p[i + 2]
      acc3 = acc3 + p[i + 3]
      acc4 = acc4 + p[i + 4]
      acc5 = acc5 + p[i + 5]
      acc6 = acc6 + p[i + 6]
      acc7 = acc7 + p[i + 7]
    i = i + 8
  while i < n:
    unsafe:
      acc0 = acc0 + p[i]
    i = i + 1
  ret acc0 + acc1 + acc2 + acc3 + acc4 + acc5 + acc6 + acc7

// Max of n floats with 8 parallel trackers
fn v8f_reduce_max(p:*float, n:int) -> float:
  if n <= 0:
    ret 0.0
  m0:float = 0.0
  unsafe:
    m0 = p[0]
  m1:float = m0
  m2:float = m0
  m3:float = m0
  m4:float = m0
  m5:float = m0
  m6:float = m0
  m7:float = m0
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      v0:float = p[i + 0]
      if v0 > m0:
        m0 = v0
      v1:float = p[i + 1]
      if v1 > m1:
        m1 = v1
      v2:float = p[i + 2]
      if v2 > m2:
        m2 = v2
      v3:float = p[i + 3]
      if v3 > m3:
        m3 = v3
      v4:float = p[i + 4]
      if v4 > m4:
        m4 = v4
      v5:float = p[i + 5]
      if v5 > m5:
        m5 = v5
      v6:float = p[i + 6]
      if v6 > m6:
        m6 = v6
      v7:float = p[i + 7]
      if v7 > m7:
        m7 = v7
    i = i + 8
  while i < n:
    unsafe:
      v:float = p[i]
      if v > m0:
        m0 = v
    i = i + 1
  if m1 > m0:
    m0 = m1
  if m2 > m0:
    m0 = m2
  if m3 > m0:
    m0 = m3
  if m4 > m0:
    m0 = m4
  if m5 > m0:
    m0 = m5
  if m6 > m0:
    m0 = m6
  if m7 > m0:
    m0 = m7
  ret m0

// Min of n floats
fn v8f_reduce_min(p:*float, n:int) -> float:
  if n <= 0:
    ret 0.0
  m0:float = 0.0
  unsafe:
    m0 = p[0]
  m1:float = m0
  m2:float = m0
  m3:float = m0
  m4:float = m0
  m5:float = m0
  m6:float = m0
  m7:float = m0
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      v0:float = p[i + 0]
      if v0 < m0:
        m0 = v0
      v1:float = p[i + 1]
      if v1 < m1:
        m1 = v1
      v2:float = p[i + 2]
      if v2 < m2:
        m2 = v2
      v3:float = p[i + 3]
      if v3 < m3:
        m3 = v3
      v4:float = p[i + 4]
      if v4 < m4:
        m4 = v4
      v5:float = p[i + 5]
      if v5 < m5:
        m5 = v5
      v6:float = p[i + 6]
      if v6 < m6:
        m6 = v6
      v7:float = p[i + 7]
      if v7 < m7:
        m7 = v7
    i = i + 8
  while i < n:
    unsafe:
      v:float = p[i]
      if v < m0:
        m0 = v
    i = i + 1
  if m1 < m0:
    m0 = m1
  if m2 < m0:
    m0 = m2
  if m3 < m0:
    m0 = m3
  if m4 < m0:
    m0 = m4
  if m5 < m0:
    m0 = m5
  if m6 < m0:
    m0 = m6
  if m7 < m0:
    m0 = m7
  ret m0

// ─────────────────────────────────────────────
// Bulk map ops (apply fn to n floats, 8 at a time)
// ─────────────────────────────────────────────

fn map_relu(dst:*float, src:*float, n:int) -> int:
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      v8f_relu(dst + i, src + i)
    i = i + 8
  while i < n:
    unsafe:
      v:float = src[i]
      if v < 0.0:
        v = 0.0
      dst[i] = v
    i = i + 1
  ret 0

fn map_relu6(dst:*float, src:*float, n:int) -> int:
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      v8f_relu6(dst + i, src + i)
    i = i + 8
  while i < n:
    unsafe:
      v:float = src[i]
      if v < 0.0:
        v = 0.0
      if v > 6.0:
        v = 6.0
      dst[i] = v
    i = i + 1
  ret 0

fn map_sigmoid(dst:*float, src:*float, n:int) -> int:
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      v8f_sigmoid_approx(dst + i, src + i)
    i = i + 8
  while i < n:
    unsafe:
      x:float = src[i]
      ax:float = x
      if ax < 0.0:
        ax = 0.0 - ax
      dst[i] = 0.5 * (x / (1.0 + ax)) + 0.5
    i = i + 1
  ret 0

fn map_gelu(dst:*float, src:*float, n:int) -> int:
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      v8f_gelu_approx(dst + i, src + i)
    i = i + 8
  while i < n:
    unsafe:
      x:float = src[i]
      s:float = 1.702 * x
      as4:float = s
      if as4 < 0.0:
        as4 = 0.0 - as4
      dst[i] = x * (0.5 * (s / (1.0 + as4)) + 0.5)
    i = i + 1
  ret 0

fn map_tanh(dst:*float, src:*float, n:int) -> int:
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      v8f_tanh_approx(dst + i, src + i)
    i = i + 8
  while i < n:
    unsafe:
      s:float = 2.0 * src[i]
      as5:float = s
      if as5 < 0.0:
        as5 = 0.0 - as5
      dst[i] = 2.0 * (0.5 * (s / (1.0 + as5)) + 0.5) - 1.0
    i = i + 1
  ret 0

// Scale n floats in-place: p[i] *= s
fn scale_inplace(p:*float, n:int, s:float) -> int:
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      p[i+0] = p[i+0] * s
      p[i+1] = p[i+1] * s
      p[i+2] = p[i+2] * s
      p[i+3] = p[i+3] * s
      p[i+4] = p[i+4] * s
      p[i+5] = p[i+5] * s
      p[i+6] = p[i+6] * s
      p[i+7] = p[i+7] * s
    i = i + 8
  while i < n:
    unsafe:
      p[i] = p[i] * s
    i = i + 1
  ret 0

// AXPBY over n elements: y[i] = alpha * x[i] + beta * y[i]
fn axpby(y:*float, x:*float, n:int, alpha:float, beta:float) -> int:
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      y[i+0] = alpha * x[i+0] + beta * y[i+0]
      y[i+1] = alpha * x[i+1] + beta * y[i+1]
      y[i+2] = alpha * x[i+2] + beta * y[i+2]
      y[i+3] = alpha * x[i+3] + beta * y[i+3]
      y[i+4] = alpha * x[i+4] + beta * y[i+4]
      y[i+5] = alpha * x[i+5] + beta * y[i+5]
      y[i+6] = alpha * x[i+6] + beta * y[i+6]
      y[i+7] = alpha * x[i+7] + beta * y[i+7]
    i = i + 8
  while i < n:
    unsafe:
      y[i] = alpha * x[i] + beta * y[i]
    i = i + 1
  ret 0

// Dot product of n floats (8 parallel accumulators reduces dependency chain)
fn dot_product(a:*float, b:*float, n:int) -> float:
  s0:float = 0.0
  s1:float = 0.0
  s2:float = 0.0
  s3:float = 0.0
  s4:float = 0.0
  s5:float = 0.0
  s6:float = 0.0
  s7:float = 0.0
  i:int = 0
  lim:int = n - (n % 8)
  while i < lim:
    unsafe:
      s0 = s0 + a[i+0] * b[i+0]
      s1 = s1 + a[i+1] * b[i+1]
      s2 = s2 + a[i+2] * b[i+2]
      s3 = s3 + a[i+3] * b[i+3]
      s4 = s4 + a[i+4] * b[i+4]
      s5 = s5 + a[i+5] * b[i+5]
      s6 = s6 + a[i+6] * b[i+6]
      s7 = s7 + a[i+7] * b[i+7]
    i = i + 8
  while i < n:
    unsafe:
      s0 = s0 + a[i] * b[i]
    i = i + 1
  ret s0 + s1 + s2 + s3 + s4 + s5 + s6 + s7
