/* Complex Data Type definition for S-Lang */
/* Copyright (c) 1997 John E. Davis
* All rights reserved.
*
* You may distribute under the terms of either the GNU General Public
* License or the Perl Artistic License.
*/
#include "config.h"
#include "sl-feat.h"
#include <stdio.h>
#include "slang.h"
#include "_slang.h"
/* The rest of the file is enclosed in this #if */
#if SLANG_HAS_COMPLEX
#if SLANG_HAS_FLOAT
# include <math.h>
#endif
#ifdef PI
# undef PI
#endif
#define PI 3.14159265358979323846
int SLang_pop_complex (double *r, double *i)
{
SLang_Object_Type obj;
double *c;
switch (SLang_peek_at_stack ())
{
default:
case SLANG_COMPLEX_TYPE:
if (-1 == _SLang_pop_object_of_type (SLANG_COMPLEX_TYPE, &obj))
return -1;
c = (double *)obj.v.p_val;
*r = c[0];
*i = c[1];
SLang_free_object (&obj);
break;
case SLANG_INT_TYPE:
case SLANG_DOUBLE_TYPE:
*i = 0.0;
if (-1 == SLang_pop_double (r, NULL, NULL))
return -1;
break;
case -1:
return -1;
}
return 0;
}
int SLang_push_complex (double r, double i)
{
double *c;
c = (double *) SLmalloc (2 * sizeof (double));
if (c == NULL)
return -1;
c[0] = r;
c[1] = i;
if (-1 == _SLang_push_void_star (SLANG_COMPLEX_TYPE, (VOID_STAR) c))
{
SLfree ((char *) c);
return -1;
}
return 0;
}
double *SLcomplex_times (double *c, double *a, double *b)
{
double a_real, b_real, a_imag, b_imag;
a_real = a[0];
b_real = b[0];
a_imag = a[1];
b_imag = b[1];
c[0] = a_real * b_real - a_imag * b_imag;
c[1] = a_imag * b_real + a_real * b_imag;
return c;
}
double *SLcomplex_divide (double *c, double *a, double *b)
{
double a_real, b_real, a_imag, b_imag;
double ratio, invden;
a_real = a[0];
b_real = b[0];
a_imag = a[1];
b_imag = b[1];
/* Do it this way to avoid overflow in the denom */
if (fabs(b_real) > fabs(b_imag))
{
ratio = b_imag / b_real;
invden = 1.0 / (b_real + b_imag * ratio);
c[0] = (a_real + ratio * a_imag) * invden;
c[1] = (a_imag - a_real * ratio) * invden;
}
else
{
ratio = b_real / b_imag;
invden = 1.0 / (b_real * ratio + b_imag);
c[0] = (a_real * ratio + a_imag) * invden;
c[1] = (a_imag * ratio - a_real) * invden;
}
return c;
}
/* a^b = exp (b log a); */
double *SLcomplex_pow (double *c, double *a, double *b)
{
return SLcomplex_exp (c, SLcomplex_times (c, b, SLcomplex_log (c, a)));
}
double SLcomplex_abs (double *z)
{
return SLmath_hypot (z[0], z[1]);
}
/* It appears that FORTRAN assumes that the branch cut for the log function
* is along the -x axis. So, use this for atan2:
*/
static double my_atan2 (double y, double x)
{
double val;
val = atan (y/x);
if (x >= 0)
return val; /* I, IV */
if (y <= 0) /* III */
return val - PI;
return PI + val; /* II */
}
static void polar_form (double *r, double *theta, double *z)
{
double x, y;
*r = SLcomplex_abs (z);
x = z[0];
y = z[1];
if (x == 0.0)
{
if (y >= 0)
*theta = 0.5 * PI;
else
*theta = 1.5 * PI;
}
else *theta = my_atan2 (y, x);
}
double *SLcomplex_sin (double *sinz, double *z)
{
double x, y;
x = z[0]; y = z[1];
sinz[0] = sin (x) * cosh (y);
sinz[1] = cos (x) * sinh (y);
return sinz;
}
double *SLcomplex_cos (double *cosz, double *z)
{
double x, y;
x = z[0]; y = z[1];
cosz[0] = cos (x) * cosh (y);
cosz[1] = -sin (x) * sinh (y);
return cosz;
}
double *SLcomplex_exp (double *expz, double *z)
{
double r, i;
r = exp (z[0]);
i = z[1];
expz[0] = r * cos (i);
expz[1] = r * sin (i);
return expz;
}
double *SLcomplex_log (double *logz, double *z)
{
double r, theta;
polar_form (&r, &theta, z); /* log R.e^(ix) = log R + ix */
logz[0] = log(r);
logz[1] = theta;
return logz;
}
double *SLcomplex_log10 (double *log10z, double *z)
{
double l10 = log (10);
(void) SLcomplex_log (log10z, z);
log10z[0] = log10z[0] / l10;
log10z[1] = log10z[1] / l10;
return log10z;
}
double *SLcomplex_sqrt (double *sqrtz, double *z)
{
double theta, r;
polar_form (&r, &theta, z);
r = sqrt (r);
theta = 0.5 * theta;
sqrtz[0] = r * cos (theta);
sqrtz[1] = r * sin (theta);
return sqrtz;
}
double *SLcomplex_tan (double *tanz, double *z)
{
double x, y, invden;
x = 2 * z[0];
y = 2 * z[1];
invden = 1.0 / (cos (x) + cosh (y));
tanz[0] = invden * sin (x);
tanz[1] = invden * sinh (y);
return tanz;
}
/* Utility Function */
static void compute_alpha_beta (double *z, double *alpha, double *beta)
{
double x, y, a, b;
x = z[0];
y = z[1];
a = 0.5 * SLmath_hypot (x + 1, y);
b = 0.5 * SLmath_hypot (x - 1, y);
*alpha = a + b;
*beta = a - b;
}
double *SLcomplex_asin (double *asinz, double *z)
{
double alpha, beta;
compute_alpha_beta (z, &alpha, &beta);
asinz[0] = asin (beta);
asinz[1] = log (alpha + sqrt (alpha * alpha - 1));
return asinz;
}
double *SLcomplex_acos (double *acosz, double *z)
{
double alpha, beta;
compute_alpha_beta (z, &alpha, &beta);
acosz[0] = acos (beta);
acosz[1] = -log (alpha + sqrt (alpha * alpha - 1));
return acosz;
}
double *SLcomplex_atan (double *atanz, double *z)
{
double x, y;
double z1[2], z2[2];
x = z[0]; y = z[1];
z1[0] = x;
z1[1] = 1 + y;
z2[0] = -x;
z2[1] = 1 - y;
SLcomplex_log (z1, SLcomplex_divide (z2, z1, z2));
atanz[0] = -0.5 * z1[1];
atanz[1] = 0.5 * z1[0];
return atanz;
}
double *SLcomplex_sinh (double *sinhz, double *z)
{
double x, y;
x = z[0]; y = z[1];
sinhz[0] = sinh (x) * cos (y);
sinhz[1] = cosh (x) * sin (y);
return sinhz;
}
double *SLcomplex_cosh (double *coshz, double *z)
{
double x, y;
x = z[0]; y = z[1];
coshz[0] = cosh (x) * cos (y);
coshz[1] = sinh (x) * sin (y);
return coshz;
}
double *SLcomplex_tanh (double *tanhz, double *z)
{
double x, y, invden;
x = 2 * z[0];
y = 2 * z[1];
invden = 1.0 / (cosh (x) + cos (y));
tanhz[0] = invden * sinh (x);
tanhz[1] = invden * sin (y);
return tanhz;
}
double *SLcomplex_asinh (double *asinhz, double *z)
{
return asinhz;
}
double *SLcomplex_acosh (double *acoshz, double *z)
{
return acoshz;
}
double *SLcomplex_atanh (double *atanhz, double *z)
{
return atanhz;
}
static int complex_binary_result (int op, unsigned char a, unsigned char b,
unsigned char *c)
{
(void) a; (void) b;
switch (op)
{
default:
case SLANG_POW:
case SLANG_PLUS:
case SLANG_MINUS:
case SLANG_TIMES:
case SLANG_DIVIDE:
*c = SLANG_COMPLEX_TYPE;
break;
case SLANG_EQ:
case SLANG_NE:
*c = SLANG_INT_TYPE;
break;
}
return 1;
}
static int complex_complex_binary (int op,
unsigned char a_type, VOID_STAR ap, unsigned int na,
unsigned char b_type, VOID_STAR bp, unsigned int nb,
VOID_STAR cp)
{
int *ic;
double *a, *b, *c;
unsigned int n, n_max;
unsigned int da, db;
(void) a_type;
(void) b_type;
a = (double *) ap;
b = (double *) bp;
c = (double *) cp;
ic = (int *) cp;
if (na == 1) da = 0; else da = 2;
if (nb == 1) db = 0; else db = 2;
if (na > nb) n_max = na; else n_max = nb;
n_max = 2 * n_max;
switch (op)
{
default:
case SLANG_PLUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] + b[0];
c[n + 1] = a[1] + b[1];
a += da; b += db;
}
break;
case SLANG_MINUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] - b[0];
c[n + 1] = a[1] - b[1];
a += da; b += db;
}
break;
case SLANG_TIMES:
for (n = 0; n < n_max; n += 2)
{
SLcomplex_times (c + n, a, b);
a += da; b += db;
}
break;
case SLANG_DIVIDE: /* / */
for (n = 0; n < n_max; n += 2)
{
if ((b[0] == 0.0) && (b[1] == 0.0))
{
SLang_Error = SL_DIVIDE_ERROR;
return -1;
}
SLcomplex_divide (c + n, a, b);
a += da; b += db;
}
break;
case SLANG_EQ: /* == */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = ((a[0] == b[0]) && (a[1] == b[1]));
a += da; b += db;
}
break;
case SLANG_NE: /* != */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = ((a[0] != b[0]) || (a[1] != b[1]));
a += da; b += db;
}
break;
case SLANG_POW:
for (n = 0; n < n_max; n += 2)
{
SLcomplex_pow (c + n, a, b);
a += da; b += db;
}
break;
}
return 1;
}
static int complex_double_binary (int op,
unsigned char a_type, VOID_STAR ap, unsigned int na,
unsigned char b_type, VOID_STAR bp, unsigned int nb,
VOID_STAR cp)
{
int *ic;
double *a, *b, *c;
unsigned int n, n_max;
unsigned int da, db;
(void) a_type;
(void) b_type;
a = (double *) ap;
b = (double *) bp;
c = (double *) cp;
ic = (int *) cp;
if (na == 1) da = 0; else da = 2;
if (nb == 1) db = 0; else db = 1;
if (na > nb) n_max = na; else n_max = nb;
n_max = 2 * n_max;
switch (op)
{
default:
case SLANG_POW:
return 0;
case SLANG_PLUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] + b[0];
c[n + 1] = a[1];
a += da; b += db;
}
break;
case SLANG_MINUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] - b[0];
c[n + 1] = a[1];
a += da; b += db;
}
break;
case SLANG_TIMES:
for (n = 0; n < n_max; n += 2)
{
double b0 = b[0];
c[n] = a[0] * b0;
c[n + 1] = a[1] * b0;
a += da; b += db;
}
break;
case SLANG_DIVIDE: /* / */
for (n = 0; n < n_max; n += 2)
{
double b0 = b[0];
if (b0 == 0.0)
{
SLang_Error = SL_DIVIDE_ERROR;
return -1;
}
c[n] = a[0] / b0;
c[n + 1] = a[1] / b0;
a += da; b += db;
}
break;
case SLANG_EQ: /* == */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = ((a[0] == b[0]) && (a[1] == 0.0));
a += da; b += db;
}
break;
case SLANG_NE: /* != */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = ((a[0] != b[0]) || (a[1] != 0.0));
a += da; b += db;
}
break;
}
return 1;
}
static int double_complex_binary (int op,
unsigned char a_type, VOID_STAR ap, unsigned int na,
unsigned char b_type, VOID_STAR bp, unsigned int nb,
VOID_STAR cp)
{
int *ic;
double *a, *b, *c;
unsigned int n, n_max;
unsigned int da, db;
(void) a_type;
(void) b_type;
a = (double *) ap;
b = (double *) bp;
c = (double *) cp;
ic = (int *) cp;
if (na == 1) da = 0; else da = 1;
if (nb == 1) db = 0; else db = 2;
if (na > nb) n_max = na; else n_max = nb;
n_max = 2 * n_max;
switch (op)
{
default:
case SLANG_POW:
return 0;
case SLANG_PLUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] + b[0];
c[n + 1] = b[1];
a += da; b += db;
}
break;
case SLANG_MINUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] - b[0];
c[n + 1] = -b[1];
a += da; b += db;
}
break;
case SLANG_TIMES:
for (n = 0; n < n_max; n += 2)
{
double a0 = a[0];
c[n] = a0 * b[0];
c[n + 1] = a0 * b[1];
a += da; b += db;
}
break;
case SLANG_DIVIDE: /* / */
for (n = 0; n < n_max; n += 2)
{
double z[2];
if ((b[0] == 0.0) && (b[1] == 0.0))
{
SLang_Error = SL_DIVIDE_ERROR;
return -1;
}
z[0] = a[0];
z[1] = 0.0;
SLcomplex_divide (c + n, z, b);
a += da; b += db;
}
break;
case SLANG_EQ: /* == */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = ((a[0] == b[0]) && (0.0 == b[1]));
a += da; b += db;
}
break;
case SLANG_NE: /* != */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = ((a[0] != b[0]) || (0.0 != b[1]));
a += da; b += db;
}
break;
}
return 1;
}
static int complex_int_binary (int op,
unsigned char a_type, VOID_STAR ap, unsigned int na,
unsigned char b_type, VOID_STAR bp, unsigned int nb,
VOID_STAR cp)
{
int *ic, *b;
double *a, *c;
unsigned int n, n_max;
unsigned int da, db;
(void) a_type;
(void) b_type;
a = (double *) ap;
b = (int *) bp;
c = (double *) cp;
ic = (int *) cp;
if (na == 1) da = 0; else da = 2;
if (nb == 1) db = 0; else db = 1;
if (na > nb) n_max = na; else n_max = nb;
n_max = 2 * n_max;
switch (op)
{
default:
case SLANG_POW:
return 0;
case SLANG_PLUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] + b[0];
c[n + 1] = a[1];
a += da; b += db;
}
break;
case SLANG_MINUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] - b[0];
c[n + 1] = a[1];
a += da; b += db;
}
break;
case SLANG_TIMES:
for (n = 0; n < n_max; n += 2)
{
double b0 = b[0];
c[n] = a[0] * b0;
c[n + 1] = a[1] * b0;
a += da; b += db;
}
break;
case SLANG_DIVIDE: /* / */
for (n = 0; n < n_max; n += 2)
{
int b0 = b[0];
if (b0 == 0)
{
SLang_Error = SL_DIVIDE_ERROR;
return -1;
}
c[n] = a[0] / b0;
c[n + 1] = a[1] / b0;
a += da; b += db;
}
break;
case SLANG_EQ: /* == */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = ((a[0] == (double) b[0]) && (a[1] == 0.0));
a += da; b += db;
}
break;
case SLANG_NE: /* != */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = ((a[0] != (double) b[0]) || (a[1] != 0.0));
a += da; b += db;
}
break;
}
return 1;
}
static int int_complex_binary (int op,
unsigned char a_type, VOID_STAR ap, unsigned int na,
unsigned char b_type, VOID_STAR bp, unsigned int nb,
VOID_STAR cp)
{
int *ic, *a;
double *b, *c;
unsigned int n, n_max;
unsigned int da, db;
(void) a_type;
(void) b_type;
a = (int *) ap;
b = (double *) bp;
c = (double *) cp;
ic = (int *) cp;
if (na == 1) da = 0; else da = 1;
if (nb == 1) db = 0; else db = 2;
if (na > nb) n_max = na; else n_max = nb;
n_max = 2 * n_max;
switch (op)
{
default:
case SLANG_POW:
return 0;
case SLANG_PLUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] + b[0];
c[n + 1] = b[1];
a += da; b += db;
}
break;
case SLANG_MINUS:
for (n = 0; n < n_max; n += 2)
{
c[n] = a[0] - b[0];
c[n + 1] = -b[1];
a += da; b += db;
}
break;
case SLANG_TIMES:
for (n = 0; n < n_max; n += 2)
{
int a0 = a[0];
c[n] = a0 * b[0];
c[n + 1] = a0 * b[1];
a += da; b += db;
}
break;
case SLANG_DIVIDE: /* / */
for (n = 0; n < n_max; n += 2)
{
double z[2];
if ((b[0] == 0.0) && (b[1] == 0.0))
{
SLang_Error = SL_DIVIDE_ERROR;
return -1;
}
z[0] = (double) a[0];
z[1] = 0.0;
SLcomplex_divide (c + n, z, b);
a += da; b += db;
}
break;
case SLANG_EQ: /* == */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = (((double)a[0] == b[0]) && (0.0 == b[1]));
a += da; b += db;
}
break;
case SLANG_NE: /* != */
for (n = 0; n < n_max; n += 2)
{
ic[n/2] = (((double)a[0] != b[0]) || (0.0 != b[1]));
a += da; b += db;
}
break;
}
return 1;
}
static int complex_unary_result (int op, unsigned char a, unsigned char *b)
{
(void) a;
switch (op)
{
default:
return 0;
case SLANG_PLUSPLUS:
case SLANG_MINUSMINUS:
case SLANG_CHS:
case SLANG_MUL2:
*b = SLANG_COMPLEX_TYPE;
break;
case SLANG_SQR: /* |Real|^2 + |Imag|^2 ==> double */
case SLANG_ABS: /* |z| ==> double */
*b = SLANG_DOUBLE_TYPE;
break;
case SLANG_SIGN:
*b = SLANG_INT_TYPE;
break;
}
return 1;
}
static int complex_unary (int op,
unsigned char a_type, VOID_STAR ap, unsigned int na,
VOID_STAR bp)
{
unsigned int n;
double *a, *b;
int *ic;
(void) a_type;
a = (double *) ap;
b = (double *) bp;
ic = (int *) bp;
na = 2 * na;
switch (op)
{
default:
return 0;
case SLANG_PLUSPLUS:
for (n = 0; n < na; n += 2) b[n] = (a[n] + 1);
break;
case SLANG_MINUSMINUS:
for (n = 0; n < na; n += 2) b[n] = (a[n] - 1);
break;
case SLANG_CHS:
for (n = 0; n < na; n += 2)
{
b[n] = -(a[n]);
b[n + 1] = -(a[n + 1]);
}
break;
case SLANG_SQR: /* |Real|^2 + |Imag|^2 ==> double */
for (n = 0; n < na; n += 2)
b[n/2] = (a[n] * a[n] + a[n + 1] * a[n + 1]);
break;
case SLANG_MUL2:
for (n = 0; n < na; n += 2)
{
b[n] = (2 * a[n]);
b[n + 1] = (2 * a[n + 1]);
}
break;
case SLANG_ABS: /* |z| ==> double */
for (n = 0; n < na; n += 2)
b[n/2] = SLcomplex_abs (a + n);
break;
case SLANG_SIGN:
/* Another creative extension. Lets return an integer which indicates
* whether the complex number is in the upperhalf plane or not.
*/
for (n = 0; n < na; n += 2)
{
if (a[n + 1] < 0.0) ic[n/2] = -1;
else if (a[n + 1] > 0.0) ic[n/2] = 1;
else ic[n/2] = 0;
}
break;
}
return 1;
}
static int
complex_typecast (unsigned char from_type, VOID_STAR from, unsigned int num,
unsigned char to_type, VOID_STAR to)
{
double *z;
double *d;
int *i;
unsigned int n;
(void) to_type;
z = (double *) to;
switch (from_type)
{
default:
return 0;
case SLANG_INT_TYPE:
i = (int *) from;
for (n = 0; n < num; n++)
{
*z++ = (double) i[n];
*z++ = 0.0;
}
break;
case SLANG_DOUBLE_TYPE:
d = (double *) from;
for (n = 0; n < num; n++)
{
*z++ = d[n];
*z++ = 0.0;
}
break;
}
return 1;
}
static void complex_destroy (unsigned char type, VOID_STAR ptr)
{
(void) type;
SLfree ((char *)*(double **) ptr);
}
static int complex_push (unsigned char type, VOID_STAR ptr)
{
double *z;
(void) type;
z = *(double **) ptr;
return SLang_push_complex (z[0], z[1]);
}
static int complex_pop (unsigned char type, VOID_STAR ptr)
{
double *z;
(void) type;
z = *(double **) ptr;
return SLang_pop_complex (&z[0], &z[1]);
}
int _SLinit_slcomplex (void)
{
SLang_Class_Type *cl;
if (NULL == (cl = SLclass_allocate_class ("Complex_Type")))
return -1;
(void) SLclass_set_destroy_function (cl, complex_destroy);
(void) SLclass_set_push_function (cl, complex_push);
(void) SLclass_set_pop_function (cl, complex_pop);
if (-1 == SLclass_register_class (cl, SLANG_COMPLEX_TYPE, 2 * sizeof (double),
SLANG_CLASS_TYPE_VECTOR))
return -1;
if ((-1 == (SLclass_add_binary_op (SLANG_COMPLEX_TYPE, SLANG_COMPLEX_TYPE, complex_complex_binary, complex_binary_result)))
|| (-1 == (SLclass_add_binary_op (SLANG_COMPLEX_TYPE, SLANG_DOUBLE_TYPE, complex_double_binary, complex_binary_result)))
|| (-1 == (SLclass_add_binary_op (SLANG_DOUBLE_TYPE, SLANG_COMPLEX_TYPE, double_complex_binary, complex_binary_result)))
|| (-1 == (SLclass_add_binary_op (SLANG_COMPLEX_TYPE, SLANG_INT_TYPE, complex_int_binary, complex_binary_result)))
|| (-1 == (SLclass_add_binary_op (SLANG_INT_TYPE, SLANG_COMPLEX_TYPE, int_complex_binary, complex_binary_result)))
|| (-1 == (SLclass_add_unary_op (SLANG_COMPLEX_TYPE, complex_unary, complex_unary_result)))
|| (-1 == (SLclass_add_typecast (SLANG_DOUBLE_TYPE, SLANG_COMPLEX_TYPE, complex_typecast, 1)))
|| (-1 == (SLclass_add_typecast (SLANG_INT_TYPE, SLANG_COMPLEX_TYPE, complex_typecast, 1))))
return -1;
return 0;
}
#endif /* if SLANG_HAS_COMPLEX */
