001: // Copyright (c) 1997 Per M.A. Bothner.
002: // This is free software; for terms and warranty disclaimer see ./COPYING.
003:
004: package gnu.math;
005:
006: import java.io.*;
007:
008: /** A complex number using rectangular (Cartesian) plain double values.
009: * @author Per Bothner
010: * @author Some algorithms were transcribed from GNU libstdc++,
011: * written by Jason Merrill.
012: * Also see below for copyrights for functions taken from fdlib and f2c.
013: */
014:
015: public class DComplex extends Complex implements Externalizable {
016: double real;
017: double imag;
018:
019: public DComplex() {
020: }
021:
022: public DComplex(double real, double imag) {
023: this .real = real;
024: this .imag = imag;
025: }
026:
027: public RealNum re() {
028: return new DFloNum(real);
029: }
030:
031: public double doubleValue() {
032: return real;
033: }
034:
035: public RealNum im() {
036: return new DFloNum(imag);
037: }
038:
039: public double doubleImagValue() {
040: return imag;
041: }
042:
043: public boolean equals(Object obj) {
044: if (obj == null || !(obj instanceof Complex))
045: return false;
046: Complex y = (Complex) obj;
047: return y.unit() == Unit.Empty
048: && (Double.doubleToLongBits(real) == Double
049: .doubleToLongBits(y.reValue()))
050: && (Double.doubleToLongBits(imag) == Double
051: .doubleToLongBits(y.imValue()));
052: }
053:
054: public String toString() {
055: String prefix = "";
056:
057: String reString;
058: if (real == 1.0 / 0.0) {
059: prefix = "#i";
060: reString = "1/0";
061: } else if (real == -1.0 / 0.0) {
062: prefix = "#i";
063: reString = "-1/0";
064: } else if (Double.isNaN(real)) {
065: prefix = "#i";
066: reString = "0/0";
067: } else
068: reString = Double.toString(real);
069:
070: if (Double.doubleToLongBits(imag) == 0) // i.e. imag is 0.0 and not -0.0
071: return prefix + reString;
072:
073: String imString;
074: if (imag == 1.0 / 0.0) {
075: prefix = "#i";
076: imString = "+1/0i";
077: } else if (imag == -1.0 / 0.0) {
078: prefix = "#i";
079: imString = "-1/0i";
080: } else if (Double.isNaN(imag)) {
081: prefix = "#i";
082: imString = "+0/0i";
083: } else {
084: imString = Double.toString(imag) + "i";
085: if (imString.charAt(0) != '-')
086: imString = "+" + imString;
087: }
088:
089: return ((Double.doubleToLongBits(real) == 0 ? prefix : prefix
090: + reString) + imString);
091: }
092:
093: public String toString(int radix) {
094: if (radix == 10)
095: return toString();
096: return "#d" + toString();
097: }
098:
099: // All transcendental complex functions return DComplex
100:
101: public final Numeric neg() {
102: return new DComplex(-real, -imag);
103: }
104:
105: public Numeric add(Object y, int k) {
106: if (y instanceof Complex) {
107: Complex yc = (Complex) y;
108: if (yc.dimensions() != Dimensions.Empty)
109: throw new ArithmeticException("units mis-match");
110: return new DComplex(real + k * yc.reValue(), imag + k
111: * yc.imValue());
112: }
113: return ((Numeric) y).addReversed(this , k);
114: }
115:
116: public Numeric mul(Object y) {
117: if (y instanceof Complex) {
118: Complex yc = (Complex) y;
119: if (yc.unit() == Unit.Empty) {
120: double y_re = yc.reValue();
121: double y_im = yc.imValue();
122: return new DComplex(real * y_re - imag * y_im, real
123: * y_im + imag * y_re);
124: }
125: return Complex.times(this , yc);
126: }
127: return ((Numeric) y).mulReversed(this );
128: }
129:
130: public Numeric div(Object y) {
131: if (y instanceof Complex) {
132: Complex yc = (Complex) y;
133: return div(real, imag, yc.doubleValue(), yc
134: .doubleImagValue());
135: }
136: return ((Numeric) y).divReversed(this );
137: }
138:
139: public static DComplex power(double x_re, double x_im, double y_re,
140: double y_im) {
141: double h;
142: /* #ifdef JAVA5 */
143: // h = Math.hypot(x_re, x_im);
144: /* #else */
145: h = DComplex.hypot(x_re, x_im);
146: /* #endif */
147: double logr = Math.log(h);
148: double t = Math.atan2(x_im, x_re);
149: double r = Math.exp(logr * y_re - y_im * t);
150: t = y_im * logr + y_re * t;
151: return Complex.polar(r, t);
152: }
153:
154: public static Complex log(double x_re, double x_im) {
155: double h;
156: /* #ifdef JAVA5 */
157: // h = Math.hypot(x_re, x_im);
158: /* #else */
159: h = DComplex.hypot(x_re, x_im);
160: /* #endif */
161: return make(Math.log(h), Math.atan2(x_im, x_re));
162: }
163:
164: // The code below is adapted from f2c's libF77, and is subject to this
165: // copyright:
166:
167: /****************************************************************
168: Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
169:
170: Permission to use, copy, modify, and distribute this software
171: and its documentation for any purpose and without fee is hereby
172: granted, provided that the above copyright notice appear in all
173: copies and that both that the copyright notice and this
174: permission notice and warranty disclaimer appear in supporting
175: documentation, and that the names of AT&T Bell Laboratories or
176: Bellcore or any of their entities not be used in advertising or
177: publicity pertaining to distribution of the software without
178: specific, written prior permission.
179:
180: AT&T and Bellcore disclaim all warranties with regard to this
181: software, including all implied warranties of merchantability
182: and fitness. In no event shall AT&T or Bellcore be liable for
183: any special, indirect or consequential damages or any damages
184: whatsoever resulting from loss of use, data or profits, whether
185: in an action of contract, negligence or other tortious action,
186: arising out of or in connection with the use or performance of
187: this software.
188: ****************************************************************/
189:
190: public static DComplex div(double x_re, double x_im, double y_re,
191: double y_im) {
192: double ar = Math.abs(y_re);
193: double ai = Math.abs(y_im);
194: double nr, ni;
195: double t, d;
196: if (ar <= ai) {
197: t = y_re / y_im;
198: d = y_im * (1 + t * t);
199: nr = x_re * t + x_im;
200: ni = x_im * t - x_re;
201: } else {
202: t = y_im / y_re;
203: d = y_re * (1 + t * t);
204: nr = x_re + x_im * t;
205: ni = x_im - x_re * t;
206: }
207: return new DComplex(nr / d, ni / d);
208: }
209:
210: public static Complex sqrt(double x_re, double x_im) {
211: /* #ifdef JAVA5 */
212: // double r = Math.hypot(x_re, x_im);
213: /* #else */
214: double r = DComplex.hypot(x_re, x_im);
215: /* #endif */
216: double nr, ni;
217: if (r == 0.0)
218: nr = ni = r;
219: else if (x_re > 0) {
220: nr = Math.sqrt(0.5 * (r + x_re));
221: ni = x_im / nr / 2;
222: } else {
223: ni = Math.sqrt(0.5 * (r - x_re));
224: if (x_im < 0)
225: ni = -ni;
226: nr = x_im / ni / 2;
227: }
228: return new DComplex(nr, ni);
229: }
230:
231: // Transcribed from:
232: // http://netlib.bell-labs.com/netlib/fdlibm/e_hypot.c.Z
233: /*
234: * ====================================================
235: * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
236: *
237: * Developed at SunSoft, a Sun Microsystems, Inc. business.
238: * Permission to use, copy, modify, and distribute this
239: * software is freely granted, provided that this notice
240: * is preserved.
241: * ====================================================
242: */
243: /* __ieee754_hypot(x,y)
244: *
245: * Method :
246: * If (assume round-to-nearest) z=x*x+y*y
247: * has error less than sqrt(2)/2 ulp, than
248: * sqrt(z) has error less than 1 ulp (exercise).
249: *
250: * So, compute sqrt(x*x+y*y) with some care as
251: * follows to get the error below 1 ulp:
252: *
253: * Assume x>y>0;
254: * (if possible, set rounding to round-to-nearest)
255: * 1. if x > 2y use
256: * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
257: * where x1 = x with lower 32 bits cleared, x2 = x-x1; else
258: * 2. if x <= 2y use
259: * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
260: * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
261: * y1= y with lower 32 bits chopped, y2 = y-y1.
262: *
263: * NOTE: scaling may be necessary if some argument is too
264: * large or too tiny
265: *
266: * Special cases:
267: * hypot(x,y) is INF if x or y is +INF or -INF; else
268: * hypot(x,y) is NAN if x or y is NAN.
269: *
270: * Accuracy:
271: * hypot(x,y) returns sqrt(x^2+y^2) with error less
272: * than 1 ulps (units in the last place)
273: */
274:
275: /* #ifndef JAVA5 */
276: static double hypot(double x, double y) {
277: double a = x, b = y, t1, t2, w;
278: int j, ha, hb;
279: long la = (Double.doubleToLongBits(x) << 1) >>> 1;
280: long lb = (Double.doubleToLongBits(y) << 1) >>> 1;
281:
282: ha = (int) (la >>> 32); // high word of x
283: hb = (int) (lb >>> 32); // high word of y
284: if (hb > ha) {
285: j = ha;
286: ha = hb;
287: hb = j;
288: long l = la;
289: la = lb;
290: lb = l;
291: }
292: a = Double.longBitsToDouble(la); // a <- |a|
293: b = Double.longBitsToDouble(lb); // b <- |b|
294: /* Now a is max (abs(x), abs(y)) and b is min(abs(x), abs(y));
295: la and lb are the long bits of a and b;
296: and ha and hb are the high order bits of la and lb. */
297: if ((ha - hb) > 0x3c00000) // x/y > 2**60
298: return a + b;
299: int k = 0;
300: j = 0; // scale as high-order of double
301: if (ha > 0x5f300000) { // a>2**500
302: if (ha >= 0x7ff00000) { // Inf or NaN
303: w = a + b; // for sNaN
304: if ((la & 0xfffffffffffffL) == 0)
305: w = a;
306: if ((lb ^ 0x7ff0000000000000L) == 0)
307: w = b;
308: return w;
309: }
310: /* scale a and b by 2**-600 */
311: j = -0x25800000;
312: k += 600;
313: }
314: if (hb < 0x20b00000) { // b < 2**-500
315: if (hb <= 0x000fffff) { // subnormal b or 0
316: if (lb == 0)
317: return a;
318: t1 = Double.longBitsToDouble(0x7fd0000000000000L); // t1=2^1022
319: b *= t1;
320: a *= t1;
321: k -= 1022;
322: } else { // scale a and b by 2^600
323: k -= 600;
324: j = 0x25800000;
325: }
326: }
327: if (j != 0) {
328: ha += j;
329: hb += j;
330: la += (j << 32);
331: lb += (j << 32);
332: a = Double.longBitsToDouble(la);
333: b = Double.longBitsToDouble(lb);
334: }
335:
336: /* medium size a and b */
337: w = a - b;
338: if (w > b) {
339: t1 = Double.longBitsToDouble((long) ha << 32);
340: t2 = a - t1;
341: w = t1 * t1 - (b * (-b) - t2 * (a + t1));
342: } else {
343: a = a + a;
344: double y1 = Double.longBitsToDouble((long) hb << 32);
345: double y2 = b - y1;
346: t1 = Double
347: .longBitsToDouble(((long) (ha + 0x00100000)) << 32);
348: t2 = a - t1;
349: w = t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b));
350: }
351: w = Math.sqrt(w);
352: if (k != 0) { // t1 = 2^k
353: t1 = Double
354: .longBitsToDouble(0x3ff0000000000000L + ((long) k << 52));
355: w *= t1;
356: }
357: return w;
358: }
359:
360: /* #endif */
361:
362: /**
363: * @serialData Writes the real part, followed by the imaginary part.
364: * Both are written as doubles (using writeDouble).
365: */
366: public void writeExternal(ObjectOutput out) throws IOException {
367: out.writeDouble(real);
368: out.writeDouble(imag);
369: }
370:
371: public void readExternal(ObjectInput in) throws IOException,
372: ClassNotFoundException {
373: real = in.readDouble();
374: imag = in.readDouble();
375: }
376: }
|