double_extensions.hpp

Go to the documentation of this file.

//This file is part of Bertini 2.
//
//bertini2/double_extensions.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/double_extensions.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/double_extensions.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_DOUBLE_EXTENSIONS_HPP
#define BERTINI_DOUBLE_EXTENSIONS_HPP
 
#pragma once
 
#include <random>
#include "bertini2/records/draw_functions.hpp"
#include <complex>
 
 
namespace bertini{
 
// Forward declaration — full definition in random.hpp / random.cpp.
// Avoids pulling mpfr_complex.hpp into the Boost.Multiprecision include chain.
std::mt19937& ThreadEngine();
 
    // The two double-precision scalar types, named to parallel real_mp / complex_mp.
    // real_dbl is an explicit alias for double so a future change lives in one place.
    using real_dbl = double;  
    using complex_dbl = std::complex<double>;  

    inline
    std::complex<double> operator*(unsigned i, std::complex<double> z)
    {
        z*=i;
        return z;
    }

    inline
    bool isnan(std::complex<double> const& z)
    {
        using std::isnan;
        return isnan(z.real()) || isnan(z.imag());
    }
 

    inline
    double RandReal()
    {
        // pinned draw (b2rand/1, ADR-0044)
        return records::DrawSymmetricDouble();
    }
 
    namespace{
        using real_dbl = double;
        using complex_dbl = std::complex<double>;
    }

    inline complex_dbl pow(const complex_dbl & z, int power)
    {
        if (power < 0) {
            return pow(1./z, -power);
        }
        else if (power==0)
            return complex_dbl(1,0);
        else if(power==1)
            return z;
        else if(power==2)
            return z*z;
        else if(power==3)
            return z*z*z;
        else
        {
            unsigned int p(static_cast<unsigned int>(power));
            complex_dbl result(1,0), z_to_the_current_power_of_two = z;
            // have copy of p in memory, can freely modify it.
            do {
                if ( (p & 1) == 1 ) { // get the lowest bit of the number
                    result *= z_to_the_current_power_of_two;
                }
                z_to_the_current_power_of_two *= z_to_the_current_power_of_two; // square z_to_the_current_power_of_two
            } while (p  >>= 1);
 
            return result;
        }
    }
 
} // namespace bertini
 
 
#endif // include guard