r2rs-nmath 0.1.1

// Translation of nmath's dgamma
//
// Copyright (C) 2020 neek-sss
//
// 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 3 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, see <https://www.gnu.org/licenses/>.

use num_traits::Float;
use strafe_type::{FloatConstraint, Positive64, Rational64, Real64};

use crate::{distribution::pois::dpois_raw, traits::DPQ};

/// Computes the density of the gamma distribution,
///
///
/// $ p(x;a,s) = \frac{1/s (x/s)^{a-1} exp(-x/s)}{(a-1)!} $
///
/// where 's' is the scale (= 1/lambda in other parametrizations)
/// and 'a' is the shape parameter ( = alpha in other contexts).
///
/// The old (R 1.1.1) version of the code is available via '#define D_non_pois'
pub fn dgamma<RE: Into<Real64>, P: Into<Positive64>, RA: Into<Rational64>>(
    x: RE,
    shape: P,
    scale: RA,
    log: bool,
) -> Real64 {
    let x = x.into().unwrap();
    let shape = shape.into().unwrap();
    let scale = scale.into().unwrap();

    let mut pr = 0.0;
    if x < 0.0 {
        return f64::d_0(log).into();
    }
    if shape == 0.0 {
        /* point mass at 0 */
        let ret = if x == 0.0 {
            f64::infinity()
        } else {
            f64::d_0(log)
        };
        return ret.into();
    }
    if x == 0.0 {
        if shape < 1.0 {
            return f64::infinity().into();
        }
        if shape > 1.0 {
            return f64::d_0(log).into();
        }
        /* else */
        return if log { -scale.ln() } else { (1.0) / scale }.into();
    }

    if shape < 1.0 {
        pr = dpois_raw(shape, x / scale, log);
        return if log {
            /* NB: currently *always*  shape/x > 0  if shape < 1:
             * -- overflow to Inf happens, but underflow to 0 does NOT : */
            pr + if (shape / x).is_finite() {
                (shape / x).ln()
            } else {
                /* shape/x overflows to +Inf */
                shape.ln() - x.ln()
            }
        } else {
            pr * shape / x
        }
        .into();
    }
    /* else  shape >= 1 */
    pr = dpois_raw(shape - 1.0, x / scale, log);
    if log { pr - scale.ln() } else { pr / scale }.into()
}