num_traits.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//num_traits.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.
//
//num_traits.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 num_traits.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_NUM_TRAITS_HPP
#define BERTINI_NUM_TRAITS_HPP
 
#include <random>
#include "bertini2/records/draw_functions.hpp"
#include <complex>
#include <cmath>
#include "bertini2/mpfr_complex.hpp"
#include "bertini2/random.hpp"
 
 
namespace bertini
{
    using NumErrorT = double;

    template<typename T>
    T RandomUnit();

    template<typename T>
    struct NumTraits
    {};
 
 

    template <> struct NumTraits<double>
    {
        inline static unsigned NumDigits()
        {
            return 16;
        }

        inline static unsigned NumFuzzyDigits()
        {
            return 14;
        }

        inline
        static unsigned TolToDigits(NumErrorT tol)
        {
            return static_cast<unsigned>(ceil(-log10(tol)));
        }

        inline static
        double FromString(std::string const& s)
        {
            return boost::lexical_cast<double>(s);
        }

        inline static
        double FromRational(mpq_rational const& n, unsigned /* precision */)
        {
            return double(n);
        }
 
        using Real = double;  
        using Complex = complex_dbl;  
    };
 

    template <> struct NumTraits<complex_dbl >
    {
        inline static unsigned NumDigits()
        {
            return 16;
        }

        inline static unsigned NumFuzzyDigits()
        {
            return 14;
        }

        inline static
        complex_dbl FromString(std::string const& s)
        {
            return boost::lexical_cast<complex_dbl>(s);
        }

        inline static
        complex_dbl FromString(std::string const& s, std::string const& t)
        {
            return complex_dbl(boost::lexical_cast<double>(s),boost::lexical_cast<double>(t));
        }

        inline static
        complex_dbl FromRational(mpq_rational const& n, unsigned /* precision */)
        {
            return complex_dbl(static_cast<double>(n),0);
        }
 
        using Real = double;  
        using Complex = complex_dbl;  
    };
 

    inline
    unsigned PrecisionIncrement()
    {
        return 10;
    }

    constexpr
    unsigned DoublePrecision()
    {
        return 16;
    }

    constexpr
    unsigned LowestMultiplePrecision()
    {
        return 20;
    }

    constexpr
    unsigned MaxPrecisionAllowed()
    {
        return 1000;
    }

    inline
    unsigned Precision(double)
    {
        return DoublePrecision();
    }

    inline
    void Precision(double, unsigned prec)
    {
        if (prec!=DoublePrecision())
        {
            std::stringstream err_msg;
            err_msg << "trying to change precision of a double to " << prec;
            throw std::runtime_error(err_msg.str());
        }
    }

    inline
    unsigned Precision(complex_dbl)
    {
        return DoublePrecision();
    }

    inline
    void Precision(complex_dbl, unsigned prec)
    {
        if (prec!=DoublePrecision())
        {
            std::stringstream err_msg;
            err_msg << "trying to change precision of a double to " << prec;
            throw std::runtime_error(err_msg.str());
        }
    }

    inline
    complex_dbl rand_complex()
    {
        using std::abs;
        using std::sqrt;
        // 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).
        double const re = records::DrawSymmetricDouble();
        double const im = records::DrawSymmetricDouble();
        complex_dbl returnme(re, im);
        return returnme / sqrt( abs(returnme));
    }

    template <> inline
    complex_dbl RandomUnit<complex_dbl >()
    {
        // draw order is a contract: real first, then imaginary (see rand_complex)
        double const re = records::DrawSymmetricDouble();
        double const im = records::DrawSymmetricDouble();
        complex_dbl returnme(re, im);
        return returnme / abs(returnme);
    }

    template <>
    inline
    complex_mp RandomUnit<complex_mp>()
    {
        return multiprecision::RandomUnit();
    }

    template <>
    inline
    real_mp RandomUnit<real_mp>()
    {
        return real_mp( records::DrawSymmetricDouble() < 0.0 ? -1 : 1 );
    }

    template <typename NumberType>
    inline
    NumberType RandomOrthonormalSeed()
    {
        return RandomUnit<NumberType>();
    }

    template <>
    inline
    real_mp RandomOrthonormalSeed<real_mp>()
    {
        return real_mp( records::DrawSymmetricDouble() );
    }
}// re: namespace bertini
 
 
 
 
 
 
 
 
 
 
 
 
namespace bertini {
 
 

    template <> struct NumTraits<real_mp>
    {
        inline static unsigned NumDigits()
        {
            return ThreadPrecision();
        }

        inline static unsigned NumFuzzyDigits()
        {
            return ThreadPrecision()-3;
        }

        inline
        static unsigned TolToDigits(real_mp tol)
        {
            real_mp b = ceil(-log10(tol));
            return b.convert_to<unsigned int>();
        }

        inline static
        real_mp FromString(std::string const& s)
        {
            return real_mp(s);
        }

        inline static
        real_mp FromRational(mpq_rational const& n, unsigned precision)
        {
            return real_mp(n,precision);
        }
 
        using Real = real_mp;  
        using Complex = complex_mp;  
    };
 
 

    template <> struct NumTraits<complex_mp>
    {
        inline static unsigned NumDigits()
        {
            return ThreadPrecision();
        }

        inline static
        complex_mp FromString(std::string const& s)
        {
            return complex_mp(s);
        }

        inline static
        complex_mp FromString(std::string const& s, std::string const& t)
        {
            return complex_mp(s,t);
        }

        inline static
        complex_mp FromRational(mpq_rational const& n, unsigned precision)
        {
            return complex_mp(n,0,precision);
        }
 
        using Real = real_mp;  
        using Complex = complex_mp;  
    };

    template <> struct NumTraits<mpq_rational>
    {
        inline static
        mpq_rational FromString(std::string const& str)
        {
            std::string s = str;
            bool negative = false;
            if (!s.empty() && s[0] == '-') { negative = true; s = s.substr(1); }
            else if (!s.empty() && s[0] == '+') { s = s.substr(1); }
 
            int exp_shift = 0;
            auto e_pos = s.find_first_of("eE");
            if (e_pos != std::string::npos) {
                exp_shift = std::stoi(s.substr(e_pos + 1));
                s = s.substr(0, e_pos);
            }
 
            auto dot_pos = s.find('.');
            int decimal_places = 0;
            if (dot_pos != std::string::npos) {
                decimal_places = static_cast<int>(s.size()) - static_cast<int>(dot_pos) - 1;
                s.erase(dot_pos, 1);
            }
 
            if (s.empty() || s.find_first_not_of('0') == std::string::npos) {
                s = "0";
            } else {
                s = s.substr(s.find_first_not_of('0'));
            }
            mpz_int numer(s);
            if (negative) numer = -numer;
 
            int net_exp = decimal_places - exp_shift;
            if (net_exp > 0) {
                mpz_int denom = 1;
                for (int i = 0; i < net_exp; ++i) denom *= 10;
                return mpq_rational(numer, denom);
            } else if (net_exp < 0) {
                mpz_int mult = 1;
                for (int i = 0; i < -net_exp; ++i) mult *= 10;
                return mpq_rational(numer * mult, 1);
            } else {
                return mpq_rational(numer, 1);
            }
        }
    };
 
}
 
#endif