# Logarithmic random distribution satisfying RandomNumberDistribution concept

I provide logarithmic_distribution to satisfy the C++ RandomNumberDistribution concept implementing the PDF

p(x) = Norm * (floor + logfac * log(x) )


with min <= x <= max. Parameters are floor, logfac, min, and max, while Norm is a normalisation constant.

This is the third, improved implementation, see also my previous post (the second version was removed in edit before any answer was received; this version is more similar to the first).

Apart from general fitness of the code, I have some specific questions (marked QUESTION in the code). Note that I'm not asking about the mathematical correctness of the implementation.

#include <iostream>
#include <random>
#include <tuple>

namespace my_random {

class logarithmic_distribution;

class parameter_for_logarithmic_distribution
: std::tuple<double,double,double,double>
{
using base = std::tuple<double,double,double,double>;
using self = parameter_for_logarithmic_distribution;
friend class logarithmic_distribution;
friend std::ostream&operator<<(std::ostream&, self const&);
friend std::istream&operator>>(std::istream&, self&);
void check_validity() const;

public:
using distribution_type = logarithmic_distribution;

// constructors: must support the same as for distribution_type
//   (including auto-generated copy constructtor and assignment)
// QUESTION: I check validity of parameters here (and in stream input)
//           Is this the correct/preferred approach? Or should we only
//           check later, when the parameters are actually used (i.e. in
//           the constructor of logarithmic_distribution::implementation)?
parameter_for_logarithmic_distribution(double floor, double logfac, double min, double max)
: base(floor,logfac,min,max)
{
check_validity();
}
parameter_for_logarithmic_distribution(double floor, double logfac, double max)
: parameter_for_logarithmic_distribution(floor,logfac,1.0,max)
{}
parameter_for_logarithmic_distribution()
: parameter_for_logarithmic_distribution(0,1,1,2)
{}

// comparison (strictly only == is required)
friend bool operator==(self const&p1, self const&p2)
{
return static_cast<base const&>(p1) == static_cast<base const&>(p2);
}
friend bool operator!=(self const&p1, self const&p2)
{
return !(p1==p2);
}

// getter methods
//   will be inherited by distribution, thus automatically satisfying the
//   corresponding requirement from concept RandomNumberDistribution
double floor() const { return std::get<0>(*this); }
double logfac() const { return std::get<1>(*this); }
double min() const { return std::get<2>(*this); }
double max() const { return std::get<3>(*this); }
};  // class my_random::parameter_for_logarithmic_distribution

// random distribution  p(x) = Norm * (floor + logfac * log(x))

//   satisfying the RandomNumberDistribution concept, see
//   http://en.cppreference.com/w/cpp/concept/RandomNumberDistribution
class logarithmic_distribution
// QUESTION: I derive logarithmic_distribution from its param_type, which
//           has the benefit that all the getters are inherited too, thus
//           satisfying the corresponding requirement. The same holds for
//           comparison and stream output.
//           However, is the alternative (having a member 'parameters')
//           perhaps preferrable?
: public parameter_for_logarithmic_distribution
{
// further data
double min_log_min;
double cdf_at_max;
std::uniform_real_distribution<double> uniform;

public:
// public types required for RandomNumberDistribution
using result_type = double;
using param_type = parameter_for_logarithmic_distribution;

// constructors and assignment as required for RandomNumberDistribution
//   (including auto-generated copy constructor and assignment)
logarithmic_distribution(param_type const&);
logarithmic_distribution(double floor, double logfac, double min, double max)
: logarithmic_distribution(param_type(floor,logfac,min,max))
{}
logarithmic_distribution(double floor, double logfac, double max)
: logarithmic_distribution(param_type(floor,logfac,max))
{}
logarithmic_distribution()
: logarithmic_distribution(param_type())
{}

// setters as required for RandomNumberDistribution
void param(param_type const&p)
{
if(static_cast<param_type const&>(*this) != p)
operator=(logarithmic_distribution(p));
}
void reset() const {}

// getters as required for RandomNumberDistribution (including inherited)
param_type param() const { return *this; }

// state input (differs from parameter input)
friend std::istream&operator>>(std::istream&is, logarithmic_distribution&d)
{
param_type p;
is >> p;
d.param(p);
return is;
}

// random number generation interface (as required)
template<typename Generator>
result_type operator()(Generator &g)
{
return sample(uniform(g));
}

template<typename Generator>
result_type operator()(Generator &g, param_type const &p)
{
return logarithmic_distribution(p)(g);
}

// miscellaneous getters, not required
double pdf(double) const;
double cdf(double) const;
double mean() const;
double variance() const;

private:
double _pdf(double) const;
double _cdf(double, double) const;
double sample(double) const;
};// class my_random::logarithmic_distribution::implementation
} // namespace my_random


