1#ifndef _GLUCAT_CLIFFORD_ALGEBRA_IMP_H
2#define _GLUCAT_CLIFFORD_ALGEBRA_IMP_H
64 template<
typename Scalar_T,
typename Index_Set_T,
typename Multivector_T>
68 {
return "clifford_algebra"; }
71 template<
typename Scalar_T,
typename Index_Set_T,
typename Multivector_T>
75 default_truncation = std::numeric_limits<Scalar_T>::epsilon();
80 template<
typename, const index_t, const index_t,
typename>
class Multivector,
81 template<
typename, const index_t, const index_t,
typename>
class RHS,
82 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
86 operator!= (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
bool
87 {
return !(lhs == rhs); }
90 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
91 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
94 operator!= (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const Scalar_T& scr) ->
bool
95 {
return !(lhs == scr); }
98 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
99 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
102 operator!= (
const Scalar_T& scr,
const Multivector<Scalar_T,LO,HI,Tune_P>& rhs) ->
bool
103 {
return !(rhs == scr); }
108 template<
typename, const index_t, const index_t,
typename>
class Multivector,
109 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
114 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
115 static const auto scalar_eps = std::numeric_limits<Scalar_T>::epsilon();
116 static const auto nbr_different_bits =
117 std::numeric_limits<Scalar_T>::digits / Tune_P::denom_different_bits + Tune_P::extra_different_bits;
118 static const auto abs_tol = scalar_eps *
120 using framed_multi_t =
typename multivector_t::framed_multi_t;
121 const auto nbr_terms = double(framed_multi_t(val).truncated(scalar_eps).nbr_terms());
122 return abs_tol * abs_tol * std::max(Scalar_T(nbr_terms), Scalar_T(1));
128 template<
typename, const index_t, const index_t,
typename>
class Multivector,
129 template<
typename, const index_t, const index_t,
typename>
class RHS,
130 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
135 const RHS<Scalar_T,LO,HI,Tune_P>& rhs,
136 const Scalar_T threshold) -> Scalar_T
138 const auto relative =
norm(rhs) > threshold;
139 const auto abs_norm_diff =
norm(rhs-lhs);
141 ? abs_norm_diff/
norm(rhs)
148 template<
typename, const index_t, const index_t,
typename>
class Multivector,
149 template<
typename, const index_t, const index_t,
typename>
class RHS,
150 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
155 const RHS<Scalar_T,LO,HI,Tune_P>& rhs,
156 const Scalar_T threshold,
157 const Scalar_T tolerance) ->
bool
163 template<
typename, const index_t, const index_t,
typename>
class Multivector,
164 template<
typename, const index_t, const index_t,
typename>
class RHS,
165 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
170 const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
bool
177 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
178 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
181 operator+ (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const Scalar_T& scr) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
184 return result += scr;
188 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
189 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
192 operator+ (
const Scalar_T& scr,
const Multivector<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
200 template<
typename, const index_t, const index_t,
typename>
class Multivector,
201 template<
typename, const index_t, const index_t,
typename>
class RHS,
202 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
206 operator+ (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
209 return result += rhs;
213 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
214 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
217 operator- (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const Scalar_T& scr) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
220 return result -= scr;
224 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
225 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
228 operator- (
const Scalar_T& scr,
const Multivector<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
229 {
return -rhs + scr; }
234 template<
typename, const index_t, const index_t,
typename>
class Multivector,
235 template<
typename, const index_t, const index_t,
typename>
class RHS,
236 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
240 operator- (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
243 return result -= rhs;
247 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
248 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
251 operator* (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const Scalar_T& scr) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
254 return result *= scr;
258 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
259 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
262 operator* (
const Scalar_T& scr,
const Multivector<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
271 template<
typename, const index_t, const index_t,
typename>
class Multivector,
272 template<
typename, const index_t, const index_t,
typename>
class RHS,
273 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
277 operator* (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
279 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
280 return lhs * multivector_t(rhs);
286 template<
typename, const index_t, const index_t,
typename>
class Multivector,
287 template<
typename, const index_t, const index_t,
typename>
class RHS,
288 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
292 operator^ (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
294 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
295 return lhs ^ multivector_t(rhs);
301 template<
typename, const index_t, const index_t,
typename>
class Multivector,
302 template<
typename, const index_t, const index_t,
typename>
class RHS,
303 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
307 operator& (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
309 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
310 return lhs & multivector_t(rhs);
316 template<
typename, const index_t, const index_t,
typename>
class Multivector,
317 template<
typename, const index_t, const index_t,
typename>
class RHS,
318 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
322 operator% (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
324 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
325 return lhs % multivector_t(rhs);
331 template<
typename, const index_t, const index_t,
typename>
class Multivector,
332 template<
typename, const index_t, const index_t,
typename>
class RHS,
333 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
337 star (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) -> Scalar_T
339 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
340 return star(lhs, multivector_t(rhs));
344 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
345 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
348 operator/ (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const Scalar_T& scr) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
351 return result /= scr;
355 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
356 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
359 operator/ (
const Scalar_T& scr,
const Multivector<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
361 Multivector<Scalar_T,LO,HI,Tune_P> result = scr;
362 return result /= rhs;
368 template<
typename, const index_t, const index_t,
typename>
class Multivector,
369 template<
typename, const index_t, const index_t,
typename>
class RHS,
370 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
374 operator/ (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
376 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
377 return lhs / multivector_t(rhs);
383 template<
typename, const index_t, const index_t,
typename>
class Multivector,
384 template<
typename, const index_t, const index_t,
typename>
class RHS,
385 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
389 operator| (
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
391 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
392 return lhs | multivector_t(rhs);
396 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
397 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
400 inv(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
401 {
return val.inv(); }
404 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
405 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
407 pow(
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
int rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
409 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
410 if (lhs == Scalar_T(0))
420 auto result = multivector_t(Scalar_T(1));
440 template<
typename, const index_t, const index_t,
typename>
class Multivector,
441 template<
typename, const index_t, const index_t,
typename>
class RHS,
442 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
446 pow(
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
const RHS<Scalar_T,LO,HI,Tune_P>& rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
450 if (lhs == Scalar_T(0))
452 const Scalar_T m = rhs.scalar();
463 return exp(
log(lhs) * rhs);
467 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
468 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
470 outer_pow(
const Multivector<Scalar_T,LO,HI,Tune_P>& lhs,
int rhs) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
471 {
return lhs.outer_pow(rhs); }
474 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
475 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
478 scalar(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) -> Scalar_T
479 {
return val.scalar(); }
482 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
483 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
486 real(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) -> Scalar_T
487 {
return val.scalar(); }
492 template<
typename, const index_t, const index_t,
typename>
class Multivector,
493 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P
497 imag(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) -> Scalar_T
498 {
return Scalar_T(0); }
501 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
502 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
505 pure(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
506 {
return val - val.scalar(); }
509 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
510 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
513 even(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
514 {
return val.even(); }
517 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
518 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
521 odd(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
522 {
return val.odd(); }
525 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
526 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
529 vector_part(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const std::vector<Scalar_T>
530 {
return val.vector_part(); }
533 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
534 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
537 involute(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
538 {
return val.involute(); }
541 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
542 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
545 reverse(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
546 {
return val.reverse(); }
549 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
550 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
553 conj(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
554 {
return val.conj(); }
557 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
558 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
561 quad(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) -> Scalar_T
562 {
return val.quad(); }
565 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
566 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
569 norm(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) -> Scalar_T
570 {
return val.norm(); }
573 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
574 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
577 abs(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) -> Scalar_T
581 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
582 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
585 max_abs(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) -> Scalar_T
586 {
return val.max_abs(); }
589 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
590 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
592 complexifier(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
594 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
597 auto frm = val.frame();
598 using array_t = std::array<index_t, 4>;
599 auto incp = array_t{0, 2, 1, 0};
600 auto incq = array_t{1, 0, 0, 0};
601 auto bott =
pos_mod((frm.count_pos() - frm.count_neg()), 4);
628 auto new_bott =
pos_mod(frm.count_pos() - frm.count_neg(), 4);
630 if ((incp[new_bott] == 0) && (incq[new_bott] == 0))
631 return multivector_t(frm, Scalar_T(1));
634 return traits_t::NaN();
639 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
640 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
643 elliptic(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
647 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
648 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
653 const Multivector<Scalar_T,LO,HI,Tune_P>& i,
const bool prechecked =
false)
657 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
658 using index_set_t =
typename multivector_t::index_set_t;
659 using error_t =
typename multivector_t::error_t;
661 const auto i_frame = i.frame();
664 (val.frame() | i_frame) != i_frame ||
666 throw error_t(
"check_complex(val, i): i is not a valid complexifier for val");
671 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
672 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
675 sqrt(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
676 {
return sqrt(val, i, prechecked); }
679 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
680 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
683 sqrt(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
687 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
688 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
690 clifford_exp(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
698 const auto scalar_val = val.scalar();
699 const auto scalar_exp = traits_t::exp(scalar_val);
700 if (traits_t::isNaN_or_isInf(scalar_exp))
701 return traits_t::NaN();
702 if (val == scalar_val)
705 using multivector_t = Multivector<Scalar_T,LO,HI,Tune_P>;
706 auto A = val - scalar_val;
707 const auto pure_scale2 = A.norm();
709 if (traits_t::isNaN_or_isInf(pure_scale2))
710 return traits_t::NaN();
711 if (pure_scale2 == Scalar_T(0))
714 const auto ilog2_scale =
715 std::max(0, traits_t::to_int(ceil((
log2(pure_scale2) + Scalar_T(A.frame().count()))/Scalar_T(2))) - 3);
716 const auto i_scale = traits_t::pow(Scalar_T(2), ilog2_scale);
717 if (traits_t::isNaN_or_isInf(i_scale))
718 return traits_t::NaN();
721 multivector_t pure_exp;
723 using limits_t = std::numeric_limits<Scalar_T>;
724 const auto nbr_even_powers = 2*(limits_t::digits / 32) + 4;
725 using nbr_t =
decltype(nbr_even_powers);
728 const auto max_power = 2*nbr_even_powers + 1;
729 static std::array<Scalar_T, max_power+1> c;
730 if (c[0] != Scalar_T(1))
734 k =
decltype(max_power)(0);
737 c[k+1] = c[k]*(max_power-k) / ((2*max_power-k)*(k+1));
741 std::array<multivector_t, nbr_even_powers> AA;
743 AA[1] = AA[0] * AA[0];
746 k != nbr_even_powers;
748 AA[k] = AA[k-2] * AA[1];
751 auto residual = multivector_t();
752 auto U = multivector_t(c[0]);
755 k != nbr_even_powers;
758 const auto& term = AA[k]*c[2*k + 2] - residual;
759 const auto& sum = U + term;
760 residual = (sum - U) - term;
763 residual = multivector_t();
764 auto AV = multivector_t(c[1]);
767 k != nbr_even_powers;
770 const auto& term = AA[k]*c[2*k + 3] - residual;
771 const auto& sum = AV + term;
772 residual = (sum - AV) - term;
776 pure_exp = (U+AV) / (U-AV);
779 k =
decltype(ilog2_scale)(0);
782 pure_exp *= pure_exp;
783 return pure_exp * scalar_exp;
787 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
788 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
791 log(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
792 {
return log(val, i, prechecked); }
795 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
796 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
799 log(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
803 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
804 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
807 cosh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
811 return traits_t::NaN();
813 const auto& s = val.scalar();
815 return traits_t::cosh(s);
816 return (
exp(val)+
exp(-val)) / Scalar_T(2);
821 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
822 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
825 acosh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
830 return traits_t::NaN();
832 const auto radical =
sqrt(val*val - Scalar_T(1), i,
true);
833 return (
norm(val + radical) >=
norm(val))
834 ?
log(val + radical, i,
true)
835 : -
log(val - radical, i,
true);
840 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
841 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
844 acosh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
848 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
849 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
851 cos(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
855 return traits_t::NaN();
857 const auto& s = val.scalar();
859 return traits_t::cos(s);
863 static const auto& twopi = Scalar_T(2) * traits_t::pi();
865 (val - s + traits_t::fmod(s, twopi));
866 return (
exp(z)+
exp(-z)) / Scalar_T(2);
870 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
871 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
874 cos(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
879 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
880 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
883 acos(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
887 return traits_t::NaN();
889 const auto& s = val.scalar();
890 if (val == s && traits_t::abs(s) <= Scalar_T(1))
891 return traits_t::acos(s);
894 return i *
acosh(val, i,
true);
899 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
900 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
903 acos(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
907 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
908 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
911 sinh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
915 return traits_t::NaN();
917 const auto& s = val.scalar();
919 return traits_t::sinh(s);
921 return (
exp(val)-
exp(-val)) / Scalar_T(2);
926 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
927 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
930 asinh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
935 return traits_t::NaN();
937 const auto radical =
sqrt(val*val + Scalar_T(1), i,
true);
938 return (
norm(val + radical) >=
norm(val))
939 ?
log( val + radical, i,
true)
940 : -
log(-val + radical, i,
true);
945 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
946 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
949 asinh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
953 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
954 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
956 sin(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
960 return traits_t::NaN();
962 const auto& s = val.scalar();
964 return traits_t::sin(s);
968 static const auto& twopi = Scalar_T(2) * traits_t::pi();
970 (val - s + traits_t::fmod(s, twopi));
971 return i * (
exp(-z)-
exp(z)) / Scalar_T(2);
975 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
976 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
979 sin(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
984 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
985 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
988 asin(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
992 return traits_t::NaN();
994 const auto& s = val.scalar();
995 if (val == s && traits_t::abs(s) <= Scalar_T(1))
996 return traits_t::asin(s);
999 return -i *
asinh(i * val, i,
true);
1004 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
1005 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
1008 asin(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
1012 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
1013 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
1016 tanh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
1020 return traits_t::NaN();
1022 const auto& s = val.scalar();
1024 return traits_t::tanh(s);
1031 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
1032 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
1035 atanh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
1041 : (
norm(val + Scalar_T(1)) >
norm(val - Scalar_T(1)))
1042 ? (
log(val + Scalar_T(1), i,
true) -
log(-val + Scalar_T(1), i,
true)) / Scalar_T(2)
1043 :
log((val + Scalar_T(1)) / (-val + Scalar_T(1)), i,
true) / Scalar_T(2);
1048 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
1049 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
1052 atanh(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
1056 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
1057 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
1060 tan(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
1064 return traits_t::NaN();
1066 const auto& s = val.scalar();
1068 return traits_t::tan(s);
1071 return sin(val, i,
true) /
cos(val, i,
true);
1075 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
1076 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
1079 tan(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
1084 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
1085 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
1088 atan(
const Multivector<Scalar_T,LO,HI,Tune_P>& val,
const Multivector<Scalar_T,LO,HI,Tune_P>& i,
bool prechecked) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
1092 return traits_t::NaN();
1094 const auto& s = val.scalar();
1096 return traits_t::atan(s);
1099 return -i *
atanh(i * val, i,
true);
1104 template<
template<
typename, const index_t, const index_t,
typename>
class Multivector,
1105 typename Scalar_T,
const index_t LO,
const index_t HI,
typename Tune_P >
1108 atan(
const Multivector<Scalar_T,LO,HI,Tune_P>& val) ->
const Multivector<Scalar_T,LO,HI,Tune_P>
clifford_algebra<> declares the operations of a Clifford algebra
static auto classname() -> const std::string
Extra traits which extend numeric limits.
static auto sqrt(const Scalar_T &val) -> Scalar_T
Square root of scalar.
static auto pow(const Scalar_T &val, int n) -> Scalar_T
Integer power.
auto operator|(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Transformation via twisted adjoint action.
auto even(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Even part.
auto operator*(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const Scalar_T &scr) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Product of multivector and scalar.
auto operator&(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inner product.
auto exp(const framed_multi< Scalar_T, LO, HI, Tune_P > &val) -> const framed_multi< Scalar_T, LO, HI, Tune_P >
Exponential of multivector.
auto cosh(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Hyperbolic cosine of multivector.
auto scalar(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar part.
auto tanh(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Hyperbolic tangent of multivector.
auto approx_equal(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs, const Scalar_T threshold, const Scalar_T tolerance) -> bool
Test for approximate equality of multivectors.
auto tan(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Tangent of multivector with specified complexifier.
auto error_squared(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs, const Scalar_T threshold) -> Scalar_T
Relative or absolute error using the quadratic norm.
int index_t
Size of index_t should be enough to represent LO, HI.
auto elliptic(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
auto sinh(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Hyperbolic sine of multivector.
auto imag(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Imaginary part: deprecated (always 0)
static void check_complex(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false)
Check that i is a valid complexifier for val.
auto involute(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Main involution, each {i} is replaced by -{i} in each term, eg. {1}*{2} -> (-{2})*(-{1})
auto operator+(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const Scalar_T &scr) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Geometric sum of multivector and scalar.
auto conj(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Conjugation, rev o invo == invo o rev.
auto operator%(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Left contraction.
auto sin(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Sine of multivector with specified complexifier.
auto pos_mod(LHS_T lhs, RHS_T rhs) -> LHS_T
Modulo function which works reliably for lhs < 0.
auto outer_pow(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, int rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Outer product power of multivector.
auto pow(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, int rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Integer power of multivector.
auto abs(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Absolute value == sqrt(norm)
auto inv(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Geometric multiplicative inverse.
auto asinh(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inverse hyperbolic sine of multivector with specified complexifier.
auto pure(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Pure part.
auto atanh(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inverse hyperbolic tangent of multivector with specified complexifier.
auto atan(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inverse tangent of multivector with specified complexifier.
auto log2(const Scalar_T &x) -> Scalar_T
Log base 2 of scalar.
auto clifford_exp(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Exponential of multivector.
auto cos(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Cosine of multivector with specified complexifier.
auto operator/(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const Scalar_T &scr) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Quotient of multivector and scalar.
auto log(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Natural logarithm of multivector with specified complexifier.
auto asin(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inverse sine of multivector with specified complexifier.
auto real(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Real part: synonym for scalar part.
auto acosh(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inverse hyperbolic cosine of multivector with specified complexifier.
auto norm(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar_T norm == sum of norm of coordinates.
auto error_squared_tol(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Quadratic norm error tolerance relative to a specific multivector.
auto star(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> Scalar_T
Hestenes scalar product.
auto operator^(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Outer product.
auto max_abs(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Maximum of absolute values of components of multivector: multivector infinity norm.
auto sqrt(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Square root of multivector with specified complexifier.
auto reverse(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Reversion, eg. {1}*{2} -> {2}*{1}.
auto complexifier(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Square root of -1 which commutes with all members of the frame of the given multivector.
auto odd(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Odd part.
auto vector_part(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> const std::vector< Scalar_T >
Vector part of multivector, as a vector_t with respect to frame()
auto quad(const Multivector< Scalar_T, LO, HI, Tune_P > &val) -> Scalar_T
Scalar_T quadratic form == (rev(x)*x)(0)
auto operator-(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const Scalar_T &scr) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Geometric difference of multivector and scalar.
auto acos(const Multivector< Scalar_T, LO, HI, Tune_P > &val, const Multivector< Scalar_T, LO, HI, Tune_P > &i, const bool prechecked=false) -> const Multivector< Scalar_T, LO, HI, Tune_P >
Inverse cosine of multivector with specified complexifier.
auto operator!=(const Multivector< Scalar_T, LO, HI, Tune_P > &lhs, const RHS< Scalar_T, LO, HI, Tune_P > &rhs) -> bool
Test for inequality of multivectors.