***********************************************************************;
* Functions to compute the confluent hypergeometric function, M(a,b,x);
* Methods are described in detail in 
* "Computing the Confluent Hypergeometric Function, M(a,b,z)",
* by K. E. Muller, 1998 (manuscript in review).;
***********************************************************************;
* For j in {1,2,3,4,5,1C,2C,4C,5C}
* function _MABXj_(A,B,X,MAXTERMS) returns   sum(a,b,x)  //
*                                            LFACTOR     //
*                                            NTERMS      //
*                                            CONVERGED   //
*                                            max(abs(term(i))) ;
* Note that M(a,b,x)=sum(a,b,x)*exp(LFACTOR) ;
* NTERMS= # of terms at which series terminated;
* CONVERGED=1 if program judges computation converged, 0 otherwise;
* Failure to converge to 1E-10 relative accuracy gives a missing value;
***********************************************************************;
* All four inputs are scalar and real numeric values;
* MAXTERMS, controls the maximum number of terms
* evaluated before the function quits and declares non-convergence;
***********************************************************************;


START _MABX1_(A,B,X,MAXTERMS);
*Method 1.  Use defining series;
*WARNING*Fails obviously or secretly for large |xa/b|;

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

A_1=A-1;
B_1=B-1;
TERMI= X# (A / B);
MABX=1+TERMI;
HOLD={1}//MABX;
NTERMS=1;
TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;
  DO I=2 TO MAXTERMS;
  NTERMS=NTERMS+1;
  IF ABS(TERMI/MABX)<=TOLERANC THEN 
                      RETURN(MABX//{0}//NTERMS//{1}//MAX(ABS(HOLD)));
  IF ABS(TERMI)>1E304 THEN GO TO DIVERGE;
  TERMI=TERMI # (A_1+I) # (X / ( (B_1+I)#I ) );
  MABX=MABX+TERMI;
  HOLD=HOLD//MABX;
  END;
DIVERGE: RETURN( {.}//{0}//NTERMS//{0}//MAX(ABS(HOLD)) ); *Failed;
FINISH _MABX1_;


START _MABX1C_(A,B,X,MAXTERMS);
*Use defining series and continued fraction with forward recurrence;
**WARNING**This may fail obviously or secretly for large |xa/b|,

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;
BIGA={1,1}//(1+A#(X/B)); *=BIGA[-1]//BIGA[0]//BIGA[1];
  DO I=2 TO MAXTERMS;
  R_I=(A+I-1)/(I#(B+I-1));
  IF R_I=0 THEN RETURN(
     BIGA[I+1]//{0}//(NROW(BIGA)-1)//{1}//MAX(ABS(BIGA)) );*converged;
  IF ABS(BIGA[I+1])>1E304 THEN GO TO DIVERGE;
  BIGA_I1=BIGA[I+1]+R_I#X#(BIGA[I+1]-BIGA[I]);
  BIGA=BIGA//BIGA_I1;
  IF ABS((BIGA_I1-BIGA[I+1])/BIGA_I1)<=TOLERANC THEN RETURN(
     BIGA_I1//{0}//(NROW(BIGA)-1)//{1}//MAX(ABS(BIGA)));*converged;
  END;
DIVERGE: RETURN({.,.}//(NROW(BIGA)-1)//{0}//MAX(ABS(BIGA))); 
FINISH _MABX1C_;


START _MABX2_(A,B,X,MAXTERMS);
*Use asymptotic approximation;
*See eqn's 16.2.5-6 p466, W. J. Thompson, 1997, which are ;
*Abramowitz and Stegun, p508, eqn 13.5.1, reduced for X real;
*Requires b>0 a>0 and |x|>>0;

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

*Requires b>a>0 and |x|>>0;
IF (A<=0)|(B<=A) THEN RETURN(J(5,1,.));
* As detailed in Muller(1998), abs(a(b-a)/x)<1 a better test;

TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;
TERMI=1;
SUM=1;
NTERMS=1;
HOLD=SUM;
IF X>0 THEN LFACTOR=LGAMMA(B) - LGAMMA(A)   - LOG(X) #(B-A) + X;
       ELSE LFACTOR=LGAMMA(B) - LGAMMA(B-A) - LOG(-X)#A;

IF X<0 THEN GO TO NEGATIVE;

****************************;
*This block of code for X>>0;
B_A_1 =B-A-1;
_1_A_1=1-A-1;
  DO I=1 TO MAXTERMS;
  LAST=TERMI;
  TERMI=TERMI # (B_A_1+I) # (_1_A_1+I) / (I#X) ;
  NTERMS=NTERMS+1;
  SUM=SUM+TERMI;
  HOLD=HOLD//SUM;
  IF ABS(TERMI/SUM)<=TOLERANC 
    THEN RETURN( SUM//LFACTOR//NTERMS//{1}//MAX(ABS(HOLD)));*Converged;
  IF (ABS(TERMI)>ABS(LAST))&(I>2) THEN
     RETURN( {.} //LFACTOR//NTERMS//{0}//MAX(ABS(HOLD))); *Diverged;
  END;
RETURN( {.} // LFACTOR //NTERMS//{0}//MAX(ABS(HOLD))); *Failed;

****************************;
*This block of code for X<<0;
NEGATIVE:;
_X=-X;
A_B =A-B;
A_1=A-1;
  DO I=1 TO MAXTERMS;
  LAST=TERMI;
  TERMI=TERMI # (A_B+I) # (A_1+I) / (I#_X) ;
  SUM=SUM+TERMI;
  HOLD=HOLD//SUM;
  NTERMS=NTERMS+1;
  IF ABS(TERMI/SUM)<=TOLERANC 
     THEN RETURN( SUM//LFACTOR//I//{1}//MAX(ABS(HOLD))) ; *Converged;
  IF (ABS(TERMI)>ABS(LAST))&(I>2) THEN
     RETURN( {.} //LFACTOR//NTERMS//{0}//MAX(ABS(HOLD))); *Diverged;
  END;
RETURN( {.} // LFACTOR //NTERMS//{0}//MAX(ABS(HOLD))); *Failed;

FINISH _MABX2_;


START _MABX2C_(A,B,X,MAXTERMS);
*Asymptotic approximation with continued fraction, forward recurrence;
*See eqn's 16.2.5-6 p466, W. J. Thompson, 1997, which are ;
*Abramowitz and Stegun, p508, eqn 13.5.1, reduced for X real;
*Requires b>a>0 and |X|>>0;

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

*Requires b>a>0 and |x|>>0;
IF (A<=0)|(B<=A) THEN RETURN(J(5,1,.));
* As detailed in Muller(1998), abs(a(b-a)/x)<1 a better test;

TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;
IF X>0 THEN LFACTOR=LGAMMA(B) - LGAMMA(A)   - LOG(X) #(B-A) + X ;
       ELSE LFACTOR=LGAMMA(B) - LGAMMA(B-A) - LOG(-X)#A;

IF X<0 THEN GO TO NEGATIVE;

****************************;
*This block of code for X>>0;
B_A=B-A;
_1_A=1-A;
BIGA=  {1}   //  {1}  //(1+(B_A)#(1-A)/X);
*** =BIGA[-1]//BIGA[0]//BIGA[1];
  DO I=2 TO MAXTERMS;
  R_I=(B_A+I-1)#(_1_A+I-1)/I;
  IF R_I=0 THEN RETURN(
     BIGA[I+1]//LFACTOR//NROW(BIGA)//{1}//MAX(ABS(BIGA)));*Converged;
  IF (ABS(R_I/X)>1) & (I>2) THEN RETURN(
      {.}//LFACTOR//NROW(BIGA)//{0}//MAX(ABS(BIGA))); *Diverged;
  BIGA_I1=BIGA[I+1]+(R_I/X)#(BIGA[I+1]-BIGA[I]);
  BIGA=BIGA//BIGA_I1;
  IF ABS((BIGA_I1-BIGA[I+1])/BIGA_I1)<=TOLERANC THEN RETURN(
     BIGA_I1//LFACTOR//(NROW(BIGA)-1)//{1}//MAX(ABS(BIGA)));*Converged;
  END;
RETURN({.}//LFACTOR//(NROW(BIGA)-1)//{0}//MAX(ABS(BIGA))); *failed;

****************************;
*This block of code for X<<0;
NEGATIVE:;
_X=-X;
A_B=A-B;
A_1=A-1;
BIGA=  {1}   // {1}   //(1+(A_B+1)#A/(_X));
*** =BIGA[-1]//BIGA[0]//BIGA[1];
  DO I=2 TO MAXTERMS;
  R_I=(A_B+I)#(A_1+I)/I;
  IF R_I=0 THEN RETURN(
      BIGA[I+1]//LFACTOR//NROW(BIGA)//{1}//MAX(ABS(BIGA))); *Converged;
  IF (ABS(R_I/_X)>1) & (I>2) THEN RETURN( 
     {.}//LFACTOR//NROW(BIGA)//{0}//MAX(ABS(BIGA)) ); *Diverged;
  BIGA_I1=BIGA[I+1]+(R_I/(_X))#(BIGA[I+1]-BIGA[I]);
  BIGA=BIGA//BIGA_I1;
  IF ABS((BIGA_I1-BIGA[I+1])/BIGA_I1)<=TOLERANC THEN RETURN(
     BIGA_I1//LFACTOR//(NROW(BIGA)-1)//{1}//MAX(ABS(BIGA))); *Converged;
  END;
RETURN({.}//LFACTOR//(NROW(BIGA)-1)//{0}//MAX(ABS(BIGA))); *Failed;

FINISH _MABX2C_;


START _MABX3_(A,B,X,MAXTERMS);
*Rational approximation in Y. L. Luke (1977), p182-183:
*"Algorithms for the Computation of Mathematical Functions,"
*New York, Academic press;

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

XA=ABS(X);
IF X>0 THEN D=B-A;
       ELSE D=A;
IF X>0 THEN LFACTOR=X;
       ELSE LFACTOR=0;

BIGQ0={1};
BIGQ1= 1 + (D+1)#XA/(2#B) ;
BIGQ2= 1 + (D+2)#XA/(2#(B+1)) + (D+1)#(D+2)#XA#XA/(12#B#(B+1)) ;
BIGQ=BIGQ0//BIGQ1//BIGQ2;

BIGP0={1};
BIGP1=BIGQ1-D#XA/B;
BIGP2=BIGQ2-(D#XA/B)#(1 +(D+2)#XA/(2#(B+1))) + D#(D+1)#XA#XA/(2#B#(B+1));
BIGP=BIGP0//BIGP1//BIGP2;

HOLD=BIGP/BIGQ;
TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;
  DO I=3 TO MAXTERMS;
  F1=(I-D-2)/(2#(2#I-3)#(I+B-1));
  F2=(I+D)#(I+D-1) / (4#(2#I-1)#(2#I-3)#(I+B-2)#(I+B-1)) ;
  F3=-(I+D-2)#(I+D-1)#(I-D-2) / (8#((2#I-3)##2)#(2#I-5)#
                                   (I+B-3)#(I+B-2)#(I+B-1) ) ;
  F4=-(I+D-1)#(I-B-1) / ( 2#(2#I-3)#(I+B-2)#(I+B-1) ) ;
  IF MAX(ABS(BIGP[I]//BIGQ[I]))>1E304 THEN GO TO DIVERGE;
  BIGP_I= (1+F1#XA)#BIGP[I] + (F4+F2#XA)#XA#BIGP[I-1] +
                                F3#XA#XA#XA#BIGP[I-2];
  BIGQ_I= (1+F1#XA)#BIGQ[I] + (F4+F2#XA)#XA#BIGQ[I-1] +
                                F3#XA#XA#XA#BIGQ[I-2];
  HOLD=HOLD//(BIGP_I/BIGQ_I);
  BIGP=BIGP//BIGP_I;
  BIGQ=BIGQ//BIGQ_I;
  IF ABS((HOLD[I+1]-HOLD[I])/HOLD[I+1])<=TOLERANC THEN
           RETURN(HOLD[I+1]//LFACTOR//NROW(BIGP)//{1}//
                  MAX(ABS(BIGP//BIGQ)) );
  END;
DIVERGE: RETURN({.}//LFACTOR//NROW(BIGP)//{0}//MAX(ABS(BIGP//BIGQ)));
FINISH _MABX3_;


START _MABX4_(A,B,X,MAXTERMS);
*Series approximation due to Muller (1998);
*Requires b>a>0;
*Works best with abs(a(b-a)/x)<<1;

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

*Requires b>a>0;
IF (A<=0)|(B<=A) THEN RETURN(J(5,1,.));

TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;

IF X>0 THEN LFACTOR=LGAMMA(B) - LGAMMA(A)   - LOG(X) #(B-A) + X;
       ELSE LFACTOR=LGAMMA(B) - LGAMMA(B-A) - LOG(-X)#A;

IF X<0 THEN GO TO NEGATIVE;

****************************;
*This block of code for X>0;
FACTORI=(B-A)#(1-A) / X ;
TERMI=PROBGAM(X,B-A+1) # FACTORI;
SUM=PROBGAM(X,B-A)+TERMI;
HOLD=PROBGAM(X,B-A)//SUM;
B_A_1 =B-A-1;
_1_A_1=1-A-1;
  DO I=2 TO MAXTERMS;
  IF ABS(TERMI/SUM)<=TOLERANC THEN 
     RETURN(SUM // LFACTOR//I// {1}//MAX(ABS(HOLD))); *Converged;
  LAST=TERMI;
  IF ABS(FACTORI)>1E304 THEN GO TO DIVERGE1;
  FACTORI=FACTORI # (B_A_1+I) # (_1_A_1+I) / (I#X) ;
  TERMI=PROBGAM(X,I+B-A) #FACTORI;
  SUM=SUM+TERMI;
  HOLD=HOLD//SUM;
  END;
DIVERGE1: RETURN( {.}//LFACTOR//I//{0}//MAX(ABS(HOLD)) );

****************************;
*This block of code for X<0;
NEGATIVE:;
_X=-X;
FACTORI= (1+A-B)#A /_X ;
TERMI= PROBGAM(_X,A+1)#FACTORI ;
SUM=PROBGAM(_X,A) + TERMI;
HOLD=PROBGAM(_X,A)//SUM;
A_B =1+A-B-1;
A_1=A-1;
  DO I=2 TO MAXTERMS;
  IF ABS(TERMI/SUM)<=TOLERANC THEN 
     RETURN(SUM // LFACTOR//I// {1}//MAX(ABS(HOLD))); *Converged;
  LAST=TERMI;
  IF ABS(FACTORI)>1E304 THEN GO TO DIVERGE2;
  FACTORI=FACTORI#(A_B+I) # (A_1+I) / (I#_X);
  TERMI= PROBGAM(_X,A+I) #  FACTORI;
  SUM=SUM+TERMI;
  HOLD=HOLD//SUM;
  END;
DIVERGE2: RETURN( {.}//LFACTOR//I//{0}//MAX(ABS(HOLD)) ); 

FINISH _MABX4_;


START _MABX4C_(A,B,X,MAXTERMS);
*Series approximation due to Muller (1998);
*Use continued fraction, forward recurrence ;
*Requires b>a>0;
*Works best if abs(a(b-a)/x)<<1;

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

*Requires b>a>0;
IF (A<=0)|(B<=A) THEN RETURN(J(5,1,.));

TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;
IF X>0 THEN LFACTOR=LGAMMA(B) - LGAMMA(A)   - LOG(X) #(B-A) + X ;
       ELSE LFACTOR=LGAMMA(B) - LGAMMA(B-A) - LOG(-X)#A;

IF X<0 THEN GO TO NEGATIVE;

****************************;
*This block of code for X>0;
B_A=B-A;
_1_A=1-A;
PROB=PROBGAM(X,B_A)//PROBGAM(X,B_A+1);
BIGA=  {1}   //PROB[1]//(PROB[1]+PROB[2]#(B_A)#(1-A)/X);
*** =BIGA[-1]//BIGA[0]//BIGA[1];
  DO I=2 TO MAXTERMS;
  PROB=PROB//PROBGAM(X,B_A+I);
  R_I1=(B_A+I-1)#(_1_A+I-1)/I;
  IF PROB[I]=0 THEN RETURN({.}//
       LFACTOR//(NROW(BIGA)-1)//{0}//MAX(ABS(BIGA))); *Did not converge;
  R_I2=PROB[I+1]/PROB[I];
  R_I=R_I1#R_I2;
  IF R_I=0 THEN RETURN(
     BIGA[I+1]//LFACTOR//NROW(BIGA)//{1}//MAX(ABS(BIGA))); *Converged;
  BIGA_I1=BIGA[I+1]+(R_I/X)#(BIGA[I+1]-BIGA[I]);
  BIGA=BIGA//BIGA_I1;
  IF ABS((BIGA_I1-BIGA[I+1])/BIGA_I1)<=TOLERANC THEN RETURN(
     BIGA_I1//LFACTOR//(NROW(BIGA)-1)//{1}//MAX(ABS(BIGA))); *Converged;
  END;
RETURN({.}//LFACTOR//(NROW(BIGA)-1)//{0}//MAX(ABS(BIGA))); *Failed;

***************************;
*This block of code for X<0;
NEGATIVE:;
_X=-X;
A_B=A-B;
A_1=A-1;
PROB=PROBGAM(_X,A)//PROBGAM(_X,A+1);
BIGA=  {1}   //PROB[1]//(PROB[1]+PROB[2]#(A_B+1)#A/_X);
*** =BIGA[-1]//BIGA[0]//BIGA[1];
  DO I=2 TO MAXTERMS;
  PROB=PROB//PROBGAM(_X,A+I);
  R_I1=(A_B+I)#(A_1+I)/I;
  IF PROB[I]=0 THEN RETURN({.}//
     LFACTOR//(NROW(BIGA)-1)//{0}//MAX(ABS(BIGA)) );*Did not converge;
  R_I2=PROB[I+1]/PROB[I];
  R_I=R_I1#R_I2;
  IF R_I=0 THEN RETURN(
    BIGA[I+1]//LFACTOR//NROW(BIGA)//{1}//MAX(ABS(BIGA))); *Converged;
  BIGA_I1=BIGA[I+1]+(R_I/(-X))#(BIGA[I+1]-BIGA[I]);
  BIGA=BIGA//BIGA_I1;
  IF ABS((BIGA_I1-BIGA[I+1])/BIGA_I1)<=TOLERANC THEN RETURN(
     BIGA_I1//LFACTOR//(NROW(BIGA)-1)//{1}//MAX(ABS(BIGA)));*Converged;
  END;
RETURN({.}//LFACTOR//(NROW(BIGA)-1)//{0}//MAX(ABS(BIGA))); *Failed;

FINISH _MABX4C_;


START _MABX5_(A,B,X,MAXTERMS);
*Beta moment approximation due to Muller (1998);
*Requires b>a>0;

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

*Requires b>a>0;
IF (A<=0)|(B<=A) THEN RETURN(J(5,1,.));

TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;
LFACTOR= A#X/B ;
B_A=B-A;
MU2= A#B_A/ ( B#B#(B+1) );
MU3= 2#A#B_A#(B-2#A) /
     ( (B##3) # (B+1) # (B+2) );
MU={1,0}//MU2//MU3;
MEAN=A/B;
LOG_X=LOG(ABS(X));
SUM=1+MU2#X#(X/2)+MU3#X#X#(X/6);
HOLD={1,1}//(1+MU2#X#(X/2))//(1+MU2#X#(X/2)+MU3#X#X#(X/6));
  DO M_1=4 TO MAXTERMS BY 1;
  SUM_J=0;
  M=M_1-1;
    DO J=1 TO M;
    KSET=(DO(0,J-1,1))`;
    KPROD=B_A#((M-KSET)/(B#(B+M-KSET)));
    PRODUCT=KPROD[#];
    SUM_J=SUM_J+PRODUCT#MU[M_1-J]/(B+M-J);
    END;
  MU_M1=A#SUM_J - MEAN#MU[M_1]#M/(B+M);
  MU=MU//MU_M1;

  LFACTRM1=M_1#LOG_X-LGAMMA(M_1+1);
  IF MU_M1=0 THEN TERM_M1=0;
             ELSE DO;
                  LOGM1=LOG(ABS(MU_M1))+LFACTRM1;
                  IF LOGM1>706 THEN GO TO DIVERGE;
                  TERM_M1=SIGN(MU_M1)#EXP(LOGM1);
                  END;
  IF (X<0)&(MOD(M_1,2)=1) THEN TERM_M1=-TERM_M1;
  SUM=SUM+TERM_M1;
  HOLD=HOLD//SUM;

  IF ABS(TERM_M1/SUM)<=TOLERANC THEN GO TO CONVERGE;
  END;
DIVERGE: RETURN({.}//LFACTOR//NROW(MU)//{0}//MAX(ABS(HOLD))); 

CONVERGE:  RETURN(SUM//LFACTOR//NROW(MU)//{1}//MAX(ABS(HOLD)));
FINISH _MABX5_;


START _MABX5C_(A,B,X,MAXTERMS);
*Beta moment approximation due to Muller (1998);
*Use continued fraction with forward recurrence;
*Requires b>a>0;

*Treat general special cases;
IF (ROUND(B,1)=B)&(B<=0) THEN RETURN(J(5,1,.));
IF (X=0)|(A=0) THEN RETURN({1,0,1,1,1});
IF A=B THEN RETURN(EXP(X)//{0,1,1,1});

*Requires b>a>0;
IF (A<=0)|(B<=A) THEN RETURN(J(5,1,.));

TOLERANC=1E-10; *numeric zero, relative inaccuracy tolerated;
LFACTOR=A#X/B;
B_A=B-A;
MU2= A#B_A/ ( B#B#(B+1) );
MU3= 2#A#B_A#(B-2#A) /  ( (B##3) # (B+1) # (B+2) );
MU={1,0}//MU2//MU3;
MEAN=A/B;

M_I=1;
FACTORI=1;
BIGA={1}//(MU2/2)//(MU2/2+MU3#X/6);
BIGA_I1=BIGA[3];
  DO I=2 TO MAXTERMS BY 1;
  M=I+1;
  M_1=M+1;

  SUM_J=0;  *compute mu[m];
    DO J=1 TO M;
    KSET=(DO(0,J-1,1))`;
    KPROD=B_A#((M-KSET)/(B#(B+M-KSET)));
    PRODUCT=KPROD[#];
    SUM_J=SUM_J+PRODUCT#MU[M_1-J]/(B+M-J);
    END;
  MU_M1=A#SUM_J - MEAN#MU[M_1]#M/(B+M);
  MU=MU//MU_M1;

  R_I1=X/(I+2);
  IF MU[I+2]=0 THEN GO TO DIVERGE;
  R_I2=MU[I+3]/MU[I+2];
  R_I=R_I1#R_I2;

  IF ABS(BIGA[I+1])>1E304 THEN GO TO DIVERGE;
  BIGA_I1=BIGA[I+1]+R_I#(BIGA[I+1]-BIGA[I]);
  BIGA=BIGA//BIGA_I1;
  IF ABS((BIGA_I1-BIGA[I+1])/BIGA[I+1])<=TOLERANC THEN GO TO CONVERGE;
  END;
DIVERGE: RETURN({.}//LFACTOR//NROW(MU)//{0}//MAX(ABS(BIGA))); *Failed;

CONVERGE: RETURN(
         (1+X#BIGA_I1#X)//LFACTOR//NROW(MU)//{1}//MAX(ABS(BIGA)));
FINISH _MABX5C_;