diff --git a/src/include/complex.h b/src/include/complex.h index 3e7ae1d6b..3a083d4f0 100644 --- a/src/include/complex.h +++ b/src/include/complex.h @@ -47,24 +47,24 @@ typedef struct { * COMPLEX NUMBER DATA STRUCTURE * * >>> Structure fields: - * Real (RealNumber) - * The real portion of the number. Real must be the first + * real (realNumber) + * The real portion of the number. real must be the first * field in this structure. - * Imag (RealNumber) + * imag (realNumber) * The imaginary portion of the number. This field must follow - * immediately after Real. + * immediately after real. */ #define spREAL double -/* Begin `RealNumber'. */ -typedef spREAL RealNumber, *RealVector; +// /* Begin `realNumber'. */ +// typedef spREAL realNumber, *realVector; -/* Begin `ComplexNumber'. */ -typedef struct -{ RealNumber Real; - RealNumber Imag; -} ComplexNumber, *ComplexVector; +// /* Begin `ComplexNumber'. */ +// typedef struct +// { RealNumber Real; +// RealNumber Imag; +// } ComplexNumber, *ComplexVector; /* Some defines used mainly in cmath.c. */ @@ -283,206 +283,208 @@ typedef struct (A).imag = (B.imag) - (C.imag); \ } +#ifndef _SPDEFS_H + /* Macro function that returns the approx absolute value of a complex number. */ -#define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag)) +#define ELEMENT_MAG(ptr) (ABS((ptr)->real) + ABS((ptr)->imag)) -#define CMPLX_ASSIGN_VALUE(cnum, vReal, vImag) \ -{ (cnum).Real = vReal; \ - (cnum).Imag = vImag; \ +#define CMPLX_ASSIGN_VALUE(cnum, vreal, vimag) \ +{ (cnum).real = vreal; \ + (cnum).imag = vimag; \ } /* Complex assignment statements. */ #define CMPLX_ASSIGN(to,from) \ -{ (to).Real = (from).Real; \ - (to).Imag = (from).Imag; \ +{ (to).real = (from).real; \ + (to).imag = (from).imag; \ } #define CMPLX_CONJ_ASSIGN(to,from) \ -{ (to).Real = (from).Real; \ - (to).Imag = -(from).Imag; \ +{ (to).real = (from).real; \ + (to).imag = -(from).imag; \ } #define CMPLX_NEGATE_ASSIGN(to,from) \ -{ (to).Real = -(from).Real; \ - (to).Imag = -(from).Imag; \ +{ (to).real = -(from).real; \ + (to).imag = -(from).imag; \ } #define CMPLX_CONJ_NEGATE_ASSIGN(to,from) \ -{ (to).Real = -(from).Real; \ - (to).Imag = (from).Imag; \ +{ (to).real = -(from).real; \ + (to).imag = (from).imag; \ } -#define CMPLX_CONJ(a) (a).Imag = -(a).Imag +#define CMPLX_CONJ(a) (a).imag = -(a).imag -#define CONJUGATE(a) (a).Imag = -(a).Imag +#define CONJUGATE(a) (a).imag = -(a).imag #define CMPLX_NEGATE(a) \ -{ (a).Real = -(a).Real; \ - (a).Imag = -(a).Imag; \ +{ (a).real = -(a).real; \ + (a).imag = -(a).imag; \ } #define CMPLX_NEGATE_SELF(cnum) \ -{ (cnum).real = -(cnum).Real; \ - (cnum).imag = -(cnum).Imag; \ +{ (cnum).real = -(cnum).real; \ + (cnum).imag = -(cnum).imag; \ } /* Macro that returns the approx magnitude (L-1 norm) of a complex number. */ -#define CMPLX_1_NORM(a) (ABS((a).Real) + ABS((a).Imag)) +#define CMPLX_1_NORM(a) (ABS((a).real) + ABS((a).imag)) /* Macro that returns the approx magnitude (L-infinity norm) of a complex. */ -#define CMPLX_INF_NORM(a) (MAX (ABS((a).Real),ABS((a).Imag))) +#define CMPLX_INF_NORM(a) (MAX (ABS((a).real),ABS((a).imag))) /* Macro function that returns the magnitude (L-2 norm) of a complex number. */ -#define CMPLX_2_NORM(a) (sqrt((a).Real*(a).Real + (a).Imag*(a).Imag)) +#define CMPLX_2_NORM(a) (sqrt((a).real*(a).real + (a).imag*(a).imag)) /* Macro function that performs complex addition. */ #define CMPLX_ADD(to,from_a,from_b) \ -{ (to).Real = (from_a).Real + (from_b).Real; \ - (to).Imag = (from_a).Imag + (from_b).Imag; \ +{ (to).real = (from_a).real + (from_b).real; \ + (to).imag = (from_a).imag + (from_b).imag; \ } /* Macro function that performs addition of a complex and a scalar. */ #define CMPLX_ADD_SELF_SCALAR(cnum, scalar) \ -{ (cnum).Real += scalar; \ +{ (cnum).real += scalar; \ } /* Macro function that performs complex subtraction. */ #define CMPLX_SUBT(to,from_a,from_b) \ -{ (to).Real = (from_a).Real - (from_b).Real; \ - (to).Imag = (from_a).Imag - (from_b).Imag; \ +{ (to).real = (from_a).real - (from_b).real; \ + (to).imag = (from_a).imag - (from_b).imag; \ } /* Macro function that is equivalent to += operator for complex numbers. */ #define CMPLX_ADD_ASSIGN(to,from) \ -{ (to).Real += (from).Real; \ - (to).Imag += (from).Imag; \ +{ (to).real += (from).real; \ + (to).imag += (from).imag; \ } /* Macro function that is equivalent to -= operator for complex numbers. */ #define CMPLX_SUBT_ASSIGN(to,from) \ -{ (to).Real -= (from).Real; \ - (to).Imag -= (from).Imag; \ +{ (to).real -= (from).real; \ + (to).imag -= (from).imag; \ } /* Macro function that multiplies a complex number by a scalar. */ #define SCLR_MULT(to,sclr,cmplx) \ -{ (to).Real = (sclr) * (cmplx).Real; \ - (to).Imag = (sclr) * (cmplx).Imag; \ +{ (to).real = (sclr) * (cmplx).real; \ + (to).imag = (sclr) * (cmplx).imag; \ } /* Macro function that multiply-assigns a complex number by a scalar. */ #define SCLR_MULT_ASSIGN(to,sclr) \ -{ (to).Real *= (sclr); \ - (to).Imag *= (sclr); \ +{ (to).real *= (sclr); \ + (to).imag *= (sclr); \ } /* Macro function that multiplies two complex numbers. */ #define CMPLX_MULT(to,from_a,from_b) \ -{ (to).real = (from_a).Real * (from_b).Real - \ - (from_a).Imag * (from_b).Imag; \ - (to).imag = (from_a).Real * (from_b).Imag + \ - (from_a).Imag * (from_b).Real; \ +{ (to).real = (from_a).real * (from_b).real - \ + (from_a).imag * (from_b).imag; \ + (to).imag = (from_a).real * (from_b).imag + \ + (from_a).imag * (from_b).real; \ } /* Macro function that multiplies a complex number and a scalar. */ #define CMPLX_MULT_SCALAR(to,from, scalar) \ -{ (to).Real = (from).Real * scalar; \ - (to).Imag = (from).Imag * scalar; \ +{ (to).real = (from).real * scalar; \ + (to).imag = (from).imag * scalar; \ } /* Macro function that implements *= for a complex and a scalar number. */ #define CMPLX_MULT_SELF_SCALAR(cnum, scalar) \ -{ (cnum).Real *= scalar; \ - (cnum).Imag *= scalar; \ +{ (cnum).real *= scalar; \ + (cnum).imag *= scalar; \ } /* Macro function that multiply-assigns a complex number by a scalar. */ #define SCLR_MULT_ASSIGN(to,sclr) \ -{ (to).Real *= (sclr); \ - (to).Imag *= (sclr); \ +{ (to).real *= (sclr); \ + (to).imag *= (sclr); \ } /* Macro function that implements to *= from for complex numbers. */ #define CMPLX_MULT_ASSIGN(to,from) \ -{ RealNumber to_real_ = (to).Real; \ - (to).Real = to_real_ * (from).Real - \ - (to).Imag * (from).Imag; \ - (to).Imag = to_real_ * (from).Imag + \ - (to).Imag * (from).Real; \ +{ realNumber to_real_ = (to).real; \ + (to).real = to_real_ * (from).real - \ + (to).imag * (from).imag; \ + (to).imag = to_real_ * (from).imag + \ + (to).imag * (from).real; \ } /* Macro function that multiplies two complex numbers, the first of which is * conjugated. */ #define CMPLX_CONJ_MULT(to,from_a,from_b) \ -{ (to).Real = (from_a).Real * (from_b).Real + \ - (from_a).Imag * (from_b).Imag; \ - (to).Imag = (from_a).Real * (from_b).Imag - \ - (from_a).Imag * (from_b).Real; \ +{ (to).real = (from_a).real * (from_b).real + \ + (from_a).imag * (from_b).imag; \ + (to).imag = (from_a).real * (from_b).imag - \ + (from_a).imag * (from_b).real; \ } /* Macro function that multiplies two complex numbers and then adds them * to another. to = add + mult_a * mult_b */ #define CMPLX_MULT_ADD(to,mult_a,mult_b,add) \ -{ (to).Real = (mult_a).Real * (mult_b).Real - \ - (mult_a).Imag * (mult_b).Imag + (add).Real; \ - (to).Imag = (mult_a).Real * (mult_b).Imag + \ - (mult_a).Imag * (mult_b).Real + (add).Imag; \ +{ (to).real = (mult_a).real * (mult_b).real - \ + (mult_a).imag * (mult_b).imag + (add).real; \ + (to).imag = (mult_a).real * (mult_b).imag + \ + (mult_a).imag * (mult_b).real + (add).imag; \ } /* Macro function that subtracts the product of two complex numbers from * another. to = subt - mult_a * mult_b */ #define CMPLX_MULT_SUBT(to,mult_a,mult_b,subt) \ -{ (to).Real = (subt).Real - (mult_a).Real * (mult_b).Real + \ - (mult_a).Imag * (mult_b).Imag; \ - (to).Imag = (subt).Imag - (mult_a).Real * (mult_b).Imag - \ - (mult_a).Imag * (mult_b).Real; \ +{ (to).real = (subt).real - (mult_a).real * (mult_b).real + \ + (mult_a).imag * (mult_b).imag; \ + (to).imag = (subt).imag - (mult_a).real * (mult_b).imag - \ + (mult_a).imag * (mult_b).real; \ } /* Macro function that multiplies two complex numbers and then adds them * to another. to = add + mult_a* * mult_b where mult_a* represents mult_a * conjugate. */ #define CMPLX_CONJ_MULT_ADD(to,mult_a,mult_b,add) \ -{ (to).Real = (mult_a).Real * (mult_b).Real + \ - (mult_a).Imag * (mult_b).Imag + (add).Real; \ - (to).Imag = (mult_a).Real * (mult_b).Imag - \ - (mult_a).Imag * (mult_b).Real + (add).Imag; \ +{ (to).real = (mult_a).real * (mult_b).real + \ + (mult_a).imag * (mult_b).imag + (add).real; \ + (to).imag = (mult_a).real * (mult_b).imag - \ + (mult_a).imag * (mult_b).real + (add).imag; \ } /* Macro function that multiplies two complex numbers and then adds them * to another. to += mult_a * mult_b */ #define CMPLX_MULT_ADD_ASSIGN(to,from_a,from_b) \ -{ (to).Real += (from_a).Real * (from_b).Real - \ - (from_a).Imag * (from_b).Imag; \ - (to).Imag += (from_a).Real * (from_b).Imag + \ - (from_a).Imag * (from_b).Real; \ +{ (to).real += (from_a).real * (from_b).real - \ + (from_a).imag * (from_b).imag; \ + (to).imag += (from_a).real * (from_b).imag + \ + (from_a).imag * (from_b).real; \ } /* Macro function that multiplies two complex numbers and then subtracts them * from another. */ #define CMPLX_MULT_SUBT_ASSIGN(to,from_a,from_b) \ -{ (to).Real -= (from_a).Real * (from_b).Real - \ - (from_a).Imag * (from_b).Imag; \ - (to).Imag -= (from_a).Real * (from_b).Imag + \ - (from_a).Imag * (from_b).Real; \ +{ (to).real -= (from_a).real * (from_b).real - \ + (from_a).imag * (from_b).imag; \ + (to).imag -= (from_a).real * (from_b).imag + \ + (from_a).imag * (from_b).real; \ } /* Macro function that multiplies two complex numbers and then adds them * to the destination. to += from_a* * from_b where from_a* represents from_a * conjugate. */ #define CMPLX_CONJ_MULT_ADD_ASSIGN(to,from_a,from_b) \ -{ (to).Real += (from_a).Real * (from_b).Real + \ - (from_a).Imag * (from_b).Imag; \ - (to).Imag += (from_a).Real * (from_b).Imag - \ - (from_a).Imag * (from_b).Real; \ +{ (to).real += (from_a).real * (from_b).real + \ + (from_a).imag * (from_b).imag; \ + (to).imag += (from_a).real * (from_b).imag - \ + (from_a).imag * (from_b).real; \ } /* Macro function that multiplies two complex numbers and then subtracts them * from the destination. to -= from_a* * from_b where from_a* represents from_a * conjugate. */ #define CMPLX_CONJ_MULT_SUBT_ASSIGN(to,from_a,from_b) \ -{ (to).Real -= (from_a).Real * (from_b).Real + \ - (from_a).Imag * (from_b).Imag; \ - (to).Imag -= (from_a).Real * (from_b).Imag - \ - (from_a).Imag * (from_b).Real; \ +{ (to).real -= (from_a).real * (from_b).real + \ + (from_a).imag * (from_b).imag; \ + (to).imag -= (from_a).real * (from_b).imag - \ + (from_a).imag * (from_b).real; \ } /* @@ -491,55 +493,56 @@ typedef struct /* Complex division: to = num / den */ #define CMPLX_DIV(to,num,den) \ -{ RealNumber r_, s_; \ - if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ - ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ - { r_ = (den).Imag / (den).Real; \ - s_ = (den).Real + r_*(den).Imag; \ - (to).Real = ((num).Real + r_*(num).Imag)/s_; \ - (to).Imag = ((num).Imag - r_*(num).Real)/s_; \ +{ realNumber r_, s_; \ + if (((den).real >= (den).imag AND (den).real > -(den).imag) OR \ + ((den).real < (den).imag AND (den).real <= -(den).imag)) \ + { r_ = (den).imag / (den).real; \ + s_ = (den).real + r_*(den).imag; \ + (to).real = ((num).real + r_*(num).imag)/s_; \ + (to).imag = ((num).imag - r_*(num).real)/s_; \ } \ else \ - { r_ = (den).Real / (den).Imag; \ - s_ = (den).Imag + r_*(den).Real; \ - (to).Real = (r_*(num).Real + (num).Imag)/s_; \ - (to).Imag = (r_*(num).Imag - (num).Real)/s_; \ + { r_ = (den).real / (den).imag; \ + s_ = (den).imag + r_*(den).real; \ + (to).real = (r_*(num).real + (num).imag)/s_; \ + (to).imag = (r_*(num).imag - (num).real)/s_; \ } \ } /* Complex division and assignment: num /= den */ #define CMPLX_DIV_ASSIGN(num,den) \ -{ RealNumber r_, s_, t_; \ - if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ - ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ - { r_ = (den).Imag / (den).Real; \ - s_ = (den).Real + r_*(den).Imag; \ - t_ = ((num).Real + r_*(num).Imag)/s_; \ - (num).Imag = ((num).Imag - r_*(num).Real)/s_; \ - (num).Real = t_; \ +{ realNumber r_, s_, t_; \ + if (((den).real >= (den).imag AND (den).real > -(den).imag) OR \ + ((den).real < (den).imag AND (den).real <= -(den).imag)) \ + { r_ = (den).imag / (den).real; \ + s_ = (den).real + r_*(den).imag; \ + t_ = ((num).real + r_*(num).imag)/s_; \ + (num).imag = ((num).imag - r_*(num).real)/s_; \ + (num).real = t_; \ } \ else \ - { r_ = (den).Real / (den).Imag; \ - s_ = (den).Imag + r_*(den).Real; \ - t_ = (r_*(num).Real + (num).Imag)/s_; \ - (num).Imag = (r_*(num).Imag - (num).Real)/s_; \ - (num).Real = t_; \ + { r_ = (den).real / (den).imag; \ + s_ = (den).imag + r_*(den).real; \ + t_ = (r_*(num).real + (num).imag)/s_; \ + (num).imag = (r_*(num).imag - (num).real)/s_; \ + (num).real = t_; \ } \ } /* Complex reciprocation: to = 1.0 / den */ #define CMPLX_RECIPROCAL(to,den) \ -{ RealNumber r_; \ - if (((den).Real >= (den).Imag && (den).Real > -(den).Imag) || \ - ((den).Real < (den).Imag && (den).Real <= -(den).Imag)) \ - { r_ = (den).Imag / (den).Real; \ - (to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \ +{ realNumber r_; \ + if (((den).real >= (den).imag && (den).real > -(den).imag) || \ + ((den).real < (den).imag && (den).real <= -(den).imag)) \ + { r_ = (den).imag / (den).real; \ + (to).imag = -r_*((to).real = 1.0/((den).real + r_*(den).imag)); \ } \ else \ - { r_ = (den).Real / (den).Imag; \ - (to).Real = -r_*((to).Imag = -1.0/((den).Imag + r_*(den).Real));\ + { r_ = (den).real / (den).imag; \ + (to).real = -r_*((to).imag = -1.0/((den).imag + r_*(den).real));\ } \ } +#endif /* _SPDEF_H */ #endif /*_COMPLEX_H */ diff --git a/src/include/memory.h b/src/include/memory.h index 5bfe00ee5..4bf871d87 100644 --- a/src/include/memory.h +++ b/src/include/memory.h @@ -25,4 +25,25 @@ extern void txfree(void *ptr); #define REALLOC(x,y) trealloc((char *)(x),(unsigned)(y)) #define ZERO(PTR,TYPE) (bzero((PTR),sizeof(TYPE))) +#ifdef CIDER + + +#define RALLOC(ptr,type,number)\ +if ((number) && (!(ptr = (type *)calloc((number), (unsigned)(sizeof(type)))))) {\ + return(E_NOMEM);\ +} + +#define XALLOC(ptr,type,number) \ +if ((number) && (!(ptr = (type *)calloc((number), (unsigned)(sizeof(type)))))) {\ + SPfrontEnd->IFerror( E_PANIC, "Out of Memory", NIL(IFuid) );\ + exit( 1 );\ +} + +#define XCALLOC(ptr,type,number) \ +if ((number) && (!(ptr = (type *)calloc((number), (unsigned)(sizeof(type)))))) {\ + fprintf( stderr, "Out of Memory\n" );\ + exit( 1 );\ +} +#endif /* CIDER */ + #endif diff --git a/src/include/smpdefs.h b/src/include/smpdefs.h index a1f7432a2..70cdafac2 100644 --- a/src/include/smpdefs.h +++ b/src/include/smpdefs.h @@ -12,7 +12,7 @@ Modified: 2000 AlansFixes #include #include -#include +#include "complex.h" int SMPaddElt( SMPmatrix *, int , int , double ); double * SMPmakeElt( SMPmatrix * , int , int ); diff --git a/src/maths/sparse/spbuild.c b/src/maths/sparse/spbuild.c index 07d3e9c69..fceff9447 100644 --- a/src/maths/sparse/spbuild.c +++ b/src/maths/sparse/spbuild.c @@ -11,6 +11,7 @@ * >>> User accessible functions contained in this file: * spClear * spGetElement + * spFindElement * spGetAdmittance * spGetQuad * spGetOnes @@ -144,10 +145,89 @@ spClear(void *eMatrix ) +/* + * SINGLE ELEMENT LOCATION IN MATRIX BY INDEX + * + * Finds element [Row,Col] and returns a pointer to it. If element is + * not found then it is created and spliced into matrix. This routine + * is only to be used after spCreate() and before spMNA_Preorder(), + * spFactor() or spOrderAndFactor(). Returns a pointer to the + * Real portion of a MatrixElement. This pointer is later used by + * spADD_xxx_ELEMENT to directly access element. + * + * >>> Returns: + * Returns a pointer to the element. This pointer is then used to directly + * access the element during successive builds. + * + * >>> Arguments: + * Matrix (char *) + * Pointer to the matrix that the element is to be added to. + * Row (int) + * Row index for element. Must be in the range of [0..Size] unless + * the options EXPANDABLE or TRANSLATE are used. Elements placed in + * row zero are discarded. In no case may Row be less than zero. + * Col (int) + * Column index for element. Must be in the range of [0..Size] unless + * the options EXPANDABLE or TRANSLATE are used. Elements placed in + * column zero are discarded. In no case may Col be less than zero. + * + * >>> Local variables: + * pElement (RealNumber *) + * Pointer to the element. + * + * >>> Possible errors: + * spNO_MEMORY + * Error is not cleared in this routine. + */ +RealNumber * +spFindElement( void *eMatrix, int Row, int Col ) +{ +MatrixPtr Matrix = (MatrixPtr)eMatrix; +RealNumber *pElement; +int index; +ElementPtr spcFindElementInCol(); +void Translate(); +/* Begin `spFindElement'. */ + assert( IS_SPARSE( Matrix ) AND Row >= 0 AND Col >= 0 ); + if ((Row == 0) OR (Col == 0)) + return &Matrix->TrashCan.Real; +#if TRANSLATE + Translate( Matrix, &Row, &Col ); + if (Matrix->Error == spNO_MEMORY) return NULL; +#endif + +#if NOT TRANSLATE + assert(Row <= Matrix->Size AND Col <= Matrix->Size); +#endif + +/* + * The condition part of the following if statement tests to see if the + * element resides along the diagonal, if it does then it tests to see + * if the element has been created yet (Diag pointer not NULL). The + * pointer to the element is then assigned to Element after it is cast + * into a pointer to a RealNumber. This casting makes the pointer into + * a pointer to Real. This statement depends on the fact that Real + * is the first record in the MatrixElement structure. + */ + + if ((Row != Col) OR ((pElement = (RealNumber *)Matrix->Diag[Row]) == NULL)) + { +/* + * Element does not exist or does not reside along diagonal. Search + * column for element. As in the if statement above, the pointer to the + * element which is returned by spcFindElementInCol is cast into a + * pointer to Real, a RealNumber. + */ + pElement = (RealNumber*)spcFindElementInCol( Matrix, + &(Matrix->FirstInCol[Col]), + Row, Col, NO ); + } + return pElement; +} diff --git a/src/maths/sparse/spdefs.h b/src/maths/sparse/spdefs.h index f13f42511..8db08028f 100644 --- a/src/maths/sparse/spdefs.h +++ b/src/maths/sparse/spdefs.h @@ -39,14 +39,11 @@ */ #include -#include "complex.h" #undef ABORT #undef MALLOC #undef FREE #undef REALLOC -#undef CMPLX_MULT -#undef CMPLX_RECIPROCAL @@ -91,8 +88,29 @@ /* Macro procedure that swaps two entities. */ #define SWAP(type, a, b) {type swapx; swapx = a; a = b; b = swapx;} -/* Macro function that returns the approx absolute value of a complex number. */ + + +/* Real and Complex numbers definition */ + +#define spREAL double + +/* Begin `realNumber'. */ +typedef spREAL RealNumber, *RealVector; + +/* Begin `ComplexNumber'. */ +typedef struct +{ RealNumber Real; + RealNumber Imag; +} ComplexNumber, *ComplexVector; + +/* Macro function that returns the approx absolute value of a complex + number. */ #define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag)) + +#define CMPLX_ASSIGN_VALUE(cnum, vReal, vImag) \ +{ (cnum).Real = vReal; \ + (cnum).Imag = vImag; \ +} /* Complex assignment statements. */ #define CMPLX_ASSIGN(to,from) \ @@ -111,12 +129,21 @@ { (to).Real = -(from).Real; \ (to).Imag = (from).Imag; \ } + #define CMPLX_CONJ(a) (a).Imag = -(a).Imag + +#define CONJUGATE(a) (a).Imag = -(a).Imag + #define CMPLX_NEGATE(a) \ { (a).Real = -(a).Real; \ (a).Imag = -(a).Imag; \ } +#define CMPLX_NEGATE_SELF(cnum) \ +{ (cnum).Real = -(cnum).Real; \ + (cnum).Imag = -(cnum).Imag; \ +} + /* Macro that returns the approx magnitude (L-1 norm) of a complex number. */ #define CMPLX_1_NORM(a) (ABS((a).Real) + ABS((a).Imag)) @@ -132,6 +159,11 @@ (to).Imag = (from_a).Imag + (from_b).Imag; \ } +/* Macro function that performs addition of a complex and a scalar. */ +#define CMPLX_ADD_SELF_SCALAR(cnum, scalar) \ +{ (cnum).Real += scalar; \ +} + /* Macro function that performs complex subtraction. */ #define CMPLX_SUBT(to,from_a,from_b) \ { (to).Real = (from_a).Real - (from_b).Real; \ @@ -149,7 +181,7 @@ { (to).Real -= (from).Real; \ (to).Imag -= (from).Imag; \ } - + /* Macro function that multiplies a complex number by a scalar. */ #define SCLR_MULT(to,sclr,cmplx) \ { (to).Real = (sclr) * (cmplx).Real; \ @@ -169,13 +201,32 @@ (to).Imag = (from_a).Real * (from_b).Imag + \ (from_a).Imag * (from_b).Real; \ } + +/* Macro function that multiplies a complex number and a scalar. */ +#define CMPLX_MULT_SCALAR(to,from, scalar) \ +{ (to).Real = (from).Real * scalar; \ + (to).Imag = (from).Imag * scalar; \ +} + +/* Macro function that implements *= for a complex and a scalar number. */ + +#define CMPLX_MULT_SELF_SCALAR(cnum, scalar) \ +{ (cnum).Real *= scalar; \ + (cnum).Imag *= scalar; \ +} + +/* Macro function that multiply-assigns a complex number by a scalar. */ +#define SCLR_MULT_ASSIGN(to,sclr) \ +{ (to).Real *= (sclr); \ + (to).Imag *= (sclr); \ +} /* Macro function that implements to *= from for complex numbers. */ #define CMPLX_MULT_ASSIGN(to,from) \ -{ RealNumber to_real_ = (to).Real; \ - (to).Real = to_real_ * (from).Real - \ +{ RealNumber to_Real_ = (to).Real; \ + (to).Real = to_Real_ * (from).Real - \ (to).Imag * (from).Imag; \ - (to).Imag = to_real_ * (from).Imag + \ + (to).Imag = to_Real_ * (from).Imag + \ (to).Imag * (from).Real; \ } @@ -254,11 +305,53 @@ (from_a).Imag * (from_b).Real; \ } +/* + * Macro functions that provide complex division. + */ + +/* Complex division: to = num / den */ +#define CMPLX_DIV(to,num,den) \ +{ RealNumber r_, s_; \ + if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ + ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ + { r_ = (den).Imag / (den).Real; \ + s_ = (den).Real + r_*(den).Imag; \ + (to).Real = ((num).Real + r_*(num).Imag)/s_; \ + (to).Imag = ((num).Imag - r_*(num).Real)/s_; \ + } \ + else \ + { r_ = (den).Real / (den).Imag; \ + s_ = (den).Imag + r_*(den).Real; \ + (to).Real = (r_*(num).Real + (num).Imag)/s_; \ + (to).Imag = (r_*(num).Imag - (num).Real)/s_; \ + } \ +} + +/* Complex division and assignment: num /= den */ +#define CMPLX_DIV_ASSIGN(num,den) \ +{ RealNumber r_, s_, t_; \ + if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ + ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ + { r_ = (den).Imag / (den).Real; \ + s_ = (den).Real + r_*(den).Imag; \ + t_ = ((num).Real + r_*(num).Imag)/s_; \ + (num).Imag = ((num).Imag - r_*(num).Real)/s_; \ + (num).Real = t_; \ + } \ + else \ + { r_ = (den).Real / (den).Imag; \ + s_ = (den).Imag + r_*(den).Real; \ + t_ = (r_*(num).Real + (num).Imag)/s_; \ + (num).Imag = (r_*(num).Imag - (num).Real)/s_; \ + (num).Real = t_; \ + } \ +} + /* Complex reciprocation: to = 1.0 / den */ #define CMPLX_RECIPROCAL(to,den) \ { RealNumber r_; \ - if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ - ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ + if (((den).Real >= (den).Imag && (den).Real > -(den).Imag) || \ + ((den).Real < (den).Imag && (den).Real <= -(den).Imag)) \ { r_ = (den).Imag / (den).Real; \ (to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \ } \ @@ -269,6 +362,12 @@ } + + + + +/* Allocation */ + #define ALLOC(type,number) ((type *)tmalloc((unsigned)(sizeof(type)*(number)))) #define REALLOC(ptr,type,number) \ ptr = (type *)trealloc((char *)ptr,(unsigned)(sizeof(type)*(number)))