Mercurial > octave

comparison liboctave/numeric/lo-specfun.cc @ 21231:5f318c8ec634

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
eliminate feature tests from lo-specfun.h * lo-specfun.h, lo-specfun.cc (xacosh, xasinh, xatanh, xerf, xerfc xexpm1, xlog1p, xcbrt): Rename to have 'x' prefix. Conditionally define in .cc file. Change all uses Move complex versions of acosh, asinh, and atanh functions here.
author John W. Eaton <jwe@octave.org>
date Tue, 09 Feb 2016 04:15:50 -0500
parents f7121e111991
children 40de9f8f23a6
comparison
equal deleted inserted replaced
21230:721539013db4 21231:5f318c8ec634
196 F77_FUNC (dpsifn, DPSIFN) (const double*, const octave_idx_type&, 196 F77_FUNC (dpsifn, DPSIFN) (const double*, const octave_idx_type&,
197 const octave_idx_type&, const octave_idx_type&, 197 const octave_idx_type&, const octave_idx_type&,
198 double*, octave_idx_type*, octave_idx_type*); 198 double*, octave_idx_type*, octave_idx_type*);
199 } 199 }
200 200
201 #if ! defined (HAVE_ACOSH)
202 double 201 double
203 acosh (double x) 202 xacosh (double x)
204 { 203 {
204 #if defined (HAVE_ACOSH)
205 return acosh (x);
206 #else
205 double retval; 207 double retval;
206 F77_XFCN (xdacosh, XDACOSH, (x, retval)); 208 F77_XFCN (xdacosh, XDACOSH, (x, retval));
207 return retval; 209 return retval;
208 }
209 #endif 210 #endif
210 211 }
211 #if ! defined (HAVE_ACOSHF) 212
212 float 213 float
213 acoshf (float x) 214 xacosh (float x)
214 { 215 {
216 #if defined (HAVE_ACOSHF)
217 return acoshf (x);
218 #else
215 float retval; 219 float retval;
216 F77_XFCN (xacosh, XACOSH, (x, retval)); 220 F77_XFCN (xacosh, XACOSH, (x, retval));
217 return retval; 221 return retval;
218 }
219 #endif 222 #endif
220 223 }
221 #if ! defined (HAVE_ASINH) 224
225 Complex
226 xacosh (const Complex& x)
227 {
228 return log (x + sqrt (x + 1.0) * sqrt (x - 1.0));
229 }
230
231 FloatComplex
232 xacosh (const FloatComplex& x)
233 {
234 return log (x + sqrt (x + 1.0f) * sqrt (x - 1.0f));
235 }
236
222 double 237 double
223 asinh (double x) 238 xasinh (double x)
224 { 239 {
240 #if defined (HAVE_ASINH)
241 return asinh (x);
242 #else
225 double retval; 243 double retval;
226 F77_XFCN (xdasinh, XDASINH, (x, retval)); 244 F77_XFCN (xdasinh, XDASINH, (x, retval));
227 return retval; 245 return retval;
228 }
229 #endif 246 #endif
230 247 }
231 #if ! defined (HAVE_ASINHF) 248
232 float 249 float
233 asinhf (float x) 250 xasinh (float x)
234 { 251 {
252 #if defined (HAVE_ASINHF)
253 return asinhf (x);
254 #else
235 float retval; 255 float retval;
236 F77_XFCN (xasinh, XASINH, (x, retval)); 256 F77_XFCN (xasinh, XASINH, (x, retval));
237 return retval; 257 return retval;
238 }
239 #endif 258 #endif
240 259 }
241 #if ! defined (HAVE_ATANH) 260
261 Complex
262 xasinh (const Complex& x)
263 {
264 return log (x + sqrt (x*x + 1.0));
265 }
266
267 FloatComplex
268 xasinh (const FloatComplex& x)
269 {
270 return log (x + sqrt (x*x + 1.0f));
271 }
272
242 double 273 double
243 atanh (double x) 274 xatanh (double x)
244 { 275 {
276 #if defined (HAVE_ATANH)
277 return atanh (x);
278 #else
245 double retval; 279 double retval;
246 F77_XFCN (xdatanh, XDATANH, (x, retval)); 280 F77_XFCN (xdatanh, XDATANH, (x, retval));
247 return retval; 281 return retval;
248 }
249 #endif 282 #endif
250 283 }
251 #if ! defined (HAVE_ATANHF) 284
252 float 285 float
253 atanhf (float x) 286 xatanh (float x)
254 { 287 {
288 #if defined (HAVE_ATANHF)
289 return atanhf (x);
290 #else
255 float retval; 291 float retval;
256 F77_XFCN (xatanh, XATANH, (x, retval)); 292 F77_XFCN (xatanh, XATANH, (x, retval));
257 return retval; 293 return retval;
258 }
259 #endif 294 #endif
260 295 }
261 #if ! defined (HAVE_ERF) 296
297 Complex
298 xatanh (const Complex& x)
299 {
300 return log ((1.0 + x) / (1.0 - x)) / 2.0;
301 }
302
303 FloatComplex
304 xatanh (const FloatComplex& x)
305 {
306 return log ((1.0f + x) / (1.0f - x)) / 2.0f;
307 }
308
262 double 309 double
263 erf (double x) 310 xerf (double x)
264 { 311 {
312 #if defined (HAVE_ERF)
313 return erf (x);
314 #else
265 double retval; 315 double retval;
266 F77_XFCN (xderf, XDERF, (x, retval)); 316 F77_XFCN (xderf, XDERF, (x, retval));
267 return retval; 317 return retval;
268 }
269 #endif 318 #endif
270 319 }
271 #if ! defined (HAVE_ERFF) 320
272 float 321 float
273 erff (float x) 322 xerf (float x)
274 { 323 {
324 #if defined (HAVE_ERFF)
325 return erff (x);
326 #else
275 float retval; 327 float retval;
276 F77_XFCN (xerf, XERF, (x, retval)); 328 F77_XFCN (xerf, XERF, (x, retval));
277 return retval; 329 return retval;
278 }
279 #endif 330 #endif
280 331 }
281 #if ! defined (HAVE_ERFC)
282 double
283 erfc (double x)
284 {
285 double retval;
286 F77_XFCN (xderfc, XDERFC, (x, retval));
287 return retval;
288 }
289 #endif
290
291 #if ! defined (HAVE_ERFCF)
292 float
293 erfcf (float x)
294 {
295 float retval;
296 F77_XFCN (xerfc, XERFC, (x, retval));
297 return retval;
298 }
299 #endif
300 332
301 // Complex error function from the Faddeeva package 333 // Complex error function from the Faddeeva package
302 Complex 334 Complex
303 erf (const Complex& x) 335 xerf (const Complex& x)
304 { 336 {
305 return Faddeeva::erf (x); 337 return Faddeeva::erf (x);
306 } 338 }
339
307 FloatComplex 340 FloatComplex
308 erf (const FloatComplex& x) 341 xerf (const FloatComplex& x)
309 { 342 {
310 Complex xd (real (x), imag (x)); 343 Complex xd (real (x), imag (x));
311 Complex ret = Faddeeva::erf (xd, std::numeric_limits<float>::epsilon ()); 344 Complex ret = Faddeeva::erf (xd, std::numeric_limits<float>::epsilon ());
312 return FloatComplex (real (ret), imag (ret)); 345 return FloatComplex (real (ret), imag (ret));
313 } 346 }
314 347
348 double
349 xerfc (double x)
350 {
351 #if defined (HAVE_ERFC)
352 return erfc (x);
353 #else
354 double retval;
355 F77_XFCN (xderfc, XDERFC, (x, retval));
356 return retval;
357 #endif
358 }
359
360 float
361 xerfc (float x)
362 {
363 #if defined (HAVE_ERFCF)
364 return erfcf (x);
365 #else
366 float retval;
367 F77_XFCN (xerfc, XERFC, (x, retval));
368 return retval;
369 #endif
370 }
371
315 // Complex complementary error function from the Faddeeva package 372 // Complex complementary error function from the Faddeeva package
316 Complex 373 Complex
317 erfc (const Complex& x) 374 xerfc (const Complex& x)
318 { 375 {
319 return Faddeeva::erfc (x); 376 return Faddeeva::erfc (x);
320 } 377 }
378
321 FloatComplex 379 FloatComplex
322 erfc (const FloatComplex& x) 380 xerfc (const FloatComplex& x)
323 { 381 {
324 Complex xd (real (x), imag (x)); 382 Complex xd (real (x), imag (x));
325 Complex ret = Faddeeva::erfc (xd, std::numeric_limits<float>::epsilon ()); 383 Complex ret = Faddeeva::erfc (xd, std::numeric_limits<float>::epsilon ());
326 return FloatComplex (real (ret), imag (ret)); 384 return FloatComplex (real (ret), imag (ret));
327 } 385 }
516 return result + FloatComplex (0., M_PI); 574 return result + FloatComplex (0., M_PI);
517 else 575 else
518 return result; 576 return result;
519 } 577 }
520 578
521 #if ! defined (HAVE_EXPM1)
522 double 579 double
523 expm1 (double x) 580 xexpm1 (double x)
524 { 581 {
582 #if defined (HAVE_EXPM1)
583 return expm1 (x);
584 #else
525 double retval; 585 double retval;
526 586
527 double ax = fabs (x); 587 double ax = fabs (x);
528 588
529 if (ax < 0.1) 589 if (ax < 0.1)
549 } 609 }
550 else 610 else
551 retval = exp (x) - 1; 611 retval = exp (x) - 1;
552 612
553 return retval; 613 return retval;
554 }
555 #endif 614 #endif
615 }
556 616
557 Complex 617 Complex
558 expm1 (const Complex& x) 618 xexpm1 (const Complex& x)
559 { 619 {
560 Complex retval; 620 Complex retval;
561 621
562 if (std:: abs (x) < 1) 622 if (std:: abs (x) < 1)
563 { 623 {
564 double im = x.imag (); 624 double im = x.imag ();
565 double u = expm1 (x.real ()); 625 double u = xexpm1 (x.real ());
566 double v = sin (im/2); 626 double v = sin (im/2);
567 v = -2*v*v; 627 v = -2*v*v;
568 retval = Complex (u*v + u + v, (u+1) * sin (im)); 628 retval = Complex (u*v + u + v, (u+1) * sin (im));
569 } 629 }
570 else 630 else
571 retval = std::exp (x) - Complex (1); 631 retval = std::exp (x) - Complex (1);
572 632
573 return retval; 633 return retval;
574 } 634 }
575 635
576 #if ! defined (HAVE_EXPM1F)
577 float 636 float
578 expm1f (float x) 637 xexpm1 (float x)
579 { 638 {
639 #if defined (HAVE_EXPM1F)
640 return expm1f (x);
641 #else
580 float retval; 642 float retval;
581 643
582 float ax = fabs (x); 644 float ax = fabs (x);
583 645
584 if (ax < 0.1) 646 if (ax < 0.1)
604 } 666 }
605 else 667 else
606 retval = exp (x) - 1; 668 retval = exp (x) - 1;
607 669
608 return retval; 670 return retval;
609 }
610 #endif 671 #endif
672 }
611 673
612 FloatComplex 674 FloatComplex
613 expm1 (const FloatComplex& x) 675 xexpm1 (const FloatComplex& x)
614 { 676 {
615 FloatComplex retval; 677 FloatComplex retval;
616 678
617 if (std:: abs (x) < 1) 679 if (std:: abs (x) < 1)
618 { 680 {
619 float im = x.imag (); 681 float im = x.imag ();
620 float u = expm1 (x.real ()); 682 float u = xexpm1 (x.real ());
621 float v = sin (im/2); 683 float v = sin (im/2);
622 v = -2*v*v; 684 v = -2*v*v;
623 retval = FloatComplex (u*v + u + v, (u+1) * sin (im)); 685 retval = FloatComplex (u*v + u + v, (u+1) * sin (im));
624 } 686 }
625 else 687 else
626 retval = std::exp (x) - FloatComplex (1); 688 retval = std::exp (x) - FloatComplex (1);
627 689
628 return retval; 690 return retval;
629 } 691 }
630 692
631 #if ! defined (HAVE_LOG1P)
632 double 693 double
633 log1p (double x) 694 xlog1p (double x)
634 { 695 {
696 #if defined (HAVE_LOG1P)
697 return log1p (x);
698 #else
635 double retval; 699 double retval;
636 700
637 double ax = fabs (x); 701 double ax = fabs (x);
638 702
639 if (ax < 0.2) 703 if (ax < 0.2)
647 } 711 }
648 else 712 else
649 retval = gnulib::log (1 + x); 713 retval = gnulib::log (1 + x);
650 714
651 return retval; 715 return retval;
652 }
653 #endif 716 #endif
717 }
654 718
655 Complex 719 Complex
656 log1p (const Complex& x) 720 xlog1p (const Complex& x)
657 { 721 {
658 Complex retval; 722 Complex retval;
659 723
660 double r = x.real (), i = x.imag (); 724 double r = x.real (), i = x.imag ();
661 725
662 if (fabs (r) < 0.5 && fabs (i) < 0.5) 726 if (fabs (r) < 0.5 && fabs (i) < 0.5)
663 { 727 {
664 double u = 2*r + r*r + i*i; 728 double u = 2*r + r*r + i*i;
665 retval = Complex (log1p (u / (1+sqrt (u+1))), 729 retval = Complex (xlog1p (u / (1+sqrt (u+1))),
666 atan2 (1 + r, i)); 730 atan2 (1 + r, i));
667 } 731 }
668 else 732 else
669 retval = std::log (Complex (1) + x); 733 retval = std::log (Complex (1) + x);
670 734
671 return retval; 735 return retval;
672 } 736 }
673 737
674 #if ! defined (HAVE_CBRT) 738 template <typename T>
675 double cbrt (double x) 739 T
676 { 740 xxcbrt (T x)
677 static const double one_third = 0.3333333333333333333; 741 {
742 static const T one_third = 0.3333333333333333333f;
678 if (xfinite (x)) 743 if (xfinite (x))
679 { 744 {
680 // Use pow. 745 // Use pow.
681 double y = std::pow (std::abs (x), one_third) * signum (x); 746 T y = std::pow (std::abs (x), one_third) * signum (x);
682 // Correct for better accuracy. 747 // Correct for better accuracy.
683 return (x / (y*y) + y + y) / 3; 748 return (x / (y*y) + y + y) / 3;
684 } 749 }
685 else 750 else
686 return x; 751 return x;
687 } 752 }
753
754 double
755 xcbrt (double x)
756 {
757 #if defined (HAVE_CBRT)
758 return cbrt (x);
759 #else
760 return xxcbrt (x);
688 #endif 761 #endif
689 762 }
690 #if ! defined (HAVE_LOG1PF) 763
691 float 764 float
692 log1pf (float x) 765 xlog1p (float x)
693 { 766 {
767 #if defined (HAVE_LOG1PF)
768 return log1pf (x);
769 #else
694 float retval; 770 float retval;
695 771
696 float ax = fabs (x); 772 float ax = fabs (x);
697 773
698 if (ax < 0.2) 774 if (ax < 0.2)
706 } 782 }
707 else 783 else
708 retval = gnulib::logf (1.0f + x); 784 retval = gnulib::logf (1.0f + x);
709 785
710 return retval; 786 return retval;
711 }
712 #endif 787 #endif
788 }
713 789
714 FloatComplex 790 FloatComplex
715 log1p (const FloatComplex& x) 791 xlog1p (const FloatComplex& x)
716 { 792 {
717 FloatComplex retval; 793 FloatComplex retval;
718 794
719 float r = x.real (), i = x.imag (); 795 float r = x.real (), i = x.imag ();
720 796
721 if (fabs (r) < 0.5 && fabs (i) < 0.5) 797 if (fabs (r) < 0.5 && fabs (i) < 0.5)
722 { 798 {
723 float u = 2*r + r*r + i*i; 799 float u = 2*r + r*r + i*i;
724 retval = FloatComplex (log1p (u / (1+sqrt (u+1))), 800 retval = FloatComplex (xlog1p (u / (1+sqrt (u+1))),
725 atan2 (1 + r, i)); 801 atan2 (1 + r, i));
726 } 802 }
727 else 803 else
728 retval = std::log (FloatComplex (1) + x); 804 retval = std::log (FloatComplex (1) + x);
729 805
730 return retval; 806 return retval;
731 } 807 }
732 808
733 #if ! defined (HAVE_CBRTF) 809 float
734 float cbrtf (float x) 810 xcbrt (float x)
735 { 811 {
736 static const float one_third = 0.3333333333333333333f; 812 #if defined (HAVE_CBRTF)
737 if (xfinite (x)) 813 return cbrtf (x);
738 { 814 #else
739 // Use pow. 815 return xxcbrt (x);
740 float y = std::pow (std::abs (x), one_third) * signum (x);
741 // Correct for better accuracy.
742 return (x / (y*y) + y + y) / 3;
743 }
744 else
745 return x;
746 }
747 #endif 816 #endif
817 }
748 818
749 static inline Complex 819 static inline Complex
750 zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr); 820 zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
751 821
752 static inline Complex 822 static inline Complex
2890 Complex rc_log1p (double x) 2960 Complex rc_log1p (double x)
2891 { 2961 {
2892 const double pi = 3.14159265358979323846; 2962 const double pi = 3.14159265358979323846;
2893 return (x < -1.0 2963 return (x < -1.0
2894 ? Complex (gnulib::log (-(1.0 + x)), pi) 2964 ? Complex (gnulib::log (-(1.0 + x)), pi)
2895 : Complex (log1p (x))); 2965 : Complex (xlog1p (x)));
2896 } 2966 }
2897 2967
2898 FloatComplex rc_log1p (float x) 2968 FloatComplex rc_log1p (float x)
2899 { 2969 {
2900 const float pi = 3.14159265358979323846f; 2970 const float pi = 3.14159265358979323846f;
2901 return (x < -1.0f 2971 return (x < -1.0f
2902 ? FloatComplex (gnulib::logf (-(1.0f + x)), pi) 2972 ? FloatComplex (gnulib::logf (-(1.0f + x)), pi)
2903 : FloatComplex (log1pf (x))); 2973 : FloatComplex (xlog1p (x)));
2904 } 2974 }
2905 2975
2906 // This algorithm is due to P. J. Acklam. 2976 // This algorithm is due to P. J. Acklam.
2907 // 2977 //
2908 // See http://home.online.no/~pjacklam/notes/invnorm/ 2978 // See http://home.online.no/~pjacklam/notes/invnorm/