pared down mod2pi pull request

base/exports.jl

@@ -382,6 +382,7 @@ export
     mod,
     mod1,
     modf,
+    mod2pi,
     nan,
     nextfloat,
     nextpow,

base/math.jl

@@ -11,7 +11,7 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
        cbrt, sqrt, erf, erfc, erfcx, erfi, dawson,
        ceil, floor, trunc, round, significand, 
        lgamma, hypot, gamma, lfact, max, min, ldexp, frexp,
-       clamp, modf, ^, 
+       clamp, modf, ^, mod2pi,
        airy, airyai, airyprime, airyaiprime, airybi, airybiprime,
        besselj0, besselj1, besselj, bessely0, bessely1, bessely,
        hankelh1, hankelh2, besseli, besselk, besselh,
@@ -644,9 +644,9 @@ hankelh2(nu, z) = besselh(nu, 2, z)
 
 
 function angle_restrict_symm(theta)
-    P1 = 4 * 7.8539812564849853515625e-01
-    P2 = 4 * 3.7748947079307981766760e-08
-    P3 = 4 * 2.6951514290790594840552e-15
+    const P1 = 4 * 7.8539812564849853515625e-01
+    const P2 = 4 * 3.7748947079307981766760e-08
+    const P3 = 4 * 2.6951514290790594840552e-15
 
     y = 2*floor(theta/(2*pi))
     r = ((theta - y*P1) - y*P2) - y*P3
@@ -664,7 +664,7 @@ const clg_coeff = [76.18009172947146,
                    -0.5395239384953e-5]
 
 function clgamma_lanczos(z)
-    sqrt2pi = 2.5066282746310005
+    const sqrt2pi = 2.5066282746310005
 
     y = x = z
     temp = x + 5.5
@@ -685,7 +685,7 @@ function lgamma(z::Complex)
     if real(z) <= 0.5
         a = clgamma_lanczos(1-z)
         b = log(sinpi(z))
-        logpi = 1.14472988584940017
+        const logpi = 1.14472988584940017
         z = logpi - b - a
     else
         z = clgamma_lanczos(z)
@@ -1337,4 +1337,90 @@ end
 erfcinv(x::Integer) = erfcinv(float(x))
 @vectorize_1arg Real erfcinv
 
+## mod2pi-related calculations ##
+
+function add22condh(xh::Float64, xl::Float64, yh::Float64, yl::Float64)
+    # as above, but only compute and return high double
+    r = xh+yh
+    s = (abs(xh) > abs(yh)) ? (xh-r+yh+yl+xl) : (yh-r+xh+xl+yl)
+    zh = r+s
+    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/openlibm/blob/master/src/e_rem_pio2.c
+
+    y = [0.0,0.0]
+    n = ccall((:__ieee754_rem_pio2, libm), 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)
+
+const pi2o2_h  = 3.141592653589793      # convert(Float64, pi * BigFloat(1))
+const pi2o2_l  = 1.2246467991473532e-16 # convert(Float64, pi * BigFloat(1) - pi2o2_h)
+
+const pi3o2_h  = 4.71238898038469       # convert(Float64, pi * BigFloat(3/2))
+const pi3o2_l  = 1.8369701987210297e-16 # convert(Float64, pi * BigFloat(3/2) - pi3o2_h)
+
+const pi4o2_h  = 6.283185307179586      # convert(Float64, pi * BigFloat(2))
+const pi4o2_l  = 2.4492935982947064e-16 # convert(Float64, pi * BigFloat(2) - pi4o2_h)
+
+function mod2pi(x::Float64) # or modtau(x)
+# with r = mod2pi(x)
+# a) 0 <= r < 2π  (note: boundary open or closed - a bit fuzzy, due to rem_pio2 implementation)
+# b) r-x = k*2π with k integer
+
+# note: mod(n,4) is 0,1,2,3; while mod(n-1,4)+1 is 1,2,3,4.
+# We use the latter to push negative y in quadrant 0 into the positive (one revolution, + 4*pi/2)
+
+    if x < pi4o2_h
+        if 0.0 <= x return x end
+        if x > -pi4o2_h 
+            return add22condh(x,0.0,pi4o2_h,pi4o2_l)
+        end
+    end
+
+    (n,y) = ieee754_rem_pio2(x)
+
+    if iseven(n)
+        if n & 2 == 2 # add pi
+            return add22condh(y[1],y[2],pi2o2_h,pi2o2_l)
+        else # add 0 or 2pi
+            if y[1] > 0.0
+                return y[1]
+            else # else add 2pi
+                return add22condh(y[1],y[2],pi4o2_h,pi4o2_l)
+            end
+        end
+    else # add pi/2 or 3pi/2
+        if n & 2 == 2 # add 3pi/2
+            return add22condh(y[1],y[2],pi3o2_h,pi3o2_l) 
+        else # add pi/2
+            return add22condh(y[1],y[2],pi1o2_h,pi1o2_l) 
+        end
+    end
+end
+
+mod2pi(x::Float32) = float32(mod2pi(float64(x)))
+mod2pi(x::Int32) = mod2pi(float64(x))
+function mod2pi(x::Int64)
+  fx = float64(x)
+  fx == x || error("Integer argument to mod2pi is too large: $x")
+  mod2pi(fx)
+end
+
 end # module

