GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
lo-mappers.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2013 John W. Eaton
4 Copyright (C) 2010 VZLU Prague
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include <cfloat>
29 
30 #include "lo-error.h"
31 #include "lo-ieee.h"
32 #include "lo-mappers.h"
33 #include "lo-math.h"
34 #include "lo-specfun.h"
35 #include "lo-utils.h"
36 #include "oct-cmplx.h"
37 
38 #include "f77-fcn.h"
39 
40 // double -> double mappers.
41 
42 // Both xtrunc and xround belong here so we can keep gnulib:: out of
43 // lo-mappers.h.
44 
45 double
46 xtrunc (double x)
47 {
48  return gnulib::trunc (x);
49 }
50 
51 double
52 xcopysign (double x, double y)
53 {
54  return gnulib::copysign (x, y);
55 }
56 
57 double xfloor (double x)
58 {
59  return gnulib::floor (x);
60 }
61 
62 double
63 xround (double x)
64 {
65  return gnulib::round (x);
66 }
67 
68 double
69 xroundb (double x)
70 {
71  double t = xround (x);
72 
73  if (fabs (x - t) == 0.5)
74  t = 2 * xtrunc (0.5 * t);
75 
76  return t;
77 }
78 
79 double
80 signum (double x)
81 {
82  double tmp = 0.0;
83 
84  if (x < 0.0)
85  tmp = -1.0;
86  else if (x > 0.0)
87  tmp = 1.0;
88 
89  return xisnan (x) ? octave_NaN : tmp;
90 }
91 
92 double
93 xlog2 (double x)
94 {
95 #if defined (HAVE_LOG2)
96  return log2 (x);
97 #else
98 #if defined (M_LN2)
99  static double ln2 = M_LN2;
100 #else
101  static double ln2 = log (2);
102 #endif
103 
104  return log (x) / ln2;
105 #endif
106 }
107 
108 Complex
109 xlog2 (const Complex& x)
110 {
111 #if defined (M_LN2)
112  static double ln2 = M_LN2;
113 #else
114  static double ln2 = log (2);
115 #endif
116 
117  return std::log (x) / ln2;
118 }
119 
120 double
121 xexp2 (double x)
122 {
123 #if defined (HAVE_EXP2)
124  return exp2 (x);
125 #else
126 #if defined (M_LN2)
127  static double ln2 = M_LN2;
128 #else
129  static double ln2 = log (2);
130 #endif
131 
132  return exp (x * ln2);
133 #endif
134 }
135 
136 double
137 xlog2 (double x, int& exp)
138 {
139  return gnulib::frexp (x, &exp);
140 }
141 
142 Complex
143 xlog2 (const Complex& x, int& exp)
144 {
145  double ax = std::abs (x);
146  double lax = xlog2 (ax, exp);
147  return (ax != lax) ? (x / ax) * lax : x;
148 }
149 
150 // double -> bool mappers.
151 
152 #if ! defined(HAVE_CMATH_ISNAN)
153 bool
154 xisnan (double x)
155 {
156  return lo_ieee_isnan (x);
157 }
158 #endif
159 
160 #if ! defined(HAVE_CMATH_ISFINITE)
161 bool
162 xfinite (double x)
163 {
164  return lo_ieee_finite (x);
165 }
166 #endif
167 
168 #if ! defined(HAVE_CMATH_ISINF)
169 bool
170 xisinf (double x)
171 {
172  return lo_ieee_isinf (x);
173 }
174 #endif
175 
176 bool
177 octave_is_NA (double x)
178 {
179  return lo_ieee_is_NA (x);
180 }
181 
182 bool
184 {
185  return lo_ieee_isnan (x);
186 }
187 
188 // (double, double) -> double mappers.
189 
190 // complex -> complex mappers.
191 
192 Complex
193 acos (const Complex& x)
194 {
195  static Complex i (0, 1);
196 
197  return -i * (log (x + i * (sqrt (1.0 - x*x))));
198 }
199 
200 Complex
201 acosh (const Complex& x)
202 {
203  return log (x + sqrt (x*x - 1.0));
204 }
205 
206 Complex
207 asin (const Complex& x)
208 {
209  static Complex i (0, 1);
210 
211  return -i * log (i*x + sqrt (1.0 - x*x));
212 }
213 
214 Complex
215 asinh (const Complex& x)
216 {
217  return log (x + sqrt (x*x + 1.0));
218 }
219 
220 Complex
221 atan (const Complex& x)
222 {
223  static Complex i (0, 1);
224 
225  return i * log ((i + x) / (i - x)) / 2.0;
226 }
227 
228 Complex
229 atanh (const Complex& x)
230 {
231  return log ((1.0 + x) / (1.0 - x)) / 2.0;
232 }
233 
234 // complex -> bool mappers.
235 
236 bool
238 {
239  return (octave_is_NA (real (x)) || octave_is_NA (imag (x)));
240 }
241 
242 bool
244 {
245  return (xisnan (real (x)) || xisnan (imag (x)));
246 }
247 
248 // (complex, complex) -> complex mappers.
249 
250 // FIXME: need to handle NA too?
251 
252 Complex
253 xmin (const Complex& x, const Complex& y)
254 {
255  return abs (x) <= abs (y) ? x : (xisnan (x) ? x : y);
256 }
257 
258 Complex
259 xmax (const Complex& x, const Complex& y)
260 {
261  return abs (x) >= abs (y) ? x : (xisnan (x) ? x : y);
262 }
263 
264 
265 // float -> float mappers.
266 
267 // Both xtrunc and xround belong here so we can keep gnulib:: out of
268 // lo-mappers.h.
269 
270 float
271 xtrunc (float x)
272 {
273  return gnulib::truncf (x);
274 }
275 
276 float
277 xcopysign (float x, float y)
278 {
279  return gnulib::copysignf (x, y);
280 }
281 
282 float xfloor (float x)
283 {
284  return gnulib::floorf (x);
285 }
286 
287 float
288 xround (float x)
289 {
290  return gnulib::roundf (x);
291 }
292 
293 float
294 xroundb (float x)
295 {
296  float t = xround (x);
297 
298  if (fabsf (x - t) == 0.5)
299  t = 2 * xtrunc (0.5 * t);
300 
301  return t;
302 }
303 
304 float
305 signum (float x)
306 {
307  float tmp = 0.0;
308 
309  if (x < 0.0)
310  tmp = -1.0;
311  else if (x > 0.0)
312  tmp = 1.0;
313 
314  return xisnan (x) ? octave_Float_NaN : tmp;
315 }
316 
317 float
318 xlog2 (float x)
319 {
320 #if defined (HAVE_LOG2F)
321  return log2f (x);
322 #elif defined (HAVE_LOG2)
323  return log2 (x);
324 #else
325 #if defined (M_LN2)
326  static float ln2 = M_LN2;
327 #else
328  static float ln2 = log2 (2);
329 #endif
330 
331  return log (x) / ln2;
332 #endif
333 }
334 
336 xlog2 (const FloatComplex& x)
337 {
338 #if defined (M_LN2)
339  static float ln2 = M_LN2;
340 #else
341  static float ln2 = log (2);
342 #endif
343 
344  return std::log (x) / ln2;
345 }
346 
347 float
348 xexp2 (float x)
349 {
350 #if defined (HAVE_EXP2F)
351  return exp2f (x);
352 #elif defined (HAVE_EXP2)
353  return exp2 (x);
354 #else
355 #if defined (M_LN2)
356  static float ln2 = M_LN2;
357 #else
358  static float ln2 = log2 (2);
359 #endif
360 
361  return exp (x * ln2);
362 #endif
363 }
364 
365 float
366 xlog2 (float x, int& exp)
367 {
368  return gnulib::frexpf (x, &exp);
369 }
370 
372 xlog2 (const FloatComplex& x, int& exp)
373 {
374  float ax = std::abs (x);
375  float lax = xlog2 (ax, exp);
376  return (ax != lax) ? (x / ax) * lax : x;
377 }
378 
379 // float -> bool mappers.
380 
381 #if ! defined(HAVE_CMATH_ISNANF)
382 bool
383 xisnan (float x)
384 {
385  return lo_ieee_isnan (x);
386 }
387 #endif
388 
389 #if ! defined(HAVE_CMATH_ISFINITEF)
390 bool
391 xfinite (float x)
392 {
393  return lo_ieee_finite (x);
394 }
395 #endif
396 
397 #if ! defined(HAVE_CMATH_ISINFF)
398 bool
399 xisinf (float x)
400 {
401  return lo_ieee_isinf (x);
402 }
403 #endif
404 
405 bool
406 octave_is_NA (float x)
407 {
408  return lo_ieee_is_NA (x);
409 }
410 
411 bool
413 {
414  return lo_ieee_isnan (x);
415 }
416 
417 // (float, float) -> float mappers.
418 
419 // complex -> complex mappers.
420 
422 acos (const FloatComplex& x)
423 {
424  static FloatComplex i (0, 1);
425 
426  return -i * (log (x + i * (sqrt (static_cast<float>(1.0) - x*x))));
427 }
428 
430 acosh (const FloatComplex& x)
431 {
432  return log (x + sqrt (x*x - static_cast<float>(1.0)));
433 }
434 
436 asin (const FloatComplex& x)
437 {
438  static FloatComplex i (0, 1);
439 
440  return -i * log (i*x + sqrt (static_cast<float>(1.0) - x*x));
441 }
442 
444 asinh (const FloatComplex& x)
445 {
446  return log (x + sqrt (x*x + static_cast<float>(1.0)));
447 }
448 
450 atan (const FloatComplex& x)
451 {
452  static FloatComplex i (0, 1);
453 
454  return i * log ((i + x) / (i - x)) / static_cast<float>(2.0);
455 }
456 
458 atanh (const FloatComplex& x)
459 {
460  return log ((static_cast<float>(1.0) + x) / (static_cast<float>
461  (1.0) - x)) / static_cast<float>(2.0);
462 }
463 
464 // complex -> bool mappers.
465 
466 bool
467 octave_is_NA (const FloatComplex& x)
468 {
469  return (octave_is_NA (real (x)) || octave_is_NA (imag (x)));
470 }
471 
472 bool
474 {
475  return (xisnan (real (x)) || xisnan (imag (x)));
476 }
477 
478 // (complex, complex) -> complex mappers.
479 
480 // FIXME: need to handle NA too?
481 
483 xmin (const FloatComplex& x, const FloatComplex& y)
484 {
485  return abs (x) <= abs (y) ? x : (xisnan (x) ? x : y);
486 }
487 
489 xmax (const FloatComplex& x, const FloatComplex& y)
490 {
491  return abs (x) >= abs (y) ? x : (xisnan (x) ? x : y);
492 }
493 
494 Complex
495 rc_acos (double x)
496 {
497  return fabs (x) > 1.0 ? acos (Complex (x)) : Complex (acos (x));
498 }
499 
501 rc_acos (float x)
502 {
503  return fabsf (x) > 1.0f ? acos (FloatComplex (x)) : FloatComplex (acosf (x));
504 }
505 
506 Complex
507 rc_acosh (double x)
508 {
509  return x < 1.0 ? acosh (Complex (x)) : Complex (acosh (x));
510 }
511 
513 rc_acosh (float x)
514 {
515  return x < 1.0f ? acosh (FloatComplex (x)) : FloatComplex (acoshf (x));
516 }
517 
518 Complex
519 rc_asin (double x)
520 {
521  return fabs (x) > 1.0 ? asin (Complex (x)) : Complex (asin (x));
522 }
523 
525 rc_asin (float x)
526 {
527  return fabsf (x) > 1.0f ? asin (FloatComplex (x)) : FloatComplex (asinf (x));
528 }
529 
530 Complex
531 rc_atanh (double x)
532 {
533  return fabs (x) > 1.0 ? atanh (Complex (x)) : Complex (atanh (x));
534 }
535 
537 rc_atanh (float x)
538 {
539  return fabsf (x) > 1.0f ? atanh (FloatComplex (x))
540  : FloatComplex (atanhf (x));
541 }
542 
543 Complex
544 rc_log (double x)
545 {
546  const double pi = 3.14159265358979323846;
547  return x < 0.0 ? Complex (log (-x), pi) : Complex (log (x));
548 }
549 
551 rc_log (float x)
552 {
553  const float pi = 3.14159265358979323846f;
554  return x < 0.0f ? FloatComplex (logf (-x), pi) : FloatComplex (logf (x));
555 }
556 
557 Complex
558 rc_log2 (double x)
559 {
560  const double pil2 = 4.53236014182719380962; // = pi / log(2)
561  return x < 0.0 ? Complex (xlog2 (-x), pil2) : Complex (xlog2 (x));
562 }
563 
565 rc_log2 (float x)
566 {
567  const float pil2 = 4.53236014182719380962f; // = pi / log(2)
568  return x < 0.0f ? FloatComplex (xlog2 (-x), pil2) : FloatComplex (xlog2 (x));
569 }
570 
571 Complex
572 rc_log10 (double x)
573 {
574  const double pil10 = 1.36437635384184134748; // = pi / log(10)
575  return x < 0.0 ? Complex (log10 (-x), pil10) : Complex (log10 (x));
576 }
577 
579 rc_log10 (float x)
580 {
581  const float pil10 = 1.36437635384184134748f; // = pi / log(10)
582  return x < 0.0f ? FloatComplex (log10 (-x), pil10)
583  : FloatComplex (log10f (x));
584 }
585 
586 Complex
587 rc_sqrt (double x)
588 {
589  return x < 0.0 ? Complex (0.0, sqrt (-x)) : Complex (sqrt (x));
590 }
591 
593 rc_sqrt (float x)
594 {
595  return x < 0.0f ? FloatComplex (0.0f, sqrtf (-x)) : FloatComplex (sqrtf (x));
596 }
597 
598 bool
599 xnegative_sign (double x)
600 {
601  return __lo_ieee_signbit (x);
602 }
603 
604 bool
605 xnegative_sign (float x)
606 {
607  return __lo_ieee_float_signbit (x);
608 }
609 
610 // Convert X to the nearest integer value. Should not pass NaN to
611 // this function.
612 
613 // Sometimes you need a large integer, but not always.
614 
616 NINTbig (double x)
617 {
622  else
623  return static_cast<octave_idx_type> ((x > 0) ? (x + 0.5) : (x - 0.5));
624 }
625 
627 NINTbig (float x)
628 {
633  else
634  return static_cast<octave_idx_type> ((x > 0) ? (x + 0.5) : (x - 0.5));
635 }
636 
637 int
638 NINT (double x)
639 {
642  else if (x < std::numeric_limits<int>::min ())
644  else
645  return static_cast<int> ((x > 0) ? (x + 0.5) : (x - 0.5));
646 }
647 
648 int
649 NINT (float x)
650 {
653  else if (x < std::numeric_limits<int>::min ())
655  else
656  return static_cast<int> ((x > 0) ? (x + 0.5) : (x - 0.5));
657 }