GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
lo-specfun.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2017 John W. Eaton
4 Copyright (C) 2007-2010 D. Martin
5 Copyright (C) 2010 Jaroslav Hajek
6 Copyright (C) 2010 VZLU Prague
7 Copyright (C) 2016 CarnĂ« Draug
8 
9 This file is part of Octave.
10 
11 Octave is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
15 
16 Octave is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19 for more details.
20 
21 You should have received a copy of the GNU General Public License
22 along with Octave; see the file COPYING. If not, see
23 <http://www.gnu.org/licenses/>.
24 
25 */
26 
27 #if defined (HAVE_CONFIG_H)
28 # include "config.h"
29 #endif
30 
31 #include "Range.h"
32 #include "CColVector.h"
33 #include "CMatrix.h"
34 #include "dRowVector.h"
35 #include "dMatrix.h"
36 #include "dNDArray.h"
37 #include "CNDArray.h"
38 #include "fCColVector.h"
39 #include "fCMatrix.h"
40 #include "fRowVector.h"
41 #include "fMatrix.h"
42 #include "fNDArray.h"
43 #include "fCNDArray.h"
44 #include "f77-fcn.h"
45 #include "lo-amos-proto.h"
46 #include "lo-error.h"
47 #include "lo-ieee.h"
48 #include "lo-slatec-proto.h"
49 #include "lo-specfun.h"
50 #include "mx-inlines.cc"
51 #include "lo-mappers.h"
52 #include "lo-math.h"
53 
54 #include "Faddeeva.hh"
55 
56 namespace octave
57 {
58  namespace math
59  {
60  double
61  acosh (double x)
62  {
63 #if defined (HAVE_ACOSH)
64  return ::acosh (x);
65 #else
66  double retval;
67  F77_XFCN (xdacosh, XDACOSH, (x, retval));
68  return retval;
69 #endif
70  }
71 
72  float
73  acosh (float x)
74  {
75 #if defined (HAVE_ACOSHF)
76  return acoshf (x);
77 #else
78  float retval;
79  F77_XFCN (xacosh, XACOSH, (x, retval));
80  return retval;
81 #endif
82  }
83 
84  Complex
85  acosh (const Complex& x)
86  {
87 #if defined (HAVE_COMPLEX_STD_ACOSH)
88  return std::acosh (x);
89 #else
90  return log (x + sqrt (x + 1.0) * sqrt (x - 1.0));
91 #endif
92  }
93 
95  acosh (const FloatComplex& x)
96  {
97 #if defined (HAVE_COMPLEX_STD_ACOSH)
98  return std::acosh (x);
99 #else
100  return log (x + sqrt (x + 1.0f) * sqrt (x - 1.0f));
101 #endif
102  }
103 
104  double
105  asinh (double x)
106  {
107 #if defined (HAVE_ASINH)
108  return ::asinh (x);
109 #else
110  double retval;
111  F77_XFCN (xdasinh, XDASINH, (x, retval));
112  return retval;
113 #endif
114  }
115 
116  float
117  asinh (float x)
118  {
119 #if defined (HAVE_ASINHF)
120  return asinhf (x);
121 #else
122  float retval;
123  F77_XFCN (xasinh, XASINH, (x, retval));
124  return retval;
125 #endif
126  }
127 
128  Complex
129  asinh (const Complex& x)
130  {
131 #if defined (HAVE_COMPLEX_STD_ASINH)
132  return std::asinh (x);
133 #else
134  return log (x + sqrt (x*x + 1.0));
135 #endif
136  }
137 
140  {
141 #if defined (HAVE_COMPLEX_STD_ASINH)
142  return std::asinh (x);
143 #else
144  return log (x + sqrt (x*x + 1.0f));
145 #endif
146  }
147 
148  double
149  atanh (double x)
150  {
151 #if defined (HAVE_ATANH)
152  return ::atanh (x);
153 #else
154  double retval;
155  F77_XFCN (xdatanh, XDATANH, (x, retval));
156  return retval;
157 #endif
158  }
159 
160  float
161  atanh (float x)
162  {
163 #if defined (HAVE_ATANHF)
164  return atanhf (x);
165 #else
166  float retval;
167  F77_XFCN (xatanh, XATANH, (x, retval));
168  return retval;
169 #endif
170  }
171 
172  Complex
173  atanh (const Complex& x)
174  {
175 #if defined (HAVE_COMPLEX_STD_ATANH)
176  return std::atanh (x);
177 #else
178  return log ((1.0 + x) / (1.0 - x)) / 2.0;
179 #endif
180  }
181 
184  {
185 #if defined (HAVE_COMPLEX_STD_ATANH)
186  return std::atanh (x);
187 #else
188  return log ((1.0f + x) / (1.0f - x)) / 2.0f;
189 #endif
190  }
191 
192  double
193  erf (double x)
194  {
195 #if defined (HAVE_ERF)
196  return ::erf (x);
197 #else
198  double retval;
199  F77_XFCN (xderf, XDERF, (x, retval));
200  return retval;
201 #endif
202  }
203 
204  float
205  erf (float x)
206  {
207 #if defined (HAVE_ERFF)
208  return erff (x);
209 #else
210  float retval;
211  F77_XFCN (xerf, XERF, (x, retval));
212  return retval;
213 #endif
214  }
215 
216  // Complex error function from the Faddeeva package
217  Complex
218  erf (const Complex& x)
219  {
220  return Faddeeva::erf (x);
221  }
222 
224  erf (const FloatComplex& x)
225  {
226  Complex xd (x.real (), x.imag ());
227  Complex ret = Faddeeva::erf (xd, std::numeric_limits<float>::epsilon ());
228  return FloatComplex (ret.real (), ret.imag ());
229  }
230 
231  double
232  erfc (double x)
233  {
234 #if defined (HAVE_ERFC)
235  return ::erfc (x);
236 #else
237  double retval;
238  F77_XFCN (xderfc, XDERFC, (x, retval));
239  return retval;
240 #endif
241  }
242 
243  float
244  erfc (float x)
245  {
246 #if defined (HAVE_ERFCF)
247  return erfcf (x);
248 #else
249  float retval;
250  F77_XFCN (xerfc, XERFC, (x, retval));
251  return retval;
252 #endif
253  }
254 
255  // Complex complementary error function from the Faddeeva package
256  Complex
257  erfc (const Complex& x)
258  {
259  return Faddeeva::erfc (x);
260  }
261 
263  erfc (const FloatComplex& x)
264  {
265  Complex xd (x.real (), x.imag ());
266  Complex ret = Faddeeva::erfc (xd, std::numeric_limits<float>::epsilon ());
267  return FloatComplex (ret.real (), ret.imag ());
268  }
269 
270  // Real and complex scaled complementary error function from Faddeeva package
271  float erfcx (float x) { return Faddeeva::erfcx(x); }
272  double erfcx (double x) { return Faddeeva::erfcx(x); }
273 
274  Complex
275  erfcx (const Complex& x)
276  {
277  return Faddeeva::erfcx (x);
278  }
279 
282  {
283  Complex xd (x.real (), x.imag ());
284  Complex ret = Faddeeva::erfcx (xd, std::numeric_limits<float>::epsilon ());
285  return FloatComplex (ret.real (), ret.imag ());
286  }
287 
288  // Real and complex imaginary error function from Faddeeva package
289  float erfi (float x) { return Faddeeva::erfi(x); }
290  double erfi (double x) { return Faddeeva::erfi(x); }
291 
292  Complex
293  erfi (const Complex& x)
294  {
295  return Faddeeva::erfi (x);
296  }
297 
299  erfi (const FloatComplex& x)
300  {
301  Complex xd (x.real (), x.imag ());
302  Complex ret = Faddeeva::erfi (xd, std::numeric_limits<float>::epsilon ());
303  return FloatComplex (ret.real (), ret.imag ());
304  }
305 
306  // Real and complex Dawson function (= scaled erfi) from Faddeeva package
307  float dawson (float x) { return Faddeeva::Dawson(x); }
308  double dawson (double x) { return Faddeeva::Dawson(x); }
309 
310  Complex
311  dawson (const Complex& x)
312  {
313  return Faddeeva::Dawson (x);
314  }
315 
318  {
319  Complex xd (x.real (), x.imag ());
320  Complex ret = Faddeeva::Dawson (xd, std::numeric_limits<float>::epsilon ());
321  return FloatComplex (ret.real (), ret.imag ());
322  }
323 
324  double
325  gamma (double x)
326  {
327  double result;
328 
329  // Special cases for (near) compatibility with Matlab instead of
330  // tgamma. Matlab does not have -0.
331 
332  if (x == 0)
333  result = (octave::math::negative_sign (x)
336  else if ((x < 0 && octave::math::x_nint (x) == x) || octave::math::isinf (x))
338  else if (octave::math::isnan (x))
340  else
341  {
342 #if defined (HAVE_TGAMMA)
343  result = tgamma (x);
344 #else
345  F77_XFCN (xdgamma, XDGAMMA, (x, result));
346 #endif
347  }
348 
349  return result;
350  }
351 
352  double
353  lgamma (double x)
354  {
355 #if defined (HAVE_LGAMMA)
356  return ::lgamma (x);
357 #else
358  double result;
359  double sgngam;
360 
361  if (octave::math::isnan (x))
362  result = x;
363  else if ((x <= 0 && octave::math::x_nint (x) == x) || octave::math::isinf (x))
365  else
366  F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));
367 
368  return result;
369 #endif
370  }
371 
372  Complex
373  rc_lgamma (double x)
374  {
375  double result;
376 
377 #if defined (HAVE_LGAMMA_R)
378  int sgngam;
379  result = lgamma_r (x, &sgngam);
380 #else
381  double sgngam = 0.0;
382 
383  if (octave::math::isnan (x))
384  result = x;
385  else if ((x <= 0 && octave::math::x_nint (x) == x) || octave::math::isinf (x))
387  else
388  F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));
389 
390 #endif
391 
392  if (sgngam < 0)
393  return result + Complex (0., M_PI);
394  else
395  return result;
396  }
397 
398  float
399  gamma (float x)
400  {
401  float result;
402 
403  // Special cases for (near) compatibility with Matlab instead of
404  // tgamma. Matlab does not have -0.
405 
406  if (x == 0)
407  result = (octave::math::negative_sign (x)
410  else if ((x < 0 && octave::math::x_nint (x) == x) || octave::math::isinf (x))
412  else if (octave::math::isnan (x))
414  else
415  {
416 #if defined (HAVE_TGAMMA)
417  result = tgammaf (x);
418 #else
419  F77_XFCN (xgamma, XGAMMA, (x, result));
420 #endif
421  }
422 
423  return result;
424  }
425 
426  float
427  lgamma (float x)
428  {
429 #if defined (HAVE_LGAMMAF)
430  return lgammaf (x);
431 #else
432  float result;
433  float sgngam;
434 
435  if (octave::math::isnan (x))
436  result = x;
437  else if ((x <= 0 && octave::math::x_nint (x) == x) || octave::math::isinf (x))
439  else
440  F77_XFCN (algams, ALGAMS, (x, result, sgngam));
441 
442  return result;
443 #endif
444  }
445 
447  rc_lgamma (float x)
448  {
449  float result;
450 
451 #if defined (HAVE_LGAMMAF_R)
452  int sgngam;
453  result = lgammaf_r (x, &sgngam);
454 #else
455  float sgngam = 0.0f;
456 
457  if (octave::math::isnan (x))
458  result = x;
459  else if ((x <= 0 && octave::math::x_nint (x) == x) || octave::math::isinf (x))
461  else
462  F77_XFCN (algams, ALGAMS, (x, result, sgngam));
463 
464 #endif
465 
466  if (sgngam < 0)
467  return result + FloatComplex (0., M_PI);
468  else
469  return result;
470  }
471 
472  double
473  expm1 (double x)
474  {
475 #if defined (HAVE_EXPM1)
476  return ::expm1 (x);
477 #else
478  double retval;
479 
480  double ax = fabs (x);
481 
482  if (ax < 0.1)
483  {
484  ax /= 16;
485 
486  // use Taylor series to calculate exp(x)-1.
487  double t = ax;
488  double s = 0;
489  for (int i = 2; i < 7; i++)
490  s += (t *= ax/i);
491  s += ax;
492 
493  // use the identity (a+1)^2-1 = a*(a+2)
494  double e = s;
495  for (int i = 0; i < 4; i++)
496  {
497  s *= e + 2;
498  e *= e + 2;
499  }
500 
501  retval = (x > 0) ? s : -s / (1+s);
502  }
503  else
504  retval = exp (x) - 1;
505 
506  return retval;
507 #endif
508  }
509 
510  Complex
511  expm1 (const Complex& x)
512  {
513  Complex retval;
514 
515  if (std:: abs (x) < 1)
516  {
517  double im = x.imag ();
518  double u = expm1 (x.real ());
519  double v = sin (im/2);
520  v = -2*v*v;
521  retval = Complex (u*v + u + v, (u+1) * sin (im));
522  }
523  else
524  retval = std::exp (x) - Complex (1);
525 
526  return retval;
527  }
528 
529  float
530  expm1 (float x)
531  {
532 #if defined (HAVE_EXPM1F)
533  return expm1f (x);
534 #else
535  float retval;
536 
537  float ax = fabs (x);
538 
539  if (ax < 0.1)
540  {
541  ax /= 16;
542 
543  // use Taylor series to calculate exp(x)-1.
544  float t = ax;
545  float s = 0;
546  for (int i = 2; i < 7; i++)
547  s += (t *= ax/i);
548  s += ax;
549 
550  // use the identity (a+1)^2-1 = a*(a+2)
551  float e = s;
552  for (int i = 0; i < 4; i++)
553  {
554  s *= e + 2;
555  e *= e + 2;
556  }
557 
558  retval = (x > 0) ? s : -s / (1+s);
559  }
560  else
561  retval = exp (x) - 1;
562 
563  return retval;
564 #endif
565  }
566 
569  {
571 
572  if (std:: abs (x) < 1)
573  {
574  float im = x.imag ();
575  float u = expm1 (x.real ());
576  float v = sin (im/2);
577  v = -2*v*v;
578  retval = FloatComplex (u*v + u + v, (u+1) * sin (im));
579  }
580  else
581  retval = std::exp (x) - FloatComplex (1);
582 
583  return retval;
584  }
585 
586  double
587  log1p (double x)
588  {
589 #if defined (HAVE_LOG1P)
590  return ::log1p (x);
591 #else
592  double retval;
593 
594  double ax = fabs (x);
595 
596  if (ax < 0.2)
597  {
598  // approximation log (1+x) ~ 2*sum ((x/(2+x)).^ii ./ ii), ii = 1:2:2n+1
599  double u = x / (2 + x), t = 1, s = 0;
600  for (int i = 2; i < 12; i += 2)
601  s += (t *= u*u) / (i+1);
602 
603  retval = 2 * (s + 1) * u;
604  }
605  else
606  retval = std::log (1 + x);
607 
608  return retval;
609 #endif
610  }
611 
612  Complex
613  log1p (const Complex& x)
614  {
615  Complex retval;
616 
617  double r = x.real (), i = x.imag ();
618 
619  if (fabs (r) < 0.5 && fabs (i) < 0.5)
620  {
621  double u = 2*r + r*r + i*i;
622  retval = Complex (log1p (u / (1+sqrt (u+1))),
623  atan2 (1 + r, i));
624  }
625  else
626  retval = std::log (Complex (1) + x);
627 
628  return retval;
629  }
630 
631  template <typename T>
632  T
633  xcbrt (T x)
634  {
635  static const T one_third = 0.3333333333333333333f;
636  if (octave::math::finite (x))
637  {
638  // Use pow.
639  T y = std::pow (std::abs (x), one_third) * signum (x);
640  // Correct for better accuracy.
641  return (x / (y*y) + y + y) / 3;
642  }
643  else
644  return x;
645  }
646 
647  double
648  cbrt (double x)
649  {
650 #if defined (HAVE_CBRT)
651  return ::cbrt (x);
652 #else
653  return xxcbrt (x);
654 #endif
655  }
656 
657  float
658  log1p (float x)
659  {
660 #if defined (HAVE_LOG1PF)
661  return log1pf (x);
662 #else
663  float retval;
664 
665  float ax = fabs (x);
666 
667  if (ax < 0.2)
668  {
669  // approximation log (1+x) ~ 2*sum ((x/(2+x)).^ii ./ ii), ii = 1:2:2n+1
670  float u = x / (2 + x), t = 1.0f, s = 0;
671  for (int i = 2; i < 12; i += 2)
672  s += (t *= u*u) / (i+1);
673 
674  retval = 2 * (s + 1.0f) * u;
675  }
676  else
677  retval = std::log (1.0f + x);
678 
679  return retval;
680 #endif
681  }
682 
685  {
687 
688  float r = x.real (), i = x.imag ();
689 
690  if (fabs (r) < 0.5 && fabs (i) < 0.5)
691  {
692  float u = 2*r + r*r + i*i;
693  retval = FloatComplex (log1p (u / (1+sqrt (u+1))),
694  atan2 (1 + r, i));
695  }
696  else
697  retval = std::log (FloatComplex (1) + x);
698 
699  return retval;
700  }
701 
702  float
703  cbrt (float x)
704  {
705 #if defined (HAVE_CBRTF)
706  return cbrtf (x);
707 #else
708  return xxcbrt (x);
709 #endif
710  }
711 
712  static inline Complex
713  zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
714 
715  static inline Complex
716  zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
717 
718  static inline Complex
719  zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
720 
721  static inline Complex
722  zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
723 
724  static inline Complex
725  zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
726 
727  static inline Complex
728  zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
729 
730  static inline Complex
732  {
733  static const Complex inf_val
736 
737  static const Complex nan_val
740 
741  Complex retval;
742 
743  switch (ierr)
744  {
745  case 0:
746  case 3:
747  retval = val;
748  break;
749 
750  case 2:
751  retval = inf_val;
752  break;
753 
754  default:
755  retval = nan_val;
756  break;
757  }
758 
759  return retval;
760  }
761 
762  static inline bool
764  {
765  return x == static_cast<double> (static_cast<long> (x));
766  }
767 
768  static inline Complex
769  zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
770  {
771  Complex retval;
772 
773  if (alpha >= 0.0)
774  {
775  double yr = 0.0;
776  double yi = 0.0;
777 
778  octave_idx_type nz;
779 
780  double zr = z.real ();
781  double zi = z.imag ();
782 
783  F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
784 
785  if (kode != 2)
786  {
787  double expz = exp (std::abs (zi));
788  yr *= expz;
789  yi *= expz;
790  }
791 
792  if (zi == 0.0 && zr >= 0.0)
793  yi = 0.0;
794 
795  retval = bessel_return_value (Complex (yr, yi), ierr);
796  }
797  else if (is_integer_value (alpha))
798  {
799  // zbesy can overflow as z->0, and cause troubles for generic case below
800  alpha = -alpha;
801  Complex tmp = zbesj (z, alpha, kode, ierr);
802  if ((static_cast<long> (alpha)) & 1)
803  tmp = - tmp;
804  retval = bessel_return_value (tmp, ierr);
805  }
806  else
807  {
808  alpha = -alpha;
809 
810  Complex tmp = cos (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
811 
812  if (ierr == 0 || ierr == 3)
813  {
814  tmp -= sin (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
815 
816  retval = bessel_return_value (tmp, ierr);
817  }
818  else
821  }
822 
823  return retval;
824  }
825 
826  static inline Complex
827  zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
828  {
829  Complex retval;
830 
831  if (alpha >= 0.0)
832  {
833  double yr = 0.0;
834  double yi = 0.0;
835 
836  octave_idx_type nz;
837 
838  double wr, wi;
839 
840  double zr = z.real ();
841  double zi = z.imag ();
842 
843  ierr = 0;
844 
845  if (zr == 0.0 && zi == 0.0)
846  {
848  yi = 0.0;
849  }
850  else
851  {
852  F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, 2, 1, &yr, &yi, nz,
853  &wr, &wi, ierr);
854 
855  if (kode != 2)
856  {
857  double expz = exp (std::abs (zi));
858  yr *= expz;
859  yi *= expz;
860  }
861 
862  if (zi == 0.0 && zr >= 0.0)
863  yi = 0.0;
864  }
865 
866  return bessel_return_value (Complex (yr, yi), ierr);
867  }
868  else if (is_integer_value (alpha - 0.5))
869  {
870  // zbesy can overflow as z->0, and cause troubles for generic case below
871  alpha = -alpha;
872  Complex tmp = zbesj (z, alpha, kode, ierr);
873  if ((static_cast<long> (alpha - 0.5)) & 1)
874  tmp = - tmp;
875  retval = bessel_return_value (tmp, ierr);
876  }
877  else
878  {
879  alpha = -alpha;
880 
881  Complex tmp = cos (M_PI * alpha) * zbesy (z, alpha, kode, ierr);
882 
883  if (ierr == 0 || ierr == 3)
884  {
885  tmp += sin (M_PI * alpha) * zbesj (z, alpha, kode, ierr);
886 
887  retval = bessel_return_value (tmp, ierr);
888  }
889  else
892  }
893 
894  return retval;
895  }
896 
897  static inline Complex
898  zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
899  {
900  Complex retval;
901 
902  if (alpha >= 0.0)
903  {
904  double yr = 0.0;
905  double yi = 0.0;
906 
907  octave_idx_type nz;
908 
909  double zr = z.real ();
910  double zi = z.imag ();
911 
912  F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
913 
914  if (kode != 2)
915  {
916  double expz = exp (std::abs (zr));
917  yr *= expz;
918  yi *= expz;
919  }
920 
921  if (zi == 0.0 && zr >= 0.0)
922  yi = 0.0;
923 
924  retval = bessel_return_value (Complex (yr, yi), ierr);
925  }
926  else if (is_integer_value (alpha))
927  {
928  // zbesi can overflow as z->0, and cause troubles for generic case below
929  alpha = -alpha;
930  Complex tmp = zbesi (z, alpha, kode, ierr);
931  retval = bessel_return_value (tmp, ierr);
932  }
933  else
934  {
935  alpha = -alpha;
936 
937  Complex tmp = zbesi (z, alpha, kode, ierr);
938 
939  if (ierr == 0 || ierr == 3)
940  {
941  Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha)
942  * zbesk (z, alpha, kode, ierr);
943 
944  if (kode == 2)
945  {
946  // Compensate for different scaling factor of besk.
947  tmp2 *= exp (-z - std::abs (z.real ()));
948  }
949 
950  tmp += tmp2;
951 
952  retval = bessel_return_value (tmp, ierr);
953  }
954  else
957  }
958 
959  return retval;
960  }
961 
962  static inline Complex
963  zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
964  {
965  Complex retval;
966 
967  if (alpha >= 0.0)
968  {
969  double yr = 0.0;
970  double yi = 0.0;
971 
972  octave_idx_type nz;
973 
974  double zr = z.real ();
975  double zi = z.imag ();
976 
977  ierr = 0;
978 
979  if (zr == 0.0 && zi == 0.0)
980  {
982  yi = 0.0;
983  }
984  else
985  {
986  F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
987 
988  if (kode != 2)
989  {
990  Complex expz = exp (-z);
991 
992  double rexpz = expz.real ();
993  double iexpz = expz.imag ();
994 
995  double tmp = yr*rexpz - yi*iexpz;
996 
997  yi = yr*iexpz + yi*rexpz;
998  yr = tmp;
999  }
1000 
1001  if (zi == 0.0 && zr >= 0.0)
1002  yi = 0.0;
1003  }
1004 
1005  retval = bessel_return_value (Complex (yr, yi), ierr);
1006  }
1007  else
1008  {
1009  Complex tmp = zbesk (z, -alpha, kode, ierr);
1010 
1011  retval = bessel_return_value (tmp, ierr);
1012  }
1013 
1014  return retval;
1015  }
1016 
1017  static inline Complex
1018  zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
1019  {
1020  Complex retval;
1021 
1022  if (alpha >= 0.0)
1023  {
1024  double yr = 0.0;
1025  double yi = 0.0;
1026 
1027  octave_idx_type nz;
1028 
1029  double zr = z.real ();
1030  double zi = z.imag ();
1031 
1032  F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz, ierr);
1033 
1034  if (kode != 2)
1035  {
1036  Complex expz = exp (Complex (0.0, 1.0) * z);
1037 
1038  double rexpz = expz.real ();
1039  double iexpz = expz.imag ();
1040 
1041  double tmp = yr*rexpz - yi*iexpz;
1042 
1043  yi = yr*iexpz + yi*rexpz;
1044  yr = tmp;
1045  }
1046 
1047  retval = bessel_return_value (Complex (yr, yi), ierr);
1048  }
1049  else
1050  {
1051  alpha = -alpha;
1052 
1053  static const Complex eye = Complex (0.0, 1.0);
1054 
1055  Complex tmp = exp (M_PI * alpha * eye) * zbesh1 (z, alpha, kode, ierr);
1056 
1057  retval = bessel_return_value (tmp, ierr);
1058  }
1059 
1060  return retval;
1061  }
1062 
1063  static inline Complex
1064  zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr)
1065  {
1066  Complex retval;
1067 
1068  if (alpha >= 0.0)
1069  {
1070  double yr = 0.0;
1071  double yi = 0.0;
1072 
1073  octave_idx_type nz;
1074 
1075  double zr = z.real ();
1076  double zi = z.imag ();
1077 
1078  F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz, ierr);
1079 
1080  if (kode != 2)
1081  {
1082  Complex expz = exp (-Complex (0.0, 1.0) * z);
1083 
1084  double rexpz = expz.real ();
1085  double iexpz = expz.imag ();
1086 
1087  double tmp = yr*rexpz - yi*iexpz;
1088 
1089  yi = yr*iexpz + yi*rexpz;
1090  yr = tmp;
1091  }
1092 
1093  retval = bessel_return_value (Complex (yr, yi), ierr);
1094  }
1095  else
1096  {
1097  alpha = -alpha;
1098 
1099  static const Complex eye = Complex (0.0, 1.0);
1100 
1101  Complex tmp = exp (-M_PI * alpha * eye) * zbesh2 (z, alpha, kode, ierr);
1102 
1103  retval = bessel_return_value (tmp, ierr);
1104  }
1105 
1106  return retval;
1107  }
1108 
1109  typedef Complex (*dptr) (const Complex&, double, int, octave_idx_type&);
1110 
1111  static inline Complex
1112  do_bessel (dptr f, const char *, double alpha, const Complex& x,
1113  bool scaled, octave_idx_type& ierr)
1114  {
1115  Complex retval;
1116 
1117  retval = f (x, alpha, (scaled ? 2 : 1), ierr);
1118 
1119  return retval;
1120  }
1121 
1122  static inline ComplexMatrix
1123  do_bessel (dptr f, const char *, double alpha, const ComplexMatrix& x,
1124  bool scaled, Array<octave_idx_type>& ierr)
1125  {
1126  octave_idx_type nr = x.rows ();
1127  octave_idx_type nc = x.cols ();
1128 
1129  ComplexMatrix retval (nr, nc);
1130 
1131  ierr.resize (dim_vector (nr, nc));
1132 
1133  for (octave_idx_type j = 0; j < nc; j++)
1134  for (octave_idx_type i = 0; i < nr; i++)
1135  retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j));
1136 
1137  return retval;
1138  }
1139 
1140  static inline ComplexMatrix
1141  do_bessel (dptr f, const char *, const Matrix& alpha, const Complex& x,
1142  bool scaled, Array<octave_idx_type>& ierr)
1143  {
1144  octave_idx_type nr = alpha.rows ();
1145  octave_idx_type nc = alpha.cols ();
1146 
1147  ComplexMatrix retval (nr, nc);
1148 
1149  ierr.resize (dim_vector (nr, nc));
1150 
1151  for (octave_idx_type j = 0; j < nc; j++)
1152  for (octave_idx_type i = 0; i < nr; i++)
1153  retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
1154 
1155  return retval;
1156  }
1157 
1158  static inline ComplexMatrix
1159  do_bessel (dptr f, const char *fn, const Matrix& alpha,
1160  const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr)
1161  {
1163 
1164  octave_idx_type x_nr = x.rows ();
1165  octave_idx_type x_nc = x.cols ();
1166 
1167  octave_idx_type alpha_nr = alpha.rows ();
1168  octave_idx_type alpha_nc = alpha.cols ();
1169 
1170  if (x_nr != alpha_nr || x_nc != alpha_nc)
1171  (*current_liboctave_error_handler)
1172  ("%s: the sizes of alpha and x must conform", fn);
1173 
1174  octave_idx_type nr = x_nr;
1175  octave_idx_type nc = x_nc;
1176 
1177  retval.resize (nr, nc);
1178 
1179  ierr.resize (dim_vector (nr, nc));
1180 
1181  for (octave_idx_type j = 0; j < nc; j++)
1182  for (octave_idx_type i = 0; i < nr; i++)
1183  retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
1184 
1185  return retval;
1186  }
1187 
1188  static inline ComplexNDArray
1189  do_bessel (dptr f, const char *, double alpha, const ComplexNDArray& x,
1190  bool scaled, Array<octave_idx_type>& ierr)
1191  {
1192  dim_vector dv = x.dims ();
1193  octave_idx_type nel = dv.numel ();
1194  ComplexNDArray retval (dv);
1195 
1196  ierr.resize (dv);
1197 
1198  for (octave_idx_type i = 0; i < nel; i++)
1199  retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i));
1200 
1201  return retval;
1202  }
1203 
1204  static inline ComplexNDArray
1205  do_bessel (dptr f, const char *, const NDArray& alpha, const Complex& x,
1206  bool scaled, Array<octave_idx_type>& ierr)
1207  {
1208  dim_vector dv = alpha.dims ();
1209  octave_idx_type nel = dv.numel ();
1210  ComplexNDArray retval (dv);
1211 
1212  ierr.resize (dv);
1213 
1214  for (octave_idx_type i = 0; i < nel; i++)
1215  retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i));
1216 
1217  return retval;
1218  }
1219 
1220  static inline ComplexNDArray
1221  do_bessel (dptr f, const char *fn, const NDArray& alpha,
1222  const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr)
1223  {
1224  dim_vector dv = x.dims ();
1226 
1227  if (dv != alpha.dims ())
1229  ("%s: the sizes of alpha and x must conform", fn);
1230 
1231  octave_idx_type nel = dv.numel ();
1232 
1233  retval.resize (dv);
1234  ierr.resize (dv);
1235 
1236  for (octave_idx_type i = 0; i < nel; i++)
1237  retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i));
1238 
1239  return retval;
1240  }
1241 
1242  static inline ComplexMatrix
1243  do_bessel (dptr f, const char *, const RowVector& alpha,
1244  const ComplexColumnVector& x, bool scaled,
1246  {
1247  octave_idx_type nr = x.numel ();
1248  octave_idx_type nc = alpha.numel ();
1249 
1250  ComplexMatrix retval (nr, nc);
1251 
1252  ierr.resize (dim_vector (nr, nc));
1253 
1254  for (octave_idx_type j = 0; j < nc; j++)
1255  for (octave_idx_type i = 0; i < nr; i++)
1256  retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j));
1257 
1258  return retval;
1259  }
1260 
1261 #define SS_BESSEL(name, fcn) \
1262  Complex \
1263  name (double alpha, const Complex& x, bool scaled, octave_idx_type& ierr) \
1264  { \
1265  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1266  }
1267 
1268 #define SM_BESSEL(name, fcn) \
1269  ComplexMatrix \
1270  name (double alpha, const ComplexMatrix& x, bool scaled, \
1271  Array<octave_idx_type>& ierr) \
1272  { \
1273  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1274  }
1275 
1276 #define MS_BESSEL(name, fcn) \
1277  ComplexMatrix \
1278  name (const Matrix& alpha, const Complex& x, bool scaled, \
1279  Array<octave_idx_type>& ierr) \
1280  { \
1281  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1282  }
1283 
1284 #define MM_BESSEL(name, fcn) \
1285  ComplexMatrix \
1286  name (const Matrix& alpha, const ComplexMatrix& x, bool scaled, \
1287  Array<octave_idx_type>& ierr) \
1288  { \
1289  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1290  }
1291 
1292 #define SN_BESSEL(name, fcn) \
1293  ComplexNDArray \
1294  name (double alpha, const ComplexNDArray& x, bool scaled, \
1295  Array<octave_idx_type>& ierr) \
1296  { \
1297  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1298  }
1299 
1300 #define NS_BESSEL(name, fcn) \
1301  ComplexNDArray \
1302  name (const NDArray& alpha, const Complex& x, bool scaled, \
1303  Array<octave_idx_type>& ierr) \
1304  { \
1305  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1306  }
1307 
1308 #define NN_BESSEL(name, fcn) \
1309  ComplexNDArray \
1310  name (const NDArray& alpha, const ComplexNDArray& x, bool scaled, \
1311  Array<octave_idx_type>& ierr) \
1312  { \
1313  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1314  }
1315 
1316 #define RC_BESSEL(name, fcn) \
1317  ComplexMatrix \
1318  name (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, \
1319  Array<octave_idx_type>& ierr) \
1320  { \
1321  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1322  }
1323 
1324 #define ALL_BESSEL(name, fcn) \
1325  SS_BESSEL (name, fcn) \
1326  SM_BESSEL (name, fcn) \
1327  MS_BESSEL (name, fcn) \
1328  MM_BESSEL (name, fcn) \
1329  SN_BESSEL (name, fcn) \
1330  NS_BESSEL (name, fcn) \
1331  NN_BESSEL (name, fcn) \
1332  RC_BESSEL (name, fcn)
1333 
1340 
1341 #undef ALL_BESSEL
1342 #undef SS_BESSEL
1343 #undef SM_BESSEL
1344 #undef MS_BESSEL
1345 #undef MM_BESSEL
1346 #undef SN_BESSEL
1347 #undef NS_BESSEL
1348 #undef NN_BESSEL
1349 #undef RC_BESSEL
1350 
1351  static inline FloatComplex
1352  cbesj (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1353 
1354  static inline FloatComplex
1355  cbesy (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1356 
1357  static inline FloatComplex
1358  cbesi (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1359 
1360  static inline FloatComplex
1361  cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1362 
1363  static inline FloatComplex
1364  cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1365 
1366  static inline FloatComplex
1367  cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1368 
1369  static inline FloatComplex
1371  {
1372  static const FloatComplex inf_val
1375 
1376  static const FloatComplex nan_val
1379 
1381 
1382  switch (ierr)
1383  {
1384  case 0:
1385  case 3:
1386  retval = val;
1387  break;
1388 
1389  case 2:
1390  retval = inf_val;
1391  break;
1392 
1393  default:
1394  retval = nan_val;
1395  break;
1396  }
1397 
1398  return retval;
1399  }
1400 
1401  static inline bool
1403  {
1404  return x == static_cast<float> (static_cast<long> (x));
1405  }
1406 
1407  static inline FloatComplex
1408  cbesj (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
1409  {
1411 
1412  if (alpha >= 0.0)
1413  {
1414  FloatComplex y = 0.0;
1415 
1416  octave_idx_type nz;
1417 
1418  F77_FUNC (cbesj, CBESJ) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1,
1419  F77_CMPLX_ARG (&y), nz, ierr);
1420 
1421  if (kode != 2)
1422  {
1423  float expz = exp (std::abs (z.imag ()));
1424  y *= expz;
1425  }
1426 
1427  if (z.imag () == 0.0 && z.real () >= 0.0)
1428  y = FloatComplex (y.real (), 0.0);
1429 
1430  retval = bessel_return_value (y, ierr);
1431  }
1432  else if (is_integer_value (alpha))
1433  {
1434  // zbesy can overflow as z->0, and cause troubles for generic case below
1435  alpha = -alpha;
1436  FloatComplex tmp = cbesj (z, alpha, kode, ierr);
1437  if ((static_cast<long> (alpha)) & 1)
1438  tmp = - tmp;
1439  retval = bessel_return_value (tmp, ierr);
1440  }
1441  else
1442  {
1443  alpha = -alpha;
1444 
1445  FloatComplex tmp = cosf (static_cast<float> (M_PI) * alpha)
1446  * cbesj (z, alpha, kode, ierr);
1447 
1448  if (ierr == 0 || ierr == 3)
1449  {
1450  tmp -= sinf (static_cast<float> (M_PI) * alpha)
1451  * cbesy (z, alpha, kode, ierr);
1452 
1453  retval = bessel_return_value (tmp, ierr);
1454  }
1455  else
1458  }
1459 
1460  return retval;
1461  }
1462 
1463  static inline FloatComplex
1464  cbesy (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
1465  {
1467 
1468  if (alpha >= 0.0)
1469  {
1470  FloatComplex y = 0.0;
1471 
1472  octave_idx_type nz;
1473 
1474  FloatComplex w;
1475 
1476  ierr = 0;
1477 
1478  if (z.real () == 0.0 && z.imag () == 0.0)
1479  {
1481  }
1482  else
1483  {
1484  F77_FUNC (cbesy, CBESY) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1,
1485  F77_CMPLX_ARG (&y), nz, F77_CMPLX_ARG (&w), ierr);
1486 
1487  if (kode != 2)
1488  {
1489  float expz = exp (std::abs (z.imag ()));
1490  y *= expz;
1491  }
1492 
1493  if (z.imag () == 0.0 && z.real () >= 0.0)
1494  y = FloatComplex (y.real (), 0.0);
1495  }
1496 
1497  return bessel_return_value (y, ierr);
1498  }
1499  else if (is_integer_value (alpha - 0.5))
1500  {
1501  // zbesy can overflow as z->0, and cause troubles for generic case below
1502  alpha = -alpha;
1503  FloatComplex tmp = cbesj (z, alpha, kode, ierr);
1504  if ((static_cast<long> (alpha - 0.5)) & 1)
1505  tmp = - tmp;
1506  retval = bessel_return_value (tmp, ierr);
1507  }
1508  else
1509  {
1510  alpha = -alpha;
1511 
1512  FloatComplex tmp = cosf (static_cast<float> (M_PI) * alpha)
1513  * cbesy (z, alpha, kode, ierr);
1514 
1515  if (ierr == 0 || ierr == 3)
1516  {
1517  tmp += sinf (static_cast<float> (M_PI) * alpha)
1518  * cbesj (z, alpha, kode, ierr);
1519 
1520  retval = bessel_return_value (tmp, ierr);
1521  }
1522  else
1525  }
1526 
1527  return retval;
1528  }
1529 
1530  static inline FloatComplex
1531  cbesi (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
1532  {
1534 
1535  if (alpha >= 0.0)
1536  {
1537  FloatComplex y = 0.0;
1538 
1539  octave_idx_type nz;
1540 
1541  F77_FUNC (cbesi, CBESI) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1,
1542  F77_CMPLX_ARG (&y), nz, ierr);
1543 
1544  if (kode != 2)
1545  {
1546  float expz = exp (std::abs (z.real ()));
1547  y *= expz;
1548  }
1549 
1550  if (z.imag () == 0.0 && z.real () >= 0.0)
1551  y = FloatComplex (y.real (), 0.0);
1552 
1553  retval = bessel_return_value (y, ierr);
1554  }
1555  else
1556  {
1557  alpha = -alpha;
1558 
1559  FloatComplex tmp = cbesi (z, alpha, kode, ierr);
1560 
1561  if (ierr == 0 || ierr == 3)
1562  {
1563  FloatComplex tmp2 = static_cast<float> (2.0 / M_PI)
1564  * sinf (static_cast<float> (M_PI) * alpha)
1565  * cbesk (z, alpha, kode, ierr);
1566 
1567  if (kode == 2)
1568  {
1569  // Compensate for different scaling factor of besk.
1570  tmp2 *= exp (-z - std::abs (z.real ()));
1571  }
1572 
1573  tmp += tmp2;
1574 
1575  retval = bessel_return_value (tmp, ierr);
1576  }
1577  else
1580  }
1581 
1582  return retval;
1583  }
1584 
1585  static inline FloatComplex
1586  cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
1587  {
1589 
1590  if (alpha >= 0.0)
1591  {
1592  FloatComplex y = 0.0;
1593 
1594  octave_idx_type nz;
1595 
1596  ierr = 0;
1597 
1598  if (z.real () == 0.0 && z.imag () == 0.0)
1599  {
1601  }
1602  else
1603  {
1604  F77_FUNC (cbesk, CBESK) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1,
1605  F77_CMPLX_ARG (&y), nz, ierr);
1606 
1607  if (kode != 2)
1608  {
1609  FloatComplex expz = exp (-z);
1610 
1611  float rexpz = expz.real ();
1612  float iexpz = expz.imag ();
1613 
1614  float tmp_r = y.real () * rexpz - y.imag () * iexpz;
1615  float tmp_i = y.real () * iexpz + y.imag () * rexpz;
1616 
1617  y = FloatComplex (tmp_r, tmp_i);
1618  }
1619 
1620  if (z.imag () == 0.0 && z.real () >= 0.0)
1621  y = FloatComplex (y.real (), 0.0);
1622  }
1623 
1624  retval = bessel_return_value (y, ierr);
1625  }
1626  else
1627  {
1628  FloatComplex tmp = cbesk (z, -alpha, kode, ierr);
1629 
1630  retval = bessel_return_value (tmp, ierr);
1631  }
1632 
1633  return retval;
1634  }
1635 
1636  static inline FloatComplex
1637  cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
1638  {
1640 
1641  if (alpha >= 0.0)
1642  {
1643  FloatComplex y = 0.0;
1644 
1645  octave_idx_type nz;
1646 
1647  F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, 1,
1648  F77_CMPLX_ARG (&y), nz, ierr);
1649 
1650  if (kode != 2)
1651  {
1652  FloatComplex expz = exp (FloatComplex (0.0, 1.0) * z);
1653 
1654  float rexpz = expz.real ();
1655  float iexpz = expz.imag ();
1656 
1657  float tmp_r = y.real () * rexpz - y.imag () * iexpz;
1658  float tmp_i = y.real () * iexpz + y.imag () * rexpz;
1659 
1660  y = FloatComplex (tmp_r, tmp_i);
1661  }
1662 
1663  retval = bessel_return_value (y, ierr);
1664  }
1665  else
1666  {
1667  alpha = -alpha;
1668 
1669  static const FloatComplex eye = FloatComplex (0.0, 1.0);
1670 
1671  FloatComplex tmp = exp (static_cast<float> (M_PI) * alpha * eye)
1672  * cbesh1 (z, alpha, kode, ierr);
1673 
1674  retval = bessel_return_value (tmp, ierr);
1675  }
1676 
1677  return retval;
1678  }
1679 
1680  static inline FloatComplex
1681  cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr)
1682  {
1684 
1685  if (alpha >= 0.0)
1686  {
1687  FloatComplex y = 0.0;
1688 
1689  octave_idx_type nz;
1690 
1691  F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 2, 1,
1692  F77_CMPLX_ARG (&y), nz, ierr);
1693 
1694  if (kode != 2)
1695  {
1696  FloatComplex expz = exp (-FloatComplex (0.0, 1.0) * z);
1697 
1698  float rexpz = expz.real ();
1699  float iexpz = expz.imag ();
1700 
1701  float tmp_r = y.real () * rexpz - y.imag () * iexpz;
1702  float tmp_i = y.real () * iexpz + y.imag () * rexpz;
1703 
1704  y = FloatComplex (tmp_r, tmp_i);
1705  }
1706 
1707  retval = bessel_return_value (y, ierr);
1708  }
1709  else
1710  {
1711  alpha = -alpha;
1712 
1713  static const FloatComplex eye = FloatComplex (0.0, 1.0);
1714 
1715  FloatComplex tmp = exp (-static_cast<float> (M_PI) * alpha * eye)
1716  * cbesh2 (z, alpha, kode, ierr);
1717 
1718  retval = bessel_return_value (tmp, ierr);
1719  }
1720 
1721  return retval;
1722  }
1723 
1724  typedef FloatComplex (*fptr) (const FloatComplex&, float, int,
1725  octave_idx_type&);
1726 
1727  static inline FloatComplex
1728  do_bessel (fptr f, const char *, float alpha, const FloatComplex& x,
1729  bool scaled, octave_idx_type& ierr)
1730  {
1731  FloatComplex retval;
1732 
1733  retval = f (x, alpha, (scaled ? 2 : 1), ierr);
1734 
1735  return retval;
1736  }
1737 
1738  static inline FloatComplexMatrix
1739  do_bessel (fptr f, const char *, float alpha, const FloatComplexMatrix& x,
1740  bool scaled, Array<octave_idx_type>& ierr)
1741  {
1742  octave_idx_type nr = x.rows ();
1743  octave_idx_type nc = x.cols ();
1744 
1745  FloatComplexMatrix retval (nr, nc);
1746 
1747  ierr.resize (dim_vector (nr, nc));
1748 
1749  for (octave_idx_type j = 0; j < nc; j++)
1750  for (octave_idx_type i = 0; i < nr; i++)
1751  retval(i,j) = f (x(i,j), alpha, (scaled ? 2 : 1), ierr(i,j));
1752 
1753  return retval;
1754  }
1755 
1756  static inline FloatComplexMatrix
1757  do_bessel (fptr f, const char *, const FloatMatrix& alpha,
1758  const FloatComplex& x,
1759  bool scaled, Array<octave_idx_type>& ierr)
1760  {
1761  octave_idx_type nr = alpha.rows ();
1762  octave_idx_type nc = alpha.cols ();
1763 
1764  FloatComplexMatrix retval (nr, nc);
1765 
1766  ierr.resize (dim_vector (nr, nc));
1767 
1768  for (octave_idx_type j = 0; j < nc; j++)
1769  for (octave_idx_type i = 0; i < nr; i++)
1770  retval(i,j) = f (x, alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
1771 
1772  return retval;
1773  }
1774 
1775  static inline FloatComplexMatrix
1776  do_bessel (fptr f, const char *fn, const FloatMatrix& alpha,
1777  const FloatComplexMatrix& x, bool scaled,
1779  {
1781 
1782  octave_idx_type x_nr = x.rows ();
1783  octave_idx_type x_nc = x.cols ();
1784 
1785  octave_idx_type alpha_nr = alpha.rows ();
1786  octave_idx_type alpha_nc = alpha.cols ();
1787 
1788  if (x_nr != alpha_nr || x_nc != alpha_nc)
1789  (*current_liboctave_error_handler)
1790  ("%s: the sizes of alpha and x must conform", fn);
1791 
1792  octave_idx_type nr = x_nr;
1793  octave_idx_type nc = x_nc;
1794 
1795  retval.resize (nr, nc);
1796 
1797  ierr.resize (dim_vector (nr, nc));
1798 
1799  for (octave_idx_type j = 0; j < nc; j++)
1800  for (octave_idx_type i = 0; i < nr; i++)
1801  retval(i,j) = f (x(i,j), alpha(i,j), (scaled ? 2 : 1), ierr(i,j));
1802 
1803  return retval;
1804  }
1805 
1806  static inline FloatComplexNDArray
1807  do_bessel (fptr f, const char *, float alpha, const FloatComplexNDArray& x,
1808  bool scaled, Array<octave_idx_type>& ierr)
1809  {
1810  dim_vector dv = x.dims ();
1811  octave_idx_type nel = dv.numel ();
1813 
1814  ierr.resize (dv);
1815 
1816  for (octave_idx_type i = 0; i < nel; i++)
1817  retval(i) = f (x(i), alpha, (scaled ? 2 : 1), ierr(i));
1818 
1819  return retval;
1820  }
1821 
1822  static inline FloatComplexNDArray
1823  do_bessel (fptr f, const char *, const FloatNDArray& alpha,
1824  const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr)
1825  {
1826  dim_vector dv = alpha.dims ();
1827  octave_idx_type nel = dv.numel ();
1829 
1830  ierr.resize (dv);
1831 
1832  for (octave_idx_type i = 0; i < nel; i++)
1833  retval(i) = f (x, alpha(i), (scaled ? 2 : 1), ierr(i));
1834 
1835  return retval;
1836  }
1837 
1838  static inline FloatComplexNDArray
1839  do_bessel (fptr f, const char *fn, const FloatNDArray& alpha,
1840  const FloatComplexNDArray& x, bool scaled,
1842  {
1843  dim_vector dv = x.dims ();
1845 
1846  if (dv != alpha.dims ())
1848  ("%s: the sizes of alpha and x must conform", fn);
1849 
1850  octave_idx_type nel = dv.numel ();
1851 
1852  retval.resize (dv);
1853  ierr.resize (dv);
1854 
1855  for (octave_idx_type i = 0; i < nel; i++)
1856  retval(i) = f (x(i), alpha(i), (scaled ? 2 : 1), ierr(i));
1857 
1858  return retval;
1859  }
1860 
1861  static inline FloatComplexMatrix
1862  do_bessel (fptr f, const char *, const FloatRowVector& alpha,
1863  const FloatComplexColumnVector& x, bool scaled,
1865  {
1866  octave_idx_type nr = x.numel ();
1867  octave_idx_type nc = alpha.numel ();
1868 
1869  FloatComplexMatrix retval (nr, nc);
1870 
1871  ierr.resize (dim_vector (nr, nc));
1872 
1873  for (octave_idx_type j = 0; j < nc; j++)
1874  for (octave_idx_type i = 0; i < nr; i++)
1875  retval(i,j) = f (x(i), alpha(j), (scaled ? 2 : 1), ierr(i,j));
1876 
1877  return retval;
1878  }
1879 
1880 #define SS_BESSEL(name, fcn) \
1881  FloatComplex \
1882  name (float alpha, const FloatComplex& x, bool scaled, octave_idx_type& ierr) \
1883  { \
1884  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1885  }
1886 
1887 #define SM_BESSEL(name, fcn) \
1888  FloatComplexMatrix \
1889  name (float alpha, const FloatComplexMatrix& x, bool scaled, \
1890  Array<octave_idx_type>& ierr) \
1891  { \
1892  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1893  }
1894 
1895 #define MS_BESSEL(name, fcn) \
1896  FloatComplexMatrix \
1897  name (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, \
1898  Array<octave_idx_type>& ierr) \
1899  { \
1900  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1901  }
1902 
1903 #define MM_BESSEL(name, fcn) \
1904  FloatComplexMatrix \
1905  name (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, \
1906  Array<octave_idx_type>& ierr) \
1907  { \
1908  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1909  }
1910 
1911 #define SN_BESSEL(name, fcn) \
1912  FloatComplexNDArray \
1913  name (float alpha, const FloatComplexNDArray& x, bool scaled, \
1914  Array<octave_idx_type>& ierr) \
1915  { \
1916  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1917  }
1918 
1919 #define NS_BESSEL(name, fcn) \
1920  FloatComplexNDArray \
1921  name (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, \
1922  Array<octave_idx_type>& ierr) \
1923  { \
1924  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1925  }
1926 
1927 #define NN_BESSEL(name, fcn) \
1928  FloatComplexNDArray \
1929  name (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, \
1930  Array<octave_idx_type>& ierr) \
1931  { \
1932  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1933  }
1934 
1935 #define RC_BESSEL(name, fcn) \
1936  FloatComplexMatrix \
1937  name (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, \
1938  Array<octave_idx_type>& ierr) \
1939  { \
1940  return do_bessel (fcn, #name, alpha, x, scaled, ierr); \
1941  }
1942 
1943 #define ALL_BESSEL(name, fcn) \
1944  SS_BESSEL (name, fcn) \
1945  SM_BESSEL (name, fcn) \
1946  MS_BESSEL (name, fcn) \
1947  MM_BESSEL (name, fcn) \
1948  SN_BESSEL (name, fcn) \
1949  NS_BESSEL (name, fcn) \
1950  NN_BESSEL (name, fcn) \
1951  RC_BESSEL (name, fcn)
1952 
1954  ALL_BESSEL (bessely, cbesy)
1955  ALL_BESSEL (besseli, cbesi)
1956  ALL_BESSEL (besselk, cbesk)
1957  ALL_BESSEL (besselh1, cbesh1)
1958  ALL_BESSEL (besselh2, cbesh2)
1959 
1960 #undef ALL_BESSEL
1961 #undef SS_BESSEL
1962 #undef SM_BESSEL
1963 #undef MS_BESSEL
1964 #undef MM_BESSEL
1965 #undef SN_BESSEL
1966 #undef NS_BESSEL
1967 #undef NN_BESSEL
1968 #undef RC_BESSEL
1969 
1970  Complex
1971  airy (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr)
1972  {
1973  double ar = 0.0;
1974  double ai = 0.0;
1975 
1976  octave_idx_type nz;
1977 
1978  double zr = z.real ();
1979  double zi = z.imag ();
1980 
1981  octave_idx_type id = deriv ? 1 : 0;
1982 
1983  F77_FUNC (zairy, ZAIRY) (zr, zi, id, 2, ar, ai, nz, ierr);
1984 
1985  if (! scaled)
1986  {
1987  Complex expz = exp (- 2.0 / 3.0 * z * sqrt (z));
1988 
1989  double rexpz = expz.real ();
1990  double iexpz = expz.imag ();
1991 
1992  double tmp = ar*rexpz - ai*iexpz;
1993 
1994  ai = ar*iexpz + ai*rexpz;
1995  ar = tmp;
1996  }
1997 
1998  if (zi == 0.0 && (! scaled || zr >= 0.0))
1999  ai = 0.0;
2000 
2001  return bessel_return_value (Complex (ar, ai), ierr);
2002  }
2003 
2004  Complex
2005  biry (const Complex& z, bool deriv, bool scaled, octave_idx_type& ierr)
2006  {
2007  double ar = 0.0;
2008  double ai = 0.0;
2009 
2010  double zr = z.real ();
2011  double zi = z.imag ();
2012 
2013  octave_idx_type id = deriv ? 1 : 0;
2014 
2015  F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 2, ar, ai, ierr);
2016 
2017  if (! scaled)
2018  {
2019  Complex expz = exp (std::abs (std::real (2.0 / 3.0 * z * sqrt (z))));
2020 
2021  double rexpz = expz.real ();
2022  double iexpz = expz.imag ();
2023 
2024  double tmp = ar*rexpz - ai*iexpz;
2025 
2026  ai = ar*iexpz + ai*rexpz;
2027  ar = tmp;
2028  }
2029 
2030  if (zi == 0.0 && (! scaled || zr >= 0.0))
2031  ai = 0.0;
2032 
2033  return bessel_return_value (Complex (ar, ai), ierr);
2034  }
2035 
2037  airy (const ComplexMatrix& z, bool deriv, bool scaled,
2039  {
2040  octave_idx_type nr = z.rows ();
2041  octave_idx_type nc = z.cols ();
2042 
2043  ComplexMatrix retval (nr, nc);
2044 
2045  ierr.resize (dim_vector (nr, nc));
2046 
2047  for (octave_idx_type j = 0; j < nc; j++)
2048  for (octave_idx_type i = 0; i < nr; i++)
2049  retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j));
2050 
2051  return retval;
2052  }
2053 
2055  biry (const ComplexMatrix& z, bool deriv, bool scaled,
2057  {
2058  octave_idx_type nr = z.rows ();
2059  octave_idx_type nc = z.cols ();
2060 
2061  ComplexMatrix retval (nr, nc);
2062 
2063  ierr.resize (dim_vector (nr, nc));
2064 
2065  for (octave_idx_type j = 0; j < nc; j++)
2066  for (octave_idx_type i = 0; i < nr; i++)
2067  retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j));
2068 
2069  return retval;
2070  }
2071 
2073  airy (const ComplexNDArray& z, bool deriv, bool scaled,
2075  {
2076  dim_vector dv = z.dims ();
2077  octave_idx_type nel = dv.numel ();
2078  ComplexNDArray retval (dv);
2079 
2080  ierr.resize (dv);
2081 
2082  for (octave_idx_type i = 0; i < nel; i++)
2083  retval(i) = airy (z(i), deriv, scaled, ierr(i));
2084 
2085  return retval;
2086  }
2087 
2089  biry (const ComplexNDArray& z, bool deriv, bool scaled,
2091  {
2092  dim_vector dv = z.dims ();
2093  octave_idx_type nel = dv.numel ();
2094  ComplexNDArray retval (dv);
2095 
2096  ierr.resize (dv);
2097 
2098  for (octave_idx_type i = 0; i < nel; i++)
2099  retval(i) = biry (z(i), deriv, scaled, ierr(i));
2100 
2101  return retval;
2102  }
2103 
2104  FloatComplex
2105  airy (const FloatComplex& z, bool deriv, bool scaled, octave_idx_type& ierr)
2106  {
2107  FloatComplex a;
2108 
2109  octave_idx_type nz;
2110 
2111  octave_idx_type id = deriv ? 1 : 0;
2112 
2113  F77_FUNC (cairy, CAIRY) (F77_CONST_CMPLX_ARG (&z), id, 2, F77_CMPLX_ARG (&a),
2114  nz, ierr);
2115 
2116  float ar = a.real ();
2117  float ai = a.imag ();
2118 
2119  if (! scaled)
2120  {
2121  FloatComplex expz = exp (- 2.0f / 3.0f * z * sqrt (z));
2122 
2123  float rexpz = expz.real ();
2124  float iexpz = expz.imag ();
2125 
2126  float tmp = ar*rexpz - ai*iexpz;
2127 
2128  ai = ar*iexpz + ai*rexpz;
2129  ar = tmp;
2130  }
2131 
2132  if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0))
2133  ai = 0.0;
2134 
2135  return bessel_return_value (FloatComplex (ar, ai), ierr);
2136  }
2137 
2138  FloatComplex
2139  biry (const FloatComplex& z, bool deriv, bool scaled, octave_idx_type& ierr)
2140  {
2141  FloatComplex a;
2142 
2143  octave_idx_type id = deriv ? 1 : 0;
2144 
2145  F77_FUNC (cbiry, CBIRY) (F77_CONST_CMPLX_ARG (&z), id, 2, F77_CMPLX_ARG (&a),
2146  ierr);
2147 
2148  float ar = a.real ();
2149  float ai = a.imag ();
2150 
2151  if (! scaled)
2152  {
2153  FloatComplex expz = exp (std::abs (std::real (2.0f / 3.0f * z * sqrt (z))));
2154 
2155  float rexpz = expz.real ();
2156  float iexpz = expz.imag ();
2157 
2158  float tmp = ar*rexpz - ai*iexpz;
2159 
2160  ai = ar*iexpz + ai*rexpz;
2161  ar = tmp;
2162  }
2163 
2164  if (z.imag () == 0.0 && (! scaled || z.real () >= 0.0))
2165  ai = 0.0;
2166 
2167  return bessel_return_value (FloatComplex (ar, ai), ierr);
2168  }
2169 
2171  airy (const FloatComplexMatrix& z, bool deriv, bool scaled,
2173  {
2174  octave_idx_type nr = z.rows ();
2175  octave_idx_type nc = z.cols ();
2176 
2177  FloatComplexMatrix retval (nr, nc);
2178 
2179  ierr.resize (dim_vector (nr, nc));
2180 
2181  for (octave_idx_type j = 0; j < nc; j++)
2182  for (octave_idx_type i = 0; i < nr; i++)
2183  retval(i,j) = airy (z(i,j), deriv, scaled, ierr(i,j));
2184 
2185  return retval;
2186  }
2187 
2189  biry (const FloatComplexMatrix& z, bool deriv, bool scaled,
2191  {
2192  octave_idx_type nr = z.rows ();
2193  octave_idx_type nc = z.cols ();
2194 
2195  FloatComplexMatrix retval (nr, nc);
2196 
2197  ierr.resize (dim_vector (nr, nc));
2198 
2199  for (octave_idx_type j = 0; j < nc; j++)
2200  for (octave_idx_type i = 0; i < nr; i++)
2201  retval(i,j) = biry (z(i,j), deriv, scaled, ierr(i,j));
2202 
2203  return retval;
2204  }
2205 
2207  airy (const FloatComplexNDArray& z, bool deriv, bool scaled,
2209  {
2210  dim_vector dv = z.dims ();
2211  octave_idx_type nel = dv.numel ();
2213 
2214  ierr.resize (dv);
2215 
2216  for (octave_idx_type i = 0; i < nel; i++)
2217  retval(i) = airy (z(i), deriv, scaled, ierr(i));
2218 
2219  return retval;
2220  }
2221 
2223  biry (const FloatComplexNDArray& z, bool deriv, bool scaled,
2225  {
2226  dim_vector dv = z.dims ();
2227  octave_idx_type nel = dv.numel ();
2229 
2230  ierr.resize (dv);
2231 
2232  for (octave_idx_type i = 0; i < nel; i++)
2233  retval(i) = biry (z(i), deriv, scaled, ierr(i));
2234 
2235  return retval;
2236  }
2237 
2238  OCTAVE_NORETURN static void
2240  const dim_vector& d3)
2241  {
2242  std::string d1_str = d1.str ();
2243  std::string d2_str = d2.str ();
2244  std::string d3_str = d3.str ();
2245 
2246  (*current_liboctave_error_handler)
2247  ("betainc: nonconformant arguments (x is %s, a is %s, b is %s)",
2248  d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
2249  }
2250 
2251  OCTAVE_NORETURN static void
2253  const dim_vector& d3)
2254  {
2255  std::string d1_str = d1.str ();
2256  std::string d2_str = d2.str ();
2257  std::string d3_str = d3.str ();
2258 
2259  (*current_liboctave_error_handler)
2260  ("betaincinv: nonconformant arguments (x is %s, a is %s, b is %s)",
2261  d1_str.c_str (), d2_str.c_str (), d3_str.c_str ());
2262  }
2263 
2264  double
2265  betainc (double x, double a, double b)
2266  {
2267  double retval;
2268  F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval));
2269  return retval;
2270  }
2271 
2273  betainc (double x, double a, const Array<double>& b)
2274  {
2275  dim_vector dv = b.dims ();
2276  octave_idx_type nel = dv.numel ();
2277 
2278  Array<double> retval (dv);
2279 
2280  double *pretval = retval.fortran_vec ();
2281 
2282  for (octave_idx_type i = 0; i < nel; i++)
2283  *pretval++ = betainc (x, a, b(i));
2284 
2285  return retval;
2286  }
2287 
2289  betainc (double x, const Array<double>& a, double b)
2290  {
2291  dim_vector dv = a.dims ();
2292  octave_idx_type nel = dv.numel ();
2293 
2294  Array<double> retval (dv);
2295 
2296  double *pretval = retval.fortran_vec ();
2297 
2298  for (octave_idx_type i = 0; i < nel; i++)
2299  *pretval++ = betainc (x, a(i), b);
2300 
2301  return retval;
2302  }
2303 
2305  betainc (double x, const Array<double>& a, const Array<double>& b)
2306  {
2308  dim_vector dv = a.dims ();
2309 
2310  if (dv != b.dims ())
2311  err_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ());
2312 
2313  octave_idx_type nel = dv.numel ();
2314 
2315  retval.resize (dv);
2316 
2317  double *pretval = retval.fortran_vec ();
2318 
2319  for (octave_idx_type i = 0; i < nel; i++)
2320  *pretval++ = betainc (x, a(i), b(i));
2321 
2322  return retval;
2323  }
2324 
2326  betainc (const Array<double>& x, double a, double b)
2327  {
2328  dim_vector dv = x.dims ();
2329  octave_idx_type nel = dv.numel ();
2330 
2331  Array<double> retval (dv);
2332 
2333  double *pretval = retval.fortran_vec ();
2334 
2335  for (octave_idx_type i = 0; i < nel; i++)
2336  *pretval++ = betainc (x(i), a, b);
2337 
2338  return retval;
2339  }
2340 
2342  betainc (const Array<double>& x, double a, const Array<double>& b)
2343  {
2345  dim_vector dv = x.dims ();
2346 
2347  if (dv != b.dims ())
2348  err_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ());
2349 
2350  octave_idx_type nel = dv.numel ();
2351 
2352  retval.resize (dv);
2353 
2354  double *pretval = retval.fortran_vec ();
2355 
2356  for (octave_idx_type i = 0; i < nel; i++)
2357  *pretval++ = betainc (x(i), a, b(i));
2358 
2359  return retval;
2360  }
2361 
2363  betainc (const Array<double>& x, const Array<double>& a, double b)
2364  {
2366  dim_vector dv = x.dims ();
2367 
2368  if (dv != a.dims ())
2369  err_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0));
2370 
2371  octave_idx_type nel = dv.numel ();
2372 
2373  retval.resize (dv);
2374 
2375  double *pretval = retval.fortran_vec ();
2376 
2377  for (octave_idx_type i = 0; i < nel; i++)
2378  *pretval++ = betainc (x(i), a(i), b);
2379 
2380  return retval;
2381  }
2382 
2385  {
2387  dim_vector dv = x.dims ();
2388 
2389  if (dv != a.dims () || dv != b.dims ())
2390  err_betainc_nonconformant (dv, a.dims (), b.dims ());
2391 
2392  octave_idx_type nel = dv.numel ();
2393 
2394  retval.resize (dv);
2395 
2396  double *pretval = retval.fortran_vec ();
2397 
2398  for (octave_idx_type i = 0; i < nel; i++)
2399  *pretval++ = betainc (x(i), a(i), b(i));
2400 
2401  return retval;
2402  }
2403 
2404  float
2405  betainc (float x, float a, float b)
2406  {
2407  float retval;
2408  F77_XFCN (xbetai, XBETAI, (x, a, b, retval));
2409  return retval;
2410  }
2411 
2412  Array<float>
2413  betainc (float x, float a, const Array<float>& b)
2414  {
2415  dim_vector dv = b.dims ();
2416  octave_idx_type nel = dv.numel ();
2417 
2418  Array<float> retval (dv);
2419 
2420  float *pretval = retval.fortran_vec ();
2421 
2422  for (octave_idx_type i = 0; i < nel; i++)
2423  *pretval++ = betainc (x, a, b(i));
2424 
2425  return retval;
2426  }
2427 
2428  Array<float>
2429  betainc (float x, const Array<float>& a, float b)
2430  {
2431  dim_vector dv = a.dims ();
2432  octave_idx_type nel = dv.numel ();
2433 
2434  Array<float> retval (dv);
2435 
2436  float *pretval = retval.fortran_vec ();
2437 
2438  for (octave_idx_type i = 0; i < nel; i++)
2439  *pretval++ = betainc (x, a(i), b);
2440 
2441  return retval;
2442  }
2443 
2444  Array<float>
2445  betainc (float x, const Array<float>& a, const Array<float>& b)
2446  {
2448  dim_vector dv = a.dims ();
2449 
2450  if (dv != b.dims ())
2451  err_betainc_nonconformant (dim_vector (0, 0), dv, b.dims ());
2452 
2453  octave_idx_type nel = dv.numel ();
2454 
2455  retval.resize (dv);
2456 
2457  float *pretval = retval.fortran_vec ();
2458 
2459  for (octave_idx_type i = 0; i < nel; i++)
2460  *pretval++ = betainc (x, a(i), b(i));
2461 
2462  return retval;
2463  }
2464 
2465  Array<float>
2466  betainc (const Array<float>& x, float a, float b)
2467  {
2468  dim_vector dv = x.dims ();
2469  octave_idx_type nel = dv.numel ();
2470 
2471  Array<float> retval (dv);
2472 
2473  float *pretval = retval.fortran_vec ();
2474 
2475  for (octave_idx_type i = 0; i < nel; i++)
2476  *pretval++ = betainc (x(i), a, b);
2477 
2478  return retval;
2479  }
2480 
2481  Array<float>
2482  betainc (const Array<float>& x, float a, const Array<float>& b)
2483  {
2485  dim_vector dv = x.dims ();
2486 
2487  if (dv != b.dims ())
2488  err_betainc_nonconformant (dv, dim_vector (0, 0), b.dims ());
2489 
2490  octave_idx_type nel = dv.numel ();
2491 
2492  retval.resize (dv);
2493 
2494  float *pretval = retval.fortran_vec ();
2495 
2496  for (octave_idx_type i = 0; i < nel; i++)
2497  *pretval++ = betainc (x(i), a, b(i));
2498 
2499  return retval;
2500  }
2501 
2502  Array<float>
2503  betainc (const Array<float>& x, const Array<float>& a, float b)
2504  {
2506  dim_vector dv = x.dims ();
2507 
2508  if (dv != a.dims ())
2509  err_betainc_nonconformant (dv, a.dims (), dim_vector (0, 0));
2510 
2511  octave_idx_type nel = dv.numel ();
2512 
2513  retval.resize (dv);
2514 
2515  float *pretval = retval.fortran_vec ();
2516 
2517  for (octave_idx_type i = 0; i < nel; i++)
2518  *pretval++ = betainc (x(i), a(i), b);
2519 
2520  return retval;
2521  }
2522 
2523  Array<float>
2524  betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b)
2525  {
2527  dim_vector dv = x.dims ();
2528 
2529  if (dv != a.dims () || dv != b.dims ())
2530  err_betainc_nonconformant (dv, a.dims (), b.dims ());
2531 
2532  octave_idx_type nel = dv.numel ();
2533 
2534  retval.resize (dv);
2535 
2536  float *pretval = retval.fortran_vec ();
2537 
2538  for (octave_idx_type i = 0; i < nel; i++)
2539  *pretval++ = betainc (x(i), a(i), b(i));
2540 
2541  return retval;
2542  }
2543 
2544  // FIXME: there is still room for improvement here...
2545 
2546  double
2547  gammainc (double x, double a, bool& err)
2548  {
2549  if (a < 0.0 || x < 0.0)
2550  (*current_liboctave_error_handler)
2551  ("gammainc: A and X must be non-negative");
2552 
2553  err = false;
2554 
2555  double retval;
2556 
2557  F77_XFCN (xgammainc, XGAMMAINC, (a, x, retval));
2558 
2559  return retval;
2560  }
2561 
2562  Matrix
2563  gammainc (double x, const Matrix& a)
2564  {
2565  octave_idx_type nr = a.rows ();
2566  octave_idx_type nc = a.cols ();
2567 
2568  Matrix retval (nr, nc);
2569 
2570  bool err;
2571 
2572  for (octave_idx_type j = 0; j < nc; j++)
2573  for (octave_idx_type i = 0; i < nr; i++)
2574  {
2575  retval(i,j) = gammainc (x, a(i,j), err);
2576 
2577  if (err)
2578  return Matrix ();
2579  }
2580 
2581  return retval;
2582  }
2583 
2584  Matrix
2585  gammainc (const Matrix& x, double a)
2586  {
2587  octave_idx_type nr = x.rows ();
2588  octave_idx_type nc = x.cols ();
2589 
2590  Matrix retval (nr, nc);
2591 
2592  bool err;
2593 
2594  for (octave_idx_type j = 0; j < nc; j++)
2595  for (octave_idx_type i = 0; i < nr; i++)
2596  {
2597  retval(i,j) = gammainc (x(i,j), a, err);
2598 
2599  if (err)
2600  return Matrix ();
2601  }
2602 
2603  return retval;
2604  }
2605 
2606  Matrix
2607  gammainc (const Matrix& x, const Matrix& a)
2608  {
2609  octave_idx_type nr = x.rows ();
2610  octave_idx_type nc = x.cols ();
2611 
2612  octave_idx_type a_nr = a.rows ();
2613  octave_idx_type a_nc = a.cols ();
2614 
2615  if (nr != a_nr || nc != a_nc)
2616  (*current_liboctave_error_handler)
2617  ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
2618  nr, nc, a_nr, a_nc);
2619 
2620  Matrix retval (nr, nc);
2621 
2622  bool err;
2623 
2624  for (octave_idx_type j = 0; j < nc; j++)
2625  for (octave_idx_type i = 0; i < nr; i++)
2626  {
2627  retval(i,j) = gammainc (x(i,j), a(i,j), err);
2628 
2629  if (err)
2630  return Matrix ();
2631  }
2632 
2633  return retval;
2634  }
2635 
2636  NDArray
2637  gammainc (double x, const NDArray& a)
2638  {
2639  dim_vector dv = a.dims ();
2640  octave_idx_type nel = dv.numel ();
2641 
2642  NDArray retval (dv);
2643 
2644  bool err;
2645 
2646  for (octave_idx_type i = 0; i < nel; i++)
2647  {
2648  retval(i) = gammainc (x, a(i), err);
2649 
2650  if (err)
2651  return NDArray ();
2652  }
2653 
2654  return retval;
2655  }
2656 
2657  NDArray
2658  gammainc (const NDArray& x, double a)
2659  {
2660  dim_vector dv = x.dims ();
2661  octave_idx_type nel = dv.numel ();
2662 
2663  NDArray retval (dv);
2664 
2665  bool err;
2666 
2667  for (octave_idx_type i = 0; i < nel; i++)
2668  {
2669  retval(i) = gammainc (x(i), a, err);
2670 
2671  if (err)
2672  return NDArray ();
2673  }
2674 
2675  return retval;
2676  }
2677 
2678  NDArray
2679  gammainc (const NDArray& x, const NDArray& a)
2680  {
2681  dim_vector dv = x.dims ();
2682  octave_idx_type nel = dv.numel ();
2683 
2684  if (dv != a.dims ())
2685  {
2686  std::string x_str = dv.str ();
2687  std::string a_str = a.dims ().str ();
2688 
2689  (*current_liboctave_error_handler)
2690  ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
2691  x_str.c_str (), a_str.c_str ());
2692  }
2693 
2694  NDArray retval (dv);
2695 
2696  bool err;
2697 
2698  for (octave_idx_type i = 0; i < nel; i++)
2699  {
2700  retval(i) = gammainc (x(i), a(i), err);
2701 
2702  if (err)
2703  return NDArray ();
2704  }
2705 
2706  return retval;
2707  }
2708 
2709  float
2710  gammainc (float x, float a, bool& err)
2711  {
2712  if (a < 0.0 || x < 0.0)
2713  (*current_liboctave_error_handler)
2714  ("gammainc: A and X must be non-negative");
2715 
2716  err = false;
2717 
2718  float retval;
2719 
2720  F77_XFCN (xsgammainc, XSGAMMAINC, (a, x, retval));
2721 
2722  return retval;
2723  }
2724 
2725  FloatMatrix
2726  gammainc (float x, const FloatMatrix& a)
2727  {
2728  octave_idx_type nr = a.rows ();
2729  octave_idx_type nc = a.cols ();
2730 
2731  FloatMatrix retval (nr, nc);
2732 
2733  bool err;
2734 
2735  for (octave_idx_type j = 0; j < nc; j++)
2736  for (octave_idx_type i = 0; i < nr; i++)
2737  {
2738  retval(i,j) = gammainc (x, a(i,j), err);
2739 
2740  if (err)
2741  return FloatMatrix ();
2742  }
2743 
2744  return retval;
2745  }
2746 
2747  FloatMatrix
2748  gammainc (const FloatMatrix& x, float a)
2749  {
2750  octave_idx_type nr = x.rows ();
2751  octave_idx_type nc = x.cols ();
2752 
2753  FloatMatrix retval (nr, nc);
2754 
2755  bool err;
2756 
2757  for (octave_idx_type j = 0; j < nc; j++)
2758  for (octave_idx_type i = 0; i < nr; i++)
2759  {
2760  retval(i,j) = gammainc (x(i,j), a, err);
2761 
2762  if (err)
2763  return FloatMatrix ();
2764  }
2765 
2766  return retval;
2767  }
2768 
2769  FloatMatrix
2771  {
2772  octave_idx_type nr = x.rows ();
2773  octave_idx_type nc = x.cols ();
2774 
2775  octave_idx_type a_nr = a.rows ();
2776  octave_idx_type a_nc = a.cols ();
2777 
2778  if (nr != a_nr || nc != a_nc)
2779  (*current_liboctave_error_handler)
2780  ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)",
2781  nr, nc, a_nr, a_nc);
2782 
2783  FloatMatrix retval (nr, nc);
2784 
2785  bool err;
2786 
2787  for (octave_idx_type j = 0; j < nc; j++)
2788  for (octave_idx_type i = 0; i < nr; i++)
2789  {
2790  retval(i,j) = gammainc (x(i,j), a(i,j), err);
2791 
2792  if (err)
2793  return FloatMatrix ();
2794  }
2795 
2796  return retval;
2797  }
2798 
2799  FloatNDArray
2800  gammainc (float x, const FloatNDArray& a)
2801  {
2802  dim_vector dv = a.dims ();
2803  octave_idx_type nel = dv.numel ();
2804 
2805  FloatNDArray retval (dv);
2806 
2807  bool err;
2808 
2809  for (octave_idx_type i = 0; i < nel; i++)
2810  {
2811  retval(i) = gammainc (x, a(i), err);
2812 
2813  if (err)
2814  return FloatNDArray ();
2815  }
2816 
2817  return retval;
2818  }
2819 
2820  FloatNDArray
2821  gammainc (const FloatNDArray& x, float a)
2822  {
2823  dim_vector dv = x.dims ();
2824  octave_idx_type nel = dv.numel ();
2825 
2826  FloatNDArray retval (dv);
2827 
2828  bool err;
2829 
2830  for (octave_idx_type i = 0; i < nel; i++)
2831  {
2832  retval(i) = gammainc (x(i), a, err);
2833 
2834  if (err)
2835  return FloatNDArray ();
2836  }
2837 
2838  return retval;
2839  }
2840 
2841  FloatNDArray
2843  {
2844  dim_vector dv = x.dims ();
2845  octave_idx_type nel = dv.numel ();
2846 
2847  if (dv != a.dims ())
2848  {
2849  std::string x_str = dv.str ();
2850  std::string a_str = a.dims ().str ();
2851 
2852  (*current_liboctave_error_handler)
2853  ("gammainc: nonconformant arguments (arg 1 is %s, arg 2 is %s)",
2854  x_str.c_str (), a_str.c_str ());
2855  }
2856 
2857  FloatNDArray retval (dv);
2858 
2859  bool err;
2860 
2861  for (octave_idx_type i = 0; i < nel; i++)
2862  {
2863  retval(i) = gammainc (x(i), a(i), err);
2864 
2865  if (err)
2866  return FloatNDArray ();
2867  }
2868 
2869  return retval;
2870  }
2871 
2872  Complex rc_log1p (double x)
2873  {
2874  const double pi = 3.14159265358979323846;
2875  return (x < -1.0
2876  ? Complex (std::log (-(1.0 + x)), pi)
2877  : Complex (log1p (x)));
2878  }
2879 
2880  FloatComplex rc_log1p (float x)
2881  {
2882  const float pi = 3.14159265358979323846f;
2883  return (x < -1.0f
2884  ? FloatComplex (std::log (-(1.0f + x)), pi)
2885  : FloatComplex (log1p (x)));
2886  }
2887 
2888  // This algorithm is due to P. J. Acklam.
2889  //
2890  // See http://home.online.no/~pjacklam/notes/invnorm/
2891  //
2892  // The rational approximation has relative accuracy 1.15e-9 in the whole
2893  // region. For doubles, it is refined by a single step of Halley's 3rd
2894  // order method. For single precision, the accuracy is already OK, so
2895  // we skip it to get faster evaluation.
2896 
2897  static double do_erfinv (double x, bool refine)
2898  {
2899  // Coefficients of rational approximation.
2900  static const double a[] =
2901  {
2902  -2.806989788730439e+01, 1.562324844726888e+02,
2903  -1.951109208597547e+02, 9.783370457507161e+01,
2904  -2.168328665628878e+01, 1.772453852905383e+00
2905  };
2906  static const double b[] =
2907  {
2908  -5.447609879822406e+01, 1.615858368580409e+02,
2909  -1.556989798598866e+02, 6.680131188771972e+01,
2910  -1.328068155288572e+01
2911  };
2912  static const double c[] =
2913  {
2914  -5.504751339936943e-03, -2.279687217114118e-01,
2915  -1.697592457770869e+00, -1.802933168781950e+00,
2916  3.093354679843505e+00, 2.077595676404383e+00
2917  };
2918  static const double d[] =
2919  {
2920  7.784695709041462e-03, 3.224671290700398e-01,
2921  2.445134137142996e+00, 3.754408661907416e+00
2922  };
2923 
2924  static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2.
2925  static const double pbreak = 0.95150;
2926  double ax = fabs (x), y;
2927 
2928  // Select case.
2929  if (ax <= pbreak)
2930  {
2931  // Middle region.
2932  const double q = 0.5 * x, r = q*q;
2933  const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
2934  const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
2935  y = yn / yd;
2936  }
2937  else if (ax < 1.0)
2938  {
2939  // Tail region.
2940  const double q = sqrt (-2*std::log (0.5*(1-ax)));
2941  const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
2942  const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
2943  y = yn / yd * octave::math::signum (-x);
2944  }
2945  else if (ax == 1.0)
2947  else
2949 
2950  if (refine)
2951  {
2952  // One iteration of Halley's method gives full precision.
2953  double u = (erf (y) - x) * spi2 * exp (y*y);
2954  y -= u / (1 + y*u);
2955  }
2956 
2957  return y;
2958  }
2959 
2960  double erfinv (double x)
2961  {
2962  return do_erfinv (x, true);
2963  }
2964 
2965  float erfinv (float x)
2966  {
2967  return do_erfinv (x, false);
2968  }
2969 
2970  // The algorthim for erfcinv is an adaptation of the erfinv algorithm
2971  // above from P. J. Acklam. It has been modified to run over the
2972  // different input domain of erfcinv. See the notes for erfinv for an
2973  // explanation.
2974 
2975  static double do_erfcinv (double x, bool refine)
2976  {
2977  // Coefficients of rational approximation.
2978  static const double a[] =
2979  {
2980  -2.806989788730439e+01, 1.562324844726888e+02,
2981  -1.951109208597547e+02, 9.783370457507161e+01,
2982  -2.168328665628878e+01, 1.772453852905383e+00
2983  };
2984  static const double b[] =
2985  {
2986  -5.447609879822406e+01, 1.615858368580409e+02,
2987  -1.556989798598866e+02, 6.680131188771972e+01,
2988  -1.328068155288572e+01
2989  };
2990  static const double c[] =
2991  {
2992  -5.504751339936943e-03, -2.279687217114118e-01,
2993  -1.697592457770869e+00, -1.802933168781950e+00,
2994  3.093354679843505e+00, 2.077595676404383e+00
2995  };
2996  static const double d[] =
2997  {
2998  7.784695709041462e-03, 3.224671290700398e-01,
2999  2.445134137142996e+00, 3.754408661907416e+00
3000  };
3001 
3002  static const double spi2 = 8.862269254527579e-01; // sqrt(pi)/2.
3003  static const double pbreak_lo = 0.04850; // 1-pbreak
3004  static const double pbreak_hi = 1.95150; // 1+pbreak
3005  double y;
3006 
3007  // Select case.
3008  if (x >= pbreak_lo && x <= pbreak_hi)
3009  {
3010  // Middle region.
3011  const double q = 0.5*(1-x), r = q*q;
3012  const double yn = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q;
3013  const double yd = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0;
3014  y = yn / yd;
3015  }
3016  else if (x > 0.0 && x < 2.0)
3017  {
3018  // Tail region.
3019  const double q = (x < 1
3020  ? sqrt (-2*std::log (0.5*x))
3021  : sqrt (-2*std::log (0.5*(2-x))));
3022 
3023  const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
3024 
3025  const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
3026 
3027  y = yn / yd;
3028 
3029  if (x < pbreak_lo)
3030  y = -y;
3031  }
3032  else if (x == 0.0)
3034  else if (x == 2.0)
3036  else
3038 
3039  if (refine)
3040  {
3041  // One iteration of Halley's method gives full precision.
3042  double u = (erf (y) - (1-x)) * spi2 * exp (y*y);
3043  y -= u / (1 + y*u);
3044  }
3045 
3046  return y;
3047  }
3048 
3049  double erfcinv (double x)
3050  {
3051  return do_erfcinv (x, true);
3052  }
3053 
3054  float erfcinv (float x)
3055  {
3056  return do_erfcinv (x, false);
3057  }
3058 
3059  //
3060  // Incomplete Beta function ratio
3061  //
3062  // Algorithm based on the one by John Burkardt.
3063  // See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html
3064  //
3065  // The original code is distributed under the GNU LGPL v3 license.
3066  //
3067  // Reference:
3068  //
3069  // KL Majumder, GP Bhattacharjee,
3070  // Algorithm AS 63:
3071  // The incomplete Beta Integral,
3072  // Applied Statistics,
3073  // Volume 22, Number 3, 1973, pages 409-411.
3074  //
3075  double
3076  betain (double x, double p, double q, double beta, bool& err)
3077  {
3078  double acu = 0.1E-14, ai, cx;
3079  bool indx;
3080  int ns;
3081  double pp, psq, qq, rx, temp, term, value, xx;
3082 
3083  value = x;
3084  err = false;
3085 
3086  // Check the input arguments.
3087 
3088  if ((p <= 0.0 || q <= 0.0) || (x < 0.0 || 1.0 < x))
3089  {
3090  err = true;
3091  return value;
3092  }
3093 
3094  // Special cases.
3095 
3096  if (x == 0.0 || x == 1.0)
3097  {
3098  return value;
3099  }
3100 
3101  // Change tail if necessary and determine S.
3102 
3103  psq = p + q;
3104  cx = 1.0 - x;
3105 
3106  if (p < psq * x)
3107  {
3108  xx = cx;
3109  cx = x;
3110  pp = q;
3111  qq = p;
3112  indx = true;
3113  }
3114  else
3115  {
3116  xx = x;
3117  pp = p;
3118  qq = q;
3119  indx = false;
3120  }
3121 
3122  term = 1.0;
3123  ai = 1.0;
3124  value = 1.0;
3125  ns = static_cast<int> (qq + cx * psq);
3126 
3127  // Use the Soper reduction formula.
3128 
3129  rx = xx / cx;
3130  temp = qq - ai;
3131  if (ns == 0)
3132  {
3133  rx = xx;
3134  }
3135 
3136  for ( ; ; )
3137  {
3138  term *= temp * rx / (pp + ai);
3139  value += term;
3140  temp = fabs (term);
3141 
3142  if (temp <= acu && temp <= acu * value)
3143  {
3144  value *= exp (pp * std::log (xx)
3145  + (qq - 1.0) * std::log (cx) - beta) / pp;
3146 
3147  if (indx)
3148  {
3149  value = 1.0 - value;
3150  }
3151  break;
3152  }
3153 
3154  ai += 1.0;
3155  ns -= 1;
3156 
3157  if (0 <= ns)
3158  {
3159  temp = qq - ai;
3160  if (ns == 0)
3161  {
3162  rx = xx;
3163  }
3164  }
3165  else
3166  {
3167  temp = psq;
3168  psq += 1.0;
3169  }
3170  }
3171 
3172  return value;
3173  }
3174 
3175  //
3176  // Inverse of the incomplete Beta function
3177  //
3178  // Algorithm based on the one by John Burkardt.
3179  // See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html
3180  //
3181  // The original code is distributed under the GNU LGPL v3 license.
3182  //
3183  // Reference:
3184  //
3185  // GW Cran, KJ Martin, GE Thomas,
3186  // Remark AS R19 and Algorithm AS 109:
3187  // A Remark on Algorithms AS 63: The Incomplete Beta Integral
3188  // and AS 64: Inverse of the Incomplete Beta Integeral,
3189  // Applied Statistics,
3190  // Volume 26, Number 1, 1977, pages 111-114.
3191  //
3192  double
3193  betaincinv (double y, double p, double q)
3194  {
3195  double a, acu, adj, fpu, g, h;
3196  int iex;
3197  bool indx;
3198  double pp, prev, qq, r, s, sae = -37.0, sq, t, tx, value, w, xin, ycur, yprev;
3199 
3200  double beta = lgamma (p) + lgamma (q) - lgamma (p + q);
3201  bool err = false;
3202  fpu = pow (10.0, sae);
3203  value = y;
3204 
3205  // Test for admissibility of parameters.
3206 
3207  if (p <= 0.0 || q <= 0.0)
3208  (*current_liboctave_error_handler) ("betaincinv: wrong parameters");
3209  if (y < 0.0 || 1.0 < y)
3210  (*current_liboctave_error_handler) ("betaincinv: wrong parameter Y");
3211 
3212  if (y == 0.0 || y == 1.0)
3213  return value;
3214 
3215  // Change tail if necessary.
3216 
3217  if (0.5 < y)
3218  {
3219  a = 1.0 - y;
3220  pp = q;
3221  qq = p;
3222  indx = true;
3223  }
3224  else
3225  {
3226  a = y;
3227  pp = p;
3228  qq = q;
3229  indx = false;
3230  }
3231 
3232  // Calculate the initial approximation.
3233 
3234  r = sqrt (- std::log (a * a));
3235 
3236  ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r);
3237 
3238  if (1.0 < pp && 1.0 < qq)
3239  {
3240  r = (ycur * ycur - 3.0) / 6.0;
3241  s = 1.0 / (pp + pp - 1.0);
3242  t = 1.0 / (qq + qq - 1.0);
3243  h = 2.0 / (s + t);
3244  w = ycur * sqrt (h + r) / h - (t - s) * (r + 5.0 / 6.0 - 2.0 / (3.0 * h));
3245  value = pp / (pp + qq * exp (w + w));
3246  }
3247  else
3248  {
3249  r = qq + qq;
3250  t = 1.0 / (9.0 * qq);
3251  t = r * pow (1.0 - t + ycur * sqrt (t), 3);
3252 
3253  if (t <= 0.0)
3254  {
3255  value = 1.0 - exp ((std::log ((1.0 - a) * qq) + beta) / qq);
3256  }
3257  else
3258  {
3259  t = (4.0 * pp + r - 2.0) / t;
3260 
3261  if (t <= 1.0)
3262  {
3263  value = exp ((std::log (a * pp) + beta) / pp);
3264  }
3265  else
3266  {
3267  value = 1.0 - 2.0 / (t + 1.0);
3268  }
3269  }
3270  }
3271 
3272  // Solve for X by a modified Newton-Raphson method,
3273  // using the function BETAIN.
3274 
3275  r = 1.0 - pp;
3276  t = 1.0 - qq;
3277  yprev = 0.0;
3278  sq = 1.0;
3279  prev = 1.0;
3280 
3281  if (value < 0.0001)
3282  {
3283  value = 0.0001;
3284  }
3285 
3286  if (0.9999 < value)
3287  {
3288  value = 0.9999;
3289  }
3290 
3291  iex = std::max (- 5.0 / pp / pp - 1.0 / pow (a, 0.2) - 13.0, sae);
3292 
3293  acu = pow (10.0, iex);
3294 
3295  for ( ; ; )
3296  {
3297  ycur = betain (value, pp, qq, beta, err);
3298 
3299  if (err)
3300  {
3301  return value;
3302  }
3303 
3304  xin = value;
3305  ycur = (ycur - a) * exp (beta + r * std::log (xin)
3306  + t * std::log (1.0 - xin));
3307 
3308  if (ycur * yprev <= 0.0)
3309  {
3310  prev = std::max (sq, fpu);
3311  }
3312 
3313  g = 1.0;
3314 
3315  for ( ; ; )
3316  {
3317  for ( ; ; )
3318  {
3319  adj = g * ycur;
3320  sq = adj * adj;
3321 
3322  if (sq < prev)
3323  {
3324  tx = value - adj;
3325 
3326  if (0.0 <= tx && tx <= 1.0)
3327  {
3328  break;
3329  }
3330  }
3331  g /= 3.0;
3332  }
3333 
3334  if (prev <= acu)
3335  {
3336  if (indx)
3337  {
3338  value = 1.0 - value;
3339  }
3340  return value;
3341  }
3342 
3343  if (ycur * ycur <= acu)
3344  {
3345  if (indx)
3346  {
3347  value = 1.0 - value;
3348  }
3349  return value;
3350  }
3351 
3352  if (tx != 0.0 && tx != 1.0)
3353  {
3354  break;
3355  }
3356 
3357  g /= 3.0;
3358  }
3359 
3360  if (tx == value)
3361  {
3362  break;
3363  }
3364 
3365  value = tx;
3366  yprev = ycur;
3367  }
3368 
3369  if (indx)
3370  {
3371  value = 1.0 - value;
3372  }
3373 
3374  return value;
3375  }
3376 
3378  betaincinv (double x, double a, const Array<double>& b)
3379  {
3380  dim_vector dv = b.dims ();
3381  octave_idx_type nel = dv.numel ();
3382 
3383  Array<double> retval (dv);
3384 
3385  double *pretval = retval.fortran_vec ();
3386 
3387  for (octave_idx_type i = 0; i < nel; i++)
3388  *pretval++ = betaincinv (x, a, b(i));
3389 
3390  return retval;
3391  }
3392 
3394  betaincinv (double x, const Array<double>& a, double b)
3395  {
3396  dim_vector dv = a.dims ();
3397  octave_idx_type nel = dv.numel ();
3398 
3399  Array<double> retval (dv);
3400 
3401  double *pretval = retval.fortran_vec ();
3402 
3403  for (octave_idx_type i = 0; i < nel; i++)
3404  *pretval++ = betaincinv (x, a(i), b);
3405 
3406  return retval;
3407  }
3408 
3410  betaincinv (double x, const Array<double>& a, const Array<double>& b)
3411  {
3413  dim_vector dv = a.dims ();
3414 
3415  if (dv != b.dims ())
3417 
3418  octave_idx_type nel = dv.numel ();
3419 
3420  retval.resize (dv);
3421 
3422  double *pretval = retval.fortran_vec ();
3423 
3424  for (octave_idx_type i = 0; i < nel; i++)
3425  *pretval++ = betaincinv (x, a(i), b(i));
3426 
3427  return retval;
3428  }
3429 
3431  betaincinv (const Array<double>& x, double a, double b)
3432  {
3433  dim_vector dv = x.dims ();
3434  octave_idx_type nel = dv.numel ();
3435 
3436  Array<double> retval (dv);
3437 
3438  double *pretval = retval.fortran_vec ();
3439 
3440  for (octave_idx_type i = 0; i < nel; i++)
3441  *pretval++ = betaincinv (x(i), a, b);
3442 
3443  return retval;
3444  }
3445 
3447  betaincinv (const Array<double>& x, double a, const Array<double>& b)
3448  {
3450  dim_vector dv = x.dims ();
3451 
3452  if (dv != b.dims ())
3453  err_betaincinv_nonconformant (dv, dim_vector (0, 0), b.dims ());
3454 
3455  octave_idx_type nel = dv.numel ();
3456 
3457  retval.resize (dv);
3458 
3459  double *pretval = retval.fortran_vec ();
3460 
3461  for (octave_idx_type i = 0; i < nel; i++)
3462  *pretval++ = betaincinv (x(i), a, b(i));
3463 
3464  return retval;
3465  }
3466 
3468  betaincinv (const Array<double>& x, const Array<double>& a, double b)
3469  {
3471  dim_vector dv = x.dims ();
3472 
3473  if (dv != a.dims ())
3474  err_betaincinv_nonconformant (dv, a.dims (), dim_vector (0, 0));
3475 
3476  octave_idx_type nel = dv.numel ();
3477 
3478  retval.resize (dv);
3479 
3480  double *pretval = retval.fortran_vec ();
3481 
3482  for (octave_idx_type i = 0; i < nel; i++)
3483  *pretval++ = betaincinv (x(i), a(i), b);
3484 
3485  return retval;
3486  }
3487 
3490  const Array<double>& b)
3491  {
3493  dim_vector dv = x.dims ();
3494 
3495  if (dv != a.dims () && dv != b.dims ())
3496  err_betaincinv_nonconformant (dv, a.dims (), b.dims ());
3497 
3498  octave_idx_type nel = dv.numel ();
3499 
3500  retval.resize (dv);
3501 
3502  double *pretval = retval.fortran_vec ();
3503 
3504  for (octave_idx_type i = 0; i < nel; i++)
3505  *pretval++ = betaincinv (x(i), a(i), b(i));
3506 
3507  return retval;
3508  }
3509 
3510  void
3511  ellipj (double u, double m, double& sn, double& cn, double& dn, double& err)
3512  {
3513  static const int Nmax = 16;
3514  double m1, t=0, si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi;
3515  int n, Nn, ii;
3516 
3517  if (m < 0 || m > 1)
3518  {
3519  (*current_liboctave_warning_with_id_handler)
3520  ("Octave:ellipj-invalid-m",
3521  "ellipj: invalid M value, required value 0 <= M <= 1");
3522 
3523  sn = cn = dn = lo_ieee_nan_value ();
3524 
3525  return;
3526  }
3527 
3528  double sqrt_eps = sqrt (std::numeric_limits<double>::epsilon ());
3529  if (m < sqrt_eps)
3530  {
3531  // For small m, (Abramowitz and Stegun, Section 16.13)
3532  si_u = sin (u);
3533  co_u = cos (u);
3534  t = 0.25*m*(u - si_u*co_u);
3535  sn = si_u - t * co_u;
3536  cn = co_u + t * si_u;
3537  dn = 1 - 0.5*m*si_u*si_u;
3538  }
3539  else if ((1 - m) < sqrt_eps)
3540  {
3541  // For m1 = (1-m) small (Abramowitz and Stegun, Section 16.15)
3542  m1 = 1 - m;
3543  si_u = sinh (u);
3544  co_u = cosh (u);
3545  ta_u = tanh (u);
3546  se_u = 1/co_u;
3547  sn = ta_u + 0.25*m1*(si_u*co_u - u)*se_u*se_u;
3548  cn = se_u - 0.25*m1*(si_u*co_u - u)*ta_u*se_u;
3549  dn = se_u + 0.25*m1*(si_u*co_u + u)*ta_u*se_u;
3550  }
3551  else
3552  {
3553  // Arithmetic-Geometric Mean (AGM) algorithm
3554  // (Abramowitz and Stegun, Section 16.4)
3555  a[0] = 1;
3556  b = sqrt (1 - m);
3557  c[0] = sqrt (m);
3558  for (n = 1; n < Nmax; ++n)
3559  {
3560  a[n] = (a[n - 1] + b)/2;
3561  c[n] = (a[n - 1] - b)/2;
3562  b = sqrt (a[n - 1]*b);
3563  if (c[n]/a[n] < std::numeric_limits<double>::epsilon ()) break;
3564  }
3565  if (n >= Nmax - 1)
3566  {
3567  err = 1;
3568  return;
3569  }
3570  Nn = n;
3571  for (ii = 1; n > 0; ii = ii*2, --n) ; // ii = pow(2,Nn)
3572  phi = ii*a[Nn]*u;
3573  for (n = Nn; n > 0; --n)
3574  {
3575  phi = (asin ((c[n]/a[n])* sin (phi)) + phi)/2;
3576  }
3577  sn = sin (phi);
3578  cn = cos (phi);
3579  dn = sqrt (1 - m*sn*sn);
3580  }
3581  }
3582 
3583  void
3584  ellipj (const Complex& u, double m, Complex& sn, Complex& cn, Complex& dn,
3585  double& err)
3586  {
3587  double m1 = 1 - m, ss1, cc1, dd1;
3588 
3589  ellipj (u.imag (), m1, ss1, cc1, dd1, err);
3590  if (u.real () == 0)
3591  {
3592  // u is pure imag: Jacoby imag. transf.
3593  sn = Complex (0, ss1/cc1);
3594  cn = 1/cc1; // cn.imag = 0;
3595  dn = dd1/cc1; // dn.imag = 0;
3596  }
3597  else
3598  {
3599  // u is generic complex
3600  double ss, cc, dd, ddd;
3601 
3602  ellipj (u.real (), m, ss, cc, dd, err);
3603  ddd = cc1*cc1 + m*ss*ss*ss1*ss1;
3604  sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd);
3605  cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd);
3606  dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
3607  }
3608  }
3609 
3610  static const double pi = 3.14159265358979323846;
3611 
3612  template <typename T>
3613  static inline T
3614  xlog (const T& x)
3615  {
3616  return log (x);
3617  }
3618 
3619  template <>
3620  inline double
3621  xlog (const double& x)
3622  {
3623  return std::log (x);
3624  }
3625 
3626  template <>
3627  inline float
3628  xlog (const float& x)
3629  {
3630  return std::log (x);
3631  }
3632 
3633  template <typename T>
3634  static T
3636  {
3637  // Coefficients for C.Lanczos expansion of psi function from XLiFE++
3638  // gammaFunctions psi_coef[k] = - (2k+1) * lg_coef[k] (see melina++
3639  // gamma functions -1/12, 3/360,-5/1260, 7/1680,-9/1188,
3640  // 11*691/360360,-13/156, 15*3617/122400, ? , ?
3641  static const T dg_coeff[10] =
3642  {
3643  -0.83333333333333333e-1, 0.83333333333333333e-2,
3644  -0.39682539682539683e-2, 0.41666666666666667e-2,
3645  -0.75757575757575758e-2, 0.21092796092796093e-1,
3646  -0.83333333333333333e-1, 0.4432598039215686,
3647  -0.3053954330270122e+1, 0.125318899521531e+2
3648  };
3649 
3650  T overz2 = T (1.0) / (zc * zc);
3651  T overz2k = overz2;
3652 
3653  T p = 0;
3654  for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2)
3655  p += dg_coeff[k] * overz2k;
3656  p += xlog (zc) - T (0.5) / zc;
3657  return p;
3658  }
3659 
3660  template <typename T>
3661  T
3662  xpsi (T z)
3663  {
3664  static const double euler_mascheroni =
3665  0.577215664901532860606512090082402431042;
3666 
3667  const bool is_int = (octave::math::floor (z) == z);
3668 
3669  T p = 0;
3670  if (z <= 0)
3671  {
3672  // limits - zeros of the gamma function
3673  if (is_int)
3674  p = -octave::numeric_limits<T>::Inf (); // Matlab returns -Inf for psi (0)
3675  else
3676  // Abramowitz and Stegun, page 259, eq 6.3.7
3677  p = psi (1 - z) - (pi / tan (pi * z));
3678  }
3679  else if (is_int)
3680  {
3681  // Abramowitz and Stegun, page 258, eq 6.3.2
3682  p = - euler_mascheroni;
3683  for (octave_idx_type k = z - 1; k > 0; k--)
3684  p += 1.0 / k;
3685  }
3686  else if (octave::math::floor (z + 0.5) == z + 0.5)
3687  {
3688  // Abramowitz and Stegun, page 258, eq 6.3.3 and 6.3.4
3689  for (octave_idx_type k = z; k > 0; k--)
3690  p += 1.0 / (2 * k - 1);
3691 
3692  p = - euler_mascheroni - 2 * std::log (2) + 2 * (p);
3693  }
3694  else
3695  {
3696  // adapted from XLiFE++ gammaFunctions
3697 
3698  T zc = z;
3699  // Use formula for derivative of LogGamma(z)
3700  if (z < 10)
3701  {
3702  const signed char n = 10 - z;
3703  for (signed char k = n - 1; k >= 0; k--)
3704  p -= 1.0 / (k + z);
3705  zc += n;
3706  }
3707  p += lanczos_approximation_psi (zc);
3708  }
3709 
3710  return p;
3711  }
3712 
3713  // explicit instantiations
3714  double psi (double z) { return xpsi (z); }
3715  float psi (float z) { return xpsi (z); }
3716 
3717  template <typename T>
3718  std::complex<T>
3719  xpsi (const std::complex<T>& z)
3720  {
3721  // adapted from XLiFE++ gammaFunctions
3722 
3723  typedef typename std::complex<T>::value_type P;
3724 
3725  P z_r = z.real ();
3726  P z_ra = z_r;
3727 
3728  std::complex<T> dgam (0.0, 0.0);
3729  if (z.imag () == 0)
3730  dgam = std::complex<T> (psi (z_r), 0.0);
3731  else if (z_r < 0)
3732  dgam = psi (P (1.0) - z)- (P (pi) / tan (P (pi) * z));
3733  else
3734  {
3735  // Use formula for derivative of LogGamma(z)
3736  std::complex<T> z_m = z;
3737  if (z_ra < 8)
3738  {
3739  unsigned char n = 8 - z_ra;
3740  z_m = z + std::complex<T> (n, 0.0);
3741 
3742  // Recurrence formula. For | Re(z) | < 8, use recursively
3743  //
3744  // DiGamma(z) = DiGamma(z+1) - 1/z
3745  std::complex<T> z_p = z + P (n - 1);
3746  for (unsigned char k = n; k > 0; k--, z_p -= 1.0)
3747  dgam -= P (1.0) / z_p;
3748  }
3749 
3750  // for | Re(z) | > 8, use derivative of C.Lanczos expansion for
3751  // LogGamma
3752  //
3753  // psi(z) = log(z) - 1/(2z) - 1/12z^2 + 3/360z^4 - 5/1260z^6
3754  // + 7/1680z^8 - 9/1188z^10 + ...
3755  //
3756  // (Abramowitz&Stegun, page 259, formula 6.3.18
3757  dgam += lanczos_approximation_psi (z_m);
3758  }
3759  return dgam;
3760  }
3761 
3762  // explicit instantiations
3763  Complex psi (const Complex& z) { return xpsi (z); }
3764  FloatComplex psi (const FloatComplex& z) { return xpsi (z); }
3765 
3766  template <typename T>
3767  static inline void
3768  fortran_psifn (const T z, const octave_idx_type n, T* ans,
3770 
3771  template <>
3772  inline void
3773  fortran_psifn<double> (const double z, const octave_idx_type n,
3774  double* ans, octave_idx_type* ierr)
3775  {
3776  octave_idx_type flag = 0;
3777  F77_XFCN (dpsifn, DPSIFN, (&z, n, 1, 1, ans, &flag, ierr));
3778  }
3779 
3780  template <>
3781  inline void
3782  fortran_psifn<float> (const float z, const octave_idx_type n,
3783  float* ans, octave_idx_type* ierr)
3784  {
3785  octave_idx_type flag = 0;
3786  F77_XFCN (psifn, PSIFN, (&z, n, 1, 1, ans, &flag, ierr));
3787  }
3788 
3789  template <typename T>
3790  T
3791  xpsi (const octave_idx_type n, T z)
3792  {
3793  T ans;
3794  octave_idx_type ierr = 0;
3795  fortran_psifn<T> (z, n, &ans, &ierr);
3796  if (ierr == 0)
3797  {
3798  // Remember that psifn and dpsifn return scales values
3799  // When n is 1: do nothing since ((-1)**(n+1)/gamma(n+1)) == 1
3800  // When n is 0: change sign since ((-1)**(n+1)/gamma(n+1)) == -1
3801  if (n > 1)
3802  // FIXME: xgamma here is a killer for our precision since it grows
3803  // way too fast.
3804  ans = ans / (pow (-1.0, n + 1) / gamma (double (n+1)));
3805  else if (n == 0)
3806  ans = -ans;
3807  }
3808  else if (ierr == 2)
3809  ans = - octave::numeric_limits<T>::Inf ();
3810  else // we probably never get here
3812 
3813  return ans;
3814  }
3815 
3816  double psi (octave_idx_type n, double z) { return xpsi (n, z); }
3817  float psi (octave_idx_type n, float z) { return xpsi (n, z); }
3818  }
3819 }
3820 
3821 #if defined (OCTAVE_USE_DEPRECATED_FUNCTIONS)
3822 
3823 ComplexMatrix besselj (double alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3824 ComplexMatrix bessely (double alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3825 ComplexMatrix besseli (double alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3826 ComplexMatrix besselk (double alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3827 ComplexMatrix besselh1 (double alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3828 ComplexMatrix besselh2 (double alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3829 
3830 ComplexMatrix besselj (const Matrix& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3831 ComplexMatrix bessely (const Matrix& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3832 ComplexMatrix besseli (const Matrix& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3833 ComplexMatrix besselk (const Matrix& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3834 ComplexMatrix besselh1 (const Matrix& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3835 ComplexMatrix besselh2 (const Matrix& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3836 
3837 ComplexMatrix besselj (const Matrix& alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3838 ComplexMatrix bessely (const Matrix& alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3839 ComplexMatrix besseli (const Matrix& alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3840 ComplexMatrix besselk (const Matrix& alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3841 ComplexMatrix besselh1 (const Matrix& alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3842 ComplexMatrix besselh2 (const Matrix& alpha, const ComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3843 
3844 ComplexNDArray besselj (double alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3845 ComplexNDArray bessely (double alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3846 ComplexNDArray besseli (double alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3847 ComplexNDArray besselk (double alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3848 ComplexNDArray besselh1 (double alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3849 ComplexNDArray besselh2 (double alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3850 
3851 ComplexNDArray besselj (const NDArray& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3852 ComplexNDArray bessely (const NDArray& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3853 ComplexNDArray besseli (const NDArray& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3854 ComplexNDArray besselk (const NDArray& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3855 ComplexNDArray besselh1 (const NDArray& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3856 ComplexNDArray besselh2 (const NDArray& alpha, const Complex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3857 
3858 ComplexNDArray besselj (const NDArray& alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3859 ComplexNDArray bessely (const NDArray& alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3860 ComplexNDArray besseli (const NDArray& alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3861 ComplexNDArray besselk (const NDArray& alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3862 ComplexNDArray besselh1 (const NDArray& alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3863 ComplexNDArray besselh2 (const NDArray& alpha, const ComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3864 
3865 ComplexMatrix besselj (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3866 ComplexMatrix bessely (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3867 ComplexMatrix besseli (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3868 ComplexMatrix besselk (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3869 ComplexMatrix besselh1 (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3870 ComplexMatrix besselh2 (const RowVector& alpha, const ComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3871 
3872 FloatComplexMatrix besselj (float alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3873 FloatComplexMatrix bessely (float alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3874 FloatComplexMatrix besseli (float alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3875 FloatComplexMatrix besselk (float alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3876 FloatComplexMatrix besselh1 (float alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3877 FloatComplexMatrix besselh2 (float alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3878 
3879 FloatComplexMatrix besselj (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3880 FloatComplexMatrix bessely (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3881 FloatComplexMatrix besseli (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3882 FloatComplexMatrix besselk (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3883 FloatComplexMatrix besselh1 (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3884 FloatComplexMatrix besselh2 (const FloatMatrix& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3885 
3886 FloatComplexMatrix besselj (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3887 FloatComplexMatrix bessely (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3888 FloatComplexMatrix besseli (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3889 FloatComplexMatrix besselk (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3890 FloatComplexMatrix besselh1 (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3891 FloatComplexMatrix besselh2 (const FloatMatrix& alpha, const FloatComplexMatrix& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3892 
3893 FloatComplexNDArray besselj (float alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3894 FloatComplexNDArray bessely (float alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3895 FloatComplexNDArray besseli (float alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3896 FloatComplexNDArray besselk (float alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3897 FloatComplexNDArray besselh1 (float alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3898 FloatComplexNDArray besselh2 (float alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3899 
3900 FloatComplexNDArray besselj (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3901 FloatComplexNDArray bessely (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3902 FloatComplexNDArray besseli (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3903 FloatComplexNDArray besselk (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3904 FloatComplexNDArray besselh1 (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3905 FloatComplexNDArray besselh2 (const FloatNDArray& alpha, const FloatComplex& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3906 
3907 FloatComplexNDArray besselj (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3908 FloatComplexNDArray bessely (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3909 FloatComplexNDArray besseli (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3910 FloatComplexNDArray besselk (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3911 FloatComplexNDArray besselh1 (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3912 FloatComplexNDArray besselh2 (const FloatNDArray& alpha, const FloatComplexNDArray& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3913 
3914 FloatComplexMatrix besselj (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselj (alpha, x, scaled, ierr); }
3915 FloatComplexMatrix bessely (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::bessely (alpha, x, scaled, ierr); }
3916 FloatComplexMatrix besseli (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besseli (alpha, x, scaled, ierr); }
3917 FloatComplexMatrix besselk (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselk (alpha, x, scaled, ierr); }
3918 FloatComplexMatrix besselh1 (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh1 (alpha, x, scaled, ierr); }
3919 FloatComplexMatrix besselh2 (const FloatRowVector& alpha, const FloatComplexColumnVector& x, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::besselh2 (alpha, x, scaled, ierr); }
3920 
3921 ComplexMatrix airy (const ComplexMatrix& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::airy (z, deriv, scaled, ierr); }
3922 ComplexMatrix biry (const ComplexMatrix& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::biry (z, deriv, scaled, ierr); }
3923 
3924 ComplexNDArray airy (const ComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::airy (z, deriv, scaled, ierr); }
3925 ComplexNDArray biry (const ComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::biry (z, deriv, scaled, ierr); }
3926 
3927 FloatComplexMatrix airy (const FloatComplexMatrix& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::airy (z, deriv, scaled, ierr); }
3928 FloatComplexMatrix biry (const FloatComplexMatrix& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::biry (z, deriv, scaled, ierr); }
3929 
3930 FloatComplexNDArray airy (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::airy (z, deriv, scaled, ierr); }
3931 FloatComplexNDArray biry (const FloatComplexNDArray& z, bool deriv, bool scaled, Array<octave_idx_type>& ierr) { return octave::math::biry (z, deriv, scaled, ierr); }
3932 
3933 Array<double> betainc (double x, double a, const Array<double>& b) { return octave::math::betainc (x, a, b); }
3934 Array<double> betainc (double x, const Array<double>& a, double b) { return octave::math::betainc (x, a, b); }
3935 Array<double> betainc (double x, const Array<double>& a, const Array<double>& b) { return octave::math::betainc (x, a, b); }
3936 
3937 Array<double> betainc (const Array<double>& x, double a, double b) { return octave::math::betainc (x, a, b); }
3938 Array<double> betainc (const Array<double>& x, double a, const Array<double>& b) { return octave::math::betainc (x, a, b); }
3939 Array<double> betainc (const Array<double>& x, const Array<double>& a, double b) { return octave::math::betainc (x, a, b); }
3940 Array<double> betainc (const Array<double>& x, const Array<double>& a, const Array<double>& b) { return octave::math::betainc (x, a, b); }
3941 
3942 Array<float> betainc (float x, float a, const Array<float>& b) { return octave::math::betainc (x, a, b); }
3943 Array<float> betainc (float x, const Array<float>& a, float b) { return octave::math::betainc (x, a, b); }
3944 Array<float> betainc (float x, const Array<float>& a, const Array<float>& b) { return octave::math::betainc (x, a, b); }
3945 
3946 Array<float> betainc (const Array<float>& x, float a, float b) { return octave::math::betainc (x, a, b); }
3947 Array<float> betainc (const Array<float>& x, float a, const Array<float>& b) { return octave::math::betainc (x, a, b); }
3948 Array<float> betainc (const Array<float>& x, const Array<float>& a, float b) { return octave::math::betainc (x, a, b); }
3949 Array<float> betainc (const Array<float>& x, const Array<float>& a, const Array<float>& b) { return octave::math::betainc (x, a, b); }
3950 
3951 Matrix gammainc (double x, const Matrix& a) { return octave::math::gammainc (x, a); }
3952 Matrix gammainc (const Matrix& x, double a) { return octave::math::gammainc (x, a); }
3953 Matrix gammainc (const Matrix& x, const Matrix& a) { return octave::math::gammainc (x, a); }
3954 
3955 NDArray gammainc (double x, const NDArray& a) { return octave::math::gammainc (x, a); }
3956 NDArray gammainc (const NDArray& x, double a) { return octave::math::gammainc (x, a); }
3957 NDArray gammainc (const NDArray& x, const NDArray& a) { return octave::math::gammainc (x, a); }
3958 
3959 FloatMatrix gammainc (float x, const FloatMatrix& a) { return octave::math::gammainc (x, a); }
3960 FloatMatrix gammainc (const FloatMatrix& x, float a) { return octave::math::gammainc (x, a); }
3961 FloatMatrix gammainc (const FloatMatrix& x, const FloatMatrix& a) { return octave::math::gammainc (x, a); }
3962 
3963 FloatNDArray gammainc (float x, const FloatNDArray& a) { return octave::math::gammainc (x, a); }
3964 FloatNDArray gammainc (const FloatNDArray& x, float a) { return octave::math::gammainc (x, a); }
3965 FloatNDArray gammainc (const FloatNDArray& x, const FloatNDArray& a) { return octave::math::gammainc (x, a); }
3966 
3967 Array<double> betaincinv (double x, double a, const Array<double>& b) { return octave::math::betaincinv (x, a, b); }
3968 Array<double> betaincinv (double x, const Array<double>& a, double b) { return octave::math::betaincinv (x, a, b); }
3969 Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b) { return octave::math::betaincinv (x, a, b); }
3970 
3971 Array<double> betaincinv (const Array<double>& x, double a, double b) { return octave::math::betaincinv (x, a, b); }
3972 Array<double> betaincinv (const Array<double>& x, double a, const Array<double>& b) { return octave::math::betaincinv (x, a, b); }
3973 Array<double> betaincinv (const Array<double>& x, const Array<double>& a, double b) { return octave::math::betaincinv (x, a, b); }
3974 Array<double> betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b) { return octave::math::betaincinv (x, a, b); }
3975 
3976 #endif
uint32_t id
Definition: graphics.cc:11587
float cbrt(float x)
Definition: lo-specfun.cc:703
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
Definition: fCMatrix.h:189
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE const F77_INT F77_INT * ierr
Octave interface to the compression and uncompression libraries.
Definition: aepbalance.cc:47
subroutine dpsifn(X, N, KODE, M, ANS, NZ, IERR)
Definition: dpsifn.f:2
subroutine xerf(x, result)
Definition: xerf.f:1
std::string str(char sep= 'x') const
Definition: dim-vector.cc:73
double psi(double z)
Definition: lo-specfun.cc:3714
static FloatComplex cbesh1(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:1637
subroutine dlgams(X, DLGAM, SGNGAM)
Definition: dlgams.f:2
subroutine xdgamma(x, result)
Definition: xdgamma.f:1
static OCTAVE_NORETURN void err_betaincinv_nonconformant(const dim_vector &d1, const dim_vector &d2, const dim_vector &d3)
Definition: lo-specfun.cc:2252
float erfcx(float x)
Definition: lo-specfun.cc:271
subroutine cairy(Z, ID, KODE, AI, NZ, IERR)
Definition: cairy.f:1
ar P
Definition: __luinc__.cc:49
subroutine zairy(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
Definition: zairy.f:1
function erf(X)
Definition: erf.f:2
std::complex< double > erfi(std::complex< double > z, double relerr=0)
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
identity matrix If supplied two scalar respectively For allows like xample val
Definition: data.cc:5068
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE * f
FloatComplexMatrix besselh1(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
Definition: lo-specfun.cc:1957
static Complex zbesk(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:963
FloatComplex log1p(const FloatComplex &x)
Definition: lo-specfun.cc:684
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
Complex besselk(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1337
double betaincinv(double y, double p, double q)
Definition: lo-specfun.cc:3193
bool isnan(double x)
Definition: lo-mappers.cc:347
for large enough k
Definition: lu.cc:606
double atanh(double x)
Definition: lo-specfun.cc:149
Complex rc_lgamma(double x)
Definition: lo-specfun.cc:373
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the base of natural logarithms The constant ex $e satisfies the equation log(e)
function atanh(X)
Definition: atanh.f:2
Complex besselh1(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1338
subroutine cbiry(Z, ID, KODE, BI, IERR)
Definition: cbiry.f:1
static T xlog(const T &x)
Definition: lo-specfun.cc:3614
std::complex< double > Dawson(std::complex< double > z, double relerr=0)
octave_value real(void) const
Definition: ov.h:1382
double lo_ieee_nan_value(void)
Definition: lo-ieee.cc:126
u
Definition: lu.cc:138
FloatComplexMatrix besselh2(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
Definition: lo-specfun.cc:1958
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function t
Definition: ov-usr-fcn.cc:935
double expm1(double x)
Definition: lo-specfun.cc:473
T xcbrt(T x)
Definition: lo-specfun.cc:633
double asinh(double x)
Definition: lo-specfun.cc:105
subroutine xderfc(x, result)
Definition: xderfc.f:1
s
Definition: file-io.cc:2682
Complex asin(const Complex &x)
Definition: lo-mappers.cc:150
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
octave_idx_type a_nc
Definition: sylvester.cc:72
subroutine xerfc(x, result)
Definition: xerfc.f:1
i e
Definition: data.cc:2724
subroutine xdacosh(x, result)
Definition: xdacosh.f:1
double cbrt(double x)
Definition: lo-specfun.cc:648
std::complex< double > erf(std::complex< double > z, double relerr=0)
FloatComplexMatrix besseli(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
Definition: lo-specfun.cc:1955
octave_idx_type rows(void) const
Definition: Array.h:401
double acosh(double x)
Definition: lo-specfun.cc:61
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition: dim-vector.h:389
octave_value imag(void) const
Definition: ov.h:1373
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
static Complex zbesj(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:769
#define ALL_BESSEL(name, fcn)
Definition: lo-specfun.cc:1943
FloatComplexNDArray biry(const FloatComplexNDArray &z, bool deriv, bool scaled, Array< octave_idx_type > &ierr)
Definition: lo-specfun.cc:2223
subroutine xgamma(x, result)
Definition: xgamma.f:1
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:439
float erfi(float x)
Definition: lo-specfun.cc:289
subroutine xdatanh(x, result)
Definition: xdatanh.f:1
FloatNDArray gammainc(const FloatNDArray &x, const FloatNDArray &a)
Definition: lo-specfun.cc:2842
double h
Definition: graphics.cc:11205
subroutine algams(X, ALGAM, SGNGAM)
Definition: algams.f:2
double gamma(double x)
Definition: lo-specfun.cc:325
static void fortran_psifn(const T z, const octave_idx_type n, T *ans, octave_idx_type *ierr)
Complex besseli(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1336
FloatComplex expm1(const FloatComplex &x)
Definition: lo-specfun.cc:568
octave_idx_type a_nr
Definition: sylvester.cc:71
FloatComplexMatrix besselk(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
Definition: lo-specfun.cc:1956
function asinh(X)
Definition: asinh.f:2
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
static double do_erfinv(double x, bool refine)
Definition: lo-specfun.cc:2897
static FloatComplex cbesk(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:1586
std::complex< double > w(std::complex< double > z, double relerr=0)
static Complex zbesi(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:898
static double do_erfcinv(double x, bool refine)
Definition: lo-specfun.cc:2975
Array< float > betainc(const Array< float > &x, const Array< float > &a, const Array< float > &b)
Definition: lo-specfun.cc:2524
double signum(double x)
Definition: lo-mappers.h:259
subroutine zbiry(ZR, ZI, ID, KODE, BIR, BII, IERR)
Definition: zbiry.f:1
float lgamma(float x)
Definition: lo-specfun.cc:427
float dawson(float x)
Definition: lo-specfun.cc:307
static Complex zbesh2(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:1064
subroutine xacosh(x, result)
Definition: xacosh.f:1
subroutine psifn(X, N, KODE, M, ANS, NZ, IERR)
Definition: psifn.f:2
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1028
static Complex bessel_return_value(const Complex &val, octave_idx_type ierr)
Definition: lo-specfun.cc:731
subroutine xgammainc(a, x, result)
Definition: xgmainc.f:1
static T lanczos_approximation_psi(const T zc)
Definition: lo-specfun.cc:3635
octave_int< T > pow(const octave_int< T > &a, const octave_int< T > &b)
double tmp
Definition: data.cc:6300
octave_value retval
Definition: data.cc:6294
F77_RET_T F77_FUNC(xstopx, XSTOPX) const
Definition: f77-fcn.c:53
bool finite(double x)
Definition: lo-mappers.cc:367
std::complex< double > erfcx(std::complex< double > z, double relerr=0)
void fortran_psifn< float >(const float z, const octave_idx_type n, float *ans, octave_idx_type *ierr)
Definition: lo-specfun.cc:3782
subroutine zbesh(ZR, ZI, FNU, KODE, M, N, CYR, CYI, NZ, IERR)
Definition: zbesh.f:1
Definition: dMatrix.h:37
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
subroutine xdasinh(x, result)
Definition: xdasinh.f:1
Array< double > betaincinv(const Array< double > &x, const Array< double > &a, const Array< double > &b)
Definition: lo-specfun.cc:3489
With real return the complex result
Definition: data.cc:3375
double erfcinv(double x)
Definition: lo-specfun.cc:3049
bool isinf(double x)
Definition: lo-mappers.cc:387
subroutine cbesh(Z, FNU, KODE, M, N, CY, NZ, IERR)
Definition: cbesh.f:1
function acosh(X)
Definition: acosh.f:2
subroutine xdbetai(x, a, b, result)
Definition: xdbetai.f:1
static FloatComplex cbesy(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:1464
double erf(double x)
Definition: lo-specfun.cc:193
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:228
#define Inf
Definition: Faddeeva.cc:247
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
SparseMatrix atan2(const double &x, const SparseMatrix &y)
Definition: dSparse.cc:606
static Complex zbesy(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:827
FloatComplex erfc(const FloatComplex &x)
Definition: lo-specfun.cc:263
static FloatComplex cbesh2(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:1681
double betainc(double x, double a, double b)
Definition: lo-specfun.cc:2265
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition: CMatrix.h:184
Complex(* dptr)(const Complex &, double, int, octave_idx_type &)
Definition: lo-specfun.cc:1109
static FloatComplex cbesi(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:1531
static double wi[256]
Definition: randmtzig.cc:440
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
double log1p(double x)
Definition: lo-specfun.cc:587
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
bool negative_sign(double x)
Definition: lo-mappers.cc:324
p
Definition: lu.cc:138
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol zero divided by nd tex zero divided by nd ifnottex and any operation involving another NaN value(5+NaN).Note that NaN always compares not equal to NaN(NaN!
subroutine xbetai(x, a, b, result)
Definition: xbetai.f:1
issues an error eealso double
Definition: ov-bool-mat.cc:594
static FloatComplex cbesj(const FloatComplex &z, float alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:1408
static Complex zbesh1(const Complex &z, double alpha, int kode, octave_idx_type &ierr)
Definition: lo-specfun.cc:1018
Complex rc_log1p(double x)
Definition: lo-specfun.cc:2872
the element is set to zero In other the statement xample y
Definition: data.cc:5342
b
Definition: cellfun.cc:398
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol NaN(Not a Number).NaN is the result of operations which do not produce a well defined 0 result.Common operations which produce a NaN are arithmetic with infinity ex($\infty-\infty $)
subroutine xsgammainc(a, x, result)
Definition: xsgmainc.f:1
Complex besselh2(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1339
double erfinv(double x)
Definition: lo-specfun.cc:2960
void ellipj(double u, double m, double &sn, double &cn, double &dn, double &err)
Definition: lo-specfun.cc:3511
#define F77_CONST_CMPLX_ARG(x)
Definition: f77-fcn.h:342
void fortran_psifn< double >(const double z, const octave_idx_type n, double *ans, octave_idx_type *ierr)
Definition: lo-specfun.cc:3773
double floor(double x)
Definition: lo-mappers.cc:330
double erfc(double x)
Definition: lo-specfun.cc:232
double gammainc(double x, double a, bool &err)
Definition: lo-specfun.cc:2547
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
OCTAVE_EXPORT octave_value_list error nd deftypefn *const octave_scalar_map err
Definition: error.cc:1036
FloatComplex(* fptr)(const FloatComplex &, float, int, octave_idx_type &)
Definition: lo-specfun.cc:1724
Complex airy(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1971
FloatComplexMatrix besselj(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
Definition: lo-specfun.cc:1953
double lgamma(double x)
Definition: lo-specfun.cc:353
FloatComplexMatrix bessely(const FloatRowVector &alpha, const FloatComplexColumnVector &x, bool scaled, Array< octave_idx_type > &ierr)
Definition: lo-specfun.cc:1954
std::complex< double > Complex
Definition: oct-cmplx.h:31
const T * fortran_vec(void) const
Definition: Array.h:584
subroutine xatanh(x, result)
Definition: xatanh.f:1
FloatComplexNDArray airy(const FloatComplexNDArray &z, bool deriv, bool scaled, Array< octave_idx_type > &ierr)
Definition: lo-specfun.cc:2207
ar
Definition: __qp__.cc:492
static bool is_integer_value(double x)
Definition: lo-specfun.cc:763
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:136
octave_idx_type cols(void) const
Definition: Array.h:409
double betain(double x, double p, double q, double beta, bool &err)
Definition: lo-specfun.cc:3076
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:854
dim_vector dv
Definition: sub2ind.cc:263
T x_nint(T x)
Definition: lo-mappers.h:299
static OCTAVE_NORETURN void err_betainc_nonconformant(const dim_vector &d1, const dim_vector &d2, const dim_vector &d3)
Definition: lo-specfun.cc:2239
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
Complex besselj(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1334
static const double pi
Definition: lo-specfun.cc:3610
Complex bessely(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1335
static Complex do_bessel(dptr f, const char *, double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1112
subroutine xasinh(x, result)
Definition: xasinh.f:1
subroutine xderf(x, result)
Definition: xderf.f:1
Complex biry(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:2005
std::complex< double > erfc(std::complex< double > z, double relerr=0)