Algorithms and software are presented for efficiently computing reference solutions of the general rate model with proven error bounds. Moreover, algorithms and software are presented for efficiently computing moments of arbitrary order. The methods are based on numerical inverse Laplace transform, and support both quasi-stationary and dynamic linear binding models. The inlet concentration profiles are treated in a most general way using piecewise cubic polynomials. Algorithmic differentiation obviates manual derivation of the required derivatives. Arbitrary precision arithmetics are applied for minimizing numerical roundoff errors, and several convergence acceleration techniques are evaluated. The implemented software package is freely available as open source on GitHub. (C) 2015 Elsevier Ltd. All rights reserved.