mbernoa.m

function mbernoa
%This program is a direct conversion of the corresponding Fortran program in
%S. Zhang & J. Jin "Computation of Special Functions" (Wiley, 1996).
%online: http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
%
%Converted by f2matlab open source project:
%online: https://sourceforge.net/projects/f2matlab/
% written by Ben Barrowes (barrowes@alum.mit.edu)
%

%       ===========================================================
%       Purpose: This program computes Bernoulli number Bn using
%                subroutine BERNOA
%       Example: Compute Bernouli number Bn for n = 0,1,...,10
%                Computed results:
%                   n            Bn
%                 --------------------------
%                   0     .100000000000D+01
%                   1    -.500000000000D+00
%                   2     .166666666667D+00
%                   4    -.333333333333D-01
%                   6     .238095238095D-01
%                   8    -.333333333333D-01
%                  10     .757575757576D-01
%       ===========================================================
n=[];b=[];
  b=0;
 b=zeros(1,200+1);
fprintf(1,'%s \n','  please enter nmax');
%        READ(*,*)N
n=10;
[n,b]=bernoa(n,b);
fprintf(1,'%s \n','   n            bn');
fprintf(1,'%s \n',' --------------------------');
fprintf(1,[repmat(' ',1,2),'%3g','%22.12g' ' \n'],0,b(0+1));
fprintf(1,[repmat(' ',1,2),'%3g','%22.12g' ' \n'],1,b(1+1));
for  k=2:2:n;
fprintf(1,[repmat(' ',1,2),'%3g','%22.12g' ' \n'],k,b(k+1));
end;  k=n+1;
%format(2x,i3,d22.12);
end
function [n,bn]=bernoa(n,bn,varargin);
%       ======================================
%       Purpose: Compute Bernoulli number Bn
%       Input :  n --- Serial number
%       Output:  BN(n)--- Bn
%       ======================================
bn(0+1)=1.0d0;
bn(1+1)=-0.5d0;
for  m=2:n;
s=-(1.0d0./(m+1.0d0)-0.5d0);
for  k=2:m-1;
r=1.0d0;
for  j=2:k;
r=r.*(j+m-k)./j;
end;  j=k+1;
s=s-r.*bn(k+1);
end;  k=m-1+1;
bn(m+1)=s;
end;  m=fix(n)+1;
for  m=3:2:n;
bn(m+1)=0.0d0;
end;  m=fix(n)+1;
return;
end