diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
index 51f801fc..2a45e38f 100644
--- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
@@ -1261,7 +1261,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 throw new ArgumentNullException("leftSide");
             }
 
-            return (DenseMatrix)leftSide.Modulus(rightSide);
+            return (DenseMatrix)leftSide.Remainder(rightSide);
         }
 
         public override Cholesky<Complex> Cholesky()
diff --git a/src/Numerics/LinearAlgebra/Complex/Matrix.cs b/src/Numerics/LinearAlgebra/Complex/Matrix.cs
index d2ec6f96..ea61d027 100644
--- a/src/Numerics/LinearAlgebra/Complex/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Matrix.cs
@@ -455,31 +455,67 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         }
 
         /// <summary>
-        /// Pointwise modulus this matrix with another matrix and stores the result into the result matrix.
+        /// Pointwise canonical modulus, where the result has the sign of the divisor,
+        /// of this matrix with another matrix and stores the result into the result matrix.
         /// </summary>
         /// <param name="divisor">The pointwise denominator matrix to use</param>
         /// <param name="result">The result of the modulus.</param>
-        protected override void DoPointwiseModulus(Matrix<Complex> divisor, Matrix<Complex> result)
+        protected override sealed void DoPointwiseModulus(Matrix<Complex> divisor, Matrix<Complex> result)
         {
             throw new NotSupportedException();
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of this matrix with another matrix and stores the result into the result matrix.
+        /// </summary>
+        /// <param name="divisor">The pointwise denominator matrix to use</param>
+        /// <param name="result">The result of the modulus.</param>
+        protected override sealed void DoPointwiseRemainder(Matrix<Complex> divisor, Matrix<Complex> result)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
-        protected override void DoModulus(Complex divisor, Matrix<Complex> result)
+        protected override sealed void DoModulus(Complex divisor, Matrix<Complex> result)
         {
             throw new NotSupportedException();
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override sealed void DoModulusByThis(Complex dividend, Matrix<Complex> result)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
-        protected override void DoModulusByThis(Complex dividend, Matrix<Complex> result)
+        protected override sealed void DoRemainder(Complex divisor, Matrix<Complex> result)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override sealed void DoRemainderByThis(Complex dividend, Matrix<Complex> result)
         {
             throw new NotSupportedException();
         }
diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
index 453b6eb4..284a9276 100644
--- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
@@ -1400,7 +1400,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 throw new ArgumentNullException("leftSide");
             }
 
-            return (SparseMatrix)leftSide.Modulus(rightSide);
+            return (SparseMatrix)leftSide.Remainder(rightSide);
         }
 
         public override string ToTypeString()
diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
index 6819fc0b..f1ef1fe6 100644
--- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
@@ -1258,7 +1258,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 throw new ArgumentNullException("leftSide");
             }
 
-            return (DenseMatrix)leftSide.Modulus(rightSide);
+            return (DenseMatrix)leftSide.Remainder(rightSide);
         }
 
         public override Cholesky<Complex32> Cholesky()
diff --git a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
index abcf6d44..c4589e0b 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Matrix.cs
@@ -449,31 +449,67 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         }
 
         /// <summary>
-        /// Pointwise modulus this matrix with another matrix and stores the result into the result matrix.
+        /// Pointwise canonical modulus, where the result has the sign of the divisor,
+        /// of this matrix with another matrix and stores the result into the result matrix.
         /// </summary>
         /// <param name="divisor">The pointwise denominator matrix to use</param>
         /// <param name="result">The result of the modulus.</param>
-        protected override void DoPointwiseModulus(Matrix<Complex32> divisor, Matrix<Complex32> result)
+        protected override sealed void DoPointwiseModulus(Matrix<Complex32> divisor, Matrix<Complex32> result)
         {
             throw new NotSupportedException();
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of this matrix with another matrix and stores the result into the result matrix.
+        /// </summary>
+        /// <param name="divisor">The pointwise denominator matrix to use</param>
+        /// <param name="result">The result of the modulus.</param>
+        protected override sealed void DoPointwiseRemainder(Matrix<Complex32> divisor, Matrix<Complex32> result)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
-        protected override void DoModulus(Complex32 divisor, Matrix<Complex32> result)
+        protected override sealed void DoModulus(Complex32 divisor, Matrix<Complex32> result)
         {
             throw new NotSupportedException();
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override sealed void DoModulusByThis(Complex32 dividend, Matrix<Complex32> result)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
-        protected override void DoModulusByThis(Complex32 dividend, Matrix<Complex32> result)
+        protected override sealed void DoRemainder(Complex32 divisor, Matrix<Complex32> result)
+        {
+            throw new NotSupportedException();
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override sealed void DoRemainderByThis(Complex32 dividend, Matrix<Complex32> result)
         {
             throw new NotSupportedException();
         }
diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
index 6d2eb1a1..9efdcdd0 100644
--- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
@@ -1394,7 +1394,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 throw new ArgumentNullException("leftSide");
             }
 
-            return (SparseMatrix)leftSide.Modulus(rightSide);
+            return (SparseMatrix)leftSide.Remainder(rightSide);
         }
 
         public override string ToTypeString()
diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
index bfcefb99..c91c3192 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
@@ -859,9 +859,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
-        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
         protected override void DoModulus(double divisor, Matrix<double> result)
         {
@@ -878,20 +879,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             CommonParallel.For(0, _values.Length, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
                 {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] %= divisor;
-                    }
-                });
+                    v[i] = Euclid.Modulus(v[i], divisor);
+                }
+            });
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
-        /// <param name="result">Matrix to store the results in.</param>
+        /// <param name="result">A vector to store the results in.</param>
         protected override void DoModulusByThis(double dividend, Matrix<double> result)
         {
             var denseResult = result as DenseMatrix;
@@ -902,13 +904,68 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
                 {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = dividend%_values[i];
-                    }
-                });
+                    v[i] = Euclid.Modulus(dividend, _values[i]);
+                }
+            });
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoRemainder(double divisor, Matrix<double> result)
+        {
+            var denseResult = result as DenseMatrix;
+            if (denseResult == null)
+            {
+                base.DoRemainder(divisor, result);
+                return;
+            }
+
+            if (!ReferenceEquals(this, result))
+            {
+                CopyTo(result);
+            }
+
+            CommonParallel.For(0, _values.Length, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
+                {
+                    v[i] %= divisor;
+                }
+            });
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override void DoRemainderByThis(double dividend, Matrix<double> result)
+        {
+            var denseResult = result as DenseMatrix;
+            if (denseResult == null)
+            {
+                base.DoRemainderByThis(dividend, result);
+                return;
+            }
+
+            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
+                {
+                    v[i] = dividend%_values[i];
+                }
+            });
         }
 
         /// <summary>
@@ -1139,7 +1196,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("leftSide");
             }
 
-            return (DenseMatrix)leftSide.Modulus(rightSide);
+            return (DenseMatrix)leftSide.Remainder(rightSide);
         }
 
         public override Cholesky<double> Cholesky()
diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs
index 991004d7..5b13743a 100644
--- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs
+++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs
@@ -504,7 +504,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             var dense = result as DenseVector;
             if (dense == null)
             {
-                base.DoModulus(divisor, result);
+                base.DoRemainder(divisor, result);
             }
             else
             {
diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
index 3f6702dd..ced251a6 100644
--- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
@@ -882,7 +882,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
@@ -896,20 +897,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             CommonParallel.For(0, _data.Length, 4096, (a, b) =>
+            {
+                var r = diagonalResult._data;
+                for (var i = a; i < b; i++)
                 {
-                    var r = diagonalResult._data;
-                    for (var i = a; i < b; i++)
-                    {
-                        r[i] = _data[i]%divisor;
-                    }
-                });
+                    r[i] = Euclid.Modulus(_data[i], divisor);
+                }
+            });
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
-        /// <param name="result">Matrix to store the results in.</param>
+        /// <param name="result">A vector to store the results in.</param>
         protected override void DoModulusByThis(double dividend, Matrix<double> result)
         {
             var diagonalResult = result as DiagonalMatrix;
@@ -920,13 +922,63 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             CommonParallel.For(0, _data.Length, 4096, (a, b) =>
+            {
+                var r = diagonalResult._data;
+                for (var i = a; i < b; i++)
                 {
-                    var r = diagonalResult._data;
-                    for (var i = a; i < b; i++)
-                    {
-                        r[i] = dividend%_data[i];
-                    }
-                });
+                    r[i] = Euclid.Modulus(dividend, _data[i]);
+                }
+            });
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoRemainder(double divisor, Matrix<double> result)
+        {
+            var diagonalResult = result as DiagonalMatrix;
+            if (diagonalResult == null)
+            {
+                base.DoRemainder(divisor, result);
+                return;
+            }
+
+            CommonParallel.For(0, _data.Length, 4096, (a, b) =>
+            {
+                var r = diagonalResult._data;
+                for (var i = a; i < b; i++)
+                {
+                    r[i] = _data[i]%divisor;
+                }
+            });
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override void DoRemainderByThis(double dividend, Matrix<double> result)
+        {
+            var diagonalResult = result as DiagonalMatrix;
+            if (diagonalResult == null)
+            {
+                base.DoRemainderByThis(dividend, result);
+                return;
+            }
+
+            CommonParallel.For(0, _data.Length, 4096, (a, b) =>
+            {
+                var r = diagonalResult._data;
+                for (var i = a; i < b; i++)
+                {
+                    r[i] = dividend%_data[i];
+                }
+            });
         }
     }
 }
diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs
index 153d3583..220d00f2 100644
--- a/src/Numerics/LinearAlgebra/Double/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs
@@ -100,7 +100,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
         /// <summary>
         /// Returns the conjugate transpose of this matrix.
