qsys_md1.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_QSYS_QSYS_MD1_H
#define LINE_API_QSYS_QSYS_MD1_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the M/D/1 queue: the M/G/1 moments with deterministic
 * service plus the queue-length distribution.
 *
 * Port of matlab/src/api/qsys/qsys_md1.m. The mean measures are qsys_mg1's
 * with cs = 0 and E[S^3] = 1/mu^3. The distribution follows from the
 * embedded-chain recursion with a_0 = exp(-rho) and the Poisson tails
 * Abar_k = P(N > k):
 *
 *   p_0 = 1 - rho,   p_n = (p_0 Abar_{n-1} + sum_{j=1}^{n-1} p_j Abar_{n-j}) / a_0,
 *
 * run until the mass left is below 1e-14 (or to nmax when given).
 *
 * ARITHMETIC: double, as the reference.
 */
 
#include <cmath>
#include <vector>
 
#include "line/api/qsys/qsys_mg1.h"
#include "line/api/qsys/qsys_res.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/**
 * @brief Exact M/D/1 queue.
 *
 * @param lambda arrival rate
 * @param mu     service rate (service time 1/mu)
 * @param nmax   truncation of the distribution; negative for MATLAB's default
 */
inline QsysRes qsys_md1(double lambda, double mu, long nmax = -1) {
    using namespace resdetail;
    if (!(lambda > 0.0) || !(mu > 0.0)) throw Error("qsys_md1: lambda and mu must be positive scalars");
    const double rho = lambda / mu;
    if (rho >= 1.0) throw_unstable("qsys_md1", "rho", rho);
    QsysRes res = qsys_mg1_res(lambda, mu, 0.0, 1.0 / (mu * mu * mu));
 
    const long cap = nmax >= 0 ? nmax : 100000;
    std::vector<double> P(static_cast<std::size_t>(cap + 1), 0.0);
    P[0] = 1.0 - rho;
    const double a0 = std::exp(-rho);
    std::vector<double> Ab;  // Ab[k] = Abar_k, up to and including the first zero
    for (long k = 0; k <= cap; ++k) {
        Ab.push_back(gammainc_lower(rho, k + 1.0));
        if (Ab.back() == 0.0) break;
    }
    auto abar = [&](long k) { return k < static_cast<long>(Ab.size()) ? Ab[static_cast<std::size_t>(k)] : 0.0; };
    double mass = P[0];
    long last = cap;
    for (long n = 1; n <= cap; ++n) {
        double s = P[0] * abar(n - 1);
        for (long j = 1; j <= n - 1; ++j) {
            const double ab = abar(n - j);
            if (ab == 0.0) continue;
            s += P[static_cast<std::size_t>(j)] * ab;
        }
        P[static_cast<std::size_t>(n)] = s / a0;
        mass += P[static_cast<std::size_t>(n)];
        if (nmax < 0 && 1.0 - mass < 1e-14) {
            last = n;
            break;
        }
    }
    P.resize(static_cast<std::size_t>(last + 1));
    res.set("queueLengthDist", P);
    res.analyzer = "exact:M/D/1";
    return res;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MD1_H