1/* 2 * zeta2 - Hurwitz Zeta function 3 * 4 * Calc is open software; you can redistribute it and/or modify it under 5 * the terms of the version 2.1 of the GNU Lesser General Public License 6 * as published by the Free Software Foundation. 7 * 8 * Calc is distributed in the hope that it will be useful, but WITHOUT 9 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 10 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 11 * Public License for more details. 12 * 13 * A copy of version 2.1 of the GNU Lesser General Public License is 14 * distributed with calc under the filename COPYING-LGPL. You should have 15 * received a copy with calc; if not, write to Free Software Foundation, Inc. 16 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 17 * 18 * Under source code control: 2013/08/11 01:31:28 19 * File existed as early as: 2013 20 */ 21 22 23/* 24 * hide internal function from resource debugging 25 */ 26static resource_debug_level; 27resource_debug_level = config("resource_debug", 0); 28 29 30define hurwitzzeta(s,a){ 31 local realpart_a imagpart_s tmp tmp1 tmp2 tmp3; 32 local sum1 sum2 sum3 i k n precision result limit; 33 local limit_function offset offset_squared rest_sum eps; 34 /* 35 According to Linas Vepstas' "An efficient algorithm for accelerating 36 the convergence of oscillatory series, useful for computing the 37 poly-logarithm and Hurwitz zeta functions" the Euler-Maclaurin series 38 is the fastest in most cases. 39 40 With a lot of help of the PARI/GP implementation by Prof. Henri Cohen, 41 hence the different license. 42 */ 43 eps=epsilon( epsilon() * 1e-3); 44 realpart_a=re(a); 45 if(realpart_a>1.5){ 46 tmp=floor(realpart_a-0.5); 47 sum1 = 0; 48 for( i = 1 ; i <= tmp ; i++){ 49 sum1 += ( a - i )^( -s ); 50 } 51 epsilon(eps); 52 return (hurwitzzeta(s,a-tmp)-sum1); 53 } 54 if(realpart_a<=0){ 55 tmp=ceil(-realpart_a+0.5); 56 for( i = 0 ; i <= tmp-1 ; i++){ 57 sum2 += ( a + i )^( -s ); 58 } 59 epsilon(eps); 60 return (hurwitzzeta(s,a+tmp)+sum2); 61 } 62 precision=digits(1/epsilon()); 63 realpart_a=re(s); 64 imagpart_s=im(s); 65 epsilon(1e-9); 66 result=s-1.; 67 if(abs(result)<0.1){ 68 result=-1; 69 } 70 else 71 result=ln(result); 72 limit=(precision*ln(10)-re((s-.5)*result)+(1.*realpart_a)*ln(2.*pi()))/2; 73 limit=max(2,ceil(max(limit,abs(s*1.)/2))); 74 limit_function=ceil(sqrt((limit+realpart_a/2-.25)^2+(imagpart_s*1.)^2/4)/ 75 pi()); 76 if (config("user_debug") > 0) { 77 print "limit_function = " limit_function; 78 print "limit = " limit; 79 print "prec = " precision; 80 } 81 /* Full precision plus 5 digits angstzuschlag*/ 82 epsilon( (10^(-precision)) * 1e-5); 83 tmp3=(a+limit_function+0.)^(-s); 84 sum3 = tmp3/2; 85 for(n=0;n<=limit_function-1;n++){ 86 sum3 += (a+n)^(-s); 87 } 88 result=sum3; 89 offset=a+limit_function; 90 offset_squared=1./(offset*offset); 91 tmp1=2*s-1; 92 tmp2=s*(s-1); 93 rest_sum=bernoulli(2*limit); 94 for(k=2*limit-2;k>=2;k-=2){ 95 rest_sum=bernoulli(k)+offset_squared* 96 (k*k+tmp1*k+tmp2)*rest_sum/((k+1)*(k+2)); 97 } 98 rest_sum=offset*(1+offset_squared*tmp2*rest_sum/2); 99 result+=rest_sum*tmp3/(s-1); 100 epsilon(eps); 101 return result; 102} 103 104 105/* 106 * restore internal function from resource debugging 107 * report important interface functions 108 */ 109config("resource_debug", resource_debug_level),; 110if (config("resource_debug") & 3) { 111 print "hurwitzzeta(s,a)"; 112} 113