c++复数类库.txt

将数学运算中的复数运算抽象为C++类的源代码

 C++复数类库

/*
Copyright (C) 1988 Free Software Foundation
    written by Doug Lea ([e-mail])

This file is part of the GNU C++ Library.  This library is free
software; you can redistribute it and/or modify it under the terms of
the GNU Library General Public License as published by the Free
Software Foundation; either version 2 of the License, or (at your
option) any later version.  This library is distributed in the hope
that it will be useful, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the GNU Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free Software
Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
*/

//# $Id: Complex.h,v 15.2 2000/12/17 22:16:34 wbrouw Exp $

#ifndef _Complex_h
#ifdef __GNUG__
#pragma interface
#endif
#define _Complex_h 1


#include <aips/aips.h>
#include <aips/Utilities/generic.h>
#include <iostream.h>
#include <aips/Mathematics/Math.h>
#include <stdlib.h>


// <summary> Complex numbers </summary>

// <synopsis>
// This documentation is lifted straight from the libg++ info page. Our classes
// are based on those classes, although we have made some modifications to them.
// The ANSI/ISO standard library has a complex<t> type. Until those classes are
// widely available, we will  use these complex classes. Although the 
// documentation below refers to Complex; DComplex (double precision) and
// IComplex (integral complex) also exist and freely interconvert.
// 
//    Class `Complex' is implemented in a way similar to that described by
// Stroustrup. In keeping with libg++ conventions, the class is named
// `Complex', not `complex'.  Complex arithmetic and relational operators
// are provided (`+, -, *, /, +=, -=, *=, /=, ==, !=').  Attempted
// division by (0, 0) triggers an exception.
// 
//    Complex numbers may be constructed and used in the following ways:
// <dl>
// <dt>Complex x;</dt>
// <dd>  Declares an uninitialized Complex. </dd>
// 
// <dt>Complex x = 2; Complex y(2.0);</dt>
// <dd>  Set x and y to the Complex value (2.0, 0.0); </dd>
// 
// <dt>Complex x(2, 3);</dt>
// <dd>  Sets x to the Complex value (2, 3); </dd>
// 
// <dt>Complex u(x); Complex v = x;</dt>
// <dd>  Set u and v to the same value as x. </dd>
// 
// <dt>double real(Complex& x);</dt>
// <dd>  returns the real part of x. </dd>
// 
// <dt>double imag(Complex& x);</dt>
// <dd>  returns the imaginary part of x. </dd>
// 
// <dt>double abs(Complex& x);</dt>
// <dd>  returns the magnitude of x. </dd>
// 
// <dt>double norm(Complex& x);</dt>
// <dd>  returns the square of the magnitude of x. </dd>
// 
// <dt>double arg(Complex& x);</dt>
// <dd>  returns the argument (amplitude) of x. </dd>
// 
// <dt>Complex polar(double r, double t = 0.0);</dt>
// <dd>  returns a Complex with abs of r and arg of t. </dd>
// 
// <dt>Complex conj(Complex& x);</dt>
// <dd>  returns the complex conjugate o </dd>f x.
// <dt>Complex cos(Complex& x);</dt>
// <dd>  returns the complex cosine of x. </dd>
// 
// <dt>Complex sin(Complex& x);</dt>
// <dd>  returns the complex sine of x. </dd>
// 
// <dt>Complex cosh(Complex& x);</dt>
// <dd>  returns the complex hyperbolic cosine of x. </dd>
// 
// <dt>Complex sinh(Complex& x);</dt>
// <dd>  returns the complex hyperbolic sine of x. </dd>
// 
// <dt>Complex exp(Complex& x);</dt>
// <dd>  returns the exponential of x. </dd>
// 
// <dt>Complex log(Complex& x);</dt>
// <dd>  returns the natural log of x. </dd>
// 
// <dt>Complex pow(Complex& x, long p);</dt>
// <dd>  returns x raised to the p power. </dd>
// 
// <dt>Complex pow(Complex& x, Complex& p);</dt>
// <dd>  returns x raised to the p power. </dd>
// 
// <dt>Complex sqrt(Complex& x);</dt>
// <dd>  returns the square root of x. </dd>
// 
// <dt> Complex min(Complex x,Complex y);
// <dd> Returns the minumum of x,y (using operator<=, i.e. the norm).
//
// <dt> Complex max(Complex x,Complex y);
// <dd> Returns the maximum of x,y (using operator>=, i.e. the norm).
//
// <dt>Bool near(Complex val1, Complex val2, Double tol = 1.0e-5);</dt>
// <dd>  returns whether val1 is relatively near val2 (see Math.h). </dd>
//
// <dt>Bool nearAbs(Complex val1, Complex val2, Double tol = 1.0e-5);</dt>
// <dd>  returns whether val1 is absolutely near val2 (see Math.h). </dd>
//
// <dt>ostream << x;</dt>
// <dd>  prints x in the form (re, im). </dd>
// 
// <dt>istream >> x;</dt>
//  <dd> reads x in the form (re, im), or just (re) or re in which case the
//      imaginary part is set to zero. </dd>
// </dl>
// 
// </synopsis>

