[Axiom-developer] special functions

[Bill Page]

Yigal,

On January 30, 2006 1:11 AM on you wrote:
> 
> I think the numerical function in SPAD for gamma(a,x) should
> have several parts.  One for x>a+1, and a few others.  Is
> this appropriate?
> 

Well, the best way to implement gamma is probably exactly
the way it is done in Recipes - by first implementing the
functions lngamma, P and Q. If you wanted to do that, I think
it would be ok. In general it would be good to be able to
use the algorithms more or less directly from Recipes since
these are so very well known and a lot of thought has gone
into their development.

But the version of gamma that I adapted for SPAD is probably
"good enough" in most cases. The only real issue is that it
does not converge are rapidly as the series approximation in
the case x<=a+1. But the test I did suggest that this is not
really a problem unless you are very worried about performance
and this is not an issue on most modern computers.

I do currently have a technical problem with the SPAD code for
the case of gamma(a,0) - see the web page.

http://wiki.axiom-developer.org/SandBoxGamma

It seems to be an issue with the SPAD compiler since I can
successfully compile the same code using the Aldor compiler.
Unfortunately Aldor is not yet easily available with the open
source version of Axiom.

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.

http://wiki.axiom-developer.org/axiom--test--1/src/interp/SfsfunBoot

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.

For example:

(1) -> sin(1)

   (1)  sin(1)
                                  Type: Expression Integer
(2) -> sin(1)::Expression Float

   (2)  0.8414709848 0789650665
                                    Type: Expression Float
(3) -> Gamma(1)

   (3)  1.0
                                         Type: DoubleFloat
(4) -> Gamma(1,2)
         _
   (4)  | (1,2)
                                  Type: Expression Integer

-----------

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.

Comments from everyone would be most welcome.

Regards,
Bill Page.

See also: Related thread on axiom-math.