Mercurial > octave
changeset 24906:451f4f291f46 stable
Modified Bessel functions to compute the output with any input
(instead of returning NaN(, giving IERR=4 (see bug #48316)
--
changed libinterp/corefcn/besselj.cc
changed liboctave/external/amos/README
changed liboctave/external/amos/cbesh.f
changed liboctave/external/amos/cbesi.f
changed liboctave/external/amos/cbesj.f
changed liboctave/external/amos/cbesk.f
changed liboctave/external/amos/zbesh.f
changed liboctave/external/amos/zbesi.f
changed liboctave/external/amos/zbesj.f
changed liboctave/external/amos/zbesk.f
changed liboctave/numeric/lo-specfun.cc
changed scripts/specfun/bessel.m
| author | Michele Ginesi <michele.ginesi@gmail.com> |
|---|---|
| date | Thu, 07 Sep 2017 15:58:26 +0200 |
| parents | 662faf9de127 |
| children | bd89440407aa |
| files | libinterp/corefcn/besselj.cc liboctave/external/amos/README liboctave/external/amos/cbesh.f liboctave/external/amos/cbesi.f liboctave/external/amos/cbesj.f liboctave/external/amos/cbesk.f liboctave/external/amos/zbesh.f liboctave/external/amos/zbesi.f liboctave/external/amos/zbesj.f liboctave/external/amos/zbesk.f liboctave/numeric/lo-specfun.cc scripts/help/bessel.m |
| diffstat | 12 files changed, 51 insertions(+), 147 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/besselj.cc Thu Sep 07 15:49:57 2017 +0200 +++ b/libinterp/corefcn/besselj.cc Thu Sep 07 15:58:26 2017 +0200 @@ -335,7 +335,7 @@ accuracy. @item -Complete loss of significance by argument reduction, return @code{NaN}. +Loss of significance by argument reduction, output may be inaccurate. @item Error---no computation, algorithm termination condition not met, return @@ -648,7 +648,7 @@ of machine accuracy. @item -Complete loss of significance by argument reduction, return @code{NaN}. +Loss of significance by argument reduction, output may be inaccurate. @item Error---no computation, algorithm termination condition not met,
--- a/liboctave/external/amos/README Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/README Thu Sep 07 15:58:26 2017 +0200 @@ -13,3 +13,21 @@ jwe@octave.org Wed Nov 11 17:29:50 1998 + +More files have been modified to recover Matlab compatibility for +Bessel functions. Now the output is always computed, independently +from the magnitude of the input + + cbesh.f + cbesi.f + cbesj.f + cbesk.f + zbesh.f + zbesi.f + zbesj.f + zbesk.f + +Michele Ginesi +michele.ginesi@gmail.com + +Sat Jul 22 11:43:40 2017
--- a/liboctave/external/amos/cbesh.f Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/cbesh.f Thu Sep 07 15:58:26 2017 +0200 @@ -219,6 +219,7 @@ C----------------------------------------------------------------------- C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- + 35 CONTINUE UFL = R1MACH(1)*1.0E+3 IF (AZ.LT.UFL) GO TO 220 IF (FNU.GT.FNUL) GO TO 90 @@ -325,7 +326,6 @@ IERR=5 RETURN 240 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/cbesi.f Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/cbesi.f Thu Sep 07 15:58:26 2017 +0200 @@ -201,6 +201,7 @@ AA=SQRT(AA) IF(AZ.GT.AA) IERR=3 IF(FN.GT.AA) IERR=3 + 35 CONTINUE ZN = Z CSGN = CONE IF (XX.GE.0.0E0) GO TO 40 @@ -252,7 +253,6 @@ IERR=5 RETURN 140 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/cbesj.f Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/cbesj.f Thu Sep 07 15:58:26 2017 +0200 @@ -199,6 +199,7 @@ C CALCULATE CSGN=EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE C WHEN FNU IS LARGE C----------------------------------------------------------------------- + 35 CONTINUE INU = INT(FNU) INUH = INU/2 IR = INU - 2*INUH @@ -247,7 +248,6 @@ IERR=5 RETURN 140 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/cbesk.f Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/cbesk.f Thu Sep 07 15:58:26 2017 +0200 @@ -202,6 +202,7 @@ C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- C UFL = EXP(-ELIM) + 35 CONTINUE UFL = R1MACH(1)*1.0E+3 IF (AZ.LT.UFL) GO TO 180 IF (FNU.GT.FNUL) GO TO 80 @@ -270,7 +271,6 @@ IERR=5 RETURN 210 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/zbesh.f Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/zbesh.f Thu Sep 07 15:58:26 2017 +0200 @@ -220,6 +220,7 @@ C----------------------------------------------------------------------- C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- + 35 CONTINUE UFL = D1MACH(1)*1.0D+3 IF (AZ.LT.UFL) GO TO 230 IF (FNU.GT.FNUL) GO TO 90 @@ -342,7 +343,6 @@ IERR=5 RETURN 260 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/zbesi.f Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/zbesi.f Thu Sep 07 15:58:26 2017 +0200 @@ -202,6 +202,7 @@ AA = DSQRT(AA) IF (AZ.GT.AA) IERR=3 IF (FN.GT.AA) IERR=3 + 35 CONTINUE ZNR = ZR ZNI = ZI CSGNR = CONER @@ -263,7 +264,6 @@ IERR=5 RETURN 260 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/external/amos/zbesj.f Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/zbesj.f Thu Sep 07 15:58:26 2017 +0200 @@ -200,6 +200,7 @@ C CALCULATE CSGN=EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE C WHEN FNU IS LARGE C----------------------------------------------------------------------- + 35 CONTINUE CII = 1.0D0 INU = INT(SNGL(FNU)) INUH = INU/2 @@ -260,7 +261,7 @@ IERR=5 RETURN 260 CONTINUE - NZ=0 IERR=4 + GO TO 35 RETURN END
--- a/liboctave/external/amos/zbesk.f Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/external/amos/zbesk.f Thu Sep 07 15:58:26 2017 +0200 @@ -205,6 +205,7 @@ C OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE C----------------------------------------------------------------------- C UFL = DEXP(-ELIM) + 35 CONTINUE UFL = D1MACH(1)*1.0D+3 IF (AZ.LT.UFL) GO TO 180 IF (FNU.GT.FNUL) GO TO 80 @@ -275,7 +276,6 @@ IERR=5 RETURN 260 CONTINUE - NZ=0 IERR=4 - RETURN + GO TO 35 END
--- a/liboctave/numeric/lo-specfun.cc Thu Sep 07 15:49:57 2017 +0200 +++ b/liboctave/numeric/lo-specfun.cc Thu Sep 07 15:58:26 2017 +0200 @@ -77,6 +77,7 @@ { case 0: case 3: + case 4: retval = val; break; @@ -109,6 +110,7 @@ { case 0: case 3: + case 4: retval = val; break; @@ -304,17 +306,10 @@ double zr = z.real (); double zi = z.imag (); - F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, t_ierr); + F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - double expz = exp (std::abs (zi)); - yr *= expz; - yi *= expz; - } - if (zi == 0.0 && zr >= 0.0) yi = 0.0; @@ -375,18 +370,11 @@ } else { - F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, 2, 1, &yr, &yi, nz, + F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, kode, 1, &yr, &yi, nz, &wr, &wi, t_ierr); ierr = t_ierr; - if (kode != 2) - { - double expz = exp (std::abs (zi)); - yr *= expz; - yi *= expz; - } - if (zi == 0.0 && zr >= 0.0) yi = 0.0; } @@ -437,17 +425,10 @@ double zr = z.real (); double zi = z.imag (); - F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 2, 1, &yr, &yi, nz, t_ierr); + F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - double expz = exp (std::abs (zr)); - yr *= expz; - yi *= expz; - } - if (zi == 0.0 && zr >= 0.0) yi = 0.0; @@ -513,24 +494,11 @@ } else { - F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, 2, 1, &yr, &yi, nz, + F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, kode, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - Complex expz = exp (-z); - - double rexpz = expz.real (); - double iexpz = expz.imag (); - - double tmp = yr*rexpz - yi*iexpz; - - yi = yr*iexpz + yi*rexpz; - yr = tmp; - } - if (zi == 0.0 && zr >= 0.0) yi = 0.0; } @@ -562,24 +530,11 @@ double zr = z.real (); double zi = z.imag (); - F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz, + F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 1, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - Complex expz = exp (Complex (0.0, 1.0) * z); - - double rexpz = expz.real (); - double iexpz = expz.imag (); - - double tmp = yr*rexpz - yi*iexpz; - - yi = yr*iexpz + yi*rexpz; - yr = tmp; - } - retval = bessel_return_value (Complex (yr, yi), ierr); } else @@ -611,24 +566,11 @@ double zr = z.real (); double zi = z.imag (); - F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz, + F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 2, 1, &yr, &yi, nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - Complex expz = exp (-Complex (0.0, 1.0) * z); - - double rexpz = expz.real (); - double iexpz = expz.imag (); - - double tmp = yr*rexpz - yi*iexpz; - - yi = yr*iexpz + yi*rexpz; - yr = tmp; - } - retval = bessel_return_value (Complex (yr, yi), ierr); } else @@ -922,17 +864,11 @@ F77_INT nz, t_ierr; - F77_FUNC (cbesj, CBESJ) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, + F77_FUNC (cbesj, CBESJ) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - float expz = exp (std::abs (z.imag ())); - y *= expz; - } - if (z.imag () == 0.0 && z.real () >= 0.0) y = FloatComplex (y.real (), 0.0); @@ -990,18 +926,12 @@ } else { - F77_FUNC (cbesy, CBESY) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, + F77_FUNC (cbesy, CBESY) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, F77_CMPLX_ARG (&y), nz, F77_CMPLX_ARG (&w), t_ierr); ierr = t_ierr; - if (kode != 2) - { - float expz = exp (std::abs (z.imag ())); - y *= expz; - } - if (z.imag () == 0.0 && z.real () >= 0.0) y = FloatComplex (y.real (), 0.0); } @@ -1050,17 +980,11 @@ F77_INT nz, t_ierr; - F77_FUNC (cbesi, CBESI) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, + F77_FUNC (cbesi, CBESI) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - float expz = exp (std::abs (z.real ())); - y *= expz; - } - if (z.imag () == 0.0 && z.real () >= 0.0) y = FloatComplex (y.real (), 0.0); @@ -1115,24 +1039,11 @@ } else { - F77_FUNC (cbesk, CBESK) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, + F77_FUNC (cbesk, CBESK) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - FloatComplex expz = exp (-z); - - float rexpz = expz.real (); - float iexpz = expz.imag (); - - float tmp_r = y.real () * rexpz - y.imag () * iexpz; - float tmp_i = y.real () * iexpz + y.imag () * rexpz; - - y = FloatComplex (tmp_r, tmp_i); - } - if (z.imag () == 0.0 && z.real () >= 0.0) y = FloatComplex (y.real (), 0.0); } @@ -1160,24 +1071,11 @@ F77_INT nz, t_ierr; - F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 1, 1, + F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 1, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - FloatComplex expz = exp (FloatComplex (0.0, 1.0) * z); - - float rexpz = expz.real (); - float iexpz = expz.imag (); - - float tmp_r = y.real () * rexpz - y.imag () * iexpz; - float tmp_i = y.real () * iexpz + y.imag () * rexpz; - - y = FloatComplex (tmp_r, tmp_i); - } - retval = bessel_return_value (y, ierr); } else @@ -1206,24 +1104,11 @@ F77_INT nz, t_ierr; - F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, 2, 2, 1, + F77_FUNC (cbesh, CBESH) (F77_CONST_CMPLX_ARG (&z), alpha, kode, 2, 1, F77_CMPLX_ARG (&y), nz, t_ierr); ierr = t_ierr; - if (kode != 2) - { - FloatComplex expz = exp (-FloatComplex (0.0, 1.0) * z); - - float rexpz = expz.real (); - float iexpz = expz.imag (); - - float tmp_r = y.real () * rexpz - y.imag () * iexpz; - float tmp_i = y.real () * iexpz + y.imag () * rexpz; - - y = FloatComplex (tmp_r, tmp_i); - } - retval = bessel_return_value (y, ierr); } else
--- a/scripts/help/bessel.m Thu Sep 07 15:49:57 2017 +0200 +++ b/scripts/help/bessel.m Thu Sep 07 15:58:26 2017 +0200 @@ -76,7 +76,7 @@ ## machine accuracy. ## ## @item -## Complete loss of significance by argument reduction, return @code{NaN}. +## Loss of significance by argument reduction, output may be inaccurate. ## ## @item ## Error---no computation, algorithm termination condition not met, return