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) 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.

 
#ifndef BERTINI_RANDOM_HPP
#define BERTINI_RANDOM_HPP
 
 
 
 
#include "bertini2/mpfr_complex.hpp"
#include "bertini2/records/derive.hpp"
#include <boost/random.hpp>
#include <random>
#include <cstdint>
 
 
namespace bertini
{

    std::mt19937& ThreadEngine();

    // 0 is a seed here and entropy in the classic format: each follows its users' convention (ADR-0072)
    void SetGlobalSeed(unsigned long seed);

    unsigned long SetGlobalSeedFromEntropy();

    unsigned long ApplyClassicRandomSeed(unsigned long randomseed);

    unsigned long GetGlobalSeed();

    void ReseedThisThread(uint64_t stream_key);

    unsigned long DerivedWorkerSeed(uint64_t worker_index);

    unsigned long DeriveSolveSeed();
 

    template <unsigned long digits = 50>
    inline
    mpz_int RandomInt()
    {
        // pinned draw (b2rand/1, ADR-0044): uniform on [-2^(digits/log10(2)), +same]
        return records::ThreadDrawStream().IntSymmetric(mpz_int(1) << digits*1000L/301L);
    }
 

    template <unsigned long digits = 50>
    mpq_rational RandomRat()
    {
        // pinned draws (b2rand/1, ADR-0044).  A zero denominator is redrawn (deterministically):
        // the legacy draw could in principle hand mpq a denominator of 0.
        mpz_int const bound = mpz_int(1) << digits*1000L/301L;
        auto& stream = records::ThreadDrawStream();
        mpz_int const num = stream.IntSymmetric(bound);
        mpz_int den = stream.IntSymmetric(bound);
        while (den == 0)
            den = stream.IntSymmetric(bound);
        return mpq_rational(num, den);
    }
 

    template <unsigned int length_in_digits>
    real_mp RandomMp()
    {
        // Pinned draw (b2rand/1, ADR-0044) from the single per-thread stream shared by
        // every random type, so SetGlobalSeed()/ReseedThisThread() control them all
        // uniformly.  Uniform on [0,1) at length_in_digits digits, bit-identical on
        // every platform.
        real_mp a{records::ThreadDrawStream().UnitRealMp(length_in_digits)};
        return a;
    }

    template <unsigned int length_in_digits>
    void RandomMpAssign(real_mp & a)
    {
        a = RandomMp<length_in_digits>();
    }

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

    real_mp RandomMp();

    real_mp RandomMp(unsigned num_digits);

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

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

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

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

    inline complex RandomReal()
    {
        return complex(RandomMp(real_mp(-1),real_mp(1))); // ,0
    }

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

    inline complex rand()
    {
        // DRAW ORDER IS A CONTRACT (cross-platform reproducibility, 2026-07-03): the two
        // component draws are sequenced EXPLICITLY -- real first, then imaginary.  Never
        // put two draws in one full-expression: C++ argument evaluation order is
        // unspecified, and gcc really does order them differently on x86_64 vs aarch64
        // (found as an architecture-split seeded-homotopy digest: every (re, im) pair of
        // every seeded coefficient was transposed between the two).
        real_mp const re = RandomMp(real_mp(-1),real_mp(1));
        real_mp const im = RandomMp(real_mp(-1),real_mp(1));
        return complex( re, im );
    }
 

    inline complex rand_unit()
    {
        // draw order is a contract: real first, then imaginary (see rand())
        real_mp const re = RandomMp(real_mp(-1),real_mp(1));
        real_mp const im = RandomMp(real_mp(-1),real_mp(1));
        complex returnme( re, im );
        return returnme / abs(returnme);   // normalize to modulus 1 (NOT sqrt(abs), which left modulus sqrt|z|)
    }

    inline complex RandomUnit()
    {
        return rand_unit();
    }
 

