Hi,
For real argument, we could easily interface std::riemann_zeta :
https://en.cppreference.com/w/cpp/numeric/special_functions/riemann_zeta
If you have a compiler (under windows you can install the minGW atoms
module), you can run the following script:
code=[
"#define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1 "
"#include <cmath>"
"#include ""double.hxx"" "
"#include ""function.hxx"" "
"extern ""C"" "
"{ "
"#include ""Scierror.h"" "
"#include ""localization.h"" "
"} "
"/* ====================================================================
*/ "
"types::Function::ReturnValue sci_zeta(types::typed_list &in, int
_iRetCount, types::typed_list &out) "
"{ "
"types::Double* pDblOut; "
"types::Double* pDblIn; "
"if (in.size() != 1) "
"{ "
"Scierror(77, _(""%s: Wrong number of input argument(s): %d
expected.\n""), ""zeta"", 1); "
"return types::Function::Error; "
"} "
"if (_iRetCount >1) "
"{ "
"Scierror(78, _(""%s: Wrong number of output argument(s): %d
expected.""), ""zeta"", 1); "
"return types::Function::Error; "
"} "
"if (in[0]->isDouble() == false ||
in[0]->getAs<types::Double>()->isComplex()) "
"{ "
"Scierror(999, _(""%s: Wrong type for input argument #%d: real
expected.\n""), ""zeta"", 1); "
"return types::Function::Error; "
"} "
"pDblIn = in[0]->getAs<types::Double>(); "
"pDblOut = pDblIn->clone(); "
"for (int i = 0 ; i <pDblIn->getSize() ; ++i) "
"{ "
"pDblOut->set(i,std::riemann_zeta(pDblIn->get(i))); "
"} "
"out.push_back(pDblOut); "
"return types::Function::OK; "
"} "
];
ulink
files = fullfile(TMPDIR,"sci_zeta.cpp");
mputl(code,files)
ilib_verbose(1)
ilib_build("zeta", ["zeta" "sci_zeta" "cppsci"],files,[])
exec loader.sce
tic
disp(zeta(0.5))
disp(toc())
Le 22/05/2022 à 08:31, Lester Anderson a écrit :
Hi all,
After a lot of trial and error, I have managed to get a set of
functions to compute the approximations of Riemann's Zeta for negative
and positive real values; values of n > 1e6 seem to give better results:
function zs=zeta_s(z, n)
// Summation loop
zs=1;
if z == 0
zs = -0.5
elseif z == 1
zs = %inf
else
for i = 2: n-1
zs = zs + i.^-z;
end
end
endfunction
function zfn=zeta_functional_eqn(s)
// Riemann's functional equation
// Analytic continuation for negative values
zfn = 2.^s .* %pi.^(s - 1) .* sin(%pi.*s./2) .* gamma(1 - s)
.* zeta_s((1 - s),n)
endfunction
For even values of s < -20 the values of Zeta(s) increase in value and
are not as close to zero as expected e.g. zeta_functional_eqn(-40)
gives 7.5221382. At small even values e.g. -10, the result is of the
order of ~1e-18 (close enough to zero). Any ideas why the even zeta
values increase or how to reduce that response?
The solution over the critical strip (zero to one) is not so efficient
unless n is very large( > 1e8), and there seems to be a performance
issue when using a for-loop compared to vectorisation. Vectorised n
speeds things up quite a bit.
function zs2=zeta_0_1(s, n)
zs2=0
for i = 1: n
zs2 = zs2 + (-1).^(i + 1)./(i.^s );
end
zs2 = 1./(1 - 2.^( 1-s )).*zs2;
endfunction
function zs1=zeta_0_1(s, n)
// Vectorised version
zs1=0
k=linspace(1,n,n);
zs1 = sum((-1).^(k+ 1)./(k.^s ));
zs1 = 1./(1 - 2.^( 1-s )).*zs1;
endfunction
For example, calculating the approximation of Zeta(0.5) using a
for-loop takes ~150s to give a value of -1.4602337981325388 (quite
close), whereas the vectorised version does the computation in under
20s, both tested using n=1e8. Can the functions be optimised to
improve speed and accuracy?
Using Scilab 6.1.1 on Windows 10 (16 Gb RAM).
Thanks
Lester
_______________________________________________
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users
--
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet
_______________________________________________
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users
