[Scilab-users] Riemann Zeta update

[Stéphane Mottelet]

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 at 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20220525/9e7c9d1c/attachment.htm>