r2rs-nmath 0.1.1

// Translation of nmath's dweibull
//
// 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, Rational64, Real64};

use crate::traits::DPQ;

/// The density function of the Weibull distribution.
pub fn dweibull<RE: Into<Real64>, RA1: Into<Rational64>, RA2: Into<Rational64>>(
    x: RE,
    shape: RA1,
    scale: RA2,
    log: bool,
) -> Real64 {
    let x = x.into().unwrap();
    let shape = shape.into().unwrap();
    let scale = scale.into().unwrap();

    let mut tmp1 = 0.0;
    let mut tmp2 = 0.0;

    if x < 0.0 {
        return f64::d_0(log).into();
    }
    if !x.is_finite() {
        return f64::d_0(log).into();
    }
    /* need to handle x == 0 separately */
    if x == 0.0 && shape < 1.0 {
        return f64::infinity().into();
    }
    tmp1 = (x / scale).powf(shape - 1.0);
    tmp2 = tmp1 * (x / scale);
    /* These are incorrect if tmp1 == 0 */
    if log {
        let x = -tmp2;
        let y = (shape * tmp1 / scale).ln();

        // Hopefully fixes some NaN's (bug in original)
        if x.is_infinite() && y.is_infinite() && x.signum() != y.signum() {
            0.0
        } else {
            x + y
        }
    } else if (-tmp2).exp() == 0.0 {
        0.0
    } else {
        (shape * tmp1 * (-tmp2).exp()) / scale
    }
    .into()
}