source: src/main/java/agents/anac/y2015/agentBuyogV2/flanagan/complex/ComplexPoly.java

Last change on this file was 127, checked in by Wouter Pasman, 6 years ago

#41 ROLL BACK of rev.126 . So this version is equal to rev. 125

File size: 52.5 KB
Line 
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
40package agents.anac.y2015.agentBuyogV2.flanagan.complex;
41
42import java.util.ArrayList;
43
44import agents.anac.y2015.agentBuyogV2.flanagan.circuits.Phasor;
45import agents.anac.y2015.agentBuyogV2.flanagan.control.BlackBox;
46import agents.anac.y2015.agentBuyogV2.flanagan.io.FileOutput;
47import agents.anac.y2015.agentBuyogV2.flanagan.math.ArrayMaths;
48import agents.anac.y2015.agentBuyogV2.flanagan.math.Conv;
49import agents.anac.y2015.agentBuyogV2.flanagan.math.Fmath;
50import agents.anac.y2015.agentBuyogV2.flanagan.math.Polynomial;
51
52import java.math.*;
53
54
55public 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
Note: See TracBrowser for help on using the repository browser.