function [f] = deta(z,k)
%DETA Calculates Dirichlet functions of the form
%
% f = sum((-1)^n/(k*n+1)^z)
%
% over the entire complex plane
% Z may be complex and any size
% Best accuracy for Abs(z) < 100
%
% Usage: f = deta(z)
% or f = deta(z,k)
%
% where k determines which Dirichlet function to sum
% For Eta (Zeta, Lambda): k=1
% For Betad: k=2
%
% Intended for use with the author's Special Function math library.
%
% This function can use a LOT of memory when size(z)
% is large. Consider using the Memory and Pack commands.
% Also, consider breaking z up into smaller chunks.
%
% Requires a complex Gamma routine.
% Tested under version 5.3.1
%
% Version 2. Less memory required for a faster and more accurate computation.
%
%see also: Zeta, Eta, Lambda, Betad
%see also: sym/zeta.m
%see also: mhelp zeta
%Andrew Odlyzko has Riemann Zeta critical line zeros listed on:
%http://www.research.att.com/~amo/zeta_tables/index.html
%Paul Godfrey
%pgodfrey@conexant.com
%March 24, 2001, Version 1
%May 19, 2006, Version 2
if nargin==1
k=1;
end
k=k(1);
if k<1 | k>2
warning('Unknown function being calculated! Results valid only for Real(z)>0.5')
% k=1 --> Eta --> Zeta or Lambda --> Bern numbers
% k=2 --> Betad --> Euler numbers
end
[sizz]=size(z);
z=z(:).'; % make z a row vector
rz=real(z);
iz=imag(z);
iszero=find(z==0);
iseven=find(z==(round(z/2)*2) & rz<0 & iz==0);
isodd= find(z==(round((z-1)/2)*2+1) & rz<0 & iz==0);
r=1/2; % reflection point
L=find(rz< r);
if ~isempty(L)
zL=z(L);
z(L)=1-zL;
end
%series coefficients were calculated using
% c(m)=1/2*(-1)^m*sum from n=m to 64 of (binomial(n,m)/2^n) for m=0:64
%coefficients are symmetrical about the 0.5 value. Each pair sums to +-1
%abs(coefficients) look like erfc(k*m)/2 due to binomial terms
%sum(cm) must = 0.5 = eta(0) = betad(0)
%worst case error occurs for z = 0.5 + i*large
cm=[ 1.00000000000000000000;
-1.00000000000000000000;
1.00000000000000000000;
-1.00000000000000000000;
1.00000000000000000000;
-1.00000000000000000000;
1.00000000000000000000;
-1.00000000000000000000;
1.00000000000000000000;
-1.00000000000000000000;
1.00000000000000000000;
-1.00000000000000000000;
1.00000000000000000000;
-1.00000000000000000000;
1.00000000000000000000;
-.99999999999999999997;
.99999999999999999982;
-.99999999999999999878;
.99999999999999999230;
-.99999999999999995447;
.99999999999999974639;
-.99999999999999866635;
.99999999999999336436;
-.99999999999996869857;
.99999999999985975801;
-.99999999999940220763;
.99999999999757200611;
-.99999999999059012624;
.99999999996515613529;
-.99999999987657568406;
.99999999958130751329;
-.99999999863835432275;
.99999999575056017674;
-.99999998726219556575;
.99999996329504842881;
-.99999989824136334285;
.99999972837896339617;
-.99999930142752569235;
.99999826775562388311;
-.99999585585451966155;
.99999042907703516303;
-.99997864899908100771;
.99995396693098658703;
-.99990402879321415451;
.99980642243302258185;
-.99962205486377183349;
.99928538365035742345;
-.99869083661390218869;
.99767515209329116264;
-.99599616339595456855;
.99330978148021601802;
-.98914852321858179268;
.98290663582613045468;
-.97383823338804077494;
.96107529662332196641;
-.94367129194415995478;
.92067314290383872514;
-.89121937308026943102;
.85465607260963168661;
-.81065616865344050265;
.75932294737121745470;
-.70125750034706351522;
.63757281651412048482;
-.56984466069178424614;
.50000000000000000000;
-.43015533930821575386;
.36242718348587951518;
-.29874249965293648478;
.24067705262878254530;
-.18934383134655949735;
.14534392739036831339;
-.10878062691973056898;
.79326857096161274865e-1;
-.56328708055840045216e-1;
.38924703376678033590e-1;
-.26161766611959225064e-1;
.17093364173869545322e-1;
-.10851476781418207318e-1;
.66902185197839819816e-2;
-.40038366040454314481e-2;
.23248479067088373647e-2;
-.13091633860978113143e-2;
.71461634964257655302e-3;
-.37794513622816650749e-3;
.19357756697741814923e-3;
-.95971206785845488970e-4;
.46033069013412965117e-4;
-.21351000918992292409e-4;
.95709229648369713430e-5;
-.41441454803384526498e-5;
.17322443761168887862e-5;
-.69857247430764713040e-6;
.27162103660382992473e-6;
-.10175863665714996119e-6;
.36704951571187421956e-7;
-.12737804434253854871e-7;
.42494398232565498619e-8;
-.13616456772471574361e-8;
.41869248671347827662e-9;
-.12342431594030601457e-9;
.34843864708354335956e-10;
-.94098737605662302149e-11;
.24279938925459658938e-11;
-.59779237374453738245e-12;
.14024199404418025462e-12;
-.31301427448857128942e-13;
.66356387857651004875e-14;
-.13336468301471878290e-14;
.25361143178057599114e-15;
-.45531217416366371003e-16;
.76984511683282582508e-17;
-.12225722610064191311e-17;
.18180600804398070112e-18;
-.25230554058481114317e-19;
.32550517447267863440e-20;
-.38868187771535226060e-21;
.42740686869144698810e-22;
-.43027767751216363886e-23;
.39383676555996902370e-24;
-.32506176440823300892e-25;
.23952940142278239919e-26;
-.15564160233229265221e-27;
.87791490932414168166e-29;
-.42112329914907008714e-30;
.16702696209117826039e-31;
-.52587209151973559528e-33;
.12322119532494628802e-34;
-.19101783200862172004e-36;
.14693679385278593850e-38];
cm=double(cm(:));
cm=flipud(cm).'; % sum from small to big
nmax=size(cm,2);
n=(1:k:k*nmax)';
n=flipud(n);
if 0==1
% z is a LR vector
% n is an UD vector
[Z,N]=meshgrid(z,n);
% this can take a LOT of memory
f=cm*(N.^-Z);
% but it's really fast to form the series expansion N.^-Z
% and then sum it by an inner product cm*() :)
else
f=0;
for in=1:length(n)
f=f+cm(in)*n(in).^-z;
end
end
%reflect across 1/2
if ~isempty(L)
zz=z(L);
if k==1
% Eta function reflection
% for test: deta(1,1) should = log(2)
t=(2-2.^(zz+1))./(2.^zz-2)./pi.^zz;
f(L)=t.*cos(pi/2*zz).*gamma(zz).*f(L);
if ~isempty(iszero)
f(iszero)=0.5;
end
if ~isempty(iseven)
f(iseven)=0;
end
end
if k==2
% Betad function reflection
%for test: deta(0,2) should = 0.5
%for test: deta(1,2) should = pi/4
f(L)=(2/pi).^zz.*sin(pi/2*zz).*gamma(zz).*f(L);
if ~isempty(isodd)
f(isodd)=0;
end
end
if k>2
% Insert reflection formula for other Dirichlet functions here
f(L)=(1/pi).^zz.*gamma(zz).*f(L);
f(L)=NaN;
end
end
f = reshape(f,sizz);
return
%a demo of this routine is
x=0:1/16:1;
y=0:1/16:32;
[X,Y]=meshgrid(x,y);
z=X+i*Y;
clear X Y
f = deta(z,1);
p=find(abs(f)>5);
if ~isempty(p)
f(p)=NaN;
end
mesh(x,y,abs(f));
view([83 3]);
axis([0 1 0 32 0 6]);
rotate3d
return
