Arithmetic.hpp 7.17 KB
Newer Older
1
2
3
4
#pragma once

#include <algorithm>

5
6
7
8
#include <amdis/common/Math.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/operations/Basic.hpp>
#include <amdis/operations/Composer.hpp>
9
10
11
12
13
14
15
16
17

namespace AMDiS
{
  namespace Operation
  {
    /** \addtogroup operations
     *  @{
     **/

18
    /// Functor that represents A+B
19
20
21
22
23
    struct Plus
    {
      template <class... Ts>
      constexpr auto operator()(Ts const&... ts) const
      {
24
        return Math::sum(ts...);
25
26
27
      }
    };

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

34
    // g + [0] => g
35
36
    template <class G>
    struct ComposerBuilder<Plus, G, Zero>
37
      : ComposerBuilder<Id, G> {};
38
39
40
41
42

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

    template <class... Int>
46
    constexpr int order(Plus const&, Int... orders)
47
48
49
50
51
    {
      return Math::max(int(orders)...);
    }

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

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

60
    /// Functor that represents A-B
61
62
63
64
65
66
67
68
    struct Minus
    {
      template <class T, class S>
      constexpr auto operator()(T const& lhs, S const& rhs) const
      {
        return lhs - rhs;
      }

69
      friend constexpr int order(Minus const&, int lhs, int rhs)
70
71
72
73
      {
        return Math::max(lhs, rhs);
      }

74
      friend constexpr auto partial(Minus const&, index_t<0>)
75
76
77
78
      {
        return One{};
      }

79
      friend constexpr auto partial(Minus const&, index_t<1>)
80
81
82
83
84
      {
        return StaticConstant<int,-1>{};
      }
    };

85
#ifndef DOXYGEN
86
87
88
    // g - [0] => g
    template <class G>
    struct ComposerBuilder<Minus, G, Zero>
89
      : ComposerBuilder<Id, G> {};
90
91
92
93
94

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

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

99
    /// Functor that represents A*B
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    struct Multiplies
    {
#ifdef AMDIS_HAS_CXX_FOLD_EXPRESSIONS
      template <class... Ts>
      constexpr auto operator()(Ts const&... ts) const
      {
        return (ts * ...);
      }
#else
      template <class T, class S>
      constexpr auto operator()(T const& lhs, S const& rhs) const
      {
        return lhs * rhs;
      }
#endif
    };

117
118
119

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

124
    // g * [0] => [0]
125
126
127
128
129
130
131
132
    template <class G>
    struct ComposerBuilder<Multiplies, G, Zero>
      : ComposerBuilder<Id, Zero> {};

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


    template <class... Int>
137
    constexpr int order(Multiplies const&, Int... orders)
138
    {
139
      return Math::sum(int(orders)...);
140
141
142
143
144
    }

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

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

153
154
155
156
157
158
159
160
161
162
    /// Functor that represents A/B
    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
163
      friend constexpr auto partial(Divides const&, index_t<0>)
164
165
166
167
168
      {
        return compose(Divides{}, One{}, Arg<1>{});
      }

      // d_1 f(x,y) = (y - x)/y^2
169
      friend constexpr auto partial(Divides const&, index_t<1>);
170
171
172
173
174
175
176
177
178
179
180
181
182
    };

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

    /// Functor that represents A-B
    struct Negate
    {
      template <class T>
      constexpr auto operator()(T const& x) const
      {
        return -x;
      }

183
      friend constexpr int order(Negate const&, int d)
184
185
186
187
      {
        return d;
      }

188
      friend constexpr auto partial(Negate const&, index_t<0>)
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
      {
        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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
213
    // forward declaration
214
    template <int p, bool positive>
Praetorius, Simon's avatar
Praetorius, Simon committed
215
216
217
218
219
    struct PowImpl;

    template <int p>
    struct PowType
    {
220
      using type = PowImpl<p, (p>0)>;
Praetorius, Simon's avatar
Praetorius, Simon committed
221
222
223
224
225
226
227
228
229
230
    };

    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>;

231
    /// Functor that represents x^p
232
    template <int p>
233
    struct PowImpl<p, true>
234
235
236
237
238
239
240
    {
      template <class T>
      constexpr auto operator()(T const& x) const
      {
        return Math::pow<p>(x);
      }

241
      friend constexpr int order(PowImpl const&, int d)
242
243
244
245
      {
        return p*d;
      }

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

252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
    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
281
    inline constexpr auto partial(Divides const&, index_t<1>)
282
283
284
285
286
    {
      return compose(Divides{}, compose(Minus{}, Arg<1>{}, Arg<0>{}),
                                compose(Pow<2>{}, Arg<1>{}));
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
287
    /// Functor that represents x^p, \see \ref Pow
288
289
    struct Pow_
    {
290
      constexpr Pow_(int p)
291
        : p_(p)
292
      {}
293
294
295
296
297
298
299

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

300
      friend constexpr int order(Pow_ const& P, int d)
301
302
303
304
      {
        return P.p_ * d;
      }

305
      friend constexpr auto partial(Pow_ const& P, index_t<0>)
306
307
308
309
310
311
312
      {
        return compose(Multiplies{}, Constant<int>{P.p_}, Pow_{P.p_-1});
      }

      int p_;
    };

313
314
315
316
    /** @} **/

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