    inline complex rand_bounded_modulus()
    {
        // draw order is a contract: real first, then imaginary (see rand())
        real_mp re = RandomMp(real_mp(-1),real_mp(1));
        real_mp im = RandomMp(real_mp(-1),real_mp(1));
        complex z( re, im );
        auto m = abs(z);
        while (m == 0)   // measure-zero, but a zero coefficient is degenerate -- redraw
        {
            re = RandomMp(real_mp(-1),real_mp(1));
            im = RandomMp(real_mp(-1),real_mp(1));
            z = complex( re, im );
            m = abs(z);
        }
        return z / sqrt(m);
    }

    inline complex RandomComplexBoundedModulus()
    {
        return rand_bounded_modulus();
    }

    inline
    void RandomComplexBoundedModulusAssign(complex & a, unsigned num_digits)
    {
        auto cached = ThreadPrecision();
        SetThreadPrecision(num_digits);
        a.precision(num_digits);
 
        // draw order is a contract: real first, then imaginary (see rand())
        real_mp re = RandomMp(real_mp(-1),real_mp(1),num_digits);
        real_mp im = RandomMp(real_mp(-1),real_mp(1),num_digits);
        complex z( re, im );
        auto m = abs(z);
        while (m == 0)
        {
            re = RandomMp(real_mp(-1),real_mp(1),num_digits);
            im = RandomMp(real_mp(-1),real_mp(1),num_digits);
            z = complex( re, im );
            m = abs(z);
        }
        a = std::move(z / sqrt(m));
        SetThreadPrecision(cached);
    }

    inline
    complex RandomComplexBoundedModulus(unsigned num_digits)
    {
        complex a;
        RandomComplexBoundedModulusAssign(a, num_digits);
        return a;
    }
 

    inline complex rand_real_bounded_modulus()
    {
        complex z( RandomMp(real_mp(-1),real_mp(1)) );   // imaginary part 0
        auto m = abs(z);
        while (m == 0)   // measure-zero, but a zero coefficient is degenerate -- redraw
        {
            z = complex( RandomMp(real_mp(-1),real_mp(1)) );
            m = abs(z);
        }
        return z / sqrt(m);   // stays real (imag 0 / real = 0)
    }

    inline complex RandomRealBoundedModulus()
    {
        return rand_real_bounded_modulus();
    }

    inline
    void RandomRealBoundedModulusAssign(complex & a, unsigned num_digits)
    {
        auto cached = ThreadPrecision();
        SetThreadPrecision(num_digits);
        a.precision(num_digits);
 
        complex z( RandomMp(real_mp(-1),real_mp(1),num_digits) );   // imaginary part 0
        auto m = abs(z);
        while (m == 0)
        {
            z = complex( RandomMp(real_mp(-1),real_mp(1),num_digits) );
            m = abs(z);
        }
        a = std::move(z / sqrt(m));   // stays real
        SetThreadPrecision(cached);
    }

    inline
    complex RandomRealBoundedModulus(unsigned num_digits)
    {
        complex a;
        RandomRealBoundedModulusAssign(a, num_digits);
        return a;
    }

    inline
    void rand_assign(complex & a, unsigned num_digits)
    {
        auto cached = ThreadPrecision();
        SetThreadPrecision(num_digits);
 
        // draw order is a contract: real first, then imaginary (see rand())
        real_mp const re = RandomMp(num_digits);
        real_mp const im = RandomMp(num_digits);
        complex_mp temp( re, im );
        a = std::move(temp);
        SetThreadPrecision(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 = ThreadPrecision();
        SetThreadPrecision(num_digits);
        a.precision(num_digits);
 
        // draw order is a contract: real first, then imaginary (see rand())
        real_mp const re = RandomMp(num_digits);
        real_mp const im = RandomMp(num_digits);
        complex temp(re, im);
        a = std::move(temp/abs(temp));   // normalize to modulus 1 (NOT sqrt(abs))
        SetThreadPrecision(cached);
    }

    inline
    complex RandomUnit(unsigned num_digits)
    {
        complex a;
        RandomUnitAssign(a,num_digits);
        return a;
    }
 
}  // namespace multiprecision
 
 
} // namespaces
// }
 
 
 
#endif