diff --git a/src/unitful.jl b/src/unitful.jl
index 338250cb..bba23f68 100644
--- a/src/unitful.jl
+++ b/src/unitful.jl
@@ -22,14 +22,35 @@ for PT in ParticleSymbols
             $PT{Quantity{S,D,U},N}(fill(y, N))
 
         end
+    end
 
-        function Base.:*(p::$PT{T,N}, y::Quantity{S,D,U}) where {S, D, U, T, N}
-            NT = promote_type(S,T)
-            $PT{Quantity{NT,D,U},N}(p.particles .* y)
-        end
-
-        function Base.:*(p::$PT, y::FreeUnits)
-            $PT(p.particles .* y)
+    for op in (*, /)
+        f = nameof(op)
+        @eval begin
+            function Base.$f(p::$PT{T,N}, y::Quantity{S,D,U}) where {S, D, U, T, N}
+                NT = promote_type(T, S)
+                $PT{Quantity{NT,D,U},N}($(op).(p.particles , y))
+            end
+
+            function Base.$f(p::$PT{T,N}, y::Quantity{S,D,U}) where {S, D, U, T <: Quantity, N}
+                QT = Base.promote_op($op, T, typeof(y))
+                $PT{QT,N}($(op).(p.particles, y))
+            end
+
+            # Below is just the reverse signature of above
+            function Base.$f(y::Quantity{S,D,U}, p::$PT{T,N}) where {S, D, U, T, N}
+                NT = promote_type(T, S)
+                $PT{Quantity{NT,D,U},N}($(op).(y, p.particles))
+            end
+
+            function Base.$f(y::Quantity{S,D,U}, p::$PT{T,N}) where {S, D, U, T <: Quantity, N}
+                QT = Base.promote_op($op, typeof(y), T)
+                $PT{QT,N}($(op).(y, p.particles))
+            end
+
+            function Base.$f(p::$PT, y::FreeUnits)
+                $PT($(op).(p.particles, y))
+            end
         end
 
     end
diff --git a/test/test_unitful.jl b/test/test_unitful.jl
index 6441d08e..39914600 100644
--- a/test/test_unitful.jl
+++ b/test/test_unitful.jl
@@ -14,6 +14,7 @@ register_primitive(unitful_testfunction) # must be outside testset
 @testset "Unitful" begin
     @info "Testing Unitful"
 
+    PT = Particles
     for PT in (Particles, StaticParticles)
         p1 = PT(100, Uniform(-0.5,1.5)) * 1u"V"
         p2 = PT(100, Uniform(-0.5,1.5)) * u"V"
@@ -24,6 +25,20 @@ register_primitive(unitful_testfunction) # must be outside testset
 
         p3 = unitful_testfunction(p1)
         @test extrema(p3) == (0.0u"V", 1.0u"V")
+
+        @test (1 ± 0.5)u"m" * (1 ± 0)u"kg" ≈ (1 ± 0.5)u"kg*m"
+        @test (1 ± 0.5)u"m" * 1u"kg" ≈ (1 ± 0.5)u"kg*m"
+
+        @test (1 ± 0.5)u"m" / (1 ± 0)u"kg" ≈ (1 ± 0.5)u"m/kg"
+        @test (1 ± 0.5)u"m" / 1u"kg" ≈ (1 ± 0.5)u"m/kg"
+
+        @test (1 ± 0.5)u"m" + (1 ± 0)u"m" ≈ (2 ± 0.5)u"m"
+        @test (1 ± 0.5)u"m" + 1u"m" ≈ (2 ± 0.5)u"m"
+
+
+        @test 1u"m" * (1 ± 0.5)u"kg" ≈ (1 ± 0.5)u"kg*m"
+        @test 1u"m" / (1 ± 0.5)u"kg" ≈ (1 ± 0.5)u"m/kg"
+        @test 1u"m" + (1 ± 0.5)u"m" ≈ (2 ± 0.5)u"m"
     end
 
 end