diff --git a/changelog b/changelog index ccd8d2f..03f1d4b 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,6 @@ +20091029 tpd src/axiom-website/patches.html 20091029.03.tpd.patch +20091029 tpd src/input/Makefile add numericgamma.input +20091029 tpd src/input/numericgamma.input added 20091029 tpd src/axiom-website/patches.html 20091029.02.tpd.patch 20091029 tpd src/input/Makefile remove lexp.input 20091029 tpd src/input/lexp.input removed, duplicate of LEXP domain diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index ccc660b..67acdfe 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -2195,5 +2195,7 @@ books/bookvol10.3 Product add makeprod example
src/input/distexpr.input added
20091029.02.tpd.patch src/input/lexp.input removed, duplicate of LEXP domain
+20091029.03.tpd.patch +src/input/numericgamma.input added
diff --git a/src/input/Makefile.pamphlet b/src/input/Makefile.pamphlet index fe9fd07..388573a 100644 --- a/src/input/Makefile.pamphlet +++ b/src/input/Makefile.pamphlet @@ -347,11 +347,12 @@ REGRES= algaggr.regress algbrbf.regress algfacob.regress alist.regress \ lpoly.regress lupfact.regress lword.regress macbug.regress \ macros.regress magma.regress mapleok.regress \ matbug.regress mathml.regress \ - matrix1.regress matrix22.regress matrix.regress \ + matrix1.regress matrix22.regress matrix.regress \ mfinfact.regress mkfunc.regress mpoly.regress mset2.regress \ mset.regress multfact.regress multiple.regress ndftip.regress \ negfloats.regress nepip.regress newlodo.regress newton.regress \ - newtonlisp.regress nlode.regress nonlinhomodiffeq.regress \ + newtonlisp.regress numericgamma.regress \ + nlode.regress nonlinhomodiffeq.regress \ none.regress noonburg.regress noptip.regress \ nqip.regress nsfip.regress numbers.regress octonion.regress \ oct.regress ode.regress odpol.regress op1.regress \ @@ -647,7 +648,7 @@ FILES= ${OUT}/algaggr.input ${OUT}/algbrbf.input ${OUT}/algfacob.input \ ${OUT}/ndftip.input ${OUT}/newlodo.input \ ${OUT}/negfloats.input \ ${OUT}/nepip.input ${OUT}/newton.input ${OUT}/newtonlisp.input \ - ${OUT}/nonlinhomodiffeq.input \ + ${OUT}/numericgamma.input ${OUT}/nonlinhomodiffeq.input \ ${OUT}/nlode.input ${OUT}/none.input ${OUT}/noonburg.input \ ${OUT}/noptip.input ${OUT}/nqip.input ${OUT}/nsfip.input \ ${OUT}/ntube.input ${OUT}/oct.input ${OUT}/ode.input \ @@ -970,7 +971,7 @@ DOCFILES= \ ${DOC}/ndftip.input.dvi ${DOC}/negfloats.input.dvi \ ${DOC}/nepip.input.dvi ${DOC}/newlodo.input.dvi \ ${DOC}/newton.input.dvi ${DOC}/newtonlisp.input.dvi \ - ${DOC}/nlode.input.dvi \ + ${DOC}/numericgamma.input.dvi ${DOC}/nlode.input.dvi \ ${DOC}/none.input.dvi ${DOC}/nonlinhomodiffeq.input.dvi \ ${DOC}/noonburg.input.dvi \ ${DOC}/noptip.input.dvi ${DOC}/nqip.input.dvi \ diff --git a/src/input/numericgamma.input.pamphlet b/src/input/numericgamma.input.pamphlet new file mode 100644 index 0000000..337822c --- /dev/null +++ b/src/input/numericgamma.input.pamphlet @@ -0,0 +1,450 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/input numericgamma.input} +\author{Bill Page} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\begin{chunk}{*} +)set break resume +)sys rm -f numericgamma.output +)spool numericgamma.output +)set message test on +)set message auto off +)clear all +)sys cp $AXIOM/../../src/input/numericgamma.input.pamphlet . +)lisp (tangle "numericgamma.input.pamphlet" "sfx.spad" "sfx.spad") +)co sfx.spad + +--S 1 of 36 +Gam(a:Float,x:Float):Float == + if x < 0.0 or a < 0.0 then error "Invalid arguments" + if x = 0.0 then return Gamma(a) + + ITMAX ==> 100 -- Maximum allowed number of iterations + FPMIN ==> 1.0e-1000 -- near the smallest representable number + -- (there is no smallest representable float) + + EPS := (10.0^(-digits()$Float+1))$Float -- Relative accuracy + + an: Float + del: Float + + b:Float:=x+1.0-a -- Set up for evaluating continued fractions + c:Float:=1.0/FPMIN -- by modified Lentz's method + d:Float:=1.0/b -- with b_0 = 0 + h:Float:=d + i:=1 + repeat -- iterate to convergence + an:=-i*(i-a) + b:=b+2.0 + d:=an*d+b + if abs(d) < FPMIN then d:=FPMIN + c:=b+an/c; + if abs(c) < FPMIN then c:=FPMIN + d:=1.0/d + del:=d*c + h:=h*del + if i > ITMAX or abs(del-1.0) < EPS then break + i:=i+1 + if i > ITMAX then error("a too large, ITMAX too small") + exp(-x)*x^a*h -- put factors in front +--R +--R Function declaration Gam : (Float,Float) -> Float has been added to +--R workspace. +--R Type: Void +--E 1 + +--S 2 of 36 +Gam(0,1) +--R +--R Compiling function Gam with type (Float,Float) -> Float +--R +--RDaly Bug +--R Error signalled from user code in function Gam: +--R a too large, ITMAX too small +--E 2 + +--S 3 of 36 +Gam(1.1.1) +--R +--R There are 1 exposed and 1 unexposed library operations named elt +--R having 1 argument(s) but none was determined to be applicable. +--R Use HyperDoc Browse, or issue +--R )display op elt +--R to learn more about the available operations. Perhaps +--R package-calling the operation or using coercions on the arguments +--R will allow you to apply the operation. +--R +--RDaly Bug +--R Cannot find application of object of type Float to argument(s) of +--R type(s) +--R Float +--R +--E 3 + +--S 4 of 36 +Gam(5,10) +--R +--R +--R (2) 0.7020645138 4706574415 +--R Type: Float +--E 4 + +--S 5 of 36 +Gam(5,11) +--R +--R +--R (3) 0.3625104156 5228203538 +--R Type: Float +--E 5 + +--S 6 of 36 +Gam(7,0) +--R +--R +--R (4) 720.0000000000 0011369 +--R Type: Float +--E 6 + +--S 7 of 36 +digits 100 +--R +--R +--R (5) 20 +--R Type: PositiveInteger +--E 7 + +--S 8 of 36 +Gam(0,1) +--R +--R +--RDaly Bug +--R Error signalled from user code in function Gam: +--R a too large, ITMAX too small +--E 8 + +--S 9 of 36 +Gam(1,1.1) +--R +--R +--R (6) +--R 0.3328710836 9807955328 8846906431 3155216124 7952156921 2491793331 386750747 +--R 0 8541284431 1612617072 7005478519 +--R Type: Float +--E 9 + +--S 10 of 36 +Gam(1,1) +--R +--R +--R (7) +--R 0.3678794411 7144232159 5523770161 4608674458 1113103176 7834507836 801697461 +--R 4 9574489980 3357147274 3459196437 +--R Type: Float +--E 10 + +--S 11 of 36 +Gam(1,1.1) +--R +--R +--R (8) +--R 0.3328710836 9807955328 8846906431 3155216124 7952156921 2491793331 386750747 +--R 0 8541284431 1612617072 7005478519 +--R Type: Float +--E 11 + +--S 12 of 36 +Gam(5,10) +--R +--R +--R (9) +--R 0.7020645138 4706574414 6387196628 3546367191 6532623256 0684622278 670587055 +--R 0 5584357048 3474646670 2985365058 +--R Type: Float +--E 12 + +--S 13 of 36 +Gam(5,11) +--R +--R +--R (10) +--R 0.3625104156 5228203538 0753904311 4079803866 4530925132 7036797697 419049037 +--R 4 2658968752 0305953551 1648548436 +--R Type: Float +--E 13 + +\end{chunk} +Note that in the case of Gam(7,0) the function Gamma(a) is called in the +above code but it is actually implemented in Axiom as Gamma(a)\$DoubleFloat +with a fixed precision so the precision of the result is not affected by +the setting of digits(). +\begin{chunk}{*} + +--S 14 of 36 +Gam(7,0) +--R +--R +--R (11) 720.0000000000 0011368683 7721616029 7393798828 125 +--R Type: Float +--E 14 + +--S 15 of 36 +Gam(7,0.1) +--R +--R +--R (12) +--R 719.9999999869 1035963050 9717349089 5137595484 2683243460 6577519316 5312727 +--R 417 6619922456 9102294155 8764196 +--R Type: Float +--E 15 + +--S 16 of 36 +Gam(7,0.2) +--R +--R +--R (13) +--R 719.9999984646 1597708521 8246915701 3222705579 4693807497 2229652513 6047137 +--R 980 7138425860 0596921944 0451807 +--R Type: Float +--E 16 + +\end{chunk} +\begin{chunk}{sfx.spad} +)abbrev package SFX SpecialFunctionExtended +SpecialFunctionExtended: Exports == Implementation where + ITMAX ==> 100.0::DoubleFloat -- Maximum allowed number of iterations + FPMIN ==> 1.0e-323::DoubleFloat -- near the smallest representable number + + Exports ==> with + NGamma:(DoubleFloat,DoubleFloat)->DoubleFloat + + Implementation ==> add + + NGamma(a:DoubleFloat,x:DoubleFloat):DoubleFloat == + if x < 0 or a < 0 then error "Invalid arguments" + if x = 0 then return Gamma(a) + EPS := (10.0::DoubleFloat**(-digits()$DoubleFloat))$DoubleFloat + b:DoubleFloat:=x+1-a -- Set up for evaluating continued fractions + c:DoubleFloat:=1/FPMIN -- by modified Lentz's method + d:DoubleFloat:=1/b -- with b_0 = 0 + h:DoubleFloat:=d + i:DoubleFloat:=1 + repeat -- iterate to convergence + an:DoubleFloat:=-i*(i-a) + b:=b+2.0::DoubleFloat + d:=an*d+b + if abs(d) < FPMIN then d:=FPMIN + c:=b+an/c + if abs(c) < FPMIN then c:=FPMIN + d:=1/d + del:DoubleFloat:=d*c + h:=h*del + if i > ITMAX or abs(del-1) < EPS then break + i:=i+1 + if i > ITMAX then error("a too large, ITMAX too small") + exp(-x+log(x)*a)*h -- put factors in front +\end{chunk} + +\begin{chunk}{*} + +--S 17 of 36 +NGamma(a,x) +--R +--R There are 1 exposed and 0 unexposed library operations named NGamma +--R having 2 argument(s) but none was determined to be applicable. +--R Use HyperDoc Browse, or issue +--R )display op NGamma +--R to learn more about the available operations. Perhaps +--R package-calling the operation or using coercions on the arguments +--R will allow you to apply the operation. +--R +--RDaly Bug +--R Cannot find a definition or applicable library operation named +--R NGamma with argument type(s) +--R Variable a +--R Variable x +--R +--R Perhaps you should use "@" to indicate the required return type, +--R or "$" to specify which version of the function you need. +--E 17 + +--S 18 of 36 +NGamma(0,1) +--R +--R +--R (14) 0.21938393439551901 +--R Type: DoubleFloat +--E 18 + +--S 19 of 36 +NGamma(0,2) +--R +--R +--R (15) 4.8900510708060993E-2 +--R Type: DoubleFloat +--E 19 + +--S 20 of 36 +NGamma(1,1) +--R +--R +--R (16) 0.36787944117144233 +--R Type: DoubleFloat +--E 20 + +--S 21 of 36 +NGamma(1,1.1) +--R +--R +--R (17) 0.33287108369807966 +--R Type: DoubleFloat +--E 21 + +--S 22 of 36 +NGamma(5,10) +--R +--R +--R (18) 0.70206451384706692 +--R Type: DoubleFloat +--E 22 + +--S 23 of 36 +NGamma(5,11) +--R +--R +--R (19) 0.36251041565228215 +--R Type: DoubleFloat +--E 23 + +--S 24 of 36 +NGamma(7,0) +--R +--R +--R (20) 720.00000000000011 +--R Type: DoubleFloat +--E 24 + +--S 25 of 36 +NGamma(7,0.1) +--R +--R +--R (21) 719.99620051670286 +--R Type: DoubleFloat +--E 25 + +--S 26 of 36 +NGamma(7,0.2) +--R +--R +--R (22) 719.99974844402846 +--R Type: DoubleFloat +--E 26 + +)set functions compile on + +--S 27 of 36 +j:=120 +--R +--R +--R (23) 120 +--R Type: PositiveInteger +--E 27 + +--S 28 of 36 +nume(a) == cons(1::Float,[((a-i)*i)::Float for i in 1..]) +--R +--R Type: Void +--E 28 + +--S 29 of 36 +dene(a,x) == [(x+2*i+1-a)::Float for i in 0..] +--R +--R Type: Void +--E 29 + +--S 30 of 36 +cfe(a,x) == continuedFraction(0,nume(a),dene(a,x)) +--R +--R Type: Void +--E 30 + +--S 31 of 36 +ccfe(a,x) == convergents cfe(a,x) +--R +--R Type: Void +--E 31 + +--S 32 of 36 +gamcfe(a,x) == exp(-x)*x^a*(ccfe(a,x).j)::Float +--R +--R Type: Void +--E 32 + +--S 33 of 36 +gamcfe(2,3) +--R +--R Compiling function nume with type PositiveInteger -> Stream Float +--R Compiling function dene with type (PositiveInteger,PositiveInteger) +--R -> Stream Float +--R Compiling function cfe with type (PositiveInteger,PositiveInteger) +--R -> ContinuedFraction Float +--R Compiling function ccfe with type (PositiveInteger,PositiveInteger) +--R -> Stream Fraction Float +--R Compiling function gamcfe with type (PositiveInteger,PositiveInteger +--R ) -> Expression Float +--R +--R (29) +--R 0.1991482734 7145577191 7369662600 2471065267 9836875369 2862270510 910424242 +--R 6 7092079820 0616216976 9465333782 +--R Type: Expression Float +--E 33 + +--S 34 of 36 +E1fun(x) == gamcfe(0,x) +--R +--R Type: Void +--E 34 + +--S 35 of 36 +E1fun(2.0) +--R +--R Compiling function nume with type NonNegativeInteger -> Stream Float +--R +--R Compiling function dene with type (NonNegativeInteger,Float) -> +--R Stream Float +--R Compiling function cfe with type (NonNegativeInteger,Float) -> +--R ContinuedFraction Float +--R Compiling function ccfe with type (NonNegativeInteger,Float) -> +--R Stream Fraction Float +--R Compiling function gamcfe with type (NonNegativeInteger,Float) -> +--R Float +--R Compiling function E1fun with type Float -> Float +--R +--R (31) +--R 0.0489005107 0806111956 7239826914 3472898212 1544510421 3277251841 716377988 +--R 0 9149832755 9949235928 1965882172 4 +--R Type: Float +--E 35 + +--S 36 of 36 +E1fun(2.0)-E1(2.0) +--R +--R +--R (32) 1.1102230246251565E-16 +--R Type: OnePointCompletion DoubleFloat +--E 36 + +)spool +)lisp (bye) + +\end{chunk} +\eject +\begin{thebibliography}{99} +\bibitem{1} http://axiom-wiki.newsynthesis.org/SandBoxGamma +\end{thebibliography} +\end{document}