Arithmetic.hpp

#pragma once

#include <algorithm>

#include <amdis/common/Math.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/operations/Basic.hpp>
#include <amdis/operations/Composer.hpp>

namespace AMDiS
{
  namespace Operation
  {
    struct Plus
    {
      template <class... Ts>
      constexpr auto operator()(Ts const&... ts) const
      {
        return Math::sum(ts...);
      }
    };

#ifndef DOXYGEN
    // [0] + g => g
    template <class G>
    struct ComposerBuilder<Plus, Zero, G>
      : ComposerBuilder<Id, G> {};

    // g + [0] => g
    template <class G>
    struct ComposerBuilder<Plus, G, Zero>
      : ComposerBuilder<Id, G> {};

    // [0] + [0] => [0]
    template <>
    struct ComposerBuilder<Plus, Zero, Zero>
      : ComposerBuilder<Id, Zero> {};
#endif

    template <class... Int>
    constexpr int order(Plus const&, Int... orders)
    {
      return Math::max(int(orders)...);
    }

    template <std::size_t I>
    constexpr auto partial(Plus const&, index_t<I>)
    {
      static_assert((I < 2), "Derivatives of `Plus` only defined for the binary case.");
      return One{};
    }

    // -------------------------------------------------------------------------

    struct Minus
    {
      template <class T, class S>
      constexpr auto operator()(T const& lhs, S const& rhs) const
      {
        return lhs - rhs;
      }

      friend constexpr int order(Minus const&, int lhs, int rhs)
      {
        return Math::max(lhs, rhs);
      }

      friend constexpr auto partial(Minus const&, index_t<0>)
      {
        return One{};
      }

      friend constexpr auto partial(Minus const&, index_t<1>)
      {
        return StaticConstant<int,-1>{};
      }
    };

#ifndef DOXYGEN
    // g - [0] => g
    template <class G>
    struct ComposerBuilder<Minus, G, Zero>
      : ComposerBuilder<Id, G> {};

    // [0] - [0] => [0]
    template <>
    struct ComposerBuilder<Minus, Zero, Zero>
      : ComposerBuilder<Id, Zero> {};
#endif

    // -------------------------------------------------------------------------

    struct Multiplies
    {
      template <class... Ts>
      constexpr auto operator()(Ts const&... ts) const
      {
        return (ts * ...);
      }
    };


#ifndef DOXYGEN
    // [0] * g => [0]
    template <class G>
    struct ComposerBuilder<Multiplies, Zero, G>
      : ComposerBuilder<Id, Zero> {};

    // g * [0] => [0]
    template <class G>
    struct ComposerBuilder<Multiplies, G, Zero>
      : ComposerBuilder<Id, Zero> {};

    // [0] * [0] => [0]
    template <>
    struct ComposerBuilder<Multiplies, Zero, Zero>
      : ComposerBuilder<Id, Zero> {};
#endif


    template <class... Int>
    constexpr int order(Multiplies const&, Int... orders)
    {
      return Math::sum(int(orders)...);
    }

    // only for binary *
    // d_0 (x * y) = y, d_1 (x * y) = x
    template <std::size_t I>
    constexpr auto partial(Multiplies const&, index_t<I>)
    {
      static_assert((I < 2), "Derivatives of `Multiplies` only defined for the binary case.");
      return Arg<1-I>{};
    }

    // -------------------------------------------------------------------------

    template <class Factor>
    struct Scale
    {
      Factor const factor_;

      template <class T>
      constexpr auto operator()(T const& value) const
      {
        return factor_ * value;
      }

      friend constexpr int order(Scale const&, int d)
      {
        return d;
      }

      // d_x f*x = f
      friend constexpr auto partial(Scale const& s, index_t<0>)
      {
        return Constant<Factor>{s.factor_};
      }
    };

#ifndef DOXYGEN
    // f*(g*A) = (f*g)*A
    template <class F, class G>
    struct ComposerBuilder<Scale<F>, Scale<G>>
    {
      using type = Scale<std::common_type_t<F,G>>;

