#ifndef SPLINE_H
#define SPLINE_H
/** @file
 *  @class SPLINE
 *  @brief Kubické splajny.
 */
#include "real_fix.h"
using namespace FIX;
typedef __SIZE_TYPE__ size_t;

// Zjištění počtu prvků statického pole.
template<class T, size_t N>constexpr size_t array_size (T (&) [N]) { return N; }
/*!
 * @class Pair
 * @brief Prvek vstupní tabulky pro výpočet kubických splajnů.
 * */
struct Pair {
  double x, y;
};
/*! @class SplineSet
    @brief Výsledné koeficienty pro výpočet iterace pomocí polynomu 3. řádu.
*/
struct SplineSet {
  real x;   //!< začátky intervalů x
  real a;   //!< záčátky intervalů y, koeficienty a
  real b;   //!< koeficienty b (x^1)
  real c;   //!< koeficienty c (x^2)
  real d;   //!< koeficienty d (x^3)
};
//! V podstatě statický vektor pro SplineSet (bez nutnosti dymanické alokace)
template<const int N> class SPLINE {
  SplineSet data [N - 1];   // zde jsou schovány koeficienty (privátně)
  public:
    //! Konstruktor
    explicit constexpr SPLINE (const Pair * const p, const bool reverse = false) noexcept {
      double x [N], y [N];  // přeskupení dat - možnost inverzní funkce při reverse = true
      if (reverse) {for (int i=0; i<N; i++) { x[i] = p[i].y; y[i] = p[i].x; }}
      else         {for (int i=0; i<N; i++) { x[i] = p[i].x; y[i] = p[i].y; }}      
      /* Tenhle příšerně složitý konstruktor je převzat ze stackoverflow.com (původní odkaz už zmizel)
       * Není to zas taková sranda - koeficienty se počítají řešením tridiagonální matice řádu N.
       * V těch indexech se snadno zabloudí, tohle kupodivu funguje (přirozené splajny) ač je to dost krátké.
       * Je zajímavé, že tohle najdete spíš ve Fortranu než v C, matematici jsou zřejmě hodně konzervativní.
       * */
      const int n = N - 1;
      double h [n];
      for (int i = 0; i < n; ++i) { h[i] = x[i+1]-x[i]; }
      double alpha [n]; alpha [0] = 0.0;
      for (int i = 1; i < n; ++i) {
        alpha[i] = 3.0 * (y[i+1] - y[i]) / h[i] - 3.0 * (y[i] - y[i-1]) / h[i-1];
      }
      double c [n+1], l [n+1], mu [n+1], z [n+1]; l [0] = 1.0;  mu[0] = 0.0;  z [0] = 0.0;
      for (int i = 1; i < n; ++i) {     // přímý chod
        l [i] = 2.0 * (x[i+1] - x[i-1]) - h[i-1] * mu[i-1];
        mu[i] = h[i] / l[i];
        z [i] = (alpha[i] - h[i-1] * z[i-1]) / l[i];
      }
      l[n] = 1.0;  z[n] = 0.0;  c[n] = 0.0;
      double b [n], d [n];
      for (int j = n-1; j >= 0; --j) {  // zpětný chod
        c[j] = z [j] - mu[j] * c[j+1];
        b[j] = (y[j+1] - y[j]) / h[j] - h[j] * (c[j+1] + 2*c[j]) / 3.0;
        d[j] = (c[j+1] - c[j]) / (3.0 * h[j]);
      }
      for (int i = 0; i < n; ++i) {     // závěrečný zápis koeficientů
        // Pro FIX je potřeba nastavit různé posuny pro koeficienty, rozsah je příliš velký a pak by to bylo nepřesné.
        data[i].x = to_real(x[i], 16);     data[i].a = to_real(y[i], 16);
        data[i].b = to_real(b[i], 20);     data[i].c = to_real(c[i], 28); data[i].d = to_real(d[i], 36);
      }
    }
    const SplineSet & operator[] (const int index) const {
      return data [index];
    }
    const size_t size () const {
      return N - 1;
    }
    const real interpolate (const real x) const {
      const unsigned    n = range (x);                     // určíme interval ve kterém budeme počítat
      const SplineSet & s = data [n];                      // koeficienty polynomu v tomto intervalu
      const real        r = x - s.x;                       // souřadnice x v tomto intervalu
      return cubic_poly (s, r);                            // vlastní výpočet
    }
    /***************************************************************/
    /** @class iterator
     *  @brief pro range-based for () */
    class iterator {
      const SplineSet * ptr;
      public:
        iterator(const SplineSet * _ptr) : ptr (_ptr) {}
        iterator  operator++ ()                             { ++ptr;   return * this;    }
        bool      operator!= (const iterator & other) const { return   ptr != other.ptr; }
        const SplineSet & operator*  ()               const { return * ptr;              }
    };
    iterator begin () const { return iterator (data        ); }
    iterator end   () const { return iterator (data + N - 1); }
  protected:
    const unsigned range (const real x) const {
      int l=0, r=size() - 1u;
      while (l <= r) {
        const int s = (l + r) >> 1;
        const real x0 = data[s + 0].x;
        const real x1 = data[s + 1].x;
        if      (x <  x0) r = s - 1;
        else if (x >= x1) l = s + 1;
        else return s;
      }
      return size() - 1u;
    }
    const real cubic_poly (const SplineSet & s, const real x) const {
      real r(s.d);      // d => 36               v konstruktoru
      mull (r, x, 24);  // 36+16-24 = 28 == c
      r += s.c;         // c => 28
      mull (r, x, 24);  // 28+16-24 = 20 == b
      r += s.b;         // b => 20
      mull (r, x, 20);  // 20+16-20 = 16 == a
      r += s.a;         // a=x => 16
      return r;         // r=y => 16
    }
};

#endif // SPLINE_H