|
3146
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 1996 John W. Eaton |
|
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 2, or (at your option) any |
|
|
10 later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, write to the Free |
|
|
19 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
|
20 |
|
|
21 */ |
|
|
22 |
|
|
23 #ifdef HAVE_CONFIG_H |
|
|
24 #include <config.h> |
|
|
25 #endif |
|
|
26 |
|
|
27 #include "Range.h" |
|
|
28 #include "dColVector.h" |
|
|
29 #include "dMatrix.h" |
|
|
30 #include "f77-fcn.h" |
|
|
31 #include "lo-error.h" |
|
|
32 #include "mx-inlines.cc" |
|
|
33 |
|
|
34 extern "C" |
|
|
35 { |
|
|
36 int F77_FCN (rjbesl, RJBESL) (const double&, const double&, |
|
|
37 const int&, double*, int&); |
|
|
38 |
|
|
39 int F77_FCN (rybesl, RYBESL) (const double&, const double&, |
|
|
40 const int&, double*, int&); |
|
|
41 |
|
|
42 int F77_FCN (ribesl, RIBESL) (const double&, const double&, |
|
|
43 const int&, const int&, double*, int&); |
|
|
44 |
|
|
45 int F77_FCN (rkbesl, RKBESL) (const double&, const double&, |
|
|
46 const int&, const int&, double*, int&); |
|
|
47 |
|
|
48 int F77_FCN (xdacosh, XDACOSH) (const double&, double&); |
|
|
49 |
|
|
50 int F77_FCN (xdasinh, XDASINH) (const double&, double&); |
|
|
51 |
|
|
52 int F77_FCN (xdatanh, XDATANH) (const double&, double&); |
|
|
53 |
|
|
54 int F77_FCN (xderf, XDERF) (const double&, double&); |
|
|
55 |
|
|
56 int F77_FCN (xderfc, XDERFC) (const double&, double&); |
|
|
57 |
|
|
58 int F77_FCN (xdbetai, XDBETAI) (const double&, const double&, |
|
|
59 const double&, double&); |
|
|
60 |
|
|
61 int F77_FCN (xdgamma, XDGAMMA) (const double&, double&); |
|
|
62 |
|
|
63 int F77_FCN (xdgami, XDGAMI) (const double&, const double&, double&); |
|
|
64 |
|
|
65 int F77_FCN (dlgams, DLGAMS) (const double&, double&, double&); |
|
|
66 } |
|
|
67 |
|
|
68 #if !defined (HAVE_ACOSH) |
|
|
69 double |
|
|
70 acosh (double x) |
|
|
71 { |
|
|
72 double retval; |
|
|
73 F77_XFCN (dacosh, DACOSH, (x, retval)); |
|
|
74 return retval; |
|
|
75 } |
|
|
76 #endif |
|
|
77 |
|
|
78 #if !defined (HAVE_ASINH) |
|
|
79 double |
|
|
80 asinh (double x) |
|
|
81 { |
|
|
82 double retval; |
|
|
83 F77_XFCN (dasinh, DASINH, (x, retval)); |
|
|
84 return retval; |
|
|
85 } |
|
|
86 #endif |
|
|
87 |
|
|
88 #if !defined (HAVE_ATANH) |
|
|
89 double |
|
|
90 atanh (double x) |
|
|
91 { |
|
|
92 double retval; |
|
|
93 F77_XFCN (datanh, DATANH, (x, retval)); |
|
|
94 return retval; |
|
|
95 } |
|
|
96 #endif |
|
|
97 |
|
|
98 #if !defined (HAVE_ERF) |
|
|
99 double |
|
|
100 erf (double x) |
|
|
101 { |
|
|
102 double retval; |
|
|
103 F77_XFCN (derf, DERF, (x, retval)); |
|
|
104 return retval; |
|
|
105 } |
|
|
106 #endif |
|
|
107 |
|
|
108 #if !defined (HAVE_ERFC) |
|
|
109 double |
|
|
110 erfc (double x) |
|
|
111 { |
|
|
112 double retval; |
|
|
113 F77_XFCN (derfc, DERFC, (x, retval)); |
|
|
114 return retval; |
|
|
115 } |
|
|
116 #endif |
|
|
117 |
|
|
118 #if !defined (HAVE_GAMMA) |
|
|
119 double |
|
|
120 gamma (double x) |
|
|
121 { |
|
|
122 double retval; |
|
|
123 F77_XFCN (xdgamma, XDGAMMA, (x, retval)); |
|
|
124 return retval; |
|
|
125 } |
|
|
126 #endif |
|
|
127 |
|
|
128 #if !defined (HAVE_LGAMMA) |
|
|
129 // If the system doesn't have lgamma, assume that it doesn't have |
|
|
130 // signgam either. |
|
|
131 |
|
|
132 int signgam; |
|
|
133 |
|
|
134 double |
|
|
135 lgamma (double x) |
|
|
136 { |
|
|
137 double retval; |
|
|
138 double sgngam; |
|
|
139 |
|
|
140 F77_XFCN (dlgams, DLGAMS, (x, retval, sgngam)); |
|
|
141 |
|
|
142 signgam = (int) sgngam; |
|
|
143 |
|
|
144 return retval; |
|
|
145 } |
|
|
146 #endif |
|
|
147 |
|
|
148 int |
|
|
149 F77_FCN (ribesl, RIBESL) (const double& x, const double& alpha, |
|
|
150 const int& nb, double *result, int& ncalc) |
|
|
151 { |
|
|
152 int retval = 0; |
|
|
153 F77_FCN (ribesl, RIBESL) (x, alpha, nb, 1, result, ncalc); |
|
|
154 return retval; |
|
|
155 } |
|
|
156 |
|
|
157 int |
|
|
158 F77_FCN (rkbesl, RKBESL) (const double& x, const double& alpha, |
|
|
159 const int& nb, double *result, int& ncalc) |
|
|
160 { |
|
|
161 int retval = 0; |
|
|
162 F77_FCN (rkbesl, RKBESL) (x, alpha, nb, 1, result, ncalc); |
|
|
163 return retval; |
|
|
164 } |
|
|
165 |
|
|
166 typedef int (*fptr) (const double&, const double&, const int&, double*, int&); |
|
|
167 |
|
|
168 static Matrix |
|
|
169 do_bessel (fptr f, const char *fn, double alpha, const Matrix& x) |
|
|
170 { |
|
|
171 Matrix retval; |
|
|
172 |
|
|
173 if (alpha >= 0.0) |
|
|
174 { |
|
|
175 int nb = 1; |
|
|
176 |
|
|
177 if (alpha >= 1.0) |
|
|
178 { |
|
|
179 double integral; |
|
|
180 alpha = modf (alpha, &integral); |
|
|
181 nb = static_cast<int> (integral) + 1; |
|
|
182 } |
|
|
183 |
|
|
184 Array<double> tmp (nb); |
|
|
185 |
|
|
186 double *ptmp = tmp.fortran_vec (); |
|
|
187 |
|
|
188 int x_nr = x.rows (); |
|
|
189 int x_nc = x.cols (); |
|
|
190 |
|
|
191 retval.resize (x_nr, x_nc); |
|
|
192 |
|
|
193 int ncalc; |
|
|
194 |
|
|
195 for (int j = 0; j < x_nc; j++) |
|
|
196 { |
|
|
197 for (int i = 0; i < x_nr; i++) |
|
|
198 { |
|
|
199 f (x(i,j), alpha, nb, ptmp, ncalc); |
|
|
200 |
|
|
201 if (ncalc < 0) |
|
|
202 { |
|
|
203 (*current_liboctave_error_handler) |
|
|
204 ("%s: error computing function values", fn); |
|
|
205 |
|
|
206 return Matrix (); |
|
|
207 } |
|
|
208 |
|
|
209 retval(i,j) = tmp(nb-1); |
|
|
210 } |
|
|
211 } |
|
|
212 } |
|
|
213 else |
|
|
214 (*current_liboctave_error_handler) |
|
|
215 ("%s: alpha must be greater than zero", fn); |
|
|
216 |
|
|
217 return retval; |
|
|
218 } |
|
|
219 |
|
|
220 static Matrix |
|
|
221 do_bessel (fptr f, const char *fn, const Range& alpha_range, |
|
|
222 const ColumnVector& x) |
|
|
223 { |
|
|
224 Matrix retval; |
|
|
225 |
|
|
226 double alpha_base = alpha_range.base (); |
|
|
227 |
|
|
228 if (alpha_base >= 0.0) |
|
|
229 { |
|
|
230 int nb_orig = alpha_range.nelem (); |
|
|
231 int nb = nb_orig; |
|
|
232 int offset = 0; |
|
|
233 |
|
|
234 if (alpha_base >= 1.0) |
|
|
235 { |
|
|
236 double integral; |
|
|
237 alpha_base = modf (alpha_base, &integral); |
|
|
238 offset = static_cast<int> (integral); |
|
|
239 nb += offset; |
|
|
240 } |
|
|
241 |
|
|
242 Array<double> tmp (nb); |
|
|
243 |
|
|
244 double *ptmp = tmp.fortran_vec (); |
|
|
245 |
|
|
246 int x_len = x.length (); |
|
|
247 |
|
|
248 retval.resize (x_len, nb_orig); |
|
|
249 |
|
|
250 int ncalc; |
|
|
251 |
|
|
252 for (int i = 0; i < x_len; i++) |
|
|
253 { |
|
|
254 f (x(i), alpha_base, nb, ptmp, ncalc); |
|
|
255 |
|
|
256 if (ncalc < 0) |
|
|
257 { |
|
|
258 (*current_liboctave_error_handler) |
|
|
259 ("%s: error computing function values", fn); |
|
|
260 |
|
|
261 return Matrix (); |
|
|
262 } |
|
|
263 |
|
|
264 for (int j = 0; j < nb_orig; j++) |
|
|
265 retval(i,j) = tmp(offset+j); |
|
|
266 } |
|
|
267 } |
|
|
268 else |
|
|
269 (*current_liboctave_error_handler) |
|
|
270 ("%s: alpha must be greater than zero", fn); |
|
|
271 |
|
|
272 return retval; |
|
|
273 } |
|
|
274 |
|
|
275 Matrix |
|
|
276 besselj (double alpha, const Matrix& x) |
|
|
277 { |
|
|
278 return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x); |
|
|
279 } |
|
|
280 |
|
|
281 Matrix |
|
|
282 bessely (double alpha, const Matrix& x) |
|
|
283 { |
|
|
284 return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x); |
|
|
285 } |
|
|
286 |
|
|
287 Matrix |
|
|
288 besseli (double alpha, const Matrix& x) |
|
|
289 { |
|
|
290 return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x); |
|
|
291 } |
|
|
292 |
|
|
293 Matrix |
|
|
294 besselk (double alpha, const Matrix& x) |
|
|
295 { |
|
|
296 return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x); |
|
|
297 } |
|
|
298 |
|
|
299 Matrix |
|
|
300 besselj (const Range& alpha, const ColumnVector& x) |
|
|
301 { |
|
|
302 return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x); |
|
|
303 } |
|
|
304 |
|
|
305 Matrix |
|
|
306 bessely (const Range& alpha, const ColumnVector& x) |
|
|
307 { |
|
|
308 return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x); |
|
|
309 } |
|
|
310 |
|
|
311 Matrix |
|
|
312 besseli (const Range& alpha, const ColumnVector& x) |
|
|
313 { |
|
|
314 return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x); |
|
|
315 } |
|
|
316 |
|
|
317 Matrix |
|
|
318 besselk (const Range& alpha, const ColumnVector& x) |
|
|
319 { |
|
|
320 return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x); |
|
|
321 } |
|
|
322 |
|
|
323 static void |
|
|
324 gripe_betainc_nonconformant (int r1, int c1, int r2, int c2, int r3, |
|
|
325 int c3) |
|
|
326 { |
|
|
327 (*current_liboctave_error_handler) |
|
|
328 ("betainc: nonconformant arguments (x is %dx%d, a is %dx%d, b is %dx%d)", |
|
|
329 r1, c1, r2, c2, r3, c3); |
|
|
330 } |
|
|
331 |
|
|
332 double |
|
|
333 betainc (double x, double a, double b) |
|
|
334 { |
|
|
335 double retval; |
|
|
336 F77_XFCN (xdbetai, XDBETAI, (x, a, b, retval)); |
|
|
337 return retval; |
|
|
338 } |
|
|
339 |
|
|
340 Matrix |
|
|
341 betainc (double x, double a, const Matrix& b) |
|
|
342 { |
|
|
343 int nr = b.rows (); |
|
|
344 int nc = b.cols (); |
|
|
345 |
|
|
346 Matrix retval (nr, nc); |
|
|
347 |
|
|
348 for (int j = 0; j < nc; j++) |
|
|
349 for (int i = 0; i < nr; i++) |
|
|
350 retval(i,j) = betainc (x, a, b(i,j)); |
|
|
351 |
|
|
352 return retval; |
|
|
353 } |
|
|
354 |
|
|
355 Matrix |
|
|
356 betainc (double x, const Matrix& a, double b) |
|
|
357 { |
|
|
358 int nr = a.rows (); |
|
|
359 int nc = a.cols (); |
|
|
360 |
|
|
361 Matrix retval (nr, nc); |
|
|
362 |
|
|
363 for (int j = 0; j < nc; j++) |
|
|
364 for (int i = 0; i < nr; i++) |
|
|
365 retval(i,j) = betainc (x, a(i,j), b); |
|
|
366 |
|
|
367 return retval; |
|
|
368 } |
|
|
369 |
|
|
370 Matrix |
|
|
371 betainc (double x, const Matrix& a, const Matrix& b) |
|
|
372 { |
|
|
373 Matrix retval; |
|
|
374 |
|
|
375 int a_nr = a.rows (); |
|
|
376 int a_nc = a.cols (); |
|
|
377 |
|
|
378 int b_nr = b.rows (); |
|
|
379 int b_nc = b.cols (); |
|
|
380 |
|
|
381 if (a_nr == b_nr && a_nc == b_nc) |
|
|
382 { |
|
|
383 retval.resize (a_nr, a_nc); |
|
|
384 |
|
|
385 for (int j = 0; j < a_nc; j++) |
|
|
386 for (int i = 0; i < a_nr; i++) |
|
|
387 retval(i,j) = betainc (x, a(i,j), b(i,j)); |
|
|
388 } |
|
|
389 else |
|
|
390 gripe_betainc_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc); |
|
|
391 |
|
|
392 return retval; |
|
|
393 } |
|
|
394 |
|
|
395 Matrix |
|
|
396 betainc (const Matrix& x, double a, double b) |
|
|
397 { |
|
|
398 int nr = x.rows (); |
|
|
399 int nc = x.cols (); |
|
|
400 |
|
|
401 Matrix retval (nr, nc); |
|
|
402 |
|
|
403 for (int j = 0; j < nc; j++) |
|
|
404 for (int i = 0; i < nr; i++) |
|
|
405 retval(i,j) = betainc (x(i,j), a, b); |
|
|
406 |
|
|
407 return retval; |
|
|
408 } |
|
|
409 |
|
|
410 Matrix |
|
|
411 betainc (const Matrix& x, double a, const Matrix& b) |
|
|
412 { |
|
|
413 Matrix retval; |
|
|
414 |
|
|
415 int nr = x.rows (); |
|
|
416 int nc = x.cols (); |
|
|
417 |
|
|
418 int b_nr = b.rows (); |
|
|
419 int b_nc = b.cols (); |
|
|
420 |
|
|
421 if (nr == b_nr && nc == b_nc) |
|
|
422 { |
|
|
423 retval.resize (nr, nc); |
|
|
424 |
|
|
425 for (int j = 0; j < nc; j++) |
|
|
426 for (int i = 0; i < nr; i++) |
|
|
427 retval(i,j) = betainc (x(i,j), a, b(i,j)); |
|
|
428 } |
|
|
429 else |
|
|
430 gripe_betainc_nonconformant (nr, nc, 1, 1, b_nr, b_nc); |
|
|
431 |
|
|
432 return retval; |
|
|
433 } |
|
|
434 |
|
|
435 Matrix |
|
|
436 betainc (const Matrix& x, const Matrix& a, double b) |
|
|
437 { |
|
|
438 Matrix retval; |
|
|
439 |
|
|
440 int nr = x.rows (); |
|
|
441 int nc = x.cols (); |
|
|
442 |
|
|
443 int a_nr = a.rows (); |
|
|
444 int a_nc = a.cols (); |
|
|
445 |
|
|
446 if (nr == a_nr && nc == a_nc) |
|
|
447 { |
|
|
448 retval.resize (nr, nc); |
|
|
449 |
|
|
450 for (int j = 0; j < nc; j++) |
|
|
451 for (int i = 0; i < nr; i++) |
|
|
452 retval(i,j) = betainc (x(i,j), a(i,j), b); |
|
|
453 } |
|
|
454 else |
|
|
455 gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, 1, 1); |
|
|
456 |
|
|
457 return retval; |
|
|
458 } |
|
|
459 |
|
|
460 Matrix |
|
|
461 betainc (const Matrix& x, const Matrix& a, const Matrix& b) |
|
|
462 { |
|
|
463 Matrix retval; |
|
|
464 |
|
|
465 int nr = x.rows (); |
|
|
466 int nc = x.cols (); |
|
|
467 |
|
|
468 int a_nr = a.rows (); |
|
|
469 int a_nc = a.cols (); |
|
|
470 |
|
|
471 int b_nr = b.rows (); |
|
|
472 int b_nc = b.cols (); |
|
|
473 |
|
|
474 if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc) |
|
|
475 { |
|
|
476 retval.resize (nr, nc); |
|
|
477 |
|
|
478 for (int j = 0; j < nc; j++) |
|
|
479 for (int i = 0; i < nr; i++) |
|
|
480 retval(i,j) = betainc (x(i,j), a(i,j), b(i,j)); |
|
|
481 } |
|
|
482 else |
|
|
483 gripe_betainc_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc); |
|
|
484 |
|
|
485 return retval; |
|
|
486 } |
|
|
487 |
|
|
488 double |
|
|
489 gammainc (double x, double a) |
|
|
490 { |
|
|
491 double retval; |
|
|
492 F77_XFCN (xdgami, XDGAMI, (x, a, retval)); |
|
|
493 return retval; |
|
|
494 } |
|
|
495 |
|
|
496 Matrix |
|
|
497 gammainc (double x, const Matrix& a) |
|
|
498 { |
|
|
499 int nr = a.rows (); |
|
|
500 int nc = a.cols (); |
|
|
501 |
|
|
502 Matrix retval (nr, nc); |
|
|
503 |
|
|
504 for (int j = 0; j < nc; j++) |
|
|
505 for (int i = 0; i < nr; i++) |
|
|
506 retval(i,j) = gammainc (x, a(i,j)); |
|
|
507 |
|
|
508 return retval; |
|
|
509 } |
|
|
510 |
|
|
511 Matrix |
|
|
512 gammainc (const Matrix& x, double a) |
|
|
513 { |
|
|
514 int nr = x.rows (); |
|
|
515 int nc = x.cols (); |
|
|
516 |
|
|
517 Matrix retval (nr, nc); |
|
|
518 |
|
|
519 for (int j = 0; j < nc; j++) |
|
|
520 for (int i = 0; i < nr; i++) |
|
|
521 retval(i,j) = gammainc (x(i,j), a); |
|
|
522 |
|
|
523 return retval; |
|
|
524 } |
|
|
525 |
|
|
526 Matrix |
|
|
527 gammainc (const Matrix& x, const Matrix& a) |
|
|
528 { |
|
|
529 Matrix retval; |
|
|
530 |
|
|
531 int nr = x.rows (); |
|
|
532 int nc = x.cols (); |
|
|
533 |
|
|
534 int a_nr = a.rows (); |
|
|
535 int a_nc = a.cols (); |
|
|
536 |
|
|
537 if (nr == a_nr && nc == a_nc) |
|
|
538 { |
|
|
539 retval.resize (nr, nc); |
|
|
540 |
|
|
541 for (int j = 0; j < nc; j++) |
|
|
542 for (int i = 0; i < nr; i++) |
|
|
543 retval(i,j) = gammainc (x(i,j), a(i,j)); |
|
|
544 } |
|
|
545 else |
|
|
546 (*current_liboctave_error_handler) |
|
|
547 ("gammainc: nonconformant arguments (arg 1 is %dx%d, arg 2 is %dx%d)", |
|
|
548 nr, nc, a_nr, a_nc); |
|
|
549 |
|
|
550 return retval; |
|
|
551 } |
|
|
552 |
|
|
553 /* |
|
|
554 ;;; Local Variables: *** |
|
|
555 ;;; mode: C++ *** |
|
|
556 ;;; End: *** |
|
|
557 */ |