// <group name="Complex numbers">

//# It seems as though the Sun Cfront compiler does not automatically
//# cast parameters to binary operators (?). By defining this, MANY
//# versions of the operators will be defined for each possibility...
#define COMPLEX_STUPID_COMPILER

#if defined(G_COMPLEX)
#undef G_COMPLEX
#endif

#define G_COMPLEX(type)  g_name2(type,G_COMPLEX)

#define G_COMPLEXdeclare2(type,CASTS)      \
class G_COMPLEX(type)        \
{          \
public:         \
          \
  type           re;        \
  type           im;        \
          \
          \
          \
  type           real() const { return re; }     \
  type           imag() const { return im; }     \
          \
  type           &real() { return re; }      \
  type           &imag() { return im; }      \
          \
                 G_COMPLEX(type)() : re((type) 0), im((type) 0) {}  \
                 G_COMPLEX(type)(const G_COMPLEX(type)& y) :re(y.re), im(y.im) {}\
                 G_COMPLEX(type)(type r, type i= (type) 0) :re(r), im(i) {} \
          \
                 CASTS        \
          \
                 ~G_COMPLEX(type)() {}      \
          \
  G_COMPLEX(type)& operator =  (G_COMPLEX(type) y) {re = y.re; im = y.im; return *this;} \
          \
          \
  G_COMPLEX(type)& operator += (const G_COMPLEX(type) y) {re += y.re;  im += y.im; return *this;}\
  G_COMPLEX(type)& operator += (type y) { re += y; return *this; }  \
          \
  G_COMPLEX(type)& operator -= (const G_COMPLEX(type) y) { re -= y.re;  im -= y.im; return *this;}\
  G_COMPLEX(type)& operator -= (type y) { re -= y; return *this;}  \
          \
  G_COMPLEX(type)& operator *= (const G_COMPLEX(type) y) {   \
   type r = re * y.re - im * y.im;   \
   im = re * y.im + im * y.re;    \
   re = r;       \
   return *this;       \
        }        \
  G_COMPLEX(type)& operator *= (type y) { re *=  y; im *=  y; return *this;} \
          \
  G_COMPLEX(type)& operator /= (const G_COMPLEX(type) y);    \
  G_COMPLEX(type)& operator /= (type y);      \
          \
  void           error(const char* msg) const;     \
};          \
          \
/*          \
* non-inline functions        \
*/          \
G_COMPLEX(type) operator /  (const G_COMPLEX(type) x, const G_COMPLEX(type) y);\
G_COMPLEX(type) operator /  (const G_COMPLEX(type) x, type y);   \
G_COMPLEX(type) operator /  (type x, const G_COMPLEX(type) y);   \
          \
G_COMPLEX(type) cos(const G_COMPLEX(type) x);    \
G_COMPLEX(type) sin(const G_COMPLEX(type) x);    \
          \
G_COMPLEX(type) cosh(const G_COMPLEX(type) x);    \
G_COMPLEX(type) sinh(const G_COMPLEX(type) x);    \
          \
