#include <boost/random/variate_generator.hpp>
#include <boost/random/uniform_on_sphere.hpp>
#include <boost/random.hpp>
#include <iostream>
#include <vector>

int main(int argc, const char *argv[])
{
    // Typedefs to the generator, container, distribution and variate generator
    typedef boost::mt19937                                              generator_t;
    typedef std::vector<double>                                          container_t;
    typedef boost::uniform_on_sphere<double, container_t>                spherical_dist_t;
    typedef boost::variate_generator<generator_t*, spherical_dist_t >   variate_generator_t;

    //Specify the dimensions to try and the the number of samples per dimension
    const unsigned int        minDim = 1u;
    const unsigned int        maxDim = 25u;
    const unsigned int        numSamples = 10000u;

    // Track the largest error found:
    double maxError = 0.0;

    // Create the generator and specify the maximum dimension to try
    generator_t         gen;

    // Iterate over the dimensions drawing samples and announcing if they're not of magnitude 1.0
    for (unsigned int n = minDim; n <= maxDim; ++n)
    {
        // Create the distribution and the variate generator
        spherical_dist_t     dist(n);
        variate_generator_t  variate(&gen, dist);

        // Draw some samples
        for (unsigned int i = 0u; i < numSamples; ++i)
        {
            // Draw a sample
            container_t sample = variate();

            // Calculate it's magnitude
            double magnitude = 0.0;
            for (unsigned int j = 0u; j < n; ++j)
            {
                magnitude = magnitude + sample.at(j)*sample.at(j);
            }
            magnitude = std::sqrt(magnitude);

            // Calculate the error in the magnitude:
            double magnitudeError = std::abs(1.0 - magnitude);

            // Update the maximum error:
            maxError = std::max(maxError, magnitudeError);

            // Check if this is above an arbitrary tolerance and is an example of the bug:
            if ( magnitudeError > 1E-12)
            {
                std::cout << std::scientific << "A " << n << " dimensional vector has a magnitude of " << magnitude << " which is " << magnitudeError << " away from 1.0" << std::endl;
            }
        }
    }

    // Statistics:
    std::cout << "Drawing " << (maxDim - minDim + 1u)*numSamples << " samples from " << minDim << "D to " << maxDim << "D found a maximum absolute magnitude error of " << maxError << std::endl;

    return 0;
}