Mercurial > octave
changeset 20176:5212cda4ac0e
maint: Periodic merge of stable to default.
| author | Rik <rik@octave.org> |
|---|---|
| date | Sun, 10 May 2015 21:37:19 -0700 |
| parents | c0f64bc26eee (diff) 014edaafa3ad (current diff) |
| children | 89616a98b02e |
| files | libinterp/corefcn/strfns.cc m4/acinclude.m4 scripts/image/gmap40.m scripts/image/imformats.m scripts/strings/cstrcat.m scripts/strings/strcat.m |
| diffstat | 14 files changed, 1300 insertions(+), 70 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Sun May 10 21:35:59 2015 -0700 +++ b/NEWS Sun May 10 21:37:19 2015 -0700 @@ -1,3 +1,11 @@ +Summary of important user-visible changes for version 4.2: +--------------------------------------------------------- + + ** Other new functions added in 4.2: + + psi + + Summary of important user-visible changes for version 4.0: ---------------------------------------------------------
--- a/doc/interpreter/arith.txi Sun May 10 21:35:59 2015 -0700 +++ b/doc/interpreter/arith.txi Sun May 10 21:37:19 2015 -0700 @@ -315,6 +315,8 @@ @anchor{XREFgammaln} @DOCSTRING(lgamma) +@DOCSTRING(psi) + @node Rational Approximations @section Rational Approximations
--- a/libinterp/corefcn/module.mk Sun May 10 21:35:59 2015 -0700 +++ b/libinterp/corefcn/module.mk Sun May 10 21:37:19 2015 -0700 @@ -233,6 +233,7 @@ corefcn/pr-output.cc \ corefcn/procstream.cc \ corefcn/profiler.cc \ + corefcn/psi.cc \ corefcn/quad.cc \ corefcn/quadcc.cc \ corefcn/qz.cc \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/psi.cc Sun May 10 21:37:19 2015 -0700 @@ -0,0 +1,240 @@ +/* + +Copyright (C) 2015 Carnë Draug + +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 3 of the License, 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, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "ov.h" +#include "defun.h" +#include "error.h" +#include "dNDArray.h" +#include "fNDArray.h" + +#include "lo-specfun.h" + +DEFUN (psi, args, , +"-*- texinfo -*-\n\ +@deftypefn {Function File} {} psi (@var{z})\n\ +@deftypefnx {Function File} {} psi (@var{k}, @var{z})\n\ +Compute the psi (polygamma) function.\n\ +\n\ +The polygamma functions are the @var{k}th derivative of the logarithm\n\ +of the gamma function. If unspecified, @var{k} defaults to zero. A value\n\ +of zero computes the digamma function, a value of 1, the trigamma function,\n\ +and so on.\n\ +\n\ +The digamma function is defined:\n\ +\n\ +@tex\n\ +$$\n\ +\\Psi (z) = {d (log (\\Gamma (z))) \\over dx}\n\ +$$\n\ +@end tex\n\ +@ifnottex\n\ +@example\n\ +@group\n\ +psi (z) = d (log (gamma (z))) / dx\n\ +@end group\n\ +@end example\n\ +@end ifnottex\n\ +\n\ +When computing the digamma function (when @var{k} equals zero), @var{z}\n\ +can have any value real or complex value. However, for polygamma functions\n\ +(@var{k} higher than 0), @var{z} must be real and non-negative.\n\ +\n\ +@seealso{gamma, gammainc, gammaln}\n\ +@end deftypefn") +{ + octave_value retval; + + const octave_idx_type nargin = args.length (); + if (nargin < 1 || nargin > 2) + { + print_usage (); + return retval; + } + + const octave_value oct_z = (nargin == 1) ? args(0) : args(1); + const octave_idx_type k = (nargin == 1) ? 0 : args(0).idx_type_value (); + if (error_state || k < 0) + { + error ("psi: K must be a non-negative integer"); + return retval; + } + else if (k == 0) + { +#define FLOAT_BRANCH(T, A, M, E) \ + if (oct_z.is_ ## T ##_type ()) \ + { \ + const A ## NDArray z = oct_z.M ## array_value (); \ + A ## NDArray psi_z (z.dims ()); \ +\ + const E* zv = z.data (); \ + 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++ = psi (*zv++); \ +\ + retval = psi_z; \ + } + + if (oct_z.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 + { + FLOAT_BRANCH(double, , , double) + else FLOAT_BRANCH(single, Float, float_, float) + else + { + error ("psi: Z must be a floating point"); + } + } + +#undef FLOAT_BRANCH + } + else + { + if (! oct_z.is_real_type ()) + { + error ("psi: Z must be real value for polygamma (K > 0)"); + return retval; + } + +#define FLOAT_BRANCH(T, A, M, E) \ + if (oct_z.is_ ## T ##_type ()) \ + { \ + const A ## NDArray z = oct_z.M ## array_value (); \ + A ## NDArray psi_z (z.dims ()); \ +\ + const E* zv = z.data (); \ + E* psi_zv = psi_z.fortran_vec (); \ + const octave_idx_type n = z.numel (); \ + for (octave_idx_type i = 0; i < n; i++) \ + { \ + if (*zv < 0) \ + { \ + error ("psi: Z must be non-negative for polygamma (K > 0)"); \ + return retval; \ + } \ + *psi_zv++ = psi (k, *zv++); \ + } \ + retval = psi_z; \ + } + + FLOAT_BRANCH(double, , , double) + else FLOAT_BRANCH(single, Float, float_, float) + else + { + error ("psi: Z must be a floating point for polygamma (K > 0)"); + } + +#undef FLOAT_BRANCH + } + + return retval; +} + +/* +%!shared em +%! em = 0.577215664901532860606512090082402431042; # Euler-Mascheroni Constant + +%!assert (psi (ones (7, 3, 5)), repmat (-em, [7 3 5])) +%!assert (psi ([0 1]), [-Inf -em]) +%!assert (psi ([-20:1]), [repmat(-Inf, [1 21]) -em]) +%!assert (psi (single ([0 1])), single ([-Inf -em])) + +## Abramowitz and Stegun, page 258, eq 6.3.5 +%!test +%! z = [-10:.1:-.1 .1:.1:20]; # drop the 0 +%! assert (psi (z + 1), psi (z) + 1 ./ z, eps*1000) + +## Abramowitz and Stegun, page 258, eq 6.3.2 +%!assert (psi (1), -em) + +## Abramowitz and Stegun, page 258, eq 6.3.3 +%!assert (psi (1/2), -em - 2 * log (2)) + +## The following tests are from Pascal Sebah and Xavier Gourdon (2002) +## "Introduction to the Gamma Function" + +## Interesting identities of the digamma function, in section of 5.1.3 +%!assert (psi (1/3), - em - (3/2) * log(3) - ((sqrt (3) / 6) * pi), eps*10) +%!assert (psi (1/4), - em -3 * log (2) - pi /2) +%!assert (psi (1/6), - em -2 * log (2) - (3/2) * log (3) - ((sqrt (3) / 2) * pi), eps*10) + +## First 6 zeros of the digamma function, in section of 5.1.5 (and also on +## Abramowitz and Stegun, page 258, eq 6.3.19) +%!assert (psi ( 1.46163214496836234126265954232572132846819620400644), 0, eps) +%!assert (psi (-0.504083008264455409258269304533302498955385182368579), 0, eps) +%!assert (psi (-1.573498473162390458778286043690434612655040859116846), 0, eps) +%!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10) +%!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps*10) +%!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) + +## 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) + +## Abramowitz and Stegun, page 260 eq 6.4.5 +%!test +%! for z = 0:20 +%! assert (psi (1, z + 0.5), 0.5 * (pi^2) - 4 * sum ((2*(1:z) -1) .^(-2)), eps*10) +%! endfor + +## Abramowitz and Stegun, page 260 eq 6.4.6 +%!test +%! z = 0.1:0.1:20; +%! for n = 0:8 +%! ## our precision goes down really quick when computing n is too high, +%! assert (psi (n, z+1), psi (n, z) + ((-1)^n) * factorial (n) * (z.^(-n-1)), 0.1) +%! endfor + +## Test input validation +%!error psi () +%!error psi (1, 2, 3) +%!error <Z must be> psi ("non numeric") +%!error <K must be a non-negative integer> psi (-5, 1) +%!error <Z must be non-negative for polygamma> psi (5, -1) +%!error <Z must be a floating point> psi (5, uint8 (-1)) +%!error <Z must be real value for polygamma> psi (5, 5i) + +*/
--- a/libinterp/corefcn/strfns.cc Sun May 10 21:35:59 2015 -0700 +++ b/libinterp/corefcn/strfns.cc Sun May 10 21:37:19 2015 -0700 @@ -210,78 +210,70 @@ octave_value retval; int nargin = args.length (); - - if (nargin > 0) - { - int n_elts = 0; - - size_t max_len = 0; - - std::queue<string_vector> args_as_strings; + int n_elts = 0; + size_t max_len = 0; + std::queue<string_vector> args_as_strings; - for (int i = 0; i < nargin; i++) - { - string_vector s = args(i).all_strings (); - - if (error_state) - { - error ("strvcat: unable to convert some args to strings"); - return retval; - } + for (int i = 0; i < nargin; i++) + { + string_vector s = args(i).all_strings (); - size_t n = s.length (); + if (error_state) + { + error ("strvcat: unable to convert some args to strings"); + return retval; + } - // do not count empty strings in calculation of number of elements - if (n > 0) + size_t n = s.length (); + + // do not count empty strings in calculation of number of elements + if (n > 0) + { + for (size_t j = 0; j < n; j++) { - for (size_t j = 0; j < n; j++) - { - if (s[j].length () > 0) - n_elts++; - } + if (s[j].length () > 0) + n_elts++; } - - size_t s_max_len = s.max_length (); - - if (s_max_len > max_len) - max_len = s_max_len; - - args_as_strings.push (s); } - string_vector result (n_elts); + size_t s_max_len = s.max_length (); + + if (s_max_len > max_len) + max_len = s_max_len; + + args_as_strings.push (s); + } + + string_vector result (n_elts); - octave_idx_type k = 0; + octave_idx_type k = 0; - for (int i = 0; i < nargin; i++) + for (int i = 0; i < nargin; i++) + { + string_vector s = args_as_strings.front (); + args_as_strings.pop (); + + size_t n = s.length (); + + if (n > 0) { - string_vector s = args_as_strings.front (); - args_as_strings.pop (); - - size_t n = s.length (); - - if (n > 0) + for (size_t j = 0; j < n; j++) { - for (size_t j = 0; j < n; j++) + std::string t = s[j]; + if (t.length () > 0) { - std::string t = s[j]; - if (t.length () > 0) - { - size_t t_len = t.length (); + size_t t_len = t.length (); - if (max_len > t_len) - t += std::string (max_len - t_len, ' '); + if (max_len > t_len) + t += std::string (max_len - t_len, ' '); - result[k++] = t; - } + result[k++] = t; } } } + } - retval = octave_value (result, '\''); - } - else - print_usage (); + retval = octave_value (result, '\''); return retval; } @@ -298,8 +290,7 @@ %!assert (strvcat ({100,{100, {""}}}), ["d";"d"]) %!assert (strvcat (["a";"be"], {"c", 100}), ["a";"be";"c";"d"]) %!assert (strvcat ("a", "bb", "ccc"), ["a "; "bb "; "ccc"]) - -%!error strvcat () +%!assert (strvcat (), "") */
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/cruft/slatec-fn/dpsifn.f Sun May 10 21:37:19 2015 -0700 @@ -0,0 +1,368 @@ +*DECK DPSIFN + SUBROUTINE DPSIFN (X, N, KODE, M, ANS, NZ, IERR) +C***BEGIN PROLOGUE DPSIFN +C***PURPOSE Compute derivatives of the Psi function. +C***LIBRARY SLATEC +C***CATEGORY C7C +C***TYPE DOUBLE PRECISION (PSIFN-S, DPSIFN-D) +C***KEYWORDS DERIVATIVES OF THE GAMMA FUNCTION, POLYGAMMA FUNCTION, +C PSI FUNCTION +C***AUTHOR Amos, D. E., (SNLA) +C***DESCRIPTION +C +C The following definitions are used in DPSIFN: +C +C Definition 1 +C PSI(X) = d/dx (ln(GAMMA(X)), the first derivative of +C the log GAMMA function. +C Definition 2 +C K K +C PSI(K,X) = d /dx (PSI(X)), the K-th derivative of PSI(X). +C ___________________________________________________________________ +C DPSIFN computes a sequence of SCALED derivatives of +C the PSI function; i.e. for fixed X and M it computes +C the M-member sequence +C +C ((-1)**(K+1)/GAMMA(K+1))*PSI(K,X) +C for K = N,...,N+M-1 +C +C where PSI(K,X) is as defined above. For KODE=1, DPSIFN returns +C the scaled derivatives as described. KODE=2 is operative only +C when K=0 and in that case DPSIFN returns -PSI(X) + LN(X). That +C is, the logarithmic behavior for large X is removed when KODE=2 +C and K=0. When sums or differences of PSI functions are computed +C the logarithmic terms can be combined analytically and computed +C separately to help retain significant digits. +C +C Note that CALL DPSIFN(X,0,1,1,ANS) results in +C ANS = -PSI(X) +C +C Input X is DOUBLE PRECISION +C X - Argument, X .gt. 0.0D0 +C N - First member of the sequence, 0 .le. N .le. 100 +C N=0 gives ANS(1) = -PSI(X) for KODE=1 +C -PSI(X)+LN(X) for KODE=2 +C KODE - Selection parameter +C KODE=1 returns scaled derivatives of the PSI +C function. +C KODE=2 returns scaled derivatives of the PSI +C function EXCEPT when N=0. In this case, +C ANS(1) = -PSI(X) + LN(X) is returned. +C M - Number of members of the sequence, M.ge.1 +C +C Output ANS is DOUBLE PRECISION +C ANS - A vector of length at least M whose first M +C components contain the sequence of derivatives +C scaled according to KODE. +C NZ - Underflow flag +C NZ.eq.0, A normal return +C NZ.ne.0, Underflow, last NZ components of ANS are +C set to zero, ANS(M-K+1)=0.0, K=1,...,NZ +C IERR - Error flag +C IERR=0, A normal return, computation completed +C IERR=1, Input error, no computation +C IERR=2, Overflow, X too small or N+M-1 too +C large or both +C IERR=3, Error, N too large. Dimensioned +C array TRMR(NMAX) is not large enough for N +C +C The nominal computational accuracy is the maximum of unit +C roundoff (=D1MACH(4)) and 1.0D-18 since critical constants +C are given to only 18 digits. +C +C PSIFN is the single precision version of DPSIFN. +C +C *Long Description: +C +C The basic method of evaluation is the asymptotic expansion +C for large X.ge.XMIN followed by backward recursion on a two +C term recursion relation +C +C W(X+1) + X**(-N-1) = W(X). +C +C This is supplemented by a series +C +C SUM( (X+K)**(-N-1) , K=0,1,2,... ) +C +C which converges rapidly for large N. Both XMIN and the +C number of terms of the series are calculated from the unit +C roundoff of the machine environment. +C +C***REFERENCES Handbook of Mathematical Functions, National Bureau +C of Standards Applied Mathematics Series 55, edited +C by M. Abramowitz and I. A. Stegun, equations 6.3.5, +C 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964. +C D. E. Amos, A portable Fortran subroutine for +C derivatives of the Psi function, Algorithm 610, ACM +C Transactions on Mathematical Software 9, 4 (1983), +C pp. 494-502. +C***ROUTINES CALLED D1MACH, I1MACH +C***REVISION HISTORY (YYMMDD) +C 820601 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890911 Removed unnecessary intrinsics. (WRB) +C 891006 Cosmetic changes to prologue. (WRB) +C 891006 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 920501 Reformatted the REFERENCES section. (WRB) +C***END PROLOGUE DPSIFN + INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ, + * FN + INTEGER I1MACH + DOUBLE PRECISION ANS, ARG, B, DEN, ELIM, EPS, FLN, + * FX, RLN, RXSQ, R1M4, R1M5, S, SLOPE, T, TA, TK, TOL, TOLS, TRM, + * TRMR, TSS, TST, TT, T1, T2, WDTOL, X, XDMLN, XDMY, XINC, XLN, + * XM, XMIN, XQ, YINT + DOUBLE PRECISION D1MACH + DIMENSION B(22), TRM(22), TRMR(100), ANS(*) + SAVE NMAX, B + DATA NMAX /100/ +C----------------------------------------------------------------------- +C BERNOULLI NUMBERS +C----------------------------------------------------------------------- + DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), + * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19), + * B(20), B(21), B(22) /1.00000000000000000D+00, + * -5.00000000000000000D-01,1.66666666666666667D-01, + * -3.33333333333333333D-02,2.38095238095238095D-02, + * -3.33333333333333333D-02,7.57575757575757576D-02, + * -2.53113553113553114D-01,1.16666666666666667D+00, + * -7.09215686274509804D+00,5.49711779448621554D+01, + * -5.29124242424242424D+02,6.19212318840579710D+03, + * -8.65802531135531136D+04,1.42551716666666667D+06, + * -2.72982310678160920D+07,6.01580873900642368D+08, + * -1.51163157670921569D+10,4.29614643061166667D+11, + * -1.37116552050883328D+13,4.88332318973593167D+14, + * -1.92965793419400681D+16/ +C +C***FIRST EXECUTABLE STATEMENT DPSIFN + IERR = 0 + NZ=0 + IF (X.LE.0.0D0) IERR=1 + IF (N.LT.0) IERR=1 + IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 + IF (M.LT.1) IERR=1 + IF (IERR.NE.0) RETURN + MM=M + NX = MIN(-I1MACH(15),I1MACH(16)) + R1M5 = D1MACH(5) + R1M4 = D1MACH(4)*0.5D0 + WDTOL = MAX(R1M4,0.5D-18) +C----------------------------------------------------------------------- +C ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT +C----------------------------------------------------------------------- + ELIM = 2.302D0*(NX*R1M5-3.0D0) + XLN = LOG(X) + 41 CONTINUE + NN = N + MM - 1 + FN = NN + T = (FN+1)*XLN +C----------------------------------------------------------------------- +C OVERFLOW AND UNDERFLOW TEST FOR SMALL AND LARGE X +C----------------------------------------------------------------------- + IF (ABS(T).GT.ELIM) GO TO 290 + IF (X.LT.WDTOL) GO TO 260 +C----------------------------------------------------------------------- +C COMPUTE XMIN AND THE NUMBER OF TERMS OF THE SERIES, FLN+1 +C----------------------------------------------------------------------- + RLN = R1M5*I1MACH(14) + RLN = MIN(RLN,18.06D0) + FLN = MAX(RLN,3.0D0) - 3.0D0 + YINT = 3.50D0 + 0.40D0*FLN + SLOPE = 0.21D0 + FLN*(0.0006038D0*FLN+0.008677D0) + XM = YINT + SLOPE*FN + MX = INT(XM) + 1 + XMIN = MX + IF (N.EQ.0) GO TO 50 + XM = -2.302D0*RLN - MIN(0.0D0,XLN) + ARG = XM/N + ARG = MIN(0.0D0,ARG) + EPS = EXP(ARG) + XM = 1.0D0 - EPS + IF (ABS(ARG).LT.1.0D-3) XM = -ARG + FLN = X*XM/EPS + XM = XMIN - X + IF (XM.GT.7.0D0 .AND. FLN.LT.15.0D0) GO TO 200 + 50 CONTINUE + XDMY = X + XDMLN = XLN + XINC = 0.0D0 + IF (X.GE.XMIN) GO TO 60 + NX = INT(X) + XINC = XMIN - NX + XDMY = X + XINC + XDMLN = LOG(XDMY) + 60 CONTINUE +C----------------------------------------------------------------------- +C GENERATE W(N+MM-1,X) BY THE ASYMPTOTIC EXPANSION +C----------------------------------------------------------------------- + T = FN*XDMLN + T1 = XDMLN + XDMLN + T2 = T + XDMLN + TK = MAX(ABS(T),ABS(T1),ABS(T2)) + IF (TK.GT.ELIM) GO TO 380 + TSS = EXP(-T) + TT = 0.5D0/XDMY + T1 = TT + TST = WDTOL*TT + IF (NN.NE.0) T1 = TT + 1.0D0/FN + RXSQ = 1.0D0/(XDMY*XDMY) + TA = 0.5D0*RXSQ + T = (FN+1)*TA + S = T*B(3) + IF (ABS(S).LT.TST) GO TO 80 + TK = 2.0D0 + DO 70 K=4,22 + T = T*((TK+FN+1)/(TK+1.0D0))*((TK+FN)/(TK+2.0D0))*RXSQ + TRM(K) = T*B(K) + IF (ABS(TRM(K)).LT.TST) GO TO 80 + S = S + TRM(K) + TK = TK + 2.0D0 + 70 CONTINUE + 80 CONTINUE + S = (S+T1)*TSS + IF (XINC.EQ.0.0D0) GO TO 100 +C----------------------------------------------------------------------- +C BACKWARD RECUR FROM XDMY TO X +C----------------------------------------------------------------------- + NX = INT(XINC) + NP = NN + 1 + IF (NX.GT.NMAX) GO TO 390 + IF (NN.EQ.0) GO TO 160 + XM = XINC - 1.0D0 + FX = X + XM +C----------------------------------------------------------------------- +C THIS LOOP SHOULD NOT BE CHANGED. FX IS ACCURATE WHEN X IS SMALL +C----------------------------------------------------------------------- + DO 90 I=1,NX + TRMR(I) = FX**(-NP) + S = S + TRMR(I) + XM = XM - 1.0D0 + FX = X + XM + 90 CONTINUE + 100 CONTINUE + ANS(MM) = S + IF (FN.EQ.0) GO TO 180 +C----------------------------------------------------------------------- +C GENERATE LOWER DERIVATIVES, J.LT.N+MM-1 +C----------------------------------------------------------------------- + IF (MM.EQ.1) RETURN + DO 150 J=2,MM + FN = FN - 1 + TSS = TSS*XDMY + T1 = TT + IF (FN.NE.0) T1 = TT + 1.0D0/FN + T = (FN+1)*TA + S = T*B(3) + IF (ABS(S).LT.TST) GO TO 120 + TK = 4 + FN + DO 110 K=4,22 + TRM(K) = TRM(K)*(FN+1)/TK + IF (ABS(TRM(K)).LT.TST) GO TO 120 + S = S + TRM(K) + TK = TK + 2.0D0 + 110 CONTINUE + 120 CONTINUE + S = (S+T1)*TSS + IF (XINC.EQ.0.0D0) GO TO 140 + IF (FN.EQ.0) GO TO 160 + XM = XINC - 1.0D0 + FX = X + XM + DO 130 I=1,NX + TRMR(I) = TRMR(I)*FX + S = S + TRMR(I) + XM = XM - 1.0D0 + FX = X + XM + 130 CONTINUE + 140 CONTINUE + MX = MM - J + 1 + ANS(MX) = S + IF (FN.EQ.0) GO TO 180 + 150 CONTINUE + RETURN +C----------------------------------------------------------------------- +C RECURSION FOR N = 0 +C----------------------------------------------------------------------- + 160 CONTINUE + DO 170 I=1,NX + S = S + 1.0D0/(X+NX-I) + 170 CONTINUE + 180 CONTINUE + IF (KODE.EQ.2) GO TO 190 + ANS(1) = S - XDMLN + RETURN + 190 CONTINUE + IF (XDMY.EQ.X) RETURN + XQ = XDMY/X + ANS(1) = S - LOG(XQ) + RETURN +C----------------------------------------------------------------------- +C COMPUTE BY SERIES (X+K)**(-(N+1)) , K=0,1,2,... +C----------------------------------------------------------------------- + 200 CONTINUE + NN = INT(FLN) + 1 + NP = N + 1 + T1 = (N+1)*XLN + T = EXP(-T1) + S = T + DEN = X + DO 210 I=1,NN + DEN = DEN + 1.0D0 + TRM(I) = DEN**(-NP) + S = S + TRM(I) + 210 CONTINUE + ANS(1) = S + IF (N.NE.0) GO TO 220 + IF (KODE.EQ.2) ANS(1) = S + XLN + 220 CONTINUE + IF (MM.EQ.1) RETURN +C----------------------------------------------------------------------- +C GENERATE HIGHER DERIVATIVES, J.GT.N +C----------------------------------------------------------------------- + TOL = WDTOL/5.0D0 + DO 250 J=2,MM + T = T/X + S = T + TOLS = T*TOL + DEN = X + DO 230 I=1,NN + DEN = DEN + 1.0D0 + TRM(I) = TRM(I)/DEN + S = S + TRM(I) + IF (TRM(I).LT.TOLS) GO TO 240 + 230 CONTINUE + 240 CONTINUE + ANS(J) = S + 250 CONTINUE + RETURN +C----------------------------------------------------------------------- +C SMALL X.LT.UNIT ROUND OFF +C----------------------------------------------------------------------- + 260 CONTINUE + ANS(1) = X**(-N-1) + IF (MM.EQ.1) GO TO 280 + K = 1 + DO 270 I=2,MM + ANS(K+1) = ANS(K)/X + K = K + 1 + 270 CONTINUE + 280 CONTINUE + IF (N.NE.0) RETURN + IF (KODE.EQ.2) ANS(1) = ANS(1) + XLN + RETURN + 290 CONTINUE + IF (T.GT.0.0D0) GO TO 380 + NZ=0 + IERR=2 + RETURN + 380 CONTINUE + NZ=NZ+1 + ANS(MM)=0.0D0 + MM=MM-1 + IF (MM.EQ.0) RETURN + GO TO 41 + 390 CONTINUE + NZ=0 + IERR=3 + RETURN + END
--- a/liboctave/cruft/slatec-fn/module.mk Sun May 10 21:35:59 2015 -0700 +++ b/liboctave/cruft/slatec-fn/module.mk Sun May 10 21:37:19 2015 -0700 @@ -34,6 +34,7 @@ cruft/slatec-fn/dlnrel.f \ cruft/slatec-fn/dpchim.f \ cruft/slatec-fn/dpchst.f \ + cruft/slatec-fn/dpsifn.f \ cruft/slatec-fn/erf.f \ cruft/slatec-fn/gami.f \ cruft/slatec-fn/gamit.f \ @@ -44,6 +45,7 @@ cruft/slatec-fn/inits.f \ cruft/slatec-fn/pchim.f \ cruft/slatec-fn/pchst.f \ + cruft/slatec-fn/psifn.f \ cruft/slatec-fn/r9lgmc.f \ cruft/slatec-fn/r9lgit.f \ cruft/slatec-fn/r9gmit.f \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/cruft/slatec-fn/psifn.f Sun May 10 21:37:19 2015 -0700 @@ -0,0 +1,368 @@ +*DECK PSIFN + SUBROUTINE PSIFN (X, N, KODE, M, ANS, NZ, IERR) +C***BEGIN PROLOGUE PSIFN +C***PURPOSE Compute derivatives of the Psi function. +C***LIBRARY SLATEC +C***CATEGORY C7C +C***TYPE SINGLE PRECISION (PSIFN-S, DPSIFN-D) +C***KEYWORDS DERIVATIVES OF THE GAMMA FUNCTION, POLYGAMMA FUNCTION, +C PSI FUNCTION +C***AUTHOR Amos, D. E., (SNLA) +C***DESCRIPTION +C +C The following definitions are used in PSIFN: +C +C Definition 1 +C PSI(X) = d/dx (ln(GAMMA(X)), the first derivative of +C the LOG GAMMA function. +C Definition 2 +C K K +C PSI(K,X) = d /dx (PSI(X)), the K-th derivative of PSI(X). +C ___________________________________________________________________ +C PSIFN computes a sequence of SCALED derivatives of +C the PSI function; i.e. for fixed X and M it computes +C the M-member sequence +C +C ((-1)**(K+1)/GAMMA(K+1))*PSI(K,X) +C for K = N,...,N+M-1 +C +C where PSI(K,X) is as defined above. For KODE=1, PSIFN returns +C the scaled derivatives as described. KODE=2 is operative only +C when K=0 and in that case PSIFN returns -PSI(X) + LN(X). That +C is, the logarithmic behavior for large X is removed when KODE=1 +C and K=0. When sums or differences of PSI functions are computed +C the logarithmic terms can be combined analytically and computed +C separately to help retain significant digits. +C +C Note that CALL PSIFN(X,0,1,1,ANS) results in +C ANS = -PSI(X) +C +C Input +C X - Argument, X .gt. 0.0E0 +C N - First member of the sequence, 0 .le. N .le. 100 +C N=0 gives ANS(1) = -PSI(X) for KODE=1 +C -PSI(X)+LN(X) for KODE=2 +C KODE - Selection parameter +C KODE=1 returns scaled derivatives of the PSI +C function. +C KODE=2 returns scaled derivatives of the PSI +C function EXCEPT when N=0. In this case, +C ANS(1) = -PSI(X) + LN(X) is returned. +C M - Number of members of the sequence, M .ge. 1 +C +C Output +C ANS - A vector of length at least M whose first M +C components contain the sequence of derivatives +C scaled according to KODE. +C NZ - Underflow flag +C NZ.eq.0, A normal return +C NZ.ne.0, Underflow, last NZ components of ANS are +C set to zero, ANS(M-K+1)=0.0, K=1,...,NZ +C IERR - Error flag +C IERR=0, A normal return, computation completed +C IERR=1, Input error, no computation +C IERR=2, Overflow, X too small or N+M-1 too +C large or both +C IERR=3, Error, N too large. Dimensioned +C array TRMR(NMAX) is not large enough for N +C +C The nominal computational accuracy is the maximum of unit +C roundoff (=R1MACH(4)) and 1.0E-18 since critical constants +C are given to only 18 digits. +C +C DPSIFN is the Double Precision version of PSIFN. +C +C *Long Description: +C +C The basic method of evaluation is the asymptotic expansion +C for large X.ge.XMIN followed by backward recursion on a two +C term recursion relation +C +C W(X+1) + X**(-N-1) = W(X). +C +C This is supplemented by a series +C +C SUM( (X+K)**(-N-1) , K=0,1,2,... ) +C +C which converges rapidly for large N. Both XMIN and the +C number of terms of the series are calculated from the unit +C roundoff of the machine environment. +C +C***REFERENCES Handbook of Mathematical Functions, National Bureau +C of Standards Applied Mathematics Series 55, edited +C by M. Abramowitz and I. A. Stegun, equations 6.3.5, +C 6.3.18, 6.4.6, 6.4.9 and 6.4.10, pp.258-260, 1964. +C D. E. Amos, A portable Fortran subroutine for +C derivatives of the Psi function, Algorithm 610, ACM +C Transactions on Mathematical Software 9, 4 (1983), +C pp. 494-502. +C***ROUTINES CALLED I1MACH, R1MACH +C***REVISION HISTORY (YYMMDD) +C 820601 DATE WRITTEN +C 890531 Changed all specific intrinsics to generic. (WRB) +C 890531 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 920501 Reformatted the REFERENCES section. (WRB) +C***END PROLOGUE PSIFN + INTEGER I, IERR, J, K, KODE, M, MM, MX, N, NMAX, NN, NP, NX, NZ + INTEGER I1MACH + REAL ANS, ARG, B, DEN, ELIM, EPS, FLN, FN, FNP, FNS, FX, RLN, + * RXSQ, R1M4, R1M5, S, SLOPE, T, TA, TK, TOL, TOLS, TRM, TRMR, + * TSS, TST, TT, T1, T2, WDTOL, X, XDMLN, XDMY, XINC, XLN, XM, + * XMIN, XQ, YINT + REAL R1MACH + DIMENSION B(22), TRM(22), TRMR(100), ANS(*) + SAVE NMAX, B + DATA NMAX /100/ +C----------------------------------------------------------------------- +C BERNOULLI NUMBERS +C----------------------------------------------------------------------- + DATA B(1), B(2), B(3), B(4), B(5), B(6), B(7), B(8), B(9), B(10), + * B(11), B(12), B(13), B(14), B(15), B(16), B(17), B(18), B(19), + * B(20), B(21), B(22) /1.00000000000000000E+00, + * -5.00000000000000000E-01,1.66666666666666667E-01, + * -3.33333333333333333E-02,2.38095238095238095E-02, + * -3.33333333333333333E-02,7.57575757575757576E-02, + * -2.53113553113553114E-01,1.16666666666666667E+00, + * -7.09215686274509804E+00,5.49711779448621554E+01, + * -5.29124242424242424E+02,6.19212318840579710E+03, + * -8.65802531135531136E+04,1.42551716666666667E+06, + * -2.72982310678160920E+07,6.01580873900642368E+08, + * -1.51163157670921569E+10,4.29614643061166667E+11, + * -1.37116552050883328E+13,4.88332318973593167E+14, + * -1.92965793419400681E+16/ +C +C***FIRST EXECUTABLE STATEMENT PSIFN + IERR = 0 + NZ=0 + IF (X.LE.0.0E0) IERR=1 + IF (N.LT.0) IERR=1 + IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1 + IF (M.LT.1) IERR=1 + IF (IERR.NE.0) RETURN + MM=M + NX = MIN(-I1MACH(12),I1MACH(13)) + R1M5 = R1MACH(5) + R1M4 = R1MACH(4)*0.5E0 + WDTOL = MAX(R1M4,0.5E-18) +C----------------------------------------------------------------------- +C ELIM = APPROXIMATE EXPONENTIAL OVER AND UNDERFLOW LIMIT +C----------------------------------------------------------------------- + ELIM = 2.302E0*(NX*R1M5-3.0E0) + XLN = LOG(X) + 41 CONTINUE + NN = N + MM - 1 + FN = NN + FNP = FN + 1.0E0 + T = FNP*XLN +C----------------------------------------------------------------------- +C OVERFLOW AND UNDERFLOW TEST FOR SMALL AND LARGE X +C----------------------------------------------------------------------- + IF (ABS(T).GT.ELIM) GO TO 290 + IF (X.LT.WDTOL) GO TO 260 +C----------------------------------------------------------------------- +C COMPUTE XMIN AND THE NUMBER OF TERMS OF THE SERIES, FLN+1 +C----------------------------------------------------------------------- + RLN = R1M5*I1MACH(11) + RLN = MIN(RLN,18.06E0) + FLN = MAX(RLN,3.0E0) - 3.0E0 + YINT = 3.50E0 + 0.40E0*FLN + SLOPE = 0.21E0 + FLN*(0.0006038E0*FLN+0.008677E0) + XM = YINT + SLOPE*FN + MX = INT(XM) + 1 + XMIN = MX + IF (N.EQ.0) GO TO 50 + XM = -2.302E0*RLN - MIN(0.0E0,XLN) + FNS = N + ARG = XM/FNS + ARG = MIN(0.0E0,ARG) + EPS = EXP(ARG) + XM = 1.0E0 - EPS + IF (ABS(ARG).LT.1.0E-3) XM = -ARG + FLN = X*XM/EPS + XM = XMIN - X + IF (XM.GT.7.0E0 .AND. FLN.LT.15.0E0) GO TO 200 + 50 CONTINUE + XDMY = X + XDMLN = XLN + XINC = 0.0E0 + IF (X.GE.XMIN) GO TO 60 + NX = INT(X) + XINC = XMIN - NX + XDMY = X + XINC + XDMLN = LOG(XDMY) + 60 CONTINUE +C----------------------------------------------------------------------- +C GENERATE W(N+MM-1,X) BY THE ASYMPTOTIC EXPANSION +C----------------------------------------------------------------------- + T = FN*XDMLN + T1 = XDMLN + XDMLN + T2 = T + XDMLN + TK = MAX(ABS(T),ABS(T1),ABS(T2)) + IF (TK.GT.ELIM) GO TO 380 + TSS = EXP(-T) + TT = 0.5E0/XDMY + T1 = TT + TST = WDTOL*TT + IF (NN.NE.0) T1 = TT + 1.0E0/FN + RXSQ = 1.0E0/(XDMY*XDMY) + TA = 0.5E0*RXSQ + T = FNP*TA + S = T*B(3) + IF (ABS(S).LT.TST) GO TO 80 + TK = 2.0E0 + DO 70 K=4,22 + T = T*((TK+FN+1.0E0)/(TK+1.0E0))*((TK+FN)/(TK+2.0E0))*RXSQ + TRM(K) = T*B(K) + IF (ABS(TRM(K)).LT.TST) GO TO 80 + S = S + TRM(K) + TK = TK + 2.0E0 + 70 CONTINUE + 80 CONTINUE + S = (S+T1)*TSS + IF (XINC.EQ.0.0E0) GO TO 100 +C----------------------------------------------------------------------- +C BACKWARD RECUR FROM XDMY TO X +C----------------------------------------------------------------------- + NX = INT(XINC) + NP = NN + 1 + IF (NX.GT.NMAX) GO TO 390 + IF (NN.EQ.0) GO TO 160 + XM = XINC - 1.0E0 + FX = X + XM +C----------------------------------------------------------------------- +C THIS LOOP SHOULD NOT BE CHANGED. FX IS ACCURATE WHEN X IS SMALL +C----------------------------------------------------------------------- + DO 90 I=1,NX + TRMR(I) = FX**(-NP) + S = S + TRMR(I) + XM = XM - 1.0E0 + FX = X + XM + 90 CONTINUE + 100 CONTINUE + ANS(MM) = S + IF (FN.EQ.0.0E0) GO TO 180 +C----------------------------------------------------------------------- +C GENERATE LOWER DERIVATIVES, J.LT.N+MM-1 +C----------------------------------------------------------------------- + IF (MM.EQ.1) RETURN + DO 150 J=2,MM + FNP = FN + FN = FN - 1.0E0 + TSS = TSS*XDMY + T1 = TT + IF (FN.NE.0.0E0) T1 = TT + 1.0E0/FN + T = FNP*TA + S = T*B(3) + IF (ABS(S).LT.TST) GO TO 120 + TK = 3.0E0 + FNP + DO 110 K=4,22 + TRM(K) = TRM(K)*FNP/TK + IF (ABS(TRM(K)).LT.TST) GO TO 120 + S = S + TRM(K) + TK = TK + 2.0E0 + 110 CONTINUE + 120 CONTINUE + S = (S+T1)*TSS + IF (XINC.EQ.0.0E0) GO TO 140 + IF (FN.EQ.0.0E0) GO TO 160 + XM = XINC - 1.0E0 + FX = X + XM + DO 130 I=1,NX + TRMR(I) = TRMR(I)*FX + S = S + TRMR(I) + XM = XM - 1.0E0 + FX = X + XM + 130 CONTINUE + 140 CONTINUE + MX = MM - J + 1 + ANS(MX) = S + IF (FN.EQ.0.0E0) GO TO 180 + 150 CONTINUE + RETURN +C----------------------------------------------------------------------- +C RECURSION FOR N = 0 +C----------------------------------------------------------------------- + 160 CONTINUE + DO 170 I=1,NX + S = S + 1.0E0/(X+NX-I) + 170 CONTINUE + 180 CONTINUE + IF (KODE.EQ.2) GO TO 190 + ANS(1) = S - XDMLN + RETURN + 190 CONTINUE + IF (XDMY.EQ.X) RETURN + XQ = XDMY/X + ANS(1) = S - LOG(XQ) + RETURN +C----------------------------------------------------------------------- +C COMPUTE BY SERIES (X+K)**(-(N+1)) , K=0,1,2,... +C----------------------------------------------------------------------- + 200 CONTINUE + NN = INT(FLN) + 1 + NP = N + 1 + T1 = (FNS+1.0E0)*XLN + T = EXP(-T1) + S = T + DEN = X + DO 210 I=1,NN + DEN = DEN + 1.0E0 + TRM(I) = DEN**(-NP) + S = S + TRM(I) + 210 CONTINUE + ANS(1) = S + IF (N.NE.0) GO TO 220 + IF (KODE.EQ.2) ANS(1) = S + XLN + 220 CONTINUE + IF (MM.EQ.1) RETURN +C----------------------------------------------------------------------- +C GENERATE HIGHER DERIVATIVES, J.GT.N +C----------------------------------------------------------------------- + TOL = WDTOL/5.0E0 + DO 250 J=2,MM + T = T/X + S = T + TOLS = T*TOL + DEN = X + DO 230 I=1,NN + DEN = DEN + 1.0E0 + TRM(I) = TRM(I)/DEN + S = S + TRM(I) + IF (TRM(I).LT.TOLS) GO TO 240 + 230 CONTINUE + 240 CONTINUE + ANS(J) = S + 250 CONTINUE + RETURN +C----------------------------------------------------------------------- +C SMALL X.LT.UNIT ROUND OFF +C----------------------------------------------------------------------- + 260 CONTINUE + ANS(1) = X**(-N-1) + IF (MM.EQ.1) GO TO 280 + K = 1 + DO 270 I=2,MM + ANS(K+1) = ANS(K)/X + K = K + 1 + 270 CONTINUE + 280 CONTINUE + IF (N.NE.0) RETURN + IF (KODE.EQ.2) ANS(1) = ANS(1) + XLN + RETURN + 290 CONTINUE + IF (T.GT.0.0E0) GO TO 380 + NZ=0 + IERR=2 + RETURN + 380 CONTINUE + NZ=NZ+1 + ANS(MM)=0.0E0 + MM=MM-1 + IF(MM.EQ.0) RETURN + GO TO 41 + 390 CONTINUE + IERR=3 + NZ=0 + RETURN + END
--- a/liboctave/numeric/lo-specfun.cc Sun May 10 21:35:59 2015 -0700 +++ b/liboctave/numeric/lo-specfun.cc Sun May 10 21:37:19 2015 -0700 @@ -1,8 +1,10 @@ /* Copyright (C) 1996-2015 John W. Eaton +Copyright (C) 2007-2010 D. Martin Copyright (C) 2010 Jaroslav Hajek Copyright (C) 2010 VZLU Prague +Copyright (C) 2015 Carnë Draug This file is part of Octave. @@ -45,6 +47,7 @@ #include "lo-specfun.h" #include "mx-inlines.cc" #include "lo-mappers.h" +#include "lo-math.h" #include "Faddeeva.hh" @@ -183,6 +186,16 @@ F77_RET_T F77_FUNC (algams, ALGAMS) (const float&, float&, float&); + + F77_RET_T + F77_FUNC (psifn, PSIFN) (const float*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + float*, octave_idx_type*, octave_idx_type*); + + F77_RET_T + F77_FUNC (dpsifn, DPSIFN) (const double*, const octave_idx_type&, + const octave_idx_type&, const octave_idx_type&, + double*, octave_idx_type*, octave_idx_type*); } #if !defined (HAVE_ACOSH) @@ -3724,3 +3737,185 @@ dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd); } } + +static const double pi = 3.14159265358979323846; + +template<class T> +static T +Lanczos_approximation_psi (const T zc) +{ + // Coefficients for C.Lanczos expansion of psi function from XLiFE++ gammaFunctions + // psi_coef[k] = - (2k+1) * lg_coef[k] (see melina++ gamma functions) + // -1/12, 3/360,-5/1260, 7/1680,-9/1188, 11*691/360360,-13/156, 15*3617/122400, ? , ? + static const T dg_coeff[10] = { + -0.83333333333333333e-1, 0.83333333333333333e-2, + -0.39682539682539683e-2, 0.41666666666666667e-2, + -0.75757575757575758e-2, 0.21092796092796093e-1, + -0.83333333333333333e-1, 0.4432598039215686, + -0.3053954330270122e+1, 0.125318899521531e+2 + }; + + T overz2 = T (1.0) / (zc * zc); + T overz2k = overz2; + + T p = 0; + for (octave_idx_type k = 0; k < 10; k++, overz2k *= overz2) + p += dg_coeff[k] * overz2k; + p += log (zc) - T (0.5) / zc; + return p; +} + +template<class T> +T +psi (const T& z) +{ + static const double euler_mascheroni = 0.577215664901532860606512090082402431042; + + const bool is_int = (xfloor (z) == z); + + T p = 0; + if (z <= 0) + { + // limits - zeros of the gamma function + if (is_int) + p = -octave_Inf; // Matlab returns -Inf for psi (0) + else + // Abramowitz and Stegun, page 259, eq 6.3.7 + p = psi (1 - z) - (pi / tan (pi * z)); + } + else if (is_int) + { + // Abramowitz and Stegun, page 258, eq 6.3.2 + p = - euler_mascheroni; + for (octave_idx_type k = z - 1; k > 0; k--) + p += 1.0 / k; + } + else if (xfloor (z + 0.5) == z + 0.5) + { + // Abramowitz and Stegun, page 258, eq 6.3.3 and 6.3.4 + for (octave_idx_type k = z; k > 0; k--) + p += 1.0 / (2 * k - 1); + + p = - euler_mascheroni - 2 * log (2) + 2 * (p); + } + else + { + // adapted from XLiFE++ gammaFunctions + + T zc = z; + // Use formula for derivative of LogGamma(z) + if (z < 10) + { + const signed char n = 10 - z; + for (signed char k = n - 1; k >= 0; k--) + p -= 1.0 / (k + z); + zc += n; + } + p += Lanczos_approximation_psi (zc); + } + + return p; +} + +// explicit instantiations +template double psi<double> (const double& z); +template float psi<float> (const float& z); + +template<class T> +std::complex<T> +psi (const std::complex<T>& z) +{ + // adapted from XLiFE++ gammaFunctions + + typedef typename std::complex<T>::value_type P; + + P z_r = z.real (); + P z_ra = z_r; + + std::complex<T> dgam (0.0, 0.0); + if (z.imag () == 0) + dgam = std::complex<T> (psi (z_r), 0.0); + else if (z_r < 0) + dgam = psi (P (1.0) - z)- (P (pi) / tan (P (pi) * z)); + else + { + // Use formula for derivative of LogGamma(z) + std::complex<T> z_m = z; + if (z_ra < 8) + { + unsigned char n = 8 - z_ra; + z_m = z + std::complex<T> (n, 0.0); + + // Recurrence formula + // for | Re(z) | < 8 , use recursively DiGamma(z) = DiGamma(z+1) - 1/z + std::complex<T> z_p = z + P (n - 1); + for (unsigned char k = n; k > 0; k--, z_p -= 1.0) + dgam -= P (1.0) / z_p; + } + + // for | Re(z) | > 8, use derivative of C.Lanczos expansion for LogGamma + // psi(z) = log(z) - 1/(2z) - 1/12z^2 + 3/360z^4 - 5/1260z^6 + 7/1680z^8 - 9/1188z^10 + ... + // (Abramowitz&Stegun, page 259, formula 6.3.18 + dgam += Lanczos_approximation_psi (z_m); + } + return dgam; +} + +// explicit instantiations +template Complex psi<double> (const Complex& z); +template FloatComplex psi<float> (const FloatComplex& z); + + +template<typename T> +static inline void +fortran_psifn (const T z, const octave_idx_type n, T* ans, + octave_idx_type* ierr); + +template<> +inline void +fortran_psifn<double> (const double z, const octave_idx_type n, + double* ans, octave_idx_type* ierr) +{ + octave_idx_type flag = 0; + F77_XFCN (dpsifn, DPSIFN, (&z, n, 1, 1, ans, &flag, ierr)); +} + +template<> +inline void +fortran_psifn<float> (const float z, const octave_idx_type n, + float* ans, octave_idx_type* ierr) +{ + octave_idx_type flag = 0; + F77_XFCN (psifn, PSIFN, (&z, n, 1, 1, ans, &flag, ierr)); +} + +template<class T> +T +psi (const octave_idx_type n, const T z) +{ + T ans; + octave_idx_type ierr = 0; + fortran_psifn<T> (z, n, &ans, &ierr); + if (ierr == 0) + { + // Remember that psifn and dpsifn return scales values + // When n is 1: do nothing since ((-1)**(n+1)/gamma(n+1)) == 1 + // When n is 0: change sign since ((-1)**(n+1)/gamma(n+1)) == -1 + if (n > 1) + // FIXME xgamma here is a killer for our precision since it grows + // way too fast + ans = ans / (pow (-1.0, n + 1) / xgamma (double (n+1))); + else if (n == 0) + ans = -ans; + } + else if (ierr == 2) + ans = - octave_Inf; + else // we probably never get here + ans = octave_NaN; + + return ans; +} + +// explicit instantiations +template double psi<double> (const octave_idx_type n, const double z); +template float psi<float> (const octave_idx_type n, const float z);
--- a/liboctave/numeric/lo-specfun.h Sun May 10 21:35:59 2015 -0700 +++ b/liboctave/numeric/lo-specfun.h Sun May 10 21:37:19 2015 -0700 @@ -663,4 +663,24 @@ ellipj (const Complex& u, double m, Complex& sn, Complex& cn, Complex& dn, double& err); +//! Digamma function. +//! +//! Only defined for double and float. +template<class T> +extern OCTAVE_API T psi (const T& z); + +//! Digamma function for complex input. +//! +//! Only defined for double and float. +template<class T> +extern OCTAVE_API std::complex<T> psi (const std::complex<T>& z); + +//! Polygamma function. +//! +//! Only defined for double and float. +//! @param n must be non-negative. If zero, the digamma function is computed. +//! @param z must be real and non-negative. +template<class T> +extern OCTAVE_API T psi (const octave_idx_type n, const T z); + #endif
--- a/m4/acinclude.m4 Sun May 10 21:35:59 2015 -0700 +++ b/m4/acinclude.m4 Sun May 10 21:37:19 2015 -0700 @@ -291,7 +291,7 @@ [AC_LANG_PUSH(C++) ac_octave_save_CPPFLAGS="$CPPFLAGS" CPPFLAGS="$QT_CPPFLAGS $CPPFLAGS" - AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ + AC_PREPROC_IFELSE([AC_LANG_PROGRAM([[ #include <Qsci/qsciglobal.h> ]], [[ #if QSCINTILLA_VERSION < 0x020600 @@ -417,7 +417,7 @@ [AC_LANG_PUSH(C++) ac_octave_save_CPPFLAGS="$CPPFLAGS" CPPFLAGS="$QT_CPPFLAGS $CPPFLAGS" - AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[ + AC_PREPROC_IFELSE([AC_LANG_PROGRAM([[ #include <Qt/qglobal.h> ]], [[ #if QT_VERSION < 0x040700
--- a/scripts/image/imformats.m Sun May 10 21:35:59 2015 -0700 +++ b/scripts/image/imformats.m Sun May 10 21:37:19 2015 -0700 @@ -79,7 +79,7 @@ persistent formats = default_formats (); if (nargin == 0 && nargout == 0) - error ("imformats: pretty print not yet implemented."); + pretty_print_formats (formats); elseif (nargin >= 1) if (isstruct (arg1)) arrayfun (@is_valid_format, arg1); @@ -281,6 +281,40 @@ end_try_catch endfunction +function pretty_print_formats (formats) + ## define header names (none should be shorter than 3 characters) + headers = {"Extension", "isa", "Info", "Read", "Write", "Alpha", "Description"}; + cols_length = cellfun (@numel, headers); + + ## Adjust the maximal length of the extensions column + extensions = cellfun (@strjoin, {formats.ext}, {", "}, + "UniformOutput", false); + cols_length(1) = max (max (cellfun (@numel, extensions)), cols_length(1)); + headers{1} = postpad (headers{1}, cols_length(1), " "); + + ## Print the headers + disp (strjoin (headers, " | ")); + under_headers = cellfun (@(x) repmat ("-", 1, numel (x)), headers, + "UniformOutput", false); + disp (strjoin (under_headers, "-+-")); + + template = strjoin (arrayfun (@(x) sprintf ("%%-%is", x), cols_length, + "UniformOutput", false), " | "); + + ## Print the function handle for this things won't be a pretty table. So + ## instead we replace them with "yes" or "no", based on the support it has. + yes_no_cols = cat (2, {formats.isa}(:), {formats.info}(:), {formats.read}(:), + {formats.write}(:), {formats.alpha}(:)); + empty = cellfun (@isempty, yes_no_cols); + yes_no_cols(empty) = "no"; + yes_no_cols(! empty) = "yes"; + + descriptions = {formats.description}; + table = cat (2, extensions(:), yes_no_cols, descriptions(:)); + printf ([template "\n"], table'{:}); + +endfunction + ## When imread or imfinfo are called, the file must exist or the ## function defined by imformats will never be called. Because ## of this, we must create a file for the tests to work.
--- a/scripts/strings/cstrcat.m Sun May 10 21:35:59 2015 -0700 +++ b/scripts/strings/cstrcat.m Sun May 10 21:37:19 2015 -0700 @@ -45,14 +45,16 @@ function st = cstrcat (varargin) - if (nargin < 1) - print_usage (); - elseif (! iscellstr (varargin)) + if (nargin == 0) + ## Special because if varargin is empty, iscellstr still returns + ## true but then "[varargin{:}]" would be of class double. + st = ""; + elseif (iscellstr (varargin)) + st = [varargin{:}]; + else error ("cstrcat: expecting arguments to character strings"); endif - st = [varargin{:}]; - endfunction @@ -65,7 +67,8 @@ %!assert (cstrcat ("foo", "bar"), "foobar") %!assert (cstrcat (["a"; "bb"], ["foo"; "bar"]), ["a foo"; "bbbar"]) +%!assert (cstrcat (), "") + ## Test input validation -%!error cstrcat () %!error cstrcat (1, 2)
--- a/scripts/strings/strcat.m Sun May 10 21:35:59 2015 -0700 +++ b/scripts/strings/strcat.m Sun May 10 21:37:19 2015 -0700 @@ -82,10 +82,8 @@ function st = strcat (varargin) if (nargin == 0) - print_usage (); - endif - - if (nargin == 1) + st = ""; + elseif (nargin == 1) st = varargin{1}; else ## Convert to cells of strings @@ -149,5 +147,5 @@ %!assert (strcat (1, 2), strcat (char (1), char (2))) %!assert (strcat ("", 2), strcat ([], char (2))) -%!error strcat () +%!assert (strcat (), "")