//  (C) Copyright John Maddock 2006.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#include <iostream>
#include <iomanip>
#include <sstream>
#include <cmath>
#include <assert.h>
#include <limits>

#include <boost/tr1/random.hpp>

template <class T>
T next_after(T v)
{
   int exponent;
   std::frexp(v, &exponent);
   return v + std::ldexp(T(1), exponent-std::numeric_limits<T>::digits);
}


template <class Real, class charT, class traits>
std::basic_ostream<charT, traits>& write_real(
   std::basic_ostream<charT, traits>& os,
   Real f)
{
   static const charT digits[] = 
   { '0', '1', '2', '3', '4', '5', '6', '7', '8', '9', 
   'a', 'b', 'c', 'd', 'e', 'f', '-', 'x', '.', 'p' };
   int exponent, digit;
   f = std::frexp(f, &exponent);
   // handle the sign:
   if(f < 0)
   {
      os.put(digits[16]);
      f = -f;
   }
   // Now the prefix:
   os.put(digits[0]).put(digits[17]);
   // extract first hex digit:
   f = ldexp(f, 4);
   digit = (int)std::floor(f);
   os.put(digits[digit]);
   f -= digit;
   if(f != 0)
   {
      // put the decimal point:
      os.put(digits[18]);
      do
      {
         // extract next hex digit:
         f = ldexp(f, 4);
         digit = (int)std::floor(f);
         os.put(digits[digit]);
         f -= digit;
      }while(f != 0);
   }
   // lastly the exponent:
   os.put(digits[19]);
   os << std::dec << exponent - 4;
   return os;
}

template <class charT, class traits>
std::basic_istream<charT, traits>& format_error(
   std::basic_istream<charT, traits>& is)
{
   is.setstate(std::ios::failbit);
   return is;
}

template <class charT>
int get_digit(charT c)
{
   if((c >= charT('0')) && (c <= charT('9')))
      return c - '0';
   if((c >= charT('a')) && (c <= charT('f')))
      return (c - 'a') + 10;
   if((c >= charT('A')) && (c <= charT('F')))
      return (c - 'a') + 10;
   return -1;
}

template <class Real, class charT, class traits>
std::basic_istream<charT, traits>& read_real(
   std::basic_istream<charT, traits>& is,
   Real& f)
{
   charT next_char;
   int sign = 1;
   Real value;
   int exponent(0), exponent_shift(0), digit;

   // skip white space:
   is >> std::ws;

   // deal with the sign if we have one:
   next_char = is.peek();
   if(next_char == charT('-'))
   {
      is.get();
      sign = -sign;
   }
   // Now we must have an "0x" prefix or it's an error:
   if(is.get() != charT('0'))
      return format_error(is);
   if(is.get() != charT('x'))
      return format_error(is);
   // now get the leading digit:
   digit = get_digit(is.get());
   if(digit < 0)
      return format_error(is);
   value = digit;
   // try for an optional decimal point:
   if(is.peek() == charT('.'))
   {
      is.get();
      // read off the digits after the point:
      while((digit = get_digit(is.peek())) >= 0)
      {
         is.get();
         value = std::ldexp(value, 4) + digit;
         exponent_shift += 4;
      }
   }
   // we must have a "p" or the format is invalid:
   if(is.get() != charT('p'))
      return format_error(is);
   // OK, just the exponent to get now:
   is >> std::dec >> exponent;

   // whew, all we have to do now is reconstruct the result:
   f = std::ldexp(value, exponent - exponent_shift);
   if(sign < 0)
      f = -f;

   return is;
}

