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