#### RemiAudio
#### Copyright (C) 2022-2024 Remilia Scarlet <remilia@posteo.jp>
####
#### This program is free software: you can redistribute it and/or modify it
#### under the terms of the GNU Affero General Public License as published by
#### the Free Software Foundation, either version 3 of the License, or (at your
#### option) any later version.
####
#### This program is distributed in the hope that it will be useful, but WITHOUT
#### ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
#### FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public
#### License for more details.
####
#### You should have received a copy of the GNU Affero General Public License
#### along with this program.  If not, see <https://www.gnu.org/licenses/>.

module RemiAudio::Interpolate
  extend self

  # An interpolation type.
  enum Mode
    Linear
    Cubic
    Hermite
    HermiteAlt
    BSpline
    Parabolic2X
    Optimal2X4P2O
    Optimal32X4P4O
  end

  # Performs linear interpolation.
  @[AlwaysInline]
  def linear(x1 : Float64, x2 : Float64, mu : Float64) : Float64
    (x1 + (mu * (x2 - x1)))
  end

  # Performs cubic interpolation.
  @[AlwaysInline]
  def cubic(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
    c1 = (x2 - x0) *  0.5
    c3 = (x1 - x2) * 1.5 + (x3 - x0) * 0.5
    c2 = c1 + x0 - x1 - c3
    ((c3 * mu + c2) * mu + c1) * mu + x1
  end

  # Performs hermite interpolation.
  #
  # This is an alternate version with slightly different coefficients.
  @[AlwaysInline]
  def hermiteAlt(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
    mu2 = mu * mu
    c0 = 3 * x1 - 3 * x2 + x3 - x0
    c1 = 2 * x0 - 5 * x1 + 4 * x2 - x3
    c2 = x2 - x0
    c3 = 2 * x1

    (c0 * mu * mu2 + c1 * mu2 + c2 * mu + c3) / 2
  end

  ### These below are taken from the "Elephant Paper".
  ###
  ### http://yehar.com/blog/wp-content/uploads/2009/08/deip.pdf

  # Performs B-Spline interpolation.
  @[AlwaysInline]
  def bspline(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
    ym1py1 = x0 + x2

    c0 = 1 / 6.0 * ym1py1 + 2 / 3.0 * x1
    c1 = 1 / 2.0 * (x2 - x0)
    c2 = 1 / 2.0 * ym1py1 - x1
    c3 = 1 / 2.0 * (x1 - x2) + 1 / 6.0 * (x3 - x0)
    ((c3 * mu + c2) * mu + c1) * mu + c0
  end

  # Performs hermite interpolation.
  @[AlwaysInline]
  def hermite(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
    # Goddess do I miss prefix notation
    c0 = x1
    c1 = 0.5 * (x2 - x0)
    c2 = x0 - (2.5 * x1) + (2 * x2) - (0.5 * x3)
    c3 = (0.5 * (x3 - x0)) + (1.5 * (x1 - x2))
    (((((c3 * mu) + c2) * mu) + c1) * mu) + c0
  end

  # Performs parabolic interpolation.  This expects 2x oversampled data.
  @[AlwaysInline]
  def parabolic2x(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
    y1mym1 : Float64 = x2 - x0
    c0 : Float64 = (0.5 * x1) + (0.25 * (x0 + x2))
    c1 : Float64 = 0.5 * y1mym1
    c2 : Float64 = 0.25 * (x1 - x1 - y1mym1)
    (((c2 * mu) + c1) * mu) + c0
  end

  # Performs optimal 4-point 2nd order interpolation.  This expects 2x
  # oversampled data.
  @[AlwaysInline]
  def optimal2x4P2O(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
    z : Float64 = mu - 0.5
    even1 : Float64 = x2 + x1
    odd1 : Float64 = x2 - x1
    even2 : Float64 = x3 + x0
    odd2 : Float64 = x3 - x0
    c0 : Float64 = (even1 * 0.42334633257225274_f64) + (even2 * 0.07668732202139628_f64)
    c1 : Float64 = (odd1 * 0.26126047291143606_f64) + (odd2 * 0.24778879018226652_f64)
    c2 : Float64 = (even1 * -0.213439787561776841_f64) + (even2 * 0.21303593243799016_f64)
    (((c2 * z) + c1) * z) + c0
  end

  # Performs optimal 4-point 4th order interpolation.  This expects 32x
  # oversampled data.
  @[AlwaysInline]
  def optimal32x4P4O(x0 : Float64, x1 : Float64, x2 : Float64, x3 : Float64, mu : Float64) : Float64
    z : Float64 = mu - 0.5
    even1 : Float64 = x2 + x1
    odd1 : Float64 = x2 - x1
    even2 : Float64 = x3 + x0
    odd2 : Float64 = x3 - x0
    c0 = (even1 * 0.46835497211269561_f64) + (even2 * 0.03164502784253309_f64)
    c1 = (odd1 * 0.56001293337091440_f64) + (odd2 * 0.14666238593949288_f64)
    c2 = (even1 * -0.250038759826233691_f64) + (even2 * 0.25003876124297131_f64)
    c3 = (odd1 * -0.49949850957839148_f64) + (odd2 * 0.16649935475113800_f64)
    c4 = (even1 * 0.00016095224137360_f64) + (even2 * -0.00016095810460478_f64)

    (((c4 * z + c3) * z + c2) * z + c1) * z + c0
  end
end
