33 lines
1,016 B
Text
33 lines
1,016 B
Text
Bern: procedure options (main); /* 4 July 2014 */
|
|
declare i fixed binary;
|
|
declare B complex fixed (31);
|
|
|
|
Bernoulli: procedure (n) returns (complex fixed (31));
|
|
declare n fixed binary;
|
|
declare anum(0:n) fixed (31), aden(0:n) fixed (31);
|
|
declare (j, m) fixed;
|
|
declare F fixed (31);
|
|
|
|
do m = 0 to n;
|
|
anum(m) = 1;
|
|
aden(m) = m+1;
|
|
do j = m to 1 by -1;
|
|
anum(j-1) = j*( aden(j)*anum(j-1) - aden(j-1)*anum(j) );
|
|
aden(j-1) = ( aden(j-1) * aden(j) );
|
|
F = gcd(abs(anum(j-1)), abs(aden(j-1)) );
|
|
if F ^= 1 then
|
|
do;
|
|
anum(j-1) = anum(j-1) / F;
|
|
aden(j-1) = aden(j-1) / F;
|
|
end;
|
|
end;
|
|
end;
|
|
return ( complex(anum(0), aden(0)) );
|
|
end Bernoulli;
|
|
|
|
do i = 0, 1, 2 to 36 by 2; /* 36 is upper limit imposed by hardware. */
|
|
B = Bernoulli(i);
|
|
put skip edit ('B(' , trim(i) , ')=' , real(B) , '/' , trim(imag(B)) )
|
|
(3 A, column(10), F(32), 2 A);
|
|
end;
|
|
end Bern;
|