octave_packages/specfun-1.1.0/zeta.m

   1 ## Copyright (C) 2006   Sylvain Pelissier   <sylvain.pelissier@gmail.com>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{z} =} zeta (@var{t})
  18 ## Compute the Riemann's Zeta function.
  19 ## @seealso{Si}
  20 ## @end deftypefn
  21
  22 function z = zeta(t)
  23   if (nargin != 1)
  24     print_usage;
  25   endif
  26   z = zeros(size(t));
  27   for j = 1:prod(size(t))
  28     if(real(t(j)) >= 0)
  29       if(imag(t(j)) == 0 && real(t(j)) > 1)
  30         F= @(x) 1./(gamma(t(j))).*x.^(t(j)-1)./(exp(x)-1);
  31         z(j) = quad(F,0,Inf);
  32       elseif(t(j) == 0)
  33         z(j) = -0.5;
  34       elseif(t(j) == 1)
  35         z(j) = Inf;
  36       else
  37         for k = 1:100
  38           z(j) += (-1).^(k-1)./(k.^t(j));
  39         endfor
  40         z(j) = 1./(1-2.^(1-t(j))).*z(j);
  41       endif
  42     else
  43       z(j) = 2.^t(j).*pi.^(t(j)-1).*sin(pi.*t(j)./2).*gamma(1-t(j)).*zeta(1-t(j));
  44     endif
  45   endfor
  46 endfunction