Mercurial > octave-nkf
annotate liboctave/lo-specfun.cc @ 11904:059fadc0cbc3 release-3-0-x
fix scaling factor for negative alpha in zbesi,cbesi
add bessel function tests
add all scaling factors to bessel documentation
| author | Brian Gough <bjg@gnu.org> |
|---|---|
| date | Mon, 12 Jan 2009 10:42:03 +0100 |
| parents | 6525eb2fba0f |
| children |
| rev | line source |
|---|---|
| 3146 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1996, 1998, 2002, 2003, 2004, 2005, 2006, 2007 |
| 4 John W. Eaton | |
| 3146 | 5 |
| 6 This file is part of Octave. | |
| 7 | |
| 8 Octave is free software; you can redistribute it and/or modify it | |
| 9 under the terms of the GNU General Public License as published by the | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 3146 | 12 |
| 13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 16 for more details. | |
| 17 | |
| 18 You should have received a copy of the GNU General Public License | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 3146 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include "Range.h" | |
| 3220 | 29 #include "CColVector.h" |
| 30 #include "CMatrix.h" | |
| 31 #include "dRowVector.h" | |
| 3146 | 32 #include "dMatrix.h" |
| 4844 | 33 #include "dNDArray.h" |
| 34 #include "CNDArray.h" | |
| 3146 | 35 #include "f77-fcn.h" |
| 36 #include "lo-error.h" | |
| 3220 | 37 #include "lo-ieee.h" |
| 38 #include "lo-specfun.h" | |
| 3146 | 39 #include "mx-inlines.cc" |
| 5701 | 40 #include "lo-mappers.h" |
| 3146 | 41 |
| 4064 | 42 #ifndef M_PI |
| 43 #define M_PI 3.14159265358979323846 | |
| 44 #endif | |
| 45 | |
| 3146 | 46 extern "C" |
| 47 { | |
| 4552 | 48 F77_RET_T |
| 49 F77_FUNC (zbesj, ZBESJ) (const double&, const double&, const double&, | |
| 5275 | 50 const octave_idx_type&, const octave_idx_type&, double*, double*, |
| 51 octave_idx_type&, octave_idx_type&); | |
| 3146 | 52 |
| 4552 | 53 F77_RET_T |
| 54 F77_FUNC (zbesy, ZBESY) (const double&, const double&, const double&, | |
| 5275 | 55 const octave_idx_type&, const octave_idx_type&, double*, double*, |
| 56 octave_idx_type&, double*, double*, octave_idx_type&); | |
| 3220 | 57 |
| 4552 | 58 F77_RET_T |
| 59 F77_FUNC (zbesi, ZBESI) (const double&, const double&, const double&, | |
| 5275 | 60 const octave_idx_type&, const octave_idx_type&, double*, double*, |
| 61 octave_idx_type&, octave_idx_type&); | |
| 3146 | 62 |
| 4552 | 63 F77_RET_T |
| 64 F77_FUNC (zbesk, ZBESK) (const double&, const double&, const double&, | |
| 5275 | 65 const octave_idx_type&, const octave_idx_type&, double*, double*, |
| 66 octave_idx_type&, octave_idx_type&); | |
| 3220 | 67 |
| 4552 | 68 F77_RET_T |
| 69 F77_FUNC (zbesh, ZBESH) (const double&, const double&, const double&, | |
| 5275 | 70 const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, double*, |
| 71 double*, octave_idx_type&, octave_idx_type&); | |
| 4552 | 72 |
| 73 F77_RET_T | |
| 5275 | 74 F77_FUNC (zairy, ZAIRY) (const double&, const double&, const octave_idx_type&, |
| 75 const octave_idx_type&, double&, double&, octave_idx_type&, octave_idx_type&); | |
| 3146 | 76 |
| 4552 | 77 F77_RET_T |
| 5275 | 78 F77_FUNC (zbiry, ZBIRY) (const double&, const double&, const octave_idx_type&, |
| 79 const octave_idx_type&, double&, double&, octave_idx_type&); | |
| 4552 | 80 |
| 81 F77_RET_T | |
| 82 F77_FUNC (xdacosh, XDACOSH) (const double&, double&); | |
| 3220 | 83 |
| 4552 | 84 F77_RET_T |
| 85 F77_FUNC (xdasinh, XDASINH) (const double&, double&); | |
| 3146 | 86 |
| 4552 | 87 F77_RET_T |
| 88 F77_FUNC (xdatanh, XDATANH) (const double&, double&); | |
| 3146 | 89 |
| 4552 | 90 F77_RET_T |
| 91 F77_FUNC (xderf, XDERF) (const double&, double&); | |
| 3146 | 92 |
| 4552 | 93 F77_RET_T |
| 94 F77_FUNC (xderfc, XDERFC) (const double&, double&); | |
| 3146 | 95 |
| 4552 | 96 F77_RET_T |
| 97 F77_FUNC (xdbetai, XDBETAI) (const double&, const double&, | |
| 98 const double&, double&); | |
| 3146 | 99 |
| 4552 | 100 F77_RET_T |
| 101 F77_FUNC (xdgamma, XDGAMMA) (const double&, double&); | |
| 3146 | 102 |
| 4552 | 103 F77_RET_T |
| 104 F77_FUNC (xgammainc, XGAMMAINC) (const double&, const double&, double&); | |
| 3146 | 105 |
| 4552 | 106 F77_RET_T |
| 107 F77_FUNC (dlgams, DLGAMS) (const double&, double&, double&); | |
| 3146 | 108 } |
| 109 | |
| 110 #if !defined (HAVE_ACOSH) | |
| 111 double | |
| 112 acosh (double x) | |
| 113 { | |
| 114 double retval; | |
| 5278 | 115 F77_XFCN (xdacosh, XDACOSH, (x, retval)); |
| 3146 | 116 return retval; |
| 117 } | |
| 118 #endif | |
| 119 | |
| 120 #if !defined (HAVE_ASINH) | |
| 121 double | |
| 122 asinh (double x) | |
| 123 { | |
| 124 double retval; | |
| 5278 | 125 F77_XFCN (xdasinh, XDASINH, (x, retval)); |
| 3146 | 126 return retval; |
| 127 } | |
| 128 #endif | |
| 129 | |
| 130 #if !defined (HAVE_ATANH) | |
| 131 double | |
| 132 atanh (double x) | |
| 133 { | |
| 134 double retval; | |
| 5278 | 135 F77_XFCN (xdatanh, XDATANH, (x, retval)); |
| 3146 | 136 return retval; |
| 137 } | |
| 138 #endif | |
| 139 | |
| 140 #if !defined (HAVE_ERF) | |
| 141 double | |
| 142 erf (double x) | |
| 143 { | |
| 144 double retval; | |
| 5278 | 145 F77_XFCN (xderf, XDERF, (x, retval)); |
| 3146 | 146 return retval; |
| 147 } | |
| 148 #endif | |
| 149 | |
| 150 #if !defined (HAVE_ERFC) | |
| 151 double | |
| 152 erfc (double x) | |
| 153 { | |
| 154 double retval; | |
| 5278 | 155 F77_XFCN (xderfc, XDERFC, (x, retval)); |
| 3146 | 156 return retval; |
| 157 } | |
| 158 #endif | |
| 159 | |
| 160 double | |
| 3156 | 161 xgamma (double x) |
| 3146 | 162 { |
| 6969 | 163 #if defined (HAVE_TGAMMA) |
| 164 return tgamma (x); | |
| 165 #else | |
| 3156 | 166 double result; |
| 5701 | 167 |
| 168 if (xisnan (x)) | |
| 169 result = x; | |
| 170 else if ((x <= 0 && D_NINT (x) == x) || xisinf (x)) | |
| 171 result = octave_Inf; | |
| 172 else | |
| 173 F77_XFCN (xdgamma, XDGAMMA, (x, result)); | |
| 6969 | 174 |
| 3156 | 175 return result; |
| 6969 | 176 #endif |
| 3146 | 177 } |
| 178 | |
| 179 double | |
| 3156 | 180 xlgamma (double x) |
| 3146 | 181 { |
| 6969 | 182 #if defined (HAVE_LGAMMA) |
| 183 return lgamma (x); | |
| 184 #else | |
| 3156 | 185 double result; |
| 3146 | 186 double sgngam; |
| 4497 | 187 |
| 5701 | 188 if (xisnan (x)) |
| 189 result = x; | |
| 6969 | 190 else if (xisinf (x)) |
| 5701 | 191 result = octave_Inf; |
| 5700 | 192 else |
| 193 F77_XFCN (dlgams, DLGAMS, (x, result, sgngam)); | |
| 4497 | 194 |
| 3156 | 195 return result; |
| 6969 | 196 #endif |
| 6961 | 197 } |
| 198 | |
| 3220 | 199 static inline Complex |
| 5275 | 200 zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr); |
| 3220 | 201 |
| 202 static inline Complex | |
| 5275 | 203 zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr); |
| 3220 | 204 |
| 205 static inline Complex | |
| 5275 | 206 zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr); |
| 3220 | 207 |
| 208 static inline Complex | |
| 5275 | 209 zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr); |
| 3220 | 210 |
| 211 static inline Complex | |
| 5275 | 212 zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr); |
| 3220 | 213 |
| 214 static inline Complex | |
| 5275 | 215 zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr); |
| 3220 | 216 |
| 217 static inline Complex | |
| 5275 | 218 bessel_return_value (const Complex& val, octave_idx_type ierr) |
| 3146 | 219 { |
| 3220 | 220 static const Complex inf_val = Complex (octave_Inf, octave_Inf); |
| 221 static const Complex nan_val = Complex (octave_NaN, octave_NaN); | |
| 222 | |
| 223 Complex retval; | |
| 224 | |
| 225 switch (ierr) | |
| 226 { | |
| 227 case 0: | |
| 228 case 3: | |
| 229 retval = val; | |
| 230 break; | |
| 231 | |
| 232 case 2: | |
| 233 retval = inf_val; | |
| 234 break; | |
| 235 | |
| 236 default: | |
| 237 retval = nan_val; | |
| 238 break; | |
| 239 } | |
| 240 | |
| 3146 | 241 return retval; |
| 242 } | |
| 243 | |
| 4911 | 244 static inline bool |
| 245 is_integer_value (double x) | |
| 246 { | |
| 247 return x == static_cast<double> (static_cast<long> (x)); | |
| 248 } | |
| 249 | |
| 3220 | 250 static inline Complex |
| 5275 | 251 zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr) |
| 3146 | 252 { |
| 3220 | 253 Complex retval; |
| 254 | |
| 255 if (alpha >= 0.0) | |
| 256 { | |
| 257 double yr = 0.0; | |
| 258 double yi = 0.0; | |
| 259 | |
| 5275 | 260 octave_idx_type nz; |
| 3220 | 261 |
| 262 double zr = z.real (); | |
| 263 double zi = z.imag (); | |
| 264 | |
| 4506 | 265 F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr); |
| 266 | |
| 267 if (kode != 2) | |
| 268 { | |
| 269 double expz = exp (std::abs (zi)); | |
| 270 yr *= expz; | |
| 271 yi *= expz; | |
| 272 } | |
| 3220 | 273 |
| 4490 | 274 if (zi == 0.0 && zr >= 0.0) |
| 3220 | 275 yi = 0.0; |
| 276 | |
| 277 retval = bessel_return_value (Complex (yr, yi), ierr); | |
| 278 } | |
| 4911 | 279 else if (is_integer_value (alpha)) |
| 280 { | |
| 281 // zbesy can overflow as z->0, and cause troubles for generic case below | |
| 282 alpha = -alpha; | |
| 283 Complex tmp = zbesj (z, alpha, kode, ierr); | |
| 284 if ((static_cast <long> (alpha)) & 1) | |
| 285 tmp = - tmp; | |
| 286 retval = bessel_return_value (tmp, ierr); | |
| 287 } | |
| 3220 | 288 else |
| 289 { | |
| 290 alpha = -alpha; | |
| 291 | |
| 292 Complex tmp = cos (M_PI * alpha) * zbesj (z, alpha, kode, ierr); | |
| 293 | |
| 294 if (ierr == 0 || ierr == 3) | |
| 295 { | |
| 296 tmp -= sin (M_PI * alpha) * zbesy (z, alpha, kode, ierr); | |
| 297 | |
| 298 retval = bessel_return_value (tmp, ierr); | |
| 299 } | |
| 300 else | |
| 301 retval = Complex (octave_NaN, octave_NaN); | |
| 302 } | |
| 303 | |
| 3146 | 304 return retval; |
| 305 } | |
| 306 | |
| 3220 | 307 static inline Complex |
| 5275 | 308 zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr) |
| 3146 | 309 { |
| 3220 | 310 Complex retval; |
| 3146 | 311 |
| 312 if (alpha >= 0.0) | |
| 313 { | |
| 3220 | 314 double yr = 0.0; |
| 315 double yi = 0.0; | |
| 316 | |
| 5275 | 317 octave_idx_type nz; |
| 3220 | 318 |
| 319 double wr, wi; | |
| 3146 | 320 |
| 3220 | 321 double zr = z.real (); |
| 322 double zi = z.imag (); | |
| 323 | |
| 324 ierr = 0; | |
| 325 | |
| 326 if (zr == 0.0 && zi == 0.0) | |
| 3146 | 327 { |
| 3220 | 328 yr = -octave_Inf; |
| 329 yi = 0.0; | |
| 330 } | |
| 331 else | |
| 332 { | |
| 4506 | 333 F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, 2, 1, &yr, &yi, nz, |
| 334 &wr, &wi, ierr); | |
| 335 | |
| 336 if (kode != 2) | |
| 337 { | |
| 338 double expz = exp (std::abs (zi)); | |
| 339 yr *= expz; | |
| 340 yi *= expz; | |
| 341 } | |
| 3220 | 342 |
| 4490 | 343 if (zi == 0.0 && zr >= 0.0) |
| 3220 | 344 yi = 0.0; |
| 3146 | 345 } |
| 346 | |
| 3220 | 347 return bessel_return_value (Complex (yr, yi), ierr); |
| 348 } | |
| 4911 | 349 else if (is_integer_value (alpha - 0.5)) |
| 350 { | |
| 351 // zbesy can overflow as z->0, and cause troubles for generic case below | |
| 352 alpha = -alpha; | |
| 353 Complex tmp = zbesj (z, alpha, kode, ierr); | |
| 354 if ((static_cast <long> (alpha - 0.5)) & 1) | |
| 355 tmp = - tmp; | |
| 356 retval = bessel_return_value (tmp, ierr); | |
| 357 } | |
| 3220 | 358 else |
| 359 { | |
| 360 alpha = -alpha; | |
| 3146 | 361 |
| 3220 | 362 Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr); |
| 3146 | 363 |
| 3220 | 364 if (ierr == 0 || ierr == 3) |
| 3146 | 365 { |
| 3220 | 366 tmp += sin (M_PI * alpha) * zbesj (z, alpha, kode, ierr); |
| 367 | |
| 368 retval = bessel_return_value (tmp, ierr); | |
| 369 } | |
| 370 else | |
| 371 retval = Complex (octave_NaN, octave_NaN); | |
| 372 } | |
| 373 | |
| 374 return retval; | |
| 375 } | |
| 376 | |
| 377 static inline Complex | |
| 5275 | 378 zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr) |
| 3220 | 379 { |
| 380 Complex retval; | |
| 3146 | 381 |
| 3220 | 382 if (alpha >= 0.0) |
| 383 { | |
| 384 double yr = 0.0; | |
| 385 double yi = 0.0; | |
| 386 | |
| 5275 | 387 octave_idx_type nz; |
| 3146 | 388 |
| 3220 | 389 double zr = z.real (); |
| 390 double zi = z.imag (); | |
| 391 | |
| 4506 | 392 F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr); |
| 393 | |
| 394 if (kode != 2) | |
| 395 { | |
| 396 double expz = exp (std::abs (zr)); | |
| 397 yr *= expz; | |
| 398 yi *= expz; | |
| 399 } | |
| 3146 | 400 |
| 4490 | 401 if (zi == 0.0 && zr >= 0.0) |
| 3220 | 402 yi = 0.0; |
| 403 | |
| 404 retval = bessel_return_value (Complex (yr, yi), ierr); | |
| 3146 | 405 } |
| 406 else | |
| 3220 | 407 { |
| 408 alpha = -alpha; | |
| 409 | |
| 410 Complex tmp = zbesi (z, alpha, kode, ierr); | |
| 411 | |
| 412 if (ierr == 0 || ierr == 3) | |
| 413 { | |
|
11904
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
414 Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha) |
| 7176 | 415 * zbesk (z, alpha, kode, ierr); |
|
11904
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
416 |
|
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
417 if (kode == 2) |
|
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
418 { |
|
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
419 // Compensate for different scaling factor of besk. |
|
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
420 tmp2 *= exp(-z - std::abs(z.real())); |
|
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
421 } |
|
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
422 |
|
059fadc0cbc3
fix scaling factor for negative alpha in zbesi,cbesi
Brian Gough <bjg@gnu.org>
parents:
7176
diff
changeset
|
423 tmp += tmp2; |
| 3220 | 424 |
| 425 retval = bessel_return_value (tmp, ierr); | |
| 426 } | |
| 427 else | |
| 428 retval = Complex (octave_NaN, octave_NaN); | |
| 429 } | |
| 430 | |
| 431 return retval; | |
| 432 } | |
| 433 | |
| 434 static inline Complex | |
| 5275 | 435 zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr) |
| 3220 | 436 { |
| 437 Complex retval; | |
| 438 | |
| 439 if (alpha >= 0.0) | |
| 440 { | |
| 441 double yr = 0.0; | |
| 442 double yi = 0.0; | |
| 443 | |
| 5275 | 444 octave_idx_type nz; |
| 3220 | 445 |
| 446 double zr = z.real (); | |
| 447 double zi = z.imag (); | |
| 448 | |
| 449 ierr = 0; | |
| 450 | |
| 451 if (zr == 0.0 && zi == 0.0) | |
| 452 { | |
| 453 yr = octave_Inf; | |
| 454 yi = 0.0; | |
| 455 } | |
| 456 else | |
| 457 { | |
| 4506 | 458 F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr); |
| 459 | |
| 460 if (kode != 2) | |
| 461 { | |
| 462 Complex expz = exp (-z); | |
| 463 | |
| 464 double rexpz = real (expz); | |
| 465 double iexpz = imag (expz); | |
| 466 | |
| 467 double tmp = yr*rexpz - yi*iexpz; | |
| 468 | |
| 469 yi = yr*iexpz + yi*rexpz; | |
| 470 yr = tmp; | |
| 471 } | |
| 3220 | 472 |
| 4490 | 473 if (zi == 0.0 && zr >= 0.0) |
| 3220 | 474 yi = 0.0; |
| 475 } | |
| 476 | |
| 477 retval = bessel_return_value (Complex (yr, yi), ierr); | |
| 478 } | |
| 479 else | |
| 480 { | |
| 481 Complex tmp = zbesk (z, -alpha, kode, ierr); | |
| 482 | |
| 483 retval = bessel_return_value (tmp, ierr); | |
| 484 } | |
| 3146 | 485 |
| 486 return retval; | |
| 487 } | |
| 488 | |
| 3220 | 489 static inline Complex |
| 5275 | 490 zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr) |
| 3146 | 491 { |
| 3220 | 492 Complex retval; |
| 3146 | 493 |
| 3220 | 494 if (alpha >= 0.0) |
| 3146 | 495 { |
| 3220 | 496 double yr = 0.0; |
| 497 double yi = 0.0; | |
| 498 | |
| 5275 | 499 octave_idx_type nz; |
| 3220 | 500 |
| 501 double zr = z.real (); | |
| 502 double zi = z.imag (); | |
| 3146 | 503 |
| 4506 | 504 F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz, ierr); |
| 505 | |
| 506 if (kode != 2) | |
| 507 { | |
| 508 Complex expz = exp (Complex (0.0, 1.0) * z); | |
| 509 | |
| 510 double rexpz = real (expz); | |
| 511 double iexpz = imag (expz); | |
| 512 | |
| 513 double tmp = yr*rexpz - yi*iexpz; | |
| 514 | |
| 515 yi = yr*iexpz + yi*rexpz; | |
| 516 yr = tmp; | |
| 517 } | |
| 3146 | 518 |
| 3220 | 519 retval = bessel_return_value (Complex (yr, yi), ierr); |
| 520 } | |
| 521 else | |
| 522 { | |
| 523 alpha = -alpha; | |
| 524 | |
| 525 static const Complex eye = Complex (0.0, 1.0); | |
| 3146 | 526 |
| 3220 | 527 Complex tmp = exp (M_PI * alpha * eye) * zbesh1 (z, alpha, kode, ierr); |
| 3146 | 528 |
| 3220 | 529 retval = bessel_return_value (tmp, ierr); |
| 530 } | |
| 3146 | 531 |
| 3220 | 532 return retval; |
| 533 } | |
| 3146 | 534 |
| 3220 | 535 static inline Complex |
| 5275 | 536 zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr) |
| 3220 | 537 { |
| 538 Complex retval; | |
| 3146 | 539 |
| 3220 | 540 if (alpha >= 0.0) |
| 541 { | |
| 542 double yr = 0.0; | |
| 543 double yi = 0.0; | |
| 544 | |
| 5275 | 545 octave_idx_type nz; |
| 3146 | 546 |
| 3220 | 547 double zr = z.real (); |
| 548 double zi = z.imag (); | |
| 3146 | 549 |
| 4506 | 550 F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz, ierr); |
| 551 | |
| 552 if (kode != 2) | |
| 553 { | |
| 554 Complex expz = exp (-Complex (0.0, 1.0) * z); | |
| 555 | |
| 556 double rexpz = real (expz); | |
| 557 double iexpz = imag (expz); | |
| 558 | |
| 559 double tmp = yr*rexpz - yi*iexpz; | |
| 560 | |
| 561 yi = yr*iexpz + yi*rexpz; | |
| 562 yr = tmp; | |
| 563 } | |
| 3220 | 564 |
| 565 retval = bessel_return_value (Complex (yr, yi), ierr); | |
| 3146 | 566 } |
| 567 else | |
| 3220 | 568 { |
| 569 alpha = -alpha; | |
| 570 | |
| 571 static const Complex eye = Complex (0.0, 1.0); | |
| 572 | |
| 573 Complex tmp = exp (-M_PI * alpha * eye) * zbesh2 (z, alpha, kode, ierr); | |
| 574 | |
| 575 retval = bessel_return_value (tmp, ierr); | |
| 576 } | |
| 577 | |
| 578 return retval; | |
| 579 } | |
| 580 | |
| 5275 | 581 typedef Complex (*fptr) (const Complex&, double, int, octave_idx_type&); |
| 3220 | 582 |
| 583 static inline Complex | |
| 584 do_bessel (fptr f, const char *, double alpha, const Complex& x, | |
| 5275 | 585 bool scaled, octave_idx_type& ierr) |
| 3220 | 586 { |
| 587 Complex retval; | |
| 588 | |
| 589 retval = f (x, alpha, (scaled ? 2 : 1), ierr); | |
| 590 | |
| 591 return retval; | |
| 592 } | |
| 593 | |
| 594 static inline ComplexMatrix | |
| 595 do_bessel (fptr f, const char *, double alpha, const ComplexMatrix& x, | |
| 5275 | 596 bool scaled, Array2<octave_idx_type>& ierr) |
| 3220 | 597 { |
| 5275 | 598 octave_idx_type nr = x.rows (); |
| 599 octave_idx_type nc = x.cols (); | |
| 3220 | 600 |
| 601 ComplexMatrix retval (nr, nc); | |
| 602 | |
| 603 ierr.resize (nr, nc); | |
| 604 | |
| 5275 | 605 for (octave_idx_type j = 0; j < nc; j++) |
| 606 for (octave_idx_type i = 0; i < nr; i++) | |
| 3220 | 607 retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j)); |
| 608 | |
| 609 return retval; | |
| 610 } | |
| 611 | |
| 612 static inline ComplexMatrix | |
| 613 do_bessel (fptr f, const char *, const Matrix& alpha, const Complex& x, | |
| 5275 | 614 bool scaled, Array2<octave_idx_type>& ierr) |
| 3220 | 615 { |
| 5275 | 616 octave_idx_type nr = alpha.rows (); |
| 617 octave_idx_type nc = alpha.cols (); | |
| 3220 | 618 |
| 619 ComplexMatrix retval (nr, nc); | |
| 620 | |
| 621 ierr.resize (nr, nc); | |
| 622 | |
| 5275 | 623 for (octave_idx_type j = 0; j < nc; j++) |
| 624 for (octave_idx_type i = 0; i < nr; i++) | |
| 3220 | 625 retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j)); |
| 3146 | 626 |
| 627 return retval; | |
| 628 } | |
| 629 | |
| 3220 | 630 static inline ComplexMatrix |
| 631 do_bessel (fptr f, const char *fn, const Matrix& alpha, | |
| 5275 | 632 const ComplexMatrix& x, bool scaled, Array2<octave_idx_type>& ierr) |
| 3146 | 633 { |
| 3220 | 634 ComplexMatrix retval; |
| 635 | |
| 5275 | 636 octave_idx_type x_nr = x.rows (); |
| 637 octave_idx_type x_nc = x.cols (); | |
| 3220 | 638 |
| 5275 | 639 octave_idx_type alpha_nr = alpha.rows (); |
| 640 octave_idx_type alpha_nc = alpha.cols (); | |
| 3220 | 641 |
| 642 if (x_nr == alpha_nr && x_nc == alpha_nc) | |
| 643 { | |
| 5275 | 644 octave_idx_type nr = x_nr; |
| 645 octave_idx_type nc = x_nc; | |
| 3220 | 646 |
| 647 retval.resize (nr, nc); | |
| 648 | |
| 649 ierr.resize (nr, nc); | |
| 650 | |
| 5275 | 651 for (octave_idx_type j = 0; j < nc; j++) |
| 652 for (octave_idx_type i = 0; i < nr; i++) | |
| 3220 | 653 retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j)); |
| 654 } | |
| 655 else | |
| 656 (*current_liboctave_error_handler) | |
| 657 ("%s: the sizes of alpha and x must conform", fn); | |
| 658 | |
| 659 return retval; | |
| 3146 | 660 } |
| 661 | |
| 4844 | 662 static inline ComplexNDArray |
| 663 do_bessel (fptr f, const char *, double alpha, const ComplexNDArray& x, | |
| 5275 | 664 bool scaled, ArrayN<octave_idx_type>& ierr) |
| 4844 | 665 { |
| 666 dim_vector dv = x.dims (); | |
| 5275 | 667 octave_idx_type nel = dv.numel (); |
| 4844 | 668 ComplexNDArray retval (dv); |
| 669 | |
| 670 ierr.resize (dv); | |
| 671 | |
| 5275 | 672 for (octave_idx_type i = 0; i < nel; i++) |
| 4844 | 673 retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i)); |
| 674 | |
| 675 return retval; | |
| 676 } | |
| 677 | |
| 678 static inline ComplexNDArray | |
| 679 do_bessel (fptr f, const char *, const NDArray& alpha, const Complex& x, | |
| 5275 | 680 bool scaled, ArrayN<octave_idx_type>& ierr) |
| 4844 | 681 { |
| 682 dim_vector dv = alpha.dims (); | |
| 5275 | 683 octave_idx_type nel = dv.numel (); |
| 4844 | 684 ComplexNDArray retval (dv); |
| 685 | |
| 686 ierr.resize (dv); | |
| 687 | |
| 5275 | 688 for (octave_idx_type i = 0; i < nel; i++) |
| 4844 | 689 retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i)); |
| 690 | |
| 691 return retval; | |
| 692 } | |
| 693 | |
| 694 static inline ComplexNDArray | |
| 695 do_bessel (fptr f, const char *fn, const NDArray& alpha, | |
| 5275 | 696 const ComplexNDArray& x, bool scaled, ArrayN<octave_idx_type>& ierr) |
| 4844 | 697 { |
| 698 dim_vector dv = x.dims (); | |
| 699 ComplexNDArray retval; | |
| 700 | |
| 701 if (dv == alpha.dims ()) | |
| 702 { | |
| 5275 | 703 octave_idx_type nel = dv.numel (); |
| 4844 | 704 |
| 705 retval.resize (dv); | |
| 706 ierr.resize (dv); | |
| 707 | |
| 5275 | 708 for (octave_idx_type i = 0; i < nel; i++) |
| 4844 | 709 retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i)); |
| 710 } | |
| 711 else | |
| 712 (*current_liboctave_error_handler) | |
| 713 ("%s: the sizes of alpha and x must conform", fn); | |
| 714 | |
| 715 return retval; | |
| 716 } | |
| 717 | |
| 3220 | 718 static inline ComplexMatrix |
| 719 do_bessel (fptr f, const char *, const RowVector& alpha, | |
| 5275 | 720 const ComplexColumnVector& x, bool scaled, Array2<octave_idx_type>& ierr) |
| 3146 | 721 { |
| 5275 | 722 octave_idx_type nr = x.length (); |
| 723 octave_idx_type nc = alpha.length (); | |
| 3220 | 724 |
| 725 ComplexMatrix retval (nr, nc); | |
| 3146 | 726 |
| 3220 | 727 ierr.resize (nr, nc); |
| 728 | |
| 5275 | 729 for (octave_idx_type j = 0; j < nc; j++) |
| 730 for (octave_idx_type i = 0; i < nr; i++) | |
| 3220 | 731 retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j)); |
| 732 | |
| 733 return retval; | |
| 3146 | 734 } |
| 735 | |
| 3220 | 736 #define SS_BESSEL(name, fcn) \ |
| 737 Complex \ | |
| 5275 | 738 name (double alpha, const Complex& x, bool scaled, octave_idx_type& ierr) \ |
| 3220 | 739 { \ |
| 740 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ | |
| 741 } | |
| 742 | |
| 743 #define SM_BESSEL(name, fcn) \ | |
| 744 ComplexMatrix \ | |
| 745 name (double alpha, const ComplexMatrix& x, bool scaled, \ | |
| 5275 | 746 Array2<octave_idx_type>& ierr) \ |
| 3220 | 747 { \ |
| 748 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ | |
| 749 } | |
| 750 | |
| 751 #define MS_BESSEL(name, fcn) \ | |
| 752 ComplexMatrix \ | |
| 753 name (const Matrix& alpha, const Complex& x, bool scaled, \ | |
| 5275 | 754 Array2<octave_idx_type>& ierr) \ |
| 3220 | 755 { \ |
| 756 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ | |
| 757 } | |
| 758 | |
| 759 #define MM_BESSEL(name, fcn) \ | |
| 760 ComplexMatrix \ | |
| 761 name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \ | |
| 5275 | 762 Array2<octave_idx_type>& ierr) \ |
| 3220 | 763 { \ |
| 764 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ | |
| 765 } | |
| 766 | |
| 4844 | 767 #define SN_BESSEL(name, fcn) \ |
| 768 ComplexNDArray \ | |
| 769 name (double alpha, const ComplexNDArray& x, bool scaled, \ | |
| 5275 | 770 ArrayN<octave_idx_type>& ierr) \ |
| 4844 | 771 { \ |
| 772 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ | |
| 773 } | |
| 774 | |
| 775 #define NS_BESSEL(name, fcn) \ | |
| 776 ComplexNDArray \ | |
| 777 name (const NDArray& alpha, const Complex& x, bool scaled, \ | |
| 5275 | 778 ArrayN<octave_idx_type>& ierr) \ |
| 4844 | 779 { \ |
| 780 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ | |
| 781 } | |
| 782 | |
| 783 #define NN_BESSEL(name, fcn) \ | |
| 784 ComplexNDArray \ | |
| 785 name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \ | |
| 5275 | 786 ArrayN<octave_idx_type>& ierr) \ |
| 4844 | 787 { \ |
| 788 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ | |
| 789 } | |
| 790 | |
| 3220 | 791 #define RC_BESSEL(name, fcn) \ |
| 792 ComplexMatrix \ | |
| 793 name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \ | |
| 5275 | 794 Array2<octave_idx_type>& ierr) \ |
| 3220 | 795 { \ |
| 796 return do_bessel (fcn, #name, alpha, x, scaled, ierr); \ | |
| 797 } | |
| 798 | |
| 799 #define ALL_BESSEL(name, fcn) \ | |
| 800 SS_BESSEL (name, fcn) \ | |
| 801 SM_BESSEL (name, fcn) \ | |
| 802 MS_BESSEL (name, fcn) \ | |
| 803 MM_BESSEL (name, fcn) \ | |
| 4844 | 804 SN_BESSEL (name, fcn) \ |
| 805 NS_BESSEL (name, fcn) \ | |
| 806 NN_BESSEL (name, fcn) \ | |
| 3220 | 807 RC_BESSEL (name, fcn) |
| 808 | |
| 809 ALL_BESSEL (besselj, zbesj) | |
| 810 ALL_BESSEL (bessely, zbesy) | |
| 811 ALL_BESSEL (besseli, zbesi) | |
| 812 ALL_BESSEL (besselk, zbesk) | |
| 813 ALL_BESSEL (besselh1, zbesh1) | |
| 814 ALL_BESSEL (besselh2, zbesh2) | |
| 815 | |
| 816 Complex | |
| 5275 | 817 airy (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr) |
| 3146 | 818 { |
| 3220 | 819 double ar = 0.0; |
| 820 double ai = 0.0; | |
| 821 | |
| 5275 | 822 octave_idx_type nz; |
| 3220 | 823 |
| 824 double zr = z.real (); | |
| 825 double zi = z.imag (); | |
| 3146 | 826 |
| 5275 | 827 octave_idx_type id = deriv ? 1 : 0; |
| 3220 | 828 |
| 4506 | 829 F77_FUNC (zairy, ZAIRY) (zr, zi, id, 2, ar, ai, nz, ierr); |
| 830 | |
| 831 if (! scaled) | |
| 832 { | |
| 833 Complex expz = exp (- 2.0 / 3.0 * z * sqrt(z)); | |
| 3220 | 834 |
| 4506 | 835 double rexpz = real (expz); |
| 836 double iexpz = imag (expz); | |
| 837 | |
| 838 double tmp = ar*rexpz - ai*iexpz; | |
| 839 | |
| 840 ai = ar*iexpz + ai*rexpz; | |
| 841 ar = tmp; | |
| 842 } | |
| 3220 | 843 |
| 4490 | 844 if (zi == 0.0 && (! scaled || zr >= 0.0)) |
| 3225 | 845 ai = 0.0; |
| 846 | |
| 3220 | 847 return bessel_return_value (Complex (ar, ai), ierr); |
| 3146 | 848 } |
| 849 | |
| 3220 | 850 Complex |
| 5275 | 851 biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr) |
| 3146 | 852 { |
| 3220 | 853 double ar = 0.0; |
| 854 double ai = 0.0; | |
| 855 | |
| 856 double zr = z.real (); | |
| 857 double zi = z.imag (); | |
| 858 | |
| 5275 | 859 octave_idx_type id = deriv ? 1 : 0; |
| 3220 | 860 |
| 4506 | 861 F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 2, ar, ai, ierr); |
| 862 | |
| 863 if (! scaled) | |
| 864 { | |
| 865 Complex expz = exp (std::abs (real (2.0 / 3.0 * z * sqrt (z)))); | |
| 3220 | 866 |
| 4506 | 867 double rexpz = real (expz); |
| 868 double iexpz = imag (expz); | |
| 869 | |
| 870 double tmp = ar*rexpz - ai*iexpz; | |
| 871 | |
| 872 ai = ar*iexpz + ai*rexpz; | |
| 873 ar = tmp; | |
| 874 } | |
| 3220 | 875 |
| 4490 | 876 if (zi == 0.0 && (! scaled || zr >= 0.0)) |
| 3225 | 877 ai = 0.0; |
| 878 | |
| 3220 | 879 return bessel_return_value (Complex (ar, ai), ierr); |
| 3146 | 880 } |
| 881 | |
| 3220 | 882 ComplexMatrix |
| 5275 | 883 airy (const ComplexMatrix& z, bool deriv, bool scaled, Array2<octave_idx_type>& ierr) |
| 3146 | 884 { |
| 5275 | 885 octave_idx_type nr = z.rows (); |
| 886 octave_idx_type nc = z.cols (); | |
| 3220 | 887 |
| 888 ComplexMatrix retval (nr, nc); | |
| 889 | |
| 890 ierr.resize (nr, nc); | |
| 891 | |
| 5275 | 892 for (octave_idx_type j = 0; j < nc; j++) |
| 893 for (octave_idx_type i = 0; i < nr; i++) | |
| 3220 | 894 retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j)); |
| 895 | |
| 896 return retval; | |
| 3146 | 897 } |
| 898 | |
| 3220 | 899 ComplexMatrix |
| 5275 | 900 biry (const ComplexMatrix& z, bool deriv, bool scaled, Array2<octave_idx_type>& ierr) |
| 3146 | 901 { |
| 5275 | 902 octave_idx_type nr = z.rows (); |
| 903 octave_idx_type nc = z.cols (); | |
| 3220 | 904 |
| 905 ComplexMatrix retval (nr, nc); | |
| 906 | |
| 907 ierr.resize (nr, nc); | |
| 908 | |
| 5275 | 909 for (octave_idx_type j = 0; j < nc; j++) |
| 910 for (octave_idx_type i = 0; i < nr; i++) | |
| 3220 | 911 retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j)); |
| 912 | |
| 913 return retval; | |
| 3146 | 914 } |
| 915 | |
| 4844 | 916 ComplexNDArray |
| 5275 | 917 airy (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<octave_idx_type>& ierr) |
| 4844 | 918 { |
| 919 dim_vector dv = z.dims (); | |
| 5275 | 920 octave_idx_type nel = dv.numel (); |
| 4844 | 921 ComplexNDArray retval (dv); |
| 922 | |
| 923 ierr.resize (dv); | |
| 924 | |
| 5275 | 925 for (octave_idx_type i = 0; i < nel; i++) |
| 4844 | 926 retval (i) = airy (z(i), deriv, scaled, ierr(i)); |
| 927 | |
| 928 return retval; | |
| 929 } | |
| 930 | |
| 931 ComplexNDArray | |
| 5275 | 932 biry (const ComplexNDArray& z, bool deriv, bool scaled, ArrayN<octave_idx_type>& ierr) |
| 4844 | 933 { |
| 934 dim_vector dv = z.dims (); | |
| 5275 | 935 octave_idx_type nel = dv.numel (); |
| 4844 | 936 ComplexNDArray retval (dv); |
| 937 | |
| 938 ierr.resize (dv); | |
| 939 | |
| 5275 | 940 for (octave_idx_type i = 0; i < nel; i++) |
| 4844 | 941 retval (i) = biry (z(i), deriv, scaled, ierr(i)); |
| 942 | |
| 943 return retval; | |
| 944 } | |
| 945 | |
| 3146 | 946 static void |
| 5275 | 947 gripe_betainc_nonconformant (octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2, octave_idx_type r3, |
| 948 octave_idx_type c3) | |
| 3146 | 949 { |
| 950 (*current_liboctave_error_handler) | |
| 951 ("betainc: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)", | |
| 952 r1, c1, r2, c2, r3, c3); | |
| 953 } | |
| 954 | |
| 4844 | 955 static dim_vector null_dims (0); |
| 956 | |
| 957 static void | |
| 958 gripe_betainc_nonconformant (const dim_vector& d1, const dim_vector& d2, | |
| 959 const dim_vector& d3) | |
| 960 { | |
| 961 std::string d1_str = d1.str (); | |
| 962 std::string d2_str = d2.str (); | |
| 963 std::string d3_str = d3.str (); | |
| 964 | |
| 965 (*current_liboctave_error_handler) | |
| 966 ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)", | |
| 967 d1_str.c_str (), d2_str.c_str (), d3_str.c_str ()); | |
| 968 } | |
| 969 | |
| 3146 | 970 double |
| 971 betainc (double x, double a, double b) | |
| 972 { | |
| 973 double retval; | |
| 5700 | 974 F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval)); |
| 3146 | 975 return retval; |
| 976 } | |
| 977 | |
| 978 Matrix | |
| 979 betainc (double x, double a, const Matrix& b) | |
| 980 { | |
| 5275 | 981 octave_idx_type nr = b.rows (); |
| 982 octave_idx_type nc = b.cols (); | |
| 3146 | 983 |
| 984 Matrix retval (nr, nc); | |
| 985 | |
| 5275 | 986 for (octave_idx_type j = 0; j < nc; j++) |
| 987 for (octave_idx_type i = 0; i < nr; i++) | |
| 3146 | 988 retval(i,j) = betainc (x, a, b(i,j)); |
| 989 | |
| 990 return retval; | |
| 991 } | |
| 992 | |
| 993 Matrix | |
| 994 betainc (double x, const Matrix& a, double b) | |
| 995 { | |
| 5275 | 996 octave_idx_type nr = a.rows (); |
| 997 octave_idx_type nc = a.cols (); | |
| 3146 | 998 |
| 999 Matrix retval (nr, nc); | |
| 1000 | |
| 5275 | 1001 for (octave_idx_type j = 0; j < nc; j++) |
| 1002 for (octave_idx_type i = 0; i < nr; i++) | |
| 3146 | 1003 retval(i,j) = betainc (x, a(i,j), b); |
| 1004 | |
| 1005 return retval; | |
| 1006 } | |
| 1007 | |
| 1008 Matrix | |
| 1009 betainc (double x, const Matrix& a, const Matrix& b) | |
| 1010 { | |
| 1011 Matrix retval; | |
| 1012 | |
| 5275 | 1013 octave_idx_type a_nr = a.rows (); |
| 1014 octave_idx_type a_nc = a.cols (); | |
| 3146 | 1015 |
| 5275 | 1016 octave_idx_type b_nr = b.rows (); |
| 1017 octave_idx_type b_nc = b.cols (); | |
| 3146 | 1018 |
| 1019 if (a_nr == b_nr && a_nc == b_nc) | |
| 1020 { | |
| 1021 retval.resize (a_nr, a_nc); | |
| 1022 | |
| 5275 | 1023 for (octave_idx_type j = 0; j < a_nc; j++) |
| 1024 for (octave_idx_type i = 0; i < a_nr; i++) | |
| 3146 | 1025 retval(i,j) = betainc (x, a(i,j), b(i,j)); |
| 1026 } | |
| 1027 else | |
| 1028 gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc); | |
| 1029 | |
| 1030 return retval; | |
| 1031 } | |
| 1032 | |
| 4844 | 1033 NDArray |
| 1034 betainc (double x, double a, const NDArray& b) | |
| 1035 { | |
| 1036 dim_vector dv = b.dims (); | |
| 1037 int nel = dv.numel (); | |
| 1038 | |
| 1039 NDArray retval (dv); | |
| 1040 | |
| 1041 for (int i = 0; i < nel; i++) | |
| 1042 retval (i) = betainc (x, a, b(i)); | |
| 1043 | |
| 1044 return retval; | |
| 1045 } | |
| 1046 | |
| 1047 NDArray | |
| 1048 betainc (double x, const NDArray& a, double b) | |
| 1049 { | |
| 1050 dim_vector dv = a.dims (); | |
| 1051 int nel = dv.numel (); | |
| 1052 | |
| 1053 NDArray retval (dv); | |
| 1054 | |
| 1055 for (int i = 0; i < nel; i++) | |
| 1056 retval (i) = betainc (x, a(i), b); | |
| 1057 | |
| 1058 return retval; | |
| 1059 } | |
| 1060 | |
| 1061 NDArray | |
| 1062 betainc (double x, const NDArray& a, const NDArray& b) | |
| 1063 { | |
| 1064 NDArray retval; | |
| 1065 dim_vector dv = a.dims (); | |
| 1066 | |
| 1067 if (dv == b.dims ()) | |
| 1068 { | |
| 1069 int nel = dv.numel (); | |
| 1070 | |
| 1071 retval.resize (dv); | |
| 1072 | |
| 1073 for (int i = 0; i < nel; i++) | |
| 1074 retval (i) = betainc (x, a(i), b(i)); | |
| 1075 } | |
| 1076 else | |
| 1077 gripe_betainc_nonconformant (dim_vector (0), dv, b.dims ()); | |
| 1078 | |
| 1079 return retval; | |
| 1080 } | |
| 1081 | |
| 1082 | |
| 3146 | 1083 Matrix |
| 1084 betainc (const Matrix& x, double a, double b) | |
| 1085 { | |
| 5275 | 1086 octave_idx_type nr = x.rows (); |
| 1087 octave_idx_type nc = x.cols (); | |
| 3146 | 1088 |
| 1089 Matrix retval (nr, nc); | |
| 1090 | |
| 5275 | 1091 for (octave_idx_type j = 0; j < nc; j++) |
| 1092 for (octave_idx_type i = 0; i < nr; i++) | |
| 3146 | 1093 retval(i,j) = betainc (x(i,j), a, b); |
| 1094 | |
| 1095 return retval; | |
| 1096 } | |
| 1097 | |
| 1098 Matrix | |
| 1099 betainc (const Matrix& x, double a, const Matrix& b) | |
| 1100 { | |
| 1101 Matrix retval; | |
| 1102 | |
| 5275 | 1103 octave_idx_type nr = x.rows (); |
| 1104 octave_idx_type nc = x.cols (); | |
| 3146 | 1105 |
| 5275 | 1106 octave_idx_type b_nr = b.rows (); |
| 1107 octave_idx_type b_nc = b.cols (); | |
| 3146 | 1108 |
| 1109 if (nr == b_nr && nc == b_nc) | |
| 1110 { | |
| 1111 retval.resize (nr, nc); | |
| 1112 | |
| 5275 | 1113 for (octave_idx_type j = 0; j < nc; j++) |
| 1114 for (octave_idx_type i = 0; i < nr; i++) | |
| 3146 | 1115 retval(i,j) = betainc (x(i,j), a, b(i,j)); |
| 1116 } | |
| 1117 else | |
| 1118 gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc); | |
| 1119 | |
| 1120 return retval; | |
| 1121 } | |
| 1122 | |
| 1123 Matrix | |
| 1124 betainc (const Matrix& x, const Matrix& a, double b) | |
| 1125 { | |
| 1126 Matrix retval; | |
| 1127 | |
| 5275 | 1128 octave_idx_type nr = x.rows (); |
| 1129 octave_idx_type nc = x.cols (); | |
| 3146 | 1130 |
| 5275 | 1131 octave_idx_type a_nr = a.rows (); |
| 1132 octave_idx_type a_nc = a.cols (); | |
| 3146 | 1133 |
| 1134 if (nr == a_nr && nc == a_nc) | |
| 1135 { | |
| 1136 retval.resize (nr, nc); | |
| 1137 | |
| 5275 | 1138 for (octave_idx_type j = 0; j < nc; j++) |
| 1139 for (octave_idx_type i = 0; i < nr; i++) | |
| 3146 | 1140 retval(i,j) = betainc (x(i,j), a(i,j), b); |
| 1141 } | |
| 1142 else | |
| 1143 gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1); | |
| 1144 | |
| 1145 return retval; | |
| 1146 } | |
| 1147 | |
| 1148 Matrix | |
| 1149 betainc (const Matrix& x, const Matrix& a, const Matrix& b) | |
| 1150 { | |
| 1151 Matrix retval; | |
| 1152 | |
| 5275 | 1153 octave_idx_type nr = x.rows (); |
| 1154 octave_idx_type nc = x.cols (); | |
| 3146 | 1155 |
| 5275 | 1156 octave_idx_type a_nr = a.rows (); |
| 1157 octave_idx_type a_nc = a.cols (); | |
| 3146 | 1158 |
| 5275 | 1159 octave_idx_type b_nr = b.rows (); |
| 1160 octave_idx_type b_nc = b.cols (); | |
| 3146 | 1161 |
| 1162 if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc) | |
| 1163 { | |
| 1164 retval.resize (nr, nc); | |
| 1165 | |
| 5275 | 1166 for (octave_idx_type j = 0; j < nc; j++) |
| 1167 for (octave_idx_type i = 0; i < nr; i++) | |
| 3146 | 1168 retval(i,j) = betainc (x(i,j), a(i,j), b(i,j)); |
| 1169 } | |
| 1170 else | |
| 1171 gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc); | |
| 1172 | |
| 1173 return retval; | |
| 1174 } | |
| 1175 | |
| 4844 | 1176 NDArray |
| 1177 betainc (const NDArray& x, double a, double b) | |
| 1178 { | |
| 1179 dim_vector dv = x.dims (); | |
| 1180 int nel = dv.numel (); | |
| 1181 | |
| 1182 NDArray retval (dv); | |
| 1183 | |
| 1184 for (int i = 0; i < nel; i++) | |
| 1185 retval (i) = betainc (x(i), a, b); | |
| 1186 | |
| 1187 return retval; | |
| 1188 } | |
| 1189 | |
| 1190 NDArray | |
| 1191 betainc (const NDArray& x, double a, const NDArray& b) | |
| 1192 { | |
| 1193 NDArray retval; | |
| 1194 dim_vector dv = x.dims (); | |
| 1195 | |
| 1196 if (dv == b.dims ()) | |
| 1197 { | |
| 1198 int nel = dv.numel (); | |
| 1199 | |
| 1200 retval.resize (dv); | |
| 1201 | |
| 1202 for (int i = 0; i < nel; i++) | |
| 1203 retval (i) = betainc (x(i), a, b(i)); | |
| 1204 } | |
| 1205 else | |
| 1206 gripe_betainc_nonconformant (dv, dim_vector (0), b.dims ()); | |
| 1207 | |
| 1208 return retval; | |
| 1209 } | |
| 1210 | |
| 1211 NDArray | |
| 1212 betainc (const NDArray& x, const NDArray& a, double b) | |
| 1213 { | |
| 1214 NDArray retval; | |
| 1215 dim_vector dv = x.dims (); | |
| 1216 | |
| 1217 if (dv == a.dims ()) | |
| 1218 { | |
| 1219 int nel = dv.numel (); | |
| 1220 | |
| 1221 retval.resize (dv); | |
| 1222 | |
| 1223 for (int i = 0; i < nel; i++) | |
| 1224 retval (i) = betainc (x(i), a(i), b); | |
| 1225 } | |
| 1226 else | |
| 1227 gripe_betainc_nonconformant (dv, a.dims (), dim_vector (0)); | |
| 1228 | |
| 1229 return retval; | |
| 1230 } | |
| 1231 | |
| 1232 NDArray | |
| 1233 betainc (const NDArray& x, const NDArray& a, const NDArray& b) | |
| 1234 { | |
| 1235 NDArray retval; | |
| 1236 dim_vector dv = x.dims (); | |
| 1237 | |
| 1238 if (dv == a.dims () && dv == b.dims ()) | |
| 1239 { | |
| 1240 int nel = dv.numel (); | |
| 1241 | |
| 1242 retval.resize (dv); | |
| 1243 | |
| 1244 for (int i = 0; i < nel; i++) | |
| 1245 retval (i) = betainc (x(i), a(i), b(i)); | |
| 1246 } | |
| 1247 else | |
| 1248 gripe_betainc_nonconformant (dv, a.dims (), b.dims ()); | |
| 1249 | |
| 1250 return retval; | |
| 1251 } | |
| 1252 | |
| 5775 | 1253 // FIXME -- there is still room for improvement here... |
| 3164 | 1254 |
| 3146 | 1255 double |
| 4004 | 1256 gammainc (double x, double a, bool& err) |
| 3146 | 1257 { |
| 1258 double retval; | |
| 3164 | 1259 |
| 4004 | 1260 err = false; |
| 3164 | 1261 |
| 4004 | 1262 if (a < 0.0 || x < 0.0) |
| 1263 { | |
| 1264 (*current_liboctave_error_handler) | |
| 1265 ("gammainc: A and X must be non-negative"); | |
| 1266 | |
| 1267 err = true; | |
| 1268 } | |
| 1269 else | |
| 5278 | 1270 F77_XFCN (xgammainc, XGAMMAINC, (a, x, retval)); |
| 3164 | 1271 |
| 3146 | 1272 return retval; |
| 1273 } | |
| 1274 | |
| 1275 Matrix | |
| 1276 gammainc (double x, const Matrix& a) | |
| 1277 { | |
| 5275 | 1278 octave_idx_type nr = a.rows (); |
| 1279 octave_idx_type nc = a.cols (); | |
| 3146 | 1280 |
| 4004 | 1281 Matrix result (nr, nc); |
| 1282 Matrix retval; | |
| 1283 | |
| 1284 bool err; | |
| 3146 | 1285 |
| 5275 | 1286 for (octave_idx_type j = 0; j < nc; j++) |
| 1287 for (octave_idx_type i = 0; i < nr; i++) | |
| 4004 | 1288 { |
| 1289 result(i,j) = gammainc (x, a(i,j), err); | |
| 1290 | |
| 1291 if (err) | |
| 1292 goto done; | |
| 1293 } | |
| 1294 | |
| 1295 retval = result; | |
| 1296 | |
| 1297 done: | |
| 3146 | 1298 |
| 1299 return retval; | |
| 1300 } | |
| 1301 | |
| 1302 Matrix | |
| 1303 gammainc (const Matrix& x, double a) | |
| 1304 { | |
| 5275 | 1305 octave_idx_type nr = x.rows (); |
| 1306 octave_idx_type nc = x.cols (); | |
| 3146 | 1307 |
| 4004 | 1308 Matrix result (nr, nc); |
| 1309 Matrix retval; | |
| 1310 | |
| 1311 bool err; | |
| 3146 | 1312 |
| 5275 | 1313 for (octave_idx_type j = 0; j < nc; j++) |
| 1314 for (octave_idx_type i = 0; i < nr; i++) | |
| 4004 | 1315 { |
| 1316 result(i,j) = gammainc (x(i,j), a, err); | |
| 1317 | |
| 1318 if (err) | |
| 1319 goto done; | |
| 1320 } | |
| 1321 | |
| 1322 retval = result; | |
| 1323 | |
| 1324 done: | |
| 3146 | 1325 |
| 1326 return retval; | |
| 1327 } | |
| 1328 | |
| 1329 Matrix | |
| 1330 gammainc (const Matrix& x, const Matrix& a) | |
| 1331 { | |
| 4004 | 1332 Matrix result; |
| 3146 | 1333 Matrix retval; |
| 1334 | |
| 5275 | 1335 octave_idx_type nr = x.rows (); |
| 1336 octave_idx_type nc = x.cols (); | |
| 3146 | 1337 |
| 5275 | 1338 octave_idx_type a_nr = a.rows (); |
| 1339 octave_idx_type a_nc = a.cols (); | |
| 3146 | 1340 |
| 1341 if (nr == a_nr && nc == a_nc) | |
| 1342 { | |
| 4004 | 1343 result.resize (nr, nc); |
| 1344 | |
| 1345 bool err; | |
| 3146 | 1346 |
| 5275 | 1347 for (octave_idx_type j = 0; j < nc; j++) |
| 1348 for (octave_idx_type i = 0; i < nr; i++) | |
| 4004 | 1349 { |
| 1350 result(i,j) = gammainc (x(i,j), a(i,j), err); | |
| 1351 | |
| 1352 if (err) | |
| 1353 goto done; | |
| 1354 } | |
| 1355 | |
| 1356 retval = result; | |
| 3146 | 1357 } |
| 1358 else | |
| 1359 (*current_liboctave_error_handler) | |
| 1360 ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)", | |
| 1361 nr, nc, a_nr, a_nc); | |
| 1362 | |
| 4004 | 1363 done: |
| 1364 | |
| 3146 | 1365 return retval; |
| 1366 } | |
| 1367 | |
| 4844 | 1368 NDArray |
| 1369 gammainc (double x, const NDArray& a) | |
| 1370 { | |
| 1371 dim_vector dv = a.dims (); | |
| 1372 int nel = dv.numel (); | |
| 1373 | |
| 1374 NDArray retval; | |
| 1375 NDArray result (dv); | |
| 1376 | |
| 1377 bool err; | |
| 1378 | |
| 1379 for (int i = 0; i < nel; i++) | |
| 1380 { | |
| 1381 result (i) = gammainc (x, a(i), err); | |
| 1382 | |
| 1383 if (err) | |
| 1384 goto done; | |
| 1385 } | |
| 1386 | |
| 1387 retval = result; | |
| 1388 | |
| 1389 done: | |
| 1390 | |
| 1391 return retval; | |
| 1392 } | |
| 1393 | |
| 1394 NDArray | |
| 1395 gammainc (const NDArray& x, double a) | |
| 1396 { | |
| 1397 dim_vector dv = x.dims (); | |
| 1398 int nel = dv.numel (); | |
| 1399 | |
| 1400 NDArray retval; | |
| 1401 NDArray result (dv); | |
| 1402 | |
| 1403 bool err; | |
| 1404 | |
| 1405 for (int i = 0; i < nel; i++) | |
| 1406 { | |
| 1407 result (i) = gammainc (x(i), a, err); | |
| 1408 | |
| 1409 if (err) | |
| 1410 goto done; | |
| 1411 } | |
| 1412 | |
| 1413 retval = result; | |
| 1414 | |
| 1415 done: | |
| 1416 | |
| 1417 return retval; | |
| 1418 } | |
| 1419 | |
| 1420 NDArray | |
| 1421 gammainc (const NDArray& x, const NDArray& a) | |
| 1422 { | |
| 1423 dim_vector dv = x.dims (); | |
| 1424 int nel = dv.numel (); | |
| 1425 | |
| 1426 NDArray retval; | |
| 1427 NDArray result; | |
| 1428 | |
| 1429 if (dv == a.dims ()) | |
| 1430 { | |
| 1431 result.resize (dv); | |
| 1432 | |
| 1433 bool err; | |
| 1434 | |
| 1435 for (int i = 0; i < nel; i++) | |
| 1436 { | |
| 1437 result (i) = gammainc (x(i), a(i), err); | |
| 1438 | |
| 1439 if (err) | |
| 1440 goto done; | |
| 1441 } | |
| 1442 | |
| 1443 retval = result; | |
| 1444 } | |
| 1445 else | |
| 1446 { | |
| 1447 std::string x_str = dv.str (); | |
| 1448 std::string a_str = a.dims ().str (); | |
| 1449 | |
| 1450 (*current_liboctave_error_handler) | |
| 1451 ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)", | |
| 1452 x_str.c_str (), a_str. c_str ()); | |
| 1453 } | |
| 1454 | |
| 1455 done: | |
| 1456 | |
| 1457 return retval; | |
| 1458 } | |
| 1459 | |
| 3146 | 1460 /* |
| 1461 ;;; Local Variables: *** | |
| 1462 ;;; mode: C++ *** | |
| 1463 ;;; End: *** | |
| 1464 */ |