/* boost random/uniform_on_sphere.hpp header file
 *
 * Copyright Jens Maurer 2000-2001
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 *
 * Revision history
 *  2001-02-18  moved to individual header files
 */

#ifndef BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP
#define BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP

#include <vector>
#include <algorithm>     // std::transform
#include <functional>    // std::bind2nd, std::divides
#include <boost/assert.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/normal_distribution.hpp>

namespace boost {
namespace random {

/**
 * Instantiations of class template uniform_on_sphere model a
 * \random_distribution. Such a distribution produces random
 * numbers uniformly distributed on the unit sphere of arbitrary
 * dimension @c dim. The @c Cont template parameter must be a STL-like
 * container type with begin and end operations returning non-const
 * ForwardIterators of type @c Cont::iterator. 
 */
template<class RealType = double, class Cont = std::vector<RealType> >
class uniform_on_sphere
{
public:
    typedef RealType input_type;
    typedef Cont result_type;

    class param_type
    {
    public:

        typedef uniform_on_sphere distribution_type;

        /**
         * Constructs the parameters of a uniform_on_sphere
         * distribution, given the dimension of the sphere.
         */
        explicit param_type(int dim_arg = 2) : _dim(dim_arg)
        {
            BOOST_ASSERT(_dim >= 0);
        }

        /** Returns the dimension of the sphere. */
        int dim() const { return _dim; }

        /** Writes the parameters to a @c std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        {
            os << parm._dim;
            return os;
        }

        /** Reads the parameters from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            is >> parm._dim;
            return is;
        }

        /** Returns true if the two sets of parameters are equal. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        { return lhs._dim == rhs._dim; }

        /** Returns true if the two sets of parameters are different. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)

    private:
        int _dim;
    };

    /**
     * Constructs a @c uniform_on_sphere distribution.
     * @c dim is the dimension of the sphere.
     *
     * Requires: dim >= 0
     */
    explicit uniform_on_sphere(int dim_arg = 2)
      : _container(dim_arg), _dim(dim_arg) { }

    /**
     * Constructs a @c uniform_on_sphere distribution from its parameters.
     */
    explicit uniform_on_sphere(const param_type& parm)
      : _container(parm.dim()), _dim(parm.dim()) { }

    // compiler-generated copy ctor and assignment operator are fine

    /** Returns the dimension of the sphere. */
    int dim() const { return _dim; }

    /** Returns the parameters of the distribution. */
    param_type param() const { return param_type(_dim); }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        _dim = parm.dim();
        _container.resize(_dim);
    }

    /**
     * Returns the smallest value that the distribution can produce.
     * Note that this is required to approximate the standard library's
     * requirements.  The behavior is defined according to lexicographical
     * comparison so that for a container type of std::vector,
     * dist.min() <= x <= dist.max() where x is any value produced
     * by the distribution.
     */
    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
    {
        result_type result(_dim);
        if(_dim != 0) {
            result.front() = RealType(-1.0);
        }
        return result;
    }
    /**
     * Returns the largest value that the distribution can produce.
     * Note that this is required to approximate the standard library's
     * requirements.  The behavior is defined according to lexicographical
     * comparison so that for a container type of std::vector,
     * dist.min() <= x <= dist.max() where x is any value produced
     * by the distribution.
     */
    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    {
        result_type result(_dim);
        if(_dim != 0) {
            result.front() = RealType(1.0);
        }
        return result;
    }

    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() {}

    /**
     * Returns a point uniformly distributed over the surface of
     * a sphere of dimension dim().
     */
    template<class Engine>
    const result_type & operator()(Engine& eng)
    {
        using std::sqrt;
        switch(_dim)
        {
        case 0: break;
        case 1:
            {
                if(uniform_01<RealType>()(eng) < 0.5) {
                    *_container.begin() = -1;
                } else {
                    *_container.begin() = 1;
                }
            }
        case 2:
            {
                uniform_01<RealType> uniform;
                RealType sqsum;
                RealType x, y;
                do {
                    x = uniform(eng) * 2 - 1;
                    y = uniform(eng) * 2 - 1;
                    sqsum = x*x + y*y;
                } while(sqsum == 0 || sqsum > 1);
                RealType mult = 1/sqrt(sqsum);
                typename Cont::iterator iter = _container.begin();
                *iter = x * mult;
                iter++;
                *iter = y * mult;
                break;
            }
        case 3:
            {
                uniform_01<RealType> uniform;
                RealType sqsum;
                RealType x, y;
                do {
                    x = uniform(eng) * 2 - 1;
                    y = uniform(eng) * 2 - 1;
                    sqsum = x*x + y*y;
                } while(sqsum > 1);
                RealType mult = 2 * sqrt(1 - sqsum);
                typename Cont::iterator iter = _container.begin();
                *iter = x * mult;
                ++iter;
                *iter = y * mult;
                ++iter;
                *iter = 2 * sqsum - 1;
                break;
            }
        default:
            {
                detail::unit_normal_distribution<RealType> normal;
                RealType sqsum;
                do {
                    sqsum = 0;
                    for(typename Cont::iterator it = _container.begin();
                        it != _container.end();
                        ++it) {
                        RealType val = normal(eng);
                        *it = val;
                        sqsum += val * val;
                    }
                } while(sqsum == 0);
                // for all i: result[i] /= sqrt(sqsum)
                std::transform(_container.begin(), _container.end(), _container.begin(),
                               std::bind2nd(std::multiplies<RealType>(), 1/sqrt(sqsum)));
            }
        }
        return _container;
    }

    /**
     * Returns a point uniformly distributed over the surface of
     * a sphere of dimension param.dim().
     */
    template<class Engine>
    result_type operator()(Engine& eng, const param_type& parm) const
    {
        return uniform_on_sphere(parm)(eng);
    }

    /** Writes the distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, uniform_on_sphere, sd)
    {
        os << sd._dim;
        return os;
    }

    /** Reads the distribution from a @c std::istream. */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, uniform_on_sphere, sd)
    {
        is >> sd._dim;
        sd._container.resize(sd._dim);
        return is;
    }

    /**
     * Returns true if the two distributions will produce identical
     * sequences of values, given equal generators.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(uniform_on_sphere, lhs, rhs)
    { return lhs._dim == rhs._dim; }

    /**
     * Returns true if the two distributions may produce different
     * sequences of values, given equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(uniform_on_sphere)

private:
    result_type _container;
    int _dim;
};

} // namespace random

using random::uniform_on_sphere;

} // namespace boost

#endif // BOOST_RANDOM_UNIFORM_ON_SPHERE_HPP