/**********
Copyright 1990 Regents of the University of California.  All rights reserved.
Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group
**********/

/*
 * 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 "cmath2.h"


extern void checkseed(void); /* seed random or set by 'set rndseed=value'*/
extern double drand(void);  /* from randnumb.c */
extern double gauss0(void);  /* from randnumb.c */
extern int poisson(double);  /* from randnumb.c */
extern double exprand(double);  /* from randnumb.c */


static double
cx_max_local(void *data, short int type, int length)
{
    double largest = 0.0;

    if (type == VF_COMPLEX) {
	ngcomplex_t *cc = (ngcomplex_t *) data;
	int i;

        for (i = 0; i < length; i++)
            if (cmag(cc[i]) > largest)
                largest = cmag(cc[i]);
    } else {
	double *dd = (double *) data;
	int i;

        for (i = 0; i < length; i++)
            if (FTEcabs(dd[i]) > largest)
                largest = FTEcabs(dd[i]);
    }
    return largest;
}

/* Normalize the data so that the magnitude of the greatest value is 1. */

void *
cx_norm(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double largest = 0.0;

    largest = cx_max_local(data, type, length);
    if (largest == 0.0) {
	fprintf(cp_err, "Error: can't normalize a 0 vector\n");
	return (NULL);
    }

    *newlength = length;
    if (type == VF_COMPLEX) {
	ngcomplex_t *c;
	ngcomplex_t *cc = (ngcomplex_t *) data;
	int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;

        for (i = 0; i < length; i++) {
            realpart(c[i]) = realpart(cc[i]) / largest;
            imagpart(c[i]) = imagpart(cc[i]) / largest;
        }
        return ((void *) c);
    } else {
	double *d;
	double *dd = (double *) data;
	int i;

        d = alloc_d(length);
        *newtype = VF_REAL;

        for (i = 0; i < length; i++)
            d[i] = dd[i] / largest;
        return ((void *) d);
    }
}

void *
cx_uminus(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
	ngcomplex_t *c;
	ngcomplex_t *cc = (ngcomplex_t *) data;
	int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = - realpart(cc[i]);
            imagpart(c[i]) = - imagpart(cc[i]);
        }
        return ((void *) c);
    } else {
	double *d;
	double *dd = (double *) data;
	int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = - dd[i];
        return ((void *) d);
    }
}


/* random integers drawn from a uniform distribution 
   *data in: integer numbers, their absolut values are used, 
       maximum is RAND_MAX (32767) 
   *data out: random integers in interval [0, data[i][ 
       standard library function rand() is used
*/
void *
cx_rnd(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    checkseed();
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            int j, k;
            j = (int)floor(realpart(cc[i]));
            k = (int)floor(imagpart(cc[i]));
            realpart(c[i]) = j ? rand() % j : 0;
            imagpart(c[i]) = k ? rand() % k : 0;
        }
        return ((void *) c);
    } else {
        double *d;
        double *dd = (double *) data;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++) {
            int j;
            j = (int)floor(dd[i]);
            d[i] = j ? rand() % j : 0;
        }
        return ((void *) d);
    }
}

/* random numbers drawn from a uniform distribution 
   *data out: random numbers in interval [-1, 1[ 
*/
void *
cx_sunif(void *data, short int type, int length, int *newlength, short int *newtype)
{
    NG_IGNORE(data);

    *newlength = length;
    checkseed();
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = drand();
            imagpart(c[i]) = drand();
        }
        return ((void *) c);
    } else {
        double *d;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++) {
            d[i] = drand();
        }
        return ((void *) d);
    }
}

/* random numbers drawn from a poisson distribution 
   *data in:  lambda
   *data out: random integers according to poisson distribution, 
       with lambda given by each vector element 
*/
void *
cx_poisson(void *data, short int type, int length, int *newlength, short int *newtype)
{
    NG_IGNORE(data);

    *newlength = length;
    checkseed();
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = poisson (realpart(cc[i]));
            imagpart(c[i]) = poisson (imagpart(cc[i]));
        }
        return ((void *) c);
    } else {
        double *d;
        double *dd = (double *) data;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++) {
            d[i] = poisson(dd[i]);
        }
        return ((void *) d);
    }
}

/* random numbers drawn from an exponential distribution 
   *data in:  Mean values
   *data out: exponentially distributed random numbers, 
       with mean given by each vector element 
*/
void *
cx_exponential(void *data, short int type, int length, int *newlength, short int *newtype)
{
    NG_IGNORE(data);

    *newlength = length;
    checkseed();
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = exprand(realpart(cc[i]));
            imagpart(c[i]) = exprand(imagpart(cc[i]));
        }
        return ((void *) c);
    } else {
        double *d;
        double *dd = (double *) data;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++) {
            d[i] = exprand(dd[i]);
        }
        return ((void *) d);
    }
}

