110 lines
No EOL
3.2 KiB
MQL5
110 lines
No EOL
3.2 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| zeta.mqh |
|
|
//| Copyright 2025, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
#property copyright "Copyright 2025, MetaQuotes Ltd."
|
|
#property link "https://www.mql5.com"
|
|
#include "error.mqh"
|
|
#include "const.mqh"
|
|
namespace special {
|
|
namespace cephes {
|
|
|
|
namespace detail {
|
|
/* Expansion coefficients
|
|
* for Euler-Maclaurin summation formula
|
|
* (2k)! / B2k
|
|
* where B2k are Bernoulli numbers
|
|
*/
|
|
double zeta_A[] = {
|
|
12.0,
|
|
-720.0,
|
|
30240.0,
|
|
-1209600.0,
|
|
47900160.0,
|
|
-1.8924375803183791606e9, /*1.307674368e12/691 */
|
|
7.47242496e10,
|
|
-2.950130727918164224e12, /*1.067062284288e16/3617 */
|
|
1.1646782814350067249e14, /*5.109094217170944e18/43867 */
|
|
-4.5979787224074726105e15, /*8.028576626982912e20/174611 */
|
|
1.8152105401943546773e17, /*1.5511210043330985984e23/854513 */
|
|
-7.1661652561756670113e18 /*1.6938241367317436694528e27/236364091 */
|
|
};
|
|
|
|
/* 30 Nov 86 -- error in third coefficient fixed */
|
|
} // namespace detail
|
|
|
|
double inline zeta(double x, double q) {
|
|
int i;
|
|
double a, b, k, s, t, w;
|
|
|
|
if (x == 1.0)
|
|
return (infinity());
|
|
|
|
if (x < 1.0) {
|
|
set_error("zeta", SF_ERROR_DOMAIN, NULL);
|
|
return (quiet_NaN());
|
|
}
|
|
|
|
if (q <= 0.0) {
|
|
if (q == floor(q)) {
|
|
set_error("zeta", SF_ERROR_SINGULAR, NULL);
|
|
return (infinity());
|
|
}
|
|
if (x != floor(x))
|
|
{
|
|
set_error("zeta", SF_ERROR_DOMAIN, NULL);
|
|
return (quiet_NaN());
|
|
}
|
|
}
|
|
|
|
/* Asymptotic expansion
|
|
* https://dlmf.nist.gov/25.11#E43
|
|
*/
|
|
if (q > 1e8) {
|
|
return (1 / (x - 1) + 1 / (2 * q)) * pow(q, 1 - x);
|
|
}
|
|
|
|
/* Euler-Maclaurin summation formula */
|
|
|
|
/* Permit negative q but continue sum until n+q > +9 .
|
|
* This case should be handled by a reflection formula.
|
|
* If q<0 and x is an integer, there is a relation to
|
|
* the polyGamma function.
|
|
*/
|
|
s = pow(q, -x);
|
|
a = q;
|
|
i = 0;
|
|
b = 0.0;
|
|
while ((i < 9) || (a <= 9.0)) {
|
|
i += 1;
|
|
a += 1.0;
|
|
b = pow(a, -x);
|
|
s += b;
|
|
if (abs(b / s) < detail::MACHEP)
|
|
return (s);
|
|
}
|
|
|
|
w = a;
|
|
s += b * w / (x - 1.0);
|
|
s -= 0.5 * b;
|
|
a = 1.0;
|
|
k = 0.0;
|
|
for (i = 0; i < 12; i++) {
|
|
a *= x + k;
|
|
b /= w;
|
|
t = a * b / detail::zeta_A[i];
|
|
s = s + t;
|
|
t = abs(t / s);
|
|
if (t < detail::MACHEP)
|
|
return (s);
|
|
k += 1.0;
|
|
a *= x + k;
|
|
b /= w;
|
|
k += 1.0;
|
|
}
|
|
return (s);
|
|
}
|
|
|
|
} // namespace cephes
|
|
} // namespace special |