diff --git a/src/Numerics/Polynomial.cs b/src/Numerics/Polynomial.cs index 20084a28..ec604a14 100644 --- a/src/Numerics/Polynomial.cs +++ b/src/Numerics/Polynomial.cs @@ -234,9 +234,9 @@ namespace MathNet.Numerics // Negate, and normalize (scale such that the polynomial becomes monic) double aN = Coefficients[n]; double[] p = new double[n]; - for (int ii = n - 1; ii >= 0; ii--) + for (int i = n - 1; i >= 0; i--) { - p[ii] = -Coefficients[ii] / aN; + p[i] = -Coefficients[i] / aN; } DenseMatrix A0 = DenseMatrix.CreateDiagonal(n - 1, n - 1, 1.0); @@ -465,88 +465,71 @@ namespace MathNet.Numerics /// A tuple holding quotient in first and remainder in second public static Tuple DivideRemainder(Polynomial a, Polynomial b) { - if (b.Degree < 0) + var bDegree = b.Degree; + if (bDegree < 0) + { throw new DivideByZeroException("b polynomial ends with zero"); + } - if (a.Degree < 0) + var aDegree = a.Degree; + if (aDegree < 0) { // zero divided by non-zero is zero without remainder return Tuple.Create(a, a); } - var c1 = a.Coefficients.ToArray(); - var c2 = b.Coefficients.ToArray(); - - var n1 = c1.Length; - var n2 = c2.Length; - - double[] quo = null; - double[] rem = null; - - if (n2 == 1) // division by scalar + if (bDegree == 0) { - var fact = c2[0]; - quo = new double[n1]; - for (int i = 0; i < n1; i++) - quo[i] = c1[i] / fact; - rem = new double[] { 0 }; + // division by scalar + return Tuple.Create(Divide(a, b.Coefficients[0]), new Polynomial()); } - else if (n1 < n2) // denominator degree higher than nominator degree + + if (aDegree < bDegree) { + // denominator degree higher than nominator degree // quotient always be 0 and return c1 as remainder - quo = new double[] { 0 }; - rem = c1.ToArray(); + return Tuple.Create(new Polynomial(), a); } - else - { - var dn = n1 - n2; - var scl = c2[n2 - 1]; - var c22 = new double[n2 - 1]; - for (int ii = 0; ii < c22.Length; ii++) - { - c22[ii] = c2[ii] / scl; - } - - int i = dn; - int j = n1 - 1; - while (i >= 0) - { - var v = c1[j]; - for (int k = i; k < j; k++) - c1[k] -= c22[k - i] * v; - i--; - j--; - } - var j1 = j + 1; - var l1 = n1 - j1; - - rem = new double[j1]; - quo = new double[l1]; + var c1 = a.Coefficients.ToArray(); + var c2 = b.Coefficients.ToArray(); - for (int k = 0; k < l1; k++) - { - quo[k] = c1[k + j1] / scl; - } + var scl = c2[bDegree]; + var c22 = new double[bDegree]; + for (int ii = 0; ii < c22.Length; ii++) + { + c22[ii] = c2[ii] / scl; + } - for (int k = 0; k < j1; k++) + int i = aDegree - bDegree; + int j = aDegree; + while (i >= 0) + { + var v = c1[j]; + for (int k = i; k < j; k++) { - rem[k] = c1[k]; + c1[k] -= c22[k - i] * v; } - + i--; + j--; } - if (rem == null) - throw new NullReferenceException("resulting remainder was null"); + var j1 = j + 1; + var l1 = aDegree - j; - if (quo == null) - throw new NullReferenceException("resulting quotient was null"); + var quo = new double[l1]; + for (int k = 0; k < l1; k++) + { + quo[k] = c1[k + j1] / scl; + } - // output mapping - var pQuo = new Polynomial(quo); - var pRem = new Polynomial(rem); + var rem = new double[j1]; + for (int k = 0; k < j1; k++) + { + rem[k] = c1[k]; + } - return new Tuple(pQuo, pRem); + return Tuple.Create(new Polynomial(quo), new Polynomial(rem)); } #endregion @@ -954,7 +937,7 @@ namespace MathNet.Numerics } #endregion - #region + #region Clone public Polynomial Clone() { int degree = EvaluateDegree(Coefficients);