Code shared from the Rust Playground

playground.rs

/// Calculates the factorial of a non-negative integer n.
/// Since the input z is the upper limit of the sum, and k goes up to z,
/// we can assume k is small enough that the factorial fits in u64.
/// However, for the division in the formula to work with f64,
/// we'll return the result as an f64.
fn factorial(n: u32) -> f64 {
    if n == 0 {
        1.0
    } else {
        (1..=n).map(|i| i as f64).product()
    }
}

/// Implements the given mathematical expression.
/// 1 - \sum_{k=0}^{z} \frac{\lambda^k e^{-\lambda}}{k!} \left(1 - \left(\frac{q}{p}\right)^{(z-k)}\right)
///
/// Parameters:
/// - lambda: The parameter \lambda (a float, must be non-negative).
/// - p: The parameter p (a float, denominator in the base of the power, must be non-zero).
/// - q: The parameter q (a float, numerator in the base of the power).
/// - z: The upper limit of the sum (an unsigned integer).
///
/// Returns:
/// The result of the expression as an f64.
fn calculate_expression(lambda: f64, p: f64, q: f64, z: u32) -> f64 {
    // Input validation (optional but good practice)
    if p == 0.0 {
        panic!("Parameter 'p' cannot be zero.");
    }
    if lambda.is_sign_negative() {
        eprintln!("Warning: lambda should typically be non-negative in this context.");
    }

    let ratio_qp = q / p;
    let exp_neg_lambda = (-lambda).exp(); // e^{-\lambda}
    let mut sum: f64 = 0.0;

    for k in 0..=z {
        // 1. Poisson probability term: P(k; \lambda) = (\lambda^k e^{-\lambda}) / k!
        let lambda_pow_k = lambda.powi(k as i32);
        let k_factorial = factorial(k);
        let poisson_term = (lambda_pow_k * exp_neg_lambda) / k_factorial;

        // 2. The second term: (1 - (q/p)^(z-k))
        let exponent = (z - k) as i32;
        let power_term = 1.0 - ratio_qp.powi(exponent);

        // 3. Add the product to the total sum
        sum += poisson_term * power_term;
    }

    // 4. Final result: 1 - sum
    1.0 - sum
}

fn main() {
    // Example usage:
    let lambda = 2.5; // Example \lambda value
    let p = 0.6; // Example p value
    let q = 0.4; // Example q value (assuming p+q is related, e.g., p+q=1 for a binomial context)
    let z = 3; // Example upper limit z

    let result = calculate_expression(lambda, p, q, z);

    println!(
        "Parameters: \u{03BB} = {}, p = {}, q = {}, z = {}",
        lambda, p, q, z
    );
    println!("The result of the expression is: {:.6}", result);
}
//