-        /// </summary>        
+        /// </summary>
         /// <returns>The conjugate transpose of this matrix.</returns>
         public override sealed Matrix<double> ConjugateTranspose()
         {
@@ -376,7 +376,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
@@ -386,17 +387,52 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             {
                 for (var column = 0; column < ColumnCount; column++)
                 {
-                    result.At(row, column, At(row, column)%divisor);
+                    result.At(row, column, Euclid.Modulus(At(row, column), divisor));
                 }
             }
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
-        /// <param name="result">Matrix to store the results in.</param>
+        /// <param name="result">A vector to store the results in.</param>
         protected override void DoModulusByThis(double dividend, Matrix<double> result)
+        {
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = 0; column < ColumnCount; column++)
+                {
+                    result.At(row, column, Euclid.Modulus(dividend, At(row, column)));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoRemainder(double divisor, Matrix<double> result)
+        {
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = 0; column < ColumnCount; column++)
+                {
+                    result.At(row, column, At(row, column)%divisor);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override void DoRemainderByThis(double dividend, Matrix<double> result)
         {
             for (var row = 0; row < RowCount; row++)
             {
@@ -440,11 +476,29 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         }
 
         /// <summary>
-        /// Pointwise modulus this matrix with another matrix and stores the result into the result matrix.
+        /// Pointwise canonical modulus, where the result has the sign of the divisor,
+        /// of this matrix with another matrix and stores the result into the result matrix.
         /// </summary>
         /// <param name="divisor">The pointwise denominator matrix to use</param>
         /// <param name="result">The result of the modulus.</param>
         protected override void DoPointwiseModulus(Matrix<double> divisor, Matrix<double> result)
+        {
+            for (var j = 0; j < ColumnCount; j++)
+            {
+                for (var i = 0; i < RowCount; i++)
+                {
+                    result.At(i, j, Euclid.Modulus(At(i, j), divisor.At(i, j)));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of this matrix with another matrix and stores the result into the result matrix.
+        /// </summary>
+        /// <param name="divisor">The pointwise denominator matrix to use</param>
+        /// <param name="result">The result of the modulus.</param>
+        protected override void DoPointwiseRemainder(Matrix<double> divisor, Matrix<double> result)
         {
             for (var j = 0; j < ColumnCount; j++)
             {
diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
index e3a9c384..49aa0bc8 100644
--- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
@@ -1129,9 +1129,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
-        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
         protected override void DoModulus(double divisor, Matrix<double> result)
         {
@@ -1147,6 +1148,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 CopyTo(result);
             }
 
+            var resultStorage = sparseResult._storage;
+            for (var index = 0; index < resultStorage.Values.Length; index++)
+            {
+                resultStorage.Values[index] = Euclid.Modulus(resultStorage.Values[index], divisor);
+            }
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoRemainder(double divisor, Matrix<double> result)
+        {
+            var sparseResult = result as SparseMatrix;
+            if (sparseResult == null)
+            {
+                base.DoRemainder(divisor, result);
+                return;
+            }
+
+            if (!ReferenceEquals(this, result))
+            {
+                CopyTo(result);
+            }
+
             var resultStorage = sparseResult._storage;
             for (var index = 0; index < resultStorage.Values.Length; index++)
             {
@@ -1420,7 +1448,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("leftSide");
             }
 
-            return (SparseMatrix)leftSide.Modulus(rightSide);
+            return (SparseMatrix)leftSide.Remainder(rightSide);
         }
 
         public override string ToTypeString()
diff --git a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs
index 1d778495..4d3d3e36 100644
--- a/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs
+++ b/src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs
@@ -178,19 +178,37 @@ namespace MathNet.Numerics.LinearAlgebra
         protected abstract void DoDivideByThis(T dividend, Matrix<T> result);
 
         /// <summary>
-        /// Computes the modulus for the given divisor each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
         protected abstract void DoModulus(T divisor, Matrix<T> result);
 
         /// <summary>
-        /// Computes the modulus for the given dividend for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
         /// <param name="result">A vector to store the results in.</param>
         protected abstract void DoModulusByThis(T dividend, Matrix<T> result);
 
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected abstract void DoRemainder(T divisor, Matrix<T> result);
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected abstract void DoRemainderByThis(T dividend, Matrix<T> result);
+
         /// <summary>
         /// Pointwise multiplies this matrix with another matrix and stores the result into the result matrix.
         /// </summary>
@@ -206,12 +224,21 @@ namespace MathNet.Numerics.LinearAlgebra
         protected abstract void DoPointwiseDivide(Matrix<T> divisor, Matrix<T> result);
 
         /// <summary>
-        /// Pointwise modulus this matrix with another matrix and stores the result into the result matrix.
+        /// Pointwise canonical modulus, where the result has the sign of the divisor,
+        /// of this matrix with another matrix and stores the result into the result matrix.
         /// </summary>
         /// <param name="divisor">The pointwise denominator matrix to use</param>
         /// <param name="result">The result of the modulus.</param>
         protected abstract void DoPointwiseModulus(Matrix<T> divisor, Matrix<T> result);
 
+        /// <summary>
+        /// Pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of this matrix with another matrix and stores the result into the result matrix.
+        /// </summary>
+        /// <param name="divisor">The pointwise denominator matrix to use</param>
+        /// <param name="result">The result of the modulus.</param>
+        protected abstract void DoPointwiseRemainder(Matrix<T> divisor, Matrix<T> result);
+
         /// <summary>
         /// Adds a scalar to each element of the matrix.
         /// </summary>
@@ -1020,7 +1047,8 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Computes the modulus (matrix % divisor) for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <returns>A matrix containing the results.</returns>
@@ -1032,7 +1060,8 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Computes the modulus (matrix % divisor) for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
@@ -1047,7 +1076,8 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Computes the modulus (dividend % matrix) for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
         /// <returns>A matrix containing the results.</returns>
@@ -1059,7 +1089,8 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Computes the modulus (dividend % matrix) for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
@@ -1073,6 +1104,64 @@ namespace MathNet.Numerics.LinearAlgebra
             DoModulusByThis(dividend, result);
         }
 
+        /// <summary>
+        /// Computes the remainder (matrix % divisor), where the result has the sign of the dividend,
+        /// for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <returns>A matrix containing the results.</returns>
+        public Matrix<T> Remainder(T divisor)
+        {
+            var result = Build.SameAs(this);
+            DoRemainder(divisor, result);
+            return result;
+        }
+
+        /// <summary>
+        /// Computes the remainder (matrix % divisor), where the result has the sign of the dividend,
+        /// for each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        public void Remainder(T divisor, Matrix<T> result)
+        {
+            if (ColumnCount != result.ColumnCount || RowCount != result.RowCount)
+            {
+                throw DimensionsDontMatch<ArgumentException>(this, result);
+            }
+
+            DoRemainder(divisor, result);
+        }
+
+        /// <summary>
+        /// Computes the remainder (dividend % matrix), where the result has the sign of the dividend,
+        /// for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <returns>A matrix containing the results.</returns>
+        public Matrix<T> RemainderByThis(T dividend)
+        {
+            var result = Build.SameAs(this);
+            DoRemainderByThis(dividend, result);
+            return result;
+        }
+
+        /// <summary>
+        /// Computes the remainder (dividend % matrix), where the result has the sign of the dividend,
+        /// for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        public void RemainderByThis(T dividend, Matrix<T> result)
+        {
+            if (ColumnCount != result.ColumnCount || RowCount != result.RowCount)
+            {
+                throw DimensionsDontMatch<ArgumentException>(this, result);
+            }
+
+            DoRemainderByThis(dividend, result);
+        }
+
         /// <summary>
         /// Pointwise multiplies this matrix with another matrix.
         /// </summary>
@@ -1144,11 +1233,11 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Pointwise modulus this matrix by another matrix.
+        /// Pointwise canonical modulus, where the result has the sign of the divisor,
+        /// of this matrix by another matrix.
         /// </summary>
         /// <param name="divisor">The pointwise denominator matrix to use.</param>
         /// <exception cref="ArgumentException">If this matrix and <paramref name="divisor"/> are not the same size.</exception>
-        /// <returns>A new matrix that is the pointwise modulus of this matrix and <paramref name="divisor"/>.</returns>
         public Matrix<T> PointwiseModulus(Matrix<T> divisor)
         {
             if (ColumnCount != divisor.ColumnCount || RowCount != divisor.RowCount)
@@ -1162,7 +1251,8 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Pointwise modulus this matrix by another matrix and stores the result into the result matrix.
+        /// Pointwise canonical modulus, where the result has the sign of the divisor,
+        /// of this matrix by another matrix and stores the result into the result matrix.
         /// </summary>
         /// <param name="divisor">The pointwise denominator matrix to use.</param>
         /// <param name="result">The matrix to store the result of the pointwise modulus.</param>
@@ -1178,6 +1268,42 @@ namespace MathNet.Numerics.LinearAlgebra
             DoPointwiseModulus(divisor, result);
         }
 
+        /// <summary>
+        /// Pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of this matrix by another matrix.
+        /// </summary>
+        /// <param name="divisor">The pointwise denominator matrix to use.</param>
+        /// <exception cref="ArgumentException">If this matrix and <paramref name="divisor"/> are not the same size.</exception>
+        public Matrix<T> PointwiseRemainder(Matrix<T> divisor)
+        {
+            if (ColumnCount != divisor.ColumnCount || RowCount != divisor.RowCount)
+            {
+                throw DimensionsDontMatch<ArgumentException>(this, divisor);
+            }
+
+            var result = Build.SameAs(this, divisor);
+            DoPointwiseRemainder(divisor, result);
+            return result;
+        }
+
+        /// <summary>
+        /// Pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of this matrix by another matrix and stores the result into the result matrix.
+        /// </summary>
+        /// <param name="divisor">The pointwise denominator matrix to use.</param>
+        /// <param name="result">The matrix to store the result of the pointwise remainder.</param>
+        /// <exception cref="ArgumentException">If this matrix and <paramref name="divisor"/> are not the same size.</exception>
+        /// <exception cref="ArgumentException">If this matrix and <paramref name="result"/> are not the same size.</exception>
+        public void PointwiseRemainder(Matrix<T> divisor, Matrix<T> result)
+        {
+            if (ColumnCount != result.ColumnCount || RowCount != result.RowCount || ColumnCount != divisor.ColumnCount || RowCount != divisor.RowCount)
+            {
+                throw DimensionsDontMatch<ArgumentException>(this, divisor, result);
+            }
+
+            DoPointwiseRemainder(divisor, result);
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>
diff --git a/src/Numerics/LinearAlgebra/Matrix.Operators.cs b/src/Numerics/LinearAlgebra/Matrix.Operators.cs
index 6b615b4b..239d9eab 100644
--- a/src/Numerics/LinearAlgebra/Matrix.Operators.cs
+++ b/src/Numerics/LinearAlgebra/Matrix.Operators.cs
@@ -239,7 +239,8 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         /// <summary>
-        /// Computes the modulus of each element of the matrix of the given divisor.
+        /// Computes the pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of each element of the matrix of the given divisor.
         /// </summary>
         /// <param name="dividend">The matrix whose elements we want to compute the modulus of.</param>
         /// <param name="divisor">The divisor to use.</param>
@@ -247,11 +248,12 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <exception cref="ArgumentNullException">If <paramref name="dividend"/> is <see langword="null" />.</exception>
         public static Matrix<T> operator %(Matrix<T> dividend, T divisor)
         {
-            return dividend.Modulus(divisor);
+            return dividend.Remainder(divisor);
         }
 
         /// <summary>
-        /// Computes the modulus of the given dividend of each element of the matrix.
+        /// Computes the pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of the given dividend of each element of the matrix.
         /// </summary>
         /// <param name="dividend">The dividend we want to compute the modulus of.</param>
         /// <param name="divisor">The matrix whose elements we want to use as divisor.</param>
@@ -259,20 +261,20 @@ namespace MathNet.Numerics.LinearAlgebra
         /// <exception cref="ArgumentNullException">If <paramref name="divisor"/> is <see langword="null" />.</exception>
         public static Matrix<T> operator %(T dividend, Matrix<T> divisor)
         {
-            return divisor.ModulusByThis(dividend);
+            return divisor.RemainderByThis(dividend);
         }
 
         /// <summary>
-        /// Computes the pointwise modulus of each element of two matrices.
+        /// Computes the pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of each element of two matrices.
         /// </summary>
-        /// <param name="dividend">The matrix whose elements we want to compute the modulus of.</param>
+        /// <param name="dividend">The matrix whose elements we want to compute the remainder of.</param>
         /// <param name="divisor">The divisor to use.</param>
-        /// <returns>The result of the calculation</returns>
         /// <exception cref="ArgumentException">If <paramref name="dividend"/> and <paramref name="divisor"/> are not the same size.</exception>
         /// <exception cref="ArgumentNullException">If <paramref name="dividend"/> is <see langword="null" />.</exception>
         public static Matrix<T> operator %(Matrix<T> dividend, Matrix<T> divisor)
         {
-            return dividend.PointwiseModulus(divisor);
+            return dividend.PointwiseRemainder(divisor);
         }
 
         [SpecialName]
@@ -288,9 +290,9 @@ namespace MathNet.Numerics.LinearAlgebra
         }
 
         [SpecialName]
-        public static Vector<T> op_DotPercent(Vector<T> dividend, Vector<T> divisor)
+        public static Matrix<T> op_DotPercent(Matrix<T> dividend, Matrix<T> divisor)
         {
-            return dividend.PointwiseModulus(divisor);
+            return dividend.PointwiseRemainder(divisor);
         }
     }
 }
diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
index 10b65a1c..0904e179 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
@@ -859,9 +859,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
-        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
         protected override void DoModulus(float divisor, Matrix<float> result)
         {
@@ -878,20 +879,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             CommonParallel.For(0, _values.Length, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
                 {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] %= divisor;
-                    }
-                });
+                    v[i] = Euclid.Modulus(v[i], divisor);
+                }
+            });
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
-        /// <param name="result">Matrix to store the results in.</param>
+        /// <param name="result">A vector to store the results in.</param>
         protected override void DoModulusByThis(float dividend, Matrix<float> result)
         {
             var denseResult = result as DenseMatrix;
@@ -902,13 +904,68 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
                 {
-                    var v = denseResult._values;
-                    for (int i = a; i < b; i++)
-                    {
-                        v[i] = dividend%_values[i];
-                    }
-                });
+                    v[i] = Euclid.Modulus(dividend, _values[i]);
+                }
+            });
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoRemainder(float divisor, Matrix<float> result)
+        {
+            var denseResult = result as DenseMatrix;
+            if (denseResult == null)
+            {
+                base.DoRemainder(divisor, result);
+                return;
+            }
+
+            if (!ReferenceEquals(this, result))
+            {
+                CopyTo(result);
+            }
+
+            CommonParallel.For(0, _values.Length, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
+                {
+                    v[i] %= divisor;
+                }
+            });
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override void DoRemainderByThis(float dividend, Matrix<float> result)
+        {
+            var denseResult = result as DenseMatrix;
+            if (denseResult == null)
+            {
+                base.DoRemainderByThis(dividend, result);
+                return;
+            }
+
+            CommonParallel.For(0, _values.Length, 4096, (a, b) =>
+            {
+                var v = denseResult._values;
+                for (int i = a; i < b; i++)
+                {
+                    v[i] = dividend%_values[i];
+                }
+            });
         }
 
         /// <summary>
