Remove ieee754_rem_pio2 in favor of a rem_pio2_kernel written in Julia.

LICENSE.md

@@ -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]
 - [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]

base/math.jl

@@ -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)
@@ -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)
@@ -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
@@ -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
@@ -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"))