Mercurial > hg > octave-max
changeset 3119:f3e1da120048
[project @ 1997-11-30 21:13:58 by jwe]
author | jwe |
---|---|
date | Sun, 30 Nov 1997 21:17:04 +0000 |
parents | 74cc8e2fe2c0 |
children | a2371c4a1d50 |
files | libcruft/ChangeLog libcruft/specfun/ribesl.f libcruft/specfun/rjbesl.f libcruft/specfun/rkbesl.f libcruft/specfun/rybesl.f liboctave/ChangeLog liboctave/Makefile.in liboctave/dbleBessel.cc liboctave/dbleBessel.h src/ChangeLog src/DLD-FUNCTIONS/bessel.cc src/Makefile.in |
diffstat | 12 files changed, 526 insertions(+), 10 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog +++ b/libcruft/ChangeLog @@ -1,3 +1,11 @@ +Sat Nov 29 13:02:14 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * specfun/ribesl.f (ribesl): Use d1mach to get machine parameters. + SAVE static data between calls. Use PARAMETERS where possible. + * specfun/rjbesl.f (rjbesl): Likewise. + * specfun/rkbesl.f (rkbesl): Likewise. + * specfun/rybesl.f (rybesl): Likewise. + Fri Nov 28 14:05:12 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * specfun: New subdirectory.
--- a/libcruft/specfun/ribesl.f +++ b/libcruft/specfun/ribesl.f @@ -182,6 +182,8 @@ C------------------------------------------------------------------- DATA FIRST /.TRUE./ C------------------------------------------------------------------- + SAVE FIRST, NSIG, ENTEN, ENSIG, RTNSIG, ENMTEN, EXPARG, XLARGE +C------------------------------------------------------------------- C Statement functions for conversion C------------------------------------------------------------------- CONV(N) = DBLE(N)
--- a/libcruft/specfun/rjbesl.f +++ b/libcruft/specfun/rjbesl.f @@ -173,6 +173,8 @@ C--------------------------------------------------------------------- DATA FIRST /.TRUE./ C--------------------------------------------------------------------- + SAVE FACT, FIRST, NSIG, ENTEN, ENSIG, RTNSIG, ENMTEN, XLARGE +C--------------------------------------------------------------------- C Statement functions for conversion and the gamma function. C--------------------------------------------------------------------- CONV(I) = DBLE(I)
--- a/libcruft/specfun/rkbesl.f +++ b/libcruft/specfun/rkbesl.f @@ -187,9 +187,13 @@ DATA ESTF/4.18341D1, 7.1075D0, 6.4306D0, 4.25110D1, 1.35633D0, 1 8.45096D1, 2.0D1/ C--------------------------------------------------------------------- + DATA FIRST /.TRUE./ +C--------------------------------------------------------------------- + SAVE P, Q, R, S, T, ESTM, ESTF + SAVE FIRST, EPS, XINF, XMIN, SQXMIN, XMAX +C--------------------------------------------------------------------- C Machine dependent parameters C--------------------------------------------------------------------- - DATA FIRST /.TRUE./ IF (FIRST) THEN EPS = D1MACH (4) XINF = D1MACH (2)
--- a/libcruft/specfun/rybesl.f +++ b/libcruft/specfun/rybesl.f @@ -176,9 +176,12 @@ 9 -0.76852840844786673690D-01,-0.28387654227602353814D+00, A 0.92187029365045265648D+00/ C---------------------------------------------------------------------- + DATA FIRST /.TRUE./ +C---------------------------------------------------------------------- + SAVE CH, FIRST, EPS, XINF, XMIN, DEL, XLARGE, THRESH +C---------------------------------------------------------------------- C Machine-dependent constants C---------------------------------------------------------------------- - DATA FIRST /.TRUE./ IF (FIRST) THEN EPS = D1MACH (4) XINF = D1MACH (2)
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,8 @@ +Sun Nov 30 14:59:12 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * dbleBessel.h, dbleBessel.cc: New files. + * Makefile.in (INCLUDES, SOURCES): Add them to the lists. + Wed Nov 26 20:02:13 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * lo-sysdep.cc (octave_getcwd): Prefer getcwd over getwd.
--- a/liboctave/Makefile.in +++ b/liboctave/Makefile.in @@ -42,9 +42,9 @@ INCLUDES := Bounds.h CollocWt.h DAE.h DAEFunc.h DASSL.h FEGrid.h \ LinConst.h LP.h LPsolve.h LSODE.h NLConst.h NLEqn.h NLFunc.h \ - NLP.h ODE.h ODEFunc.h Objective.h QP.h Quad.h \ - Range.h base-de.h base-min.h byte-swap.h cmd-edit.h cmd-hist.h \ - data-conv.h dir-ops.h file-ops.h file-stat.h getopt.h \ + NLP.h ODE.h ODEFunc.h Objective.h QP.h Quad.h Range.h base-de.h \ + base-min.h byte-swap.h cmd-edit.h cmd-hist.h data-conv.h \ + dbleBessel.h dir-ops.h file-ops.h file-stat.h getopt.h \ glob-match.h idx-vector.h lo-ieee.h lo-mappers.h lo-sysdep.h \ lo-utils.h mach-info.h oct-alloc.h oct-cmplx.h oct-env.h \ oct-math.h oct-group.h oct-passwd.h oct-syscalls.h pathlen.h \ @@ -76,9 +76,9 @@ mx-m-cs.cc mx-m-dm.cc mx-s-cdm.cc mx-s-cm.cc mx-s-dm.cc SOURCES := Bounds.cc CollocWt.cc DAE.cc DASSL.cc FEGrid.cc \ - LinConst.cc LPsolve.cc LSODE.cc NLEqn.cc \ - Quad.cc Range.cc acosh.c asinh.c atanh.c cmd-edit.cc \ - cmd-hist.cc data-conv.cc dir-ops.cc erf.c erfc.c f2c-main.c \ + LinConst.cc LPsolve.cc LSODE.cc NLEqn.cc Quad.cc Range.cc \ + acosh.c asinh.c atanh.c cmd-edit.cc cmd-hist.cc data-conv.cc \ + dbleBessel.cc dir-ops.cc erf.c erfc.c f2c-main.c \ file-ops.cc file-stat.cc filemode.c gamma.c getopt.c getopt1.c \ glob-match.cc idx-vector.cc lgamma.c lo-ieee.cc lo-mappers.cc \ lo-sysdep.cc lo-utils.cc mach-info.cc mkdir.c oct-alloc.cc \
new file mode 100644 --- /dev/null +++ b/liboctave/dbleBessel.cc @@ -0,0 +1,221 @@ +/* + +Copyright (C) 1996 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if defined (__GNUG__) +#pragma implementation +#endif + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "Range.h" +#include "dColVector.h" +#include "dMatrix.h" +#include "dbleBessel.h" +#include "f77-fcn.h" +#include "lo-error.h" +#include "mx-inlines.cc" + +extern "C" +{ + int F77_FCN (rjbesl, RJBESL) (const double&, const double&, + const int&, double*, int&); + + int F77_FCN (rybesl, RYBESL) (const double&, const double&, + const int&, double*, int&); + + int F77_FCN (ribesl, RIBESL) (const double&, const double&, + const int&, const int&, double*, int&); + + int F77_FCN (rkbesl, RKBESL) (const double&, const double&, + const int&, const int&, double*, int&); +} + +int +F77_FCN (ribesl, RIBESL) (const double& x, const double& alpha, + const int& nb, double *result, int& ncalc) +{ + int retval = 0; + F77_FCN (ribesl, RIBESL) (x, alpha, nb, 1, result, ncalc); + return retval; +} + +int +F77_FCN (rkbesl, RKBESL) (const double& x, const double& alpha, + const int& nb, double *result, int& ncalc) +{ + int retval = 0; + F77_FCN (rkbesl, RKBESL) (x, alpha, nb, 1, result, ncalc); + return retval; +} + +typedef int (*fptr) (const double&, const double&, const int&, double*, int&); + +static Matrix +do_bessel (fptr f, const char *fn, double alpha, const Matrix& x) +{ + Matrix retval; + + if (alpha >= 0.0) + { + int nb = 1; + + if (alpha >= 1.0) + { + double integral; + alpha = modf (alpha, &integral); + nb = static_cast<int> (integral) + 1; + } + + Array<double> tmp (nb); + + double *ptmp = tmp.fortran_vec (); + + int x_nr = x.rows (); + int x_nc = x.cols (); + + retval.resize (x_nr, x_nc); + + int ncalc; + + for (int j = 0; j < x_nc; j++) + { + for (int i = 0; i < x_nr; i++) + { + f (x(i,j), alpha, nb, ptmp, ncalc); + + // XXX FIXME XXX -- check ncalc. + + retval(i,j) = tmp(nb-1); + } + } + } + else + (*current_liboctave_error_handler) + ("%s: alpha must be greater than zero", fn); + + return retval; +} + +static Matrix +do_bessel (fptr f, const char *fn, const Range& alpha_range, + const ColumnVector& x) +{ + Matrix retval; + + double alpha_base = alpha_range.base (); + + if (alpha_base >= 0.0) + { + int nb_orig = alpha_range.nelem (); + int nb = nb_orig; + int offset = 0; + + if (alpha_base >= 1.0) + { + double integral; + alpha_base = modf (alpha_base, &integral); + offset = static_cast<int> (integral); + nb += offset; + } + + Array<double> tmp (nb); + + double *ptmp = tmp.fortran_vec (); + + int x_len = x.length (); + + retval.resize (x_len, nb_orig); + + int ncalc; + + for (int i = 0; i < x_len; i++) + { + f (x(i), alpha_base, nb, ptmp, ncalc); + + // XXX FIXME XXX -- check ncalc. + + for (int j = 0; j < nb_orig; j++) + retval(i,j) = tmp(offset+j); + } + } + else + (*current_liboctave_error_handler) + ("%s: alpha must be greater than zero", fn); + + return retval; +} + +Matrix +besselj (double alpha, const Matrix& x) +{ + return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x); +} + +Matrix +bessely (double alpha, const Matrix& x) +{ + return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x); +} + +Matrix +besseli (double alpha, const Matrix& x) +{ + return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x); +} + +Matrix +besselk (double alpha, const Matrix& x) +{ + return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x); +} + +Matrix +besselj (const Range& alpha, const ColumnVector& x) +{ + return do_bessel (F77_FCN (rjbesl, RJBESL), "besselj", alpha, x); +} + +Matrix +bessely (const Range& alpha, const ColumnVector& x) +{ + return do_bessel (F77_FCN (rybesl, RYBESL), "bessely", alpha, x); +} + +Matrix +besseli (const Range& alpha, const ColumnVector& x) +{ + return do_bessel (F77_FCN (ribesl, RIBESL), "besseli", alpha, x); +} + +Matrix +besselk (const Range& alpha, const ColumnVector& x) +{ + return do_bessel (F77_FCN (rkbesl, RKBESL), "besselk", alpha, x); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
new file mode 100644 --- /dev/null +++ b/liboctave/dbleBessel.h @@ -0,0 +1,50 @@ +/* + +Copyright (C) 1996 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#if !defined (octave_BESSEL_h) +#define octave_BESSEL_h 1 + +#if defined (__GNUG__) +#pragma interface +#endif + +class ColumnVector; +class Matrix; +class Range; + +extern Matrix besselj (double alpha, const Matrix& x); +extern Matrix bessely (double alpha, const Matrix& x); +extern Matrix besseli (double alpha, const Matrix& x); +extern Matrix besselk (double alpha, const Matrix& x); + +extern Matrix besselj (const Range& alpha, const ColumnVector& x); +extern Matrix bessely (const Range& alpha, const ColumnVector& x); +extern Matrix besseli (const Range& alpha, const ColumnVector& x); +extern Matrix besselk (const Range& alpha, const ColumnVector& x); + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +Sun Nov 30 14:58:56 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * DLD-FUNCTIONS/bessel.cc: New file. + * Makefile.in (DLD_XSRC): Add it to the list. + Thu Nov 27 23:28:59 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * lex.l (handle_string): Constructor for string class takes
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/bessel.cc @@ -0,0 +1,216 @@ +/* + +Copyright (C) 1997 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "dbleBessel.h" + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "help.h" +#include "oct-obj.h" +#include "utils.h" + +#define DO_BESSEL(type, alpha, x) \ + do \ + { \ + switch (type) \ + { \ + case 'j': \ + retval = besselj (alpha, x); \ + break; \ + \ + case 'y': \ + retval = bessely (alpha, x); \ + break; \ + \ + case 'i': \ + retval = besseli (alpha, x); \ + break; \ + \ + case 'k': \ + retval = besselk (alpha, x); \ + break; \ + \ + default: \ + break; \ + } \ + } \ + while (0) + +static void +gripe_bessel_arg_1 (const char *fn) +{ + error ("%s: alpha must be scalar or range with increment equal to 1", fn); +} + +octave_value_list +do_bessel (char type, const char *fn, const octave_value_list& args) +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 2) + { + octave_value alpha_arg = args(0); + + if (alpha_arg.is_scalar_type ()) + { + double alpha = alpha_arg.double_value (); + + if (! error_state) + { + Matrix x = args(1).matrix_value (); + + if (! error_state) + DO_BESSEL (type, alpha, x); + else + error ("%s: expecting matrix as second argument", fn); + } + else + gripe_bessel_arg_1 (fn); + } + else + { + Range alpha; + + if (! alpha_arg.is_range ()) + { + ColumnVector tmp = alpha_arg.vector_value (); + + if (! error_state) + { + int len = tmp.length (); + + double base = tmp(0); + + for (int i = 1; i < len; i++) + { + if (tmp(i) != base + i) + { + gripe_bessel_arg_1 (fn); + break; + } + } + + if (! error_state) + alpha = Range (tmp(0), tmp(len-1)); + } + } + else + alpha = alpha_arg.range_value (); + + if (! error_state) + { + ColumnVector x = args(1).vector_value (); + + if (! error_state) + DO_BESSEL (type, alpha, x); + else + error ("%s: expecting vector as second argument", fn); + } + } + } + else + print_usage (fn); + + return retval; +} + +DEFUN_DLD (besselj, args, , + "besselj (alpha, x)\n\ +\n\ +Compute Bessel functions of the first kind.\n\ +\n\ +X must be a real matrix, vector or scalar.\n\ +\n\ +If ALPHA is a scalar, the result is the same size as X. If ALPHA is a\n\ +range, X must be a vector or scalar, and the result is a matrix with\n\ +length(X) rows and length(ALPHA) columns.\n\ +\n\ +ALPHA must be greater than or equal to zero. If ALPHA is a range, it\n\ +must have an increment equal to one.") +{ + return do_bessel ('j', "besselj", args); +} + +DEFUN_DLD (bessely, args, , + "bessely (alpha, x)\n\ +\n\ +Compute Bessel functions of the second kind.\n\ +\n\ +X must be a real matrix, vector or scalar.\n\ +\n\ +If ALPHA is a scalar, the result is the same size as X. If ALPHA is a\n\ +range, X must be a vector or scalar, and the result is a matrix with\n\ +length(X) rows and length(ALPHA) columns.\n\ +\n\ +ALPHA must be greater than or equal to zero. If ALPHA is a range, it\n\ +must have an increment equal to one.") +{ + return do_bessel ('y', "bessely", args); +} + +DEFUN_DLD (besseli, args, , + "besseli (alpha, x)\n\ +\n\ +Compute modified Bessel functions of the first kind.\n\ +\n\ +X must be a real matrix, vector or scalar.\n\ +\n\ +If ALPHA is a scalar, the result is the same size as X. If ALPHA is a\n\ +range, X must be a vector or scalar, and the result is a matrix with\n\ +length(X) rows and length(ALPHA) columns.\n\ +\n\ +ALPHA must be greater than or equal to zero. If ALPHA is a range, it\n\ +must have an increment equal to one.") +{ + return do_bessel ('i', "besseli", args); +} + +DEFUN_DLD (besselk, args, , + "besselk (alpha, x)\n\ +\n\ +Compute modified Bessel functions of the second kind.\n\ +\n\ +X must be a real matrix, vector or scalar.\n\ +\n\ +If ALPHA is a scalar, the result is the same size as X. If ALPHA is a\n\ +range, X must be a vector or scalar, and the result is a matrix with\n\ +length(X) rows and length(ALPHA) columns.\n\ +\n\ +ALPHA must be greater than or equal to zero. If ALPHA is a range, it\n\ +must have an increment equal to one.") +{ + return do_bessel ('k', "besselk", args); +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -39,8 +39,8 @@ endif endif -DLD_XSRC := balance.cc chol.cc colloc.cc dassl.cc det.cc eig.cc \ - expm.cc fft.cc fft2.cc filter.cc find.cc fsolve.cc \ +DLD_XSRC := balance.cc bessel.cc chol.cc colloc.cc dassl.cc det.cc \ + eig.cc expm.cc fft.cc fft2.cc filter.cc find.cc fsolve.cc \ getgrent.cc getpwent.cc getrusage.cc givens.cc hess.cc \ ifft.cc ifft2.cc inv.cc log.cc lpsolve.cc lsode.cc \ lu.cc minmax.cc pinv.cc qr.cc quad.cc qzval.cc rand.cc \