From Odwiki
/****************************************************************************
File: wikiComplex.h
Authors: The od[wiki] community at
http://www.odforce.net/wiki
Description: VEX package for complex number arithmetic.
Complex numbers in this VEX implementation are
represented using the 'vector' data type, where
the x-component is used for the real part, and
the y-component for the imaginary part. However; the
module also defines the data type 'complex' so as to
make application code a little clearer and make the
underlying vector type transparent -- i.e: application
code should use the type 'complex' after including
this module.
Wherever it makes sense, functions are given alternative
versions using the suffixes 'R' (Real) and 'C' (complex)
to denote various combinations of data types for their
arguments. This lets us do complex arithmetic with mixed
data types... but due to the fact that we can't overload,
the function names start getting a little complicated, but
the overall pattern is consistent.
-----------------------------------------------------------------------------
Part of the small but growing od[wiki] library of VEX functions...
This software is placed in the public domain and is provided as is
without express or implied warranty.
*****************************************************************************/
#ifndef wikiComplex_h__GUARD
#define wikiComplex_h__GUARD
#include <wikiMathConst.h>
//------------------------------------------------------------------
// Define the data type -- obviously this makes 'complex' a
// reserved word after including this module.
//------------------------------------------------------------------
#define complex vector
//------------------------------------------------------------------
// Create a complex number
//------------------------------------------------------------------
complex cx_set(float re, im) {
return set(re,im,0);
}
complex cx_setR(float re) {
return set(re,0,0);
}
complex cx_setV(vector z) {
return set(z.x,z.y,0);
}
//------------------------------------------------------------------
// Retrieve the real part of the complex number z
//------------------------------------------------------------------
float cx_re(complex z) {
return z.x;
}
//------------------------------------------------------------------
// Retrieve the imaginary part of the complex number z
//------------------------------------------------------------------
float cx_im(complex z) {
return z.y;
}
//------------------------------------------------------------------
// The distance from origin in polar coords -- a.k.a
// 'complex modulus', or 'norm', or 'length'.
//------------------------------------------------------------------
float cx_mod(complex z) {
float l=0.0;
if(z.x!=0.0 || z.y!=0.0) {
l = length(set(z.x, z.y, 0.0));
}
return l;
}
float cx_modR(float x) {
return x<0.0 ? -x : x;
}
float cx_mod2(complex z) {
float l=0.0;
if(z.x!=0.0 || z.y!=0.0) {
l = length2(set(z.x, z.y, 0.0));
}
return l;
}
float cx_mod2R(float x) {
return x*x;
}
#define cx_norm cx_mod
#define cx_normR cx_modR
#define cx_norm2 cx_mod2
#define cx_norm2R cx_mod2R
#define cx_length cx_mod
#define cx_lengthR cx_modR
#define cx_length2 cx_mod2
#define cx_length2R cx_mod2R
//------------------------------------------------------------------
// Absolute Value (the magnitude of the complex number)
// An alternative implementation of cx_mod adapted from
// Numerical Recipes.
//------------------------------------------------------------------
float cx_abs(complex z) {
float result;
float x=abs(z.x), y=abs(z.y), d;
if(x==0.0 && y==0) {
result=0.0;
} else if(x>=y) {
d = y/x;
result = x*sqrt(1.0+d*d);
} else {
d = x/y;
result = y*sqrt(1.0+d*d);
}
return result;
}
//------------------------------------------------------------------
// 'Argument' (angle in radians with the x-axis in polar coords)
//------------------------------------------------------------------
float cx_arg(complex z) {
return atan(z.y, z.x);
}
//------------------------------------------------------------------
// Conjugate (mult of conjugate pairs produce a real)
// Conjugate appears "flipped" across the real axis.
//------------------------------------------------------------------
complex cx_conj(complex z) {
return cx_set(z.x, -z.y);
}
//------------------------------------------------------------------
// Inverse
//------------------------------------------------------------------
complex cx_inv(complex z) {
float re,im,s;
if(abs(z.x)>=abs(z.y)) {
s = 1.0 / (z.x + z.y*(z.y/z.x));
re = s;
im = s * (-z.y/z.x);
} else {
s = 1.0 / (z.x*(z.x/z.y) + z.y);
re = s * (z.x/z.y);
im = -s;
}
return cx_set(re,im);
}
//------------------------------------------------------------------
// Negation (-z)
//------------------------------------------------------------------
complex cx_neg(complex z) {
return cx_set(-z.x,-z.y);
}
complex cx_negR(float x) {
return cx_setR(-x);
}
//------------------------------------------------------------------
// Addition (a+
//------------------------------------------------------------------
complex cx_add(complex a, {
return cx_set(a.x+b.x, a.y+b.y);
}
complex cx_addR(complex a; float {
return cx_set(a.x+b, a.y);
}
complex cx_addRC(float a; complex {
return cx_set(a+b.x, b.y);
}
//------------------------------------------------------------------
// Subtraction (a-
//------------------------------------------------------------------
complex cx_sub(complex a, {
return cx_set(a.x-b.x, a.y-b.y);
}
complex cx_subR(complex a; float {
return cx_set(a.x-b, a.y);
}
complex cx_subRC(float a; complex {
return cx_set(a-b.x, -b.y);
}
//------------------------------------------------------------------
// Multiplication (a*
//------------------------------------------------------------------
complex cx_mul(complex a, {
return cx_set(a.x*b.x - a.y*b.y, a.x*b.y + a.y*b.x);
}
complex cx_mulR(complex a; float {
return cx_set(a.x*b, a.y*;
}
complex cx_mulRC(float a; complex {
return cx_set(a*b.x, a*b.y);
}
//------------------------------------------------------------------
// Division (a/b)
// Adapted from numerical Recipes
// This does not check for division by zero!
//------------------------------------------------------------------
complex cx_div(complex a, {
float re,im,s;
if(abs(b.x) >= abs(b.y)) {
float byox = b.y/b.x;
s = 1.0 / (b.x + b.y*byox);
re = s * (a.x + a.y*byox);
im = s * (a.y - a.x*byox);
} else {
float bxoy = b.x/b.y;
s = 1.0 / (b.x*bxoy + b.y);
re = s * (a.x*bxoy + a.y);
im = s * (a.y*bxoy - a.x);
}
return cx_set(re,im);
}
complex cx_divR(complex a; float {
float l2 = b*b;
return l2==0.0 ? 0.0 : cx_set((a.x*/l2, (a.y*/l2);
}
complex cx_divRC(float a; complex {
return cx_div(cx_setR(a),;
}
//------------------------------------------------------------------
// Exponential
//------------------------------------------------------------------
complex cx_exp(complex z) {
return cx_set(exp(z.x)*cos(z.y), exp(z.x)*sin(z.y));
}
complex cx_expR(float x) {
return cx_setR(exp(x));
}
//------------------------------------------------------------------
// Power
//------------------------------------------------------------------
// complex raised to complex power
complex cx_pow(complex base, exp) {
float re = log(cx_abs(base)),
im = cx_arg(base),
re2 = (re*exp.x) - (im*exp.y),
im2 = (re*exp.y) + (im*exp.x),
ere = exp(re2);
return cx_set(ere*cos(im2), ere*sin(im2));
}
// complex raised to a real power
complex cx_powR(complex base; float exp) {
float re = log(cx_abs(base)),
im = cx_arg(base),
ere = exp(re);
return cx_set(ere*cos(im), ere*sin(im));
}
// real raised to complex power
complex cx_powRC(float base; complex exp) {
float re = log(cx_abs(base)),
im = WM_PI_2, // atan2(0,base)
re2 = (re*exp.x) - (im*exp.y),
im2 = (re*exp.y) + (im*exp.x),
ere = exp(re2);
return cx_set(ere*cos(im2), ere*sin(im2));
}
//------------------------------------------------------------------
// Natural Logarithm (principal branch, where {-pi < arg <= pi} )
//------------------------------------------------------------------
complex cx_log(complex z) {
return cx_set(log(cx_mod(z)), cx_arg(z));
}
complex cx_logR(float x) {
return cx_setR(log(x));
}
//------------------------------------------------------------------
// Square Root (one of the two roots)
//------------------------------------------------------------------
complex cx_sqrt(complex z) {
float re,im,tmp;
float mag = cx_abs(z);
if(mag>0.0) {
if(z.x>0.0) {
tmp = sqrt(0.5*(mag+z.x));
re = tmp;
im = 0.5*z.y / tmp;
} else {
tmp = sqrt(0.5*(mag-z.x));
if(z.y<0.0) tmp=-tmp;
re = 0.5*z.y / tmp;
im = tmp;
}
} else {
re = im = 0.0;
}
return cx_set(re,im);
}
complex cx_sqrtR(float x) {
complex z;
if(x>=0) {
z = x==0.0 ? 0.0 : cx_setR(sqrt(x));
} else {
z = cx_sqrt(cx_setR(x));
}
return z;
}
//------------------------------------------------------------------
// Nth Root
//------------------------------------------------------------------
complex cx_root(complex z; float n) {
float p = 1.0/n;
float th_n = cx_arg(z)*p;
return cx_set(pow(cx_mod(z),p)*cos(th_n), sin(th_n));
}
complex cx_rootR(float x; float n) {
return cx_root(cx_setR(x),n);
}
//------------------------------------------------------------------
// Sine
//------------------------------------------------------------------
complex cx_sin(complex z) {
return cx_set(cosh(z.y)*sin(z.x), sinh(z.y)*cos(z.x));
}
complex cx_sinR(float x) {
return cx_setR(sin(x));
}
//------------------------------------------------------------------
// Cosine
//------------------------------------------------------------------
complex cx_cos(complex z) {
return cx_set(cosh(z.y)*cos(z.x), -sinh(z.y)*sin(z.x));
}
complex cx_cosR(float x) {
return cx_setR(cos(x));
}
//------------------------------------------------------------------
// Tangent
//------------------------------------------------------------------
complex cx_tan(complex z) {
return cx_div(cx_sin(z),cx_cos(z));
}
complex cx_tanR(float x) {
return cx_setR(tan(x));
}
//------------------------------------------------------------------
// Hyperbolic Sine
//------------------------------------------------------------------
complex cx_sinh(complex z) {
return cx_set(sinh(z.x)*cos(z.y),cosh(z.x)*sin(z.y));
}
complex cx_sinhR(float x) {
return cx_setR(sinh(x));
}
//------------------------------------------------------------------
// Hyperbolic Cosine
//------------------------------------------------------------------
complex cx_cosh(complex z) {
return cx_set(cosh(z.x)*cos(z.y),sinh(z.x)*sin(z.y));
}
complex cx_coshR(float z) {
return cx_setR(cosh(x));
}
#endif // End wikiComplex_h__GUARD
