module sfbern; % Procedures for computing Bernoulli numbers.
%
% Author: Chris Cannam, Oct 1992.
% Module for Euler numbers added by Kerry Gaskell, Sep 1993
%
% Note there is currently no Bernoulli polynomial function.
% There was one in an older version but it won't convert directly.
% This is Something To Be Done.
fluid '(compute!-bernoulli);
imports complex!*on!*switch, complex!*off!*switch,
complex!*restore!*switch;
exports nearest!-int!-to!-bf, bernoulli!*calc, multi!*bern,
single!*bern, retrieve!*bern;
algebraic operator bernoulli;
symbolic operator bernoulli!*calc;
algebraic (bernoullirules := {
bernoulli(~n) => 1 when numberp n and n = 0,
bernoulli(~n) => -1/2 when numberp n and n = 1,
bernoulli(~n) => 0 when numberp n and impart n = 0
and n = floor n and n/2 neq floor (n/2) and n > 0,
bernoulli(~n) => bernoulli!*calc n when numberp n
and impart n = 0 and n = floor n and n > 0
})$
algebraic (let bernoullirules);
algebraic procedure bernoulli!*calc n;
begin scalar precom, result, prepre;
% Loading the SPECFAUX module will do two things. First it will
% set compute!-bernoulli to true, so that there is no future attempt
% to load it. Then it will set up a table of values in the variable
% bernoulli!-alist, where the table is computed at compile time rather
% than load or run-time. This will make compiling specfaux.red a fairly
% slow process. It also has bad consequences for any attempt to run
% this code interpreted.
% Note: ACN find the "algebraic symbolic" stuff here pretty heavy
% and confusing, but without it REDUCE sticks in calls to aeval (etc)
% in places where that is not wanted. Maybe a future version of the
% language will make mixed algebraic/symbolic mode code less delicate.
if (lisp null compute!-bernoulli) then
symbolic <<errorset!*('(load_package '(specfaux)), nil); nil>>;
precom := complex!*off!*switch();
if (prepre := precision(0)) < !!nfpd
then precision (!!nfpd + 1);
result := algebraic symbolic retrieve!*bern(n);
precision prepre;
complex!*restore!*switch(precom);
return result;
end;
symbolic procedure retrieve!*bern n;
begin scalar info, result;
integer heldpre;
info := assoc(n, bernoulli!-alist);
if not info then result := bern!*calc (n, '(() () ()))
else
<< info := cdr info;
if !*rounded then
if (heldpre := cadr info) and heldpre >= c!:prec!:() then
result := mk!*sq !*f2q rd!:prep caddr info
else if (result := car info) then
result := mk!*sq !*f2q mkround
divbf(i2bf!: caadr result,
i2bf!: cdadr result)
else result := bern!*calc(n, info)
else if not (result := car info) then
result := bern!*calc(n,info) >>;
return result;
end;
symbolic procedure bern!*calc(n, info);
begin scalar result;
result := single!*bern(n/2);
if !*rounded then
info := list (car info, c!:prec!:(), result)
else info := list (result, cadr info, caddr info);
bernoulli!-alist := (n . info) . bernoulli!-alist;
return result;
end;
%
% Computation of Bernoulli numbers using the algorithms of
% one Herbert S. Wilf, presented by Sandra Fillebrown in the
% Journal of Algorithms 13 (1992)
%
% Chris Cannam, October 1992
%
%
% Useful auxiliary fn.
%
symbolic procedure nearest!-int!-to!-bf(x);
(conv!:bf2i rb)
where rb = (if lessp!:(x,bfz!*)
then difference!:(x,bfhalf!*)
else plus!:(x,bfhalf!*));
%
% Procedure to compute B(2k) for k = 2 ... n
%
% Returns list of even bernoullis from B(4) to B(2n),
% in reverse order; only works when compiled, owing
% to reliance upon msd!:, which is a compiled inline
% macro.
%
% If called with rounded mode off, it computes the
% exact quotient; otherwise it will usually approximate
% (to the correct precision) if it saves time to do so.
%
symbolic procedure multi!*bern(n);
begin scalar results, primes, tprimes, r0, rk, rkm1, b2k,
tpi, pie, tk, n2k;
integer thisp, gn, prepre, prernd, p2k, k2, plim, d2k;
results := nil;
prernd := !*rounded;
if not prernd then on rounded;
prepre := precision 0;
if new!*bfs then
<< gn := 2 * n * msd!: n;
if gn < (log2of10*!!nfpd) then precision (!!nfpd + 2)
else if prepre > (gn/3) or not prernd then
precision (gn/3 + 1)
else precision (prepre + 2) >>
else
<< gn := 2 * n * length explode n;
if gn < !!nfpd then precision (!!nfpd + 2)
else if prepre > gn or not prernd then
precision (gn + 2)
else precision (prepre + 2) >>;
tpi := pi!*(); pie := divbf(bfone!*, timbf(tpi, e!*()));
if n < 1786 then primes := !*primelist!*
else
<< primes := nil;
for thisp := 3573 step 2 until (2*n + 1) do
if primep thisp then primes := thisp . primes;
primes := append(!*primelist!*, reverse primes) >>;
r0 := sq2bf!* algebraic ((2*pi)**(-2));
rkm1 := timbf(i2bf!: 4, r0);
for k := 2:n do
<< k2 := 2*k;
rk := timbf(i2bf!:(k2*(k2 - 1)), timbf(r0, rkm1));
rkm1 := rk;
tk := bfone!*; d2k := 1;
plim := 1 + conv!:bf2i timbf(i2bf!: k2, pie);
tprimes := cdr primes; thisp := car primes;
while thisp <= plim do
<< p2k := thisp ** k2;
tk := timbf(tk, divbf(i2bf!: p2k, i2bf!: (p2k - 1)));
thisp := car tprimes;
tprimes := cdr tprimes >>;
tprimes := cdr primes; thisp := car primes;
while thisp <= k+1 do
<< if cdr divide (k2, thisp - 1) = 0 then
d2k := d2k * thisp;
thisp := car tprimes;
tprimes := cdr tprimes >>;
if primep (k2 + 1) then d2k := d2k * (k2 + 1);
n2k := timbf(timbf(rk, tk), i2bf!: d2k);
if prernd then
b2k := mk!*sq !*f2q mkround
divbf (i2bf!: (((-1)**(k+1)) *
nearest!-int!-to!-bf n2k),
i2bf!: d2k)
else b2k := list ('!*sq, (((-1)**(k+1)) *
nearest!-int!-to!-bf n2k) . d2k, t);
results := b2k . results >>;
precision prepre;
if not prernd then off rounded;
return results;
end;
%
% Procedure to compute B(2n). If it is called with rounded
% mode off, it computes the exact quotient; otherwise it
% will approximate (to the correct precision) whenever it
% saves time to do so.
%
symbolic procedure single!*bern(n);
begin scalar result, primes, tprimes, rn, tn, n2n, pie;
integer d2n, thisp, gn, prepre, prernd, p2n, n2, plim;
prernd := !*rounded;
if not prernd then on rounded;
prepre := precision 0;
if new!*bfs then
<< gn := 2 * n * msd!: n;
if gn < (log2of10*!!nfpd) then precision (!!nfpd + 2)
else if prepre > (gn/3) or not prernd then
precision (gn/3 + 1)
else precision (prepre + 2) >>
else
<< gn := 2 * n * length explode n;
if gn < !!nfpd then precision (!!nfpd + 2)
else if prepre > gn or not prernd then
precision (gn + 1)
else precision (prepre + 2) >>;
pie := divbf(bfone!*, timbf(pi!*(), e!*()));
if n < 1786 then primes := !*primelist!*
else
<< primes := nil;
for thisp := 3573 step 2 until (2*n + 1) do
if primep thisp then primes := thisp . primes;
primes := append(!*primelist!*, reverse primes) >>;
n2 := 2*n;
rn := divbf(i2bf!: (2 * factorial n2),
sq2bf!* algebraic ((2*pi)**(n2)));
tn := bfone!*; d2n := 1;
plim := 1 + conv!:bf2i timbf(i2bf!: n2, pie);
tprimes := cdr primes; thisp := car primes;
while thisp <= plim do
<< p2n := thisp ** n2;
tn := timbf(tn, divbf(i2bf!: p2n, i2bf!: (p2n - 1)));
thisp := car tprimes;
tprimes := cdr tprimes >>;
tprimes := cdr primes; thisp := car primes;
while thisp <= n+1 do
<< if cdr divide (n2, thisp - 1) = 0 then
d2n := d2n * thisp;
thisp := car tprimes;
tprimes := cdr tprimes >>;
if primep (n2 + 1) then d2n := d2n * (n2 + 1);
n2n := timbf(timbf(rn, tn), i2bf!: d2n);
precision prepre;
if prernd then
result := mkround divbf(i2bf!: (((-1)**(n+1)) *
nearest!-int!-to!-bf n2n),
i2bf!: d2n)
else
<< off rounded;
result := list ('!*sq, (((-1)**(n+1)) *
nearest!-int!-to!-bf n2n) . d2n, t) >>;
return result;
end;
% Euler numbers module by Kerry Gaskell
algebraic operator Euler;
algebraic
let { Euler(0) => 1,
Euler(~n) => Euler_aux(n) when fixp n and n > 0};
flag('(euler_aux),'opfn);
symbolic procedure Euler_aux(n);
if not evenp n then 0 else
begin scalar nn,list_of_bincoeff, newlist, old, curr,eulers,sum;
list_of_bincoeff := { 1 };
eulers :={ -1,1};
nn := -2;
while N > 0 do
<< nn := nn + 1;
old := 0;
newlist := {};
while not(list_of_bincoeff = {}) do
<< curr := first list_of_bincoeff;
newlist := (old + curr) . newlist;
old := curr;
list_of_bincoeff := rest list_of_bincoeff;
>>;
list_of_bincoeff := 1 . newlist;
n := n -1 ;
% now that we have got the row of Pascal's triangle
% we can use it and compute the Next Euler number.
if nn > 0 and evenp nn then <<
curr := list_of_bincoeff;
old := eulers; sum := 0;
while old do <<
curr := cddr curr;
sum := sum - first old * first curr;
old := cdr old;
>>;
eulers := sum . eulers;
>>;
>>;
return first eulers;
end;
endmodule;
end;