complex_arith.h

Go to the documentation of this file.

#ifndef _COMPLEX_ARITH_H_
#define _COMPLEX_ARITH_H_

#include <math.h>

typedef struct _complex_type
{
    double real;    
    double imag;    
} _complex_type;

#define complex     _complex_type

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex make_complex(double a, double b)
{
    complex z;
    z.real = a;
    z.imag = b;
    return z;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_conj(complex z)
{
    complex result;
    result.real = z.real;
    result.imag = -z.imag;
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_add(complex w, complex z)
{
    complex result;
    result.real = w.real + z.real;
    result.imag = w.imag + z.imag;
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_neg(complex z)
{
    complex result;
    result.real = -z.real;
    result.imag = -z.imag;
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_sub(complex w, complex z)
{
    complex result;
    result.real = w.real - z.real;
    result.imag = w.imag - z.imag;
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_mul(complex w, complex z)
{
    complex result;
    result.real = w.real * z.real - w.imag * z.imag;
    result.imag = w.real * z.imag + w.imag * z.real;
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_inv(complex z)
{
    complex result;
    
    /* There are two mathematically-equivalent formulas for the inverse. For
     * accuracy, choose the formula with the smaller value of |ratio|.
     */
    if (fabs(z.real) >= fabs(z.imag))
    {
        double ratio = z.imag / z.real;
        double denom = z.real + z.imag * ratio;
        result.real = 1 / denom;
        result.imag = -ratio / denom;
    }
    else
    {
        double ratio = z.real / z.imag;
        double denom = z.real * ratio + z.imag;
        result.real = ratio / denom;
        result.imag = -1 / denom;
    }
    
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_div(complex w, complex z)
{
    complex result;
    
    /* For accuracy, choose the formula with the smaller value of |ratio|. */
    if (fabs(z.real) >= fabs(z.imag))
    {
        double ratio = z.imag / z.real;
        double denom = z.real + z.imag * ratio;
        result.real = (w.real + w.imag * ratio) / denom;
        result.imag = (w.imag - w.real * ratio) / denom;
    }
    else
    {
        double ratio = z.real / z.imag;
        double denom = z.real * ratio + z.imag;
        result.real = (w.real * ratio + w.imag) / denom;
        result.imag = (w.imag * ratio - w.real) / denom;
    }
    
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static double c_mag(complex z)
{
    z.real = fabs(z.real);
    z.imag = fabs(z.imag);
    
    /* For accuracy, choose the formula with the smaller value of |ratio|. */
    if (z.real >= z.imag)
    {
        double ratio = z.imag / z.real;
        return z.real * sqrt(1.0 + ratio * ratio);
    }
    else
    {
        double ratio = z.real / z.imag;
        return z.imag * sqrt(1.0 + ratio * ratio);
    }
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static double c_arg(complex z)
{
    return atan2(z.imag, z.real);
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_pow(complex w, complex z)
{
    complex result;
    double mag_w = c_mag(w);
    double arg_w = c_arg(w);
    double mag = pow(mag_w, z.real) * exp(-z.imag * arg_w);
    double arg = z.real * arg_w + z.imag * log(mag_w);
    result.real = mag * cos(arg);
    result.imag = mag * sin(arg);
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_real_pow(complex w, double x)
{
    complex result;
    double mag = pow(c_mag(w), x);
    double arg = c_arg(w) * x;
    result.real = mag * cos(arg);
    result.imag = mag * sin(arg);
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_sqrt(complex z)
{
    double r = c_mag(z);
    complex result;
    result.real = sqrt((r + z.real) / 2);
    result.imag = sqrt((r - z.real) / 2);
    
    if (z.imag < 0)
        result.imag = -result.imag;
    
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_exp(complex z)
{
    double r = exp(z.real);
    complex result;
    result.real = r * cos(z.imag);
    result.imag = r * sin(z.imag);
    return result;
}

#ifdef __GNUC__
__attribute__((pure,unused))
#endif

static complex c_log(complex z)
{
    complex result;
    result.real = log(c_mag(z));
    result.imag = c_arg(z);
    return result;
}

#endif /* _COMPLEX_ARITH_H_ */