00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026
00027 #include "Range.h"
00028 #include "CColVector.h"
00029 #include "CMatrix.h"
00030 #include "dRowVector.h"
00031 #include "dMatrix.h"
00032 #include "dNDArray.h"
00033 #include "CNDArray.h"
00034 #include "f77-fcn.h"
00035 #include "lo-error.h"
00036 #include "lo-ieee.h"
00037 #include "lo-specfun.h"
00038 #include "mx-inlines.cc"
00039
00040 #ifndef M_PI
00041 #define M_PI 3.14159265358979323846
00042 #endif
00043
00044 extern "C"
00045 {
00046 F77_RET_T
00047 F77_FUNC (zbesj, ZBESJ) (const double&, const double&, const double&,
00048 const int&, const int&, double*, double*,
00049 int&, int&);
00050
00051 F77_RET_T
00052 F77_FUNC (zbesy, ZBESY) (const double&, const double&, const double&,
00053 const int&, const int&, double*, double*,
00054 int&, double*, double*, int&);
00055
00056 F77_RET_T
00057 F77_FUNC (zbesi, ZBESI) (const double&, const double&, const double&,
00058 const int&, const int&, double*, double*,
00059 int&, int&);
00060
00061 F77_RET_T
00062 F77_FUNC (zbesk, ZBESK) (const double&, const double&, const double&,
00063 const int&, const int&, double*, double*,
00064 int&, int&);
00065
00066 F77_RET_T
00067 F77_FUNC (zbesh, ZBESH) (const double&, const double&, const double&,
00068 const int&, const int&, const int&, double*,
00069 double*, int&, int&);
00070
00071 F77_RET_T
00072 F77_FUNC (zairy, ZAIRY) (const double&, const double&, const int&,
00073 const int&, double&, double&, int&, int&);
00074
00075 F77_RET_T
00076 F77_FUNC (zbiry, ZBIRY) (const double&, const double&, const int&,
00077 const int&, double&, double&, int&);
00078
00079 F77_RET_T
00080 F77_FUNC (xdacosh, XDACOSH) (const double&, double&);
00081
00082 F77_RET_T
00083 F77_FUNC (xdasinh, XDASINH) (const double&, double&);
00084
00085 F77_RET_T
00086 F77_FUNC (xdatanh, XDATANH) (const double&, double&);
00087
00088 F77_RET_T
00089 F77_FUNC (xderf, XDERF) (const double&, double&);
00090
00091 F77_RET_T
00092 F77_FUNC (xderfc, XDERFC) (const double&, double&);
00093
00094 F77_RET_T
00095 F77_FUNC (xdbetai, XDBETAI) (const double&, const double&,
00096 const double&, double&);
00097
00098 F77_RET_T
00099 F77_FUNC (xdgamma, XDGAMMA) (const double&, double&);
00100
00101 F77_RET_T
00102 F77_FUNC (xgammainc, XGAMMAINC) (const double&, const double&, double&);
00103
00104 F77_RET_T
00105 F77_FUNC (dlgams, DLGAMS) (const double&, double&, double&);
00106 }
00107
00108 #if !defined (HAVE_ACOSH)
00109 double
00110 acosh (double x)
00111 {
00112 double retval;
00113 F77_FUNC (xdacosh, XDACOSH) (x, retval);
00114 return retval;
00115 }
00116 #endif
00117
00118 #if !defined (HAVE_ASINH)
00119 double
00120 asinh (double x)
00121 {
00122 double retval;
00123 F77_FUNC (xdasinh, XDASINH) (x, retval);
00124 return retval;
00125 }
00126 #endif
00127
00128 #if !defined (HAVE_ATANH)
00129 double
00130 atanh (double x)
00131 {
00132 double retval;
00133 F77_FUNC (xdatanh, XDATANH) (x, retval);
00134 return retval;
00135 }
00136 #endif
00137
00138 #if !defined (HAVE_ERF)
00139 double
00140 erf (double x)
00141 {
00142 double retval;
00143 F77_FUNC (xderf, XDERF) (x, retval);
00144 return retval;
00145 }
00146 #endif
00147
00148 #if !defined (HAVE_ERFC)
00149 double
00150 erfc (double x)
00151 {
00152 double retval;
00153 F77_FUNC (xderfc, XDERFC) (x, retval);
00154 return retval;
00155 }
00156 #endif
00157
00158 double
00159 xgamma (double x)
00160 {
00161 double result;
00162 F77_FUNC (xdgamma, XDGAMMA) (x, result);
00163 return result;
00164 }
00165
00166 double
00167 xlgamma (double x)
00168 {
00169 double result;
00170 double sgngam;
00171
00172 if (x < 0)
00173 (*current_liboctave_error_handler)
00174 ("xlgamma: argument must be nonnegative");
00175
00176 F77_FUNC (dlgams, DLGAMS) (x, result, sgngam);
00177
00178 return result;
00179 }
00180
00181 static inline Complex
00182 zbesj (const Complex& z, double alpha, int kode, int& ierr);
00183
00184 static inline Complex
00185 zbesy (const Complex& z, double alpha, int kode, int& ierr);
00186
00187 static inline Complex
00188 zbesi (const Complex& z, double alpha, int kode, int& ierr);
00189
00190 static inline Complex
00191 zbesk (const Complex& z, double alpha, int kode, int& ierr);
00192
00193 static inline Complex
00194 zbesh1 (const Complex& z, double alpha, int kode, int& ierr);
00195
00196 static inline Complex
00197 zbesh2 (const Complex& z, double alpha, int kode, int& ierr);
00198
00199 static inline Complex
00200 bessel_return_value (const Complex& val, int ierr)
00201 {
00202 static const Complex inf_val = Complex (octave_Inf, octave_Inf);
00203 static const Complex nan_val = Complex (octave_NaN, octave_NaN);
00204
00205 Complex retval;
00206
00207 switch (ierr)
00208 {
00209 case 0:
00210 case 3:
00211 retval = val;
00212 break;
00213
00214 case 2:
00215 retval = inf_val;
00216 break;
00217
00218 default:
00219 retval = nan_val;
00220 break;
00221 }
00222
00223 return retval;
00224 }
00225
00226 static inline bool
00227 is_integer_value (double x)
00228 {
00229 return x == static_cast<double> (static_cast<long> (x));
00230 }
00231
00232 static inline Complex
00233 zbesj (const Complex& z, double alpha, int kode, int& ierr)
00234 {
00235 Complex retval;
00236
00237 if (alpha >= 0.0)
00238 {
00239 double yr = 0.0;
00240 double yi = 0.0;
00241
00242 int nz;
00243
00244 double zr = z.real ();
00245 double zi = z.imag ();
00246
00247 F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
00248
00249 if (kode != 2)
00250 {
00251 double expz = exp (std::abs (zi));
00252 yr *= expz;
00253 yi *= expz;
00254 }
00255
00256 if (zi == 0.0 && zr >= 0.0)
00257 yi = 0.0;
00258
00259 retval = bessel_return_value (Complex (yr, yi), ierr);
00260 }
00261 else if (is_integer_value (alpha))
00262 {
00263
00264 alpha = -alpha;
00265 Complex tmp = zbesj (z, alpha, kode, ierr);
00266 if ((static_cast <long> (alpha)) & 1)
00267 tmp = - tmp;
00268 retval = bessel_return_value (tmp, ierr);
00269 }
00270 else
00271 {
00272 alpha = -alpha;
00273
00274 Complex tmp = cos (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
00275
00276 if (ierr == 0 || ierr == 3)
00277 {
00278 tmp -= sin (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
00279
00280 retval = bessel_return_value (tmp, ierr);
00281 }
00282 else
00283 retval = Complex (octave_NaN, octave_NaN);
00284 }
00285
00286 return retval;
00287 }
00288
00289 static inline Complex
00290 zbesy (const Complex& z, double alpha, int kode, int& ierr)
00291 {
00292 Complex retval;
00293
00294 if (alpha >= 0.0)
00295 {
00296 double yr = 0.0;
00297 double yi = 0.0;
00298
00299 int nz;
00300
00301 double wr, wi;
00302
00303 double zr = z.real ();
00304 double zi = z.imag ();
00305
00306 ierr = 0;
00307
00308 if (zr == 0.0 && zi == 0.0)
00309 {
00310 yr = -octave_Inf;
00311 yi = 0.0;
00312 }
00313 else
00314 {
00315 F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, 2, 1, &yr, &yi, nz,
00316 &wr, &wi, ierr);
00317
00318 if (kode != 2)
00319 {
00320 double expz = exp (std::abs (zi));
00321 yr *= expz;
00322 yi *= expz;
00323 }
00324
00325 if (zi == 0.0 && zr >= 0.0)
00326 yi = 0.0;
00327 }
00328
00329 return bessel_return_value (Complex (yr, yi), ierr);
00330 }
00331 else if (is_integer_value (alpha - 0.5))
00332 {
00333
00334 alpha = -alpha;
00335 Complex tmp = zbesj (z, alpha, kode, ierr);
00336 if ((static_cast <long> (alpha - 0.5)) & 1)
00337 tmp = - tmp;
00338 retval = bessel_return_value (tmp, ierr);
00339 }
00340 else
00341 {
00342 alpha = -alpha;
00343
00344 Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
00345
00346 if (ierr == 0 || ierr == 3)
00347 {
00348 tmp += sin (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
00349
00350 retval = bessel_return_value (tmp, ierr);
00351 }
00352 else
00353 retval = Complex (octave_NaN, octave_NaN);
00354 }
00355
00356 return retval;
00357 }
00358
00359 static inline Complex
00360 zbesi (const Complex& z, double alpha, int kode, int& ierr)
00361 {
00362 Complex retval;
00363
00364 if (alpha >= 0.0)
00365 {
00366 double yr = 0.0;
00367 double yi = 0.0;
00368
00369 int nz;
00370
00371 double zr = z.real ();
00372 double zi = z.imag ();
00373
00374 F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
00375
00376 if (kode != 2)
00377 {
00378 double expz = exp (std::abs (zr));
00379 yr *= expz;
00380 yi *= expz;
00381 }
00382
00383 if (zi == 0.0 && zr >= 0.0)
00384 yi = 0.0;
00385
00386 retval = bessel_return_value (Complex (yr, yi), ierr);
00387 }
00388 else
00389 {
00390 alpha = -alpha;
00391
00392 Complex tmp = zbesi (z, alpha, kode, ierr);
00393
00394 if (ierr == 0 || ierr == 3)
00395 {
00396 if (! is_integer_value (alpha - 0.5))
00397 tmp += (2.0 / M_PI) * sin (M_PI * alpha)
00398 * zbesk (z, alpha, kode, ierr);
00399
00400 retval = bessel_return_value (tmp, ierr);
00401 }
00402 else
00403 retval = Complex (octave_NaN, octave_NaN);
00404 }
00405
00406 return retval;
00407 }
00408
00409 static inline Complex
00410 zbesk (const Complex& z, double alpha, int kode, int& ierr)
00411 {
00412 Complex retval;
00413
00414 if (alpha >= 0.0)
00415 {
00416 double yr = 0.0;
00417 double yi = 0.0;
00418
00419 int nz;
00420
00421 double zr = z.real ();
00422 double zi = z.imag ();
00423
00424 ierr = 0;
00425
00426 if (zr == 0.0 && zi == 0.0)
00427 {
00428 yr = octave_Inf;
00429 yi = 0.0;
00430 }
00431 else
00432 {
00433 F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
00434
00435 if (kode != 2)
00436 {
00437 Complex expz = exp (-z);
00438
00439 double rexpz = real (expz);
00440 double iexpz = imag (expz);
00441
00442 double tmp = yr*rexpz - yi*iexpz;
00443
00444 yi = yr*iexpz + yi*rexpz;
00445 yr = tmp;
00446 }
00447
00448 if (zi == 0.0 && zr >= 0.0)
00449 yi = 0.0;
00450 }
00451
00452 retval = bessel_return_value (Complex (yr, yi), ierr);
00453 }
00454 else
00455 {
00456 Complex tmp = zbesk (z, -alpha, kode, ierr);
00457
00458 retval = bessel_return_value (tmp, ierr);
00459 }
00460
00461 return retval;
00462 }
00463
00464 static inline Complex
00465 zbesh1 (const Complex& z, double alpha, int kode, int& ierr)
00466 {
00467 Complex retval;
00468
00469 if (alpha >= 0.0)
00470 {
00471 double yr = 0.0;
00472 double yi = 0.0;
00473
00474 int nz;
00475
00476 double zr = z.real ();
00477 double zi = z.imag ();
00478
00479 F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz, ierr);
00480
00481 if (kode != 2)
00482 {
00483 Complex expz = exp (Complex (0.0, 1.0) * z);
00484
00485 double rexpz = real (expz);
00486 double iexpz = imag (expz);
00487
00488 double tmp = yr*rexpz - yi*iexpz;
00489
00490 yi = yr*iexpz + yi*rexpz;
00491 yr = tmp;
00492 }
00493
00494 retval = bessel_return_value (Complex (yr, yi), ierr);
00495 }
00496 else
00497 {
00498 alpha = -alpha;
00499
00500 static const Complex eye = Complex (0.0, 1.0);
00501
00502 Complex tmp = exp (M_PI * alpha * eye) * zbesh1 (z, alpha, kode, ierr);
00503
00504 retval = bessel_return_value (tmp, ierr);
00505 }
00506
00507 return retval;
00508 }
00509
00510 static inline Complex
00511 zbesh2 (const Complex& z, double alpha, int kode, int& ierr)
00512 {
00513 Complex retval;
00514
00515 if (alpha >= 0.0)
00516 {
00517 double yr = 0.0;
00518 double yi = 0.0;
00519
00520 int nz;
00521
00522 double zr = z.real ();
00523 double zi = z.imag ();
00524
00525 F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz, ierr);
00526
00527 if (kode != 2)
00528 {
00529 Complex expz = exp (-Complex (0.0, 1.0) * z);
00530
00531 double rexpz = real (expz);
00532 double iexpz = imag (expz);
00533
00534 double tmp = yr*rexpz - yi*iexpz;
00535
00536 yi = yr*iexpz + yi*rexpz;
00537 yr = tmp;
00538 }
00539
00540 retval = bessel_return_value (Complex (yr, yi), ierr);
00541 }
00542 else
00543 {
00544 alpha = -alpha;
00545
00546 static const Complex eye = Complex (0.0, 1.0);
00547
00548 Complex tmp = exp (-M_PI * alpha * eye) * zbesh2 (z, alpha, kode, ierr);
00549
00550 retval = bessel_return_value (tmp, ierr);
00551 }
00552
00553 return retval;
00554 }
00555
00556 typedef Complex (*fptr) (const Complex&, double, int, int&);
00557
00558 static inline Complex
00559 do_bessel (fptr f, const char *, double alpha, const Complex& x,
00560 bool scaled, int& ierr)
00561 {
00562 Complex retval;
00563
00564 retval = f (x, alpha, (scaled ? 2 : 1), ierr);
00565
00566 return retval;
00567 }
00568
00569 static inline ComplexMatrix
00570 do_bessel (fptr f, const char *, double alpha, const ComplexMatrix& x,
00571 bool scaled, Array2<int>& ierr)
00572 {
00573 int nr = x.rows ();
00574 int nc = x.cols ();
00575
00576 ComplexMatrix retval (nr, nc);
00577
00578 ierr.resize (nr, nc);
00579
00580 for (int j = 0; j < nc; j++)
00581 for (int i = 0; i < nr; i++)
00582 retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j));
00583
00584 return retval;
00585 }
00586
00587 static inline ComplexMatrix
00588 do_bessel (fptr f, const char *, const Matrix& alpha, const Complex& x,
00589 bool scaled, Array2<int>& ierr)
00590 {
00591 int nr = alpha.rows ();
00592 int nc = alpha.cols ();
00593
00594 ComplexMatrix retval (nr, nc);
00595
00596 ierr.resize (nr, nc);
00597
00598 for (int j = 0; j < nc; j++)
00599 for (int i = 0; i < nr; i++)
00600 retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
00601
00602 return retval;
00603 }
00604
00605 static inline ComplexMatrix
00606 do_bessel (fptr f, const char *fn, const Matrix& alpha,
00607 const ComplexMatrix& x, bool scaled, Array2<int>& ierr)
00608 {
00609 ComplexMatrix retval;
00610
00611 int x_nr = x.rows ();
00612 int x_nc = x.cols ();
00613
00614 int alpha_nr = alpha.rows ();
00615 int alpha_nc = alpha.cols ();
00616
00617 if (x_nr == alpha_nr && x_nc == alpha_nc)
00618 {
00619 int nr = x_nr;
00620 int nc = x_nc;
00621
00622 retval.resize (nr, nc);
00623
00624 ierr.resize (nr, nc);
00625
00626 for (int j = 0; j < nc; j++)
00627 for (int i = 0; i < nr; i++)
00628 retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
00629 }
00630 else
00631 (*current_liboctave_error_handler)
00632 ("%s: the sizes of alpha and x must conform", fn);
00633
00634 return retval;
00635 }
00636
00637 static inline ComplexNDArray
00638 do_bessel (fptr f, const char *, double alpha, const ComplexNDArray& x,
00639 bool scaled, ArrayN<int>& ierr)
00640 {
00641 dim_vector dv = x.dims ();
00642 int nel = dv.numel ();
00643 ComplexNDArray retval (dv);
00644
00645 ierr.resize (dv);
00646
00647 for (int i = 0; i < nel; i++)
00648 retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i));
00649
00650 return retval;
00651 }
00652
00653 static inline ComplexNDArray
00654 do_bessel (fptr f, const char *, const NDArray& alpha, const Complex& x,
00655 bool scaled, ArrayN<int>& ierr)
00656 {
00657 dim_vector dv = alpha.dims ();
00658 int nel = dv.numel ();
00659 ComplexNDArray retval (dv);
00660
00661 ierr.resize (dv);
00662
00663 for (int i = 0; i < nel; i++)
00664 retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i));
00665
00666 return retval;
00667 }
00668
00669 static inline ComplexNDArray
00670 do_bessel (fptr f, const char *fn, const NDArray& alpha,
00671 const ComplexNDArray& x, bool scaled, ArrayN<int>& ierr)
00672 {
00673 dim_vector dv = x.dims ();
00674 ComplexNDArray retval;
00675
00676 if (dv == alpha.dims ())
00677 {
00678 int nel = dv.numel ();
00679
00680 retval.resize (dv);
00681 ierr.resize (dv);
00682
00683 for (int i = 0; i < nel; i++)
00684 retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i));
00685 }
00686 else
00687 (*current_liboctave_error_handler)
00688 ("%s: the sizes of alpha and x must conform", fn);
00689
00690 return retval;
00691 }
00692
00693 static inline ComplexMatrix
00694 do_bessel (fptr f, const char *, const RowVector& alpha,
00695 const ComplexColumnVector& x, bool scaled, Array2<int>& ierr)
00696 {
00697 int nr = x.length ();
00698 int nc = alpha.length ();
00699
00700 ComplexMatrix retval (nr, nc);
00701
00702 ierr.resize (nr, nc);
00703
00704 for (int j = 0; j < nc; j++)
00705 for (int i = 0; i < nr; i++)
00706 retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j));
00707
00708 return retval;
00709 }
00710
00711 #define SS_BESSEL(name, fcn) \
00712 Complex \
00713 name (double alpha, const Complex& x, bool scaled, int& ierr) \
00714 { \
00715 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
00716 }
00717
00718 #define SM_BESSEL(name, fcn) \
00719 ComplexMatrix \
00720 name (double alpha, const ComplexMatrix& x, bool scaled, \
00721 Array2<int>& ierr) \
00722 { \
00723 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
00724 }
00725
00726 #define MS_BESSEL(name, fcn) \
00727 ComplexMatrix \
00728 name (const Matrix& alpha, const Complex& x, bool scaled, \
00729 Array2<int>& ierr) \
00730 { \
00731 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
00732 }
00733
00734 #define MM_BESSEL(name, fcn) \
00735 ComplexMatrix \
00736 name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \
00737 Array2<int>& ierr) \
00738 { \
00739 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
00740 }
00741
00742 #define SN_BESSEL(name, fcn) \
00743 ComplexNDArray \
00744 name (double alpha, const ComplexNDArray& x, bool scaled, \
00745 ArrayN<int>& ierr) \
00746 { \
00747 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
00748 }
00749
00750 #define NS_BESSEL(name, fcn) \
00751 ComplexNDArray \
00752 name (const NDArray& alpha, const Complex& x, bool scaled, \
00753 ArrayN<int>& ierr) \
00754 { \
00755 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
00756 }
00757
00758 #define NN_BESSEL(name, fcn) \
00759 ComplexNDArray \
00760 name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \
00761 ArrayN<int>& ierr) \
00762 { \
00763 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
00764 }
00765
00766 #define RC_BESSEL(name, fcn) \
00767 ComplexMatrix \
00768 name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
00769 Array2<int>& ierr) \
00770 { \
00771 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
00772 }
00773
00774 #define ALL_BESSEL(name, fcn) \
00775 SS_BESSEL (name, fcn) \
00776 SM_BESSEL (name, fcn) \
00777 MS_BESSEL (name, fcn) \
00778 MM_BESSEL (name, fcn) \
00779 SN_BESSEL (name, fcn) \
00780 NS_BESSEL (name, fcn) \
00781 NN_BESSEL (name, fcn) \
00782 RC_BESSEL (name, fcn)
00783
00784 ALL_BESSEL (besselj, zbesj)
00785 ALL_BESSEL (bessely, zbesy)
00786 ALL_BESSEL (besseli, zbesi)
00787 ALL_BESSEL (besselk, zbesk)
00788 ALL_BESSEL (besselh1, zbesh1)
00789 ALL_BESSEL (besselh2, zbesh2)
00790
00791 Complex
00792 airy (const Complex& z, bool deriv, bool scaled, int& ierr)
00793 {
00794 double ar = 0.0;
00795 double ai = 0.0;
00796
00797 int nz;
00798
00799 double zr = z.real ();
00800 double zi = z.imag ();
00801
00802 int id = deriv ? 1 : 0;
00803
00804 F77_FUNC (zairy, ZAIRY) (zr, zi, id, 2, ar, ai, nz, ierr);
00805
00806 if (! scaled)
00807 {
00808 Complex expz = exp (- 2.0 / 3.0 * z * sqrt(z));
00809
00810 double rexpz = real (expz);
00811 double iexpz = imag (expz);
00812
00813 double tmp = ar*rexpz - ai*iexpz;
00814
00815 ai = ar*iexpz + ai*rexpz;
00816 ar = tmp;
00817 }
00818
00819 if (zi == 0.0 && (! scaled || zr >= 0.0))
00820 ai = 0.0;
00821
00822 return bessel_return_value (Complex (ar, ai), ierr);
00823 }
00824
00825 Complex
00826 biry (const Complex& z, bool deriv, bool scaled, int& ierr)
00827 {
00828 double ar = 0.0;
00829 double ai = 0.0;
00830
00831 double zr = z.real ();
00832 double zi = z.imag ();
00833
00834 int id = deriv ? 1 : 0;
00835
00836 F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 2, ar, ai, ierr);
00837
00838 if (! scaled)
00839 {
00840 Complex expz = exp (std::abs (real (2.0 / 3.0 * z * sqrt (z))));
00841
00842 double rexpz = real (expz);
00843 double iexpz = imag (expz);
00844
00845 double tmp = ar*rexpz - ai*iexpz;
00846
00847 ai = ar*iexpz + ai*rexpz;
00848 ar = tmp;
00849 }
00850
00851 if (zi == 0.0 && (! scaled || zr >= 0.0))
00852 ai = 0.0;
00853
00854 return bessel_return_value (Complex (ar, ai), ierr);
00855 }
00856
00857 ComplexMatrix
00858 airy (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr)
00859 {
00860 int nr = z.rows ();
00861 int nc = z.cols ();
00862
00863 ComplexMatrix retval (nr, nc);
00864
00865 ierr.resize (nr, nc);
00866
00867 for (int j = 0; j < nc; j++)
00868 for (int i = 0; i < nr; i++)
00869 retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j));
00870
00871 return retval;
00872 }
00873
00874 ComplexMatrix
00875 biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<int>& ierr)
00876 {
00877 int nr = z.rows ();
00878 int nc = z.cols ();
00879
00880 ComplexMatrix retval (nr, nc);
00881
00882 ierr.resize (nr, nc);
00883
00884 for (int j = 0; j < nc; j++)
00885 for (int i = 0; i < nr; i++)
00886 retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j));
00887
00888 return retval;
00889 }
00890
00891 ComplexNDArray
00892 airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr)
00893 {
00894 dim_vector dv = z.dims ();
00895 int nel = dv.numel ();
00896 ComplexNDArray retval (dv);
00897
00898 ierr.resize (dv);
00899
00900 for (int i = 0; i < nel; i++)
00901 retval (i) = airy (z(i), deriv, scaled, ierr(i));
00902
00903 return retval;
00904 }
00905
00906 ComplexNDArray
00907 biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<int>& ierr)
00908 {
00909 dim_vector dv = z.dims ();
00910 int nel = dv.numel ();
00911 ComplexNDArray retval (dv);
00912
00913 ierr.resize (dv);
00914
00915 for (int i = 0; i < nel; i++)
00916 retval (i) = biry (z(i), deriv, scaled, ierr(i));
00917
00918 return retval;
00919 }
00920
00921 static void
00922 gripe_betainc_nonconformant (int r1, int c1, int r2, int c2, int r3,
00923 int c3)
00924 {
00925 (*current_liboctave_error_handler)
00926 ("betainc: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)",
00927 r1, c1, r2, c2, r3, c3);
00928 }
00929
00930 static dim_vector null_dims (0);
00931
00932 static void
00933 gripe_betainc_nonconformant (const dim_vector& d1, const dim_vector& d2,
00934 const dim_vector& d3)
00935 {
00936 std::string d1_str = d1.str ();
00937 std::string d2_str = d2.str ();
00938 std::string d3_str = d3.str ();
00939
00940 (*current_liboctave_error_handler)
00941 ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)",
00942 d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
00943 }
00944
00945 double
00946 betainc (double x, double a, double b)
00947 {
00948 double retval;
00949 F77_FUNC (xdbetai, XDBETAI) (x, a, b, retval);
00950 return retval;
00951 }
00952
00953 Matrix
00954 betainc (double x, double a, const Matrix& b)
00955 {
00956 int nr = b.rows ();
00957 int nc = b.cols ();
00958
00959 Matrix retval (nr, nc);
00960
00961 for (int j = 0; j < nc; j++)
00962 for (int i = 0; i < nr; i++)
00963 retval(i,j) = betainc (x, a, b(i,j));
00964
00965 return retval;
00966 }
00967
00968 Matrix
00969 betainc (double x, const Matrix& a, double b)
00970 {
00971 int nr = a.rows ();
00972 int nc = a.cols ();
00973
00974 Matrix retval (nr, nc);
00975
00976 for (int j = 0; j < nc; j++)
00977 for (int i = 0; i < nr; i++)
00978 retval(i,j) = betainc (x, a(i,j), b);
00979
00980 return retval;
00981 }
00982
00983 Matrix
00984 betainc (double x, const Matrix& a, const Matrix& b)
00985 {
00986 Matrix retval;
00987
00988 int a_nr = a.rows ();
00989 int a_nc = a.cols ();
00990
00991 int b_nr = b.rows ();
00992 int b_nc = b.cols ();
00993
00994 if (a_nr == b_nr && a_nc == b_nc)
00995 {
00996 retval.resize (a_nr, a_nc);
00997
00998 for (int j = 0; j < a_nc; j++)
00999 for (int i = 0; i < a_nr; i++)
01000 retval(i,j) = betainc (x, a(i,j), b(i,j));
01001 }
01002 else
01003 gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc);
01004
01005 return retval;
01006 }
01007
01008 NDArray
01009 betainc (double x, double a, const NDArray& b)
01010 {
01011 dim_vector dv = b.dims ();
01012 int nel = dv.numel ();
01013
01014 NDArray retval (dv);
01015
01016 for (int i = 0; i < nel; i++)
01017 retval (i) = betainc (x, a, b(i));
01018
01019 return retval;
01020 }
01021
01022 NDArray
01023 betainc (double x, const NDArray& a, double b)
01024 {
01025 dim_vector dv = a.dims ();
01026 int nel = dv.numel ();
01027
01028 NDArray retval (dv);
01029
01030 for (int i = 0; i < nel; i++)
01031 retval (i) = betainc (x, a(i), b);
01032
01033 return retval;
01034 }
01035
01036 NDArray
01037 betainc (double x, const NDArray& a, const NDArray& b)
01038 {
01039 NDArray retval;
01040 dim_vector dv = a.dims ();
01041
01042 if (dv == b.dims ())
01043 {
01044 int nel = dv.numel ();
01045
01046 retval.resize (dv);
01047
01048 for (int i = 0; i < nel; i++)
01049 retval (i) = betainc (x, a(i), b(i));
01050 }
01051 else
01052 gripe_betainc_nonconformant (dim_vector (0), dv, b.dims ());
01053
01054 return retval;
01055 }
01056
01057
01058 Matrix
01059 betainc (const Matrix& x, double a, double b)
01060 {
01061 int nr = x.rows ();
01062 int nc = x.cols ();
01063
01064 Matrix retval (nr, nc);
01065
01066 for (int j = 0; j < nc; j++)
01067 for (int i = 0; i < nr; i++)
01068 retval(i,j) = betainc (x(i,j), a, b);
01069
01070 return retval;
01071 }
01072
01073 Matrix
01074 betainc (const Matrix& x, double a, const Matrix& b)
01075 {
01076 Matrix retval;
01077
01078 int nr = x.rows ();
01079 int nc = x.cols ();
01080
01081 int b_nr = b.rows ();
01082 int b_nc = b.cols ();
01083
01084 if (nr == b_nr && nc == b_nc)
01085 {
01086 retval.resize (nr, nc);
01087
01088 for (int j = 0; j < nc; j++)
01089 for (int i = 0; i < nr; i++)
01090 retval(i,j) = betainc (x(i,j), a, b(i,j));
01091 }
01092 else
01093 gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc);
01094
01095 return retval;
01096 }
01097
01098 Matrix
01099 betainc (const Matrix& x, const Matrix& a, double b)
01100 {
01101 Matrix retval;
01102
01103 int nr = x.rows ();
01104 int nc = x.cols ();
01105
01106 int a_nr = a.rows ();
01107 int a_nc = a.cols ();
01108
01109 if (nr == a_nr && nc == a_nc)
01110 {
01111 retval.resize (nr, nc);
01112
01113 for (int j = 0; j < nc; j++)
01114 for (int i = 0; i < nr; i++)
01115 retval(i,j) = betainc (x(i,j), a(i,j), b);
01116 }
01117 else
01118 gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1);
01119
01120 return retval;
01121 }
01122
01123 Matrix
01124 betainc (const Matrix& x, const Matrix& a, const Matrix& b)
01125 {
01126 Matrix retval;
01127
01128 int nr = x.rows ();
01129 int nc = x.cols ();
01130
01131 int a_nr = a.rows ();
01132 int a_nc = a.cols ();
01133
01134 int b_nr = b.rows ();
01135 int b_nc = b.cols ();
01136
01137 if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc)
01138 {
01139 retval.resize (nr, nc);
01140
01141 for (int j = 0; j < nc; j++)
01142 for (int i = 0; i < nr; i++)
01143 retval(i,j) = betainc (x(i,j), a(i,j), b(i,j));
01144 }
01145 else
01146 gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc);
01147
01148 return retval;
01149 }
01150
01151 NDArray
01152 betainc (const NDArray& x, double a, double b)
01153 {
01154 dim_vector dv = x.dims ();
01155 int nel = dv.numel ();
01156
01157 NDArray retval (dv);
01158
01159 for (int i = 0; i < nel; i++)
01160 retval (i) = betainc (x(i), a, b);
01161
01162 return retval;
01163 }
01164
01165 NDArray
01166 betainc (const NDArray& x, double a, const NDArray& b)
01167 {
01168 NDArray retval;
01169 dim_vector dv = x.dims ();
01170
01171 if (dv == b.dims ())
01172 {
01173 int nel = dv.numel ();
01174
01175 retval.resize (dv);
01176
01177 for (int i = 0; i < nel; i++)
01178 retval (i) = betainc (x(i), a, b(i));
01179 }
01180 else
01181 gripe_betainc_nonconformant (dv, dim_vector (0), b.dims ());
01182
01183 return retval;
01184 }
01185
01186 NDArray
01187 betainc (const NDArray& x, const NDArray& a, double b)
01188 {
01189 NDArray retval;
01190 dim_vector dv = x.dims ();
01191
01192 if (dv == a.dims ())
01193 {
01194 int nel = dv.numel ();
01195
01196 retval.resize (dv);
01197
01198 for (int i = 0; i < nel; i++)
01199 retval (i) = betainc (x(i), a(i), b);
01200 }
01201 else
01202 gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0));
01203
01204 return retval;
01205 }
01206
01207 NDArray
01208 betainc (const NDArray& x, const NDArray& a, const NDArray& b)
01209 {
01210 NDArray retval;
01211 dim_vector dv = x.dims ();
01212
01213 if (dv == a.dims () && dv == b.dims ())
01214 {
01215 int nel = dv.numel ();
01216
01217 retval.resize (dv);
01218
01219 for (int i = 0; i < nel; i++)
01220 retval (i) = betainc (x(i), a(i), b(i));
01221 }
01222 else
01223 gripe_betainc_nonconformant (dv, a.dims (), b.dims ());
01224
01225 return retval;
01226 }
01227
01228
01229
01230 double
01231 gammainc (double x, double a, bool& err)
01232 {
01233 double retval;
01234
01235 err = false;
01236
01237 if (a < 0.0 || x < 0.0)
01238 {
01239 (*current_liboctave_error_handler)
01240 ("gammainc: A and X must be non-negative");
01241
01242 err = true;
01243 }
01244 else
01245 F77_FUNC (xgammainc, XGAMMAINC) (a, x, retval);
01246
01247 return retval;
01248 }
01249
01250 Matrix
01251 gammainc (double x, const Matrix& a)
01252 {
01253 int nr = a.rows ();
01254 int nc = a.cols ();
01255
01256 Matrix result (nr, nc);
01257 Matrix retval;
01258
01259 bool err;
01260
01261 for (int j = 0; j < nc; j++)
01262 for (int i = 0; i < nr; i++)
01263 {
01264 result(i,j) = gammainc (x, a(i,j), err);
01265
01266 if (err)
01267 goto done;
01268 }
01269
01270 retval = result;
01271
01272 done:
01273
01274 return retval;
01275 }
01276
01277 Matrix
01278 gammainc (const Matrix& x, double a)
01279 {
01280 int nr = x.rows ();
01281 int nc = x.cols ();
01282
01283 Matrix result (nr, nc);
01284 Matrix retval;
01285
01286 bool err;
01287
01288 for (int j = 0; j < nc; j++)
01289 for (int i = 0; i < nr; i++)
01290 {
01291 result(i,j) = gammainc (x(i,j), a, err);
01292
01293 if (err)
01294 goto done;
01295 }
01296
01297 retval = result;
01298
01299 done:
01300
01301 return retval;
01302 }
01303
01304 Matrix
01305 gammainc (const Matrix& x, const Matrix& a)
01306 {
01307 Matrix result;
01308 Matrix retval;
01309
01310 int nr = x.rows ();
01311 int nc = x.cols ();
01312
01313 int a_nr = a.rows ();
01314 int a_nc = a.cols ();
01315
01316 if (nr == a_nr && nc == a_nc)
01317 {
01318 result.resize (nr, nc);
01319
01320 bool err;
01321
01322 for (int j = 0; j < nc; j++)
01323 for (int i = 0; i < nr; i++)
01324 {
01325 result(i,j) = gammainc (x(i,j), a(i,j), err);
01326
01327 if (err)
01328 goto done;
01329 }
01330
01331 retval = result;
01332 }
01333 else
01334 (*current_liboctave_error_handler)
01335 ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
01336 nr, nc, a_nr, a_nc);
01337
01338 done:
01339
01340 return retval;
01341 }
01342
01343 NDArray
01344 gammainc (double x, const NDArray& a)
01345 {
01346 dim_vector dv = a.dims ();
01347 int nel = dv.numel ();
01348
01349 NDArray retval;
01350 NDArray result (dv);
01351
01352 bool err;
01353
01354 for (int i = 0; i < nel; i++)
01355 {
01356 result (i) = gammainc (x, a(i), err);
01357
01358 if (err)
01359 goto done;
01360 }
01361
01362 retval = result;
01363
01364 done:
01365
01366 return retval;
01367 }
01368
01369 NDArray
01370 gammainc (const NDArray& x, double a)
01371 {
01372 dim_vector dv = x.dims ();
01373 int nel = dv.numel ();
01374
01375 NDArray retval;
01376 NDArray result (dv);
01377
01378 bool err;
01379
01380 for (int i = 0; i < nel; i++)
01381 {
01382 result (i) = gammainc (x(i), a, err);
01383
01384 if (err)
01385 goto done;
01386 }
01387
01388 retval = result;
01389
01390 done:
01391
01392 return retval;
01393 }
01394
01395 NDArray
01396 gammainc (const NDArray& x, const NDArray& a)
01397 {
01398 dim_vector dv = x.dims ();
01399 int nel = dv.numel ();
01400
01401 NDArray retval;
01402 NDArray result;
01403
01404 if (dv == a.dims ())
01405 {
01406 result.resize (dv);
01407
01408 bool err;
01409
01410 for (int i = 0; i < nel; i++)
01411 {
01412 result (i) = gammainc (x(i), a(i), err);
01413
01414 if (err)
01415 goto done;
01416 }
01417
01418 retval = result;
01419 }
01420 else
01421 {
01422 std::string x_str = dv.str ();
01423 std::string a_str = a.dims ().str ();
01424
01425 (*current_liboctave_error_handler)
01426 ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
01427 x_str.c_str (), a_str. c_str ());
01428 }
01429
01430 done:
01431
01432 return retval;
01433 }
01434
01435
01436
01437
01438
01439