Add Julia ports of hyperbolic (including inverse) functions from openlibm. #24159

pkofod · 2017-10-15T20:52:36Z

Still need to fix the license stuff and add benchmarks in BaseBenchmarks.

Edit: since I added some inverse hyperbolics as well, their tests are missing.

giordano · 2017-10-16T11:54:52Z

base/special/hyperbolic.jl

+    #    2. Find the the branch and the expression to calculate and return it
+    #      a)   x <= COSH_SMALL_X
+    #               return T(1)
+    #      b)   0 <= x <= ln2/2


Maybe COSH_SMALL_X <= x <= ln2/2? ;-)

That would make sense :)

giordano · 2017-10-16T11:56:12Z

If these functions are faster than the fastmath ones, you may also want to change them :-)

pkofod · 2017-10-16T12:14:22Z

If these functions are faster than the fastmath ones, you may also want to change them :-)

Yes I will! I know I missed it with some of the others. Let's benchmark first :)

musm · 2017-10-16T13:21:29Z

base/special/hyperbolic.jl

+# is preserved.
+# ====================================================
+
+_ldexp_exp(x::Float64, i::T) where T = ccall(("__ldexp_exp", libm), Float64, (Float64, T), x, i)


Just a quick comment, this should just use the julia ldexp function, if you don't need the error handling, you can write a _ldexp function and remove the error/exception handling in ldexp

Oh! I looked for a Julia binding/version.. but I must have mistyped the name when searching.

Is this worth changing now?

Unfortunately, this function is actually not in Julia yet, but yes, it is worth changing, and it should be a simple function to port.

edit: to clarify this requires porting 7-8 lines of code if I recall correctly. It should be in the k_exp.c file, and I might actually have port locally already.

edit2: yeah, I had them locally. I'll make a pr tomorrow.

is this not the function @musm was suggesting?

no, that computes x*2^n, what is needed here is exp(n)*2^n

edit: as I said, I'll make a pr tomorrow, as I'm off to bed now, but it goes something like the following, although it will not run because I have not included a few helper functions in the gist.

https://gist.github.com/pkofod/75a237b43cbcb82ea26410aecd8dc4ca

edit2: I know the doc string is wrong for the second function... I'm still mid-PR-preparation, but I'm too tired to continue today

KristofferC · 2017-10-17T10:50:58Z

base/special/hyperbolic.jl

 # acosh methods
-ACOSH_LN2(T) = 6.93147180559945286227e-01
-ACOSH_LN2(T) = 6.9314718246e-01
+AH_LN2(T::Type{Float64}) = 6.93147180559945286227e-01


No need for T

KristofferC · 2017-10-17T10:51:49Z

base/special/hyperbolic.jl

@@ -172,9 +176,41 @@ end
 tanh(x::Real) = tanh(float(x))

 # Inverse hyperbolic functions
+# asinh methods
+AH_LN2(T) =  6.93147180559945286227e-01


Can remove?

stevengj · 2017-10-17T18:02:49Z

base/special/hyperbolic.jl

+
+H_LARGE_X(::Type{Float32}) = 88.72283f0
+H_OVERFLOW_X(::Type{Float32}) = 89.415985f0
+function sinh(x::T) where T <: Union{Float32, Float64}


~~Wouldn't it be cleaner to use the syntax sinh(x::Union{Float32, Float64}) for functions like this?~~ Oh, nevermind, you use T in the function body.

yes I use these to pick the branching cutoffs. Do you think it would be cleaner to write sinh(x::Union{Float32, Float64}) and then just define H_LARGE_X(::Float32) = 88.72283f0? on the actual x instead of the type of x? Should compile down to the same thing, right?

I think it's fine as-is.

stevengj · 2017-10-17T18:06:11Z

test/math.jl

+        by = tanh(big(x))
+        @test Float32((tanh(x) - by))/eps(abs(Float32(by))) <= 1
+        bym = tanh(big(-x))
+        @test Float32(abs(tanh(-x) - bym))/eps(abs(Float32(bym))) <= 1


