Qucs-core  0.0.19
vector.cpp
Go to the documentation of this file.
00001 /*
00002  * vector.cpp - vector class implementation
00003  *
00004  * Copyright (C) 2003, 2004, 2005, 2006, 2007 Stefan Jahn <stefan@lkcc.org>
00005  * Copyright (C) 2006, 2007 Gunther Kraut <gn.kraut@t-online.de>
00006  *
00007  * This is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2, or (at your option)
00010  * any later version.
00011  *
00012  * This software is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this package; see the file COPYING.  If not, write to
00019  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00020  * Boston, MA 02110-1301, USA.
00021  *
00022  * $Id$
00023  *
00024  */
00025 
00026 #if HAVE_CONFIG_H
00027 # include <config.h>
00028 #endif
00029 
00030 #include <limits>
00031 
00032 #include <stdio.h>
00033 #include <stdlib.h>
00034 #include <string.h>
00035 #include <cmath>
00036 #include <float.h>
00037 #include <assert.h>
00038 
00039 #include "real.h"
00040 #include "complex.h"
00041 #include "object.h"
00042 #include "logging.h"
00043 #include "strlist.h"
00044 #include "vector.h"
00045 #include "consts.h"
00046 
00047 namespace qucs {
00048 
00049 // Constructor creates an unnamed instance of the vector class.
00050 vector::vector () : object () {
00051   capacity = size = 0;
00052   data = NULL;
00053   dependencies = NULL;
00054   origin = NULL;
00055   requested = 0;
00056   next = prev = nullptr;
00057 }
00058 
00059 /* Constructor creates an unnamed instance of the vector class with a
00060    given initial size. */
00061 vector::vector (int s) : object () {
00062   assert (s >= 0);
00063   capacity = size = s;
00064   data = s > 0 ? (nr_complex_t *)
00065     calloc (capacity, sizeof (nr_complex_t)) : NULL;
00066   dependencies = NULL;
00067   origin = NULL;
00068   requested = 0;
00069   next = prev = nullptr;
00070 }
00071 
00072 /* Constructor creates an unnamed instance of the vector class with a
00073    given initial size and content. */
00074 vector::vector (int s, nr_complex_t val) : object () {
00075   assert (s >= 0);
00076   capacity = size = s;
00077   data = s > 0 ? (nr_complex_t *)
00078     calloc (capacity, sizeof (nr_complex_t)) : NULL;
00079   for (int i = 0; i < s; i++) data[i] = val;
00080   dependencies = NULL;
00081   origin = NULL;
00082   requested = 0;
00083   next = prev = nullptr;
00084 }
00085 
00086 // Constructor creates an named instance of the vector class.
00087 vector::vector (const std::string &n) : object (n) {
00088   capacity = size = 0;
00089   data = NULL;
00090   dependencies = NULL;
00091   origin = NULL;
00092   requested = 0;
00093   next = prev = nullptr;
00094 }
00095 
00096 /* This constructor creates a named instance of the vector class with
00097    a given initial size. */
00098   vector::vector (const std::string &n, int s) : object (n) {
00099   assert (s >= 0);
00100   capacity = size = s;
00101   data = s > 0 ? (nr_complex_t *)
00102     calloc (capacity, sizeof (nr_complex_t)) : NULL;
00103   dependencies = NULL;
00104   origin = NULL;
00105   requested = 0;
00106   next = prev = nullptr;
00107 }
00108 
00109 /* The copy constructor creates a new instance based on the given
00110    vector object. */
00111 vector::vector (const vector & v) : object (v) {
00112   size = v.size;
00113   capacity = v.capacity;
00114   data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity);
00115   memcpy (data, v.data, sizeof (nr_complex_t) * size);
00116   dependencies = v.dependencies ? new strlist (*v.dependencies) : NULL;
00117   origin = v.origin ? strdup (v.origin) : NULL;
00118   requested = v.requested;
00119   next = v.next;
00120   prev = v.prev;
00121 }
00122 
00123 /* The assignment copy constructor creates a new instance based on the
00124    given vector object.  It copies the data only and leaves any other
00125    properties untouched. */
00126 const vector& vector::operator=(const vector & v) {
00127   if (&v != this) {
00128     size = v.size;
00129     capacity = v.capacity;
00130     if (data) { free (data); data = NULL; }
00131     if (capacity > 0) {
00132       data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity);
00133       if (size > 0) memcpy (data, v.data, sizeof (nr_complex_t) * size);
00134     }
00135   }
00136   return *this;
00137 }
00138 
00139 // Destructor deletes a vector object.
00140 vector::~vector () {
00141   free (data);
00142   delete dependencies;
00143   free (origin);
00144 }
00145 
00146 // Returns data dependencies.
00147 strlist * vector::getDependencies (void) {
00148   return dependencies;
00149 }
00150 
00151 // Sets the data dependencies.
00152 void vector::setDependencies (strlist * s) {
00153   delete dependencies;
00154   dependencies = s;
00155 }
00156 
00157 /* The function appends a new complex data item to the end of the
00158    vector and ensures that the vector can hold the increasing number
00159    of data items. */
00160 void vector::add (nr_complex_t c) {
00161   if (data == NULL) {
00162     size = 0; capacity = 64;
00163     data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity);
00164   }
00165   else if (size >= capacity) {
00166     capacity *= 2;
00167     data = (nr_complex_t *) realloc (data, sizeof (nr_complex_t) * capacity);
00168   }
00169   data[size++] = c;
00170 }
00171 
00172 /* This function appends the given vector to the vector. */
00173 void vector::add (vector * v) {
00174   if (v != NULL) {
00175     if (data == NULL) {
00176       size = 0; capacity = v->getSize ();
00177       data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity);
00178     }
00179     else if (size + v->getSize () > capacity) {
00180       capacity += v->getSize ();
00181       data = (nr_complex_t *) realloc (data, sizeof (nr_complex_t) * capacity);
00182     }
00183     for (int i = 0; i < v->getSize (); i++) data[size++] = v->get (i);
00184   }
00185 }
00186 
00187 // Returns the complex data item at the given position.
00188 nr_complex_t vector::get (int i) {
00189   return data[i];
00190 }
00191 
00192 void vector::set (nr_double_t d, int i) {
00193   data[i] = nr_complex_t (d);
00194 }
00195 
00196 void vector::set (const nr_complex_t z, int i) {
00197   data[i] = nr_complex_t (z);
00198 }
00199 
00200 // The function returns the current size of the vector.
00201 int vector::getSize (void) const {
00202   return size;
00203 }
00204 
00205 int vector::checkSizes (vector v1, vector v2) {
00206   if (v1.getSize () != v2.getSize ()) {
00207     logprint (LOG_ERROR, "vector '%s' and '%s' have different sizes\n",
00208               v1.getName (), v2.getName ());
00209     return 0;
00210   }
00211   return 1;
00212 }
00213 
00214 // searches the maximum value of the vector elements.
00215 // complex numbers in the 1. and 4. quadrant are counted as "abs(c)".
00216 // complex numbers in the 2. and 3. quadrant are counted as "-abs(c)".
00217 nr_double_t vector::maximum (void) {
00218   nr_complex_t c;
00219   nr_double_t d, max_D = -std::numeric_limits<nr_double_t>::max();
00220   for (int i = 0; i < getSize (); i++) {
00221     c = data[i];
00222     d = fabs (arg (c)) < pi_over_2 ? abs (c) : -abs (c);
00223     if (d > max_D) max_D = d;
00224   }
00225   return max_D;
00226 }
00227 
00228 // searches the minimum value of the vector elements.
00229 // complex numbers in the 1. and 4. quadrant are counted as "abs(c)".
00230 // complex numbers in the 2. and 3. quadrant are counted as "-abs(c)".
00231 nr_double_t vector::minimum (void) {
00232   nr_complex_t c;
00233   nr_double_t d, min_D = +std::numeric_limits<nr_double_t>::max();
00234   for (int i = 0; i < getSize (); i++) {
00235     c = data[i];
00236     d = fabs (arg (c)) < pi_over_2 ? abs (c) : -abs (c);
00237     if (d < min_D) min_D = d;
00238   }
00239   return min_D;
00240 }
00241 
00242 /* Unwraps a phase vector in radians.  Adds +/- 2*Pi if consecutive
00243    values jump about |Pi|. */
00244 vector unwrap (vector v, nr_double_t tol, nr_double_t step) {
00245   vector result (v.getSize ());
00246   nr_double_t add = 0;
00247   result (0) = v (0);
00248   for (int i = 1; i < v.getSize (); i++) {
00249     nr_double_t diff = real (v (i) - v (i-1));
00250     if (diff > +tol) {
00251       add -= step;
00252     } else if (diff < -tol) {
00253       add += step;
00254     }
00255     result (i) = v (i) + add;
00256   }
00257   return result;
00258 }
00259 
00260 nr_complex_t sum (vector v) {
00261   nr_complex_t result (0.0);
00262   for (int i = 0; i < v.getSize (); i++) result += v.get (i);
00263   return result;
00264 }
00265 
00266 nr_complex_t prod (vector v) {
00267   nr_complex_t result (1.0);
00268   for (int i = 0; i < v.getSize (); i++) result *= v.get (i);
00269   return result;
00270 }
00271 
00272 nr_complex_t avg (vector v) {
00273   nr_complex_t result (0.0);
00274   for (int i = 0; i < v.getSize (); i++) result += v.get (i);
00275   return result / (nr_double_t) v.getSize ();
00276 }
00277 
00278 vector signum (vector v) {
00279   vector result (v);
00280   for (int i = 0; i < v.getSize (); i++) result.set (signum (v.get (i)), i);
00281   return result;
00282 }
00283 
00284 vector sign (vector v) {
00285   vector result (v);
00286   for (int i = 0; i < v.getSize (); i++) result.set (sign (v.get (i)), i);
00287   return result;
00288 }
00289 
00290 vector xhypot (vector v, const nr_complex_t z) {
00291   vector result (v);
00292   for (int i = 0; i < v.getSize (); i++) result.set (xhypot (v.get(i), z), i);
00293   return result;
00294 }
00295 
00296 vector xhypot (vector v, const nr_double_t d) {
00297   vector result (v);
00298   for (int i = 0; i < v.getSize (); i++) result.set (xhypot (v.get(i), d), i);
00299   return result;
00300 }
00301 
00302 vector xhypot (const nr_complex_t z, vector v) {
00303   vector result (v);
00304   for (int i = 0; i < v.getSize (); i++) result.set (xhypot (z, v.get (i)), i);
00305   return result;
00306 }
00307 
00308 vector xhypot (const nr_double_t d, vector v) {
00309   vector result (v);
00310   for (int i = 0; i < v.getSize (); i++) result.set (xhypot (d, v.get (i)), i);
00311   return result;
00312 }
00313 
00314 vector xhypot (vector v1, vector v2) {
00315   int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize ();
00316   if (len1 >= len2) {
00317     assert (len1 % len2 == 0);
00318     len = len1;
00319   } else {
00320     assert (len2 % len1 == 0);
00321     len = len2;
00322   }
00323   vector res (len);
00324   for (j = i = n = 0; n < len; n++) {
00325     res (n) = xhypot (v1 (i), v2 (j));
00326     if (++i >= len1) i = 0; if (++j >= len2) j = 0;
00327   }
00328   return res;
00329 }
00330 
00331 vector sinc (vector v) {
00332   vector result (v);
00333   for (int i = 0; i < v.getSize (); i++) result.set (sinc (v.get (i)), i);
00334   return result;
00335 }
00336 
00337 vector abs (vector v) {
00338   vector result (v);
00339   for (int i = 0; i < v.getSize (); i++) result.set (abs (v.get (i)), i);
00340   return result;
00341 }
00342 
00343 vector norm (vector v) {
00344   vector result (v);
00345   for (int i = 0; i < v.getSize (); i++) result.set (norm (v.get (i)), i);
00346   return result;
00347 }
00348 
00349 vector arg (vector v) {
00350   vector result (v);
00351   for (int i = 0; i < v.getSize (); i++) result.set (arg (v.get (i)), i);
00352   return result;
00353 }
00354 
00355 vector real (vector v) {
00356   vector result (v);
00357   for (int i = 0; i < v.getSize (); i++) result.set (real (v.get (i)), i);
00358   return result;
00359 }
00360 
00361 vector imag (vector v) {
00362   vector result (v);
00363   for (int i = 0; i < v.getSize (); i++) result.set (imag (v.get (i)), i);
00364   return result;
00365 }
00366 
00367 vector conj (vector v) {
00368   vector result (v);
00369   for (int i = 0; i < v.getSize (); i++) result.set (conj (v.get (i)), i);
00370   return result;
00371 }
00372 
00373 vector dB (vector v) {
00374   vector result (v);
00375   for (int i = 0; i < v.getSize (); i++)
00376     result.set (10.0 * std::log10 (norm (v.get (i))), i);
00377   return result;
00378 }
00379 
00380 vector sqrt (vector v) {
00381   vector result (v);
00382   for (int i = 0; i < v.getSize (); i++) result.set (sqrt (v.get (i)), i);
00383   return result;
00384 }
00385 
00386 vector exp (vector v) {
00387   vector result (v);
00388   for (int i = 0; i < v.getSize (); i++) result.set (exp (v.get (i)), i);
00389   return result;
00390 }
00391 
00392 vector limexp (vector v) {
00393   vector result (v);
00394   for (int i = 0; i < v.getSize (); i++) result.set (limexp (v.get (i)), i);
00395   return result;
00396 }
00397 
00398 vector log (vector v) {
00399   vector result (v);
00400   for (int i = 0; i < v.getSize (); i++) result.set (log (v.get (i)), i);
00401   return result;
00402 }
00403 
00404 vector log10 (vector v) {
00405   vector result (v);
00406   for (int i = 0; i < v.getSize (); i++) result.set (log10 (v.get (i)), i);
00407   return result;
00408 }
00409 
00410 vector log2 (vector v) {
00411   vector result (v);
00412   for (int i = 0; i < v.getSize (); i++) result.set (log2 (v.get (i)), i);
00413   return result;
00414 }
00415 
00416 vector pow (vector v, const nr_complex_t z) {
00417   vector result (v);
00418   for (int i = 0; i < v.getSize (); i++) result.set (pow (v.get(i), z), i);
00419   return result;
00420 }
00421 
00422 vector pow (vector v, const nr_double_t d) {
00423   vector result (v);
00424   for (int i = 0; i < v.getSize (); i++) result.set (pow (v.get(i), d), i);
00425   return result;
00426 }
00427 
00428 vector pow (const nr_complex_t z, vector v) {
00429   vector result (v);
00430   for (int i = 0; i < v.getSize (); i++) result.set (pow (z, v.get (i)), i);
00431   return result;
00432 }
00433 
00434 vector pow (const nr_double_t d, vector v) {
00435   vector result (v);
00436   for (int i = 0; i < v.getSize (); i++) result.set (pow (d, v.get (i)), i);
00437   return result;
00438 }
00439 
00440 vector pow (vector v1, vector v2) {
00441   int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize ();
00442   if (len1 >= len2) {
00443     assert (len1 % len2 == 0);
00444     len = len1;
00445   } else {
00446     assert (len2 % len1 == 0);
00447     len = len2;
00448   }
00449   vector res (len);
00450   for (j = i = n = 0; n < len; n++) {
00451     res (n) = pow (v1 (i), v2 (j));
00452     if (++i >= len1) i = 0; if (++j >= len2) j = 0;
00453   }
00454   return res;
00455 }
00456 
00457 vector sin (vector v) {
00458   vector result (v);
00459   for (int i = 0; i < v.getSize (); i++) result.set (sin (v.get (i)), i);
00460   return result;
00461 }
00462 
00463 vector asin (vector v) {
00464   vector result (v);
00465   for (int i = 0; i < v.getSize (); i++) result.set (asin (v.get (i)), i);
00466   return result;
00467 }
00468 
00469 vector acos (vector v) {
00470   vector result (v);
00471   for (int i = 0; i < v.getSize (); i++) result.set (acos (v.get (i)), i);
00472   return result;
00473 }
00474 
00475 vector cos (vector v) {
00476   vector result (v);
00477   for (int i = 0; i < v.getSize (); i++) result.set (cos (v.get (i)), i);
00478   return result;
00479 }
00480 
00481 vector tan (vector v) {
00482   vector result (v);
00483   for (int i = 0; i < v.getSize (); i++) result.set (tan (v.get (i)), i);
00484   return result;
00485 }
00486 
00487 vector atan (vector v) {
00488   vector result (v);
00489   for (int i = 0; i < v.getSize (); i++) result.set (atan (v.get (i)), i);
00490   return result;
00491 }
00492 
00493 vector cot (vector v) {
00494   vector result (v);
00495   for (int i = 0; i < v.getSize (); i++) result.set (cot (v.get (i)), i);
00496   return result;
00497 }
00498 
00499 vector acot (vector v) {
00500   vector result (v);
00501   for (int i = 0; i < v.getSize (); i++) result.set (acot (v.get (i)), i);
00502   return result;
00503 }
00504 
00505 vector sinh (vector v) {
00506   vector result (v);
00507   for (int i = 0; i < v.getSize (); i++) result.set (sinh (v.get (i)), i);
00508   return result;
00509 }
00510 
00511 vector asinh (vector v) {
00512   vector result (v);
00513   for (int i = 0; i < v.getSize (); i++) result.set (asinh (v.get (i)), i);
00514   return result;
00515 }
00516 
00517 vector cosh (vector v) {
00518   vector result (v);
00519   for (int i = 0; i < v.getSize (); i++) result.set (cosh (v.get (i)), i);
00520   return result;
00521 }
00522 
00523 vector sech (vector v) {
00524   vector result (v);
00525   for (int i = 0; i < v.getSize (); i++) result.set (sech (v.get (i)), i);
00526   return result;
00527 }
00528 
00529 vector cosech (vector v) {
00530   vector result (v);
00531   for (int i = 0; i < v.getSize (); i++) result.set (cosech (v.get (i)), i);
00532   return result;
00533 }
00534 
00535 vector acosh (vector v) {
00536   vector result (v);
00537   for (int i = 0; i < v.getSize (); i++) result.set (acosh (v.get (i)), i);
00538   return result;
00539 }
00540 
00541 vector asech (vector v) {
00542   vector result (v);
00543   for (int i = 0; i < v.getSize (); i++) result.set (asech (v.get (i)), i);
00544   return result;
00545 }
00546 
00547 vector tanh (vector v) {
00548   vector result (v);
00549   for (int i = 0; i < v.getSize (); i++) result.set (tanh (v.get (i)), i);
00550   return result;
00551 }
00552 
00553 vector atanh (vector v) {
00554   vector result (v);
00555   for (int i = 0; i < v.getSize (); i++) result.set (atanh (v.get (i)), i);
00556   return result;
00557 }
00558 
00559 vector coth (vector v) {
00560   vector result (v);
00561   for (int i = 0; i < v.getSize (); i++) result.set (coth (v.get (i)), i);
00562   return result;
00563 }
00564 
00565 vector acoth (vector v) {
00566   vector result (v);
00567   for (int i = 0; i < v.getSize (); i++) result.set (acoth (v.get (i)), i);
00568   return result;
00569 }
00570 
00571 // converts impedance to reflexion coefficient
00572 vector ztor (vector v, nr_complex_t zref) {
00573   vector result (v);
00574   for (int i = 0; i < v.getSize (); i++) result (i) = ztor (v (i), zref);
00575   return result;
00576 }
00577 
00578 // converts admittance to reflexion coefficient
00579 vector ytor (vector v, nr_complex_t zref) {
00580   vector result (v);
00581   for (int i = 0; i < v.getSize (); i++) result (i) = ytor (v (i), zref);
00582   return result;
00583 }
00584 
00585 // converts reflexion coefficient to impedance
00586 vector rtoz (vector v, nr_complex_t zref) {
00587   vector result (v);
00588   for (int i = 0; i < v.getSize (); i++) result (i) = rtoz (v (i), zref);
00589   return result;
00590 }
00591 
00592 // converts reflexion coefficient to admittance
00593 vector rtoy (vector v, nr_complex_t zref) {
00594   vector result (v);
00595   for (int i = 0; i < v.getSize (); i++) result (i) = rtoy (v (i), zref);
00596   return result;
00597 }
00598 
00599 // differentiates 'var' with respect to 'dep' exactly 'n' times
00600 vector diff (vector var, vector dep, int n) {
00601   int k, xi, yi, exchange = 0;
00602   vector x, y;
00603   // exchange dependent and independent variable if necessary
00604   if (var.getSize () < dep.getSize ()) {
00605     x = vector (var);
00606     y = vector (dep);
00607     exchange++;
00608   }
00609   else {
00610     x = vector (dep);
00611     y = vector (var);
00612   }
00613   assert (y.getSize () % x.getSize () == 0 && x.getSize () >= 2);
00614   vector result (y);
00615   nr_complex_t c;
00616 
00617   for (k = 0; k < n; k++) {  // differentiate n times
00618     for (yi = xi = 0; yi < y.getSize (); yi++, xi++) {
00619       if (xi == x.getSize ()) xi = 0; // roll through independent vector
00620       if (xi == 0) {
00621         c = (y.get (yi + 1) - y.get (yi)) / (x.get (xi + 1) - x.get (xi));
00622       } else if (xi == x.getSize () - 1) {
00623         c = (y.get (yi) - y.get (yi - 1)) / (x.get (xi) - x.get (xi - 1));
00624       }
00625       else {
00626         c =
00627           ((y.get (yi) - y.get (yi - 1)) / (x.get (xi) - x.get (xi - 1)) +
00628            (y.get (yi + 1) - y.get (yi)) / (x.get (xi + 1) - x.get (xi))) /
00629           2.0;
00630       }
00631       result.set (exchange ? 1.0 / c : c, yi);
00632     }
00633     y = result;
00634   }
00635   return result;
00636 }
00637 
00638 vector vector::operator=(const nr_complex_t c) {
00639   for (int i = 0; i < size; i++) data[i] = c;
00640   return *this;
00641 }
00642 
00643 vector vector::operator=(const nr_double_t d) {
00644   for (int i = 0; i < size; i++) data[i] = d;
00645   return *this;
00646 }
00647 
00648 vector vector::operator+=(vector v) {
00649   int i, n, len = v.getSize ();
00650   assert (size % len == 0);
00651   for (i = n = 0; i < size; i++) { data[i] += v (n); if (++n >= len) n = 0; }
00652   return *this;
00653 }
00654 
00655 vector vector::operator+=(const nr_complex_t c) {
00656   for (int i = 0; i < size; i++) data[i] += c;
00657   return *this;
00658 }
00659 
00660 vector vector::operator+=(const nr_double_t d) {
00661   for (int i = 0; i < size; i++) data[i] += d;
00662   return *this;
00663 }
00664 
00665 vector operator+(vector v1, vector v2) {
00666   int len1 = v1.getSize (), len2 = v2.getSize ();
00667   vector result;
00668   if (len1 >= len2) {
00669     result  = v1;
00670     result += v2;
00671   } else {
00672     result  = v2;
00673     result += v1;
00674   }
00675   return result;
00676 }
00677 
00678 vector operator+(vector v, const nr_complex_t c) {
00679   vector result (v);
00680   result += c;
00681   return result;
00682 }
00683 
00684 vector operator+(const nr_complex_t c, vector v) {
00685   return v + c;
00686 }
00687 
00688 vector operator+(vector v, const nr_double_t d) {
00689   vector result (v);
00690   result += d;
00691   return result;
00692 }
00693 
00694 vector operator+(const nr_double_t d, vector v) {
00695   return v + d;
00696 }
00697 
00698 vector vector::operator-() {
00699   vector result (size);
00700   for (int i = 0; i < size; i++) result (i) = -data[i];
00701   return result;
00702 }
00703 
00704 vector vector::operator-=(vector v) {
00705   int i, n, len = v.getSize ();
00706   assert (size % len == 0);
00707   for (i = n = 0; i < size; i++) { data[i] -= v (n); if (++n >= len) n = 0; }
00708   return *this;
00709 }
00710 
00711 vector vector::operator-=(const nr_complex_t c) {
00712   for (int i = 0; i < size; i++) data[i] -= c;
00713   return *this;
00714 }
00715 
00716 vector vector::operator-=(const nr_double_t d) {
00717   for (int i = 0; i < size; i++) data[i] -= d;
00718   return *this;
00719 }
00720 
00721 vector operator-(vector v1, vector v2) {
00722   int len1 = v1.getSize (), len2 = v2.getSize ();
00723   vector result;
00724   if (len1 >= len2) {
00725     result  = v1;
00726     result -= v2;
00727   } else {
00728     result  = -v2;
00729     result += v1;
00730   }
00731   return result;
00732 }
00733 
00734 vector operator-(vector v, const nr_complex_t c) {
00735   vector result (v);
00736   result -= c;
00737   return result;
00738 }
00739 
00740 vector operator-(vector v, const nr_double_t d) {
00741   vector result (v);
00742   result -= d;
00743   return result;
00744 }
00745 
00746 vector operator-(const nr_complex_t c, vector v) {
00747   vector result (-v);
00748   result += c;
00749   return result;
00750 }
00751 
00752 vector operator-(const nr_double_t d, vector v) {
00753   vector result (-v);
00754   result += d;
00755   return result;
00756 }
00757 
00758 vector vector::operator*=(vector v) {
00759   int i, n, len = v.getSize ();
00760   assert (size % len == 0);
00761   for (i = n = 0; i < size; i++) { data[i] *= v (n); if (++n >= len) n = 0; }
00762   return *this;
00763 }
00764 
00765 vector vector::operator*=(const nr_complex_t c) {
00766   for (int i = 0; i < size; i++) data[i] *= c;
00767   return *this;
00768 }
00769 
00770 vector vector::operator*=(const nr_double_t d) {
00771   for (int i = 0; i < size; i++) data[i] *= d;
00772   return *this;
00773 }
00774 
00775 vector operator*(vector v1, vector v2) {
00776   int len1 = v1.getSize (), len2 = v2.getSize ();
00777   vector result;
00778   if (len1 >= len2) {
00779     result  = v1;
00780     result *= v2;
00781   } else {
00782     result  = v2;
00783     result *= v1;
00784   }
00785   return result;
00786 }
00787 
00788 vector operator*(vector v, const nr_complex_t c) {
00789   vector result (v);
00790   result *= c;
00791   return result;
00792 }
00793 
00794 vector operator*(vector v, const nr_double_t d) {
00795   vector result (v);
00796   result *= d;
00797   return result;
00798 }
00799 
00800 vector operator*(const nr_complex_t c, vector v) {
00801   return v * c;
00802 }
00803 
00804 vector operator*(const nr_double_t d, vector v) {
00805   return v * d;
00806 }
00807 
00808 vector vector::operator/=(vector v) {
00809   int i, n, len = v.getSize ();
00810   assert (size % len == 0);
00811   for (i = n = 0; i < size; i++) { data[i] /= v (n); if (++n >= len) n = 0; }
00812   return *this;
00813 }
00814 
00815 vector vector::operator/=(const nr_complex_t c) {
00816   for (int i = 0; i < size; i++) data[i] /= c;
00817   return *this;
00818 }
00819 
00820 vector vector::operator/=(const nr_double_t d) {
00821   for (int i = 0; i < size; i++) data[i] /= d;
00822   return *this;
00823 }
00824 
00825 vector operator/(vector v1, vector v2) {
00826   int len1 = v1.getSize (), len2 = v2.getSize ();
00827   vector result;
00828   if (len1 >= len2) {
00829     assert (len1 % len2 == 0);
00830     result  = v1;
00831     result /= v2;
00832   } else {
00833     assert (len2 % len1 == 0);
00834     result  = 1 / v2;
00835     result *= v1;
00836   }
00837   return result;
00838 }
00839 
00840 vector operator/(vector v, const nr_complex_t c) {
00841   vector result (v);
00842   result /= c;
00843   return result;
00844 }
00845 
00846 vector operator/(vector v, const nr_double_t d) {
00847   vector result (v);
00848   result /= d;
00849   return result;
00850 }
00851 
00852 vector operator/(const nr_complex_t c, vector v) {
00853   vector result (v);
00854   result  = c;
00855   result /= v;
00856   return result;
00857 }
00858 
00859 vector operator/(const nr_double_t d, vector v) {
00860   vector result (v);
00861   result  = d;
00862   result /= v;
00863   return result;
00864 }
00865 
00866 vector operator%(vector v, const nr_complex_t z) {
00867   int len = v.getSize ();
00868   vector result (len);
00869   for (int i = 0; i < len; i++) result (i) = v (i) % z;
00870   return result;
00871 }
00872 
00873 vector operator%(vector v, const nr_double_t d) {
00874   int len = v.getSize ();
00875   vector result (len);
00876   for (int i = 0; i < len; i++) result (i) = v (i) % d;
00877   return result;
00878 }
00879 
00880 vector operator%(const nr_complex_t z, vector v) {
00881   int len = v.getSize ();
00882   vector result (len);
00883   for (int i = 0; i < len; i++) result (i) = z % v (i);
00884   return result;
00885 }
00886 
00887 vector operator%(const nr_double_t d, vector v) {
00888   int len = v.getSize ();
00889   vector result (len);
00890   for (int i = 0; i < len; i++) result (i) = d % v (i);
00891   return result;
00892 }
00893 
00894 vector operator%(vector v1, vector v2) {
00895   int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize ();
00896   if (len1 >= len2) {
00897     assert (len1 % len2 == 0);
00898     len = len1;
00899   } else {
00900     assert (len2 % len1 == 0);
00901     len = len2;
00902   }
00903   vector res (len);
00904   for (j = i = n = 0; n < len; n++) {
00905     res (n) = v1 (i) % v2 (j);
00906     if (++i >= len1) i = 0;  if (++j >= len2) j = 0;
00907   }
00908   return res;
00909 }
00910 
00911 /* This function reverses the order of the data list. */
00912 void vector::reverse (void) {
00913   nr_complex_t * buffer = (nr_complex_t *)
00914     malloc (sizeof (nr_complex_t) * size);
00915   for (int i = 0; i < size; i++) buffer[i] = data[size - 1 - i];
00916   free (data);
00917   data = buffer;
00918   capacity = size;
00919 }
00920 
00921 // Sets the origin (the analysis) of the vector.
00922 void vector::setOrigin (const char * n) {
00923   free (origin);
00924   origin = n ? strdup (n) : NULL;
00925 }
00926 
00927 // Returns the origin (the analysis) of the vector.
00928 char * vector::getOrigin (void) {
00929   return origin;
00930 }
00931 
00932 /* The function returns the number of entries with the given value
00933    deviating no more than the given epsilon. */
00934 int vector::contains (nr_complex_t val, nr_double_t eps) {
00935   int count = 0;
00936   for (int i = 0; i < size; i++) {
00937     if (abs (data[i] - val) <= eps) count++;
00938   }
00939   return count;
00940 }
00941 
00942 // Sorts the vector either in ascending or descending order.
00943 void vector::sort (bool ascending) {
00944   nr_complex_t t;
00945   for (int i = 0; i < size; i++) {
00946     for (int n = 0; n < size - 1; n++) {
00947       if (ascending ? data[n] > data[n+1] : data[n] < data[n+1]) {
00948         t = data[n];
00949         data[n] = data[n+1];
00950         data[n+1] = t;
00951       }
00952     }
00953   }
00954 }
00955 
00956 /* The function creates a linear stepped vector of values starting at
00957    the given start value, ending with the given stop value and
00958    containing points elements. */
00959 vector linspace (nr_double_t start, nr_double_t stop, int points) {
00960   vector result (points);
00961   nr_double_t val, step = (stop - start) / (points - 1);
00962   for (int i = 0; i < points; i++) {
00963     val = start + (i * step);
00964     if (i != 0 && fabs (val) < fabs (step) / 4 && fabs (val) < std::numeric_limits<nr_double_t>::epsilon())
00965       val = 0.0;
00966     result.set (val, i);
00967   }
00968   return result;
00969 }
00970 
00971 /* The function creates a logarithmic stepped vector of values
00972    starting at the given start value, ending with the given stop value
00973    and containing points elements. */
00974 vector logspace (nr_double_t start, nr_double_t stop, int points) {
00975   assert (start * stop > 0);
00976   vector result (points);
00977   nr_double_t step, first, last, d;
00978 
00979   // ensure the last value being larger than the first
00980   if (fabs (start) > fabs (stop)) {
00981     first = fabs (stop);
00982     last = fabs (start);
00983   }
00984   else {
00985     first = fabs (start);
00986     last = fabs (stop);
00987   }
00988   // check direction and sign of values
00989   d = fabs (start) > fabs (stop) ? -1 : 1;
00990   // compute logarithmic step size
00991   step = (::log (last) - ::log (first)) / (points - 1);
00992   for (int i = 0, j = points - 1; i < points; i++, j--) {
00993     if (d > 0)
00994       result.set (start * ::exp (step * i), i);
00995     else
00996       result.set (stop * ::exp (step * i), j);
00997   }
00998   return result;
00999 }
01000 
01001 vector cumsum (vector v) {
01002   vector result (v);
01003   nr_complex_t val (0.0);
01004   for (int i = 0; i < v.getSize (); i++) {
01005     val += v.get (i);
01006     result.set (val, i);
01007   }
01008   return result;
01009 }
01010 
01011 vector cumavg (vector v) {
01012   vector result (v);
01013   nr_complex_t val (0.0);
01014   for (int i = 0; i < v.getSize (); i++) {
01015     val = (val * (nr_double_t) i + v.get (i)) / (i + 1.0);
01016     result.set (val, i);
01017   }
01018   return result;
01019 }
01020 
01021 vector cumprod (vector v) {
01022   vector result (v);
01023   nr_complex_t val (1.0);
01024   for (int i = 0; i < v.getSize (); i++) {
01025     val *= v.get (i);
01026     result.set (val, i);
01027   }
01028   return result;
01029 }
01030 
01031 vector ceil (vector v) {
01032   vector result (v);
01033   for (int i = 0; i < v.getSize (); i++) result.set (ceil (v.get (i)), i);
01034   return result;
01035 }
01036 
01037 vector fix (vector v) {
01038   vector result (v);
01039   for (int i = 0; i < v.getSize (); i++) result.set (fix (v.get (i)), i);
01040   return result;
01041 }
01042 
01043 vector floor (vector v) {
01044   vector result (v);
01045   for (int i = 0; i < v.getSize (); i++) result.set (floor (v.get (i)), i);
01046   return result;
01047 }
01048 
01049 vector round (vector v) {
01050   vector result (v);
01051   for (int i = 0; i < v.getSize (); i++) result.set (round (v.get (i)), i);
01052   return result;
01053 }
01054 
01055 vector sqr (vector v) {
01056   vector result (v);
01057   for (int i = 0; i < v.getSize (); i++) result.set (sqr (v.get (i)), i);
01058   return result;
01059 }
01060 
01061 vector step (vector v) {
01062   vector result (v);
01063   for (int i = 0; i < v.getSize (); i++) result.set (step (v.get (i)), i);
01064   return result;
01065 }
01066 
01067 static nr_double_t integrate_n (vector v) { /* using trapezoidal rule */
01068   nr_double_t result = 0.0;
01069   for (int i = 1; i < v.getSize () - 1; i++) result += norm (v.get (i));
01070   result += 0.5 * norm (v.get (0));
01071   result += 0.5 * norm (v.get (v.getSize () - 1));
01072   return result;
01073 }
01074 
01075 nr_double_t vector::rms (void) {
01076   nr_double_t result = std::sqrt (integrate_n (*this) / getSize ());
01077   return result;
01078 }
01079 
01080 nr_double_t vector::variance (void) {
01081   nr_double_t result = 0.0;
01082   nr_complex_t average = avg (*this);
01083   for (int i = 0; i < getSize (); i++) result += norm (get (i) - average);
01084   if (getSize () > 1)
01085     return result / (getSize () - 1);
01086   return 0.0;
01087 }
01088 
01089 nr_double_t vector::stddev (void) {
01090   return std::sqrt (variance ());
01091 }
01092 
01093 vector erf (vector v) {
01094   vector result (v);
01095   for (int i = 0; i < v.getSize (); i++) result.set (erf (v.get (i)), i);
01096   return result;
01097 }
01098 
01099 vector erfc (vector v) {
01100   vector result (v);
01101   for (int i = 0; i < v.getSize (); i++) result.set (erfc (v.get (i)), i);
01102   return result;
01103 }
01104 
01105 vector erfinv (vector v) {
01106   vector result (v);
01107   for (int i = 0; i < v.getSize (); i++) result.set (erfinv (v.get (i)), i);
01108   return result;
01109 }
01110 
01111 vector erfcinv (vector v) {
01112   vector result (v);
01113   for (int i = 0; i < v.getSize (); i++) result.set (erfcinv (v.get (i)), i);
01114   return result;
01115 }
01116 
01117 vector rad2deg (vector v) {
01118   vector result (v);
01119   for (int i = 0; i < v.getSize (); i++) result.set (rad2deg (v.get (i)), i);
01120   return result;
01121 }
01122 
01123 vector deg2rad (vector v) {
01124   vector result (v);
01125   for (int i = 0; i < v.getSize (); i++) result.set (deg2rad (v.get (i)), i);
01126   return result;
01127 }
01128 
01129 vector i0 (vector v) {
01130   vector result (v);
01131   for (int i = 0; i < v.getSize (); i++) result.set (i0 (v.get (i)), i);
01132   return result;
01133 }
01134 
01135 vector jn (const int n, vector v) {
01136   vector result (v);
01137   for (int i = 0; i < v.getSize (); i++) result.set (jn (n, v.get (i)), i);
01138   return result;
01139 }
01140 
01141 vector yn (const int n, vector v) {
01142   vector result (v);
01143   for (int i = 0; i < v.getSize (); i++) result.set (yn (n, v.get (i)), i);
01144   return result;
01145 }
01146 
01147 vector polar (const nr_complex_t a, vector v) {
01148   vector result (v);
01149   for (int i = 0; i < v.getSize (); i++) result.set (qucs::polar (a, v.get (i)), i);
01150   return result;
01151 }
01152 
01153 vector polar (vector v, const nr_complex_t p) {
01154   vector result (v);
01155   for (int i = 0; i < v.getSize (); i++) result.set (qucs::polar (v.get (i), p), i);
01156   return result;
01157 }
01158 
01159 vector polar (vector a, vector p) {
01160   int j, i, n, len, len1 = a.getSize (), len2 = p.getSize ();
01161   if (len1 >= len2) {
01162     assert (len1 % len2 == 0);
01163     len = len1;
01164   } else {
01165     assert (len2 % len1 == 0);
01166     len = len2;
01167   }
01168   vector res (len);
01169   for (j = i = n = 0; n < len; n++) {
01170     res (n) = qucs::polar (a (i), p (j));
01171     if (++i >= len1) i = 0;  if (++j >= len2) j = 0;
01172   }
01173   return res;
01174 }
01175 
01176 vector atan2 (const nr_double_t y, vector v) {
01177   vector result (v);
01178   for (int i = 0; i < v.getSize (); i++)
01179     result.set (atan2 (y, v.get (i)), i);
01180   return result;
01181 }
01182 
01183 vector atan2 (vector v, const nr_double_t x) {
01184   vector result (v);
01185   for (int i = 0; i < v.getSize (); i++)
01186     result.set (atan2 (v.get (i), x) , i);
01187   return result;
01188 }
01189 
01190 vector atan2 (vector y, vector x) {
01191   int j, i, n, len, len1 = y.getSize (), len2 = x.getSize ();
01192   if (len1 >= len2) {
01193     assert (len1 % len2 == 0);
01194     len = len1;
01195   } else {
01196     assert (len2 % len1 == 0);
01197     len = len2;
01198   }
01199   vector res (len);
01200   for (j = i = n = 0; n < len; n++) {
01201     res (n) = atan2 (y (i), x (j));
01202     if (++i >= len1) i = 0; if (++j >= len2) j = 0;
01203   }
01204   return res;
01205 }
01206 
01207 vector w2dbm (vector v) {
01208   vector result (v);
01209   for (int i = 0; i < v.getSize (); i++)
01210     result.set (10.0 * log10 (v.get (i) / 0.001), i);
01211   return result;
01212 }
01213 
01214 vector dbm2w (vector v) {
01215   vector result (v);
01216   for (int i = 0; i < v.getSize (); i++)
01217     result.set (0.001 * pow (10.0 , v.get (i) / 10.0), i);
01218   return result;
01219 }
01220 
01221 nr_double_t integrate (vector v, const nr_double_t h) {
01222   nr_double_t s = real (v.get (0) ) / 2;
01223   for (int i = 1; i < v.getSize () - 1; i++)
01224     s += real (v.get (i));
01225   return (s + real (v.get (v.getSize () - 1) ) / 2) * h;
01226 }
01227 
01228 nr_complex_t integrate (vector v, const nr_complex_t h) {
01229   nr_complex_t s;
01230   s = v.get (0) / 2.0;
01231   for (int i = 1; i < v.getSize () - 1; i++)
01232     s += v.get (i);
01233   return (s + v.get (v.getSize () - 1) / 2.0) * h;
01234 }
01235 
01236 vector dbm (vector v, const nr_complex_t z) {
01237   vector result (v);
01238   for (int i = 0; i < v.getSize (); i++)
01239     result.set (10.0 * log10 (norm (v.get (i)) / conj (z) / 0.001), i);
01240   return result;
01241 }
01242 
01243 vector runavg (const nr_complex_t x, const int n) {
01244   vector result (n);
01245   for (int i = 0; i < n; i++) result.set (x, i);
01246   return result;
01247 }
01248 
01249 vector runavg (vector v, const int n) {
01250   nr_complex_t s (0.0), y;
01251   int len = v.getSize () - n + 1, i;
01252   vector result (len);
01253   for (i = 0; i < n; i++) s += v.get (i);
01254   y = s / (nr_double_t) n; // first running average value
01255   result.set (y, 0);
01256   for (i = 0; i < len - 1; i++) {
01257     y += (v.get (i + n) - v.get (i)) / (nr_double_t) n;
01258     result.set (y, i + 1);
01259   }
01260   return result;
01261 }
01262 
01263 #ifdef DEBUG
01264 // Debug function: Prints the vector object.
01265 void vector::print (void) {
01266   for (int r = 0; r < size; r++) {
01267     fprintf (stderr, "%+.2e%+.2ei\n", (double) real (get (r)),
01268              (double) imag (get (r)));
01269   }
01270 }
01271 #endif /* DEBUG */
01272 
01273 } // namespace qucs