Remove ieee754_rem_pio2 in favor of a rem_pio2_kernel written in Juli…

…a. (#22603)

* Remove ieee754_rem_pio2 in favor of a rem_pio2_kernel written in Julia.

* Add missing begin key.

* Remove _approx.

* Move to separate files.

* Fix LICENSE.md to mention FDLIBM instead of Openlibm.

* Address comments.

* Strengthen test to faithfully rounded.

* Fix LICENSE.md message for rem_pio2.

* Fix style in LICENSE.md entry.

* Remove semicolons.

* Move highword up, and remove duplicate unsafe_trunc.

* Fix LICENSE.md by removing a bullet and changing license of base/special/exp.jl.

* Change license info for base/special/exp.jl.

* Small changes.

* Get and reset precision for BigFloats, and space before rem in -rem.

* setprecision do

* Add comments, move test, and switch to muladd in some places.

* Fix y1 branches of rem2pi.

* Small changes.

* Move comment in rem_pio2.jl and add test for fast branch of mod2pi.

* rint docstring fix and make it clear what the constant is.

* Update comment for INV2PI.

* Fix wrong test set name.

* Tests against ieee754_rem_pio2 output.

* Inline cody_waite functions.

* rint -> round, remove rint, remove one argument cody waite, replace Int(x) with trunc(Int, x).

* Add some tests.

* Inline rem_pio2_kernel, and rearrange code slightly.

* fix xhp

* Use DoubleFloat64.

* Move constants into functions.

* Fix escaping of mod

* Fix tests and remove specific variables.

* Fix tests

* Fix issues raised in comments.

* More appropriate ulp test (test against eps of reference number).

* Change link to a stable github link.

LICENSE.md

@@ -73,7 +73,8 @@ The following components of Julia's standard library have separate licenses:
 - base/grisu/* (see [double-conversion](https://github.com/google/double-conversion/blob/master/LICENSE))
 - 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/exp.jl (see [FDLIBM](http://www.netlib.org/fdlibm/e_exp.c) [Freely distributable with preserved copyright notice])
+- base/special/rem_pio2.jl (see [FDLIBM](http://www.netlib.org/fdlibm/e_rem_pio2.c) [Freely distributable with preserved copyright notice])
 
 
 Julia's build process uses the following external tools:

base/math.jl

@@ -781,25 +781,6 @@ function add22condh(xh::Float64, xl::Float64, yh::Float64, yl::Float64)
     return zh
 end
 
-@inline 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 = Ref{NTuple{2,Float64}}()
-    n = ccall((:__ieee754_rem_pio2, openspecfun), Cint, (Float64, Ptr{Void}), 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)
@@ -847,47 +828,47 @@ function rem2pi end
 function rem2pi(x::Float64, ::RoundingMode{:Nearest})
     abs(x) < pi && return x
 
-    (n,y) = ieee754_rem_pio2(x)
+    n,y = 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 y.hi <= 0
+                return add22condh(y.hi,y.lo,pi2o2_h,pi2o2_l)
             else
-                return add22condh(y[1],y[2],-pi2o2_h,-pi2o2_l)
+                return add22condh(y.hi,y.lo,-pi2o2_h,-pi2o2_l)
             end
         else          # n % 4 == 0: add 0
-            return y[1]
+            return y.hi+y.lo
         end
     else
         if n & 2 == 2 # n % 4 == 3: subtract pi/2
-            return add22condh(y[1],y[2],-pi1o2_h,-pi1o2_l)
+            return add22condh(y.hi,y.lo,-pi1o2_h,-pi1o2_l)
         else          # n % 4 == 1: add pi/2
-            return add22condh(y[1],y[2],pi1o2_h,pi1o2_l)
+            return add22condh(y.hi,y.lo,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,y = rem_pio2_kernel(x)
 
     if iseven(n)
         if n & 2 == 2 # n % 4 == 2: add pi
-            z = add22condh(y[1],y[2],pi2o2_h,pi2o2_l)
+            z = add22condh(y.hi,y.lo,pi2o2_h,pi2o2_l)
         else          # n % 4 == 0: add 0 or 2pi
-            if y[1] > 0
-                z = y[1]
+            if y.hi > 0
+                z = y.hi+y.lo
             else      # negative: add 2pi
-                z = add22condh(y[1],y[2],pi4o2_h,pi4o2_l)
+                z = add22condh(y.hi,y.lo,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(y.hi,y.lo,pi3o2_h,pi3o2_l)
         else          # n % 4 == 1: add pi/2
-            z = add22condh(y[1],y[2],pi1o2_h,pi1o2_l)
+            z = add22condh(y.hi,y.lo,pi1o2_h,pi1o2_l)
         end
     end
     copysign(z,x)
@@ -901,23 +882,23 @@ function rem2pi(x::Float64, ::RoundingMode{:Down})
         end
     end
 
-    (n,y) = ieee754_rem_pio2(x)
+    n,y = 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(y.hi,y.lo,pi2o2_h,pi2o2_l)
         else          # n % 4 == 0: add 0 or 2pi
-            if y[1] > 0
-                return y[1]
+            if y.hi > 0
+                return y.hi+y.lo
             else      # negative: add 2pi
-                return add22condh(y[1],y[2],pi4o2_h,pi4o2_l)
+                return add22condh(y.hi,y.lo,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(y.hi,y.lo,pi3o2_h,pi3o2_l)
         else          # n % 4 == 1: add pi/2
-            return add22condh(y[1],y[2],pi1o2_h,pi1o2_l)
+            return add22condh(y.hi,y.lo,pi1o2_h,pi1o2_l)
         end
     end
 end
@@ -930,23 +911,23 @@ function rem2pi(x::Float64, ::RoundingMode{:Up})
         end
     end
 
-    (n,y) = ieee754_rem_pio2(x)
+    n,y = 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(y.hi,y.lo,-pi2o2_h,-pi2o2_l)
         else          # n % 4 == 0: sub 0 or 2pi
-            if y[1] < 0
-                return y[1]
+            if y.hi < 0
+                return y.hi+y.lo
             else      # positive: sub 2pi
-                return add22condh(y[1],y[2],-pi4o2_h,-pi4o2_l)
+                return add22condh(y.hi,y.lo,-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(y.hi,y.lo,-pi1o2_h,-pi1o2_l)
         else          # n % 4 == 1: sub 3pi/2
-            return add22condh(y[1],y[2],-pi3o2_h,-pi3o2_l)
+            return add22condh(y.hi,y.lo,-pi3o2_h,-pi3o2_l)
         end
     end
 end
@@ -1023,6 +1004,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"))

base/special/exp.jl

@@ -1,4 +1,4 @@
-# Based on FreeBSD lib/msun/src/e_exp.c
+# Based on FDLIBM http://www.netlib.org/fdlibm/e_exp.c
 # which is made available under the following licence
 
 ## Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. Permission