No tests for e.g. acosh?

In general, it would be good to make sure you have test coverage for all of the branches.

Maybe I should have put a wip. All these will get tests and benchmarks, sorry for not making that more clear.

The inverse functions were not there originally, but I decided to add them here and just haven't gotten the tests in yet.

pkofod · 2017-11-23T21:21:28Z

I'm not using the __ldexp_exp from openlibm again, I will implement these in another PR.

All functions now have tests. Need benchmarks and license-thing.

pkofod · 2017-11-24T08:48:51Z

benchmark PR JuliaCI/BaseBenchmarks.jl#141

StefanKarpinski · 2018-02-09T17:39:52Z

If either @simonbyrne or @stevengj have a look at this and approve, it's good to squash merge then!

stevengj · 2018-02-12T21:53:38Z

test/math.jl

+    end
+    for x in (pcnfloat(2.7755602085408512e-17)..., pcnfloat(22.0)..., pcnfloat(709.7822265633563)...)
+        by = cosh(big(x))
+        @test Float64((cosh(x) - by))/eps(abs(Float64(by))) <= 1


Would be clearer to use @test cosh(x) ≈ cosh(big(y)) rtol=eps(Float64) here and similarly in the other tests.

big(y) is meant to be big(x), right?

stevengj · 2018-02-12T21:56:49Z

test/math.jl

+        by = acosh(big(x))
+        @test Float64((acosh(x) - by))/eps(abs(Float64(by))) <= 1
+    end
+    for x in (nextfloat(1f0), pcnfloat(2f0)..., pcnfloat(2f0^28)...)


Couldn't you merge these two loops by putting them into the for T loop as

