Merge pull request #122 from klinemichael/modpi_fix

Fixed mod2pi, mod1pi, modhalfpi, and modqrtrpi

Author: JeffreySarnoff

.DS_Store

@@ -1,11 +1,13 @@
+## Edited by Michael Kline, MIT License
+## This work was supported in part by the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists (WDTS) under the Science Undergraduate Laboratory Internships Program (SULI) program at Oak Ridge National Laboratory, administered by the Oak Ridge Institute for Science and Education.
+
 function mod2pi(x::DoubleFloat{T}) where {T<:IEEEFloat}
 	s = signbit(x)
-	if s
-		x = -x
-	end
-	x < pi2o1(DoubleFloat{T}) && return x
-	x < pi4o1(DoubleFloat{T}) && return mod2pi(x - pi2o1(DoubleFloat{T}))
-	himdlo = mul323(inv_pi_2o1_t64, HILO(x))
+	c = s ? -1 : 1
+    y = c * x
+	y < pi2o1(DoubleFloat{T}) && return s ? pi2o1(DoubleFloat{T}) - y : y
+	y < pi4o1(DoubleFloat{T}) && return mod2pi(x - c*pi2o1(DoubleFloat{T}))
+	himdlo = mul323(inv_pi_2o1_t64, HILO(y))
 	hi, md, lo = three_sum([modf(x)[1] for x in himdlo]...,)
 	if hi >= 1.0
 	    hi = hi - 1.0
@@ -23,12 +25,11 @@ end
 
 function mod1pi(x::DoubleFloat{T}) where {T<:IEEEFloat}
 	s = signbit(x)
-	if s
-		x = -x
-	end
-	x < pi1o1(DoubleFloat{T}) && return x
-	x < pi2o1(DoubleFloat{T}) && return mod1pi(x - pi1o1(DoubleFloat{T}))
-	himdlo = mul323(inv_pi_1o1_t64, HILO(x))
+	c = s ? -1 : 1
+    y = c * x
+	x < pi1o1(DoubleFloat{T}) && return s ? pi1o1(DoubleFloat{T}) - y : y
+	x < pi2o1(DoubleFloat{T}) && return mod1pi(x - c*pi1o1(DoubleFloat{T}))
+	himdlo = mul323(inv_pi_1o1_t64, HILO(y))
 	hi, md, lo = three_sum([modf(x)[1] for x in himdlo]...,)
 	if hi >= 1.0
 	    hi = hi - 1.0
@@ -46,12 +47,11 @@ end
 
 function modhalfpi(x::DoubleFloat{T}) where {T<:IEEEFloat}
 	s = signbit(x)
-	if s
-		x = -x
-	end
-	x < pi1o2(DoubleFloat{T}) && return x
-	x < pi1o1(DoubleFloat{T}) && return modhalfpi(x - pi1o2(DoubleFloat{T}))
-	himdlo = mul323(inv_pi_1o2_t64, HILO(x))
+	c = s ? -1 : 1
+    y = c * x
+	x < pi1o2(DoubleFloat{T}) && return s ? pi1o2(DoubleFloat{T}) - y : y
+	x < pi1o1(DoubleFloat{T}) && return modhalfpi(x - c*pi1o2(DoubleFloat{T}))
+	himdlo = mul323(inv_pi_1o2_t64, HILO(y))
 	hi, md, lo = three_sum([modf(x)[1] for x in himdlo]...,)
 	if hi >= 1.0
 	    hi = hi - 1.0
@@ -69,12 +69,11 @@ end
 
 function modqrtrpi(x::DoubleFloat{T}) where {T<:IEEEFloat}
 	s = signbit(x)
-	if s
-		x = -x
-	end
-	x < pi1o4(DoubleFloat{T}) && return x
-	x < pi1o2(DoubleFloat{T}) && return modqrtrpi(x - pi1o4(DoubleFloat{T}))
-	himdlo = mul323(inv_pi_1o4_t64, HILO(x))
+	c = s ? -1 : 1
+    y = c * x
+	x < pi1o4(DoubleFloat{T}) && return s ? pi1o4(DoubleFloat{T}) - y : y
+	x < pi1o2(DoubleFloat{T}) && return modqrtrpi(x - c*pi1o4(DoubleFloat{T}))
+	himdlo = mul323(inv_pi_1o4_t64, HILO(y))
 	hi, md, lo = three_sum([modf(x)[1] for x in himdlo]...,)
 	if hi >= 1.0
 	    hi = hi - 1.0
@@ -132,7 +131,7 @@ end
 # within (pi/2) quadrant determine nearest multiple of pi/16
 
 function whichquadrant(x::DoubleFloat{T}) where {T<:IEEEFloat}
-    x = negpi_pospi(x)	
+    x = negpi_pospi(x)
     if signbit(x) # quadrant -2 or -1
         quadrant = -x < DoubleFloat{T}(pi)/2 ? -2 : -1
     else          # quadrant  1 or  2

src/.DS_Store

@@ -3,7 +3,7 @@
 
     @test DoubleFloats.mod2pi(x) == Double64(6.111805926475162, -1.667563077572613e-16)
     @test DoubleFloats.mod1pi(x) == Double64(2.9702132728853683, 1.54868222178066e-16)
-    # these tests are unusable, sometimes they pass and othertimes they fail 
+    # these tests are unusable, sometimes they pass and othertimes they fail
     # @test DoubleFloats.modhalfpi(x) == Double64(1.3994169460904717, 9.362948122802731e-17)
     # @test DoubleFloats.modqrtrpi(x) == Double64(0.6140187826930236, -4.8008991213172176e-17)
 end
@@ -13,6 +13,6 @@ end
 
     @test DoubleFloats.rem2pi(x) == Double64(6.111805926475162, -1.667563077572613e-16)
     @test DoubleFloats.rem1pi(x) == Double64(2.9702132728853683, 1.54868222178066e-16)
-    @test DoubleFloats.rem2pi(-x) == Double64(-6.111805926475162, 1.667563077572613e-16)
-    @test DoubleFloats.rem1pi(-x) == Double64(-2.9702132728853683, -1.54868222178066e-16)
+    @test DoubleFloats.rem2pi(-x) == Double64(-0.1713793807044248, 4.647966647701768e-18)
+    @test DoubleFloats.rem1pi(-x) == Double64(2.9702132728853683, 1.54868222178066e-16)
 end