wehub-resource-sync
108 lines
4.2 KiB
C++
108 lines
4.2 KiB
C++
/*
|
|
* ******************************************************************************
|
|
* *
|
|
* *
|
|
* * This program and the accompanying materials are made available under the
|
|
* * terms of the Apache License, Version 2.0 which is available at
|
|
* * https://www.apache.org/licenses/LICENSE-2.0.
|
|
* *
|
|
* * See the NOTICE file distributed with this work for additional
|
|
* * information regarding copyright ownership.
|
|
* * Unless required by applicable law or agreed to in writing, software
|
|
* * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
|
|
* * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
|
|
* * License for the specific language governing permissions and limitations
|
|
* * under the License.
|
|
* *
|
|
* * SPDX-License-Identifier: Apache-2.0
|
|
* *****************************************************************************
|
|
*/
|
|
|
|
//
|
|
// @author Yurii Shyrma (iuriish@yahoo.com)
|
|
//
|
|
|
|
#ifndef LIBND4J_GAMMAMATHFUNC_H
|
|
#define LIBND4J_GAMMAMATHFUNC_H
|
|
#include <ops/declarable/helpers/helpers.h>
|
|
|
|
#include "array/NDArray.h"
|
|
|
|
namespace sd {
|
|
namespace ops {
|
|
namespace helpers {
|
|
|
|
// calculate the digamma function for each element for array
|
|
SD_LIB_HIDDEN void diGamma(LaunchContext* context, NDArray& x, NDArray& z);
|
|
|
|
// calculate the polygamma function
|
|
SD_LIB_HIDDEN void polyGamma(LaunchContext* context, NDArray& n, NDArray& x, NDArray& z);
|
|
|
|
// calculate the digamma function for one element
|
|
// implementation is based on serial representation written in terms of the Hurwitz zeta function as polygamma =
|
|
// (-1)^{n+1} * n! * zeta(n+1, x)
|
|
template <typename T>
|
|
SD_HOST_DEVICE T diGammaScalar(T x) {
|
|
const int xInt = static_cast<int>(x);
|
|
|
|
// negative and zero
|
|
if (x <= 0) {
|
|
if (x == xInt) // integer
|
|
return DataTypeUtils::infOrMax<T>();
|
|
else
|
|
return diGammaScalar<T>(1 - x) -
|
|
M_PI / math::sd_tan<T, T>(M_PI * x); // use reflection formula psi(1-x) = psi(x) + pi*cot(pi*x)
|
|
}
|
|
|
|
// positive integer
|
|
if (x == xInt &&
|
|
xInt <= 20) { // psi(n) = -Euler_Mascheroni_const + sum_from_k=1_to_n-1( 1/k ), for n = 1,2,3,...inf, we use this
|
|
// formula only for n <= 20 to avoid time consuming sum calculation for bigger n
|
|
T result = static_cast<T>(-0.577215664901532);
|
|
for (LongType i = 1; i <= xInt - 1; ++i) {
|
|
result += static_cast<T>(1) / i;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
// positive half-integer
|
|
if (x - xInt == 0.5 && xInt <= 20) { // psi(n+0.5) = -Euler_Mascheroni_const - 2*ln(2) + sum_from_k=1_to_n( 2/(2*k-1)
|
|
// ) , for n = 1,2,3,...inf, we use this formula only for n <= 20 to avoid
|
|
// time consuming sum calculation for bigger n
|
|
T result = static_cast<T>(-0.577215664901532 - static_cast<T>(2) * math::sd_log<T, T>(static_cast<T>(2)));
|
|
for (LongType i = 1; i <= xInt; ++i) {
|
|
result += static_cast<T>(2) / (2 * i - 1);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
// positive, smaller then 5; we should use number > 5 in order to have satisfactory accuracy in asymptotic expansion
|
|
if (x < 5) return diGammaScalar<T>(1 + x) - static_cast<T>(1) / x; // recurrence formula psi(x) = psi(x+1) - 1/x.
|
|
|
|
// *** other positive **** //
|
|
|
|
// truncated expansion formula (from wiki)
|
|
// psi(x) = log(x) - 1/(2*x) - 1/(12*x^2) + 1/(120*x^4) - 1/(252*x^6) + 1/(240*x^8) - 5/(660*x^10) + 691/(32760*x^12)
|
|
// - 1/(12*x^14) + ...
|
|
|
|
if (x >= (sizeof(T) > 4 ? 1.e16 : 1.e8)) // if x is too big take into account only log(x)
|
|
return math::sd_log<T, T>(x);
|
|
|
|
// coefficients used in truncated asymptotic expansion formula
|
|
const T coeffs[7] = {-(T)1 / 12, (T)1 / 120, -(T)1 / 252, (T)1 / 240, -(T)5 / 660, (T)691 / 32760, -(T)1 / 12};
|
|
// const T coeffs[7] = {-0.0833333333333333, 0.00833333333333333, -0.00396825396825397, 0.00416666666666667,
|
|
// -0.00757575757575758, 0.0210927960927961, -0.0833333333333333};
|
|
|
|
const T x2Inv = static_cast<T>(1) / (x * x);
|
|
T result = static_cast<T>(0);
|
|
|
|
for (int i = 6; i >= 0; --i) result = (result + coeffs[i]) * x2Inv;
|
|
return result + math::sd_log<T, T>(x) - static_cast<T>(0.5) / x;
|
|
}
|
|
|
|
} // namespace helpers
|
|
} // namespace ops
|
|
} // namespace sd
|
|
|
|
#endif // LIBND4J_GAMMAMATHFUNC_H
|