doc/manual/mathematical-operations.rst

@@ -314,6 +314,9 @@ Function        Description
 ``fld(x,y)``    floored division; quotient rounded towards ``-Inf``
 ``rem(x,y)``    remainder; satisfies ``x == div(x,y)*y + rem(x,y)``; sign matches ``x``
 ``mod(x,y)``    modulus; satisfies ``x == fld(x,y)*y + mod(x,y)``; sign matches ``y``
+``mod2pi(x)``   modulus with respect to 2pi;  ``0 <= mod2pi(x)  < 2pi``
+``modpi(x)``    modulus with respect to pi;   ``0 <= modpi(x)   < pi``
+``modpio2(x)``  modulus with respect to pi/2; ``0 <= modpio2(x) < pi/2``
 ``gcd(x,y...)`` greatest common divisor of ``x``, ``y``,...; sign matches ``x``
 ``lcm(x,y...)`` least common multiple of ``x``, ``y``,...; sign matches ``x``
 =============== =======================================================================

doc/stdlib/base.rst

@@ -1972,6 +1972,18 @@ Mathematical Operators
 
    Modulus after division, returning in the range [0,m)
 
+.. function:: modpi(x)
+
+   Modulus after division by pi, returning in the range [0,pi). More accurate than mod(x,pi).
+
+.. function:: mod2pi(x)
+
+   Modulus after division by 2pi, returning in the range [0,2pi). More accurate than mod(x,2pi).
+
+.. function:: modpio2(x)
+
+   Modulus after division by pi/2, returning in the range [0,pi/2). More accurate than mod(x,pi/2).
+
 .. function:: rem(x, m)
 
    Remainder after division
@@ -2468,7 +2480,7 @@ Mathematical Functions
 
 .. function:: round(x, [digits, [base]])
 
-   ``round(x)`` returns the nearest integral value of the same type as ``x`` to ``x``. ``round(x, digits)`` rounds to the specified number of digits after the decimal place, or before if negative, e.g., ``round(pi,2)`` is ``3.14``. ``round(x, digits, base)`` rounds using a different base, defaulting to 10, e.g., ``round(pi, 3, 2)`` is ``3.125``.
+   ``round(x)`` returns the nearest integral value of the same type as ``x`` to ``x``. ``round(x, digits)`` rounds to the specified number of digits after the decimal place, or before if negative, e.g., ``round(pi,2)`` is ``3.14``. ``round(x, digits, base)`` rounds using a different base, defaulting to 10, e.g., ``round(pi, 1, 8)`` is ``3.125``.
 
 .. function:: ceil(x, [digits, [base]])
 

test/Makefile

@@ -6,7 +6,7 @@ TESTS = all core keywordargs numbers strings unicode collections hashing	\
 	math functional bigint sorting statistics spawn parallel arpack file	\
 	git pkg resolve suitesparse complex version pollfd mpfr	broadcast       \
 	socket floatapprox priorityqueue readdlm regex float16 combinatorics    \
-	sysinfo rounding ranges
+	sysinfo rounding ranges mod2pi
 
 default: all
 

test/mod2pi.jl