source file

#include "logarithmic_distribution.h"
#include <exception>
#include <cmath>

namespace {
// definite integral of PDF(x) from min to x
inline double log_dist_cdf(double x, double logx, double floor, double logfac, double min, double min_log_min)
{
return (floor-logfac) * (x-min) + logfac * (x*logx - min_log_min);
}
// indefinite integral of  2 * x * PDF(x), not normalised
inline double log_dist_for_mean(double x, double logx, double floor, double logfac)
{
return x*x*(floor+logfac*(logx-0.5));
}
// indefinite integral of  3 * x^2 * PDF(x), not normalised
inline double log_dist_for_var(double x, double logx, double floor, double logfac)
{
return x*x*x*(floor+logfac*(logx-0.3333333333333333333));
}
}   // namespace {

namespace my_random {

void parameter_for_logarithmic_distribution::check_validity() const
{
if(min() <= 0.)
throw std::runtime_error("logarithmic pdf: min <= 0");
if(max() <= min())
throw std::runtime_error("logarithmic pdf: min >= max");
if((floor()+logfac()*std::log(min())) < 0.0)
throw std::runtime_error("logarithmic pdf <0 at x=min");
}

// I/O of parameters according to RandomNumberDistribution concept
//
// actually, the requirement is for the distribution_type, but we delegate
// to the param_type.
//
// note:  must not change flags of I/O streams
std::ostream&operator<<(std::ostream&os, parameter_for_logarithmic_distribution const&p)
{
using ios_base = std::ostream::ios_base;
const auto flags = os.flags();
const auto fill = os.fill();
const auto prec = os.precision();
const auto space = os.widen(' ');
os.flags(ios_base::scientific | ios_base::left);
os.fill(space);
os.precision(std::numeric_limits<double>::max_digits10);
os << p.floor() << space
<< p.logfac() << space
<< p.min() << space
<< p.max();
os.flags(flags);
os.fill(fill);
os.precision(prec);
return os;
}

std::istream&operator>>(std::istream&is, parameter_for_logarithmic_distribution &p)
{
using ios_base = std::istream::ios_base;
const auto flags = is.flags();
is.flags(ios_base::dec | ios_base::skipws);
is >> std::get<0>(p)
>> std::get<1>(p)
>> std::get<2>(p)
>> std::get<3>(p);
is.flags(flags);
p.check_validity();
return is;
}

inline double logarithmic_distribution::_pdf(double logx) const
{
return floor() + logfac() * logx;
}

inline double logarithmic_distribution::_cdf(double x, double logx) const
{
return log_dist_cdf(x,logx,floor(),logfac(),min(),min_log_min);
}

logarithmic_distribution::logarithmic_distribution(param_type const&p)
: param_type  (p)
, min_log_min (p.min() * std::log(p.min()))
, cdf_at_max  (log_dist_cdf(p.max(),std::log(p.max()),
p.floor(),p.logfac(),p.min(),min_log_min))
, uniform     (0.0,cdf_at_max)
{}

double logarithmic_distribution::pdf(double x) const
{
return _pdf(std::log(x))/cdf_at_max;
}

double logarithmic_distribution::cdf(double x) const
{
return _cdf(x,std::log(x))/cdf_at_max;
}

double logarithmic_distribution::mean() const
{
return 0.5*(log_dist_for_mean(max(),std::log(max()),floor(),logfac()) -
log_dist_for_mean(min(),std::log(min()),floor(),logfac())
)/cdf_at_max;
}

double logarithmic_distribution::variance() const
{
auto logmin = std::log(min());
auto logmax = std::log(max());
auto mom1 = 0.5*
(log_dist_for_mean(max(),logmax,floor(),logfac()) -
log_dist_for_mean(min(),logmin,floor(),logfac()))/cdf_at_max;
auto mom2 = 0.33333333333333333*
(log_dist_for_var (max(),logmax,floor(),logfac()) -
log_dist_for_var (min(),logmin,floor(),logfac()))/cdf_at_max;
return mom2-mom1*mom1;
}

// given a random value C for the non-normalised CDF(x), find x
double logarithmic_distribution::sample(const double C) const
{
using std::abs;
if(C <= 0)
return min();
if(C >= cdf_at_max)
return max();
// bracketed Newton-Raphson for f(x) = CDF(x)-C  (not normalised)
const unsigned max_iterations = 100;
const double eps = 10*std::numeric_limits<double>::epsilon();
auto xl = min();
auto xh = max();
auto dxo= xh-xl;
auto dx = dxo;
auto x  = xl+C*dx/cdf_at_max;                   // linear interpolation
auto lx = std::log(x);
auto df = _pdf(lx);                             // df/dx = PDF(x)
auto f  = _cdf(x,lx) - C;                       // f(x)  = CDF(x)-C
for(unsigned it=0; it!=max_iterations; ++it) {
auto xo = x;
if((f>0? ((x-xh)*f>df) : ((x-xl)*f<df))     // NR outside of bracket
|| (abs(f+f) > abs(dxo*df))) {           // NR convergence is slow
dxo= dx;                                // bi-section
dx = 0.5*(xh-xl);
x  = xl+dx;
} else {
dxo= dx;                                // Newton-Raphson
dx = f/df;
x -= dx;
}
if((xo<=x && xo>=x) || abs(dx)<=eps*abs(x)) // converged?
return x;
lx = std::log(x);
df = _pdf(lx);
f  = _cdf(x,lx) - C;
if(f<0) xl=x; else xh=x;                    // maintain bracket
}
throw
std::runtime_error("logarithmic_distribution: exceeded 100 iterations");
}
} // namespace my_random