@@ -1139,7 +1196,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 throw new ArgumentNullException("leftSide");
             }
 
-            return (DenseMatrix)leftSide.Modulus(rightSide);
+            return (DenseMatrix)leftSide.Remainder(rightSide);
         }
 
         public override Cholesky<float> Cholesky()
diff --git a/src/Numerics/LinearAlgebra/Single/DenseVector.cs b/src/Numerics/LinearAlgebra/Single/DenseVector.cs
index 6527b00f..e57131ad 100644
--- a/src/Numerics/LinearAlgebra/Single/DenseVector.cs
+++ b/src/Numerics/LinearAlgebra/Single/DenseVector.cs
@@ -493,7 +493,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             var dense = result as DenseVector;
             if (dense == null)
             {
-                base.DoModulus(divisor, result);
+                base.DoRemainder(divisor, result);
             }
             else
             {
diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
index 0c365220..88a4fe33 100644
--- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
@@ -882,7 +882,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
@@ -896,20 +897,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             CommonParallel.For(0, _data.Length, 4096, (a, b) =>
+            {
+                var r = diagonalResult._data;
+                for (var i = a; i < b; i++)
                 {
-                    var r = diagonalResult._data;
-                    for (var i = a; i < b; i++)
-                    {
-                        r[i] = _data[i]%divisor;
-                    }
-                });
+                    r[i] = Euclid.Modulus(_data[i], divisor);
+                }
+            });
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
-        /// <param name="result">Matrix to store the results in.</param>
+        /// <param name="result">A vector to store the results in.</param>
         protected override void DoModulusByThis(float dividend, Matrix<float> result)
         {
             var diagonalResult = result as DiagonalMatrix;
@@ -920,13 +922,63 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             CommonParallel.For(0, _data.Length, 4096, (a, b) =>
+            {
+                var r = diagonalResult._data;
+                for (var i = a; i < b; i++)
                 {
-                    var r = diagonalResult._data;
-                    for (var i = a; i < b; i++)
-                    {
-                        r[i] = dividend%_data[i];
-                    }
-                });
+                    r[i] = Euclid.Modulus(dividend, _data[i]);
+                }
+            });
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoRemainder(float divisor, Matrix<float> result)
+        {
+            var diagonalResult = result as DiagonalMatrix;
+            if (diagonalResult == null)
+            {
+                base.DoRemainder(divisor, result);
+                return;
+            }
+
+            CommonParallel.For(0, _data.Length, 4096, (a, b) =>
+            {
+                var r = diagonalResult._data;
+                for (var i = a; i < b; i++)
+                {
+                    r[i] = _data[i]%divisor;
+                }
+            });
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override void DoRemainderByThis(float dividend, Matrix<float> result)
+        {
+            var diagonalResult = result as DiagonalMatrix;
+            if (diagonalResult == null)
+            {
+                base.DoRemainderByThis(dividend, result);
+                return;
+            }
+
+            CommonParallel.For(0, _data.Length, 4096, (a, b) =>
+            {
+                var r = diagonalResult._data;
+                for (var i = a; i < b; i++)
+                {
+                    r[i] = dividend%_data[i];
+                }
+            });
         }
     }
 }
diff --git a/src/Numerics/LinearAlgebra/Single/Matrix.cs b/src/Numerics/LinearAlgebra/Single/Matrix.cs
index de324d48..bc3b8cd5 100644
--- a/src/Numerics/LinearAlgebra/Single/Matrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/Matrix.cs
@@ -100,7 +100,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
         /// <summary>
         /// Returns the conjugate transpose of this matrix.