// Random tests of hex format reals:
template <class T>
void test_io_1(T)
{
   std::tr1::mt19937 rnd;
   std::tr1::uniform_real<T> ur_a(0.5, 1 /*-(std::numeric_limits<T>::max)(), (std::numeric_limits<T>::max)()*/);
   std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<T> > gen(rnd, ur_a);
   std::tr1::uniform_int<int> ui(std::numeric_limits<T>::min_exponent, std::numeric_limits<T>::max_exponent);
   std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_int<int> > geni(rnd, ui);

   T value = 0.0019075645054089487L;
   T read_value;
   for(int i = 0; i <= INT_MAX; ++i)
   {
      std::stringstream ss1, ss2;
      write_real(ss1, value);
      read_real(ss1, read_value);
      if(read_value != value)
      {
         std::cerr << "Bad read value found: " << write_real(std::cerr, value);
      }
      write_real(ss2, -value);
      read_real(ss2, read_value);
      if(read_value != -value)
      {
         std::cerr << "Bad read value found: " << write_real(std::cerr, value);
      }
      if(i && ((i % 1000) == 0))
         std::cout << "Tested " << i << " values...\n";
      value = std::ldexp(gen(), geni());
   }
}

// Sequential tests of hex format reals:
template <class T>
void test_io_2(T)
{
   T value = 0;
   T read_value;
   // this doesn't test every floating point value, but it does test a lot of them!
   for(int i = 0; i <= INT_MAX; ++i)
   {
      std::stringstream ss1, ss2;
      write_real(ss1, value);
      read_real(ss1, read_value);
      if(read_value != value)
      {
         std::cerr << "Bad read value found: " << write_real(std::cerr, value);
      }
      write_real(ss2, -value);
      read_real(ss2, read_value);
      if(read_value != -value)
      {
         std::cerr << "Bad read value found: " << write_real(std::cerr, value);
      }
      if(i && ((i % 1000) == 0))
         std::cout << "Tested " << i << " values...\n";
      value = next_after(value);
   }
}

// Random tests of native (decimal) format reals:
template <class T>
void test_native_io(T)
{
   std::tr1::mt19937 rnd;
   std::tr1::uniform_real<T> ur_a(0.5, 1 /*-(std::numeric_limits<T>::max)(), (std::numeric_limits<T>::max)()*/);
   std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<T> > gen(rnd, ur_a);
   std::tr1::uniform_int<int> ui(std::numeric_limits<T>::min_exponent, std::numeric_limits<T>::max_exponent);
   std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_int<int> > geni(rnd, ui);

   T value = 0.0019075645054089487L; //std::ldexp(gen(), geni());
   T read_value;
   for(int i = 0; i <= INT_MAX; ++i)
   {
      std::stringstream ss1, ss2;
      ss1 << std::setprecision(std::numeric_limits<T>::digits10+2);
      ss2 << std::setprecision(std::numeric_limits<T>::digits10+2);
      ss1 << value;
      ss1 >> read_value;
      if(read_value != value)
      {
         std::cerr << "Bad read value found: " << write_real(std::cerr, value);
      }
      ss2 << -value;
      ss2 >> read_value;
      if(read_value != -value)
      {
         std::cerr << "Bad read value found: " << write_real(std::cerr, value);
      }
      if(i && ((i % 1000) == 0))
         std::cout << "Tested " << i << " values...\n";
      value = std::ldexp(gen(), geni());
   }
}

// Sequential tests of native (decimal) format reals:
template <class T>
void test_native_io_2(T)
{
   T value = 1;
   T read_value;
   // this doesn't test every floating point value, but it does test a lot of them!
   for(int i = 0; i <= INT_MAX; ++i)
   {
      std::stringstream ss1, ss2;
      ss1 << std::setprecision(std::numeric_limits<T>::digits10+2);
      ss2 << std::setprecision(std::numeric_limits<T>::digits10+2);
      ss1 << value;
      ss1 >> read_value;
      if(read_value != value)
      {
         std::cerr << "Bad read value found: " << write_real(std::cerr, value);
      }
      ss2 << -value;
      ss2 >> read_value;
      if(read_value != -value)
      {
         std::cerr << "Bad read value found: " << write_real(std::cerr, value);
      }
      if(i && ((i % 1000) == 0))
         std::cout << "Tested " << i << " values...\n";
      value = next_after(value);
   }
}


int main()
{
   // Uncomment this test native io operators:
   //test_native_io((double)(0));
   //test_native_io_2((float)(0));

   // random tests, adjust real type as required:
   test_io_1((long double)(0));

   // sequential tests, adjust real type as required:
   test_io_2((long double)(0));
   return 0;
}