main.cpp

斯坦福Energy211/CME211课《c++编程——地球科学科学家和工程师》的课件

// This is an application file, which uses functionality of
// the Complex class to implement new functionality

#include <iostream>
#include <cmath>

// This is needed to use Complex objects
#include "complex.h"

// Declarations of functions defined in this file
Complex exp( Complex z );
Complex pow( Complex z, double x );
Complex sqrt( Complex z );

using namespace std;

int main()
{
	// ipi = 0 + i*pi, where pi = 3.14159...
	Complex ipi(0, M_PI);
	// e^(i*pi) = -1
	cout << "exp(i*pi)=" << exp(ipi) << endl;
	
	// z = i, the imaginary unit
	Complex z(0.0, 1.0);
	// Print out the value of z
	cout << "z=" << z << endl;
	// Compute the square root of i, which is
	// sqrt(2)/2 + sqrt(2)/2 i
	cout << "sqrt(z)=" << sqrt(z) << endl;
	
	// A complex number z has n nth roots, which lie on
	// the circle in the complex plane of radius r^(1/n), 
	// where r = |z|.  These roots can be found by noting
	// that if 
	//
	// z = r*exp(i*t) = r*(cos(t) + i*sin(t)), 
	//
	// where (r,t) are the polar coordinates of z, then
	//
	// z1 = r^(1/n)*(cos(t/n)+i*sin(t/n))
	//
	// is an nth root of z, and so are z1*exp(2pi*i*j/n), 
	// for j=1,...,n-1. Here we compute the four 4th roots 
	// of z=i.
	int n = 4;
	// pow(z,1/n) computes z1 above
	Complex z1( pow(z, 1.0 / n) );
	// rootj = z1*e^(2pi*i*j/n), for j=0,1,...,n-1
	Complex rootj(z1);
	// shift = 2pi*i/n simplifies computation of rootj
	Complex shift(0, 2.0*M_PI/n);
	for ( int j = 0; j < n; j++ )
	{
		// Now output z1*e^(2pi*i*j/n), for j=0,1,...,n-1
		cout << "root " << j+1 << "=" << rootj << endl;
		rootj *= exp(shift);
	}
	return 0;
}

Complex exp( Complex z )
{
	// z = u + iv => exp(z) = exp(u)(cos(v) + i cos(v))
	Complex y(cos(z.Im()), sin(z.Im()));
	y *= std::exp(z.Re());
	return y;
}

Complex pow( Complex z, double x )
{
	// z = r*exp(it) => z^x = r^x*exp(itx)
	double arg = z.Arg();
	if ( arg == 0.0 )
		arg = arg + 2.0 * M_PI;
	Complex y(cos(arg*x), sin(arg*x));
	y *= std::pow(z.Abs(), x);
	return y;
}

Complex sqrt( Complex z )
{
	// Just use pow with exponent 1/2
	return pow( z, 0.5 );
}