[Axiom-math] Re: [Axiom-developer] special functions

[Francois Maltey]

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