/* random numbers drawn from a Gaussian distribution 
   mean 0, std dev 1
*/
void *
cx_sgauss(void *data, short int type, int length, int *newlength, short int *newtype)
{
    NG_IGNORE(data);

    *newlength = length;
    checkseed();
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = gauss0();
            imagpart(c[i]) = gauss0();
        }
        return ((void *) c);
    } else {
        double *d;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++) {
            d[i] = gauss0();
        }
        return ((void *) d);
    }
}




/* Compute the avg of a vector.
   Created by A.M.Roldan 2005-05-21  */
   
void 
*cx_avg(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double sum_real = 0.0, sum_imag = 0.0;
    int i;

    if (type == VF_REAL) {

        double *d = alloc_d(length);
        double *dd = (double *) data;

        *newtype = VF_REAL;
        *newlength = length;

        for (i = 0; i < length; i++) {
            sum_real += dd[i];
            d[i] = sum_real / (double)(i+1);
        }

        return ((void *) d);

    } else {

        ngcomplex_t *c = alloc_c(length);
        ngcomplex_t *cc = (ngcomplex_t *) data;

        *newtype = VF_COMPLEX;
        *newlength = length;

        for (i = 0; i < length; i++) {
            sum_real += realpart(cc[i]);
            realpart(c[i]) = sum_real / (double)(i+1);

            sum_imag += imagpart(cc[i]);
            imagpart(c[i]) = sum_imag / (double)(i+1);
        }

       return ((void *) c);

    }
}


/* Compute the mean of a vector. */

void *
cx_mean(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = 1;
    rcheck(length > 0, "mean");
    if (type == VF_REAL) {
	double *d;
	double *dd = (double *) data;
	int i;

        d = alloc_d(1);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            *d += dd[i];
        *d /= length;
        return ((void *) d);
    } else {
	ngcomplex_t *c;
	ngcomplex_t *cc = (ngcomplex_t *) data;
	int i;

        c = alloc_c(1);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(*c) += realpart(cc[i]);
            imagpart(*c) += imagpart(cc[i]);
        }
        realpart(*c) /= length;
        imagpart(*c) /= length;
        return ((void *) c);
    }
}


void *
cx_length(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d;

    NG_IGNORE(data);
    NG_IGNORE(type);

    *newlength = 1;
    *newtype = VF_REAL;
    d = alloc_d(1);
    *d = length;
    return ((void *) d);
}

/* Return a vector from 0 to the magnitude of the argument. Length of the
 * argument is irrelevent.
 */


void *
cx_vector(void *data, short int type, int length, int *newlength, short int *newtype)
{
    ngcomplex_t *cc = (ngcomplex_t *) data;
    double *dd = (double *) data;
    int i, len;
    double *d;

    NG_IGNORE(length);

    if (type == VF_REAL)
        len = (int)FTEcabs(*dd);
    else
        len = (int)cmag(*cc);
    if (len == 0)
        len = 1;
    d = alloc_d(len);
    *newlength = len;
    *newtype = VF_REAL;
    for (i = 0; i < len; i++)
        d[i] = i;
    return ((void *) d);
}

/* Create a vector of the given length composed of all ones. */


void *
cx_unitvec(void *data, short int type, int length, int *newlength, short int *newtype)
{
    ngcomplex_t *cc = (ngcomplex_t *) data;
    double *dd = (double *) data;
    int i, len;
    double *d;

    NG_IGNORE(length);

    if (type == VF_REAL)
        len = (int)FTEcabs(*dd);
    else
        len = (int)cmag(*cc);
    if (len == 0)
        len = 1;
    d = alloc_d(len);
    *newlength = len;
    *newtype = VF_REAL;
    for (i = 0; i < len; i++)
        d[i] = 1;
    return ((void *) d);
}

/* Calling methods for these functions are:
 *  cx_something(data1, data2, datatype1, datatype2, length)
 *
 * The length of the two data vectors is always the same, and is the length
 * of the result. The result type is complex iff one of the args is
 * complex.
 */

void *
cx_plus(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 *c, c1, c2;
    int i;

    if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
        d = alloc_d(length);
        for (i = 0; i < length; i++)
            d[i] = dd1[i] + dd2[i];
        return ((void *) d);
    } else {
        c = alloc_c(length);
        for (i = 0; i < length; i++) {
            if (datatype1 == VF_REAL) {
                realpart(c1) = dd1[i];
                imagpart(c1) = 0.0;
            } else {
                realpart(c1) = realpart(cc1[i]);
                imagpart(c1) = imagpart(cc1[i]);
            }
            if (datatype2 == VF_REAL) {
                realpart(c2) = dd2[i];
                imagpart(c2) = 0.0;
            } else {
                realpart(c2) = realpart(cc2[i]);
                imagpart(c2) = imagpart(cc2[i]);
            }
            realpart(c[i]) = realpart(c1) + realpart(c2);
            imagpart(c[i]) = imagpart(c1) + imagpart(c2);
        }
        return ((void *) c);
    }
}