-        /// </summary>        
+        /// </summary>
         /// <returns>The conjugate transpose of this matrix.</returns>
         public override sealed Matrix<float> ConjugateTranspose()
         {
@@ -347,7 +347,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
         /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
@@ -357,17 +358,52 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             {
                 for (var column = 0; column < ColumnCount; column++)
                 {
-                    result.At(row, column, At(row, column)%divisor);
+                    result.At(row, column, Euclid.Modulus(At(row, column), divisor));
                 }
             }
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given dividend for each element of the matrix.
         /// </summary>
         /// <param name="dividend">The scalar numerator to use.</param>
-        /// <param name="result">Matrix to store the results in.</param>
+        /// <param name="result">A vector to store the results in.</param>
         protected override void DoModulusByThis(float dividend, Matrix<float> result)
+        {
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = 0; column < ColumnCount; column++)
+                {
+                    result.At(row, column, Euclid.Modulus(dividend, At(row, column)));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoRemainder(float divisor, Matrix<float> result)
+        {
+            for (var row = 0; row < RowCount; row++)
+            {
+                for (var column = 0; column < ColumnCount; column++)
+                {
+                    result.At(row, column, At(row, column)%divisor);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given dividend for each element of the matrix.
+        /// </summary>
+        /// <param name="dividend">The scalar numerator to use.</param>
+        /// <param name="result">A vector to store the results in.</param>
+        protected override void DoRemainderByThis(float dividend, Matrix<float> result)
         {
             for (var row = 0; row < RowCount; row++)
             {
@@ -440,11 +476,29 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         }
 
         /// <summary>
-        /// Pointwise modulus this matrix with another matrix and stores the result into the result matrix.
+        /// Pointwise canonical modulus, where the result has the sign of the divisor,
+        /// of this matrix with another matrix and stores the result into the result matrix.
         /// </summary>
         /// <param name="divisor">The pointwise denominator matrix to use</param>
         /// <param name="result">The result of the modulus.</param>
         protected override void DoPointwiseModulus(Matrix<float> divisor, Matrix<float> result)
+        {
+            for (var j = 0; j < ColumnCount; j++)
+            {
+                for (var i = 0; i < RowCount; i++)
+                {
+                    result.At(i, j, Euclid.Modulus(At(i, j), divisor.At(i, j)));
+                }
+            }
+        }
+
+        /// <summary>
+        /// Pointwise remainder (% operator), where the result has the sign of the dividend,
+        /// of this matrix with another matrix and stores the result into the result matrix.
+        /// </summary>
+        /// <param name="divisor">The pointwise denominator matrix to use</param>
+        /// <param name="result">The result of the modulus.</param>
+        protected override void DoPointwiseRemainder(Matrix<float> divisor, Matrix<float> result)
         {
             for (var j = 0; j < ColumnCount; j++)
             {
diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
index 403fd6f7..97b19ec4 100644
--- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
+++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
@@ -1136,9 +1136,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         }
 
         /// <summary>
-        /// Computes the modulus for each element of the matrix.
+        /// Computes the canonical modulus, where the result has the sign of the divisor,
+        /// for the given divisor each element of the matrix.
         /// </summary>
-        /// <param name="divisor">The divisor to use.</param>
+        /// <param name="divisor">The scalar denominator to use.</param>
         /// <param name="result">Matrix to store the results in.</param>
         protected override void DoModulus(float divisor, Matrix<float> result)
         {
@@ -1154,6 +1155,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 CopyTo(result);
             }
 
+            var resultStorage = sparseResult._storage;
+            for (var index = 0; index < resultStorage.Values.Length; index++)
+            {
+                resultStorage.Values[index] = Euclid.Modulus(resultStorage.Values[index], divisor);
+            }
+        }
+
+        /// <summary>
+        /// Computes the remainder (% operator), where the result has the sign of the dividend,
+        /// for the given divisor each element of the matrix.
+        /// </summary>
+        /// <param name="divisor">The scalar denominator to use.</param>
+        /// <param name="result">Matrix to store the results in.</param>
+        protected override void DoRemainder(float divisor, Matrix<float> result)
+        {
+            var sparseResult = result as SparseMatrix;
+            if (sparseResult == null)
+            {
+                base.DoRemainder(divisor, result);
+                return;
+            }
+
+            if (!ReferenceEquals(this, result))
+            {
+                CopyTo(result);
+            }
+
             var resultStorage = sparseResult._storage;
             for (var index = 0; index < resultStorage.Values.Length; index++)
             {
@@ -1427,7 +1455,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 throw new ArgumentNullException("leftSide");
             }
 
-            return (SparseMatrix)leftSide.Modulus(rightSide);
+            return (SparseMatrix)leftSide.Remainder(rightSide);
         }
 
         public override string ToTypeString()
diff --git a/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs b/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs
index 840d8e23..e5c16f69 100644
--- a/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs
+++ b/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs
@@ -767,7 +767,6 @@ namespace MathNet.Numerics.LinearAlgebra
         /// of this vector with another vector.
         /// </summary>
         /// <param name="divisor">The pointwise denominator vector to use.</param>
-        /// <returns>A new vector which is the pointwise modulus of the two vectors.</returns>
         /// <exception cref="ArgumentException">If this vector and <paramref name="divisor"/> are not the same size.</exception>
         public Vector<T> PointwiseModulus(Vector<T> divisor)
         {
@@ -808,7 +807,6 @@ namespace MathNet.Numerics.LinearAlgebra
         /// of this vector with another vector.
         /// </summary>
         /// <param name="divisor">The pointwise denominator vector to use.</param>
-        /// <returns>A new vector which is the pointwise remainder of the two vectors.</returns>
         /// <exception cref="ArgumentException">If this vector and <paramref name="divisor"/> are not the same size.</exception>
         public Vector<T> PointwiseRemainder(Vector<T> divisor)
         {
diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs
index 80fb7e8c..caa58e90 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs
@@ -54,7 +54,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrix.ColumnCount; j++)
                 {
-                    Assert.AreEqual(matrix[i, j] * scalar, clone[i, j]);
+                    Assert.AreEqual(matrix[i, j]*scalar, clone[i, j]);
                 }
             }
         }
@@ -67,14 +67,14 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         {
             var matrix = TestMatrices["Singular3x3"];
             var x = new DenseVector(new[] { 1.0, 2.0, 3.0 });
-            var y = matrix * x;
+            var y = matrix*x;
 
             Assert.AreEqual(matrix.RowCount, y.Count);
 
             for (var i = 0; i < matrix.RowCount; i++)
             {
                 var ar = matrix.Row(i);
-                var dot = ar * x;
+                var dot = ar*x;
                 Assert.AreEqual(dot, y[i]);
             }
         }
@@ -93,7 +93,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             for (var i = 0; i < matrix.RowCount; i++)
             {
                 var ar = matrix.Row(i);
-                var dot = ar * x;
+                var dot = ar*x;
                 Assert.AreEqual(dot, y[i]);
             }
         }
@@ -115,7 +115,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             for (var i = 0; i < matrix.RowCount; i++)
             {
                 var ar = matrix.Row(i);
-                var dot = ar * y;
+                var dot = ar*y;
                 Assert.AreEqual(dot, x[i]);
             }
         }
@@ -142,13 +142,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         public void CanOperatorLeftMultiplyWithScalar(double scalar)
         {
             var matrix = TestMatrices["Singular3x3"];
-            var clone = scalar * matrix;
+            var clone = scalar*matrix;
 
             for (var i = 0; i < matrix.RowCount; i++)
             {
                 for (var j = 0; j < matrix.ColumnCount; j++)
                 {
-                    Assert.AreEqual(scalar * matrix[i, j], clone[i, j]);
+                    Assert.AreEqual(scalar*matrix[i, j], clone[i, j]);
                 }
             }
         }
@@ -163,13 +163,13 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         public void CanOperatorRightMultiplyWithScalar(double scalar)
         {
             var matrix = TestMatrices["Singular3x3"];
-            var clone = matrix * scalar;
+            var clone = matrix*scalar;
 
             for (var i = 0; i < matrix.RowCount; i++)
             {
                 for (var j = 0; j < matrix.ColumnCount; j++)
                 {
-                    Assert.AreEqual(matrix[i, j] * scalar, clone[i, j]);
+                    Assert.AreEqual(matrix[i, j]*scalar, clone[i, j]);
                 }
             }
         }
@@ -191,7 +191,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrix.ColumnCount; j++)
                 {
-                    Assert.AreEqual(matrix[i, j] * scalar, result[i, j]);
+                    Assert.AreEqual(matrix[i, j]*scalar, result[i, j]);
                 }
             }
         }
