用逼近法计算不完全贝塔函数的值

function bp = betap（xab）
%自变量的值：x
%第一个参数：a
%第二个参数：b
%自变量取x值时的不完全贝塔函数值：gp

format long;
if  a <=0 || b<=0 
    disp（‘a和b不能小于0!‘）;
    return;
end

if  x == 0 || x == 1
    bp = 0;
end
    
c1 = power（xa）;
c2 = power（1-xb）;
c3 = Beta（ab）;

if x<（（a+1）/（a+b+2））
    n = 1;
    while 1
        f1 = fi（2*nxab）;
        f2 = fi（2*n+2xab）;
        if abs（f2-f1）<1.0e-30
            bp  = f2*c1*c2/a/c3;
            break;
        else
            n = n + 1;
        end
    end       
else
    if x ==0.5 && a == b
        bp = 0.5;