void *
cx_minus(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 *c, c1, c2;
    int i;

    if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
        d = alloc_d(length);
        for (i = 0; i < length; i++)
            d[i] = dd1[i] - dd2[i];
        return ((void *) d);
    } else {
        c = alloc_c(length);
        for (i = 0; i < length; i++) {
            if (datatype1 == VF_REAL) {
                realpart(c1) = dd1[i];
                imagpart(c1) = 0.0;
            } else {
                realpart(c1) = realpart(cc1[i]);
                imagpart(c1) = imagpart(cc1[i]);
            }
            if (datatype2 == VF_REAL) {
                realpart(c2) = dd2[i];
                imagpart(c2) = 0.0;
            } else {
                realpart(c2) = realpart(cc2[i]);
                imagpart(c2) = imagpart(cc2[i]);
            }
            realpart(c[i]) = realpart(c1) - realpart(c2);
            imagpart(c[i]) = imagpart(c1) - imagpart(c2);
        }
        return ((void *) c);
    }
}

void *
cx_times(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 *c, c1, c2;
    int i;

    if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
        d = alloc_d(length);
        for (i = 0; i < length; i++)
            d[i] = dd1[i] * dd2[i];
        return ((void *) d);
    } else {
        c = alloc_c(length);
        for (i = 0; i < length; i++) {
            if (datatype1 == VF_REAL) {
                realpart(c1) = dd1[i];
                imagpart(c1) = 0.0;
            } else {
                realpart(c1) = realpart(cc1[i]);
                imagpart(c1) = imagpart(cc1[i]);
            }
            if (datatype2 == VF_REAL) {
                realpart(c2) = dd2[i];
                imagpart(c2) = 0.0;
            } else {
                realpart(c2) = realpart(cc2[i]);
                imagpart(c2) = imagpart(cc2[i]);
            }
            realpart(c[i]) = realpart(c1) * realpart(c2)
                - imagpart(c1) * imagpart(c2);
            imagpart(c[i]) = imagpart(c1) * realpart(c2)
                + realpart(c1) * imagpart(c2);
        }
        return ((void *) c);
    }
}

void *
cx_mod(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 *c, c1, c2;
    int i, r1, r2, i1, i2, r3, i3;

    if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
        d = alloc_d(length);
        for (i = 0; i < length; i++) {
            r1 = (int)floor(FTEcabs(dd1[i]));
            rcheck(r1 > 0, "mod");
            r2 = (int)floor(FTEcabs(dd2[i]));
            rcheck(r2 > 0, "mod");
            r3 = r1 % r2;
            d[i] = (double) r3;
        }
        return ((void *) d);
    } else {
        c = alloc_c(length);
        for (i = 0; i < length; i++) {
            if (datatype1 == VF_REAL) {
                realpart(c1) = dd1[i];
                imagpart(c1) = 0.0;
            } else {
                realpart(c1) = realpart(cc1[i]);
                imagpart(c1) = imagpart(cc1[i]);
            }
            if (datatype2 == VF_REAL) {
                realpart(c2) = dd2[i];
                imagpart(c2) = 0.0;
            } else {
                realpart(c2) = realpart(cc2[i]);
                imagpart(c2) = imagpart(cc2[i]);
            }
            r1 = (int)floor(FTEcabs(realpart(c1)));
            rcheck(r1 > 0, "mod");
            r2 = (int)floor(FTEcabs(realpart(c2)));
            rcheck(r2 > 0, "mod");
            i1 = (int)floor(FTEcabs(imagpart(c1)));
            rcheck(i1 > 0, "mod");
            i2 = (int)floor(FTEcabs(imagpart(c2)));
            rcheck(i2 > 0, "mod");
            r3 = r1 % r2;
            i3 = i1 % i2;
            realpart(c[i]) = (double) r3;
            imagpart(c[i]) = (double) i3;
        }
        return ((void *) c);
    }
}


/* Routoure JM : Compute the max of a vector. */

