function ibetaImp( x, a, b, regularized, upper, out, stride, offset ) {
	var lambda;
	var prefix;
	var fract;
	var bbar;
	var div;
	var tmp;
	var i0;
	var i1;
	var k;
	var n;
	var p;
	var y;

	y = 1.0 - x;
	i0 = offset;
	i1 = offset + stride;

	// Derivative not set...
	out[ i1 ] = -1;
	if ( isnan( x ) || x < 0.0 || x > 1.0 ) {
		out[ i0 ] = NaN;
		out[ i1 ] = NaN;
		return out;
	}
	if ( regularized ) {
		if ( a < 0.0 || b < 0.0 ) {
			out[ i0 ] = NaN;
			out[ i1 ] = NaN;
			return out;
		}
		// Extend to a few very special cases...
		if ( a === 0.0 ) {
			if ( b === 0.0 ) {
				out[ i0 ] = NaN;
				out[ i1 ] = NaN;
				return out;
			}
			if ( b > 0.0 ) {
				out[ i0 ] = ( upper ) ? 0.0 : 1.0;
				return out;
			}
		} else if ( b === 0.0 ) {
			if ( a > 0.0 ) {
				out[ i0 ] = ( upper ) ? 1.0 : 0.0;
				return out;
			}
		}
	} else if ( a <= 0.0 || b <= 0.0 ) {
		out[ i0 ] = NaN;
		out[ i1 ] = NaN;
		return out;
	}
	if ( x === 0.0 ) {
		if ( a === 1.0 ) {
			out[ i1 ] = 1.0;
		} else {
			out[ i1 ] = ( a < 1.0 ) ? MAX_FLOAT64 / 2.0 : MIN_FLOAT64 * 2.0;
		}
		if ( upper ) {
			out[ i0 ] = ( regularized ) ? 1.0 : beta( a, b );
			return out;
		}
		out[ i0 ] = 0.0;
		return out;
	}
	if ( x === 1.0 ) {
		if ( b === 1.0 ) {
			out[ i1 ] = 1.0;
		} else {
			out[ i1 ] = ( b < 1.0 ) ? MAX_FLOAT64 / 2.0 : MIN_FLOAT64 * 2.0;
		}
		if ( upper ) {
			out[ i0 ] = 0.0;
		} else {
			out[ i0 ] = ( regularized ) ? 1.0 : beta( a, b );
		}
		return out;
	}
	if ( a === 0.5 && b === 0.5 ) {
		out[ i1 ] = ONE_OVER_PI * sqrt( y * x );

		// We have an arcsine distribution:
		p = ( upper ) ? asin( sqrt(y) ) : asin( sqrt(x) );
		p /= HALF_PI;
		if ( !regularized ) {
			p *= PI;
		}
		out[ i0 ] = p;
		return out;
	}
	if ( a === 1.0 ) {
		tmp = b;
		b = a;
		a = tmp;

		tmp = y;
		y = x;
		x = tmp;

		upper = !upper;
	}
	if ( b === 1.0 ) {
		// Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
		if ( a === 1.0 ) {
			out[ i0 ] = ( upper ) ? y : x;
			out[ i1 ] = 1.0;
			return out;
		}
		out[ i1 ] = a * pow( x, a - 1.0 );
		if ( y < 0.5 ) {
			p = ( upper ) ? -expm1( a * log1p(-y) ) : exp( a * log1p(-y) );
		} else {
			p = ( upper ) ? -( pow( x, a ) - 1.0 ) : pow( x, a );
		}
		if ( !regularized ) {
			p /= a;
		}
		out[ i0 ] = p;
		return out;
	}
	if ( min( a, b ) <= 1.0 ) {
		if ( x > 0.5 ) {
			tmp = b;
			b = a;
			a = tmp;

			tmp = y;
			y = x;
			x = tmp;

			upper = !upper;
		}
		if ( max( a, b ) <= 1.0 ) {
			// Both a,b < 1:
			if ( (a >= min( 0.2, b ) ) || ( pow(x, a) <= 0.9 ) ) {
				if ( upper ) {
					fract = -( ( regularized ) ? 1.0 : beta( a, b ) );
					upper = false;
					fract = -ibetaSeries( a, b, x, fract, regularized, out, y );
				} else {
					fract = ibetaSeries( a, b, x, 0, regularized, out, y );
				}
			} else {
				tmp = b;
				b = a;
				a = tmp;

				tmp = y;
				y = x;
				x = tmp;

				upper = !upper;
				if ( y >= 0.3 ) {
					if ( upper ) {
						fract = -( ( regularized ) ? 1.0 : beta( a, b ) );
						upper = false;
						fract = -ibetaSeries( a, b, x, fract, regularized, out, y ); // eslint-disable-line max-len
					} else {
						fract = ibetaSeries( a, b, x, 0, regularized, out, y );
					}
				} else {
					// Sidestep on a, and then use the series representation:
					if ( regularized ) {
						prefix = 1;
					} else {
						prefix = risingFactorialRatio( a + b, a, 20 );
					}
					fract = ibetaAStep( a, b, x, y, 20, regularized, out );
					if ( upper ) {
						fract -= ( ( regularized ) ? 1 : beta( a, b ) );
						upper = false;
						fract = -betaSmallBLargeASeries( a + 20.0, b, x, y, fract, prefix, regularized ); // eslint-disable-line max-len
					} else {
						fract = betaSmallBLargeASeries( a + 20.0, b, x, y, fract, prefix, regularized ); // eslint-disable-line max-len
					}
				}
			}
		} else if ( b <= 1.0 || ( x < 0.1 && ( pow( b * x, a ) <= 0.7 ) ) ) {
			if ( upper ) {
				fract = -( ( regularized ) ? 1 : beta( a, b ) );
				upper = false;
				fract = -ibetaSeries( a, b, x, fract, regularized, out, y );
			} else {
				fract = ibetaSeries( a, b, x, 0.0, regularized, out, y );
			}
		} else {
			tmp = b;
			b = a;
			a = tmp;

			tmp = y;
			y = x;
			x = tmp;
			upper = !upper;

			if ( y >= 0.3 ) {
				if ( upper ) {
					fract = -(( regularized ) ? 1.0 : beta( a, b ));
					upper = false;
					fract = -ibetaSeries( a, b, x, fract, regularized, out, y );
				} else {
					fract = ibetaSeries( a, b, x, 0.0, regularized, out, y );
				}
			}
			else if ( a >= 15.0 ) {
				if ( upper ) {
					fract = -(( regularized ) ? 1.0 : beta( a, b ));
					upper = false;
					fract = -betaSmallBLargeASeries( a, b, x, y, fract, 1.0, regularized ); // eslint-disable-line max-len
				} else {
					fract = betaSmallBLargeASeries( a, b, x, y, 0.0, 1.0, regularized ); // eslint-disable-line max-len
				}
			}
			else {
				if ( regularized ) {
					prefix = 1;
				} else {
					// Sidestep to improve errors:
					prefix = risingFactorialRatio( a + b, a, 20.0 );
				}
				fract = ibetaAStep( a, b, x, y, 20.0, regularized, out );
				if ( upper ) {
					fract -= ( ( regularized ) ? 1.0 : beta( a, b ) );
					upper = false;
					fract = -betaSmallBLargeASeries( a + 20.0, b, x, y, fract, prefix, regularized ); // eslint-disable-line max-len
				} else {
					fract = betaSmallBLargeASeries( a + 20.0, b, x, y, fract, prefix, regularized ); // eslint-disable-line max-len
				}
			}
		}
	} else {
		// Both a,b >= 1:
		if ( a < b ) {
			lambda = a - ( (a + b) * x );
		} else {
			lambda = ( (a + b) * y ) - b;
		}
		if ( lambda < 0.0 ) {
			tmp = b;
			b = a;
			a = tmp;

			tmp = y;
			y = x;
			x = tmp;
			upper = !upper;
		}
		if ( b < 40.0 ) {
			if (
				floor(a) === a &&
				floor(b) === b &&
				a < MAX_INT32 - 100
			) {
				// Relate to the binomial distribution and use a finite sum:
				k = a - 1.0;
				n = b + k;
				fract = binomialCCDF( n, k, x, y );
				if ( !regularized ) {
					fract *= beta( a, b );
				}
			}
			else if ( b * x <= 0.7 ) {
				if ( upper ) {
					fract = -( ( regularized ) ? 1.0 : beta( a, b ) );
					upper = false;
					fract = -ibetaSeries( a, b, x, fract, regularized, out, y );
				} else {
					fract = ibetaSeries( a, b, x, 0.0, regularized, out, y );
				}
			}
			else if ( a > 15.0 ) {
				// Sidestep so we can use the series representation:
				n = floor( b );
				if ( n === b ) {
					n -= 1;
				}
				bbar = b - n;
				if ( regularized ) {
					prefix = 1;
				} else {
					prefix = risingFactorialRatio( a + bbar, bbar, n );
				}
				fract = ibetaAStep( bbar, a, y, x, n, regularized );
				fract = betaSmallBLargeASeries( a, bbar, x, y, fract, 1.0, regularized ); // eslint-disable-line max-len
				fract /= prefix;
			}
			else if ( regularized ) {
				n = floor( b );
				bbar = b - n;
				if ( bbar <= 0 ) {
					n -= 1;
					bbar += 1;
				}
				fract = ibetaAStep( bbar, a, y, x, n, regularized );
				fract += ibetaAStep( a, bbar, x, y, 20.0, regularized );
				if ( upper ) {
					fract -= 1;
				}
				fract = betaSmallBLargeASeries( a + 20.0, bbar, x, y, fract, 1, regularized ); // eslint-disable-line max-len
				if ( upper ) {
					fract = -fract;
					upper = false;
				}
			}
			else {
				fract = ibetaFraction2( a, b, x, y, regularized, out );
			}
		} else {
			fract = ibetaFraction2( a, b, x, y, regularized, out );
		}
	}
	if ( out[ i1 ] < 0.0 ) {
		out[ i1 ] = ibetaPowerTerms( a, b, x, y, true );
	}
	div = y * x;
	if ( out[ i1 ] !== 0.0 ) {
		if ( ( MAX_FLOAT64 * div < out[ i1 ] ) ) {
			// Overflow, return an arbitrarily large value:
			out[ i1 ] = MAX_FLOAT64 / 2.0;
		} else {
			out[ i1 ] /= div;
		}
	}
	out[ i0 ] = ( upper ) ? ( ( regularized ) ? 1.0 : beta( a, b ) ) - fract : fract; // eslint-disable-line max-len
	return out;
}
Evaluates the incomplete beta function and its first derivative and assigns results to a provided output array.

## Notes

-   This function divides up the input range and selects the right implementation method for each domain.

@param {Probability} x - function input
@param {NonNegativeNumber} a - function parameter
@param {NonNegativeNumber} b - function parameter
@param {boolean} regularized - boolean indicating if the function should evaluate the regularized boolean beta function
@param {boolean} upper - boolean indicating if the function should return the upper tail of the incomplete beta function instead
@param {(Array|TypedArray|Object)} out - output array
@param {integer} stride - output array stride
@param {NonNegativeInteger} offset - output array index offset
@returns {(Array|TypedArray|Object)} function value and first derivative

@example
var out = ibetaImp( 0.5, 2.0, 2.0, false, false, [ 0.0, 0.0 ], 1, 0 );
// returns [ ~0.083, ~1.5 ]

@example
var out = ibetaImp( 0.2, 1.0, 2.0, false, true, [ 0.0, 0.0 ], 1, 0 );
// returns [ 0.32, 1.6 ]

@example
var out = ibetaImp( 0.2, 1.0, 2.0, true, true, [ 0.0, 0.0 ], 1, 0 );
// returns [ 0.64, 1.6 ]