DComplex


1   // Copyright (c) 1997  Per M.A. Bothner.
2   // This is free software;  for terms and warranty disclaimer see ./COPYING.
3   
4   package gnu.math;
5   import java.io.*;
6   
7   /** A complex number using rectangular (Cartesian) plain double values.
8    * @author Per Bothner
9    * @author Some algorithms were transcribed from GNU libstdc++,
10   * written by Jason Merrill.
11   * Also see below for copyrights for functions taken from fdlib and f2c.
12   */
13  
14  public class DComplex extends Complex implements Externalizable
15  { 
16    double real;
17    double imag;
18  
19    public DComplex ()
20    {
21    }
22  
23    public DComplex (double real, double imag)
24    {
25      this.real = real;
26      this.imag = imag;
27    }
28  
29    public RealNum re () { return new DFloNum (real); }
30    public double doubleValue() { return real; }
31    public RealNum im () { return new DFloNum (imag); }
32    public double doubleImagValue () { return imag; }
33  
34    public boolean equals (Object   obj)
35    {
36      if (obj == null || ! (obj instanceof Complex))
37        return false;
38      Complex y = (Complex)obj;
39      return y.unit() == Unit.Empty
40        && (Double.doubleToLongBits(real)
41        == Double.doubleToLongBits(y.reValue()))
42        && (Double.doubleToLongBits(imag)
43        == Double.doubleToLongBits(y.imValue()));
44    }
45  
46    public String   toString ()
47    {
48      String   prefix = "";
49  
50      String   reString;
51      if (real == 1.0/0.0)
52        {
53      prefix = "#i"; reString = "1/0";
54        }
55      else if (real == -1.0/0.0)
56        {
57      prefix = "#i"; reString = "-1/0";
58        }
59      else if (Double.isNaN (real))
60        {
61      prefix = "#i"; reString = "0/0";
62        }
63      else
64        reString = Double.toString (real);
65  
66      if (Double.doubleToLongBits (imag) == 0)  // i.e. imag is 0.0 and not -0.0
67        return prefix + reString;
68  
69      String   imString;
70      if (imag == 1.0/0.0)
71        {
72      prefix = "#i"; imString = "+1/0i";
73        }
74      else if (imag == -1.0/0.0)
75        {
76      prefix = "#i"; imString = "-1/0i";
77        }
78      else if (Double.isNaN (imag))
79        {
80      prefix = "#i"; imString = "+0/0i";
81        }
82      else
83        {
84      imString = Double.toString (imag) + "i";
85      if (imString.charAt (0) != '-')
86        imString = "+" + imString;
87        }
88  
89      return ((Double.doubleToLongBits (real) == 0 ? prefix : prefix + reString)
90              + imString);
91    }
92  
93    public String   toString (int radix)
94    {
95      if (radix == 10)
96        return toString ();
97      return "#d" + toString ();
98    }
99  
100   // All transcendental complex functions return DComplex
101 
102   public final Numeric neg () { return new DComplex (-real, -imag); }
103 
104   public Numeric add (Object   y, int k)
105   {
106     if (y instanceof Complex)
107       {
108         Complex yc = (Complex)y;
109     if (yc.dimensions() != Dimensions.Empty)
110       throw new ArithmeticException   ("units mis-match");
111     return new DComplex (real + k * yc.reValue(),
112                  imag + k * yc.imValue());
113       }
114     return ((Numeric)y).addReversed(this, k);
115   }
116 
117   public Numeric mul (Object   y)
118   {
119     if (y instanceof Complex)
120       {
121         Complex yc = (Complex)y;
122     if (yc.unit() == Unit.Empty)
123       {
124         double y_re = yc.reValue();
125         double y_im = yc.imValue();
126         return new DComplex (real * y_re - imag * y_im,
127                  real * y_im + imag * y_re);
128       }
129     return Complex.times(this, yc);
130       }
131     return ((Numeric)y).mulReversed(this);
132   }
133 
134   public Numeric div (Object   y)
135   {
136     if (y instanceof Complex)
137       {
138     Complex yc = (Complex) y;
139     return div (real, imag,
140             yc.doubleValue(), yc.doubleImagValue());
141       }
142     return ((Numeric)y).divReversed(this);
143   }
144 
145   public static DComplex power (double x_re, double x_im,
146                 double y_re, double y_im)
147   {
148     double h;
149     /* #ifdef JAVA5 */
150     // h = Math.hypot(x_re, x_im);
151     /* #else */
152     h = DComplex.hypot(x_re, x_im);
153     /* #endif */
154     double logr = Math.log (h);
155     double t = Math.atan2 (x_im, x_re);
156     double r = Math.exp (logr * y_re - y_im * t);
157     t = y_im * logr + y_re * t;
158     return Complex.polar (r, t);
159   }
160 
161   public static Complex log (double x_re, double x_im)
162   {
163     double h;
164     /* #ifdef JAVA5 */
165     // h = Math.hypot(x_re, x_im);
166     /* #else */
167     h = DComplex.hypot(x_re, x_im);
168     /* #endif */
169     return make(Math.log(h), Math.atan2(x_im, x_re));
170   }
171 
172   // The code below is adapted from f2c's libF77, and is subject to this
173   // copyright:
174  
175   /****************************************************************
176     Copyright 1990, 1991, 1992, 1993 by AT&T Bell Laboratories and Bellcore.
177  
178     Permission to use, copy, modify, and distribute this software
179     and its documentation for any purpose and without fee is hereby
180     granted, provided that the above copyright notice appear in all
181     copies and that both that the copyright notice and this
182     permission notice and warranty disclaimer appear in supporting
183     documentation, and that the names of AT&T Bell Laboratories or
184     Bellcore or any of their entities not be used in advertising or
185     publicity pertaining to distribution of the software without
186     specific, written prior permission.
187  
188     AT&T and Bellcore disclaim all warranties with regard to this
189     software, including all implied warranties of merchantability
190     and fitness.  In no event shall AT&T or Bellcore be liable for
191     any special, indirect or consequential damages or any damages
192     whatsoever resulting from loss of use, data or profits, whether
193     in an action of contract, negligence or other tortious action,
194     arising out of or in connection with the use or performance of
195     this software.
196     ****************************************************************/
197 
198   public static DComplex div (double x_re, double x_im,
199                   double y_re, double y_im)
200   {
201     double ar = Math.abs (y_re);
202     double ai = Math.abs (y_im);
203     double nr, ni;
204     double t, d;
205     if (ar <= ai)
206       {
207     t = y_re / y_im;
208     d = y_im * (1 + t*t);
209     nr = x_re * t + x_im;
210     ni = x_im * t - x_re;
211       }
212     else
213       {
214     t = y_im / y_re;
215     d = y_re * (1 + t*t);
216     nr = x_re + x_im * t;
217     ni = x_im - x_re * t;
218       }
219     return new DComplex (nr / d, ni / d);
220   }
221   
222   public static Complex sqrt (double x_re, double x_im)
223   {
224     /* #ifdef JAVA5 */
225     // double r = Math.hypot(x_re, x_im);
226     /* #else */
227     double r = DComplex.hypot(x_re, x_im);
228     /* #endif */
229     double nr, ni;
230     if (r == 0.0)
231       nr = ni = r;
232     else if (x_re > 0)
233       {
234     nr = Math.sqrt (0.5 * (r + x_re));
235     ni = x_im / nr / 2;
236       }
237     else
238       {
239     ni = Math.sqrt (0.5 * (r - x_re));
240     if (x_im < 0)
241       ni = - ni;
242     nr = x_im / ni / 2;
243       }
244     return new DComplex (nr, ni);
245   }
246 
247   // Transcribed from:
248   // http://netlib.bell-labs.com/netlib/fdlibm/e_hypot.c.Z
249   /*
250    * ====================================================
251    * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
252    *
253    * Developed at SunSoft, a Sun Microsystems, Inc. business.
254    * Permission to use, copy, modify, and distribute this
255    * software is freely granted, provided that this notice 
256    * is preserved.
257    * ====================================================
258    */
259   /* __ieee754_hypot(x,y)
260    *
261    * Method :                  
262    *      If (assume round-to-nearest) z=x*x+y*y 
263    *      has error less than sqrt(2)/2 ulp, than 
264    *      sqrt(z) has error less than 1 ulp (exercise).
265    *
266    *      So, compute sqrt(x*x+y*y) with some care as 
267    *      follows to get the error below 1 ulp:
268    *
269    *      Assume x>y>0;
270    *      (if possible, set rounding to round-to-nearest)
271    *      1. if x > 2y  use
272    *              x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
273    *      where x1 = x with lower 32 bits cleared, x2 = x-x1; else
274    *      2. if x <= 2y use
275    *              t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
276    *      where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, 
277    *      y1= y with lower 32 bits chopped, y2 = y-y1.
278    *              
279    *      NOTE: scaling may be necessary if some argument is too 
280    *            large or too tiny
281    *
282    * Special cases:
283    *      hypot(x,y) is INF if x or y is +INF or -INF; else
284    *      hypot(x,y) is NAN if x or y is NAN.
285    *
286    * Accuracy:
287    *      hypot(x,y) returns sqrt(x^2+y^2) with error less 
288    *      than 1 ulps (units in the last place) 
289    */
290 
291   /* #ifndef JAVA5 */
292   static double hypot (double x, double y)
293   {
294     double a=x,b=y,t1,t2,w;
295     int j,ha,hb;
296     long la = (Double.doubleToLongBits(x) << 1) >>> 1;
297     long lb = (Double.doubleToLongBits(y) << 1) >>> 1;
298 
299     ha = (int)(la >>> 32);        // high word of  x
300     hb = (int)(lb >>> 32);        // high word of  y
301     if (hb > ha)
302       {
303         j=ha; ha=hb; hb=j;
304         long l=la; la=lb; lb=l;
305       }
306     a = Double.longBitsToDouble(la);   // a <- |a|
307     b = Double.longBitsToDouble(lb);   // b <- |b|
308     /* Now a is max (abs(x), abs(y)) and b is min(abs(x), abs(y));
309        la and lb are the long bits of a and b;
310        and ha and hb are the high order bits of la and lb. */
311     if ((ha-hb) > 0x3c00000) // x/y > 2**60
312       return a+b;
313     int k=0;
314     j = 0;  // scale as high-order of double
315     if (ha > 0x5f300000)
316       {   // a>2**500
317         if (ha >= 0x7ff00000)
318           {       // Inf or NaN
319             w = a+b;                 // for sNaN
320             if ((la & 0xfffffffffffffL) == 0)
321               w = a;
322             if ((lb^0x7ff0000000000000L) == 0)
323               w = b;
324             return w;
325           }
326         /* scale a and b by 2**-600 */
327         j = -0x25800000;  k += 600;
328       }
329     if (hb < 0x20b00000)
330       {   // b < 2**-500
331         if (hb <= 0x000fffff)
332           {      // subnormal b or 0  
333             if (lb == 0)
334               return a;
335             t1 = Double.longBitsToDouble(0x7fd0000000000000L); // t1=2^1022
336             b *= t1;
337             a *= t1;
338             k -= 1022;
339           }
340         else
341           {            // scale a and b by 2^600
342             k -= 600;
343             j = 0x25800000;
344           }
345       }
346     if (j != 0)
347       {
348         ha += j; hb += j;
349         la += (j << 32);  lb += (j << 32);
350         a = Double.longBitsToDouble(la);
351         b = Double.longBitsToDouble(lb);
352       }
353 
354     /* medium size a and b */
355     w = a-b;
356     if (w>b)
357       {
358         t1 = Double.longBitsToDouble ((long) ha << 32);
359         t2 = a-t1;
360         w  = t1*t1-(b*(-b)-t2*(a+t1)); 
361       }
362     else
363       {
364         a = a+a;
365         double y1 = Double.longBitsToDouble ((long) hb << 32);
366         double y2 = b - y1;
367         t1 = Double.longBitsToDouble (((long) (ha+0x00100000)) << 32);
368         t2 = a - t1;
369         w  = t1*y1-(w*(-w)-(t1*y2+t2*b));
370       }
371     w = Math.sqrt(w);
372     if(k!=0)
373       { // t1 = 2^k
374         t1 = Double.longBitsToDouble (0x3ff0000000000000L + ((long)k << 52));
375         w *= t1;
376       }
377     return w;
378   }
379   /* #endif */
380 
381   /**
382    * @serialData Writes the real part, followed by the imaginary part.
383    *   Both are written as doubles (using writeDouble).
384    */
385   public void writeExternal(ObjectOutput out) throws IOException
386   {
387     out.writeDouble(real);
388     out.writeDouble(imag);
389   }
390 
391   public void readExternal(ObjectInput in)
392     throws IOException, ClassNotFoundException  
393   {
394     real = in.readDouble();
395     imag = in.readDouble();
396   }
397 }
398
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags