Complex


1   /*****************************************************************************
2   
3   JEP - Java Math Expression Parser 2.3.0
4         October 3 2004
5         (c) Copyright 2004, Nathan Funk and Richard Morris
6         See LICENSE.txt for license information.
7   
8   *****************************************************************************/
9   
10  
11  package  org.nfunk.jep.type;
12  
13  /**
14   * Represents a complex number with double precision real and imaginary
15   * components. Includes complex arithmetic functions.<p>
16   * The two main sources of reference used for creating this class were:<br>
17   * - "Numerical Recipes in C - The Art of Scientific Computing"
18   *    (ISBN 0-521-43108-5) http://www.nr.com and <br>
19   * - The org.netlib.math.complex package (http://www.netlib.org) which was 
20   *   developed by Sandy Anderson and Priyantha Jayanetti (published under
21   *   GPL).<p>
22   * Some of the arithmetic functions in this class are based on the mathematical
23   * equations given in the source of the netlib package. The functions were
24   * validated by comparing results with the netlib complex class.<p>
25   * It is important to note that the netlib complex package is more
26   * extensive and efficient (e.g. Garbage collector friendly) than this
27   * implementation. If high precision and efficiency if of necessity it is 
28   * recommended to use the netlib package.
29   *
30   * @author Nathan Funk
31   * @version 2.3.0 alpha now extends Number, has add and sub methods.
32   * @version 2.3.0 beta 1 now overrides equals and hashCode.
33  /* @version 2.3.0 beta 2 does not implement Number anymore, as caused too many problems.
34   */
35  
36  public class Complex
37  {
38      /** the real component */
39      private double re;
40      
41      /** the imaginary component */
42      private double im;
43      
44  
45  //------------------------------------------------------------------------
46  // Constructors
47  
48      /**
49       * Default constructor.
50       */
51      public Complex() {
52          re = 0;
53          im = 0;
54      }
55  
56      /**
57       * Constructor from a single double value. The complex number is
58       * initialized with the real component equal to the parameter, and
59       * the imaginary component equal to zero.
60       */
61      public Complex(double re_in) {
62          re = re_in;
63          im = 0;
64      }
65  
66      /**
67       * Construct from a Number. This constructor uses the doubleValue()
68       * method of the parameter to initialize the real component of the
69       * complex number. The imaginary component is initialized to zero.
70       */
71      public Complex(Number   re_in) {
72          re = re_in.doubleValue();
73          im = 0;
74      }
75      
76      /**
77       * Copy constructor
78       */
79      public Complex(Complex z) {
80          re = z.re;
81          im = z.im;
82      }
83  
84      /**
85       * Initialize the real and imaginary components to the values given
86       * by the parameters.
87       */
88      public Complex(double re_in, double im_in) {
89          re = re_in;
90          im = im_in;
91      }
92  
93      /**
94       * Returns the real component of this object
95       */
96      public double re() {
97          return re;  
98      }
99  
100     /**
101      * Returns the imaginary component of this object
102      */
103     public double im() {
104         return im;
105     }
106     
107     /**
108      * Copies the values from the parameter object to this object
109      */
110     public void set(Complex z) {
111         re = z.re;
112         im = z.im;
113     }
114     
115     /**
116      * Sets the real and imaginary values of the object.
117      */
118     public void set(double re_in, double im_in) {
119         re = re_in;
120         im = im_in;
121     }
122 
123     /**
124      * Sets the real component of the object
125      */
126     public void setRe(double re_in) {
127         re = re_in;
128     }
129 
130     /**
131      * Sets the imaginary component of the object
132      */
133     public void setIm(double im_in) {
134         im = im_in;
135     }
136 
137 //------------------------------------------------------------------------
138 // Various functions
139 
140     /**
141      * Compares this object with the Complex number given as parameter
142      * <pre>b</pre>. The <pre>tolerance</pre> parameter is the radius
143      * within which the <pre>b</pre> number must lie for the two
144      * complex numbers to be considered equal.
145      *
146      * @return <pre>true</pre> if the complex number are considered equal,
147      * <pre>false</pre> otherwise.
148      */
149     public boolean equals(Complex b, double tolerance) {
150         double temp1 = (re - b.re);
151         double temp2 = (im - b.im);
152         
153         return (temp1*temp1 + temp2*temp2) <= tolerance*tolerance;
154     }
155     /**
156      * Compares this object against the specified object. 
157      * The result is true if and only if the argument is not null 
158      * and is a Complex object that represents the same complex number. 
159      * Equality follows the same pattern as Double aplies to each field:
160      * <ul>
161      * <li>If d1 and d2 both represent Double.NaN, then the equals method returns true, even though Double.NaN==Double.NaN has the value false.
162      * <li>If d1 represents +0.0 while d2 represents -0.0, or vice versa, the equal test has the value false, even though +0.0==-0.0 has the value true.
163      * </ul>
164      * This definition allows hash tables to operate properly.
165 
166      * @since 2.3.0.2
167      */
168     public boolean equals(Object   o) {
169         if(!(o instanceof Complex)) return false;
170         Complex c = (Complex) o;
171         return(Double.doubleToLongBits(this.re) == Double.doubleToLongBits(c.re) 
172             && Double.doubleToLongBits(this.im) == Double.doubleToLongBits(c.im));
173     }
174     /**
175      * Always override hashCode when you override equals.
176      * Efective Java, Joshua Bloch, Sun Press
177      */
178     public int hashCode() {
179         int result = 17;
180         long xl = Double.doubleToLongBits(this.re);
181         long yl = Double.doubleToLongBits(this.im);
182         int xi = (int)(xl^(xl>>32));
183         int yi = (int)(yl^(yl>>32));
184         result = 37*result+xi;
185         result = 37*result+yi;
186         return result;
187     }
188     /**
189      * Returns the value of this complex number as a string in the format:
190      * <pre>(real, imaginary)</pre>.
191      */
192     public String   toString() {
193         return "(" + re + ", " + im + ")";
194     }
195     
196     /**
197      * Returns <tt>true</tt> if either the real or imaginary component of this
198      * <tt>Complex</tt> is an infinite value.
199      *
200      * <p>
201      * @return  <tt>true</tt> if either component of the <tt>Complex</tt> object is infinite; <tt>false</tt>, otherwise.
202      * <p>
203      **/
204     public boolean isInfinite() {
205         return (Double.isInfinite(re) || Double.isInfinite(im));
206     }
207 
208     /**
209      * Returns <tt>true</tt> if either the real or imaginary component of this
210      * <tt>Complex</tt> is a Not-a-Number (<tt>NaN</tt>) value.
211      *
212      * <p>
213      * @return  <tt>true</tt> if either component of the <tt>Complex</tt> object is <tt>NaN</tt>; <tt>false</tt>, otherwise.
214      * <p>
215      **/
216     public boolean isNaN() {
217         return (Double.isNaN(re) || Double.isNaN(im));
218     }
219     
220     /**
221      * Returns the absolute value of the complex number.
222      * <p>
223      * Adapted from Numerical Recipes in C -
224      * The Art of Scientific Computing<br>
225      * ISBN 0-521-43108-5
226      */
227     public double abs() {
228         double absRe = Math.abs(re);
229         double absIm = Math.abs(im);
230         
231         if (absRe == 0 && absIm == 0) {
232             return 0;
233         } else if (absRe>absIm) {
234             double temp = absIm/absRe;
235             return absRe*Math.sqrt(1 + temp*temp);  
236         } else {
237             double temp = absRe/absIm;
238             return absIm*Math.sqrt(1 + temp*temp);
239         }
240     }
241 
242     /**
243      * Returns the square of the absolute value (re*re+im*im).
244      */
245     public double abs2() {
246         return re*re+im*im; 
247     }
248 
249     /**
250      * Returns the argument of this complex number (Math.atan2(re,im))
251      */
252     public double arg() {
253         return Math.atan2(im,re);
254     }
255 
256     /**
257      * Returns the negative value of this complex number.
258      */
259     public Complex neg() {
260         return new Complex(-re,-im);
261     }
262 
263     /**
264      * Multiply the complex number with a double value.
265      * @return The result of the multiplication
266      */
267     public Complex mul(double b) {
268         return new Complex(re*b, im*b);
269     }
270 
271     /**
272      * Adds the complex number with another complex value.
273      * @return The result of the addition
274      * @since 2.3.0.1
275      */
276     public Complex add(Complex b) {
277         return new Complex(re+b.re,im+b.im);
278     }
279 
280     /**
281      * Adds the complex number with another complex value.
282      * @return The result of the addition
283      * @since 2.3.0.1
284      */
285     public Complex sub(Complex b) {
286         return new Complex(re-b.re,im-b.im);
287     }
288     /**
289      * Multiply the complex number with another complex value.
290      * @return The result of the multiplication
291      */
292     public Complex mul(Complex b) {
293         return new Complex(re*b.re - im*b.im,
294                            im*b.re + re*b.im);
295     }
296     
297     /**
298      * Returns the result of dividing this complex number by the parameter.
299      */
300     public Complex div(Complex b) {
301         // Adapted from Numerical Recipes in C - The Art of Scientific Computing
302         // ISBN 0-521-43108-5
303         double resRe, resIm;
304         double r, den;
305         
306         if (Math.abs(b.re) >= Math.abs(b.im)) {
307             r = b.im/b.re;
308             den = b.re + r*b.im;
309             resRe = (re+r*im)/den;
310             resIm = (im-r*re)/den;
311         } else {
312             r = b.re/b.im;
313             den = b.im + r*b.re;
314             resRe = (re*r+im)/den;
315             resIm = (im*r-re)/den;
316         }
317         
318         return new Complex(resRe, resIm);
319     }
320 
321     /**
322       * Returns the value of this complex number raised to the power
323       * of a real component (in double precision).<p>
324       * This method considers special cases where a simpler algorithm
325       * would return "ugly" results.<br>
326       * For example when the expression (-1e40)^0.5 is evaluated without
327       * considering the special case, the argument of the base is the
328       * double number closest to pi. When sin and cos are used for the
329       * final evaluation of the result, the slight difference of the
330       * argument from pi causes a non-zero value for the real component
331       * of the result. Because the value of the base is so high, the error
332       * is magnified.Although the error is normal for floating 
333       * point calculations, the consideration of commonly occuring special
334       * cases improves the accuracy and aesthetics of the results.<p>
335       * If you know a more elegant way to solve this problem, please let
336       * me know at nathanfunk@hotmail.com .
337       */
338     public Complex power(double exponent) {
339         // z^exp = abs(z)^exp * (cos(exp*arg(z)) + i*sin(exp*arg(z)))
340         double scalar = Math.pow(abs(),exponent);
341         boolean specialCase = false;
342         int factor = 0;
343 
344         // consider special cases to avoid floating point errors
345         // for power expressions such as (-1e20)^2
346         if (im==0 && re<0) {specialCase = true; factor = 2;}
347         if (re==0 && im>0) {specialCase = true; factor = 1;}
348         if (re==0 && im<0) {specialCase = true; factor = -1;}
349         
350         if (specialCase && factor*exponent == (int)(factor*exponent)) {
351             short[] cSin = {0,1,0,-1}; //sin of 0, pi/2, pi, 3pi/2
352             short[] cCos = {1,0,-1,0}; //cos of 0, pi/2, pi, 3pi/2
353             
354             int x = ((int)(factor*exponent))%4;
355             if (x<0) x = 4+x;
356             
357             return new Complex(scalar*cCos[x], scalar*cSin[x]);
358         }
359         
360         double temp = exponent * arg();
361         
362         return new Complex(scalar*Math.cos(temp), scalar*Math.sin(temp));
363     }
364     
365     /**
366      * Returns the value of this complex number raised to the power of
367      * a complex exponent
368      */
369     public Complex power(Complex exponent) {
370         if (exponent.im == 0) return power(exponent.re);
371         
372         double temp1Re = Math.log(abs());
373         double temp1Im = arg();
374         
375         double temp2Re = (temp1Re*exponent.re) - (temp1Im*exponent.im);
376         double temp2Im = (temp1Re*exponent.im) + (temp1Im*exponent.re);
377 
378         double scalar = Math.exp(temp2Re);
379         
380         return new Complex(scalar*Math.cos(temp2Im), scalar*Math.sin(temp2Im));
381     }
382 
383     /**
384      * Returns the logarithm of this complex number.
385      */ 
386     public Complex log() {
387         return new Complex(Math.log(abs()), arg()); 
388     }
389     
390     /**
391      * Calculates the square root of this object.
392      * Adapted from Numerical Recipes in C - The Art of Scientific Computing
393      * (ISBN 0-521-43108-5)
394      */
395     public Complex sqrt() {
396         Complex c;
397         double absRe,absIm,w,r;
398         
399         if (re == 0 && im == 0) {
400             c = new Complex(0,0);
401         } else {
402             absRe = Math.abs(re);
403             absIm = Math.abs(im);
404             
405             if (absRe>=absIm) {
406                 r = absIm/absRe;
407                 w = Math.sqrt(absRe)*Math.sqrt(0.5*(1.0+Math.sqrt(1.0+r*r)));
408             } else {
409                 r = absRe/absIm;
410                 w = Math.sqrt(absIm)*Math.sqrt(0.5*(r  +Math.sqrt(1.0+r*r)));
411             }
412             
413             if (re>=0) {
414                 c = new Complex(w, im/(2.0*w));
415             } else {
416                 if (im<0) w = -w;
417                 c = new Complex(im/(2.0*w), w);
418             }
419         }
420         
421         return c;
422     }
423 
424 //------------------------------------------------------------------------
425 // Trigonometric functions
426 
427     /**
428      * Returns the sine of this complex number.
429      */
430     public Complex sin() {
431         double izRe, izIm;
432         double temp1Re, temp1Im;
433         double temp2Re, temp2Im;
434         double scalar;
435         
436         //  sin(z)  =  ( exp(i*z) - exp(-i*z) ) / (2*i)
437         izRe = -im;
438         izIm =  re;
439         
440         // first exp
441         scalar = Math.exp(izRe);
442         temp1Re = scalar * Math.cos(izIm);
443         temp1Im = scalar * Math.sin(izIm);
444         
445         // second exp
446         scalar = Math.exp(-izRe);
447         temp2Re = scalar * Math.cos(-izIm);
448         temp2Im = scalar * Math.sin(-izIm);
449         
450         temp1Re -= temp2Re;
451         temp1Im -= temp2Im;
452         
453         return new Complex(0.5*temp1Im, -0.5*temp1Re);
454     }
455 
456     /**
457      * Returns the cosine of this complex number.
458      */
459     public Complex cos() {
460         double izRe, izIm;
461         double temp1Re, temp1Im;
462         double temp2Re, temp2Im;
463         double scalar;
464         
465         //  cos(z)  =  ( exp(i*z) + exp(-i*z) ) / 2
466         izRe = -im;
467         izIm =  re;
468         
469         // first exp
470         scalar = Math.exp(izRe);
471         temp1Re = scalar * Math.cos(izIm);
472         temp1Im = scalar * Math.sin(izIm);
473         
474         // second exp
475         scalar = Math.exp(-izRe);
476         temp2Re = scalar * Math.cos(-izIm);
477         temp2Im = scalar * Math.sin(-izIm);
478         
479         temp1Re += temp2Re;
480         temp1Im += temp2Im;
481         
482         return new Complex(0.5*temp1Re, 0.5*temp1Im);
483     }
484 
485 
486     /**
487      * Returns the tangent of this complex number.
488      */
489     public Complex tan() {
490         // tan(z) = sin(z)/cos(z)
491         double izRe, izIm;
492         double temp1Re, temp1Im;
493         double temp2Re, temp2Im;
494         double scalar;
495         Complex sinResult, cosResult;
496         
497         //  sin(z)  =  ( exp(i*z) - exp(-i*z) ) / (2*i)
498         izRe = -im;
499         izIm =  re;
500         
501         // first exp
502         scalar = Math.exp(izRe);
503         temp1Re = scalar * Math.cos(izIm);
504         temp1Im = scalar * Math.sin(izIm);
505         
506         // second exp
507         scalar = Math.exp(-izRe);
508         temp2Re = scalar * Math.cos(-izIm);
509         temp2Im = scalar * Math.sin(-izIm);
510         
511         temp1Re -= temp2Re;
512         temp1Im -= temp2Im;
513         
514         sinResult = new Complex(0.5*temp1Re, 0.5*temp1Im);
515 
516         //  cos(z)  =  ( exp(i*z) + exp(-i*z) ) / 2
517         izRe = -im;
518         izIm =  re;
519         
520         // first exp
521         scalar = Math.exp(izRe);
522         temp1Re = scalar * Math.cos(izIm);
523         temp1Im = scalar * Math.sin(izIm);
524         
525         // second exp
526         scalar = Math.exp(-izRe);
527         temp2Re = scalar * Math.cos(-izIm);
528         temp2Im = scalar * Math.sin(-izIm);
529         
530         temp1Re += temp2Re;
531         temp1Im += temp2Im;
532         
533         cosResult = new Complex(0.5*temp1Re, 0.5*temp1Im);
534         
535         return sinResult.div(cosResult);
536     }
537     
538 
539 //------------------------------------------------------------------------
540 // Inverse trigonometric functions
541     
542     public Complex asin() {
543         Complex result;
544         double tempRe, tempIm;
545         
546         //  asin(z)  =  -i * log(i*z + sqrt(1 - z*z))
547 
548         tempRe =  1.0 - ( (re*re) - (im*im) );
549         tempIm =  0.0 - ( (re*im) + (im*re) );
550 
551         result =  new Complex(tempRe, tempIm);
552         result = result.sqrt();
553         
554         result.re += -im;
555         result.im +=  re;
556         
557         tempRe = Math.log(result.abs());
558         tempIm = result.arg();
559         
560         result.re =   tempIm;
561         result.im = - tempRe;
562         
563         return result;
564     }
565 
566     public Complex acos() {
567         Complex result;
568         double tempRe, tempIm;
569         
570         //  acos(z)  =  -i * log( z + i * sqrt(1 - z*z) )
571 
572         tempRe =  1.0 - ( (re*re) - (im*im) );
573         tempIm =  0.0 - ( (re*im) + (im*re) );
574 
575         result =  new Complex(tempRe, tempIm);
576         result = result.sqrt();
577         
578         tempRe = -result.im;
579         tempIm =  result.re;
580         
581         result.re = re + tempRe;
582         result.im = im + tempIm;
583         
584         tempRe = Math.log(result.abs());
585         tempIm = result.arg();
586         
587         result.re =   tempIm;
588         result.im = - tempRe;
589         
590         return result;
591     }
592 
593     public Complex atan() {
594         // atan(z) = -i/2 * log((i-z)/(i+z))
595         
596         double tempRe, tempIm;
597         Complex result = new Complex(-re, 1.0 - im);
598         
599         tempRe = re;
600         tempIm = 1.0 + im;
601         
602         result = result.div(new Complex(tempRe, tempIm));
603         
604         tempRe = Math.log(result.abs());
605         tempIm = result.arg();
606         
607         result.re =  0.5*tempIm;
608         result.im = -0.5*tempRe;
609         
610         return result;
611     }
612 
613 
614 //------------------------------------------------------------------------
615 // Hyperbolic trigonometric functions
616 
617     public Complex sinh() {
618         double scalar;
619         double temp1Re, temp1Im;
620         double temp2Re, temp2Im;
621         //  sinh(z)  =  ( exp(z) - exp(-z) ) / 2
622         
623         // first exp
624         scalar = Math.exp(re);
625         temp1Re = scalar * Math.cos(im);
626         temp1Im = scalar * Math.sin(im);
627         
628         // second exp
629         scalar = Math.exp(-re);
630         temp2Re = scalar * Math.cos(-im);
631         temp2Im = scalar * Math.sin(-im);
632         
633         temp1Re -= temp2Re;
634         temp1Im -= temp2Im;
635         
636         return new Complex(0.5*temp1Re, 0.5*temp1Im);
637     }
638 
639 
640     public Complex cosh() {
641         double scalar;
642         double temp1Re, temp1Im;
643         double temp2Re, temp2Im;
644         //  cosh(z)  =  ( exp(z) + exp(-z) ) / 2
645         
646         // first exp
647         scalar = Math.exp(re);
648         temp1Re = scalar * Math.cos(im);
649         temp1Im = scalar * Math.sin(im);
650         
651         // second exp
652         scalar = Math.exp(-re);
653         temp2Re = scalar * Math.cos(-im);
654         temp2Im = scalar * Math.sin(-im);
655         
656         temp1Re += temp2Re;
657         temp1Im += temp2Im;
658         
659         return new Complex(0.5*temp1Re, 0.5*temp1Im);
660     }
661 
662     public Complex tanh() {
663         double scalar;
664         double temp1Re, temp1Im;
665         double temp2Re, temp2Im;
666         Complex sinRes, cosRes;
667         //  tanh(z)  =  sinh(z) / cosh(z)
668 
669         scalar = Math.exp(re);
670         temp1Re = scalar * Math.cos(im);
671         temp1Im = scalar * Math.sin(im);
672         
673         scalar = Math.exp(-re);
674         temp2Re = scalar * Math.cos(-im);
675         temp2Im = scalar * Math.sin(-im);
676         
677         temp1Re -= temp2Re;
678         temp1Im -= temp2Im;
679         
680         sinRes = new Complex(0.5*temp1Re, 0.5*temp1Im);
681 
682         scalar = Math.exp(re);
683         temp1Re = scalar * Math.cos(im);
684         temp1Im = scalar * Math.sin(im);
685         
686         scalar = Math.exp(-re);
687         temp2Re = scalar * Math.cos(-im);
688         temp2Im = scalar * Math.sin(-im);
689         
690         temp1Re += temp2Re;
691         temp1Im += temp2Im;
692         
693         cosRes = new Complex(0.5*temp1Re, 0.5*temp1Im);
694         
695         return sinRes.div(cosRes);
696     }
697 
698     
699 //------------------------------------------------------------------------
700 // Inverse hyperbolic trigonometric functions
701 
702     public Complex asinh() {
703         Complex result;
704         //  asinh(z)  =  log(z + sqrt(z*z + 1))
705 
706         result = new Complex(
707             ((re*re) - (im*im)) + 1,
708             (re*im) + (im*re));
709 
710         result = result.sqrt();
711         
712         result.re += re;
713         result.im += im;
714         
715         double temp = result.arg();
716         result.re = Math.log(result.abs());
717         result.im = temp;
718         
719         return result;
720     }
721 
722     public Complex acosh() {
723         Complex result;
724         
725         //  acosh(z)  =  log(z + sqrt(z*z - 1))
726 
727         result = new Complex(
728             ((re*re) - (im*im)) - 1,
729             (re*im) + (im*re));
730 
731         result = result.sqrt();
732         
733         result.re += re;
734         result.im += im;
735         
736         double temp = result.arg();
737         result.re = Math.log(result.abs());
738         result.im = temp;
739         
740         return result;
741     }
742 
743     public Complex atanh() {
744         //  atanh(z)  =  1/2 * log( (1+z)/(1-z) )
745         
746         double tempRe, tempIm;
747         Complex result = new Complex(1.0 + re, im);
748         
749         tempRe = 1.0 - re;
750         tempIm =     - im;
751         
752         result = result.div(new Complex(tempRe, tempIm));
753         
754         tempRe = Math.log(result.abs());
755         tempIm = result.arg();
756         
757         result.re = 0.5*tempRe;
758         result.im = 0.5*tempIm;
759         
760         return result;
761     }
762 
763     /**
764      * Converts an [r,theta] pair to a complex number r * e^(i theta).
765      * @param r The radius
766      * @param theta The angle
767      * @return The complex result.
768      * @since 2.3.0.1
769      */
770     public static Complex polarValueOf(Number   r,Number   theta)
771     {
772         double rad = r.doubleValue();
773         double ang = theta.doubleValue();
774         return new Complex(rad*Math.cos(ang), rad*Math.sin(ang));
775         
776     }
777     /** Returns real part.
778      * @since 2.3.0.1
779      */
780     public double doubleValue() {
781         return re;
782     }
783 
784     /** Returns real part.
785      * @since 2.3.0.1
786      */
787     public float floatValue() {
788         return (float) re;
789     }
790 
791     /** Returns real part.
792      * @since 2.3.0.1
793      */
794     public int intValue() {
795         return (int) re;
796     }
797 
798     /** Returns real part.
799      * @since 2.3.0.1
800      */
801     public long longValue() {
802         return (long) re;
803     }
804 
805 }
806
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags