Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove ieee754_rem_pio2 in favor of a rem_pio2_kernel written in Julia. #22603

Merged
merged 39 commits into from
Aug 2, 2017
Merged
Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
39 commits
Select commit Hold shift + click to select a range
3f951d1
Remove ieee754_rem_pio2 in favor of a rem_pio2_kernel written in Julia.
pkofod Jun 28, 2017
62a91ed
Add missing begin key.
pkofod Jun 29, 2017
16ec5f0
Remove _approx.
pkofod Jun 29, 2017
6f8e725
Move to separate files.
pkofod Jun 29, 2017
c6f593e
Fix LICENSE.md to mention FDLIBM instead of Openlibm.
pkofod Jun 29, 2017
5fd35f4
Address comments.
pkofod Jun 30, 2017
1d2faa6
Strengthen test to faithfully rounded.
pkofod Jun 30, 2017
7878d2c
Fix LICENSE.md message for rem_pio2.
pkofod Jun 30, 2017
0017edd
Fix style in LICENSE.md entry.
pkofod Jun 30, 2017
8258c48
Remove semicolons.
pkofod Jun 30, 2017
83da1a2
Move highword up, and remove duplicate unsafe_trunc.
pkofod Jun 30, 2017
27455bf
Fix LICENSE.md by removing a bullet and changing license of base/spec…
pkofod Jun 30, 2017
bb626b3
Change license info for base/special/exp.jl.
pkofod Jul 1, 2017
45d4ba1
Small changes.
pkofod Jul 1, 2017
ff19e40
Get and reset precision for BigFloats, and space before rem in -rem.
pkofod Jul 1, 2017
375b8af
setprecision do
pkofod Jul 1, 2017
8ca63c1
Add comments, move test, and switch to muladd in some places.
pkofod Jul 5, 2017
1ef4918
Fix y1 branches of rem2pi.
pkofod Jul 5, 2017
a1f1c34
Small changes.
pkofod Jul 5, 2017
ecc6963
Move comment in rem_pio2.jl and add test for fast branch of mod2pi.
pkofod Jul 6, 2017
f8090f5
rint docstring fix and make it clear what the constant is.
pkofod Jul 6, 2017
eed7e62
Update comment for INV2PI.
pkofod Jul 6, 2017
dd2c152
Fix wrong test set name.
pkofod Jul 6, 2017
c06e2f3
Tests against ieee754_rem_pio2 output.
pkofod Jul 8, 2017
6a5dc4a
Inline cody_waite functions.
pkofod Jul 8, 2017
50e52e3
rint -> round, remove rint, remove one argument cody waite, replace I…
pkofod Jul 9, 2017
7b6213e
Add some tests.
pkofod Jul 9, 2017
ae8577c
Inline rem_pio2_kernel, and rearrange code slightly.
pkofod Jul 13, 2017
ba62caf
fix xhp
pkofod Jul 14, 2017
f209c0a
Use DoubleFloat64.
pkofod Jul 21, 2017
994e158
Move constants into functions.
pkofod Jul 22, 2017
aba4579
Fix escaping of mod
pkofod Jul 22, 2017
4865a41
Fix tests and remove specific variables.
pkofod Jul 22, 2017
6dda42c
Fix tests
pkofod Jul 23, 2017
8adb120
Fix issues raised in comments.
pkofod Jul 26, 2017
0545479
More appropriate ulp test (test against eps of reference number).
pkofod Jul 26, 2017
8ee70a1
Merge branch 'master' into rempio2gsoc
pkofod Jul 30, 2017
765c566
Change link to a stable github link.
pkofod Jul 31, 2017
32d2839
Merge remote-tracking branch 'pkofod/rempio2gsoc' into rempio2gsoc
pkofod Jul 31, 2017
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions LICENSE.md
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@ Julia includes code from the following projects, which have their own licenses:
- [LLVM](http://releases.llvm.org/3.9.0/LICENSE.TXT) (for parts of src/jitlayers.cpp and src/disasm.cpp) [BSD-3, effectively]
- [MUSL](http://git.musl-libc.org/cgit/musl/tree/COPYRIGHT) (for getopt implementation on Windows) [MIT]
- [MINGW](https://sourceforge.net/p/mingw/mingw-org-wsl/ci/legacy/tree/mingwrt/mingwex/dirname.c) (for dirname implementation on Windows) [MIT]
- [FDLIBM](http://www.netlib.org/fdlibm/readme) [Freely distributable]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you look down this file I think you also need to update and add a line below

base/special/exp.jl (see FREEBSD MSUN [FreeBSD/2-clause BSD/Simplified BSD License])

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, alright, thanks.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the difference between this section and the one below? It's not obvious to me.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

projects we use vs how we use them? I think this section might actually be for core language/compiler, not base

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, yes, it's not very clear. It looks like this section only applies to core features, in which case FDLIB shouldn't be added here (just like MSUN is listed below, but not here).

BTW, if it's kept here, the new line should mention what it's used for, like existing lines.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this section is the right place for non-core julia code under base

- [NetBSD](http://www.netbsd.org/about/redistribution.html) (for setjmp, longjmp, and strptime implementations on Windows) [BSD-3]
- [Python](https://docs.python.org/2/license.html) (for strtod implementation on Windows) [BSD-3, effectively]
- [randmtzig.c](https://github.com/JuliaLang/julia/blob/master/test/perf/micro/randmtzig.c) for Gaussian random number generation (for C benchmarks only) [BSD-3]
Expand Down Expand Up @@ -74,6 +75,7 @@ The following components of Julia's standard library have separate licenses:
- base/sparse/umfpack.jl (see [SUITESPARSE](http://faculty.cse.tamu.edu/davis/suitesparse.html))
- base/sparse/cholmod.jl (see [SUITESPARSE](http://faculty.cse.tamu.edu/davis/suitesparse.html))
- base/special/exp.jl (see [FREEBSD MSUN](https://github.com/freebsd/freebsd) [FreeBSD/2-clause BSD/Simplified BSD License])
- base/special/rem_pio2.jl [FDLIBM](http://www.netlib.org/fdlibm/readme) [Freely distributable]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should point to a copy of the license, as "freely distributable" is a bit vague (and in particular it doesn't say that the copyright attribution should be preserved when distributing).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can point to http://www.netlib.org/fdlibm/e_rem_pio2.c ? There isn't really a license file in fdlibm as far as I can see. I could write "Freely distributable with preserved copyright notice."?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doing both sounds like a good idea to me.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The license is the same as the exp code line

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why? I mean where is it stated?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I know next to nothing about licensing, but is it necessary to go with the msun license if they simply took it from fdlibm?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hmmm, good point, maybe someone more knowledgeable here can comment, if you do end up changing it best to also modify the one for exp.jl since they should be the same

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually what you have probably makes the most sense



Julia's build process uses the following external tools:
Expand Down
76 changes: 29 additions & 47 deletions base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -738,25 +738,6 @@ function add22condh(xh::Float64, xl::Float64, yh::Float64, yl::Float64)
return zh
end

function ieee754_rem_pio2(x::Float64)
# rem_pio2 essentially computes x mod pi/2 (ie within a quarter circle)
# and returns the result as
# y between + and - pi/4 (for maximal accuracy (as the sign bit is exploited)), and
# n, where n specifies the integer part of the division, or, at any rate,
# in which quadrant we are.
# The invariant fulfilled by the returned values seems to be
# x = y + n*pi/2 (where y = y1+y2 is a double-double and y2 is the "tail" of y).
# Note: for very large x (thus n), the invariant might hold only modulo 2pi
# (in other words, n might be off by a multiple of 4, or a multiple of 100)

# this is just wrapping up
# https://github.com/JuliaLang/openspecfun/blob/master/rem_pio2/e_rem_pio2.c

y = [0.0,0.0]
n = ccall((:__ieee754_rem_pio2, openspecfun), Cint, (Float64,Ptr{Float64}), x, y)
return (n,y)
end

# multiples of pi/2, as double-double (ie with "tail")
const pi1o2_h = 1.5707963267948966 # convert(Float64, pi * BigFloat(1/2))
const pi1o2_l = 6.123233995736766e-17 # convert(Float64, pi * BigFloat(1/2) - pi1o2_h)
Expand Down Expand Up @@ -804,47 +785,47 @@ function rem2pi end
function rem2pi(x::Float64, ::RoundingMode{:Nearest})
abs(x) < pi && return x

(n,y) = ieee754_rem_pio2(x)
n,y1,y2 = rem_pio2_kernel(x)

if iseven(n)
if n & 2 == 2 # n % 4 == 2: add/subtract pi
if y[1] <= 0
return add22condh(y[1],y[2],pi2o2_h,pi2o2_l)
if y1 <= 0
return add22condh(y1,y2,pi2o2_h,pi2o2_l)
else
return add22condh(y[1],y[2],-pi2o2_h,-pi2o2_l)
return add22condh(y1,y2,-pi2o2_h,-pi2o2_l)
end
else # n % 4 == 0: add 0
return y[1]
return y1
end
else
if n & 2 == 2 # n % 4 == 3: subtract pi/2
return add22condh(y[1],y[2],-pi1o2_h,-pi1o2_l)
return add22condh(y1,y2,-pi1o2_h,-pi1o2_l)
else # n % 4 == 1: add pi/2
return add22condh(y[1],y[2],pi1o2_h,pi1o2_l)
return add22condh(y1,y2,pi1o2_h,pi1o2_l)
end
end
end
function rem2pi(x::Float64, ::RoundingMode{:ToZero})
ax = abs(x)
ax <= 2*Float64(pi,RoundDown) && return x

(n,y) = ieee754_rem_pio2(ax)
n,y1,y2 = rem_pio2_kernel(ax)

if iseven(n)
if n & 2 == 2 # n % 4 == 2: add pi
z = add22condh(y[1],y[2],pi2o2_h,pi2o2_l)
z = add22condh(y1,y2,pi2o2_h,pi2o2_l)
else # n % 4 == 0: add 0 or 2pi
if y[1] > 0
z = y[1]
if y1 > 0
z = y1
else # negative: add 2pi
z = add22condh(y[1],y[2],pi4o2_h,pi4o2_l)
z = add22condh(y1,y2,pi4o2_h,pi4o2_l)
end
end
else
if n & 2 == 2 # n % 4 == 3: add 3pi/2
z = add22condh(y[1],y[2],pi3o2_h,pi3o2_l)
z = add22condh(y1,y2,pi3o2_h,pi3o2_l)
else # n % 4 == 1: add pi/2
z = add22condh(y[1],y[2],pi1o2_h,pi1o2_l)
z = add22condh(y1,y2,pi1o2_h,pi1o2_l)
end
end
copysign(z,x)
Expand All @@ -858,23 +839,23 @@ function rem2pi(x::Float64, ::RoundingMode{:Down})
end
end

(n,y) = ieee754_rem_pio2(x)
n,y1,y2 = rem_pio2_kernel(x)

if iseven(n)
if n & 2 == 2 # n % 4 == 2: add pi
return add22condh(y[1],y[2],pi2o2_h,pi2o2_l)
return add22condh(y1,y2,pi2o2_h,pi2o2_l)
else # n % 4 == 0: add 0 or 2pi
if y[1] > 0
return y[1]
if y1 > 0
return y1
else # negative: add 2pi
return add22condh(y[1],y[2],pi4o2_h,pi4o2_l)
return add22condh(y1,y2,pi4o2_h,pi4o2_l)
end
end
else
if n & 2 == 2 # n % 4 == 3: add 3pi/2
return add22condh(y[1],y[2],pi3o2_h,pi3o2_l)
return add22condh(y1,y2,pi3o2_h,pi3o2_l)
else # n % 4 == 1: add pi/2
return add22condh(y[1],y[2],pi1o2_h,pi1o2_l)
return add22condh(y1,y2,pi1o2_h,pi1o2_l)
end
end
end
Expand All @@ -887,23 +868,23 @@ function rem2pi(x::Float64, ::RoundingMode{:Up})
end
end

(n,y) = ieee754_rem_pio2(x)
n,y1,y2 = rem_pio2_kernel(x)

if iseven(n)
if n & 2 == 2 # n % 4 == 2: sub pi
return add22condh(y[1],y[2],-pi2o2_h,-pi2o2_l)
return add22condh(y1,y2,-pi2o2_h,-pi2o2_l)
else # n % 4 == 0: sub 0 or 2pi
if y[1] < 0
return y[1]
if y1 < 0
return y1
else # positive: sub 2pi
return add22condh(y[1],y[2],-pi4o2_h,-pi4o2_l)
return add22condh(y1,y2,-pi4o2_h,-pi4o2_l)
end
end
else
if n & 2 == 2 # n % 4 == 3: sub pi/2
return add22condh(y[1],y[2],-pi1o2_h,-pi1o2_l)
return add22condh(y1,y2,-pi1o2_h,-pi1o2_l)
else # n % 4 == 1: sub 3pi/2
return add22condh(y[1],y[2],-pi3o2_h,-pi3o2_l)
return add22condh(y1,y2,-pi3o2_h,-pi3o2_l)
end
end
end
Expand Down Expand Up @@ -977,6 +958,7 @@ include(joinpath("special", "exp.jl"))
include(joinpath("special", "exp10.jl"))
include(joinpath("special", "trig.jl"))
include(joinpath("special", "gamma.jl"))
include(joinpath("special", "rem_pio2.jl"))

module JuliaLibm
include(joinpath("special", "log.jl"))
Expand Down
Loading