/*
 *  Mathlib : A C Library of Special Functions
 *  Copyright (C) 1998 Ross Ihaka
 *  Copyright (C) 2000-12 The R Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.R-project.org/Licenses/
 *
 *  SYNOPSIS
 *
 *    #include <Rmath.h>
 *    double dnbeta(double x, double a, double b, double ncp, int give_log);
 *
 *  DESCRIPTION
 *
 *    Computes the density of the noncentral beta distribution with
 *    noncentrality parameter ncp.  The noncentral beta distribution
 *    has density:
 *
 *		       Inf
 *	f(x|a,b,ncp) = SUM  p(i) *  x^(a+i-1) * (1-x)^(b-1) / B(a+i,b)
 *		       i=0
 *
 *    where:
 *
 *		p(k) = exp(-ncp/2) (ncp/2)^k / k!
 *
 *	      B(a,b) = Gamma(a) * Gamma(b) / Gamma(a+b)
 *
 *
 *    This can be computed efficiently by using the recursions:
 *
 *	      p(k+1) = ncp/2 / (k+1) * p(k)
 *
 *      B(a+k+1,b)   = (a+k)/(a+b+k) * B(a+k,b)
 *
 * The new algorithm first determines for which k the k-th term is maximal,
 * and then sums outwards to both sides from the 'mid'.
 */

#include "nmath.h"
#include "dpq.h"

double dnbeta(double x, double a, double b, double ncp, int give_log)
{
    const static double eps = 1.e-15;

    int kMax;
    double k, ncp2, dx2, d, D;
    LDOUBLE sum, term, p_k, q;

#ifdef IEEE_754
    if (ISNAN(x) || ISNAN(a) || ISNAN(b) || ISNAN(ncp))
	return x + a + b + ncp;
#endif
    if (ncp < 0 || a <= 0 || b <= 0)
	ML_ERR_return_NAN;

    if (!R_FINITE(a) || !R_FINITE(b) || !R_FINITE(ncp))
	ML_ERR_return_NAN;

    if (x < 0 || x > 1) return(R_D__0);
    if(ncp == 0)
	return dbeta(x, a, b, give_log);

    /* New algorithm, starting with *largest* term : */
    ncp2 = 0.5 * ncp;
    dx2 = ncp2*x;
    d = (dx2 - a - 1)/2;
    D = d*d + dx2 * (a + b) - a;
    if(D <= 0) {
	kMax = 0;
    } else {
	D = ceil(d + sqrt(D));
	kMax = (D > 0) ? (int)D : 0;
    }

    /* The starting "middle term" --- first look at it's log scale: */
    term = dbeta(x, a + kMax, b, /* log = */ TRUE);
    p_k = dpois_raw(kMax, ncp2,              TRUE);
    if(x == 0. || !R_FINITE(term) || !R_FINITE((double)p_k)) /* if term = +Inf */
	return R_D_exp((double)(p_k + term));

    /* Now if s_k := p_k * t_k  {here = exp(p_k + term)} would underflow,
     * we should rather scale everything and re-scale at the end:*/

    p_k += term; /* = log(p_k) + log(t_k) == log(s_k) -- used at end to rescale */
    /* mid = 1 = the rescaled value, instead of  mid = exp(p_k); */

    /* Now sum from the inside out */
    sum = term = 1. /* = mid term */;
    /* middle to the left */
    k = kMax;
    while(k > 0 && term > sum * eps) {
	k--;
	q = /* 1 / r_k = */ (k+1)*(k+a) / (k+a+b) / dx2;
	term *= q;
	sum += term;
    }
    /* middle to the right */
    term = 1.;
    k = kMax;
    do {
	q = /* r_{old k} = */ dx2 * (k+a+b) / (k+a) / (k+1);
	k++;
	term *= q;
	sum += term;
    } while (term > sum * eps);

#ifdef HAVE_LONG_DOUBLE
    return R_D_exp((double)(p_k + logl(sum)));
#else
    return R_D_exp((double)(p_k + log(sum)));
#endif
}