Matrix functions and scalar powers of gr_mat

[fredrik-johansson]

• Adds gr_mat_func_jordan for computing arbitrary analytic functions of square matrices using Jordan decomposition (generalized out from the previous code for exp and log)

• Improves that method to fall back on diagonalization over inexact fields in the case where one has n distinct eigenvalues (would previously attempt to compute ranks and fail)

• Adds gr_mat_pow_si, gr_mat_pow_fmpq, gr_mat_pow_scalar, gr_mat_sqrt, gr_mat_rsqrt

• Also fixes a couple of other minor issues caught while testing (including a memory leak in gr_mat_fflu).

So, amusingly one can now do this:

>>> Mat(CC)([[1, -1], [1, 1]]) ** RR.pi()
[[([-2.32073556181001 +/- 5.94e-15] + [+/- 4.77e-15]*I), ([-1.85449839019252 +/- 6.53e-15] + [+/- 5.20e-15]*I)],
[([1.85449839019252 +/- 6.74e-15] + [+/- 5.31e-15]*I), ([-2.32073556181001 +/- 5.59e-15] + [+/- 4.81e-15]*I)]]
>>> _ ** (1 / RR.pi())
[[([1.00000000000 +/- 6.60e-13] + [+/- 6.60e-13]*I), ([-1.0000000000 +/- 9.82e-13] + [+/- 9.80e-13]*I)],
[([1.00000000000 +/- 6.61e-13] + [+/- 6.59e-13]*I), ([1.0000000000 +/- 9.81e-13] + [+/- 9.80e-13]*I)]]

Note that this code has many limitations. For example, it is not very useful for the purpose of detecting whether matrices over non-closed fields or rings are perfect powers:

>>> (Mat(ZZ)([[1,2],[3,4]]) ** 2)
[[7, 10],
[15, 22]]
>>> (Mat(ZZ)([[1,2],[3,4]]) ** 2).sqrt()
Traceback (most recent call last):
  ...
FlintUnableError: failed to compute sqrt(x) in {Matrices (any shape) over Integer ring (fmpz)} for {x = [[7, 10], [15, 22]]}

This does work though:

>>> (Mat(QQbar)([[1,2],[3,4]]) ** 2).sqrt()
[[Root a = 1.56670 of 11*a^2-27, Root a = 1.74078 of 33*a^2-100],
[Root a = 2.61116 of 11*a^2-75, Root a = 4.17786 of 11*a^2-192]]
>>> _ ** 2
[[7, 10],
[15, 22]]
>>> (Mat(CC)([[1,2],[3,4]]) ** 2).sqrt()
[[([1.566698903601 +/- 3.84e-13] + [+/- 9.85e-14]*I), ([1.740776559557 +/- 1.35e-13] + [+/- 1.08e-13]*I)],
[([2.611164839335 +/- 5.97e-13] + [+/- 1.22e-13]*I), ([4.177863742937 +/- 3.96e-13] + [+/- 1.34e-13]*I)]]
>>> _ ** 2
[[([7.000000000000 +/- 8.54e-13] + [+/- 8.03e-13]*I), ([10.0000000000 +/- 1.09e-12] + [+/- 1.03e-12]*I)],
[([15.00000000000 +/- 1.41e-12] + [+/- 1.31e-12]*I), ([22.00000000000 +/- 1.73e-12] + [+/- 1.61e-12]*I)]]

In the following example, there is a negative real eigenvalue and the error bounds blow up on the branch cut. Here we would need a non-principal internal square root:

>>> Mat(CC)([[1,2],[3,4]]).sqrt()
[[([0.5536886 +/- 3.29e-8] + [+/- 0.465]*I), ([0.8069607 +/- 2.71e-8] + [+/- 0.213]*I)],
[([1.2104411 +/- 3.97e-8] + [+/- 0.319]*I), ([1.7641297 +/- 4.24e-8] + [+/- 0.146]*I)]]