diff --git a/src/ConcreteOperations/linear_combination.jl b/src/ConcreteOperations/linear_combination.jl
index 2a637a35bc..16a9c4b0c0 100644
--- a/src/ConcreteOperations/linear_combination.jl
+++ b/src/ConcreteOperations/linear_combination.jl
@@ -1,9 +1,10 @@
 export linear_combination
 
 """
-    linear_combination(P1::SimpleSparsePolynomialZonotope, P2::SimpleSparsePolynomialZonotope)
+    linear_combination(P1::SimpleSparsePolynomialZonotope,
+                       P2::SimpleSparsePolynomialZonotope)
 
-Compute the linear combination of simple sparse polynomial zonotopes `P1` and `P2`.
+Compute the linear combination of two simple sparse polynomial zonotopes.
 
 ### Input
 
@@ -21,14 +22,20 @@ The linear combination of two sets ``P₁`` and ``P₂`` is defined as
 ```math
 \\{1/2(1+λ)p₁ + 1/2(1-λ)p₂ | p₁ ∈ P₁, p₂ ∈ P₂, λ ∈ [-1, 1]\\}.
 ```
+
+This method implements the algorithm described in Proposition 3.1.25 of [1].
+
+[1] N. Kochdumper. *Extensions of polynomial zonotopes and their application to
+verification of cyber-physical systems*. 2021.
 """
-function linear_combination(P1::SimpleSparsePolynomialZonotope, P2::SimpleSparsePolynomialZonotope)
+function linear_combination(P1::SimpleSparsePolynomialZonotope,
+                            P2::SimpleSparsePolynomialZonotope)
     c1, c2 = center(P1), center(P2)
     G1, G2 = genmat(P1), genmat(P2)
     E1, E2 = expmat(P1), expmat(P2)
 
-    c = 0.5 * (c1 + c2)
-    G = 0.5 * hcat(c1 - c2, G1, G1, G2, -G2)
+    c = (c1 + c2) / 2
+    G = hcat(c1 - c2, G1, G1, G2, -G2) / 2
 
     N = promote_type(eltype(E1), eltype(E2))
     E = hcat(vcat(zeros(N, nparams(P1) + nparams(P2)), one(N)),