namespace my_random {


Oh, dear! A branding crisis. Be bold! Pick some crazy name that google has zero hits for, and run with it.

Also, many identifiers mention "logarithmic_distribution", which is a lovely and standard name. Maybe consider abbreviating to "log_distribution"?

tiny nit, typo: "constructtor"

// QUESTION: I check validity of parameters [in ctor]

Yes! This is definitely the correct and preferred approach. Bugs are inevitable, but we wish to write shallow bugs, which are immediately called out and are trivial to diagnose / repair. DbC asks us to check the pre-conditions early, so we can immediately pin the blame on the caller, before backtrace shows we're a few convoluted levels down in library calls.

parameter_for_logarithmic_distribution(floor, logfac, 1.0, max)


It seems we're setting min to one here. Maybe a comment to that effect? Maybe enforce that min shall be >= 1 ?

Anyway, the subsequent defaulting looks perfectly nice.

friend bool operator==


Ok, I'm going to be super picky here, because I think that's the level of feedback you were soliciting.

This is a perfectly nice implementation. But I will note that the ctor admits of several different NaNs, potentially. And comparing NaNs is trouble, it doesn't always offer intuitive results. In particular, operator!= might not do what someone expects.

Consider prohibiting inf / NaN in the ctor. There are some lovely validity checks. But asking whether min <= 0 won't flush out a NaN min, nor a positive infinity. I guess it's a matter of taste, depending on your level of paranoia.

You have a deadlink. Apparently the concept was elevated to a named requirement. New page is at

https://en.cppreference.com/w/cpp/named_req/RandomNumberDistribution

You ask if an alternative of having member 'parameters' might be preferable.

No, I have not yet seen anything in the code that would make me think that.

I appreciate that you call out "required" versus "not required". It's fine the way it is. But here is a small suggestion.

Maybe offer a pair of classes, one inheriting from the other? With the base being documented as offering "only what the spec requires", and the other being the one you anticipate an app developer would typically link against.

That way the code tells us what is required, rather than comments. In a way that automated unit tests can exploit.

The 0.3333333333333333333 constant is perhaps less convenient than 1. / 3. ? I mean, it's gonna be a single compile-time constant, either way.

Same constant appears in variance mom2. If you keep it as-is, consider making it a MANIFEST_CONSTANT.

    , min_log_min (p.min() * std::log(p.min()))
, cdf_at_max  (log_dist_cdf(p.max(),std::log(p.max()),


style nit: Consider preferring , comma at end-of-line.

if ((xo<=x && xo>=x) || ...


That seems, ummmm, slightly tricky. Why not equality? Are we looking for NaN, there? Would a stricter initial validity check help us to simplify this?

Overall?

Looks great, ship it!

• "I mean, it's gonna be a single compile-time constant, either way." - OP states C++11 so, no, it is not mandated to be. Commented Jan 14, 2023 at 3:12
• Sure, a compliant compiler is free to keep taking reciprocal of 3 every time if it wants. I'm just not aware of one that people would choose to use that actually does that. Constant folding is among the most trivial of optimization passes. Is there a compiler / arch / -O setting on godbolt.org you'd care to highlight? My concern was "how shall author clearly convey intent to collaborators?" rather than being focused on communicating details to the machine. When I read the source and don't count decimal digits, it isn't obvious to me we have the proper 53 bits of mantissa.
– J_H
Commented Jan 14, 2023 at 3:41
• oh, wow. Had completely forgotten about this post and now an answer after almost 5 years. I think the NaN tests are really missing in check_validity. The pseudo-equality test was to avoid compiler warnings about floating-point equality tests. You said "ship it" -- where though? Commented Feb 5, 2023 at 15:18