diff --git a/base/exports.jl b/base/exports.jl index e029518b47f44..bb3d933994c83 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -381,7 +381,10 @@ export maxintfloat, mod, mod1, + mod2pi, modf, + modpi, + modpio2, nan, nextfloat, nextpow, diff --git a/base/math.jl b/base/math.jl index a92550213d104..a3b84fff58923 100644 --- a/base/math.jl +++ b/base/math.jl @@ -12,6 +12,7 @@ export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, ceil, floor, trunc, round, significand, lgamma, hypot, gamma, lfact, max, min, ldexp, frexp, clamp, modf, ^, + mod2pi, modpi, modpio2, airy, airyai, airyprime, airyaiprime, airybi, airybiprime, besselj0, besselj1, besselj, bessely0, bessely1, bessely, hankelh1, hankelh2, besseli, besselk, besselh, @@ -644,9 +645,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 +665,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 +686,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 +1338,169 @@ end erfcinv(x::Integer) = erfcinv(float(x)) @vectorize_1arg Real erfcinv + +function add22Cond(xh::Float64, xl::Float64, yh::Float64, yl::Float64) + # This algorithm, due to Dekker [1], computes the sum of + # two double-double numbers as a double-double, with a relative error smaller than 2^−103 + # [1] http://gdz.sub.uni-goettingen.de/dms/load/img/?PPN=PPN362160546_0018&DMDID=DMDLOG_0023&LOGID=LOG_0023&PHYSID=PHYS_0232 + # [2] http://ftp.nluug.nl/pub/os/BSD/FreeBSD/distfiles/crlibm/crlibm-1.0beta3.pdf + r = xh+yh + s = (abs(xh) > abs(yh)) ? (xh-r+yh+yl+xl) : (yh-r+xh+xl+yl) + zh = r+s + zl = r-zh+s + return (zh,zl) +end + +function add22Cond_h(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?source=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 add22Cond_h(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 add22Cond_h(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 add22Cond_h(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 add22Cond_h(y[1],y[2],pi3o2_h,pi3o2_l) + else # add pi/2 + return add22Cond_h(y[1],y[2],pi1o2_h,pi1o2_l) + end + end +end + +function modpi(x::Float64) +# with r = modpi(x) +# a) 0 <= r < π (note: boundary open or closed - a bit fuzzy, due to rem_pio2 implementation) +# b) r-x = k*π with k integer + if x < pi2o2_h + if 0.0 <= x return x end + if x > -pi2o2_h + return add22Cond_h(x,0.0,pi2o2_h,pi2o2_l) + end + end + (n,y) = ieee754_rem_pio2(x) + if iseven(n) + if y[1] > 0.0 + return y[1] + else # else add pi + return add22Cond_h(y[1],y[2],pi2o2_h,pi2o2_l) + end + else # add pi/2 + return add22Cond_h(y[1],y[2],pi1o2_h,pi1o2_l) + end +end + + +function modpio2(x::Float64) +# with r = modpio2(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: we explicitly test for 0 <= values < pi/2, because +# ieee754_rem_pio2 behaves weirdly for arguments that are already +# within -pi/4, pi/4 e.g. +# ieee754_rem_pio2(0.19633954084936206) returns +# (1,[-1.3744567859455346,-6.12323399538461e-17]) +# which does not add up to 0.19633954084936206 + if x < pi1o2_h + if x >= 0.0 return x end + if x > -pi1o2_h + zh = add22Cond_h(x,0.0,pi1o2_h,pi1o2_l) + return zh + end + end + (n,y) = ieee754_rem_pio2(x) + if y[1] > 0.0 + return y[1] + else + zh = add22Cond_h(y[1],y[2],pi1o2_h,pi1o2_l) + return zh + 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 + +modpi(x::Float32) = float32(modpi(float64(x))) +modpi(x::Int32) = modpi(float64(x)) +function modpi(x::Int64) + fx = float64(x) + fx == x || error("Integer argument to modpi is too large: $x") + modpi(fx) +end + +modpio2(x::Float32)= float32(modpio2(float64(x))) +modpio2(x::Int32) = modpio2(float64(x)) +function modpio2(x::Int64) + fx = float64(x) + fx == x || error("Integer argument to modpio2 is too large: $x") + modpio2(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/math-modpi.jl b/test/math-modpi.jl new file mode 100644 index 0000000000000..51cfefc690592 --- /dev/null +++ b/test/math-modpi.jl @@ -0,0 +1,233 @@ + + +# 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, x mod 2pi, x mod pi, x mod pi/2 (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.5707963267948968 6.123233995736766e-17; + 1.5707963267948966 1.5707963267948966 1.5707963267948966 1.5707963267948966; + -3.141592653589793 3.1415926535897936 1.2246467991473532e-16 1.2246467991473532e-16; + 3.141592653589793 3.141592653589793 3.141592653589793 1.5707963267948966; + 6.283185307179586 6.283185307179586 3.141592653589793 1.5707963267948963; + -6.283185307179586 2.4492935982947064e-16 2.4492935982947064e-16 2.4492935982947064e-16; + 45.553093477052 1.5707963267948966 1.5707963267948966 6.189806365883577e-19; + 3.14159265359 3.14159265359 2.0682310711021444e-13 2.0682310711021444e-13; + -3.14159265359 3.1415926535895866 3.1415926535895866 1.5707963267946898; + -3.9269808169872418 2.356204490192345 2.356204490192345 0.7854081633974481; + -3.73063127613788 2.5525540310417068 2.5525540310417068 0.98175770424681; + -3.5342817352885176 2.748903571891069 2.748903571891069 1.1781072450961723; + -3.337932194439156 2.945253112740431 2.945253112740431 1.374456785945534; + -3.1415826535897935 3.141602653589793 9.999999999743887e-6 9.999999999743887e-6; + -2.9452331127404316 3.337952194439155 0.19635954084936158 0.19635954084936158; + -2.7488835718910694 3.5343017352885173 0.39270908169872387 0.39270908169872387; + -2.5525340310417075 3.730651276137879 0.5890586225480857 0.5890586225480857; + -2.356184490192345 3.9270008169872415 0.785408163397448 0.785408163397448; + -2.1598349493429834 4.123350357836603 0.9817577042468099 0.9817577042468099; + -1.9634854084936209 4.319699898685966 1.1781072450961725 1.1781072450961725; + -1.7671358676442588 4.516049439535328 1.3744567859455346 1.3744567859455346; + -1.5707863267948967 4.71239898038469 1.5708063267948966 9.9999999999047e-6; + -1.3744367859455346 4.908748521234052 1.7671558676442587 0.19635954084936197; + -1.1780872450961726 5.105098062083414 1.9635054084936208 0.39270908169872404; + -0.9817377042468104 5.301447602932776 2.159854949342983 0.5890586225480863; + -0.7853881633974483 5.4977971437821385 2.356204490192345 0.7854081633974483; + -0.5890386225480863 5.6941466846315 2.552554031041707 0.9817577042468104; + -0.3926890816987242 5.890496225480862 2.748903571891069 1.1781072450961725; + -0.1963395408493621 6.0868457663302244 2.9452531127404313 1.3744567859455346; + 1.0e-5 1.0e-5 1.0e-5 1.0e-5; + 0.19635954084936205 0.19635954084936205 0.19635954084936205 0.19635954084936205; + 0.3927090816987241 0.3927090816987241 0.3927090816987241 0.3927090816987241; + 0.5890586225480862 0.5890586225480862 0.5890586225480862 0.5890586225480862; + 0.7854081633974482 0.7854081633974482 0.7854081633974482 0.7854081633974482; + 0.9817577042468103 0.9817577042468103 0.9817577042468103 0.9817577042468103; + 1.1781072450961723 1.1781072450961723 1.1781072450961723 1.1781072450961723; + 1.3744567859455343 1.3744567859455343 1.3744567859455343 1.3744567859455343; + 1.5708063267948964 1.5708063267948964 1.5708063267948964 9.999999999782235e-6; + 1.7671558676442585 1.7671558676442585 1.7671558676442585 0.19635954084936186; + 1.9635054084936205 1.9635054084936205 1.9635054084936205 0.3927090816987239; + 2.159854949342982 2.159854949342982 2.159854949342982 0.5890586225480855; + 2.3562044901923445 2.3562044901923445 2.3562044901923445 0.7854081633974478; + 2.5525540310417063 2.5525540310417063 2.5525540310417063 0.9817577042468096; + 2.7489035718910686 2.7489035718910686 2.7489035718910686 1.178107245096172; + 2.9452531127404304 2.9452531127404304 2.9452531127404304 1.3744567859455339; + 3.1416026535897927 3.1416026535897927 9.999999999498959e-6 9.999999999498959e-6; + 3.3379521944391546 3.3379521944391546 0.19635954084936136 0.19635954084936136; + 3.534301735288517 3.534301735288517 0.39270908169872365 0.39270908169872365; + 3.7306512761378787 3.7306512761378787 0.5890586225480855 0.5890586225480855; + 3.927000816987241 3.927000816987241 0.7854081633974478 0.7854081633974478; + -3.9270008169872415 2.356184490192345 2.356184490192345 0.7853881633974484; + -3.7306512761378796 2.552534031041707 2.552534031041707 0.9817377042468103; + -3.5343017352885173 2.7488835718910694 2.7488835718910694 1.1780872450961726; + -3.3379521944391555 2.945233112740431 2.945233112740431 1.3744367859455344; + -3.141602653589793 3.1415826535897935 3.1415826535897935 1.5707863267948967; + -2.9452531127404313 3.3379321944391553 0.1963395408493619 0.1963395408493619; + -2.748903571891069 3.5342817352885176 0.3926890816987242 0.3926890816987242; + -2.552554031041707 3.7306312761378795 0.589038622548086 0.589038622548086; + -2.356204490192345 3.9269808169872418 0.7853881633974483 0.7853881633974483; + -2.159854949342983 4.123330357836603 0.9817377042468102 0.9817377042468102; + -1.9635054084936208 4.3196798986859655 1.1780872450961726 1.1780872450961726; + -1.7671558676442587 4.516029439535328 1.3744367859455346 1.3744367859455346; + -1.5708063267948966 4.71237898038469 1.5707863267948967 1.5707863267948967; + -1.3744567859455346 4.908728521234052 1.7671358676442588 0.19633954084936206; + -1.1781072450961725 5.105078062083414 1.9634854084936209 0.39268908169872413; + -0.9817577042468104 5.301427602932776 2.159834949342983 0.5890386225480863; + -0.7854081633974483 5.497777143782138 2.3561844901923448 0.7853881633974483; + -0.5890586225480863 5.694126684631501 2.552534031041707 0.9817377042468104; + -0.39270908169872415 5.890476225480862 2.748883571891069 1.1780872450961726; + -0.19635954084936208 6.086825766330224 2.945233112740431 1.3744367859455346; + -1.0e-5 6.283175307179587 3.141582653589793 1.5707863267948967; + 0.19633954084936206 0.19633954084936206 0.19633954084936206 0.19633954084936206; + 0.39268908169872413 0.39268908169872413 0.39268908169872413 0.39268908169872413; + 0.5890386225480861 0.5890386225480861 0.5890386225480861 0.5890386225480861; + 0.7853881633974482 0.7853881633974482 0.7853881633974482 0.7853881633974482; + 0.9817377042468103 0.9817377042468103 0.9817377042468103 0.9817377042468103; + 1.1780872450961724 1.1780872450961724 1.1780872450961724 1.1780872450961724; + 1.3744367859455344 1.3744367859455344 1.3744367859455344 1.3744367859455344; + 1.5707863267948965 1.5707863267948965 1.5707863267948965 1.5707863267948965; + 1.7671358676442586 1.7671358676442586 1.7671358676442586 0.19633954084936195; + 1.9634854084936206 1.9634854084936206 1.9634854084936206 0.392689081698724; + 2.1598349493429825 2.1598349493429825 2.1598349493429825 0.5890386225480858; + 2.3561844901923448 2.3561844901923448 2.3561844901923448 0.7853881633974481; + 2.5525340310417066 2.5525340310417066 2.5525340310417066 0.98173770424681; + 2.748883571891069 2.748883571891069 2.748883571891069 1.1780872450961724; + 2.9452331127404308 2.9452331127404308 2.9452331127404308 1.3744367859455342; + 3.141582653589793 3.141582653589793 3.141582653589793 1.5707863267948965; + 3.337932194439155 3.337932194439155 0.19633954084936167 0.19633954084936167; + 3.534281735288517 3.534281735288517 0.39268908169872396 0.39268908169872396; + 3.730631276137879 3.730631276137879 0.5890386225480858 0.5890386225480858; + 3.9269808169872413 3.9269808169872413 0.7853881633974481 0.7853881633974481; + 22.0 3.1504440784612404 0.008851424871447331 0.008851424871447331; + 333.0 6.2743640266615035 3.13277137307171 1.5619750462768134; + 355.0 3.1416227979431572 3.014435336405372e-5 3.014435336405372e-5; + 103993.0 6.283166177843807 3.1415735242540137 1.570777197459117; + 104348.0 3.141603668607378 1.10150175844633e-5 1.10150175844633e-5; + 208341.0 3.141584539271598 3.141584539271598 1.5707882124767014; + 312689.0 2.9006993893361787e-6 2.9006993893361787e-6 2.9006993893361787e-6; + 833719.0 3.1415903406703767 3.1415903406703767 1.5707940138754801; + 1.146408e6 3.1415932413697663 5.877799728814151e-7 5.877799728814151e-7; + 4.272943e6 6.283184757600089 3.1415921040102956 1.5707957772153989; + 5.419351e6 3.1415926917902683 3.820047507089661e-8 3.820047507089661e-8; + 8.0143857e7 6.283185292406739 3.1415926388169466 1.5707963120220498; + 1.65707065e8 3.1415926622445745 8.654781434964792e-9 8.654781434964792e-9; + 2.45850922e8 3.141592647471728 3.141592647471728 1.5707963206768312; + 4.11557987e8 2.5367160519636766e-9 2.5367160519636766e-9 2.5367160519636766e-9; + 1.068966896e9 3.14159265254516 3.14159265254516 1.5707963257502633; + 2.549491779e9 4.474494938161497e-10 4.474494938161497e-10 4.474494938161497e-10; + 6.167950454e9 3.141592653440059 3.141592653440059 1.5707963266451623; + 1.4885392687e10 1.4798091093322177e-10 1.4798091093322177e-10 1.4798091093322177e-10; + 2.1053343141e10 3.14159265358804 3.14159265358804 1.5707963267931433; + 1.783366216531e12 6.969482408757582e-13 6.969482408757582e-13 6.969482408757582e-13; + 3.587785776203e12 3.141592653589434 3.141592653589434 1.570796326794537; + 5.371151992734e12 3.1415926535901306 3.374642143850602e-13 3.374642143850602e-13; + 8.958937768937e12 6.283185307179564 3.1415926535897714 1.5707963267948746; + 1.39755218526789e14 3.1415926535898 7.167032800493559e-15 7.167032800493559e-15; + 4.28224593349304e14 3.1415926535897927 3.1415926535897927 1.5707963267948961; + 5.706674932067741e15 4.237546464512562e-16 4.237546464512562e-16 4.237546464512562e-16; + 6.134899525417045e15 3.141592653589793 3.141592653589793 1.5707963267948966; +] + +function testModPi() + verbose = false + numTestCases = size(testCases,1) + println("Testing mod2pi, modpi, modpio2 with $numTestCases instances... ") + + modFns = [mod2pi,modpi,modpio2] + xDivisors = [2pi,pi,pi/2] + errsNew, errsOld = Array(Float64,0), Array(Float64,0) + for rowIdx in 1:numTestCases + xExact = testCases[rowIdx,1] + verbose && print(lpad(string(xExact),25," ")) + for colIdx in 1:3 + 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 + if verbose + if sum(errsNew[end-2:]) == 0 + print(" ok ") + else print(" ERR ") end + print("diffs new $(errsNew[end-2]) $(errsNew[end-1]) $(errsNew[end]) ") + print("diffs old $(errsOld[end-2]) $(errsOld[end-1]) $(errsOld[end])\n") + end + end + sort!(errsNew) + sort!(errsOld) + totalErrNew = sum(errsNew) + totalErrOld = sum(errsOld) + @test_approx_eq totalErrNew 0.0 + println("Total err = $totalErrNew (new), $totalErrOld (old).") +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) + +# pi +@test_approx_eq modpi(10) mod(10,pi) +@test_approx_eq modpi(-10) mod(-10,pi) +@test_approx_eq modpi(355) 3.014435336405372e-5 +@test_approx_eq modpi(int32(355)) 3.014435336405372e-5 +@test_approx_eq modpi(355.0) 3.014435336405372e-5 +@test_approx_eq modpi(355.0f0) 3.0144354f-5 +@test modpi(2^60) == modpi(2.0^60) +@test_throws modpi(2^60-1) + + +# pi/2 +@test_approx_eq modpio2(10) mod(10,pi/2) +@test_approx_eq modpio2(-10) mod(-10,pi/2) +@test_approx_eq modpio2(355) 3.014435336405372e-5 +@test_approx_eq modpio2(int32(355)) 3.014435336405372e-5 +@test_approx_eq modpio2(355.0) 3.014435336405372e-5 +@test_approx_eq modpio2(355.0f0) 3.0144354f-5 +@test modpio2(2^60) == modpio2(2.0^60) +@test_throws modpio2(2^60-1) + diff --git a/test/runtests.jl b/test/runtests.jl index 297762c186572..8ea8272cbe652 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,8 @@ 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", + "math-modpi"] tests = ARGS==["all"] ? testnames : ARGS