diff --git a/src/maths/cmaths/cmath2.c b/src/maths/cmaths/cmath2.c index 7fbcdd62b..c27b3a5b6 100644 --- a/src/maths/cmaths/cmath2.c +++ b/src/maths/cmaths/cmath2.c @@ -30,15 +30,15 @@ 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; + ngcomplex_t *cc = (ngcomplex_t *) data; + int i; for (i = 0; i < length; i++) if (largest < cmag(cc[i])) largest = cmag(cc[i]); } else { - double *dd = (double *) data; - int i; + double *dd = (double *) data; + int i; for (i = 0; i < length; i++) if (largest < FTEcabs(dd[i])) @@ -47,6 +47,7 @@ cx_max_local(void *data, short int type, int length) return largest; } + /* Normalize the data so that the magnitude of the greatest value is 1. */ void * @@ -56,15 +57,15 @@ cx_norm(void *data, short int type, int length, int *newlength, short int *newty largest = cx_max_local(data, type, length); if (largest == 0.0) { - fprintf(cp_err, "Error: can't normalize a 0 vector\n"); - return (NULL); + 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; + ngcomplex_t *c; + ngcomplex_t *cc = (ngcomplex_t *) data; + int i; c = alloc_c(length); *newtype = VF_COMPLEX; @@ -75,9 +76,9 @@ cx_norm(void *data, short int type, int length, int *newlength, short int *newty } return ((void *) c); } else { - double *d; - double *dd = (double *) data; - int i; + double *d; + double *dd = (double *) data; + int i; d = alloc_d(length); *newtype = VF_REAL; @@ -88,14 +89,15 @@ cx_norm(void *data, short int type, int length, int *newlength, short int *newty } } + 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; + ngcomplex_t *c; + ngcomplex_t *cc = (ngcomplex_t *) data; + int i; c = alloc_c(length); *newtype = VF_COMPLEX; @@ -105,9 +107,9 @@ cx_uminus(void *data, short int type, int length, int *newlength, short int *new } return ((void *) c); } else { - double *d; - double *dd = (double *) data; - int i; + double *d; + double *dd = (double *) data; + int i; d = alloc_d(length); *newtype = VF_REAL; @@ -124,6 +126,7 @@ cx_uminus(void *data, short int type, int length, int *newlength, short int *new * 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) { @@ -160,9 +163,11 @@ cx_rnd(void *data, short int type, int length, int *newlength, short int *newtyp } } + /* 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) { @@ -194,11 +199,13 @@ cx_sunif(void *data, short int type, int length, int *newlength, short int *newt } } + /* 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) { @@ -230,11 +237,13 @@ cx_poisson(void *data, short int type, int length, int *newlength, short int *ne } } + /* 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) { @@ -266,9 +275,11 @@ cx_exponential(void *data, short int type, int length, int *newlength, short int } } -/* random numbers drawn from a Gaussian distribution + +/* 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) { @@ -301,12 +312,10 @@ cx_sgauss(void *data, short int type, int length, int *newlength, short int *new } - - /* 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) { @@ -344,7 +353,7 @@ cx_avg(void *data, short int type, int length, int *newlength, short int *newtyp imagpart(c[i]) = sum_imag / (double)(i+1); } - return ((void *) c); + return ((void *) c); } } @@ -358,9 +367,9 @@ cx_mean(void *data, short int type, int length, int *newlength, short int *newty *newlength = 1; rcheck(length > 0, "mean"); if (type == VF_REAL) { - double *d; - double *dd = (double *) data; - int i; + double *d; + double *dd = (double *) data; + int i; d = alloc_d(1); *newtype = VF_REAL; @@ -369,9 +378,9 @@ cx_mean(void *data, short int type, int length, int *newlength, short int *newty *d /= length; return ((void *) d); } else { - ngcomplex_t *c; - ngcomplex_t *cc = (ngcomplex_t *) data; - int i; + ngcomplex_t *c; + ngcomplex_t *cc = (ngcomplex_t *) data; + int i; c = alloc_c(1); *newtype = VF_COMPLEX; @@ -401,11 +410,11 @@ cx_length(void *data, short int type, int length, int *newlength, short int *new 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) { @@ -430,8 +439,8 @@ cx_vector(void *data, short int type, int length, int *newlength, short int *new return ((void *) d); } -/* Create a vector of the given length composed of all ones. */ +/* 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) @@ -457,6 +466,7 @@ cx_unitvec(void *data, short int type, int length, int *newlength, short int *ne return ((void *) d); } + /* Calling methods for these functions are: * cx_something(data1, data2, datatype1, datatype2, length) * @@ -503,6 +513,7 @@ cx_plus(void *data1, void *data2, short int datatype1, short int datatype2, int } } + void * cx_minus(void *data1, void *data2, short int datatype1, short int datatype2, int length) { @@ -541,6 +552,7 @@ cx_minus(void *data1, void *data2, short int datatype1, short int datatype2, int } } + void * cx_times(void *data1, void *data2, short int datatype1, short int datatype2, int length) { @@ -581,6 +593,7 @@ cx_times(void *data1, void *data2, short int datatype1, short int datatype2, int } } + void * cx_mod(void *data1, void *data2, short int datatype1, short int datatype2, int length) { @@ -645,41 +658,43 @@ cx_max(void *data, short int type, int length, int *newlength, short int *newtyp /* 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 (largest < dd[i]) - 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 (largest_real < realpart(cc[i])) - largest_real = realpart(cc[i]); - if (largest_complex < imagpart(cc[i])) - largest_complex = imagpart(cc[i]); + 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 (largest < dd[i]) + 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 (largest_real < realpart(cc[i])) + largest_real = realpart(cc[i]); + if (largest_complex < imagpart(cc[i])) + 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 * @@ -689,35 +704,35 @@ cx_min(void *data, short int type, int length, int *newlength, short int *newtyp /* 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 (smallest > dd[i]) - 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 (smallest_real > realpart(cc[i])) - smallest_real = realpart(cc[i]); - if (smallest_complex > imagpart(cc[i])) - smallest_complex = imagpart(cc[i]); + 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 (smallest > dd[i]) + 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 (smallest_real > realpart(cc[i])) + smallest_real = realpart(cc[i]); + if (smallest_complex > imagpart(cc[i])) + smallest_complex = imagpart(cc[i]); } realpart(*c) = smallest_real; imagpart(*c) = smallest_complex; @@ -735,49 +750,50 @@ cx_d(void *data, short int type, int length, int *newlength, short int *newtype) /* 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]); + 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; + ngcomplex_t *c; + ngcomplex_t *cc = (ngcomplex_t *) data; + int i; c = alloc_c(length); *newtype = VF_COMPLEX; @@ -787,9 +803,9 @@ cx_floor(void *data, short int type, int length, int *newlength, short int *newt } return ((void *) c); } else { - double *d; - double *dd = (double *) data; - int i; + double *d; + double *dd = (double *) data; + int i; d = alloc_d(length); *newtype = VF_REAL; @@ -799,14 +815,15 @@ cx_floor(void *data, short int type, int length, int *newlength, short int *newt } } + 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; + ngcomplex_t *c; + ngcomplex_t *cc = (ngcomplex_t *) data; + int i; c = alloc_c(length); *newtype = VF_COMPLEX; @@ -816,9 +833,9 @@ cx_ceil(void *data, short int type, int length, int *newlength, short int *newty } return ((void *) c); } else { - double *d; - double *dd = (double *) data; - int i; + double *d; + double *dd = (double *) data; + int i; d = alloc_d(length); *newtype = VF_REAL; @@ -828,14 +845,15 @@ cx_ceil(void *data, short int type, int length, int *newlength, short int *newty } } + 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; + ngcomplex_t *c; + ngcomplex_t *cc = (ngcomplex_t *) data; + int i; c = alloc_c(length); *newtype = VF_COMPLEX; @@ -845,9 +863,9 @@ cx_nint(void *data, short int type, int length, int *newlength, short int *newty } return ((void *) c); } else { - double *d; - double *dd = (double *) data; - int i; + double *d; + double *dd = (double *) data; + int i; d = alloc_d(length); *newtype = VF_REAL;