This is complex.c in view mode; [Download] [Up]
/*
complex - functions and definitions for complex types
Copyright (C) 1995, Adam Fedor
*/
#include <math.h>
#include "MathArray/complex.h"
/* Addition, multiplication of complex numbers */
complex_double
c_add (complex_double x, complex_double y)
{
x.real += y.real;
x.imag += y.imag;
return x;
}
complex_double
c_sub (complex_double x, complex_double y)
{
x.real -= y.real;
x.imag -= y.imag;
return x;
}
complex_double
c_mult (complex_double x, complex_double y)
{
complex_double t;
t.real = (x.real*y.real - x.imag*y.imag);
t.imag = (x.real*y.imag + x.imag*y.real);
return t;
}
complex_double
c_div (complex_double x, complex_double y)
{
complex_double t;
double r;
double denom;
/* FIXME: raise exception if denom == 0? */
if (fabs(y.real) >= fabs(y.imag)) {
r = y.imag / y.real;
denom = y.real + r * y.imag;
t.real = (x.real + r * x.imag) / denom;
t.imag = (x.imag - r * x.real) / denom;
} else {
r = y.real / y.imag;
denom = y.imag + r * y.real;
t.real = (x.real * r + x.imag) / denom;
t.imag = (x.imag * r - x.real) / denom;
}
return t;
}
/* Functions you typically find with complex numbers */
double
c_abs (complex_double z)
{
double x, y, t, temp;
x = fabs(z.real);
y = fabs(z.imag);
if (x == 0.0)
t = y;
else if (y == 0.0)
t = x;
else if (x > y) {
temp = y / x;
t = x * sqrt(1.0+temp*temp);
} else {
temp = x / y;
t = y * sqrt(1.0+temp*temp);
}
return t;
}
double
c_imag (complex_double z)
{
return z.imag;
}
double
c_real (complex_double z)
{
return z.real;
}
complex_double
c_set (double real, double imag)
{
complex_double d;
d.real = real;
d.imag = imag;
return d;
}
complex_double
c_f2d (complex_float f)
{
complex_double d;
d.real = f.real;
d.imag = f.imag;
return d;
}
complex_float
c_d2f (complex_double d)
{
complex_float f;
f.real = d.real;
f.imag = d.imag;
return f;
}
complex_double
ccos (complex_double z)
{
complex_double t;
t.real = cos(z.real) * cosh(z.imag);
t.imag = -sin(z.real) * sinh(z.imag);
return t;
}
complex_double
cexp (complex_double z)
{
complex_double t;
t.real = exp(z.real) * cos(z.imag);
t.imag = exp(z.real) * sin(z.imag);
return t;
}
complex_double
clog (complex_double z)
{
complex_double t;
t.real = log(c_abs(z));
t.imag = atan2(z.imag, z.real);
return t;
}
complex_double
clog10 (complex_double z)
{
complex_double t;
/* FIXME: correct implementation of clog10? */
t.real = log10(c_abs(z));
t.imag = atan2(z.imag, z.real);
return t;
}
complex_double
conj (complex_double z)
{
z.imag = -z.imag;
return z;
}
complex_double
cpow (complex_double x, complex_double y)
{
return cexp(c_mult(y, clog(x)));
}
complex_double
csin (complex_double z)
{
complex_double t;
t.real = sin(z.real) * cosh(z.imag);
t.imag = cos(z.real) * sinh(z.imag);
return t;
}
complex_double
csqrt (complex_double z)
{
complex_double t;
if ((z.real == 0.0) && (z.imag == 0.0)) {
t.real = 0.0;
t.imag = 0.0;
} else {
double x, y, w, r;
x = fabs(z.real);
y = fabs(z.imag);
if (x >= y) {
r = y/x;
w = sqrt(x)*sqrt(0.5*(1.0+sqrt(1.0+r*r)));
} else {
r = x/y;
w = sqrt(y)*sqrt(0.5*(r+sqrt(1.0+r*r)));
}
if (z.real >= 0.0) {
t.real = w;
t.imag = z.imag/(2.0*w);
} else {
t.imag = (z.imag >= 0) ? w : -w;
t.real = z.imag/(2.0*t.imag);
}
}
return t;
}
/* Not so common functions of complex numbers */
complex_double
cacos( complex_double x )
{
complex_double one;
complex_double nimag;
one.real = 1.0;
one.imag = 0.0;
nimag.real = 0.0;
nimag.imag = -1.0;
return c_mult(nimag, clog( c_add(x, csqrt( c_sub( c_mult(x, x), one)))));
}
complex_double
casin( complex_double x )
{
complex_double one;
complex_double nimag;
complex_double ix;
one.real = 1.0;
one.imag = 0.0;
nimag.real = 0.0;
nimag.imag = 1.0;
ix = c_mult(nimag, x);
nimag.real = 0.0;
nimag.imag = -1.0;
return c_mult(nimag, clog( c_add(ix, csqrt( c_sub( c_mult(x, x), one)))));
}
complex_double
catan( complex_double x )
{
complex_double imag;
complex_double ihalf;
imag.real = 0.0;
imag.imag = 1.0;
ihalf.real = 0.0;
ihalf.imag = 0.5;
return c_mult( ihalf, clog( c_div( c_add(imag, x), c_sub(imag,x ))));
}
complex_double
catan2( complex_double x, complex_double __y )
{
complex_double t;
abort();
return t;
}
complex_double
ccosh( complex_double x )
{
complex_double two;
complex_double nx;
two.real = 2.0;
two.imag = 0.0;
nx.real = -x.real;
nx.imag = -x.imag;
return c_div( c_add( cexp(x), cexp(nx)), two);
}
complex_double
cmod( complex_double x, complex_double __y )
{
complex_double t;
abort();
return t;
}
complex_double
csinh( complex_double x )
{
complex_double two;
complex_double nx;
two.real = 2.0;
two.imag = 0.0;
nx.real = -x.real;
nx.imag = -x.imag;
return c_div( c_sub( cexp(x), cexp(nx)), two);
}
complex_double
ctan( complex_double x )
{
return c_div(csin(x), ccos(x));
}
complex_double
ctanh( complex_double x )
{
return c_div(csinh(x), ccosh(x));
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.