@@ -0,0 +1,195 @@
+# NOTES on range reduction
+# [1] compute numbers near pi: http://www.cs.berkeley.edu/~wkahan/testpi/nearpi.c
+# [2] range reduction: http://hal-ujm.ccsd.cnrs.fr/docs/00/08/69/04/PDF/RangeReductionIEEETC0305.pdf
+# [3] precise addition, see Add22: http://ftp.nluug.nl/pub/os/BSD/FreeBSD/distfiles/crlibm/crlibm-1.0beta3.pdf
+
+# Examples:
+# ΓΓ = 6411027962775774 / 2^47  # see [2] above, section 1.2
+# julia> mod(ΓΓ, pi/2)   # "naive" way - easily wrong
+# 1.7763568394002505e-15
+# julia> modpio2( ΓΓ )   # using function provided here
+# 6.189806365883577e-19
+# Wolfram Alpha: mod(6411027962775774 / 2^47, pi/2)
+# 6.189806365883577000150671465609655958633034115366621088... × 10^-19
+
+
+# Test Cases. Each row contains: x and x mod 2pi (as from Wolfram Alpha)
+# The values x are:
+# -pi/2, pi/2, -pi, pi, 2pi, -2pi 
+#   (or rather, the Float64 approx to those numbers. 
+#   Thus, x mod pi will result in a small, but positive number)
+# ΓΓ = 6411027962775774 / 2^47  
+#   from [2], section 1.2: 
+#   the Float64 greater than 8, and less than 2**63 − 1 closest to a multiple of π/4 is
+#   Γ = 6411027962775774 / 2^48. We take ΓΓ = 2*Γ to get cancellation with pi/2 already
+# 3.14159265359, -3.14159265359
+# pi/16*k +/- 0.00001 for k in [-20:20] # to cover all quadrants
+# numerators of continuous fraction approximations to pi
+#   see http://oeis.org/A002485
+#   (reason: for max cancellation, we want x = k*pi + eps for small eps, so x/k ≈ pi)
+
+testCases = [
+      -1.5707963267948966         4.71238898038469
+       1.5707963267948966       1.5707963267948966
+       -3.141592653589793       3.1415926535897936
+        3.141592653589793        3.141592653589793
+        6.283185307179586        6.283185307179586
+       -6.283185307179586   2.4492935982947064e-16
+          45.553093477052       1.5707963267948966
+            3.14159265359            3.14159265359
+           -3.14159265359       3.1415926535895866
+      -3.9269808169872418        2.356204490192345
+        -3.73063127613788       2.5525540310417068
+      -3.5342817352885176        2.748903571891069
+       -3.337932194439156        2.945253112740431
+      -3.1415826535897935        3.141602653589793
+      -2.9452331127404316        3.337952194439155
+      -2.7488835718910694       3.5343017352885173
+      -2.5525340310417075        3.730651276137879
+       -2.356184490192345       3.9270008169872415
+      -2.1598349493429834        4.123350357836603
+      -1.9634854084936209        4.319699898685966
+      -1.7671358676442588        4.516049439535328
+      -1.5707863267948967         4.71239898038469
+      -1.3744367859455346        4.908748521234052
+      -1.1780872450961726        5.105098062083414
+      -0.9817377042468104        5.301447602932776
+      -0.7853881633974483       5.4977971437821385
+      -0.5890386225480863          5.6941466846315
+      -0.3926890816987242        5.890496225480862
+      -0.1963395408493621       6.0868457663302244
+                   1.0e-5                   1.0e-5
+      0.19635954084936205      0.19635954084936205
+       0.3927090816987241       0.3927090816987241
+       0.5890586225480862       0.5890586225480862
+       0.7854081633974482       0.7854081633974482
+       0.9817577042468103       0.9817577042468103
+       1.1781072450961723       1.1781072450961723
+       1.3744567859455343       1.3744567859455343
+       1.5708063267948964       1.5708063267948964
+       1.7671558676442585       1.7671558676442585
+       1.9635054084936205       1.9635054084936205
+        2.159854949342982        2.159854949342982
+       2.3562044901923445       2.3562044901923445
+       2.5525540310417063       2.5525540310417063
+       2.7489035718910686       2.7489035718910686
+       2.9452531127404304       2.9452531127404304
+       3.1416026535897927       3.1416026535897927
+       3.3379521944391546       3.3379521944391546
+        3.534301735288517        3.534301735288517
+       3.7306512761378787       3.7306512761378787
+        3.927000816987241        3.927000816987241
+      -3.9270008169872415        2.356184490192345
+      -3.7306512761378796        2.552534031041707
+      -3.5343017352885173       2.7488835718910694
+      -3.3379521944391555        2.945233112740431
+       -3.141602653589793       3.1415826535897935
+      -2.9452531127404313       3.3379321944391553
+       -2.748903571891069       3.5342817352885176
+       -2.552554031041707       3.7306312761378795
+       -2.356204490192345       3.9269808169872418
+       -2.159854949342983        4.123330357836603
+      -1.9635054084936208       4.3196798986859655
+      -1.7671558676442587        4.516029439535328
+      -1.5708063267948966         4.71237898038469
+      -1.3744567859455346        4.908728521234052
+      -1.1781072450961725        5.105078062083414
+      -0.9817577042468104        5.301427602932776
+      -0.7854081633974483        5.497777143782138
+      -0.5890586225480863        5.694126684631501
+     -0.39270908169872415        5.890476225480862
+     -0.19635954084936208        6.086825766330224
+                  -1.0e-5        6.283175307179587
+      0.19633954084936206      0.19633954084936206
+      0.39268908169872413      0.39268908169872413
+       0.5890386225480861       0.5890386225480861
+       0.7853881633974482       0.7853881633974482
+       0.9817377042468103       0.9817377042468103
+       1.1780872450961724       1.1780872450961724
+       1.3744367859455344       1.3744367859455344
+       1.5707863267948965       1.5707863267948965
+       1.7671358676442586       1.7671358676442586
+       1.9634854084936206       1.9634854084936206
+       2.1598349493429825       2.1598349493429825
+       2.3561844901923448       2.3561844901923448
+       2.5525340310417066       2.5525340310417066
+        2.748883571891069        2.748883571891069
+       2.9452331127404308       2.9452331127404308
+        3.141582653589793        3.141582653589793
+        3.337932194439155        3.337932194439155
+        3.534281735288517        3.534281735288517
+        3.730631276137879        3.730631276137879
+       3.9269808169872413       3.9269808169872413
+                     22.0       3.1504440784612404
+                    333.0       6.2743640266615035
+                    355.0       3.1416227979431572
+                 103993.0        6.283166177843807
+                 104348.0        3.141603668607378
+                 208341.0        3.141584539271598
+                 312689.0    2.9006993893361787e-6
+                 833719.0       3.1415903406703767
+               1.146408e6       3.1415932413697663
+               4.272943e6        6.283184757600089
+               5.419351e6       3.1415926917902683
+              8.0143857e7        6.283185292406739
+             1.65707065e8       3.1415926622445745
+             2.45850922e8        3.141592647471728
+             4.11557987e8    2.5367160519636766e-9
+            1.068966896e9         3.14159265254516
+            2.549491779e9    4.474494938161497e-10
+            6.167950454e9        3.141592653440059
+          1.4885392687e10   1.4798091093322177e-10
+          2.1053343141e10         3.14159265358804
+        1.783366216531e12    6.969482408757582e-13
+        3.587785776203e12        3.141592653589434
+        5.371151992734e12       3.1415926535901306
+        8.958937768937e12        6.283185307179564
+      1.39755218526789e14          3.1415926535898
+      4.28224593349304e14       3.1415926535897927
+     5.706674932067741e15    4.237546464512562e-16
+     6.134899525417045e15        3.141592653589793
+]
+
+function testModPi()
+    numTestCases = size(testCases,1)
+    modFns = [mod2pi]
+    xDivisors = [2pi]
+    errsNew, errsOld = Array(Float64,0), Array(Float64,0)
+    for rowIdx in 1:numTestCases
+        xExact = testCases[rowIdx,1]
+        for colIdx in 1:1
+            xSoln = testCases[rowIdx,colIdx+1]
+            xDivisor = xDivisors[colIdx]
+            modFn = modFns[colIdx]
+            # 2. want: xNew := modFn(xExact)  ≈  xSoln  <--- this is the crucial bit, xNew close to xSoln
+            # 3. know: xOld := mod(xExact,xDivisor) might be quite a bit off from xSoln - that's expected
+            xNew = modFn(xExact)
+            xOld = mod(xExact,xDivisor)
+
+            newDiff  = abs(xNew - xSoln)  # should be zero, ideally (our new function)
+            oldDiff  = abs(xOld - xSoln)  # should be zero in a perfect world, but often bigger due to cancellation
+            oldDiff  = min(oldDiff, abs(xDivisor - oldDiff)) # we are being generous here:
+            # if xOld happens to end up "on the wrong side of 0", eg 
+            # if xSoln = 3.14 (correct), but xOld reports 0.01, 
+            # we don't take the long way around the circle of 3.14 - 0.01, but the short way of 3.1415.. - (3.14 - 0.1) 
+            push!(errsNew,abs(newDiff))
+            push!(errsOld,abs(oldDiff))
+        end
+    end
+    sort!(errsNew)
+    sort!(errsOld)
+    totalErrNew = sum(errsNew)
+    totalErrOld = sum(errsOld)
+    @test_approx_eq totalErrNew 0.0
+end
+testModPi()
+
+# 2pi
+@test_approx_eq mod2pi(10)          mod(10,2pi)
+@test_approx_eq mod2pi(-10)         mod(-10,2pi)
+@test_approx_eq mod2pi(355)         3.1416227979431572
+@test_approx_eq mod2pi(int32(355))  3.1416227979431572
+@test_approx_eq mod2pi(355.0)       3.1416227979431572
+@test_approx_eq mod2pi(355.0f0)     3.1416228f0
+@test mod2pi(2^60) == mod2pi(2.0^60)
+@test_throws mod2pi(2^60-1)

test/runtests.jl

@@ -6,7 +6,7 @@ testnames = ["core", "keywordargs", "numbers", "strings", "unicode",
              "arpack", "file", "suitesparse", "version",
              "resolve", "pollfd", "mpfr", "broadcast", "complex",
              "socket", "floatapprox", "readdlm", "regex", "float16",
-             "combinatorics", "sysinfo", "rounding", "ranges"]
+             "combinatorics", "sysinfo", "rounding", "ranges", "mod2pi"]
 
 tests = ARGS==["all"] ? testnames : ARGS