diff --git a/base/exports.jl b/base/exports.jl index e029518b47f44..d8a83cd7f4a1e 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -382,6 +382,7 @@ export mod, mod1, modf, + mod2pi, nan, nextfloat, nextpow, diff --git a/base/math.jl b/base/math.jl index a92550213d104..6a5476624f730 100644 --- a/base/math.jl +++ b/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 diff --git a/doc/manual/mathematical-operations.rst b/doc/manual/mathematical-operations.rst index b2d4c55fb11ae..7d6e9df38b851 100644 --- a/doc/manual/mathematical-operations.rst +++ b/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`` =============== ======================================================================= diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst index 38ffce7d0e6fa..c83a326362d82 100644 --- a/doc/stdlib/base.rst +++ b/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]]) diff --git a/test/Makefile b/test/Makefile index bad86d6c6d872..2c503fce875b5 100644 --- a/test/Makefile +++ b/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 diff --git a/test/mod2pi.jl b/test/mod2pi.jl new file mode 100644 index 0000000000000..79dfa5c1c4302 --- /dev/null +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index 297762c186572..68397af2f1790 100644 --- a/test/runtests.jl +++ b/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