Re: binocdf inaccuracy in Octave

[Daniel J Sebald]

On 07/08/2013 05:03 PM, Pascal Dupuis wrote:
> 2013/7/8 Daniel J Sebald<address@hidden>:
> > I couldn't find that.  All I see in xdbetai.f is:
> >
> >        subroutine xdbetai (x, a, b, result)
> >        external dbetai
> >        double precision x, a, b, result, dbetai
> >        result = dbetai (x, a, b)
> >        return
> >        end
>
> Yeah, this is the C-to-fortran interface. The computation is performed
> inside liboctave/cruft/slatec-fn/dbetai.f

Oh, I thought it would be in some library somewhere.  Thanks...

Well, it does appear that internally the routine address numerical
aspects somewhat:

      Y = X
      P = PIN
      Q = QIN
      IF (Q.LE.P .AND. X.LT.0.8D0) GO TO 20
      IF (X.LT.0.2D0) GO TO 20
      Y = 1.0D0 - Y
      P = QIN
      Q = PIN
C
 20   IF ((P+Q)*Y/(P+1.D0).LT.EPS) GO TO 80

So that when the routine gets to the end, the line that Rik pointed out:

      IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI

and the up-front modifications it entails, hopefully do something
beneficial rather than harmful.

Dan