// Copyright 2023 Google LLC
//
// This source code is licensed under the BSD-style license found in the
// LICENSE file in the root directory of this source tree.

$assert (P, H) == (19, 9)
$assert FMA in [0, 3]
$assert BATCH_TILE % 8 == 0
$assert BATCH_TILE >= 8
$SIMD_TILE = BATCH_TILE // 8
#include <assert.h>
#include <stddef.h>
#include <math.h>

#include <immintrin.h>

#include <immintrin.h>

#include "xnnpack/common.h"
#include "xnnpack/intrinsics-polyfill.h"
#include "xnnpack/microparams.h"
#include "xnnpack/vunary.h"

$COEFS = {
$  19: "-0x1.1D841Cp-32f",
$  17: "0x1.C4FC88p-26f",
$  15: "-0x1.332066p-20f",
$  13: "0x1.D1AEA2p-16f",
$  11: "-0x1.B2782Ep-12f",
$  9: "0x1.03CAEAp-8f",
$  7: "-0x1.967628p-6f",
$  5: "0x1.ABC35Cp-4f",
$  3: "-0x1.499D08p-2f",
$}

$POLY_SUFFIX = "p%dh%dt2" % (P, H)
$PARAMS_STRUCT = "avx_polynomial_" + POLY_SUFFIX
$ISA = "fma3" if FMA else "f16c"
void xnn_f16_vtanh_ukernel__${ISA}_polynomial_${POLY_SUFFIX}_u${BATCH_TILE}(
    size_t batch,
    const xnn_float16* input,
    xnn_float16* output,
    const struct xnn_f16_default_params params[restrict XNN_MIN_ELEMENTS(1)]) XNN_OOB_READS
{
  assert(batch != 0);
  assert(batch % sizeof(uint16_t) == 0);
  assert(input != NULL);
  assert(output != NULL);

  const __m256 vneg_sat_cutoff = _mm256_set1_ps(-0x1.1F0000p+2f);
  const __m256 vpos_sat_cutoff = _mm256_set1_ps(0x1.1F0000p+2f);
  XNN_FORCE_REALIZATION(vneg_sat_cutoff);
  XNN_FORCE_REALIZATION(vpos_sat_cutoff);
  $for i in reversed(range(3, P+1, 2)):
    const __m256 vc${i} = _mm256_set1_ps(${COEFS[i]});
    XNN_FORCE_REALIZATION(vc${i});

  const uint16_t* i = (const uint16_t*) input;
  uint16_t* o = (uint16_t*) output;
  $if BATCH_TILE > 8:
    for (; batch >= ${BATCH_TILE} * sizeof(uint16_t); batch -= ${BATCH_TILE} * sizeof(uint16_t)) {
      __m256 vx0 = _mm256_cvtph_ps(_mm_loadu_si128((const __m128i*) i));
      $for N in range(1, SIMD_TILE):
        __m256 vx${N} = _mm256_cvtph_ps(_mm_loadu_si128((const __m128i*) (i + ${N * 8})));
      i += ${BATCH_TILE};

      $for N in range(SIMD_TILE):
        vx${N} = _mm256_max_ps(vneg_sat_cutoff, vx${N});
      $for N in range(SIMD_TILE):
        vx${N} = _mm256_min_ps(vpos_sat_cutoff, vx${N});

      $for N in range(SIMD_TILE):
        const __m256 vt${N} = _mm256_mul_ps(vx${N}, vx${N});

      $if FMA == 3:
        $for N in range(SIMD_TILE):
          __m256 vp${N} = vc${P};
        $for i in reversed(range(3, P, 2)):
          $for N in range(SIMD_TILE):
            vp${N} = _mm256_fmadd_ps(vp${N}, vt${N}, vc${i});
      $else:
        $for N in range(SIMD_TILE):
          __m256 vp${N} = _mm256_add_ps(_mm256_mul_ps(vc${P}, vt${N}), vc${P-2});
        $for i in reversed(range(3, P-2, 2)):
          $for N in range(SIMD_TILE):
            vp${N} = _mm256_add_ps(_mm256_mul_ps(vp${N}, vt${N}), vc${i});

      $for N in range(SIMD_TILE):
        const __m256 vxt${N} = _mm256_mul_ps(vx${N}, vt${N});
      $for N in range(SIMD_TILE):
        $if FMA == 3:
          const __m256 vy${N} = _mm256_fmadd_ps(vp${N}, vxt${N}, vx${N});
        $else:
          const __m256 vy${N} = _mm256_add_ps(_mm256_mul_ps(vp${N}, vxt${N}), vx${N});

      _mm_storeu_si128((__m128i*) o, _mm256_cvtps_ph(vy0, _MM_FROUND_TO_NEAREST_INT));
      $for N in range(1, SIMD_TILE):
        _mm_storeu_si128((__m128i*) (o + ${N * 8}), _mm256_cvtps_ph(vy${N}, _MM_FROUND_TO_NEAREST_INT));
      o += ${BATCH_TILE};
    }
  for (; batch >= 8 * sizeof(uint16_t); batch -= 8 * sizeof(uint16_t)) {
    __m256 vx = _mm256_cvtph_ps(_mm_loadu_si128((const __m128i*) i));
    i += 8;

    vx = _mm256_max_ps(vneg_sat_cutoff, vx);
    vx = _mm256_min_ps(vpos_sat_cutoff, vx);

    const __m256 vt = _mm256_mul_ps(vx, vx);

    $if FMA == 3:
      __m256 vp = vc${P};
      $for i in reversed(range(3, P, 2)):
        vp = _mm256_fmadd_ps(vp, vt, vc${i});
    $else:
      __m256 vp = _mm256_add_ps(_mm256_mul_ps(vc${P}, vt), vc${P-2});
      $for i in reversed(range(3, P-2, 2)):
        vp = _mm256_add_ps(_mm256_mul_ps(vp, vt), vc${i});

    const __m256 vxt = _mm256_mul_ps(vx, vt);
    $if FMA == 3:
      const __m256 vy = _mm256_fmadd_ps(vp, vxt, vx);
    $else:
      const __m256 vy = _mm256_add_ps(_mm256_mul_ps(vp, vxt), vx);

    _mm_storeu_si128((__m128i*) o, _mm256_cvtps_ph(vy, _MM_FROUND_TO_NEAREST_INT));
    o += 8;
  }
  if (batch != 0) {
    __m256 vx = _mm256_cvtph_ps(_mm_loadu_si128((const __m128i*) i));

    vx = _mm256_max_ps(vneg_sat_cutoff, vx);
    vx = _mm256_min_ps(vpos_sat_cutoff, vx);

    const __m256 vt = _mm256_mul_ps(vx, vx);

    $if FMA == 3:
      __m256 vp = vc${P};
      $for i in reversed(range(3, P, 2)):
        vp = _mm256_fmadd_ps(vp, vt, vc${i});
    $else:
      __m256 vp = _mm256_add_ps(_mm256_mul_ps(vc${P}, vt), vc${P-2});
      $for i in reversed(range(3, P-2, 2)):
        vp = _mm256_add_ps(_mm256_mul_ps(vp, vt), vc${i});

    const __m256 vxt = _mm256_mul_ps(vx, vt);
    $if FMA == 3:
      const __m256 vy = _mm256_fmadd_ps(vp, vxt, vx);
    $else:
      const __m256 vy = _mm256_add_ps(_mm256_mul_ps(vp, vxt), vx);

    __m128i vh = _mm256_cvtps_ph(vy, _MM_FROUND_TO_NEAREST_INT);

    if (batch & (4 * sizeof(uint16_t))) {
      _mm_storel_epi64((__m128i*) o, vh);
      vh = _mm_unpackhi_epi64(vh, vh);
      o += 4;
    }
    if (batch & (2 * sizeof(uint16_t))) {
      _mm_storeu_si32(o, vh);
      vh = _mm_srli_epi64(vh, 32);
      o += 2;
    }
    if (batch & (1 * sizeof(uint16_t))) {
      *o = (uint16_t) _mm_extract_epi16(vh, 0);
    }
  }
}