package scientific.linear
import std.math.*
import std.unittest.*
import std.unittest.testmacro.*
import scientific.numbers.*
public func expm(a: Matrix<Float64>): Matrix<Float64> {
if (!a.isSquare()) {
throw IllegalArgumentException("expm: matrix must be square")
}
var A_L1 = norm(a, r'1')
var e: Matrix<Float64> = eye(a.getRows(), a.getCols())
var ucopy: Matrix<Float64> = e
var vcopy: Matrix<Float64> = e
var n_squarings: Int64 = 0
if (A_L1 < 0.01495585217958292) {
var (u,v) = pade3(a)
ucopy = u
vcopy = v
} else if (A_L1 < 0.253939833006323) {
var (u,v) = pade5(a)
ucopy = u
vcopy = v
} else if (A_L1 < 0.9504178996162932) {
var (u,v) = pade7(a)
ucopy = u
vcopy = v
} else if (A_L1 < 2.097847961257068) {
var (u,v) = pade9(a)
ucopy = u
vcopy = v
} else {
var maxnorm = 5.371920351148152
n_squarings = max(0, Int64(ceil(log2(A_L1 / maxnorm))))
var aa = a * (1.0 / exp2(Float64(n_squarings)))
var (u,v) = pade13(aa)
ucopy = u
vcopy = v
}
var p = ucopy + vcopy
var q = -1.0 * ucopy + vcopy
var r = solve(q,p)
// Squaring step to undo scaling
for (i in 0..n_squarings) {
r = cj_dgemm(r, r)
}
return r
}
public func expm(a: Matrix<Complex64>): Matrix<Complex64> {
if (!a.isSquare()) {
throw IllegalArgumentException("expm: matrix must be square")
}
var A_L1 = norm(a, r'1')
var e: Matrix<Complex64> = eye(a.getRows(), a.getCols())
var ucopy: Matrix<Complex64> = e
var vcopy: Matrix<Complex64> = e
var n_squarings: Int64 = 0
if (A_L1 < 0.01495585217958292) {
var (u,v) = pade3(a)
ucopy = u
vcopy = v
} else if (A_L1 < 0.253939833006323) {
var (u,v) = pade5(a)
ucopy = u
vcopy = v
} else if (A_L1 < 0.9504178996162932) {
var (u,v) = pade7(a)
ucopy = u
vcopy = v
} else if (A_L1 < 2.097847961257068) {
var (u,v) = pade9(a)
ucopy = u
vcopy = v
} else {
var maxnorm = 5.371920351148152
n_squarings = max(0, Int64(ceil(log2(A_L1 / maxnorm))))
var temp = 1.0 / exp2(Float64(n_squarings))
var complexTemp = Complex64(temp, 0.0)
var aa = a * complexTemp
var (u,v) = pade13(aa)
ucopy = u
vcopy = v
}
var p = ucopy + vcopy
var temp2 = Complex64(-1.0, 0.0)
var q = temp2 * ucopy + vcopy
var r = solve(q,p)
// Squaring step to undo scaling
for (i in 0..n_squarings) {
r = cj_zgemm(r, r)
}
return r
}
public func pade3(a: Matrix<Float64>): (Matrix<Float64>, Matrix<Float64>) {
var b = vector<Float64> (
[120.0, 60.0, 12.0, 1.0]
)
var ident: Matrix<Float64> = eye(a.getRows(), a.getCols())
var a2 = cj_dgemm(a, a)
var u = cj_dgemm(a, (a2*b[3] + ident*b[1]))
var v = a2*b[2] + ident*b[0]
return (u,v)
}
public func pade3(a: Matrix<Complex64>): (Matrix<Complex64>, Matrix<Complex64>) {
var b = vector<Complex64>([
Complex64(120.0, 0.0), Complex64(60.0, 0.0), Complex64(12.0, 0.0), Complex64(1.0, 0.0)
])
var ident: Matrix<Complex64> = eye(a.getRows(), a.getCols())
var a2 = cj_zgemm(a, a)
var u = cj_zgemm(a, (a2*b[3] + ident*b[1]))
var v = a2*b[2] + ident*b[0]
return (u,v)
}
public func pade5(a: Matrix<Float64>): (Matrix<Float64>, Matrix<Float64>) {
var b = vector<Float64> (
[30240.0, 15120.0, 3360.0, 420.0, 30.0, 1.0]
)
var ident: Matrix<Float64> = eye(a.getRows(), a.getCols())
var a2 = cj_dgemm(a, a)
var a4 = cj_dgemm(a2, a2)
var u = cj_dgemm(a, (b[5]*a4 + b[3]*a2 + b[1]*ident))
var v = b[4]*a4 + b[2]*a2 + b[0]*ident
return (u,v)
}
public func pade5(a: Matrix<Complex64>): (Matrix<Complex64>, Matrix<Complex64>) {
var b = vector<Complex64>([
Complex64(30240.0, 0.0), Complex64(15120.0, 0.0), Complex64(3360.0, 0.0), Complex64(420.0, 0.0),
Complex64(30.0, 0.0), Complex64(1.0, 0.0)
])
var ident: Matrix<Complex64> = eye(a.getRows(), a.getCols())
var a2 = cj_zgemm(a, a)
var a4 = cj_zgemm(a2, a2)
var u = cj_zgemm(a, (b[5]*a4 + b[3]*a2 + b[1]*ident))
var v = b[4]*a4 + b[2]*a2 + b[0]*ident
return (u,v)
}
public func pade7(a: Matrix<Float64>): (Matrix<Float64>, Matrix<Float64>) {
var b = vector<Float64> (
[17297280.0, 8648640.0, 1995840.0, 277200.0, 25200.0, 1512.0, 56.0, 1.0]
)
var ident: Matrix<Float64> = eye(a.getRows(), a.getCols())
var a2 = cj_dgemm(a, a)
var a4 = cj_dgemm(a2, a2)
var a6 = cj_dgemm(a4, a2)
var u = cj_dgemm(a, (b[7]*a6 + b[5]*a4 + b[3]*a2 + b[1]*ident))
var v = b[6]*a6 + b[4]*a4 + b[2]*a2 + b[0]*ident
return (u,v)
}
public func pade7(a: Matrix<Complex64>): (Matrix<Complex64>, Matrix<Complex64>) {
var b = vector<Complex64>([
Complex64(17297280.0, 0.0), Complex64(8648640.0, 0.0), Complex64(1995840.0, 0.0), Complex64(277200.0, 0.0),
Complex64(25200.0, 0.0), Complex64(1512.0, 0.0), Complex64(56.0, 0.0), Complex64(1.0, 0.0)
])
var ident: Matrix<Complex64> = eye(a.getRows(), a.getCols())
var a2 = cj_zgemm(a, a)
var a4 = cj_zgemm(a2, a2)
var a6 = cj_zgemm(a4, a2)
var u = cj_zgemm(a, (b[7]*a6 + b[5]*a4 + b[3]*a2 + b[1]*ident))
var v = b[6]*a6 + b[4]*a4 + b[2]*a2 + b[0]*ident
return (u,v)
}
public func pade9(a: Matrix<Float64>): (Matrix<Float64>, Matrix<Float64>) {
var b = vector<Float64> (
[17643225600.0, 8821612800.0, 2075673600.0, 302702400.0, 30270240.0,
2162160.0, 110880.0, 3690.0, 90.0, 1.0]
)
var ident: Matrix<Float64> = eye(a.getRows(), a.getCols())
var a2 = cj_dgemm(a, a)
var a4 = cj_dgemm(a2, a2)
var a6 = cj_dgemm(a4, a2)
var a8 = cj_dgemm(a6, a2)
var u = cj_dgemm(a, (b[9]*a8 + b[7]*a6 + b[5]*a4 + b[3]*a2 + b[1]*ident))
var v = b[8]*a8 + b[6]*a6 + b[4]*a4 + b[2]*a2 + b[0]*ident
return (u,v)
}
public func pade9(a: Matrix<Complex64>): (Matrix<Complex64>, Matrix<Complex64>) {
var b = vector<Complex64>([
Complex64(17643225600.0, 0.0), Complex64(8821612800.0, 0.0), Complex64(2075673600.0, 0.0),
Complex64(302702400.0, 0.0), Complex64(30270240.0, 0.0), Complex64(2162160.0, 0.0),
Complex64(110880.0, 0.0), Complex64(3690.0, 0.0), Complex64(90.0, 0.0), Complex64(1.0, 0.0)
])
var ident: Matrix<Complex64> = eye(a.getRows(), a.getCols())
var a2 = cj_zgemm(a, a)
var a4 = cj_zgemm(a2, a2)
var a6 = cj_zgemm(a4, a2)
var a8 = cj_zgemm(a6, a2)
var u = cj_zgemm(a, (b[9]*a8 + b[7]*a6 + b[5]*a4 + b[3]*a2 + b[1]*ident))
var v = b[8]*a8 + b[6]*a6 + b[4]*a4 + b[2]*a2 + b[0]*ident
return (u,v)
}
public func pade13(a: Matrix<Float64>): (Matrix<Float64>, Matrix<Float64>) {
var b = vector<Float64> (
[64764752532480000.0, 32382376266240000.0, 7771770303897600.0,
1187353796428800.0, 129060195264000.0, 10559470521600.0, 670442572800.0,
33522128640.0, 1323241920.0, 40840800.0, 960960.0, 16380.0, 182.0, 1.0]
)
var ident: Matrix<Float64> = eye(a.getRows(), a.getCols())
var a2 = cj_dgemm(a, a)
var a4 = cj_dgemm(a2, a2)
var a6 = cj_dgemm(a4, a2)
var u = cj_dgemm(a, cj_dgemm(a6, b[13]*a6 + b[11]*a4 + b[9]*a2) + b[7]*a6 + b[5]*a4 + b[3]*a2 + b[1]*ident)
var v = cj_dgemm(a6, b[12]*a6 + b[10]*a4 + b[8]*a2) + b[6]*a6 + b[4]*a4 + b[2]*a2 + b[0]*ident
return (u,v)
}
public func pade13(a: Matrix<Complex64>): (Matrix<Complex64>, Matrix<Complex64>) {
var b = vector<Complex64>([
Complex64(64764752532480000.0, 0.0), Complex64(32382376266240000.0, 0.0), Complex64(7771770303897600.0, 0.0),
Complex64(1187353796428800.0, 0.0), Complex64(129060195264000.0, 0.0), Complex64(10559470521600.0, 0.0),
Complex64(670442572800.0, 0.0), Complex64(33522128640.0, 0.0), Complex64(1323241920.0, 0.0),
Complex64(40840800.0, 0.0), Complex64(960960.0, 0.0), Complex64(16380.0, 0.0),
Complex64(182.0, 0.0), Complex64(1.0, 0.0)
])
var ident: Matrix<Complex64> = eye(a.getRows(), a.getCols())
var a2 = cj_zgemm(a, a)
var a4 = cj_zgemm(a2, a2)
var a6 = cj_zgemm(a4, a2)
var u = cj_zgemm(a, cj_zgemm(a6, b[13]*a6 + b[11]*a4 + b[9]*a2) + b[7]*a6 + b[5]*a4 + b[3]*a2 + b[1]*ident)
var v = cj_zgemm(a6, b[12]*a6 + b[10]*a4 + b[8]*a2) + b[6]*a6 + b[4]*a4 + b[2]*a2 + b[0]*ident
return (u,v)
}
public func cosm(a: Matrix<Complex64>): Matrix<Complex64> {
return Complex64(0.5, 0.0) * (expm(Complex64(0.0, 1.0) * a ) + expm(Complex64(0.0, -1.0) * a))
}
public func cosm(a: Matrix<Float64>): Matrix<Float64> {
var temp = expm(Complex64(0.0, 1.0) * toComplex64(a))
return getReal(temp)
}
public func sinm(a: Matrix<Complex64>): Matrix<Complex64> {
return Complex64(0.0, -0.5) * (expm(Complex64(0.0, 1.0) * a ) - expm(Complex64(0.0, -1.0) * a))
}
public func sinm(a: Matrix<Float64>): Matrix<Float64> {
var temp = expm(Complex64(0.0, 1.0) * toComplex64(a))
return getImag(temp)
}
public func tanm(a: Matrix<Float64>): Matrix<Float64> {
var result = solve(cosm(a), sinm(a))
return result
}
public func tanm(a: Matrix<Complex64>): Matrix<Complex64> {
var result = solve(cosm(a), sinm(a))
return result
}
public func coshm(a: Matrix<Float64>): Matrix<Float64> {
var result = 0.5 * (expm(a) + expm(-1.0 * a))
return result
}
public func coshm(a: Matrix<Complex64>): Matrix<Complex64> {
var result = Complex64(0.5, 0.0) * (expm(a) + expm(Complex64(-1.0, 0.0) * a))
return result
}
public func sinhm(a: Matrix<Float64>): Matrix<Float64> {
var result = 0.5 * (expm(a) - expm(-1.0 * a))
return result
}
public func sinhm(a: Matrix<Complex64>): Matrix<Complex64> {
var result = Complex64(0.5, 0.0) * (expm(a) - expm(Complex64(-1.0, 0.0) * a))
return result
}
public func tanhm(a: Matrix<Float64>): Matrix<Float64> {
var result = solve(coshm(a), sinhm(a))
return result
}
public func tanhm(a: Matrix<Complex64>): Matrix<Complex64> {
var result = solve(coshm(a), sinhm(a))
return result
}
public func toComplex64(a: Matrix<Float64>): Matrix<Complex64> {
let row = a.getRows()
let col = a.getCols()
var res = Matrix<Complex64>(row, col)
for (i in 0..row) {
for (j in 0..col) {
res[i,j] = Complex64(a[i,j], 0.0)
}
}
return res
}
public func getReal(a: Matrix<Complex64>): Matrix<Float64> {
let row = a.getRows()
let col = a.getCols()
var res = Matrix<Float64>(row, col)
for (i in 0..row) {
for (j in 0..col) {
res[i,j] = a[i,j].real
}
}
return res
}
public func getImag(a: Matrix<Complex64>): Matrix<Float64> {
let row = a.getRows()
let col = a.getCols()
var res = Matrix<Float64>(row, col)
for (i in 0..row) {
for (j in 0..col) {
res[i,j] = a[i,j].imag
}
}
return res
}
public func matExp(a: Matrix<Complex64>): Matrix<Complex64> {
let row = a.getRows()
let col = a.getCols()
var res = Matrix<Complex64>(row, col)
for (i in 0..row) {
for (j in 0..col) {
res[i,j] = Complex64(cos(a[i,j].imag) * exp(a[i,j].real), sin(a[i,j].imag) * exp(a[i,j].real))
}
}
return res
}
public func getExp(a: Matrix<Complex64>): Matrix<Complex64> {
let row = a.getRows()
let col = a.getCols()
var res = Matrix<Complex64>(row, row)
for (i in 0..row) {
for (j in 0..row) {
if (i==0 || j==0) {
res[i,j] = Complex64(1.0, 0.0)
} else {
var index = (i * j) % row
res[i,j] = a[index, 0]
}
}
}
return res
}
public func signm(A: Matrix<Float64>): Matrix<Float64> {
if (A.getRows() != A.getCols()) {
throw IllegalArgumentException("signm: input must be a square matrix.")
}
// Count number of steps
var i: Int64 = 0
var curA: Matrix<Float64> = A
while (i < 100) {
let nextA = 0.5 * (curA + inv(curA))
if (norm(nextA - curA, r'F') < 1e-9) {
return nextA
} else {
i += 1
curA = nextA
}
}
throw IllegalArgumentException("signm: convergence failed.")
}
public func sqrtm(A: Matrix<Float64>): Matrix<Float64> {
if (A.getRows() != A.getCols()) {
throw IllegalArgumentException("sqrtm: input must be a square matrix.")
}
// Count number of steps
var i: Int64 = 0
let sz = A.getRows()
var X: Matrix<Float64> = A
var Y: Matrix<Float64> = eye(sz)
while (i < 100) {
let nextX = 0.5 * (X + inv(Y))
let nextY = 0.5 * (Y + inv(X))
if (norm(nextX - X, r'F') < 1e-9 && norm(nextY - Y, r'F') < 1e-9) {
return nextX
} else {
i += 1
X = nextX
Y = nextY
}
}
throw IllegalArgumentException("sqrtm: convergence failed.")
}
public func logm(A: Matrix<Float64>): Matrix<Float64> {
if (A.getRows() != A.getCols()) {
throw IllegalArgumentException("sqrtm: input must be a square matrix.")
}
var k: Int64 = 0
var curA: Matrix<Float64> = A
let I: Matrix<Float64> = eye(A.getRows())
while (k < 100) {
curA = sqrtm(curA)
k += 1
if (norm(curA - I, r'F') < 1e-2) {
break
}
}
// Use Pade approximation to compute log(curA)
let diff = curA - I
let diff2 = diff * diff
let diff3 = diff2 * diff
let D = 60.0 * I + 90.0 * diff + 36.0 * diff2 + 3.0 * diff3
let N = 60.0 * diff + 60.0 * diff2 + 11.0 * diff3
let logA_pade = inv(D) * N
return pow(2.0, Float64(k)) * logA_pade
}
@Test
public class TestMatFuncs {
@TestCase
func testExpm(): Unit {
var a1 = matrix<Float64>(
[[1.0, 2.0],
[-1.0, 3.0]]
)
var a2 = matrix<Float64>(
[[1.0, 1.0, 0.0],
[ 0.0, 0.0, 2.0],
[ 0.0, 0.0, -1.0]]
)
let r1 = expm(a1)
let r2 = expm(a2)
var expected_r1 = matrix<Float64>(
[[-2.225352, 12.435353],
[-6.217676, 10.210000]]
)
var expected_r2 = matrix<Float64>(
[[2.718282, 1.718282, 1.086161],
[0.000000, 1.000000, 1.264241],
[0.000000, 0.000000, 0.367879]]
)
@Assert(approxEqual(r1, expected_r1, atol:1e-6))
@Assert(approxEqual(r2, expected_r2, atol:1e-6))
}
@TestCase
func testComplex64Expm(): Unit {
var a = matrix<Complex64>([
[Complex64(0.0, 0.0), Complex64(0.0, 1.0), Complex64(0.0, 0.0)],
[Complex64(0.0, 0.0), Complex64(0.0, 0.0), Complex64(0.0, 2.0)],
[Complex64(0.0, 0.0), Complex64(0.0, 0.0), Complex64(0.0, -1.0)]
])
let r = expm(a)
var expected_r = matrix<Complex64>([
[Complex64(1.0, 0.0), Complex64(0.0, 1.0), Complex64(-0.91939539, 0.31705803)],
[Complex64(0.0, 0.0), Complex64(1.0, 0.0), Complex64( 0.91939539, 1.68294197)],
[Complex64(0.0, 0.0), Complex64(0.0, 0.0), Complex64( 0.54030231, -0.841470985)]
])
@Assert(approxEqualC(r, expected_r, atol:1e-6))
}
@TestCase
func testCosm(): Unit {
var a = matrix<Complex64>([
[Complex64( 1.0, 0.0), Complex64(2.0, 0.0)],
[Complex64(-1.0, 0.0), Complex64(3.0, 0.0)]
])
let r = cosm(Complex64(0.0, 1.0) * a)
var expected_r = matrix<Complex64>([
[Complex64(-1.01917479, 0.0), Complex64(6.1037956, 0.0)],
[Complex64(-3.0518978, 0.0), Complex64(5.08462081, 0.0)]
])
@Assert(approxEqualC(r, expected_r, atol:1e-6))
}
@TestCase
func testCosmReal(): Unit {
var a = matrix<Float64>(
[[1.0, 2.0],
[-1.0, 3.0]]
)
let r = cosm(a)
var expected_r = matrix<Float64>(
[[0.4264593, -2.13721484],
[1.06860742, -1.71075555]]
)
@Assert(approxEqual(r, expected_r, atol:1e-6))
}
@TestCase
func testSinm(): Unit {
var a = matrix<Complex64>([
[Complex64( 1.0, 0.0), Complex64(2.0, 0.0)],
[Complex64(-1.0, 0.0), Complex64(3.0, 0.0)]
])
let r = sinm(Complex64(0.0, 1.0) * a)
var expected_r = matrix<Complex64>([
[Complex64(0.0, -1.20617747), Complex64(0.0, 6.33155703)],
[Complex64(0.0, -3.16577851), Complex64(0.0, 5.12537955)]
])
@Assert(approxEqualC(r, expected_r, atol:1e-6))
}
@TestCase
func testSinmReal(): Unit {
var a = matrix<Float64>(
[[1.0, 2.0],
[-1.0, 3.0]]
)
let r = sinm(a)
var expected_r = matrix<Float64>(
[[1.89217551, -0.97811252],
[0.48905626, 0.91406299]]
)
@Assert(approxEqual(r, expected_r, atol:1e-6))
}
@TestCase
func testExpmSinmCosm(): Unit {
var a = matrix<Float64>(
[[1.0, 2.0],
[-1.0, 3.0]]
)
var r1 = expm(Complex64(0.0, 1.0) * toComplex64(a))
var r2 = toComplex64(cosm(a)) + Complex64(0.0, 1.0) * toComplex64(sinm(a))
var expected_r = matrix<Complex64>([
[Complex64(0.42645930, 1.89217551), Complex64(-2.13721484, -0.97811252)],
[Complex64(1.06860742, 0.48905626), Complex64(-1.71075555, 0.91406299)]
])
@Assert(approxEqualC(r1, expected_r, atol:1e-6))
@Assert(approxEqualC(r2, expected_r, atol:1e-6))
}
@TestCase
func testTanm(): Unit {
var a = matrix<Float64>(
[[1.0, 3.0],
[1.0, 4.0]]
)
let s = sinm(a)
let c = cosm(a)
let t = tanm(a)
let sc = cj_dgemm(s, inv(c))
var expected_t = matrix<Float64>(
[[-2.00876993, -8.41880636],
[-2.80626879, -10.42757629]]
)
@Assert(approxEqual(t, expected_t, atol:1e-6))
@Assert(approxEqual(sc, expected_t, atol:1e-6))
}
@TestCase
func testComplex64Tanm(): Unit {
var a = matrix<Complex64>([
[Complex64(1.0, 0.0), Complex64(0.0, 2.0), Complex64(0.0, 3.0)],
[Complex64(0.0, 4.0), Complex64(5.0, 0.0), Complex64(0.0, 6.0)],
[Complex64(0.0, 7.0), Complex64(0.0, 8.0), Complex64(0.0, 0.0)]
])
let t: Matrix<Complex64> = tanm(a)
var expected_t = matrix<Complex64>([
[Complex64(-0.10065095, -0.52460407), Complex64(0.04550901, 0.43519723), Complex64(-0.08926418, 0.41064267)],
[Complex64(0.08990665, 0.80916157), Complex64(0.28488148, -0.33118316), Complex64(0.0641805, 0.6943741)],
[Complex64(-0.20754218, 0.99898816), Complex64(0.08520355, 0.90542117), Complex64(-0.18930083, -0.1468804)]
])
@Assert(approxEqualC(t, expected_t, atol:1e-6))
}
@TestCase
func testCoshm(): Unit {
var a = matrix<Float64>(
[[1.0, 3.0],
[1.0, 4.0]]
)
let r = coshm(a)
var expected_r = matrix<Float64>(
[[11.24592233, 38.76236492],
[12.92078831, 50.00828725]]
)
@Assert(approxEqual(r, expected_r, atol:1e-6))
}
@TestCase
func testSinhm(): Unit {
var a = matrix<Float64>(
[[1.0, 3.0],
[1.0, 4.0]]
)
let r = sinhm(a)
var expected_r = matrix<Float64>(
[[10.57300653, 39.28826594],
[13.09608865, 49.86127247]]
)
@Assert(approxEqual(r, expected_r, atol:1e-6))
}
@TestCase
func testTanhm(): Unit {
var a = matrix<Float64>(
[[1.0, 3.0],
[1.0, 4.0]]
)
let t = tanhm(a)
let s = sinhm(a)
let c = coshm(a)
let sc = cj_dgemm(s, inv(c))
var expected_t = matrix<Float64>(
[[0.3428582, 0.51987926],
[0.17329309, 0.86273746]]
)
@Assert(approxEqual(t, expected_t, atol:1e-6))
@Assert(approxEqual(sc, expected_t, atol:1e-6))
}
@TestCase
func testSignm(): Unit {
var a = matrix<Float64>([
[1.0, 2.0, 3.0],
[1.0, 2.0, 1.0],
[1.0, 1.0, 1.0]
])
var expected_res = matrix<Float64>([
[-0.13127464, 0.15312833, 1.83967663],
[ 0.22292636, 0.96982489, -0.36252242],
[ 0.51565075, -0.06979803, 0.16144975]
])
var res = signm(a)
@Assert(approxEqual(res, expected_res, atol:1e-7))
}
@TestCase
func testSqrtm(): Unit {
var a = matrix<Float64>([
[1.0, 3.0],
[1.0, 4.0]
])
var expected_res = matrix<Float64>([
[ 0.75592895, 1.13389342],
[ 0.37796447, 1.88982237]
])
var res = sqrtm(a)
@Assert(approxEqual(res, expected_res, atol:1e-7))
}
@TestCase
func testLogm(): Unit {
var a = matrix<Float64>([
[1.0, 3.0],
[1.0, 4.0]
])
var expected_res = matrix<Float64>([
[-1.02571087, 2.05142174],
[ 0.68380725, 1.02571087]
])
var res = logm(a)
@Assert(approxEqual(res, expected_res, atol:1e-7))
}
}