491 lines
13 KiB
C
491 lines
13 KiB
C
/**********
|
|
Copyright 1990 Regents of the University of California. All rights reserved.
|
|
Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group
|
|
**********/
|
|
|
|
/** \file cmath3.c
|
|
\brief functions for the control language parser
|
|
|
|
Routines to do complex mathematical functions. These routines require
|
|
the -lm libraries. We sacrifice a lot of space to be able
|
|
to avoid having to do a seperate call for every vector element,
|
|
but it pays off in time savings. These routines should never
|
|
allow FPE's to happen.
|
|
|
|
Complex functions are called as follows:
|
|
cx_something(data, type, length, &newlength, &newtype),
|
|
and return a char * that is cast to complex or double.
|
|
*/
|
|
|
|
|
|
#include "ngspice/ngspice.h"
|
|
#include "ngspice/cpdefs.h"
|
|
#include "ngspice/dvec.h"
|
|
|
|
#include "cmath.h"
|
|
#include "cmath3.h"
|
|
|
|
|
|
static ngcomplex_t *cexp_sp3(ngcomplex_t *c); /* cexp exist's in some newer compiler */
|
|
static int cln(ngcomplex_t *c, ngcomplex_t *rv);
|
|
static void ctimes(ngcomplex_t *c1, ngcomplex_t *c2, ngcomplex_t *rv);
|
|
|
|
void *cx_divide(void *data1, void *data2,
|
|
short int datatype1, short int datatype2, int length)
|
|
{
|
|
int xrc = 0;
|
|
void *rv;
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t *c, c1, c2;
|
|
int i;
|
|
|
|
if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
|
|
double *d;
|
|
rv = d = alloc_d(length);
|
|
for (i = 0; i < length; i++) {
|
|
rcheck(dd2[i] != 0, "divide");
|
|
d[i] = dd1[i] / dd2[i];
|
|
}
|
|
}
|
|
else {
|
|
rv = c = alloc_c(length);
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
}
|
|
else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
}
|
|
else {
|
|
c2 = cc2[i];
|
|
}
|
|
rcheck((realpart(c2) != 0) || (imagpart(c2) != 0), "divide");
|
|
#define xx5 realpart(c1)
|
|
#define xx6 imagpart(c1)
|
|
cdiv(xx5, xx6, realpart(c2), imagpart(c2), realpart(c[i]),
|
|
imagpart(c[i]));
|
|
}
|
|
}
|
|
|
|
EXITPOINT:
|
|
if (xrc != 0) { /* Free resources on error */
|
|
tfree(rv);
|
|
rv = NULL;
|
|
}
|
|
|
|
return rv;
|
|
} /* end of function cx_divide */
|
|
|
|
|
|
|
|
/* Should just use "j( )" */
|
|
/* The comma operator. What this does (unless it is part of the argument
|
|
* list of a user-defined function) is arg1 + j(arg2).
|
|
*/
|
|
|
|
void *
|
|
cx_comma(void *data1, void *data2, short int datatype1, short int datatype2, int length)
|
|
{
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t *c, c1, c2;
|
|
int i;
|
|
|
|
c = alloc_c(length);
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
} else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
} else {
|
|
c2 = cc2[i];
|
|
}
|
|
|
|
realpart(c[i]) = realpart(c1) + imagpart(c2);
|
|
imagpart(c[i]) = imagpart(c1) + realpart(c2);
|
|
}
|
|
return ((void *) c);
|
|
}
|
|
|
|
void *cx_power(void *data1, void *data2,
|
|
short int datatype1, short int datatype2, int length)
|
|
{
|
|
int xrc = 0;
|
|
void *rv;
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
|
|
if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
|
|
double *d;
|
|
rv = d = alloc_d(length);
|
|
|
|
int i;
|
|
for (i = 0; i < length; i++) {
|
|
rcheck((dd1[i] >= 0) || (floor(dd2[i]) == ceil(dd2[i])), "power");
|
|
d[i] = pow(dd1[i], dd2[i]);
|
|
}
|
|
}
|
|
else {
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t *c, c1, c2, *t;
|
|
rv = c = alloc_c(length);
|
|
|
|
int i;
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
}
|
|
else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
}
|
|
else {
|
|
c2 = cc2[i];
|
|
}
|
|
|
|
if ((realpart(c1) == 0.0) && (imagpart(c1) == 0.0)) {
|
|
realpart(c[i]) = 0.0;
|
|
imagpart(c[i]) = 0.0;
|
|
}
|
|
else { /* if ((imagpart(c1) != 0.0) &&
|
|
(imagpart(c2) != 0.0)) */
|
|
ngcomplex_t tmp, tmp2;
|
|
if (cln(&c1, &tmp) != 0) {
|
|
(void) fprintf(cp_err, "power of 0 + i 0 not allowed.\n");
|
|
xrc = -1;
|
|
goto EXITPOINT;
|
|
}
|
|
ctimes(&c2, &tmp, &tmp2);
|
|
t = cexp_sp3(&tmp2);
|
|
c[i] = *t;
|
|
/*
|
|
} else {
|
|
realpart(c[i]) = pow(realpart(c1),
|
|
realpart(c2));
|
|
imagpart(c[i]) = 0.0;
|
|
*/
|
|
}
|
|
}
|
|
}
|
|
|
|
EXITPOINT:
|
|
if (xrc != 0) { /* Free resources on error */
|
|
txfree(rv);
|
|
rv = NULL;
|
|
}
|
|
|
|
return rv;
|
|
} /* end of function cx_power */
|
|
|
|
|
|
|
|
/* These are unnecessary... Only cx_power uses them... */
|
|
|
|
static ngcomplex_t *
|
|
cexp_sp3(ngcomplex_t *c)
|
|
{
|
|
static ngcomplex_t r;
|
|
double d;
|
|
|
|
d = exp(realpart(*c));
|
|
realpart(r) = d * cos(imagpart(*c));
|
|
if (imagpart(*c) != 0.0)
|
|
imagpart(r) = d * sin(imagpart(*c));
|
|
else
|
|
imagpart(r) = 0.0;
|
|
return (&r);
|
|
}
|
|
|
|
static int cln(ngcomplex_t *c, ngcomplex_t *rv)
|
|
{
|
|
double c_r = c->cx_real;
|
|
double c_i = c->cx_imag;
|
|
|
|
if (c_r == 0 && c_i == 0) {
|
|
(void) fprintf(cp_err, "Complex log of 0 + i0 is undefined.\n");
|
|
return -1;
|
|
}
|
|
|
|
rv->cx_real = log(cmag(*c));
|
|
if (c_i != 0.0) {
|
|
rv->cx_imag = atan2(c_i, c_r);
|
|
}
|
|
else {
|
|
rv->cx_imag = 0.0;
|
|
}
|
|
return 0;
|
|
} /* end of functon cln */
|
|
|
|
|
|
|
|
static void ctimes(ngcomplex_t *c1, ngcomplex_t *c2, ngcomplex_t *rv)
|
|
{
|
|
rv->cx_real = realpart(*c1) * realpart(*c2) -
|
|
imagpart(*c1) * imagpart(*c2);
|
|
rv->cx_imag = imagpart(*c1) * realpart(*c2) +
|
|
realpart(*c1) * imagpart(*c2);
|
|
return;
|
|
}
|
|
|
|
|
|
|
|
/* Now come all the relational and logical functions. It's overkill to put
|
|
* them here, but... Note that they always return a real value, with the
|
|
* result the same length as the arguments.
|
|
*/
|
|
|
|
void *
|
|
cx_eq(void *data1, void *data2, short int datatype1, short int datatype2, int length)
|
|
{
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
double *d;
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t c1, c2;
|
|
int i;
|
|
|
|
d = alloc_d(length);
|
|
if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
|
|
for (i = 0; i < length; i++)
|
|
if (dd1[i] == dd2[i])
|
|
d[i] = 1.0;
|
|
else
|
|
d[i] = 0.0;
|
|
} else {
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
} else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
} else {
|
|
c2 = cc2[i];
|
|
}
|
|
d[i] = ((realpart(c1) == realpart(c2)) &&
|
|
(imagpart(c1) == imagpart(c2)));
|
|
}
|
|
}
|
|
return ((void *) d);
|
|
}
|
|
|
|
void *
|
|
cx_gt(void *data1, void *data2, short int datatype1, short int datatype2, int length)
|
|
{
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
double *d;
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t c1, c2;
|
|
int i;
|
|
|
|
d = alloc_d(length);
|
|
if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
|
|
for (i = 0; i < length; i++)
|
|
if (dd1[i] > dd2[i])
|
|
d[i] = 1.0;
|
|
else
|
|
d[i] = 0.0;
|
|
} else {
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
} else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
} else {
|
|
c2 = cc2[i];
|
|
}
|
|
d[i] = ((realpart(c1) > realpart(c2)) &&
|
|
(imagpart(c1) > imagpart(c2)));
|
|
}
|
|
}
|
|
return ((void *) d);
|
|
}
|
|
|
|
void *
|
|
cx_lt(void *data1, void *data2, short int datatype1, short int datatype2, int length)
|
|
{
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
double *d;
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t c1, c2;
|
|
int i;
|
|
|
|
d = alloc_d(length);
|
|
if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
|
|
for (i = 0; i < length; i++)
|
|
if (dd1[i] < dd2[i])
|
|
d[i] = 1.0;
|
|
else
|
|
d[i] = 0.0;
|
|
} else {
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
} else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
} else {
|
|
c2 = cc2[i];
|
|
}
|
|
d[i] = ((realpart(c1) < realpart(c2)) &&
|
|
(imagpart(c1) < imagpart(c2)));
|
|
}
|
|
}
|
|
return ((void *) d);
|
|
}
|
|
|
|
void *
|
|
cx_ge(void *data1, void *data2, short int datatype1, short int datatype2, int length)
|
|
{
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
double *d;
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t c1, c2;
|
|
int i;
|
|
|
|
d = alloc_d(length);
|
|
if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
|
|
for (i = 0; i < length; i++)
|
|
if (dd1[i] >= dd2[i])
|
|
d[i] = 1.0;
|
|
else
|
|
d[i] = 0.0;
|
|
} else {
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
} else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
} else {
|
|
c2 = cc2[i];
|
|
}
|
|
d[i] = ((realpart(c1) >= realpart(c2)) &&
|
|
(imagpart(c1) >= imagpart(c2)));
|
|
}
|
|
}
|
|
return ((void *) d);
|
|
}
|
|
|
|
void *
|
|
cx_le(void *data1, void *data2, short int datatype1, short int datatype2, int length)
|
|
{
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
double *d;
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t c1, c2;
|
|
int i;
|
|
|
|
d = alloc_d(length);
|
|
if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
|
|
for (i = 0; i < length; i++)
|
|
if (dd1[i] <= dd2[i])
|
|
d[i] = 1.0;
|
|
else
|
|
d[i] = 0.0;
|
|
} else {
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
} else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
} else {
|
|
c2 = cc2[i];
|
|
}
|
|
d[i] = ((realpart(c1) <= realpart(c2)) &&
|
|
(imagpart(c1) <= imagpart(c2)));
|
|
}
|
|
}
|
|
return ((void *) d);
|
|
}
|
|
|
|
void *
|
|
cx_ne(void *data1, void *data2, short int datatype1, short int datatype2, int length)
|
|
{
|
|
double *dd1 = (double *) data1;
|
|
double *dd2 = (double *) data2;
|
|
double *d;
|
|
ngcomplex_t *cc1 = (ngcomplex_t *) data1;
|
|
ngcomplex_t *cc2 = (ngcomplex_t *) data2;
|
|
ngcomplex_t c1, c2;
|
|
int i;
|
|
|
|
d = alloc_d(length);
|
|
if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
|
|
for (i = 0; i < length; i++)
|
|
if (dd1[i] != dd2[i])
|
|
d[i] = 1.0;
|
|
else
|
|
d[i] = 0.0;
|
|
} else {
|
|
for (i = 0; i < length; i++) {
|
|
if (datatype1 == VF_REAL) {
|
|
realpart(c1) = dd1[i];
|
|
imagpart(c1) = 0.0;
|
|
} else {
|
|
c1 = cc1[i];
|
|
}
|
|
if (datatype2 == VF_REAL) {
|
|
realpart(c2) = dd2[i];
|
|
imagpart(c2) = 0.0;
|
|
} else {
|
|
c2 = cc2[i];
|
|
}
|
|
d[i] = ((realpart(c1) != realpart(c2)) &&
|
|
(imagpart(c1) != imagpart(c2)));
|
|
}
|
|
}
|
|
return ((void *) d);
|
|
}
|
|
|