[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Axiom-developer] special functions
From: |
Francois Maltey |
Subject: |
Re: [Axiom-developer] special functions |
Date: |
31 Jan 2006 14:10:41 +0100 |
User-agent: |
Gnus/5.09 (Gnus v5.9.0) Emacs/21.4 |
Hello Bill, and others !
>
> Comments from everyone would be most welcome.
>
So I take some time to explain my point of view as a former mupad user.
You show 2 problems which causes axiom-source difficult to read :
// A //
> A bigger issue seems to me to be the way Gamma(x) is currently
> impletement in Axiom. I was a little horrified to see that it
> is actually implemented in BOOT code. Horrified because while
> BOOT is (nearly) as readable as SPAD, it is much closer to the
> machine level and therefore correspondingly harder to debug.
> In my opinion BOOT has it's place at the deeper levels of the
> Axiom interpreter and the SPAD compiler itself, but I think this
> is completely unnecessary when it comes to Algebra code.
When I want to understand a function I make a lot of grep on *.spad files
in order to find some interessing lines.
I really do this for the sin function.
I don't imagine that I must read the boot-loader for the gamma function.
The mupad file are in trees ; all files have the same structure.
I send the skeleton of the main igamma function in mupad lib.
There are sin.mu, gamma.mu, igamma.mu, Ei.mu files in the SPECFUNC directory.
There are also some igamma implentation in
STDLIB/DIFF/igamma.mu for differential computation
STDLIB/CONJ/igamma.mu for conjugate
SERIES/igamma.mu for series
...
// B //
> The second thing is that I do not like the way that the special
> functions have been implemented as part of the DoubleFloat
> domain. This is quite in contrast to the way that the trignometric
> functions have been implemented. One consequence is that it is
> hard to control when one wants to treat Gamma(x) symbolically
> (as one might do in the case of symbolic integration etc.) and
> when to treat it numerically.
When I use axiom I ask (to me) the question :
<< how can I explain this idea (to my students) with very view words >>
And I see it's not so easy to do with Float, Expression Float, DoubleFloat :
> For example:
>
sin(1) : sin(1) Type: Expression Integer
sin(1)::Expression Float Right
(sin 1)::Float Wrong, Why ?
(2/3)::Float
Gamma(1) Type: DoubleFloat
I prefer Integer
Gamma(1,2) Expression Integer is right.
> I don't think there should be this difference in behaviour
> between Gamma(x), Gamma(a,x) and sin(x). Fixing this will
> likely require some deeper surgery in Axiom.
This is nasty.
Computer Algebra is so difficult
that elementary principes must be very sharp.
If-no debug becomes impossible as maple with structural type mistake as 0=0.0
// C // The skeleton of the igamma function :
igamma :=
proc(a, x)
option noDebug;
local fa, ceila, floora, fraca, t1, t2, am1, fx, k, absx, result;
save DIGITS;
begin
if args(0) < 2 then error("2 arguments expected"); end_if;
....
case type(x) /* mupad does map */
of DOM_SET do
of "_union" do
return(map(x, u-> igamma(a,u)))
end_case;
....
if testargs() then /* arithmetical type */
....
if x = infinity then /* for limits */
if x=0 then return(gamma(a)) end_if; /* easy properties */
if a=0 then return(Ei(x)) end_if;
if a=1 then return(exp(-x)) end_if;
if x=float(0) then return(gamma(float(a))) end_if; /* first float cases */
if a=float(0) then return(Ei(float(x))) end_if;
if a=float(1) then return(exp(-float(x))) end_if;
if a = float(1/2) or x = float (x) then /* numerical computation */
...
if a = 1/2 then
return(float(PI^(1/2))*erfc(float(x)^(1/2))) end_if;
....
if domtype(x) = DOM_FLOAT and x >= 0 then /* igamma numerical call */
if domtype(a) = DOM_FLOAT and a >= 0 then
....
//-----------------------------------------------
// normalization: use recursion formula to reduce
// index a to strip 0<= a < 1. Do not do this for
// complex a, because in this case float evaluation
// is not possible, anyway.
// Restrict this to |a| <= 500, because this
// would be much too expensive for huge |a|
if testtype(a, Type::Real) and specfunc::abs(a) <= 500 then
....
// Warning: dangerous float evaluation via recursion!
// This is likely to be subject to cancellation, so
// we need to boost DIGITS !
DIGITS := ....
if 1 < a then // normalize a by shifting the first argument
// to the region 0 < a < 1:
// use recursion:
// igamma(a, x) = (a-1)*igamma(a-1, x) + x^(a-1)*exp(-x)
// = (a-1)*(a-2)*..*(a-k)*igmma(a - k, x) + exp(-x)*(x^(a-1)
// + (a-1)*x^(a-2) + ..)
....
if a < 0 then // normalize a by shifting the first argument
// to the region 0 < a < 1:
....
end_proc:
// some others functions as :
igamma:= slot(igamma, "info",
"igamma -- the incomplete Gamma function [try ?igamma for details]"):
igamma:= slot(igamma, "conjugate",
loadproc(slot(igamma,"conjugate"),pathname("STDLIB","CONJ"),"igamma")):
igamma:= slot(igamma, "diff",
loadproc(slot(igamma,"diff"),pathname("STDLIB","DIFF"),"igamma")):
....
// end of file
I hope you understand the point of view of an axiom end-user.
François