Some prime number research, providing the ability to select precision and see the tradeoffs on speed. Do everything to make the computation fast. Yes, I feel dirty about the sighandler stuff.

README.md

Installation

# let's install the dependencies
brew install cmake primesieve boost mpfr

# download this gist as a zip file and set it up
# unpack the zip file to this directory
# open your web browser to https://gist.github.com/pstoll/cf0f18a720c88baae8a865acca74ec14
# and "Download Zip" 
cp ~/Downloads/cf0f18a720c88baae8a865acca74ec14-*.zip ./cf0f18a720c88baae8a865acca74ec14.zip
unzip cf0f18a720c88baae8a865acca74ec14.zip
mkdir cf0f18a720c88baae8a865acca74ec14/* .
rm -rf cf0f18a720c88baae8a865acca74ec14/ cf0f18a720c88baae8a865acca74ec14.zip
mkdir cmake
mv cmake-FindMPFR.cmake cmake/FindMPFR.cmake

# get it ready to build
cmake .

# now actually build the executable
make

cmake-FindMPFR.cmake

# Try to find the MPFR library
# See http://www.mpfr.org/
#
# This module supports requiring a minimum version, e.g. you can do
#   find_package(MPFR 2.3.0)
# to require version 2.3.0 to newer of MPFR.
#
# Once done this will define
#
#  MPFR_FOUND - system has MPFR lib with correct version
#  MPFR_INCLUDES - the MPFR include directory
#  MPFR_LIBRARIES - the MPFR library
#  MPFR_VERSION - MPFR version

# Copyright (c) 2006, 2007 Montel Laurent, <montel@kde.org>
# Copyright (c) 2008, 2009 Gael Guennebaud, <g.gael@free.fr>
# Copyright (c) 2010 Jitse Niesen, <jitse@maths.leeds.ac.uk>
# Redistribution and use is allowed according to the terms of the BSD license.

# Set MPFR_INCLUDES

find_path(MPFR_INCLUDES
  NAMES
  mpfr.h
  PATHS
  $ENV{GMPDIR}
  ${INCLUDE_INSTALL_DIR}
)

# Set MPFR_FIND_VERSION to 1.0.0 if no minimum version is specified

if(NOT MPFR_FIND_VERSION)
  if(NOT MPFR_FIND_VERSION_MAJOR)
    set(MPFR_FIND_VERSION_MAJOR 1)
  endif(NOT MPFR_FIND_VERSION_MAJOR)
  if(NOT MPFR_FIND_VERSION_MINOR)
    set(MPFR_FIND_VERSION_MINOR 0)
  endif(NOT MPFR_FIND_VERSION_MINOR)
  if(NOT MPFR_FIND_VERSION_PATCH)
    set(MPFR_FIND_VERSION_PATCH 0)
  endif(NOT MPFR_FIND_VERSION_PATCH)

  set(MPFR_FIND_VERSION "${MPFR_FIND_VERSION_MAJOR}.${MPFR_FIND_VERSION_MINOR}.${MPFR_FIND_VERSION_PATCH}")
endif(NOT MPFR_FIND_VERSION)


if(MPFR_INCLUDES)

  # Set MPFR_VERSION
  
  file(READ "${MPFR_INCLUDES}/mpfr.h" _mpfr_version_header)
  
  string(REGEX MATCH "define[ \t]+MPFR_VERSION_MAJOR[ \t]+([0-9]+)" _mpfr_major_version_match "${_mpfr_version_header}")
  set(MPFR_MAJOR_VERSION "${CMAKE_MATCH_1}")
  string(REGEX MATCH "define[ \t]+MPFR_VERSION_MINOR[ \t]+([0-9]+)" _mpfr_minor_version_match "${_mpfr_version_header}")
  set(MPFR_MINOR_VERSION "${CMAKE_MATCH_1}")
  string(REGEX MATCH "define[ \t]+MPFR_VERSION_PATCHLEVEL[ \t]+([0-9]+)" _mpfr_patchlevel_version_match "${_mpfr_version_header}")
  set(MPFR_PATCHLEVEL_VERSION "${CMAKE_MATCH_1}")
  
  set(MPFR_VERSION ${MPFR_MAJOR_VERSION}.${MPFR_MINOR_VERSION}.${MPFR_PATCHLEVEL_VERSION})
  
  # Check whether found version exceeds minimum version
  
  if(${MPFR_VERSION} VERSION_LESS ${MPFR_FIND_VERSION})
    set(MPFR_VERSION_OK FALSE)
    message(STATUS "MPFR version ${MPFR_VERSION} found in ${MPFR_INCLUDES}, "
                   "but at least version ${MPFR_FIND_VERSION} is required")
  else(${MPFR_VERSION} VERSION_LESS ${MPFR_FIND_VERSION})
    set(MPFR_VERSION_OK TRUE)
  endif(${MPFR_VERSION} VERSION_LESS ${MPFR_FIND_VERSION})

endif(MPFR_INCLUDES)

# Set MPFR_LIBRARIES

find_library(MPFR_LIBRARIES mpfr PATHS $ENV{GMPDIR} ${LIB_INSTALL_DIR})

# Epilogue

include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(MPFR DEFAULT_MSG
                                  MPFR_INCLUDES MPFR_LIBRARIES MPFR_VERSION_OK)
mark_as_advanced(MPFR_INCLUDES MPFR_LIBRARIES)

CMakeLists.txt

cmake_minimum_required(VERSION 3.5)
set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/cmake)

find_package(primesieve "7.0" REQUIRED )

#if (${CMAKE_SYSTEM_NAME} MATCHES "Emscripten")
#    set(CMAKE_C_COMPILER "emcc")
#    set(CMAKE_CXX_COMPILER "emcc")
#endif ()

if(USE_QUAD_FLOAT)
  add_library(quadmath STATIC IMPORTED)
  find_library(quadmath_lib NAMES quadmath.a)
endif()

find_package(MPFR REQUIRED)
include_directories(${MPFR_INCLUDES})
if(NOT MPFR_FOUND)
    message(FATAL_ERROR "Could not find MPFR!")
endif(NOT MPFR_FOUND)

find_package(Boost 1.65.0)
if(Boost_FOUND)
  include_directories(${Boost_INCLUDE_DIRS})
endif()
set(Boost_USE_STATIC_LIBS        ON) # only find static libs
set(Boost_USE_MULTITHREADED      ON)
set(Boost_USE_STATIC_RUNTIME    ON)


add_executable (hex_sums hex_sums.cpp)
target_link_libraries(hex_sums primesieve::primesieve ${QUAD_FLOAT_LIBS} ${MPFR_LIBRARIES})

if (${CMAKE_SYSTEM_NAME} MATCHES "Emscripten")
    set_target_properties(client PROPERTIES LINK_FLAGS "-o hexsum.js -s WASM=1 -s DEMANGLE_SUPPORT=1  --bind")


set_target_properties(hex_sums PROPERTIES LINK_FLAGS "-o hexsum.js -s WASM=1 -O2 -std=c++11")
else ()
    #target_link_libraries(hex_sums)
endif ()

hex_sums.cpp

//
// cmake -DCMAKE_C_COMPILER=/usr/local/bin/gcc-7 -DCMAKE_CXX_COMPILER=/usr/local/bin/g++-7  -DBUILD_EXAMPLES=ON ..
// 
// using prime iterator to do prime hexagon calculations
// 

#include <primesieve.hpp>
#include <iostream>

#if USE_QUAD_FLOAT
#include <boost/multiprecision/float128.hpp>
#endif
#include <boost/multiprecision/mpfr.hpp>
#include <boost/math/constants/constants.hpp>

namespace bmp = boost::multiprecision;

#include <chrono>
typedef std::chrono::high_resolution_clock Clock;
typedef uint64_t prime_value_t;

#include <string>
#include <cstdlib>
#include <cxxabi.h>
#include <map>

#include <iostream>
#include <csignal>

static const char my_version[] = "1.0.2";

template<typename T>
std::string type_name()
{
    int status;
    std::string tname = typeid(T).name();
    char *demangled_name = abi::__cxa_demangle(tname.c_str(), NULL, NULL, &status);
    if(status == 0) {
        tname = demangled_name;
        std::free(demangled_name);
    }   
    return tname;
}

#define HEXDEBUG(x) if (0) do x while(0)


// type is to hold a sum of prime fractions, which we are varying precision on
template <class TS, class TV=uint64_t> // template class sum (TS) and template class value (TV)
class prime_result_t {
 public:
  TV start_, end_, cur_num_;
  int32_t start_sign_;
  TS start_int_sum_, start_prime_sum_;
  TS neg_int_sum_, pos_int_sum_;
  TS neg_prime_sum_, pos_prime_sum_;
  TS sum_exponent_;
  
  prime_result_t( TV start, TV end, int32_t start_sign, TS sum_exponent):
    start_(start), end_(end), cur_num_(0),start_sign_(start_sign),
    start_int_sum_(0), start_prime_sum_(0), sum_exponent_(sum_exponent),
    neg_int_sum_(0), pos_int_sum_(0), neg_prime_sum_(0), pos_prime_sum_(0)
  {
    instance_ = this;
    init_sums();
  };

  prime_result_t( TV end, TS sum_exponent):
    start_(4), end_(end), cur_num_(0),start_sign_(1),
    start_int_sum_(0), start_prime_sum_(0), sum_exponent_(sum_exponent),
    neg_int_sum_(0), pos_int_sum_(0), neg_prime_sum_(0), pos_prime_sum_(0)
  {
    instance_ = this;
    init_sums();
  };

  void init_sums() {
    const TS one   = 1,
      two   = 2,
      three = 3 ,
      four  = 4;

    const TS exp = sum_exponent_;
    start_int_sum_ = one + pow(one/two, exp) + pow(one/three, exp) + pow(one/four, exp);
    start_prime_sum_ = pow(one/two, exp) + pow(one/three, exp);
  }
  
  const void report_computation()
  {
    const TS total_int_sum   = start_int_sum_   + pos_int_sum_   - neg_int_sum_;
    const TS total_prime_sum = start_prime_sum_ + pos_prime_sum_ - neg_prime_sum_;
    const TS total_ratio = total_int_sum / total_prime_sum;
    
    //std::cout.precision(std::numeric_limits<T>::max_digits10);
    std::cout.precision(50);

    long precision = std::numeric_limits<TS>::digits10;
    long bytes = sizeof(TS);
    // dynamic types have zero here. so ask it...
    if (precision == 0) {
      long bits = mpfr_get_default_prec();
      long digits = bmp::detail::digits2_2_10(bits);
      precision = digits;
      bytes = mpfr_custom_get_size(bits);
    }
    
    // type_name<TS>() 
    std::cout << "Sum Range primes " << "(" << start_ << "-" << end_ << ", @" << cur_num_ << ")"
	      << " with sum_exponent= " << sum_exponent_ << std::endl
	      << " using floating point with " << precision << " digits"
	      << " and " << bytes << " bytes" << std::endl
	      << "Net ∑ Integers: " << total_int_sum << std::endl
	      << "Net ∑ Primes  : " << total_prime_sum << std::endl
	      << "Ratio         : " << total_ratio << std::endl
      
	      << "Neg ∑ Ints    : " << neg_int_sum_ << std::endl
	      << "Pos ∑ Ints    : " << (start_int_sum_ + pos_int_sum_) << std::endl
	      << "Delta ∑ Ints  : " << (start_int_sum_ + pos_int_sum_ - neg_int_sum_) << std::endl
	      << "Pos ∑/Neg ∑   : " << (start_int_sum_ + pos_int_sum_)/(neg_int_sum_) << std::endl
      
	      << "Neg ∑ Prime   : " << neg_prime_sum_ << std::endl
	      << "Pos ∑ Prime   : " << (start_prime_sum_ + pos_prime_sum_) << std::endl
	      << "Delta ∑ Prime : " << (start_prime_sum_ + pos_prime_sum_ - neg_prime_sum_) << std::endl
	      << "Pos ∑/Neg ∑   : " << (start_prime_sum_ + pos_prime_sum_)/(neg_prime_sum_) << std::endl;
  }

  void sig_handler(int signum) {
    std::cout << "Instance signal called for (" << signum << "). \n";
    report_computation();
    
  }

  // static void static_sig_handler(int signum) {
  //   std::cout << "Static signal handler called for(" << signum << "). \n";
  //   instance_->sig_handler(signum);
  // }

  prime_result_t<TS,TV> *instance_;
};

template <class T> int compute_prime_sums(uint64_t start, uint64_t end, int32_t start_sign, prime_result_t<T>& res );

typedef double my_double;
#ifdef USE_QUAD_FLOAT
typedef bmp::number<bmp::arithmetic_backend<bmp::float128>, bmp::et_off> my_quad;
typedef bmp::float128 my_float128;
#endif
typedef bmp::number<bmp::mpfr_float_backend<64, bmp::allocate_stack> > my_f64;
typedef bmp::number<bmp::mpfr_float_backend<128, bmp::allocate_stack> > my_f128;
typedef bmp::number<bmp::mpfr_float_backend<256, bmp::allocate_stack> > my_f256;
typedef bmp::number<bmp::mpfr_float_backend<512, bmp::allocate_stack> > my_f512;
typedef bmp::number<bmp::mpfr_float_backend<1024, bmp::allocate_stack> > my_f1024;
typedef bmp::mpfr_float my_f_any;

template<typename T> void my_signal_wrapper(int signum) {};

// Wrapper function uses a global to remember the object:
prime_result_t<double>* res_double_ptr;
template<> void my_signal_wrapper<double>(int signum)
{
  if (res_double_ptr) { res_double_ptr->sig_handler(signum); }
}

#ifdef USE_QUAD_FLOAT
prime_result_t<my_float128>* res_float128_ptr;
template<> void my_signal_wrapper<my_float128>(int signum)
{
  if (res_float128_ptr) { res_float128_ptr->sig_handler(signum); }
}
#endif

prime_result_t<my_f64>* res_f64_ptr;
template<> void my_signal_wrapper<my_f64>(int signum)
{
  if (res_f64_ptr) { res_f64_ptr->sig_handler(signum); }
}

prime_result_t<my_f128>* res_f128_ptr;
template<> void my_signal_wrapper<my_f128>(int signum)
{
  if (res_f128_ptr) { res_f128_ptr->sig_handler(signum); }
}

prime_result_t<my_f256>* res_f256_ptr;
template<> void my_signal_wrapper<my_f256>(int signum)
{
  if (res_f256_ptr) { res_f256_ptr->sig_handler(signum); }
}

prime_result_t<my_f512>* res_f512_ptr;
template<> void my_signal_wrapper<my_f512>(int signum)
{
  if (res_f512_ptr) { res_f512_ptr->sig_handler(signum); }
}

prime_result_t<my_f1024>* res_f1024_ptr;
template<> void my_signal_wrapper<my_f1024>(int signum)
{
  if (res_f1024_ptr) { res_f1024_ptr->sig_handler(signum); }
}

prime_result_t<my_f_any>* res_f_any_ptr;
template<> void my_signal_wrapper<my_f_any>(int signum)
{
  if (res_f_any_ptr) { res_f_any_ptr->sig_handler(signum); }
}

template <typename T>
void time_one_sum(prime_result_t<T> &res)
{
  auto t1 = Clock::now();
  compute_prime_sums( res.start_, res.end_, res.start_sign_, res.sum_exponent_, res );
  auto t2 = Clock::now();

  res.report_computation();
  std::cout << " (computed in "
	    << std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count()
	    << " ms)" << std::endl;
  
}


template <class T>
int compute_prime_sums(const uint64_t prime_start, const uint64_t end, const int32_t start_sign, const T sum_exponent, prime_result_t<T>& res ) {
    
  primesieve::iterator it(prime_start); // Generate primes > start
  bool positive = (start_sign > 0);
  
  T neg_int_sum = 0.0,
    pos_int_sum = 0.0,
    *int_sum_ptr;

  T neg_prime_sum = 0.0,
    pos_prime_sum = 0.0,
    *prime_sum_ptr;
  
  prime_value_t cur_prime, next_prime;
  prime_value_t ii = 0;

  prime_value_t start = it.next_prime();

  HEXDEBUG( std::cout << "start_prime = " << prime_start
	    << ", start=" << start
	    << ", end=" << end << std::endl; );
  
  // iterate over the primes
  for (cur_prime = start, next_prime = it.next_prime(); next_prime < end; next_prime = it.next_prime()) {
    
    HEXDEBUG( std::cout << "cur_prime = " << cur_prime << std::endl; );

    positive = !positive;	// switch flag to other sign
    
    prime_sum_ptr = (positive ? &res.pos_prime_sum_ : &res.neg_prime_sum_);

    {
      T v1 = (T)1/(T)cur_prime;
      T v2 = pow(v1, sum_exponent);
      *prime_sum_ptr += v2;
    }
    //res.report_computation();
    
    int_sum_ptr = (positive ? &res.pos_int_sum_ : &res.neg_int_sum_);
    // iterate over ints from current prime to next prime
    for (ii = cur_prime; ii < next_prime; ++ii) {
      res.cur_num_ = ii;
      T v1 = (T)1/(T)ii;
      T v2 = pow(v1, sum_exponent);
      *int_sum_ptr += v2;
      HEXDEBUG( std::cout << "cur_int = " << ii << std::endl; );
      //res.report_computation();
    }
    
    // save prime for next loop
    cur_prime = next_prime;
  }

  // these values are computed and stored in the loops above...
  //res.neg_int_sum_ = neg_int_sum;
  //res.pos_int_sum_ = pos_int_sum;

  //res.neg_prime_sum_ = neg_prime_sum;
  //res.pos_prime_sum_ = pos_prime_sum;

  return 0;
}

std::map<std::string, my_f_any> float_const_map;


void usage(char** argv, std::string message){
  if (message.length() > 0) {
    std::cerr << message << std::endl;
  }
  std::cerr << "usage: " << argv[0] << " prime_range_power_of_10 sum_exponent float_digits" << std::endl;
  std::cerr << "  version: " <<  my_version << std::endl;
  std::cerr << "  example: " << argv[0] << " 8 1.5 256" << std::endl;
  std::cerr << "  Any number of float digit precision is supported but optimized version are available for:" << std::endl;
  std::cerr << "    double: "
	    << std::numeric_limits<double>::digits10 << " digits of precision"
	    << " using " << sizeof(double) << " bytes" << std::endl;
#ifdef USE_QUAD_FLOAT
  std::cerr << "    quad-float: "
	    << std::numeric_limits<my_float128>::digits10 << " digits of precision"
	    << " using " << sizeof(my_float128) << " bytes" << std::endl;
#endif
  std::cerr << "    mpfr 64: "
	    << std::numeric_limits<my_f64>::digits10 << " digits of precision"
	    << " using " << sizeof(my_f64) << " bytes" << std::endl;
  std::cerr << "    mpfr 128: "
	    << std::numeric_limits<my_f128>::digits10 << " digits of precision"
	    << " using " << sizeof(my_f128) << " bytes" << std::endl;
  std::cerr << "    mpfr 256: "
	    << std::numeric_limits<my_f256>::digits10 << " digits of precision"
	    << " using " << sizeof(my_f256) << " bytes" << std::endl;
  std::cerr << "    mpfr 512: "
	    << std::numeric_limits<my_f512>::digits10 << " digits of precision"
	    << " using " << sizeof(my_f512) << " bytes" << std::endl;
  std::cerr << "    mpfr 1024: "
	    << std::numeric_limits<my_f1024>::digits10 << " digits of precision"
	    << " using " << sizeof(my_f1024) << " bytes" << std::endl;
  std::cerr << "  Known constants for sum_exponet are:" << std::endl;
  for(auto ii = float_const_map.begin(); ii != float_const_map.end(); ++ii) {
    std::cerr << "     " << ii->first << " = " << ii->second << std::endl;
  }
  exit(1);
}

void build_map(std::map<std::string, my_f_any> & m) {
  m["half"] = boost::math::constants::half<my_f_any>();
  m["third"] = boost::math::constants::third<my_f_any>();
  m["two_thirds"] = boost::math::constants::two_thirds<my_f_any>();
  m["three_quarters"] = boost::math::constants::three_quarters<my_f_any>();
  m["root_two"] = boost::math::constants::root_two<my_f_any>();
  m["root_three"] = boost::math::constants::root_three<my_f_any>();
  m["half_root_two"] = boost::math::constants::half_root_two<my_f_any>();
  m["ln_two"] = boost::math::constants::ln_two<my_f_any>();
  m["ln_ten"] = boost::math::constants::ln_ten<my_f_any>();
  m["ln_ln_two"] = boost::math::constants::ln_ln_two<my_f_any>();
  m["root_ln_four"] = boost::math::constants::root_ln_four<my_f_any>();
  m["one_div_root_two"] = boost::math::constants::one_div_root_two<my_f_any>();
  m["pi"] = boost::math::constants::pi<my_f_any>();
  m["half_pi"] = boost::math::constants::half_pi<my_f_any>();
  m["third_pi"] = boost::math::constants::third_pi<my_f_any>();
  m["sixth_pi"] = boost::math::constants::sixth_pi<my_f_any>();
  m["two_pi"] = boost::math::constants::two_pi<my_f_any>();
  m["two_thirds_pi"] = boost::math::constants::two_thirds_pi<my_f_any>();
  m["three_quarters_pi"] = boost::math::constants::three_quarters_pi<my_f_any>();
  m["four_thirds_pi"] = boost::math::constants::four_thirds_pi<my_f_any>();
  m["one_div_two_pi"] = boost::math::constants::one_div_two_pi<my_f_any>();
  m["root_pi"] = boost::math::constants::root_pi<my_f_any>();
  m["root_half_pi"] = boost::math::constants::root_half_pi<my_f_any>();
  m["root_two_pi"] = boost::math::constants::root_two_pi<my_f_any>();
  m["one_div_root_pi"] = boost::math::constants::one_div_root_pi<my_f_any>();
  m["one_div_root_two_pi"] = boost::math::constants::one_div_root_two_pi<my_f_any>();
  m["root_one_div_pi"] = boost::math::constants::root_one_div_pi<my_f_any>();
  m["pi_minus_three"] = boost::math::constants::pi_minus_three<my_f_any>();
  m["four_minus_pi"] = boost::math::constants::four_minus_pi<my_f_any>();
  m["pi_pow_e"] = boost::math::constants::pi_pow_e<my_f_any>();
  m["pi_sqr"] = boost::math::constants::pi_sqr<my_f_any>();
  m["pi_sqr_div_six"] = boost::math::constants::pi_sqr_div_six<my_f_any>();
  m["pi_cubed"] = boost::math::constants::pi_cubed<my_f_any>();
  m["cbrt_pi"] = boost::math::constants::cbrt_pi<my_f_any>();
  m["one_div_cbrt_pi"] = boost::math::constants::one_div_cbrt_pi<my_f_any>();
  m["e"] = boost::math::constants::e<my_f_any>();
  m["exp_minus_half"] = boost::math::constants::exp_minus_half<my_f_any>();
  m["e_pow_pi"] = boost::math::constants::e_pow_pi<my_f_any>();
  m["root_e"] = boost::math::constants::root_e<my_f_any>();
  m["log10_e"] = boost::math::constants::log10_e<my_f_any>();
  m["one_div_log10_e"] = boost::math::constants::one_div_log10_e<my_f_any>();
  m["degree"] = boost::math::constants::degree<my_f_any>();
  m["radian"] = boost::math::constants::radian<my_f_any>();
  m["sin_one"] = boost::math::constants::sin_one<my_f_any>();
  m["cos_one"] = boost::math::constants::cos_one<my_f_any>();
  m["sinh_one"] = boost::math::constants::sinh_one<my_f_any>();
  m["cosh_one"] = boost::math::constants::cosh_one<my_f_any>();
  m["phi"] = boost::math::constants::phi<my_f_any>();
  m["ln_phi"] = boost::math::constants::ln_phi<my_f_any>();
  m["one_div_ln_phi"] = boost::math::constants::one_div_ln_phi<my_f_any>();
  m["euler"] = boost::math::constants::euler<my_f_any>();
  m["one_div_euler"] = boost::math::constants::one_div_euler<my_f_any>();
  m["euler_sqr"] = boost::math::constants::euler_sqr<my_f_any>();
  m["zeta_two"] = boost::math::constants::zeta_two<my_f_any>();
  m["zeta_three"] = boost::math::constants::zeta_three<my_f_any>();
  m["catalan"] = boost::math::constants::catalan<my_f_any>();
  m["glaisher"] = boost::math::constants::glaisher<my_f_any>();
#if 0
  // holy cow - this takes *forever* to run at runtime. you are not welcome, sir khinchin <cough> constant.
  m["khinchin"] = boost::math::constants::khinchin<my_f_any>();
#endif
  m["extreme_value_skewness"] = boost::math::constants::extreme_value_skewness<my_f_any>();
  m["rayleigh_skewness"] = boost::math::constants::rayleigh_skewness<my_f_any>();
  m["rayleigh_kurtosis_excess"] = boost::math::constants::rayleigh_kurtosis_excess<my_f_any>();
  m["rayleigh_kurtosis"] = boost::math::constants::rayleigh_kurtosis<my_f_any>();
  
}

template <class T>
T get_float_const_value(std::map<std::string, my_f_any> & m, std::string s){
  T val;
  auto ii = m.find(s);
  if (ii != m.end()) {
    val =  static_cast<T>(ii->second);
  } else {
    val = static_cast<T>(s);
    if (isnan(val)) {
      std::cerr << "could not convert " << s << " into a number" << std::endl;
    }
  }
  return val;
}

int main(int argc, char** argv)
{
  double max_prime_power = 10;
  int float_digits = 128;
  std::string sum_exp_str;

  build_map(float_const_map);
  
  if (argc < 3) {
    usage(argv, "");
  }
  
  if (argc > 1) {
    max_prime_power = std::atof(argv[1]);
    if (isnan(max_prime_power)) {
      usage(argv, "Could not convert power to a number");
    }
  }

  if (argc > 2) {
    sum_exp_str = argv[2];
    //double sum_exponent = 1.0;
    //sum_exponent = std::atof(argv[2]);
    //if (isnan(sum_exponent)){ 
    //  usage(argv, "Could not convert power to an floating point number");
    //}

  }

  if (argc > 3) {
    float_digits = std::atoi(argv[3]);
    if (float_digits == 0) {
      usage(argv, "Could not convert FLOAT_DIGITS input to an int");
    }
    
  }

  const uint64_t nn = pow(10, max_prime_power);

  long bits = bmp::detail::digits10_2_2((float_digits));
  //std::cout << "converting " << float_digits << " floating point digits to " << bits << " bits " << std::endl;

  switch (float_digits) {
  case std::numeric_limits<double>::digits10:
    {
      std::cout << "using optimized double data type..." << std::endl;
      my_f_any v = get_float_const_value<my_f_any>(float_const_map, sum_exp_str);
      double sum_exp = static_cast<double>(v);
      // std::atof(sum_exp_str.c_str());
      prime_result_t<double> res_double(nn, sum_exp);
      res_double_ptr = &res_double;
      signal(SIGINT, &my_signal_wrapper<double> );  
      time_one_sum<double>(res_double);
    }
    break;
#ifdef USE_QUAD_FLOAT
  case std::numeric_limits<my_float128>::digits10:
    {
      my_float128 sum_exp = get_float_const_value<my_float128>(float_const_map, sum_exp_str);
      //my_float128 sum_exp = static_cast<my_float128>(sum_exp_str);
      prime_result_t<my_float128> res_float128(nn, sum_exp);
      res_float128_ptr = &res_float128;
      signal(SIGINT, &my_signal_wrapper<my_float128> );  
      time_one_sum<my_float128>(float128);
    }
    break;
#endif
  case std::numeric_limits<my_f256>::digits10:
    {
      my_f256 sum_exp = get_float_const_value<my_f256>(float_const_map, sum_exp_str);
      //my_f256 sum_exp = static_cast<my_f256>(sum_exp_str);
      prime_result_t<my_f256> res_f256(nn, sum_exp);
      res_f256_ptr = &res_f256;
      signal(SIGINT, &my_signal_wrapper<my_f256> );  
      time_one_sum<my_f256>(res_f256);
    }
    break;
  case std::numeric_limits<my_f64>::digits10:
    {
      my_f64 sum_exp = get_float_const_value<my_f64>(float_const_map, sum_exp_str);
      //my_f64 sum_exp = static_cast<my_f64>(sum_exp_str);
      prime_result_t<my_f64> res_f64(nn, sum_exp);
      res_f64_ptr = &res_f64;
      signal(SIGINT, &my_signal_wrapper<my_f64> );  
      time_one_sum<my_f64>(res_f64);
    }
    break;
  case std::numeric_limits<my_f128>::digits10:
    {
      my_f128 sum_exp = get_float_const_value<my_f128>(float_const_map, sum_exp_str);
      //my_f128 sum_exp = static_cast<my_f128>(sum_exp_str);
      prime_result_t<my_f128> res_f128(nn, sum_exp);
      res_f128_ptr = &res_f128;
      signal(SIGINT, &my_signal_wrapper<my_f128> );  
      time_one_sum<my_f128>(res_f128);
    }
    break;
  case std::numeric_limits<my_f512>::digits10:
    {
      my_f512 sum_exp = get_float_const_value<my_f512>(float_const_map, sum_exp_str);
      prime_result_t<my_f512> res_f512(nn, sum_exp);
      res_f512_ptr = &res_f512;
      signal(SIGINT, &my_signal_wrapper<my_f512> );  
      time_one_sum<my_f512>(res_f512);
    }
    break;
  case std::numeric_limits<my_f1024>::digits10:
    { 
      my_f1024 sum_exp = get_float_const_value<my_f1024>(float_const_map, sum_exp_str);
      prime_result_t<my_f1024> res_f1024(nn, sum_exp);
      res_f1024_ptr = &res_f1024;
      signal(SIGINT, &my_signal_wrapper<my_f1024> );  
      time_one_sum<my_f1024>(res_f1024);
    }
    break;
  default:
    {
      my_f_any sum_exp = get_float_const_value<my_f_any>(float_const_map, sum_exp_str);
      prime_result_t<my_f_any> res_f_any(nn, sum_exp);
      mpfr_set_default_prec(bits);
      bmp::mpfr_float::default_precision(float_digits);
      res_f_any_ptr = &res_f_any;
      signal(SIGINT, &my_signal_wrapper<my_f_any> );
      time_one_sum<my_f_any>(res_f_any);
    }
    break;

  }
  
  return 0;
}