patches to matrix log (#33245)

* patches to matrix log

Avoid integer overflow if `s > 63`.
Correct logic for `s == 0`.
Only use fancy divided difference formulae if eigenvalues are close - avoids dangerous roundoff error if they are in opposite sectors.

* add tests

Author: RalphAS

stdlib/LinearAlgebra/src/triangular.jl

@@ -2183,13 +2183,14 @@ function log(A0::UpperTriangular{T}) where T<:BlasFloat
     end
 
     # Compute accurate superdiagonal of T
-    blockpower!(A, A0, 1 / 2^s)
+    blockpower!(A, A0, 0.5^s)
 
     # Compute accurate diagonal of T
     for i = 1:n
         a = A0[i,i]
         if s == 0
-            r = a - 1
+            A[i,i] = a - 1
+            continue
         end
         s0 = s
         if angle(a) >= pi / 2
@@ -2226,7 +2227,7 @@ function log(A0::UpperTriangular{T}) where T<:BlasFloat
     end
 
     # Scale back
-    lmul!(2^s, Y)
+    lmul!(2.0^s, Y)
 
     # Compute accurate diagonal and superdiagonal of log(T)
     for k = 1:n-1
@@ -2242,8 +2243,12 @@ function log(A0::UpperTriangular{T}) where T<:BlasFloat
             Y[k,k+1] = A0[k,k+1] * (logAkp1 - logAk) / (Akp1 - Ak)
         else
             z = (Akp1 - Ak)/(Akp1 + Ak)
-            w = atanh(z) + im * pi * (unw(logAkp1-logAk) + unw(log1p(z)-log1p(-z)))
-            Y[k,k+1] = 2 * A0[k,k+1] * w / (Akp1 - Ak)
+            if abs(z) > 1
+                Y[k,k+1] = A0[k,k+1] * (logAkp1 - logAk) / (Akp1 - Ak)
+            else
+                w = atanh(z) + im * pi * (unw(logAkp1-logAk) - unw(log1p(z)-log1p(-z)))
+                Y[k,k+1] = 2 * A0[k,k+1] * w / (Akp1 - Ak)
+            end
         end
     end
 
@@ -2396,9 +2401,13 @@ function blockpower!(A::UpperTriangular, A0::UpperTriangular, p)
             logAk = log(Ak)
             logAkp1 = log(Akp1)
             z = (Akp1 - Ak)/(Akp1 + Ak)
-            w = atanh(z) + im * pi * (unw(logAkp1-logAk) + unw(log1p(z)-log1p(-z)))
-            dd = 2 * exp(p*(logAk+logAkp1)/2) * sinh(p*w) / (Akp1 - Ak);
-            A[k,k+1] = A0[k,k+1] * dd
+            if abs(z) > 1
+                A[k,k+1] = A0[k,k+1] * (Akp1p - Akp) / (Akp1 - Ak)
+            else
+                w = atanh(z) + im * pi * (unw(logAkp1-logAk) - unw(log1p(z)-log1p(-z)))
+                dd = 2 * exp(p*(logAk+logAkp1)/2) * sinh(p*w) / (Akp1 - Ak);
+                A[k,k+1] = A0[k,k+1] * dd
+            end
         end
     end
 end

stdlib/LinearAlgebra/test/dense.jl

@@ -909,4 +909,13 @@ end
     end
 end
 
+@testset "Matrix log PR #33245" begin
+    # edge case for divided difference
+    A1 = triu(ones(3,3),1) + diagm([1.0, -2eps()-1im, -eps()+0.75im])
+    @test exp(log(A1)) ≈ A1
+    # case where no sqrt is needed (s=0)
+    A2 = [1.01 0.01 0.01; 0 1.01 0.01; 0 0 1.01]
+    @test exp(log(A2)) ≈ A2
+end
+
 end # module TestDense