PageRenderTime 53ms CodeModel.GetById 27ms RepoModel.GetById 0ms app.codeStats 0ms

/Numerics/Numerics/Core/Polynomial.cs

#
C# | 285 lines | 150 code | 34 blank | 101 comment | 58 complexity | 771ad0ea8e2b3d80c481aa9e75aea86a MD5 | raw file
  1. using System;
  2. using System.Collections.Generic;
  3. using System.Text;
  4. namespace Meta.Numerics {
  5. /// <summary>
  6. /// Represents a polynomial with real coefficients.
  7. /// </summary>
  8. public class Polynomial {
  9. /// <summary>
  10. /// Initializes a new polynomial with the given coefficients.
  11. /// </summary>
  12. /// <param name="coefficients">The coefficients of the polynomial.</param>
  13. /// <returns>The specified polynomial.</returns>
  14. /// <remarks>
  15. /// <para>Coefficients should be arranged from low to high order, so that the kth entry is the coefficient of x<sup>k</sup>. For example,
  16. /// to specify the polynomial 5 - 6 x + 7 x<sup>2</sup>, give the values 5, -6, 7.</para>
  17. /// </remarks>
  18. public static Polynomial FromCoefficients (params double[] coefficients) {
  19. if (coefficients == null) throw new ArgumentNullException("coefficients");
  20. if (coefficients.Length == 0) throw new InvalidOperationException();
  21. return (new Polynomial(coefficients));
  22. }
  23. /// <summary>
  24. /// Initializes a new polynomial that passes through the given points.
  25. /// </summary>
  26. /// <param name="points">An N X 2 array whose first column contains the x values of points and whose second column contains the corresponding y values.</param>
  27. /// <returns>A polynomial of degree N-1 that passes through all the given points.</returns>
  28. public static Polynomial FromPoints (double[,] points) {
  29. if (points == null) throw new ArgumentNullException("points");
  30. if (points.GetLength(0) == 0) throw new ArgumentException("The first dimension of the points array must have length at least one.", "points");
  31. if (points.GetLength(1) != 2) throw new ArgumentException("The second dimension of the points array must have length two.", "points");
  32. double[] x = new double[points.GetLength(0)];
  33. double[] y = new double[points.GetLength(0)];
  34. for (int i = 0; i < points.GetLength(0); i++) {
  35. x[i] = points[i, 0];
  36. y[i] = points[i, 1];
  37. }
  38. PolynomialInterpolator interpolator = new PolynomialInterpolator(x, y);
  39. return (new InterpolatingPolynomial(interpolator));
  40. }
  41. /// <summary>
  42. /// Initializes a new polynomial that passes through the given points.
  43. /// </summary>
  44. /// <param name="points">A collection of points.</param>
  45. /// <returns>A polynomial that passes through all the given points.</returns>
  46. public static Polynomial FromPoints (ICollection<XY> points) {
  47. if (points == null) throw new ArgumentNullException("points");
  48. if (points.Count == 0) throw new ArgumentException("There must be at least one point in the points collection.", "points");
  49. double[] x = new double[points.Count];
  50. double[] y = new double[points.Count];
  51. int i = 0;
  52. foreach (XY point in points) {
  53. x[i] = point.X;
  54. y[i] = point.Y;
  55. i++;
  56. }
  57. PolynomialInterpolator interpolator = new PolynomialInterpolator(x, y);
  58. return (new InterpolatingPolynomial(interpolator));
  59. }
  60. internal Polynomial (double[] coefficients) {
  61. this.coefficients = coefficients;
  62. // set order, accounting for any trailing zeros in the array of coefficients
  63. order = coefficients.Length - 1;
  64. while ((order > 0) && (coefficients[order] == 0.0)) order--;
  65. }
  66. private double[] coefficients;
  67. private int order;
  68. /// <summary>
  69. /// Gets the degree of the polynomial.
  70. /// </summary>
  71. /// <remarks>
  72. /// <para>The degree of a polynomial is the highest power of the variable that appears. For example, the degree of 5 + 6 x + 7 x<sup>2</sup> is 2.</para>
  73. /// </remarks>
  74. public virtual int Degree {
  75. get {
  76. return (order);
  77. }
  78. }
  79. /// <summary>
  80. /// Gets the specificed coefficient.
  81. /// </summary>
  82. /// <param name="n">The power of the variable for which the coefficient is desired.</param>
  83. /// <returns>The coefficient of x<sup>n</sup>.</returns>
  84. /// <exception cref="ArgumentOutOfRangeException"><paramref name="n"/> is negative.</exception>
  85. public virtual double Coefficient (int n) {
  86. if (n < 0) {
  87. throw new ArgumentOutOfRangeException("n");
  88. } else if (n >= coefficients.Length) {
  89. return (0.0);
  90. } else {
  91. return (coefficients[n]);
  92. }
  93. }
  94. /// <summary>
  95. /// Evaluates the polynomial for the given input value.
  96. /// </summary>
  97. /// <param name="x">The value of the variable.</param>
  98. /// <returns>The value of the polynomial.</returns>
  99. public virtual double Evaluate (double x) {
  100. // This is horner's scheme, which is well-known but non-obvious, at least to me.
  101. // a_0 + a_1 x + a_2 x^2 + a_3 x^3 = ((a_3 x + a_2) x + a_1) x + a_0
  102. // It requires 2d flops, i.e. only one multiply and one add for each degree added, so it is more efficient than naive evaluation.
  103. // It is also ostensibly less subject to rounding errors than naive evaluation.
  104. double t = 0.0;
  105. for (int i = coefficients.Length - 1; i >= 0; i--) {
  106. t = t * x + coefficients[i];
  107. }
  108. return (t);
  109. }
  110. /// <summary>
  111. /// Differentiates the polynomial.
  112. /// </summary>
  113. /// <returns>The derivative of the polynomail.</returns>
  114. public virtual Polynomial Differentiate () {
  115. double[] coefficients = new double[this.Degree];
  116. for (int i = 0; i < coefficients.Length; i++) coefficients[i] = (i + 1) * this.Coefficient(i + 1);
  117. return (new Polynomial(coefficients));
  118. }
  119. /// <summary>
  120. /// Integrates the polynomail.
  121. /// </summary>
  122. /// <param name="C">The integration constant.</param>
  123. /// <returns>The integral of the polynomial.</returns>
  124. public virtual Polynomial Integrate (double C) {
  125. double[] coefficients = new double[this.Degree + 2];
  126. coefficients[0] = C;
  127. for (int i = 1; i < coefficients.Length; i++) coefficients[i] = this.Coefficient(i - 1) / i;
  128. return (new Polynomial(coefficients));
  129. }
  130. /// <summary>
  131. /// Generates a string representation of the polynomial.
  132. /// </summary>
  133. /// <returns>A string representation of the polynomial.</returns>
  134. public override string ToString () {
  135. StringBuilder text = new StringBuilder(Coefficient(0).ToString());
  136. for (int i = 1; i <= this.Degree; i++) {
  137. double coefficient = Coefficient(i);
  138. if (coefficient < 0.0) {
  139. text.AppendFormat(" - {0}", -coefficient);
  140. } else {
  141. text.AppendFormat(" + {0}", coefficient);
  142. }
  143. if (i == 1) {
  144. text.Append(" x");
  145. } else {
  146. text.AppendFormat(" x^{0}", i);
  147. }
  148. }
  149. return text.ToString();
  150. }
  151. /// <summary>
  152. /// Negates a polynomial.
  153. /// </summary>
  154. /// <param name="p">The polynomial.</param>
  155. /// <returns>The addative inverse of the polynomial.</returns>
  156. /// <exception cref="ArgumentNullException"><paramref name="p"/> is null.</exception>
  157. public static Polynomial operator - (Polynomial p) {
  158. if (p == null) throw new ArgumentNullException("p");
  159. double[] coefficients = new double[p.Degree + 1];
  160. for (int i = 0; i < coefficients.Length; i++) {
  161. coefficients[i] = -p.Coefficient(i);
  162. }
  163. return (new Polynomial(coefficients));
  164. }
  165. /// <summary>
  166. /// Computes the sum of two polynomials.
  167. /// </summary>
  168. /// <param name="p1">The first polynomial.</param>
  169. /// <param name="p2">The second polynomial.</param>
  170. /// <returns>The sum polynomial.</returns>
  171. /// <exception cref="ArgumentNullException"><paramref name="p1"/> or <paramref name="p2"/> is null.</exception>
  172. public static Polynomial operator + (Polynomial p1, Polynomial p2) {
  173. if (p1 == null) throw new ArgumentNullException("p1");
  174. if (p2 == null) throw new ArgumentNullException("p2");
  175. double[] coefficients = new double[Math.Max(p1.Degree, p2.Degree) + 1];
  176. for (int i = 0; i < coefficients.Length; i++) coefficients[i] = p1.Coefficient(i) + p2.Coefficient(i);
  177. return (new Polynomial(coefficients));
  178. }
  179. /// <summary>
  180. /// Computes the difference of two polynomials.
  181. /// </summary>
  182. /// <param name="p1">The first polynomial.</param>
  183. /// <param name="p2">The second polynomial.</param>
  184. /// <returns>The difference polynomial.</returns>
  185. /// <exception cref="ArgumentNullException"><paramref name="p1"/> or <paramref name="p2"/> is null.</exception>
  186. public static Polynomial operator - (Polynomial p1, Polynomial p2) {
  187. if (p1 == null) throw new ArgumentNullException("p1");
  188. if (p2 == null) throw new ArgumentNullException("p2");
  189. double[] coefficients = new double[Math.Max(p1.Degree, p2.Degree) + 1];
  190. for (int i = 0; i < coefficients.Length; i++) coefficients[i] = p1.Coefficient(i) - p2.Coefficient(i);
  191. return (new Polynomial(coefficients));
  192. }
  193. /// <summary>
  194. /// Computes the product of two polynomials.
  195. /// </summary>
  196. /// <param name="p1">The first polynomial.</param>
  197. /// <param name="p2">The second polynomial.</param>
  198. /// <returns>The product polynomial.</returns>
  199. /// <exception cref="ArgumentNullException"><paramref name="p1"/> or <paramref name="p2"/> is null.</exception>
  200. public static Polynomial operator * (Polynomial p1, Polynomial p2) {
  201. if (p1 == null) throw new ArgumentNullException("p1");
  202. if (p2 == null) throw new ArgumentNullException("p2");
  203. double[] coefficients = new double[p1.Degree + p2.Degree + 1];
  204. for (int i1 = 0; i1 <= p1.Degree; i1++) {
  205. for (int i2 = 0; i2 <= p2.Degree; i2++) {
  206. coefficients[i1 + i2] += p1.Coefficient(i1) * p2.Coefficient(i2);
  207. }
  208. }
  209. return (new Polynomial(coefficients));
  210. }
  211. /// <summary>
  212. /// Computes the quotient of two polynomials.
  213. /// </summary>
  214. /// <param name="p1">The dividend polynomial.</param>
  215. /// <param name="p2">The divisor polynomial.</param>
  216. /// <param name="remainder">The remainder polynomial.</param>
  217. /// <returns>The quotient polynomial.</returns>
  218. /// <remarks>
  219. /// <para>p<sub>1</sub> = q p<sub>2</sub> + r</para>
  220. /// </remarks>
  221. public static Polynomial Divide (Polynomial p1, Polynomial p2, out Polynomial remainder) {
  222. if (p1 == null) throw new ArgumentNullException("p1");
  223. if (p2 == null) throw new ArgumentNullException("p2");
  224. if (p2.Degree >= p1.Degree) {
  225. throw new InvalidOperationException();
  226. }
  227. // compute and store some value we will use repeatedly
  228. int d1 = p1.Degree;
  229. int d2 = p2.Degree;
  230. double a = p2.Coefficient(d2);
  231. // a is the leading coefficient of p2; it is the only number we ever divide by
  232. // (so if p2 is monic, polynomial division does not involve division at all!)
  233. // copy the coefficients of p1 into an working array; this will become our remainder when we are done
  234. double[] r = new double[d1 + 1];
  235. for (int i = 0; i < r.Length; i++) r[i] = p1.Coefficient(i);
  236. // create space for the coefficients of the quotient polynomial, which has order p1.Order - p2.Order
  237. double[] q = new double[d1 - d2 + 1];
  238. // do long division
  239. for (int k = q.Length - 1; k >= 0; k--) {
  240. q[k] = r[d2 + k] / a;
  241. r[d2 + k] = 0.0;
  242. for (int j = d2 + k - 1; j >= k; j--) {
  243. r[j] = r[j] - q[k] * p2.Coefficient(j - k);
  244. }
  245. }
  246. // form the remainder and quotient polynomials from the arrays
  247. remainder = new Polynomial(r);
  248. return (new Polynomial(q));
  249. }
  250. }
  251. }