void *
cx_max(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = 1;
    /* test if length >0 et affiche un message d'erreur */
    rcheck(length > 0, "mean");
    if (type == VF_REAL) {
      double largest=0.0;
      double *d;
      double *dd = (double *) data;
      int i;
      
      d = alloc_d(1);
      *newtype = VF_REAL;
      largest=dd[0];
      for (i = 1; i < length; i++)
        if (dd[i]>largest) largest=dd[i];
      *d=largest;
      return ((void *) d);
    } else { 
      double largest_real=0.0;
      double largest_complex=0.0;
      ngcomplex_t *c;
      ngcomplex_t *cc = (ngcomplex_t *) data;
      int i;
      
      c = alloc_c(1);
      *newtype = VF_COMPLEX;
      largest_real=realpart(*cc);
      largest_complex=imagpart(*cc);
      for (i = 0; i < length; i++) {
        if (realpart(cc[i])>largest_real) largest_real=realpart(cc[i]);
        if (imagpart(cc[i])>largest_complex) largest_complex=imagpart(cc[i]);
        }
        realpart(*c) = largest_real;
        imagpart(*c) = largest_complex;
        return ((void *) c);
    }
}
/* Routoure JM : Compute the min of a vector. */

void *
cx_min(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = 1;
    /* test if length >0 et affiche un message d'erreur */
    rcheck(length > 0, "mean");
    if (type == VF_REAL) {
      double smallest;
      double *d;
      double *dd = (double *) data;
      int i;
      
      d = alloc_d(1);
      *newtype = VF_REAL;
      smallest=dd[0];
      for (i = 1; i < length; i++)
        if (dd[i]<smallest) smallest=dd[i];
      *d=smallest;
      return ((void *) d);
    } else { 
      double smallest_real;
      double smallest_complex;
      ngcomplex_t *c;
      ngcomplex_t *cc = (ngcomplex_t *) data;
      int i;
      
      c = alloc_c(1);
      *newtype = VF_COMPLEX;
      smallest_real=realpart(*cc);
      smallest_complex=imagpart(*cc);
      for (i = 1; i < length; i++) {
        if (realpart(cc[i])<smallest_real) smallest_real=realpart(cc[i]);
        if (imagpart(cc[i])<smallest_complex) smallest_complex=imagpart(cc[i]);
        }
        realpart(*c) = smallest_real;
        imagpart(*c) = smallest_complex;
        return ((void *) c);
    }
}


/* Routoure JM : Compute the differential  of a vector. */

void *
cx_d(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    /* test if length >0 et affiche un message d'erreur */
    rcheck(length > 0, "deriv");
    if (type == VF_REAL) {
      double *d;
      double *dd = (double *) data;
      int i;
      
      d = alloc_d(length);
      *newtype = VF_REAL;
      d[0]=dd[1]-dd[0];
      d[length-1]=dd[length-1]-dd[length-2];
      for (i = 1; i < length-1; i++)
        d[i]=dd[i+1]-dd[i-1];

      return ((void *) d);
    } else { 

      ngcomplex_t *c;
      ngcomplex_t *cc = (ngcomplex_t *) data;
      int i;
      
      c = alloc_c(length);
      *newtype = VF_COMPLEX;
      realpart(*c)=realpart(cc[1])-realpart(cc[0]);
      imagpart(*c)=imagpart(cc[1])-imagpart(cc[0]);
      realpart(c[length-1])=realpart(cc[length-1])-realpart(cc[length-2]);
      imagpart(c[length-1])=imagpart(cc[length-1])-imagpart(cc[length-2]);
      
      
      for (i = 1; i < (length-1); i++) {
        realpart(c[i])=realpart(cc[i+1])-realpart(cc[i-1]);
        imagpart(c[i])=imagpart(cc[i+1])-imagpart(cc[i-1]);

        }
        return ((void *) c);
    }
}

void *
cx_floor(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
	ngcomplex_t *c;
	ngcomplex_t *cc = (ngcomplex_t *) data;
	int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = floor(realpart(cc[i]));
            imagpart(c[i]) = floor(imagpart(cc[i]));
        }
        return ((void *) c);
    } else {
	double *d;
	double *dd = (double *) data;
	int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = floor(dd[i]);
        return ((void *) d);
    }
}

void *
cx_ceil(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
	ngcomplex_t *c;
	ngcomplex_t *cc = (ngcomplex_t *) data;
	int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = ceil(realpart(cc[i]));
            imagpart(c[i]) = ceil(imagpart(cc[i]));
        }
        return ((void *) c);
    } else {
	double *d;
	double *dd = (double *) data;
	int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = ceil(dd[i]);
        return ((void *) d);
    }
}

void *
cx_nint(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
	ngcomplex_t *c;
	ngcomplex_t *cc = (ngcomplex_t *) data;
	int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = nearbyint(realpart(cc[i]));
            imagpart(c[i]) = nearbyint(imagpart(cc[i]));
        }
        return ((void *) c);
    } else {
	double *d;
	double *dd = (double *) data;
	int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = nearbyint(dd[i]);
        return ((void *) d);
    }
}