Mercurial > octave

diff liboctave/numeric/lo-specfun.cc @ 24906:451f4f291f46 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 740df3fce3fb
children bd89440407aa
line wrap: on
line diff
--- 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