Mercurial > octave-antonio

diff libinterp/corefcn/psi.cc @ 20155:1fae49e34a1a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
psi: add support for complex numbers. * libinterp/corefcn/psi.cc: add logic and input check to support complex numbers (implementation is in lo-specfun.cc). Add tests. * liboctave/numeric/lo-specfun.cc, liboctave/numeric/lo-specfun.h: add template specialization to psi() for std::complex. It is mostly taken from the implementation in XLiFE++ (also under GPLv3+, see b03c7cccadc2 commit message for more details).
author Carnë Draug <carandraug@octave.org>
date Sun, 15 Mar 2015 06:30:09 +0000
parents 45565ecec019
children bd565f3e0ecb
line wrap: on
line diff
--- a/libinterp/corefcn/psi.cc	Sun Mar 15 03:31:16 2015 +0000
+++ b/libinterp/corefcn/psi.cc	Sun Mar 15 06:30:09 2015 +0000
@@ -73,16 +73,28 @@
       E* psi_zv = psi_z.fortran_vec (); \
       const octave_idx_type n = z.numel (); \
       for (octave_idx_type i = 0; i < n; i++) \
-        psi_zv[i] = psi<E> (zv[i]); \
+        psi_zv[i] = psi (zv[i]); \
 \
       retval = psi_z; \
     } 
 
-  FLOAT_BRANCH(double, , , double)
-  else FLOAT_BRANCH(single, Float, float_, float)
+  if (args(0).is_complex_type ())
+    {
+      FLOAT_BRANCH(double, Complex, complex_, Complex)
+      else FLOAT_BRANCH(single, FloatComplex, float_complex_, FloatComplex)
+      else
+        {
+          error ("psi: Z must be a floating point");
+        }
+    }
   else
     {
-      error ("psi: Z must be a floating point");
+      FLOAT_BRANCH(double, , , double)
+      else FLOAT_BRANCH(single, Float, float_, float)
+      else
+        {
+          error ("psi: Z must be a floating point");
+        }
     }
 
 #undef FLOAT_BRANCH
@@ -126,5 +138,22 @@
 %!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10)
 %!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps)
 %!assert (psi (-4.653237761743142441714598151148207363719069416133868), 0, eps*100)
+
+## Tests for complex values
+%!shared z
+%! z = [-10:.1:-.1 .1:.1:20]; # drop the 0
+
+## Abramowitz and Stegun, page 259 eq 6.3.10
+%!assert (real (psi (i*z)), real (psi (1 - i*z)))
+
+## Abramowitz and Stegun, page 259 eq 6.3.11
+%!assert (imag (psi (i*z)), 1/2 .* 1./z + 1/2 * pi * coth (pi * z), eps *10)
+
+## Abramowitz and Stegun, page 259 eq 6.3.12
+%!assert (imag (psi (1/2 + i*z)), 1/2 * pi * tanh (pi * z), eps*10)
+
+## Abramowitz and Stegun, page 259 eq 6.3.13
+%!assert (imag (psi (1 + i*z)), - 1./(2*z) + 1/2 * pi * coth (pi * z), eps*10)
+
 */