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