DUNE-ACFEM (unstable)

fractionconstant.hh
1 #ifndef __DUNE_ACFEM_COMMON_FRACTIONCONSTANT_HH__
2 #define __DUNE_ACFEM_COMMON_FRACTIONCONSTANT_HH__
3 
4 #include <string>
5 #include <numeric>
6 #include <ratio>
7 
8 #include <dune/acfem/common/types.hh>
9 namespace Dune {
10 
11  namespace ACFem {
12 
13  namespace TypedValue {
14 
25  template<class Int, Int N, Int D>
26  struct FractionConstant
27  {
28  using IntegerType = Int;
29  static constexpr IntegerType numerator = N;
30  static constexpr IntegerType denominator = D;
31 
32  static_assert(D > 0 && (std::gcd(N, D) == 1 || std::gcd(N, D) == -1),
33  "GCD should be 1 here ... oops");
34 
35  constexpr FractionConstant()
36  {}
37 
38  template<class Int2, Int2 N2, Int2 D2>
39  constexpr FractionConstant(Constant<Int2, N2>, Constant<Int2, D2>)
40  {
41  static_assert((N == (Int)N2 / std::gcd((Int)N2, (Int)D2))
42  &&
43  (D == (Int)D2 / std::gcd((Int)N2, (Int)D2)),
44  "Assignment from integral constants only possible if their quotient matches the value of the fraction");
45  }
46 
47  template<class Int2, Int2 N2>
48  constexpr FractionConstant(Constant<Int2, N2>)
49  : FractionConstant(Constant<Int, (Int)N2>{}, Constant<Int, (Int)1>{})
50  {}
51 
52  constexpr FractionConstant(const FractionConstant&)
53  {}
54 
55  constexpr FractionConstant(FractionConstant&&)
56  {}
57 
58  constexpr FractionConstant& operator=(const FractionConstant&)
59  {
60  return *this;
61  }
62 
63  template<class T, std::enable_if_t<std::numeric_limits<std::decay_t<T> >::is_specialized, int> = 0>
64  constexpr operator T() const
65  {
66  return static_cast<T>(N) / static_cast<T>(D);
67  }
68 
69  };
70 
72  template<class T, bool RequireDecay = false, class SFINAE = void>
73  struct IsFractionConstant
74  : FalseType
75  {};
76 
77  template<class T>
78  struct IsFractionConstant<T, false, std::enable_if_t<!IsDecay<T>::value> >
79  : IsFractionConstant<std::decay_t<T> >
80  {};
81 
82  template<class I, I N, I D, bool RequireDecay>
83  struct IsFractionConstant<FractionConstant<I, N, D>, RequireDecay>
84  : TrueType
85  {};
86 
87  template<class T>
88  constexpr T sign(T val) {
89  return (T(0) < val) - (val < T(0));
90  }
91 
92  template<class T>
93  constexpr T abs(T val) {
94  return sign(val) * val;
95  }
96 
100  template<class I, I N, I D = 1>
101  using Fraction =
102  FractionConstant<I,
103  sign(D) * N / std::gcd(N, D),
104  abs(D) / std::gcd(N, D)>;
105 
107  template<std::ptrdiff_t N, std::ptrdiff_t D = 1>
108  using IntFraction = Fraction<std::ptrdiff_t, N, D>;
109 
123  template<std::ptrdiff_t N, std::ptrdiff_t D = 1>
124  constexpr auto intFraction(IntConstant<N> = IntConstant<N>{}, IntConstant<D> = IntConstant<D>{})
125  {
126  return IntFraction<N, D>{};
127  }
128 
129  template<std::size_t N, std::size_t D = 1>
130  constexpr auto intFraction(IndexConstant<N>, IndexConstant<D> = IndexConstant<D>{})
131  {
132  return IntFraction<(std::ptrdiff_t)N, (std::ptrdiff_t)D>{};
133  }
134 
135  template<class T>
136  constexpr std::ptrdiff_t numerator(T&& t = T{})
137  {
138  static_assert(IsFractionConstant<T>::value, "T must fulfil IsFractionConstant<T>");
139  return std::decay_t<T>::numerator;
140  }
141 
142  template<class T>
143  constexpr std::ptrdiff_t denominator(T&& t = T{})
144  {
145  static_assert(IsFractionConstant<T>::value, "T must fulfil IsFractionConstant<T>");
146  return std::decay_t<T>::denominator;
147  }
148 
150  template<std::ptrdiff_t S>
151  using Sign = IntFraction<(S > 0 ? 1 : (S < 0 ? -1 : 0))>;
152 
154  template<std::ptrdiff_t S>
155  constexpr auto sign(IntConstant<S> = IntConstant<S>{})
156  {
157  return Sign<S>{};
158  }
159 
161  template<class T>
162  using SignOf = Sign<T{}>;
163 
165  template<class T>
166  constexpr auto signOf(T&&)
167  {
168  return SignOf<T>{};
169  }
170 
172  template<class T, class SFINAE = void>
173  struct IsSign;
174 
176  template<class T>
177  struct IsSign<T, std::enable_if_t<!IsFractionConstant<T>::value> >
178  : FalseType
179  {};
180 
184  template<class T>
185  struct IsSign<T, std::enable_if_t<IsFractionConstant<T>::value> >
186  : BoolConstant<(numerator<T>()*numerator<T>() <= 1 && denominator<T>() == 1)>
187  {};
188 
190  template<class T1, class T2>
191  using FieldPromotionType = std::decay_t<decltype(std::declval<T1>() * std::declval<T2>())>;
192 
193  template<class I, I N, I D>
194  constexpr auto operator-(FractionConstant<I, N, D>)
195  {
196  return Fraction<I, -N, D>{};
197  }
198 
199  template<class I, I N, I D>
200  constexpr auto reciprocal(FractionConstant<I, N, D>)
201  {
202  static_assert(N != 0, "Reciprocal of zero attempted.");
203  return Fraction<I, D, N>{};
204  }
205 
206  template<class I1, I1 N1, I1 D1,
207  class I2, I2 N2, I2 D2>
208  constexpr auto operator*(FractionConstant<I1, N1, D1>, FractionConstant<I2, N2, D2>)
209  {
210  // normalize
211  using T1 = Fraction<I1, N1, D1>;
212  using T2 = Fraction<I2, N2, D2>;
213 
214  // cross normalize to reduce the chance of overflow
215  using C1 = Fraction<FieldPromotionType<I1, I2>, T1::numerator, T2::denominator>;
216  using C2 = Fraction<FieldPromotionType<I1, I2>, T2::numerator, T1::denominator>;
217 
218  // calculate
219  return Fraction<FieldPromotionType<I1, I2>,
220  C1::numerator * C2::numerator,
221  C1::denominator * C2::denominator>{};
222  }
223 
224  template<class I1, I1 N1, I1 D1,
225  class I2, I2 N2, I2 D2>
226  constexpr auto operator+(FractionConstant<I1, N1, D1>, FractionConstant<I2, N2, D2>)
227  {
228  // normalize
229  using Int = FieldPromotionType<I1, I2>;
230  using T1 = Fraction<Int, N1, D1>;
231  using T2 = Fraction<Int, N2, D2>;
232 
233  // calculate
234 
235  // reduce the chance of overflow
236  constexpr Int denomGcd = std::gcd(T1::denominator, T2::denominator);
237  constexpr Int denom1 = T1::denominator/denomGcd;
238  constexpr Int denom2 = T2::denominator/denomGcd;
239 
240  return Fraction<Int,
241  T1::numerator*denom2 + T2::numerator*denom1,
242  denom1*denom2*denomGcd>{};
243  }
244 
245  template<class I1, I1 N1, I1 D1,
246  class I2, I2 N2, I2 D2>
247  constexpr auto operator-(FractionConstant<I1, N1, D1>, FractionConstant<I2, N2, D2>)
248  {
249  // normalize
250  using Int = FieldPromotionType<I1, I2>;
251  using T1 = Fraction<Int, N1, D1>;
252  using T2 = Fraction<Int, N2, D2>;
253 
254  // calculate
255 
256  // reduce the chance of overflow
257  constexpr Int denomGcd = std::gcd(T1::denominator, T2::denominator);
258  constexpr Int denom1 = T1::denominator/denomGcd;
259  constexpr Int denom2 = T2::denominator/denomGcd;
260 
261  return Fraction<Int,
262  T1::numerator*denom2 - T2::numerator*denom1,
263  denom1*denom2*denomGcd>{};
264  }
265 
266  template<class I1, I1 N1, I1 D1,
267  class I2, I2 N2, I2 D2>
268  constexpr auto operator/(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
269  {
270  static_assert(N2 != 0, "Division by Zero.");
271 
272  return f1 * reciprocal(f2);
273  }
274 
275  template<class I1, I1 N1, class I2, I2 N2>
276  constexpr auto operator%(FractionConstant<I1, N1, 1>, FractionConstant<I2, N2, 1>)
277  {
278  return Fraction<FieldPromotionType<I1, I2>, (N1 % N2), 1>{};
279  }
280 
281  template<class I1, I1 N1, class I2, I2 N2>
282  constexpr auto div(FractionConstant<I1, N1, 1>, FractionConstant<I2, N2, 1>)
283  {
284  return Fraction<FieldPromotionType<I1, I2>, N1 / N2, 1>{};
285  }
286 
288  template<class I1, I1 N1, I1 D1, class I2>
289  constexpr auto pow(FractionConstant<I1, N1, D1>, FractionConstant<I2, 0, 1>)
290  {
291  return FractionConstant<I1, 1, 1>{};
292  }
293 
295  template<class I1, I1 N1, I1 D1, class I2>
296  constexpr auto pow(FractionConstant<I1, N1, D1>, FractionConstant<I2, 1, 1>)
297  {
298  return FractionConstant<I1, N1, D1>{};
299  }
300 
302  template<class I1, I1 N1, I1 D1, class I2, I2 N2, std::enable_if_t<(N2 > 1), int> = 0>
303  constexpr auto pow(FractionConstant<I1, N1, D1>, FractionConstant<I2, N2, 1>)
304  {
305  return FractionConstant<I1, N1, D1>{}*pow(FractionConstant<I1, N1, D1>{}, FractionConstant<I2, N2-1, 1>{});
306  }
307 
308  template<class I, I N, I D>
309  constexpr auto sqr(FractionConstant<I, N, D>)
310  {
311  return pow(FractionConstant<I, N, D>{}, FractionConstant<I, 2, 1>{});
312  }
313 
315  template<class I1, I1 N1, I1 D1, class I2, I2 N2, std::enable_if_t<(N2 < 0 && N1 != 0), int> = 0>
316  constexpr auto pow(FractionConstant<I1, N1, D1>, FractionConstant<I2, N2, 1>)
317  {
318  return pow(Fraction<I1, D1, N1>{}, FractionConstant<I2, N2, 1>{});
319  }
320 
321  template<class I, I N, I D>
322  constexpr auto floor(FractionConstant<I, N, D>)
323  {
324  if constexpr (D == 1) {
325  return FractionConstant<I, N, 1>{};
326  } else if constexpr (N < 0) {
327  // -1.5 = -3/2 -> -2
328  // -0.5 = -1/2 -> -1
329  return FractionConstant<I, N/D-1, 1>{};
330  } else {
331  // 1.5 = 3/2 -> 1
332  return FractionConstant<I, N/D, 1>{};
333  }
334  }
335 
336  template<class I, I N, I D>
337  constexpr auto ceil(FractionConstant<I, N, D>)
338  {
339  if constexpr (D == 1) {
340  return FractionConstant<I, N, 1>{};
341  } else {
342  // -1.5 = -3/2 -> (-3+1)/2 = -1
343  // 1.5 = 3/2 -> (3+1)/2 = 2
344  return FractionConstant<I, (N+1)/D, 1>{};
345  }
346  }
347 
349  template<class I, I N, I D>
350  constexpr auto round(FractionConstant<I, N, D>)
351  {
352  // -9/14, (-9+14)/14 = 5/14 = 0
353  return FractionConstant<I, (N < 0 ? (N-D/2)/D : (N-D/2)/D) , 1>{};
354  }
355 
356  template<class I, I N, I D>
357  constexpr auto abs(FractionConstant<I, N, D>)
358  {
359  if constexpr (N < 0) {
360  return FractionConstant<I, -N, D>{};
361  } else {
362  return FractionConstant<I, N, D>{};
363  }
364  }
365 
366  template<class I1, I1 N1, I1 D1,
367  class I2, I2 N2, I2 D2>
368  constexpr bool operator==(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
369  {
370  // closure type
371  using Int = FieldPromotionType<I1, I2>;
372 
373  // normalize
374  using T1 = Fraction<Int, N1, D1>;
375  using T2 = Fraction<Int, N2, D2>;
376 
377  // reduce the chance of overflow
378  constexpr Int denomGcd = std::gcd(T1::denominator, T2::denominator);
379  constexpr Int denom1 = T1::denominator/denomGcd;
380  constexpr Int denom2 = T2::denominator/denomGcd;
381 
382  return T1::numerator*denom2 == T2::numerator*denom1;
383  }
384 
385  template<class I1, I1 N1, I1 D1,
386  class I2, I2 N2, I2 D2>
387  constexpr bool operator!=(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
388  {
389  return !(f1 == f2);
390  }
391 
392  template<class I1, I1 N1, I1 D1,
393  class I2, I2 N2, I2 D2>
394  constexpr bool operator<(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
395  {
396  // closure type
397  using Int = FieldPromotionType<I1, I2>;
398 
399  // normalize
400  using T1 = Fraction<Int, N1, D1>;
401  using T2 = Fraction<Int, N2, D2>;
402 
403  // reduce the chance of overflow
404  constexpr Int denomGcd = std::gcd(T1::denominator, T2::denominator);
405  constexpr Int denom1 = T1::denominator/denomGcd;
406  constexpr Int denom2 = T2::denominator/denomGcd;
407 
408  return T1::numerator*denom2 < T2::numerator*denom1;
409  }
410 
411  template<class I1, I1 N1, I1 D1,
412  class I2, I2 N2, I2 D2>
413  constexpr bool operator<=(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
414  {
415  return f1 < f2 || f1 == f2;
416  }
417 
418  template<class I1, I1 N1, I1 D1,
419  class I2, I2 N2, I2 D2>
420  constexpr bool operator>(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
421  {
422  return f2 < f1;
423  }
424 
425  template<class I1, I1 N1, I1 D1,
426  class I2, I2 N2, I2 D2>
427  constexpr bool operator>=(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
428  {
429  return f2 <= f1;
430  }
431 
432  using std::min;
433 
434  template<class I1, I1 N1, I1 D1,
435  class I2, I2 N2, I2 D2,
436  std::enable_if_t<!std::is_same<FractionConstant<I1, N1, D1>,
437  FractionConstant<I2, N2, D2> >::value, int> = 0>
438  constexpr auto min(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
439  {
440  if constexpr (f1 <= f2) {
441  return f1;
442  } else {
443  return f2;
444  }
445  }
446 
447  using std::max;
448 
449  template<class I1, I1 N1, I1 D1,
450  class I2, I2 N2, I2 D2,
451  std::enable_if_t<!std::is_same<FractionConstant<I1, N1, D1>,
452  FractionConstant<I2, N2, D2> >::value, int> = 0>
453  constexpr auto max(FractionConstant<I1, N1, D1> f1, FractionConstant<I2, N2, D2> f2)
454  {
455  if constexpr (f1 >= f2) {
456  return f1;
457  } else {
458  return f2;
459  }
460  }
461 
462  template<class I, I N, I D>
463  std::string toString(FractionConstant<I, N, D>)
464  {
465  // normalize
466  using T = Fraction<I, N, D>;
467  if constexpr (T::denominator == 1) {
468  return std::to_string(T::numerator) + "_c";
469  } else {
470  return "(" + std::to_string(T::numerator) + "/" + std::to_string(T::denominator) + ")_c";
471  }
472  }
473 
475  template<class I, I N, I D>
476  std::ostream& operator<<(std::ostream& out, FractionConstant<I, N, D> frac)
477  {
478  return out << toString(frac);
479  }
480 
481  } // NS TypedValue
482 
483  template<class T>
484  using IsFractionConstant = TypedValue::IsFractionConstant<T>;
485 
486  } // namespace ACFem
487 
488 } // Namespace Dune
489 
490 #endif // __DUNE_ACFEM_COMMON_FRACTIONCONSTANT_HH__
Dune::ACFem::operator<<
std::ostream & operator<<(std::ostream &out, TypeString< T > &&t)
Output operator for TypePrint tag-structure.
Definition: ostream.hh:39
Dune::ACFem::Tensor::pow
auto pow(T1 &&t1, T2 &&t2)
Power operations with scalar exponents are promoted to component-wise power operations by multiplying...
Definition: expressions.hh:525
Dune::ACFem::Constant
integral_constant< T, V > Constant
Short-cut for any integral constant.
Definition: types.hh:40
Dune::ACFem::BoolConstant
Constant< bool, V > BoolConstant
Short-cut for integral constant of type bool.
Definition: types.hh:48
Dune::ACFem::FalseType
BoolConstant< false > FalseType
Alias for std::false_type.
Definition: types.hh:110
Dune::ACFem::TrueType
BoolConstant< true > TrueType
Alias for std::true_type.
Definition: types.hh:107
Dune::ACFem::TypedValue::FractionConstant
A class implementing compile-time constant fractions of integers.
Definition: fractionconstant.hh:27
Dune::ACFem::TypedValue::IsFractionConstant
Definition: fractionconstant.hh:75
Dune::ACFem::TypedValue::IsSign
Definition: fractionconstant.hh:173
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 6, 22:30, 2024)