/*=============================================================================
   Copyright (c) 2018 Joel de Guzman. All rights reserved.
   Copyright (c) 2024 AshRa al-Baschir, Shr3D Industries. All rights reserved.

   Distributed under the Boost Software License, Version 1.0. (See accompanying
   file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
=============================================================================*/

#include "amalgamated_q.h"

#include <assert.h>
#include <array>
#include <algorithm>
#include <chrono>
#include <limits.h>
#include <math.h>
#include <vector>

//#define CYCFI_Q_DECIBEL_HPP_FEBRUARY_21_2018

#define CYCFI_ASSERT(x, msg) assert(((x) && msg)); (void)(x), void(msg)

typedef double decibel;

namespace cycfi
{
  ////////////////////////////////////////////////////////////////////////////
  // Ignore utility for ignoring unused params.
  ////////////////////////////////////////////////////////////////////////////
  template <typename... T>
  inline void ignore(T&&...) {}

  ////////////////////////////////////////////////////////////////////////////
  // Constants
  ////////////////////////////////////////////////////////////////////////////
  constexpr auto pi = 3.14159265358979323846264338327f;

  ////////////////////////////////////////////////////////////////////////////
  // Non-copyable base class
  ////////////////////////////////////////////////////////////////////////////
  struct non_copyable
  {
    non_copyable() = default;
    non_copyable(non_copyable const& rhs) = delete;
    ~non_copyable() = default;
    non_copyable& operator=(non_copyable const&) = delete;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Unused (parameter and variable) utility
  ////////////////////////////////////////////////////////////////////////////
  template <typename... T>
  void unused(T&&...) {}

  ////////////////////////////////////////////////////////////////////////////
  // Time
  ////////////////////////////////////////////////////////////////////////////
  using duration = std::chrono::duration<double>;
  using microseconds = std::chrono::duration<double, std::micro>;
  using milliseconds = std::chrono::duration<double, std::milli>;
  using seconds = std::chrono::duration<double>;
  using minutes = std::chrono::duration<double, std::ratio<60>>;
  using hours = std::chrono::duration<double, std::ratio<60 * 60>>;
  using time_point = std::chrono::time_point<std::chrono::steady_clock, duration>;

  ////////////////////////////////////////////////////////////////////////////
  // constexpr utilities
  ////////////////////////////////////////////////////////////////////////////
  constexpr bool equal(char const* lhs, char const* rhs)
  {
    if (lhs == rhs)
      return true;
    if (!lhs || !rhs)
      return false;
    while (*lhs || *rhs)
      if (*lhs++ != *rhs++)
        return false;
    return true;
  }

  ////////////////////////////////////////////////////////////////////////////
  // Metaprogramming utilities
  ////////////////////////////////////////////////////////////////////////////
  template <typename T>
  struct remove_cvref
  {
    using type = typename std::remove_cv<typename std::remove_reference<T>::type>::type;
  };

  template <typename T>
  using remove_cvref_t = typename remove_cvref<T>::type;

  template <typename T, typename... Rest>
  struct is_arithmetic
  {
    static constexpr bool value
      = std::is_arithmetic<T>::value && is_arithmetic<Rest...>::value;
  };

  template <typename T>
  struct is_arithmetic<T>
  {
    static constexpr bool value = std::is_arithmetic<T>::value;
  };

  namespace detail
  {
    template <std::size_t bits>
    struct int_that_fits_impl { using type = void; };

    template <>
    struct int_that_fits_impl<8> { using type = std::int8_t; };

    template <>
    struct int_that_fits_impl<16> { using type = std::int16_t; };

    template <>
    struct int_that_fits_impl<32> { using type = std::int32_t; };

    template <>
    struct int_that_fits_impl<64> { using type = std::int64_t; };

    template <std::size_t bits>
    struct uint_that_fits_impl { using type = void; };

    template <>
    struct uint_that_fits_impl<8> { using type = std::uint8_t; };

    template <>
    struct uint_that_fits_impl<16> { using type = std::uint16_t; };

    template <>
    struct uint_that_fits_impl<32> { using type = std::uint32_t; };

    template <>
    struct uint_that_fits_impl<64> { using type = uint64_t; };

    constexpr std::size_t size_that_fits_int(std::size_t bits)
    {
      if (bits <= 8)
        return 8;
      else if (bits <= 16)
        return 16;
      else if (bits <= 32)
        return 32;
      else if (bits <= 64)
        return 64;
      return 0;
    }
  }

  template <std::size_t bits>
  struct int_that_fits
    : detail::int_that_fits_impl<detail::size_that_fits_int(bits)>
  {
    using type = typename
      detail::int_that_fits_impl<detail::size_that_fits_int(bits)>::type;
#if !defined(_MSC_VER)
    static_assert(!std::is_same<type, void>::value,
      "Error: No int type fits specified number of bits."
      );
#endif
  };

  template <std::size_t bits>
  struct uint_that_fits
    : detail::uint_that_fits_impl<detail::size_that_fits_int(bits)>
  {
    using type = typename
      detail::uint_that_fits_impl<detail::size_that_fits_int(bits)>::type;
#if !defined(_MSC_VER)
    static_assert(!std::is_same<type, void>::value,
      "Error: No int type fits specified number of bits."
      );
#endif
  };

  using natural_int = typename int_that_fits<sizeof(void*)* CHAR_BIT>::type;
  using natural_uint = typename uint_that_fits<sizeof(void*)* CHAR_BIT>::type;

  ////////////////////////////////////////////////////////////////////////////
  // Constants
 ////////////////////////////////////////////////////////////////////////////
  template <typename T>
  struct int_traits
  {
    static constexpr T max = std::numeric_limits<T>::max();
    static constexpr T min = std::numeric_limits<T>::min();
  };

  template <typename T>
  constexpr T int_max()
  {
    return int_traits<T>::max;
  }

  template <typename T>
  constexpr T int_min()
  {
    return int_traits<T>::min;
  }

  ////////////////////////////////////////////////////////////////////////////
  // integer and binary functions
  ////////////////////////////////////////////////////////////////////////////
  constexpr int16_t promote(int8_t i)
  {
    return i;
  }

  constexpr uint16_t promote(uint8_t i)
  {
    return i;
  }

  constexpr int32_t promote(int16_t i)
  {
    return i;
  }

  constexpr std::uint32_t promote(uint16_t i)
  {
    return i;
  }

  constexpr int64_t promote(int32_t i)
  {
    return i;
  }

  constexpr uint64_t promote(uint32_t i)
  {
    return i;
  }

  constexpr double promote(float i)
  {
    return i;
  }

  constexpr long double promote(double i)
  {
    return i;
  }

  template <typename T>
  constexpr T pow2(std::size_t n)
  {
    return (n == 0) ? T(1) : T(2) * pow2<T>(n - 1);
  }

  // This is needed to force compile-time evaluation
  template <typename T, size_t n>
  struct static_pow2
  {
    constexpr static T val = pow2<T>(n);
  };

  // smallest power of 2 that fits n
  template <typename T>
  constexpr T smallest_pow2(T n, T m = 1)
  {
    return (m < n) ? smallest_pow2(n, m << 1) : m;
  }

  template <typename T>
  constexpr bool is_pow2(T n)
  {
    return (n & (n - 1)) == 0;
  }

  ////////////////////////////////////////////////////////////////////////////
  // static int types
 ////////////////////////////////////////////////////////////////////////////
  template <typename T, T value_>
  struct static_int
  {
    static constexpr T value = value_;
  };

  template <int i>
  using int_ = static_int<int, i>;

  template <std::size_t i>
  using uint_ = static_int<std::size_t, i>;

  template <int8_t i>
  using int8_ = static_int<int8_t, i>;

  template <uint8_t i>
  using uint8_ = static_int<uint8_t, i>;

  template <int16_t i>
  using int16_ = static_int<int16_t, i>;

  template <uint16_t i>
  using uint16_ = static_int<uint16_t, i>;

  template <int32_t i>
  using int32_ = static_int<int32_t, i>;

  template <uint32_t i>
  using uint32_ = static_int<uint32_t, i>;

  template <int64_t i>
  using int64_ = static_int<int64_t, i>;

  template <uint64_t i>
  using uint64_ = static_int<uint64_t, i>;

  ////////////////////////////////////////////////////////////////////////////
  // Utilities
  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename U>
  inline T& clamp_max(T& val, U const& max)
  {
    if (val > max)
      val = max;
    return val;
  }

  template <typename T, typename U>
  inline T& clamp_min(T& val, U const& min)
  {
    if (val < min)
      val = min;
    return val;
  }

  template <typename T, typename U, typename V>
  inline T& clamp(T& val, U const& min, V const& max)
  {
    assert(min <= max);
    clamp_min(val, min);
    clamp_max(val, max);
    return val;
  }

  template <typename T, typename U, typename V>
  inline bool within(T const& val, U const& min, V const& max)
  {
    return (val >= min) && (val <= max);
  }

  template <typename T>
  constexpr T abs(T i)
  {
    return (i >= 0) ? i : -i;
  }

  ////////////////////////////////////////////////////////////////////////////
  // deleter: generic custom deleter for, e.g. unique_ptr.
  ////////////////////////////////////////////////////////////////////////////
  template <typename T, void(&delete_)(T*)>
  struct deleter
  {
    void operator()(T* p) { delete_(p); }
  };

  ////////////////////////////////////////////////////////////////////////////
  // Return true if little endian
  ////////////////////////////////////////////////////////////////////////////
  inline bool is_little_endian()
  {
    static_assert(sizeof(char) != sizeof(short), "Error: not usable on this machine");
    short number = 0x1;
    char* p = reinterpret_cast<char*>(&number);
    return (p[0] == 1);
  }
}


/*=====================================================================*
 *                   Copyright (C) 2012 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifdef __cplusplus
#define cast_uint32_t static_cast<uint32_t>
#else
#define cast_uint32_t (uint32_t)
#endif

 /*=====================================================================*
  *                   Copyright (C) 2011 Paul Mineiro                   *
  * All rights reserved.                                                *
  *                                                                     *
  * Redistribution and use in source and binary forms, with             *
  * or without modification, are permitted provided that the            *
  * following conditions are met:                                       *
  *                                                                     *
  *     * Redistributions of source code must retain the                *
  *     above copyright notice, this list of conditions and             *
  *     the following disclaimer.                                       *
  *                                                                     *
  *     * Redistributions in binary form must reproduce the             *
  *     above copyright notice, this list of conditions and             *
  *     the following disclaimer in the documentation and/or            *
  *     other materials provided with the distribution.                 *
  *                                                                     *
  *     * Neither the name of Paul Mineiro nor the names                *
  *     of other contributors may be used to endorse or promote         *
  *     products derived from this software without specific            *
  *     prior written permission.                                       *
  *                                                                     *
  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
  * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
  * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
  * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
  * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
  * POSSIBILITY OF SUCH DAMAGE.                                         *
  *                                                                     *
  * Contact: Paul Mineiro <paul@mineiro.com>                            *
  *=====================================================================*/

#ifndef __SSE_H_
#define __SSE_H_

#ifdef __SSE2__

#include <emmintrin.h>

#ifdef __cplusplus
namespace {
#endif // __cplusplus

  typedef __m128 v4sf;
  typedef __m128i v4si;

#define v4si_to_v4sf _mm_cvtepi32_ps
#define v4sf_to_v4si _mm_cvttps_epi32

#if _MSC_VER && !__INTEL_COMPILER
  template <class T>
  __forceinline char GetChar(T value, size_t index) { return ((char*)&value)[index]; }

#define AS_4CHARS(a) \
      GetChar(int32_t(a), 0), GetChar(int32_t(a), 1), \
      GetChar(int32_t(a), 2), GetChar(int32_t(a), 3)

#define _MM_SETR_EPI32(a0, a1, a2, a3) \
      { AS_4CHARS(a0), AS_4CHARS(a1), AS_4CHARS(a2), AS_4CHARS(a3) }

#define v4sfl(x) (const v4sf { (x), (x), (x), (x) })
#define v4sil(x) (const v4si _MM_SETR_EPI32(x, x, x, x))

  __forceinline const v4sf operator+(const v4sf& a, const v4sf& b) { return _mm_add_ps(a, b); }
  __forceinline const v4sf operator-(const v4sf& a, const v4sf& b) { return _mm_sub_ps(a, b); }
  __forceinline const v4sf operator/(const v4sf& a, const v4sf& b) { return _mm_div_ps(a, b); }
  __forceinline const v4sf operator*(const v4sf& a, const v4sf& b) { return _mm_mul_ps(a, b); }

  __forceinline const v4sf operator+(const v4sf& a) { return a; }
  __forceinline const v4sf operator-(const v4sf& a) { return _mm_xor_ps(a, _mm_castsi128_ps(_mm_set1_epi32(0x80000000))); }

  __forceinline const v4sf operator&(const v4sf& a, const v4sf& b) { return _mm_and_ps(a, b); }
  __forceinline const v4sf operator|(const v4sf& a, const v4sf& b) { return _mm_or_ps(a, b); }
  __forceinline const v4sf operator^(const v4sf& a, const v4sf& b) { return _mm_xor_ps(a, b); }

  __forceinline const v4si operator&(const v4si& a, const v4si& b) { return _mm_and_si128(a, b); }
  __forceinline const v4si operator|(const v4si& a, const v4si& b) { return _mm_or_si128(a, b); }
  __forceinline const v4si operator^(const v4si& a, const v4si& b) { return _mm_xor_si128(a, b); }

  __forceinline const v4sf operator+=(v4sf& a, const v4sf& b) { return a = a + b; }
  __forceinline const v4sf operator-=(v4sf& a, const v4sf& b) { return a = a - b; }
  __forceinline const v4sf operator*=(v4sf& a, const v4sf& b) { return a = a * b; }
  __forceinline const v4sf operator/=(v4sf& a, const v4sf& b) { return a = a / b; }

  __forceinline const v4si operator|=(v4si& a, const v4si& b) { return a = a | b; }
  __forceinline const v4si operator&=(v4si& a, const v4si& b) { return a = a & b; }
  __forceinline const v4si operator^=(v4si& a, const v4si& b) { return a = a ^ b; }
#else
#define v4sfl(x) ((const v4sf) { (x), (x), (x), (x) })
#define v2dil(x) ((const v4si) { (x), (x) })
#define v4sil(x) v2dil((((long long) (x)) << 32) | (long long) (x))
#endif

  typedef union { v4sf f; float array[4]; } v4sfindexer;
#define v4sf_index(_findx, _findi)      \
  ({                                    \
     v4sfindexer _findvx = { _findx } ; \
     _findvx.array[_findi];             \
  })
  typedef union { v4si i; int array[4]; } v4siindexer;
#define v4si_index(_iindx, _iindi)      \
  ({                                    \
     v4siindexer _iindvx = { _iindx } ; \
     _iindvx.array[_iindi];             \
  })

  typedef union { v4sf f; v4si i; } v4sfv4sipun;
#if _MSC_VER && !__INTEL_COMPILER
#define v4sf_fabs(x) _mm_and_ps(x, _mm_castsi128_ps(_mm_set1_epi32(0x7fffffff)))
#else
#define v4sf_fabs(x)                  \
  ({                                    \
     v4sfv4sipun vx;                    \
     vx.f = x;                          \
     vx.i &= v4sil (0x7FFFFFFF);        \
     vx.f;                              \
  })
#endif

#ifdef __cplusplus
} // end namespace
#endif // __cplusplus

#endif // __SSE2__

#endif // __SSE_H_
/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __FAST_EXP_H_
#define __FAST_EXP_H_

 // Underflow of exponential is common practice in numerical routines,
 // so handle it here.

static inline float
fastpow2(float p)
{
  float offset = (p < 0) ? 1.0f : 0.0f;
  float clipp = (p < -126) ? -126.0f : p;
  int w = clipp;
  float z = clipp - w + offset;
  union { uint32_t i; float f; } v = { cast_uint32_t((1 << 23) * (clipp + 121.2740575f + 27.7280233f / (4.84252568f - z) - 1.49012907f * z)) };

  return v.f;
}

static inline float
fastexp(float p)
{
  return fastpow2(1.442695040f * p);
}

static inline float
fasterpow2(float p)
{
  float clipp = (p < -126) ? -126.0f : p;
  union { uint32_t i; float f; } v = { cast_uint32_t((1 << 23) * (clipp + 126.94269504f)) };
  return v.f;
}

static inline float
fasterexp(float p)
{
  return fasterpow2(1.442695040f * p);
}

#ifdef __SSE2__

static inline v4sf
vfastpow2(const v4sf p)
{
  v4sf ltzero = _mm_cmplt_ps(p, v4sfl(0.0f));
  v4sf offset = _mm_and_ps(ltzero, v4sfl(1.0f));
  v4sf lt126 = _mm_cmplt_ps(p, v4sfl(-126.0f));
  v4sf clipp = _mm_or_ps(_mm_andnot_ps(lt126, p), _mm_and_ps(lt126, v4sfl(-126.0f)));
  v4si w = v4sf_to_v4si(clipp);
  v4sf z = clipp - v4si_to_v4sf(w) + offset;

  const v4sf c_121_2740838 = v4sfl(121.2740575f);
  const v4sf c_27_7280233 = v4sfl(27.7280233f);
  const v4sf c_4_84252568 = v4sfl(4.84252568f);
  const v4sf c_1_49012907 = v4sfl(1.49012907f);
  union { v4si i; v4sf f; } v = {
    v4sf_to_v4si(
      v4sfl(1 << 23) *
      (clipp + c_121_2740838 + c_27_7280233 / (c_4_84252568 - z) - c_1_49012907 * z)
    )
  };

  return v.f;
}

static inline v4sf
vfastexp(const v4sf p)
{
  const v4sf c_invlog_2 = v4sfl(1.442695040f);

  return vfastpow2(c_invlog_2 * p);
}

static inline v4sf
vfasterpow2(const v4sf p)
{
  const v4sf c_126_94269504 = v4sfl(126.94269504f);
  v4sf lt126 = _mm_cmplt_ps(p, v4sfl(-126.0f));
  v4sf clipp = _mm_or_ps(_mm_andnot_ps(lt126, p), _mm_and_ps(lt126, v4sfl(-126.0f)));
  union { v4si i; v4sf f; } v = { v4sf_to_v4si(v4sfl(1 << 23) * (clipp + c_126_94269504)) };
  return v.f;
}

static inline v4sf
vfasterexp(const v4sf p)
{
  const v4sf c_invlog_2 = v4sfl(1.442695040f);

  return vfasterpow2(c_invlog_2 * p);
}

#endif //__SSE2__

#endif // __FAST_EXP_H_
/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __FAST_LOG_H_
#define __FAST_LOG_H_

static inline float
fastlog2(float x)
{
  union { float f; uint32_t i; } vx = { x };
  union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | 0x3f000000 };
  float y = vx.i;
  y *= 1.1920928955078125e-7f;

  return y - 124.22551499f
    - 1.498030302f * mx.f
    - 1.72587999f / (0.3520887068f + mx.f);
}

static inline float
fastlog(float x)
{
  return 0.69314718f * fastlog2(x);
}

static inline float
fasterlog2(float x)
{
  union { float f; uint32_t i; } vx = { x };
  float y = vx.i;
  y *= 1.1920928955078125e-7f;
  return y - 126.94269504f;
}

static inline float
fasterlog(float x)
{
  //  return 0.69314718f * fasterlog2 (x);

  union { float f; uint32_t i; } vx = { x };
  float y = vx.i;
  y *= 8.2629582881927490e-8f;
  return y - 87.989971088f;
}

#ifdef __SSE2__

static inline v4sf
vfastlog2(v4sf x)
{
  union { v4sf f; v4si i; } vx = { x };
  union { v4si i; v4sf f; } mx; mx.i = (vx.i & v4sil(0x007FFFFF)) | v4sil(0x3f000000);
  v4sf y = v4si_to_v4sf(vx.i);
  y *= v4sfl(1.1920928955078125e-7f);

  const v4sf c_124_22551499 = v4sfl(124.22551499f);
  const v4sf c_1_498030302 = v4sfl(1.498030302f);
  const v4sf c_1_725877999 = v4sfl(1.72587999f);
  const v4sf c_0_3520087068 = v4sfl(0.3520887068f);

  return y - c_124_22551499
    - c_1_498030302 * mx.f
    - c_1_725877999 / (c_0_3520087068 + mx.f);
}

static inline v4sf
vfastlog(v4sf x)
{
  const v4sf c_0_69314718 = v4sfl(0.69314718f);

  return c_0_69314718 * vfastlog2(x);
}

static inline v4sf
vfasterlog2(v4sf x)
{
  union { v4sf f; v4si i; } vx = { x };
  v4sf y = v4si_to_v4sf(vx.i);
  y *= v4sfl(1.1920928955078125e-7f);

  const v4sf c_126_94269504 = v4sfl(126.94269504f);

  return y - c_126_94269504;
}

static inline v4sf
vfasterlog(v4sf x)
{
  //  const v4sf c_0_69314718 = v4sfl (0.69314718f);
  //
  //  return c_0_69314718 * vfasterlog2 (x);

  union { v4sf f; v4si i; } vx = { x };
  v4sf y = v4si_to_v4sf(vx.i);
  y *= v4sfl(8.2629582881927490e-8f);

  const v4sf c_87_989971088 = v4sfl(87.989971088f);

  return y - c_87_989971088;
}

#endif // __SSE2__

#endif // __FAST_LOG_H_
/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

 // fasterfc: not actually faster than erfcf(3) on newer machines!
 // ... although vectorized version is interesting
 //     and fastererfc is very fast

