Skip to content

Commit 7da6f59

Browse files
chore: fix comments & follow stdlib's conventions
--- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: passed - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed ---
1 parent 5031684 commit 7da6f59

3 files changed

Lines changed: 198 additions & 76 deletions

File tree

lib/node_modules/@stdlib/math/base/special/digammaf/lib/main.js

Lines changed: 91 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -14,27 +14,79 @@
1414
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
1515
* See the License for the specific language governing permissions and
1616
* limitations under the License.
17+
*
18+
*
19+
* ## Notice
20+
*
21+
* The original C code, long comment, copyright, license, and constants are from [Cephes]{@link http://www.netlib.org/cephes}. The implementation follows the original, but has been modified for JavaScript.
22+
*
23+
* ```text
24+
* Cephes Math Library Release 2.2: June, 1992
25+
* Copyright 1984, 1987, 1992 by Stephen L. Moshier
26+
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
27+
* ```
1728
*/
1829

1930
'use strict';
2031

2132
// MODULES //
2233

23-
var float64ToFloat32 = require( '@stdlib/number/float64/base/to-float32' );
34+
var f32 = require( '@stdlib/number/float64/base/to-float32' );
2435
var floorf = require( '@stdlib/math/base/special/floorf' );
2536
var lnf = require( '@stdlib/math/base/special/lnf' );
26-
var cotf = require( '@stdlib/math/base/special/cotf' );
37+
var tanf = require( '@stdlib/math/base/special/tanf' );
2738
var FLOAT32_PI = require( '@stdlib/constants/float32/pi' );
2839
var EULERGAMMA = require( '@stdlib/constants/float32/eulergamma' );
2940
var isnanf = require( '@stdlib/math/base/assert/is-nanf' );
3041
var polyvalP = require( './polyval_p.js' );
3142

3243

44+
// VARIABLES //
45+
46+
// Threshold above which the asymptotic expansion's correction term is below float32 epsilon:
47+
var ASYMPTOTIC_THRESHOLD = 1.0e8;
48+
49+
// Minimum value of x for which we use the asymptotic expansion directly:
50+
var MIN_SAFE_ASYMPTOTIC = 10.0;
51+
52+
3353
// MAIN //
3454

3555
/**
3656
* Evaluates the digamma function for a single-precision floating-point number.
3757
*
58+
* ## Method
59+
*
60+
* 1. For \\\\(x \u003c= 0\\\\), we use the reflection formula
61+
*
62+
* ```tex
63+
* \\psi(1-x) = \\psi(x) + \\frac{\\pi}{\\tan(\\pi \\cdot x)}
64+
* ```
65+
*
66+
* to make \\\\(x\\\\) positive.
67+
*
68+
* 2. For positive integer \\\\(x \u003c= 10\\\\), we use the recurrence
69+
*
70+
* ```tex
71+
* \\psi(n) = -\\gamma + \\sum_{k=1}^{n-1} \\frac{1}{k}
72+
* ```
73+
*
74+
* where \\\\(\\gamma\\\\) is the Euler-Mascheroni constant.
75+
*
76+
* 3. For general \\\\(x\\\\), we apply the recurrence relation
77+
*
78+
* ```tex
79+
* \\psi(x+1) = \\psi(x) + \\frac{1}{x}
80+
* ```
81+
*
82+
* until \\\\(x \u003e= 10\\\\), then use the asymptotic expansion
83+
*
84+
* ```tex
85+
* \\psi(x) = \\ln(x) - \\frac{1}{2x} - \\sum_{k=1}^{\\infty} \\frac{B_{2k}}{2k x^{2k}}
86+
* ```
87+
*
88+
* where \\\\(B_{2k}\\\\) are Bernoulli numbers.
89+
*
3890
* @param {number} x - input value
3991
* @returns {number} function value
4092
*
@@ -64,19 +116,17 @@ var polyvalP = require( './polyval_p.js' );
64116
*/
65117
function digammaf(x) {
66118
var negative;
67-
var nz;
119+
var tmp;
120+
var rem;
68121
var xx;
69-
var p;
70-
var q;
71122
var s;
72123
var w;
73124
var y;
74125
var z;
75126
var i;
76127
var n;
77128

78-
xx = float64ToFloat32(x);
79-
nz = 0.0;
129+
xx = f32(x);
80130
negative = false;
81131

82132
if (isnanf(xx)) {
@@ -86,58 +136,57 @@ function digammaf(x) {
86136
// Handle negative arguments using reflection formula...
87137
if (xx <= 0.0) {
88138
negative = true;
89-
q = xx;
90-
p = floorf(q);
139+
140+
// Argument reduction for tan:
141+
rem = f32(xx - floorf(xx));
142+
143+
// Reflect:
144+
xx = f32(1.0 - xx);
91145

92146
// Check for singularity (negative integer or zero)...
93-
if (p === q) {
147+
if (rem === 0.0) {
94148
return NaN;
95149
}
96-
97-
// Compute fractional part...
98-
nz = float64ToFloat32(q - p);
99-
if (nz === 0.5) {
100-
nz = 0.0;
101-
} else {
102-
if (nz > 0.5) {
103-
p += 1.0;
104-
nz = float64ToFloat32(q - p);
105-
}
106-
nz = float64ToFloat32(FLOAT32_PI * cotf(float64ToFloat32(FLOAT32_PI * nz))); // eslint-disable-line max-len
150+
if (rem > 0.5) {
151+
rem = f32(rem - 1.0);
107152
}
108-
xx = float64ToFloat32(1.0 - xx);
153+
tmp = f32(FLOAT32_PI / tanf(f32(FLOAT32_PI * rem)));
154+
} else {
155+
tmp = 0.0;
109156
}
110157

111158
// Use direct formula for small positive integers...
112-
if (xx <= 10.0 && xx === floorf(xx)) {
159+
if (xx <= MIN_SAFE_ASYMPTOTIC && xx === floorf(xx)) {
113160
y = 0.0;
114161
n = xx;
115162
for (i = 1; i < n; i++) {
116-
w = i;
117-
y = float64ToFloat32(y + float64ToFloat32(1.0 / w));
163+
y = f32(y + f32(1.0 / i));
118164
}
119-
y = float64ToFloat32(y - EULERGAMMA);
120-
} else {
121-
// Use recurrence to make x large enough for asymptotic expansion...
122-
s = xx;
123-
w = 0.0;
124-
while (s < 10.0) {
125-
w = float64ToFloat32(w + float64ToFloat32(1.0 / s));
126-
s = float64ToFloat32(s + 1.0);
165+
y = f32(y - EULERGAMMA);
166+
if (negative) {
167+
y = f32(y - tmp);
127168
}
169+
return y;
170+
}
128171

129-
// Asymptotic expansion...
130-
if (s < 1.0e8) {
131-
z = float64ToFloat32(1.0 / float64ToFloat32(s * s));
132-
y = float64ToFloat32(z * polyvalP(z));
133-
} else {
134-
y = 0.0;
135-
}
136-
y = float64ToFloat32(lnf(s) - float64ToFloat32(0.5 / s) - y - w);
172+
// Use recurrence to make x large enough for asymptotic expansion...
173+
s = xx;
174+
w = 0.0;
175+
while (s < MIN_SAFE_ASYMPTOTIC) {
176+
w = f32(w + f32(1.0 / s));
177+
s = f32(s + 1.0);
137178
}
138179

180+
// Asymptotic expansion...
181+
if (s < ASYMPTOTIC_THRESHOLD) {
182+
z = f32(1.0 / f32(s * s));
183+
y = f32(z * polyvalP(z));
184+
} else {
185+
y = 0.0;
186+
}
187+
y = f32(lnf(s) - f32(0.5 / s) - y - w);
139188
if (negative) {
140-
y = float64ToFloat32(y - nz);
189+
y = f32(y - tmp);
141190
}
142191
return y;
143192
}
Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
#/
2+
# @license Apache-2.0
3+
#
4+
# Copyright (c) 2026 The Stdlib Authors.
5+
#
6+
# Licensed under the Apache License, Version 2.0 (the "License");
7+
# you may not use this file except in compliance with the License.
8+
# You may obtain a copy of the License at
9+
#
10+
# http://www.apache.org/licenses/LICENSE-2.0
11+
#
12+
# Unless required by applicable law or agreed to in writing, software
13+
# distributed under the License is distributed on an "AS IS" BASIS,
14+
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15+
# See the License for the specific language governing permissions and
16+
# limitations under the License.
17+
#/
18+
19+
# VARIABLES #
20+
21+
ifndef VERBOSE
22+
QUIET := @
23+
else
24+
QUIET :=
25+
endif
26+
27+
# Determine the OS ([1][1], [2][2]).
28+
#
29+
# [1]: https://en.wikipedia.org/wiki/Uname#Examples
30+
# [2]: http://stackoverflow.com/a/27776822/2225624
31+
OS ?= $(shell uname)
32+
ifneq (, $(findstring MINGW,$(OS)))
33+
OS := WINNT
34+
else
35+
ifneq (, $(findstring MSYS,$(OS)))
36+
OS := WINNT
37+
else
38+
ifneq (, $(findstring CYGWIN,$(OS)))
39+
OS := WINNT
40+
else
41+
ifneq (, $(findstring Windows_NT,$(OS)))
42+
OS := WINNT
43+
endif
44+
endif
45+
endif
46+
endif
47+
48+
49+
# RULES #
50+
51+
#/
52+
# Removes generated files for building an add-on.
53+
#
54+
# @example
55+
# make clean-addon
56+
#/
57+
clean-addon:
58+
$(QUIET) -rm -f *.o *.node
59+
60+
.PHONY: clean-addon
61+
62+
#/
63+
# Removes generated files.
64+
#
65+
# @example
66+
# make clean
67+
#/
68+
clean: clean-addon
69+
70+
.PHONY: clean

0 commit comments

Comments
 (0)