1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 1996-2021 The Octave Project Developers 4 // 5 // See the file COPYRIGHT.md in the top-level directory of this 6 // distribution or <https://octave.org/copyright/>. 7 // 8 // This file is part of Octave. 9 // 10 // Octave is free software: you can redistribute it and/or modify it 11 // under the terms of the GNU General Public License as published by 12 // the Free Software Foundation, either version 3 of the License, or 13 // (at your option) any later version. 14 // 15 // Octave is distributed in the hope that it will be useful, but 16 // WITHOUT ANY WARRANTY; without even the implied warranty of 17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 // GNU General Public License for more details. 19 // 20 // You should have received a copy of the GNU General Public License 21 // along with Octave; see the file COPYING. If not, see 22 // <https://www.gnu.org/licenses/>. 23 // 24 //////////////////////////////////////////////////////////////////////// 25 26 #if defined (HAVE_CONFIG_H) 27 # include "config.h" 28 #endif 29 30 #include <cmath> 31 32 #include <algorithm> 33 #include <limits> 34 #include <string> 35 36 #include "CColVector.h" 37 #include "CMatrix.h" 38 #include "CNDArray.h" 39 #include "Faddeeva.hh" 40 #include "dMatrix.h" 41 #include "dNDArray.h" 42 #include "dRowVector.h" 43 #include "f77-fcn.h" 44 #include "fCColVector.h" 45 #include "fCMatrix.h" 46 #include "fCNDArray.h" 47 #include "fMatrix.h" 48 #include "fNDArray.h" 49 #include "fRowVector.h" 50 #include "lo-amos-proto.h" 51 #include "lo-error.h" 52 #include "lo-ieee.h" 53 #include "lo-mappers.h" 54 #include "lo-slatec-proto.h" 55 #include "lo-specfun.h" 56 #include "mx-inlines.cc" 57 58 namespace octave 59 { 60 namespace math 61 { 62 static inline Complex bessel_return_value(const Complex & val,octave_idx_type ierr)63 bessel_return_value (const Complex& val, octave_idx_type ierr) 64 { 65 static const Complex inf_val 66 = Complex (numeric_limits<double>::Inf (), 67 numeric_limits<double>::Inf ()); 68 69 static const Complex nan_val 70 = Complex (numeric_limits<double>::NaN (), 71 numeric_limits<double>::NaN ()); 72 73 Complex retval; 74 75 switch (ierr) 76 { 77 case 0: 78 case 3: 79 case 4: 80 retval = val; 81 break; 82 83 case 2: 84 retval = inf_val; 85 break; 86 87 default: 88 retval = nan_val; 89 break; 90 } 91 92 return retval; 93 } 94 95 static inline FloatComplex bessel_return_value(const FloatComplex & val,octave_idx_type ierr)96 bessel_return_value (const FloatComplex& val, octave_idx_type ierr) 97 { 98 static const FloatComplex inf_val 99 = FloatComplex (numeric_limits<float>::Inf (), 100 numeric_limits<float>::Inf ()); 101 102 static const FloatComplex nan_val 103 = FloatComplex (numeric_limits<float>::NaN (), 104 numeric_limits<float>::NaN ()); 105 106 FloatComplex retval; 107 108 switch (ierr) 109 { 110 case 0: 111 case 3: 112 case 4: 113 retval = val; 114 break; 115 116 case 2: 117 retval = inf_val; 118 break; 119 120 default: 121 retval = nan_val; 122 break; 123 } 124 125 return retval; 126 } 127 128 129 130 Complex airy(const Complex & z,bool deriv,bool scaled,octave_idx_type & ierr)131 airy (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr) 132 { 133 double ar = 0.0; 134 double ai = 0.0; 135 136 double zr = z.real (); 137 double zi = z.imag (); 138 139 F77_INT id = (deriv ? 1 : 0); 140 F77_INT nz, t_ierr; 141 142 F77_FUNC (zairy, ZAIRY) (zr, zi, id, 2, ar, ai, nz, t_ierr); 143 144 ierr = t_ierr; 145 146 if (! scaled) 147 { 148 Complex expz = exp (- 2.0 / 3.0 * z * sqrt (z)); 149 150 double rexpz = expz.real (); 151 double iexpz = expz.imag (); 152 153 double tmp = ar*rexpz - ai*iexpz; 154 155 ai = ar*iexpz + ai*rexpz; 156 ar = tmp; 157 } 158 159 if (zi == 0.0 && (! scaled || zr >= 0.0)) 160 ai = 0.0; 161 162 return bessel_return_value (Complex (ar, ai), ierr); 163 } 164 165 ComplexMatrix airy(const ComplexMatrix & z,bool deriv,bool scaled,Array<octave_idx_type> & ierr)166 airy (const ComplexMatrix& z, bool deriv, bool scaled, 167 Array<octave_idx_type>& ierr) 168 { 169 octave_idx_type nr = z.rows (); 170 octave_idx_type nc = z.cols (); 171 172 ComplexMatrix retval (nr, nc); 173 174 ierr.resize (dim_vector (nr, nc)); 175 176 for (octave_idx_type j = 0; j < nc; j++) 177 for (octave_idx_type i = 0; i < nr; i++) 178 retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j)); 179 180 return retval; 181 } 182 183 ComplexNDArray airy(const ComplexNDArray & z,bool deriv,bool scaled,Array<octave_idx_type> & ierr)184 airy (const ComplexNDArray& z, bool deriv, bool scaled, 185 Array<octave_idx_type>& ierr) 186 { 187 dim_vector dv = z.dims (); 188 octave_idx_type nel = dv.numel (); 189 ComplexNDArray retval (dv); 190 191 ierr.resize (dv); 192 193 for (octave_idx_type i = 0; i < nel; i++) 194 retval(i) = airy (z(i), deriv, scaled, ierr(i)); 195 196 return retval; 197 } 198 199 FloatComplex airy(const FloatComplex & z,bool deriv,bool scaled,octave_idx_type & ierr)200 airy (const FloatComplex& z, bool deriv, bool scaled, 201 octave_idx_type& ierr) 202 { 203 FloatComplex a; 204 205 F77_INT id = (deriv ? 1 : 0); 206 F77_INT nz, t_ierr; 207 208 F77_FUNC (cairy, CAIRY) (F77_CONST_CMPLX_ARG (&z), id, 2, 209 F77_CMPLX_ARG (&a), nz, t_ierr); 210 211 ierr = t_ierr; 212 213 float ar = a.real (); 214 float ai = a.imag (); 215 216 if (! scaled) 217 { 218 FloatComplex expz = exp (- 2.0f / 3.0f * z * sqrt (z)); 219 220 float rexpz = expz.real (); 221 float iexpz = expz.imag (); 222 223 float tmp = ar*rexpz - ai*iexpz; 224 225 ai = ar*iexpz + ai*rexpz; 226 ar = tmp; 227 } 228 229 if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0)) 230 ai = 0.0; 231 232 return bessel_return_value (FloatComplex (ar, ai), ierr); 233 } 234 235 FloatComplexMatrix airy(const FloatComplexMatrix & z,bool deriv,bool scaled,Array<octave_idx_type> & ierr)236 airy (const FloatComplexMatrix& z, bool deriv, bool scaled, 237 Array<octave_idx_type>& ierr) 238 { 239 octave_idx_type nr = z.rows (); 240 octave_idx_type nc = z.cols (); 241 242 FloatComplexMatrix retval (nr, nc); 243 244 ierr.resize (dim_vector (nr, nc)); 245 246 for (octave_idx_type j = 0; j < nc; j++) 247 for (octave_idx_type i = 0; i < nr; i++) 248 retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j)); 249 250 return retval; 251 } 252 253 FloatComplexNDArray airy(const FloatComplexNDArray & z,bool deriv,bool scaled,Array<octave_idx_type> & ierr)254 airy (const FloatComplexNDArray& z, bool deriv, bool scaled, 255 Array<octave_idx_type>& ierr) 256 { 257 dim_vector dv = z.dims (); 258 octave_idx_type nel = dv.numel (); 259 FloatComplexNDArray retval (dv); 260 261 ierr.resize (dv); 262 263 for (octave_idx_type i = 0; i < nel; i++) 264 retval(i) = airy (z(i), deriv, scaled, ierr(i)); 265 266 return retval; 267 } 268 269 static inline bool is_integer_value(double x)270 is_integer_value (double x) 271 { 272 return x == static_cast<double> (static_cast<long> (x)); 273 } 274 275 static inline Complex 276 zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr); 277 278 static inline Complex 279 zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr); 280 281 static inline Complex 282 zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr); 283 284 static inline Complex 285 zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr); 286 287 static inline Complex 288 zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr); 289 290 static inline Complex 291 zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr); 292 293 static inline Complex zbesj(const Complex & z,double alpha,int kode,octave_idx_type & ierr)294 zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr) 295 { 296 Complex retval; 297 298 if (alpha >= 0.0) 299 { 300 double yr = 0.0; 301 double yi = 0.0; 302 303 F77_INT nz, t_ierr; 304 305 double zr = z.real (); 306 double zi = z.imag (); 307 308 F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr); 309 310 ierr = t_ierr; 311 312 if (zi == 0.0 && zr >= 0.0) 313 yi = 0.0; 314 315 retval = bessel_return_value (Complex (yr, yi), ierr); 316 } 317 else if (is_integer_value (alpha)) 318 { 319 // zbesy can overflow as z->0, and cause troubles for generic case below 320 alpha = -alpha; 321 Complex tmp = zbesj (z, alpha, kode, ierr); 322 if ((static_cast<long> (alpha)) & 1) 323 tmp = - tmp; 324 retval = bessel_return_value (tmp, ierr); 325 } 326 else 327 { 328 alpha = -alpha; 329 330 Complex tmp = cos (M_PI * alpha) * zbesj (z, alpha, kode, ierr); 331 332 if (ierr == 0 || ierr == 3) 333 { 334 tmp -= sin (M_PI * alpha) * zbesy (z, alpha, kode, ierr); 335 336 retval = bessel_return_value (tmp, ierr); 337 } 338 else 339 retval = Complex (numeric_limits<double>::NaN (), 340 numeric_limits<double>::NaN ()); 341 } 342 343 return retval; 344 } 345 346 static inline Complex zbesy(const Complex & z,double alpha,int kode,octave_idx_type & ierr)347 zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr) 348 { 349 Complex retval; 350 351 if (alpha >= 0.0) 352 { 353 double yr = 0.0; 354 double yi = 0.0; 355 356 F77_INT nz, t_ierr; 357 358 double wr, wi; 359 360 double zr = z.real (); 361 double zi = z.imag (); 362 363 ierr = 0; 364 365 if (zr == 0.0 && zi == 0.0) 366 { 367 yr = -numeric_limits<double>::Inf (); 368 yi = 0.0; 369 } 370 else 371 { 372 F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, kode, 1, &yr, &yi, nz, 373 &wr, &wi, t_ierr); 374 375 ierr = t_ierr; 376 377 if (zi == 0.0 && zr >= 0.0) 378 yi = 0.0; 379 } 380 381 return bessel_return_value (Complex (yr, yi), ierr); 382 } 383 else if (is_integer_value (alpha - 0.5)) 384 { 385 // zbesy can overflow as z->0, and cause troubles for generic case below 386 alpha = -alpha; 387 Complex tmp = zbesj (z, alpha, kode, ierr); 388 if ((static_cast<long> (alpha - 0.5)) & 1) 389 tmp = - tmp; 390 retval = bessel_return_value (tmp, ierr); 391 } 392 else 393 { 394 alpha = -alpha; 395 396 Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr); 397 398 if (ierr == 0 || ierr == 3) 399 { 400 tmp += sin (M_PI * alpha) * zbesj (z, alpha, kode, ierr); 401 402 retval = bessel_return_value (tmp, ierr); 403 } 404 else 405 retval = Complex (numeric_limits<double>::NaN (), 406 numeric_limits<double>::NaN ()); 407 } 408 409 return retval; 410 } 411 412 static inline Complex zbesi(const Complex & z,double alpha,int kode,octave_idx_type & ierr)413 zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr) 414 { 415 Complex retval; 416 417 if (alpha >= 0.0) 418 { 419 double yr = 0.0; 420 double yi = 0.0; 421 422 F77_INT nz, t_ierr; 423 424 double zr = z.real (); 425 double zi = z.imag (); 426 427 F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr); 428 429 ierr = t_ierr; 430 431 if (zi == 0.0 && zr >= 0.0) 432 yi = 0.0; 433 434 retval = bessel_return_value (Complex (yr, yi), ierr); 435 } 436 else if (is_integer_value (alpha)) 437 { 438 // zbesi can overflow as z->0, and cause troubles for generic case below 439 alpha = -alpha; 440 Complex tmp = zbesi (z, alpha, kode, ierr); 441 retval = bessel_return_value (tmp, ierr); 442 } 443 else 444 { 445 alpha = -alpha; 446 447 Complex tmp = zbesi (z, alpha, kode, ierr); 448 449 if (ierr == 0 || ierr == 3) 450 { 451 Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha) 452 * zbesk (z, alpha, kode, ierr); 453 454 if (kode == 2) 455 { 456 // Compensate for different scaling factor of besk. 457 tmp2 *= exp (-z - std::abs (z.real ())); 458 } 459 460 tmp += tmp2; 461 462 retval = bessel_return_value (tmp, ierr); 463 } 464 else 465 retval = Complex (numeric_limits<double>::NaN (), 466 numeric_limits<double>::NaN ()); 467 } 468 469 return retval; 470 } 471 472 static inline Complex zbesk(const Complex & z,double alpha,int kode,octave_idx_type & ierr)473 zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr) 474 { 475 Complex retval; 476 477 if (alpha >= 0.0) 478 { 479 double yr = 0.0; 480 double yi = 0.0; 481 482 F77_INT nz, t_ierr; 483 484 double zr = z.real (); 485 double zi = z.imag (); 486 487 ierr = 0; 488 489 if (zr == 0.0 && zi == 0.0) 490 { 491 yr = numeric_limits<double>::Inf (); 492 yi = 0.0; 493 } 494 else 495 { 496 F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, kode, 1, &yr, &yi, nz, 497 t_ierr); 498 499 ierr = t_ierr; 500 501 if (zi == 0.0 && zr >= 0.0) 502 yi = 0.0; 503 } 504 505 retval = bessel_return_value (Complex (yr, yi), ierr); 506 } 507 else 508 { 509 Complex tmp = zbesk (z, -alpha, kode, ierr); 510 511 retval = bessel_return_value (tmp, ierr); 512 } 513 514 return retval; 515 } 516 517 static inline Complex zbesh1(const Complex & z,double alpha,int kode,octave_idx_type & ierr)518 zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr) 519 { 520 Complex retval; 521 522 if (alpha >= 0.0) 523 { 524 double yr = 0.0; 525 double yi = 0.0; 526 527 F77_INT nz, t_ierr; 528 529 double zr = z.real (); 530 double zi = z.imag (); 531 532 F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 1, 1, &yr, &yi, nz, 533 t_ierr); 534 535 ierr = t_ierr; 536 537 retval = bessel_return_value (Complex (yr, yi), ierr); 538 } 539 else 540 { 541 alpha = -alpha; 542 543 static const Complex eye = Complex (0.0, 1.0); 544 545 Complex tmp = exp (M_PI * alpha * eye) * zbesh1 (z, alpha, kode, ierr); 546 547 retval = bessel_return_value (tmp, ierr); 548 } 549 550 return retval; 551 } 552 553 static inline Complex zbesh2(const Complex & z,double alpha,int kode,octave_idx_type & ierr)554 zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr) 555 { 556 Complex retval; 557 558 if (alpha >= 0.0) 559 { 560 double yr = 0.0; 561 double yi = 0.0; 562 563 F77_INT nz, t_ierr; 564 565 double zr = z.real (); 566 double zi = z.imag (); 567 568 F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 2, 1, &yr, &yi, nz, 569 t_ierr); 570 571 ierr = t_ierr; 572 573 retval = bessel_return_value (Complex (yr, yi), ierr); 574 } 575 else 576 { 577 alpha = -alpha; 578 579 static const Complex eye = Complex (0.0, 1.0); 580 581 Complex tmp = exp (-M_PI * alpha * eye) * zbesh2 (z, alpha, kode, ierr); 582 583 retval = bessel_return_value (tmp, ierr); 584 } 585 586 return retval; 587 } 588 589 typedef Complex (*dptr) (const Complex&, double, int, octave_idx_type&); 590 591 static inline Complex do_bessel(dptr f,const char *,double alpha,const Complex & x,bool scaled,octave_idx_type & ierr)592 do_bessel (dptr f, const char *, double alpha, const Complex& x, 593 bool scaled, octave_idx_type& ierr) 594 { 595 Complex retval; 596 597 retval = f (x, alpha, (scaled ? 2 : 1), ierr); 598 599 return retval; 600 } 601 602 static inline ComplexMatrix do_bessel(dptr f,const char *,double alpha,const ComplexMatrix & x,bool scaled,Array<octave_idx_type> & ierr)603 do_bessel (dptr f, const char *, double alpha, const ComplexMatrix& x, 604 bool scaled, Array<octave_idx_type>& ierr) 605 { 606 octave_idx_type nr = x.rows (); 607 octave_idx_type nc = x.cols (); 608 609 ComplexMatrix retval (nr, nc); 610 611 ierr.resize (dim_vector (nr, nc)); 612 613 for (octave_idx_type j = 0; j < nc; j++) 614 for (octave_idx_type i = 0; i < nr; i++) 615 retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j)); 616 617 return retval; 618 } 619 620 static inline ComplexMatrix do_bessel(dptr f,const char *,const Matrix & alpha,const Complex & x,bool scaled,Array<octave_idx_type> & ierr)621 do_bessel (dptr f, const char *, const Matrix& alpha, const Complex& x, 622 bool scaled, Array<octave_idx_type>& ierr) 623 { 624 octave_idx_type nr = alpha.rows (); 625 octave_idx_type nc = alpha.cols (); 626 627 ComplexMatrix retval (nr, nc); 628 629 ierr.resize (dim_vector (nr, nc)); 630 631 for (octave_idx_type j = 0; j < nc; j++) 632 for (octave_idx_type i = 0; i < nr; i++) 633 retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j)); 634 635 return retval; 636 } 637 638 static inline ComplexMatrix do_bessel(dptr f,const char * fn,const Matrix & alpha,const ComplexMatrix & x,bool scaled,Array<octave_idx_type> & ierr)639 do_bessel (dptr f, const char *fn, const Matrix& alpha, 640 const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) 641 { 642 ComplexMatrix retval; 643 644 octave_idx_type x_nr = x.rows (); 645 octave_idx_type x_nc = x.cols (); 646 647 octave_idx_type alpha_nr = alpha.rows (); 648 octave_idx_type alpha_nc = alpha.cols (); 649 650 if (x_nr != alpha_nr || x_nc != alpha_nc) 651 (*current_liboctave_error_handler) 652 ("%s: the sizes of alpha and x must conform", fn); 653 654 octave_idx_type nr = x_nr; 655 octave_idx_type nc = x_nc; 656 657 retval.resize (nr, nc); 658 659 ierr.resize (dim_vector (nr, nc)); 660 661 for (octave_idx_type j = 0; j < nc; j++) 662 for (octave_idx_type i = 0; i < nr; i++) 663 retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j)); 664 665 return retval; 666 } 667 668 static inline ComplexNDArray do_bessel(dptr f,const char *,double alpha,const ComplexNDArray & x,bool scaled,Array<octave_idx_type> & ierr)669 do_bessel (dptr f, const char *, double alpha, const ComplexNDArray& x, 670 bool scaled, Array<octave_idx_type>& ierr) 671 { 672 dim_vector dv = x.dims (); 673 octave_idx_type nel = dv.numel (); 674 ComplexNDArray retval (dv); 675 676 ierr.resize (dv); 677 678 for (octave_idx_type i = 0; i < nel; i++) 679 retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i)); 680 681 return retval; 682 } 683 684 static inline ComplexNDArray do_bessel(dptr f,const char *,const NDArray & alpha,const Complex & x,bool scaled,Array<octave_idx_type> & ierr)685 do_bessel (dptr f, const char *, const NDArray& alpha, const Complex& x, 686 bool scaled, Array<octave_idx_type>& ierr) 687 { 688 dim_vector dv = alpha.dims (); 689 octave_idx_type nel = dv.numel (); 690 ComplexNDArray retval (dv); 691 692 ierr.resize (dv); 693 694 for (octave_idx_type i = 0; i < nel; i++) 695 retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i)); 696 697 return retval; 698 } 699 700 static inline ComplexNDArray do_bessel(dptr f,const char * fn,const NDArray & alpha,const ComplexNDArray & x,bool scaled,Array<octave_idx_type> & ierr)701 do_bessel (dptr f, const char *fn, const NDArray& alpha, 702 const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) 703 { 704 dim_vector dv = x.dims (); 705 ComplexNDArray retval; 706 707 if (dv != alpha.dims ()) 708 (*current_liboctave_error_handler) 709 ("%s: the sizes of alpha and x must conform", fn); 710 711 octave_idx_type nel = dv.numel (); 712 713 retval.resize (dv); 714 ierr.resize (dv); 715 716 for (octave_idx_type i = 0; i < nel; i++) 717 retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i)); 718 719 return retval; 720 } 721 722 static inline ComplexMatrix do_bessel(dptr f,const char *,const RowVector & alpha,const ComplexColumnVector & x,bool scaled,Array<octave_idx_type> & ierr)723 do_bessel (dptr f, const char *, const RowVector& alpha, 724 const ComplexColumnVector& x, bool scaled, 725 Array<octave_idx_type>& ierr) 726 { 727 octave_idx_type nr = x.numel (); 728 octave_idx_type nc = alpha.numel (); 729 730 ComplexMatrix retval (nr, nc); 731 732 ierr.resize (dim_vector (nr, nc)); 733 734 for (octave_idx_type j = 0; j < nc; j++) 735 for (octave_idx_type i = 0; i < nr; i++) 736 retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j)); 737 738 return retval; 739 } 740 741 #define SS_BESSEL(name, fcn) \ 742 Complex \ 743 name (double alpha, const Complex& x, bool scaled, octave_idx_type& ierr) \ 744 { \ 745 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 746 } 747 748 #define SM_BESSEL(name, fcn) \ 749 ComplexMatrix \ 750 name (double alpha, const ComplexMatrix& x, bool scaled, \ 751 Array<octave_idx_type>& ierr) \ 752 { \ 753 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 754 } 755 756 #define MS_BESSEL(name, fcn) \ 757 ComplexMatrix \ 758 name (const Matrix& alpha, const Complex& x, bool scaled, \ 759 Array<octave_idx_type>& ierr) \ 760 { \ 761 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 762 } 763 764 #define MM_BESSEL(name, fcn) \ 765 ComplexMatrix \ 766 name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \ 767 Array<octave_idx_type>& ierr) \ 768 { \ 769 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 770 } 771 772 #define SN_BESSEL(name, fcn) \ 773 ComplexNDArray \ 774 name (double alpha, const ComplexNDArray& x, bool scaled, \ 775 Array<octave_idx_type>& ierr) \ 776 { \ 777 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 778 } 779 780 #define NS_BESSEL(name, fcn) \ 781 ComplexNDArray \ 782 name (const NDArray& alpha, const Complex& x, bool scaled, \ 783 Array<octave_idx_type>& ierr) \ 784 { \ 785 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 786 } 787 788 #define NN_BESSEL(name, fcn) \ 789 ComplexNDArray \ 790 name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \ 791 Array<octave_idx_type>& ierr) \ 792 { \ 793 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 794 } 795 796 #define RC_BESSEL(name, fcn) \ 797 ComplexMatrix \ 798 name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \ 799 Array<octave_idx_type>& ierr) \ 800 { \ 801 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 802 } 803 804 #define ALL_BESSEL(name, fcn) \ 805 SS_BESSEL (name, fcn) \ 806 SM_BESSEL (name, fcn) \ 807 MS_BESSEL (name, fcn) \ 808 MM_BESSEL (name, fcn) \ 809 SN_BESSEL (name, fcn) \ 810 NS_BESSEL (name, fcn) \ 811 NN_BESSEL (name, fcn) \ 812 RC_BESSEL (name, fcn) 813 814 ALL_BESSEL (besselj, zbesj) 815 ALL_BESSEL (bessely, zbesy) 816 ALL_BESSEL (besseli, zbesi) 817 ALL_BESSEL (besselk, zbesk) 818 ALL_BESSEL (besselh1, zbesh1) 819 ALL_BESSEL (besselh2, zbesh2) 820 821 #undef ALL_BESSEL 822 #undef SS_BESSEL 823 #undef SM_BESSEL 824 #undef MS_BESSEL 825 #undef MM_BESSEL 826 #undef SN_BESSEL 827 #undef NS_BESSEL 828 #undef NN_BESSEL 829 #undef RC_BESSEL 830 831 static inline FloatComplex 832 cbesj (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); 833 834 static inline FloatComplex 835 cbesy (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); 836 837 static inline FloatComplex 838 cbesi (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); 839 840 static inline FloatComplex 841 cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); 842 843 static inline FloatComplex 844 cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); 845 846 static inline FloatComplex 847 cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr); 848 849 static inline bool is_integer_value(float x)850 is_integer_value (float x) 851 { 852 return x == static_cast<float> (static_cast<long> (x)); 853 } 854 855 static inline FloatComplex cbesj(const FloatComplex & z,float alpha,int kode,octave_idx_type & ierr)856 cbesj (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) 857 { 858 FloatComplex retval; 859 860 if (alpha >= 0.0) 861 { 862 FloatComplex y = 0.0; 863 864 F77_INT nz, t_ierr; 865 866 F77_FUNC (cbesj, CBESJ) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 867 F77_CMPLX_ARG (&y), nz, t_ierr); 868 869 ierr = t_ierr; 870 871 if (z.imag () == 0.0 && z.real () >= 0.0) 872 y = FloatComplex (y.real (), 0.0); 873 874 retval = bessel_return_value (y, ierr); 875 } 876 else if (is_integer_value (alpha)) 877 { 878 // zbesy can overflow as z->0, and cause troubles for generic case below 879 alpha = -alpha; 880 FloatComplex tmp = cbesj (z, alpha, kode, ierr); 881 if ((static_cast<long> (alpha)) & 1) 882 tmp = - tmp; 883 retval = bessel_return_value (tmp, ierr); 884 } 885 else 886 { 887 alpha = -alpha; 888 889 FloatComplex tmp = cosf (static_cast<float> (M_PI) * alpha) 890 * cbesj (z, alpha, kode, ierr); 891 892 if (ierr == 0 || ierr == 3) 893 { 894 tmp -= sinf (static_cast<float> (M_PI) * alpha) 895 * cbesy (z, alpha, kode, ierr); 896 897 retval = bessel_return_value (tmp, ierr); 898 } 899 else 900 retval = FloatComplex (numeric_limits<float>::NaN (), 901 numeric_limits<float>::NaN ()); 902 } 903 904 return retval; 905 } 906 907 static inline FloatComplex cbesy(const FloatComplex & z,float alpha,int kode,octave_idx_type & ierr)908 cbesy (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) 909 { 910 FloatComplex retval; 911 912 if (alpha >= 0.0) 913 { 914 FloatComplex y = 0.0; 915 916 F77_INT nz, t_ierr; 917 918 FloatComplex w; 919 920 ierr = 0; 921 922 if (z.real () == 0.0 && z.imag () == 0.0) 923 { 924 y = FloatComplex (-numeric_limits<float>::Inf (), 0.0); 925 } 926 else 927 { 928 F77_FUNC (cbesy, CBESY) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 929 F77_CMPLX_ARG (&y), nz, 930 F77_CMPLX_ARG (&w), t_ierr); 931 932 ierr = t_ierr; 933 934 if (z.imag () == 0.0 && z.real () >= 0.0) 935 y = FloatComplex (y.real (), 0.0); 936 } 937 938 return bessel_return_value (y, ierr); 939 } 940 else if (is_integer_value (alpha - 0.5)) 941 { 942 // zbesy can overflow as z->0, and cause troubles for generic case below 943 alpha = -alpha; 944 FloatComplex tmp = cbesj (z, alpha, kode, ierr); 945 if ((static_cast<long> (alpha - 0.5)) & 1) 946 tmp = - tmp; 947 retval = bessel_return_value (tmp, ierr); 948 } 949 else 950 { 951 alpha = -alpha; 952 953 FloatComplex tmp = cosf (static_cast<float> (M_PI) * alpha) 954 * cbesy (z, alpha, kode, ierr); 955 956 if (ierr == 0 || ierr == 3) 957 { 958 tmp += sinf (static_cast<float> (M_PI) * alpha) 959 * cbesj (z, alpha, kode, ierr); 960 961 retval = bessel_return_value (tmp, ierr); 962 } 963 else 964 retval = FloatComplex (numeric_limits<float>::NaN (), 965 numeric_limits<float>::NaN ()); 966 } 967 968 return retval; 969 } 970 971 static inline FloatComplex cbesi(const FloatComplex & z,float alpha,int kode,octave_idx_type & ierr)972 cbesi (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) 973 { 974 FloatComplex retval; 975 976 if (alpha >= 0.0) 977 { 978 FloatComplex y = 0.0; 979 980 F77_INT nz, t_ierr; 981 982 F77_FUNC (cbesi, CBESI) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 983 F77_CMPLX_ARG (&y), nz, t_ierr); 984 985 ierr = t_ierr; 986 987 if (z.imag () == 0.0 && z.real () >= 0.0) 988 y = FloatComplex (y.real (), 0.0); 989 990 retval = bessel_return_value (y, ierr); 991 } 992 else 993 { 994 alpha = -alpha; 995 996 FloatComplex tmp = cbesi (z, alpha, kode, ierr); 997 998 if (ierr == 0 || ierr == 3) 999 { 1000 FloatComplex tmp2 = static_cast<float> (2.0 / M_PI) 1001 * sinf (static_cast<float> (M_PI) * alpha) 1002 * cbesk (z, alpha, kode, ierr); 1003 1004 if (kode == 2) 1005 { 1006 // Compensate for different scaling factor of besk. 1007 tmp2 *= exp (-z - std::abs (z.real ())); 1008 } 1009 1010 tmp += tmp2; 1011 1012 retval = bessel_return_value (tmp, ierr); 1013 } 1014 else 1015 retval = FloatComplex (numeric_limits<float>::NaN (), 1016 numeric_limits<float>::NaN ()); 1017 } 1018 1019 return retval; 1020 } 1021 1022 static inline FloatComplex cbesk(const FloatComplex & z,float alpha,int kode,octave_idx_type & ierr)1023 cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) 1024 { 1025 FloatComplex retval; 1026 1027 if (alpha >= 0.0) 1028 { 1029 FloatComplex y = 0.0; 1030 1031 F77_INT nz, t_ierr; 1032 1033 ierr = 0; 1034 1035 if (z.real () == 0.0 && z.imag () == 0.0) 1036 { 1037 y = FloatComplex (numeric_limits<float>::Inf (), 0.0); 1038 } 1039 else 1040 { 1041 F77_FUNC (cbesk, CBESK) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 1042 F77_CMPLX_ARG (&y), nz, t_ierr); 1043 1044 ierr = t_ierr; 1045 1046 if (z.imag () == 0.0 && z.real () >= 0.0) 1047 y = FloatComplex (y.real (), 0.0); 1048 } 1049 1050 retval = bessel_return_value (y, ierr); 1051 } 1052 else 1053 { 1054 FloatComplex tmp = cbesk (z, -alpha, kode, ierr); 1055 1056 retval = bessel_return_value (tmp, ierr); 1057 } 1058 1059 return retval; 1060 } 1061 1062 static inline FloatComplex cbesh1(const FloatComplex & z,float alpha,int kode,octave_idx_type & ierr)1063 cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) 1064 { 1065 FloatComplex retval; 1066 1067 if (alpha >= 0.0) 1068 { 1069 FloatComplex y = 0.0; 1070 1071 F77_INT nz, t_ierr; 1072 1073 F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 1, 1074 F77_CMPLX_ARG (&y), nz, t_ierr); 1075 1076 ierr = t_ierr; 1077 1078 retval = bessel_return_value (y, ierr); 1079 } 1080 else 1081 { 1082 alpha = -alpha; 1083 1084 static const FloatComplex eye = FloatComplex (0.0, 1.0); 1085 1086 FloatComplex tmp = exp (static_cast<float> (M_PI) * alpha * eye) 1087 * cbesh1 (z, alpha, kode, ierr); 1088 1089 retval = bessel_return_value (tmp, ierr); 1090 } 1091 1092 return retval; 1093 } 1094 1095 static inline FloatComplex cbesh2(const FloatComplex & z,float alpha,int kode,octave_idx_type & ierr)1096 cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr) 1097 { 1098 FloatComplex retval; 1099 1100 if (alpha >= 0.0) 1101 { 1102 FloatComplex y = 0.0;; 1103 1104 F77_INT nz, t_ierr; 1105 1106 F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 2, 1, 1107 F77_CMPLX_ARG (&y), nz, t_ierr); 1108 1109 ierr = t_ierr; 1110 1111 retval = bessel_return_value (y, ierr); 1112 } 1113 else 1114 { 1115 alpha = -alpha; 1116 1117 static const FloatComplex eye = FloatComplex (0.0, 1.0); 1118 1119 FloatComplex tmp = exp (-static_cast<float> (M_PI) * alpha * eye) 1120 * cbesh2 (z, alpha, kode, ierr); 1121 1122 retval = bessel_return_value (tmp, ierr); 1123 } 1124 1125 return retval; 1126 } 1127 1128 typedef FloatComplex (*fptr) (const FloatComplex&, float, int, 1129 octave_idx_type&); 1130 1131 static inline FloatComplex do_bessel(fptr f,const char *,float alpha,const FloatComplex & x,bool scaled,octave_idx_type & ierr)1132 do_bessel (fptr f, const char *, float alpha, const FloatComplex& x, 1133 bool scaled, octave_idx_type& ierr) 1134 { 1135 FloatComplex retval; 1136 1137 retval = f (x, alpha, (scaled ? 2 : 1), ierr); 1138 1139 return retval; 1140 } 1141 1142 static inline FloatComplexMatrix do_bessel(fptr f,const char *,float alpha,const FloatComplexMatrix & x,bool scaled,Array<octave_idx_type> & ierr)1143 do_bessel (fptr f, const char *, float alpha, const FloatComplexMatrix& x, 1144 bool scaled, Array<octave_idx_type>& ierr) 1145 { 1146 octave_idx_type nr = x.rows (); 1147 octave_idx_type nc = x.cols (); 1148 1149 FloatComplexMatrix retval (nr, nc); 1150 1151 ierr.resize (dim_vector (nr, nc)); 1152 1153 for (octave_idx_type j = 0; j < nc; j++) 1154 for (octave_idx_type i = 0; i < nr; i++) 1155 retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j)); 1156 1157 return retval; 1158 } 1159 1160 static inline FloatComplexMatrix do_bessel(fptr f,const char *,const FloatMatrix & alpha,const FloatComplex & x,bool scaled,Array<octave_idx_type> & ierr)1161 do_bessel (fptr f, const char *, const FloatMatrix& alpha, 1162 const FloatComplex& x, 1163 bool scaled, Array<octave_idx_type>& ierr) 1164 { 1165 octave_idx_type nr = alpha.rows (); 1166 octave_idx_type nc = alpha.cols (); 1167 1168 FloatComplexMatrix retval (nr, nc); 1169 1170 ierr.resize (dim_vector (nr, nc)); 1171 1172 for (octave_idx_type j = 0; j < nc; j++) 1173 for (octave_idx_type i = 0; i < nr; i++) 1174 retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j)); 1175 1176 return retval; 1177 } 1178 1179 static inline FloatComplexMatrix do_bessel(fptr f,const char * fn,const FloatMatrix & alpha,const FloatComplexMatrix & x,bool scaled,Array<octave_idx_type> & ierr)1180 do_bessel (fptr f, const char *fn, const FloatMatrix& alpha, 1181 const FloatComplexMatrix& x, bool scaled, 1182 Array<octave_idx_type>& ierr) 1183 { 1184 FloatComplexMatrix retval; 1185 1186 octave_idx_type x_nr = x.rows (); 1187 octave_idx_type x_nc = x.cols (); 1188 1189 octave_idx_type alpha_nr = alpha.rows (); 1190 octave_idx_type alpha_nc = alpha.cols (); 1191 1192 if (x_nr != alpha_nr || x_nc != alpha_nc) 1193 (*current_liboctave_error_handler) 1194 ("%s: the sizes of alpha and x must conform", fn); 1195 1196 octave_idx_type nr = x_nr; 1197 octave_idx_type nc = x_nc; 1198 1199 retval.resize (nr, nc); 1200 1201 ierr.resize (dim_vector (nr, nc)); 1202 1203 for (octave_idx_type j = 0; j < nc; j++) 1204 for (octave_idx_type i = 0; i < nr; i++) 1205 retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j)); 1206 1207 return retval; 1208 } 1209 1210 static inline FloatComplexNDArray do_bessel(fptr f,const char *,float alpha,const FloatComplexNDArray & x,bool scaled,Array<octave_idx_type> & ierr)1211 do_bessel (fptr f, const char *, float alpha, const FloatComplexNDArray& x, 1212 bool scaled, Array<octave_idx_type>& ierr) 1213 { 1214 dim_vector dv = x.dims (); 1215 octave_idx_type nel = dv.numel (); 1216 FloatComplexNDArray retval (dv); 1217 1218 ierr.resize (dv); 1219 1220 for (octave_idx_type i = 0; i < nel; i++) 1221 retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i)); 1222 1223 return retval; 1224 } 1225 1226 static inline FloatComplexNDArray do_bessel(fptr f,const char *,const FloatNDArray & alpha,const FloatComplex & x,bool scaled,Array<octave_idx_type> & ierr)1227 do_bessel (fptr f, const char *, const FloatNDArray& alpha, 1228 const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) 1229 { 1230 dim_vector dv = alpha.dims (); 1231 octave_idx_type nel = dv.numel (); 1232 FloatComplexNDArray retval (dv); 1233 1234 ierr.resize (dv); 1235 1236 for (octave_idx_type i = 0; i < nel; i++) 1237 retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i)); 1238 1239 return retval; 1240 } 1241 1242 static inline FloatComplexNDArray do_bessel(fptr f,const char * fn,const FloatNDArray & alpha,const FloatComplexNDArray & x,bool scaled,Array<octave_idx_type> & ierr)1243 do_bessel (fptr f, const char *fn, const FloatNDArray& alpha, 1244 const FloatComplexNDArray& x, bool scaled, 1245 Array<octave_idx_type>& ierr) 1246 { 1247 dim_vector dv = x.dims (); 1248 FloatComplexNDArray retval; 1249 1250 if (dv != alpha.dims ()) 1251 (*current_liboctave_error_handler) 1252 ("%s: the sizes of alpha and x must conform", fn); 1253 1254 octave_idx_type nel = dv.numel (); 1255 1256 retval.resize (dv); 1257 ierr.resize (dv); 1258 1259 for (octave_idx_type i = 0; i < nel; i++) 1260 retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i)); 1261 1262 return retval; 1263 } 1264 1265 static inline FloatComplexMatrix do_bessel(fptr f,const char *,const FloatRowVector & alpha,const FloatComplexColumnVector & x,bool scaled,Array<octave_idx_type> & ierr)1266 do_bessel (fptr f, const char *, const FloatRowVector& alpha, 1267 const FloatComplexColumnVector& x, bool scaled, 1268 Array<octave_idx_type>& ierr) 1269 { 1270 octave_idx_type nr = x.numel (); 1271 octave_idx_type nc = alpha.numel (); 1272 1273 FloatComplexMatrix retval (nr, nc); 1274 1275 ierr.resize (dim_vector (nr, nc)); 1276 1277 for (octave_idx_type j = 0; j < nc; j++) 1278 for (octave_idx_type i = 0; i < nr; i++) 1279 retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j)); 1280 1281 return retval; 1282 } 1283 1284 #define SS_BESSEL(name, fcn) \ 1285 FloatComplex \ 1286 name (float alpha, const FloatComplex& x, bool scaled, \ 1287 octave_idx_type& ierr) \ 1288 { \ 1289 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 1290 } 1291 1292 #define SM_BESSEL(name, fcn) \ 1293 FloatComplexMatrix \ 1294 name (float alpha, const FloatComplexMatrix& x, bool scaled, \ 1295 Array<octave_idx_type>& ierr) \ 1296 { \ 1297 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 1298 } 1299 1300 #define MS_BESSEL(name, fcn) \ 1301 FloatComplexMatrix \ 1302 name (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, \ 1303 Array<octave_idx_type>& ierr) \ 1304 { \ 1305 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 1306 } 1307 1308 #define MM_BESSEL(name, fcn) \ 1309 FloatComplexMatrix \ 1310 name (const FloatMatrix& alpha, const FloatComplexMatrix& x, \ 1311 bool scaled, Array<octave_idx_type>& ierr) \ 1312 { \ 1313 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 1314 } 1315 1316 #define SN_BESSEL(name, fcn) \ 1317 FloatComplexNDArray \ 1318 name (float alpha, const FloatComplexNDArray& x, bool scaled, \ 1319 Array<octave_idx_type>& ierr) \ 1320 { \ 1321 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 1322 } 1323 1324 #define NS_BESSEL(name, fcn) \ 1325 FloatComplexNDArray \ 1326 name (const FloatNDArray& alpha, const FloatComplex& x, \ 1327 bool scaled, Array<octave_idx_type>& ierr) \ 1328 { \ 1329 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 1330 } 1331 1332 #define NN_BESSEL(name, fcn) \ 1333 FloatComplexNDArray \ 1334 name (const FloatNDArray& alpha, const FloatComplexNDArray& x, \ 1335 bool scaled, Array<octave_idx_type>& ierr) \ 1336 { \ 1337 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 1338 } 1339 1340 #define RC_BESSEL(name, fcn) \ 1341 FloatComplexMatrix \ 1342 name (const FloatRowVector& alpha, \ 1343 const FloatComplexColumnVector& x, bool scaled, \ 1344 Array<octave_idx_type>& ierr) \ 1345 { \ 1346 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ 1347 } 1348 1349 #define ALL_BESSEL(name, fcn) \ 1350 SS_BESSEL (name, fcn) \ 1351 SM_BESSEL (name, fcn) \ 1352 MS_BESSEL (name, fcn) \ 1353 MM_BESSEL (name, fcn) \ 1354 SN_BESSEL (name, fcn) \ 1355 NS_BESSEL (name, fcn) \ 1356 NN_BESSEL (name, fcn) \ 1357 RC_BESSEL (name, fcn) 1358 ALL_BESSEL(besselj,cbesj)1359 ALL_BESSEL (besselj, cbesj) 1360 ALL_BESSEL (bessely, cbesy) 1361 ALL_BESSEL (besseli, cbesi) 1362 ALL_BESSEL (besselk, cbesk) 1363 ALL_BESSEL (besselh1, cbesh1) 1364 ALL_BESSEL (besselh2, cbesh2) 1365 1366 #undef ALL_BESSEL 1367 #undef SS_BESSEL 1368 #undef SM_BESSEL 1369 #undef MS_BESSEL 1370 #undef MM_BESSEL 1371 #undef SN_BESSEL 1372 #undef NS_BESSEL 1373 #undef NN_BESSEL 1374 #undef RC_BESSEL 1375 1376 Complex 1377 biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr) 1378 { 1379 double ar = 0.0; 1380 double ai = 0.0; 1381 1382 double zr = z.real (); 1383 double zi = z.imag (); 1384 1385 F77_INT id = (deriv ? 1 : 0); 1386 F77_INT t_ierr; 1387 1388 F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 2, ar, ai, t_ierr); 1389 1390 ierr = t_ierr; 1391 1392 if (! scaled) 1393 { 1394 Complex expz = exp (std::abs (std::real (2.0 / 3.0 * z * sqrt (z)))); 1395 1396 double rexpz = expz.real (); 1397 double iexpz = expz.imag (); 1398 1399 double tmp = ar*rexpz - ai*iexpz; 1400 1401 ai = ar*iexpz + ai*rexpz; 1402 ar = tmp; 1403 } 1404 1405 if (zi == 0.0 && (! scaled || zr >= 0.0)) 1406 ai = 0.0; 1407 1408 return bessel_return_value (Complex (ar, ai), ierr); 1409 } 1410 1411 ComplexMatrix biry(const ComplexMatrix & z,bool deriv,bool scaled,Array<octave_idx_type> & ierr)1412 biry (const ComplexMatrix& z, bool deriv, bool scaled, 1413 Array<octave_idx_type>& ierr) 1414 { 1415 octave_idx_type nr = z.rows (); 1416 octave_idx_type nc = z.cols (); 1417 1418 ComplexMatrix retval (nr, nc); 1419 1420 ierr.resize (dim_vector (nr, nc)); 1421 1422 for (octave_idx_type j = 0; j < nc; j++) 1423 for (octave_idx_type i = 0; i < nr; i++) 1424 retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j)); 1425 1426 return retval; 1427 } 1428 1429 ComplexNDArray biry(const ComplexNDArray & z,bool deriv,bool scaled,Array<octave_idx_type> & ierr)1430 biry (const ComplexNDArray& z, bool deriv, bool scaled, 1431 Array<octave_idx_type>& ierr) 1432 { 1433 dim_vector dv = z.dims (); 1434 octave_idx_type nel = dv.numel (); 1435 ComplexNDArray retval (dv); 1436 1437 ierr.resize (dv); 1438 1439 for (octave_idx_type i = 0; i < nel; i++) 1440 retval(i) = biry (z(i), deriv, scaled, ierr(i)); 1441 1442 return retval; 1443 } 1444 1445 FloatComplex biry(const FloatComplex & z,bool deriv,bool scaled,octave_idx_type & ierr)1446 biry (const FloatComplex& z, bool deriv, bool scaled, 1447 octave_idx_type& ierr) 1448 { 1449 FloatComplex a; 1450 1451 F77_INT id = (deriv ? 1 : 0); 1452 F77_INT t_ierr; 1453 1454 F77_FUNC (cbiry, CBIRY) (F77_CONST_CMPLX_ARG (&z), id, 2, 1455 F77_CMPLX_ARG (&a), t_ierr); 1456 1457 ierr = t_ierr; 1458 1459 float ar = a.real (); 1460 float ai = a.imag (); 1461 1462 if (! scaled) 1463 { 1464 FloatComplex expz 1465 = exp (std::abs (std::real (2.0f / 3.0f * z * sqrt (z)))); 1466 1467 float rexpz = expz.real (); 1468 float iexpz = expz.imag (); 1469 1470 float tmp = ar*rexpz - ai*iexpz; 1471 1472 ai = ar*iexpz + ai*rexpz; 1473 ar = tmp; 1474 } 1475 1476 if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0)) 1477 ai = 0.0; 1478 1479 return bessel_return_value (FloatComplex (ar, ai), ierr); 1480 } 1481 1482 FloatComplexMatrix biry(const FloatComplexMatrix & z,bool deriv,bool scaled,Array<octave_idx_type> & ierr)1483 biry (const FloatComplexMatrix& z, bool deriv, bool scaled, 1484 Array<octave_idx_type>& ierr) 1485 { 1486 octave_idx_type nr = z.rows (); 1487 octave_idx_type nc = z.cols (); 1488 1489 FloatComplexMatrix retval (nr, nc); 1490 1491 ierr.resize (dim_vector (nr, nc)); 1492 1493 for (octave_idx_type j = 0; j < nc; j++) 1494 for (octave_idx_type i = 0; i < nr; i++) 1495 retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j)); 1496 1497 return retval; 1498 } 1499 1500 FloatComplexNDArray biry(const FloatComplexNDArray & z,bool deriv,bool scaled,Array<octave_idx_type> & ierr)1501 biry (const FloatComplexNDArray& z, bool deriv, bool scaled, 1502 Array<octave_idx_type>& ierr) 1503 { 1504 dim_vector dv = z.dims (); 1505 octave_idx_type nel = dv.numel (); 1506 FloatComplexNDArray retval (dv); 1507 1508 ierr.resize (dv); 1509 1510 for (octave_idx_type i = 0; i < nel; i++) 1511 retval(i) = biry (z(i), deriv, scaled, ierr(i)); 1512 1513 return retval; 1514 } 1515 1516 // Real and complex Dawson function (= scaled erfi) from Faddeeva package dawson(double x)1517 double dawson (double x) { return Faddeeva::Dawson (x); } dawson(float x)1518 float dawson (float x) { return Faddeeva::Dawson (x); } 1519 1520 Complex dawson(const Complex & x)1521 dawson (const Complex& x) 1522 { 1523 return Faddeeva::Dawson (x); 1524 } 1525 1526 FloatComplex dawson(const FloatComplex & x)1527 dawson (const FloatComplex& x) 1528 { 1529 Complex xd (x.real (), x.imag ()); 1530 Complex ret; 1531 ret = Faddeeva::Dawson (xd, std::numeric_limits<float>::epsilon ()); 1532 return FloatComplex (ret.real (), ret.imag ()); 1533 } 1534 1535 void ellipj(double u,double m,double & sn,double & cn,double & dn,double & err)1536 ellipj (double u, double m, double& sn, double& cn, double& dn, double& err) 1537 { 1538 static const int Nmax = 16; 1539 double m1, t=0, si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi; 1540 int n, Nn, ii; 1541 1542 if (m < 0 || m > 1) 1543 { 1544 (*current_liboctave_warning_with_id_handler) 1545 ("Octave:ellipj-invalid-m", 1546 "ellipj: invalid M value, required value 0 <= M <= 1"); 1547 1548 sn = cn = dn = lo_ieee_nan_value (); 1549 1550 return; 1551 } 1552 1553 double sqrt_eps = std::sqrt (std::numeric_limits<double>::epsilon ()); 1554 if (m < sqrt_eps) 1555 { 1556 // For small m, (Abramowitz and Stegun, Section 16.13) 1557 si_u = sin (u); 1558 co_u = cos (u); 1559 t = 0.25*m*(u - si_u*co_u); 1560 sn = si_u - t * co_u; 1561 cn = co_u + t * si_u; 1562 dn = 1 - 0.5*m*si_u*si_u; 1563 } 1564 else if ((1 - m) < sqrt_eps) 1565 { 1566 // For m1 = (1-m) small (Abramowitz and Stegun, Section 16.15) 1567 m1 = 1 - m; 1568 si_u = sinh (u); 1569 co_u = cosh (u); 1570 ta_u = tanh (u); 1571 se_u = 1/co_u; 1572 sn = ta_u + 0.25*m1*(si_u*co_u - u)*se_u*se_u; 1573 cn = se_u - 0.25*m1*(si_u*co_u - u)*ta_u*se_u; 1574 dn = se_u + 0.25*m1*(si_u*co_u + u)*ta_u*se_u; 1575 } 1576 else 1577 { 1578 // Arithmetic-Geometric Mean (AGM) algorithm 1579 // (Abramowitz and Stegun, Section 16.4) 1580 a[0] = 1; 1581 b = std::sqrt (1 - m); 1582 c[0] = std::sqrt (m); 1583 for (n = 1; n < Nmax; ++n) 1584 { 1585 a[n] = (a[n - 1] + b)/2; 1586 c[n] = (a[n - 1] - b)/2; 1587 b = std::sqrt (a[n - 1]*b); 1588 if (c[n]/a[n] < std::numeric_limits<double>::epsilon ()) break; 1589 } 1590 if (n >= Nmax - 1) 1591 { 1592 err = 1; 1593 return; 1594 } 1595 Nn = n; 1596 for (ii = 1; n > 0; ii *= 2, --n) {} // ii = pow(2,Nn) 1597 phi = ii*a[Nn]*u; 1598 for (n = Nn; n > 0; --n) 1599 { 1600 phi = (std::asin ((c[n]/a[n])* sin (phi)) + phi)/2; 1601 } 1602 sn = sin (phi); 1603 cn = cos (phi); 1604 dn = std::sqrt (1 - m*sn*sn); 1605 } 1606 } 1607 1608 void ellipj(const Complex & u,double m,Complex & sn,Complex & cn,Complex & dn,double & err)1609 ellipj (const Complex& u, double m, Complex& sn, Complex& cn, Complex& dn, 1610 double& err) 1611 { 1612 double m1 = 1 - m, ss1, cc1, dd1; 1613 1614 ellipj (u.imag (), m1, ss1, cc1, dd1, err); 1615 if (u.real () == 0) 1616 { 1617 // u is pure imag: Jacoby imag. transf. 1618 sn = Complex (0, ss1/cc1); 1619 cn = 1/cc1; // cn.imag = 0; 1620 dn = dd1/cc1; // dn.imag = 0; 1621 } 1622 else 1623 { 1624 // u is generic complex 1625 double ss, cc, dd, ddd; 1626 1627 ellipj (u.real (), m, ss, cc, dd, err); 1628 ddd = cc1*cc1 + m*ss*ss*ss1*ss1; 1629 sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd); 1630 cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd); 1631 dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd); 1632 } 1633 } 1634 1635 // Complex error function from the Faddeeva package 1636 Complex erf(const Complex & x)1637 erf (const Complex& x) 1638 { 1639 return Faddeeva::erf (x); 1640 } 1641 1642 FloatComplex erf(const FloatComplex & x)1643 erf (const FloatComplex& x) 1644 { 1645 Complex xd (x.real (), x.imag ()); 1646 Complex ret = Faddeeva::erf (xd, std::numeric_limits<float>::epsilon ()); 1647 return FloatComplex (ret.real (), ret.imag ()); 1648 } 1649 1650 // Complex complementary error function from the Faddeeva package 1651 Complex erfc(const Complex & x)1652 erfc (const Complex& x) 1653 { 1654 return Faddeeva::erfc (x); 1655 } 1656 1657 FloatComplex erfc(const FloatComplex & x)1658 erfc (const FloatComplex& x) 1659 { 1660 Complex xd (x.real (), x.imag ()); 1661 Complex ret = Faddeeva::erfc (xd, std::numeric_limits<float>::epsilon ()); 1662 return FloatComplex (ret.real (), ret.imag ()); 1663 } 1664 1665 // The algorithm for erfcinv is an adaptation of the erfinv algorithm 1666 // above from P. J. Acklam. It has been modified to run over the 1667 // different input domain of erfcinv. See the notes for erfinv for an 1668 // explanation. 1669 do_erfcinv(double x,bool refine)1670 static double do_erfcinv (double x, bool refine) 1671 { 1672 // Coefficients of rational approximation. 1673 static const double a[] = 1674 { 1675 -2.806989788730439e+01, 1.562324844726888e+02, 1676 -1.951109208597547e+02, 9.783370457507161e+01, 1677 -2.168328665628878e+01, 1.772453852905383e+00 1678 }; 1679 static const double b[] = 1680 { 1681 -5.447609879822406e+01, 1.615858368580409e+02, 1682 -1.556989798598866e+02, 6.680131188771972e+01, 1683 -1.328068155288572e+01 1684 }; 1685 static const double c[] = 1686 { 1687 -5.504751339936943e-03, -2.279687217114118e-01, 1688 -1.697592457770869e+00, -1.802933168781950e+00, 1689 3.093354679843505e+00, 2.077595676404383e+00 1690 }; 1691 static const double d[] = 1692 { 1693 7.784695709041462e-03, 3.224671290700398e-01, 1694 2.445134137142996e+00, 3.754408661907416e+00 1695 }; 1696 1697 static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. 1698 static const double pbreak_lo = 0.04850; // 1-pbreak 1699 static const double pbreak_hi = 1.95150; // 1+pbreak 1700 double y; 1701 1702 // Select case. 1703 if (x >= pbreak_lo && x <= pbreak_hi) 1704 { 1705 // Middle region. 1706 const double q = 0.5*(1-x), r = q*q; 1707 const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q; 1708 const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0; 1709 y = yn / yd; 1710 } 1711 else if (x > 0.0 && x < 2.0) 1712 { 1713 // Tail region. 1714 const double q = (x < 1 1715 ? std::sqrt (-2*std::log (0.5*x)) 1716 : std::sqrt (-2*std::log (0.5*(2-x)))); 1717 1718 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]; 1719 1720 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0; 1721 1722 y = yn / yd; 1723 1724 if (x < pbreak_lo) 1725 y = -y; 1726 } 1727 else if (x == 0.0) 1728 return numeric_limits<double>::Inf (); 1729 else if (x == 2.0) 1730 return -numeric_limits<double>::Inf (); 1731 else 1732 return numeric_limits<double>::NaN (); 1733 1734 if (refine) 1735 { 1736 // One iteration of Halley's method gives full precision. 1737 double u = (erf (y) - (1-x)) * spi2 * exp (y*y); 1738 y -= u / (1 + y*u); 1739 } 1740 1741 return y; 1742 } 1743 erfcinv(double x)1744 double erfcinv (double x) 1745 { 1746 return do_erfcinv (x, true); 1747 } 1748 erfcinv(float x)1749 float erfcinv (float x) 1750 { 1751 return do_erfcinv (x, false); 1752 } 1753 1754 // Real and complex scaled complementary error function from Faddeeva pkg. erfcx(double x)1755 double erfcx (double x) { return Faddeeva::erfcx (x); } erfcx(float x)1756 float erfcx (float x) { return Faddeeva::erfcx (x); } 1757 1758 Complex erfcx(const Complex & x)1759 erfcx (const Complex& x) 1760 { 1761 return Faddeeva::erfcx (x); 1762 } 1763 1764 FloatComplex erfcx(const FloatComplex & x)1765 erfcx (const FloatComplex& x) 1766 { 1767 Complex xd (x.real (), x.imag ()); 1768 Complex ret; 1769 ret = Faddeeva::erfcx (xd, std::numeric_limits<float>::epsilon ()); 1770 return FloatComplex (ret.real (), ret.imag ()); 1771 } 1772 1773 // Real and complex imaginary error function from Faddeeva package erfi(double x)1774 double erfi (double x) { return Faddeeva::erfi (x); } erfi(float x)1775 float erfi (float x) { return Faddeeva::erfi (x); } 1776 1777 Complex erfi(const Complex & x)1778 erfi (const Complex& x) 1779 { 1780 return Faddeeva::erfi (x); 1781 } 1782 1783 FloatComplex erfi(const FloatComplex & x)1784 erfi (const FloatComplex& x) 1785 { 1786 Complex xd (x.real (), x.imag ()); 1787 Complex ret = Faddeeva::erfi (xd, std::numeric_limits<float>::epsilon ()); 1788 return FloatComplex (ret.real (), ret.imag ()); 1789 } 1790 1791 // This algorithm is due to P. J. Acklam. 1792 // 1793 // See http://home.online.no/~pjacklam/notes/invnorm/ 1794 // 1795 // The rational approximation has relative accuracy 1.15e-9 in the whole 1796 // region. For doubles, it is refined by a single step of Halley's 3rd 1797 // order method. For single precision, the accuracy is already OK, so 1798 // we skip it to get faster evaluation. 1799 do_erfinv(double x,bool refine)1800 static double do_erfinv (double x, bool refine) 1801 { 1802 // Coefficients of rational approximation. 1803 static const double a[] = 1804 { 1805 -2.806989788730439e+01, 1.562324844726888e+02, 1806 -1.951109208597547e+02, 9.783370457507161e+01, 1807 -2.168328665628878e+01, 1.772453852905383e+00 1808 }; 1809 static const double b[] = 1810 { 1811 -5.447609879822406e+01, 1.615858368580409e+02, 1812 -1.556989798598866e+02, 6.680131188771972e+01, 1813 -1.328068155288572e+01 1814 }; 1815 static const double c[] = 1816 { 1817 -5.504751339936943e-03, -2.279687217114118e-01, 1818 -1.697592457770869e+00, -1.802933168781950e+00, 1819 3.093354679843505e+00, 2.077595676404383e+00 1820 }; 1821 static const double d[] = 1822 { 1823 7.784695709041462e-03, 3.224671290700398e-01, 1824 2.445134137142996e+00, 3.754408661907416e+00 1825 }; 1826 1827 static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2. 1828 static const double pbreak = 0.95150; 1829 double ax = fabs (x), y; 1830 1831 // Select case. 1832 if (ax <= pbreak) 1833 { 1834 // Middle region. 1835 const double q = 0.5 * x, r = q*q; 1836 const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q; 1837 const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0; 1838 y = yn / yd; 1839 } 1840 else if (ax < 1.0) 1841 { 1842 // Tail region. 1843 const double q = std::sqrt (-2*std::log (0.5*(1-ax))); 1844 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]; 1845 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0; 1846 y = yn / yd * math::signum (-x); 1847 } 1848 else if (ax == 1.0) 1849 return numeric_limits<double>::Inf () * math::signum (x); 1850 else 1851 return numeric_limits<double>::NaN (); 1852 1853 if (refine) 1854 { 1855 // One iteration of Halley's method gives full precision. 1856 double u = (erf (y) - x) * spi2 * exp (y*y); 1857 y -= u / (1 + y*u); 1858 } 1859 1860 return y; 1861 } 1862 erfinv(double x)1863 double erfinv (double x) 1864 { 1865 return do_erfinv (x, true); 1866 } 1867 erfinv(float x)1868 float erfinv (float x) 1869 { 1870 return do_erfinv (x, false); 1871 } 1872 1873 Complex expm1(const Complex & x)1874 expm1 (const Complex& x) 1875 { 1876 Complex retval; 1877 1878 if (std::abs (x) < 1) 1879 { 1880 double im = x.imag (); 1881 double u = expm1 (x.real ()); 1882 double v = sin (im/2); 1883 v = -2*v*v; 1884 retval = Complex (u*v + u + v, (u+1) * sin (im)); 1885 } 1886 else 1887 retval = std::exp (x) - Complex (1); 1888 1889 return retval; 1890 } 1891 1892 FloatComplex expm1(const FloatComplex & x)1893 expm1 (const FloatComplex& x) 1894 { 1895 FloatComplex retval; 1896 1897 if (std::abs (x) < 1) 1898 { 1899 float im = x.imag (); 1900 float u = expm1 (x.real ()); 1901 float v = sin (im/2); 1902 v = -2*v*v; 1903 retval = FloatComplex (u*v + u + v, (u+1) * sin (im)); 1904 } 1905 else 1906 retval = std::exp (x) - FloatComplex (1); 1907 1908 return retval; 1909 } 1910 1911 double gamma(double x)1912 gamma (double x) 1913 { 1914 double result; 1915 1916 // Special cases for (near) compatibility with Matlab instead of tgamma. 1917 // Matlab does not have -0. 1918 1919 if (x == 0) 1920 result = (math::negative_sign (x) 1921 ? -numeric_limits<double>::Inf () 1922 : numeric_limits<double>::Inf ()); 1923 else if ((x < 0 && math::x_nint (x) == x) 1924 || math::isinf (x)) 1925 result = numeric_limits<double>::Inf (); 1926 else if (math::isnan (x)) 1927 result = numeric_limits<double>::NaN (); 1928 else 1929 result = std::tgamma (x); 1930 1931 return result; 1932 } 1933 1934 float gamma(float x)1935 gamma (float x) 1936 { 1937 float result; 1938 1939 // Special cases for (near) compatibility with Matlab instead of tgamma. 1940 // Matlab does not have -0. 1941 1942 if (x == 0) 1943 result = (math::negative_sign (x) 1944 ? -numeric_limits<float>::Inf () 1945 : numeric_limits<float>::Inf ()); 1946 else if ((x < 0 && math::x_nint (x) == x) 1947 || math::isinf (x)) 1948 result = numeric_limits<float>::Inf (); 1949 else if (math::isnan (x)) 1950 result = numeric_limits<float>::NaN (); 1951 else 1952 result = std::tgammaf (x); 1953 1954 return result; 1955 } 1956 1957 Complex log1p(const Complex & x)1958 log1p (const Complex& x) 1959 { 1960 Complex retval; 1961 1962 double r = x.real (), i = x.imag (); 1963 1964 if (fabs (r) < 0.5 && fabs (i) < 0.5) 1965 { 1966 double u = 2*r + r*r + i*i; 1967 retval = Complex (log1p (u / (1+std::sqrt (u+1))), 1968 atan2 (1 + r, i)); 1969 } 1970 else 1971 retval = std::log (Complex (1) + x); 1972 1973 return retval; 1974 } 1975 1976 FloatComplex log1p(const FloatComplex & x)1977 log1p (const FloatComplex& x) 1978 { 1979 FloatComplex retval; 1980 1981 float r = x.real (), i = x.imag (); 1982 1983 if (fabs (r) < 0.5 && fabs (i) < 0.5) 1984 { 1985 float u = 2*r + r*r + i*i; 1986 retval = FloatComplex (log1p (u / (1+std::sqrt (u+1))), 1987 atan2 (1 + r, i)); 1988 } 1989 else 1990 retval = std::log (FloatComplex (1) + x); 1991 1992 return retval; 1993 } 1994 1995 static const double pi = 3.14159265358979323846; 1996 1997 template <typename T> 1998 static inline T xlog(const T & x)1999 xlog (const T& x) 2000 { 2001 return log (x); 2002 } 2003 2004 template <> 2005 inline double xlog(const double & x)2006 xlog (const double& x) 2007 { 2008 return std::log (x); 2009 } 2010 2011 template <> 2012 inline float xlog(const float & x)2013 xlog (const float& x) 2014 { 2015 return std::log (x); 2016 } 2017 2018 template <typename T> 2019 static T lanczos_approximation_psi(const T zc)2020 lanczos_approximation_psi (const T zc) 2021 { 2022 // Coefficients for C.Lanczos expansion of psi function from XLiFE++ 2023 // gammaFunctions psi_coef[k] = - (2k+1) * lg_coef[k] (see melina++ 2024 // gamma functions -1/12, 3/360,-5/1260, 7/1680,-9/1188, 2025 // 11*691/360360,-13/156, 15*3617/122400, ? , ? 2026 static const T dg_coeff[10] = 2027 { 2028 -0.83333333333333333e-1, 0.83333333333333333e-2, 2029 -0.39682539682539683e-2, 0.41666666666666667e-2, 2030 -0.75757575757575758e-2, 0.21092796092796093e-1, 2031 -0.83333333333333333e-1, 0.4432598039215686, 2032 -0.3053954330270122e+1, 0.125318899521531e+2 2033 }; 2034 2035 T overz2 = T (1.0) / (zc * zc); 2036 T overz2k = overz2; 2037 2038 T p = 0; 2039 for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2) 2040 p += dg_coeff[k] * overz2k; 2041 p += xlog (zc) - T (0.5) / zc; 2042 return p; 2043 } 2044 2045 template <typename T> 2046 T xpsi(T z)2047 xpsi (T z) 2048 { 2049 static const double euler_mascheroni 2050 = 0.577215664901532860606512090082402431042; 2051 2052 const bool is_int = (std::floor (z) == z); 2053 2054 T p = 0; 2055 if (z <= 0) 2056 { 2057 // limits - zeros of the gamma function 2058 if (is_int) 2059 p = -numeric_limits<T>::Inf (); // Matlab returns -Inf for psi (0) 2060 else 2061 // Abramowitz and Stegun, page 259, eq 6.3.7 2062 p = psi (1 - z) - (pi / tan (pi * z)); 2063 } 2064 else if (is_int) 2065 { 2066 // Abramowitz and Stegun, page 258, eq 6.3.2 2067 p = - euler_mascheroni; 2068 for (octave_idx_type k = z - 1; k > 0; k--) 2069 p += 1.0 / k; 2070 } 2071 else if (std::floor (z + 0.5) == z + 0.5) 2072 { 2073 // Abramowitz and Stegun, page 258, eq 6.3.3 and 6.3.4 2074 for (octave_idx_type k = z; k > 0; k--) 2075 p += 1.0 / (2 * k - 1); 2076 2077 p = - euler_mascheroni - 2 * std::log (2) + 2 * (p); 2078 } 2079 else 2080 { 2081 // adapted from XLiFE++ gammaFunctions 2082 2083 T zc = z; 2084 // Use formula for derivative of LogGamma(z) 2085 if (z < 10) 2086 { 2087 const signed char n = 10 - z; 2088 for (signed char k = n - 1; k >= 0; k--) 2089 p -= 1.0 / (k + z); 2090 zc += n; 2091 } 2092 p += lanczos_approximation_psi (zc); 2093 } 2094 2095 return p; 2096 } 2097 2098 // explicit instantiations psi(double z)2099 double psi (double z) { return xpsi (z); } psi(float z)2100 float psi (float z) { return xpsi (z); } 2101 2102 template <typename T> 2103 std::complex<T> xpsi(const std::complex<T> & z)2104 xpsi (const std::complex<T>& z) 2105 { 2106 // adapted from XLiFE++ gammaFunctions 2107 2108 typedef typename std::complex<T>::value_type P; 2109 2110 P z_r = z.real (); 2111 P z_ra = z_r; 2112 2113 std::complex<T> dgam (0.0, 0.0); 2114 if (z.imag () == 0) 2115 dgam = std::complex<T> (psi (z_r), 0.0); 2116 else if (z_r < 0) 2117 dgam = psi (P (1.0) - z)- (P (pi) / tan (P (pi) * z)); 2118 else 2119 { 2120 // Use formula for derivative of LogGamma(z) 2121 std::complex<T> z_m = z; 2122 if (z_ra < 8) 2123 { 2124 unsigned char n = 8 - z_ra; 2125 z_m = z + std::complex<T> (n, 0.0); 2126 2127 // Recurrence formula. For | Re(z) | < 8, use recursively 2128 // 2129 // DiGamma(z) = DiGamma(z+1) - 1/z 2130 std::complex<T> z_p = z + P (n - 1); 2131 for (unsigned char k = n; k > 0; k--, z_p -= 1.0) 2132 dgam -= P (1.0) / z_p; 2133 } 2134 2135 // for | Re(z) | > 8, use derivative of C.Lanczos expansion for 2136 // LogGamma 2137 // 2138 // psi(z) = log(z) - 1/(2z) - 1/12z^2 + 3/360z^4 - 5/1260z^6 2139 // + 7/1680z^8 - 9/1188z^10 + ... 2140 // 2141 // (Abramowitz&Stegun, page 259, formula 6.3.18 2142 dgam += lanczos_approximation_psi (z_m); 2143 } 2144 return dgam; 2145 } 2146 2147 // explicit instantiations psi(const Complex & z)2148 Complex psi (const Complex& z) { return xpsi (z); } psi(const FloatComplex & z)2149 FloatComplex psi (const FloatComplex& z) { return xpsi (z); } 2150 2151 template <typename T> 2152 static inline void 2153 fortran_psifn (T z, octave_idx_type n, T& ans, octave_idx_type& ierr); 2154 2155 template <> 2156 inline void fortran_psifn(double z,octave_idx_type n_arg,double & ans,octave_idx_type & ierr)2157 fortran_psifn<double> (double z, octave_idx_type n_arg, 2158 double& ans, octave_idx_type& ierr) 2159 { 2160 F77_INT n = to_f77_int (n_arg); 2161 F77_INT flag = 0; 2162 F77_INT t_ierr; 2163 F77_XFCN (dpsifn, DPSIFN, (z, n, 1, 1, ans, flag, t_ierr)); 2164 ierr = t_ierr; 2165 } 2166 2167 template <> 2168 inline void fortran_psifn(float z,octave_idx_type n_arg,float & ans,octave_idx_type & ierr)2169 fortran_psifn<float> (float z, octave_idx_type n_arg, 2170 float& ans, octave_idx_type& ierr) 2171 { 2172 F77_INT n = to_f77_int (n_arg); 2173 F77_INT flag = 0; 2174 F77_INT t_ierr; 2175 F77_XFCN (psifn, PSIFN, (z, n, 1, 1, ans, flag, t_ierr)); 2176 ierr = t_ierr; 2177 } 2178 2179 template <typename T> 2180 T xpsi(octave_idx_type n,T z)2181 xpsi (octave_idx_type n, T z) 2182 { 2183 T ans; 2184 octave_idx_type ierr = 0; 2185 fortran_psifn<T> (z, n, ans, ierr); 2186 if (ierr == 0) 2187 { 2188 // Remember that psifn and dpsifn return scales values 2189 // When n is 1: do nothing since ((-1)**(n+1)/gamma(n+1)) == 1 2190 // When n is 0: change sign since ((-1)**(n+1)/gamma(n+1)) == -1 2191 if (n > 1) 2192 // FIXME: xgamma here is a killer for our precision since it grows 2193 // way too fast. 2194 ans = ans / (std::pow (-1.0, n + 1) / gamma (double (n+1))); 2195 else if (n == 0) 2196 ans = -ans; 2197 } 2198 else if (ierr == 2) 2199 ans = - numeric_limits<T>::Inf (); 2200 else // we probably never get here 2201 ans = numeric_limits<T>::NaN (); 2202 2203 return ans; 2204 } 2205 psi(octave_idx_type n,double z)2206 double psi (octave_idx_type n, double z) { return xpsi (n, z); } psi(octave_idx_type n,float z)2207 float psi (octave_idx_type n, float z) { return xpsi (n, z); } 2208 2209 Complex rc_lgamma(double x)2210 rc_lgamma (double x) 2211 { 2212 double result; 2213 2214 #if defined (HAVE_LGAMMA_R) 2215 int sgngam; 2216 result = lgamma_r (x, &sgngam); 2217 #else 2218 result = std::lgamma (x); 2219 int sgngam = signgam; 2220 #endif 2221 2222 if (sgngam < 0) 2223 return result + Complex (0., M_PI); 2224 else 2225 return result; 2226 } 2227 2228 FloatComplex rc_lgamma(float x)2229 rc_lgamma (float x) 2230 { 2231 float result; 2232 2233 #if defined (HAVE_LGAMMAF_R) 2234 int sgngam; 2235 result = lgammaf_r (x, &sgngam); 2236 #else 2237 result = std::lgammaf (x); 2238 int sgngam = signgam; 2239 #endif 2240 2241 if (sgngam < 0) 2242 return result + FloatComplex (0., M_PI); 2243 else 2244 return result; 2245 } 2246 rc_log1p(double x)2247 Complex rc_log1p (double x) 2248 { 2249 return (x < -1.0 2250 ? Complex (std::log (-(1.0 + x)), M_PI) 2251 : Complex (log1p (x))); 2252 } 2253 rc_log1p(float x)2254 FloatComplex rc_log1p (float x) 2255 { 2256 return (x < -1.0f 2257 ? FloatComplex (std::log (-(1.0f + x)), M_PI) 2258 : FloatComplex (log1p (x))); 2259 } 2260 } 2261 } 2262