random.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
// bertini2/random.hpp 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 of the License, or
//(at your option) any later version.
//
// bertini2/random.hpp 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.
//
//You should have received a copy of the GNU General Public License
//along with  bertini2/random.hpp.  If not, see <http://www.gnu.org/licenses/>.
//
// Copyright(C) 2015 - 2021 by Bertini2 Development Team
//
// See <http://www.gnu.org/licenses/> for a copy of the license, 
// as well as COPYING.  Bertini2 is provided with permitted 
// additional terms in the b2/licenses/ directory.
 
// individual authors of this file include:
// silviana amethyst, university of wisconsin eau claire

 
#ifndef BERTINI_RANDOM_HPP
#define BERTINI_RANDOM_HPP
 
 
 
 
#include "bertini2/mpfr_complex.hpp"
#include <boost/random.hpp>
 
 
namespace bertini
{
    template <unsigned long digits = 50>
    inline
    mpz_int RandomInt()
    {
        using namespace boost::random;
        static mt19937 mt;
        static uniform_int_distribution<mpz_int> ui(-(mpz_int(1) << digits*1000L/301L), mpz_int(1) << digits*1000L/301L);
        return ui(mt);
    }
    
    
    template <unsigned long digits = 50>
    mpq_rational RandomRat()
    {
        using namespace boost::random;
        static mt19937 mt;
        static uniform_int_distribution<mpz_int> ui(-(mpz_int(1) << digits*1000L/301L), mpz_int(1) << digits*1000L/301L);
        return mpq_rational(ui(mt),ui(mt));
    }
 

    template <unsigned int length_in_digits>
    mpfr_float RandomMp()
    {   
 
        using namespace boost::multiprecision;
        using namespace boost::random;
 
        static uniform_real_distribution<number<mpfr_float_backend<length_in_digits>, et_on> > distribution(0,1);
        static independent_bits_engine<mt19937, length_in_digits*1000L/301L, mpz_int> bit_generator;
 
        mpfr_float a{distribution(bit_generator)};
        return a;
    }
    
    template <unsigned int length_in_digits>
    void RandomMpAssign(mpfr_float & a)
    {   
        a = RandomMp<length_in_digits>();
    }

    template <unsigned int length_in_digits>
    mpfr_float RandomMp(const mpfr_float & left, const mpfr_float & right)
    {
        return (right-left)*RandomMp<length_in_digits>()+left;
    }
 
 

    mpfr_float RandomMp();

    mpfr_float RandomMp(unsigned num_digits);

    mpfr_float RandomMp(const mpfr_float & a, const mpfr_float & b);

    mpfr_float RandomMp(const mpfr_float & a, const mpfr_float & b, unsigned num_digits);

    void RandomMpAssign(mpfr_float & a, unsigned num_digits);
 
    
 
    
} // re: namespace bertini
 
 
 
 
namespace bertini{
 
namespace multiprecision{
 
 
using complex = bertini::mpfr_complex;
using bertini::RandomMp;
 
 

    inline 
    void RandomRealAssign(complex & a, unsigned num_digits)
    {
        auto cached = DefaultPrecision();
        DefaultPrecision(num_digits);
        complex temp(RandomMp(mpfr_float(-1),mpfr_float(1),num_digits)); // ,0
        a.swap(temp);
        DefaultPrecision(cached);
    }

    inline complex RandomReal()
    {
        return complex(RandomMp(mpfr_float(-1),mpfr_float(1))); // ,0
    }
    
    inline complex RandomReal(unsigned num_digits)
    {
        auto cached = DefaultPrecision();
        DefaultPrecision(num_digits);
        auto result = complex(RandomMp(mpfr_float(-1),mpfr_float(1),num_digits));// ,0
        DefaultPrecision(cached);
        return result;
    }
 
 
 
    
        
 

    inline complex rand()
    {
        return complex( RandomMp(mpfr_float(-1),mpfr_float(1)), RandomMp(mpfr_float(-1),mpfr_float(1)) );
    }
 

    inline complex rand_unit()
    {
        complex returnme( RandomMp(mpfr_float(-1),mpfr_float(1)), RandomMp(mpfr_float(-1),mpfr_float(1)) );
        return returnme / sqrt( abs(returnme));
    }
 
    inline complex RandomUnit()
    {
        return rand_unit();
    }
 
    inline 
    void rand_assign(complex & a, unsigned num_digits)
    {
        auto cached = DefaultPrecision();
        DefaultPrecision(num_digits);
        
        mpfr_complex temp( RandomMp(num_digits), RandomMp(num_digits) );
        a = std::move(temp);
        DefaultPrecision(cached);
    }
 
    inline 
    void RandomComplexAssign(complex & a, unsigned num_digits)
    {
        rand_assign(a,num_digits);
    }
 
    inline 
    complex RandomComplex(unsigned num_digits)
    {
        complex z;
        RandomComplexAssign(z, num_digits);
        return z;
    }
 
    inline 
    void RandomUnitAssign(complex & a, unsigned num_digits)
    {
        auto cached = DefaultPrecision();
        DefaultPrecision(num_digits);
        a.precision(num_digits);
        
        complex temp(RandomMp(num_digits),RandomMp(num_digits));
        a = std::move(temp/sqrt(abs(temp)));
        DefaultPrecision(cached);
    }
 
    inline 
    complex RandomUnit(unsigned num_digits)
    {
        complex a;
        RandomUnitAssign(a,num_digits);
        return a;
    }
 
}  // namespace multiprecision
 
 
} // namespaces
// }
 
 
 
#endif