samer@4: function [f] = deta(z,k) samer@4: % deta - Calculates Dirichlet functions of the form samer@4: % samer@4: % f = sum((-1)^n/(k*n+1)^z) samer@4: % samer@4: % over the entire complex plane samer@4: % Z may be complex and any size samer@4: % Best accuracy for Abs(z) < 100 samer@4: % samer@4: % Usage: f = deta(z) samer@4: % or f = deta(z,k) samer@4: % samer@4: % where k determines which Dirichlet function to sum samer@4: % For Eta (Zeta, Lambda): k=1 samer@4: % For Betad: k=2 samer@4: % samer@4: % This function can use a LOT of memory when size(z) samer@4: % is large. Consider using the Memory and Pack commands. samer@4: % Also, consider breaking z up into smaller chunks. samer@4: % samer@4: % Requires a complex Gamma routine. samer@4: % Tested under version 5.3.1 samer@4: % samer@4: %see also: Zeta, Eta, Lambda, Betad samer@4: %see also: sym/zeta.m samer@4: %see also: mhelp zeta samer@4: samer@4: %Andrew Odlyzko has Riemann Zeta critical line zeros listed on: samer@4: %http://www.research.att.com/~amo/zeta_tables/index.html samer@4: samer@4: %Paul Godfrey samer@4: %pgodfrey@conexant.com samer@4: %March 24, 2001 samer@4: samer@4: if nargin==1 samer@4: k=1; samer@4: end samer@4: k=k(1); samer@4: if k<1 | k>2 samer@4: warning('Unknown function being calculated! Results valid only for Real(z)>0.5') samer@4: % k=1 --> Eta --> Zeta or Lambda --> Bern numbers samer@4: % k=2 --> Betad --> Euler numbers samer@4: end samer@4: samer@4: [sizz]=size(z); samer@4: z=z(:).'; % make z a row vector samer@4: rz=real(z); samer@4: iz=imag(z); samer@4: iszero=find(z==0); samer@4: iseven=find(z==(round(z/2)*2) & rz<0 & iz==0); samer@4: isodd= find(z==(round((z-1)/2)*2+1) & rz<0 & iz==0); samer@4: samer@4: r=1/2; % reflection point samer@4: L=find(rz< r); samer@4: if ~isempty(L) samer@4: zL=z(L); samer@4: z(L)=1-zL; samer@4: end samer@4: samer@4: %series coefficients were calculated using samer@4: % c(m)=sum from n=m to 64 of (binomial(n,m)/2^n) for m=0:64 samer@4: samer@4: %coefficients are symmetrical about the 0.5 value. Each pair sums to +-1 samer@4: %abs(coefficients) look like erfc(k*m)/2 due to binomial terms samer@4: %sum(cm) must = 0.5 = eta(0) = betad(0) samer@4: %worst case error occurs for z = 0.5 + i*large samer@4: samer@4: cm= [ .99999999999999999997; samer@4: -.99999999999999999821; samer@4: .99999999999999994183; samer@4: -.99999999999999875788; samer@4: .99999999999998040668; samer@4: -.99999999999975652196; samer@4: .99999999999751767484; samer@4: -.99999999997864739190; samer@4: .99999999984183784058; samer@4: -.99999999897537734890; samer@4: .99999999412319859549; samer@4: -.99999996986230482845; samer@4: .99999986068828287678; samer@4: -.99999941559419338151; samer@4: .99999776238757525623; samer@4: -.99999214148507363026; samer@4: .99997457616475604912; samer@4: -.99992394671207596228; samer@4: .99978893483826239739; samer@4: -.99945495809777621055; samer@4: .99868681159465798081; samer@4: -.99704078337369034566; samer@4: .99374872693175507536; samer@4: -.98759401271422391785; samer@4: .97682326283354439220; samer@4: -.95915923302922997013; samer@4: .93198380256105393618; samer@4: -.89273040299591077603; samer@4: .83945793215750220154; samer@4: -.77148960729470505477; samer@4: .68992761745934847866; samer@4: -.59784149990330073143; samer@4: .50000000000000000000; samer@4: -.40215850009669926857; samer@4: .31007238254065152134; samer@4: -.22851039270529494523; samer@4: .16054206784249779846; samer@4: -.10726959700408922397; samer@4: .68016197438946063823e-1; samer@4: -.40840766970770029873e-1; samer@4: .23176737166455607805e-1; samer@4: -.12405987285776082154e-1; samer@4: .62512730682449246388e-2; samer@4: -.29592166263096543401e-2; samer@4: .13131884053420191908e-2; samer@4: -.54504190222378945440e-3; samer@4: .21106516173760261250e-3; samer@4: -.76053287924037718971e-4; samer@4: .25423835243950883896e-4; samer@4: -.78585149263697370338e-5; samer@4: .22376124247437700378e-5; samer@4: -.58440580661848562719e-6; samer@4: .13931171712321674741e-6; samer@4: -.30137695171547022183e-7; samer@4: .58768014045093054654e-8; samer@4: -.10246226511017621219e-8; samer@4: .15816215942184366772e-9; samer@4: -.21352608103961806529e-10; samer@4: .24823251635643084345e-11; samer@4: -.24347803504257137241e-12; samer@4: .19593322190397666205e-13; samer@4: -.12421162189080181548e-14; samer@4: .58167446553847312884e-16; samer@4: -.17889335846010823161e-17; samer@4: .27105054312137610850e-19]; samer@4: samer@4: cm=flipud(cm).'; % sum from small to big samer@4: nmax=size(cm,2); samer@4: n=(1:k:k*nmax)'; samer@4: n=flipud(n); samer@4: % z is a LR vector samer@4: % n is an UD vector samer@4: [Z,N]=meshgrid(z,n); samer@4: samer@4: % this can take a LOT of memory samer@4: f=cm*(N.^-Z); samer@4: % but it's really fast to form the series expansion N.^-Z samer@4: % and then sum it by an inner product cm*() :) samer@4: samer@4: %reflect across 1/2 samer@4: if ~isempty(L) samer@4: zz=z(L); samer@4: if k==1 samer@4: % Eta function reflection samer@4: % for test: deta(1,1) should = log(2) samer@4: t=(2-2.^(zz+1))./(2.^zz-2)./pi.^zz; samer@4: f(L)=t.*cos(pi/2*zz).*gamma(zz).*f(L); samer@4: if ~isempty(iszero) samer@4: f(iszero)=0.5; samer@4: end samer@4: if ~isempty(iseven) samer@4: f(iseven)=0; samer@4: end samer@4: end samer@4: if k==2 samer@4: % Betad function reflection samer@4: %for test: deta(0,2) should = 0.5 samer@4: %for test: deta(1,2) should = pi/4 samer@4: f(L)=(2/pi).^zz.*sin(pi/2*zz).*gamma(zz).*f(L); samer@4: if ~isempty(isodd) samer@4: f(isodd)=0; samer@4: end samer@4: end samer@4: if k>2 samer@4: % Insert reflection formula for other Dirichlet functions here samer@4: f(L)=(1/pi).^zz.*gamma(zz).*f(L); samer@4: f(L)=NaN; samer@4: end samer@4: end samer@4: samer@4: f = reshape(f,sizz); samer@4: return samer@4: