Mercurial > octave-antonio

comparison liboctave/numeric/lo-specfun.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
comparison
equal deleted inserted replaced
20154:45565ecec019 20155:1fae49e34a1a
3726 cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd); 3726 cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd);
3727 dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd); 3727 dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
3728 } 3728 }
3729 } 3729 }
3730 3730
3731 static const double euler_mascheroni = 0.577215664901532860606512090082402431042;
3732 static const double pi = 3.14159265358979323846;
3733 // Coefficients for C.Lanczos expansion of psi function from XLiFE++ gammaFunctions
3734 // psi_coef[k] = - (2k+1) * lg_coef[k] (see melina++ gamma functions)
3735 // -1/12, 3/360,-5/1260, 7/1680,-9/1188, 11*691/360360,-13/156, 15*3617/122400, ? , ?
3736 static const double psi_coeff[10] = {
3737 -0.83333333333333333e-1, 0.83333333333333333e-2,
3738 -0.39682539682539683e-2, 0.41666666666666667e-2,
3739 -0.75757575757575758e-2, 0.21092796092796093e-1,
3740 -0.83333333333333333e-1, 0.4432598039215686,
3741 -0.3053954330270122e+1, 0.125318899521531e+2
3742 };
3743
3731 template<class T> 3744 template<class T>
3732 T 3745 T
3733 psi (const T& z) 3746 psi (const T& z)
3734 { 3747 {
3735 const T euler_mascheroni = 0.577215664901532860606512090082402431042;
3736 const T pi = 3.14159265358979323846;
3737
3738 const bool is_int = (xfloor (z) == z); 3748 const bool is_int = (xfloor (z) == z);
3739 3749
3740 T p = 0; 3750 T p = 0;
3741 if (z <= 0) 3751 if (z <= 0)
3742 { 3752 {
3766 } 3776 }
3767 else 3777 else
3768 { 3778 {
3769 // adapted from XLiFE++ gammaFunctions 3779 // adapted from XLiFE++ gammaFunctions
3770 3780
3771 // Coefficients for C.Lanczos expansion of DiGamma function
3772 // dg_coef[k] = - (2k+1) * lg_coef[k] (see melina++ gamma functions)
3773 // -1/12, 3/360,-5/1260, 7/1680,-9/1188, 11*691/360360,-13/156, 15*3617/122400, ? , ?
3774 const T dg_coeff[10] ={-0.83333333333333333e-1, 0.83333333333333333e-2,
3775 -0.39682539682539683e-2, 0.41666666666666667e-2,
3776 -0.75757575757575758e-2, 0.21092796092796093e-1,
3777 -0.83333333333333333e-1, 0.4432598039215686,
3778 -0.3053954330270122e+1, 0.125318899521531e+2};
3779
3780 T zc = z; 3781 T zc = z;
3781 // Use formula for derivative of LogGamma(z) 3782 // Use formula for derivative of LogGamma(z)
3782 if (z < 10) 3783 if (z < 10)
3783 { 3784 {
3784 const octave_idx_type n = 10 - z; 3785 const octave_idx_type n = 10 - z;
3789 T overz2 = 1.0 / (zc*zc); 3790 T overz2 = 1.0 / (zc*zc);
3790 T overz2k = overz2; 3791 T overz2k = overz2;
3791 3792
3792 p += log (zc) - 0.5 / zc; 3793 p += log (zc) - 0.5 / zc;
3793 for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2) 3794 for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2)
3794 p += dg_coeff[k] * overz2k; 3795 p += psi_coeff[k] * overz2k;
3795 } 3796 }
3796 3797
3797 return p; 3798 return p;
3798 } 3799 }
3799 3800
3800 // explicit instantiations 3801 // explicit instantiations
3801 template double psi<double> (const double& z); 3802 template double psi<double> (const double& z);
3802 template float psi<float> (const float& z); 3803 template float psi<float> (const float& z);
3803 3804
3805 template<class T>
3806 std::complex<T>
3807 psi (const std::complex<T>& z)
3808 {
3809 // adapted from XLiFE++ gammaFunctions
3810
3811 typedef typename std::complex<T>::value_type P;
3812
3813 P z_r = z.real ();
3814 P z_ra = z_r;
3815
3816 if (z.imag () == 0)
3817 return std::complex<T> (psi (z_r), 0.0);
3818 else if (z_r < 0)
3819 return psi (P (1.0) - z)- (P (pi) / tan (P (pi) * z));
3820 else
3821 {
3822 // Use formula for derivative of LogGamma(z)
3823
3824 std::complex<T> dgam = 0.0;
3825 std::complex<T> z_p = z;
3826
3827 octave_idx_type n = 0;
3828 std::complex<T> z_m = z_p;
3829 if (z_ra < 8)
3830 {
3831 n = 8 - octave_idx_type (z_ra);
3832 z_m = z_p + std::complex<T> (n, 0.0);
3833 }
3834
3835 // for | Re(z) | > 8, use derivative of C.Lanczos expansion for LogGamma
3836 // psi(z) = log(z) - 1/(2z) - 1/12z^2 + 3/360z^4 - 5/1260z^6 + 7/1680z^8 - 9/1188z^10 + ...
3837 // (Abramowitz&Stegun, page 259, formula 6.3.18
3838 std::complex<T> overz = P (1.0) / z_m;
3839 std::complex<T> overz2 = overz * overz;
3840 std::complex<T> overz2k = overz2;
3841
3842 dgam += log (z_m) - P (0.5) * overz;
3843 for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2)
3844 dgam += P (psi_coeff[k]) * overz2k;
3845
3846 // Recurrence formula
3847 // for | Re(z) | < 8 , use recursively DiGamma(z) = DiGamma(z+1) - 1/z
3848 for (octave_idx_type k = 0; k < n; k++, z_p += 1.0)
3849 dgam -= P (1.0) / z_p;
3850
3851 return dgam;
3852 }
3853 }
3854
3855 // explicit instantiations
3856 template Complex psi<double> (const Complex& z);
3857 template FloatComplex psi<float> (const FloatComplex& z);
3858