Rack/plugins/community/repos/LindenbergResearch/src/dsp/Horner.h

/*
  Implementation of the Horner method of polynomial evaluation

  Copyright (C) 2011 Darko Veberic, darko.veberic@ijs.si

  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU 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 General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once


namespace dsp {

    template<typename Float, class Tag, unsigned int order>
    struct Polynomial {
        static inline Float Coeff();
        static inline Float Coeff(const Float x);
    };


    template<typename Float, class Tag, unsigned int order>
    struct Horner {
        static Float Recurse(const Float term, const Float x) {
            return Horner<Float, Tag, order - 1>::Recurse(term * x + Polynomial<Float, Tag, order>::Coeff(), x);
        }


        static Float RecurseAlt(const Float term, const Float x) {
            return Horner<Float, Tag, order - 1>::RecurseAlt(Polynomial<Float, Tag, order>::Coeff() + x * term, x);
        }


        static Float Eval(const Float x) { return Horner<Float, Tag, order - 1>::Recurse(Polynomial<Float, Tag, order>::Coeff(), x); }


        static Float EvalAlt(const Float x) { return Horner<Float, Tag, order - 1>::RecurseAlt(Polynomial<Float, Tag, order>::Coeff(), x); }


        //
        static Float Recurse(const Float term, const Float x, const Float y) {
            return Horner<Float, Tag, order - 1>::Recurse(term * x + Polynomial<Float, Tag, order>::Coeff(y), x, y);
        }


        static Float Eval(const Float x, const Float y) {
            return Horner<Float, Tag, order - 1>::Recurse(Polynomial<Float, Tag, order>::Coeff(y), x, y);
        }
    };


    template<typename Float, class Tag>
    struct Horner<Float, Tag, 0> {
        static Float Recurse(const Float term, const Float x) { return term * x + Polynomial<Float, Tag, 0>::Coeff(); }


        static Float Eval(const Float x) { return Polynomial<Float, Tag, 0>::Coeff(); }


        //
        static Float Recurse(const Float term, const Float x, const Float y) { return term * x + Polynomial<Float, Tag, 0>::Coeff(y); }


        static Float Eval(const Float x, const Float y) { return Polynomial<Float, Tag, 0>::Coeff(y); }
    };


#define HORNER_COEFF(_Tag_, _i_, _c_)  \
  template<typename Float> struct Polynomial<Float, _Tag_, _i_> { static Float Coeff() { return Float(_c_); } }
#define HORNER_COEFF2(_Tag_, _i_, _c_y_)  \
  template<typename Float> struct Polynomial<Float, _Tag_, _i_> { static Float Coeff(const Float y) { return Float(_c_y_); } }


#define HORNER0(F, x, c0)                                                (F)(c0)
#define HORNER1(F, x, c1, c0)                                 HORNER0(F, x, (F)(c1)*(x) + (F)(c0)                                )
#define HORNER2(F, x, c2, c1, c0)                             HORNER1(F, x, (F)(c2)*(x) + (F)(c1),                             c0)
#define HORNER3(F, x, c3, c2, c1, c0)                         HORNER2(F, x, (F)(c3)*(x) + (F)(c2),                         c1, c0)
#define HORNER4(F, x, c4, c3, c2, c1, c0)                     HORNER3(F, x, (F)(c4)*(x) + (F)(c3),                     c2, c1, c0)
#define HORNER5(F, x, c5, c4, c3, c2, c1, c0)                 HORNER4(F, x, (F)(c5)*(x) + (F)(c4),                 c3, c2, c1, c0)
#define HORNER6(F, x, c6, c5, c4, c3, c2, c1, c0)             HORNER5(F, x, (F)(c6)*(x) + (F)(c5),             c4, c3, c2, c1, c0)
#define HORNER7(F, x, c7, c6, c5, c4, c3, c2, c1, c0)         HORNER6(F, x, (F)(c7)*(x) + (F)(c6),         c5, c4, c3, c2, c1, c0)
#define HORNER8(F, x, c8, c7, c6, c5, c4, c3, c2, c1, c0)     HORNER7(F, x, (F)(c8)*(x) + (F)(c7),     c6, c5, c4, c3, c2, c1, c0)
#define HORNER9(F, x, c9, c8, c7, c6, c5, c4, c3, c2, c1, c0) HORNER8(F, x, (F)(c9)*(x) + (F)(c8), c7, c6, c5, c4, c3, c2, c1, c0)

}