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