for x in (nextfloat(T(1)), pcnfloat(T(2))..., pcnfloat(T(2^28))...) @test acosh(x) ≈ acosh(big(x) rtol=eps(T) end

stevengj · 2018-02-12T21:58:43Z

test/math.jl

+        by = tanh(big(x))
+        @test Float64((tanh(x) - by))/eps(abs(Float64(by))) <= 1
+        bym = tanh(big(-x))
+        @test Float64(abs(tanh(-x) - bym))/eps(abs(Float64(bym))) <= 1


@test tanh(x) ≈ tanh(big(x)) ≈ -tanh(x) rtol=eps(Float64)

You could merge the for x loops into

for x in Iterators.flatten(pcnfloat.(T == Float64 ? [2.0^-28,1.0,22.0] : [2f0^-12,1f0,9f0])) @test tanh(x) ≈ tanh(big(x)) ≈ -tanh(x) rtol=eps(T) end

stevengj · 2018-02-12T22:07:10Z

base/special/hyperbolic.jl

+acosh(x::Real) = acosh(float(x))
+
+# atanh methods
+@noinline atanh_domain_error(x) = throw(DomainError(x, "atanh(x) is only defined for x <= 1."))


for |x| ≤ 1

stevengj · 2018-02-12T22:08:22Z

base/special/hyperbolic.jl

+asinh(x::Real) = asinh(float(x))
+
+# acosh methods
+@noinline acosh_domain_error(x) = throw(DomainError(x, "acosh(x) is only defined for x >= 1."))


Unicode ≥ is a bit nicer

pkofod · 2018-02-12T22:58:57Z

Thanks @stevengj ! I hope I've addressed your comments.

stevengj · 2018-02-13T00:05:16Z

test/math.jl

+        @test sinh(-T(1000)) === -T(Inf)
+        @test isnan_type(T, sinh(T(NaN)))
+    end
+    for x in (pcnfloat(2.0^-28)..., pcnfloat(22.0)..., pcnfloat(709.7822265633563)...)


Seems like it would be better to use SINH_SMALL_X(T), H_MEDIUM_X(T), etcetera (importing those symbols into the test). That way it will unify the tests for both T, make it clearer where these thresholds come from, and ensure the tests are updated if the thresholds change (or if the sinh code is rewritten).

stevengj · 2018-02-13T00:07:26Z

test/math.jl

+        @test isnan_type(T, sinh(T(NaN)))
+    end
+    for x in (pcnfloat(2.0^-28)..., pcnfloat(22.0)..., pcnfloat(709.7822265633563)...)
+        @test sinh(x) ≈ sinh(big(x)) rtol=eps(Float64)


Seems simpler to do @test sinh(x) ≈ sinh(big(x)) ≈ -sinh(x) rtol=eps(T), but it looks like this syntax isn't supported by @test. 😢

Yeah, you suggested that above and I tried it but without luck so I just split it out into two.

stevengj · 2018-02-13T00:16:20Z

test/math.jl

+        @test asinh(nextfloat(zero(T))) === nextfloat(zero(T))
+        @test asinh(prevfloat(zero(T))) === prevfloat(zero(T))
+        @test isnan_type(T, asinh(T(NaN)))
+        for x in Iterators.flatten(pcnfloat.(T == Float64 ? [2.0^-28,2.0,2.0^28] : [2f0^-12,2f0,22f28]))


I'm confused about this conditional. Looking back at the asinh code, it seems like it uses 2^-28 and 2^28 for both single and double precision. 2^-12 is only used for tanh.

Again, this would be clarified by using TANH_SMALL_X(T) etcetera here for precision-dependent thresholds.

This is a good idea. Originally it was unsure if I there were going to be tests in separate commits before the actual change (just to make sure nothing changed), but I guess this will be squashed anyway, so might as well just use the thresholds from the actual code.

stevengj · 2018-02-13T16:22:46Z

test/math.jl

@@ -907,7 +901,7 @@ end
        @test asinh(nextfloat(zero(T))) === nextfloat(zero(T))
        @test asinh(prevfloat(zero(T))) === prevfloat(zero(T))
        @test isnan_type(T, asinh(T(NaN)))
-        for x in Iterators.flatten(pcnfloat.(T == Float64 ? [2.0^-28,2.0,2.0^28] : [2f0^-12,2f0,22f28]))
+        for x in Iterators.flatten(pcnfloat.([T(2)^-28,T(2),T(2)]))


The second T(2) doesn't look right. I assume you mean T(2)^28?

stevengj · 2018-02-13T16:24:08Z

test/math.jl

@@ -919,7 +913,7 @@ end
        @test_throws DomainError acosh(T(0.1))
        @test acosh(one(T)) === zero(T)
        @test isnan_type(T, acosh(T(NaN)))
-        for x in (nextfloat(T(1)), pcnfloat(T(2))..., pcnfloat(T(2^28))...)
+        for x in Iterators.flatten(pcnfloat.([T(1) T(2), T(2^28)]))


Shouldn't it be Iterators.flatten([nextfloat(T(1)); pcnfloat.([T(2), T(2^28)])])? Or [nextfloat(T(1)), pcnfloat(T(2))..., pcnfloat(T(2^28))...]

I seem to have pushed without this fix.

pkofod · 2018-02-13T19:13:46Z

stevengj approved these changes an hour ago

Thanks for great comments @stevengj ! The tests are so much cleaner now.

…libm. (JuliaLang/julia#24159)

…libm. (#24159)

Add Julia ports of hyperbolic functions from openlibm.

fd37f3d

ararslan added maths Mathematical functions potential benchmark Could make a good benchmark in BaseBenchmarks labels Oct 15, 2017

Actually save the correct tests...

0bead47

giordano reviewed Oct 16, 2017

View reviewed changes

musm reviewed Oct 16, 2017

View reviewed changes

pkofod added 4 commits October 16, 2017 22:52

Fix comments in hyperbolic functions and add acosh.

53886c5

Use Julia version of ldexp instead of openlibm.

9cba2fa

Add asinh

423d46a

Fix missing end and improve comments.

dd8c892

KristofferC reviewed Oct 17, 2017

View reviewed changes

pkofod added 2 commits October 17, 2017 13:00

Spelling.

4daaa8b

Changes

6bd838c

stevengj reviewed Oct 17, 2017

View reviewed changes

Add atanh and minor fixes in general.

5f8f975

pkofod changed the title ~~Add Julia ports of hyperbolic functions from openlibm.~~ Add Julia ports of hyperbolic (including inverse) functions from openlibm. Oct 17, 2017

pkofod added 5 commits October 17, 2017 22:46

Remove atanh from loop in math.jl.

ecbc401

_ldexp_exp to ldexp

2aec485

Readd ldexp_exp

7bf4cdf

Remove semi-colon.

90b5c45

Add tests for inverse hyperbolics.

ea5c860

pkofod added 3 commits November 23, 2017 22:28

Fix licences.

79d7477

Fix a missing f0.

ee9b458

Clean up tested values.

1fbca61

StefanKarpinski requested review from simonbyrne and stevengj February 9, 2018 17:39

Merge branch 'master' into hyperbolic

2b396ad

stevengj reviewed Feb 12, 2018

View reviewed changes

pkofod added 3 commits February 12, 2018 23:27

Improve domain error messages.

c15347b

Make tests neater.

79ae7d1

Merge branch 'hyperbolic' of github.com:pkofod/julia into hyperbolic

eb1390a

stevengj reviewed Feb 13, 2018

View reviewed changes

stevengj mentioned this pull request Feb 13, 2018

@test a ≈ b atol=ε #19901

Merged

stevengj reviewed Feb 13, 2018

View reviewed changes

pkofod added 3 commits February 13, 2018 15:31

Fix wrong merge conflict resolution, and make tests prettier.

c3bedbd

More testing.

d92b7bc

Fix end

b751c37

stevengj reviewed Feb 13, 2018

View reviewed changes

More test changes.

5c22157

stevengj approved these changes Feb 13, 2018

View reviewed changes

StefanKarpinski merged commit 76e5232 into JuliaLang:master Feb 14, 2018

pkofod deleted the hyperbolic branch February 14, 2018 20:16

KristofferC removed the potential benchmark Could make a good benchmark in BaseBenchmarks label Oct 8, 2018

vchuravy pushed a commit to JuliaPackaging/LazyArtifacts.jl that referenced this pull request Oct 2, 2023

Add Julia ports of hyperbolic (including inverse) functions from open…

8f87f03

…libm. (JuliaLang/julia#24159)

Keno pushed a commit that referenced this pull request Jun 5, 2024

Add Julia ports of hyperbolic (including inverse) functions from open…

560bd96

…libm. (#24159)

Add Julia ports of hyperbolic (including inverse) functions from openlibm. #24159

Add Julia ports of hyperbolic (including inverse) functions from openlibm. #24159

Conversation

pkofod commented Oct 15, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

giordano commented Oct 16, 2017

pkofod commented Oct 16, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pkofod Feb 14, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pkofod Feb 14, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Oct 17, 2017 • edited Loading

Choose a reason for hiding this comment

pkofod Oct 17, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pkofod commented Nov 23, 2017

pkofod commented Nov 24, 2017

StefanKarpinski commented Feb 9, 2018

stevengj Feb 12, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pkofod commented Feb 12, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stevengj Feb 13, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pkofod commented Feb 13, 2018

pkofod commented Oct 15, 2017 •

edited

Loading

pkofod Feb 14, 2018 •

edited

Loading

pkofod Feb 14, 2018 •

edited

Loading

stevengj Oct 17, 2017 •

edited

Loading

pkofod Oct 17, 2017 •

edited

Loading

stevengj Feb 12, 2018 •

edited

Loading

stevengj Feb 13, 2018 •

edited

Loading