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