1 #ifndef _GLUCAT_MATRIX_MULTI_IMP_H 2 #define _GLUCAT_MATRIX_MULTI_IMP_H 39 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS) 40 # pragma GCC diagnostic push 41 # pragma GCC diagnostic ignored "-Wunused-local-typedefs" 43 #include <boost/numeric/ublas/matrix.hpp> 44 #include <boost/numeric/ublas/matrix_expression.hpp> 45 #include <boost/numeric/ublas/matrix_proxy.hpp> 46 #include <boost/numeric/ublas/matrix_sparse.hpp> 47 #include <boost/numeric/ublas/operation.hpp> 48 #include <boost/numeric/ublas/operation_sparse.hpp> 49 #include <boost/numeric/ublas/triangular.hpp> 50 #include <boost/numeric/ublas/lu.hpp> 51 # if defined(_GLUCAT_GCC_IGNORE_UNUSED_LOCAL_TYPEDEFS) 52 # pragma GCC diagnostic pop 66 template<
typename Scalar_T, const index_t LO, const index_t HI >
70 {
return "matrix_multi"; }
80 static const int offset_log2_dim[] = {0, 1, 0, 1, 1, 2, 1, 1};
82 return (p+q)/2 + offset_log2_dim[bott];
87 template<
typename Matrix_Index_T, const index_t LO, const index_t HI >
95 template<
typename Scalar_T, const index_t LO, const index_t HI >
103 template<
typename Scalar_T, const index_t LO, const index_t HI >
104 template<
typename Other_Scalar_T >
111 typedef typename other_matrix_t::const_iterator1 other_const_iterator1;
112 typedef typename other_matrix_t::const_iterator2 other_const_iterator2;
113 for (other_const_iterator1
117 for (other_const_iterator2
118 val_it2 = val_it1.begin();
119 val_it2 != val_it1.end();
125 template<
typename Scalar_T, const index_t LO, const index_t HI >
126 template<
typename Other_Scalar_T >
136 this->
m_matrix.resize(dim, dim,
false);
139 typedef typename other_matrix_t::const_iterator1 other_const_iterator1;
140 typedef typename other_matrix_t::const_iterator2 other_const_iterator2;
141 for (other_const_iterator1
145 for (other_const_iterator2
146 val_it2 = val_it1.begin();
147 val_it2 != val_it1.end();
154 template<
typename Scalar_T, const index_t LO, const index_t HI >
166 template<
typename Scalar_T, const index_t LO, const index_t HI >
172 this->
m_matrix.resize(dim, dim,
false);
174 *
this +=
term_t(ist, crd);
178 template<
typename Scalar_T, const index_t LO, const index_t HI >
183 if (!prechecked && (ist | frm) != frm)
184 throw error_t(
"multivector_t(ist,crd,frm): cannot initialize with value outside of frame");
186 this->
m_matrix.resize(dim, dim,
false);
188 *
this +=
term_t(ist, crd);
192 template<
typename Scalar_T, const index_t LO, const index_t HI >
198 this->
m_matrix.resize(dim, dim,
false);
204 template<
typename Scalar_T, const index_t LO, const index_t HI >
210 template<
typename Scalar_T, const index_t LO, const index_t HI >
217 throw error_t(
"multivector_t(vec,frm): cannot initialize with vector not matching frame");
219 this->
m_matrix.resize(dim, dim,
false);
221 typename vector_t::const_iterator vec_it = vec.begin();
237 template<
typename Scalar_T, const index_t LO, const index_t HI >
243 template<
typename Scalar_T, const index_t LO, const index_t HI >
249 template<
typename Scalar_T, const index_t LO, const index_t HI >
250 template<
typename Other_Scalar_T >
258 *
this = val.template fast_matrix_multi<Scalar_T>(this->
m_frame);
264 this->
m_matrix.resize(dim, dim,
false);
269 val_it = val.begin();
276 template<
typename Scalar_T, const index_t LO, const index_t HI >
277 template<
typename Other_Scalar_T >
285 *
this = val.template fast_matrix_multi<Scalar_T>(our_frame);
292 this->
m_matrix.resize(dim, dim,
false);
297 val_it = val.begin();
304 template<
typename Scalar_T, const index_t LO, const index_t HI >
305 template<
typename Matrix_T >
312 typedef typename Matrix_T::const_iterator1 const_iterator1;
313 typedef typename Matrix_T::const_iterator2 const_iterator2;
315 mtx_it1 = mtx.begin1();
316 mtx_it1 != mtx.end1();
319 mtx_it2 = mtx_it1.begin();
320 mtx_it2 != mtx_it1.end();
326 template<
typename Scalar_T, const index_t LO, const index_t HI >
333 template<
typename Scalar_T, const index_t LO, const index_t HI >
347 template<
typename Scalar_T, const index_t LO, const index_t HI >
358 framed_multi_t framed_lhs;
359 framed_multi_t framed_rhs;
365 our_frame |= framed_lhs.frame() | framed_rhs.frame();
376 template<
typename Scalar_T, const index_t LO, const index_t HI >
396 #if defined(_GLUCAT_USE_DENSE_MATRICES) 397 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
399 typedef typename matrix_t::const_iterator1 const_iterator1;
400 typedef typename matrix_t::const_iterator2 const_iterator2;
404 it1 = lhs_ref.m_matrix.begin1();
405 it1 != lhs_ref.m_matrix.end1();
412 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
414 it1 = rhs_ref.m_matrix.begin1();
415 it1 != rhs_ref.m_matrix.end1();
422 return ublas::norm_inf(lhs_ref.m_matrix - rhs_ref.m_matrix) == 0;
425 const_iterator1 this_it1 = lhs_ref.
m_matrix.begin1();
426 const_iterator1 it1 = rhs_ref.m_matrix.begin1();
428 (this_it1 != lhs_ref.m_matrix.end1()) && (it1 != rhs_ref.m_matrix.end1());
431 if ( this_it1.index1() != it1.index1() )
433 const_iterator2 this_it2 = this_it1.begin();
434 const_iterator2 it2 = it1.begin();
436 (this_it2 != this_it1.end()) && (it2 != it1.end());
438 if ( (this_it2.index2() != it2.index2()) || (*this_it2 != *it2) )
440 if ( (this_it2 != this_it1.end()) || (it2 != it1.end()) )
443 return (this_it1 == lhs_ref.m_matrix.end1()) && (it1 == rhs_ref.m_matrix.end1());
448 template<
typename Scalar_T, const index_t LO, const index_t HI >
454 if (scr != Scalar_T(0))
456 else if (ublas::norm_inf(this->
m_matrix) != 0)
461 return !(dim == 1 && this->
isnan());
466 template<
typename Scalar_T, const index_t LO, const index_t HI >
474 template<
typename Scalar_T, const index_t LO, const index_t HI >
482 return *
this *= Scalar_T(2);
496 template<
typename Scalar_T, const index_t LO, const index_t HI >
504 return *
this = Scalar_T(0);
518 template<
typename Scalar_T, const index_t LO, const index_t HI >
526 template<
typename Scalar_T, const index_t LO, const index_t HI >
534 if (traits_t::isNaN_or_isInf(scr) || this->
isnan())
535 return *
this = traits_t::NaN();
536 if (scr == Scalar_T(0))
544 template<
typename Scalar_T, const index_t LO, const index_t HI >
552 #if defined(_GLUCAT_CHECK_ISNAN) 558 multivector_t lhs_reframed;
559 multivector_t rhs_reframed;
560 const index_set_t our_frame =
reframe(lhs, rhs, lhs_reframed, rhs_reframed);
561 const multivector_t& lhs_ref = (lhs.
m_frame == our_frame)
564 const multivector_t& rhs_ref = (rhs.
m_frame == our_frame)
569 #if defined(_GLUCAT_USE_DENSE_MATRICES) 572 const matrix_index_t dim = lhs_ref.m_matrix.size1();
574 result.m_matrix.clear();
575 ublas::axpy_prod(lhs_ref.m_matrix, rhs_ref.m_matrix, result.m_matrix,
true);
578 typedef typename matrix_t::expression_type expression_t;
581 multivector_t(ublas::sparse_prod<expression_t>(lhs_ref.m_matrix, rhs_ref.m_matrix), our_frame);
586 template<
typename Scalar_T, const index_t LO, const index_t HI >
591 {
return *
this = *
this * rhs; }
594 template<
typename Scalar_T, const index_t LO, const index_t HI >
605 template<
typename Scalar_T, const index_t LO, const index_t HI >
610 {
return *
this = *
this ^ rhs; }
613 template<
typename Scalar_T, const index_t LO, const index_t HI >
624 template<
typename Scalar_T, const index_t LO, const index_t HI >
629 {
return *
this = *
this & rhs; }
632 template<
typename Scalar_T, const index_t LO, const index_t HI >
643 template<
typename Scalar_T, const index_t LO, const index_t HI >
648 {
return *
this = *
this % rhs; }
651 template<
typename Scalar_T, const index_t LO, const index_t HI >
655 {
return (lhs * rhs).scalar(); }
658 template<
typename Scalar_T, const index_t LO, const index_t HI >
663 {
return *
this *= Scalar_T(1)/scr; }
666 template<
typename Scalar_T, const index_t LO, const index_t HI >
672 #if defined(_GLUCAT_CHECK_ISNAN) 674 return traits_t::NaN();
677 if (rhs == Scalar_T(0))
678 return traits_t::NaN();
684 multivector_t lhs_reframed;
685 multivector_t rhs_reframed;
686 const index_set_t our_frame =
reframe(lhs, rhs, lhs_reframed, rhs_reframed);
687 const multivector_t& lhs_ref = (lhs.
m_frame == our_frame)
690 const multivector_t& rhs_ref = (rhs.
m_frame == our_frame)
703 const matrix_t& AT = ublas::trans(rhs_ref.m_matrix);
706 typedef ublas::permutation_matrix<matrix_index_t> permutation_t;
708 permutation_t pvector(AT.size1());
709 if (! ublas::lu_factorize(LU, pvector))
711 const matrix_t& BT = ublas::trans(lhs_ref.m_matrix);
713 ublas::lu_substitute(LU, pvector, XT);
714 #if defined(_GLUCAT_CHECK_ISNAN) 716 return traits_t::NaN();
726 ublas::axpy_prod(AT, XT, R,
false);
727 #if defined(_GLUCAT_CHECK_ISNAN) 729 return traits_t::NaN();
732 Scalar_T nr = ublas::norm_inf(R);
733 if ( nr != Scalar_T(0) && !traits_t::isNaN_or_isInf(nr) )
736 Scalar_T nrold = nr + Scalar_T(1);
749 ublas::lu_substitute(LU, pvector, D);
753 ublas::axpy_prod(AT, XTnew, R,
false);
754 nr = ublas::norm_inf(R);
762 return traits_t::NaN();
766 template<
typename Scalar_T, const index_t LO, const index_t HI >
771 {
return *
this = *
this / rhs; }
774 template<
typename Scalar_T, const index_t LO, const index_t HI >
778 {
return rhs * lhs / rhs.
involute(); }
781 template<
typename Scalar_T, const index_t LO, const index_t HI >
786 {
return *
this = rhs * *
this / rhs.
involute(); }
789 template<
typename Scalar_T, const index_t LO, const index_t HI >
797 template<
typename Scalar_T, const index_t LO, const index_t HI >
805 template<
typename Scalar_T, const index_t LO, const index_t HI >
811 throw error_t(
"outer_pow(m): negative exponent");
826 template<
typename Scalar_T, const index_t LO, const index_t HI >
834 template<
typename Scalar_T, const index_t LO, const index_t HI >
842 template<
typename Scalar_T, const index_t LO, const index_t HI >
856 template<
typename Scalar_T, const index_t LO, const index_t HI >
862 if ((grade < 0) || (grade > HI-LO))
869 template<
typename Scalar_T, const index_t LO, const index_t HI >
880 template<
typename Scalar_T, const index_t LO, const index_t HI >
885 {
return *
this - this->
scalar(); }
888 template<
typename Scalar_T, const index_t LO, const index_t HI >
896 template<
typename Scalar_T, const index_t LO, const index_t HI >
904 template<
typename Scalar_T, const index_t LO, const index_t HI >
911 template<
typename Scalar_T, const index_t LO, const index_t HI >
916 if (!prechecked && (this->
frame() | frm) != frm)
917 throw error_t(
"vector_part(frm): value is outside of requested frame");
920 if (this->
frame() != frm)
923 const index_t begin_index = frm.min();
924 const index_t end_index = frm.max()+1;
939 template<
typename Scalar_T, const index_t LO, const index_t HI >
947 template<
typename Scalar_T, const index_t LO, const index_t HI >
955 template<
typename Scalar_T, const index_t LO, const index_t HI >
963 template<
typename Scalar_T, const index_t LO, const index_t HI >
974 template<
typename Scalar_T, const index_t LO, const index_t HI >
985 template<
typename Scalar_T, const index_t LO, const index_t HI >
993 template<
typename Scalar_T, const index_t LO, const index_t HI >
1002 template<
typename Scalar_T, const index_t LO, const index_t HI >
1006 write(
const std::string& msg)
const 1010 template<
typename Scalar_T, const index_t LO, const index_t HI >
1014 write(std::ofstream& ofile,
const std::string& msg)
const 1017 throw error_t(
"write(ofile,msg): cannot write to output file");
1022 template<
typename Scalar_T, const index_t LO, const index_t HI >
1025 operator<< (std::ostream& os, const matrix_multi<Scalar_T,LO,HI>& val)
1027 os << typename matrix_multi<Scalar_T,LO,HI>::framed_multi_t(val);
1032 template<
typename Scalar_T, const index_t LO, const index_t HI >
1047 template<
typename Scalar_T, const index_t LO, const index_t HI >
1053 if (std::numeric_limits<Scalar_T>::has_quiet_NaN)
1060 template<
typename Scalar_T, const index_t LO, const index_t HI >
1068 template<
typename Scalar_T, const index_t LO, const index_t HI >
1074 if (term.second != Scalar_T(0))
1080 template<
typename Multivector_T,
typename Matrix_T,
typename Basis_Matrix_T >
1091 typedef typename basis_matrix_t::value_type basis_scalar_t;
1097 if (ublas::norm_inf(X) == 0)
1100 const basis_matrix_t& I = matrix::unit<basis_matrix_t>(2);
1101 basis_matrix_t J(2,2,2);
1103 J(0,1) = basis_scalar_t(-1);
1104 J(1,0) = basis_scalar_t( 1);
1105 basis_matrix_t K = J;
1106 K(0,1) = basis_scalar_t( 1);
1107 basis_matrix_t JK = I;
1108 JK(0,0) = basis_scalar_t(-1);
1113 const index_set_t ist_mnpn = ist_mn | ist_pn;
1123 result +=
term_t(ist_mn, j_x);
1124 result +=
term_t(ist_pn, k_x);
1125 return result +=
term_t(ist_mnpn, jk_x);
1132 const framed_multi_t& i_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1134 const framed_multi_t& j_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1136 const framed_multi_t& k_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1138 const framed_multi_t& jk_x = fast<framed_multi_t, matrix_t, basis_matrix_t>
1141 result = i_x.even() - jk_x.odd();
1142 result += (j_x.even() - k_x.odd()) * mn;
1143 result += (k_x.even() - j_x.odd()) * pn;
1144 return result += (jk_x.even() - i_x.odd()) * mnpn;
1149 template<
typename Scalar_T, const index_t LO, const index_t HI >
1158 return (this->
template fast_framed_multi<Scalar_T>()).template fast_matrix_multi<Scalar_T>(frm);
1162 template<
typename Scalar_T, const index_t LO, const index_t HI >
1163 template <
typename Other_Scalar_T>
1188 const index_t level = (p+q)/2;
1192 framed_multi_t val = fast<framed_multi_t, matrix_t, basis_matrix_t>(this->
m_matrix, level);
1195 switch (
pos_mod(orig_p-orig_q, 8))
1205 if (orig_p-orig_q > 4)
1208 if (orig_p-orig_q < -3)
1217 template<
typename Scalar_T, const index_t LO, const index_t HI,
typename Matrix_T >
1219 public std::map< std::pair< const index_set<LO,HI>, const index_set<LO,HI> >,
1236 friend class friend_for_private_destructor;
1240 template<
typename Scalar_T, const index_t LO, const index_t HI >
1245 typedef std::pair<const index_set_t, const index_set_t> index_set_pair_t;
1246 const index_set_pair_t& unfolded_pair = index_set_pair_t(ist, this->
m_frame);
1249 typedef typename basis_table_t::const_iterator basis_const_iterator_t;
1250 basis_table_t& basis_cache = basis_table_t::basis();
1257 const basis_const_iterator_t basis_it = basis_cache.find(unfolded_pair);
1258 if (basis_it != basis_cache.end())
1259 return *(basis_it->second);
1263 const index_set_pair_t& folded_pair = index_set_pair_t(folded_set, folded_frame);
1264 typedef std::pair<const index_set_pair_t, basis_matrix_t*> basis_pair_t;
1267 const basis_const_iterator_t basis_it = basis_cache.find(folded_pair);
1268 if (basis_it != basis_cache.end())
1271 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1275 const index_t folded_max = folded_frame.
max();
1276 const index_t folded_min = folded_frame.
min();
1291 basis_cache.insert(basis_pair_t(folded_pair, result_ptr));
1292 basis_cache.insert(basis_pair_t(unfolded_pair, result_ptr));
1298 template<
typename Scalar_T, const index_t LO, const index_t HI >
1312 return traits_t::NaN();
1315 const int nbr_even_powers = array_size/2 - 1;
1318 std::vector<multivector_t> XX(nbr_even_powers);
1320 XX[1] = XX[0] * XX[0];
1323 k != nbr_even_powers;
1325 XX[k] = XX[k-2] * XX[1];
1328 multivector_t N = a[1];
1331 k != nbr_even_powers;
1333 N += XX[k] * a[2*k + 3];
1338 k != nbr_even_powers;
1340 N += XX[k] * a[2*k + 2];
1341 multivector_t D = b[1];
1344 k != nbr_even_powers;
1346 D += XX[k] * b[2*k + 3];
1351 k != nbr_even_powers;
1353 D += XX[k] * b[2*k + 2];
1358 template<
typename Scalar_T, const index_t LO, const index_t HI >
1366 const multivector_t& invM =
inv(M);
1367 M = ((M + invM)/Scalar_T(2) + Scalar_T(1)) / Scalar_T(2);
1368 Y *= (invM + Scalar_T(1)) / Scalar_T(2);
1372 template<
typename Scalar_T, const index_t LO, const index_t HI >
1380 if (val == Scalar_T(0))
1383 typedef std::numeric_limits<Scalar_T> limits_t;
1384 static const Scalar_T tol =
std::pow(limits_t::epsilon(), 2);
1385 static const Scalar_T tol2 = tol * tol;
1387 multivector_t M = val;
1388 multivector_t Y = val;
1389 Scalar_T norm_M_1 =
norm(M - Scalar_T(1));
1393 step != sqrt_max_steps && norm_M_1 > tol2;
1397 return traits_t::NaN();
1399 norm_M_1 =
norm(M - Scalar_T(1));
1401 if (norm_M_1 > tol2)
1402 return traits_t::NaN();
1411 template<
typename Scalar_T >
1414 static const int array_size = 14;
1415 static const Scalar_T array[array_size];
1420 template<
typename Scalar_T >
1423 static const int array_size = 14;
1424 static const Scalar_T array[array_size];
1427 template<
typename Scalar_T >
1428 const Scalar_T pade_sqrt_a<Scalar_T>::array[pade_sqrt_a<Scalar_T>::array_size] =
1430 1.0, 27.0/4.0, 81.0/4.0, 2277.0/64.0,
1431 10395.0/256.0, 32319.0/1024.0, 8721.0/512.0, 26163.0/4096.0,
1432 53703.0/32768.0, 36465.0/131072.0, 3861.0/131072.0, 7371.0/4194304.0,
1433 819.0/16777216.0, 27.0/67108864.0
1435 template<
typename Scalar_T >
1436 const Scalar_T pade_sqrt_b<Scalar_T>::array[pade_sqrt_b<Scalar_T>::array_size] =
1438 1.0, 25.0/4.0, 69.0/4.0, 1771.0/64.0,
1439 7315.0/256.0, 20349.0/1024.0, 4845.0/512.0, 12597.0/4096.0,
1440 21879.0/32768.0, 12155.0/131072.0, 1001.0/131072.0, 1365.0/4194304.0,
1441 91.0/16777216.0, 1.0/67108864.0
1445 struct pade_sqrt_a<float>
1447 static const int array_size = 10;
1448 static const float array[array_size];
1451 struct pade_sqrt_b<float>
1453 static const int array_size = 10;
1454 static const float array[array_size];
1456 const float pade_sqrt_a<float>::array[pade_sqrt_a<float>::array_size] =
1458 1.0, 19.0/4.0, 19.0/2.0, 665.0/64.0,
1459 1729.0/256.0, 2717.0/1024.0, 627.0/1024.0, 627.0/8192.0,
1460 285.0/65536.0, 19.0/262144.0
1462 const float pade_sqrt_b<float>::array[pade_sqrt_a<float>::array_size] =
1464 1.0, 17.0/4.0, 15.0/2.0, 455.0/64.0,
1465 1001.0/256.0, 1287.0/1024.0, 231.0/1024.0, 165.0/8192.0,
1466 45.0/65536, 1.0/262144.0
1470 struct pade_sqrt_a<long double>
1472 static const int array_size = 18;
1473 static const long double array[array_size];
1476 struct pade_sqrt_b<long double>
1478 static const int array_size = 18;
1479 static const long double array[array_size];
1481 const long double pade_sqrt_a<long double>::array[pade_sqrt_a<long double>::array_size] =
1483 1.0L, 35.0L/4.0L, 35.0L, 5425.0L/64.0L,
1484 35525.0L/256.0L, 166257.0L/1024.0L, 143325.0L/1024.0L, 740025.0L/8192.0L,
1485 2877875.0L/65536.0L, 4206125.0L/262144.0L, 572033.0L/131072.0L, 1820105.0L/2097152.0L,
1486 1028755.0L/8388608.0L, 395675.0L/33554432.0L, 24225.0L/33554432.0L, 6783.0L/268435456.0L,
1487 1785.0L/4294967296.0L, 35.0L/17179869184.0L
1489 const long double pade_sqrt_b<long double>::array[pade_sqrt_a<long double>::array_size] =
1491 1.0L, 33.0L/4.0L, 31.0L, 4495.0L/64.0L,
1492 27405.0L/256.0L, 118755.0L/1024.0L, 94185.0L/1024.0L, 444015.0L/8192.0L,
1493 1562275.0L/65536.0L, 2042975.0L/262144.0L, 245157.0L/131072.0L, 676039.0L/2097152.0L,
1494 323323.0L/8388608.0L, 101745.0L/33554432.0L, 4845.0L/33554432.0L, 969.0L/268435456.0L,
1495 153.0L/4294967296.0L, 1.0L/17179869184.0L
1498 #if defined(_GLUCAT_USE_QD) 1500 struct pade_sqrt_a<dd_real>
1502 static const int array_size = 22;
1503 static const dd_real array[array_size];
1506 struct pade_sqrt_b<dd_real>
1508 static const int array_size = 22;
1509 static const dd_real array[array_size];
1511 const dd_real pade_sqrt_a<dd_real>::array[pade_sqrt_a<dd_real>::array_size] =
1513 dd_real(
"1"), dd_real(
"43")/dd_real(
"4"),
1514 dd_real(
"215")/dd_real(
"4"), dd_real(
"10621")/dd_real(
"64"),
1515 dd_real(
"90687")/dd_real(
"256"), dd_real(
"567987")/dd_real(
"1024"),
1516 dd_real(
"168861")/dd_real(
"256"), dd_real(
"1246355")/dd_real(
"2048"),
1517 dd_real(
"7228859")/dd_real(
"16384"), dd_real(
"16583853")/dd_real(
"65536"),
1518 dd_real(
"7538115")/dd_real(
"65536"), dd_real(
"173376645")/dd_real(
"4194304"),
1519 dd_real(
"195747825")/dd_real(
"16777216"), dd_real(
"171655785")/dd_real(
"67108864"),
1520 dd_real(
"14375115")/dd_real(
"33554432"), dd_real(
"14375115")/dd_real(
"268435456"),
1521 dd_real(
"20764055")/dd_real(
"4294967296"), dd_real(
"5167525")/dd_real(
"17179869184"),
1522 dd_real(
"206701")/dd_real(
"17179869184"), dd_real(
"76153")/dd_real(
"274877906944"),
1523 dd_real(
"3311")/dd_real(
"1099511627776") , dd_real(
"43")/dd_real(
"4398046511104")
1525 const dd_real pade_sqrt_b<dd_real>::array[pade_sqrt_a<dd_real>::array_size] =
1527 dd_real(
"1"), dd_real(
"41")/dd_real(
"4"),
1528 dd_real(
"195")/dd_real(
"4"), dd_real(
"9139")/dd_real(
"64"),
1529 dd_real(
"73815")/dd_real(
"256"), dd_real(
"435897")/dd_real(
"1024"),
1530 dd_real(
"121737")/dd_real(
"256"), dd_real(
"840565")/dd_real(
"2048"),
1531 dd_real(
"4539051")/dd_real(
"16384"), dd_real(
"9641775")/dd_real(
"65536"),
1532 dd_real(
"4032015")/dd_real(
"65536"), dd_real(
"84672315")/dd_real(
"4194304"),
1533 dd_real(
"86493225")/dd_real(
"16777216"), dd_real(
"67863915")/dd_real(
"67108864"),
1534 dd_real(
"5014575")/dd_real(
"33554432"), dd_real(
"4345965")/dd_real(
"268435456"),
1535 dd_real(
"5311735")/dd_real(
"4294967296"), dd_real(
"1081575")/dd_real(
"17179869184"),
1536 dd_real(
"33649")/dd_real(
"17179869184"), dd_real(
"8855")/dd_real(
"274877906944"),
1537 dd_real(
"231")/dd_real(
"1099511627776"), dd_real(
"1")/dd_real(
"4398046511104")
1541 struct pade_sqrt_a<qd_real>
1543 static const int array_size = 34;
1544 static const qd_real array[array_size];
1547 struct pade_sqrt_b<qd_real>
1549 static const int array_size = 34;
1550 static const qd_real array[array_size];
1552 const qd_real pade_sqrt_a<qd_real>::array[pade_sqrt_a<qd_real>::array_size] =
1554 qd_real(
"1"), qd_real(
"67")/qd_real(
"4"),
1555 qd_real(
"134"), qd_real(
"43617")/qd_real(
"64"),
1556 qd_real(
"633485")/qd_real(
"256"), qd_real(
"6992857")/qd_real(
"1024"),
1557 qd_real(
"15246721")/qd_real(
"1024"), qd_real(
"215632197")/qd_real(
"8192"),
1558 qd_real(
"2518145487")/qd_real(
"65536"), qd_real(
"12301285425")/qd_real(
"262144"),
1559 qd_real(
"6344873535")/qd_real(
"131072"), qd_real(
"89075432355")/qd_real(
"2097152"),
1560 qd_real(
"267226297065")/qd_real(
"8388608"), qd_real(
"687479618945")/qd_real(
"33554432"),
1561 qd_real(
"379874182975")/qd_real(
"33554432"), qd_real(
"1443521895305")/qd_real(
"268435456"),
1562 qd_real(
"9425348845815")/qd_real(
"4294967296"), qd_real(
"13195488384141")/qd_real(
"17179869184"),
1563 qd_real(
"987417498133")/qd_real(
"4294967296"), qd_real(
"8055248011085")/qd_real(
"137438953472"),
1564 qd_real(
"6958363175533")/qd_real(
"549755813888"), qd_real(
"5056698705201")/qd_real(
"2199023255552"),
1565 qd_real(
"766166470485")/qd_real(
"2199023255552"), qd_real(
"766166470485")/qd_real(
"17592186044416"),
1566 qd_real(
"623623871325")/qd_real(
"140737488355328"), qd_real(
"203123203803")/qd_real(
"562949953421312"),
1567 qd_real(
"6478601247")/qd_real(
"281474976710656"), qd_real(
"5038912081")/qd_real(
"4503599627370496"),
1568 qd_real(
"719844583")/qd_real(
"18014398509481984"), qd_real(
"71853815")/qd_real(
"72057594037927936"),
1569 qd_real(
"1165197")/qd_real(
"72057594037927936"), qd_real(
"87703")/qd_real(
"576460752303423488"),
1570 qd_real(
"12529")/qd_real(
"18446744073709551616"), qd_real(
"67")/qd_real(
"73786976294838206464")
1572 const qd_real pade_sqrt_b<qd_real>::array[pade_sqrt_a<qd_real>::array_size] =
1574 qd_real(
"1"), qd_real(
"65")/qd_real(
"4"),
1575 qd_real(
"126"), qd_real(
"39711")/qd_real(
"64"),
1576 qd_real(
"557845")/qd_real(
"256"), qd_real(
"5949147")/qd_real(
"1024"),
1577 qd_real(
"12515965")/qd_real(
"1024"), qd_real(
"170574723")/qd_real(
"8192"),
1578 qd_real(
"1916797311")/qd_real(
"65536"), qd_real(
"8996462475")/qd_real(
"262144"),
1579 qd_real(
"4450881435")/qd_real(
"131072"), qd_real(
"59826782925")/qd_real(
"2097152"),
1580 qd_real(
"171503444385")/qd_real(
"8388608"), qd_real(
"420696483235")/qd_real(
"33554432"),
1581 qd_real(
"221120793075")/qd_real(
"33554432"), qd_real(
"797168807855")/qd_real(
"268435456"),
1582 qd_real(
"4923689695575")/qd_real(
"4294967296"), qd_real(
"6499270398159")/qd_real(
"17179869184"),
1583 qd_real(
"456864812569")/qd_real(
"4294967296"), qd_real(
"3486599885395")/qd_real(
"137438953472"),
1584 qd_real(
"2804116503573")/qd_real(
"549755813888"), qd_real(
"1886827875075")/qd_real(
"2199023255552"),
1585 qd_real(
"263012370465")/qd_real(
"2199023255552"), qd_real(
"240141729555")/qd_real(
"17592186044416"),
1586 qd_real(
"176848560525")/qd_real(
"140737488355328"), qd_real(
"51538723353")/qd_real(
"562949953421312"),
1587 qd_real(
"1450433115")/qd_real(
"281474976710656"), qd_real(
"977699359")/qd_real(
"4503599627370496"),
1588 qd_real(
"118183439")/qd_real(
"18014398509481984"), qd_real(
"9652005")/qd_real(
"72057594037927936"),
1589 qd_real(
"121737")/qd_real(
"72057594037927936"), qd_real(
"6545")/qd_real(
"576460752303423488"),
1590 qd_real(
"561")/qd_real(
"18446744073709551616"), qd_real(
"1")/qd_real(
"73786976294838206464")
1598 template<
typename Scalar_T, const index_t LO, const index_t HI >
1609 return traits_t::NaN();
1617 typedef typename traits_t::demoted::type demoted_scalar_t;
1620 const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
1621 const demoted_multivector_t& demoted_i = demoted_multivector_t(i);
1628 typedef typename traits_t::promoted::type promoted_scalar_t;
1631 const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
1632 const promoted_multivector_t& promoted_i = promoted_multivector_t(i);
1643 template<
typename Scalar_T, const index_t LO, const index_t HI >
1654 return traits_t::NaN();
1658 const Scalar_T realval = val.
scalar();
1661 if (realval < Scalar_T(0))
1669 #if !defined(_GLUCAT_USE_EIGENVALUES) 1670 const multivector_t val2 = val*val;
1671 const Scalar_T real_val2 = val2.
scalar();
1672 if (val2 == real_val2 && real_val2 > Scalar_T(0))
1673 return matrix_sqrt(-i * val, i) * (i + Scalar_T(1)) / sqrt_2;
1677 const Scalar_T scale =
1678 (realval != Scalar_T(0) &&
norm(val/realval - Scalar_T(1)) < Scalar_T(1))
1680 : (realval < Scalar_T(0))
1684 if (traits_t::isNaN_or_isInf(sqrt_scale))
1685 return traits_t::NaN();
1688 multivector_t rescale = sqrt_scale;
1689 if (scale < Scalar_T(0))
1690 rescale = i * sqrt_scale;
1692 const multivector_t& unitval = val / scale;
1693 const Scalar_T max_norm = Scalar_T(1.0/4.0);
1695 #if defined(_GLUCAT_USE_EIGENVALUES) 1696 multivector_t scaled_result;
1704 scaled_result =
matrix_sqrt(-i * unitval, i) * (i + Scalar_T(1)) / sqrt_2;
1709 scaled_result =
matrix_sqrt(
exp(i*safe_arg) * unitval, i) *
exp(-i*safe_arg/Scalar_T(2));
1714 (
norm(unitval - Scalar_T(1)) < max_norm)
1717 pade_sqrt_a<Scalar_T>::array,
1718 pade_sqrt_b<Scalar_T>::array,
1719 unitval - Scalar_T(1))
1724 if (scaled_result.isnan())
1725 return traits_t::NaN();
1727 return scaled_result * rescale;
1729 const multivector_t& scaled_result =
1730 (
norm(unitval - Scalar_T(1)) < max_norm)
1733 pade_sqrt_a<Scalar_T>::array,
1734 pade_sqrt_b<Scalar_T>::array,
1735 unitval - Scalar_T(1))
1738 if (scaled_result.isnan())
1741 return traits_t::NaN();
1743 const multivector_t& mi_unitval = -i * unitval;
1745 const multivector_t& scaled_mi_result =
1746 (
norm(mi_unitval - Scalar_T(1)) < max_norm)
1749 pade_sqrt_a<Scalar_T>::array,
1750 pade_sqrt_b<Scalar_T>::array,
1751 mi_unitval - Scalar_T(1))
1754 if (scaled_mi_result.isnan())
1755 return traits_t::NaN();
1757 return scaled_mi_result * rescale * (i + Scalar_T(1)) / sqrt_2;
1760 return scaled_result * rescale;
1768 template<
typename Scalar_T >
1771 static const int array_size = 14;
1772 static const Scalar_T array[array_size];
1777 template<
typename Scalar_T >
1780 static const int array_size = 14;
1781 static const Scalar_T array[array_size];
1783 template<
typename Scalar_T >
1784 const Scalar_T pade_log_a<Scalar_T>::array[pade_log_a<Scalar_T>::array_size] =
1786 0.0, 1.0, 6.0, 4741.0/300.0,
1787 1441.0/60.0, 107091.0/4600.0, 8638.0/575.0, 263111.0/40250.0,
1788 153081.0/80500.0, 395243.0/1101240.0, 28549.0/688275.0, 605453.0/228813200.0,
1789 785633.0/10296594000.0, 1145993.0/1873980108000.0
1791 template<
typename Scalar_T >
1792 const Scalar_T pade_log_b<Scalar_T>::array[pade_log_b<Scalar_T>::array_size] =
1794 1.0, 13.0/2.0, 468.0/25.0, 1573.0/50.0,
1795 1573.0/46.0, 11583.0/460.0, 10296.0/805.0, 2574.0/575.0,
1796 11583.0/10925.0, 143.0/874.0, 572.0/37145.0, 117.0/148580.0,
1797 13.0/742900.0, 1.0/10400600.0
1801 struct pade_log_a<float>
1803 static const int array_size = 10;
1804 static const float array[array_size];
1807 struct pade_log_b<float>
1809 static const int array_size = 10;
1810 static const float array[array_size];
1812 const float pade_log_a<float>::array[pade_log_a<float>::array_size] =
1814 0.0, 1.0, 4.0, 1337.0/204.0,
1815 385.0/68.0, 1879.0/680.0, 193.0/255.0, 197.0/1820.0,
1816 419.0/61880.0, 7129.0/61261200.0
1818 const float pade_log_b<float>::array[pade_log_a<float>::array_size] =
1820 1.0, 9.0/2.0, 144.0/17.0, 147.0/17.0,
1821 441.0/85.0, 63.0/34.0, 84.0/221.0, 9.0/221.0,
1822 9.0/4862.0, 1.0/48620.0
1826 struct pade_log_a<long double>
1828 static const int array_size = 18;
1829 static const long double array[array_size];
1832 struct pade_log_b<long double>
1834 static const int array_size = 18;
1835 static const long double array[array_size];
1837 const long double pade_log_a<long double>::array[pade_log_a<long double>::array_size] =
1839 0.0L, 1.0L, 8.0L, 3835.0L/132.0L,
1840 8365.0L/132.0L, 11363807.0L/122760.0L, 162981.0L/1705.0L, 9036157.0L/125860.0L,
1841 18009875.0L/453096.0L, 44211925.0L/2718576.0L, 4149566.0L/849555.0L, 16973929.0L/16020180.0L,
1842 172459.0L/1068012.0L, 116317061.0L/7025382936.0L, 19679783.0L/18441630207.0L, 23763863.0L/614721006900.0L,
1843 50747.0L/79318839600.0L, 42142223.0L/14295951736466400.0L
1845 const long double pade_log_b<long double>::array[pade_log_a<long double>::array_size] =
1847 1.0L, 17.0L/2.0L, 1088.0L/33.0L, 850.0L/11.0L,
1848 41650.0L/341.0L, 140777.0L/1023.0L, 1126216.0L/9889.0L, 63206.0L/899.0L,
1849 790075.0L/24273.0L, 60775.0L/5394.0L, 38896.0L/13485.0L, 21658.0L/40455.0L,
1850 21658.0L/310155.0L, 4165.0L/682341.0L, 680.0L/2047023.0L, 34.0L/3411705.0L,
1851 17.0L/129644790.0L, 1.0L/2333606220
1853 #if defined(_GLUCAT_USE_QD) 1855 struct pade_log_a<dd_real>
1857 static const int array_size = 22;
1858 static const dd_real array[array_size];
1861 struct pade_log_b<dd_real>
1863 static const int array_size = 22;
1864 static const dd_real array[array_size];
1866 const dd_real pade_log_a<dd_real>::array[pade_log_a<dd_real>::array_size] =
1868 dd_real(
"0"), dd_real(
"1"),
1869 dd_real(
"10"), dd_real(
"22781")/dd_real(
"492"),
1870 dd_real(
"21603")/dd_real(
"164"), dd_real(
"5492649")/dd_real(
"21320"),
1871 dd_real(
"978724")/dd_real(
"2665"), dd_real(
"4191605")/dd_real(
"10619"),
1872 dd_real(
"12874933")/dd_real(
"39442"), dd_real(
"11473457")/dd_real(
"54612"),
1873 dd_real(
"2406734")/dd_real(
"22755"), dd_real(
"166770367")/dd_real(
"4004880"),
1874 dd_real(
"30653165")/dd_real(
"2402928"), dd_real(
"647746389")/dd_real(
"215195552"),
1875 dd_real(
"25346331")/dd_real(
"47074027"), dd_real(
"278270613")/dd_real(
"3900419380"),
1876 dd_real(
"105689791")/dd_real(
"15601677520"), dd_real(
"606046475")/dd_real(
"1379188292768"),
1877 dd_real(
"969715")/dd_real(
"53502994116"), dd_real(
"11098301")/dd_real(
"26204577562592"),
1878 dd_real(
"118999")/dd_real(
"26204577562592"), dd_real(
"18858053")/dd_real(
"1392249205900512960")
1880 const dd_real pade_log_b<dd_real>::array[pade_log_a<dd_real>::array_size] =
1882 dd_real(
"1"), dd_real(
"21")/dd_real(
"2"),
1883 dd_real(
"2100")/dd_real(
"41"), dd_real(
"12635")/dd_real(
"82"),
1884 dd_real(
"341145")/dd_real(
"1066"), dd_real(
"1037799")/dd_real(
"2132"),
1885 dd_real(
"11069856")/dd_real(
"19721"), dd_real(
"9883800")/dd_real(
"19721"),
1886 dd_real(
"6918660")/dd_real(
"19721"), dd_real(
"293930")/dd_real(
"1517"),
1887 dd_real(
"1410864")/dd_real(
"16687"), dd_real(
"88179")/dd_real(
"3034"),
1888 dd_real(
"734825")/dd_real(
"94054"), dd_real(
"305235")/dd_real(
"188108"),
1889 dd_real(
"348840")/dd_real(
"1363783"), dd_real(
"40698")/dd_real(
"1363783"),
1890 dd_real(
"6783")/dd_real(
"2727566"), dd_real(
"9975")/dd_real(
"70916716"),
1891 dd_real(
"266")/dd_real(
"53187537"), dd_real(
"7")/dd_real(
"70916716"),
1892 dd_real(
"7")/dd_real(
"8155422340"), dd_real(
"1")/dd_real(
"538257874440")
1896 struct pade_log_a<qd_real>
1898 static const int array_size = 34;
1899 static const qd_real array[array_size];
1902 struct pade_log_b<qd_real>
1904 static const int array_size = 34;
1905 static const qd_real array[array_size];
1907 const qd_real pade_log_a<qd_real>::array[pade_log_a<qd_real>::array_size] =
1909 qd_real(
"0"), qd_real(
"1"),
1910 qd_real(
"16"), qd_real(
"95201")/qd_real(
"780"),
1911 qd_real(
"30721")/qd_real(
"52"), qd_real(
"7416257")/qd_real(
"3640"),
1912 qd_real(
"1039099")/qd_real(
"195"), qd_real(
"6097772319")/qd_real(
"555100"),
1913 qd_real(
"1564058073")/qd_real(
"85400"), qd_real(
"30404640205")/qd_real(
"1209264"),
1914 qd_real(
"725351278")/qd_real(
"25193"), qd_real(
"4092322670789")/qd_real(
"147429436"),
1915 qd_real(
"4559713849589")/qd_real(
"201040140"), qd_real(
"5049361751189")/qd_real(
"320023080"),
1916 qd_real(
"74979677195")/qd_real(
"8000577"), qd_real(
"16569850691873")/qd_real(
"3481514244"),
1917 qd_real(
"1065906022369")/qd_real(
"515779888"), qd_real(
"335956770855841")/qd_real(
"438412904800"),
1918 qd_real(
"1462444287585964")/qd_real(
"6041877844275"), qd_real(
"397242326339851")/qd_real(
"6122436215532"),
1919 qd_real(
"64211291334131")/qd_real(
"4373168725380"), qd_real(
"142322343550859")/qd_real(
"51080680851480"),
1920 qd_real(
"154355972958659")/qd_real(
"351179680853925"), qd_real(
"167483568676259")/qd_real(
"2937139148960100"),
1921 qd_real(
"4230788929433")/qd_real(
"704913395750424"), qd_real(
"197968763176019")/qd_real(
"392923948371995600"),
1922 qd_real(
"10537522306718")/qd_real(
"319250708052246425"), qd_real(
"236648286272519")/qd_real(
"144249197475035425500"),
1923 qd_real(
"260715545088119")/qd_real(
"4375558990076074573500"), qd_real(
"289596255666839")/qd_real(
"192874640282553367199880"),
1924 qd_real(
"8802625510547")/qd_real(
"361639950529787563499775"), qd_real(
"373831661521439")/qd_real(
"1659204093030665341336967700"),
1925 qd_real(
"446033437968239")/qd_real(
"464577146048586295574350956000"), qd_real(
"53676090078349")/qd_real(
"47386868896955802148583797512000")
1927 const qd_real pade_log_b<qd_real>::array[pade_log_a<qd_real>::array_size] =
1929 qd_real(
"1"), qd_real(
"33")/qd_real(
"2"),
1930 qd_real(
"8448")/qd_real(
"65"), qd_real(
"42284")/qd_real(
"65"),
1931 qd_real(
"211420")/qd_real(
"91"), qd_real(
"573562")/qd_real(
"91"),
1932 qd_real(
"32119472")/qd_real(
"2379"), qd_real(
"92917044")/qd_real(
"3965"),
1933 qd_real(
"603960786")/qd_real(
"17995"), qd_real(
"144626625")/qd_real(
"3599"),
1934 qd_real(
"2776831200")/qd_real(
"68381"), qd_real(
"16692542100")/qd_real(
"478667"),
1935 qd_real(
"12241197540")/qd_real(
"478667"), qd_real(
"1098569010")/qd_real(
"68381"),
1936 qd_real(
"31387686000")/qd_real(
"3624193"), qd_real(
"9939433900")/qd_real(
"2479711"),
1937 qd_real(
"67091178825")/qd_real(
"42155087"), qd_real(
"2683647153")/qd_real(
"4959422"),
1938 qd_real(
"19083713088")/qd_real(
"121505839"), qd_real(
"4708152900")/qd_real(
"121505839"),
1939 qd_real(
"941630580")/qd_real(
"116546417"), qd_real(
"88704330")/qd_real(
"62755763"),
1940 qd_real(
"12902448")/qd_real(
"62755763"), qd_real(
"1542684")/qd_real(
"62755763"),
1941 qd_real(
"6427850")/qd_real(
"2698497809"), qd_real(
"3471039")/qd_real(
"18889484663"),
1942 qd_real(
"8544096")/qd_real(
"774468871183"), qd_real(
"39556")/qd_real(
"79027435835"),
1943 qd_real(
"118668")/qd_real(
"7191496660985"), qd_real(
"10230")/qd_real(
"27327687311743"),
1944 qd_real(
"5456")/qd_real(
"1011124430534491"), qd_real(
"44")/qd_real(
"1011124430534491"),
1945 qd_real(
"11")/qd_real(
"70778710137414370"), qd_real(
"1")/qd_real(
"7219428434016265740")
1952 template<
typename Scalar_T, const index_t LO, const index_t HI >
1963 if (val == Scalar_T(0) || val.
isnan())
1964 return traits_t::NaN();
1966 return pade_approx(pade_log_a<Scalar_T>::array_size,
1967 pade_log_a<Scalar_T>::array,
1968 pade_log_b<Scalar_T>::array,
1973 template<
typename Scalar_T, const index_t LO, const index_t HI >
1981 if (val == Scalar_T(0) || val.
isnan())
1982 return traits_t::NaN();
1984 typedef std::numeric_limits<Scalar_T> limits_t;
1985 static const Scalar_T epsilon = limits_t::epsilon();
1986 static const Scalar_T max_inner_norm =
traits_t::pow(epsilon, 2);
1987 static const Scalar_T max_outer_norm = Scalar_T(6.0/limits_t::digits);
1988 multivector_t Y = val;
1989 multivector_t E = Scalar_T(0);
1992 Scalar_T pow_2_outer_step = Scalar_T(1);
1993 Scalar_T pow_4_outer_step = Scalar_T(1);
1994 for (outer_step = 0, norm_Y_1 =
norm(Y - Scalar_T(1));
1996 ++outer_step, norm_Y_1 =
norm(Y - Scalar_T(1)))
1998 if (Y == Scalar_T(0) || Y.isnan())
1999 return traits_t::NaN();
2002 multivector_t M = Y;
2006 norm(M - Scalar_T(1)) * pow_4_outer_step > max_inner_norm;
2010 E += (M - Scalar_T(1)) * pow_2_outer_step;
2011 pow_2_outer_step *= Scalar_T(2);
2012 pow_4_outer_step *= Scalar_T(4);
2015 return traits_t::NaN();
2017 return pade_log(Y) * pow_2_outer_step - E;
2021 template<
typename Scalar_T, const index_t LO, const index_t HI >
2027 if (val == Scalar_T(0) || val.
isnan())
2028 return traits_t::NaN();
2036 typedef typename traits_t::demoted::type demoted_scalar_t;
2039 const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
2040 const demoted_multivector_t& demoted_i = demoted_multivector_t(i);
2047 typedef typename traits_t::promoted::type promoted_scalar_t;
2050 const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
2051 const promoted_multivector_t& promoted_i = promoted_multivector_t(i);
2062 template<
typename Scalar_T, const index_t LO, const index_t HI >
2070 if (val == Scalar_T(0) || val.
isnan())
2071 return traits_t::NaN();
2074 const Scalar_T realval = val.
scalar();
2077 if (realval < Scalar_T(0))
2083 #if !defined(_GLUCAT_USE_EIGENVALUES) 2084 const multivector_t val2 = val*val;
2085 const Scalar_T real_val2 = val2.
scalar();
2086 if (val2 == real_val2 && real_val2 > 0)
2087 return matrix_log(-i * val, i) + i * pi/Scalar_T(2);
2090 const Scalar_T max_norm = Scalar_T(1.0/9.0);
2091 const Scalar_T scale =
2092 (realval != Scalar_T(0) &&
norm(val/realval - Scalar_T(1)) < max_norm)
2094 : (realval < Scalar_T(0))
2097 if (scale == Scalar_T(0))
2098 return traits_t::NaN();
2101 multivector_t rescale = log_scale;
2102 if (scale < Scalar_T(0))
2103 rescale = i * pi + log_scale;
2104 const multivector_t unitval = val/scale;
2106 return traits_t::NaN();
2107 #if defined(_GLUCAT_USE_EIGENVALUES) 2108 multivector_t scaled_result;
2116 scaled_result =
matrix_log(-i * unitval, i) + i * pi/Scalar_T(2);
2121 scaled_result =
matrix_log(
exp(i*safe_arg) * unitval, i) - i * safe_arg;
2129 multivector_t scaled_result =
cascade_log(unitval);
2131 if (scaled_result.isnan())
2132 return traits_t::NaN();
2134 return scaled_result + rescale;
2138 template<
typename Scalar_T, const index_t LO, const index_t HI >
2144 return traits_t::NaN();
2146 const Scalar_T s =
scalar(val);
2154 typedef typename traits_t::demoted::type demoted_scalar_t;
2157 const demoted_multivector_t& demoted_val = demoted_multivector_t(val);
2163 typedef typename traits_t::promoted::type promoted_scalar_t;
2166 const promoted_multivector_t& promoted_val = promoted_multivector_t(val);
2175 #endif // _GLUCAT_MATRIX_MULTI_IMP_H multivector_t & centre_pm4_qp4(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{p-4,q+4}.
index_set< LO, HI > index_set_t
friend const matrix_multi_t operator|(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
static Multivector_T fast(const Matrix_T &X, index_t level)
Inverse generalized Fast Fourier Transform.
std::vector< Scalar_T > vector_t
multivector_t & centre_qp1_pm1(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{q+1,p-1}.
virtual index_t grade() const =0
Maximum of the grades of each term.
static const matrix_multi_t random(const index_set_t frm, Scalar_T fill=Scalar_T(1))
Random multivector within a frame.
const Multivector< Scalar_T, LO, HI > sqrt(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Square root of multivector with specified complexifier.
index_t count_pos() const
Number of positive indices included in set.
const index_set_t fold() const
Fold this index set within itself as a frame.
static const matrix_multi< Scalar_T, LO, HI > pade_log(const matrix_multi< Scalar_T, LO, HI > &val)
Pade' approximation of log.
Scalar_T abs(const Multivector< Scalar_T, LO, HI > &val)
Absolute value == sqrt(norm)
A framed_multi<Scalar_T,LO,HI> is a framed approximation to a multivector.
const Multivector< Scalar_T, LO, HI > log(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Natural logarithm of multivector with specified complexifier.
static generator_table< Matrix_T > & generator()
Single instance of generator table.
virtual Scalar_T scalar() const =0
Scalar part.
map_t::const_iterator const_iterator
virtual multivector_t & operator|=(const multivector_t &rhs)=0
Transformation via twisted adjoint action.
friend const matrix_multi_t operator&(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
index_t offset_level(const index_t p, const index_t q)
Determine the log2 dim corresponding to signature p, q.
virtual multivector_t & operator*=(const Scalar_T &scr)=0
Product of multivector and scalar.
static const matrix_multi< Scalar_T, LO, HI > pade_approx(const int array_size, const Scalar_T a[], const Scalar_T b[], const matrix_multi< Scalar_T, LO, HI > &X)
Pade' approximation.
virtual const vector_t vector_part() const =0
Vector part of multivector, as a vector_t with respect to frame()
static void db_step(matrix_multi< Scalar_T, LO, HI > &M, matrix_multi< Scalar_T, LO, HI > &Y)
Single step of product form of Denman-Beavers square root iteration.
static const index_t offset_to_super[]
Offsets between the current signature and that of the real superalgebra.
virtual multivector_t & operator/=(const Scalar_T &scr)=0
Quotient of multivector and scalar.
static const matrix_multi< Scalar_T, LO, HI > cascade_log(const matrix_multi< Scalar_T, LO, HI > &val)
Incomplete square root cascade and Pade' approximation of log.
friend const matrix_multi_t operator^(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
index_t max() const
Maximum member.
friend std::istream & operator>>(std::istream &s, multivector_t &val)
virtual const multivector_t inv() const =0
Geometric multiplicative inverse.
_GLUCAT_CLIFFORD_ALGEBRA_OPERATIONS multivector_t & operator=(const multivector_t &rhs)
Assignment operator.
virtual Scalar_T operator[](const index_set_t ist) const =0
Subscripting: map from index set to scalar coordinate.
friend const matrix_multi_t operator%(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
Extra traits which extend numeric limits.
Table of basis elements used as a cache by basis_element()
virtual const multivector_t reverse() const =0
Reversion, eg. {1}*{2} -> {2}*{1}.
friend Scalar_T star(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
const Multivector< Scalar_T, LO, HI > clifford_exp(const Multivector< Scalar_T, LO, HI > &val)
Exponential of multivector.
Matrix_T::value_type norm_frob2(const Matrix_T &val)
Square of Frobenius norm.
const Multivector< Scalar_T, LO, HI > pow(const Multivector< Scalar_T, LO, HI > &lhs, int rhs)
Integer power of multivector.
virtual const multivector_t even() const =0
Even part of multivector, sum of even grade terms.
bool isnan(const Matrix_T &m)
Not a Number.
index_t count() const
Cardinality: Number of indices included in set.
static void check_complex(const Multivector< Scalar_T, LO, HI > &val, const Multivector< Scalar_T, LO, HI > &i, const bool prechecked=false)
Check that i is a valid complexifier for val.
matrix_multi multivector_t
index_set_t m_frame
Index set representing the frame for the subalgebra which contains the multivector.
const framed_multi< Other_Scalar_T, LO, HI > fast_framed_multi() const
Use inverse generalized FFT to construct a framed_multi_t.
A matrix_multi<Scalar_T,LO,HI> is a matrix approximation to a multivector.
virtual Scalar_T norm() const =0
Scalar_T norm == sum of norm of coordinates.
Structure containing classification of eigenvalues.
index_t min() const
Minimum member.
index_set< LO, HI > index_set_t
virtual const multivector_t operator-() const =0
Unary -.
std::pair< const index_set_t, Scalar_T > term_t
virtual Scalar_T max_abs() const =0
Maximum of absolute values of components of multivector: multivector infinity norm.
const framed_multi< Scalar_T, LO, HI > exp(const framed_multi< Scalar_T, LO, HI > &val)
Exponential of multivector.
virtual const multivector_t operator()(index_t grade) const =0
Pure grade-vector part.
error< multivector_t > error_t
const matrix_multi_t fast_matrix_multi(const index_set_t frm) const
Use generalized FFT to construct a matrix_multi_t.
eig_case_t m_eig_case
What kind of eigenvalues does the matrix contain?
friend const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > matrix_log(const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &val, const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &i)
Abstract exception class.
std::pair< const index_set_t, Scalar_T > term_t
virtual const multivector_t involute() const =0
Main involution, each {i} is replaced by -{i} in each term, eg. {1} -> -{1}.
multivector_t & centre_pp4_qm4(index_t &p, index_t &q)
Subalgebra isomorphism: R_{p,q} to R_{p+4,q-4}.
static const precision_t function_precision
Precision used for exp, log and sqrt functions.
virtual const multivector_t pow(int m) const =0
*this to the m
virtual const multivector_t outer_pow(int m) const =0
Outer product power.
virtual void write(const std::string &msg="") const =0
Write formatted multivector to output.
friend const matrix_multi_t operator/(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
ublas::compressed_matrix< Scalar_T, orientation_t > matrix_t
multivector_t unfold(const index_set_t frm) const
Subalgebra isomorphism: unfold each term within the given frame.
const basis_matrix_t basis_element(const index_set< LO, HI > &ist) const
Create a basis element matrix within the current frame.
friend const index_set< Other_LO, Other_HI > reframe(const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &lhs, const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &rhs, matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &lhs_reframed, matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &rhs_reframed)
virtual const multivector_t odd() const =0
Odd part of multivector, sum of odd grade terms.
static basis_table & basis()
Single instance of basis table.
framed_multi< Scalar_T, LO, HI > framed_multi_t
matrix_t::size_type matrix_index_t
virtual multivector_t & operator-=(const multivector_t &rhs)=0
Geometric difference.
friend const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > matrix_sqrt(const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &val, const matrix_multi< Other_Scalar_T, Other_LO, Other_HI > &i)
static const std::string classname()
Class name used in messages.
matrix_t m_matrix
Matrix value representing the multivector within the folded frame.
static const matrix_multi< Scalar_T, LO, HI > db_sqrt(const matrix_multi< Scalar_T, LO, HI > &val)
Product form of Denman-Beavers square root iteration.
Index set class based on std::bitset<> in Gnu standard C++ library.
virtual const multivector_t conj() const =0
Conjugation, reverse o involute == involute o reverse.
virtual multivector_t & operator^=(const multivector_t &rhs)=0
Outer product.
const RHS_T::expression_type mono_prod(const ublas::matrix_expression< LHS_T > &lhs, const ublas::matrix_expression< RHS_T > &rhs)
Product of monomial matrices.
multivector_t & operator+=(const term_t &rhs)
Add a term, if non-zero.
virtual const multivector_t truncated(const Scalar_T &limit=Scalar_T(DEFAULT_TRUNCATION)) const =0
Remove all terms with relative size smaller than limit.
virtual bool operator==(const multivector_t &val) const =0
Test for equality of multivectors.
Matrix_T::value_type trace(const Matrix_T &val)
Matrix trace.
Scalar_T m_safe_arg
Argument such that exp(pi-m_safe_arg) lies between arguments of eigenvalues.
ublas::compressed_matrix< int, orientation_t > basis_matrix_t
int index_t
Size of index_t should be enough to represent LO, HI.
const RHS_T signed_perm_nork(const LHS_T &lhs, const RHS_T &rhs)
Left inverse of Kronecker product where lhs is a signed permutation matrix.
virtual const multivector_t pure() const =0
Pure part.
eig_genus< Matrix_T > classify_eigenvalues(const Matrix_T &val)
Classify the eigenvalues of a matrix.
static Matrix_Index_T folded_dim(const index_set< LO, HI > &sub)
Determine the matrix dimension of the fold of a subalegbra.
LHS_T pos_mod(LHS_T lhs, RHS_T rhs)
Modulo function which works reliably for lhs < 0.
virtual multivector_t & operator&=(const multivector_t &rhs)=0
Inner product.
static const framed_multi_t random(const index_set_t frm, Scalar_T fill=Scalar_T(1))
Random multivector within a frame.
virtual Scalar_T quad() const =0
Scalar_T quadratic form == (rev(x)*x)(0)
index_t count_neg() const
Number of negative indices included in set.
virtual bool isnan() const =0
Check if a multivector contains any IEEE NaN values.
friend const matrix_multi_t operator*(const matrix_multi_t &lhs, const matrix_multi_t &rhs)
virtual const index_set_t frame() const =0
Subalgebra generated by all generators of terms of given multivector.
matrix_multi()
Default constructor.
virtual multivector_t & operator%=(const multivector_t &rhs)=0
Contraction.