Branch data Line data Source code
1 : : /*
2 : : * real.cpp - some real valued function implementations
3 : : *
4 : : * Copyright (C) 2008 Stefan Jahn <stefan@lkcc.org>
5 : : * Copyright (C) 2014 Guilheme Brondani Torri <guitorri@gmail.com>
6 : : *
7 : : * This is free software; you can redistribute it and/or modify
8 : : * it under the terms of the GNU General Public License as published by
9 : : * the Free Software Foundation; either version 2, or (at your option)
10 : : * any later version.
11 : : *
12 : : * This software is distributed in the hope that it will be useful,
13 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : : * GNU General Public License for more details.
16 : : *
17 : : * You should have received a copy of the GNU General Public License
18 : : * along with this package; see the file COPYING. If not, write to
19 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
20 : : * Boston, MA 02110-1301, USA.
21 : : *
22 : : * $Id$
23 : : *
24 : : */
25 : :
26 : : #if HAVE_CONFIG_H
27 : : # include <config.h>
28 : : #endif
29 : :
30 : : //#include <cstdlib>
31 : : #include <cmath>
32 : : #include <cassert>
33 : :
34 : : #include "consts.h"
35 : : #include "real.h"
36 : :
37 : : namespace qucs {
38 : :
39 : : //
40 : : // trigonometric
41 : : //
42 : :
43 : : /*! \brief Compute cosine of an angle
44 : : \param[in] z angle in radians
45 : : \return cosine of z
46 : : */
47 : 8880 : nr_double_t cos (const nr_double_t arg) {
48 : 8880 : return std::cos (arg);
49 : : }
50 : :
51 : : /*! \brief Compute sine of an angle
52 : : \param[in] z angle in radians
53 : : \return sine of z
54 : : */
55 : 9472 : nr_double_t sin (const nr_double_t arg) {
56 : 9472 : return std::sin (arg);
57 : : }
58 : :
59 : : /*! \brief Compute tangent of an angle
60 : : \param[in] z angle in radians
61 : : \return tangent of z
62 : : */
63 : 9386 : nr_double_t tan (const nr_double_t arg) {
64 : 9386 : return std::tan (arg);
65 : : }
66 : :
67 : : /*! \brief Compute arc cosine
68 : : \param[in] z arc
69 : : \return arc cosine of z
70 : : */
71 : 0 : nr_double_t acos (const nr_double_t arg) {
72 : 0 : return std::acos (arg);
73 : : }
74 : :
75 : : /*! \brief Compute arc sine
76 : : \param[in] z arc
77 : : \return arc sine of z
78 : : */
79 : 0 : nr_double_t asin (const nr_double_t arg) {
80 : 0 : return std::asin (arg);
81 : : }
82 : :
83 : : /*! \brief Compute arc tangent
84 : : \param[in] z arc
85 : : \return arc tangent of z
86 : : */
87 : 83020 : nr_double_t atan (const nr_double_t arg) {
88 : 83020 : return std::atan (arg);
89 : : }
90 : :
91 : : /*! \brief Compute arc tangent with two parameters (fortran like function)
92 : : \param[in] x proportion of x-coordinate
93 : : \param[in] y proportion of y-coordinate
94 : : \return principal value of the arc tangent of y/x, expressed in radians.
95 : : */
96 : 0 : nr_double_t atan2 (const nr_double_t x, const nr_double_t y) {
97 : 0 : return std::atan2 (x,y);
98 : : }
99 : :
100 : : //
101 : : // hyperbolic
102 : : //
103 : :
104 : : /*! \brief Compute hyperbolic cosine
105 : : \param[in] z arc
106 : : \return hyperbolic cosine of z
107 : : */
108 : 0 : nr_double_t cosh (const nr_double_t arg) {
109 : 0 : return std::cosh (arg);
110 : : }
111 : :
112 : : /*! \brief Compute hyperbolic sine
113 : : \param[in] z arc
114 : : \return hyperbolic sine of z
115 : : */
116 : 0 : nr_double_t sinh (const nr_double_t arg) {
117 : 0 : return std::sinh (arg);
118 : : }
119 : :
120 : : /*! \brief Compute hyperbolic tangent
121 : : \param[in] z arc
122 : : \return hyperbolic tangent of z
123 : : */
124 : 0 : nr_double_t tanh (const nr_double_t arg) {
125 : 0 : return std::tanh (arg);
126 : : }
127 : :
128 : : /*! \brief Compute arc hyperbolic cosine
129 : : \param[in] z arc
130 : : \return arc hyperbolic cosine of z
131 : : */
132 : 0 : nr_double_t acosh (const nr_double_t arg) {
133 : : #ifdef HAVE_STD_ACOSH
134 : : // c++11
135 : 0 : return std::acosh (arg);
136 : : #elif HAVE_ACOSH
137 : : return ::acosh (arg);
138 : : #else
139 : : return log (arg + sqrt (arg * arg - 1.0));
140 : : #endif
141 : : }
142 : :
143 : : /*! \brief Compute arc hyperbolic sine
144 : : \param[in] z arc
145 : : \return arc hyperbolic sine of z
146 : : */
147 : 0 : nr_double_t asinh (const nr_double_t arg)
148 : : {
149 : : #ifdef HAVE_STD_ASINH
150 : : // c++11
151 : 0 : return std::asinh (arg);
152 : : #elif HAVE_ASINH
153 : : return ::asinh (arg);
154 : : #else
155 : : return log (arg + sqrt (arg * arg + 1.0));
156 : : #endif
157 : : }
158 : :
159 : : /*! \brief Compute arc hyperbolic tangent
160 : : \param[in] z arc
161 : : \return arc hyperbolic tangent of z
162 : : */
163 : 0 : nr_double_t atanh (const nr_double_t arg)
164 : : {
165 : : #ifdef HAVE_STD_ATANH
166 : : // c++11
167 : 0 : return std::atanh (arg);
168 : : #elif HAVE_ATANH
169 : : return ::atanh (arg);
170 : : #else
171 : : return 0.5 * log ( 2.0 / (1.0 - arg) - 1.0);
172 : : #endif
173 : : }
174 : :
175 : :
176 : : //
177 : : // exponential and logarithmic functions
178 : : //
179 : 4572028 : nr_double_t exp (const nr_double_t arg) {
180 : 4572028 : return std::exp (arg);
181 : : }
182 : 3935996 : nr_double_t log (const nr_double_t arg) {
183 : 3935996 : return std::log (arg);
184 : : }
185 : 140 : nr_double_t log10 (const nr_double_t arg) {
186 : 140 : return std::log10 (arg);
187 : : }
188 : :
189 : : //
190 : : // power functions
191 : : //
192 : :
193 : 220513 : nr_double_t pow (const nr_double_t a, const nr_double_t b)
194 : : {
195 : 220513 : return std::pow (a,b);
196 : : }
197 : :
198 : 406833 : nr_double_t sqrt (const nr_double_t d) {
199 : 406833 : return std::sqrt (d);
200 : : }
201 : :
202 : : /*!\brief Euclidean distance function
203 : :
204 : : The xhypot() function returns \f$\sqrt{a^2+b^2}\f$.
205 : : This is the length of the hypotenuse of a right-angle triangle with sides
206 : : of length a and b, or the distance
207 : : of the point (a,b) from the origin.
208 : :
209 : : \param[in] a first length
210 : : \param[in] b second length
211 : : \return Euclidean distance from (0,0) to (a,b): \f$\sqrt{a^2+b^2}\f$
212 : : */
213 : 0 : nr_double_t xhypot (const nr_double_t a, const nr_double_t b) {
214 : : #ifdef HAVE_STD_HYPOT
215 : : return std::hypot(a,b) // c++11
216 : : #else
217 : 0 : nr_double_t c = fabs (a);
218 : 0 : nr_double_t d = fabs (b);
219 [ # # ]: 0 : if (c > d) {
220 : 0 : nr_double_t e = d / c;
221 : 0 : return c * sqrt (1 + e * e);
222 : : }
223 [ # # ]: 0 : else if (d == 0)
224 : 0 : return 0;
225 : : else {
226 : 0 : nr_double_t e = c / d;
227 : 0 : return d * sqrt (1 + e * e);
228 : : }
229 : : #endif
230 : : }
231 : :
232 : : //
233 : : // error functions
234 : : //
235 : :
236 : 0 : nr_double_t erf( nr_double_t arg) {
237 : : #ifdef HAVE_STD_ERF
238 : 0 : return std::erf (arg); // c++11
239 : : #elif HAVE_ERF
240 : : return ::erf (arg);
241 : : #endif
242 : : }
243 : :
244 : :
245 : : //
246 : : // rounding and remainder functions
247 : : //
248 : 0 : nr_double_t ceil( nr_double_t arg) {
249 : 0 : return std::ceil(arg);
250 : : }
251 : :
252 : 126104 : nr_double_t floor( nr_double_t arg) {
253 : 126104 : return std::floor(arg);
254 : : }
255 : :
256 : 0 : nr_double_t fmod( nr_double_t arg) {
257 : : #ifdef HAVE_STD_TRUNC
258 : : return std::fmod(arg);
259 : : #else
260 : 0 : return fmod(arg);
261 : : #endif
262 : : }
263 : :
264 : 0 : nr_double_t trunc( nr_double_t arg) {
265 : : #ifdef HAVE_STD_TRUNC
266 : : return qucs::trunc(arg);
267 : : #elif HAVE_TRUNC
268 : 0 : return ::trunc (arg);
269 : : #else
270 : : return arg > 0 ? floor (arg) : floor (arg + 1);
271 : : #endif
272 : : }
273 : 0 : nr_double_t round( nr_double_t arg) {
274 : : #ifdef HAVE_STD_ROUND
275 : 0 : return qucs::round(arg);
276 : : #elif HAVE_ROUND
277 : : return ::round (arg);
278 : : #else
279 : : return (arg > 0) ? floor (arg + 0.5) : ceil (arg - 0.5);
280 : : #endif
281 : : }
282 : :
283 : :
284 : : //
285 : : // Qucs extra trigonometric helper
286 : : //
287 : 3168 : nr_double_t coth (const nr_double_t d) {
288 : 3168 : return 1.0 / std::tanh (d);
289 : : }
290 : :
291 : 3168 : nr_double_t sech (const nr_double_t d) {
292 : 3168 : return (1.0 / std::cosh (d));
293 : : }
294 : :
295 : 0 : nr_double_t cosech (const nr_double_t d) {
296 : 0 : return (1.0 / std::sinh (d));
297 : : }
298 : :
299 : :
300 : : //
301 : : // Qucs extra math functions
302 : : //
303 : :
304 : : /*!\brief Square a value
305 : :
306 : : \param[in] r Real number
307 : : \return \f$x^2\f$
308 : : */
309 : 473605 : nr_double_t sqr (const nr_double_t r) {
310 : 473605 : return r * r;
311 : : }
312 : :
313 : 0 : unsigned int sqr (unsigned int r) {
314 : 0 : return r * r;
315 : : }
316 : :
317 : :
318 : :
319 : : /*!\brief Quartic function
320 : :
321 : : \param[in] r Real number
322 : : \return \f$x^4\f$
323 : : */
324 : 11442 : nr_double_t quadr (const nr_double_t r) {
325 : 11442 : return r * r * r * r;
326 : : }
327 : :
328 : :
329 : : //
330 : : // extra math functions
331 : : //
332 : :
333 : : /*!\brief Compute limited exponential
334 : :
335 : : Compute limited exponential:
336 : : \f[
337 : : \begin{cases}
338 : : \exp r & \text{if } r < \text{M\_LIMEXP} \\
339 : : \exp (\text{M\_LIMEXP})\left[1.0 + (r - \text{M\_LIMEXP})\right] &
340 : : \text{else}
341 : : \end{cases}
342 : : \f]
343 : :
344 : : #M_LIMEXP is a constant
345 : : \param[in] r real number
346 : : \return limited exponential of r
347 : : \todo Change limexp(real) limexp(complex) file order
348 : : \todo Document #M_LIMEXP
349 : : */
350 : 0 : nr_double_t limexp (const nr_double_t r) {
351 [ # # ]: 0 : return r < M_LIMEXP ? exp (r) : exp (M_LIMEXP) * (1.0 + (r - M_LIMEXP));
352 : : }
353 : :
354 : : /*!\brief real signum function
355 : :
356 : : compute \f[
357 : : \mathrm{signum}\;d=
358 : : = \begin{cases}
359 : : O & \text{if } d=0 \\
360 : : 1 & \text{if } d>0 \\
361 : : -1 & \text{if } d<0
362 : : \end{cases}
363 : : \f]
364 : : \param[in] d real number
365 : : \return signum of d
366 : : \todo Move near complex signum
367 : : */
368 : 0 : nr_double_t signum (const nr_double_t d) {
369 [ # # ]: 0 : if (d == 0) return 0;
370 [ # # ]: 0 : return d < 0 ? -1 : 1;
371 : : }
372 : :
373 : : /*!\brief real sign function
374 : :
375 : : compute \f[
376 : : \mathrm{sign}\;d=
377 : : = \begin{cases}
378 : : 1 & \text{if } d\ge 0 \\
379 : : -1 & \text{if } d<0
380 : : \end{cases}
381 : : \f]
382 : : \param[in] d real number
383 : : \return sign of d
384 : : \todo Move near complex sign
385 : : */
386 : 0 : nr_double_t sign (const nr_double_t d) {
387 [ # # ]: 0 : return d < 0 ? -1 : 1;
388 : : }
389 : :
390 : : /*!\brief Real cardinal sinus
391 : :
392 : : Compute \f$\mathrm{sinc}\;d=\frac{\sin d}{d}\f$
393 : : \param[in] d real number
394 : : \return cardianal sinus of s
395 : : \todo Why not inline
396 : : */
397 : 0 : nr_double_t sinc (const nr_double_t d) {
398 [ # # ]: 0 : if (d == 0) return 1;
399 : 0 : return sin (d) / d;
400 : : }
401 : :
402 : : /*!\brief Fix function
403 : :
404 : : Fix is nearest integral value in direction of 0,
405 : : \f[
406 : : \operatorname{fix} d=\begin{cases}
407 : : \operatorname{floor} d & \text{if } d > 0 \\
408 : : \operatorname{ceil} d & \text{else}
409 : : \end{cases}
410 : : \f]
411 : :
412 : : \param[in] d real number
413 : : \return fixed complex number
414 : : \todo Why not inline?
415 : : */
416 : 0 : nr_double_t fix (const nr_double_t d) {
417 [ # # ]: 0 : return (d > 0) ? floor (d) : ceil (d);
418 : : }
419 : :
420 : : /*!\brief Heaviside step function
421 : :
422 : : The Heaviside step function, H, also called unit step function,
423 : : is a discontinuous function whose value is zero for negative argument and
424 : : one for positive argument. For zero by convention, H(0)=0.5
425 : : \param[in] d Heaviside argument
426 : : \return Heaviside step
427 : : \todo Create Heaviside alias
428 : : */
429 : 0 : nr_double_t step (const nr_double_t d) {
430 : 0 : nr_double_t x = d;
431 [ # # ]: 0 : if (x < 0.0)
432 : 0 : x = 0.0;
433 [ # # ]: 0 : else if (x > 0.0)
434 : 0 : x = 1.0;
435 : : else
436 : 0 : x = 0.5;
437 : 0 : return x;
438 : : }
439 : :
440 : : /*!\brief Compute factorial n ie \$n!\$
441 : :
442 : : */
443 : : unsigned int
444 : 0 : factorial (unsigned int n) {
445 : 0 : unsigned int result = 1;
446 : :
447 : : /* 13! > 2^32 */
448 [ # # ]: 0 : assert (n < 13);
449 : :
450 [ # # ]: 0 : if (n == 0)
451 : 0 : return 1;
452 : :
453 [ # # ]: 0 : for (; n > 1; n--)
454 : 0 : result = result * n;
455 : :
456 : 0 : return result;
457 : : }
458 : :
459 : :
460 : : //
461 : : // overload complex manipulations on reals
462 : : //
463 : :
464 : :
465 : : /*!\brief Real part of real number
466 : :
467 : : \param[in] r Real number
468 : : \return Real part of r ie r
469 : : */
470 : 20554106 : nr_double_t real (const nr_double_t r) {
471 : 20554106 : return r;
472 : : }
473 : :
474 : : /*!\brief Imaginary part of complex number
475 : :
476 : : \param[in] r Real number
477 : : \return Imaginary part of r
478 : : */
479 : 136130 : nr_double_t imag (const nr_double_t r) {
480 : 136130 : return 0.0;
481 : : }
482 : :
483 : : /*!\brief Compute euclidian norm of real number
484 : :
485 : : Compute \f$r^2\f$
486 : : \param[in] r Real number
487 : : \return Euclidian norm of r
488 : : */
489 : 0 : nr_double_t norm (const nr_double_t r) {
490 : 0 : return r * r;
491 : : }
492 : :
493 : : /*!\brief Compute complex modulus of real number
494 : :
495 : : \param[in] r Real number
496 : : \return Modulus of r
497 : : */
498 : 99883153 : nr_double_t abs (const nr_double_t r) {
499 : 99883153 : return std::abs (r);
500 : : }
501 : :
502 : : /*!\brief Conjugate of real number
503 : :
504 : : \param[in] r Real number
505 : : \return Conjugate of real r ie r
506 : : */
507 : 0 : nr_double_t conj (const nr_double_t r) {
508 : 0 : return r;
509 : : }
510 : :
511 : : } // namespace qucs
|