static inline float
fasterfc(float x)
{
  static const float k = 3.3509633149424609f;
  static const float a = 0.07219054755431126f;
  static const float b = 15.418191568719577f;
  static const float c = 5.609846028328545f;

  union { float f; uint32_t i; } vc = { c * x };
  float xsq = x * x;
  float xquad = xsq * xsq;

  vc.i |= 0x80000000;

  return 2.0f / (1.0f + fastpow2(k * x)) - a * x * (b * xquad - 1.0f) * fasterpow2(vc.f);
}

static inline float
fastererfc(float x)
{
  static const float k = 3.3509633149424609f;

  return 2.0f / (1.0f + fasterpow2(k * x));
}

// fasterf: not actually faster than erff(3) on newer machines!
// ... although vectorized version is interesting
//     and fastererf is very fast

static inline float
fasterf(float x)
{
  return 1.0f - fasterfc(x);
}

static inline float
fastererf(float x)
{
  return 1.0f - fastererfc(x);
}

static inline float
fastinverseerf(float x)
{
  static const float invk = 0.30004578719350504f;
  static const float a = 0.020287853348211326f;
  static const float b = 0.07236892874789555f;
  static const float c = 0.9913030456864257f;
  static const float d = 0.8059775923760193f;

  float xsq = x * x;

  return invk * fastlog2((1.0f + x) / (1.0f - x))
    + x * (a - b * xsq) / (c - d * xsq);
}

static inline float
fasterinverseerf(float x)
{
  static const float invk = 0.30004578719350504f;

  return invk * fasterlog2((1.0f + x) / (1.0f - x));
}

#ifdef __SSE2__

static inline v4sf
vfasterfc(v4sf x)
{
  const v4sf k = v4sfl(3.3509633149424609f);
  const v4sf a = v4sfl(0.07219054755431126f);
  const v4sf b = v4sfl(15.418191568719577f);
  const v4sf c = v4sfl(5.609846028328545f);

  union { v4sf f; v4si i; } vc; vc.f = c * x;
  vc.i |= v4sil(0x80000000);

  v4sf xsq = x * x;
  v4sf xquad = xsq * xsq;

  return v4sfl(2.0f) / (v4sfl(1.0f) + vfastpow2(k * x)) - a * x * (b * xquad - v4sfl(1.0f)) * vfasterpow2(vc.f);
}

static inline v4sf
vfastererfc(const v4sf x)
{
  const v4sf k = v4sfl(3.3509633149424609f);

  return v4sfl(2.0f) / (v4sfl(1.0f) + vfasterpow2(k * x));
}

static inline v4sf
vfasterf(v4sf x)
{
  return v4sfl(1.0f) - vfasterfc(x);
}

static inline v4sf
vfastererf(const v4sf x)
{
  return v4sfl(1.0f) - vfastererfc(x);
}

static inline v4sf
vfastinverseerf(v4sf x)
{
  const v4sf invk = v4sfl(0.30004578719350504f);
  const v4sf a = v4sfl(0.020287853348211326f);
  const v4sf b = v4sfl(0.07236892874789555f);
  const v4sf c = v4sfl(0.9913030456864257f);
  const v4sf d = v4sfl(0.8059775923760193f);

  v4sf xsq = x * x;

  return invk * vfastlog2((v4sfl(1.0f) + x) / (v4sfl(1.0f) - x))
    + x * (a - b * xsq) / (c - d * xsq);
}

static inline v4sf
vfasterinverseerf(v4sf x)
{
  const v4sf invk = v4sfl(0.30004578719350504f);

  return invk * vfasterlog2((v4sfl(1.0f) + x) / (v4sfl(1.0f) - x));
}

#endif //__SSE2__

/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __FAST_GAMMA_H_
#define __FAST_GAMMA_H_

 /* gamma/digamma functions only work for positive inputs */

static inline float
fastlgamma(float x)
{
  float logterm = fastlog(x * (1.0f + x) * (2.0f + x));
  float xp3 = 3.0f + x;

  return -2.081061466f
    - x
    + 0.0833333f / xp3
    - logterm
    + (2.5f + x) * fastlog(xp3);
}

static inline float
fasterlgamma(float x)
{
  return -0.0810614667f
    - x
    - fasterlog(x)
    + (0.5f + x) * fasterlog(1.0f + x);
}

static inline float
fastdigamma(float x)
{
  float twopx = 2.0f + x;
  float logterm = fastlog(twopx);

  return (-48.0f + x * (-157.0f + x * (-127.0f - 30.0f * x))) /
    (12.0f * x * (1.0f + x) * twopx * twopx)
    + logterm;
}

static inline float
fasterdigamma(float x)
{
  float onepx = 1.0f + x;

  return -1.0f / x - 1.0f / (2 * onepx) + fasterlog(onepx);
}

#ifdef __SSE2__

static inline v4sf
vfastlgamma(v4sf x)
{
  const v4sf c_1_0 = v4sfl(1.0f);
  const v4sf c_2_0 = v4sfl(2.0f);
  const v4sf c_3_0 = v4sfl(3.0f);
  const v4sf c_2_081061466 = v4sfl(2.081061466f);
  const v4sf c_0_0833333 = v4sfl(0.0833333f);
  const v4sf c_2_5 = v4sfl(2.5f);

  v4sf logterm = vfastlog(x * (c_1_0 + x) * (c_2_0 + x));
  v4sf xp3 = c_3_0 + x;

  return -c_2_081061466
    - x
    + c_0_0833333 / xp3
    - logterm
    + (c_2_5 + x) * vfastlog(xp3);
}

static inline v4sf
vfasterlgamma(v4sf x)
{
  const v4sf c_0_0810614667 = v4sfl(0.0810614667f);
  const v4sf c_0_5 = v4sfl(0.5f);
  const v4sf c_1 = v4sfl(1.0f);

  return -c_0_0810614667
    - x
    - vfasterlog(x)
    + (c_0_5 + x) * vfasterlog(c_1 + x);
}

static inline v4sf
vfastdigamma(v4sf x)
{
  v4sf twopx = v4sfl(2.0f) + x;
  v4sf logterm = vfastlog(twopx);

  return (v4sfl(-48.0f) + x * (v4sfl(-157.0f) + x * (v4sfl(-127.0f) - v4sfl(30.0f) * x))) /
    (v4sfl(12.0f) * x * (v4sfl(1.0f) + x) * twopx * twopx)
    + logterm;
}

static inline v4sf
vfasterdigamma(v4sf x)
{
  const v4sf c_1_0 = v4sfl(1.0f);
  const v4sf c_2_0 = v4sfl(2.0f);
  v4sf onepx = c_1_0 + x;

  return -c_1_0 / x - c_1_0 / (c_2_0 * onepx) + vfasterlog(onepx);
}

#endif //__SSE2__

#endif // __FAST_GAMMA_H_
/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __FAST_HYPERBOLIC_H_
#define __FAST_HYPERBOLIC_H_

static inline float
fastsinh(float p)
{
  return 0.5f * (fastexp(p) - fastexp(-p));
}

static inline float
fastersinh(float p)
{
  return 0.5f * (fasterexp(p) - fasterexp(-p));
}

static inline float
fastcosh(float p)
{
  return 0.5f * (fastexp(p) + fastexp(-p));
}

static inline float
fastercosh(float p)
{
  return 0.5f * (fasterexp(p) + fasterexp(-p));
}

static inline float
fasttanh(float p)
{
  return -1.0f + 2.0f / (1.0f + fastexp(-2.0f * p));
}

static inline float
fastertanh(float p)
{
  return -1.0f + 2.0f / (1.0f + fasterexp(-2.0f * p));
}

#ifdef __SSE2__

static inline v4sf
vfastsinh(const v4sf p)
{
  const v4sf c_0_5 = v4sfl(0.5f);

  return c_0_5 * (vfastexp(p) - vfastexp(-p));
}

static inline v4sf
vfastersinh(const v4sf p)
{
  const v4sf c_0_5 = v4sfl(0.5f);

  return c_0_5 * (vfasterexp(p) - vfasterexp(-p));
}

static inline v4sf
vfastcosh(const v4sf p)
{
  const v4sf c_0_5 = v4sfl(0.5f);

  return c_0_5 * (vfastexp(p) + vfastexp(-p));
}

static inline v4sf
vfastercosh(const v4sf p)
{
  const v4sf c_0_5 = v4sfl(0.5f);

  return c_0_5 * (vfasterexp(p) + vfasterexp(-p));
}

static inline v4sf
vfasttanh(const v4sf p)
{
  const v4sf c_1 = v4sfl(1.0f);
  const v4sf c_2 = v4sfl(2.0f);

  return -c_1 + c_2 / (c_1 + vfastexp(-c_2 * p));
}

static inline v4sf
vfastertanh(const v4sf p)
{
  const v4sf c_1 = v4sfl(1.0f);
  const v4sf c_2 = v4sfl(2.0f);

  return -c_1 + c_2 / (c_1 + vfasterexp(-c_2 * p));
}

#endif //__SSE2__

#endif // __FAST_HYPERBOLIC_H_
/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __FAST_LAMBERT_W_H_
#define __FAST_LAMBERT_W_H_

 // these functions compute the upper branch aka W_0

static inline float
fastlambertw(float x)
{
  static const float threshold = 2.26445f;

  float c = (x < threshold) ? 1.546865557f : 1.0f;
  float d = (x < threshold) ? 2.250366841f : 0.0f;
  float a = (x < threshold) ? -0.737769969f : 0.0f;

  float logterm = fastlog(c * x + d);
  float loglogterm = fastlog(logterm);

  float minusw = -a - logterm + loglogterm - loglogterm / logterm;
  float expminusw = fastexp(minusw);
  float xexpminusw = x * expminusw;
  float pexpminusw = xexpminusw - minusw;

  return (2.0f * xexpminusw - minusw * (4.0f * xexpminusw - minusw * pexpminusw)) /
    (2.0f + pexpminusw * (2.0f - minusw));
}

static inline float
fasterlambertw(float x)
{
  static const float threshold = 2.26445f;

  float c = (x < threshold) ? 1.546865557f : 1.0f;
  float d = (x < threshold) ? 2.250366841f : 0.0f;
  float a = (x < threshold) ? -0.737769969f : 0.0f;

  float logterm = fasterlog(c * x + d);
  float loglogterm = fasterlog(logterm);

  float w = a + logterm - loglogterm + loglogterm / logterm;
  float expw = fasterexp(-w);

  return (w * w + expw * x) / (1.0f + w);
}

static inline float
fastlambertwexpx(float x)
{
  static const float k = 1.1765631309f;
  static const float a = 0.94537622168f;

  float logarg = fmaxf(x, k);
  float powarg = (x < k) ? a * (x - k) : 0;

  float logterm = fastlog(logarg);
  float powterm = fasterpow2(powarg);  // don't need accuracy here

  float w = powterm * (logarg - logterm + logterm / logarg);
  float logw = fastlog(w);
  float p = x - logw;

  return w * (2.0f + p + w * (3.0f + 2.0f * p)) /
    (2.0f - p + w * (5.0f + 2.0f * w));
}

static inline float
fasterlambertwexpx(float x)
{
  static const float k = 1.1765631309f;
  static const float a = 0.94537622168f;

  float logarg = fmaxf(x, k);
  float powarg = (x < k) ? a * (x - k) : 0;

  float logterm = fasterlog(logarg);
  float powterm = fasterpow2(powarg);

  float w = powterm * (logarg - logterm + logterm / logarg);
  float logw = fasterlog(w);

  return w * (1.0f + x - logw) / (1.0f + w);
}

#ifdef __SSE2__

static inline v4sf
vfastlambertw(v4sf x)
{
  const v4sf threshold = v4sfl(2.26445f);

  v4sf under = _mm_cmplt_ps(x, threshold);
  v4sf c = _mm_or_ps(_mm_and_ps(under, v4sfl(1.546865557f)),
    _mm_andnot_ps(under, v4sfl(1.0f)));
  v4sf d = _mm_and_ps(under, v4sfl(2.250366841f));
  v4sf a = _mm_and_ps(under, v4sfl(-0.737769969f));

  v4sf logterm = vfastlog(c * x + d);
  v4sf loglogterm = vfastlog(logterm);

  v4sf minusw = -a - logterm + loglogterm - loglogterm / logterm;
  v4sf expminusw = vfastexp(minusw);
  v4sf xexpminusw = x * expminusw;
  v4sf pexpminusw = xexpminusw - minusw;

  return (v4sfl(2.0f) * xexpminusw - minusw * (v4sfl(4.0f) * xexpminusw - minusw * pexpminusw)) /
    (v4sfl(2.0f) + pexpminusw * (v4sfl(2.0f) - minusw));
}

static inline v4sf
vfasterlambertw(v4sf x)
{
  const v4sf threshold = v4sfl(2.26445f);

  v4sf under = _mm_cmplt_ps(x, threshold);
  v4sf c = _mm_or_ps(_mm_and_ps(under, v4sfl(1.546865557f)),
    _mm_andnot_ps(under, v4sfl(1.0f)));
  v4sf d = _mm_and_ps(under, v4sfl(2.250366841f));
  v4sf a = _mm_and_ps(under, v4sfl(-0.737769969f));

  v4sf logterm = vfasterlog(c * x + d);
  v4sf loglogterm = vfasterlog(logterm);

  v4sf w = a + logterm - loglogterm + loglogterm / logterm;
  v4sf expw = vfasterexp(-w);

  return (w * w + expw * x) / (v4sfl(1.0f) + w);
}

static inline v4sf
vfastlambertwexpx(v4sf x)
{
  const v4sf k = v4sfl(1.1765631309f);
  const v4sf a = v4sfl(0.94537622168f);
  const v4sf two = v4sfl(2.0f);
  const v4sf three = v4sfl(3.0f);
  const v4sf five = v4sfl(5.0f);

  v4sf logarg = _mm_max_ps(x, k);
  v4sf powarg = _mm_and_ps(_mm_cmplt_ps(x, k), a * (x - k));

  v4sf logterm = vfastlog(logarg);
  v4sf powterm = vfasterpow2(powarg);  // don't need accuracy here

  v4sf w = powterm * (logarg - logterm + logterm / logarg);
  v4sf logw = vfastlog(w);
  v4sf p = x - logw;

  return w * (two + p + w * (three + two * p)) /
    (two - p + w * (five + two * w));
}

static inline v4sf
vfasterlambertwexpx(v4sf x)
{
  const v4sf k = v4sfl(1.1765631309f);
  const v4sf a = v4sfl(0.94537622168f);

  v4sf logarg = _mm_max_ps(x, k);
  v4sf powarg = _mm_and_ps(_mm_cmplt_ps(x, k), a * (x - k));

  v4sf logterm = vfasterlog(logarg);
  v4sf powterm = vfasterpow2(powarg);

  v4sf w = powterm * (logarg - logterm + logterm / logarg);
  v4sf logw = vfasterlog(w);

  return w * (v4sfl(1.0f) + x - logw) / (v4sfl(1.0f) + w);
}

#endif // __SSE2__

#endif // __FAST_LAMBERT_W_H_

/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __FAST_POW_H_
#define __FAST_POW_H_

static inline float
fastpow(float x,
  float p)
{
  return fastpow2(p * fastlog2(x));
}

static inline float
fasterpow(float x,
  float p)
{
  return fasterpow2(p * fasterlog2(x));
}

#ifdef __SSE2__

static inline v4sf
vfastpow(const v4sf x,
  const v4sf p)
{
  return vfastpow2(p * vfastlog2(x));
}

static inline v4sf
vfasterpow(const v4sf x,
  const v4sf p)
{
  return vfasterpow2(p * vfasterlog2(x));
}

#endif //__SSE2__

#endif // __FAST_POW_H_
/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __FAST_SIGMOID_H_
#define __FAST_SIGMOID_H_

static inline float
fastsigmoid(float x)
{
  return 1.0f / (1.0f + fastexp(-x));
}

static inline float
fastersigmoid(float x)
{
  return 1.0f / (1.0f + fasterexp(-x));
}

#ifdef __SSE2__

static inline v4sf
vfastsigmoid(const v4sf x)
{
  const v4sf c_1 = v4sfl(1.0f);

  return c_1 / (c_1 + vfastexp(-x));
}

static inline v4sf
vfastersigmoid(const v4sf x)
{
  const v4sf c_1 = v4sfl(1.0f);

  return c_1 / (c_1 + vfasterexp(-x));
}

#endif //__SSE2__

#endif // __FAST_SIGMOID_H_
/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/


 // http://www.devmaster.net/forums/showthread.php?t=5784
 // fast sine variants are for x \in [ -\pi, pi ]
 // fast cosine variants are for x \in [ -\pi, pi ]
 // fast tangent variants are for x \in [ -\pi / 2, pi / 2 ]
 // "full" versions of functions handle the entire range of inputs
 // although the range reduction technique used here will be hopelessly
 // inaccurate for |x| >> 1000
 //
 // WARNING: fastsinfull, fastcosfull, and fasttanfull can be slower than
 // libc calls on older machines (!) and on newer machines are only
 // slighly faster.  however:
 //   * vectorized versions are competitive
 //   * faster full versions are competitive

static inline float
fastsin(float x)
{
  static const float fouroverpi = 1.2732395447351627f;
  static const float fouroverpisq = 0.40528473456935109f;
  static const float q = 0.78444488374548933f;
  union { float f; uint32_t i; } p = { 0.20363937680730309f };
  union { float f; uint32_t i; } r = { 0.015124940802184233f };
  union { float f; uint32_t i; } s = { -0.0032225901625579573f };

  union { float f; uint32_t i; } vx = { x };
  uint32_t sign = vx.i & 0x80000000;
  vx.i = vx.i & 0x7FFFFFFF;

  float qpprox = fouroverpi * x - fouroverpisq * x * vx.f;
  float qpproxsq = qpprox * qpprox;

  p.i |= sign;
  r.i |= sign;
  s.i ^= sign;

  return q * qpprox + qpproxsq * (p.f + qpproxsq * (r.f + qpproxsq * s.f));
}

static inline float
fastersin(float x)
{
  static const float fouroverpi = 1.2732395447351627f;
  static const float fouroverpisq = 0.40528473456935109f;
  static const float q = 0.77633023248007499f;
  union { float f; uint32_t i; } p = { 0.22308510060189463f };

  union { float f; uint32_t i; } vx = { x };
  uint32_t sign = vx.i & 0x80000000;
  vx.i &= 0x7FFFFFFF;

  float qpprox = fouroverpi * x - fouroverpisq * x * vx.f;

  p.i |= sign;

  return qpprox * (q + p.f * qpprox);
}

static inline float
fastsinfull(float x)
{
  static const float twopi = 6.2831853071795865f;
  static const float invtwopi = 0.15915494309189534f;

  int k = x * invtwopi;
  float half = (x < 0) ? -0.5f : 0.5f;
  return fastsin((half + k) * twopi - x);
}

static inline float
fastersinfull(float x)
{
  static const float twopi = 6.2831853071795865f;
  static const float invtwopi = 0.15915494309189534f;

  int k = x * invtwopi;
  float half = (x < 0) ? -0.5f : 0.5f;
  return fastersin((half + k) * twopi - x);
}

static inline float
fastcos(float x)
{
  static const float halfpi = 1.5707963267948966f;
  static const float halfpiminustwopi = -4.7123889803846899f;
  float offset = (x > halfpi) ? halfpiminustwopi : halfpi;
  return fastsin(x + offset);
}

static inline float
fastercos(float x)
{
  static const float twooverpi = 0.63661977236758134f;
  static const float p = 0.54641335845679634f;

  union { float f; uint32_t i; } vx = { x };
  vx.i &= 0x7FFFFFFF;

  float qpprox = 1.0f - twooverpi * vx.f;

  return qpprox + p * qpprox * (1.0f - qpprox * qpprox);
}

static inline float
fastcosfull(float x)
{
  static const float halfpi = 1.5707963267948966f;
  return fastsinfull(x + halfpi);
}

static inline float
fastercosfull(float x)
{
  static const float halfpi = 1.5707963267948966f;
  return fastersinfull(x + halfpi);
}

static inline float
fasttan(float x)
{
  static const float halfpi = 1.5707963267948966f;
  return fastsin(x) / fastsin(x + halfpi);
}

static inline float
fastertan(float x)
{
  return fastersin(x) / fastercos(x);
}

static inline float
fasttanfull(float x)
{
  static const float twopi = 6.2831853071795865f;
  static const float invtwopi = 0.15915494309189534f;

  int k = x * invtwopi;
  float half = (x < 0) ? -0.5f : 0.5f;
  float xnew = x - (half + k) * twopi;

  return fastsin(xnew) / fastcos(xnew);
}

static inline float
fastertanfull(float x)
{
  static const float twopi = 6.2831853071795865f;
  static const float invtwopi = 0.15915494309189534f;

  int k = x * invtwopi;
  float half = (x < 0) ? -0.5f : 0.5f;
  float xnew = x - (half + k) * twopi;

  return fastersin(xnew) / fastercos(xnew);
}

#ifdef __SSE2__

static inline v4sf
vfastsin(const v4sf x)
{
  const v4sf fouroverpi = v4sfl(1.2732395447351627f);
  const v4sf fouroverpisq = v4sfl(0.40528473456935109f);
  const v4sf q = v4sfl(0.78444488374548933f);
  const v4sf p = v4sfl(0.20363937680730309f);
  const v4sf r = v4sfl(0.015124940802184233f);
  const v4sf s = v4sfl(-0.0032225901625579573f);

  union { v4sf f; v4si i; } vx = { x };
  v4si sign = vx.i & v4sil(0x80000000);
  vx.i &= v4sil(0x7FFFFFFF);

  v4sf qpprox = fouroverpi * x - fouroverpisq * x * vx.f;
  v4sf qpproxsq = qpprox * qpprox;
  union { v4sf f; v4si i; } vy; vy.f = qpproxsq * (p + qpproxsq * (r + qpproxsq * s));
  vy.i ^= sign;

  return q * qpprox + vy.f;
}