@@ -410,7 +410,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         {
             var matrixA = TestMatrices[nameA];
             var matrixB = TestMatrices[nameB];
-            var matrixC = matrixA * matrixB;
+            var matrixC = matrixA*matrixB;
 
             Assert.AreEqual(matrixC.RowCount, matrixA.RowCount);
             Assert.AreEqual(matrixC.ColumnCount, matrixB.ColumnCount);
@@ -419,7 +419,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrixC.ColumnCount; j++)
                 {
-                    AssertHelpers.AlmostEqual(matrixA.Row(i) * matrixB.Column(j), matrixC[i, j], 15);
+                    AssertHelpers.AlmostEqual(matrixA.Row(i)*matrixB.Column(j), matrixC[i, j], 15);
                 }
             }
         }
@@ -445,7 +445,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrixC.ColumnCount; j++)
                 {
-                    AssertHelpers.AlmostEqual(matrixA.Row(i) * matrixB.Row(j), matrixC[i, j], 15);
+                    AssertHelpers.AlmostEqual(matrixA.Row(i)*matrixB.Row(j), matrixC[i, j], 15);
                 }
             }
         }
@@ -479,7 +479,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrixC.ColumnCount; j++)
                 {
-                    AssertHelpers.AlmostEqual(matrixA.Row(i) * matrixB.Row(j), matrixC[i, j], 15);
+                    AssertHelpers.AlmostEqual(matrixA.Row(i)*matrixB.Row(j), matrixC[i, j], 15);
                 }
             }
         }
@@ -517,7 +517,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrixC.ColumnCount; j++)
                 {
-                    AssertHelpers.AlmostEqual(matrixA.Row(i) * matrixB.Row(j), matrixC[i, j], 15);
+                    AssertHelpers.AlmostEqual(matrixA.Row(i)*matrixB.Row(j), matrixC[i, j], 15);
                 }
             }
         }
@@ -530,7 +530,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         {
             var matrix = TestMatrices["Singular3x3"];
             var other = TestMatrices["Wide2x3"];
-            Assert.Throws<ArgumentException>(() => { var result = matrix * other; });
+            Assert.Throws<ArgumentException>(() => { var result = matrix*other; });
         }
 
         /// <summary>
@@ -557,7 +557,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrixC.ColumnCount; j++)
                 {
-                    AssertHelpers.AlmostEqual(matrixA.Row(i) * matrixB.Column(j), matrixC[i, j], 15);
+                    AssertHelpers.AlmostEqual(matrixA.Row(i)*matrixB.Column(j), matrixC[i, j], 15);
                 }
             }
         }
@@ -577,7 +577,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             for (var j = 0; j < matrix.ColumnCount; j++)
             {
                 var ar = matrix.Column(j);
-                var dot = ar * x;
+                var dot = ar*x;
                 Assert.AreEqual(dot, y[j]);
             }
         }
@@ -596,7 +596,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             for (var j = 0; j < matrix.ColumnCount; j++)
             {
                 var ar = matrix.Column(j);
-                var dot = ar * x;
+                var dot = ar*x;
                 Assert.AreEqual(dot, y[j]);
             }
         }
@@ -618,7 +618,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             for (var j = 0; j < matrix.ColumnCount; j++)
             {
                 var ar = matrix.Column(j);
-                var dot = ar * y;
+                var dot = ar*y;
                 Assert.AreEqual(dot, x[j]);
             }
         }
@@ -656,7 +656,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrixC.ColumnCount; j++)
                 {
-                    AssertHelpers.AlmostEqual(matrixA.Column(i) * matrixB.Column(j), matrixC[i, j], 15);
+                    AssertHelpers.AlmostEqual(matrixA.Column(i)*matrixB.Column(j), matrixC[i, j], 15);
                 }
             }
         }
@@ -694,7 +694,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < matrixC.ColumnCount; j++)
                 {
-                    AssertHelpers.AlmostEqual(matrixA.Column(i) * matrixB.Column(j), matrixC[i, j], 15);
+                    AssertHelpers.AlmostEqual(matrixA.Column(i)*matrixB.Column(j), matrixC[i, j], 15);
                 }
             }
         }
@@ -779,7 +779,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
         {
             var matrixA = TestMatrices["Wide2x3"];
             var matrixB = TestMatrices["Square3x3"];
-            var result = CreateMatrix(matrixA.RowCount * matrixB.RowCount, matrixA.ColumnCount * matrixB.ColumnCount);
+            var result = CreateMatrix(matrixA.RowCount*matrixB.RowCount, matrixA.ColumnCount*matrixB.ColumnCount);
             matrixA.KroneckerProduct(matrixB, result);
             for (var i = 0; i < matrixA.RowCount; i++)
             {
@@ -789,7 +789,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
                     {
                         for (var jj = 0; jj < matrixB.ColumnCount; jj++)
                         {
-                            Assert.AreEqual(result[(i * matrixB.RowCount) + ii, (j * matrixB.ColumnCount) + jj], matrixA[i, j] * matrixB[ii, jj]);
+                            Assert.AreEqual(result[(i*matrixB.RowCount) + ii, (j*matrixB.ColumnCount) + jj], matrixA[i, j]*matrixB[ii, jj]);
                         }
                     }
                 }
@@ -813,7 +813,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
                     {
                         for (var jj = 0; jj < matrixB.ColumnCount; jj++)
                         {
-                            Assert.AreEqual(result[(i * matrixB.RowCount) + ii, (j * matrixB.ColumnCount) + jj], matrixA[i, j] * matrixB[ii, jj]);
+                            Assert.AreEqual(result[(i*matrixB.RowCount) + ii, (j*matrixB.ColumnCount) + jj], matrixA[i, j]*matrixB[ii, jj]);
                         }
                     }
                 }
@@ -886,7 +886,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
                 {
                     for (var j = 0; j < data.ColumnCount; j++)
                     {
-                        Assert.AreEqual(data[i, j] * other[i, j], result[i, j]);
+                        Assert.AreEqual(data[i, j]*other[i, j], result[i, j]);
                     }
                 }
 
@@ -895,7 +895,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
                 {
                     for (var j = 0; j < data.ColumnCount; j++)
                     {
-                        Assert.AreEqual(data[i, j] * other[i, j], result[i, j]);
+                        Assert.AreEqual(data[i, j]*other[i, j], result[i, j]);
                     }
                 }
             }
@@ -939,7 +939,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < data.ColumnCount; j++)
                 {
-                    Assert.AreEqual(data[i, j] / other[i, j], result[i, j]);
+                    Assert.AreEqual(data[i, j]/other[i, j], result[i, j]);
                 }
             }
 
@@ -948,7 +948,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             {
                 for (var j = 0; j < data.ColumnCount; j++)
                 {
-                    Assert.AreEqual(data[i, j] / other[i, j], result[i, j]);
+                    Assert.AreEqual(data[i, j]/other[i, j], result[i, j]);
                 }
             }
         }
@@ -1009,74 +1009,104 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             Assert.Throws<ArgumentException>(() => matrix.Trace());
         }
 
-        /// <summary>
-        /// Can compute the modules of each element of vector.
-        /// </summary>
         [Test]
-        public void CanComputeModulus()
+        public void CanComputeRemainderUsingOperator()
         {
             var matrix = TestMatrices["Square3x3"];
-            var mod = matrix.Modulus(3.2);
+            var mod = matrix%(-3.2);
             for (var row = 0; row < matrix.RowCount; row++)
             {
                 for (var column = 0; column < matrix.ColumnCount; column++)
                 {
-                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2, mod[row, column], 14);
+                    AssertHelpers.AlmostEqual(Euclid.Remainder(matrix[row, column], -3.2), mod[row, column], 14);
                 }
             }
         }
 
-        /// <summary>
-        /// Can compute the modules of each element of vector using a result vector.
-        /// </summary>
         [Test]
