1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package edu.sc.seis.TauP;
24
25 /***
26 * This class implements complex numbers. It provides the basic operations
27 * (addition, subtraction, multiplication, division) as well as a set of
28 * complex functions.
29 *
30 * The binary operations have the form, where op is <code>plus</code>,
31 * <code>minus</code>, <code>times</code> or <code>over</code>.
32 * <pre>
33 * public static Complex op(Complex x, Complex y) // x op y
34 * public static Complex op(Complex x, double y) // x op y
35 * public static Complex op(double x, Complex y) // x op y
36 * public Complex op(Complex y) // this op y
37 * public Complex op(double y) // this op y
38 * public Complex opReverse(double x) // x op this
39 * public Complex opEquals(Complex y) // this op= y
40 * public Complex opEquals(double y) // this op= y
41 * </pre>
42 *
43 * The functions in this class follow the rules for complex arithmetic
44 * as defined C9x Annex G:"IEC 559-compatible complex arithmetic."
45 * The API is not the same, but handling of infinities, NaNs, and positive
46 * and negative zeros is intended to follow the same rules.
47 *
48 * This class depends on the standard java.lang.Math class following
49 * certain rules, as defined in the C9x Annex F, for the handling of
50 * infinities, NaNs, and positive and negative zeros. Sun's specification
51 * is that java.lang.Math should reproduce the results in the Sun's fdlibm
52 * C library. This library appears to follow the Annex F specification.
53 * At least on Windows, Sun's JDK 1.0 and 1.1 do NOT follow this specification.
54 * Sun's JDK 1.2(RC2) does follow the Annex F specification. Thesefore,
55 * this class will not give the expected results for edge cases with
56 * JDK 1.0 and 1.1.
57 */
58 public class Complex implements java.io.Serializable, Cloneable
59 {
60 /***
61 * @serial Real part of the Complex.
62 */
63 public double re;
64
65 /***
66 * @serial Imaginary part of the Complex.
67 */
68 public double im;
69
70 /***
71 * Serialization ID
72 */
73 static final long serialVersionUID = -633126172485117692L;
74
75
76 /***
77 * String used in converting Complex to String.
78 * Default is "i", but sometimes "j" is desired.
79 * Note that this is set for the class, not for
80 * a particular instance of a Complex.
81 */
82 public static String suffix = "i";
83
84
85 private final static long negZeroBits =
86 Double.doubleToLongBits(1.0/Double.NEGATIVE_INFINITY);
87
88
89
90 /***
91 * Constructs a Complex equal to the argument.
92 * @param z A Complex object
93 * If z is null then a NullPointerException is thrown.
94 */
95 public Complex(Complex z)
96 {
97 re = z.re;
98 im = z.im;
99 }
100
101
102 /***
103 * Constructs a Complex with real and imaginary parts given
104 * by the input arguments.
105 * @param re A double value equal to the real part of the Complex object.
106 * @param im A double value equal to the imaginary part of the Complex object.
107 */
108 public Complex(double re, double im)
109 {
110 this.re = re;
111 this.im = im;
112 }
113
114
115 /***
116 * Constructs a Complex with a zero imaginary part.
117 * @param re A double value equal to the real part of the Complex object.
118 */
119 public Complex(double re)
120 {
121 this.re = re;
122 this.im = 0.0;
123 }
124
125
126 /***
127 * Constructs a Complex equal to zero.
128 */
129 public Complex()
130 {
131 re = 0.0;
132 im = 0.0;
133 }
134
135
136 /***
137 * Tests if this is a complex Not-a-Number (NaN) value.
138 * @return True if either component of the Complex object is NaN;
139 * false, otherwise.
140 */
141 public boolean isNaN()
142 {
143 return (Double.isNaN(re) || Double.isNaN(im));
144 }
145
146
147 /***
148 * Tests if this is an infinite complex number.
149 * @return True if either component of the Complex object is NaN; false, otherwise.
150 */
151 public boolean isInfinite()
152 {
153 return Double.isInfinite(re) || Double.isInfinite(im);
154 }
155
156
157 /***
158 * Compares with another Complex.
159 * <p><em>Note: To be useful in hashtables this method
160 * considers two NaN double values to be equal. This
161 * is not according to IEEE specification.</em>
162 * @param z A Complex object.
163 * @return True if the real and imaginary parts of this object
164 * are equal to their counterparts in the argument; false, otherwise.
165 */
166 public boolean equals(Complex z)
167 {
168 if (isNaN() && z.isNaN()) {
169 return true;
170 } else {
171 return (re == z.re && im == z.im);
172 }
173 }
174
175
176 /***
177 * Compares this object against the specified object.
178 * <p><em>Note: To be useful in hashtables this method
179 * considers two NaN double values to be equal. This
180 * is not according to IEEE specification</em>
181 * @param obj The object to compare with.
182 * @return True if the objects are the same; false otherwise.
183 */
184 public boolean equals(Object obj)
185 {
186 if (obj == null) {
187 return false;
188 } else if (obj instanceof Complex) {
189 return equals((Complex)obj);
190 } else {
191 return false;
192 }
193 }
194
195 /***
196 * Returns a hashcode for this Complex.
197 * @return A hash code value for this object.
198 */
199 public int hashCode()
200 {
201 long re_bits = Double.doubleToLongBits(re);
202 long im_bits = Double.doubleToLongBits(im);
203 return (int)((re_bits^im_bits)^((re_bits^im_bits)>>32));
204 }
205
206 /***
207 * Returns the real part of a Complex object.
208 * @param z A Complex object.
209 * @return The real part of z.
210 */
211 public static double real(Complex z)
212 {
213 return z.re;
214 }
215
216 /***
217 * Returns the imaginary part of a Complex object.
218 * @param z A Complex object.
219 * @return The imaginary part of z.
220 */
221 public static double imag(Complex z)
222 {
223 return z.im;
224 }
225
226
227 /***
228 * Copies the contents of a double into this Complex.
229 * @param x The real part.
230 * @return The modified Complex number. The real part is
231 * set to x and the imaginary part is set to zero.
232 */
233 public Complex assign(double x)
234 {
235 re = x;
236 im = 0.0;
237 return this;
238 }
239
240 /***
241 * Copies the contents of a Complex into this Complex.
242 * @param z A Complex object.
243 * @return The modified Complex number.
244 * If the argument is null then a NullPointerException is thrown.
245 */
246 public Complex assign(Complex z)
247 {
248 re = z.re;
249 im = z.im;
250 return this;
251 }
252
253 /***
254 * Returns the negative of a Complex object, -z.
255 * @param z A Complex object.
256 * @return A newly constructed Complex initialized to
257 * the negative of the argument.
258 */
259 public static Complex negative(Complex z)
260 {
261 return new Complex(-z.re, -z.im);
262 }
263
264
265 /***
266 * Returns the complex conjugate of a Complex object.
267 * @param z A Complex object.
268 * @return A newly constructed Complex initialized to complex conjugate of z.
269 */
270 public static Complex conjugate(Complex z)
271 {
272 return new Complex(z.re, -z.im);
273 }
274
275
276 /***
277 * Returns the sum of two Complex objects, x+y.
278 * @param x A Complex object.
279 * @param y A Complex object.
280 * @return A newly constructed Complex initialized to x+y.
281 */
282 public static Complex plus(Complex x, Complex y)
283 {
284 return new Complex(x.re+y.re, x.im+y.im);
285 }
286
287 /***
288 * Returns the sum of a Complex and a double, x+y.
289 * @param x A Complex object.
290 * @param y A double value.
291 * @return A newly constructed Complex initialized to x+y.
292 */
293 public static Complex plus(Complex x, double y)
294 {
295 return new Complex(x.re+y, x.im);
296 }
297
298 /***
299 * Returns the sum of a double and a Complex, x+y.
300 * @param x A double value.
301 * @param y A Complex object.
302 * @return A newly constructed Complex initialized to x+y.
303 */
304 public static Complex plus(double x, Complex y)
305 {
306 return new Complex(x+y.re, y.im);
307 }
308
309 /***
310 * Returns the sum of this Complex and another Complex, this+y.
311 * @param y A Complex object.
312 * @return A newly constructed Complex initialized to this+y.
313 */
314 public Complex plus(Complex y)
315 {
316 return new Complex(re+y.re, im+y.im);
317 }
318
319 /***
320 * Returns the sum of this Complex a double, this+y.
321 * @param y A double value.
322 * @return A newly constructed Complex initialized to this+y.
323 */
324 public Complex plus(double y)
325 {
326 return new Complex(re+y, im);
327 }
328
329 /***
330 * Returns the sum of this Complex and a double, x+this.
331 * @param x A double value.
332 * @return A newly constructed Complex initialized to x+this.
333 */
334 public Complex plusReverse(double x)
335 {
336 return new Complex(re+x, im);
337 }
338
339 /***
340 * Adds a Complex to this Complex and returns the sum, this += y.
341 * @param y A Complex object.
342 * @return This object augmented by the argument.
343 */
344 public Complex plusEquals(Complex y)
345 {
346 re += y.re;
347 im += y.im;
348 return this;
349 }
350
351 /***
352 * Adds a double into this Complex and returns the sum, this += y.
353 * @param y A double value.
354 * @return This object augmented by the argument.
355 */
356 public Complex plusEquals(double y)
357 {
358 re += y;
359 return this;
360 }
361
362
363 /***
364 * Returns the difference of two Complex objects, x-y.
365 * @param x A Complex object.
366 * @param y A Complex object.
367 * @return A newly constructed Complex initialized to x-y.
368 */
369 public static Complex minus(Complex x, Complex y)
370 {
371 return new Complex(x.re-y.re, x.im-y.im);
372 }
373
374 /***
375 * Returns the difference of a Complex object and a double, x-y.
376 * @param x A Complex object.
377 * @param y A double value.
378 * @return A newly constructed Complex initialized to x-y.
379 */
380 public static Complex minus(Complex x, double y)
381 {
382 return new Complex(x.re-y, x.im);
383 }
384
385 /***
386 * Returns the difference of a double and a Complex object, x-y.
387 * @param x A double value.
388 * @param y A Complex object.
389 * @return A newly constructed Complex initialized to x-y..
390 */
391 public static Complex minus(double x, Complex y)
392 {
393 return new Complex(x-y.re, -y.im);
394 }
395
396 /***
397 * Returns the difference of this Complex object and
398 * another Complex object, this-y.
399 * @param y A Complex object.
400 * @return A newly constructed Complex initialized to this-y.
401 */
402 public Complex minus(Complex y)
403 {
404 return new Complex(re-y.re, im-y.im);
405 }
406
407 /***
408 * Subtracts a double from this Complex and returns the difference, this-y.
409 * @param y A double value.
410 * @return A newly constructed Complex initialized to this-y.
411 */
412 public Complex minus(double y)
413 {
414 return new Complex(re-y, im);
415 }
416
417 /***
418 * Returns the difference of this Complex object and a double, this-y.
419 * @param y A double value.
420 * @return A newly constructed Complex initialized to x-this.
421 */
422 public Complex minusReverse(double x)
423 {
424 return new Complex(x-re, -im);
425 }
426
427 /***
428 * Subtracts a Complex from this Complex and returns the difference, this -= y.
429 * @param y A Complex object.
430 * @return This object less the input argument.
431 */
432 public Complex minusEquals(Complex y)
433 {
434 re -= y.re;
435 im -= y.im;
436 return this;
437 }
438
439 /***
440 * Subtracts a double from this Complex and returns the difference, this -= y.
441 * @param y A double value.
442 * @return This object less the input argument.
443 */
444 public Complex minusEquals(double y)
445 {
446 re -= y;
447 return this;
448 }
449
450
451 /***
452 * Returns the product of two Complex objects, x*y.
453 * @param x A Complex object.
454 * @param y A Complex object.
455 * @return A newly constructed Complex initialized to x*y.
456 */
457 public static Complex times(Complex x, Complex y)
458 {
459 Complex t = new Complex(x.re*y.re-x.im*y.im, x.re*y.im+x.im*y.re);
460 if (Double.isNaN(t.re) && Double.isNaN(t.im))
461 timesNaN(x, y, t);
462 return t;
463 }
464
465
466
467
468 private static double copysign(double a, double b)
469 {
470 double abs = Math.abs(a);
471 return ((b < 0) ? -abs : abs);
472 }
473
474 /***
475 * Recovers infinities when computed x*y = NaN+i*NaN.
476 * This code is not part of times(), so that times
477 * could be inlined by an optimizing compiler.
478 * <p>
479 * This algorithm is adapted from the C9x Annex G:
480 * "IEC 559-compatible complex arithmetic."
481 * @param x First Complex operand.
482 * @param y Second Complex operand.
483 * @param t The product x*y, computed without regard to NaN.
484 * The real and/or the imaginary part of t is
485 * expected to be NaN.
486 * @return The corrected product of x*y.
487 */
488 private static void timesNaN(Complex x, Complex y, Complex t)
489 {
490 boolean recalc = false;
491 double a = x.re;
492 double b = x.im;
493 double c = y.re;
494 double d = y.im;
495
496 if (Double.isInfinite(a) || Double.isInfinite(b)) {
497
498 a = copysign(Double.isInfinite(a)?1.0:0.0, a);
499 b = copysign(Double.isInfinite(b)?1.0:0.0, b);
500 if (Double.isNaN(c)) c = copysign(0.0, c);
501 if (Double.isNaN(d)) d = copysign(0.0, d);
502 recalc = true;
503 }
504
505 if (Double.isInfinite(c) || Double.isInfinite(d)) {
506
507 a = copysign(Double.isInfinite(c)?1.0:0.0, c);
508 b = copysign(Double.isInfinite(d)?1.0:0.0, d);
509 if (Double.isNaN(a)) a = copysign(0.0, a);
510 if (Double.isNaN(b)) b = copysign(0.0, b);
511 recalc = true;
512 }
513
514 if (!recalc) {
515 if (Double.isInfinite(a*c) || Double.isInfinite(b*d) ||
516 Double.isInfinite(a*d) || Double.isInfinite(b*c)) {
517
518 if (Double.isNaN(a)) a = copysign(0.0, a);
519 if (Double.isNaN(b)) b = copysign(0.0, b);
520 if (Double.isNaN(c)) c = copysign(0.0, c);
521 if (Double.isNaN(d)) d = copysign(0.0, d);
522 recalc = true;
523 }
524 }
525
526 if (recalc) {
527 t.re = Double.POSITIVE_INFINITY * (a*c - b*d);
528 t.im = Double.POSITIVE_INFINITY * (a*d + b*c);
529 }
530 }
531
532
533 /***
534 * Returns the product of a Complex object and a double, x*y.
535 * @param x A Complex object.
536 * @param y A double value.
537 * @return A newly constructed Complex initialized to x*y.
538 */
539 public static Complex times(Complex x, double y)
540 {
541 return new Complex(x.re*y, x.im*y);
542 }
543
544 /***
545 * Returns the product of a double and a Complex object, x*y.
546 * @param x A double value.
547 * @param y A Complex object.
548 * @return A newly constructed Complex initialized to x*y.
549 */
550 public static Complex times(double x, Complex y)
551 {
552 return new Complex(x*y.re, x*y.im);
553 }
554
555 /***
556 * Returns the product of this Complex object and another Complex object, this*y.
557 * @param y A Complex object.
558 * @return A newly constructed Complex initialized to this*y.
559 */
560 public Complex times(Complex y)
561 {
562 return times(this,y);
563 }
564
565 /***
566 * Returns the product of this Complex object and a double, this*y.
567 * @param y A double value.
568 * @return A newly constructed Complex initialized to this*y.
569 */
570 public Complex times(double y)
571 {
572 return new Complex(re*y, im*y);
573 }
574
575 /***
576 * Returns the product of a double and this Complex, x*this.
577 * @param y A double value.
578 * @return A newly constructed Complex initialized to x*this.
579 */
580 public Complex timesReverse(double x)
581 {
582 return new Complex(x*re, x*im);
583 }
584
585 /***
586 * Multiplies this Complex object by another Complex and returns the product, this *= y.
587 * @param y A Complex object.
588 * @return This object multiplied by the input argument.
589 */
590 public Complex timesEquals(Complex y)
591 {
592 double new_re = re*y.re - im*y.im;
593 im = re*y.im + im*y.re;
594 re = new_re;
595 return this;
596 }
597
598 /***
599 * Multiplies this Complex by a double and returns the product, this *= y.
600 * @param y A double value.
601 * @return This object multiplied by the input argument.
602 */
603 public Complex timesEquals(double y)
604 {
605 re *= y;
606 im *= y;
607 return this;
608 }
609
610 private static boolean isFinite(double x)
611 {
612 return !(Double.isInfinite(x) || Double.isNaN(x));
613 }
614
615
616 /***
617 * Returns Complex object divided by a Complex object, x/y.
618 * @param x The numerator, a Complex object.
619 * @param y The denominator, a Complex object.
620 * @return A newly constructed Complex initialized to x/y.
621 */
622 public static Complex over(Complex x, Complex y)
623 {
624 double a = x.re;
625 double b = x.im;
626 double c = y.re;
627 double d = y.im;
628
629 double scale = Math.max(Math.abs(c), Math.abs(d));
630 boolean isScaleFinite = isFinite(scale);
631 if (isScaleFinite) {
632 c /= scale;
633 d /= scale;
634 }
635
636 double den = c*c + d*d;
637 Complex z = new Complex((a*c+b*d)/den, (b*c-a*d)/den);
638
639 if (isScaleFinite) {
640 z.re /= scale;
641 z.im /= scale;
642 }
643
644
645 if (Double.isNaN(z.re) && Double.isNaN(z.im)) {
646 if (den == 0.0 && (!Double.isNaN(a) || !Double.isNaN(b))) {
647 double s = copysign(Double.POSITIVE_INFINITY, c);
648 z.re = s * a;
649 z.im = s * b;
650
651 } else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
652 isFinite(c) && isFinite(d)) {
653 a = copysign(Double.isInfinite(a)?1.0:0.0, a);
654 b = copysign(Double.isInfinite(b)?1.0:0.0, b);
655 z.re = Double.POSITIVE_INFINITY * (a*c + b*d);
656 z.im = Double.POSITIVE_INFINITY * (b*c - a*d);
657
658 } else if (Double.isInfinite(scale) &&
659 isFinite(a) && isFinite(b)) {
660 c = copysign(Double.isInfinite(c)?1.0:0.0, c);
661 d = copysign(Double.isInfinite(d)?1.0:0.0, d);
662 z.re = 0.0 * (a*c + b*d);
663 z.im = 0.0 * (b*c - a*d);
664 }
665 }
666 return z;
667 }
668
669
670 /***
671 * Returns Complex object divided by a double, x/y.
672 * @param x The numerator, a Complex object.
673 * @param y The denominator, a double.
674 * @return A newly constructed Complex initialized to x/y.
675 */
676 public static Complex over(Complex x, double y)
677 {
678 return new Complex(x.re/y, x.im/y);
679 }
680
681 /***
682 * Returns a double divided by a Complex object, x/y.
683 * @param x A double value.
684 * @param y The denominator, a Complex object.
685 * @return A newly constructed Complex initialized to x/y.
686 */
687 public static Complex over(double x, Complex y)
688 {
689 return y.overReverse(x);
690 }
691
692 /***
693 * Returns this Complex object divided by another Complex object, this/y.
694 * @param y The denominator, a Complex object.
695 * @return A newly constructed Complex initialized to x/y.
696 */
697 public Complex over(Complex y)
698 {
699 return over(this, y);
700 }
701
702 /***
703 * Returns this Complex object divided by double, this/y.
704 * @param y The denominator, a double.
705 * @return A newly constructed Complex initialized to x/y.
706 */
707 public Complex over(double y)
708 {
709 return over(this, y);
710 }
711
712 /***
713 * Returns a double dividied by this Complex object, x/this.
714 * @param x The numerator, a double.
715 * @return A newly constructed Complex initialized to x/this.
716 */
717 public Complex overReverse(double x)
718 {
719 double den, t;
720 Complex z;
721 if (Math.abs(re) > Math.abs(im)) {
722 t = im / re;
723 den = re + im*t;
724 z = new Complex(x/den, -x*t/den);
725 } else {
726 t = re / im;
727 den = im + re*t;
728 z = new Complex(x*t/den, -x/den);
729 }
730 return z;
731 }
732
733 /***
734 * Divides this Complex by a Complex and returns the result, this /= y.
735 * @param y The denominator, a Complex object.
736 * @return This object divided by the input argument.
737 */
738 public Complex overEquals(Complex y)
739 {
740 assign(over(this,y));
741 return this;
742 }
743
744 /***
745 * Divides this Complex by a double and returns the result, this /= y.
746 * @param y The denominator, a double.
747 * @return This object divided by the input argument.
748 */
749 public Complex overEquals(double y)
750 {
751 re /= y;
752 im /= y;
753 return this;
754 }
755
756
757 /***
758 * Returns the absolute value (modulus) of a Complex, |z|.
759 * @param z A Complex object.
760 * @return A double value equal to the absolute value of the argument.
761 */
762 public static double abs(Complex z)
763 {
764 if (z.isInfinite()) return Double.POSITIVE_INFINITY;
765
766 double x = Math.abs(z.re);
767 double y = Math.abs(z.im);
768 if (x + y == 0.0) {
769 return 0.0;
770 } else if (x > y) {
771 y /= x;
772 return x*Math.sqrt(1.0+y*y);
773 } else {
774 x /= y;
775 return y*Math.sqrt(x*x+1.0);
776 }
777 }
778
779
780 /***
781 * Returns the argument (phase) of a Complex, in radians,
782 * with a branch cut along the negative real axis.
783 * @param z A Complex object.
784 * @return A double value equal to the argument (or phase) of a Complex.
785 * It is in the interval [-pi,pi].
786 */
787 public static double argument(Complex z)
788 {
789 return Math.atan2(z.im, z.re);
790 }
791
792
793 /***
794 * Returns the square root of a Complex,
795 * with a branch cut along the negative real axis.
796 * @param z A Complex object.
797 * @return A newly constructed Complex initialized
798 * to square root of z. Its real part is
799 * non-negative.
800 */
801 public static Complex sqrt(Complex z)
802 {
803 Complex result = new Complex();
804
805 if (Double.isInfinite(z.im)) {
806 result.re = Double.POSITIVE_INFINITY;
807 result.im = z.im;
808 } else if (Double.isNaN(z.re)) {
809 result.re = result.im = Double.NaN;
810 } else if (Double.isNaN(z.im)) {
811 if (Double.isInfinite(z.re)) {
812 if (z.re > 0) {
813 result.re = z.re;
814 result.im = z.im;
815 } else {
816 result.re = z.im;
817 result.im = Double.POSITIVE_INFINITY;
818 }
819 } else {
820 result.re = result.im = Double.NaN;
821 }
822 } else {
823
824
825 double t = abs(z);
826
827 if (Math.abs(z.re) <= Math.abs(z.im)) {
828
829 result.re = Math.sqrt(0.5*(t+z.re));
830 result.im = Math.sqrt(0.5*(t-z.re));
831 } else {
832
833 if (z.re > 0) {
834 result.re = t + z.re;
835 result.im = Math.abs(z.im)*Math.sqrt(0.5/result.re);
836 result.re = Math.sqrt(0.5*result.re);
837 } else {
838 result.im = t - z.re;
839 result.re = Math.abs(z.im)*Math.sqrt(0.5/result.im);
840 result.im = Math.sqrt(0.5*result.im);
841 }
842 }
843 if (z.im < 0)
844 result.im = -result.im;
845 }
846 return result;
847 }
848
849
850 /***
851 * Returns the exponential of a Complex z, exp(z).
852 * @param z A Complex object.
853 * @return A newly constructed Complex initialized to exponential
854 * of the argument.
855 */
856 public static Complex exp(Complex z)
857 {
858 Complex result = new Complex();
859
860 double r = Math.exp(z.re);
861
862 double cosa = Math.cos(z.im);
863 double sina = Math.sin(z.im);
864 if (Double.isInfinite(z.im) || Double.isNaN(z.im) || Math.abs(cosa)>1) {
865 cosa = sina = Double.NaN;
866
867 }
868
869 if (Double.isInfinite(z.re) || Double.isInfinite(r)) {
870 if (z.re < 0) {
871 r = 0;
872 if (Double.isInfinite(z.im) || Double.isNaN(z.im)) {
873 cosa = sina = 0;
874 } else {
875 cosa /= Double.POSITIVE_INFINITY;
876 sina /= Double.POSITIVE_INFINITY;
877 }
878 } else {
879 r = z.re;
880 if (Double.isNaN(z.im)) cosa = 1;
881 }
882 }
883
884 if (z.im == 0.0) {
885 result.re = r;
886 result.im = z.im;
887 } else {
888 result.re = r*cosa;
889 result.im = r*sina;
890 }
891 return result;
892 }
893
894
895 /***
896 * Returns the logarithm of a Complex z,
897 * with a branch cut along the negative real axis.
898 * @param z A Complex object.
899 * @return A newly constructed Complex initialized to logarithm
900 * of the argument. Its imaginary part is in the
901 * interval [-i*pi,i*pi].
902 */
903 public static Complex log(Complex z)
904 {
905 Complex result = new Complex();
906
907 if (Double.isNaN(z.re)) {
908 result.re = result.im = z.re;
909 if (Double.isInfinite(z.im))
910 result.re = Double.POSITIVE_INFINITY;
911 } else if (Double.isNaN(z.im)) {
912 result.re = result.im = z.im;
913 if (Double.isInfinite(z.re))
914 result.re = Double.POSITIVE_INFINITY;
915 } else {
916 result.re = Math.log(abs(z));
917 result.im = argument(z);
918 }
919 return result;
920 }
921
922
923 /***
924 * Returns the sine of a Complex.
925 * @param z A Complex object.
926 * @return A newly constructed Complex initialized to sine of the argument.
927 */
928 public static Complex sin(Complex z)
929 {
930
931 Complex iz = new Complex(-z.im,z.re);
932 Complex s = sinh(iz);
933 double re = s.im;
934 s.im = -s.re;
935 s.re = re;
936 return s;
937 }
938
939
940 /***
941 * Returns the cosine of a Complex.
942 * @param z A Complex object.
943 * @return A newly constructed Complex initialized to cosine of the argument.
944 */
945 public static Complex cos(Complex z)
946 {
947
948 return cosh(new Complex(-z.im,z.re));
949 }
950
951 /***
952 * Returns the tangent of a Complex.
953 * @param z A Complex object.
954 * @return A newly constructed Complex initialized
955 * to tangent of the argument.
956 */
957 public static Complex tan(Complex z)
958 {
959
960 Complex iz = new Complex(-z.im,z.re);
961 Complex s = tanh(iz);
962 double re = s.im;
963 s.im = -s.re;
964 s.re = re;
965 return s;
966 }
967
968 /***
969 * Returns the inverse sine (arc sine) of a Complex,
970 * with branch cuts outside the interval [-1,1] along the
971 * real axis.
972 * @param z A Complex object.
973 * @return A newly constructed Complex initialized to inverse
974 * (arc) sine of the argument. The real part of the
975 * result is in the interval [-pi/2,+pi/2].
976 */
977 public static Complex asin(Complex z)
978 {
979 Complex result = new Complex();
980
981 double r = abs(z);
982
983 if (Double.isInfinite(r)) {
984 boolean infiniteX = Double.isInfinite(z.re);
985 boolean infiniteY = Double.isInfinite(z.im);
986 if (infiniteX) {
987 double pi2 = 0.5*Math.PI;
988 result.re = (z.re>0 ? pi2 : -pi2);
989 if (infiniteY) result.re /= 2;
990 } else if (infiniteY) {
991 result.re = z.re/Double.POSITIVE_INFINITY;
992 }
993 if (Double.isNaN(z.im)) {
994 result.im = -z.re;
995 result.re = z.im;
996 } else {
997 result.im = z.im*Double.POSITIVE_INFINITY;
998 }
999 return result;
1000 } else if (Double.isNaN(r)) {
1001 result.re = result.im = Double.NaN;
1002 if (z.re == 0) result.re = z.re;
1003 } else if (r < 2.58095e-08) {
1004
1005 result.re = z.re;
1006 result.im = z.im;
1007 } else if (z.re == 0) {
1008 result.re = 0;
1009 result.im = Sfun.asinh(z.im);
1010 } else if (r <= 0.1) {
1011 Complex z2 = times(z,z);
1012
1013 for (int i = 1; i <= 8; i++) {
1014 double twoi = 2*(8-i) + 1;
1015 result = times(times(result,z2),twoi/(twoi+1.0));
1016 result.re += 1.0/twoi;
1017 }
1018 result.timesEquals(z);
1019 } else {
1020
1021
1022
1023
1024 Complex w = ((z.im < 0) ? negative(z) : z);
1025 Complex sqzp1 = sqrt(plus(w,1.0));
1026 if (sqzp1.im < 0.0)
1027 sqzp1 = negative(sqzp1);
1028 Complex sqzm1 = sqrt(minus(w,1.0));
1029 result = log(plus(w,times(sqzp1,sqzm1)));
1030
1031 double rx = result.re;
1032 result.re = 0.5*Math.PI + result.im;
1033 result.im = -rx;
1034 }
1035
1036 if (result.re > 0.5*Math.PI) {
1037 result.re = Math.PI - result.re;
1038 result.im = -result.im;
1039 }
1040 if (result.re < -0.5*Math.PI) {
1041 result.re = -Math.PI - result.re;
1042 result.im = -result.im;
1043 }
1044 if (z.im < 0) {
1045 result.re = -result.re;
1046 result.im = -result.im;
1047 }
1048 return result;
1049 }
1050
1051
1052 /***
1053 * Returns the inverse cosine (arc cosine) of a Complex,
1054 * with branch cuts outside the interval [-1,1] along the
1055 * real axis.
1056 * @param z A Complex object.
1057 * @return A newly constructed Complex initialized to
1058 * inverse (arc) cosine of the argument.
1059 * The real part of the result is in the interval [0,pi].
1060 */
1061 public static Complex acos(Complex z)
1062 {
1063 Complex result = new Complex();
1064 double r = abs(z);
1065
1066 if (Double.isInfinite(z.re) && Double.isNaN(z.im)) {
1067 result.re = Double.NaN;
1068 result.im = Double.NEGATIVE_INFINITY;
1069 } else if (Double.isInfinite(r)) {
1070 result.re = Math.atan2(Math.abs(z.im),z.re);
1071 result.im = z.im*Double.NEGATIVE_INFINITY;
1072 } else if (r == 0) {
1073 result.re = Math.PI/2;
1074 result.im = -z.im;
1075 } else {
1076 result.assign(minus(Math.PI/2,asin(z)));
1077 }
1078 return result;
1079 }
1080
1081 /***
1082 * Returns the inverse tangent (arc tangent) of a Complex,
1083 * with branch cuts outside the interval [-i,i] along the
1084 * imaginary axis.
1085 * @param z A Complex object.
1086 * @return A newly constructed Complex initialized to
1087 * inverse (arc) tangent of the argument.
1088 * Its real part is in the interval [-pi/2,pi/2].
1089 */
1090 public static Complex atan(Complex z)
1091 {
1092 Complex result = new Complex();
1093 double r = abs(z);
1094
1095 if (Double.isInfinite(r)) {
1096 double pi2 = 0.5*Math.PI;
1097 double im = (Double.isNaN(z.im) ? 0 : z.im);
1098 result.re = (z.re<0 ? -pi2 : pi2);
1099 result.im = (im<0 ? -1 : 1)/Double.POSITIVE_INFINITY;
1100 if (Double.isNaN(z.re)) result.re = z.re;
1101 } else if (Double.isNaN(r)) {
1102 result.re = result.im = Double.NaN;
1103 if (z.im == 0) result.im = z.im;
1104 } else if (r < 1.82501e-08) {
1105
1106 result.re = z.re;
1107 result.im = z.im;
1108 } else if (r < 0.1) {
1109 Complex z2 = times(z,z);
1110
1111 for (int k = 0; k < 17; k++) {
1112 Complex temp = times(z2,result);
1113 int twoi = 2*(17-k) - 1;
1114 result.re = 1.0/twoi - temp.re;
1115 result.im = -temp.im;
1116 }
1117 result.timesEquals(z);
1118 } else if (r < 9.0072e+15) {
1119
1120 double r2 = r*r;
1121 result.re = 0.5*Math.atan2(2*z.re,1.0-r2);
1122 result.im = 0.25*Math.log((r2+2*z.im+1)/(r2-2*z.im+1));
1123 } else {
1124 result.re = ((z.re < 0.0) ? -0.5*Math.PI : 0.5*Math.PI);
1125 }
1126 return result;
1127 }
1128
1129 /***
1130 * Returns the hyperbolic sine of a Complex.
1131 * @param z A Complex object.
1132 * @return A newly constructed Complex initialized to hyperbolic
1133 * sine of the argument.
1134 */
1135 public static Complex sinh(Complex z)
1136 {
1137 double coshx = Sfun.cosh(z.re);
1138 double sinhx = Sfun.sinh(z.re);
1139 double cosy = Math.cos(z.im);
1140 double siny = Math.sin(z.im);
1141 boolean infiniteX = Double.isInfinite(coshx);
1142 boolean infiniteY = Double.isInfinite(z.im);
1143 Complex result;
1144
1145 if (z.im == 0) {
1146 result = new Complex(Sfun.sinh(z.re));
1147 } else {
1148
1149 result = new Complex(sinhx*cosy, coshx*siny);
1150 if (infiniteY) {
1151 result.im = Double.NaN;
1152 if (z.re == 0) result.re = 0;
1153 }
1154 if (infiniteX) {
1155 result.re = z.re*cosy;
1156 result.im = z.re*siny;
1157 if (z.im == 0) result.im = 0;
1158 if (infiniteY) result.re = z.im;
1159 }
1160 }
1161 return result;
1162 }
1163
1164 /***
1165 * Returns the hyperbolic cosh of a Complex.
1166 * @param z A Complex object.
1167 * @return A newly constructed Complex initialized to
1168 * the hyperbolic cosine of the argument.
1169 */
1170 public static Complex cosh(Complex z)
1171 {
1172 if (z.im == 0) {
1173 return new Complex(Sfun.cosh(z.re));
1174 }
1175
1176 double coshx = Sfun.cosh(z.re);
1177 double sinhx = Sfun.sinh(z.re);
1178 double cosy = Math.cos(z.im);
1179 double siny = Math.sin(z.im);
1180 boolean infiniteX = Double.isInfinite(coshx);
1181 boolean infiniteY = Double.isInfinite(z.im);
1182
1183
1184 Complex result = new Complex(coshx*cosy, sinhx*siny);
1185 if (infiniteY) result.re = Double.NaN;
1186 if (z.re == 0) {
1187 result.im = 0;
1188 } else if (infiniteX) {
1189 result.re = z.re*cosy;
1190 result.im = z.re*siny;
1191 if (z.im == 0) result.im = 0;
1192 if (Double.isNaN(z.im)) {
1193 result.re = z.re;
1194 } else if (infiniteY) {
1195 result.re = z.im;
1196 }
1197 }
1198 return result;
1199 }
1200
1201 /***
1202 * Returns the hyperbolic tanh of a Complex.
1203 * @param z A Complex object.
1204 * @return A newly constructed Complex initialized to
1205 * the hyperbolic tangent of the argument.
1206 */
1207 public static Complex tanh(Complex z)
1208 {
1209 double sinh2x = Sfun.sinh(2*z.re);
1210
1211 if (z.im == 0) {
1212 return new Complex(Sfun.tanh(z.re));
1213 } else if (sinh2x == 0) {
1214 return new Complex(0,Math.tan(z.im));
1215 }
1216
1217 double cosh2x = Sfun.cosh(2*z.re);
1218 double cos2y = Math.cos(2*z.im);
1219 double sin2y = Math.sin(2*z.im);
1220 boolean infiniteX = Double.isInfinite(cosh2x);
1221
1222
1223 if (Double.isInfinite(z.im) || Double.isNaN(z.im)) {
1224 cos2y = sin2y = Double.NaN;
1225 }
1226
1227 if (infiniteX)
1228 return new Complex(z.re > 0 ? 1 : -1);
1229
1230
1231 double den = (cosh2x + cos2y);
1232 return new Complex(sinh2x/den, sin2y/den);
1233 }
1234
1235 /***
1236 * Returns the Complex z raised to the x power,
1237 * with a branch cut for the first parameter (z) along the
1238 * negative real axis.
1239 * @param z A Complex object.
1240 * @param x A double value.
1241 * @return A newly constructed Complex initialized to z to the power x.
1242 */
1243 public static Complex pow(Complex z, double x)
1244 {
1245 double absz = abs(z);
1246 Complex result = new Complex();
1247
1248 if (absz == 0.0) {
1249 result.assign(z);
1250 } else {
1251 double a = argument(z);
1252 double e = Math.pow(absz, x);
1253 result.re = e*Math.cos(x*a);
1254 result.im = e*Math.sin(x*a);
1255 }
1256 return result;
1257 }
1258
1259 /***
1260 * Returns the inverse hyperbolic sine (arc sinh) of a Complex,
1261 * with a branch cuts outside the interval [-i,i].
1262 * @param z A Complex object.
1263 * @return A newly constructed Complex initialized to
1264 * inverse (arc) hyperbolic sine of the argument.
1265 * Its imaginary part is in the interval [-i*pi/2,i*pi/2].
1266 */
1267 public static Complex asinh(Complex z)
1268 {
1269
1270 Complex miz = new Complex(z.im,-z.re);
1271 Complex result = asin(miz);
1272 double rx = result.im;
1273 result.im = result.re;
1274 result.re = -rx;
1275 return result;
1276 }
1277
1278 /***
1279 * Returns the inverse hyperbolic cosine (arc cosh) of a Complex,
1280 * with a branch cut at values less than one along the real axis.
1281 * @param z A Complex object.
1282 * @return A newly constructed Complex initialized to
1283 * inverse (arc) hyperbolic cosine of the argument.
1284 * The real part of the result is non-negative and its
1285 * imaginary part is in the interval [-i*pi,i*pi].
1286 */
1287 public static Complex acosh(Complex z)
1288 {
1289 Complex result = acos(z);
1290 double rx = -result.im;
1291 result.im = result.re;
1292 result.re = rx;
1293 if (result.re < 0 || isNegZero(result.re)) {
1294 result.re = -result.re;
1295 result.im = -result.im;
1296 }
1297 return result;
1298 }
1299
1300
1301 /***
1302 * Returns true is x is a negative zero.
1303 */
1304 private static boolean isNegZero(double x)
1305 {
1306 return (Double.doubleToLongBits(x) == negZeroBits);
1307 }
1308
1309 /***
1310 * Returns the inverse hyperbolic tangent (arc tanh) of a Complex,
1311 * with a branch cuts outside the interval [-1,1] on the real axis.
1312 * @param z A Complex object.
1313 * @return A newly constructed Complex initialized to
1314 * inverse (arc) hyperbolic tangent of the argument.
1315 * The imaginary part of the result is in the interval
1316 * [-i*pi/2,i*pi/2].
1317 */
1318 public static Complex atanh(Complex z)
1319 {
1320
1321 Complex miz = new Complex(z.im,-z.re);
1322 Complex result = atan(miz);
1323 double rx = result.im;
1324 result.im = result.re;
1325 result.re = -rx;
1326 return result;
1327
1328 }
1329
1330
1331 /***
1332 * Returns the Complex x raised to the Complex y power.
1333 * @param x A Complex object.
1334 * @param y A Complex object.
1335 * @return A newly constructed Complex initialized
1336 * to x<SUP><FONT SIZE="1">y</FONT></SUP><FONT SIZE="3">.
1337 */
1338 public static Complex pow(Complex x, Complex y)
1339 {
1340 return exp(times(y,log(x)));
1341 }
1342
1343
1344
1345
1346 /***
1347 * Returns a String representation for the specified Complex.
1348 * @return A String representation for this object.
1349 */
1350 public String toString()
1351 {
1352 if (im == 0.0)
1353 return String.valueOf(re);
1354
1355 if (re == 0.0)
1356 return String.valueOf(im) + suffix;
1357
1358 String sign = (im < 0.0) ? "" : "+";
1359 return (String.valueOf(re) + sign + String.valueOf(im) + suffix);
1360 }
1361
1362
1363 /***
1364 * Parses a string into a Complex.
1365 * @param s The string to be parsed.
1366 * @return A newly constructed Complex initialized to the value represented
1367 * by the string argument.
1368 * @exception NumberFormatException If the string does not contain a parsable Complex number.
1369 * @exception NullPointerException If the input argument is null.
1370 */
1371 public static Complex valueOf(String s) throws NumberFormatException
1372 {
1373 String input = s.trim();
1374 int iBeginNumber = 0;
1375 Complex z = new Complex();
1376 int state = 0;
1377 int sign = 1;
1378 boolean haveRealPart = false;
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389 for (int k = 0; k < input.length(); k++) {
1390
1391 char ch = input.charAt(k);
1392
1393 switch (ch) {
1394
1395 case '0': case '1': case '2': case '3': case '4':
1396 case '5': case '6': case '7': case '8': case '9':
1397 if (state == 0 || state == 1) {
1398 state = 2;
1399 } else if (state == 4) {
1400 state = 5;
1401 }
1402 break;
1403
1404 case '-':
1405 case '+':
1406 sign = ((ch=='+') ? 1 : -1);
1407 if (state == 0) {
1408 state = 1;
1409 } else if (state == 4) {
1410 state = 5;
1411 } else {
1412 if (!haveRealPart) {
1413
1414 z.re = Double.valueOf(input.substring(iBeginNumber,k)).doubleValue();
1415 haveRealPart = true;
1416
1417 iBeginNumber = k;
1418 state = 1;
1419 } else {
1420 throw new NumberFormatException(input);
1421 }
1422 }
1423 break;
1424
1425 case '.':
1426 if (state == 0 || state == 1 || state == 2)
1427 state = 3;
1428 else
1429 throw new NumberFormatException(input);
1430 break;
1431
1432 case 'i': case 'I':
1433 case 'j': case 'J':
1434 if (k+1 != input.length()) {
1435 throw new NumberFormatException(input);
1436 } else if (state == 0 || state == 1) {
1437 z.im = sign;
1438 return z;
1439 } else if (state == 2 || state == 3 || state == 5) {
1440 z.im = Double.valueOf(input.substring(iBeginNumber,k)).doubleValue();
1441 return z;
1442 } else {
1443 throw new NumberFormatException(input);
1444 }
1445
1446
1447 case 'e': case 'E': case 'd': case 'D':
1448 if (state == 2 || state == 3) {
1449 state = 4;
1450 } else {
1451 throw new NumberFormatException(input);
1452 }
1453 break;
1454
1455 default:
1456 throw new NumberFormatException(input);
1457 }
1458
1459 }
1460
1461 if (!haveRealPart) {
1462 z.re = Double.valueOf(input).doubleValue();
1463 return z;
1464 } else {
1465 throw new NumberFormatException(input);
1466 }
1467 }
1468
1469 }