Implementation of the Incomplete Beta function

Am Sonntag, den 25.01.2009, 11:19 -0500 schrieb Raymond Toy:
> Dieter Kaiser wrote:
> > (%i10) beta_incomplete(0.5,%i,1.5);
> > BIGFLOAT:  Unable to convert %i to a CL or BIGFLOAT number
> >  -- an error.  To debug this try debugmode(true);
> >   
> Dang.  I fixed it so that bigfloat:to converts %i to #c(0 1).  But
> beta_incomplete still fails here:
> 
> (BIGFLOAT:* (BIGFLOAT:TO ($GAMMA (TO A))) (BIGFLOAT:TO ($GAMMA (TO B))))
> 
> Here, B = #c(0 1).  Then (to b) is '$%i.  ($gamma $%i) returns '((%gamma
> simp) $%i).  bigfloat:to can't handle that.
> 
> I think the calls to $gamma need convert either the args or the results
> to some kind of float/bfloat before calling bigfloat:to.
> 
> Ray

Hello Ray,

I see the problem and I am thinking about it.

The first part of the beta-incomplete function uses a reflection
formula. Perhaps it is necessary to redesign this part completely
without using the package bigfloat.

I have checked the handling of constants. The results are interesting:

The flag numer gives a numerical result for %i, but not the function
float():

(%i30) gamma(%i),numer;
(%o30) - .4980156681183566 %i - .1549498283018101
(%i31) float(gamma(%i));
(%o31) gamma(%i)

That is different for other constants which give always a numerical
result:

(%i32) gamma(%gamma),numer;
(%o32) 1.543881741815662
(%i33) float(gamma(%gamma));
(%o33) 1.543881741815662
(%i34) gamma(%phi),numer;
(%o34) .8956731517052878
(%i35) float(gamma(%phi));
(%o35) .8956731517052878

But for %e we have a bug in the gamma function:

(%i36) gamma(%e),numer;
Maxima encountered a Lisp error:
CDR: $%E is not a list
Automatically continuing.
To reenable the Lisp debugger set *debugger-hook* to nil.

(%i37) float(gamma(%e));
(%o37) 1.567468255774053

I think the problem is the following line in the function simpgamma:

($numer (gammafloat (fpcofrat j)))

This is called without any further check. Furthermore I think simpgamma
has to be redesign a bit to be more correct in handling the flag $numer
and constants.

Some results for the sin function. The function float() works for all
constants:

(%i38) sin(%e);
(%o38) sin(%e)
(%i39) sin(%e),numer;
(%o39) sin(%e)
(%i40) float(sin(%e));
(%o40) .4107812905029088
(%i41) float(sin(%gamma));
(%o41) .5456928232039928
(%i42) sin(%gamma),numer;
(%o42) .5456928232039928
(%i43) sin(%i),numer;
(%o43) 1.175201193643801 %i
(%i44) float(sin(%i));
(%o44) 1.175201193643801 %i

Dieter Kaiser