#include <random>
#include <ctime>
#include <algorithm>
#include <cmath>
#include <tuple>
#include <iostream>
#include <iomanip>

#define WARM_UP_ENGINE

using vector = std::tuple<double,double,double> ;

const auto x = [] ( vector v ) { return std::get<0>(v) ; } ;
const auto y = [] ( vector v ) { return std::get<1>(v) ; } ;
const auto z = [] ( vector v ) { return std::get<2>(v) ; } ;

double norm( vector v ) { return std::sqrt( x(v)*x(v) + y(v)*y(v) + z(v)*z(v) ) ; }

vector random_vector()
{
    #ifdef WARM_UP_ENGINE // if high quality pseudo-randomness is criticsl

        static int sseq[ std::mt19937::state_size ] ;
        const static bool once = ( std::srand( std::time(nullptr) ),
                                   std::generate( std::begin(sseq), std::end(sseq), std::rand ),
                                   true ) ;
        static std::seed_seq seed_seq( std::begin(sseq), std::end(sseq) ) ;
        static std::mt19937 twister(seed_seq) ;

    #else // no warm up is required, should be adequate for normal use

        static std::mt19937 twister( std::time(0) ) ;

    #endif // WARM_UP_ENGINE

    static std::uniform_real_distribution<double> distr( -1000, 1000 ) ;

    return vector( distr(twister), distr(twister), distr(twister) ) ;
}

vector random_unit_vector()
{
    auto v = random_vector() ;
    constexpr double epsilon = 0.01 ;
    double m = norm(v) ;
    if( m > epsilon ) return vector( x(v)/m, y(v)/m, z(v)/m ) ;
    else return random_unit_vector() ;
}

int main()
{
    std::cout << std::fixed << std::setprecision(3) ;

    for( int i = 0 ; i < 10 ; ++i )
    {
        auto v = random_unit_vector() ;
        std::cout << std::showpos << "(" << x(v) << ',' << y(v) << ',' << z(v)
                  << ")  " ;

        // move the random unit vector to a random location
        static std::mt19937 twister( std::time(0) ) ;
        static std::uniform_real_distribution<double> distr( -10, 10 ) ;

        double x1 = distr(twister), y1 = distr(twister), z1 = distr(twister) ;
        double x2 = x1 + x(v), y2 = y1 + y(v), z2 = z1+z(v) ;

        std::cout << std::showpos << "    (" << x1 << ',' << y1 << ',' << z1 << ')'
                  << " ----> (" << x2 << ',' << y2 << ',' << z2 << ")\n\n" ;
    }
}