1 | /*
|
---|
2 | * Class Polynomial
|
---|
3 | *
|
---|
4 | * Defines a polynomial
|
---|
5 | * y = a[0] + a[1].x + a[2].x^2 + a[3].3 + . . . + a[n].x^n
|
---|
6 | * where x and all a[i] are real
|
---|
7 | * and deg is the degree of the polynomial, i.e. n,
|
---|
8 | * and includes the methods associated with polynomials,
|
---|
9 | *
|
---|
10 | * See clas ComplexPoly for Complex Polynomials
|
---|
11 | *
|
---|
12 | * WRITTEN BY: Dr Michael Thomas Flanagan
|
---|
13 | *
|
---|
14 | * See class Complex for standard complex arithmetic
|
---|
15 | *
|
---|
16 | * DATE: 6 June 2010 (adapted from ComplexPoly [February 2002])
|
---|
17 | * UPDATED: 21 January 2011
|
---|
18 | *
|
---|
19 | * DOCUMENTATION:
|
---|
20 | * See Michael Thomas Flanagan's Java library on-line web pages:
|
---|
21 | * http://www.ee.ucl.ac.uk/~mflanaga/java/Polynomial.html
|
---|
22 | * http://www.ee.ucl.ac.uk/~mflanaga/java/
|
---|
23 | *
|
---|
24 | *
|
---|
25 | * Copyright (c) 2002 - 2011
|
---|
26 | *
|
---|
27 | * PERMISSION TO COPY:
|
---|
28 | * Permission to use, copy and modify this software and its documentation for
|
---|
29 | * NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement
|
---|
30 | * to the author, Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies.
|
---|
31 | *
|
---|
32 | * Dr Michael Thomas Flanagan makes no representations about the suitability
|
---|
33 | * or fitness of the software for any or for a particular purpose.
|
---|
34 | * Michael Thomas Flanagan shall not be liable for any damages suffered
|
---|
35 | * as a result of using, modifying or distributing this software or its derivatives.
|
---|
36 | *
|
---|
37 | ***************************************************************************************/
|
---|
38 |
|
---|
39 | package agents.anac.y2015.agentBuyogV2.flanagan.math;
|
---|
40 |
|
---|
41 | import java.util.ArrayList;
|
---|
42 |
|
---|
43 | import agents.anac.y2015.agentBuyogV2.flanagan.analysis.Stat;
|
---|
44 | import agents.anac.y2015.agentBuyogV2.flanagan.complex.Complex;
|
---|
45 | import agents.anac.y2015.agentBuyogV2.flanagan.complex.ComplexPoly;
|
---|
46 | import agents.anac.y2015.agentBuyogV2.flanagan.io.FileOutput;
|
---|
47 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
|
---|
48 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
|
---|
49 |
|
---|
50 | import java.math.*;
|
---|
51 |
|
---|
52 |
|
---|
53 | public class Polynomial{
|
---|
54 |
|
---|
55 | private int deg = 0; // Degree of the polynomial
|
---|
56 | private int degwz = 0; // Degree of the polynomial with zero roots removed
|
---|
57 | private double[] coeff; // Array of polynomial coefficients
|
---|
58 | private double[] coeffwz; // Array of polynomial coefficients with zero roots removed
|
---|
59 |
|
---|
60 | private boolean suppressRootsErrorMessages = false; // = true if suppression of 'null returned' error messages in roots is required
|
---|
61 |
|
---|
62 | // CONSTRUCTORS
|
---|
63 | public Polynomial(int n){
|
---|
64 | this.deg = n;
|
---|
65 | this.coeff = new double[n+1];
|
---|
66 | }
|
---|
67 |
|
---|
68 | // Coefficients are double
|
---|
69 | public Polynomial(double[] aa){
|
---|
70 | this.deg =aa.length-1;
|
---|
71 | this.coeff = new double[this.deg+1];
|
---|
72 | for(int i=0; i<=this.deg; i++){
|
---|
73 | this.coeff[i]=aa[i];
|
---|
74 | }
|
---|
75 | }
|
---|
76 |
|
---|
77 | // Coefficients are real (float)
|
---|
78 | public Polynomial(float[] aa){
|
---|
79 | this.deg =aa.length-1;
|
---|
80 | coeff = new double[this.deg+1];
|
---|
81 | for(int i=0; i<=deg; i++){
|
---|
82 | this.coeff[i] = (double)aa[i];
|
---|
83 | }
|
---|
84 | }
|
---|
85 |
|
---|
86 | // Coefficients are long
|
---|
87 | public Polynomial(long[] aa){
|
---|
88 | this.deg =aa.length-1;
|
---|
89 | coeff = new double[this.deg+1];
|
---|
90 | for(int i=0; i<=deg; i++){
|
---|
91 | this.coeff[i] = (double)aa[i];
|
---|
92 | }
|
---|
93 | }
|
---|
94 |
|
---|
95 | // Coefficients are int
|
---|
96 | public Polynomial(int[] aa){
|
---|
97 | this.deg =aa.length-1;
|
---|
98 | coeff = new double[this.deg+1];
|
---|
99 | for(int i=0; i<=deg; i++){
|
---|
100 | this.coeff[i] = (double)aa[i];
|
---|
101 | }
|
---|
102 | }
|
---|
103 |
|
---|
104 | // Coefficients are in an a ArrayList - each element may be any relevant type -
|
---|
105 | // The ArrayList must contain the individual coefficients as irs elements and not as an array
|
---|
106 | public Polynomial(ArrayList<Object> aa){
|
---|
107 | this.deg = aa.size()-1;
|
---|
108 | coeff = new double[this.deg+1];
|
---|
109 | for(int i=0; i<=deg; i++){
|
---|
110 | int code = this.getTypeCode((Object)aa.get(i));
|
---|
111 | switch(code){
|
---|
112 | case 1: // Byte
|
---|
113 | this.coeff[i] = (double)((Byte)aa.get(i));
|
---|
114 | break;
|
---|
115 | case 2: // Short
|
---|
116 | this.coeff[i] = (double)((Short)aa.get(i));
|
---|
117 | break;
|
---|
118 | case 3: // Integer
|
---|
119 | this.coeff[i] = (double)((Integer)aa.get(i));
|
---|
120 | break;
|
---|
121 | case 4: // Long
|
---|
122 | this.coeff[i] = (double)((Long)aa.get(i));
|
---|
123 | break;
|
---|
124 | case 5: // Float
|
---|
125 | this.coeff[i] = (double)((Float)aa.get(i));
|
---|
126 | break;
|
---|
127 | case 6: // Double
|
---|
128 | this.coeff[i] = (double)((Double)aa.get(i));
|
---|
129 | break;
|
---|
130 | case 7: // BigInteger
|
---|
131 | this.coeff[i] = (double)((BigInteger)aa.get(i)).doubleValue();
|
---|
132 | break;
|
---|
133 | case 8:// BigDecimal
|
---|
134 | this.coeff[i] = (double)((BigDecimal)aa.get(i)).doubleValue();
|
---|
135 | break;
|
---|
136 | default: throw new IllegalArgumentException("Type code, " + code + ", not recognised");
|
---|
137 | }
|
---|
138 | }
|
---|
139 | }
|
---|
140 |
|
---|
141 | // Returns a code indicating the type of oObject passed as the argument
|
---|
142 | // Called by ArrayList constructor
|
---|
143 | private int getTypeCode(Object obj){
|
---|
144 |
|
---|
145 | int code = 0;
|
---|
146 |
|
---|
147 | if(obj instanceof Byte){
|
---|
148 | code = 1;
|
---|
149 | }
|
---|
150 | else{
|
---|
151 | if(obj instanceof Short){
|
---|
152 | code = 2;
|
---|
153 | }
|
---|
154 | else{
|
---|
155 | if(obj instanceof Integer){
|
---|
156 | code = 3;
|
---|
157 | }
|
---|
158 | else{
|
---|
159 | if(obj instanceof Long){
|
---|
160 | code = 4;
|
---|
161 | }
|
---|
162 | else{
|
---|
163 | if(obj instanceof Float){
|
---|
164 | code = 5;
|
---|
165 | }
|
---|
166 | else{
|
---|
167 | if(obj instanceof Double){
|
---|
168 | code = 6;
|
---|
169 | }
|
---|
170 | else{
|
---|
171 | if(obj instanceof BigInteger){
|
---|
172 | code = 7;
|
---|
173 | }
|
---|
174 | else{
|
---|
175 | if(obj instanceof BigDecimal){
|
---|
176 | code = 8;
|
---|
177 | }
|
---|
178 | }
|
---|
179 | }
|
---|
180 | }
|
---|
181 | }
|
---|
182 | }
|
---|
183 | }
|
---|
184 | }
|
---|
185 |
|
---|
186 | return code;
|
---|
187 | }
|
---|
188 |
|
---|
189 |
|
---|
190 | // Single constant - double
|
---|
191 | // y = aa
|
---|
192 | // needed in class Loop
|
---|
193 | public Polynomial(double aa){
|
---|
194 | this.deg = 0;
|
---|
195 | coeff = new double[1];
|
---|
196 | this.coeff[0] = aa;
|
---|
197 | }
|
---|
198 |
|
---|
199 | // Straight line - coefficients are double
|
---|
200 | // y = aa + bb.x
|
---|
201 | public Polynomial(double aa, double bb){
|
---|
202 | this.deg = 1;
|
---|
203 | coeff = new double[2];
|
---|
204 | this.coeff[0] = aa;
|
---|
205 | this.coeff[1] = bb;
|
---|
206 | }
|
---|
207 |
|
---|
208 | // Quadratic - coefficients are double
|
---|
209 | // y = aa + bb.x + cc.x^2
|
---|
210 | public Polynomial(double aa, double bb, double cc){
|
---|
211 | this.deg = 2;
|
---|
212 | coeff = new double[3];
|
---|
213 | this.coeff[0] = aa;
|
---|
214 | this.coeff[1] = bb;
|
---|
215 | this.coeff[2] = cc;
|
---|
216 | }
|
---|
217 |
|
---|
218 |
|
---|
219 | // Cubic - coefficients are double
|
---|
220 | // y = aa + bb.x + cc.x^2 + dd.x^3
|
---|
221 | public Polynomial(double aa, double bb, double cc, double dd){
|
---|
222 | this.deg = 3;
|
---|
223 | coeff = new double[4];
|
---|
224 | this.coeff[0] = aa;
|
---|
225 | this.coeff[1] = bb;
|
---|
226 | this.coeff[2] = cc;
|
---|
227 | this.coeff[3] = dd;
|
---|
228 | }
|
---|
229 |
|
---|
230 |
|
---|
231 | // METHODS
|
---|
232 |
|
---|
233 | // Returns a Polynomial given the polynomial's roots
|
---|
234 | public static Polynomial rootsToPoly(double[] roots){
|
---|
235 | if(roots==null)return null;
|
---|
236 |
|
---|
237 | int pdeg = roots.length;
|
---|
238 |
|
---|
239 | double[] rootCoeff = new double[2];
|
---|
240 | rootCoeff[0] = -roots[0];
|
---|
241 | rootCoeff[1] = 1.0;
|
---|
242 | Polynomial rPoly = new Polynomial(rootCoeff);
|
---|
243 | for(int i=1; i<pdeg; i++){
|
---|
244 | rootCoeff[0] = -roots[i];
|
---|
245 | Polynomial cRoot = new Polynomial(rootCoeff);
|
---|
246 | rPoly = rPoly.times(cRoot);
|
---|
247 | }
|
---|
248 | return rPoly;
|
---|
249 | }
|
---|
250 |
|
---|
251 | // Reset the polynomial with double[]
|
---|
252 | public void resetPoly(double [] aa){
|
---|
253 | if((this.deg+1)!=aa.length)throw new IllegalArgumentException("array lengths do not match");
|
---|
254 | for(int i=0; i<this.deg; i++){
|
---|
255 | this.coeff[i] = aa[i];
|
---|
256 | }
|
---|
257 | }
|
---|
258 |
|
---|
259 | // Reset the polynomial with ArrayList
|
---|
260 | public void resetPoly(ArrayList<Object> aa){
|
---|
261 | if((this.deg+1)!=aa.size())throw new IllegalArgumentException("array lengths do not match");
|
---|
262 | for(int i=0; i<=deg; i++){
|
---|
263 | int code = this.getTypeCode((Object)aa.get(i));
|
---|
264 | switch(code){
|
---|
265 | case 1: // Byte
|
---|
266 | this.coeff[i] = (double)((Byte)aa.get(i));
|
---|
267 | break;
|
---|
268 | case 2: // Short
|
---|
269 | this.coeff[i] = (double)((Short)aa.get(i));
|
---|
270 | break;
|
---|
271 | case 3: // Integer
|
---|
272 | this.coeff[i] = (double)((Integer)aa.get(i));
|
---|
273 | break;
|
---|
274 | case 4: // Long
|
---|
275 | this.coeff[i] = (double)((Long)aa.get(i));
|
---|
276 | break;
|
---|
277 | case 5: // Float
|
---|
278 | this.coeff[i] = (double)((Float)aa.get(i));
|
---|
279 | break;
|
---|
280 | case 6: // Double
|
---|
281 | this.coeff[i] = (double)((Double)aa.get(i));
|
---|
282 | break;
|
---|
283 | case 7: // BigInteger
|
---|
284 | this.coeff[i] = ((BigInteger)aa.get(i)).doubleValue();
|
---|
285 | break;
|
---|
286 | case 8: // BigDecimal
|
---|
287 | this.coeff[i] = ((BigDecimal)aa.get(i)).doubleValue();
|
---|
288 | break;
|
---|
289 | default: throw new IllegalArgumentException("Type code, " + code + ", not recognised");
|
---|
290 | }
|
---|
291 | }
|
---|
292 | }
|
---|
293 |
|
---|
294 |
|
---|
295 | // Reset a coefficient
|
---|
296 | public void resetCoeff(int i, double aa){
|
---|
297 | this.coeff[i] = aa;
|
---|
298 | }
|
---|
299 |
|
---|
300 | // Return a copy of this Polynomial [instance method]
|
---|
301 | public Polynomial copy(){
|
---|
302 | if(this==null){
|
---|
303 | return null;
|
---|
304 | }
|
---|
305 | else{
|
---|
306 | Polynomial aa = new Polynomial(this.deg);
|
---|
307 | for(int i=0; i<=this.deg; i++){
|
---|
308 | aa.coeff[i] = this.coeff[i];
|
---|
309 | }
|
---|
310 | aa.deg = this.deg;
|
---|
311 | aa.degwz = this.degwz;
|
---|
312 | aa.coeffwz = Conv.copy(this.coeffwz);
|
---|
313 | return aa;
|
---|
314 | }
|
---|
315 | }
|
---|
316 |
|
---|
317 | // Return a copy of this Polynomial [static]
|
---|
318 | public static Polynomial copy(Polynomial bb){
|
---|
319 | if(bb==null){
|
---|
320 | return null;
|
---|
321 | }
|
---|
322 | else{
|
---|
323 | Polynomial aa = new Polynomial(bb.deg);
|
---|
324 | for(int i=0; i<=bb.deg; i++){
|
---|
325 | aa.coeff[i] = bb.coeff[i];
|
---|
326 | }
|
---|
327 | aa.deg = bb.deg;
|
---|
328 | aa.degwz = bb.degwz;
|
---|
329 | aa.coeffwz = Conv.copy(bb.coeffwz);
|
---|
330 | return aa;
|
---|
331 | }
|
---|
332 | }
|
---|
333 |
|
---|
334 | // Clone a Polynomial
|
---|
335 | public Object clone(){
|
---|
336 | return (Object) this.copy();
|
---|
337 | }
|
---|
338 |
|
---|
339 | // Return a copy of the polynomial coefficients
|
---|
340 | public double[] coefficientsCopy(){
|
---|
341 | double[] aa = new double[this.deg+1];
|
---|
342 | for(int i=0; i<=this.deg; i++){
|
---|
343 | aa[i] = this.coeff[i];
|
---|
344 | }
|
---|
345 | return aa;
|
---|
346 | }
|
---|
347 |
|
---|
348 | // Return a reference to the polynomial coefficients array
|
---|
349 | public double[] coefficientsReference(){
|
---|
350 | return this.coeff;
|
---|
351 | }
|
---|
352 |
|
---|
353 | // Return a coefficient
|
---|
354 | public double getCoefficient(int i){
|
---|
355 | return this.coeff[i];
|
---|
356 | }
|
---|
357 |
|
---|
358 | // Return the degree
|
---|
359 | public int getDeg(){
|
---|
360 | return this.deg;
|
---|
361 | }
|
---|
362 |
|
---|
363 |
|
---|
364 | // Convert to a String of the form a[0] + a[1].x + a[2].x^2 etc.
|
---|
365 | public String toString(){
|
---|
366 | String ss = "";
|
---|
367 | ss = ss + this.coeff[0];
|
---|
368 | if(this.deg>0)ss = ss + " + (" + this.coeff[1] + ").x";
|
---|
369 | for(int i=2; i<=this.deg; i++){
|
---|
370 | ss = ss + " + (" + this.coeff[i] + ").x^" + i;
|
---|
371 | }
|
---|
372 | return ss;
|
---|
373 | }
|
---|
374 |
|
---|
375 | // Convert to a ComplexPoly
|
---|
376 | public ComplexPoly toComplexPoly(){
|
---|
377 | ComplexPoly cp = new ComplexPoly(this);
|
---|
378 | return cp;
|
---|
379 | }
|
---|
380 |
|
---|
381 |
|
---|
382 | // Print the polynomial to screen
|
---|
383 | public void print(){
|
---|
384 | System.out.print(this.toString());
|
---|
385 | }
|
---|
386 |
|
---|
387 | // Print the polynomial to screen with line return
|
---|
388 | public void println(){
|
---|
389 | System.out.println(this.toString());
|
---|
390 | }
|
---|
391 |
|
---|
392 | // Print the polynomial to a text file with title
|
---|
393 | public void printToText(String title){
|
---|
394 | title = title + ".txt";
|
---|
395 | FileOutput fout = new FileOutput(title, 'n');
|
---|
396 |
|
---|
397 | fout.println("Output File for a Polynomial");
|
---|
398 | fout.dateAndTimeln();
|
---|
399 | fout.println();
|
---|
400 | fout.print("Polynomial degree is ");
|
---|
401 | fout.println(this.deg);
|
---|
402 | fout.println();
|
---|
403 | fout.println("The coefficients are ");
|
---|
404 |
|
---|
405 | for(int i=0;i<=this.deg;i++){
|
---|
406 | fout.println(this.coeff[i]);
|
---|
407 | }
|
---|
408 | fout.println();
|
---|
409 | fout.println("End of file.");
|
---|
410 | fout.close();
|
---|
411 | }
|
---|
412 |
|
---|
413 | // Print the polynomial to a text file without a given title
|
---|
414 | public void printToText(){
|
---|
415 | String title = "PolynomialOut";
|
---|
416 | printToText(title);
|
---|
417 | }
|
---|
418 |
|
---|
419 | // TRANSFORM
|
---|
420 | // Transform a real polynomial to an s-Domain polynomial ratio
|
---|
421 | // Returns ArrayList:
|
---|
422 | // First element: numerator as ComplexPoly
|
---|
423 | // Second element: denominator as ComplexPoly
|
---|
424 | // Instance method
|
---|
425 | public ArrayList<ComplexPoly> sTransform(){
|
---|
426 | return Polynomial.sTransform(this.coeff);
|
---|
427 | }
|
---|
428 |
|
---|
429 | // TRANSFORM
|
---|
430 | // Transform a real polynomial to an s-Domain polynomial ratio
|
---|
431 | // Argument: Polynomial a0 + a1.x + a2.x^2 + ....
|
---|
432 | // Returns ArrayList:
|
---|
433 | // First element: numerator as ComplexPoly
|
---|
434 | // Second element: denominator as ComplexPoly
|
---|
435 | // Static method
|
---|
436 | public static ArrayList<ComplexPoly> sTransform(Polynomial poly){
|
---|
437 | return Polynomial.sTransform(poly.coefficientsCopy());
|
---|
438 | }
|
---|
439 |
|
---|
440 |
|
---|
441 | // TRANSFORM
|
---|
442 | // Transform a real polynomial to an s-Domain polynomial ratio
|
---|
443 | // Argument: coefficients of the real polynomial a0 + a1.x + a2.x^2 + ....
|
---|
444 | // Returns ArrayList:
|
---|
445 | // First element: numerator as ComplexPoly
|
---|
446 | // Second element: denominator as ComplexPoly
|
---|
447 | // Static method
|
---|
448 | public static ArrayList<ComplexPoly> sTransform(double[] coeff){
|
---|
449 | int n = coeff.length;
|
---|
450 | ComplexPoly[] sNum = new ComplexPoly[n]; // numerator of each transformed term
|
---|
451 | ComplexPoly[] sDen = new ComplexPoly[n]; // denomenator of each transformed term
|
---|
452 | ComplexPoly sNumer = null; // numerator of the completely transformed polynomial
|
---|
453 | ComplexPoly sDenom = new ComplexPoly(Complex.plusOne()); // denomenator of the completely transformed polynomial
|
---|
454 |
|
---|
455 | // s-Transform of each term of the polynomial
|
---|
456 | for(int i=0; i<n; i++){
|
---|
457 | sNum[i] = new ComplexPoly(new Complex(coeff[i]*Fmath.factorial(i), 0));
|
---|
458 | sDen[i] = new ComplexPoly(i+1);
|
---|
459 | sDen[i].resetCoeff(i+1, Complex.plusOne());
|
---|
460 | }
|
---|
461 |
|
---|
462 | // create a common denomenator
|
---|
463 | sDenom = sDen[n-1];
|
---|
464 |
|
---|
465 | // create a common numerator
|
---|
466 | for(int i=0; i<n-1; i++){
|
---|
467 | sNum[i] = sNum[i].times(sDen[n-i-2]);
|
---|
468 | }
|
---|
469 | sNumer = sNum[0];
|
---|
470 | for(int i=1; i<n; i++)sNumer = sNumer.plus(sNum[i]);
|
---|
471 |
|
---|
472 | // Output arrayList
|
---|
473 | ArrayList<ComplexPoly> al = new ArrayList<ComplexPoly>();
|
---|
474 | al.add(sNumer);
|
---|
475 | al.add(sDenom);
|
---|
476 |
|
---|
477 | return al;
|
---|
478 | }
|
---|
479 |
|
---|
480 |
|
---|
481 |
|
---|
482 | // LOGICAL TESTS
|
---|
483 | // Check if two polynomials are identical
|
---|
484 | public boolean equals(Polynomial cp){
|
---|
485 | return isEqual(cp);
|
---|
486 | }
|
---|
487 |
|
---|
488 | public boolean isEqual(Polynomial cp){
|
---|
489 | boolean ret = false;
|
---|
490 | int nDegThis = this.getDeg();
|
---|
491 | int nDegCp = cp.getDeg();
|
---|
492 | if(nDegThis==nDegCp){
|
---|
493 | boolean test = true;
|
---|
494 | int i=0;
|
---|
495 | while(test){
|
---|
496 | if(this.coeff[i]!=cp.getCoefficient(i)){
|
---|
497 | test = false;
|
---|
498 | }
|
---|
499 | else{
|
---|
500 | i++;
|
---|
501 | if(i>nDegCp){
|
---|
502 | test = false;
|
---|
503 | ret = true;
|
---|
504 | }
|
---|
505 | }
|
---|
506 | }
|
---|
507 | }
|
---|
508 | return ret;
|
---|
509 | }
|
---|
510 |
|
---|
511 | // Check if two polynomials are identical (static)
|
---|
512 | public static boolean isEqual(Polynomial cp1, Polynomial cp2){
|
---|
513 | boolean ret = false;
|
---|
514 | int nDegCp1 = cp1.getDeg();
|
---|
515 | int nDegCp2 = cp2.getDeg();
|
---|
516 | if(nDegCp1==nDegCp2){
|
---|
517 | boolean test = true;
|
---|
518 | int i=0;
|
---|
519 | while(test){
|
---|
520 | if(cp1.getCoefficient(i)!=cp2.getCoefficient(i)){
|
---|
521 | test = false;
|
---|
522 | }
|
---|
523 | else{
|
---|
524 | i++;
|
---|
525 | if(i>nDegCp1){
|
---|
526 | test = false;
|
---|
527 | ret = true;
|
---|
528 | }
|
---|
529 | }
|
---|
530 | }
|
---|
531 | }
|
---|
532 | return ret;
|
---|
533 | }
|
---|
534 |
|
---|
535 | // ADDITION OF TWO POLYNOMIALS
|
---|
536 | // Addition, instance method
|
---|
537 | public Polynomial plus(Polynomial b){
|
---|
538 |
|
---|
539 | Polynomial c = null;
|
---|
540 | if(b.deg<=this.deg){
|
---|
541 | c = new Polynomial(this.deg);
|
---|
542 | for(int i=b.deg+1; i<=this.deg; i++)c.coeff[i] = this.coeff[i];
|
---|
543 | for(int i=0; i<=b.deg; i++)c.coeff[i] = this.coeff[i] + b.coeff[i];
|
---|
544 | }
|
---|
545 | else{
|
---|
546 | c = new Polynomial(b.deg);
|
---|
547 | for(int i=this.deg+1; i<=b.deg; i++)c.coeff[i] = b.coeff[i];
|
---|
548 | for(int i=0; i<=this.deg; i++)c.coeff[i] = this.coeff[i] + b.coeff[i];
|
---|
549 | }
|
---|
550 | return c;
|
---|
551 | }
|
---|
552 |
|
---|
553 | // Addition, static method
|
---|
554 | public static Polynomial plus(Polynomial a, Polynomial b){
|
---|
555 | Polynomial c = null;
|
---|
556 | if(b.deg<=a.deg){
|
---|
557 | c = new Polynomial(a.deg);
|
---|
558 | for(int i=b.deg+1; i<=a.deg; i++)c.coeff[i] = a.coeff[i];
|
---|
559 | for(int i=0; i<=b.deg; i++)c.coeff[i] = a.coeff[i] + b.coeff[i];
|
---|
560 | }
|
---|
561 | else{
|
---|
562 | c = new Polynomial(b.deg);
|
---|
563 | for(int i=a.deg+1; i<=b.deg; i++)c.coeff[i] = b.coeff[i];
|
---|
564 | for(int i=0; i<=a.deg; i++)c.coeff[i] = a.coeff[i] + b.coeff[i];
|
---|
565 | }
|
---|
566 | return c;
|
---|
567 | }
|
---|
568 |
|
---|
569 | // Addition of a double, instance method
|
---|
570 | public Polynomial plus(double bb){
|
---|
571 | Polynomial b = new Polynomial(bb);
|
---|
572 | return this.plus(b);
|
---|
573 | }
|
---|
574 |
|
---|
575 | // Addition of a double, static method
|
---|
576 | public static Polynomial plus(Polynomial a, double bb){
|
---|
577 | Polynomial b = new Polynomial(bb);
|
---|
578 | return Polynomial.plus(a, b);
|
---|
579 | }
|
---|
580 |
|
---|
581 | // Addition of an int, instance method
|
---|
582 | public Polynomial plus(int bb){
|
---|
583 | Polynomial b = new Polynomial((double) bb);
|
---|
584 | return this.plus(b);
|
---|
585 | }
|
---|
586 |
|
---|
587 | // Addition of an int, static method
|
---|
588 | public static Polynomial plus(Polynomial a, int bb){
|
---|
589 | Polynomial b = new Polynomial((double)bb);
|
---|
590 | return Polynomial.plus(a, b);
|
---|
591 | }
|
---|
592 |
|
---|
593 | // SUBTRACTION OF TWO POLYNOMIALS
|
---|
594 | // Subtraction, instance method
|
---|
595 | public Polynomial minus(Polynomial b){
|
---|
596 |
|
---|
597 | Polynomial c = null;
|
---|
598 | if(b.deg<=this.deg){
|
---|
599 | c = new Polynomial(this.deg);
|
---|
600 | for(int i=b.deg+1; i<=this.deg; i++)c.coeff[i] = this.coeff[i];
|
---|
601 | for(int i=0; i<=b.deg; i++)c.coeff[i] = this.coeff[i] - b.coeff[i];
|
---|
602 | }
|
---|
603 | else{
|
---|
604 | c = new Polynomial(b.deg);
|
---|
605 | for(int i=this.deg+1; i<=b.deg; i++)c.coeff[i] = -b.coeff[i];
|
---|
606 | for(int i=0; i<=this.deg; i++)c.coeff[i] = this.coeff[i] - b.coeff[i];
|
---|
607 | }
|
---|
608 | return c;
|
---|
609 | }
|
---|
610 |
|
---|
611 | // Subtraction of a double, instance method
|
---|
612 | public Polynomial minus(double bb){
|
---|
613 | Polynomial b = new Polynomial(bb);
|
---|
614 | return this.minus(b);
|
---|
615 | }
|
---|
616 |
|
---|
617 | // Subtraction of an int, instance method
|
---|
618 | public Polynomial minus(int bb){
|
---|
619 | Polynomial b = new Polynomial((double)bb);
|
---|
620 | return this.minus(b);
|
---|
621 | }
|
---|
622 |
|
---|
623 | // Subtraction, static method
|
---|
624 | public static Polynomial minus(Polynomial a, Polynomial b){
|
---|
625 | Polynomial c = null;
|
---|
626 | if(b.deg<=a.deg){
|
---|
627 | c = new Polynomial(a.deg);
|
---|
628 | for(int i=b.deg+1; i<=a.deg; i++)c.coeff[i] = a.coeff[i];
|
---|
629 | for(int i=0; i<=b.deg; i++)c.coeff[i] = a.coeff[i] - b.coeff[i];
|
---|
630 | }
|
---|
631 | else{
|
---|
632 | c = new Polynomial(b.deg);
|
---|
633 | for(int i=a.deg+1; i<=b.deg; i++)c.coeff[i] = -b.coeff[i];
|
---|
634 | for(int i=0; i<=a.deg; i++)c.coeff[i] = a.coeff[i] - b.coeff[i];
|
---|
635 | }
|
---|
636 | return c;
|
---|
637 | }
|
---|
638 |
|
---|
639 | // Subtraction of a double, static method
|
---|
640 | public static Polynomial minus(Polynomial a, double bb){
|
---|
641 | Polynomial b = new Polynomial(bb);
|
---|
642 | return Polynomial.minus(a, b);
|
---|
643 | }
|
---|
644 |
|
---|
645 | // Subtraction of a int, static method
|
---|
646 | public static Polynomial minus(Polynomial a, int bb){
|
---|
647 | Polynomial b = new Polynomial((double)bb);
|
---|
648 | return Polynomial.minus(a, b);
|
---|
649 | }
|
---|
650 |
|
---|
651 |
|
---|
652 | // MULTIPLICATION OF TWO POLYNOMIALS
|
---|
653 | // Multiplication, instance method
|
---|
654 | public Polynomial times(Polynomial b){
|
---|
655 | int n = this.deg + b.deg;
|
---|
656 | Polynomial c = new Polynomial(n);
|
---|
657 | for(int i=0; i<=this.deg; i++){
|
---|
658 | for(int j=0; j<=b.deg; j++){
|
---|
659 | c.coeff[i+j] += (this.coeff[i]*b.coeff[j]);
|
---|
660 | }
|
---|
661 | }
|
---|
662 | return c;
|
---|
663 | }
|
---|
664 |
|
---|
665 | // Multiplication, static method
|
---|
666 | public static Polynomial times(Polynomial a, Polynomial b){
|
---|
667 | int n = a.deg + b.deg;
|
---|
668 | Polynomial c = new Polynomial(n);
|
---|
669 | for(int i=0; i<=a.deg; i++){
|
---|
670 | for(int j=0; j<=b.deg; j++){
|
---|
671 | c.coeff[i+j] += (a.coeff[i]*b.coeff[j]);
|
---|
672 | }
|
---|
673 | }
|
---|
674 | return c;
|
---|
675 | }
|
---|
676 |
|
---|
677 | // Multiplication by a double, instance method
|
---|
678 | public Polynomial times(double bb){
|
---|
679 | Polynomial c = new Polynomial(this.deg);
|
---|
680 | for(int i=0; i<=this.deg; i++){
|
---|
681 | c.coeff[i] = this.coeff[i]*bb;
|
---|
682 | }
|
---|
683 | return c;
|
---|
684 | }
|
---|
685 |
|
---|
686 | // Multiplication by a double, static method
|
---|
687 | public static Polynomial times(Polynomial a, double bb){
|
---|
688 | Polynomial c = new Polynomial(a.deg);
|
---|
689 | for(int i=0; i<=a.deg; i++){
|
---|
690 | c.coeff[i] = a.coeff[i]*bb;
|
---|
691 | }
|
---|
692 | return c;
|
---|
693 | }
|
---|
694 |
|
---|
695 | // Multiplication by an int, instance method
|
---|
696 | public Polynomial times(int bb){
|
---|
697 | Polynomial c = new Polynomial(this.deg);
|
---|
698 | for(int i=0; i<=this.deg; i++){
|
---|
699 | c.coeff[i] = this.coeff[i]*bb;
|
---|
700 | }
|
---|
701 | return c;
|
---|
702 | }
|
---|
703 |
|
---|
704 | // Multiplication by an int, static method
|
---|
705 | public static Polynomial times(Polynomial a, int bb){
|
---|
706 | Polynomial c = new Polynomial(a.deg);
|
---|
707 | for(int i=0; i<=a.deg; i++){
|
---|
708 | c.coeff[i] = a.coeff[i]*bb;
|
---|
709 | }
|
---|
710 | return c;
|
---|
711 | }
|
---|
712 |
|
---|
713 |
|
---|
714 | // DERIVATIVES
|
---|
715 | // Return the coefficients, as a new Polynomial, of the nth derivative
|
---|
716 | public Polynomial nthDerivative(int n){
|
---|
717 | Polynomial dnydxn;
|
---|
718 |
|
---|
719 | if(n>this.deg){
|
---|
720 | dnydxn = new Polynomial(0.0);
|
---|
721 | }
|
---|
722 | else{
|
---|
723 | dnydxn = new Polynomial(this.deg-n);
|
---|
724 | double[] nc = new double[this.deg - n + 1];
|
---|
725 |
|
---|
726 | int k = this.deg - n;
|
---|
727 | for(int i=this.deg; i>n-1; i--){
|
---|
728 | nc[k] = this.coeff[i];
|
---|
729 | for(int j=0; j<n; j++){
|
---|
730 | nc[k] = nc[k]*(i-j);
|
---|
731 | }
|
---|
732 | k--;
|
---|
733 | }
|
---|
734 | dnydxn = new Polynomial(nc);
|
---|
735 | }
|
---|
736 | return dnydxn;
|
---|
737 | }
|
---|
738 |
|
---|
739 | // EVALUATION OF A POLYNOMIAL AND ITS DERIVATIVES
|
---|
740 | // Evaluate the polynomial
|
---|
741 | public double evaluate(double x){
|
---|
742 | double y = 0.0;
|
---|
743 | if(this.deg==0){
|
---|
744 | y = this.coeff[0];
|
---|
745 | }
|
---|
746 | else{
|
---|
747 | y = this.coeff[this.deg];
|
---|
748 | for(int i=this.deg-1; i>=0; i--){
|
---|
749 | y = y*x + this.coeff[i];
|
---|
750 | }
|
---|
751 | }
|
---|
752 | return y;
|
---|
753 | }
|
---|
754 |
|
---|
755 | // Evaluate the nth derivative of the polynomial
|
---|
756 | public double nthDerivEvaluate(int n, double x){
|
---|
757 | double dnydxn = 0.0;
|
---|
758 | double[] nc = new double[this.deg+1];
|
---|
759 |
|
---|
760 | if(n==0)
|
---|
761 | {
|
---|
762 | dnydxn=this.evaluate(x);
|
---|
763 | System.out.println("n = 0 in Polynomial.nthDerivative");
|
---|
764 | System.out.println("polynomial itself evaluated and returned");
|
---|
765 | }
|
---|
766 | else{
|
---|
767 | Polynomial nthderiv = this.nthDerivative(n);
|
---|
768 | dnydxn = nthderiv.evaluate(x);
|
---|
769 | }
|
---|
770 | return dnydxn;
|
---|
771 | }
|
---|
772 |
|
---|
773 |
|
---|
774 |
|
---|
775 | // ROOTS OF POLYNOMIALS
|
---|
776 | // For general details of root searching and a discussion of the rounding errors
|
---|
777 | // see Numerical Recipes, The Art of Scientific Computing
|
---|
778 | // by W H Press, S A Teukolsky, W T Vetterling & B P Flannery
|
---|
779 | // Cambridge University Press, http://www.nr.com/
|
---|
780 |
|
---|
781 | // Calculate the roots (real or double) of a polynomial (real or double)
|
---|
782 | // polish = true ([for deg>3 see laguerreAll(...)]
|
---|
783 | // initial root estimates are all zero [for deg>3 see laguerreAll(...)]
|
---|
784 | public Complex[] roots(){
|
---|
785 | ComplexPoly cp = new ComplexPoly(this);
|
---|
786 | return cp.roots();
|
---|
787 | }
|
---|
788 |
|
---|
789 | // Calculate the roots - as above with the exception that the error messages are suppressed
|
---|
790 | // Required by BlackBox
|
---|
791 | public Complex[] rootsNoMessages(){
|
---|
792 | ComplexPoly cp = new ComplexPoly(this);
|
---|
793 | return cp.rootsNoMessages();
|
---|
794 | }
|
---|
795 |
|
---|
796 | // Calculate the roots (real or double) of a polynomial (real or double)
|
---|
797 | // initial root estimates are all zero [for deg>3 see laguerreAll(...)]
|
---|
798 | // for polish see laguerreAll(...)[for deg>3]
|
---|
799 | public Complex[] roots(boolean polish){
|
---|
800 | ComplexPoly cp = new ComplexPoly(this);
|
---|
801 | return cp.roots(polish);
|
---|
802 | }
|
---|
803 |
|
---|
804 | // Calculate the roots (real or double) of a polynomial (real or double)
|
---|
805 | // for estx see laguerreAll(...)[for deg>3] - initial estimate of first root
|
---|
806 | // polish = true see laguerreAll(...)[for deg>3]
|
---|
807 | public Complex[] roots(double estx){
|
---|
808 | ComplexPoly cp = new ComplexPoly(this);
|
---|
809 | return cp.roots(new Complex(estx, 0.0));
|
---|
810 | }
|
---|
811 |
|
---|
812 | public Complex[] roots(Complex estx){
|
---|
813 | ComplexPoly cp = new ComplexPoly(this);
|
---|
814 | return cp.roots(estx);
|
---|
815 | }
|
---|
816 |
|
---|
817 | // Calculate the roots (real or complex) of a polynomial (real or complex)
|
---|
818 | public Complex[] roots(boolean polish, double estx){
|
---|
819 | ComplexPoly cp = new ComplexPoly(this);
|
---|
820 | return cp.roots(new Complex(estx, 0.0));
|
---|
821 | }
|
---|
822 |
|
---|
823 | // ROOTS OF A QUADRATIC EQUATION
|
---|
824 | // ax^2 + bx + c = 0
|
---|
825 | // roots returned in root[]
|
---|
826 | // 4ac << b*b accomodated by these methods
|
---|
827 | public static Complex[] quadratic(double c, double b, double a){
|
---|
828 | return ComplexPoly.quadratic(new Complex(c, 0.0), new Complex(b, 0.0), new Complex(a, 0.0));
|
---|
829 | }
|
---|
830 |
|
---|
831 | // ROOTS OF A CUBIC EQUATION
|
---|
832 | // ddx^3 + ccx^2 + bbx + aa = 0
|
---|
833 | // roots returned in root[]
|
---|
834 | public static Complex[] cubic(double aa, double bb, double cc, double dd){
|
---|
835 | return ComplexPoly.cubic(new Complex(aa, 0.0), new Complex(bb, 0.0), new Complex(cc, 0.0), new Complex(dd, 0.0));
|
---|
836 | }
|
---|
837 |
|
---|
838 |
|
---|
839 | // LAGUERRE'S METHOD FOR double ROOTS OF A double POLYNOMIAL
|
---|
840 |
|
---|
841 | // Laguerre method for one of the roots
|
---|
842 | // Following the procedure in Numerical Recipes for C [Reference above]
|
---|
843 | // estx estimate of the root
|
---|
844 | // coeff[] coefficients of the polynomial
|
---|
845 | // m degree of the polynomial
|
---|
846 | public static Complex laguerre(double estx, double[] pcoeff, int m){
|
---|
847 | ArrayMaths am = new ArrayMaths(pcoeff);
|
---|
848 | return ComplexPoly.laguerre(new Complex(estx, 0.0), am.array_as_Complex(), m);
|
---|
849 | }
|
---|
850 |
|
---|
851 | public static Complex laguerre(Complex estx, double[] pcoeff, int m){
|
---|
852 | ArrayMaths am = new ArrayMaths(pcoeff);
|
---|
853 | return ComplexPoly.laguerre(estx, am.array_as_Complex(), m);
|
---|
854 | }
|
---|
855 |
|
---|
856 | // Finds all roots of a double polynomial by successive calls to laguerre
|
---|
857 | // Following the procedure in Numerical Recipes for C [Reference above]
|
---|
858 | // Initial estimates are all zero, polish=true
|
---|
859 | public Complex[] laguerreAll(){
|
---|
860 | ComplexPoly cp = new ComplexPoly(this);
|
---|
861 | return cp.laguerreAll();
|
---|
862 | }
|
---|
863 |
|
---|
864 | // Initial estimates estx, polish=true
|
---|
865 | public Complex[] laguerreAll(double estx){
|
---|
866 | ComplexPoly cp = new ComplexPoly(this);
|
---|
867 | return cp.laguerreAll(new Complex(estx, 0.0));
|
---|
868 | }
|
---|
869 |
|
---|
870 | public Complex[] laguerreAll(Complex estx){
|
---|
871 | ComplexPoly cp = new ComplexPoly(this);
|
---|
872 | return cp.laguerreAll(estx);
|
---|
873 | }
|
---|
874 |
|
---|
875 | // Initial estimates are all zero.
|
---|
876 | public Complex[] laguerreAll(boolean polish){
|
---|
877 | ComplexPoly cp = new ComplexPoly(this);
|
---|
878 | return cp.laguerreAll(polish);
|
---|
879 | }
|
---|
880 |
|
---|
881 | // Finds all roots of a double polynomial by successive calls to laguerre
|
---|
882 | // Initial estimates are estx
|
---|
883 | public Complex[] laguerreAll(boolean polish, double estx){
|
---|
884 | ComplexPoly cp = new ComplexPoly(this);
|
---|
885 | return cp.laguerreAll(polish, new Complex(estx, 0.0));
|
---|
886 | }
|
---|
887 |
|
---|
888 | public Complex[] laguerreAll(boolean polish, Complex estx){
|
---|
889 | ComplexPoly cp = new ComplexPoly(this);
|
---|
890 | return cp.laguerreAll(polish, estx);
|
---|
891 | }
|
---|
892 | }
|
---|
893 |
|
---|