opt_random.h

Go to the documentation of this file.

// Optimizations for random number functions, x86 version -*- C++ -*-

// Copyright (C) 2012-2015 Free Software Foundation, Inc.
//
// This file is part of the GNU ISO C++ Library.  This library is free
// software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the
// Free Software Foundation; either version 3, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

// Under Section 7 of GPL version 3, you are granted additional
// permissions described in the GCC Runtime Library Exception, version
// 3.1, as published by the Free Software Foundation.

// You should have received a copy of the GNU General Public License and
// a copy of the GCC Runtime Library Exception along with this program;
// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
// <http://www.gnu.org/licenses/>.

/** @file bits/opt_random.h
 *  This is an internal header file, included by other library headers.
 *  Do not attempt to use it directly. @headername{random}
 */

#ifndef _BITS_OPT_RANDOM_H
#define _BITS_OPT_RANDOM_H 1

#include <x86intrin.h>


#pragma GCC system_header


namespace std _GLIBCXX_VISIBILITY(default)
{
_GLIBCXX_BEGIN_NAMESPACE_VERSION

#ifdef __SSE3__
  template<>
    template<typename _UniformRandomNumberGenerator>
      void
      normal_distribution<double>::
      __generate(typename normal_distribution<double>::result_type* __f,
                 typename normal_distribution<double>::result_type* __t,
                 _UniformRandomNumberGenerator& __urng,
                 const param_type& __param)
      {
        typedef uint64_t __uctype;

        if (__f == __t)
          return;

        if (_M_saved_available)
          {
            _M_saved_available = false;
            *__f++ = _M_saved * __param.stddev() + __param.mean();

            if (__f == __t)
              return;
          }

        constexpr uint64_t __maskval = 0xfffffffffffffull;
        static const __m128i __mask = _mm_set1_epi64x(__maskval);
        static const __m128i __two = _mm_set1_epi64x(0x4000000000000000ull);
        static const __m128d __three = _mm_set1_pd(3.0);
        const __m128d __av = _mm_set1_pd(__param.mean());

        const __uctype __urngmin = __urng.min();
        const __uctype __urngmax = __urng.max();
        const __uctype __urngrange = __urngmax - __urngmin;
        const __uctype __uerngrange = __urngrange + 1;

        while (__f + 1 < __t)
          {
            double __le;
            __m128d __x;
            do
              {
                union
                {
                  __m128i __i;
                  __m128d __d;
                } __v;

                if (__urngrange > __maskval)
                  {
                    if (__detail::_Power_of_2(__uerngrange))
                      __v.__i = _mm_and_si128(_mm_set_epi64x(__urng(),
                                                             __urng()),
                                              __mask);
                    else
                      {
                        const __uctype __uerange = __maskval + 1;
                        const __uctype __scaling = __urngrange / __uerange;
                        const __uctype __past = __uerange * __scaling;
                        uint64_t __v1;
                        do
                          __v1 = __uctype(__urng()) - __urngmin;
                        while (__v1 >= __past);
                        __v1 /= __scaling;
                        uint64_t __v2;
                        do
                          __v2 = __uctype(__urng()) - __urngmin;
                        while (__v2 >= __past);
                        __v2 /= __scaling;

                        __v.__i = _mm_set_epi64x(__v1, __v2);
                      }
                  }
                else if (__urngrange == __maskval)
                  __v.__i = _mm_set_epi64x(__urng(), __urng());
                else if ((__urngrange + 2) * __urngrange >= __maskval
                         && __detail::_Power_of_2(__uerngrange))
                  {
                    uint64_t __v1 = __urng() * __uerngrange + __urng();
                    uint64_t __v2 = __urng() * __uerngrange + __urng();

                    __v.__i = _mm_and_si128(_mm_set_epi64x(__v1, __v2),
                                            __mask);
                  }
                else
                  {
                    size_t __nrng = 2;
                    __uctype __high = __maskval / __uerngrange / __uerngrange;
                    while (__high > __uerngrange)
                      {
                        ++__nrng;
                        __high /= __uerngrange;
                      }
                    const __uctype __highrange = __high + 1;
                    const __uctype __scaling = __urngrange / __highrange;
                    const __uctype __past = __highrange * __scaling;
                    __uctype __tmp;

                    uint64_t __v1;
                    do
                      {
                        do
                          __tmp = __uctype(__urng()) - __urngmin;
                        while (__tmp >= __past);
                        __v1 = __tmp / __scaling;
                        for (size_t __cnt = 0; __cnt < __nrng; ++__cnt)
                          {
                            __tmp = __v1;
                            __v1 *= __uerngrange;
                            __v1 += __uctype(__urng()) - __urngmin;
                          }
                      }
                    while (__v1 > __maskval || __v1 < __tmp);

                    uint64_t __v2;
                    do
                      {
                        do
                          __tmp = __uctype(__urng()) - __urngmin;
                        while (__tmp >= __past);
                        __v2 = __tmp / __scaling;
                        for (size_t __cnt = 0; __cnt < __nrng; ++__cnt)
                          {
                            __tmp = __v2;
                            __v2 *= __uerngrange;
                            __v2 += __uctype(__urng()) - __urngmin;
                          }
                      }
                    while (__v2 > __maskval || __v2 < __tmp);
                    
                    __v.__i = _mm_set_epi64x(__v1, __v2);
                  }

                __v.__i = _mm_or_si128(__v.__i, __two);
                __x = _mm_sub_pd(__v.__d, __three);
                __m128d __m = _mm_mul_pd(__x, __x);
                __le = _mm_cvtsd_f64(_mm_hadd_pd (__m, __m));
              }
            while (__le == 0.0 || __le >= 1.0);

            double __mult = (std::sqrt(-2.0 * std::log(__le) / __le)
                             * __param.stddev());

            __x = _mm_add_pd(_mm_mul_pd(__x, _mm_set1_pd(__mult)), __av);

            _mm_storeu_pd(__f, __x);
            __f += 2;
          }

        if (__f != __t)
          {
            result_type __x, __y, __r2;

            __detail::_Adaptor<_UniformRandomNumberGenerator, result_type>
              __aurng(__urng);

            do
              {
                __x = result_type(2.0) * __aurng() - 1.0;
                __y = result_type(2.0) * __aurng() - 1.0;
                __r2 = __x * __x + __y * __y;
              }
            while (__r2 > 1.0 || __r2 == 0.0);

            const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
            _M_saved = __x * __mult;
            _M_saved_available = true;
            *__f = __y * __mult * __param.stddev() + __param.mean();
          }
      }
#endif


_GLIBCXX_END_NAMESPACE_VERSION
} // namespace


#endif // _BITS_OPT_RANDOM_H