public static double powerSeries(double a,double b,double x){
  double s, t, u, v, n, t1, z, ai;
  ai=1.0 / a;
  u=(1.0 - b) * x;
  v=u / (a + 1.0);
  t1=v;
  t=u;
  n=2.0;
  s=0.0;
  z=MACHEP * ai;
  while (Math.abs(v) > z) {
    u=(n - b) * x / n;
    t*=u;
    v=t / (a + n);
    s+=v;
    n+=1.0;
  }
  s+=t1;
  s+=ai;
  u=a * Math.log(x);
  if ((a + b) < MAXGAM && Math.abs(u) < MAXLOG) {
    t=gamma(a + b) / (gamma(a) * gamma(b));
    s=s * t * Math.pow(x,a);
  }
 else {
    t=lnGamma(a + b) - lnGamma(a) - lnGamma(b) + u + Math.log(s);
    if (t < MINLOG) {
      s=0.0;
    }
 else {
      s=Math.exp(t);
    }
  }
  return s;
}

public static double powerSeries ( double a , double b , double x ) { double s , t , u , v , n , t1 , z , ai ; ai = 1.0 / a ; u = ( 1.0 - b ) * x ; v = u / ( a + 1.0 ) ; t1 = v ; t = u ; n = 2.0 ; s = 0.0 ; z = MACHEP * ai ; while ( Math . abs ( v ) > z ) { u = ( n - b ) * x / n ; t *= u ; v = t / ( a + n ) ; s += v ; n += 1.0 ; } s += t1 ; s += ai ; u = a * Math . log ( x ) ; if ( ( a + b ) < MAXGAM && Math . abs ( u ) < MAXLOG ) { t = gamma ( a + b ) / ( gamma ( a ) * gamma ( b ) ) ; s = s * t * Math . pow ( x , a ) ; } else { t = lnGamma ( a + b ) - lnGamma ( a ) - lnGamma ( b ) + u + Math . log ( s ) ; if ( t < MINLOG ) { s = 0.0 ; } else { s = Math . exp ( t ) ; } } return s ; }
power series for incomplete beta integral . use when b * x is small and x not too close to 1 .
