1 | /*
|
---|
2 | * Class ComplexPoly
|
---|
3 | *
|
---|
4 | * Defines a complex 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] may be real or complex
|
---|
7 | * and deg is the degree of the polynomial, i.e. n,
|
---|
8 | * and includes the methods associated with polynomials,
|
---|
9 | * e.g. complex root searches
|
---|
10 | *
|
---|
11 | * WRITTEN BY: Dr Michael Thomas Flanagan
|
---|
12 | *
|
---|
13 | * See class Complex for standard complex arithmetic
|
---|
14 | *
|
---|
15 | * DATE: February 2002
|
---|
16 | * UPDATED: 22 June 2003, 19 January 2005, 12 May 2005, 11 October 2005, 30 April 2007
|
---|
17 | * 22 November 2007, 10-11 August 2008, 22 August 2008, 31 October 2009, 6 November 2009
|
---|
18 | * 29 May 2010, 5-6 June 2010, 20 October 2010, 18-21 January 2011
|
---|
19 | *
|
---|
20 | * DOCUMENTATION:
|
---|
21 | * See Michael Thomas Flanagan's Java library on-line web pages:
|
---|
22 | * http://www.ee.ucl.ac.uk/~mflanaga/java/ComplexPoly.html
|
---|
23 | * http://www.ee.ucl.ac.uk/~mflanaga/java/
|
---|
24 | *
|
---|
25 | *
|
---|
26 | * Copyright (c) 2002 - 2009
|
---|
27 | *
|
---|
28 | * PERMISSION TO COPY:
|
---|
29 | * Permission to use, copy and modify this software and its documentation for
|
---|
30 | * NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement
|
---|
31 | * to the author, Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies.
|
---|
32 | *
|
---|
33 | * Dr Michael Thomas Flanagan makes no representations about the suitability
|
---|
34 | * or fitness of the software for any or for a particular purpose.
|
---|
35 | * Michael Thomas Flanagan shall not be liable for any damages suffered
|
---|
36 | * as a result of using, modifying or distributing this software or its derivatives.
|
---|
37 | *
|
---|
38 | ***************************************************************************************/
|
---|
39 |
|
---|
40 | package agents.anac.y2015.agentBuyogV2.flanagan.complex;
|
---|
41 |
|
---|
42 | import java.util.ArrayList;
|
---|
43 |
|
---|
44 | import agents.anac.y2015.agentBuyogV2.flanagan.circuits.Phasor;
|
---|
45 | import agents.anac.y2015.agentBuyogV2.flanagan.control.BlackBox;
|
---|
46 | import agents.anac.y2015.agentBuyogV2.flanagan.io.FileOutput;
|
---|
47 | import agents.anac.y2015.agentBuyogV2.flanagan.math.ArrayMaths;
|
---|
48 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
|
---|
49 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
|
---|
50 | import agents.anac.y2015.agentBuyogV2.flanagan.math.Polynomial;
|
---|
51 |
|
---|
52 | import java.math.*;
|
---|
53 |
|
---|
54 |
|
---|
55 | public class ComplexPoly{
|
---|
56 |
|
---|
57 | private int deg = 0; // Degree of the polynomial
|
---|
58 | private int degwz = 0; // Degree of the polynomial with zero roots removed
|
---|
59 | private Complex[] coeff; // Array of polynomial coefficients
|
---|
60 | private Complex[] coeffwz; // Array of polynomial coefficients with zero roots removed
|
---|
61 |
|
---|
62 | private boolean suppressRootsErrorMessages = false; // = true if suppression of 'null returned' error messages in roots is required
|
---|
63 |
|
---|
64 | // CONSTRUCTORS
|
---|
65 | public ComplexPoly(int n){
|
---|
66 | this.deg = n;
|
---|
67 | this.coeff = Complex.oneDarray(n+1);
|
---|
68 | }
|
---|
69 |
|
---|
70 | // Coefficients are complex
|
---|
71 | public ComplexPoly(Complex[] aa){
|
---|
72 | this.deg =aa.length-1;
|
---|
73 | this.coeff = Complex.oneDarray(this.deg+1);
|
---|
74 | for(int i=0; i<=this.deg; i++){
|
---|
75 | this.coeff[i]=Complex.copy(aa[i]);
|
---|
76 | }
|
---|
77 | }
|
---|
78 |
|
---|
79 | // Coefficients are real (double)
|
---|
80 | public ComplexPoly(double[] aa){
|
---|
81 | this.deg =aa.length-1;
|
---|
82 | coeff = Complex.oneDarray(this.deg+1);
|
---|
83 | for(int i=0; i<=deg; i++){
|
---|
84 | this.coeff[i].reset(aa[i], 0.0);
|
---|
85 | }
|
---|
86 | }
|
---|
87 |
|
---|
88 | // Coefficients are real (float)
|
---|
89 | public ComplexPoly(float[] aa){
|
---|
90 | this.deg =aa.length-1;
|
---|
91 | coeff = Complex.oneDarray(this.deg+1);
|
---|
92 | for(int i=0; i<=deg; i++){
|
---|
93 | this.coeff[i].reset((double)aa[i], 0.0);
|
---|
94 | }
|
---|
95 | }
|
---|
96 |
|
---|
97 | // Coefficients are int
|
---|
98 | public ComplexPoly(int[] aa){
|
---|
99 | this.deg =aa.length-1;
|
---|
100 | coeff = Complex.oneDarray(this.deg+1);
|
---|
101 | for(int i=0; i<=deg; i++){
|
---|
102 | this.coeff[i].reset((double)aa[i], 0.0);
|
---|
103 | }
|
---|
104 | }
|
---|
105 |
|
---|
106 | // Argument is a Polynolial (real polynomial)
|
---|
107 | public ComplexPoly(Polynomial aa){
|
---|
108 | this.deg = aa.getDeg();
|
---|
109 | coeff = Complex.oneDarray(this.deg+1);
|
---|
110 | for(int i=0; i<=deg; i++){
|
---|
111 | this.coeff[i].reset(aa.getCoefficient(i), 0.0);
|
---|
112 | }
|
---|
113 | }
|
---|
114 |
|
---|
115 | // Coefficients are in an a ArrayList - each element may be any relevant type -
|
---|
116 | // The ArrayList must contain the individual coefficients as irs elements and not as an array
|
---|
117 | public ComplexPoly(ArrayList<Object> aa){
|
---|
118 | this.deg =aa.size()-1;
|
---|
119 | coeff = Complex.oneDarray(this.deg+1);
|
---|
120 | for(int i=0; i<=deg; i++){
|
---|
121 | int code = this.getTypeCode((Object)aa.get(i));
|
---|
122 | switch(code){
|
---|
123 | case 1: // Byte
|
---|
124 | this.coeff[i].reset((double)((Byte)aa.get(i)), 0.0);
|
---|
125 | break;
|
---|
126 | case 2: // Short
|
---|
127 | this.coeff[i].reset((double)((Short)aa.get(i)), 0.0);
|
---|
128 | break;
|
---|
129 | case 3: // Integer
|
---|
130 | this.coeff[i].reset((double)((Integer)aa.get(i)), 0.0);
|
---|
131 | break;
|
---|
132 | case 4: // Long
|
---|
133 | this.coeff[i].reset((double)((Long)aa.get(i)), 0.0);
|
---|
134 | break;
|
---|
135 | case 5: // Float
|
---|
136 | this.coeff[i].reset((double)((Float)aa.get(i)), 0.0);
|
---|
137 | break;
|
---|
138 | case 6: // Double
|
---|
139 | this.coeff[i].reset((double)((Double)aa.get(i)), 0.0);
|
---|
140 | break;
|
---|
141 | case 7: // Complex
|
---|
142 | this.coeff[i] = (Complex)aa.get(i);
|
---|
143 | break;
|
---|
144 | case 8: // Phasor
|
---|
145 | this.coeff[i] = ((Phasor)aa.get(i)).toComplex();
|
---|
146 | break;
|
---|
147 | case 9: // BigInteger
|
---|
148 | this.coeff[i].reset(((BigInteger)aa.get(i)).doubleValue(), 0.0);
|
---|
149 | break;
|
---|
150 | case 10:// BigDecimal
|
---|
151 | this.coeff[i].reset(((BigDecimal)aa.get(i)).doubleValue(), 0.0);
|
---|
152 | break;
|
---|
153 | default: throw new IllegalArgumentException("Type code, " + code + ", not recognised");
|
---|
154 | }
|
---|
155 | }
|
---|
156 | }
|
---|
157 |
|
---|
158 | // Returns a code indicating the type of oObject passed as the argument
|
---|
159 | // Called by ArrayList constructor
|
---|
160 | private int getTypeCode(Object obj){
|
---|
161 |
|
---|
162 | int code = 0;
|
---|
163 |
|
---|
164 | if(obj instanceof Byte){
|
---|
165 | code = 1;
|
---|
166 | }
|
---|
167 | else{
|
---|
168 | if(obj instanceof Short){
|
---|
169 | code = 2;
|
---|
170 | }
|
---|
171 | else{
|
---|
172 | if(obj instanceof Integer){
|
---|
173 | code = 3;
|
---|
174 | }
|
---|
175 | else{
|
---|
176 | if(obj instanceof Long){
|
---|
177 | code = 4;
|
---|
178 | }
|
---|
179 | else{
|
---|
180 | if(obj instanceof Float){
|
---|
181 | code = 5;
|
---|
182 | }
|
---|
183 | else{
|
---|
184 | if(obj instanceof Double){
|
---|
185 | code = 6;
|
---|
186 | }
|
---|
187 | else{
|
---|
188 | if(obj instanceof Complex){
|
---|
189 | code = 7;
|
---|
190 | }
|
---|
191 | else{
|
---|
192 | if(obj instanceof Phasor){
|
---|
193 | code = 8;
|
---|
194 | }
|
---|
195 | else{
|
---|
196 | if(obj instanceof BigInteger){
|
---|
197 | code = 9;
|
---|
198 | }
|
---|
199 | else{
|
---|
200 | if(obj instanceof BigDecimal){
|
---|
201 | code = 10;
|
---|
202 | }
|
---|
203 | }
|
---|
204 | }
|
---|
205 | }
|
---|
206 | }
|
---|
207 | }
|
---|
208 | }
|
---|
209 | }
|
---|
210 | }
|
---|
211 | }
|
---|
212 |
|
---|
213 | return code;
|
---|
214 | }
|
---|
215 |
|
---|
216 |
|
---|
217 | // Single constant - complex
|
---|
218 | // y = aa
|
---|
219 | // needed in class Loop
|
---|
220 | public ComplexPoly(Complex aa){
|
---|
221 | this.deg = 0;
|
---|
222 | coeff = Complex.oneDarray(1);
|
---|
223 | this.coeff[0]=Complex.copy(aa);
|
---|
224 | }
|
---|
225 |
|
---|
226 | // Single constant - double
|
---|
227 | // y = aa
|
---|
228 | // needed in class Loop
|
---|
229 | public ComplexPoly(double aa){
|
---|
230 | this.deg = 0;
|
---|
231 | coeff = Complex.oneDarray(1);
|
---|
232 | this.coeff[0].reset(aa, 0.0);
|
---|
233 | }
|
---|
234 |
|
---|
235 | // Straight line - coefficients are complex
|
---|
236 | // y = aa + bb.x
|
---|
237 | public ComplexPoly(Complex aa, Complex bb){
|
---|
238 | this.deg = 1;
|
---|
239 | coeff = Complex.oneDarray(2);
|
---|
240 | this.coeff[0]=Complex.copy(aa);
|
---|
241 | this.coeff[1]=Complex.copy(bb);
|
---|
242 | }
|
---|
243 |
|
---|
244 | // Straight line - coefficients are real
|
---|
245 | // y = aa + bb.x
|
---|
246 | public ComplexPoly(double aa, double bb){
|
---|
247 | this.deg = 1;
|
---|
248 | coeff = Complex.oneDarray(2);
|
---|
249 | this.coeff[0].reset(aa, 0.0);
|
---|
250 | this.coeff[1].reset(bb, 0.0);
|
---|
251 | }
|
---|
252 |
|
---|
253 | // Quadratic - coefficients are complex
|
---|
254 | // y = aa + bb.x + cc.x^2
|
---|
255 | public ComplexPoly(Complex aa, Complex bb, Complex cc){
|
---|
256 | this.deg = 2;
|
---|
257 | coeff = Complex.oneDarray(3);
|
---|
258 | this.coeff[0]=Complex.copy(aa);
|
---|
259 | this.coeff[1]=Complex.copy(bb);
|
---|
260 | this.coeff[2]=Complex.copy(cc);
|
---|
261 | }
|
---|
262 |
|
---|
263 | // Quadratic - coefficients are real
|
---|
264 | // y = aa + bb.x + cc.x^2
|
---|
265 | public ComplexPoly(double aa, double bb, double cc){
|
---|
266 | this.deg = 2;
|
---|
267 | coeff = Complex.oneDarray(3);
|
---|
268 | this.coeff[0].reset(aa, 0.0);
|
---|
269 | this.coeff[1].reset(bb, 0.0);
|
---|
270 | this.coeff[2].reset(cc, 0.0);
|
---|
271 | }
|
---|
272 |
|
---|
273 | // Cubic - coefficients are complex
|
---|
274 | // y = aa + bb.x + cc.x^2 + dd.x^3
|
---|
275 | public ComplexPoly(Complex aa, Complex bb, Complex cc, Complex dd){
|
---|
276 | this.deg = 3;
|
---|
277 | coeff = Complex.oneDarray(4);
|
---|
278 | this.coeff[0]=Complex.copy(aa);
|
---|
279 | this.coeff[1]=Complex.copy(bb);
|
---|
280 | this.coeff[2]=Complex.copy(cc);
|
---|
281 | this.coeff[3]=Complex.copy(dd);
|
---|
282 | }
|
---|
283 |
|
---|
284 | // Cubic - coefficients are real
|
---|
285 | // y = aa + bb.x + cc.x^2 + dd.x^3
|
---|
286 | public ComplexPoly(double aa, double bb, double cc, double dd){
|
---|
287 | this.deg = 3;
|
---|
288 | coeff = Complex.oneDarray(4);
|
---|
289 | this.coeff[0].reset(aa, 0.0);
|
---|
290 | this.coeff[1].reset(bb, 0.0);
|
---|
291 | this.coeff[2].reset(cc, 0.0);
|
---|
292 | this.coeff[3].reset(dd, 0.0);
|
---|
293 | }
|
---|
294 |
|
---|
295 | // METHODS
|
---|
296 |
|
---|
297 | // Returns a ComplexPoly given the polynomial's roots
|
---|
298 | public static ComplexPoly rootsToPoly(Complex[] roots){
|
---|
299 | if(roots==null)return null;
|
---|
300 |
|
---|
301 | int pdeg = roots.length;
|
---|
302 |
|
---|
303 | Complex[] rootCoeff = Complex.oneDarray(2);
|
---|
304 | rootCoeff[0] = roots[0].times(Complex.minusOne());
|
---|
305 | rootCoeff[1] = Complex.plusOne();
|
---|
306 | ComplexPoly rPoly = new ComplexPoly(rootCoeff);
|
---|
307 | for(int i=1; i<pdeg; i++){
|
---|
308 | rootCoeff[0] = roots[i].times(Complex.minusOne());
|
---|
309 | ComplexPoly cRoot = new ComplexPoly(rootCoeff);
|
---|
310 | rPoly = rPoly.times(cRoot);
|
---|
311 | }
|
---|
312 | return rPoly;
|
---|
313 | }
|
---|
314 |
|
---|
315 | // Reset the polynomial with Complex[]
|
---|
316 | public void resetPoly(Complex [] aa){
|
---|
317 | if((this.deg+1)!=aa.length)throw new IllegalArgumentException("array lengths do not match");
|
---|
318 | for(int i=0; i<this.deg; i++){
|
---|
319 | this.coeff[i] = Complex.copy(aa[i]);
|
---|
320 | }
|
---|
321 | }
|
---|
322 |
|
---|
323 | // Reset the polynomial with ArrayList
|
---|
324 | public void resetPoly(ArrayList<Object> aa){
|
---|
325 | if((this.deg+1)!=aa.size())throw new IllegalArgumentException("array lengths do not match");
|
---|
326 | for(int i=0; i<=deg; i++){
|
---|
327 | int code = this.getTypeCode((Object)aa.get(i));
|
---|
328 | switch(code){
|
---|
329 | case 1: // Byte
|
---|
330 | this.coeff[i].reset((double)((Byte)aa.get(i)), 0.0);
|
---|
331 | break;
|
---|
332 | case 2: // Short
|
---|
333 | this.coeff[i].reset((double)((Short)aa.get(i)), 0.0);
|
---|
334 | break;
|
---|
335 | case 3: // Integer
|
---|
336 | this.coeff[i].reset((double)((Integer)aa.get(i)), 0.0);
|
---|
337 | break;
|
---|
338 | case 4: // Long
|
---|
339 | this.coeff[i].reset((double)((Long)aa.get(i)), 0.0);
|
---|
340 | break;
|
---|
341 | case 5: // Float
|
---|
342 | this.coeff[i].reset((double)((Float)aa.get(i)), 0.0);
|
---|
343 | break;
|
---|
344 | case 6: // Double
|
---|
345 | this.coeff[i].reset((double)((Double)aa.get(i)), 0.0);
|
---|
346 | break;
|
---|
347 | case 7: // Complex
|
---|
348 | this.coeff[i] = (Complex)aa.get(i);
|
---|
349 | break;
|
---|
350 | case 8: // Phasor
|
---|
351 | this.coeff[i] = ((Phasor)aa.get(i)).toComplex();
|
---|
352 | break;
|
---|
353 | case 9: // BigInteger
|
---|
354 | this.coeff[i].reset(((BigInteger)aa.get(i)).doubleValue(), 0.0);
|
---|
355 | break;
|
---|
356 | case 10:// BigDecimal
|
---|
357 | this.coeff[i].reset(((BigDecimal)aa.get(i)).doubleValue(), 0.0);
|
---|
358 | break;
|
---|
359 | default: throw new IllegalArgumentException("Type code, " + code + ", not recognised");
|
---|
360 | }
|
---|
361 | }
|
---|
362 | }
|
---|
363 |
|
---|
364 |
|
---|
365 | // Reset a coefficient
|
---|
366 | public void resetCoeff(int i, Complex aa){
|
---|
367 | this.coeff[i] = Complex.copy(aa);
|
---|
368 | }
|
---|
369 |
|
---|
370 | // Return a copy of this ComplexPoly [instance method]
|
---|
371 | public ComplexPoly copy(){
|
---|
372 | if(this==null){
|
---|
373 | return null;
|
---|
374 | }
|
---|
375 | else{
|
---|
376 | ComplexPoly aa = new ComplexPoly(this.deg);
|
---|
377 | for(int i=0; i<=this.deg; i++){
|
---|
378 | aa.coeff[i] = Complex.copy(this.coeff[i]);
|
---|
379 | }
|
---|
380 | aa.deg=this.deg;
|
---|
381 | aa.degwz = this.degwz;
|
---|
382 | aa.coeffwz = Conv.copy(this.coeffwz);
|
---|
383 |
|
---|
384 | return aa;
|
---|
385 | }
|
---|
386 | }
|
---|
387 |
|
---|
388 | // Return a copy of this ComplexPoly [static]
|
---|
389 | public static ComplexPoly copy(ComplexPoly bb){
|
---|
390 | if(bb==null){
|
---|
391 | return null;
|
---|
392 | }
|
---|
393 | else{
|
---|
394 | ComplexPoly aa = new ComplexPoly(bb.deg);
|
---|
395 | for(int i=0; i<=bb.deg; i++){
|
---|
396 | aa.coeff[i] = Complex.copy(bb.coeff[i]);
|
---|
397 | }
|
---|
398 | aa.deg=bb.deg;
|
---|
399 | aa.degwz = bb.degwz;
|
---|
400 | aa.coeffwz = Conv.copy(bb.coeffwz);
|
---|
401 |
|
---|
402 | return aa;
|
---|
403 | }
|
---|
404 | }
|
---|
405 |
|
---|
406 | // Clone a ComplexPoly
|
---|
407 | public Object clone(){
|
---|
408 | return (Object) this.copy();
|
---|
409 | }
|
---|
410 |
|
---|
411 | // Return a copy of the polynomial coefficients
|
---|
412 | public Complex[] polyNomCopy(){
|
---|
413 | Complex[] aa = Complex.oneDarray(this.deg+1);
|
---|
414 | for(int i=0; i<=this.deg; i++){
|
---|
415 | aa[i] = Complex.copy(this.coeff[i]);
|
---|
416 | }
|
---|
417 | return aa;
|
---|
418 | }
|
---|
419 |
|
---|
420 | // Return a reference to the polynomial
|
---|
421 | public Complex[] polyNomReference(){
|
---|
422 | return this.coeff;
|
---|
423 | }
|
---|
424 | // Return a reference to the polynomial
|
---|
425 | public Complex[] polyNomPointer(){
|
---|
426 | return this.coeff;
|
---|
427 | }
|
---|
428 |
|
---|
429 | // Return a copy of a coefficient
|
---|
430 | public Complex coeffCopy(int i){
|
---|
431 | return Complex.copy(this.coeff[i]);
|
---|
432 | }
|
---|
433 |
|
---|
434 | // Return a reference to a coefficient
|
---|
435 | public Complex coeffReference(int i){
|
---|
436 | return this.coeff[i];
|
---|
437 | }
|
---|
438 |
|
---|
439 | // Return a reference to a coefficient
|
---|
440 | public Complex coeffPointer(int i){
|
---|
441 | return this.coeff[i];
|
---|
442 | }
|
---|
443 |
|
---|
444 | // Return the degree
|
---|
445 | public int getDeg(){
|
---|
446 | return this.deg;
|
---|
447 | }
|
---|
448 |
|
---|
449 | // Sets the representation of the square root of minus one to j in Strings
|
---|
450 | public void setj(){
|
---|
451 | Complex.setj();
|
---|
452 | }
|
---|
453 |
|
---|
454 | // Sets the representation of the square root of minus one to i in Strings
|
---|
455 | public void seti(){
|
---|
456 | Complex.seti();
|
---|
457 | }
|
---|
458 |
|
---|
459 | // Convert to a String of the form (a+jb)[0] + (a+jb)[1].x + (a+jb)[2].x^2 etc.
|
---|
460 | public String toString(){
|
---|
461 | String ss = "";
|
---|
462 | ss = ss + this.coeffCopy(0).toString();
|
---|
463 | if(this.deg>0)ss = ss + " + (" + this.coeffCopy(1).toString() + ").x";
|
---|
464 | for(int i=2; i<=this.deg; i++){
|
---|
465 | ss = ss + " + (" + this.coeffCopy(i).toString() + ").x^" + i;
|
---|
466 | }
|
---|
467 | return ss;
|
---|
468 | }
|
---|
469 |
|
---|
470 | // Print the polynomial to screen
|
---|
471 | public void print(){
|
---|
472 | System.out.print(this.toString());
|
---|
473 | }
|
---|
474 |
|
---|
475 | // Print the polynomial to screen with line return
|
---|
476 | public void println(){
|
---|
477 | System.out.println(this.toString());
|
---|
478 | }
|
---|
479 |
|
---|
480 | // Print the polynomial to a text file with title
|
---|
481 | public void printToText(String title){
|
---|
482 | title = title + ".txt";
|
---|
483 | FileOutput fout = new FileOutput(title, 'n');
|
---|
484 |
|
---|
485 | fout.println("Output File for a ComplexPoly");
|
---|
486 | fout.dateAndTimeln();
|
---|
487 | fout.println();
|
---|
488 | fout.print("Polynomial degree is ");
|
---|
489 | fout.println(this.deg);
|
---|
490 | fout.println();
|
---|
491 | fout.println("The coefficients are ");
|
---|
492 |
|
---|
493 | for(int i=0;i<=this.deg;i++){
|
---|
494 | fout.println(this.coeff[i]);
|
---|
495 | }
|
---|
496 | fout.println();
|
---|
497 | fout.println("End of file.");
|
---|
498 | fout.close();
|
---|
499 | }
|
---|
500 |
|
---|
501 | // Print the polynomial to a text file without a given title
|
---|
502 | public void printToText(){
|
---|
503 | String title = "ComplexPolyOut";
|
---|
504 | printToText(title);
|
---|
505 | }
|
---|
506 |
|
---|
507 | // TRANSFORM
|
---|
508 | // Transform a Complex polynomial to an s-Domain polynomial ratio
|
---|
509 | // Argument: coefficients of the real polynomial a0 + a1.x + a2.x^2 + ....
|
---|
510 | // Returns ArrayList:
|
---|
511 | // First element: numertor as ComplexPoly
|
---|
512 | // Second element: denominator as ComplexPoly
|
---|
513 | public ArrayList<ComplexPoly> sTransform(){
|
---|
514 | return ComplexPoly.sTransform(this.coeff);
|
---|
515 | }
|
---|
516 |
|
---|
517 | // TRANSFORM
|
---|
518 | // Transform a Complex polynomial to an s-Domain polynomial ratio
|
---|
519 | // Argument: coefficients of the real polynomial a0 + a1.x + a2.x^2 + ....
|
---|
520 | // Returns ArrayList:
|
---|
521 | // First element: numertor as ComplexPoly
|
---|
522 | // Second element: denominator as ComplexPoly
|
---|
523 | public static ArrayList<ComplexPoly> sTransform(double[] coeff){
|
---|
524 | ArrayMaths am = new ArrayMaths(coeff);
|
---|
525 | Complex[] coeffc = am.getArray_as_Complex();
|
---|
526 | return ComplexPoly.sTransform(coeffc);
|
---|
527 | }
|
---|
528 |
|
---|
529 | // TRANSFORM
|
---|
530 | // Transform a real polynomial to an s-Domain polynomial ratio
|
---|
531 | // Argument: coefficients of the real polynomial a0 + a1.x + a2.x^2 + ....
|
---|
532 | // Returns ArrayList:
|
---|
533 | // First element: numertor as ComplexPoly
|
---|
534 | // Second element: denominator as ComplexPoly
|
---|
535 | public static ArrayList<ComplexPoly> sTransform(Complex[] coeff){
|
---|
536 | int n = coeff.length;
|
---|
537 | ComplexPoly[] sNum = new ComplexPoly[n]; // numerator of each transformed term
|
---|
538 | ComplexPoly[] sDen = new ComplexPoly[n]; // denomenator of each transformed term
|
---|
539 | ComplexPoly sNumer = null; // numerator of the completely transformed polynomial
|
---|
540 | ComplexPoly sDenom = new ComplexPoly(Complex.plusOne()); // denomenator of the completely transformed polynomial
|
---|
541 |
|
---|
542 | // s-Transform of each term of the polynomial
|
---|
543 | for(int i=0; i<n; i++){
|
---|
544 | sNum[i] = new ComplexPoly(coeff[i].times(new Complex(Fmath.factorial(i), 0)));
|
---|
545 | sDen[i] = new ComplexPoly(i+1);
|
---|
546 | sDen[i].resetCoeff(i+1, Complex.plusOne());
|
---|
547 | }
|
---|
548 |
|
---|
549 | // create a common denomenator
|
---|
550 | sDenom = sDen[n-1];
|
---|
551 |
|
---|
552 | // create a common numerator
|
---|
553 | for(int i=0; i<n-1; i++){
|
---|
554 | sNum[i] = sNum[i].times(sDen[n-i-2]);
|
---|
555 | }
|
---|
556 | sNumer = sNum[0];
|
---|
557 | for(int i=1; i<n; i++)sNumer = sNumer.plus(sNum[i]);
|
---|
558 |
|
---|
559 | // Output arrayList
|
---|
560 | ArrayList<ComplexPoly> al = new ArrayList<ComplexPoly>();
|
---|
561 | al.add(sNumer);
|
---|
562 | al.add(sDenom);
|
---|
563 |
|
---|
564 | return al;
|
---|
565 | }
|
---|
566 |
|
---|
567 | // LOGICAL TESTS
|
---|
568 | // Check if two polynomials are identical
|
---|
569 | public boolean equals(ComplexPoly cp){
|
---|
570 | return isEqual(cp);
|
---|
571 | }
|
---|
572 |
|
---|
573 | public boolean isEqual(ComplexPoly cp){
|
---|
574 | boolean ret = false;
|
---|
575 | int nDegThis = this.getDeg();
|
---|
576 | int nDegCp = cp.getDeg();
|
---|
577 | if(nDegThis==nDegCp){
|
---|
578 | boolean test = true;
|
---|
579 | int i=0;
|
---|
580 | while(test){
|
---|
581 | if(!this.coeff[i].isEqual(cp.coeffReference(i))){
|
---|
582 | test = false;
|
---|
583 | }
|
---|
584 | else{
|
---|
585 | i++;
|
---|
586 | if(i>nDegCp){
|
---|
587 | test = false;
|
---|
588 | ret = true;
|
---|
589 | }
|
---|
590 | }
|
---|
591 | }
|
---|
592 | }
|
---|
593 | return ret;
|
---|
594 | }
|
---|
595 |
|
---|
596 | // Check if two polynomials are identical (static)
|
---|
597 | public static boolean isEqual(ComplexPoly cp1, ComplexPoly cp2){
|
---|
598 | boolean ret = false;
|
---|
599 | int nDegCp1 = cp1.getDeg();
|
---|
600 | int nDegCp2 = cp2.getDeg();
|
---|
601 | if(nDegCp1==nDegCp2){
|
---|
602 | boolean test = true;
|
---|
603 | int i=0;
|
---|
604 | while(test){
|
---|
605 | if(!cp1.coeffReference(i).isEqual(cp2.coeffReference(i))){
|
---|
606 | test = false;
|
---|
607 | }
|
---|
608 | else{
|
---|
609 | i++;
|
---|
610 | if(i>nDegCp1){
|
---|
611 | test = false;
|
---|
612 | ret = true;
|
---|
613 | }
|
---|
614 | }
|
---|
615 | }
|
---|
616 | }
|
---|
617 | return ret;
|
---|
618 | }
|
---|
619 |
|
---|
620 | // ADDITION OF TWO POLYNOMIALS
|
---|
621 | // Addition, instance method
|
---|
622 | public ComplexPoly plus(ComplexPoly b){
|
---|
623 |
|
---|
624 | ComplexPoly c = null;
|
---|
625 | if(b.deg<=this.deg){
|
---|
626 | c = new ComplexPoly(this.deg);
|
---|
627 | for(int i=b.deg+1; i<=this.deg; i++)c.coeff[i] = Complex.copy(this.coeff[i]);
|
---|
628 | for(int i=0; i<=b.deg; i++)c.coeff[i] = this.coeff[i].plus(b.coeff[i]);
|
---|
629 | }
|
---|
630 | else{
|
---|
631 | c = new ComplexPoly(b.deg);
|
---|
632 | for(int i=this.deg+1; i<=b.deg; i++)c.coeff[i] = Complex.copy(b.coeff[i]);
|
---|
633 | for(int i=0; i<=this.deg; i++)c.coeff[i] = this.coeff[i].plus(b.coeff[i]);
|
---|
634 | }
|
---|
635 | return c;
|
---|
636 | }
|
---|
637 |
|
---|
638 | // Addition, static method
|
---|
639 | public static ComplexPoly plus(ComplexPoly a, ComplexPoly b){
|
---|
640 | ComplexPoly c = null;
|
---|
641 | if(b.deg<=a.deg){
|
---|
642 | c = new ComplexPoly(a.deg);
|
---|
643 | for(int i=b.deg+1; i<=a.deg; i++)c.coeff[i] = Complex.copy(a.coeff[i]);
|
---|
644 | for(int i=0; i<=b.deg; i++)c.coeff[i] = a.coeff[i].plus(b.coeff[i]);
|
---|
645 | }
|
---|
646 | else{
|
---|
647 | c = new ComplexPoly(b.deg);
|
---|
648 | for(int i=a.deg+1; i<=b.deg; i++)c.coeff[i] = Complex.copy(b.coeff[i]);
|
---|
649 | for(int i=0; i<=a.deg; i++)c.coeff[i] = a.coeff[i].plus(b.coeff[i]);
|
---|
650 | }
|
---|
651 | return c;
|
---|
652 | }
|
---|
653 |
|
---|
654 | // Addition of a Complex, instance method
|
---|
655 | public ComplexPoly plus(Complex bb){
|
---|
656 | ComplexPoly b = new ComplexPoly(bb);
|
---|
657 | return this.plus(b);
|
---|
658 | }
|
---|
659 |
|
---|
660 | // Addition of a Complex, static method
|
---|
661 | public static ComplexPoly plus(ComplexPoly a, Complex bb){
|
---|
662 | ComplexPoly b = new ComplexPoly(bb);
|
---|
663 | return ComplexPoly.plus(a, b);
|
---|
664 | }
|
---|
665 |
|
---|
666 | // Addition of a double, instance method
|
---|
667 | public ComplexPoly plus(double bb){
|
---|
668 | ComplexPoly b = new ComplexPoly(new Complex(bb, 0.0));
|
---|
669 | return this.plus(b);
|
---|
670 | }
|
---|
671 |
|
---|
672 | // Addition of a double, static method
|
---|
673 | public static ComplexPoly plus(ComplexPoly a, double bb){
|
---|
674 | ComplexPoly b = new ComplexPoly(new Complex(bb, 0.0));
|
---|
675 | return ComplexPoly.plus(a, b);
|
---|
676 | }
|
---|
677 |
|
---|
678 | // Addition of an int, instance method
|
---|
679 | public ComplexPoly plus(int bb){
|
---|
680 | ComplexPoly b = new ComplexPoly(new Complex((double) bb, 0.0));
|
---|
681 | return this.plus(b);
|
---|
682 | }
|
---|
683 |
|
---|
684 | // Addition of an int, static method
|
---|
685 | public static ComplexPoly plus(ComplexPoly a, int bb){
|
---|
686 | ComplexPoly b = new ComplexPoly(new Complex((double)bb, 0.0));
|
---|
687 | return ComplexPoly.plus(a, b);
|
---|
688 | }
|
---|
689 |
|
---|
690 | // SUBTRACTION OF TWO POLYNOMIALS
|
---|
691 | // Subtraction, instance method
|
---|
692 | public ComplexPoly minus(ComplexPoly b){
|
---|
693 |
|
---|
694 | ComplexPoly c = null;
|
---|
695 | if(b.deg<=this.deg){
|
---|
696 | c = new ComplexPoly(this.deg);
|
---|
697 | for(int i=b.deg+1; i<=this.deg; i++)c.coeff[i] = Complex.copy(this.coeff[i]);
|
---|
698 | for(int i=0; i<=b.deg; i++)c.coeff[i] = this.coeff[i].minus(b.coeff[i]);
|
---|
699 | }
|
---|
700 | else{
|
---|
701 | c = new ComplexPoly(b.deg);
|
---|
702 | for(int i=this.deg+1; i<=b.deg; i++)c.coeff[i] = b.coeff[i].times(Complex.minusOne());
|
---|
703 | for(int i=0; i<=this.deg; i++)c.coeff[i] = this.coeff[i].minus(b.coeff[i]);
|
---|
704 | }
|
---|
705 | return c;
|
---|
706 | }
|
---|
707 |
|
---|
708 | // Subtraction of a Complex, instance method
|
---|
709 | public ComplexPoly minus(Complex bb){
|
---|
710 | ComplexPoly b = new ComplexPoly(bb);
|
---|
711 | return this.minus(b);
|
---|
712 | }
|
---|
713 |
|
---|
714 | // Subtraction of a double, instance method
|
---|
715 | public ComplexPoly minus(double bb){
|
---|
716 | ComplexPoly b = new ComplexPoly(new Complex(bb, 0.0));
|
---|
717 | return this.minus(b);
|
---|
718 | }
|
---|
719 |
|
---|
720 | // Subtraction of an int, instance method
|
---|
721 | public ComplexPoly minus(int bb){
|
---|
722 | ComplexPoly b = new ComplexPoly(new Complex((double)bb, 0.0));
|
---|
723 | return this.minus(b);
|
---|
724 | }
|
---|
725 |
|
---|
726 | // Subtraction, static method
|
---|
727 | public static ComplexPoly minus(ComplexPoly a, ComplexPoly b){
|
---|
728 | ComplexPoly c = null;
|
---|
729 | if(b.deg<=a.deg){
|
---|
730 | c = new ComplexPoly(a.deg);
|
---|
731 | for(int i=b.deg+1; i<=a.deg; i++)c.coeff[i] = Complex.copy(a.coeff[i]);
|
---|
732 | for(int i=0; i<=b.deg; i++)c.coeff[i] = a.coeff[i].minus(b.coeff[i]);
|
---|
733 | }
|
---|
734 | else{
|
---|
735 | c = new ComplexPoly(b.deg);
|
---|
736 | for(int i=a.deg+1; i<=b.deg; i++)c.coeff[i] = b.coeff[i].times(Complex.minusOne());
|
---|
737 | for(int i=0; i<=a.deg; i++)c.coeff[i] = a.coeff[i].minus(b.coeff[i]);
|
---|
738 | }
|
---|
739 | return c;
|
---|
740 | }
|
---|
741 |
|
---|
742 | // Subtraction of a Complex, static method
|
---|
743 | public static ComplexPoly minus(ComplexPoly a, Complex bb){
|
---|
744 | ComplexPoly b = new ComplexPoly(bb);
|
---|
745 | return ComplexPoly.minus(a, b);
|
---|
746 | }
|
---|
747 |
|
---|
748 |
|
---|
749 | // Subtraction of a double, static method
|
---|
750 | public static ComplexPoly minus(ComplexPoly a, double bb){
|
---|
751 | ComplexPoly b = new ComplexPoly(new Complex(bb, 0.0));
|
---|
752 | return ComplexPoly.minus(a, b);
|
---|
753 | }
|
---|
754 |
|
---|
755 | // Subtraction of an int, static method
|
---|
756 | public static ComplexPoly minus(ComplexPoly a, int bb){
|
---|
757 | ComplexPoly b = new ComplexPoly(new Complex((double)bb, 0.0));
|
---|
758 | return ComplexPoly.minus(a, b);
|
---|
759 | }
|
---|
760 |
|
---|
761 |
|
---|
762 | // MULTIPLICATION OF TWO POLYNOMIALS
|
---|
763 | // Multiplication, instance method
|
---|
764 | public ComplexPoly times(ComplexPoly b){
|
---|
765 | int n = this.deg + b.deg;
|
---|
766 | ComplexPoly c = new ComplexPoly(n);
|
---|
767 | for(int i=0; i<=this.deg; i++){
|
---|
768 | for(int j=0; j<=b.deg; j++){
|
---|
769 | c.coeff[i+j].plusEquals(this.coeff[i].times(b.coeff[j]));
|
---|
770 | }
|
---|
771 | }
|
---|
772 | return c;
|
---|
773 | }
|
---|
774 |
|
---|
775 | // Multiplication, static method
|
---|
776 | public static ComplexPoly times(ComplexPoly a, ComplexPoly b){
|
---|
777 | int n = a.deg + b.deg;
|
---|
778 | ComplexPoly c = new ComplexPoly(n);
|
---|
779 | for(int i=0; i<=a.deg; i++){
|
---|
780 | for(int j=0; j<=b.deg; j++){
|
---|
781 | c.coeff[i+j].plusEquals(a.coeff[i].times(b.coeff[j]));
|
---|
782 | }
|
---|
783 | }
|
---|
784 | return c;
|
---|
785 | }
|
---|
786 |
|
---|
787 | // Multiplication by a Complex, instance method
|
---|
788 | public ComplexPoly times(Complex bb){
|
---|
789 | ComplexPoly c = new ComplexPoly(this.deg);
|
---|
790 | for(int i=0; i<=this.deg; i++){
|
---|
791 | c.coeff[i] = this.coeff[i].times(bb);
|
---|
792 | }
|
---|
793 | return c;
|
---|
794 | }
|
---|
795 |
|
---|
796 | // Multiplication by a Complex, static method
|
---|
797 | public static ComplexPoly times(ComplexPoly a, Complex bb){
|
---|
798 | ComplexPoly c = new ComplexPoly(a.deg);
|
---|
799 | for(int i=0; i<=a.deg; i++){
|
---|
800 | c.coeff[i] = a.coeff[i].times(bb);
|
---|
801 | }
|
---|
802 | return c;
|
---|
803 | }
|
---|
804 |
|
---|
805 | // Multiplication by a double, instance method
|
---|
806 | public ComplexPoly times(double bb){
|
---|
807 | ComplexPoly c = new ComplexPoly(this.deg);
|
---|
808 | for(int i=0; i<=this.deg; i++){
|
---|
809 | c.coeff[i] = this.coeff[i].times(new Complex(bb, 0.0));
|
---|
810 | }
|
---|
811 | return c;
|
---|
812 | }
|
---|
813 |
|
---|
814 | // Multiplication by a double, static method
|
---|
815 | public static ComplexPoly times(ComplexPoly a, double bb){
|
---|
816 | ComplexPoly c = new ComplexPoly(a.deg);
|
---|
817 | for(int i=0; i<=a.deg; i++){
|
---|
818 | c.coeff[i] = a.coeff[i].times(new Complex(bb, 0.0));
|
---|
819 | }
|
---|
820 | return c;
|
---|
821 | }
|
---|
822 |
|
---|
823 | // DERIVATIVES
|
---|
824 | // Return the coefficients, as a new ComplexPoly, of the nth derivative
|
---|
825 | public ComplexPoly nthDerivative(int n){
|
---|
826 | ComplexPoly dnydxn;
|
---|
827 |
|
---|
828 | if(n>this.deg){
|
---|
829 | dnydxn = new ComplexPoly(0.0);
|
---|
830 | }
|
---|
831 | else{
|
---|
832 | dnydxn = new ComplexPoly(this.deg-n);
|
---|
833 | Complex[] nc = Complex.oneDarray(this.deg - n + 1);
|
---|
834 |
|
---|
835 | int k = this.deg - n;
|
---|
836 | for(int i=this.deg; i>n-1; i--){
|
---|
837 | nc[k]=Complex.copy(this.coeff[i]);
|
---|
838 | for(int j=0; j<n; j++){
|
---|
839 | nc[k]=Complex.times(nc[k], i-j);
|
---|
840 | }
|
---|
841 | k--;
|
---|
842 | }
|
---|
843 | dnydxn = new ComplexPoly(nc);
|
---|
844 | }
|
---|
845 | return dnydxn;
|
---|
846 | }
|
---|
847 |
|
---|
848 | // EVALUATION OF A POLYNOMIAL AND ITS DERIVATIVES
|
---|
849 | // Evaluate the polynomial
|
---|
850 | public Complex evaluate(Complex x){
|
---|
851 | Complex y = new Complex();
|
---|
852 | if(this.deg==0){
|
---|
853 | y=Complex.copy(this.coeff[0]);
|
---|
854 | }
|
---|
855 | else{
|
---|
856 | y=Complex.copy(this.coeff[deg]);
|
---|
857 | for(int i=deg-1; i>=0; i--){
|
---|
858 | y=Complex.plus(Complex.times(y, x),this.coeff[i]);
|
---|
859 | }
|
---|
860 | }
|
---|
861 | return y;
|
---|
862 | }
|
---|
863 |
|
---|
864 | public Complex evaluate(double xx){
|
---|
865 | Complex x =new Complex(xx,0.0);
|
---|
866 | Complex y = new Complex();
|
---|
867 | if(deg==0){
|
---|
868 | y=Complex.copy(this.coeff[0]);
|
---|
869 | }
|
---|
870 | else{
|
---|
871 | y=Complex.copy(this.coeff[deg]);
|
---|
872 | for(int i=deg-1; i>=0; i--){
|
---|
873 | y=Complex.plus(Complex.times(y, x),this.coeff[i]);
|
---|
874 | }
|
---|
875 | }
|
---|
876 | return y;
|
---|
877 | }
|
---|
878 |
|
---|
879 | // Evaluate the nth derivative of the polynomial
|
---|
880 | public Complex nthDerivEvaluate(int n, Complex x){
|
---|
881 | Complex dnydxn = new Complex();
|
---|
882 | Complex[] nc = Complex.oneDarray(this.deg+1);
|
---|
883 |
|
---|
884 | if(n==0)
|
---|
885 | {
|
---|
886 | dnydxn=this.evaluate(x);
|
---|
887 | System.out.println("n = 0 in ComplexPoly.nthDerivative");
|
---|
888 | System.out.println("polynomial itself evaluated and returned");
|
---|
889 | }
|
---|
890 | else{
|
---|
891 | ComplexPoly nthderiv = this.nthDerivative(n);
|
---|
892 | dnydxn = nthderiv.evaluate(x);
|
---|
893 | }
|
---|
894 | return dnydxn;
|
---|
895 | }
|
---|
896 |
|
---|
897 | public Complex nthDerivEvaluate(int n, double xx){
|
---|
898 | Complex x = new Complex(xx, 0.0);
|
---|
899 | return nthDerivEvaluate(n, x);
|
---|
900 | }
|
---|
901 |
|
---|
902 | // ROOTS OF POLYNOMIALS
|
---|
903 | // For general details of root searching and a discussion of the rounding errors
|
---|
904 | // see Numerical Recipes, The Art of Scientific Computing
|
---|
905 | // by W H Press, S A Teukolsky, W T Vetterling & B P Flannery
|
---|
906 | // Cambridge University Press, http://www.nr.com/
|
---|
907 |
|
---|
908 | // Calculate the roots (real or complex) of a polynomial (real or complex)
|
---|
909 | // polish = true ([for deg>3 see laguerreAll(...)]
|
---|
910 | // initial root estimates are all zero [for deg>3 see laguerreAll(...)]
|
---|
911 | public Complex[] roots(){
|
---|
912 | boolean polish=true;
|
---|
913 | Complex estx = new Complex(0.0, 0.0);
|
---|
914 | return roots(polish, estx);
|
---|
915 | }
|
---|
916 |
|
---|
917 | // Calculate the roots - as above with the exception that the error messages are suppressed
|
---|
918 | // Required by BlackBox
|
---|
919 | public Complex[] rootsNoMessages(){
|
---|
920 | this.suppressRootsErrorMessages = true;
|
---|
921 | return roots();
|
---|
922 | }
|
---|
923 |
|
---|
924 | // Calculate the roots (real or complex) of a polynomial (real or complex)
|
---|
925 | // initial root estimates are all zero [for deg>3 see laguerreAll(...)]
|
---|
926 | // for polish see laguerreAll(...)[for deg>3]
|
---|
927 | public Complex[] roots(boolean polish){
|
---|
928 | Complex estx = new Complex(0.0, 0.0);
|
---|
929 | return roots(polish, estx);
|
---|
930 | }
|
---|
931 |
|
---|
932 | // Calculate the roots (real or complex) of a polynomial (real or complex)
|
---|
933 | // for estx see laguerreAll(...)[for deg>3]
|
---|
934 | // polish = true see laguerreAll(...)[for deg>3]
|
---|
935 | public Complex[] roots(Complex estx){
|
---|
936 | boolean polish=true;
|
---|
937 | return roots(polish, estx);
|
---|
938 | }
|
---|
939 |
|
---|
940 | // Calculate the roots (real or complex) of a polynomial (real or complex)
|
---|
941 | public Complex[] roots(boolean polish, Complex estx){
|
---|
942 |
|
---|
943 | Complex[] roots = Complex.oneDarray(this.deg);
|
---|
944 |
|
---|
945 | // degree = 0 - no roots
|
---|
946 | if(this.deg==0 && !this.suppressRootsErrorMessages){
|
---|
947 | System.out.println("degree of the polynomial is zero in the method ComplexPoly.roots");
|
---|
948 | System.out.println("null returned");
|
---|
949 | return null;
|
---|
950 | }
|
---|
951 |
|
---|
952 | // Check for no roots, i.e all coefficients but the first = 0
|
---|
953 | int counter = 0;
|
---|
954 | for(int i=1; i<=this.deg; i++){
|
---|
955 | if(this.coeff[i].isZero())counter++;
|
---|
956 | }
|
---|
957 | if(counter==this.deg && !this.suppressRootsErrorMessages){
|
---|
958 | System.out.println("polynomial coefficients above the zeroth are all zero in the method ComplexPoly.roots");
|
---|
959 | System.out.println("null returned");
|
---|
960 | return null;
|
---|
961 | }
|
---|
962 |
|
---|
963 | // Check for only one non-zero coefficient
|
---|
964 | int nonzerocoeff = 0;
|
---|
965 | int nonzeroindex = 0;
|
---|
966 | for(int i=0; i<=this.deg; i++){
|
---|
967 | if(!this.coeff[i].isZero()){
|
---|
968 | nonzerocoeff++;
|
---|
969 | nonzeroindex = i;
|
---|
970 | }
|
---|
971 | }
|
---|
972 | if(nonzerocoeff==1){
|
---|
973 | System.out.println("all polynomial coefficients except a[" + nonzeroindex + "] are zero in the method ComplexPoly.roots");
|
---|
974 | System.out.println("all roots returned as zero");
|
---|
975 | for(int i=0; i<this.deg; i++){
|
---|
976 | roots[i] = Complex.zero();
|
---|
977 | }
|
---|
978 | return roots;
|
---|
979 | }
|
---|
980 |
|
---|
981 | // check for leading zeros
|
---|
982 | boolean testzero=true;
|
---|
983 | int ii=0, nzeros=0;
|
---|
984 | while(testzero){
|
---|
985 | if(this.coeff[ii].isZero()){
|
---|
986 | nzeros++;
|
---|
987 | ii++;
|
---|
988 | }
|
---|
989 | else{
|
---|
990 | testzero=false;
|
---|
991 | }
|
---|
992 | }
|
---|
993 | if(nzeros>0){
|
---|
994 | this.degwz = this.deg - nzeros;
|
---|
995 | this.coeffwz = Complex.oneDarray(this.degwz+1);
|
---|
996 | for(int i=0; i<=this.degwz; i++)this.coeffwz[i] = this.coeff[i + nzeros].copy();
|
---|
997 | }
|
---|
998 | else{
|
---|
999 | this.degwz = this.deg;
|
---|
1000 | this.coeffwz = Complex.oneDarray(this.degwz+1);
|
---|
1001 | for(int i=0; i<=this.degwz; i++)this.coeffwz[i] = this.coeff[i].copy();
|
---|
1002 | }
|
---|
1003 |
|
---|
1004 | // calculate non-zero roots
|
---|
1005 | Complex[] root = Complex.oneDarray(this.degwz);
|
---|
1006 |
|
---|
1007 | switch(this.degwz){
|
---|
1008 | case 1: root[0]=Complex.negate(this.coeffwz[0].over(this.coeffwz[1]));
|
---|
1009 | break;
|
---|
1010 | case 2: root=quadratic(this.coeffwz[0],this.coeffwz[1],this.coeffwz[2]);
|
---|
1011 | break;
|
---|
1012 | case 3: root=cubic(this.coeffwz[0],this.coeffwz[1],this.coeffwz[2], this.coeffwz[3]);
|
---|
1013 | break;
|
---|
1014 | default: root=laguerreAll(polish, estx);
|
---|
1015 | }
|
---|
1016 |
|
---|
1017 | for(int i=0; i<this.degwz; i++){
|
---|
1018 | roots[i]=root[i].copy();
|
---|
1019 | }
|
---|
1020 | if(nzeros>0){
|
---|
1021 | for(int i=this.degwz; i<this.deg; i++){
|
---|
1022 | roots[i]=Complex.zero();
|
---|
1023 | }
|
---|
1024 | }
|
---|
1025 | return roots;
|
---|
1026 | }
|
---|
1027 |
|
---|
1028 | // ROOTS OF A QUADRATIC EQUATION
|
---|
1029 | // ax^2 + bx + c = 0
|
---|
1030 | // roots returned in root[]
|
---|
1031 | // 4ac << b*b accomodated by these methods
|
---|
1032 | public static Complex[] quadratic(Complex c, Complex b, Complex a){
|
---|
1033 | double qsign=1.0;
|
---|
1034 | Complex qsqrt = new Complex();
|
---|
1035 | Complex qtest = new Complex();
|
---|
1036 | Complex bconj = new Complex();
|
---|
1037 | Complex[] root = Complex.oneDarray(2);
|
---|
1038 |
|
---|
1039 | bconj = b.conjugate();
|
---|
1040 | qsqrt = Complex.sqrt((Complex.square(b)).minus((a.times(c)).times(4)));
|
---|
1041 |
|
---|
1042 | qtest = bconj.times(qsqrt);
|
---|
1043 |
|
---|
1044 | if( qtest.getReal() < 0.0 ) qsign = -1.0;
|
---|
1045 |
|
---|
1046 | qsqrt = ((qsqrt.over(qsign)).plus(b)).over(-2.0);
|
---|
1047 | root[0] = Complex.over(qsqrt, a);
|
---|
1048 | root[1] = Complex.over(c, qsqrt);
|
---|
1049 |
|
---|
1050 | return root;
|
---|
1051 | }
|
---|
1052 |
|
---|
1053 | public static Complex[] quadratic(double c, double b, double a){
|
---|
1054 | Complex aa = new Complex(a, 0.0);
|
---|
1055 | Complex bb = new Complex(b, 0.0);
|
---|
1056 | Complex cc = new Complex(c, 0.0);
|
---|
1057 |
|
---|
1058 | return quadratic(cc, bb, aa);
|
---|
1059 | }
|
---|
1060 |
|
---|
1061 | // ROOTS OF A CUBIC EQUATION
|
---|
1062 | // ddx^3 + ccx^2 + bbx + aa = 0
|
---|
1063 | // roots returned in root[]
|
---|
1064 |
|
---|
1065 | // ROOTS OF A CUBIC EQUATION
|
---|
1066 | // ddx^3 + ccx^2 + bbx + aa = 0
|
---|
1067 | // roots returned in root[]
|
---|
1068 | public static Complex[] cubic(Complex aa, Complex bb, Complex cc, Complex dd){
|
---|
1069 |
|
---|
1070 | Complex a = cc.over(dd);
|
---|
1071 | Complex b = bb.over(dd);
|
---|
1072 | Complex c = aa.over(dd);
|
---|
1073 |
|
---|
1074 | Complex[] roots = Complex.oneDarray(3);
|
---|
1075 |
|
---|
1076 | Complex bigQ = ((a.times(a)).minus(b.times(3.0))).over(9.0);
|
---|
1077 | Complex bigR = (((((a.times(a)).times(a)).times(2.0)).minus((a.times(b)).times(9.0))).plus(c.times(27.0))).over(54.0);
|
---|
1078 |
|
---|
1079 | Complex sign = Complex.plusOne();
|
---|
1080 | Complex bigAsqrtTerm = Complex.sqrt((bigR.times(bigR)).minus((bigQ.times(bigQ)).times(bigQ)));
|
---|
1081 | Complex bigRconjugate = bigR.conjugate();
|
---|
1082 | if((bigRconjugate.times(bigAsqrtTerm)).getReal()<0.0)sign = Complex.minusOne();
|
---|
1083 | Complex bigA = (Complex.pow(bigR.plus(sign.times(bigAsqrtTerm)), 1.0/3.0)).times(Complex.minusOne());
|
---|
1084 | Complex bigB = null;
|
---|
1085 | if(bigA.isZero()){
|
---|
1086 | bigB = Complex.zero();
|
---|
1087 | }
|
---|
1088 | else{
|
---|
1089 | bigB = bigQ.over(bigA);
|
---|
1090 | }
|
---|
1091 | Complex aPlusB = bigA.plus(bigB);
|
---|
1092 | Complex aMinusB = bigA.minus(bigB);
|
---|
1093 | Complex minusAplusB = aPlusB.times(Complex.minusOne());
|
---|
1094 | Complex aOver3 = a.over(3.0);
|
---|
1095 | Complex isqrt3over2 = new Complex(0.0, Math.sqrt(3.0)/2.0);
|
---|
1096 | roots[0] = aPlusB.minus(aOver3);
|
---|
1097 | roots[1] = ((minusAplusB.over(2.0)).minus(aOver3)).plus(isqrt3over2.times(aMinusB));
|
---|
1098 | roots[2] = ((minusAplusB.over(2.0)).minus(aOver3)).minus(isqrt3over2.times(aMinusB));
|
---|
1099 |
|
---|
1100 | return roots;
|
---|
1101 | }
|
---|
1102 |
|
---|
1103 | public static Complex[] cubic(double d, double c, double b, double a){
|
---|
1104 | Complex aa = new Complex(a, 0.0);
|
---|
1105 | Complex bb = new Complex(b, 0.0);
|
---|
1106 | Complex cc = new Complex(c, 0.0);
|
---|
1107 | Complex dd = new Complex(d, 0.0);
|
---|
1108 |
|
---|
1109 | return cubic(dd, cc, bb, aa);
|
---|
1110 | }
|
---|
1111 |
|
---|
1112 | // LAGUERRE'S METHOD FOR COMPLEX ROOTS OF A COMPLEX POLYNOMIAL
|
---|
1113 |
|
---|
1114 | // Laguerre method for one of the roots
|
---|
1115 | // Following the procedure in Numerical Recipes for C [Reference above]
|
---|
1116 | // estx estimate of the root
|
---|
1117 | // coeff[] coefficients of the polynomial
|
---|
1118 | // m degree of the polynomial
|
---|
1119 | public static Complex laguerre(Complex estx, Complex[] pcoeff, int m){
|
---|
1120 | double eps = 1e-7; // estimated fractional round-off error
|
---|
1121 | int mr = 8; // number of fractional values in Adam's method of breaking a limit cycle
|
---|
1122 | int mt = 1000; // number of steps in breaking a limit cycle
|
---|
1123 | int maxit = mr*mt; // maximum number of iterations allowed
|
---|
1124 | int niter = 0; // number of iterations taken
|
---|
1125 |
|
---|
1126 | // fractions used to break a limit cycle
|
---|
1127 | double frac[]={0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88, 1.0};
|
---|
1128 |
|
---|
1129 | Complex root = new Complex(); // root
|
---|
1130 | Complex b = new Complex();
|
---|
1131 | Complex d = new Complex();
|
---|
1132 | Complex f = new Complex();
|
---|
1133 | Complex g = new Complex();
|
---|
1134 | Complex g2 = new Complex();
|
---|
1135 | Complex h = new Complex();
|
---|
1136 | Complex sq = new Complex();
|
---|
1137 | Complex gp = new Complex();
|
---|
1138 | Complex gm = new Complex();
|
---|
1139 | Complex dx = new Complex();
|
---|
1140 | Complex x1 = new Complex();
|
---|
1141 | Complex temp1 = new Complex();
|
---|
1142 | Complex temp2 = new Complex();
|
---|
1143 |
|
---|
1144 | double abp = 0.0D, abm = 0.0D;
|
---|
1145 | double err = 0.0D, abx = 0.0D;
|
---|
1146 |
|
---|
1147 | for(int i=1; i<=maxit; i++){
|
---|
1148 | niter=i;
|
---|
1149 | b=Complex.copy(pcoeff[m]);
|
---|
1150 | err=Complex.abs(b);
|
---|
1151 | d=f=Complex.zero();
|
---|
1152 | abx=Complex.abs(estx);
|
---|
1153 | for(int j=m-1; j>=0;j--)
|
---|
1154 | {
|
---|
1155 | // Efficient computation of the polynomial and its first two derivatives
|
---|
1156 | f=Complex.plus(Complex.times(estx, f), d);
|
---|
1157 | d=Complex.plus(Complex.times(estx, d), b);
|
---|
1158 | b=Complex.plus(Complex.times(estx, b), pcoeff[j]);
|
---|
1159 | err=Complex.abs(b)+abx*err;
|
---|
1160 | }
|
---|
1161 | err*=eps;
|
---|
1162 |
|
---|
1163 | // Estimate of round-off error in evaluating polynomial
|
---|
1164 | if(Complex.abs(b)<=err)
|
---|
1165 | {
|
---|
1166 | root=Complex.copy(estx);
|
---|
1167 | niter=i;
|
---|
1168 | return root;
|
---|
1169 | }
|
---|
1170 | // Laguerre formula
|
---|
1171 | g=Complex.over(d, b);
|
---|
1172 | g2=Complex.square(g);
|
---|
1173 | h=Complex.minus(g2, Complex.times(2.0, Complex.over(f, b)));
|
---|
1174 | sq=Complex.sqrt(Complex.times((double)(m-1), Complex.minus(Complex.times((double)m, h), g2)));
|
---|
1175 | gp=Complex.plus(g, sq);
|
---|
1176 | gm=Complex.minus(g, sq);
|
---|
1177 | abp=Complex.abs(gp);
|
---|
1178 | abm=Complex.abs(gm);
|
---|
1179 | if( abp < abm ) gp = gm;
|
---|
1180 | temp1.setReal((double)m);
|
---|
1181 | temp2.setReal(Math.cos((double)i));
|
---|
1182 | temp2.setImag(Math.sin((double)i));
|
---|
1183 | dx=((Math.max(abp, abm) > 0.0 ? Complex.over(temp1, gp):Complex.times(Math.exp(1.0+abx),temp2)));
|
---|
1184 | x1=Complex.minus(estx, dx);
|
---|
1185 | if(Complex.isEqual(estx, x1))
|
---|
1186 | {
|
---|
1187 | root=Complex.copy(estx);
|
---|
1188 | niter=i;
|
---|
1189 | return root; // converged
|
---|
1190 | }
|
---|
1191 | if ((i % mt)!= 0){
|
---|
1192 | estx=Complex.copy(x1);
|
---|
1193 | }
|
---|
1194 | else{
|
---|
1195 | // Every so often we take a fractional step to break any limit cycle
|
---|
1196 | // (rare occurence)
|
---|
1197 | estx=Complex.minus(estx, Complex.times(frac[i/mt-1], dx));
|
---|
1198 | }
|
---|
1199 | niter=i;
|
---|
1200 | }
|
---|
1201 | // exceeded maximum allowed iterations
|
---|
1202 | root=Complex.copy(estx);
|
---|
1203 | System.out.println("Maximum number of iterations exceeded in laguerre");
|
---|
1204 | System.out.println("root returned at this point");
|
---|
1205 | return root;
|
---|
1206 | }
|
---|
1207 |
|
---|
1208 | // Finds all roots of a complex polynomial by successive calls to laguerre
|
---|
1209 | // Following the procedure in Numerical Recipes for C [Reference above]
|
---|
1210 | // Initial estimates are all zero, polish=true
|
---|
1211 | public Complex[] laguerreAll(){
|
---|
1212 | Complex estx = new Complex(0.0, 0.0);
|
---|
1213 | boolean polish = true;
|
---|
1214 | return laguerreAll(polish, estx);
|
---|
1215 | }
|
---|
1216 |
|
---|
1217 | // Initial estimates estx, polish=true
|
---|
1218 | public Complex[] laguerreAll(Complex estx){
|
---|
1219 | boolean polish = true;
|
---|
1220 | return laguerreAll(polish, estx);
|
---|
1221 | }
|
---|
1222 |
|
---|
1223 | // Initial estimates are all zero.
|
---|
1224 | public Complex[] laguerreAll(boolean polish){
|
---|
1225 | Complex estx = new Complex(0.0, 0.0);
|
---|
1226 | return laguerreAll(polish, estx);
|
---|
1227 | }
|
---|
1228 |
|
---|
1229 | // Finds all roots of a complex polynomial by successive calls to laguerre
|
---|
1230 | // Initial estimates are estx
|
---|
1231 | public Complex[] laguerreAll(boolean polish, Complex estx){
|
---|
1232 | // polish boolean variable
|
---|
1233 | // if true roots polished also by Laguerre
|
---|
1234 | // if false roots returned to be polished by another method elsewhere.
|
---|
1235 | // estx estimate of root - Preferred default value is zero to favour convergence
|
---|
1236 | // to smallest remaining root
|
---|
1237 |
|
---|
1238 | int m = this.degwz;
|
---|
1239 | double eps = 2.0e-6; // tolerance in determining round off in imaginary part
|
---|
1240 |
|
---|
1241 | Complex x = new Complex();
|
---|
1242 | Complex b = new Complex();
|
---|
1243 | Complex c = new Complex();
|
---|
1244 | Complex[] ad = new Complex[m+1];
|
---|
1245 | Complex[] roots = new Complex[m+1];
|
---|
1246 |
|
---|
1247 | // Copy polynomial for successive deflation
|
---|
1248 | for(int j=0; j<=m; j++) ad[j]=Complex.copy(this.coeffwz[j]);
|
---|
1249 |
|
---|
1250 | // Loop over each root found
|
---|
1251 | for(int j=m; j>=1; j--){
|
---|
1252 | x=Complex.copy(estx); // Preferred default value is zero to favour convergence to smallest remaining root
|
---|
1253 | // and find the root
|
---|
1254 | x=laguerre(x, ad, j);
|
---|
1255 | if(Math.abs(x.getImag())<=2.0*eps*Math.abs(x.getReal())) x.setImag(0.0);
|
---|
1256 | roots[j]=Complex.copy(x);
|
---|
1257 | b=Complex.copy(ad[j]);
|
---|
1258 | for(int jj=j-1; jj>=0; jj--){
|
---|
1259 | c=Complex.copy(ad[jj]);
|
---|
1260 | ad[jj]=Complex.copy(b);
|
---|
1261 | b=(x.times(b)).plus(c);
|
---|
1262 | }
|
---|
1263 | }
|
---|
1264 |
|
---|
1265 | if(polish){
|
---|
1266 | // polish roots using the undeflated coefficients
|
---|
1267 | for(int j=1; j<=m; j++){
|
---|
1268 | roots[j]=laguerre(roots[j], this.coeffwz, m);
|
---|
1269 | }
|
---|
1270 | }
|
---|
1271 |
|
---|
1272 | // Sort roots by their real parts by straight insertion
|
---|
1273 | for(int j=2; j<=m; j++){
|
---|
1274 | x=Complex.copy(roots[j]);
|
---|
1275 | int i=0;
|
---|
1276 | for(i=j-1; i>=1; i--){
|
---|
1277 | if(roots[i].getReal() <= x.getReal()) break;
|
---|
1278 | roots[i+1]=Complex.copy(roots[i]);
|
---|
1279 | }
|
---|
1280 | roots[i+1]=Complex.copy(x);
|
---|
1281 | }
|
---|
1282 | // shift roots to zero initial index
|
---|
1283 | for(int i=0; i<m; i++)roots[i]=Complex.copy(roots[i+1]);
|
---|
1284 | return roots;
|
---|
1285 | }
|
---|
1286 | }
|
---|
1287 |
|
---|