Artifact 4e34d732187401f362d7188acbc82aef89e8f2e7885b0993b3f4b0370a449c8d:
- Executable file
r37/packages/specfn/sfbern.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 10570) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/specfn/sfbern.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 10570) [annotate] [blame] [check-ins using]
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;