Artifact 9bfc6e6b14e631a289e735d80e7421c2678167c48a3542aea5c65096f63ebd74:
- Executable file
r37/packages/specfn/sfpsi.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: 22489) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/specfn/sfpsi.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: 22489) [annotate] [blame] [check-ins using]
module sfpsi; % Procedures relevant to the digamma, polygamma & zeta % functions. % Author: Chris Cannam, Sept/Oct '92. % Added: PSI_SIMP.RED F.J.Wright, 2 July 1993 % The polygamma rules are added by Y.K. Man on 9 July 1993 % Yiu K. Man's email is: myk@maths.qmw.ac.uk imports sq2bf!*, sf!*eval; exports do!*psi, do!*polygamma, do!*trigamma!*halves, do!*zeta, do!*zeta!*pos!*intcalc; % % A couple of global values are used (from specfns.red) which can speed % up psi calculations (a bit) when repeated calculations are made at the % same level of precision. fluid '(compute!-bernoulli); % % Here's an approximation sufficiently good for most purposes % (assuming it's right, that is); if it isn't good enough, it % won't be used. This approximation is to 506 dec. places. % algebraic (old!*precision := precision(0)); precision 510; algebraic procedure get!-eulers!-constant; begin scalar a; a := 577215664901532860606512090082402431 * 10^40 + 0421593359399235988057672348848677267776; a := a * 10^40 + 6467093694706329174674951463144724980708; a := a * 10^40 + 2480960504014486542836224173997644923536; a := a * 10^40 + 2535003337429373377376739427925952582470; a := a * 10^40 + 9491600873520394816567085323315177661152; a := a * 10^40 + 8621199501507984793745085705740029921354; a := a * 10^40 + 7861466940296043254215190587755352673313; a := a * 10^40 + 9925401296742051375413954911168510280798; a := a * 10^40 + 4234877587205038431093997361372553060889; a := a * 10^40 + 3312676001724795378367592713515772261027; a := a * 10^40 + 3492913940798430103417771778088154957066; a := a * 10^30 + 107501016191663340152278935868; a := a * (10**(-506)); return a end; algebraic (euler!*constant := get!-eulers!-constant()); algebraic precision old!*precision; algebraic clear old!*precision; % % Define some suitable rules for initial simplification of psi % (digamma) function expressions. % % Comments: % % When rounded mode is on, psi(number) is always computed % directly unless it simplifies to an expression in psi(1/2) or % psi(1), in which case it is simplified. Expressions in psi(1/2) % and psi(1) are expanded into expressions in euler!*constant. % If, however, the precision is greater than 500, then % euler!*constant is not stored sufficiently precisely, and all % such expressions will be computed without simplification. % % When rounded mode is off, psi(number) will _never_ be expanded % into an expression involving euler!*constant, but will always % be expanded into some form involving psi(p), where 0<p<1. % It should be borne in mind that computations which will need % numerical results could do without such expansion, and there- % fore such computations should be performed in rounded mode % as soon as possible. % % Expressions for the derivative and integral of psi are included. % algebraic operator psi, polygamma; symbolic operator psi!*calc; algebraic (psi!*rules := { psi(~x,~xx) => polygamma(x,xx), psi(~z) => infinity when repart z = floor repart z and impart z = 0 and z < 1, psi(~z) => -euler!*constant when numberp z and z = 1 and symbolic !*rounded and precision(0) < 501, psi(~z) => -euler!*constant - 2 * log(2) when numberp z and z = (1/2) and symbolic !*rounded and precision(0) < 501, psi(~z) => do!*psi(z) when numberp z and impart z = 0 and symbolic !*rounded, psi(~z) => (psi(z/2) + psi((z+1)/2) + 2 * log(2)) / 2 when numberp z and impart z = 0 and (z/2) = floor (z/2) and z > 0 and not symbolic !*rounded, psi(~z) => psi(z-1) + (1 / (z-1)) when numberp z and impart z = 0 and z > 1 and not symbolic !*rounded, psi(~z) => psi(1-z) + pi*cot(pi*(1-z)) when numberp z and impart z = 0 and z < 0 and not symbolic !*rounded, psi(~z) => psi(1-z) + pi*cot(pi*(1-z)) when numberp z and impart z = 0 and z > 1/2 and z < 1 and not symbolic !*rounded, df(psi(~z),z) => polygamma(1, z), int(psi(~z),z) => log gamma(~z) })$ algebraic (let psi!*rules); % PSI_SIMP.RED F.J.Wright, 2 July 1993 % The polygamma rules are added by Y.K. Man on 9 July 1993 % Support for the psi operator. % ============================= % psi(x) = df(log Gamma(x), x) as in specfn package, etc. % The specfn package does not currently provide the required % simplifications. algebraic; % Simplify to "standard form" in which argument is allowed a numeric % shift in the range 0 <= shift < 1: psi_rules := { % Rule for integer shifts (x + 3), and non-integer shifts (x + 3/2)in % a non-integer number domain (on rational) or with "on intstr, div": psi(~x+~n) => psi(x+n-1) + 1/(x+n-1) when numberp n and n >= 1, psi(~x+~n) => psi(x+n+1) - 1/(x+n) when numberp n and n < 0, polygamma(~m,~x+~n) => polygamma(m,x+n-1)+(-1)^m*factorial(m) /(x+n-1)^(m+1) when numberp n and fixp m and n >= 1, polygamma(~m,~x+~n) => polygamma(m,x+n+1)-(-1)^(m)*factorial(m) /(x+n)^(m+1) when numberp n and fixp m and n < 0, % Rule for rational shifts (x + 3/2) in the default (integer) number % domain and rational arguments (x/y + 3): psi((~x+~n)/~d) => psi((x+n-d)/d) + d/(x+n-d) when numberp(n/d) and n/d >= 1, psi((~x+~n)/~d) => psi((x+n+d)/d) - d/(x+n) when numberp(n/d) and n/d < 0, polygamma(~m,(~x+~n)/~d) => polygamma(m,(x+n-d)/d) + (-1)^m*factorial(m)*d^(m+1)/(x+n-d)^(m+1) when fixp m and numberp(n/d) and n/d >= 1, polygamma(~m,(~x+~n)/~d) => polygamma(m,(x+n+d)/d) - (-1)^m*factorial(m)*d^(m+1)/(x+n)^(m+1) when fixp m and numberp(n/d) and n/d < 0 }; % NOTE: The rational-shift rule does not work with "on intstr, div". let psi_rules; symbolic; % % Rules for initial manipulation of polygamma functions. % symbolic (operator polygamma!*calc, trigamma!*halves, printerr, polygamma_aux); symbolic procedure printerr(x); rederr x; algebraic procedure polygamma_aux(n,m); for ii:=1:(n-1) sum (1/ii**(m+1)); algebraic (polygamma!*rules := { polygamma(~n,~x) => printerr "Index of Polygamma must be an integer >= 0" when numberp n and (not fixp n or n < -1), polygamma(~n,~x) => psi(x) when numberp n and n = 0, polygamma(~n,~x) => infinity when numberp x and impart x = 0 and x = floor x and x < 1, polygamma(~n,~x) => do!*trigamma!*halves(x) when numberp n and n = 1 and numberp x and impart x = 0 and (not (x = floor x) and ((2*x) = floor (2*x))) and x > 1, polygamma(~n,~x) => ((-1) ** (n)) * (factorial n) * (- zeta(n+1) + polygamma_aux(x,n)) when fixp x and x >= 1 and not symbolic !*rounded, polygamma(~n,~x) => ((-1)**n) * factorial n * (-2 * (2**n) * zeta(n+1) + 2 * (2**n) + zeta(n+1)) when numberp x and x = (3/2) and not symbolic !*rounded, polygamma(~n,~x) => do!*polygamma(n,x) when numberp x and symbolic !*rounded and numberp n and impart n = 0 and n = floor n, df(polygamma(~n,~x), ~x) => polygamma(n+1, x), int(polygamma(~n,~x),~x) => polygamma(n-1,x) })$ algebraic (let polygamma!*rules); % % Set up rules for the initial manipulation of zeta. % % Comments: % % Zeta of positive even numbers and negative odd numbers % is evaluated (in terms of pi) always when its argument % has magnitude less than 31, and only in rounded mode % otherwise. (This is because the coefficients get a bit % big when the argument is over about 30.) % algebraic operator zeta; symbolic (operator zeta!*calc, zeta!*pos!*intcalc); algebraic (zeta!*rules := { zeta(~x) => (- (1/2)) when numberp x and x = 0, zeta(~x) => (pi ** 2) / 6 when numberp x and x = 2, zeta(~x) => (pi ** 4) / 90 when numberp x and x = 4, zeta(~x) => infinity when numberp x and x = 1, zeta(~x) => 0 when numberp x and impart x = 0 and x < 0 and (x/2) = floor(x/2), zeta(~x) => ((2*pi)**x) / (2*factorial x)*(abs bernoulli!*calc x) when numberp x and impart x = 0 and x > 0 and (x/2) = floor (x/2) and x < 31, zeta(~x) => - (bernoulli!*calc (1-x)) / (2*x) when numberp x and impart x = 0 and x < 0 and x = floor x and x > -31, zeta(~x) => ((2*pi)**x)/(2 * factorial x)*(abs bernoulli!*calc x) when numberp x and impart x = 0 and x > 0 and (x/2) = floor(x/2) and x < 201 and symbolic !*rounded, zeta(~x) => - (bernoulli!*calc (1-x)) / (1-x) when numberp x and impart x = 0 and x < 0 and x = floor x and x > -201 and symbolic !*rounded, zeta(~x) => (2**x)*(pi**(x-1))*sin(pi*x/2)*gamma(1-x)*zeta(1-x) when numberp x and impart x = 0 and x < 0 and (x neq floor x or x < -200) and symbolic !*rounded, zeta(~x) => do!*zeta!*pos!*intcalc(fix x) when symbolic !*rounded and numberp x and impart(x) = 0 and x > 1 and x = floor x and (x <= 15 or precision 0 > 100 or 2*x < precision 0), zeta(~x) => do!*zeta(x) when numberp x and impart x = 0% and x > 1 and symbolic !*rounded, df(zeta(~x),x) => -(1/2)*log(2*pi) when numberp x and x = 0 })$ algebraic (let zeta!*rules); algebraic procedure do!*psi(z); algebraic sf!*eval('psi!*calc,{z}); algebraic procedure do!*polygamma(n,z); algebraic sf!*eval('polygamma!*calc,{n,z}); algebraic procedure do!*trigamma!*halves(z); algebraic sf!*eval('trigamma!*halves,{z}); algebraic procedure do!*zeta(z); (if z <= 1.5 and precision(0) <= floor(4+3*z) then raw!*zeta(z) else if (3*z) > (10*precision(0)) then 1.0 else if z > 100 then algebraic sf!*eval('zeta!*calc,{z}) else algebraic sf!*eval('zeta!*general!*calc,{z})); algebraic procedure do!*zeta!*pos!*intcalc(z); algebraic sf!*eval('zeta!*pos!*intcalc,{z}); % % algebraic procedure psi!*calc(z); % % Compute a value of psi. Works by first computing the % smallest positive integral x at which psi(x) is easily % computable to the current precision using no more % than the first 200 bernoulli numbers, then scaling up % the given argument (if necessary) so that it can be % computed, scaling down again afterwards. % % Does not work for complex arguments. % algebraic procedure psi!*calc(z); begin scalar result, admissable, bern300, alglist!*, precom; integer prepre, scale, lowest; precom := complex!*off!*switch(); prepre := precision 0; if prepre < !!nfpd then precision (!!nfpd + 1); admissable := (1 / (10 ** prepre)) / 2; if prepre = psi!*ld(0) then lowest := psi!*ld(1) else << bern300 := abs bernoulli!*calc 300; lowest := 1 + symbolic conv!:bf2i exp!: (divbf(log!:(divbf(sq2bf!* bern300, timbf(i2bf!: 150, sq2bf!* admissable)), 4), i2bf!: 300), 3); % Use symbolic mode so as to % force less accuracy for more speed psi!*ld(0) := prepre; psi!*ld(1) := lowest >> ; if lowest>repart z then scale := ceiling(lowest - repart z) + 20; z := z + scale; result := algebraic symbolic psi!*calc!*sub(z, scale, admissable); precision prepre; complex!*restore!*switch(precom); return result; end; symbolic procedure psi!*calc!*sub(z, scale, admissable); begin scalar result, zsq, zsqp, this, bk; integer k, orda, rp; k := 2; z := sq2bf!* z; admissable := sq2bf!* admissable; zsq := timbf(z,z); zsqp := zsq; this := plubf(admissable, bfone!*); result := difbf (log!: (z,c!:prec!:()), divbf(bfone!*, timbf(bftwo!*, z))); orda := order!: admissable - 5; rp := c!:prec!:(); while greaterp!: (abs!: this, admissable) do << bk := sq2bf!* symbolic algebraic bernoulli!*calc k; this := divide!:(bk, timbf(i2bf!: k, zsqp), rp); result := difbf(result, this); k := k + 2; rp := order!: this - orda; zsqp := timbf(zsqp, zsq) >>; for n := 1:scale do result := difbf(result, divbf(bfone!*, difbf(z, i2bf!: n))); return mk!*sq !*f2q mkround result; end; % % algebraic procedure polygamma!*aux(n,z); % % Used by the procedure below, to implement the Reflection % Formula. This obtains an expression for % n % d % --- ( cot ( pi * x ) ) % n % dx % and substitutes z for x into it, returning the result. % algebraic procedure polygamma!*aux(n,z); begin scalar poly; clear dummy!*arg; poly := cot(pi * dummy!*arg); for k := 1:n do poly := df(poly, dummy!*arg); dummy!*arg := z; return poly; end; % % algebraic procedure polygamma!*calc(n,z); % % Computes a value of the polygamma function, order n, % at z. N must be an integer, and z must be real. If % z is negative, the Reflection Formula is applied by % a call to polygamma!*aux (above); then the positive % argument is fed to polygamma!*calc!*s which does the % real work. % algebraic procedure polygamma!*calc(n,z); begin scalar result, z0, prepre, precom; precom := complex!*off!*switch(); prepre := precision 0; if prepre < !!nfpd then precision (!!nfpd + 3) else precision (prepre + 3 + floor(prepre/50)); if z > 0 then << z0 := z; result := algebraic symbolic polygamma!*calc!*s(n,z0) >> else << z0 := 1-z; result := ((-1)**n)*(pi*polygamma!*aux(n,z0) + algebraic symbolic polygamma!*calc!*s(n,z0)) >>; precision prepre; complex!*restore!*switch(precom); return result; end; % % symbolic procedure polygamma!*calc!*s(n,z); % % Implementation of an asymptotic series for the poly- % gamma functions. Computes a scale factor which should % (hopefully) provide a minimum argument for which this % series is valid at the given order and precision; then % computes the series for that argument and scales down % again using the Recurrence Formula. % symbolic procedure polygamma!*calc!*s(n,z); begin scalar result, this, admissable, partial, zexp, zexp1, zsq, nfac, nfac1, kfac, rescale, signer, z0; integer k, nm1, nm2, rp, orda, min, scale; z := sq2bf!* z; signer := i2bf!:((-1)**(n-1)); admissable := divide!:(bfone!*,i2bf!:(bf!*base**c!:prec!:()),8); min := 10 + conv!:bf2i exp!:(times!:(divide!:(bfone!*,i2bf!:(300+n),8), log!:(divide!:(timbf(round!:mt(i2bf!: factorial(300+n),8), abs!: sq2bf!* symbolic algebraic bernoulli 300), times!:(admissable,round!:mt(i2bf!: factorial 300,8)), 8),8)),8); % In which Chris approximates to 8 bits % and hopes to get away with it... scale := min - (1 + conv!:bf2i z); if scale < 0 then scale := 0; z0 := plubf(z,i2bf!: scale); nfac := round!:mt(i2bf!: factorial(n-1),c!:prec!:()); zexp := texpt!:any(z0,n); result := plubf(divbf(nfac,zexp), divbf((nfac1 := timbf(i2bf!: n,nfac)), timbf(bftwo!*,(zexp1 := timbf(zexp,z0))))); nfac := nfac1; zexp := zexp1; nm1 := n-1; nm2 := n-2; rp := c!:prec!:(); nfac := timbf(nfac, i2bf!: (n+1)); kfac := bftwo!*; zexp := timbf(zexp,z0); zsq := timbf(z0,z0); partial := divbf(nfac,timbf(kfac,zexp)); k := 2; orda := order!: admissable - 5; this := bfone!*; if null compute!-bernoulli then <<errorset!*('(load_package '(specfaux)), nil); nil>>; while greaterp!:(abs!: this, admissable) do << result := plubf(result, (this := timbf(sq2bf!* retrieve!*bern k,partial))); k := k + 2; partial := divide!:(timbf(partial,i2bf!:((nm2+k)*(nm1+k))), timbf(zsq,i2bf!:((k-1)*k)),rp); rp := order!: this - orda >>; result := times!:(signer,result); if scale neq 0 then << rescale := bfz!*; nfac := round!:mt(i2bf!: factorial n,c!:prec!:()); for k := 1:scale do <<rescale := plus!:(rescale,timbf(nfac,texpt!:(z,-n-1))); z := plubf(z,bfone!*) >>; result := plubf(result,times!:(signer,rescale)) >>; return mk!*sq !*f2q mkround result; end; % % algebraic procedure trigamma!*halves(x); % % Applies a formula to derive the exact value of the trigamma % function at x where x = n+(1/2) for n = 1, 2, ... % algebraic procedure trigamma!*halves(x); begin integer prepre; scalar result, alglist!*; result := (1/2) * (pi ** 2) - (4 * (for k := 1:(round (x-(1/2))) sum ((2*k - 1) ** (-2)))); return result; end; % % algebraic procedure zeta!*calc(s); % % Calculate zeta(s). Only valid for repart(s) > 1. % % This function uses the system !*primelist!* of the first % 500 primes. If the system variable disappears or changes, % this function is helpless. % algebraic procedure zeta!*calc(z); begin scalar result, admissable, primelist, partialpl, this, modify, spl, alglist!*; integer prepre, j, rflag, thisprime, nexti; share spl; prepre := precision(0); precision prepre + 3; admissable := (1 / (10 ** (prepre + 2))); symbolic (spl := !*primelist!*); primelist := {}; result := 1; modify := 1; for k := 1:10 do << j := symbolic car spl; symbolic (spl := cdr spl); primelist := (j . primelist); modify := modify * (1 - (1 / (j**z))) >>; modify := 1 / modify; this := admissable + 1; if not symbolic cdr divide (j, 3) then j := j + 2; nexti := (if not symbolic cdr divide (j+1, 3) then 2 else 4); while ((abs this) > admissable) do << rflag := 1; partialpl := primelist; while ((partialpl neq {}) and rflag) do << thisprime := first partialpl; rflag := symbolic cdr divide(j, thisprime); partialpl := rest partialpl >>; if rflag then result := result + (this := (1 / (j**z))); j := j + nexti; nexti := 6 - nexti >>; result := result * modify; precision prepre; return result; end; algebraic procedure zeta!*pos!*intcalc(m); (((-1)**m)*polygamma(m-1,3)/factorial(m-1) + 1 + (1/(2**m))); algebraic procedure zeta!*error(z,terms); (((-1) ** (terms+2)) / ((terms+1) ** z)); algebraic procedure zeta!*general!*calc(z); begin scalar result, zp, admissable, z0; integer pre, k; pre := precision(0); admissable := algebraic symbolic (mk!*sq !*f2q mkround divide!:(bfone!*,i2bf!:(10 ** pre),8)); if (z**2) < admissable then result := ((-1/2) - ((log(2*pi))*z)/2) else if pre < !!nfpd then begin scalar sstt, stt; sstt := (for k := 2:(pre-1) sum (k**(-z))); precision (!!nfpd + 2); z0 := z; zp := pre**(-z); stt := sstt + 1; result := algebraic symbolic zeta!*general!*calc!*sub(z0,zp,admissable,pre,stt); end else <<z0 := z; zp := pre**(-z); result := algebraic symbolic zeta!*general!*calc!*sub(z0,zp,admissable,pre,'())>>; precision pre; return result; end; symbolic procedure zeta!*general!*calc!*sub(z,zp,admissable,pre,stt); begin scalar result, prere, this, fac, pre, zk1, zk2, logz; integer k; z := sq2bf!* z; zp := sq2bf!* zp; admissable := sq2bf!* admissable; if stt = nil then << result := bfone!*; k := 1; this := plus!:(admissable,bfone!*); while greaterp!: (abs!: this,admissable) and k < pre-1 do << k := k + 1; this := texpt!:any(i2bf!: k, minus!: z); result := plubf(result, this) >> >> else result := sq2bf!* stt; pre := i2bf!: pre; zk1 := plubf(z,bftwo!*); zk2 := plubf(z,bfone!*); result := plubf(result, timbf(zp,plubf(bfhalf!*,divbf(pre,difbf(z,bfone!*))))); fac := divbf(bfone!*,timbf(pre,pre)); this := timbf(divbf(z,bftwo!*),divbf(zp,pre)); result := plubf(result,divbf(this,i2bf!: 6)); k := 4; prere := plubf(result,bfone!*); while greaterp!: (abs!: difbf(prere,result), admissable) do << this := divbf(timbf(this,timbf(fac,timbf(zk1,zk2))), i2bf!:(k*(k-1))); prere := result; result := plubf(result,timbf( sq2bf!* symbolic algebraic bernoulli!*calc k, this)); zk1 := plus!:(zk1,bftwo!*); zk2 := plus!:(zk2,bftwo!*); k := k + 2; >>; return mk!*sq !*f2q mkround result; end; algebraic array stieltjes (5); % for use in raw zeta computations algebraic array stf (5); algebraic array psi!*ld (1); algebraic (psi!*ld(0) := -1); % precision at which last psi was calc'd algebraic (psi!*ld(1) := 0); % lowest post-scale value acceptable at % that precision Stieltjes (0) := 0.577215664901532860606512$ % Euler's constant Stieltjes (1) := -0.0728158233766$ Stieltjes (2) := -0.00968992187973$ Stieltjes (3) := 0.00206262773281$ Stieltjes (4) := 0.00250054826029$ Stieltjes (5) := 0.00427794456482$ Stf (0) := 1$ Stf (1) := 1$ Stf (2) := 2$ Stf (3) := 6$ Stf (4) := 24$ Stf (5) := 120$ algebraic procedure raw!*zeta(z); << z := z-1; 1/z + (for m := 0:5 sum ((-1)**m * Stieltjes(m) * z**m / Stf(m))) >>; endmodule; end;