c_cpp/etc/calc/cal/zeta2.cal

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