001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2006, GeoTools Project Managment Committee (PMC)
005: * (C) 2001, Institut de Recherche pour le Développement
006: *
007: * This library is free software; you can redistribute it and/or
008: * modify it under the terms of the GNU Lesser General Public
009: * License as published by the Free Software Foundation;
010: * version 2.1 of the License.
011: *
012: * This library is distributed in the hope that it will be useful,
013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
015: * Lesser General Public License for more details.
016: */
017: package org.geotools.math;
018:
019: // J2SE dependencies
020: import java.io.Serializable;
021: import java.util.Arrays;
022:
023: import org.geotools.resources.Utilities;
024: import org.geotools.resources.XMath;
025:
026: /**
027: * The coefficients of a polynomial equation. The equation must be in the form
028: *
029: * <code>y = c<sub>0</sub> +
030: * c<sub>1</sub>×<var>x</var> +
031: * c<sub>2</sub>×<var>x</var><sup>2</sup> +
032: * c<sub>3</sub>×<var>x</var><sup>3</sup> + ... +
033: * c<sub>n</sub>×<var>x</var><sup>n</sup></code>.
034: *
035: * The static method {@link #roots(double[])} can be used for computing the root of a polynomial
036: * equation without creating a {@code Polygon} object.
037: *
038: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/metadata/src/main/java/org/geotools/math/Polynom.java $
039: * @version $Id: Polynom.java 22443 2006-10-27 20:47:22Z desruisseaux $
040: * @author Ken Turkiwski
041: * @author Martin Desruisseaux
042: *
043: * @since 2.0
044: */
045: public class Polynom implements Serializable {
046: /**
047: * Serial version UID for compatibility with different versions.
048: */
049: private static final long serialVersionUID = 6825019711186108990L;
050:
051: /**
052: * The array when no real roots can be computed.
053: */
054: private static final double[] NO_REAL_ROOT = new double[0];
055:
056: /**
057: * The coefficients for this polynom.
058: */
059: private final double[] c;
060:
061: /**
062: * The roots of this polynom. Will be computed only when first requested.
063: */
064: private transient double[] roots;
065:
066: /**
067: * Construct a polynom with the specified coefficients.
068: *
069: * @param c The coefficients. This array is copied.
070: */
071: public Polynom(final double[] c) {
072: int n = c.length;
073: while (n != 0 && c[--n] == 0)
074: ;
075: if (n == 0) {
076: this .c = NO_REAL_ROOT;
077: } else {
078: this .c = new double[n];
079: System.arraycopy(c, 0, this .c, 0, n);
080: }
081: }
082:
083: /**
084: * Evaluate this polynomial equation for the specified <var>x</var> value.
085: * More specifically, this method compute
086: * <code>c<sub>0</sub> +
087: * c<sub>1</sub>×<var>x</var> +
088: * c<sub>2</sub>×<var>x</var><sup>2</sup> +
089: * c<sub>3</sub>×<var>x</var><sup>3</sup> + ... +
090: * c<sub>n</sub>×<var>x</var><sup>n</sup></code>.
091: */
092: public final double y(final double x) {
093: double sum = 0;
094: for (int i = c.length; --i >= 0;) {
095: sum = sum * x + c[i];
096: }
097: return sum;
098: }
099:
100: /**
101: * Find the roots of a quadratic equation.
102: * More specifically, this method solves the following equation:
103: *
104: * <blockquote><code>
105: * c0 +
106: * c1*<var>x</var> +
107: * c2*<var>x</var><sup>2</sup> == 0
108: * </code></blockquote>
109: *
110: * @return The roots. The length may be 1 or 2.
111: */
112: private static double[] quadraticRoots(double c0, double c1,
113: double c2) {
114: double d = c1 * c1 - 4 * c2 * c0;
115: if (d > 0) {
116: // Two real, distinct roots
117: d = Math.sqrt(d);
118: if (c1 < 0) {
119: d = -d;
120: }
121: final double q = 0.5 * (d - c1);
122: return new double[] { q / c2,
123: (q != 0) ? c0 / q : -0.5 * (d + c1) / c2 };
124: } else if (d == 0) {
125: // One real double root
126: return new double[] { -c1 / (2 * c2) };
127: } else {
128: // Two complex conjugate roots
129: return NO_REAL_ROOT;
130: }
131: }
132:
133: /**
134: * Find the roots of a cubic equation.
135: * More specifically, this method solves the following equation:
136: *
137: * <blockquote><code>
138: * c0 +
139: * c1*<var>x</var> +
140: * c2*<var>x</var><sup>2</sup> +
141: * c3*<var>x</var><sup>3</sup> == 0
142: * </code></blockquote>
143: *
144: * @return The roots. The length may be 1 or 3.
145: */
146: private static double[] cubicRoots(double c0, double c1, double c2,
147: double c3) {
148: c2 /= c3;
149: c1 /= c3;
150: c0 /= c3;
151: final double Q = (c2 * c2 - 3 * c1) / 9;
152: final double R = (2 * c2 * c2 * c2 - 9 * c2 * c1 + 27 * c0) / 54;
153: final double Qcubed = Q * Q * Q;
154: final double d = Qcubed - R * R;
155:
156: c2 /= 3;
157: if (d >= 0) {
158: final double theta = Math.acos(R / Math.sqrt(Qcubed)) / 3;
159: final double scale = -2 * Math.sqrt(Q);
160: final double[] roots = new double[] {
161: scale * Math.cos(theta) - c2,
162: scale * Math.cos(theta + Math.PI * 2 / 3) - c2,
163: scale * Math.cos(theta + Math.PI * 4 / 3) - c2 };
164: assert Math.abs(roots[0] * roots[1] * roots[2] + c0) < 1E-6;
165: assert Math.abs(roots[0] + roots[1] + roots[2] + c2 * 3) < 1E-6;
166: assert Math.abs(roots[0] * roots[1] + roots[0] * roots[2]
167: + roots[1] * roots[2] - c1) < 1E-6;
168: return roots;
169: } else {
170: double e = XMath.cbrt(Math.sqrt(-d) + Math.abs(R));
171: if (R > 0) {
172: e = -e;
173: }
174: return new double[] { (e + Q / e) - c2 };
175: }
176: }
177:
178: /**
179: * Find the roots of this polynome.
180: *
181: * @return The roots.
182: */
183: public double[] roots() {
184: if (roots == null) {
185: roots = roots(c);
186: }
187: return (double[]) roots.clone();
188: }
189:
190: /**
191: * Find the roots of a polynomial equation. More specifically,
192: * this method solve the following equation:
193: *
194: * <blockquote><code>
195: * c[0] +
196: * c[1]*<var>x</var> +
197: * c[2]*<var>x</var><sup>2</sup> +
198: * c[3]*<var>x</var><sup>3</sup> +
199: * ... +
200: * c[n]*<var>x</var><sup>n</sup> == 0
201: * </code></blockquote>
202: *
203: * where <var>n</var> is the array length minus 1.
204: *
205: * @param c The coefficients for the polynomial equation.
206: * @return The roots. This array may have any length up to {@code n-1}.
207: * @throws UnsupportedOperationException if there is more coefficients than this method
208: * can handle.
209: */
210: public static double[] roots(final double[] c) {
211: int n = c.length;
212: while (n != 0 && c[--n] == 0)
213: ;
214: switch (n) {
215: case 0:
216: return NO_REAL_ROOT;
217: case 1:
218: return new double[] { -c[0] / c[1] };
219: case 2:
220: return quadraticRoots(c[0], c[1], c[2]);
221: case 3:
222: return cubicRoots(c[0], c[1], c[2], c[3]);
223: default:
224: throw new UnsupportedOperationException(String.valueOf(n));
225: }
226: }
227:
228: /**
229: * Display to the standard output the roots of a polynomial equation.
230: * More specifically, this method solve the following equation:
231: *
232: * <blockquote><code>
233: * c[0] +
234: * c[1]*<var>x</var> +
235: * c[2]*<var>x</var><sup>2</sup> +
236: * c[3]*<var>x</var><sup>3</sup> +
237: * ... +
238: * c[n]*<var>x</var><sup>n</sup> == 0
239: * </code></blockquote>
240: *
241: * where <var>n</var> is the array length minus 1.
242: *
243: * @param c The coefficients for the polynomial equation.
244: */
245: public static void main(final String[] c) {
246: final double[] r = new double[c.length];
247: for (int i = 0; i < c.length; i++) {
248: r[i] = Double.parseDouble(c[i]);
249: }
250: final double[] roots = roots(r);
251: for (int i = 0; i < roots.length; i++) {
252: System.out.println(roots[i]);
253: }
254: }
255:
256: /**
257: * Returns a hash value for this polynom.
258: */
259: public int hashCode() {
260: long code = c.length;
261: for (int i = c.length; --i >= 0;) {
262: code = code * 37 + Double.doubleToLongBits(c[i]);
263: }
264: return (int) code ^ (int) (code >>> 32);
265: }
266:
267: /**
268: * Compare this polynom with the specified object for equality.
269: */
270: public boolean equals(final Object object) {
271: if (object != null && object.getClass().equals(getClass())) {
272: final Polynom that = (Polynom) object;
273: return Arrays.equals(this .c, that.c);
274: }
275: return false;
276: }
277:
278: /**
279: * Returns a string representation of this polynom.
280: */
281: public String toString() {
282: final StringBuffer buffer = new StringBuffer(Utilities
283: .getShortClassName(this ));
284: buffer.append('[');
285: for (int i = 0; i < c.length; i++) {
286: if (i != 0) {
287: buffer.append(", ");
288: }
289: buffer.append(c[i]);
290: }
291: buffer.append(']');
292: return buffer.toString();
293: }
294: }
|