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