      template <class F_, class G_>
      static constexpr type build(F_&& f, G_&& g)
      {
        return type{f.factor_ * g.factor_};
      }
    };
#endif

    // -------------------------------------------------------------------------

    struct Divides
    {
      template <class T, class S>
      constexpr auto operator()(T const& lhs, S const& rhs) const
      {
        return lhs / rhs;
      }

      // d_0 f(x,y) = 1 / y
      friend constexpr auto partial(Divides const&, index_t<0>)
      {
        return compose(Divides{}, One{}, Arg<1>{});
      }

      // d_1 f(x,y) = (y - x)/y^2
      friend constexpr auto partial(Divides const&, index_t<1>);
    };

    // -------------------------------------------------------------------------

    struct Negate
    {
      template <class T>
      constexpr auto operator()(T const& x) const
      {
        return -x;
      }

      friend constexpr int order(Negate const&, int d)
      {
        return d;
      }

      friend constexpr auto partial(Negate const&, index_t<0>)
      {
        return StaticConstant<int,-1>{};
      }
    };

#ifndef DOXYGEN
    // g + -h  => g - h
    template <class G>
    struct ComposerBuilder<Plus, G, Negate>
      : ComposerBuilder<Minus, G, Id> {};

    // [0] - g => -g
    template <class G>
    struct ComposerBuilder<Minus, Zero, G>
      : ComposerBuilder<Id, Negate> {};

    // -(g - h) => (h - g)
    template <>
    struct ComposerBuilder<Negate, Minus>
      : ComposerBuilder<Minus, Arg<1>, Arg<0>> {};
#endif

    // -------------------------------------------------------------------------

    // forward declaration
    template <int p, bool positive>
    struct PowImpl;

    template <int p>
    struct PowType
    {
      using type = PowImpl<p, (p>0)>;
    };

    template <> struct PowType<1> { using type = Id; };
    template <> struct PowType<0> { using type = One; };

    template <int p>
    using Pow = typename PowType<p>::type;

    using Sqr = Pow<2>;

    template <int p>
    struct PowImpl<p, true>
    {
      template <class T>
      constexpr auto operator()(T const& x) const
      {
        return Math::pow<p>(x);
      }

      friend constexpr int order(PowImpl const&, int d)
      {
        return p*d;
      }

      friend constexpr auto partial(PowImpl const&, index_t<0>)
      {
        return compose(Multiplies{}, StaticConstant<int,p>{}, Pow<p-1>{});
      }
    };

    template <int p>
    struct PowImpl<p, false>
      : public Composer<Divides, One, Pow<-p>>
    {
      constexpr PowImpl()
        : Composer<Divides, One, Pow<-p>>{}
      {}

      template <class T>
      constexpr auto operator()(T const& x) const
      {
        return T(1) / Math::pow<-p>(x);
      }
    };


#ifndef DOXYGEN
    // (x^p1)^p2  => x^(p1*p2)
    template <int p1, bool pos1, int p2, bool pos2>
    struct ComposerBuilder<PowImpl<p1,pos1>, PowImpl<p2,pos2>>
      : ComposerBuilder<Id, Pow<p1*p2>> {};

    // x^p * y^p  => (x*y)^p
    template <int p, bool pos>
    struct ComposerBuilder<Multiplies, PowImpl<p,pos>, PowImpl<p,pos>>
      : ComposerBuilder<Pow<p>, Multiplies> {};
#endif

    // d_1 f(x,y) = (y - x)/y^2
    inline constexpr auto partial(Divides const&, index_t<1>)
    {
      return compose(Divides{}, compose(Minus{}, Arg<1>{}, Arg<0>{}),
                                compose(Pow<2>{}, Arg<1>{}));
    }

    struct Power
    {
      constexpr Power(double p)
        : p_(p)
      {}

      template <class T>
      auto operator()(T const& x) const
      {
        return std::pow(x, p_);
      }

      friend constexpr auto partial(Power const& P, index_t<0>)
      {
        return compose(Multiplies{}, Constant<double>{P.p_}, Power{P.p_-1});
      }

      double p_;
    };

  } // end namespace Operation
} // end namespace AMDiS