diff --git a/Project.toml b/Project.toml
index 4806b6d..10f66ab 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,7 +8,7 @@ LorentzVectorBase = "592a752d-6533-4762-a71b-738712ea30ec"
 LorentzVectors = "3f54b04b-17fc-5cd4-9758-90c048d965e3"
 
 [compat]
-LorentzVectors = "^0.4"
+LorentzVectorBase = "^0.1"
 julia = "^1.6"
 
 [extras]
diff --git a/src/LorentzVectorHEP.jl b/src/LorentzVectorHEP.jl
index 104ed48..a2c6691 100644
--- a/src/LorentzVectorHEP.jl
+++ b/src/LorentzVectorHEP.jl
@@ -2,8 +2,7 @@ module LorentzVectorHEP
 
 using LorentzVectorBase
 import LorentzVectorBase: px, py, pz, energy, pt, rapidity, eta, phi, mt, mt2, pt, pt2, mass, mass2, azimuthal_angle
-
-using LorentzVectors # provides x, y, z, t
+import Base: +, -, *, /, ==, isapprox, ≈
 
 export LorentzVectorCyl, LorentzVector
 
@@ -12,6 +11,79 @@ export deltaphi, deltar, deltaeta
 export ΔR, Δϕ, Δη
 export fromPtEtaPhiE
 
+
+"""
+    LorentzVector(t, x, y, z)
+
+Lorentz 4-vector, as used in Special Relativity.
+
+The metric convention is g = diag(+1,-1,-1,-1). No distinction is made between
+co- and contravariant vectors.
+"""
+struct LorentzVector{T <: AbstractFloat}
+    t :: T
+    x :: T
+    y :: T
+    z :: T
+end
+
+"""
+    LorentzVector(t, x, y, z)
+
+Promoting constructors for LorentzVector{T}.
+"""
+LorentzVector(t, x, y, z) = LorentzVector(promote(t, x, y, z)...)
+LorentzVector(t::T, x::T, y::T, z::T) where {T <: Union{Integer, Rational, Irrational}} =
+    LorentzVector(float(t), x, y, z)
+    
+
+"""
+    dot(u, v)
+    u⋅v
+
+Inner product of 4-vectors, in the Minkowsky metric (+,-,-,-).
+"""
+function dot(u::LorentzVector, v::LorentzVector)
+    @fastmath u.t*v.t - u.x*v.x - u.y*v.y - u.z*v.z
+end
+
+
+function +(u::LorentzVector, v::LorentzVector)
+    @fastmath LorentzVector(u.t + v.t, u.x + v.x, u.y + v.y, u.z + v.z)
+end
+
+function -(u::LorentzVector)
+    @fastmath LorentzVector(-u.t, -u.x, -u.y, -u.z)
+end
+
+function -(u::LorentzVector, v::LorentzVector)
+    @fastmath u + (-v)
+end
+
+function *(λ::Number, u::LorentzVector)
+    @fastmath LorentzVector(λ*u.t, λ*u.x, λ*u.y, λ*u.z)
+end
+
+function *(u::LorentzVector, λ::Number)
+    @fastmath λ * u
+end
+
+function /(u::LorentzVector, λ::Number)
+    @fastmath u * (one(λ) / λ)
+end
+
+function ==(u::LorentzVector, v::LorentzVector)
+    u.t == v.t && u.x == v.x && u.y == v.y && u.z == v.z
+end
+
+function isapprox(u::LorentzVector, v::LorentzVector;
+                  atol::Real=0,
+                  rtol::Real=atol>0 ? 0 : √min(eps(typeof(u.x)), eps(typeof(v.x))))
+
+    err = max(abs(u.t - v.t), sqrt((u.x - v.x)^2 + (u.y - v.y)^2 + (u.z - v.z)^2))
+    err <= max(atol, rtol*max(abs(u.t), abs(v.t), sqrt(u.x^2 + u.y^2 + u.z^2), sqrt(v.x^2 + v.y^2 + v.z^2)))
+end
+
 include("cartesian.jl")
 
 include("cylindrical.jl")
diff --git a/test/runtests.jl b/test/runtests.jl
index 3bd8b43..b5d28e7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -19,7 +19,7 @@ using Test
     v1 = LorentzVectorCyl(43.71242f0, 1.4733887f0, 1.6855469f0, 0.10571289f0)
     v2 = LorentzVectorCyl(36.994347f0, 0.38684082f0, -1.3935547f0, 0.10571289f0)
     @test (v1+v2).mass == 92.55651f0
-    @test fast_mass(v1,v2) == 92.55651f0
+    @test fast_mass(v1,v2) ≈ 92.55651f0
     @test isapprox(deltar(v1,v2), 3.265188f0, atol=1e-6)
     @test isapprox(deltaphi(v1,v2), -3.0791016f0, atol=1e-6)
 
@@ -40,6 +40,24 @@ using Test
     @test phi(vcart2) ≈ -0.9884433806509134 atol=1e-9
     @test LorentzVectorHEP.phi02pi(vcart2) ≈ 5.294741926528673 atol=1e-9
 
+    @test vcart1 + vcart2 ≈ LorentzVector(vcart1.t + vcart2.t, vcart1.x + vcart2.x, vcart1.y + vcart2.y, vcart1.z + vcart2.z)
+    @test vcart1 + vcart2 == vcart2 + vcart1
+    @test -vcart1 == LorentzVector(-(vcart1.t), -(vcart1.x), -(vcart1.y), -(vcart1.z))
+    @test vcart1 - vcart2 ≈ LorentzVector(vcart1.t - vcart2.t, vcart1.x - vcart2.x, vcart1.y - vcart2.y, vcart1.z - vcart2.z)
+    @test vcart1 - vcart2 == -(vcart2 - vcart1)
+    
+    c = rand()
+    @test c*vcart1 ≈ LorentzVector(c*vcart1.t, c*vcart1.x, c*vcart1.y, c*vcart1.z)
+    @test -c*vcart1 ≈ LorentzVector(-c*vcart1.t, -c*vcart1.x, -c*vcart1.y, -c*vcart1.z)
+    @test vcart1*c == c*vcart1
+
+    if c == 0
+        c += 0.1
+    end
+
+    @test vcart1 / c ≈ vcart1 * (1/c)
+    
+
     @test deltaeta(vcart1, vcart2) ≈ 0.08825941862546584 atol=1e-9
     @test deltaphi(vcart1, vcart2) ≈ 3.0317366429628736 atol=1e-9