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;