-        public void CanComputeModulusUsingResultVector()
+        public void CanComputeRemainder()
         {
             var matrix = TestMatrices["Square3x3"];
-            var mod = CreateMatrix(matrix.RowCount, matrix.ColumnCount);
-            matrix.Modulus(3.2, mod);
+            var mod = matrix.Remainder(-3.2);
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(Euclid.Remainder(matrix[row, column], -3.2), mod[row, column], 14);
+                }
+            }
+        }
 
+        [Test]
+        public void CanComputeRemainderUsingResultVector()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = CreateMatrix(matrix.RowCount, matrix.ColumnCount);
+            matrix.Remainder(-3.2, mod);
             for (var row = 0; row < matrix.RowCount; row++)
             {
                 for (var column = 0; column < matrix.ColumnCount; column++)
                 {
-                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2, mod[row, column], 14);
+                    AssertHelpers.AlmostEqual(Euclid.Remainder(matrix[row, column], -3.2), mod[row, column], 14);
                 }
             }
         }
 
-        /// <summary>
-        /// Can compute the modules of each element of vector using a result vector.
-        /// </summary>
         [Test]
-        public void CanComputeModulusUsingSameResultVector()
+        public void CanComputeRemainderUsingSameResultVector()
         {
             var matrix = TestMatrices["Square3x3"].Clone();
-            matrix.Modulus(3.2, matrix);
+            matrix.Remainder(-3.2, matrix);
             var data = TestMatrices["Square3x3"];
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(Euclid.Remainder(data[row, column], -3.2), matrix[row, column], 14);
+                }
+            }
+        }
 
+        [Test]
+        public void CanComputeModulus()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = matrix.Modulus(-3.2);
             for (var row = 0; row < matrix.RowCount; row++)
             {
                 for (var column = 0; column < matrix.ColumnCount; column++)
                 {
-                    AssertHelpers.AlmostEqual(data[row, column] % 3.2, matrix[row, column], 14);
+                    AssertHelpers.AlmostEqual(Euclid.Modulus(matrix[row, column], -3.2), mod[row, column], 14);
                 }
             }
         }
 
-        /// <summary>
-        /// Can compute the modules of each element of vector using the operator %.
-        /// </summary>
         [Test]
-        public void CanComputeModulusUsingOperator()
+        public void CanComputeModulusUsingResultVector()
         {
             var matrix = TestMatrices["Square3x3"];
-            var mod = matrix % 3.2;
+            var mod = CreateMatrix(matrix.RowCount, matrix.ColumnCount);
+            matrix.Modulus(-3.2, mod);
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(Euclid.Modulus(matrix[row, column], -3.2), mod[row, column], 14);
+                }
+            }
+        }
+
+        [Test]
+        public void CanComputeModulusUsingSameResultVector()
+        {
+            var matrix = TestMatrices["Square3x3"].Clone();
+            matrix.Modulus(-3.2, matrix);
+            var data = TestMatrices["Square3x3"];
             for (var row = 0; row < matrix.RowCount; row++)
             {
                 for (var column = 0; column < matrix.ColumnCount; column++)
                 {
-                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2, mod[row, column], 14);
+                    AssertHelpers.AlmostEqual(Euclid.Modulus(data[row, column], -3.2), matrix[row, column], 14);
                 }
             }
         }
diff --git a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs
index 7dd60986..86e94502 100644
--- a/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs
+++ b/src/UnitTests/LinearAlgebraTests/Double/VectorTests.Arithmetic.cs
@@ -790,60 +790,82 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double
             }
         }
 
-        /// <summary>
-        /// Can compute the remainder of each element of vector.
-        /// </summary>
+        [Test]
+        public void CanComputeRemainderUsingOperator()
+        {
+            var vector = CreateVector(Data);
+            var mod = vector % (-4.5);
+            for (var index = 0; index < Data.Length; index++)
+            {
+                AssertHelpers.AlmostEqualRelative(Euclid.Remainder(Data[index], -4.5), mod[index], 14);
+            }
+        }
+
         [Test]
         public void CanComputeRemainder()
         {
             var vector = CreateVector(Data);
-            var mod = vector.Remainder(3.2);
+            var mod = vector.Remainder(-3.2);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqualRelative(Data[index] % 3.2, mod[index], 14);
+                AssertHelpers.AlmostEqualRelative(Euclid.Remainder(Data[index], -3.2), mod[index], 14);
             }
         }
 
-        /// <summary>
-        /// Can compute the remainder of each element of vector using a result vector.
-        /// </summary>
         [Test]
         public void CanComputeRemainderUsingResultVector()
         {
             var vector = CreateVector(Data);
             var mod = CreateVector(vector.Count);
-            vector.Remainder(3.2, mod);
+            vector.Remainder(-3.2, mod);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqualRelative(Data[index] % 3.2, mod[index], 14);
+                AssertHelpers.AlmostEqualRelative(Euclid.Remainder(Data[index], -3.2), mod[index], 14);
             }
         }
 
-        /// <summary>
-        /// Can compute the remainder of each element of vector using a result vector.
-        /// </summary>
         [Test]
         public void CanComputeRemainderUsingSameResultVector()
         {
             var vector = CreateVector(Data);
-            vector.Remainder(3.2, vector);
+            vector.Remainder(-3.2, vector);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqualRelative(Data[index] % 3.2, vector[index], 14);
+                AssertHelpers.AlmostEqualRelative(Euclid.Remainder(Data[index], -3.2), vector[index], 14);
             }
         }
 
-        /// <summary>
-        /// Can compute the remainder of each element of vector using the operator %.
-        /// </summary>
         [Test]
-        public void CanComputeRemainderUsingOperator()
+        public void CanComputeModulus()
+        {
+            var vector = CreateVector(Data);
+            var mod = vector.Modulus(-3.2);
+            for (var index = 0; index < Data.Length; index++)
+            {
+                AssertHelpers.AlmostEqualRelative(Euclid.Modulus(Data[index], -3.2), mod[index], 14);
+            }
+        }
+
+        [Test]
+        public void CanComputeModulusUsingResultVector()
+        {
+            var vector = CreateVector(Data);
+            var mod = CreateVector(vector.Count);
+            vector.Modulus(-3.2, mod);
+            for (var index = 0; index < Data.Length; index++)
+            {
+                AssertHelpers.AlmostEqualRelative(Euclid.Modulus(Data[index], -3.2), mod[index], 14);
+            }
+        }
+
+        [Test]
+        public void CanComputeModulusUsingSameResultVector()
         {
             var vector = CreateVector(Data);
-            var mod = vector % 4.5;
+            vector.Modulus(-3.2, vector);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqualRelative(Data[index] % 4.5, mod[index], 14);
+                AssertHelpers.AlmostEqualRelative(Euclid.Modulus(Data[index], -3.2), vector[index], 14);
             }
         }
     }
diff --git a/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs
index 00ba92b8..eafb1d8c 100644
--- a/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs
+++ b/src/UnitTests/LinearAlgebraTests/Single/MatrixTests.Arithmetic.cs
@@ -1003,74 +1003,104 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
             Assert.Throws<ArgumentException>(() => matrix.Trace());
         }
 
-        /// <summary>
-        /// Can compute the modules of each element of vector.
-        /// </summary>
         [Test]
-        public void CanComputeModulus()
+        public void CanComputeRemainderUsingOperator()
         {
             var matrix = TestMatrices["Square3x3"];
-            var mod = matrix.Modulus(3.2f);
+            var mod = matrix%(-3.2f);
             for (var row = 0; row < matrix.RowCount; row++)
             {
                 for (var column = 0; column < matrix.ColumnCount; column++)
                 {
-                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2f, mod[row, column], 6);
+                    AssertHelpers.AlmostEqual(Euclid.Remainder(matrix[row, column], -3.2f), mod[row, column], 14);
                 }
             }
         }
 