G_COMPLEX(type) exp(const G_COMPLEX(type) x);    \
G_COMPLEX(type) log(const G_COMPLEX(type) x);    \
          \
G_COMPLEX(type) sqrt(const G_COMPLEX(type) x);    \
          \
G_COMPLEX(type) log10(const G_COMPLEX(type) x);     \
          \
G_COMPLEX(type) pow(const G_COMPLEX(type) x, int p);    \
G_COMPLEX(type) pow(const G_COMPLEX(type) x, const G_COMPLEX(type) p); \
G_COMPLEX(type) pow(const G_COMPLEX(type) x, double y);   \
             \
istream&  operator >> (istream& s, G_COMPLEX(type)& x);    \
ostream&  operator << (ostream& s, const G_COMPLEX(type) x);   \
          \
          \
/*          \
* ACCESSORS         \
*/          \
inline type real(const G_COMPLEX(type) x)     \
{          \
  return x.re;        \
}          \
          \
inline type imag(const G_COMPLEX(type) x)     \
{          \
  return x.im;        \
}          \
          \
/*          \
* math.h FUNCTIONS        \
*/          \
inline double abs(const G_COMPLEX(type) x)     \
{          \
  return hypot((double)x.re, (double)x.im);    \
}          \
          \
inline double fabs(const G_COMPLEX(type) x)     \
{          \
  return abs(x);        \
}          \
          \
inline double norm(const G_COMPLEX(type) x)     \
{          \
  return (x.re * x.re + x.im * x.im);    \
}          \
          \
inline double arg(const G_COMPLEX(type) x)     \
{          \
  return atan2((double)x.im, (double)x.re);    \
}          \
          \
G_COMPLEX(type) cube(const G_COMPLEX(type) val);    \
          \
inline G_COMPLEX(type) square(const G_COMPLEX(type) val)   \
{          \
    type r = val.re;       \
    type i = val.im;       \
    return G_COMPLEX(type)( r * r - i * i, r * i * 2 );    \
}          \
          \
/*          \
* EQUALITY OPERATORS        \
*/          \
inline int  operator == (const G_COMPLEX(type) x, const G_COMPLEX(type) y) \
{          \
  return x.re == y.re && x.im == y.im;    \
}          \
          \
inline int  operator == (const G_COMPLEX(type) x, type y)   \
{          \
  return x.im == 0.0 && x.re == y;     \
}          \
          \
inline int  operator != (const G_COMPLEX(type) x, const G_COMPLEX(type) y) \
{          \
  return x.re != y.re || x.im != y.im;    \
}          \
          \
inline int  operator != (const G_COMPLEX(type) x, type y)   \
{          \
  return x.im != 0.0 || x.re != y;     \
}          \
          \
/*          \
* COMPARISON OPERATORS (using norm)      \
*/          \
inline int operator >= (const G_COMPLEX(type) x, const G_COMPLEX(type) y) { \
  return norm(x) >= norm(y);       \
}          \
          \
inline int operator >= (const G_COMPLEX(type) x, const type y) {  \
  return norm(x) >= y*y;       \
}          \
          \
inline int operator >= (const type x, const G_COMPLEX(type) y) {  \
  return x*x >= norm(y);       \
}          \
          \
inline int operator > (const G_COMPLEX(type) x, const G_COMPLEX(type) y) { \
  return norm(x) > norm(y);       \
}          \
          \
inline int operator > (const G_COMPLEX(type) x, const type y) {   \
  return norm(x) > y*y;        \
}          \
          \
inline int operator > (const type x, const G_COMPLEX(type) y) {   \
  return x*x > norm(y);        \
}          \
          \
inline int operator <= (const G_COMPLEX(type) x, const G_COMPLEX(type) y) { \
  return norm(x) <= norm(y);       \
}          \
          \
inline int operator <= (const G_COMPLEX(type) x, const type y) {  \
  return norm(x) <= y*y;       \
}          \
          \
inline int operator <= (const type x, const G_COMPLEX(type) y) {  \
  return x*x <= norm(y);       \
}          \