atan(1.0+im) is incorrectly signed #31054
Comments
|
Note:
Lines 875 to 878 in a0bb006 we have
There must be something wrong with function Also: A simple way C# using this formula |
With the sign of the imaginary part of |
|
Quick test for some of the branch cuts: We maybe should have test-suites like this for more of the complex functions, possibly also including |
|
some history: this version (oldest one, 9 year ago) is right Line 220 in 44e0594 After pr #2891, next version is wrong `atanh_60d266`function atanh_60d266(z::Complex{T}) where T<:AbstractFloat
Ω=prevfloat(typemax(T))
θ=sqrt(Ω)/4
ρ=1/θ
x=real(z)
y=imag(z)
if x > θ || abs(y) > θ #Prevent overflow
return complex(copysign(pi/2, y), real(1/z))
elseif x==one(x)
ym=abs(y)+ρ
ξ=log(sqrt(sqrt(4+y^2))/sqrt(ym))
η=copysign(pi/2+atan(ym/2), y)/2
else #Normal case
ysq=(abs(y)+ρ)^2
ξ=log1p(4x/((1-x)^2 + ysq))/4
η=angle(complex(((1-x)*(1+x)-ysq)/2, y))
end
complex(ξ, η)
endand then all those later versions are wrong. pr #2891 have mentioned a paper: Branch Cuts for Complex Elementary Functions (or: Much Ado About Nothing's Sign Bit) We may need to take a closer look at this paper to fix this bug. |
|
good digging -- I remember that we had it right at some point. That is the seminal paper. |
|
Here is some function specific information on branch cuts generated from Maple. |
|
We definitely need to get this right. |
|
Wolfram Research has this well curated resource. use the I have added a pdf with Maple's plots of the principal branches here |
|
It seems like Line 955 in 795740a is incorrect. If x == -1 and y != 0, thenLines 960 to 962 in 795740a is executed, but that code is only valid for x == 1.(Reference: https://people.freebsd.org/~das/kahan86branch.pdf) |
I found this state of affairs obliquely. It is necessary to re-vet all the complex arc functions. |
|
Thanks for looking into it @cafaxo. Care to open a pull request? |
Patchchange those judgment statement did fix the Sign Test Resultfunction test_atanh()
float_data_list = [-1.0, 0.0, 1.0]
op_list = [identity, nextfloat, prevfloat]
for op in op_list
for real in float_data_list
for imag in float_data_list
z = op(real) + op(imag)*im
@show z
try
if tanh(atanh(z)) == z
nothing
elseif tanh(atanh(z)) ≈ z
nothing
else
@show tanh(atanh(z)) atanh(z)
end
catch ex
@show ex
end
end
end
end
endBUT it cause test about base/complex.jl | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/base/complex.jl b/base/complex.jl
index dc736ebb..6f0006dc 100644
--- a/base/complex.jl
+++ b/base/complex.jl
@@ -905,7 +905,7 @@ function atanh(z::Complex{T}) where T<:AbstractFloat
x, y = reim(z)
ax = abs(x)
ay = abs(y)
- if ax > θ || ay > θ #Prevent overflow
+ if x > θ || ay > θ #Prevent overflow
if isnan(y)
if isinf(x)
return Complex(copysign(zero(x),x), y)
@@ -917,7 +917,7 @@ function atanh(z::Complex{T}) where T<:AbstractFloat
return Complex(copysign(zero(x),x), copysign(oftype(y,pi)/2, y))
end
return Complex(real(1/z), copysign(oftype(y,pi)/2, y))
- elseif ax==1
+ elseif x == 1
if y == 0
ξ = copysign(oftype(x,Inf),x)
η = zero(y)After apply this patch, all those if-else statement are the same as in the paper.
Testinclude("..\\test\\complex.jl") ResultTODO
|
|
is there a ready explanation of why the returned expressions differ |
|
This 2017 paper gives values used for testing the implementation of complex elementary functions |
|
I suggest function atanh(z::Complex{T}) where T<:AbstractFloat
Ω = prevfloat(typemax(T))
θ = sqrt(Ω)/4
ρ = 1/θ
x, y = reim(z)
ax = abs(x)
ay = abs(y)
if ax > θ || ay > θ #Prevent overflow
if isnan(y)
if isinf(x)
return Complex(copysign(zero(x),x), y)
else
return Complex(real(1/z), y)
end
end
if isinf(y)
return Complex(copysign(zero(x),x), copysign(oftype(y,pi)/2, y))
end
return Complex(real(1/z), copysign(oftype(y,pi)/2, y))
end
β = copysign(1, x)
z = β * conj(z)
x, y = reim(z)
if x == 1
if y == 0
ξ = oftype(x, Inf)
η = y
else
ym = ay+ρ
ξ = log(sqrt(sqrt(4+y*y))/sqrt(ym))
η = copysign(oftype(y,pi)/2 + atan(ym/2), y)/2
end
else #Normal case
ysq = (ay+ρ)^2
if x == 0
ξ = x
else
ξ = log1p(4x/((1-x)^2 + ysq))/4
end
η = angle(Complex((1-x)*(1+x)-ysq, 2y))/2
end
return β * complex(ξ, -η)
endAs a patch@@ -952,10 +952,16 @@
return Complex(copysign(zero(x),x), copysign(oftype(y,pi)/2, y))
end
return Complex(real(1/z), copysign(oftype(y,pi)/2, y))
- elseif ax==1
+ end
+
+ β = copysign(1, x)
+ z = β * conj(z)
+ x, y = reim(z)
+
+ if x == 1
if y == 0
- ξ = copysign(oftype(x,Inf),x)
- η = zero(y)
+ ξ = oftype(x, Inf)
+ η = y
else
ym = ay+ρ
ξ = log(sqrt(sqrt(4+y*y))/sqrt(ym))
@@ -970,7 +976,8 @@
end
η = angle(Complex((1-x)*(1+x)-ysq, 2y))/2
end
- Complex(ξ, η)
+
+ return β * complex(ξ, -η)
end
atanh(z::Complex) = atanh(float(z))This fixes the sign issue, passes all tests in Base, passes test_atanh from @inkydragon, and does not throw on |
|
@cafaxo thank you for the diligence and attention |
|
I do not know why the paper conjugates z. This should probably be dropped from my suggestion. |
|
hold off on that
conj(complex) has IEEE754 specific details as I recall
check what conj does with re aor im parts +0, -0
Jeffrey Sarnoff
… On Feb 14, 2019, at 10:57 AM, cafaxo ***@***.***> wrote:
I do not know why the paper conjugates z. This should probably be dropped from my suggestion.
—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.
|
|
conj only negates the imaginary part: Line 259 in bd07ef4 I am pretty sure that we do not need that here. |
exerpted
|
|
Fixed by #31061 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
The sign of the imaginary part should be positive.
The text was updated successfully, but these errors were encountered: