Bug Summary

File:liboctave/numeric/lo-specfun.cc
Location:line 2428, column 3
Description:Undefined or garbage value returned to caller

Annotated Source Code

1/*
2
3Copyright (C) 1996-2013 John W. Eaton
4Copyright (C) 2010 Jaroslav Hajek
5Copyright (C) 2010 VZLU Prague
6
7This file is part of Octave.
8
9Octave is free software; you can redistribute it and/or modify it
10under the terms of the GNU General Public License as published by the
11Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14Octave is distributed in the hope that it will be useful, but WITHOUT
15ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17for more details.
18
19You should have received a copy of the GNU General Public License
20along with Octave; see the file COPYING. If not, see
21<http://www.gnu.org/licenses/>.
22
23*/
24
25#ifdef HAVE_CONFIG_H1
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
51extern "C"
52{
53 F77_RET_Tint
54 F77_FUNC (zbesj, 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_Tint
60 F77_FUNC (zbesy, 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_Tint
66 F77_FUNC (zbesi, 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_Tint
72 F77_FUNC (zbesk, 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_Tint
78 F77_FUNC (zbesh, 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_Tint
84 F77_FUNC (cbesj, 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_Tint
89 F77_FUNC (cbesy, 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_Tint
95 F77_FUNC (cbesi, 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_Tint
100 F77_FUNC (cbesk, 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_Tint
105 F77_FUNC (cbesh, 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_Tint
111 F77_FUNC (zairy, 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_Tint
117 F77_FUNC (cairy, 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_Tint
122 F77_FUNC (zbiry, 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_Tint
127 F77_FUNC (cbiry, 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_Tint
132 F77_FUNC (xdacosh, XDACOSH)xdacosh_ (const double&, double&);
133
134 F77_RET_Tint
135 F77_FUNC (xacosh, XACOSH)xacosh_ (const float&, float&);
136
137 F77_RET_Tint
138 F77_FUNC (xdasinh, XDASINH)xdasinh_ (const double&, double&);
139
140 F77_RET_Tint
141 F77_FUNC (xasinh, XASINH)xasinh_ (const float&, float&);
142
143 F77_RET_Tint
144 F77_FUNC (xdatanh, XDATANH)xdatanh_ (const double&, double&);
145
146 F77_RET_Tint
147 F77_FUNC (xatanh, XATANH)xatanh_ (const float&, float&);
148
149 F77_RET_Tint
150 F77_FUNC (xderf, XDERF)xderf_ (const double&, double&);
151
152 F77_RET_Tint
153 F77_FUNC (xerf, XERF)xerf_ (const float&, float&);
154
155 F77_RET_Tint
156 F77_FUNC (xderfc, XDERFC)xderfc_ (const double&, double&);
157
158 F77_RET_Tint
159 F77_FUNC (xerfc, XERFC)xerfc_ (const float&, float&);
160
161 F77_RET_Tint
162 F77_FUNC (xdbetai, XDBETAI)xdbetai_ (const double&, const double&,
163 const double&, double&);
164
165 F77_RET_Tint
166 F77_FUNC (xbetai, XBETAI)xbetai_ (const float&, const float&,
167 const float&, float&);
168
169 F77_RET_Tint
170 F77_FUNC (xdgamma, XDGAMMA)xdgamma_ (const double&, double&);
171
172 F77_RET_Tint
173 F77_FUNC (xgamma, XGAMMA)xgamma_ (const float&, float&);
174
175 F77_RET_Tint
176 F77_FUNC (xgammainc, XGAMMAINC)xgammainc_ (const double&, const double&, double&);
177
178 F77_RET_Tint
179 F77_FUNC (xsgammainc, XSGAMMAINC)xsgammainc_ (const float&, const float&, float&);
180
181 F77_RET_Tint
182 F77_FUNC (dlgams, DLGAMS)dlgams_ (const double&, double&, double&);
183
184 F77_RET_Tint
185 F77_FUNC (algams, ALGAMS)algams_ (const float&, float&, float&);
186}
187
188#if !defined (HAVE_ACOSH1)
189double
190acosh (double x)
191{
192 double retval;
193 F77_XFCN (xdacosh, XDACOSH, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xdacosh_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xdacosh_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
194 return retval;
195}
196#endif
197
198#if !defined (HAVE_ACOSHF1)
199float
200acoshf (float x)
201{
202 float retval;
203 F77_XFCN (xacosh, XACOSH, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xacosh_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xacosh_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
204 return retval;
205}
206#endif
207
208#if !defined (HAVE_ASINH1)
209double
210asinh (double x)
211{
212 double retval;
213 F77_XFCN (xdasinh, XDASINH, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xdasinh_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xdasinh_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
214 return retval;
215}
216#endif
217
218#if !defined (HAVE_ASINHF1)
219float
220asinhf (float x)
221{
222 float retval;
223 F77_XFCN (xasinh, XASINH, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xasinh_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xasinh_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
224 return retval;
225}
226#endif
227
228#if !defined (HAVE_ATANH1)
229double
230atanh (double x)
231{
232 double retval;
233 F77_XFCN (xdatanh, XDATANH, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xdatanh_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xdatanh_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
234 return retval;
235}
236#endif
237
238#if !defined (HAVE_ATANHF1)
239float
240atanhf (float x)
241{
242 float retval;
243 F77_XFCN (xatanh, XATANH, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xatanh_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xatanh_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
244 return retval;
245}
246#endif
247
248#if !defined (HAVE_ERF1)
249double
250erf (double x)
251{
252 double retval;
253 F77_XFCN (xderf, XDERF, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xderf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xderf_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
254 return retval;
255}
256#endif
257
258#if !defined (HAVE_ERFF1)
259float
260erff (float x)
261{
262 float retval;
263 F77_XFCN (xerf, XERF, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xerf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xerf_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
264 return retval;
265}
266#endif
267
268#if !defined (HAVE_ERFC1)
269double
270erfc (double x)
271{
272 double retval;
273 F77_XFCN (xderfc, XDERFC, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xderfc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xderfc_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
274 return retval;
275}
276#endif
277
278#if !defined (HAVE_ERFCF1)
279float
280erfcf (float x)
281{
282 float retval;
283 F77_XFCN (xerfc, XERFC, (x, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xerfc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xerfc_ (x, retval); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
284 return retval;
285}
286#endif
287
288// Complex error function from the Faddeeva package
289Complex
290erf (const Complex& x)
291{
292 return Faddeeva::erf (x);
293}
294FloatComplex
295erf (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
303Complex
304erfc (const Complex& x)
305{
306 return Faddeeva::erfc (x);
307}
308FloatComplex
309erfc (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
317float erfcx (float x) { return Faddeeva::erfcx(x); }
318double erfcx (double x) { return Faddeeva::erfcx(x); }
319Complex
320erfcx (const Complex& x)
321{
322 return Faddeeva::erfcx (x);
323}
324FloatComplex
325erfcx (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
333float erfi (float x) { return Faddeeva::erfi(x); }
334double erfi (double x) { return Faddeeva::erfi(x); }
335Complex
336erfi (const Complex& x)
337{
338 return Faddeeva::erfi (x);
339}
340FloatComplex
341erfi (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
349float dawson (float x) { return Faddeeva::Dawson(x); }
350double dawson (double x) { return Faddeeva::Dawson(x); }
351Complex
352dawson (const Complex& x)
353{
354 return Faddeeva::Dawson (x);
355}
356FloatComplex
357dawson (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
364double
365xgamma (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_TGAMMA1)
378 result = tgamma (x);
379#else
380 F77_XFCN (xdgamma, XDGAMMA, (x, result))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xdgamma_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xdgamma_ (x, result); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
381#endif
382 if (xisinf (result) && (static_cast<int> (gnulib::floor (x)) % 2))
383 result = -octave_Inf;
384 }
385
386 return result;
387}
388
389double
390xlgamma (double x)
391{
392#if defined (HAVE_LGAMMA1)
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))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dlgams_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dlgams_ (x, result, sgngam); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
404
405 return result;
406#endif
407}
408
409Complex
410rc_lgamma (double x)
411{
412 double result;
413
414#if defined (HAVE_LGAMMA_R1)
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))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "dlgams_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dlgams_ (x, result, sgngam); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
426
427#endif
428
429 if (sgngam < 0)
430 return result + Complex (0., M_PI3.14159265358979323846);
431 else
432 return result;
433}
434
435float
436xgamma (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_TGAMMA1)
449 result = tgammaf (x);
450#else
451 F77_XFCN (xgamma, XGAMMA, (x, result))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xgamma_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xgamma_ (x, result); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
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
460float
461xlgamma (float x)
462{
463#if defined (HAVE_LGAMMAF1)
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))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "algams_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; algams_ (x, result, sgngam); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
475
476 return result;
477#endif
478}
479
480FloatComplex
481rc_lgamma (float x)
482{
483 float result;
484
485#if defined (HAVE_LGAMMAF_R1)
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))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "algams_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; algams_ (x, result, sgngam); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
497
498#endif
499
500 if (sgngam < 0)
501 return result + FloatComplex (0., M_PI3.14159265358979323846);
502 else
503 return result;
504}
505
506#if !defined (HAVE_EXPM11)
507double
508expm1 (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
542Complex
543expm1 (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_EXPM1F1)
562float
563expm1f (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
597FloatComplex
598expm1 (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_LOG1P1)
617double
618log1p (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
640Complex
641log1p (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_CBRT1)
660double 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_LOG1PF1)
676float
677log1pf (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
699FloatComplex
700log1p (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_CBRTF1)
719float 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
734static inline Complex
735zbesj (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
736
737static inline Complex
738zbesy (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
739
740static inline Complex
741zbesi (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
742
743static inline Complex
744zbesk (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
745
746static inline Complex
747zbesh1 (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
748
749static inline Complex
750zbesh2 (const Complex& z, double alpha, int kode, octave_idx_type& ierr);
751
752static inline Complex
753bessel_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
779static inline bool
780is_integer_value (double x)
781{
782 return x == static_cast<double> (static_cast<long> (x));
783}
784
785static inline Complex
786zbesj (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)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_PI3.14159265358979323846 * alpha) * zbesj (z, alpha, kode, ierr);
828
829 if (ierr == 0 || ierr == 3)
830 {
831 tmp -= sin (M_PI3.14159265358979323846 * 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
842static inline Complex
843zbesy (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)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_PI3.14159265358979323846 * alpha) * zbesy (z, alpha, kode, ierr);
898
899 if (ierr == 0 || ierr == 3)
900 {
901 tmp += sin (M_PI3.14159265358979323846 * 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
912static inline Complex
913zbesi (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)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_PI3.14159265358979323846) * sin (M_PI3.14159265358979323846 * 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
976static inline Complex
977zbesk (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)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
1031static inline Complex
1032zbesh1 (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)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_PI3.14159265358979323846 * alpha * eye) * zbesh1 (z, alpha, kode, ierr);
1070
1071 retval = bessel_return_value (tmp, ierr);
1072 }
1073
1074 return retval;
1075}
1076
1077static inline Complex
1078zbesh2 (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)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_PI3.14159265358979323846 * alpha * eye) * zbesh2 (z, alpha, kode, ierr);
1116
1117 retval = bessel_return_value (tmp, ierr);
1118 }
1119
1120 return retval;
1121}
1122
1123typedef Complex (*dptr) (const Complex&, double, int, octave_idx_type&);
1124
1125static inline Complex
1126do_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
1136static inline ComplexMatrix
1137do_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
1154static inline ComplexMatrix
1155do_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
1172static inline ComplexMatrix
1173do_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
1204static inline ComplexNDArray
1205do_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
1220static inline ComplexNDArray
1221do_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
1236static inline ComplexNDArray
1237do_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
1260static inline ComplexMatrix
1261do_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
1352ALL_BESSEL (besselj, zbesj)
1353ALL_BESSEL (bessely, zbesy)
1354ALL_BESSEL (besseli, zbesi)
1355ALL_BESSEL (besselk, zbesk)
1356ALL_BESSEL (besselh1, zbesh1)
1357ALL_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
1369static inline FloatComplex
1370cbesj (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1371
1372static inline FloatComplex
1373cbesy (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1374
1375static inline FloatComplex
1376cbesi (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1377
1378static inline FloatComplex
1379cbesk (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1380
1381static inline FloatComplex
1382cbesh1 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1383
1384static inline FloatComplex
1385cbesh2 (const FloatComplex& z, float alpha, int kode, octave_idx_type& ierr);
1386
1387static inline FloatComplex
1388bessel_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
1416static inline bool
1417is_integer_value (float x)
1418{
1419 return x == static_cast<float> (static_cast<long> (x));
1420}
1421
1422static inline FloatComplex
1423cbesj (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)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_PI3.14159265358979323846) * alpha)
1460 * cbesj (z, alpha, kode, ierr);
1461
1462 if (ierr == 0 || ierr == 3)
1463 {
1464 tmp -= sinf (static_cast<float> (M_PI3.14159265358979323846) * 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
1476static inline FloatComplex
1477cbesy (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)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_PI3.14159265358979323846) * alpha)
1525 * cbesy (z, alpha, kode, ierr);
1526
1527 if (ierr == 0 || ierr == 3)
1528 {
1529 tmp += sinf (static_cast<float> (M_PI3.14159265358979323846) * 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
1541static inline FloatComplex
1542cbesi (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)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_PI3.14159265358979323846)
1574 * sinf (static_cast<float> (M_PI3.14159265358979323846) * 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
1594static inline FloatComplex
1595cbesk (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)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
1644static inline FloatComplex
1645cbesh1 (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)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_PI3.14159265358979323846) * alpha * eye)
1679 * cbesh1 (z, alpha, kode, ierr);
1680
1681 retval = bessel_return_value (tmp, ierr);
1682 }
1683
1684 return retval;
1685}
1686
1687static inline FloatComplex
1688cbesh2 (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)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_PI3.14159265358979323846) * alpha * eye)
1722 * cbesh2 (z, alpha, kode, ierr);
1723
1724 retval = bessel_return_value (tmp, ierr);
1725 }
1726
1727 return retval;
1728}
1729
1730typedef FloatComplex (*fptr) (const FloatComplex&, float, int,
1731 octave_idx_type&);
1732
1733static inline FloatComplex
1734do_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
1744static inline FloatComplexMatrix
1745do_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
1762static inline FloatComplexMatrix
1763do_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
1781static inline FloatComplexMatrix
1782do_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
1814static inline FloatComplexNDArray
1815do_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
1830static inline FloatComplexNDArray
1831do_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
1846static inline FloatComplexNDArray
1847do_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
1871static inline FloatComplexMatrix
1872do_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
1963ALL_BESSEL (besselj, cbesj)
1964ALL_BESSEL (bessely, cbesy)
1965ALL_BESSEL (besseli, cbesi)
1966ALL_BESSEL (besselk, cbesk)
1967ALL_BESSEL (besselh1, cbesh1)
1968ALL_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
1980Complex
1981airy (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)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
2014Complex
2015biry (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)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
2046ComplexMatrix
2047airy (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
2064ComplexMatrix
2065biry (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
2082ComplexNDArray
2083airy (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
2098ComplexNDArray
2099biry (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
2114FloatComplex
2115airy (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)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
2148FloatComplex
2149biry (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)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
2181FloatComplexMatrix
2182airy (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
2199FloatComplexMatrix
2200biry (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
2217FloatComplexNDArray
2218airy (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
2233FloatComplexNDArray
2234biry (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
2249static void
2250gripe_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
2262static void
2263gripe_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
2275double
2276betainc (double x, double a, double b)
2277{
2278 double retval;
2279 F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xdbetai_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xdbetai_ (x, a, b, retval); octave_interrupt_immediately--
; octave_restore_current_context (saved_context); } } while (
0)
;
2280 return retval;
2281}
2282
2283Array<double>
2284betainc (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
2299Array<double>
2300betainc (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
2315Array<double>
2316betainc (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
2338Array<double>
2339betainc (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
2354Array<double>
2355betainc (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
2377Array<double>
2378betainc (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
2400Array<double>
2401betainc (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
2423float
2424betainc (float x, float a, float b)
2425{
2426 float retval;
4
'retval' declared without an initial value
2427 F77_XFCN (xbetai, XBETAI, (x, a, b, retval))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xbetai_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xbetai_ (x, a, b, retval); octave_interrupt_immediately--
; octave_restore_current_context (saved_context); } } while (
0)
;
2428 return retval;
5
Undefined or garbage value returned to caller
2429}
2430
2431Array<float>
2432betainc (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
2447Array<float>
2448betainc (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
2463Array<float>
2464betainc (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
2486Array<float>
2487betainc (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++)
1
Assuming 'i' is < 'nel'
2
Loop condition is true. Entering loop body
2497 *pretval++ = betainc (x(i), a, b);
3
Calling 'betainc'
2498
2499 return retval;
2500}
2501
2502Array<float>
2503betainc (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
2525Array<float>
2526betainc (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
2548Array<float>
2549betainc (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
2573double
2574gammainc (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))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xgammainc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xgammainc_ (a, x, retval); octave_interrupt_immediately--
; octave_restore_current_context (saved_context); } } while (
0)
;
2589
2590 return retval;
2591}
2592
2593Matrix
2594gammainc (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
2615done:
2616
2617 return retval;
2618}
2619
2620Matrix
2621gammainc (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
2642done:
2643
2644 return retval;
2645}
2646
2647Matrix
2648gammainc (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
2681done:
2682
2683 return retval;
2684}
2685
2686NDArray
2687gammainc (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
2707done:
2708
2709 return retval;
2710}
2711
2712NDArray
2713gammainc (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
2733done:
2734
2735 return retval;
2736}
2737
2738NDArray
2739gammainc (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
2773done:
2774
2775 return retval;
2776}
2777
2778float
2779gammainc (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))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "xsgammainc_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; xsgammainc_ (a, x, retval); octave_interrupt_immediately--
; octave_restore_current_context (saved_context); } } while (
0)
;
2794
2795 return retval;
2796}
2797
2798FloatMatrix
2799gammainc (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
2820done:
2821
2822 return retval;
2823}
2824
2825FloatMatrix
2826gammainc (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
2847done:
2848
2849 return retval;
2850}
2851
2852FloatMatrix
2853gammainc (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
2886done:
2887
2888 return retval;
2889}
2890
2891FloatNDArray
2892gammainc (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
2912done:
2913
2914 return retval;
2915}
2916
2917FloatNDArray
2918gammainc (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
2938done:
2939
2940 return retval;
2941}
2942
2943FloatNDArray
2944gammainc (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
2978done:
2979
2980 return retval;
2981}
2982
2983
2984Complex 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
2990FloatComplex 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
3004static 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
3067double erfinv (double x)
3068{
3069 return do_erfinv (x, true);
3070}
3071
3072float 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
3081static 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
3149double erfcinv (double x)
3150{
3151 return do_erfcinv (x, true);
3152}
3153
3154float 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//
3175double
3176betain (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//
3292double
3293betaincinv (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
3485Array<double>
3486betaincinv (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
3501Array<double>
3502betaincinv (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
3517Array<double>
3518betaincinv (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
3540Array<double>
3541betaincinv (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
3556Array<double>
3557betaincinv (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
3579Array<double>
3580betaincinv (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
3602Array<double>
3603betaincinv (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
3626void
3627ellipj (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
3697void
3698ellipj (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}