static inline v4sf
vfastersin(const v4sf x)
{
  const v4sf fouroverpi = v4sfl(1.2732395447351627f);
  const v4sf fouroverpisq = v4sfl(0.40528473456935109f);
  const v4sf q = v4sfl(0.77633023248007499f);
  const v4sf plit = v4sfl(0.22308510060189463f);
  union { v4sf f; v4si i; } p = { plit };

  union { v4sf f; v4si i; } vx = { x };
  v4si sign = vx.i & v4sil(0x80000000);
  vx.i &= v4sil(0x7FFFFFFF);

  v4sf qpprox = fouroverpi * x - fouroverpisq * x * vx.f;

  p.i |= sign;

  return qpprox * (q + p.f * qpprox);
}

static inline v4sf
vfastsinfull(const v4sf x)
{
  const v4sf twopi = v4sfl(6.2831853071795865f);
  const v4sf invtwopi = v4sfl(0.15915494309189534f);

  v4si k = v4sf_to_v4si(x * invtwopi);

  v4sf ltzero = _mm_cmplt_ps(x, v4sfl(0.0f));
  v4sf half = _mm_or_ps(_mm_and_ps(ltzero, v4sfl(-0.5f)),
    _mm_andnot_ps(ltzero, v4sfl(0.5f)));

  return vfastsin((half + v4si_to_v4sf(k)) * twopi - x);
}

static inline v4sf
vfastersinfull(const v4sf x)
{
  const v4sf twopi = v4sfl(6.2831853071795865f);
  const v4sf invtwopi = v4sfl(0.15915494309189534f);

  v4si k = v4sf_to_v4si(x * invtwopi);

  v4sf ltzero = _mm_cmplt_ps(x, v4sfl(0.0f));
  v4sf half = _mm_or_ps(_mm_and_ps(ltzero, v4sfl(-0.5f)),
    _mm_andnot_ps(ltzero, v4sfl(0.5f)));

  return vfastersin((half + v4si_to_v4sf(k)) * twopi - x);
}

static inline v4sf
vfastcos(const v4sf x)
{
  const v4sf halfpi = v4sfl(1.5707963267948966f);
  const v4sf halfpiminustwopi = v4sfl(-4.7123889803846899f);
  v4sf lthalfpi = _mm_cmpnlt_ps(x, halfpi);
  v4sf offset = _mm_or_ps(_mm_and_ps(lthalfpi, halfpiminustwopi),
    _mm_andnot_ps(lthalfpi, halfpi));
  return vfastsin(x + offset);
}

static inline v4sf
vfastercos(v4sf x)
{
  const v4sf twooverpi = v4sfl(0.63661977236758134f);
  const v4sf p = v4sfl(0.54641335845679634);

  v4sf vx = v4sf_fabs(x);
  v4sf qpprox = v4sfl(1.0f) - twooverpi * vx;

  return qpprox + p * qpprox * (v4sfl(1.0f) - qpprox * qpprox);
}

static inline v4sf
vfastcosfull(const v4sf x)
{
  const v4sf halfpi = v4sfl(1.5707963267948966f);
  return vfastsinfull(x + halfpi);
}

static inline v4sf
vfastercosfull(const v4sf x)
{
  const v4sf halfpi = v4sfl(1.5707963267948966f);
  return vfastersinfull(x + halfpi);
}

static inline v4sf
vfasttan(const v4sf x)
{
  const v4sf halfpi = v4sfl(1.5707963267948966f);
  return vfastsin(x) / vfastsin(x + halfpi);
}

static inline v4sf
vfastertan(const v4sf x)
{
  return vfastersin(x) / vfastercos(x);
}

static inline v4sf
vfasttanfull(const v4sf x)
{
  const v4sf twopi = v4sfl(6.2831853071795865f);
  const v4sf invtwopi = v4sfl(0.15915494309189534f);

  v4si k = v4sf_to_v4si(x * invtwopi);

  v4sf ltzero = _mm_cmplt_ps(x, v4sfl(0.0f));
  v4sf half = _mm_or_ps(_mm_and_ps(ltzero, v4sfl(-0.5f)),
    _mm_andnot_ps(ltzero, v4sfl(0.5f)));
  v4sf xnew = x - (half + v4si_to_v4sf(k)) * twopi;

  return vfastsin(xnew) / vfastcos(xnew);
}

static inline v4sf
vfastertanfull(const v4sf x)
{
  const v4sf twopi = v4sfl(6.2831853071795865f);
  const v4sf invtwopi = v4sfl(0.15915494309189534f);

  v4si k = v4sf_to_v4si(x * invtwopi);

  v4sf ltzero = _mm_cmplt_ps(x, v4sfl(0.0f));
  v4sf half = _mm_or_ps(_mm_and_ps(ltzero, v4sfl(-0.5f)),
    _mm_andnot_ps(ltzero, v4sfl(0.5f)));
  v4sf xnew = x - (half + v4si_to_v4sf(k)) * twopi;

  return vfastersin(xnew) / vfastercos(xnew);
}

#endif //__SSE2__


namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // min_max_range
  ////////////////////////////////////////////////////////////////////////////
  template <typename T>
  struct min_max_range
  {
    T min, max;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Some macros
  ////////////////////////////////////////////////////////////////////////////

  ////////////////////////////////////////////////////////////////////////////
  // pi
  ////////////////////////////////////////////////////////////////////////////
  constexpr double pi = 3.1415926535897932384626433832795;

  ////////////////////////////////////////////////////////////////////////////
  // abs (we need ot here because std::abs may not be constexpr)
  ////////////////////////////////////////////////////////////////////////////
  template <typename T>
  constexpr T abs(T x)
  {
    return (x < T(0)) ? -x : x;
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast trigonometric functions
  //
  // fast sine variants: x in [ -pi, pi ]
  // fast cosine variants: x in [ -pi, pi ]
  // fast tangent variants: x in [ -pi/2, pi/2 ]
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_tan(float x)
  {
    return fasttan(x);
  }

  inline float faster_tan(float x)
  {
    return fastertan(x);
  }

  inline float fast_sin(float x)
  {
    return fastsin(x);
  }

  inline float faster_sin(float x)
  {
    return fastersin(x);
  }

  inline float fast_cos(float x)
  {
    return fastcos(x);
  }

  inline float faster_cos(float x)
  {
    return fastercos(x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // fast pade-approximation of the tanh function (x should be: -3 <= x <= 3)
  ////////////////////////////////////////////////////////////////////////////
  constexpr float fast_rational_tanh(float x)
  {
    return x * (27 + x * x) / (27 + 9 * x * x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast exp
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_exp(float x)
  {
    return fastexp(x);
  }

  inline float faster_exp(float x)
  {
    return fasterexp(x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast exp Taylor series approximations (from http://www.musicdsp.org/)
  ////////////////////////////////////////////////////////////////////////////
  constexpr float fast_exp3(float x)
  {
    return (6 + x * (6 + x * (3 + x))) * 0.16666666f;
  }

  constexpr float fast_exp4(float x)
  {
    return (24 + x * (24 + x * (12 + x * (4 + x)))) * 0.041666666f;
  }

  constexpr float fast_exp5(float x)
  {
    return (120 + x * (120 + x * (60 + x * (20 + x * (5 + x))))) * 0.0083333333f;
  }

  constexpr float fast_exp6(float x)
  {
    return (720 + x * (720 + x * (360 + x * (120 + x * (30 + x * (6 + x)))))) * 0.0013888888f;
  }

  constexpr float fast_exp7(float x)
  {
    return (5040 + x * (5040 + x * (2520 + x *
      (840 + x * (210 + x * (42 + x * (7 + x))))))) * 0.00019841269f;
  }

  constexpr float fast_exp8(float x)
  {
    return (40320 + x * (40320 + x * (20160 + x * (6720 + x *
      (1680 + x * (336 + x * (56 + x * (8 + x)))))))) * 2.4801587301e-5f;
  }

  constexpr float fast_exp9(float x)
  {
    return (362880 + x * (362880 + x * (181440 + x * (60480 + x *
      (15120 + x * (3024 + x * (504 + x * (72 + x * (9 + x))))))))) * 2.75573192e-6f;
  }

  ////////////////////////////////////////////////////////////////////////////
  // linear interpolation: Interpolates a value linearly between y1 and y2
  // given mu. If mu is 0, the result is y1. If mu is 1, then the result is
  // y2.
  ////////////////////////////////////////////////////////////////////////////
  constexpr float linear_interpolate(float y1, float y2, float mu)
  {
    return y1 + mu * (y2 - y1);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast reciprocal. See http://tinyurl.com/lgmnsyg. We want to avoid
  // multiplication. So instead of 1.0f/val, we use this function which
  // is inaccurate, but fast substitute. It works by negating the exponent
  // which is assumed to be IEEE754.
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_inverse(float val)
  {
    auto x = reinterpret_cast<std::int32_t&>(val);
    x = 0x7EF311C2 - x;
    return reinterpret_cast<float&>(x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast division using multiplication and fast_inverse
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_div(float a, float b)
  {
    return a * fast_inverse(b);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast log
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_log(float x)
  {
    return fastlog(x);
  }

  inline float faster_log(float x)
  {
    return fasterlog(x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast log2
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_log2(float x)
  {
    return fastlog2(x);
  }

  inline float faster_log2(float x)
  {
    return fasterlog2(x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast pow2
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_pow2(float x)
  {
    return fastpow2(x);
  }

  inline float faster_pow2(float x)
  {
    return fasterpow2(x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast sqrt
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_sqrt(float x)
  {
    return fast_pow2(fast_log2(x) / 2);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast log10
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_log10(float x)
  {
    return 0.301029995663981f * fast_log2(x);
  }

  inline float faster_log10(float x)
  {
    return 0.301029995663981f * faster_log2(x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast pow10
  ////////////////////////////////////////////////////////////////////////////
  inline float fast_pow10(float x)
  {
    return fastpow(10, x);
  }

  inline float faster_pow10(float x)
  {
    return fasterpow(10, x);
  }

  ////////////////////////////////////////////////////////////////////////////
  // Fast random number generator
  ////////////////////////////////////////////////////////////////////////////
  inline int fast_rand()
  {
    static unsigned seed = 87263876;
    seed = (214013 * seed + 2531011);
    return (seed >> 16) & 0x7FFF;
  }

  ////////////////////////////////////////////////////////////////////////////
  // abs_within
  ////////////////////////////////////////////////////////////////////////////
  inline bool abs_within(float a, float b, float eps)
  {
    return abs(a - b) <= eps;
  }

  inline bool abs_within(int a, int b, int eps)
  {
    return abs(a - b) <= eps;
  }

  ////////////////////////////////////////////////////////////////////////////
  // rel_within
  ////////////////////////////////////////////////////////////////////////////
  inline bool rel_within(float a, float b, float eps)
  {
    return abs(a - b) <= eps * std::max(abs(a), abs(b));
  }
}

namespace cycfi::q::concepts
{
  template <typename T>
  concept Arithmetic = std::integral<T> || std::floating_point<T>;

  template <typename T>
  concept IndexableContainer = requires(T & x, std::size_t i)
  {
    { x[i] } -> std::convertible_to<typename T::value_type>;
    { x.size() } -> std::convertible_to<std::size_t>;
  };

  template <typename T>
  concept RandomAccessIteratable =
    std::random_access_iterator<typename T::iterator> &&
    requires(T & c)
  {
    { c.begin() } -> std::same_as<typename T::iterator>;
    { c.end() } -> std::same_as<typename T::iterator>;
  };
}

namespace cycfi::q
{
  namespace concepts
  {
    template <typename A, typename B>
    concept SameUnit = std::same_as<typename A::unit_type, typename B::unit_type>;
  }

  ////////////////////////////////////////////////////////////////////////////
  // On binary operations `a + b` and `a - b`, where `a` and `b` conform to
  // the `SameUnit` concept (see above), the resuling type will be whichever
  // has the `value_type` of `decltype(a.ref + b.rep)`, else if both
  // operands are promoted, then whichever has the larger `value_type` will
  // be chosen.
  //
  // Promotion logic:
  //
  // If decltype(a.rep + b.rep) is the same as a.rep choose A. Else if
  // decltype(a.rep + b.rep) is the same as b.rep choose B. Else if the
  // sizeof(a.rep) >= sizeof(b.rep) choose A. Else, choose B.
  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
  using promote_unit =
    std::conditional_t <
    std::is_same_v<
    decltype(typename A::value_type{} + typename B::value_type{})
    , typename A::value_type >
    , A
    , std::conditional_t <
    std::is_same_v<
    decltype(typename B::value_type{} + typename A::value_type{})
    , typename B::value_type >
    , B
    , std::conditional_t<
    sizeof(typename B::value_type) >= sizeof(typename B::value_type)
    , A
    , B
    >
    >
    > ;

  struct direct_unit_type {};
  constexpr static direct_unit_type direct_unit = {};

  ////////////////////////////////////////////////////////////////////////////
  // unit: Unit abstraction and encapsulation
  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Derived>
  struct unit
  {
    using derived_type = Derived;
    using value_type = T;
    // Temporary constructor. This is not
    // marked deprecated because we will use
    // this for now.
    constexpr                     unit(T val, direct_unit_type) : rep(val) {}

    constexpr                     unit(T val) : rep(val) {}
    constexpr                     unit(unit const&) = default;
    constexpr                     unit(unit&&) = default;

    constexpr unit& operator=(unit const&) = default;
    constexpr unit& operator=(unit&&) = default;

    constexpr derived_type        operator+() const;
    constexpr derived_type        operator-() const;

    constexpr derived_type& operator+=(unit rhs);
    constexpr derived_type& operator+=(concepts::Arithmetic auto rhs);

    constexpr derived_type& operator-=(unit rhs);
    constexpr derived_type& operator-=(concepts::Arithmetic auto rhs);

    constexpr derived_type& operator*=(concepts::Arithmetic auto b);
    constexpr derived_type& operator/=(concepts::Arithmetic auto b);

    constexpr derived_type const& derived() const;
    constexpr derived_type& derived();

    template <typename U> requires concepts::SameUnit<unit, U>
    constexpr unit& operator=(U rhs) { rep = rhs.rep; return derived(); }

    template <typename U> requires concepts::SameUnit<unit, U>
    constexpr unit& operator+=(U rhs) { rep += rhs.rep; return derived(); }

    template <typename U> requires concepts::SameUnit<unit, U>
    constexpr unit& operator-=(U rhs) { rep -= rhs.rep; return derived(); }

    T rep;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Free functions
  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator==(A a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator==(A a, unit<B, Derived> b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator==(unit<A, Derived> a, B b);

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator!=(A a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator!=(A a, unit<B, Derived> b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator!=(unit<A, Derived> a, B b);

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator<(A a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator<(A a, unit<B, Derived> b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator<(unit<A, Derived> a, B b);

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator<=(A a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator<=(A a, unit<B, Derived> b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator<=(unit<A, Derived> a, B b);

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator>(A a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator>(A a, unit<B, Derived> b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator>(unit<A, Derived> a, B b);

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator>=(A a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator>=(A a, unit<B, Derived> b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator>=(unit<A, Derived> a, B b);

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr promote_unit<A, B> operator+(A a, B b);

  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr promote_unit<A, B> operator-(A a, B b);

  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr typename promote_unit<A, B>::value_type operator/(A a, B b);

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr Derived operator+(A a, unit<B, Derived> b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr Derived operator-(A a, unit<B, Derived> b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr Derived operator*(A a, unit<B, Derived> b);

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr Derived operator+(unit<A, Derived> a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr Derived operator-(unit<A, Derived> a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr Derived operator*(unit<A, Derived> a, B b);

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr Derived operator/(unit<A, Derived> a, B b);

  ////////////////////////////////////////////////////////////////////////////
  // Inlines
  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Derived>
  constexpr typename unit<T, Derived>::derived_type const&
    unit<T, Derived>::derived() const
  {
    return *static_cast<derived_type const*>(this);
  }

  template <typename T, typename Derived>
  constexpr typename unit<T, Derived>::derived_type&
    unit<T, Derived>::derived()
  {
    return *static_cast<derived_type*>(this);
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Derived>
  constexpr Derived unit<T, Derived>::operator+() const
  {
    return derived();
  }

  template <typename T, typename Derived>
  constexpr Derived unit<T, Derived>::operator-() const
  {
    return derived_type{ -rep, direct_unit };
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator==(A a, B b)
  {
    return a.rep == b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator==(A a, unit<B, Derived> b)
  {
    return a == b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator==(unit<A, Derived> a, B b)
  {
    return a.rep == b;
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator!=(A a, B b)
  {
    return a.rep != b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator!=(A a, unit<B, Derived> b)
  {
    return a != b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator!=(unit<A, Derived> a, B b)
  {
    return a.rep != b;
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator<(A a, B b)
  {
    return a.rep < b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator<(A a, unit<B, Derived> b)
  {
    return a < b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator<(unit<A, Derived> a, B b)
  {
    return a.rep < b;
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator<=(A a, B b)
  {
    return a.rep <= b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator<=(A a, unit<B, Derived> b)
  {
    return a <= b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator<=(unit<A, Derived> a, B b)
  {
    return a.rep <= b;
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator>(A a, B b)
  {
    return a.rep > b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator>(A a, unit<B, Derived> b)
  {
    return a > b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator>(unit<A, Derived> a, B b)
  {
    return a.rep > b;
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr bool operator>=(A a, B b)
  {
    return a.rep >= b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr bool operator>=(A a, unit<B, Derived> b)
  {
    return a >= b.rep;
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr bool operator>=(unit<A, Derived> a, B b)
  {
    return a.rep >= b;
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Derived>
  constexpr Derived& unit<T, Derived>::operator+=(unit<T, Derived> rhs)
  {
    rep += rhs.rep;
    return derived();
  }

  template <typename T, typename Derived>
  constexpr Derived& unit<T, Derived>::operator+=(concepts::Arithmetic auto rhs)
  {
    rep += rhs;
    return derived();
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Derived>
  constexpr Derived& unit<T, Derived>::operator-=(unit<T, Derived> rhs)
  {
    rep -= rhs.rep;
    return derived();
  }

  template <typename T, typename Derived>
  constexpr Derived& unit<T, Derived>::operator-=(concepts::Arithmetic auto rhs)
  {
    rep -= rhs;
    return derived();
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Derived>
  constexpr Derived& unit<T, Derived>::operator*=(concepts::Arithmetic auto rhs)
  {
    rep *= rhs;
    return derived();
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Derived>
  constexpr Derived& unit<T, Derived>::operator/=(concepts::Arithmetic auto rhs)
  {
    rep /= rhs;
    return derived();
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr promote_unit<A, B> operator+(A a, B b)
  {
    return promote_unit<A, B>{a.rep + b.rep, direct_unit};
  }

  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr promote_unit<A, B> operator-(A a, B b)
  {
    return promote_unit<A, B>{a.rep - b.rep, direct_unit};
  }

  template <typename A, typename B>
    requires concepts::SameUnit<A, B>
  constexpr typename promote_unit<A, B>::value_type
    operator/(A a, B b)
  {
    return a.rep / b.rep;
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr Derived operator+(A a, unit<B, Derived> b)
  {
    return Derived{ a + b.rep, direct_unit };
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr Derived operator-(A a, unit<B, Derived> b)
  {
    return Derived{ a - b.rep, direct_unit };
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr Derived operator*(A a, unit<B, Derived> b)
  {
    return Derived{ a * b.rep, direct_unit };
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<A>
  constexpr Derived operator/(A a, unit<B, Derived> b)
  {
    return Derived{ a / b.rep, direct_unit };
  }

  ////////////////////////////////////////////////////////////////////////////
  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr Derived operator+(unit<A, Derived> a, B b)
  {
    return Derived{ a.rep + b, direct_unit };
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr Derived operator-(unit<A, Derived> a, B b)
  {
    return Derived{ a.rep - b, direct_unit };
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr Derived operator*(unit<A, Derived> a, B b)
  {
    return Derived{ a.rep * b, direct_unit };
  }

  template <typename A, typename B, typename Derived>
    requires concepts::Arithmetic<B>
  constexpr Derived operator/(unit<A, Derived> a, B b)
  {
    return Derived{ a.rep / b, direct_unit };
  }
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  struct duration_unit;

  struct duration : unit<double, duration>
  {
    using base_type = unit<double, duration>;
    using base_type::base_type;
    using unit_type = duration_unit;
  };

  // Free functions
  constexpr double  as_double(duration d);
  constexpr float   as_float(duration d);

  ////////////////////////////////////////////////////////////////////////////
  constexpr double as_double(duration d)
  {
    return d.rep;
  }

  constexpr float as_float(duration d)
  {
    return d.rep;
  }
}

namespace cycfi::q
{
  struct period;

  ////////////////////////////////////////////////////////////////////////////
  struct frequency_unit;

  struct frequency : unit<double, frequency>
  {
    using base_type = unit<double, frequency>;
    using base_type::base_type;
    using unit_type = frequency_unit;

    constexpr q::period  period() const;
  };

  // Free functions
  constexpr double  as_double(frequency f);
  constexpr float   as_float(frequency f);
}

namespace cycfi::q
{
  struct frequency;

  ////////////////////////////////////////////////////////////////////////////
  struct period : duration
  {
    using duration::duration;

    constexpr explicit   period(duration d);
    constexpr explicit   period(frequency f);
  };

  ////////////////////////////////////////////////////////////////////////////
  constexpr period::period(duration d)
    : duration(d)
  {
  }

  constexpr period::period(frequency f)
    : duration(1.0 / f.rep)
  {
  }
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  constexpr q::period frequency::period() const
  {
    return q::period{ 1.0 / rep };
  }

  constexpr double as_double(frequency f)
  {
    return f.rep;
  }

  constexpr float as_float(frequency f)
  {
    return f.rep;
  }
}

namespace cycfi::q::detail
{
  ////////////////////////////////////////////////////////////////////////////
  // dB lookup table
  ////////////////////////////////////////////////////////////////////////////
  constexpr float inv_db_table[] =
  {
     1.000000000000000f, 1.011579454259899f, 1.023292992280754f, 1.035142166679344f, 1.047128548050900f, 1.059253725177289f, 1.071519305237606f, 1.083926914021204f,
     1.096478196143185f, 1.109174815262401f, 1.122018454301963f, 1.135010815672315f, 1.148153621496883f, 1.161448613840343f, 1.174897554939530f, 1.188502227437018f,
     1.202264434617413f, 1.216186000646368f, 1.230268770812382f, 1.244514611771385f, 1.258925411794167f, 1.273503081016662f, 1.288249551693134f, 1.303166778452299f,
     1.318256738556407f, 1.333521432163324f, 1.348962882591654f, 1.364583136588924f, 1.380384264602885f, 1.396368361055938f, 1.412537544622754f, 1.428893958511103f,
     1.445439770745927f, 1.462177174456718f, 1.479108388168207f, 1.496235656094433f, 1.513561248436208f, 1.531087461682030f, 1.548816618912481f, 1.566751070108149f,
     1.584893192461114f, 1.603245390690041f, 1.621810097358930f, 1.640589773199539f, 1.659586907437561f, 1.678804018122560f, 1.698243652461744f, 1.717908387157588f,
     1.737800828749375f, 1.757923613958693f, 1.778279410038923f, 1.798870915128788f, 1.819700858609983f, 1.840772001468956f, 1.862087136662868f, 1.883649089489801f,
     1.905460717963247f, 1.927524913190936f, 1.949844599758045f, 1.972422736114854f, 1.995262314968880f, 2.018366363681561f, 2.041737944669529f, 2.065380155810530f,
     2.089296130854040f, 2.113489039836647f, 2.137962089502232f, 2.162718523727020f, 2.187761623949553f, 2.213094709605638f, 2.238721138568339f, 2.264644307593060f,
     2.290867652767773f, 2.317394649968478f, 2.344228815319922f, 2.371373705661655f, 2.398832919019490f, 2.426610095082415f, 2.454708915685031f, 2.483133105295571f,
     2.511886431509580f, 2.540972705549305f, 2.570395782768863f, 2.600159563165272f, 2.630267991895382f, 2.660725059798810f, 2.691534803926916f, 2.722701308077912f,
     2.754228703338166f, 2.786121168629770f, 2.818382931264454f, 2.851018267503910f, 2.884031503126606f, 2.917427014001167f, 2.951209226666386f, 2.985382618917960f,
     3.019951720402016f, 3.054921113215513f, 3.090295432513591f, 3.126079367123955f, 3.162277660168380f, 3.198895109691398f, 3.235936569296283f, 3.273406948788382f,
     3.311311214825911f, 3.349654391578277f, 3.388441561392026f, 3.427677865464504f, 3.467368504525317f, 3.507518739525680f, 3.548133892335755f, 3.589219346450053f,
     3.630780547701014f, 3.672823004980846f, 3.715352290971725f, 3.758374042884441f, 3.801893963205611f, 3.845917820453535f, 3.890451449942806f, 3.935500754557774f,
     3.981071705534972f, 4.027170343254591f, 4.073802778041127f, 4.120975190973302f, 4.168693834703354f, 4.216965034285822f, 4.265795188015927f, 4.315190768277652f,
     4.365158322401660f, 4.415704473533125f, 4.466835921509632f, 4.518559443749224f, 4.570881896148751f, 4.623810213992603f, 4.677351412871983f, 4.731512589614805f,
     4.786300923226384f, 4.841723675840994f, 4.897788193684462f, 4.954501908047901f, 5.011872336272722f, 5.069907082747043f, 5.128613839913648f, 5.188000389289611f,
     5.248074602497725f, 5.308844442309884f, 5.370317963702527f, 5.432503314924332f, 5.495408738576246f, 5.559042572704036f, 5.623413251903491f, 5.688529308438414f,
     5.754399373371569f, 5.821032177708714f, 5.888436553555890f, 5.956621435290105f, 6.025595860743578f, 6.095368972401692f, 6.165950018614822f, 6.237348354824193f,
     6.309573444801933f, 6.382634861905488f, 6.456542290346556f, 6.531305526474723f, 6.606934480075959f, 6.683439175686146f, 6.760829753919817f, 6.839116472814292f,
     6.918309709189364f, 6.998419960022734f, 7.079457843841379f, 7.161434102129020f, 7.244359600749901f, 7.328245331389041f, 7.413102413009175f, 7.498942093324558f,
     7.585775750291837f, 7.673614893618190f, 7.762471166286917f, 7.852356346100718f, 7.943282347242816f, 8.035261221856173f, 8.128305161640993f, 8.222426499470712f,
     8.317637711026711f, 8.413951416451951f, 8.511380382023766f, 8.609937521846007f, 8.709635899560805f, 8.810488730080140f, 8.912509381337454f, 9.015711376059569f,
     9.120108393559097f, 9.225714271547631f, 9.332543007969910f, 9.440608762859233f, 9.549925860214358f, 9.660508789898133f, 9.772372209558107f, 9.885530946569387f,
     10.000000000000000f, 10.115794542598982f, 10.232929922807541f, 10.351421666793437f, 10.471285480508996f, 10.592537251772887f, 10.715193052376065f, 10.839269140212034f,
     10.964781961431852f, 11.091748152624008f, 11.220184543019636f, 11.350108156723149f, 11.481536214968829f, 11.614486138403427f, 11.748975549395297f, 11.885022274370183f,
     12.022644346174131f, 12.161860006463678f, 12.302687708123818f, 12.445146117713850f, 12.589254117941675f, 12.735030810166617f, 12.882495516931343f, 13.031667784522993f,
     13.182567385564074f, 13.335214321633240f, 13.489628825916533f, 13.645831365889245f, 13.803842646028846f, 13.963683610559377f, 14.125375446227540f, 14.288939585111029f,
     14.454397707459272f, 14.621771744567184f, 14.791083881682072f, 14.962356560944336f, 15.135612484362079f, 15.310874616820303f, 15.488166189124811f, 15.667510701081493f,
     15.848931924611133f, 16.032453906900418f, 16.218100973589298f, 16.405897731995395f, 16.595869074375607f, 16.788040181225607f, 16.982436524617444f, 17.179083871575884f,
     17.378008287493753f, 17.579236139586932f, 17.782794100389228f, 17.988709151287875f, 18.197008586099834f, 18.407720014689556f, 18.620871366628677f, 18.836490894898002f,
     19.054607179632473f, 19.275249131909355f, 19.498445997580454f, 19.724227361148536f, 19.952623149688797f, 20.183663636815606f, 20.417379446695296f, 20.653801558105290f,
     20.892961308540396f, 21.134890398366466f, 21.379620895022324f, 21.627185237270201f, 21.877616239495531f, 22.130947096056374f, 22.387211385683401f, 22.646443075930595f,
     22.908676527677734f, 23.173946499684785f, 23.442288153199225f, 23.713737056616552f, 23.988329190194900f, 24.266100950824157f, 24.547089156850298f, 24.831331052955704f,
     25.118864315095795f, 25.409727055493050f, 25.703957827688633f, 26.001595631652723f, 26.302679918953814f, 26.607250597988099f, 26.915348039269155f, 27.227013080779127f,
     27.542287033381662f, 27.861211686297708f, 28.183829312644534f, 28.510182675039097f, 28.840315031266059f, 29.174270140011675f, 29.512092266663856f, 29.853826189179603f,
     30.199517204020161f, 30.549211132155140f, 30.902954325135905f, 31.260793671239558f, 31.622776601683793f, 31.988951096913972f, 32.359365692962825f, 32.734069487883815f,
     33.113112148259113f, 33.496543915782759f, 33.884415613920254f, 34.276778654645028f, 34.673685045253166f, 35.075187395256798f, 35.481338923357548f, 35.892193464500515f,
     36.307805477010142f, 36.728230049808460f, 37.153522909717260f, 37.583740428844415f, 38.018939632056124f, 38.459178204535355f, 38.904514499428068f, 39.355007545577742f,
     39.810717055349734f, 40.271703432545905f, 40.738027780411279f, 41.209751909733022f, 41.686938347033553f, 42.169650342858226f, 42.657951880159253f, 43.151907682776525f,
     43.651583224016584f, 44.157044735331255f, 44.668359215096302f, 45.185594437492242f, 45.708818961487495f, 46.238102139926035f, 46.773514128719810f, 47.315125896148054f,
     47.863009232263828f, 48.417236758409942f, 48.977881936844611f, 49.545019080479030f, 50.118723362727224f, 50.699070827470443f, 51.286138399136483f, 51.880003892896120f,
     52.480746024977257f, 53.088444423098849f, 53.703179637025272f, 54.325033149243332f, 54.954087385762456f, 55.590425727040369f, 56.234132519034908f, 56.885293084384131f,
     57.543993733715695f, 58.210321777087131f, 58.884365535558899f, 59.566214352901035f, 60.255958607435780f, 60.953689724016904f, 61.659500186148222f, 62.373483548241914f,
     63.095734448019329f, 63.826348619054862f, 64.565422903465560f, 65.313055264747234f, 66.069344800759609f, 66.834391756861450f, 67.608297539198190f, 68.391164728142925f,
     69.183097091893657f, 69.984199600227342f, 70.794578438413808f, 71.614341021290201f, 72.443596007499025f, 73.282453313890400f, 74.131024130091774f, 74.989420933245583f,
     75.857757502918361f, 76.736148936181891f, 77.624711662869160f, 78.523563461007186f, 79.432823472428140f, 80.352612218561731f, 81.283051616409907f, 82.224264994707127f,
     83.176377110267083f, 84.139514164519511f, 85.113803820237635f, 86.099375218460082f, 87.096358995608057f, 88.104887300801423f, 89.125093813374548f, 90.157113760595706f,
     91.201083935590972f, 92.257142715476334f, 93.325430079699103f, 94.406087628592360f, 95.499258602143584f, 96.605087898981353f, 97.723722095581067f, 98.855309465693907f,
     100.000000000000000f, 101.157945425989823f, 102.329299228075357f, 103.514216667934406f, 104.712854805089961f, 105.925372517728860f, 107.151930523760598f, 108.392691402120391f,
     109.647819614318507f, 110.917481526240095f, 112.201845430196300f, 113.501081567231552f, 114.815362149688283f, 116.144861384034257f, 117.489755493952913f, 118.850222743701892f,
     120.226443461741312f, 121.618600064636794f, 123.026877081238112f, 124.451461177138569f, 125.892541179416753f, 127.350308101666158f, 128.824955169313370f, 130.316677845230004f,
     131.825673855640730f, 133.352143216332394f, 134.896288259165317f, 136.458313658892394f, 138.038426460288520f, 139.636836105593773f, 141.253754462275396f, 142.889395851110208f,
     144.543977074592789f, 146.217717445671838f, 147.910838816820728f, 149.623565609443290f, 151.356124843620876f, 153.108746168203027f, 154.881661891248115f, 156.675107010814855f,
     158.489319246111421f, 160.324539069004175f, 162.181009735892985f, 164.058977319953868f, 165.958690743756136f, 167.880401812256054f, 169.824365246174438f, 171.790838715758753f,
     173.780082874937619f, 175.792361395869307f, 177.827941003892278f, 179.887091512878726f, 181.970085860998267f, 184.077200146895620f, 186.208713666286741f, 188.364908948980002f,
     190.546071796324640f, 192.752491319093650f, 194.984459975804555f, 197.242273611485331f, 199.526231496887874f, 201.836636368156178f, 204.173794466952955f, 206.538015581052917f,
     208.929613085403872f, 211.348903983664769f, 213.796208950223246f, 216.271852372702000f, 218.776162394955179f, 221.309470960563885f, 223.872113856834005f, 226.464430759305969f,
     229.086765276777243f, 231.739464996847971f, 234.422881531992260f, 237.137370566165515f, 239.883291901948979f, 242.661009508241420f, 245.470891568503106f, 248.313310529557043f,
     251.188643150957972f, 254.097270554930390f, 257.039578276886459f, 260.015956316527195f, 263.026799189538167f, 266.072505979880873f, 269.153480392691677f, 272.270130807791247f,
     275.422870333816604f, 278.612116862976961f, 281.838293126445478f, 285.101826750390956f, 288.403150312660557f, 291.742701400116573f, 295.120922666638705f, 298.538261891796026f,
     301.995172040201624f, 305.492111321551249f, 309.029543251359200f, 312.607936712395599f, 316.227766016837961f, 319.889510969139735f, 323.593656929628082f, 327.340694878838292f,
     331.131121482591084f, 334.965439157827575f, 338.844156139202369f, 342.767786546450452f, 346.736850452531655f, 350.751873952567962f, 354.813389233575322f, 358.921934645005365f,
     363.078054770101403f, 367.282300498084624f, 371.535229097172419f, 375.837404288444304f, 380.189396320561286f, 384.591782045353511f, 389.045144994280463f, 393.550075455777630f,
     398.107170553497326f, 402.717034325459053f, 407.380277804112609f, 412.097519097330405f, 416.869383470335492f, 421.696503428582275f, 426.579518801592542f, 431.519076827764991f,
     436.515832240166105f, 441.570447353312545f, 446.683592150963023f, 451.855944374922160f, 457.088189614875205f, 462.381021399260362f, 467.735141287198132f, 473.151258961480266f,
     478.630092322638518f, 484.172367584099391f, 489.778819368446136f, 495.450190804790054f, 501.187233627272462f, 506.990708274704446f, 512.861383991364846f, 518.800038928960930f,
     524.807460249772816f, 530.884444230988493f, 537.031796370252664f, 543.250331492433020f, 549.540873857624774f, 555.904257270403718f, 562.341325190349039f, 568.852930843841364f,
     575.439937337156607f, 582.103217770871538f, 588.843655355589021f, 595.662143529010336f, 602.559586074357526f, 609.536897240169310f, 616.595001861482160f, 623.734835482419157f,
     630.957344480193001f, 638.263486190548974f, 645.654229034655600f, 653.130552647472314f, 660.693448007595748f, 668.343917568614870f, 676.082975391981904f, 683.911647281429282f,
     691.830970918936259f, 699.841996002273845f, 707.945784384138051f, 716.143410212902040f, 724.435960074989907f, 732.824533138904485f, 741.310241300917710f, 749.894209332455830f,
     758.577575029183549f, 767.361489361818599f, 776.247116628691970f, 785.235634610071770f, 794.328234724281288f, 803.526122185616941f, 812.830516164099549f, 822.242649947071186f,
     831.763771102670830f, 841.395141645194713f, 851.138038202376833f, 860.993752184600794f, 870.963589956080568f, 881.048873008013743f, 891.250938133745876f, 901.571137605957006f,
     912.010839355909638f, 922.571427154762887f, 933.254300796991515f, 944.060876285923541f, 954.992586021435955f, 966.050878989813100f, 977.237220955811154f, 988.553094656939038f,
     1000.000000000000000f, 1011.579454259898284f, 1023.292992280753651f, 1035.142166679344200f, 1047.128548050899553f, 1059.253725177288743f, 1071.519305237605977f, 1083.926914021203856f,
     1096.478196143185187f, 1109.174815262400898f, 1122.018454301962947f, 1135.010815672315402f, 1148.153621496882806f, 1161.448613840342659f, 1174.897554939529073f, 1188.502227437018973f,
     1202.264434617413144f, 1216.186000646367802f, 1230.268770812381035f, 1244.514611771385717f, 1258.925411794167530f, 1273.503081016661554f, 1288.249551693133526f, 1303.166778452299923f,
     1318.256738556407527f, 1333.521432163324107f, 1348.962882591653170f, 1364.583136588923935f, 1380.384264602885196f, 1396.368361055937612f, 1412.537544622754012f, 1428.893958511102255f,
     1445.439770745927945f, 1462.177174456718376f, 1479.108388168207284f, 1496.235656094432898f, 1513.561248436208643f, 1531.087461682030380f, 1548.816618912481090f, 1566.751070108148497f,
     1584.893192461114040f, 1603.245390690041859f, 1621.810097358929852f, 1640.589773199538740f, 1659.586907437561422f, 1678.804018122560592f, 1698.243652461744205f, 1717.908387157587640f,
     1737.800828749376251f, 1757.923613958693068f, 1778.279410038922833f, 1798.870915128787374f, 1819.700858609982561f, 1840.772001468956432f, 1862.087136662867579f, 1883.649089489800190f,
     1905.460717963246225f, 1927.524913190936559f, 1949.844599758045433f, 1972.422736114853478f, 1995.262314968878854f, 2018.366363681561779f, 2041.737944669529497f, 2065.380155810529232f,
     2089.296130854038893f, 2113.489039836647862f, 2137.962089502232629f, 2162.718523727020056f, 2187.761623949551904f, 2213.094709605638855f, 2238.721138568339939f, 2264.644307593059693f,
     2290.867652767772597f, 2317.394649968479825f, 2344.228815319922887f, 2371.373705661655094f, 2398.832919019489964f, 2426.610095082414318f, 2454.708915685031116f, 2483.133105295570203f,
     2511.886431509579779f, 2540.972705549304010f, 2570.395782768864592f, 2600.159563165272175f, 2630.267991895381329f, 2660.725059798808616f, 2691.534803926916538f, 2722.701308077912927f,
     2754.228703338166270f, 2786.121168629769272f, 2818.382931264454783f, 2851.018267503909556f, 2884.031503126605912f, 2917.427014001165844f, 2951.209226666387167f, 2985.382618917960372f,
     3019.951720402016235f, 3054.921113215512378f, 3090.295432513592004f, 3126.079367123955763f, 3162.277660168379498f, 3198.895109691397010f, 3235.936569296281050f, 3273.406948788383033f,
     3311.311214825911065f, 3349.654391578275863f, 3388.441561392024141f, 3427.677865464504521f, 3467.368504525316894f, 3507.518739525679393f, 3548.133892335752989f, 3589.219346450053308f,
     3630.780547701013802f, 3672.823004980846235f, 3715.352290971723960f, 3758.374042884443043f, 3801.893963205612636f, 3845.917820453535569f, 3890.451449942804629f, 3935.500754557776418f,
     3981.071705534973262f, 4027.170343254590989f, 4073.802778041126203f, 4120.975190973304052f, 4168.693834703354696f, 4216.965034285822185f, 4265.795188015925305f, 4315.190768277650022f,
     4365.158322401661280f, 4415.704473533125565f, 4466.835921509630680f, 4518.559443749221828f, 4570.881896148751366f, 4623.810213992603167f, 4677.351412871980756f, 4731.512589614802891f,
     4786.300923226384839f, 4841.723675840994474f, 4897.788193684461476f, 4954.501908047900542f, 5011.872336272725079f, 5069.907082747044115f, 5128.613839913648007f, 5188.000389289609302f,
     5248.074602497727938f, 5308.844442309884471f, 5370.317963702526868f, 5432.503314924330880f, 5495.408738576248652f, 5559.042572704037411f, 5623.413251903491073f, 5688.529308438412954f,
     5754.399373371566071f, 5821.032177708715608f, 5888.436553555889986f, 5956.621435290103364f, 6025.595860743575031f, 6095.368972401694009f, 6165.950018614822511f, 6237.348354824191119f,
     6309.573444801930236f, 6382.634861905488833f, 6456.542290346556001f, 6531.305526474722683f, 6606.934480075957254f, 6683.439175686148701f, 6760.829753919818359f, 6839.116472814293047f,
     6918.309709189362366f, 6998.419960022738451f, 7079.457843841380964f, 7161.434102129020175f, 7244.359600749898163f, 7328.245331389044622f, 7413.102413009177326f, 7498.942093324558300f,
     7585.775750291835720f, 7673.614893618185306f, 7762.471166286919470f, 7852.356346100717928f, 7943.282347242813557f, 8035.261221856168959f, 8128.305161640994811f, 8222.426499470711860f,
     8317.637711026709439f, 8413.951416451947807f, 8511.380382023768107f, 8609.937521846008167f, 8709.635899560806138f, 8810.488730080138339f, 8912.509381337458763f, 9015.711376059571194f,
     9120.108393559095930f, 9225.714271547629323f, 9332.543007969914470f, 9440.608762859235867f, 9549.925860214360000f, 9660.508789898130999f, 9772.372209558112445f, 9885.530946569389926f,
     10000.000000000000000f, 10115.794542598981934f, 10232.929922807536059f, 10351.421666793430632f, 10471.285480508984620f, 10592.537251772897434f, 10715.193052376071137f, 10839.269140212038110f,
     10964.781961431850505f, 11091.748152624008071f, 11220.184543019629928f, 11350.108156723143111f, 11481.536214968817148f, 11614.486138403437508f, 11748.975549395303460f, 11885.022274370188825f,
     12022.644346174130987f, 12161.860006463679383f, 12302.687708123810808f, 12445.146117713844433f, 12589.254117941662116f, 12735.030810166628726f, 12882.495516931348902f, 13031.667784523000591f,
     13182.567385564074357f, 13335.214321633240615f, 13489.628825916532151f, 13645.831365889238441f, 13803.842646028837407f, 13963.683610559362933f, 14125.375446227553766f, 14288.939585111036649f,
     14454.397707459280355f, 14621.771744567182395f, 14791.083881682072388f, 14962.356560944328521f, 15135.612484362070973f, 15310.874616820286974f, 15488.166189124827724f, 15667.510701081500883f,
     15848.931924611140857f, 16032.453906900416769f, 16218.100973589298519f, 16405.897731995388312f, 16595.869074375597847f, 16788.040181225587730f, 16982.436524617460236f, 17179.083871575894591f,
     17378.008287493761600f, 17579.236139586930221f, 17782.794100389226514f, 17988.709151287872373f, 18197.008586099826061f, 18407.720014689544769f, 18620.871366628656688f, 18836.490894898019178f,
     19054.607179632483167f, 19275.249131909367861f, 19498.445997580456606f, 19724.227361148532509f, 19952.623149688788544f, 20183.663636815595964f, 20417.379446695274964f, 20653.801558105311415f,
     20892.961308540408936f, 21134.890398366474983f, 21379.620895022322657f, 21627.185237270201469f, 21877.616239495517220f, 22130.947096056363080f, 22387.211385683378467f, 22646.443075930619671f,
     22908.676527677747799f, 23173.946499684796436f, 23442.288153199227963f, 23713.737056616551854f, 23988.329190194897819f, 24266.100950824144093f, 24547.089156850284780f, 24831.331052955676569f,
     25118.864315095823258f, 25409.727055493065563f, 25703.957827688645921f, 26001.595631652722659f, 26302.679918953814195f, 26607.250597988084337f, 26915.348039269138098f, 27227.013080779099255f,
     27542.287033381689980f, 27861.211686297723645f, 28183.829312644549645f, 28510.182675039097376f, 28840.315031266058213f, 29174.270140011660260f, 29512.092266663839837f, 29853.826189179570065f,
     30199.517204020190547f, 30549.211132155152882f, 30902.954325135920953f, 31260.793671239556716f, 31622.776601683792251f, 31988.951096913973743f, 32359.365692962808680f, 32734.069487883793045f,
     33113.112148259075184f, 33496.543915782793192f, 33884.415613920275064f, 34276.778654645044298f, 34673.685045253165299f, 35075.187395256791206f, 35481.338923357528984f, 35892.193464500494883f,
     36307.805477010100731f, 36728.230049808502372f, 37153.522909717277798f, 37583.740428844430426f, 38018.939632056128175f, 38459.178204535353871f, 38904.514499428049021f, 39355.007545577725978f,
     39810.717055349690781f, 40271.703432545946271f, 40738.027780411306594f, 41209.751909733044158f, 41686.938347033552418f, 42169.650342858221848f, 42657.951880159256689f, 43151.907682776502043f,
     43651.583224016561871f, 44157.044735331204720f, 44668.359215096344997f, 45185.594437492261932f, 45708.818961487515480f, 46238.102139926035306f, 46773.514128719813016f, 47315.125896148027095f,
     47863.009232263801096f, 48417.236758409890172f, 48977.881936844663869f, 49545.019080479054537f, 50118.723362727250787f, 50699.070827470444783f, 51286.138399136478256f, 51880.003892896093021f,
     52480.746024977226625f, 53088.444423098793777f, 53703.179637025328702f, 54325.033149243361549f, 54954.087385762482882f, 55590.425727040368656f, 56234.132519034908910f, 56885.293084384131362f,
     57543.993733715666167f, 58210.321777087097871f, 58884.365535558841657f, 59566.214352901093662f, 60255.958607435808517f, 60953.689724016934633f, 61659.500186148223293f, 62373.483548241914832f,
     63095.734448019298725f, 63826.348619054828305f, 64565.422903465492709f, 65313.055264747294132f, 66069.344800759645295f, 66834.391756861485192f, 67608.297539198189043f, 68391.164728142932290f,
     69183.097091893621837f, 69984.199600227308110f, 70794.578438413736876f, 71614.341021290267236f, 72443.596007499058032f, 73282.453313890437130f, 74131.024130091769621f, 74989.420933245579363f,
     75857.757502918364480f, 76736.148936181853060f, 77624.711662869114662f, 78523.563461007099249f, 79432.823472428208333f, 80352.612218561771442f, 81283.051616409953567f, 82224.264994707118603f,
     83176.377110267087119f, 84139.514164519481710f, 85113.803820237590116f, 86099.375218459987082f, 87096.358995608141413f, 88104.887300801463425f, 89125.093813374594902f, 90157.113760595704662f,
     91201.083935590970214f, 92257.142715476278681f, 93325.430079699057387f, 94406.087628592256806f, 95499.258602143687312f, 96605.087898981408216f, 97723.722095581120811f, 98855.309465693906532f,
     100000.000000000000000f, 101157.945425989833893f, 102329.299228075367864f, 103514.216667934306315f, 104712.854805089853471f, 105925.372517728974344f, 107151.930523760704091f, 108392.691402120384737f,
     109647.819614318505046f, 110917.481526240095263f, 112201.845430196292000f, 113501.081567231434747f, 114815.362149688167847f, 116144.861384034375078f, 117489.755493953038240f, 118850.222743701888248f,
     120226.443461741306237f, 121618.600064636790194f, 123026.877081238111714f, 124451.461177138437051f, 125892.541179416613886f, 127350.308101666290895f, 128824.955169313485385f, 130316.677845230005914f,
     131825.673855640750844f, 133352.143216332391603f, 134896.288259165332420f, 136458.313658892380772f, 138038.426460288377712f, 139636.836105593625689f, 141253.754462275537662f, 142889.395851110370131f,
     144543.977074592810823f, 146217.717445671820315f, 147910.838816820731154f, 149623.565609443292487f, 151356.124843620724278f, 153108.746168202866102f, 154881.661891248280881f, 156675.107010815001559f,
     158489.319246111408575f, 160324.539069004182238f, 162181.009735892992467f, 164058.977319953875849f, 165958.690743755956646f, 167880.401812255906407f, 169824.365246174595086f, 171790.838715758931357f,
     173780.082874937623274f, 175792.361395869316766f, 177827.941003892279696f, 179887.091512878745561f, 181970.085860998253338f, 184077.200146895454964f, 186208.713666286552325f, 188364.908948980213609f,
     190546.071796324831666f, 192752.491319093649508f, 194984.459975804551505f, 197242.273611485346919f, 199526.231496887892717f, 201836.636368155974196f, 204173.794466952735092f, 206538.015581053128699f,
     208929.613085404096637f, 211348.903983664757106f, 213796.208950223255670f, 216271.852372702007415f, 218776.162394955172203f, 221309.470960563630797f, 223872.113856833777390f, 226464.430759306182154f,
     229086.765276777470717f, 231739.464996847964358f, 234422.881531992257806f, 237137.370566165540367f, 239883.291901948978193f, 242661.009508241433650f, 245470.891568502847804f, 248313.310529556794791f,
     251188.643150958203478f, 254097.270554930641083f, 257039.578276886459207f, 260015.956316527212039f, 263026.799189538171049f, 266072.505979880865198f, 269153.480392691388261f, 272270.130807791021653f,
     275422.870333816914354f, 278612.116862977214623f, 281838.293126445496455f, 285101.826750390988309f, 288403.150312660553027f, 291742.701400116609875f, 295120.922666638391092f, 298538.261891795729753f,
     301995.172040201898199f, 305492.111321551550645f, 309029.543251359194983f, 312607.936712395574432f, 316227.766016837907955f, 319889.510969139693771f, 323593.656929628108628f, 327340.694878837966826f,
     331131.121482590795495f, 334965.439157827931922f, 338844.156139202765189f, 342767.786546450457536f, 346736.850452531652991f, 350751.873952567926608f, 354813.389233575318940f, 358921.934645004977938f,
     363078.054770101036411f, 367282.300498084980063f, 371535.229097172792535f, 375837.404288444318809f, 380189.396320561238099f, 384591.782045353553258f, 389045.144994280475657f, 393550.075455777230673f,
     398107.170553496922366f, 402717.034325459506363f, 407380.277804113051388f, 412097.519097330397926f, 416869.383470335509628f, 421696.503428582276683f, 426579.518801592581440f, 431519.076827765034977f,
     436515.832240165618714f, 441570.447353312047198f, 446683.592150963493623f, 451855.944374922604766f, 457088.189614875183906f, 462381.021399260323960f, 467735.141287198115606f, 473151.258961480285507f,
     478630.092322638025507f, 484172.367584098887164f, 489778.819368446653243f, 495450.190804790530819f, 501187.233627272478770f, 506990.708274704462383f, 512861.383991364797112f, 518800.038928960973863f,
     524807.460249772295356f, 530884.444230987923220f, 537031.796370253316127f, 543250.331492433557287f, 549540.873857624828815f, 555904.257270403672010f, 562341.325190349132754f, 568852.930843841284513f,
     575439.937337156618014f, 582103.217770871007815f, 588843.655355588416569f, 595662.143529010936618f, 602559.586074358085170f, 609536.897240169346333f, 616595.001861482160166f, 623734.835482419119217f,
     630957.344480192987248f, 638263.486190548283048f, 645654.229034654912539f, 653130.552647472941317f, 660693.448007596423849f, 668343.917568614939228f, 676082.975391981890425f, 683911.647281429264694f,
     691830.970918936305679f, 699841.996002273052000f, 707945.784384137368761f, 716143.410212902701460f, 724435.960074990638532f, 732824.533138904487714f, 741310.241300917696208f, 749894.209332455880940f,
     758577.575029183528386f, 767361.489361818530597f, 776247.116628691204824f, 785235.634610070963390f, 794328.234724282170646f, 803526.122185617685318f, 812830.516164099564776f, 822242.649947071215138f,
     831763.771102670812979f, 841395.141645194729790f, 851138.038202375872061f, 860993.752184599870816f, 870963.589956081472337f, 881048.873008014634252f, 891250.938133745919913f, 901571.137605956988409f,
     912010.839355909614824f, 922571.427154762903228f, 933254.300796990515664f, 944060.876285922597162f, 954992.586021436844021f, 966050.878989814082161f, 977237.220955811208114f, 988553.094656939036213f,
     1000000.000000000000000f
  };

  constexpr float db2a(float db)
  {
    if (db < 0)
      return 1.0f / db2a(-db);

    if (db < 120.0f)
    {
      auto db_10 = db * 10.0f;
      auto index = std::uint32_t(db_10);
      auto v1 = inv_db_table[index];
      auto v2 = inv_db_table[index + 1];
      return linear_interpolate(v1, v2, db_10 - index);
    }
    return 1000000.0f;
  }
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // decibel is a highly optimized class for dealing with decibels. The
  // class provides fast conversion from linear to decibel and back. The
  // decibel value is non-linear. It operates on the logarithmic domain. The
  // decibel class is perfectly suitable for dynamics processing (e.g.
  // compressors and limiters and envelopes). Q provides fast decibel
  // computations using lookup tables and fast math computations for
  // converting to and from scalars.
  ////////////////////////////////////////////////////////////////////////////
  struct decibel_unit;

  //struct decibel : unit<double, decibel>
  //{
  //  using base_type = unit<double, decibel>;
  //  using base_type::base_type;
  //  using unit_type = decibel_unit;

  //  // Please use lin_to_db(s) instead of decibel(s).
  //  //
  //  decibel(double) = delete;
  //  //
  //  // This constructor was used in previous versions of the library to
  //  // convert linear to decibels, which can be confusing. This version
  //  // ought to have corrected this nonintuitive semantics, but this is a
  //  // disruptive change that will alter the semantics of all existing code
  //  // without warning. In order to avoid further confusion, we will mark
  //  // this constructor as deleted for the time being, making it a hard
  //  // error to alert users upgrading to this library version.
  //};

  // Free functions
  double               lin_double(double db);
  constexpr float      lin_float(double db);
  inline double       approx_db(float val);
  double              lin_to_db(double val);

  ////////////////////////////////////////////////////////////////////////////
  // Inlines
  ////////////////////////////////////////////////////////////////////////////
  inline double lin_double(double db)
  {
    return std::pow(10, db / 20);
  }

  constexpr float lin_float(double db)
  {
    return detail::db2a(db);
  }

  inline double approx_db(float val)
  {
    return  20.0f * faster_log10(val);
  }

  inline double lin_to_db(double val)
  {
    return 20.0f * fast_log10(val);
  }
}



namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  inline namespace literals
  {
    constexpr frequency operator ""_Hz(long double val)
    {
      return frequency{ double(val) };
    }

    constexpr frequency operator ""_Hz(unsigned long long int val)
    {
      return frequency{ double(val) };
    }

    constexpr frequency operator ""_KHz(long double val)
    {
      return frequency{ double(val * 1e3) };
    }

    constexpr frequency operator ""_KHz(unsigned long long int val)
    {
      return frequency{ double(val * 1e3) };
    }

    constexpr frequency operator ""_kHz(long double val)
    {
      return frequency{ double(val * 1e3) };
    }

    constexpr frequency operator ""_kHz(unsigned long long int val)
    {
      return frequency{ double(val * 1e3) };
    }

    constexpr frequency operator ""_MHz(long double val)
    {
      return frequency{ double(val * 1e6) };
    }

    constexpr frequency operator ""_MHz(unsigned long long int val)
    {
      return frequency{ double(val * 1e6) };
    }

    constexpr duration operator ""_s(long double val)
    {
      return duration{ double(val) };
    }

    constexpr duration operator ""_s(unsigned long long int val)
    {
      return duration{ double(val) };
    }

    constexpr duration operator ""_ms_(long double val)
    {
      return duration{ double(val * 1e-3) };
    }

    constexpr duration operator ""_ms_(unsigned long long int val)
    {
      return duration{ double(val * 1e-3) };
    }

    constexpr duration operator ""_us(long double val)
    {
      return duration{ double(val * 1e-6) };
    }

    constexpr duration operator ""_us(unsigned long long int val)
    {
      return duration{ double(val * 1e-6) };
    }

    constexpr double operator ""_dB(unsigned long long int val)
    {
      return val;
    }

    constexpr double operator ""_dB(long double val)
    {
      return val;
    }

    constexpr long double operator ""_pi(long double val)
    {
      return val * pi;
    }

    constexpr long double operator ""_pi(unsigned long long int val)
    {
      return val * pi;
    }
  }
}


namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // 3-point 1D median filter. Returns the median of 3 latest samples. The
  // median filter is a nonlinear digital filter often used to remove noise
  // from signals. The median filter performs better than a moving average
  // filter in reducing noise, especially outliers (sudden jumps in the
  // signal (impulsive noise)). The median filter, however, requires more
  // processing cycles.
  ////////////////////////////////////////////////////////////////////////////
  inline float median3f(float a, float b, float c)
  {
    return std::max(std::min(a, b), std::min(std::max(a, b), c));
  }

  struct median3
  {
    median3(float median_ = 0.0f)
      : _median(median_)
      , b(median_)
      , c(median_)
    {}

    float operator()(float a)
    {
      _median = median3f(a, b, c);
      c = b;
      b = a;
      return _median;
    }

    float operator()() const
    {
      return _median;
    }

    median3& operator=(float median_)
    {
      b = c = _median = median_;
      return *this;
    }

    float _median = 0.0f;
    float b = 0.0f;
    float c = 0.0f;
  };
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // The bitset class stores bits efficiently using integers <T>. Data is
  // stored in a std::vector with a size that is fixed at construction time,
  // given the number of bits required.
  //
  // Member functions are provided for:
  //
  //    1. Setting individual bits and ranges of bits
  //    2. Geting each bit at position i
  //    3. Clearing all bits
  //    4. Getting the actual integers that stores the bits.
  ////////////////////////////////////////////////////////////////////////////
  template <typename T = natural_uint>
  class bitset
  {
  public:

    using value_type = T;
    using vector_type = std::vector<T>;

    static_assert(std::is_unsigned<T>::value, "T must be unsigned");
    static constexpr auto value_size = CHAR_BIT * sizeof(T);
    static constexpr auto one = T{ 1 };

    bitset(std::size_t num_bits);
    bitset(bitset const& rhs) = default;
    bitset(bitset&& rhs) = default;

    bitset& operator=(bitset const& rhs) = default;
    bitset& operator=(bitset&& rhs) = default;

    std::size_t    size() const;
    void           clear();
    void           set(std::size_t i, bool val);
    void           set(std::size_t i, std::size_t n, bool val);
    bool           get(std::size_t i) const;

    T* data();
    T const* data() const;

  private:

    vector_type    _bits;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Implementation
  ////////////////////////////////////////////////////////////////////////////
  template <typename T>
  inline bitset<T>::bitset(std::size_t num_bits)
  {
    auto array_size = (num_bits + value_size - 1) / value_size;
    _bits.resize(array_size, 0);
  }

  template <typename T>
  inline std::size_t bitset<T>::size() const
  {
    return _bits.size() * value_size;
  }

  template <typename T>
  inline void bitset<T>::clear()
  {
    std::fill(_bits.begin(), _bits.end(), 0);
  }

  template <typename T>
  inline void bitset<T>::set(std::size_t i, bool val)
  {
    // Check that we don't get past the storage
    if (i > size())
      return;

    auto mask = 1 << (i % value_size);
    auto& ref = _bits[i / value_size];
    ref ^= (-T(val) ^ ref) & mask;
  }

  template <typename T>
  inline bool bitset<T>::get(std::size_t i) const
  {
    // Check we don't get past the storage
    if (i > size())
      return 0;

    auto mask = one << (i % value_size);
    return (_bits[i / value_size] & mask) != 0;
  }

  template <typename T>
  inline void bitset<T>::set(std::size_t i, std::size_t n, bool val)
  {
    // Check that the index (i) does not get past size
    auto size_ = size();
    if (i > size_)
      return;

    // Check that the n does not get past the size
    if ((i + n) > size_)
      n = size_ - i;

    constexpr auto all_ones = int_traits<T>::max;
    auto* p = _bits.data();
    p += i / value_size;    // Adjust the buffer pointer for the current index (i)

    // Do the first partial int
    auto mod = i & (value_size - 1);
    if (mod)
    {
      // mask off the high n bits we want to set
      mod = value_size - mod;

      // Calculate the mask
      T mask = ~(all_ones >> mod);

      // Adjust the mask if we're not going to reach the end of this int
      if (n < mod)
        mask &= (all_ones >> (mod - n));

      if (val)
        *p |= mask;
      else
        *p &= ~mask;

      // Fast exit if we're done here!
      if (n < mod)
        return;

      n -= mod;
      ++p;
    }

    // Write full ints while we can - effectively doing value_size bits at a time
    if (n >= value_size)
    {
      // Store a local value to work with
      T val_ = val ? all_ones : 0;

      do
      {
        *p++ = val_;
        n -= value_size;
      } while (n >= value_size);
    }

    // Now do the final partial int, if necessary
    if (n)
    {
      mod = n & (value_size - 1);

      // Calculate the mask
      T mask = (one << mod) - 1;

      if (val)
        *p |= mask;
      else
        *p &= ~mask;
    }
  }

  template <typename T>
  inline T* bitset<T>::data()
  {
    return _bits.data();
  }

  template <typename T>
  inline T const* bitset<T>::data() const
  {
    return _bits.data();
  }
}

namespace cycfi::q::detail
{
  template <typename C, typename = int>
  struct resizable_container
    : std::false_type
  {};

  template <typename C>
  struct resizable_container<C, decltype(std::declval<C>().resize(1), 0)>
    : std::true_type
  {};

  template <typename T>
  void init_store(std::size_t size, std::vector<T>& _data, std::size_t& _mask)
  {
    // allocate the data in a size that is a power of two for efficient indexing
    std::size_t capacity = smallest_pow2(size);
    _mask = capacity - 1;
    _data.resize(capacity, T{});
  }

  template <typename T, std::size_t N>
  void init_store(std::array<T, N>& _data, std::size_t& _mask)
  {
    static_assert(is_pow2(N),
      "Error: Storage must have a size that is a power of two");
    _mask = _data.size() - 1;
  }
}

namespace cycfi::q
{
  namespace sample_interpolation
  {
    struct none
    {
      template <typename Storage, typename T>
      T operator()(Storage const& buffer, T index) const
      {
        return buffer[std::size_t(index)];
      }
    };

    struct linear
    {
      template <typename Storage, typename T>
      T operator()(Storage const& buffer, T index) const
      {
        auto y1 = buffer[std::size_t(index)];
        auto y2 = buffer[std::size_t(index) + 1];
        auto mu = index - std::floor(index);
        return linear_interpolate(y1, y2, mu);
      }
    };
  }
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // ring_buffer
  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Storage = std::vector<T>>
  class ring_buffer
  {
  public:

    using value_type = T;
    using storage_type = Storage;
    using index_type = std::size_t;
    using interpolation_type = sample_interpolation::none;

    explicit ring_buffer();
    explicit ring_buffer(std::size_t size);
    ring_buffer(ring_buffer const& rhs) = default;
    ring_buffer(ring_buffer&& rhs) = default;

    ring_buffer& operator=(ring_buffer const& rhs) = default;
    ring_buffer& operator=(ring_buffer&& rhs) = default;

    std::size_t       size() const;
    void              push(T val);
    T const& front() const;
    T& front();
    T const& back() const;
    T& back();
    T const& operator[](std::size_t index) const;
    T& operator[](std::size_t index);
    void              clear();
    void              pop_front();

    Storage& store();
    const Storage& store() const;

  private:

    std::size_t	      _mask;
    std::size_t       _pos;
    Storage           _data;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Implementation
  ////////////////////////////////////////////////////////////////////////////
  template <typename T, typename Storage>
  inline ring_buffer<T, Storage>::ring_buffer()
    : _pos(0)
  {
    static_assert(!detail::resizable_container<Storage>::value,
      "Error: Not default constructible for resizable buffers");
    detail::init_store(_data, _mask);
  }

  template <typename T, typename Storage>
  inline ring_buffer<T, Storage>::ring_buffer(std::size_t size)
    : _pos(0)
  {
    static_assert(detail::resizable_container<Storage>::value,
      "Error: Can't be constructed with size. Storage has fixed size.");
    detail::init_store(size, _data, _mask);
  }

  // Get the size of buffer.
  template <typename T, typename Storage>
  inline std::size_t ring_buffer<T, Storage>::size() const
  {
    return _data.size();
  }

  // Push the latest element, overwriting the oldest element.
  template <typename T, typename Storage>
  inline void ring_buffer<T, Storage>::push(T val)
  {
    --_pos &= _mask;
    _data[_pos] = val;
  }

  // Get the latest element.
  template <typename T, typename Storage>
  inline T const& ring_buffer<T, Storage>::front() const
  {
    return (*this)[0];
  }

  // Get the latest element.
  template <typename T, typename Storage>
  inline T& ring_buffer<T, Storage>::front()
  {
    return (*this)[0];
  }

  // Get the oldest element.
  template <typename T, typename Storage>
  inline T const& ring_buffer<T, Storage>::back() const
  {
    return (*this)[size() - 1];
  }

  // Get the oldest element.
  template <typename T, typename Storage>
  inline T& ring_buffer<T, Storage>::back()
  {
    return (*this)[size() - 1];
  }

  // Get the nth latest element (b[0] is latest element, b[1] is the second
  // and b[size()-1] is the oldest.
  template <typename T, typename Storage>
  inline T const& ring_buffer<T, Storage>::operator[](std::size_t index) const
  {
    return _data[(_pos + index) & _mask];
  }

  // Get the nth latest element (b[0] is latest element, b[1] is the second
  // and b[size()-1] is the oldest.
  template <typename T, typename Storage>
  inline T& ring_buffer<T, Storage>::operator[](std::size_t index)
  {
    return _data[(_pos + index) & _mask];
  }

  // Clear the ring_buffer
  template <typename T, typename Storage>
  inline void ring_buffer<T, Storage>::clear()
  {
    for (auto& e : _data)
      e = T();
  }

  // Remove the front element
  template <typename T, typename Storage>
  inline void ring_buffer<T, Storage>::pop_front()
  {
    ++_pos;
  }

  // Raw access to the data storage
  template <typename T, typename Storage>
  inline Storage& ring_buffer<T, Storage>::store()
  {
    return _data;
  }

  // Raw access to the data storage
  template <typename T, typename Storage>
  inline const Storage& ring_buffer<T, Storage>::store() const
  {
    return _data;
  }
}



namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // The zero_crossing class saves zero-crossing information necessary to
  // extract accurate timing information such as periods between pulses for
  // performing analysis such as bitstream autocorrelation.
  //
  // Each zero crossing pulse is saved in a ring buffer of info elements.
  // Data include the maximum height of the waveform bounded by the pulse,
  // the pulse width, as well as the leading edge and trailing edge frame
  // positions (number of samples from the start) and y coordinates (the
  // sample values before and after each zero crossing) of the zero
  // crossings.
  //
  // Only the latest few (finite amount) of zero crossing information is
  // saved, given by the window constructor parameter. The window is the
  // number of frames (samples) of information held by the zero_crossing
  // data structure.
  //
  // Each call to the function operator, given a sample s, returns the
  // zero-crossing state (bool). is_ready() returns true when we have
  // sufficient info to perform analysis. is_ready() returns true after
  // every window/2 frames. Information about each zero crossing can be
  // obtained using the index operator[]. The leftmost edge (oldest) is at
  // the 0th index while the rightmost edge (latest) is at index
  // num_edges()-1.
  //
  // After window/2 frames, the leading edge and trailing edge frame
  // positions are shifted by -window/2 such that an edge at frame index N
  // will be shifted to N-window/2. For example, if the window size is 100
  // and the leading edge is at frame 45, it will be shifted to -5 (45-50).
  //
  // This procedure is done to ensure seamless operation from one window to
  // the next. In the example above, frame index -5 is already past the left
  // side of the window, but will still be kept as long as the trailing edge
  // is still within the window. Take note that it is also possible for the
  // latest edge to have a trailing edge that goes past the right side of
  // the window. If for example, with the same window size 100, there can be
  // an edge with a leading edge at 95 and trailing edge at 120.
  ////////////////////////////////////////////////////////////////////////////
  class zero_crossing_collector
  {
  public:

    static constexpr float pulse_height_diff = 0.8f;
    static constexpr float pulse_width_diff = 0.85f;
    static constexpr auto undefined_edge = int_min<int>();

    struct info
    {
      using crossing_data = std::pair<float, float>;

      void              update_peak(float s, std::size_t frame);
      std::size_t       period(info const& next) const;
      float             fractional_period(info const& next) const;
      int               width() const;
      bool              similar(info const& next) const;

      crossing_data     _crossing;
      float             _peak;
      int               _leading_edge = undefined_edge;
      int               _trailing_edge = undefined_edge;
      float             _width = 0.0f;
    };


    zero_crossing_collector& operator=(const zero_crossing_collector& rhs) {
      if (this != &rhs) {
        // Perform deep copy of data members
        _prev = rhs._prev;
        _hysteresis = rhs._hysteresis;
        _state = rhs._state;
        _num_edges = rhs._num_edges;
        _window_size = rhs._window_size;
        _info = rhs._info;
        _frame = rhs._frame;
        _ready = rhs._ready;
        _peak_update = rhs._peak_update;
        _peak = rhs._peak;
      }
      return *this;
    }

    zero_crossing_collector(decibel hysteresis, duration window, float sps);
    zero_crossing_collector(decibel hysteresis, std::uint32_t window);
    zero_crossing_collector(zero_crossing_collector const& rhs) = default;
    zero_crossing_collector(zero_crossing_collector&& rhs) = default;

    std::size_t          num_edges() const;
    std::size_t          capacity() const;
    std::size_t          frame() const;
    std::size_t          window_size() const;
    bool                 is_ready() const;
    float                peak_pulse() const;
    bool                 is_reset() const;

    bool                 operator()(float s);
    bool                 operator()() const;
    info const& operator[](std::size_t index) const;
    info& operator[](std::size_t index);

  private:

    void                 update_state(float s);
    void                 shift(std::size_t n);
    void                 reset();

    using info_storage = ring_buffer<info>;

    float                _prev = 0.0f;
    float          _hysteresis;
    bool                 _state = false;
    std::size_t          _num_edges = 0;
    std::size_t    _window_size;
    info_storage         _info;
    std::size_t          _frame = 0;
    bool                 _ready = false;
    float                _peak_update = 0.0f;
    float                _peak = 0.0f;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Implementation
  ////////////////////////////////////////////////////////////////////////////
  namespace detail
  {
    inline std::size_t adjust_window_size(std::size_t window)
    {
      constexpr auto bits = bitset<>::value_size;
      return std::max<std::size_t>(2, (window + bits - 1) / bits);
    }
  }

  inline zero_crossing_collector::zero_crossing_collector(decibel hysteresis, duration window, float sps)
    : zero_crossing_collector{ hysteresis, std::uint32_t(as_double(window) * sps) }
  {
  }

  inline zero_crossing_collector::zero_crossing_collector(decibel hysteresis, std::uint32_t window)
    : _hysteresis(-lin_float(hysteresis))
    , _window_size(detail::adjust_window_size(window)* bitset<>::value_size)
    , _info(_window_size / 2)
  {}

  inline void zero_crossing_collector::info::update_peak(float s, std::size_t frame)
  {
    _peak = std::max(s, _peak);
    if ((_width == 0.0f) && (s < (_peak * 0.3)))
      _width = frame - _leading_edge;
  }

  inline std::size_t zero_crossing_collector::info::period(info const& next) const
  {
    CYCFI_ASSERT(_leading_edge <= next._leading_edge, "Invalid order.");
    return next._leading_edge - _leading_edge;
  }

  inline bool zero_crossing_collector::info::similar(info const& next) const
  {
    return rel_within(_peak, next._peak, 1.0f - pulse_height_diff) &&
      rel_within(_width, next._width, 1.0f - pulse_width_diff);
  }

  inline float zero_crossing_collector::info::fractional_period(info const& next) const
  {
    CYCFI_ASSERT(_leading_edge <= next._leading_edge, "Invalid order.");

    // Get the start edge
    auto prev1 = _crossing.first;
    auto curr1 = _crossing.second;
    auto dy1 = curr1 - prev1;
    auto dx1 = -prev1 / dy1;

    // Get the next edge
    auto prev2 = next._crossing.first;
    auto curr2 = next._crossing.second;
    auto dy2 = curr2 - prev2;
    auto dx2 = -prev2 / dy2;

    // Calculate the fractional period
    auto result = next._leading_edge - _leading_edge;
    return result + (dx2 - dx1);
  }

  inline std::size_t zero_crossing_collector::num_edges() const
  {
    return _num_edges;
  }

  inline std::size_t zero_crossing_collector::capacity() const
  {
    return _info.size();
  }

  inline std::size_t zero_crossing_collector::frame() const
  {
    return _frame;
  }

  inline std::size_t zero_crossing_collector::window_size() const
  {
    return _window_size;
  }

  inline void zero_crossing_collector::reset()
  {
    _num_edges = 0;
    _state = false;
    _frame = 0;
  }

  inline bool zero_crossing_collector::is_reset() const
  {
    return _frame == 0;
  }

  inline bool zero_crossing_collector::is_ready() const
  {
    return _ready;
  }

  inline float zero_crossing_collector::peak_pulse() const
  {
    return std::max(_peak, _peak_update);
  }

  inline void zero_crossing_collector::update_state(float s)
  {
    if (_ready)
    {
      shift(_window_size / 2);
      _ready = false;
      _peak = _peak_update;
      _peak_update = 0.0f;
    }

    if (num_edges() >= capacity())
      reset();

    if (s > 0.0f)
    {
      if (!_state)
      {
        CYCFI_ASSERT(_num_edges < _info.size(), "Bad _size");
        _info.push({ { _prev, s }, s, int(_frame) });
        ++_num_edges;
        _state = 1;
      }
      else
      {
        _info[0].update_peak(s, _frame);
      }
      if (s > _peak_update)
      {
        _peak_update = s;
      }
    }
    else if (_state && s < _hysteresis)
    {
      _state = 0;
      auto& info = _info[0];
      info._trailing_edge = _frame;
      if (_peak == 0.0f)
        _peak = _peak_update;
    }

    _prev = s;
  }

  inline bool zero_crossing_collector::operator()(float s)
  {
    // Offset s by half of hysteresis, so that zero cross detection is
    // centered on the actual zero.
    s += _hysteresis / 2;

    if (num_edges() >= capacity())
      reset();

    if ((_frame == _window_size / 2) && num_edges() == 0)
      reset();

    update_state(s);

    if (++_frame >= _window_size && !_state)
    {
      // Remove half the size from _frame, so we can continue seamlessly
      _frame -= _window_size / 2;

      // We need at least two rising edges.
      if (num_edges() > 1)
        _ready = true;
      else
        reset();
    }

    return _state;
  };

  inline bool zero_crossing_collector::operator()() const
  {
    return _state;
  }

  inline zero_crossing_collector::info const&
    zero_crossing_collector::operator[](std::size_t index) const
  {
    return _info[(_num_edges - 1) - index];
  }

  inline zero_crossing_collector::info&
    zero_crossing_collector::operator[](std::size_t index)
  {
    return _info[(_num_edges - 1) - index];
  }

  inline void zero_crossing_collector::shift(std::size_t n)
  {
    _info[0]._leading_edge -= n;
    if (!_state)
      _info[0]._trailing_edge -= n;
    std::size_t i = 1;
    for (; i != _num_edges; ++i)
    {
      _info[i]._leading_edge -= n;
      int edge = (_info[i]._trailing_edge -= n);
      if (edge < 0.0f)
        break;
    }
    _num_edges = i;
  }
}

#ifdef _MSC_VER
# include <intrin.h>
# include <nmmintrin.h>
#endif

namespace cycfi::q::detail
{
  inline std::uint32_t count_bits(std::uint32_t i)
  {
#if defined(_MSC_VER)
    return __popcnt(i);
#elif defined(__GNUC__)
    return __builtin_popcount(i);
#else
# error Unsupported compiler
#endif
  }

#if (!defined(_MSC_VER) || defined(_WIN64))
  inline std::uint64_t count_bits(std::uint64_t i)
  {
#if defined(_MSC_VER)
    return _mm_popcnt_u64(i);
#elif defined(__GNUC__)
    return __builtin_popcountll(i);
#else
# error Unsupported compiler
#endif
  }
#endif // (!defined(_MSC_VER) || defined(_WIN64))
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // The bitstream_acf correlates class a bit stream (stored in a bitset) by
  // itself shifted by position, pos.
  //
  // In standard ACF (autocorrelation function) the signal is multiplied by
  // a shifted version of itself, delayed by position, pos, with the result
  // accumulated for each point, half the window of interest. The higher the
  // sum, the higher the periodicity.
  //
  // With bitstream auto correlation, the more efficient XOR operation is
  // used instead. A single XOR operation works on N bits of an integer. We
  // get a speedup factor of N (e.g. 64 for a 64 bit machine) compared to
  // standard ACF, and not to mention that integer bit operations are a lot
  // faster than floating point multiplications).
  //
  // With XOR you get a one when there’s a mismatch:
  //
  //    0 ^ 0 = 0
  //    0 ^ 1 = 1
  //    1 ^ 0 = 1
  //    1 ^ 1 = 0
  //
  // After XOR, the number of bits (set to 1) is counted. The lower the
  // count, the higher the periodicity. A count of zero gives perfect
  // correlation: there is no mismatch.
  ////////////////////////////////////////////////////////////////////////////
  template <typename T = natural_uint>
  struct bitstream_acf
  {
    static constexpr auto value_size = bitset<>::value_size;

    bitstream_acf(bitset<T> const& bits)
      : _bits(bits)
      , _mid_array(std::max<std::size_t>(((bits.size() / value_size) / 2) - 1, 1))
    {}

    std::size_t operator()(std::size_t pos) const
    {
      auto const index = pos / value_size;
      auto const shift = pos % value_size;

      auto const* p1 = _bits.data();
      auto const* p2 = _bits.data() + index;
      auto count = 0;

      if (shift == 0)
      {
        for (std::size_t i = 0; i != _mid_array; ++i)
          count += detail::count_bits(*p1++ ^ *p2++);
      }
      else
      {
        auto shift2 = value_size - shift;
        for (std::size_t i = 0; i != _mid_array; ++i)
        {
          auto v = *p2++ >> shift;
          v |= *p2 << shift2;
          count += detail::count_bits(*p1++ ^ v);
        }
      }
      return count;
    };

    bitset<T> const& _bits;
    std::size_t const    _mid_array;
  };
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // basic_moving_sum computes the moving sum of consecutive samples in a
  // window specified by max_size samples or duration d and float sps.
  //
  // basic_moving_sum can be resized as long as the new size does not exceed
  // the original size (at construction time). When resizing with
  // update=true, when downsizing, the oldest elements are subtracted from
  // the sum. When upsizing, the older elements are added to the sum,
  // otherwise, if update=false, the contents are cleared.
  ////////////////////////////////////////////////////////////////////////////
  template <typename T>
  struct basic_moving_sum
  {
    using value_type = T;

    basic_moving_sum(std::size_t max_size)
      : _buff(max_size)
      , _size(max_size)
      , _sum{ 0 }
    {
      _buff.clear();
    }

    basic_moving_sum(duration d, float sps)
      : basic_moving_sum(std::size_t(sps* as_float(d)))
    {}

    T operator()(value_type s)
    {
      _sum += s;              // Add the latest sample to the sum
      _sum -= _buff[_size - 1]; // Subtract the oldest sample from the sum
      _buff.push(s);          // Push the latest sample, erasing the oldest
      return _sum;
    }

    value_type operator()() const
    {
      return _sum;            // Return the sum
    }

    value_type sum() const
    {
      return _sum;            // Return the sum
    }

    std::size_t size() const
    {
      return _size;
    }

    void resize(std::size_t size, bool update = false)
    {
      // We cannot exceed the original size
      auto new_size = std::min(size, _buff.size());

      if (update)
      {
        if (new_size > _size) // expand
        {
          for (auto i = _size; i != new_size; ++i)
            _sum += _buff[i];
        }
        else // contract
        {
          for (auto i = new_size; i != _size; ++i)
            _sum -= _buff[i];
        }
      }
      else
      {
        clear();
      }
      _size = new_size;
    }

    void resize(duration d, float sps, bool update = false)
    {
      resize(std::size_t(sps * as_float(d)), update);
    }

    void clear()
    {
      _buff.clear();
      _sum = 0;
    }

    void fill(T val)
    {
      _buff.fill(val);
      _sum = val * _size;
    }

  private:

    using buffer = ring_buffer<T>;
    using accumulator = decltype(promote(T()));

    buffer      _buff = buffer{};
    std::size_t _size;
    accumulator _sum;
  };

  using moving_sum = basic_moving_sum<float>;
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // The moving_average is the simplest and most efficient FIR filter. It is
  // also the most common filter in DSP primarily due to its simplicity. But
  // while it is technically a low pass FIR filter, it performs poorly in
  // the frequency domain with very slow roll-off and dreadful stopband
  // attenuation. On the other hand, it performs admirably in the time
  // domain. The moving average filter is optimal in reducing random noise
  // while retaining a sharp step response.
  //
  // Averaging N samples (the moving average window size) increases the SNR
  // by the square root of N. For example, N=16 improves SNR by 4 (12dB).
  // The filter delay is exactly (N−1)/2.
  //
  // The data type, T, is a template parameter, allowing both floating point
  // as well as integer computations. Integers are typically faster than
  // floating point and are not prone to round-off errors.
  //
  // moving_average is a subclass of the moving_sum.
  ////////////////////////////////////////////////////////////////////////////
  template <typename T>
  struct basic_moving_average : basic_moving_sum<T>
  {
    using basic_moving_sum<T>::basic_moving_sum;
    using value_type = T;

    T operator()(T s)
    {
      basic_moving_sum<T>::operator()(s);
      return (*this)();
    }

    T operator()() const
    {
      // Return the average
      return this->sum() / this->size();
    }
  };

  using moving_average = basic_moving_average<float>;

  ////////////////////////////////////////////////////////////////////////////
  // Exponential moving average approximates an arithmetic moving average by
  // multiplying the last result by some factor, and adding it to the next
  // sample multiplied by some other factor.
  //
  // If b = 2/(n+1), where n is the number of samples you would have used in
  // an arithmetic average, the exponential moving average will approximate
  // the arithmetic average pretty well.
  //
  //    n: the number of samples.
  //    y: current value.
  //
  // See: https://www.dsprelated.com/showthread/comp.dsp/47981-1.php
  //
  // The exp_moving_average<n> template computes b at compile-time, where n
  // is supplied as a compile-time parameter.
  ////////////////////////////////////////////////////////////////////////////
  template <std::size_t n>
  struct exp_moving_average
  {
    static constexpr float b = 2.0f / (n + 1);
    static constexpr float b_ = 1.0f - b;

    exp_moving_average(float y_ = 0.0f)
      : y(y_)
    {}

    float operator()(float s)
    {
      return y = b * s + b_ * y;
    }

    float operator()() const
    {
      return y;
    }

    exp_moving_average& operator=(float y_)
    {
      y = y_;
      return *this;
    }

    float y = 0.0f;
  };

  ////////////////////////////////////////////////////////////////////////////
  // The rt_exp_moving_average class computes b at run time, where n is
  // supplied as a runtime parameter.
  ////////////////////////////////////////////////////////////////////////////
  struct rt_exp_moving_average
  {
    rt_exp_moving_average(float n, float y_ = 0.0f)
      : y(y_)
      , b(2.0f / (n + 1))
      , b_(1.0f - b)
    {}

    rt_exp_moving_average(duration d, float sps, float y_ = 0.0f)
      : rt_exp_moving_average(std::size_t(sps* as_float(d)), y_)
    {}

    float operator()(float s)
    {
      return y = b * s + b_ * y;
    }

    float operator()() const
    {
      return y;
    }

    rt_exp_moving_average& operator=(float y_)
    {
      y = y_;
      return *this;
    }

    void width(float n)
    {
      b = 2.0f / (n + 1);
    }

    float y = 0.0f;
    float b, b_;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Simple 2-point average filter
  ////////////////////////////////////////////////////////////////////////////
  struct moving_average2
  {
    float operator()(float s)
    {
      y = (s + _prev) / 2.0f;
      _prev = s;
      return y;
    }

    float operator()() const
    {
      return y;
    }

    float y = 0.0f;
    float _prev = 0.0f;
  };
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // fixed_pt_leaky_integrator: If you want a fast filter for integers, use
  // a fixed point leaky-integrator. k will determine the effect of the
  // filter. Choose k to be a power of 2 for efficiency (the compiler will
  // optimize the computation using shifts). k = 16 is a good starting
  // point.
  //
  // This simulates the RC filter in digital form. The equation is:
  //
  //    y[i] = rho * y[i-1] + s
  //
  // where rho < 1. To avoid floating point, we use k instead, which
  // allows for integer operations. In terms of k, rho = 1 - (1 / k).
  // So the actual formula is:
  //
  //    y[i] += s - (y[i-1] / k);
  //
  // k will also be the filter gain, so the final result should be divided
  // by k. If you need to initialize the filter (y member) to a certain
  // state, you will also need to multiply the initial value by k.
  //
  ////////////////////////////////////////////////////////////////////////////
  template <int k, typename T = int>
  struct fixed_pt_leaky_integrator
  {
    typedef T result_type;
    static constexpr int gain = k;

    T operator()(T s)
    {
      y += s - (y / k);
      return y;
    }

    T operator()() const
    {
      return y;
    }

    fixed_pt_leaky_integrator& operator=(float y_)
    {
      y = y_;
      return *this;
    }

    T y = 0;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Leaky Integrator
  ////////////////////////////////////////////////////////////////////////////
  struct leaky_integrator
  {
    leaky_integrator(float a = 0.995)
      : a(a)
    {}

    leaky_integrator(frequency f, float sps)
      : a(1.0f - (2_pi * as_double(f) / sps))
    {}

    float operator()(float s)
    {
      return y = s + a * (y - s);
    }

    float operator()() const
    {
      return y;
    }

    leaky_integrator& operator=(float y_)
    {
      y = y_;
      return *this;
    }

    void cutoff(frequency f, float sps)
    {
      a = 1.0f - (2_pi * as_double(f) / sps);
    }

    float y = 0.0f, a;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Basic one pole low-pass filter (6dB/Oct)
  ////////////////////////////////////////////////////////////////////////////
  struct one_pole_lowpass
  {
    one_pole_lowpass(float a)
      : a(a)
    {}

    one_pole_lowpass(frequency freq, float sps)
      : a(1.0 - fast_exp3(-2_pi * as_double(freq) / sps))
    {}

    float operator()(float s)
    {
      return y += a * (s - y);
    }

    float operator()() const
    {
      return y;
    }

    one_pole_lowpass& operator=(float y_)
    {
      y = y_;
      return *this;
    }

    void cutoff(frequency freq, float sps)
    {
      a = 1.0 - fast_exp3(-2_pi * as_double(freq) / sps);
    }

    float y = 0.0f, a;
  };

  ////////////////////////////////////////////////////////////////////////////
  // This filter consists of two first order low-pass filters in series,
  // with some of the difference between the two filter outputs fed back to
  // give a resonant peak.
  //
  // See: https://www.w3.org/2011/audio/audio-eq-cookbook.html
  //
  // This is probably faster than the RBJ biquds (see biquad.hpp),
  // especially when computing the coefficients (e.g. when sweeping the
  // frequency) but this filter is rather limited, quirky and inacurate. Use
  // this as a quick and dirty 'musical' filter. E.g. when you don't care
  // about the actual frequency and Q, this one can be swept quickly using a
  // normalized values from 0.0 to less than 1.0, possibly using some custom
  // response curves (e.g. synth filters). Note that _f can't be == 1.0,
  // otherwise, the _fb computation will have a divide by zero. It is also
  // possible to supply the actual frequency, given the sps (samples per
  // second), but take note that the limit is around 7kHz given a sampling
  // rate of 44100, otherwise, the divide by zero.
  //
  // Unless you are in a pinch, I'd rather just use the RBJ biquds.
  ////////////////////////////////////////////////////////////////////////////
  struct reso_filter
  {
    reso_filter(frequency f, float reso, float sps)
      : _f(2.0f * fastsin(pi * as_float(f) / sps))
      , _fb(reso + reso / (1.0f - _f))
      , _reso(reso)
    {}

    reso_filter(float f, float reso)
      : _f(f)
      , _fb(reso + reso / (1.0f - _f))
      , _reso(reso)
    {}

    float operator()(float s)
    {
      _y0 += _f * (s - _y0 + _fb * (_y0 - _y1));
      _y1 += _f * (_y0 - _y1);
      return _y1;
    }

    float operator()() const
    {
      return _y1;
    }

    void cutoff(frequency f, float sps)
    {
      _f = 2.0f * fastsin(pi * as_float(f) / sps);
      _fb = _reso + _reso / (1.0f - _f);
    }

    void cutoff(float f)
    {
      _f = f;
      _fb = _reso + _reso / (1.0f - _f);
    }

    void resonance(float reso)
    {
      _reso = reso;
      _fb = reso + reso / (1.0f - _f);
    }

    float _f, _fb, _reso;
    float _y0 = 0, _y1 = 0;
  };

  ////////////////////////////////////////////////////////////////////////////
  // dynamic_smoother based on Dynamic Smoothing Using Self Modulating Filter
  // by Andrew Simper, Cytomic, 2014, andy@cytomic.com
  //
  //    https://cytomic.com/files/dsp/DynamicSmoothing.pdf
  //
  // A robust and inexpensive dynamic smoothing algorithm based on using the
  // bandpass output of a 2 pole multimode filter to modulate its own cutoff
  // frequency. The bandpass signal is a meaure of how much the signal is
  // "changing" so is useful to increase the cutoff frequency dynamically
  // and allow for faster tracking when the input signal is changing more.
  // The absolute value of the bandpass signal is used since either a change
  // upwards or downwards should increase the cutoff.
  //
  ////////////////////////////////////////////////////////////////////////////
  struct dynamic_smoother
  {
    dynamic_smoother(frequency base, float sps)
      : dynamic_smoother(base, 0.5, sps)
    {}

    dynamic_smoother(frequency base, float sensitivity, float sps)
      : sense(sensitivity * 4.0f)  // efficient linear cutoff mapping
      , wc(as_double(base) / sps)
    {
      auto gc = std::tan(pi * wc);
      g0 = 2.0f * gc / (1.0f + gc);
    }

    float operator()(float s)
    {
      auto lowlz = low1;
      auto low2z = low2;
      auto bandz = lowlz - low2z;
      auto g = std::min(g0 + sense * std::abs(bandz), 1.0f);
      low1 = lowlz + g * (s - lowlz);
      low2 = low2z + g * (low1 - low2z);
      return low2z;
    }

    void base_frequency(frequency base, float sps)
    {
      wc = as_double(base) / sps;
      auto gc = std::tan(pi * wc);
      g0 = 2.0f * gc / (1.0f + gc);
    }

    float sense, wc, g0;
    float low1 = 0.0f;
    float low2 = 0.0f;
  };
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // During the attack phase of an audio signal, the peak envelope follower
  // closely tracks the maximum peak level. When the signal level drops
  // below the peak, the follower gradually releases the peak level with an
  // exponential decay.
  ////////////////////////////////////////////////////////////////////////////
  struct peak_envelope_follower
  {
    peak_envelope_follower(duration release, float sps);

    float                   operator()(float s);
    float                   operator()() const;
    peak_envelope_follower& operator=(float y_);
    void                    release(duration release_, float sps);

    float y = 0.0f, _release;
  };

  ////////////////////////////////////////////////////////////////////////////
  // The ar_envelope_follower follower tracks the envelope of a signal with
  // the given attack parameter and with gradual release given by the
  // release parameter. The signal decays exponentially if the signal is
  // below the peak.
  ////////////////////////////////////////////////////////////////////////////
  struct ar_envelope_follower
  {
    ar_envelope_follower(
      duration attack
      , duration release
      , float sps
    );

    float                   operator()(float s);
    float                   operator()() const;
    ar_envelope_follower& operator=(float y_);
    void                    config(duration attack, duration release, float sps);
    void                    attack(float attack_, float sps);
    void                    release(float release_, float sps);

    float y = 0.0f, _attack, _release;
  };

  using envelope_follower [[deprecated("Use ar_envelope_follower instead.")]]
    = ar_envelope_follower;

    ////////////////////////////////////////////////////////////////////////////
    // This envelope follower combines fast response, low ripple.
    //
    // Based on https://bit.ly/3zmGLIf, https://bit.ly/40OhYJf
    //
    // There is no filtering. The output is a jagged, staircase-like envelope.
    // That way, this can be useful for analysis such as onset detection. For
    // monophonic signals, the hold duration should be equal to or slightly
    // longer than 1/div the period of the lowest frequency of the signal we
    // wish to track, where `div` is the sole template parameter. The hold
    // parameter determines the staircase step duration. This staircase-like
    // envelope can be effectively smoothed out using a moving average filter
    // with the same duration as the hold parameter.
    //
    // fast_envelope_follower is provided, which has div = 2.
    ////////////////////////////////////////////////////////////////////////////
    template <std::size_t div>
    struct basic_fast_envelope_follower
    {
      static_assert(div >= 1, "div must be >= 1");
      static constexpr std::size_t size = div + 1;

      basic_fast_envelope_follower(duration hold, float sps);
      basic_fast_envelope_follower(std::size_t hold_samples);

      float    operator()(float s);
      float    operator()() const;

      std::array<float, size> _y;
      float _peak = 0;
      std::uint16_t _tick = 0, _i = 0;
      std::uint16_t const _reset;
    };

    using fast_envelope_follower = basic_fast_envelope_follower<2>;

    ////////////////////////////////////////////////////////////////////////////
    // This is a fast_envelope_follower followed by a moving average filter to
    // smooth out the staircase ripples as mentioned in the
    // fast_envelope_follower notes.
    ////////////////////////////////////////////////////////////////////////////
    template <std::size_t div>
    struct basic_fast_ave_envelope_follower
    {
      basic_fast_ave_envelope_follower(duration hold, float sps);
      basic_fast_ave_envelope_follower(std::size_t hold_samples);

      float    operator()(float s);
      float    operator()() const;

      basic_fast_envelope_follower<div> _fenv;
      moving_average _ma;
    };

    using fast_ave_envelope_follower = basic_fast_ave_envelope_follower<2>;

    ////////////////////////////////////////////////////////////////////////////
    // This rms envelope follower combines fast response, low ripple using
    // moving RMS detection and the fast_ave_envelope_follower for
    // tracking the moving RMS.
    //
    // The signal path is as follows:
    //    1. Square signal
    //    2. Fast averaging envelope follower
    //    3. Square root.
    //
    // The `fast_rms_envelope_follower_db` variant works in the dB domain,
    // which makes it easy to use as an envelope follower for dynamic range
    // effects (compressor, expander, and agc) which already work in the dB
    // domain, so we eliminate a linear to decibel conversion and optimize
    // computation by using division by 2 instead of sqrt as an added bonus.
    //
    // Designed by Joel de Guzman (June 2020)
    ////////////////////////////////////////////////////////////////////////////
    struct fast_rms_envelope_follower
    {
      constexpr static auto threshold = lin_float(-120_dB);

      fast_rms_envelope_follower(duration hold, float sps);
      float    operator()(float s);

      fast_ave_envelope_follower  _fenv;
    };

    struct fast_rms_envelope_follower_db : fast_rms_envelope_follower
    {
      using fast_rms_envelope_follower::fast_rms_envelope_follower;

      decibel operator()(float s);
    };

    ////////////////////////////////////////////////////////////////////////////
    // Implementation
    ////////////////////////////////////////////////////////////////////////////

    ////////////////////////////////////////////////////////////////////////////
    // peak_envelope_follower
    inline peak_envelope_follower::peak_envelope_follower(duration release, float sps)
      : _release(fast_exp3(-2.0f / (sps * as_double(release))))
    {}

    inline float peak_envelope_follower::operator()(float s)
    {
      if (s > y)
        y = s;
      else
        y = s + _release * (y - s);
      return y;
    }

    inline float peak_envelope_follower::operator()() const
    {
      return y;
    }

    inline peak_envelope_follower& peak_envelope_follower::operator=(float y_)
    {
      y = y_;
      return *this;
    }

    inline void peak_envelope_follower::release(duration release_, float sps)
    {
      _release = fast_exp3(-2.0f / (sps * as_double(release_)));
    }


    ////////////////////////////////////////////////////////////////////////////
    // ar_envelope_follower
    inline ar_envelope_follower::ar_envelope_follower(
      duration attack, duration release, float sps)
      : _attack(fast_exp3(-2.0f / (sps * as_double(attack))))
      , _release(fast_exp3(-2.0f / (sps * as_double(release))))
    {}

    inline float ar_envelope_follower::operator()(float s)
    {
      return y = s + ((s > y) ? _attack : _release) * (y - s);
    }

    inline float ar_envelope_follower::operator()() const
    {
      return y;
    }

    inline ar_envelope_follower& ar_envelope_follower::operator=(float y_)
    {
      y = y_;
      return *this;
    }

    inline void ar_envelope_follower::config(duration attack, duration release, float sps)
    {
      _attack = fast_exp3(-2.0f / (sps * as_double(attack)));
      _release = fast_exp3(-2.0f / (sps * as_double(release)));
    }

    inline void ar_envelope_follower::attack(float attack_, float sps)
    {
      _attack = fast_exp3(-2.0f / (sps * attack_));
    }

    inline void ar_envelope_follower::release(float release_, float sps)
    {
      _release = fast_exp3(-2.0f / (sps * release_));
    }

    ////////////////////////////////////////////////////////////////////////////
    // basic_fast_envelope_follower<div>
    template <std::size_t div>
    inline basic_fast_envelope_follower<div>
      ::basic_fast_envelope_follower(duration hold, float sps)
      : basic_fast_envelope_follower((as_float(hold)* sps))
    {}

    template <std::size_t div>
    inline basic_fast_envelope_follower<div>::
      basic_fast_envelope_follower(std::size_t hold_samples)
      : _reset(hold_samples)
    {
      _y.fill(0);
    }

    template <std::size_t div>
    inline float basic_fast_envelope_follower<div>::operator()(float s)
    {
      // Update _y
      for (auto& y : _y)
        y = std::max(s, y);

      // Reset _y in a round-robin fashion every so often (the hold parameter)
      if (_tick++ == _reset)
      {
        _tick = 0;
        _y[_i++ % size] = 0;
      }

      // The peak is the maximum of _y
      _peak = *std::max_element(_y.begin(), _y.end());
      return _peak;
    }

    template <std::size_t div>
    inline float basic_fast_envelope_follower<div>::operator()() const
    {
      return _peak;
    }

    ////////////////////////////////////////////////////////////////////////////
    // basic_fast_ave_envelope_follower<div>
    template <std::size_t div>
    inline basic_fast_ave_envelope_follower<div>::
      basic_fast_ave_envelope_follower(duration hold, float sps)
      : _fenv(hold, sps)
      , _ma(hold, sps)
    {}

    template <std::size_t div>
    inline basic_fast_ave_envelope_follower<div>::
      basic_fast_ave_envelope_follower(std::size_t hold_samples)
      : _fenv(hold_samples)
      , _ma(hold_samples)
    {}

    template <std::size_t div>
    inline float basic_fast_ave_envelope_follower<div>::operator()(float s)
    {
      return _ma(_fenv(s));
    }

    template <std::size_t div>
    inline float basic_fast_ave_envelope_follower<div>::operator()() const
    {
      return _ma();
    }

    ////////////////////////////////////////////////////////////////////////////
    // fast_rms_envelope_follower
    inline fast_rms_envelope_follower::fast_rms_envelope_follower(
      duration hold, float sps)
      : _fenv(hold, sps)
    {
    }

    inline float fast_rms_envelope_follower::operator()(float s)
    {
      auto e = _fenv(s * s);
      if (e < threshold)
        e = 0;
      return fast_sqrt(e);
    }

    ////////////////////////////////////////////////////////////////////////////
    // fast_rms_envelope_follower
    inline decibel fast_rms_envelope_follower_db::operator()(float s)
    {
      auto e = _fenv(s * s);
      if (e < threshold)
        e = 0;

      // Perform square-root in the dB domain:
      return lin_to_db(e) / 2.0f;
    }
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  class period_detector
  {
  public:

    static constexpr float pulse_threshold = 0.6f;
    static constexpr float harmonic_periodicity_factor = 16.0f;
    static constexpr float periodicity_diff_factor = 0.008f;

    struct info
    {
      float                _period = -1.0f;
      float                _periodicity = 0.0f;
    };

    period_detector(
      frequency lowest_freq
      , frequency highest_freq
      , float sps
      , decibel hysteresis
    );

    period_detector(period_detector const& rhs) = default;
    period_detector(period_detector&& rhs) = default;

    period_detector& operator=(const period_detector& rhs) {
      if (this != &rhs) {
        // Perform deep copy of data members
        _zc = rhs._zc;  // Assuming _zc supports assignment
        _fundamental = rhs._fundamental;
        _min_period = rhs._min_period;
        _range = rhs._range;
        _bits = rhs._bits;
        _weight = rhs._weight;
        _mid_point = rhs._mid_point;
        _period_diff_threshold = rhs._period_diff_threshold;
        _predicted_period = rhs._predicted_period;
        _edge_mark = rhs._edge_mark;
        _predict_edge = rhs._predict_edge;
      }
      return *this;
    }


    bool                    operator()(float s);
    bool                    operator()() const;

    bool                    is_ready() const { return _zc.is_ready(); }
    std::size_t const       minimum_period() const { return _min_period; }
    bitset<> const& bits() const { return _bits; }
    zero_crossing_collector const& edges() const { return _zc; }
    float                   predict_period() const;

    info const& fundamental() const { return _fundamental; }
    float                   harmonic(std::size_t index) const;

  private:

    void                    set_bitstream();
    void                    autocorrelate();
    int                     autocorrelate(bitstream_acf<> const& ac, std::size_t& period, bool first) const;

    zero_crossing_collector           _zc;
    info                    _fundamental;
    std::size_t        _min_period;
    int                     _range;
    bitset<>                _bits;
    float              _weight;
    std::size_t        _mid_point;
    float              _period_diff_threshold;
    mutable float           _predicted_period = -1.0f;
    std::size_t             _edge_mark = 0;
    mutable std::size_t     _predict_edge = 0;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Implementation
  ////////////////////////////////////////////////////////////////////////////
  inline period_detector::period_detector(
    frequency lowest_freq
    , frequency highest_freq
    , float sps
    , decibel hysteresis
  )
    : _zc(hysteresis, as_float(lowest_freq.period() * 2)* sps)
    , _min_period(as_float(highest_freq.period())* sps)
    , _range(as_float(highest_freq) / as_float(lowest_freq))
    , _bits(_zc.window_size())
    , _weight(2.0 / _zc.window_size())
    , _mid_point(_zc.window_size() / 2)
    , _period_diff_threshold(_mid_point* periodicity_diff_factor)
  {
#if __cpp_exceptions == 199711
    if (highest_freq <= lowest_freq)
      throw std::runtime_error(
        "Error: highest_freq <= lowest_freq."
      );
#endif
  }

  inline void period_detector::set_bitstream()
  {
    auto threshold = _zc.peak_pulse() * pulse_threshold;

    _bits.clear();
    for (std::size_t i = 0; i != _zc.num_edges(); ++i)
    {
      auto const& info = _zc[i];
      if (info._peak >= threshold)
      {
        auto pos = std::max<int>(info._leading_edge, 0);
        auto n = info._trailing_edge - pos;
        _bits.set(pos, n, 1);
      }
    }
  }

  namespace detail
  {
    struct sub_collector
    {
      // Intermediate data structure for collecting autocorrelation results
      struct info
      {
        int               _i1 = -1;
        int               _i2 = -1;
        int               _period = -1;
        float             _periodicity = 0.0f;
        std::size_t       _harmonic;
      };

      sub_collector(zero_crossing_collector const& zc, float period_diff_threshold, int range_)
        : _zc(zc)
        , _harmonic_threshold(
          period_detector::harmonic_periodicity_factor * 2 / zc.window_size())
        , _period_diff_threshold(period_diff_threshold)
        , _range(range_)
      {}

      bool empty() const
      {
        return _fundamental._period == -1;
      }

      float period_of(info const& x) const
      {
        auto const& first = _zc[x._i1];
        auto const& next = _zc[x._i2];
        return first.fractional_period(next);
      }

      void save(info const& incoming)
      {
        _fundamental = incoming;
        _fundamental._harmonic = 1;
        _first_period = period_of(_fundamental);
      };

      bool try_sub_harmonic(std::size_t harmonic, info const& incoming)
      {
        int incoming_period = incoming._period / harmonic;
        int current_period = _fundamental._period;
        if (std::abs(incoming_period - current_period) < _period_diff_threshold)
        {
          // If incoming is a different harmonic and has better
          // periodicity ...
          if (incoming._periodicity > _fundamental._periodicity &&
            harmonic != _fundamental._harmonic)
          {
            auto period_diff = std::abs(
              incoming._periodicity - _fundamental._periodicity);

            // If incoming periodicity is within the harmonic
            // periodicity threshold, then replace _fundamental with
            // incoming. Take note of the harmonic for later.
            if (period_diff <= _harmonic_threshold)
            {
              _fundamental._i1 = incoming._i1;
              _fundamental._i2 = incoming._i2;
              _fundamental._periodicity = incoming._periodicity;
              _fundamental._harmonic = harmonic;
            }

            // If not, then we save incoming (replacing the current
            // _fundamental).
            else
            {
              save(incoming);
            }
          }
          return true;
        }
        return false;
      };

      bool process_harmonics(info const& incoming)
      {
        if (incoming._period < _first_period)
          return false;

        int multiple = std::max(1.0f, std::round(period_of(incoming) / _first_period));
        return try_sub_harmonic(std::min(_range, multiple), incoming);
      }

      void operator()(info const& incoming)
      {
        if (_fundamental._period == -1.0f)
          save(incoming);

        else if (process_harmonics(incoming))
          return;

        else if (incoming._periodicity > _fundamental._periodicity)
          save(incoming);
      };

      void get(info const& info, period_detector::info& result)
      {
        if (info._period != -1.0f)
        {
          result =
          {
             period_of(info) / info._harmonic
           , info._periodicity
          };
        }
        else
        {
          result = period_detector::info{};
        }
      }

      float                   _first_period;
      info                    _fundamental;
      zero_crossing_collector const& _zc;
      float const             _harmonic_threshold;
      float const             _period_diff_threshold;
      int const               _range;
    };
  }

  inline int period_detector::autocorrelate(bitstream_acf<> const& ac, std::size_t& period, bool first) const
  {
    auto count = ac(period);
    auto mid = ac._mid_array * bitset<>::value_size;
    auto start = period;

    if (first && count == 0)   // make sure this is not a false correlation
    {
      if (ac(period / 2) == 0)  // oops false correlation!
        return -1;           // flag the return as a false correlation
    }
    else if (period < 32) // Search minimum if the resolution is low
    {
      // Search upwards for the minimum autocorrelation count
      for (auto p = start + 1; p < mid; ++p)
      {
        auto c = ac(p);
        if (c > count)
          break;
        count = c;
        period = p;
      }
      // Search downwards for the minimum autocorrelation count
      for (auto p = start - 1; p > _min_period; --p)
      {
        auto c = ac(p);
        if (c > count)
          break;
        count = c;
        period = p;
      }
    }
    return count;
  }

  inline void period_detector::autocorrelate()
  {
    auto threshold = _zc.peak_pulse() * pulse_threshold;

    CYCFI_ASSERT(_zc.num_edges() > 1, "Not enough edges.");

    bitstream_acf<> ac{ _bits };
    detail::sub_collector collect{ _zc, _period_diff_threshold, _range };

    [&]()
      {
        for (std::size_t i = 0; i != _zc.num_edges() - 1; ++i)
        {
          auto const& first = _zc[i];
          if (first._peak >= threshold)
          {
            for (std::size_t j = i + 1; j != _zc.num_edges(); ++j)
            {
              auto const& next = _zc[j];
              if (next._peak >= threshold)
              {
                auto period = first.period(next);
                if (period > _mid_point)
                  break;
                if (period >= _min_period)
                {
                  auto count = autocorrelate(ac, period, collect.empty());
                  if (count == -1)
                    return; // Return early if we have a false correlation
                  float periodicity = 1.0f - (count * _weight);
                  collect({ int(i), int(j), int(period), periodicity });
                  if (count == 0)
                    return; // Return early if we have perfect correlation
                }
              }
            }
          }
        }
      }();

      // Get the final resuts
      collect.get(collect._fundamental, _fundamental);
  }

  inline bool period_detector::operator()(float s)
  {
    // Zero crossing
    bool prev = _zc();
    bool zc = _zc(s);

    if (!zc && prev != zc)
    {
      ++_edge_mark;
      _predicted_period = -1.0f;
    }

    if (_zc.is_reset())
      _fundamental = info{};

    if (_zc.is_ready())
    {
      set_bitstream();
      autocorrelate();
      return true;
    }
    return false;
  }

  inline float period_detector::harmonic(std::size_t index) const
  {
    if (index > 0)
    {
      if (index == 1)
        return _fundamental._periodicity;

      auto target_period = _fundamental._period / index;
      if (target_period >= _min_period && target_period < _mid_point)
      {
        bitstream_acf<> ac{ _bits };
        auto count = ac(std::round(target_period));
        float periodicity = 1.0f - (count * _weight);
        return periodicity;
      }
    }
    return 0.0f;
  }

  inline bool period_detector::operator()() const
  {
    return _zc();
  }

  inline float period_detector::predict_period() const
  {
    if (_predicted_period == -1.0f && _edge_mark != _predict_edge)
    {
      _predict_edge = _edge_mark;
      if (_zc.num_edges() > 1)
      {
        auto threshold = _zc.peak_pulse() * pulse_threshold;
        for (int i = int(_zc.num_edges()) - 1; i > 0; --i)
        {
          auto const& edge2 = _zc[i];
          if (edge2._peak >= threshold)
          {
            for (int j = i - 1; j >= 0; --j)
            {
              auto const& edge1 = _zc[j];
              if (edge1.similar(edge2))
              {
                _predicted_period = edge1.fractional_period(edge2);
                return _predicted_period;
              }
            }
          }
        }
      }
    }
    return _predicted_period;
  }
}

namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  class pitch_detector
  {
  public:

    static constexpr float  max_deviation = 0.90f;
    static constexpr float  min_periodicity = 0.8f;

    pitch_detector(
      frequency lowest_freq
      , frequency highest_freq
      , float sps
      , decibel hysteresis
    );


    pitch_detector(pitch_detector const& rhs) = default;
    pitch_detector(pitch_detector&& rhs) = default;

    pitch_detector& operator=(const pitch_detector& rhs) {
      if (this != &rhs) {
        // Perform deep copy of data members
        _pd = rhs._pd;
        _frequency = rhs._frequency;
        _median = rhs._median;
        _predict_median = rhs._predict_median;
        _sps = rhs._sps;
        _frames_after_shift = rhs._frames_after_shift;
      }
      return *this;
    }


    bool                    operator()(float s);
    float                   get_frequency() const { return _frequency; }
    float                   predict_frequency(bool init = false);
    bool                    is_note_shift() const;
    std::size_t             frames_after_shift() const { return _frames_after_shift; }
    float                   periodicity() const;
    void                    reset() { _frequency = 0.0f; }

    bitset<> const& bits() const { return _pd.bits(); }
    zero_crossing_collector const& edges() const { return _pd.edges(); }
    period_detector const& get_period_detector() const { return _pd; }

  private:

    float                   calculate_frequency() const;
    float                   bias(float current, float incoming, bool& shift);
    void                    bias(float incoming);

    using exp_moving_average_type = exp_moving_average<2>;

    period_detector         _pd;
    float                   _frequency;
    median3                 _median;
    median3                 _predict_median;
    float                   _sps;
    std::size_t             _frames_after_shift = 0;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Implementation
  ////////////////////////////////////////////////////////////////////////////
  inline pitch_detector::pitch_detector(
    q::frequency lowest_freq
    , q::frequency highest_freq
    , float sps
    , decibel hysteresis
  )
    : _pd{ lowest_freq, highest_freq, sps, hysteresis }
    , _frequency{ 0.0f }
    , _sps{ sps }
  {}

  inline float pitch_detector::bias(float current, float incoming, bool& shift)
  {
    auto error = current / 32; // approx 1/2 semitone
    auto diff = std::abs(current - incoming);

    // Try fundamental
    if (diff < error)
      return incoming;

    // Try harmonics and sub-harmonics
    if (_frames_after_shift > 2)
    {
      if (current > incoming)
      {
        if (int multiple = std::round(current / incoming); multiple > 1)
        {
          auto f = incoming * multiple;
          if (std::abs(current - f) < error)
            return f;
        }
      }
      else
      {
        if (int multiple = std::round(incoming / current); multiple > 1)
        {
          auto f = incoming / multiple;
          if (std::abs(current - f) < error)
            return f;
        }
      }
    }

    // Don't do anything if the latest autocorrelation is not periodic
    // enough. Note that we only do this check on frequency shifts (i.e. at
    // this point, we are looking at a potential frequency shift, after
    // passing through the code above, checking for fundamental and
    // harmonic matches).
    if (_pd.fundamental()._periodicity > min_periodicity)
    {
      // Now we have a frequency shift
      shift = true;
      return incoming;
    }
    return current;
  }

  inline void pitch_detector::bias(float incoming)
  {
    auto current = _frequency;
    ++_frames_after_shift;
    bool shift = false;
    auto f = bias(current, incoming, shift);

    // Don't do anything if incoming is not periodic enough
    // Note that we only do this check on frequency shifts
    if (shift)
    {
      if (_pd.fundamental()._periodicity < max_deviation)
      {
        // If we don't have enough confidence in the autocorrelation
        // result, we'll try the zero-crossing edges to extract the
        // frequency and the one closest to the current frequency wins.
        bool shift2 = false;
        auto predicted = predict_frequency();
        if (predicted > 0.0f)
        {
          float f2 = bias(current, predicted, shift2);

          // If there's no shift, the edges wins
          if (!shift2)
          {
            _frequency = _median(f2);
          }
          else // else, whichever is closest to the current frequency wins.
          {
            bool predicted = std::abs(current - f) >= std::abs(current - f2);
            _frequency = _median(!predicted ? f : f2);
          }
        }
        else
        {
          _frequency = _median(f);
        }
      }
      else
      {
        // Now we have a frequency shift. Get the median of 3 (incoming
        // frequency and last two frequency shifts) to eliminate abrupt
        // changes. This will minimize potentially unwanted shifts.
        // See https://en.wikipedia.org/wiki/Median_filter
        auto f = _median(incoming);
        if (f == incoming)
          _frames_after_shift = 0;
        _frequency = f;
      }
    }
    else
    {
      _frequency = _median(f);
    }
  }

  inline bool pitch_detector::operator()(float s)
  {
    _pd(s);

    if (_pd.is_ready())
    {
      if (_frequency == 0.0f)
      {
        // Disregard if we are not periodic enough
        if (_pd.fundamental()._periodicity >= max_deviation)
        {
          auto f = calculate_frequency();
          if (f > 0.0f)
          {
            _median(f);       // Apply the median for the future
            _frequency = f;   // But assign outright now
            _frames_after_shift = 0;
          }
        }
      }
      else
      {
        if (_pd.fundamental()._periodicity < min_periodicity)
          _frames_after_shift = 0;
        auto f = calculate_frequency();
        if (f > 0.0f)
          bias(f);
      }
    }
    return _pd.is_ready();
  }

  inline float pitch_detector::calculate_frequency() const
  {
    if (_pd.fundamental()._period != -1)
      return _sps / _pd.fundamental()._period;
    return 0.0f;
  }

  inline float pitch_detector::periodicity() const
  {
    return _pd.fundamental()._periodicity;
  }

  inline bool pitch_detector::is_note_shift() const
  {
    return _frames_after_shift == 0;
  }

  inline float pitch_detector::predict_frequency(bool init)
  {
    auto period = _pd.predict_period();
    if (period < _pd.minimum_period())
      return 0.0f;
    auto f = _sps / period;
    if (_frequency != f)
    {
      f = _predict_median(f);
      if (init)
        _frequency = _median(f);
    }
    return f;
  }
}





namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // compressor dynamically modulates the gain when the signal envelope
  // rises above a specified threshold. Envelope tracking is done using an
  // external envelope follower to make it possible to use different types
  // of envelope tracking schemes, the output of which is the supplied 'env'
  // argument to the function call operator operator()(decibel env) where
  // env is the envelope of the signal in decibels obtained (e.g) using the
  // envelope_follower.
  //
  // Note that these are envelope processors that work in the logarithmic
  // domain (decibels) and process envelopes instead of the actual signal.
  // The output is the compressed envelope, also in decibels. Simply
  // multiply the signal by the result converted to float (or double). For
  // example:
  //
  //    auto gain = as_float(comp(env));
  //    auto left_out = left_signal * gain;
  //    auto right_out = right_signal * gain;
  //
  // where left_signal and right_signal are stereo input signals and
  // envelope is the computed envelope (e.g) using an envelope follower.
  //
  // The ratio parameter specifies the amount of gain applied. With the
  // typical "n:1" notation for compressors, the ratio parameter is 1/n,
  // thereby the ratio for compressors is normally from 0.0...1.0. (e.g. 4:
  // 1 compression is 1/4 or 0.25). Signal rising above the threshold is
  // attenuated, compressing the signal. For every dB above the threshold,
  // the signal is attenuated by n dB. For example, with a ratio of 4:1
  // (0.25), 1dB above the threshold is attenuated by 4dB.
  //
  // Typically, you add some makeup gain after compression to compensate for
  // the gain reduction.
  ////////////////////////////////////////////////////////////////////////////
  struct compressor
  {
    compressor(decibel threshold, float ratio);

    decibel     operator()(decibel env) const;
    void        threshold(decibel val);
    void        ratio(float ratio);
    decibel     threshold() const;
    float       ratio() const;

    decibel     _threshold;
    float       _slope;
  };

  ////////////////////////////////////////////////////////////////////////////
  // The soft_knee_compressor compressor variant provides a more gradual
  // "soft knee" gain transition around the threshold, given a knee width.
  // See compressor above.
  ////////////////////////////////////////////////////////////////////////////
  struct soft_knee_compressor
  {
    soft_knee_compressor(decibel threshold, decibel width, float ratio);

    decibel     operator()(decibel env) const;
    void        threshold(decibel val);
    void        width(decibel val);
    void        ratio(float ratio);
    decibel     threshold() const;
    decibel     width() const;
    float       ratio() const;

    decibel     _threshold, _width, _lower, _upper;
    float       _slope;
  };

  ////////////////////////////////////////////////////////////////////////////
  // The expander is the inverse of the compressor. The expander adjusts the
  // gain when the signal falls below the threshold, attenuating the signal.
  // With the typical "1:n" notation for expanders, the ratio parameter is
  // n, thereby the ratio for expanders is normally from 0.0...inf. (e.g.
  // 1 : 4 expansion is 4). A ratio of 1 : inf is a hard gate where no
  // signal passes below the threshold.
  //
  // For every dB below the threshold, the signal is attenuated by n dB. For
  // example, with a ratio of 4 : 1 (4), 1dB below the threshold is
  // attenuated by 4dB.
  ////////////////////////////////////////////////////////////////////////////
  struct expander
  {
    expander(decibel threshold, float ratio);

    decibel     operator()(decibel env) const;
    void        threshold(decibel val);
    void        ratio(float ratio);
    decibel     threshold() const;
    float       ratio() const;

    decibel     _threshold;
    float       _slope;
  };

  ////////////////////////////////////////////////////////////////////////////
  // The agc (automatic gain control) compares the envelope, env, to a
  // reference, ref, and increases or decreases the gain to maintain a
  // constant output level. A Maximum gain that can be applied when the
  // signal falls below the reference. The max constructor parameter
  // specified this "maximum" gain.
  ////////////////////////////////////////////////////////////////////////////
  struct agc
  {
    agc(decibel max);

    decibel     operator()(decibel env, decibel ref) const;
    void        max(decibel max_);
    decibel     max() const;

    decibel     _max;
  };

  ////////////////////////////////////////////////////////////////////////////
  // Implementation
  ////////////////////////////////////////////////////////////////////////////

  ////////////////////////////////////////////////////////////////////////////
  // compressor
  inline compressor::compressor(decibel threshold, float ratio)
    : _threshold(threshold)
    , _slope(1.0f - ratio)
  {}

  inline decibel compressor::operator()(decibel env) const
  {
    if (env <= _threshold)
      return 0; //db
    return _slope * (_threshold - env);
  }

  inline void compressor::threshold(decibel val)
  {
    _threshold = val;
  }

  inline void compressor::ratio(float ratio)
  {
    _slope = 1.0f - ratio;
  }

  inline decibel compressor::threshold() const
  {
    return _threshold;
  }

  inline float compressor::ratio() const
  {
    return 1.0f - _slope;
  }

  ////////////////////////////////////////////////////////////////////////////
  // soft_knee_compressor
  inline soft_knee_compressor::soft_knee_compressor(
    decibel threshold, decibel width, float ratio)
    : _threshold(threshold)
    , _width(width)
    , _lower(threshold - (_width * 0.5))
    , _upper(threshold + (_width * 0.5))
    , _slope(1.0f - ratio)
  {}

  inline decibel soft_knee_compressor::operator()(decibel env) const
  {
    if (env <= _lower)
    {
      return 0_dB;
    }
    else if (env <= _upper)
    {
      auto soft_slope = _slope * ((env - _lower) / _width) * 0.5;
      return soft_slope * (_lower - env);
    }
    else
    {
      return _slope * (_threshold - env);
    }
  }

  inline void soft_knee_compressor::threshold(decibel val)
  {
    _threshold = val;
    _lower = _threshold - (_width * 0.5);
    _upper = _threshold + (_width * 0.5);
  }

  inline void soft_knee_compressor::width(decibel val)
  {
    _width = val;
    _lower = _threshold - (_width * 0.5);
    _upper = _threshold + (_width * 0.5);
  }

  inline void soft_knee_compressor::ratio(float ratio)
  {
    _slope = 1.0f - ratio;
  }

  inline decibel soft_knee_compressor::threshold() const
  {
    return _threshold;
  }

  inline decibel soft_knee_compressor::width() const
  {
    return _width;
  }

  inline float soft_knee_compressor::ratio() const
  {
    return 1.0f - _slope;
  }

  ////////////////////////////////////////////////////////////////////////////
  // expander
  inline expander::expander(decibel threshold, float ratio)
    : _threshold(threshold)
    , _slope(ratio)
  {}

  inline decibel expander::operator()(decibel env) const
  {
    if (env >= _threshold)
      return 0_dB;
    return _slope * (env - _threshold);
  }

  inline void expander::threshold(decibel val)
  {
    _threshold = val;
  }

  inline void expander::ratio(float ratio)
  {
    _slope = ratio;
  }

  inline decibel expander::threshold() const
  {
    return _threshold;
  }

  inline float expander::ratio() const
  {
    return _slope;
  }

  ////////////////////////////////////////////////////////////////////////////
  // agc
  inline agc::agc(decibel max)
    : _max(max)
  {}

  inline decibel agc::operator()(decibel env, decibel ref) const
  {
    auto g = ref - env;
    return std::min<decibel>(g, _max);
  }

  inline void agc::max(decibel max_)
  {
    _max = max_;
  }

  inline decibel agc::max() const
  {
    return _max;
  }
}


namespace cycfi::q
{
  ////////////////////////////////////////////////////////////////////////////
  // clip a signal to range -_max...+_max
  //
  //    max: maximum value
  ////////////////////////////////////////////////////////////////////////////
  struct clip
  {
    constexpr clip(decibel max)
      : _max(lin_float(max))
    {}

    constexpr clip(float max = 1.0f)
      : _max(max)
    {}

    constexpr float operator()(float s) const
    {
      return (s > _max) ? _max : (s < -_max) ? -_max : s;
    }

    float _max;
  };

  ////////////////////////////////////////////////////////////////////////////
  // soft_clip a signal to range -1.0 to 1.0.
  ////////////////////////////////////////////////////////////////////////////
  struct soft_clip : clip
  {
    constexpr float operator()(float s) const
    {
      s = clip::operator()(s);
      return 1.5 * s - 0.5 * s * s * s;
    }
  };
}

namespace q = cycfi::q;
using namespace q::literals;

float amalgamated_q_processBlock(const float* inBlock, int blockSize, float sampleRate, float lowestFreq, float highestFreq)
{
  q::frequency lowest_freq = lowestFreq;
  q::frequency highest_freq = highestFreq;

  ////////////////////////////////////////////////////////////////////////////
  // Process
  thread_local q::pitch_detector          pd{ lowest_freq, highest_freq, sampleRate, -40.0 }; // db

  { // workaround to set changed frequency. not sure if there is a better way
    thread_local float last_lowestFreq = lowestFreq;
    thread_local float last_highestFreq = highestFreq;
    if (last_lowestFreq != lowestFreq || last_highestFreq != highestFreq)
    {
      last_lowestFreq = lowestFreq;
      last_highestFreq = highestFreq;
      pd = q::pitch_detector { lowest_freq, highest_freq, sampleRate, -40.0 }; // db
    }
  }


  auto const& bits = pd.bits();
  auto const& edges = pd.edges();
  q::bitstream_acf<>         bacf{ bits };
  auto                       min_period = as_float(highest_freq.period()) * sampleRate;

  thread_local q::peak_envelope_follower  env{ 30_ms_, sampleRate };
  thread_local q::one_pole_lowpass        lp{ highest_freq, sampleRate };
  thread_local q::one_pole_lowpass        lp2{ lowest_freq, sampleRate };

  constexpr float            slope = 1.0f / 4;
  constexpr float            makeup_gain = 4;
  thread_local q::compressor              comp{ -18.0, slope };// db
  thread_local q::clip                    clip;

  thread_local float                      onset_threshold = q::lin_float(-28.0); // db
  thread_local float                      release_threshold = q::lin_float(-60.0); // db
  thread_local float                      threshold = onset_threshold;

  thread_local bool loudEnough = false;
  thread_local float frequency = 0.0f;

  for (auto i = 0; i < blockSize; ++i)
  {
    float time = i / float(sampleRate);

    auto s = inBlock[i];

    // Bandpass filter
    s = lp(s);
    s -= lp2(s);

    // Envelope
    auto e = env(std::abs(s));
    auto e_db = q::lin_to_db(e);

    if (e > threshold)
    {
      // Compressor + makeup-gain + hard clip
      auto gain = q::lin_float(comp(e_db)) * makeup_gain;
      s = clip(s * gain);
      threshold = release_threshold;
      loudEnough = true;
    }
    else
    {
      s = 0.0f;
      threshold = onset_threshold;
      loudEnough = false;
    }

    // Pitch Detect
    const bool ready = pd(s);

    if (ready)
    {
      frequency = pd.get_frequency();
    }
  }

  if (loudEnough)
    return frequency;
  return {};
}
