Mercurial > octave-antonio
annotate libinterp/corefcn/psi.cc @ 20154:45565ecec019
New function psi to compute the digamma function.
* libinterp/corefcn/psi.cc: file for the new function file (implementation
is actually in lo-specfun.cc).
* liboctave/numeric/lo-specfun.cc, liboctave/numeric/lo-specfun.h: added
implementation of the digamma (psi )function. Partly based on diGamma()
from XLiFE++ 1.1 (file gammaFunctions.cpp which was previously melina++)
which is under GPL 3 or later and why D. Martin is also added to copyright.
* doc/interpreter/arith.txi: add function entry to the manual.
* libinterp/corefcn/module.mk: add file to the build system.
* NEWS: note new function.
author | Carnë Draug <carandraug@octave.org> |
---|---|
date | Sun, 15 Mar 2015 03:31:16 +0000 |
parents | |
children | 1fae49e34a1a |
rev | line source |
---|---|
20154
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
1 /* |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
2 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
3 Copyright (C) 2015 Carnë Draug |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
4 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
5 This file is part of Octave. |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
6 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
7 Octave is free software; you can redistribute it and/or modify it |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
8 under the terms of the GNU General Public License as published by the |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
9 Free Software Foundation; either version 3 of the License, or (at your |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
10 option) any later version. |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
11 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
15 for more details. |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
16 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
17 You should have received a copy of the GNU General Public License |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
18 along with Octave; see the file COPYING. If not, see |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
19 <http://www.gnu.org/licenses/>. |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
20 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
21 */ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
22 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
23 #ifdef HAVE_CONFIG_H |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
24 #include <config.h> |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
25 #endif |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
26 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
27 #include "ov.h" |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
28 #include "defun.h" |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
29 #include "error.h" |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
30 #include "dNDArray.h" |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
31 #include "fNDArray.h" |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
32 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
33 #include "lo-specfun.h" |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
34 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
35 DEFUN (psi, args, , |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
36 "-*- texinfo -*-\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
37 @deftypefn {Function File} {} psi (@var{z})\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
38 Compute the psi (digamma) function.\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
39 \n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
40 @tex\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
41 $$\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
42 \\Psi (z) = {d (log (\\Gamma (z))) \\over dx}\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
43 $$\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
44 @end tex\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
45 @ifnottex\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
46 @example\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
47 @group\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
48 psi (z) = d (log (gamma (z))) / dx\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
49 @end group\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
50 @end example\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
51 @end ifnottex\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
52 \n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
53 @seealso{gamma, gammainc, gammaln}\n\ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
54 @end deftypefn") |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
55 { |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
56 octave_value retval; |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
57 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
58 const octave_idx_type nargin = args.length (); |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
59 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
60 if (nargin != 1) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
61 { |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
62 print_usage (); |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
63 return retval; |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
64 } |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
65 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
66 #define FLOAT_BRANCH(T, A, M, E) \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
67 if (args(0).is_ ## T ##_type ()) \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
68 { \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
69 const A ## NDArray z = args(0).M ## array_value (); \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
70 A ## NDArray psi_z (z.dims ()); \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
71 \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
72 const E* zv = z.data (); \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
73 E* psi_zv = psi_z.fortran_vec (); \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
74 const octave_idx_type n = z.numel (); \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
75 for (octave_idx_type i = 0; i < n; i++) \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
76 psi_zv[i] = psi<E> (zv[i]); \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
77 \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
78 retval = psi_z; \ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
79 } |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
80 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
81 FLOAT_BRANCH(double, , , double) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
82 else FLOAT_BRANCH(single, Float, float_, float) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
83 else |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
84 { |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
85 error ("psi: Z must be a floating point"); |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
86 } |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
87 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
88 #undef FLOAT_BRANCH |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
89 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
90 return retval; |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
91 } |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
92 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
93 /* |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
94 %!shared em |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
95 %! em = 0.577215664901532860606512090082402431042; # Euler-Mascheroni Constant |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
96 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
97 %!assert (psi (ones (7, 3, 5)), repmat (-em, [7 3 5])) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
98 %!assert (psi ([0 1]), [-Inf -em]) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
99 %!assert (psi ([-20:1]), [repmat(-Inf, [1 21]) -em]) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
100 %!assert (psi (single ([0 1])), single ([-Inf -em])) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
101 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
102 ## Abramowitz and Stegun, page 258, eq 6.3.5 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
103 %!test |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
104 %! z = [-10:.1:-.1 .1:.1:20]; # drop the 0 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
105 %! assert (psi (z + 1), psi (z) + 1 ./ z, eps*1000) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
106 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
107 ## Abramowitz and Stegun, page 258, eq 6.3.2 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
108 %!assert (psi (1), -em) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
109 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
110 ## Abramowitz and Stegun, page 258, eq 6.3.3 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
111 %!assert (psi (1/2), -em - 2 * log (2)) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
112 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
113 ## The following tests are from Pascal Sebah and Xavier Gourdon (2002) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
114 ## "Introduction to the Gamma Function" |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
115 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
116 ## Interesting identities of the digamma function, in section of 5.1.3 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
117 %!assert (psi (1/3), - em - (3/2) * log(3) - ((sqrt (3) / 6) * pi), eps*10) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
118 %!assert (psi (1/4), - em -3 * log (2) - pi /2, eps*10) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
119 %!assert (psi (1/6), - em -2 * log (2) - (3/2) * log (3) - ((sqrt (3) / 2) * pi), eps*10) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
120 |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
121 ## First 6 zeros of the digamma function, in section of 5.1.5 (and also on |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
122 ## Abramowitz and Stegun, page 258, eq 6.3.19) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
123 %!assert (psi ( 1.46163214496836234126265954232572132846819620400644), 0, eps) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
124 %!assert (psi (-0.504083008264455409258269304533302498955385182368579), 0, eps*10) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
125 %!assert (psi (-1.573498473162390458778286043690434612655040859116846), 0, eps) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
126 %!assert (psi (-2.610720868444144650001537715718724207951074010873480), 0, eps*10) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
127 %!assert (psi (-3.635293366436901097839181566946017713948423861193530), 0, eps) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
128 %!assert (psi (-4.653237761743142441714598151148207363719069416133868), 0, eps*100) |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
129 */ |
45565ecec019
New function psi to compute the digamma function.
Carnë Draug <carandraug@octave.org>
parents:
diff
changeset
|
130 |