longfloat library HP PRIME

[robmio]

This is the algorithm:

Code:


#cas
CDF_Roy(s,m,n,theta):=
BEGIN
LOCAL A, ii, j, b, a, adzero, aa, cc;
A:=MAKEMAT(0,s,s);
b:=MAKELIST(0,x,1,s);
cc(x):=sqrt(π)*Gamma((2*m+2*n+s+x+2)/
2)/(Gamma((2*m+x+1)/2)*Gamma((2*n+x+
1)/2)*Gamma(x/2));
FOR ii FROM 1 TO s DO
b:=REPLACE(b,ii,{(Beta(m+ii,n+1,theta)
^2)/2});
FOR j FROM ii TO s-1 DO
b:=REPLACE(b,j+1,{(m+j)/(m+j+n+1)*
b(j)-Beta(2*m+ii+j,2*n+2,theta)/
(m+j+n+1)});
a:=(Beta(m+ii,n+1,theta)*Beta(m+j+1,
n+1,theta)-2*b(j+1))*cc(ii)*
cc(j+1);
A:=REPLACE(A,{ii,j+1},[[a]]);
END;
END;
aa:={};
FOR ii FROM 1 TO s DO
aa:=CONCAT(aa,{Beta(m+ii,n+1,theta)*
cc(ii)});
END;
aa:=ListToMat(aa);
adzero:=MAKELIST(0,x,1,s+1,1);
adzero:=ListToMat(adzero);
IF odd(s)==1 THEN
A:=ADDCOL(A,aa,s+1);
A:=ADDROW(A,adzero,s+1);
END;
A:=A-TRN(A);
RETURN √(DET(A));
END;
#end

[/quote]