-        /// <summary>
-        /// Can compute the modules of each element of vector using a result vector.
-        /// </summary>
         [Test]
-        public void CanComputeModulusUsingResultVector()
+        public void CanComputeRemainder()
         {
             var matrix = TestMatrices["Square3x3"];
-            var mod = CreateMatrix(matrix.RowCount, matrix.ColumnCount);
-            matrix.Modulus(3.2f, mod);
+            var mod = matrix.Remainder(-3.2f);
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(Euclid.Remainder(matrix[row, column], -3.2f), mod[row, column], 14);
+                }
+            }
+        }
 
+        [Test]
+        public void CanComputeRemainderUsingResultVector()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = CreateMatrix(matrix.RowCount, matrix.ColumnCount);
+            matrix.Remainder(-3.2f, mod);
             for (var row = 0; row < matrix.RowCount; row++)
             {
                 for (var column = 0; column < matrix.ColumnCount; column++)
                 {
-                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2f, mod[row, column], 6);
+                    AssertHelpers.AlmostEqual(Euclid.Remainder(matrix[row, column], -3.2f), mod[row, column], 14);
                 }
             }
         }
 
-        /// <summary>
-        /// Can compute the modules of each element of vector using a result vector.
-        /// </summary>
         [Test]
-        public void CanComputeModulusUsingSameResultVector()
+        public void CanComputeRemainderUsingSameResultVector()
         {
             var matrix = TestMatrices["Square3x3"].Clone();
-            matrix.Modulus(3.2f, matrix);
+            matrix.Remainder(-3.2f, matrix);
             var data = TestMatrices["Square3x3"];
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(Euclid.Remainder(data[row, column], -3.2f), matrix[row, column], 14);
+                }
+            }
+        }
 
+        [Test]
+        public void CanComputeModulus()
+        {
+            var matrix = TestMatrices["Square3x3"];
+            var mod = matrix.Modulus(-3.2f);
             for (var row = 0; row < matrix.RowCount; row++)
             {
                 for (var column = 0; column < matrix.ColumnCount; column++)
                 {
-                    AssertHelpers.AlmostEqual(data[row, column] % 3.2f, matrix[row, column], 6);
+                    AssertHelpers.AlmostEqual(Euclid.Modulus(matrix[row, column], -3.2f), mod[row, column], 14);
                 }
             }
         }
 
-        /// <summary>
-        /// Can compute the modules of each element of vector using the operator %.
-        /// </summary>
         [Test]
-        public void CanComputeModulusUsingOperator()
+        public void CanComputeModulusUsingResultVector()
         {
             var matrix = TestMatrices["Square3x3"];
-            var mod = matrix % 3.2f;
+            var mod = CreateMatrix(matrix.RowCount, matrix.ColumnCount);
+            matrix.Modulus(-3.2f, mod);
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var column = 0; column < matrix.ColumnCount; column++)
+                {
+                    AssertHelpers.AlmostEqual(Euclid.Modulus(matrix[row, column], -3.2f), mod[row, column], 14);
+                }
+            }
+        }
+
+        [Test]
+        public void CanComputeModulusUsingSameResultVector()
+        {
+            var matrix = TestMatrices["Square3x3"].Clone();
+            matrix.Modulus(-3.2f, matrix);
+            var data = TestMatrices["Square3x3"];
             for (var row = 0; row < matrix.RowCount; row++)
             {
                 for (var column = 0; column < matrix.ColumnCount; column++)
                 {
-                    AssertHelpers.AlmostEqual(matrix[row, column] % 3.2f, mod[row, column], 6);
+                    AssertHelpers.AlmostEqual(Euclid.Modulus(data[row, column], -3.2f), matrix[row, column], 14);
                 }
             }
         }
diff --git a/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Arithmetic.cs
index 418899fb..ddfeec60 100644
--- a/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Arithmetic.cs
+++ b/src/UnitTests/LinearAlgebraTests/Single/VectorTests.Arithmetic.cs
@@ -790,60 +790,82 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single
             }
         }
 
-        /// <summary>
-        /// Can compute the remainder of each element of vector.
-        /// </summary>
+        [Test]
+        public void CanComputeRemainderUsingOperator()
+        {
+            var vector = CreateVector(Data);
+            var mod = vector % (-4.5f);
+            for (var index = 0; index < Data.Length; index++)
+            {
+                AssertHelpers.AlmostEqualRelative(Euclid.Remainder(Data[index], -4.5f), mod[index], 14);
+            }
+        }
+
         [Test]
         public void CanComputeRemainder()
         {
             var vector = CreateVector(Data);
-            var mod = vector.Remainder(3.2f);
+            var mod = vector.Remainder(-3.2f);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqualRelative(Data[index] % 3.2f, mod[index], 14);
+                AssertHelpers.AlmostEqualRelative(Euclid.Remainder(Data[index], -3.2f), mod[index], 14);
             }
         }
 
-        /// <summary>
-        /// Can compute the remainder of each element of vector using a result vector.
-        /// </summary>
         [Test]
         public void CanComputeRemainderUsingResultVector()
         {
             var vector = CreateVector(Data);
             var mod = CreateVector(vector.Count);
-            vector.Remainder(3.2f, mod);
+            vector.Remainder(-3.2f, mod);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqualRelative(Data[index] % 3.2f, mod[index], 14);
+                AssertHelpers.AlmostEqualRelative(Euclid.Remainder(Data[index], -3.2f), mod[index], 14);
             }
         }
 
-        /// <summary>
-        /// Can compute the remainder of each element of vector using a result vector.
-        /// </summary>
         [Test]
         public void CanComputeRemainderUsingSameResultVector()
         {
             var vector = CreateVector(Data);
-            vector.Remainder(3.2f, vector);
+            vector.Remainder(-3.2f, vector);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqualRelative(Data[index] % 3.2f, vector[index], 14);
+                AssertHelpers.AlmostEqualRelative(Euclid.Remainder(Data[index], -3.2f), vector[index], 14);
             }
         }
 
-        /// <summary>
-        /// Can compute the remainder of each element of vector using the operator %.
-        /// </summary>
         [Test]
-        public void CanComputeRemainderUsingOperator()
+        public void CanComputeModulus()
+        {
+            var vector = CreateVector(Data);
+            var mod = vector.Modulus(-3.2f);
+            for (var index = 0; index < Data.Length; index++)
+            {
+                AssertHelpers.AlmostEqualRelative(Euclid.Modulus(Data[index], -3.2f), mod[index], 14);
+            }
+        }
+
+        [Test]
+        public void CanComputeModulusUsingResultVector()
+        {
+            var vector = CreateVector(Data);
+            var mod = CreateVector(vector.Count);
+            vector.Modulus(-3.2f, mod);
+            for (var index = 0; index < Data.Length; index++)
+            {
+                AssertHelpers.AlmostEqualRelative(Euclid.Modulus(Data[index], -3.2f), mod[index], 14);
+            }
+        }
+
+        [Test]
+        public void CanComputeModulusUsingSameResultVector()
         {
             var vector = CreateVector(Data);
-            var mod = vector % 4.5f;
+            vector.Modulus(-3.2f, vector);
             for (var index = 0; index < Data.Length; index++)
             {
-                AssertHelpers.AlmostEqualRelative(Data[index] % 4.5f, mod[index], 14);
+                AssertHelpers.AlmostEqualRelative(Euclid.Modulus(Data[index], -3.2f), vector[index], 14);
             }
         }
     }