Correct lyap definition for complex scalar arguments (#39901)

Author: Regan Koopmans

stdlib/LinearAlgebra/src/dense.jl

@@ -1568,4 +1568,4 @@ function lyap(A::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:BlasFloat}
     rmul!(Q*(Y * adjoint(Q)), inv(scale))
 end
 lyap(A::StridedMatrix{T}, C::StridedMatrix{T}) where {T<:Integer} = lyap(float(A), float(C))
-lyap(a::T, c::T) where {T<:Number} = -c/(2a)
+lyap(a::Union{Real,Complex}, c::Union{Real,Complex}) = -c/(2real(a))

stdlib/LinearAlgebra/test/dense.jl

@@ -767,6 +767,10 @@ end
     @test diag(zeros(0,1),2) == []
 end
 
+@testset "issue #39857" begin
+    @test lyap(1.0+2.0im, 3.0+4.0im) == -1.5 - 2.0im
+end
+
 @testset "Matrix to real power" for elty in (Float64, ComplexF64)
 # Tests proposed at Higham, Deadman: Testing Matrix Function Algorithms Using Identities, March 2014
     #Aa : only positive real eigenvalues