| 307 | |
| 308 | |
| 309 | qmatrix getPowerOfDiagonalMatrix(qmatrix m, qcomp p) { |
| 310 | DEMAND( isDiagonal(m) ); |
| 311 | |
| 312 | qmatrix out = getZeroMatrix(m.size()); |
| 313 | |
| 314 | // pow(qcomp,qcomp) introduces wildly erroneous |
| 315 | // imaginary components when both base is real |
| 316 | // and negative, and exponent is real and integer |
| 317 | // (so ergo does not produce complex numbers). |
| 318 | // We divert to real-pow in that scenario! |
| 319 | |
| 320 | for (size_t i=0; i<m.size(); i++) { |
| 321 | bool mIsRe = std::imag(m[i][i]) == 0; |
| 322 | bool mIsNeg = std::real(m[i][i]) < 0; |
| 323 | bool pIsRe = std::imag(p) == 0; |
| 324 | bool pIsInt = std::trunc(std::real(p)) == std::real(p); |
| 325 | |
| 326 | // use pow(qreal,qreal) or pow(qcomp,qcomp) |
| 327 | out[i][i] = (mIsRe && mIsNeg && pIsRe && pIsInt)? |
| 328 | qcomp(std::pow(std::real(m[i][i]), std::real(p)),0): |
| 329 | std::pow(m[i][i], p); |
| 330 | } |
| 331 | |
| 332 | return out; |
| 333 | } |
| 334 | |
| 335 | |
| 336 | qmatrix getExponentialOfDiagonalMatrix(qmatrix m) { |
no test coverage detected