Subject: [Boost-commit] svn:boost r73708 - in sandbox/e_float: boost/e_float libs/e_float/src/e_float/efx libs/e_float/src/e_float/gmp libs/e_float/src/e_float/mpfr libs/e_float/src/functions/elementary libs/e_float/src/functions/integer libs/e_float/src/utility libs/e_float/test
From: e_float_at_[hidden]
Date: 2011-08-12 15:50:55
Author: christopher_kormanyos
Date: 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
New Revision: 73708
URL: http://svn.boost.org/trac/boost/changeset/73708
Log:
- Included guard digits in the definition of max_digits10. (Moving toward C++ conformance.)
- General code cleanup.
- Added (but no build for) Linpack benchmark and F2C header.
Added:
sandbox/e_float/libs/e_float/test/f2c.h (contents, props changed)
sandbox/e_float/libs/e_float/test/linpack-benchmark.cpp (contents, props changed)
Text files modified:
sandbox/e_float/boost/e_float/e_float_base.hpp | 16 -
sandbox/e_float/boost/e_float/e_float_efx.hpp | 12
sandbox/e_float/boost/e_float/e_float_global_math.hpp | 2
sandbox/e_float/boost/e_float/e_float_gmp.hpp | 2
sandbox/e_float/boost/e_float/e_float_limits.hpp | 14
sandbox/e_float/boost/e_float/e_float_mpfr.hpp | 2
sandbox/e_float/libs/e_float/src/e_float/efx/e_float_efx.cpp | 321 ++++++++++++++++-----------------------
sandbox/e_float/libs/e_float/src/e_float/gmp/e_float_gmp.cpp | 82 ++++------
sandbox/e_float/libs/e_float/src/e_float/mpfr/e_float_mpfr.cpp | 38 +---
sandbox/e_float/libs/e_float/src/functions/elementary/elementary_math.cpp | 3
sandbox/e_float/libs/e_float/src/functions/integer/prime_factors.h | 10 +
sandbox/e_float/libs/e_float/src/utility/util_numeric_cast.h | 9 +
12 files changed, 217 insertions(+), 294 deletions(-)
Modified: sandbox/e_float/boost/e_float/e_float_base.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_base.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_base.hpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -29,8 +29,8 @@
// The supported range is 30-300.
// Note: This is a compile-time constant.
- #define E_FLOAT_DIGITS10 100
- #define E_FLOAT_DIGITS10_MAX 300
+ #define E_FLOAT_DIGITS10 100
+ #define E_FLOAT_DIGITS10_LIMIT 300
#if defined(E_FLOAT_TYPE_EFX)
namespace efx { class e_float; }
@@ -48,17 +48,11 @@
class e_float_base
{
public:
-
- // The value of ef_digits10_setting is the desired number of decimal digits
- // of precision in the e_float implementation. It is limited to the range of
- // 30...300 decimal digits. The digit tolerance is invariant and it is set
- // to be 15 percent (or maximally 150 digits) larger than ef_digits10.
- // All of these quantities are invariant and they are set at compile time.
static const INT32 ef_digits10_setting = E_FLOAT_DIGITS10;
- static const INT32 ef_digits10_max = E_FLOAT_DIGITS10_MAX;
- static const INT32 ef_digits10 = ((ef_digits10_setting < 30) ? 30 : ((ef_digits10_setting > ef_digits10_max) ? ef_digits10_max : ef_digits10_setting));
+ static const INT32 ef_digits10_limit = E_FLOAT_DIGITS10_LIMIT;
+ static const INT32 ef_digits10 = ((ef_digits10_setting < static_cast<INT32>(30)) ? static_cast<INT32>(30) : ((ef_digits10_setting > ef_digits10_limit) ? ef_digits10_limit : ef_digits10_setting));
static const INT32 ef_digits10_extra = static_cast<INT32>(((static_cast<INT64>(ef_digits10) * 15LL) + 50LL) / 100LL);
- static const INT32 ef_digits10_tol = static_cast<INT32>(ef_digits10 + ((ef_digits10_extra < 15) ? 15 : ((ef_digits10_extra > 150) ? 150 : ef_digits10_extra)));
+ static const INT32 ef_max_digits10 = static_cast<INT32>(ef_digits10 + ((ef_digits10_extra < static_cast<INT32>(5)) ? static_cast<INT32>(5) : ((ef_digits10_extra > static_cast<INT32>(30)) ? static_cast<INT32>(30) : ef_digits10_extra)));
virtual ~e_float_base() { }
Modified: sandbox/e_float/boost/e_float/e_float_efx.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_efx.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_efx.hpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -25,7 +25,6 @@
class e_float : public ::e_float_base
{
public:
- static const INT32 ef_elem_digits10 = static_cast<INT32>(8);
static const INT32 ef_radix = static_cast<INT32>(10);
static const INT32 ef_digits = ef_digits10;
@@ -34,8 +33,10 @@
static const INT64 ef_max_exp10 = static_cast<INT64>(+3063937869882635616LL); // Approx. [ef_max_exp / log10(2)], also an even multiple of 8
static const INT64 ef_min_exp10 = static_cast<INT64>(-3063937869882635616LL);
+ static const INT32 ef_elem_digits10 = static_cast<INT32>(8);
+
private:
- static const INT32 ef_digits10_num_base = static_cast<INT32>((ef_digits10_tol / ef_elem_digits10) + (((ef_digits10_tol % ef_elem_digits10) != 0) ? 1 : 0));
+ static const INT32 ef_digits10_num_base = static_cast<INT32>((ef_max_digits10 / ef_elem_digits10) + (((ef_max_digits10 % ef_elem_digits10) != 0) ? 1 : 0));
static const INT32 ef_elem_number = static_cast<INT32>(ef_digits10_num_base + 2);
typedef enum enum_fpclass
@@ -118,7 +119,7 @@
// Elementary primitives.
virtual e_float& calculate_inv (void);
virtual e_float& calculate_sqrt(void);
- virtual e_float& negate(void) { if(!iszero()) { neg = !neg; } return *this; }
+ virtual e_float& negate(void) { if(!iszero()) { neg = (!neg); } return *this; }
// Comparison functions
virtual bool isnan (void) const { return (fpclass == ef_NaN); }
@@ -144,7 +145,8 @@
virtual e_float extract_decimal_part (void) const;
private:
- static bool data_elem_is_nonzero_predicate(const UINT32& d) { return (d != static_cast<UINT32>(0u)); }
+ static bool data_elem_is_non_zero_predicate(const UINT32& d) { return (d != static_cast<UINT32>(0u)); }
+ static bool data_elem_is_non_nine_predicate(const UINT32& d) { return (d != static_cast<UINT32>(e_float::ef_elem_mask - 1)); }
void from_unsigned_long_long(const unsigned long long u);
void from_unsigned_long(const unsigned long u);
@@ -155,7 +157,7 @@
static UINT32 mul_loop_n (UINT32* const u, UINT32 n, const INT32 p);
static UINT32 div_loop_n (UINT32* const u, UINT32 n, const INT32 p);
- virtual INT64 get_order_exact(void) const;
+ virtual INT64 get_order_exact(void) const { return get_order_fast(); }
virtual INT64 get_order_fast(void) const;
virtual void get_output_string(std::string& str, INT64& my_exp, const std::size_t number_of_digits) const;
Modified: sandbox/e_float/boost/e_float/e_float_global_math.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_global_math.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_global_math.hpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -13,7 +13,7 @@
namespace ef
{
- inline INT64 tol(void) { return static_cast<INT64>(e_float::ef_digits10_tol); }
+ inline INT64 tol(void) { return static_cast<INT64>(e_float::ef_max_digits10); }
inline e_float fabs(const e_float& x) { return (x.isneg() ? e_float(x).negate() : x); }
}
Modified: sandbox/e_float/boost/e_float/e_float_gmp.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_gmp.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_gmp.hpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -66,7 +66,7 @@
static const INT64 ef_min_exp10 = static_cast<INT64>(-646456990LL);
private:
- static const INT32 ef_digits2 = static_cast<INT32>(((static_cast<INT64>(ef_digits10_tol) * 3322LL) + 500LL) / 1000LL);
+ static const INT32 ef_digits2 = static_cast<INT32>(((static_cast<INT64>(ef_max_digits10) * 3322LL) + 500LL) / 1000LL);
typedef enum enum_fpclass
{
Modified: sandbox/e_float/boost/e_float/e_float_limits.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_limits.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_limits.hpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -38,13 +38,13 @@
static const bool is_bounded = true;
static const bool is_modulo = false;
static const bool is_iec559 = false;
- static const INT32 digits = e_float::ef_digits; // Differs from int.
- static const INT32 digits10 = e_float::ef_digits10; // Differs from int.
- static const INT32 max_digits10 = e_float::ef_digits10 + 1; // Differs from int.
- static const INT64 min_exponent = e_float::ef_min_exp; // Differs from int.
- static const INT64 min_exponent10 = e_float::ef_min_exp10; // Differs from int.
- static const INT64 max_exponent = e_float::ef_max_exp; // Differs from int.
- static const INT64 max_exponent10 = e_float::ef_max_exp10; // Differs from int.
+ static const int digits = e_float::ef_digits;
+ static const int digits10 = e_float::ef_digits10;
+ static const int max_digits10 = e_float::ef_max_digits10;
+ static const INT64 min_exponent = e_float::ef_min_exp; // Type differs from int.
+ static const INT64 min_exponent10 = e_float::ef_min_exp10; // Type differs from int.
+ static const INT64 max_exponent = e_float::ef_max_exp; // Type differs from int.
+ static const INT64 max_exponent10 = e_float::ef_max_exp10; // Type differs from int.
static const int radix = e_float::ef_radix;
static const std::float_round_style round_style = std::round_to_nearest;
static const bool has_infinity = true;
Modified: sandbox/e_float/boost/e_float/e_float_mpfr.hpp
==============================================================================
--- sandbox/e_float/boost/e_float/e_float_mpfr.hpp (original)
+++ sandbox/e_float/boost/e_float/e_float_mpfr.hpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -74,7 +74,7 @@
static const INT64 ef_min_exp10 = static_cast<INT64>(-323228496LL);
private:
- static const INT32 ef_digits2 = static_cast<INT32>(((static_cast<INT64>(ef_digits10_tol) * 3322LL) + 500LL) / 1000LL);
+ static const INT32 ef_digits2 = static_cast<INT32>(((static_cast<INT64>(ef_max_digits10) * 3322LL) + 500LL) / 1000LL);
::mpfr_t rop;
public:
Modified: sandbox/e_float/libs/e_float/src/e_float/efx/e_float_efx.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/e_float/efx/e_float_efx.cpp (original)
+++ sandbox/e_float/libs/e_float/src/e_float/efx/e_float_efx.cpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -166,8 +166,8 @@
// mantissa expressed as an unsigned long long.
from_unsigned_long_long(fb.get_mantissa());
- // Scale the UINT64 representation to the fractional part of
- // the double and multiply with the base-2 exponent.
+ // Scale the unsigned long long representation to the fractional
+ // part of the float and multiply with the base-2 exponent.
const int p2 = fb.get_exponent() - (std::numeric_limits<float>::digits - 1);
if(p2 != 0) { operator*=(ef::pow2(static_cast<INT64>(p2))); }
@@ -201,8 +201,8 @@
// mantissa expressed as an unsigned long long.
from_unsigned_long_long(db.get_mantissa());
- // Scale the UINT64 representation to the fractional part of
- // the double and multiply with the base-2 exponent.
+ // Scale the unsigned long long representation to the fractional
+ // part of the double and multiply with the base-2 exponent.
const int p2 = db.get_exponent() - (std::numeric_limits<double>::digits - 1);
if(p2 != 0) { operator*=(ef::pow2(static_cast<INT64>(p2))); }
@@ -236,8 +236,8 @@
// mantissa expressed as an unsigned long long.
from_unsigned_long_long(ldb.get_mantissa());
- // Scale the UINT64 representation to the fractional part of
- // the double and multiply with the base-2 exponent.
+ // Scale the unsigned long long representation to the fractional
+ // part of the long double and multiply with the base-2 exponent.
const int p2 = ldb.get_exponent() - (std::numeric_limits<long double>::digits - 1);
if(p2 != 0) { operator*=(ef::pow2(static_cast<INT64>(p2))); }
@@ -282,66 +282,49 @@
fpclass (ef_finite),
prec_elem(ef_elem_number)
{
- // Create an e_float from mantissa and exponent. No, this ctor
- // does not maintain the full precision of double.
+ // Create an e_float from mantissa and exponent.
+ // This ctor does not maintain the full precision of double.
const bool mantissa_is_iszero = (::fabs(mantissa) < ((std::numeric_limits<double>::min)() * 2.0));
if(mantissa_is_iszero)
{
- if(exponent == static_cast<INT64>(0))
- {
- operator=(ef::one());
- }
- else
- {
- operator=(ef::zero());
- }
+ operator=((exponent == static_cast<INT64>(0)) ? ef::one() : ef::zero());
+ return;
}
- else
- {
- const bool b_neg = (mantissa < 0.0);
- neg = b_neg;
+ const bool b_neg = (mantissa < 0.0);
- double d = ((!b_neg) ? mantissa : -mantissa);
- INT64 e = exponent;
+ neg = b_neg;
- while(d > 1.0)
- {
- d /= 10.0;
- ++e;
- }
+ double d = ((!b_neg) ? mantissa : -mantissa);
+ INT64 e = exponent;
- while(d < 1.0)
- {
- d *= 10.0;
- --e;
- }
+ while(d > 10.0) { d /= 10.0; ++e; }
+ while(d < 1.0) { d *= 10.0; --e; }
- INT32 shift = static_cast<INT32>(e % static_cast<INT32>(ef_elem_digits10));
+ INT32 shift = static_cast<INT32>(e % static_cast<INT32>(ef_elem_digits10));
- while(static_cast<INT32>(shift-- % ef_elem_digits10) != static_cast<INT32>(0))
- {
- d *= 10.0;
- --e;
- }
+ while(static_cast<INT32>(shift-- % ef_elem_digits10) != static_cast<INT32>(0))
+ {
+ d *= 10.0;
+ --e;
+ }
- exp = e;
- neg = b_neg;
+ exp = e;
+ neg = b_neg;
- std::fill(data.begin(), data.end(), static_cast<UINT32>(0u));
+ std::fill(data.begin(), data.end(), static_cast<UINT32>(0u));
- static const INT32 digit_ratio = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) / static_cast<INT32>(ef_elem_digits10));
- static const INT32 digit_loops = static_cast<INT32>(digit_ratio + static_cast<INT32>(2));
+ static const INT32 digit_ratio = static_cast<INT32>(static_cast<INT32>(std::numeric_limits<double>::digits10) / static_cast<INT32>(ef_elem_digits10));
+ static const INT32 digit_loops = static_cast<INT32>(digit_ratio + static_cast<INT32>(2));
- for(INT32 i = static_cast<INT32>(0); i < digit_loops; i++)
- {
- UINT32 n = static_cast<UINT32>(static_cast<UINT64>(d));
- data[i] = static_cast<UINT32>(n);
- d -= static_cast<double>(n);
- d *= static_cast<double>(ef_elem_mask);
- }
+ for(INT32 i = static_cast<INT32>(0); i < digit_loops; i++)
+ {
+ UINT32 n = static_cast<UINT32>(static_cast<UINT64>(d));
+ data[i] = static_cast<UINT32>(n);
+ d -= static_cast<double>(n);
+ d *= static_cast<double>(ef_elem_mask);
}
}
@@ -418,7 +401,7 @@
carry = static_cast<UINT64>(sum / static_cast<UINT32>(ef_elem_mask));
}
- w[0] = static_cast<UINT32>(carry);
+ w[0u] = static_cast<UINT32>(carry);
}
UINT32 efx::e_float::mul_loop_n(UINT32* const u, UINT32 n, const INT32 p)
@@ -452,12 +435,17 @@
void efx::e_float::precision(const INT32 prec_digits)
{
- const INT32 elems_x = static_cast<INT32>( static_cast<INT32>( (prec_digits + (ef_elem_digits10 / 2)) / ef_elem_digits10)
+ if(prec_digits >= ef_digits10)
+ {
+ prec_elem = ef_elem_number;
+ }
+ else
+ {
+ const INT32 elems = static_cast<INT32>( static_cast<INT32>( (prec_digits + (ef_elem_digits10 / 2)) / ef_elem_digits10)
+ static_cast<INT32>(((prec_digits % ef_elem_digits10) != 0) ? 1 : 0));
- const INT32 elems = ((elems_x < static_cast<INT32>(2)) ? static_cast<INT32>(2) : elems_x);
-
- prec_elem = ((elems > ef_elem_number) ? ef_elem_number : elems);
+ prec_elem = (std::min)(ef_elem_number, (std::max)(elems, static_cast<INT32>(2)));
+ }
}
efx::e_float& efx::e_float::operator=(const e_float& v)
@@ -496,7 +484,7 @@
// Get the offset for the add/sub operation.
static const INT64 max_delta_exp = static_cast<INT64>((ef_elem_number - 1) * ef_elem_digits10);
- const INT64 ofs_exp = exp - v.exp;
+ const INT64 ofs_exp = static_cast<INT64>(exp - v.exp);
// Check if the operation is out of range, requiring special handling.
if(v.iszero() || (ofs_exp > max_delta_exp))
@@ -659,9 +647,7 @@
}
// Is it necessary to justify the data?
- const array_type::const_iterator first_nonzero_elem = std::find_if(data.begin(),
- data.end(),
- data_elem_is_nonzero_predicate);
+ const array_type::const_iterator first_nonzero_elem = std::find_if(data.begin(), data.end(), data_elem_is_non_zero_predicate);
if(first_nonzero_elem != data.begin())
{
@@ -677,13 +663,8 @@
// Justify the data
const std::size_t sj = static_cast<std::size_t>(std::distance<array_type::const_iterator>(data.begin(), first_nonzero_elem));
- std::copy(data.begin() + static_cast<std::size_t>(sj),
- data.end(),
- data.begin());
-
- std::fill(data.end() - sj,
- data.end(),
- static_cast<UINT32>(0u));
+ std::copy(data.begin() + static_cast<std::size_t>(sj), data.end(), data.begin());
+ std::fill(data.end() - sj, data.end(), static_cast<UINT32>(0u));
exp -= static_cast<INT64>(sj * static_cast<std::size_t>(ef_elem_digits10));
}
@@ -868,10 +849,7 @@
}
// Set up the multiplication loop.
- const INT32 jmax = static_cast<INT32>(data.rend() - std::find_if(data.rbegin(),
- data.rend(),
- data_elem_is_nonzero_predicate));
-
+ const INT32 jmax = static_cast<INT32>(data.rend() - std::find_if(data.rbegin(), data.rend(), data_elem_is_non_zero_predicate));
const INT32 jm1 = static_cast<INT32>(jmax + static_cast<INT32>(1));
const INT32 prec = static_cast<INT32>(prec_elem);
const INT32 jm = (std::min)(jm1, prec);
@@ -1099,9 +1077,9 @@
// Setup the iteration.
// Estimate the square root using simple manipulations.
const double sqd = ::sqrt(dd);
-
+
*this = e_float(sqd, static_cast<INT64>(ne / static_cast<INT64>(2)));
-
+
// Estimate 1.0 / (2.0 * x0) using simple manipulations.
e_float vi(0.5 / sqd, static_cast<INT64>(-ne / static_cast<INT64>(2)));
@@ -1129,7 +1107,7 @@
// Next iteration of *this
*this += vi * (-(*this * *this) + x);
}
-
+
prec_elem = ef_elem_number;
return *this;
@@ -1141,8 +1119,13 @@
// Return +1 for u > v
// 0 for u = v
// -1 for u < v
- return ((data != vd) ? ((data > vd) ? static_cast<INT32>(1) : static_cast<INT32>(-1))
- : static_cast<INT32>(0));
+
+ array_type::size_type i = static_cast<array_type::size_type>(0u);
+
+ while((i < data.size()) && (data[i] == vd[i])) { ++i; }
+
+ return ((i == data.size()) ? static_cast<INT32>(0)
+ : ((data[i] > vd[i]) ? static_cast<INT32>(1) : static_cast<INT32>(-1)));
}
INT32 efx::e_float::cmp(const e_float& v) const
@@ -1156,12 +1139,12 @@
{
// The value of *this is zero and v is either zero or non-zero.
return (v.iszero() ? static_cast<INT32>(0)
- : (v.isneg() ? static_cast<INT32>(1) : static_cast<INT32>(-1)));
+ : (v.neg ? static_cast<INT32>(1) : static_cast<INT32>(-1)));
}
else if(v.iszero())
{
// The value of v is zero and *this is non-zero.
- return (isneg() ? static_cast<INT32>(-1) : static_cast<INT32>(1));
+ return (neg ? static_cast<INT32>(-1) : static_cast<INT32>(1));
}
else
{
@@ -1185,75 +1168,66 @@
// Compare the data.
const INT32 val_cmp_data = cmp_data(v.data);
- return (neg ? static_cast<INT32>(-val_cmp_data) : val_cmp_data);
+ return ((!neg) ? val_cmp_data : static_cast<INT32>(-val_cmp_data));
}
}
}
bool efx::e_float::iszero(void) const
{
- return ( isfinite()
- && ( (data[0u] == static_cast<UINT32>(0u))
- || (exp < std::numeric_limits<e_float>::min_exponent10)));
+ if(fpclass == ef_finite)
+ {
+ return ((data[0u] == static_cast<UINT32>(0u)) || (exp < std::numeric_limits<e_float>::min_exponent10));
+ }
+ else
+ {
+ return false;
+ }
}
bool efx::e_float::isone(void) const
{
// Check if the value of *this is identically 1 or very close to 1.
- if( (!neg)
- && isfinite()
- && (data[0u] == static_cast<UINT64>(1u))
- && (exp == static_cast<INT64>(0))
- )
- {
- const array_type::const_iterator pos_first_nonzero_elem =
- std::find_if(data.begin(),
- data.end(),
- data_elem_is_nonzero_predicate);
+ const bool not_negative_and_is_finite = ((!neg) && isfinite());
- return (pos_first_nonzero_elem == data.end());
- }
- else
+ if(not_negative_and_is_finite)
{
- return false;
+ if((data[0u] == static_cast<UINT32>(1u)) && (exp == static_cast<INT64>(0)))
+ {
+ const array_type::const_iterator it_non_zero = std::find_if(data.begin(), data.end(), data_elem_is_non_zero_predicate);
+ return (it_non_zero == data.end());
+ }
+ else if((data[0u] == static_cast<UINT32>(ef_elem_mask - 1)) && (exp == static_cast<INT64>(-ef_elem_digits10)))
+ {
+ const array_type::const_iterator it_non_nine = std::find_if(data.begin(), data.end(), data_elem_is_non_nine_predicate);
+ return (it_non_nine == data.end());
+ }
}
+
+ return false;
}
bool efx::e_float::isint(void) const
{
- // Check if the value of *this is pure integer.
- if(!isfinite())
- {
- return false;
- }
+ if(fpclass != ef_finite) { return false; }
- if(iszero())
- {
- return true;
- }
+ if(iszero()) { return true; }
- if(exp < static_cast<INT64>(0))
- {
- // The absolute value of the number is smaller than 1.
- // It can not be an integer.
- return false;
- }
+ if(exp < static_cast<INT64>(0)) { return false; } // |*this| < 1.
const array_type::size_type offset_decimal_part = static_cast<array_type::size_type>(exp / ef_elem_digits10) + 1u;
- if(offset_decimal_part >= data.size())
+ if(offset_decimal_part >= static_cast<array_type::size_type>(ef_elem_number))
{
// The number is too large to resolve the integer part.
// It considered to be a pure integer.
return true;
}
- array_type::const_iterator pos_first_nonzero_elem = std::find_if(data.begin() + offset_decimal_part,
- data.end(),
- data_elem_is_nonzero_predicate);
+ array_type::const_iterator it_first_non_zero = std::find_if(data.begin() + offset_decimal_part, data.end(), data_elem_is_non_zero_predicate);
- return (pos_first_nonzero_elem == data.end());
+ return (it_first_non_zero == data.end());
}
efx::e_float& efx::e_float::operator++(void) { return *this += ef::one(); }
@@ -1284,16 +1258,13 @@
static const double d_mask = static_cast<double>(ef_elem_mask);
- mantissa = static_cast<double>(data[0])
- + (static_cast<double>(data[1]) / d_mask
- + static_cast<double>(data[2]) / (d_mask * d_mask));
+ mantissa = static_cast<double>(data[0])
+ + (static_cast<double>(data[1]) / d_mask)
+ + ((static_cast<double>(data[2]) / d_mask) / d_mask);
mantissa /= static_cast<double>(p10);
- if(neg)
- {
- mantissa = -mantissa;
- }
+ if(neg) { mantissa = -mantissa; }
}
double efx::e_float::extract_double(void) const
@@ -1494,9 +1465,7 @@
const size_t first_clear = (static_cast<size_t>(x.exp) / static_cast<size_t>(ef_elem_digits10)) + 1u;
const size_t last_clear = static_cast<size_t>(ef_elem_number);
- std::fill(x.data.begin() + first_clear,
- x.data.begin() + last_clear,
- static_cast<UINT32>(0u));
+ std::fill(x.data.begin() + first_clear, x.data.begin() + last_clear, static_cast<UINT32>(0u));
return x;
}
@@ -1540,14 +1509,12 @@
const size_t first_clear = static_cast<size_t>(ef_elem_number - first_copy);
const size_t last_clear = static_cast<size_t>(ef_elem_number);
- std::fill(x.data.begin() + first_clear,
- x.data.begin() + last_clear,
- static_cast<UINT32>(0u));
+ std::fill(x.data.begin() + first_clear, x.data.begin() + last_clear, static_cast<UINT32>(0u));
// Is it necessary to justify the data?
const array_type::const_iterator first_nonzero_elem = std::find_if(x.data.begin(),
x.data.end(),
- data_elem_is_nonzero_predicate);
+ data_elem_is_non_zero_predicate);
std::size_t sj = static_cast<std::size_t>(0u);
@@ -1583,57 +1550,33 @@
const efx::e_float& efx::e_float::my_value_nan(void) const
{
static e_float val = ef::zero();
-
val.fpclass = ef_NaN;
-
static const e_float qnan(val);
-
return qnan;
}
const efx::e_float& efx::e_float::my_value_inf(void) const
{
static e_float val = ef::zero();
-
val.fpclass = ef_inf;
-
static const e_float inf(val);
-
return inf;
}
const efx::e_float& efx::e_float::my_value_max(void) const
{
static const INT64 exp10_max = std::numeric_limits<e_float>::max_exponent10;
-
static const e_float val("1E" + Util::lexical_cast(exp10_max));
-
return val;
}
const efx::e_float& efx::e_float::my_value_min(void) const
{
static const INT64 exp10_min = std::numeric_limits<e_float>::min_exponent10;
-
static const e_float val("1E" + Util::lexical_cast(exp10_min));
-
return val;
}
-INT64 efx::e_float::get_order_exact(void) const
-{
- if(iszero())
- {
- return static_cast<INT64>(0);
- }
- else
- {
- const std::string str = Util::lexical_cast(data[0]);
- const INT64 str_len_minus_one = static_cast<INT64>(static_cast<INT64>(str.length()) - static_cast<INT64>(1));
- return static_cast<INT64>(exp + str_len_minus_one);
- }
-}
-
INT64 efx::e_float::get_order_fast(void) const
{
if(iszero())
@@ -1642,8 +1585,8 @@
}
else
{
- const double dx = ::log10(static_cast<double>(data[0])) + 0.5;
- return static_cast<INT64>(exp + static_cast<INT64>(dx));
+ const double dx = ::log10(static_cast<double>(data[0])) + (std::numeric_limits<double>::epsilon() * 0.9);
+ return static_cast<INT64>(exp + static_cast<INT64>(static_cast<INT32>(dx)));
}
}
@@ -1699,13 +1642,13 @@
}
else
{
- // Round up this digit
+ // Round up this digit.
++str.at(ix);
}
}
else
{
- // Round up last digit.
+ // Round up the last digit.
++str[ix];
}
}
@@ -1716,7 +1659,7 @@
{
std::string str(s);
- // Get possible exponent and remove it
+ // Get a possible exponent and remove it.
exp = static_cast<INT64>(0);
std::size_t pos;
@@ -1725,12 +1668,12 @@
|| ((pos = str.find('E')) != std::string::npos)
)
{
- exp = Util::numeric_cast<INT64>(std::string(str.c_str() + (pos + 1u)));
-
+ // Remove the exponent part from the string.
+ exp = Util::numeric_cast<INT64>(static_cast<const char* const>(str.c_str() + (pos + 1u)));
str = str.substr(static_cast<std::size_t>(0u), pos);
}
- // Get possible +/- sign and remove it
+ // Get a possible +/- sign and remove it.
neg = false;
if((pos = str.find(static_cast<char>('-'))) != std::string::npos)
@@ -1744,14 +1687,15 @@
str.erase(pos, static_cast<std::size_t>(1u));
}
- // Remove leading zeros for all input types
+ // Remove leading zeros for all input types.
const std::string::iterator fwd_it_leading_zero = std::find_if(str.begin(), str.end(), char_is_nonzero_predicate);
if(fwd_it_leading_zero != str.begin())
{
if(fwd_it_leading_zero == str.end())
{
- // The string contains nothing but leading zeros. The string is zero.
+ // The string contains nothing but leading zeros.
+ // This string represents zero.
operator=(ef::zero());
return true;
}
@@ -1761,53 +1705,50 @@
}
}
- // Put the input string into the standard e_float input
- // form aaa.bbbbE+/-n, where aa has 1...unit digits, bbbb
- // has an even multiple of unit digits which are possibly
- // zero padded on the right-end, and n is a signed 32-bit integer
- // which is an even multiple of unit.
+ // Put the input string into the standard e_float input form
+ // aaa.bbbbE+/-n, where aa has 1...ef_elem_digits10, bbbb has an
+ // even multiple of ef_elem_digits10 which are possibly zero padded
+ // on the right-end, and n is a signed 32-bit integer which is an
+ // even multiple of ef_elem_digits10.
- // Find possible decimal point
+ // Find a possible decimal point.
pos = str.find(static_cast<char>('.'));
if(pos != std::string::npos)
{
// Remove all trailing insignificant zeros.
- const std::string::const_reverse_iterator rev_it_non_zero_elem =
- std::find_if(str.rbegin(), str.rend(), char_is_nonzero_predicate);
+ const std::string::const_reverse_iterator rit_non_zero = std::find_if(str.rbegin(), str.rend(), char_is_nonzero_predicate);
- if(rev_it_non_zero_elem != str.rbegin())
+ if(rit_non_zero != str.rbegin())
{
- const std::string::size_type ofs = str.length() - std::distance<std::string::const_reverse_iterator>(str.rbegin(), rev_it_non_zero_elem);
+ const std::string::size_type ofs = str.length() - std::distance<std::string::const_reverse_iterator>(str.rbegin(), rit_non_zero);
str.erase(str.begin() + ofs, str.end());
}
- // Check if input is identically zero
+ // Check if the input is identically zero.
if(str == std::string("."))
{
operator=(ef::zero());
return true;
}
-
+
// Remove leading significant zeros just after the decimal point
// and adjust the exponent accordingly.
// Note that the while-loop operates only on strings of the form ".000abcd..."
// and peels away the zeros just after the decimal point.
if(str.at(static_cast<std::size_t>(0u)) == static_cast<char>('.'))
{
- const std::string::iterator fwd_it_first_non_zero_decimal =
- std::find_if(str.begin() + 1u, str.end(), char_is_nonzero_predicate);
+ const std::string::iterator it_non_zero = std::find_if(str.begin() + 1u, str.end(), char_is_nonzero_predicate);
- std::string::size_type delta_exp = static_cast<std::size_t>(0u);
+ std::size_t delta_exp = static_cast<std::size_t>(0u);
if(str.at(static_cast<std::size_t>(1u)) == static_cast<char>('0'))
{
- delta_exp = std::distance<std::string::const_iterator>(str.begin() + 1u,
- fwd_it_first_non_zero_decimal);
+ delta_exp = std::distance<std::string::const_iterator>(str.begin() + 1u, it_non_zero);
}
- // Bring one single digit into the mantissa and adjust exponent accordingly
- str.erase(str.begin(), fwd_it_first_non_zero_decimal);
+ // Bring one single digit into the mantissa and adjust exponent accordingly.
+ str.erase(str.begin(), it_non_zero);
str.insert(static_cast<std::size_t>(1u), ".");
exp -= static_cast<INT64>(delta_exp + 1u);
}
@@ -1818,8 +1759,7 @@
str.append(".");
}
- // Shift the decimal point such that the exponent is an even
- // multiple of std::numeric_limits<e_float>::elem_digits10
+ // Shift the decimal point such that the exponent is an even multiple of ef_elem_digits10.
std::size_t n_shift = static_cast<std::size_t>(0u);
const std::size_t n_exp_rem = static_cast<std::size_t>(exp % static_cast<INT64>(ef_elem_digits10));
@@ -1852,7 +1792,7 @@
exp -= static_cast<INT64>(n_shift);
}
- // Cut the size of the mantissa to <= ef_elem_digits10
+ // Cut the size of the mantissa to <= ef_elem_digits10.
pos = str.find(static_cast<char>('.'));
pos_plus_one = static_cast<std::size_t>(pos + 1u);
@@ -1869,8 +1809,8 @@
exp += static_cast<INT64>(static_cast<INT64>(n) * static_cast<INT64>(ef_elem_digits10));
}
- // Pad the decimal part such that its value
- // is an even multiple of ef_elem_digits10
+ // Pad the decimal part such that its value is an even
+ // multiple of ef_elem_digits10.
pos = str.find(static_cast<char>('.'));
pos_plus_one = static_cast<std::size_t>(pos + 1u);
@@ -1884,7 +1824,7 @@
str.append(static_cast<std::size_t>(n_cnt), static_cast<char>('0'));
}
- // Truncate decimal part if it is too long
+ // Truncate decimal part if it is too long.
const std::size_t max_dec = static_cast<std::size_t>((ef_elem_number - 1) * ef_elem_digits10);
if(static_cast<std::size_t>(str.length() - pos) > max_dec)
@@ -1893,9 +1833,10 @@
static_cast<std::size_t>(pos_plus_one + max_dec));
}
- // Now the input string has the standard e_float input form (see comment above).
+ // Now the input string has the standard e_float input form.
+ // (See the comment above.)
- // Set the data to 0.
+ // Set all the data elements to 0.
std::fill(data.begin(), data.end(), static_cast<UINT32>(0u));
// Extract the data.
Modified: sandbox/e_float/libs/e_float/src/e_float/gmp/e_float_gmp.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/e_float/gmp/e_float_gmp.cpp (original)
+++ sandbox/e_float/libs/e_float/src/e_float/gmp/e_float_gmp.cpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -20,6 +20,7 @@
#include "e_float_gmp_protos.h"
#include "../../utility/util_lexical_cast.h"
+#include "../../utility/util_numeric_cast.h"
#if defined(__GNUC__)
static inline INT32 _isnan (float x) { return static_cast<INT32>(std::isnan <float>(x)); }
@@ -44,13 +45,12 @@
void gmp::e_float::init(void)
{
- static bool precision_is_initialized = false;
+ static bool precision_is_initialized;
- if(!precision_is_initialized)
+ if(precision_is_initialized == false)
{
precision_is_initialized = true;
-
- ::mpf_set_default_prec(static_cast<unsigned long int>(ef_digits2));
+ ::mpf_set_default_prec(static_cast<unsigned long>(ef_digits2 + static_cast<INT32>(4)));
}
}
@@ -69,14 +69,14 @@
gmp::e_float::e_float() : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
::mpf_init(rop);
}
gmp::e_float::e_float(const char n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
const bool b_neg = (std::numeric_limits<char>::is_signed ? (n < static_cast<char>(0)) : false);
@@ -85,7 +85,7 @@
}
gmp::e_float::e_float(const wchar_t n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
const bool b_neg = (std::numeric_limits<wchar_t>::is_signed ? (n < static_cast<wchar_t>(0)) : false);
@@ -94,7 +94,7 @@
}
gmp::e_float::e_float(const signed char n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
const bool b_neg = (n < static_cast<signed char>(0));
@@ -103,7 +103,7 @@
}
gmp::e_float::e_float(const signed short n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
const bool b_neg = (n < static_cast<signed short>(0));
@@ -112,7 +112,7 @@
}
gmp::e_float::e_float(const signed int n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
const bool b_neg = (n < static_cast<signed int>(0));
@@ -121,7 +121,7 @@
}
gmp::e_float::e_float(const signed long n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
const bool b_neg = (n < static_cast<signed long>(0));
@@ -130,7 +130,7 @@
}
gmp::e_float::e_float(const signed long long n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
const bool b_neg = (n < static_cast<signed long long>(0));
@@ -139,42 +139,42 @@
}
gmp::e_float::e_float(const unsigned char n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
from_unsigned_long(static_cast<unsigned long>(n));
}
gmp::e_float::e_float(const unsigned short n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
from_unsigned_long(static_cast<unsigned long>(n));
}
gmp::e_float::e_float(const unsigned int n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
from_unsigned_long(static_cast<unsigned long>(n));
}
gmp::e_float::e_float(const unsigned long n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
from_unsigned_long(static_cast<unsigned long>(n));
}
gmp::e_float::e_float(const unsigned long long n) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
from_unsigned_long_long(static_cast<unsigned long long>(n));
}
gmp::e_float::e_float(const float f) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
@@ -192,8 +192,8 @@
// mantissa expressed as an unsigned long long.
from_unsigned_long_long(fb.get_mantissa());
- // Scale the UINT64 representation to the fractional part of
- // the double and multiply with the base-2 exponent.
+ // Scale the unsigned long long representation to the fractional
+ // part of the float and multiply with the base-2 exponent.
const int p2 = fb.get_exponent() - (std::numeric_limits<float>::digits - 1);
if(p2 != 0) { operator*=(ef::pow2(static_cast<INT64>(p2))); }
@@ -205,14 +205,14 @@
}
gmp::e_float::e_float(const double d) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
::mpf_init_set_d(rop, d);
}
gmp::e_float::e_float(const long double ld) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
@@ -230,8 +230,8 @@
// mantissa expressed as an unsigned long long.
from_unsigned_long_long(fb.get_mantissa());
- // Scale the UINT64 representation to the fractional part of
- // the double and multiply with the base-2 exponent.
+ // Scale the unsigned long long representation to the fractional
+ // part of the long double and multiply with the base-2 exponent.
const int p2 = fb.get_exponent() - (std::numeric_limits<long double>::digits - 1);
if(p2 != 0) { operator*=(ef::pow2(static_cast<INT64>(p2))); }
@@ -243,14 +243,14 @@
}
gmp::e_float::e_float(const char* const s) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
static_cast<void>(rd_string(s));
}
gmp::e_float::e_float(const std::string& str) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
static_cast<void>(rd_string(str.c_str()));
@@ -264,7 +264,7 @@
}
gmp::e_float::e_float(const double mantissa, const INT64 exponent) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
@@ -290,7 +290,7 @@
}
gmp::e_float::e_float(const ::mpf_t& op) : fpclass (ef_finite),
- prec_elem(ef_digits10_tol)
+ prec_elem(ef_max_digits10)
{
init();
::mpf_init_set(rop, op);
@@ -877,10 +877,7 @@
std::copy(str.begin(), std::find(str.begin(), str.end(), c0), str_sll.begin());
// Get the signed long long result.
- std::stringstream ss;
- ss << str_sll;
- signed long long n;
- ss >> n;
+ const signed long long n = Util::numeric_cast<signed long long>(str_sll);
return ((!b_neg) ? n : -n);
}
@@ -928,10 +925,7 @@
std::copy(str.begin(), std::find(str.begin(), str.end(), c0), str_ull.begin());
// Get the unsigned long long result.
- std::stringstream ss;
- ss << str_ull;
- unsigned long long n;
- ss >> n;
+ const unsigned long long n = Util::numeric_cast<unsigned long long>(str_ull);
return n;
}
@@ -975,20 +969,10 @@
const std::string str = std::string(buf.data());
// Extract the base-10 exponent.
- INT64 my_exp;
-
const std::size_t pos_letter_e = str.rfind(static_cast<char>('e'));
- if(pos_letter_e != std::string::npos)
- {
- std::stringstream ss;
- ss << static_cast<const char*>(str.c_str() + (pos_letter_e + 1u));
- ss >> my_exp;
- }
- else
- {
- my_exp = static_cast<INT64>(0);
- }
+ const INT64 my_exp = ((pos_letter_e != std::string::npos) ? Util::numeric_cast<INT64>(static_cast<const char* const>(str.c_str() + (pos_letter_e + 1u)))
+ : static_cast<INT64>(0));
return my_exp;
}
@@ -1023,7 +1007,7 @@
const std::string str_fmt = std::string("%.") + (Util::lexical_cast(the_number_of_digits_scientific) + "Fe");
// Get the string representation of the e_float in scientific notation (lowercase, noshowpos).
- std::tr1::array<char, static_cast<std::size_t>(e_float::ef_digits10_tol + 32)> buf = {{ static_cast<char>('0') }};
+ std::tr1::array<char, static_cast<std::size_t>(e_float::ef_max_digits10 + 32)> buf = {{ static_cast<char>('0') }};
static_cast<void>(gmp_sprintf(buf.data(), str_fmt.c_str(), rop));
Modified: sandbox/e_float/libs/e_float/src/e_float/mpfr/e_float_mpfr.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/e_float/mpfr/e_float_mpfr.cpp (original)
+++ sandbox/e_float/libs/e_float/src/e_float/mpfr/e_float_mpfr.cpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -19,6 +19,7 @@
#include "e_float_mpfr_protos.h"
#include "../../utility/util_lexical_cast.h"
+#include "../../utility/util_numeric_cast.h"
namespace
{
@@ -31,13 +32,12 @@
void mpfr::e_float::init(void)
{
- static bool precision_is_initialized = false;
+ static bool precision_is_initialized;
- if(!precision_is_initialized)
+ if(precision_is_initialized == false)
{
precision_is_initialized = true;
-
- ::mpfr_set_default_prec(static_cast<mp_prec_t>(ef_digits2));
+ ::mpfr_set_default_prec(static_cast<mp_prec_t>(ef_digits2 + static_cast<INT32>(4)));
}
}
@@ -82,8 +82,8 @@
// mantissa expressed as an unsigned long long.
from_unsigned_long_long(fb.get_mantissa());
- // Scale the UINT64 representation to the fractional part of
- // the double and multiply with the base-2 exponent.
+ // Scale the unsigned long long representation to the fractional
+ // part of the float and multiply with the base-2 exponent.
const int p2 = fb.get_exponent() - (std::numeric_limits<float>::digits - 1);
if(p2 != 0) { mpfr_mul_2si(rop, rop, static_cast<signed long>(p2), GMP_RNDN); }
@@ -430,10 +430,7 @@
std::copy(str.begin(), str.begin() + str_sll.size(), str_sll.begin());
- std::stringstream ss;
- ss << str_sll;
- signed long long n;
- ss >> n;
+ const signed long long n = Util::numeric_cast<signed long long>(str_sll);
return ((!b_neg) ? n : -n);
}
@@ -470,10 +467,7 @@
std::copy(str.begin(), str.begin() + str_ull.size(), str_ull.begin());
- std::stringstream ss;
- ss << str_ull;
- unsigned long long n;
- ss >> n;
+ const unsigned long long n = Util::numeric_cast<unsigned long long>(str_ull);
return n;
}
@@ -516,20 +510,10 @@
const std::string str = std::string(buf.data());
// Extract the base-10 exponent.
- INT64 my_exp;
-
const std::size_t pos_letter_e = str.rfind(static_cast<char>('e'));
- if(pos_letter_e != std::string::npos)
- {
- std::stringstream ss;
- ss << static_cast<const char*>(str.c_str() + (pos_letter_e + 1u));
- ss >> my_exp;
- }
- else
- {
- my_exp = static_cast<INT64>(0);
- }
+ const INT64 my_exp = ((pos_letter_e != std::string::npos) ? Util::numeric_cast<INT64>(static_cast<const char* const>(str.c_str() + (pos_letter_e + 1u)))
+ : static_cast<INT64>(0));
return my_exp;
}
@@ -564,7 +548,7 @@
const std::string str_fmt = std::string("%.") + (Util::lexical_cast(the_number_of_digits_scientific) + "RNe");
// Get the string representation of the e_float in scientific notation (lowercase, noshowpos).
- std::tr1::array<char, static_cast<std::size_t>(e_float::ef_digits10_tol + 32)> buf = {{ static_cast<char>(0) }};
+ std::tr1::array<char, static_cast<std::size_t>(e_float::ef_max_digits10 + 32)> buf = {{ static_cast<char>(0) }};
::mpfr_sprintf(buf.data(), str_fmt.c_str(), rop);
Modified: sandbox/e_float/libs/e_float/src/functions/elementary/elementary_math.cpp
==============================================================================
--- sandbox/e_float/libs/e_float/src/functions/elementary/elementary_math.cpp (original)
+++ sandbox/e_float/libs/e_float/src/functions/elementary/elementary_math.cpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -155,8 +155,7 @@
bool ef::small_arg(const e_float& x)
{
static const double lim_d = static_cast<double>(static_cast<INT32>(ef::tol())) / 10.0;
- static const INT64 lim_n = static_cast<INT64>(lim_d);
- static const INT64 lim = (lim_n < 6 ? 6 : lim_n);
+ static const INT64 lim = (std::max)(static_cast<INT64>(lim_d), static_cast<INT64>(6));
return (x.order() < static_cast<INT64>(-lim));
}
Modified: sandbox/e_float/libs/e_float/src/functions/integer/prime_factors.h
==============================================================================
--- sandbox/e_float/libs/e_float/src/functions/integer/prime_factors.h (original)
+++ sandbox/e_float/libs/e_float/src/functions/integer/prime_factors.h 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -1,3 +1,13 @@
+
+// Copyright Christopher Kormanyos 2002 - 2011.
+// Distributed under 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)
+
+// This work is based on an earlier work:
+// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
+// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+
#ifndef _PRIME_FACTORS_2011_06_23_H_
#define _PRIME_FACTORS_2011_06_23_H_
Modified: sandbox/e_float/libs/e_float/src/utility/util_numeric_cast.h
==============================================================================
--- sandbox/e_float/libs/e_float/src/utility/util_numeric_cast.h (original)
+++ sandbox/e_float/libs/e_float/src/utility/util_numeric_cast.h 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -24,6 +24,15 @@
ss >> t;
return t;
}
+
+ template<typename T> inline T numeric_cast(const char* const s)
+ {
+ std::stringstream ss;
+ ss << s;
+ T t;
+ ss >> t;
+ return t;
+ }
}
#endif // _UTIL_NUMERIC_CAST_2009_11_24_H_
Added: sandbox/e_float/libs/e_float/test/f2c.h
==============================================================================
--- (empty file)
+++ sandbox/e_float/libs/e_float/test/f2c.h 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -0,0 +1,227 @@
+/* f2c.h -- Standard Fortran to C header file */
+
+/** barf [ba:rf] 2. "He suggested using FORTRAN, and everybody barfed."
+
+ - From The Shogakukan DICTIONARY OF NEW ENGLISH (Second edition) */
+
+#ifndef F2C_INCLUDE
+#define F2C_INCLUDE
+
+typedef long int integer;
+typedef unsigned long int uinteger;
+typedef char *address;
+typedef short int shortint;
+typedef float real;
+typedef double doublereal;
+typedef struct { real r, i; } complex;
+typedef struct { doublereal r, i; } doublecomplex;
+typedef long int logical;
+typedef short int shortlogical;
+typedef char logical1;
+typedef char integer1;
+//#ifdef INTEGER_STAR_8 /* Adjust for integer*8. */
+typedef long long longint; /* system-dependent */
+typedef unsigned long long ulongint; /* system-dependent */
+#define qbit_clear(a,b) ((a) & ~((ulongint)1 << (b)))
+#define qbit_set(a,b) ((a) | ((ulongint)1 << (b)))
+//#endif
+
+#define TRUE_ (1)
+#define FALSE_ (0)
+
+/* Extern is for use with -E */
+#ifndef Extern
+#define Extern extern
+#endif
+
+/* I/O stuff */
+
+#ifdef f2c_i2
+/* for -i2 */
+typedef short flag;
+typedef short ftnlen;
+typedef short ftnint;
+#else
+typedef long int flag;
+typedef long int ftnlen;
+typedef long int ftnint;
+#endif
+
+/*external read, write*/
+typedef struct
+{ flag cierr;
+ ftnint ciunit;
+ flag ciend;
+ char *cifmt;
+ ftnint cirec;
+} cilist;
+
+/*internal read, write*/
+typedef struct
+{ flag icierr;
+ char *iciunit;
+ flag iciend;
+ char *icifmt;
+ ftnint icirlen;
+ ftnint icirnum;
+} icilist;
+
+/*open*/
+typedef struct
+{ flag oerr;
+ ftnint ounit;
+ char *ofnm;
+ ftnlen ofnmlen;
+ char *osta;
+ char *oacc;
+ char *ofm;
+ ftnint orl;
+ char *oblnk;
+} olist;
+
+/*close*/
+typedef struct
+{ flag cerr;
+ ftnint cunit;
+ char *csta;
+} cllist;
+
+/*rewind, backspace, endfile*/
+typedef struct
+{ flag aerr;
+ ftnint aunit;
+} alist;
+
+/* inquire */
+typedef struct
+{ flag inerr;
+ ftnint inunit;
+ char *infile;
+ ftnlen infilen;
+ ftnint *inex; /*parameters in standard's order*/
+ ftnint *inopen;
+ ftnint *innum;
+ ftnint *innamed;
+ char *inname;
+ ftnlen innamlen;
+ char *inacc;
+ ftnlen inacclen;
+ char *inseq;
+ ftnlen inseqlen;
+ char *indir;
+ ftnlen indirlen;
+ char *infmt;
+ ftnlen infmtlen;
+ char *inform;
+ ftnint informlen;
+ char *inunf;
+ ftnlen inunflen;
+ ftnint *inrecl;
+ ftnint *innrec;
+ char *inblank;
+ ftnlen inblanklen;
+} inlist;
+
+#define VOID void
+
+union Multitype { /* for multiple entry points */
+ integer1 g;
+ shortint h;
+ integer i;
+ /* longint j; */
+ real r;
+ doublereal d;
+ complex c;
+ doublecomplex z;
+ };
+
+typedef union Multitype Multitype;
+
+/*typedef long int Long;*/ /* No longer used; formerly in Namelist */
+
+struct Vardesc { /* for Namelist */
+ char *name;
+ char *addr;
+ ftnlen *dims;
+ int type;
+ };
+typedef struct Vardesc Vardesc;
+
+struct Namelist {
+ char *name;
+ Vardesc **vars;
+ int nvars;
+ };
+typedef struct Namelist Namelist;
+
+#define abs(x) ((x) >= 0 ? (x) : -(x))
+#define dabs(x) (doublereal)abs(x)
+#ifndef min
+#define min(a,b) ((a) <= (b) ? (a) : (b))
+#endif
+#ifndef max
+#define max(a,b) ((a) >= (b) ? (a) : (b))
+#endif
+#define dmin(a,b) (doublereal)min(a,b)
+#define dmax(a,b) (doublereal)max(a,b)
+#define bit_test(a,b) ((a) >> (b) & 1)
+#define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
+#define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
+
+/* procedure parameter types for -A and -C++ */
+
+#define F2C_proc_par_types 1
+#ifdef __cplusplus
+typedef int /* Unknown procedure type */ (*U_fp)(...);
+typedef shortint (*J_fp)(...);
+typedef integer (*I_fp)(...);
+typedef real (*R_fp)(...);
+typedef doublereal (*D_fp)(...), (*E_fp)(...);
+typedef /* Complex */ VOID (*C_fp)(...);
+typedef /* Double Complex */ VOID (*Z_fp)(...);
+typedef logical (*L_fp)(...);
+typedef shortlogical (*K_fp)(...);
+typedef /* Character */ VOID (*H_fp)(...);
+typedef /* Subroutine */ int (*S_fp)(...);
+#else
+typedef int /* Unknown procedure type */ (*U_fp)();
+typedef shortint (*J_fp)();
+typedef integer (*I_fp)();
+typedef real (*R_fp)();
+typedef doublereal (*D_fp)(), (*E_fp)();
+typedef /* Complex */ VOID (*C_fp)();
+typedef /* Double Complex */ VOID (*Z_fp)();
+typedef logical (*L_fp)();
+typedef shortlogical (*K_fp)();
+typedef /* Character */ VOID (*H_fp)();
+typedef /* Subroutine */ int (*S_fp)();
+#endif
+/* E_fp is for real functions when -R is not specified */
+typedef VOID C_f; /* complex function */
+typedef VOID H_f; /* character function */
+typedef VOID Z_f; /* double complex function */
+typedef doublereal E_f; /* real function with -R not specified */
+
+/* undef any lower-case symbols that your C compiler predefines, e.g.: */
+
+#ifndef Skip_f2c_Undefs
+#undef cray
+#undef gcos
+#undef mc68010
+#undef mc68020
+#undef mips
+#undef pdp11
+#undef sgi
+#undef sparc
+#undef sun
+#undef sun2
+#undef sun3
+#undef sun4
+#undef u370
+#undef u3b
+#undef u3b2
+#undef u3b5
+#undef unix
+#undef vax
+#endif
+#endif
Added: sandbox/e_float/libs/e_float/test/linpack-benchmark.cpp
==============================================================================
--- (empty file)
+++ sandbox/e_float/libs/e_float/test/linpack-benchmark.cpp 2011-08-12 15:50:52 EDT (Fri, 12 Aug 2011)
@@ -0,0 +1,1252 @@
+///////////////////////////////////////////////////////////////////////////////
+// Copyright 2011 John Maddock. Distributed under 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)
+
+/* 1000d.f -- translated by f2c (version 20050501).
+You must link the resulting object file with libf2c:
+on Microsoft Windows system, link with libf2c.lib;
+on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+or, if you install libf2c.a in a standard place, with -lf2c -lm
+-- in that order, at the end of the command line, as in
+cc *.o -lf2c -lm
+Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+http://www.netlib.org/f2c/libf2c.zip
+*/
+
+/* http://netlib.sandia.gov/f2c/index.html */
+
+#define TEST_EF_NATIVE
+
+#include <boost/lexical_cast.hpp>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#ifdef TEST_BIG_NUMBER
+#include <boost/math/big_number/gmp.hpp>
+typedef boost::math::mpf_real_100 real_type;
+#elif defined(TEST_GMPXX)
+#include <gmpxx.h>
+typedef mpf_class real_type;
+#elif defined(TEST_EF_NATIVE)
+#include <e_float/e_float.hpp>
+typedef e_float real_type;
+#define CAST_TO_RT(x) real_type(x)
+#elif defined(TEST_MATH_EF)
+#define E_FLOAT_TYPE_GMP
+#include <boost/math/bindings/e_float.hpp>
+typedef boost::math::ef::e_float real_type;
+//#define CAST_TO_RT(x) real_type(x)
+#else
+typedef double real_type;
+#endif
+
+#ifndef CAST_TO_RT
+#define CAST_TO_RT(x) x
+#endif
+
+#include "f2c.h"
+
+extern "C" integer s_wsfe(cilist *);
+extern "C" integer e_wsfe(void);
+extern "C" integer do_fio(integer *, char *, ftnlen);
+extern "C" integer s_wsle(cilist *);
+extern "C" integer do_lio(integer *, integer *, char *, ftnlen);
+extern "C" integer e_wsle(void);
+extern "C" int s_stop(char *, ftnlen);
+
+#undef abs
+#undef dabs
+#define dabs abs
+#undef min
+#undef max
+#undef dmin
+#undef dmax
+#define dmin min
+#define dmax max
+
+#include <time.h>
+
+
+#if defined(TEST_EF_NATIVE)
+#include <e_float/e_float_functions.hpp>
+using namespace ef;
+#endif
+
+using std::min;
+using std::max;
+
+/* Table of constant values */
+
+static integer c__0 = 0;
+static real_type c_b7 = CAST_TO_RT(1);
+static integer c__1 = 1;
+static integer c__9 = 9;
+
+inline double second_(void)
+{
+ return ((double)(clock())) / CLOCKS_PER_SEC;
+}
+
+int dgefa_(real_type *, integer *, integer *, integer *, integer *), dgesl_(real_type *, integer *, integer *, integer *, real_type *, integer *);
+int dmxpy_(integer *, real_type *, integer *, integer *, real_type *, real_type *);
+int matgen_(real_type *, integer *, integer *, real_type *, real_type *);
+real_type epslon_(real_type *);
+real_type ran_(integer *);
+int dscal_(integer *, real_type *, real_type *, integer *);
+int daxpy_(integer *, real_type *, real_type *, integer *, real_type *, integer *);
+integer idamax_(integer *, real_type *, integer *);
+real_type ddot_(integer *, real_type *, integer *, real_type *, integer *);
+int daxpy_(integer *, real_type *, real_type *, integer *, real_type *, integer *);
+int dmxpy_(integer *, real_type *, integer *, integer *, real_type *, real_type *);
+
+
+extern "C" int MAIN__()
+{
+#ifdef TEST_BIG_NUMBER
+ std::cout << "Testing big_number<mpf_real<100> >" << std::endl;
+#elif defined(TEST_GMPXX)
+ std::cout << "Testing mpfr_class at 100 decimal degits" << std::endl;
+ mpf_set_default_prec(((100 + 1) * 1000L) / 301L);
+#elif defined(TEST_EF_NATIVE)
+ std::cout << "Testing xxx::e_float" << std::endl;
+#elif defined(TEST_MATH_EF)
+ std::cout << "Testing boost::math::ef::e_float" << std::endl;
+#else
+ std::cout << "Testing double" << std::endl;
+#endif
+
+
+ /* Format strings */
+ static char fmt_1[] = "(\002 Please send the results of this run to:\002"
+ "//\002 Jack J. Dongarra\002/\002 Computer Science Department\002/"
+ "\002 University of Tennessee\002/\002 Knoxville, Tennessee 37996"
+ "-1300\002//\002 Fax: 615-974-8296\002//\002 Internet: dongarra_at_c"
+ "s.utk.edu\002/)";
+ static char fmt_40[] = "(\002 norm. resid resid mac"
+ "hep\002,\002 x(1) x(n)\002)";
+ static char fmt_50[] = "(1p5e16.8)";
+ static char fmt_60[] = "(//\002 times are reported for matrices of or"
+ "der \002,i5)";
+ static char fmt_70[] = "(6x,\002factor\002,5x,\002solve\002,6x,\002tota"
+ "l\002,5x,\002mflops\002,7x,\002unit\002,6x,\002ratio\002)";
+ static char fmt_80[] = "(\002 times for array with leading dimension o"
+ "f\002,i4)";
+ static char fmt_110[] = "(6(1pe11.3))";
+
+ /* System generated locals */
+ integer i__1;
+ real_type d__1, d__2, d__3;
+
+ /* Builtin functions */
+
+ /* Local variables */
+ static real_type a[1001000] /* was [1001][1000] */, b[1000];
+ static integer i__, n;
+ static real_type x[1000];
+ static double t1;
+ static integer lda;
+ static double ops;
+ static real_type eps;
+ static integer info;
+ static double time[6], cray, total;
+ static integer ipvt[1000];
+ static real_type resid, norma;
+ static real_type normx;
+ static real_type residn;
+
+ /* Fortran I/O blocks */
+ static cilist io___4 = { 0, 6, 0, fmt_1, 0 };
+ static cilist io___20 = { 0, 6, 0, fmt_40, 0 };
+ static cilist io___21 = { 0, 6, 0, fmt_50, 0 };
+ static cilist io___22 = { 0, 6, 0, fmt_60, 0 };
+ static cilist io___23 = { 0, 6, 0, fmt_70, 0 };
+ static cilist io___24 = { 0, 6, 0, fmt_80, 0 };
+ static cilist io___25 = { 0, 6, 0, fmt_110, 0 };
+ static cilist io___26 = { 0, 6, 0, 0, 0 };
+
+
+ lda = 1001;
+
+ /* this program was updated on 10/12/92 to correct a */
+ /* problem with the random number generator. The previous */
+ /* random number generator had a short period and produced */
+ /* singular matrices occasionally. */
+
+ n = 1000;
+ cray = .056f;
+ s_wsfe(&io___4);
+ e_wsfe();
+ /* Computing 3rd power */
+ d__1 = (real_type) n;
+ /* Computing 2nd power */
+ d__2 = (real_type) n;
+// ops = boost::lexical_cast<double>(real_type(d__1 * (d__1 * d__1) * 2. / 3. + d__2 * d__2 * 2.));
+ ops = static_cast<double>(real_type(d__1 * (d__1 * d__1) * 2. / 3. + d__2 * d__2 * 2.));
+
+ matgen_(a, &lda, &n, b, &norma);
+
+ /* ****************************************************************** */
+ /* ****************************************************************** */
+ /* you should replace the call to dgefa and dgesl */
+ /* by calls to your linear equation solver. */
+ /* ****************************************************************** */
+ /* ****************************************************************** */
+
+ t1 = second_();
+ dgefa_(a, &lda, &n, ipvt, &info);
+ time[0] = second_() - t1;
+ t1 = second_();
+ dgesl_(a, &lda, &n, ipvt, b, &c__0);
+ time[1] = second_() - t1;
+ total = time[0] + time[1];
+ /* ****************************************************************** */
+ /* ****************************************************************** */
+
+ /* compute a residual to verify results. */
+
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__ - 1] = b[i__ - 1];
+ /* L10: */
+ }
+ matgen_(a, &lda, &n, b, &norma);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ b[i__ - 1] = -b[i__ - 1];
+ /* L20: */
+ }
+ dmxpy_(&n, b, &n, &lda, x, a);
+ resid = CAST_TO_RT(0);
+ normx = CAST_TO_RT(0);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ /* Computing MAX */
+ d__2 = resid, d__3 = (d__1 = b[i__ - 1], abs(d__1));
+ resid = max(d__2,d__3);
+ /* Computing MAX */
+ d__2 = normx, d__3 = (d__1 = x[i__ - 1], abs(d__1));
+ normx = max(d__2,d__3);
+ /* L30: */
+ }
+ eps = epslon_(&c_b7);
+ residn = resid / (n * norma * normx * eps);
+ s_wsfe(&io___20);
+ e_wsfe();
+ s_wsfe(&io___21);
+ /*
+ do_fio(&c__1, (char *)&residn, (ftnlen)sizeof(real_type));
+ do_fio(&c__1, (char *)&resid, (ftnlen)sizeof(real_type));
+ do_fio(&c__1, (char *)&eps, (ftnlen)sizeof(real_type));
+ do_fio(&c__1, (char *)&x[0], (ftnlen)sizeof(real_type));
+ do_fio(&c__1, (char *)&x[n - 1], (ftnlen)sizeof(real_type));
+ */
+ std::cout << std::setw(12) << std::setprecision(5) << residn << " " << resid << " " << eps << " " << x[0] << " " << x[n-1] << std::endl;
+ e_wsfe();
+
+ s_wsfe(&io___22);
+ do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
+ e_wsfe();
+ s_wsfe(&io___23);
+ e_wsfe();
+
+ time[2] = total;
+ time[3] = ops / (total * 1e6);
+ time[4] = 2. / time[3];
+ time[5] = total / cray;
+ s_wsfe(&io___24);
+ do_fio(&c__1, (char *)&lda, (ftnlen)sizeof(integer));
+ e_wsfe();
+ s_wsfe(&io___25);
+ for (i__ = 1; i__ <= 6; ++i__) {
+ // do_fio(&c__1, (char *)&time[i__ - 1], (ftnlen)sizeof(real_type));
+ std::cout << std::setw(12) << std::setprecision(5) << time[i__ - 1];
+ }
+ e_wsfe();
+ s_wsle(&io___26);
+ do_lio(&c__9, &c__1, " end of tests -- this version dated 10/12/92", (
+ ftnlen)44);
+ e_wsle();
+
+ s_stop("", (ftnlen)0);
+ return 0;
+} /* MAIN__ */
+
+/* Subroutine */ int matgen_(real_type *a, integer *lda, integer *n,
+ real_type *b, real_type *norma)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2;
+ real_type d__1, d__2;
+
+ /* Local variables */
+ static integer i__, j;
+ static integer init[4];
+
+
+ /* Parameter adjustments */
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ --b;
+
+ /* Function Body */
+ init[0] = 1;
+ init[1] = 2;
+ init[2] = 3;
+ init[3] = 1325;
+ *norma = CAST_TO_RT(0);
+ i__1 = *n;
+ for (j = 1; j <= i__1; ++j) {
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ a[i__ + j * a_dim1] = ran_(init) - .5f;
+ /* Computing MAX */
+ d__2 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
+ *norma = max(d__2,*norma);
+ /* L20: */
+ }
+ /* L30: */
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ b[i__] = CAST_TO_RT(0);
+ /* L35: */
+ }
+ i__1 = *n;
+ for (j = 1; j <= i__1; ++j) {
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ b[i__] += a[i__ + j * a_dim1];
+ /* L40: */
+ }
+ /* L50: */
+ }
+ return 0;
+} /* matgen_ */
+
+/* Subroutine */ int dgefa_(real_type *a, integer *lda, integer *n, integer *
+ ipvt, integer *info)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2, i__3;
+
+ /* Local variables */
+ static integer j, k, l;
+ static real_type t;
+ static integer kp1, nm1;
+
+
+ /* dgefa factors a double precision matrix by gaussian elimination. */
+
+ /* dgefa is usually called by dgeco, but it can be called */
+ /* directly with a saving in time if rcond is not needed. */
+ /* (time for dgeco) = (1 + 9/n)*(time for dgefa) . */
+
+ /* on entry */
+
+ /* a double precision(lda, n) */
+ /* the matrix to be factored. */
+
+ /* lda integer */
+ /* the leading dimension of the array a . */
+
+ /* n integer */
+ /* the order of the matrix a . */
+
+ /* on return */
+
+ /* a an upper triangular matrix and the multipliers */
+ /* which were used to obtain it. */
+ /* the factorization can be written a = l*u where */
+ /* l is a product of permutation and unit lower */
+ /* triangular matrices and u is upper triangular. */
+
+ /* ipvt integer(n) */
+ /* an integer vector of pivot indices. */
+
+ /* info integer */
+ /* = 0 normal value. */
+ /* = k if u(k,k) .eq. 0.0 . this is not an error */
+ /* condition for this subroutine, but it does */
+ /* indicate that dgesl or dgedi will divide by zero */
+ /* if called. use rcond in dgeco for a reliable */
+ /* indication of singularity. */
+
+ /* linpack. this version dated 08/14/78 . */
+ /* cleve moler, university of new mexico, argonne national lab. */
+
+ /* subroutines and functions */
+
+ /* blas daxpy,dscal,idamax */
+
+ /* internal variables */
+
+
+
+ /* gaussian elimination with partial pivoting */
+
+ /* Parameter adjustments */
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ --ipvt;
+
+ /* Function Body */
+ *info = 0;
+ nm1 = *n - 1;
+ if (nm1 < 1) {
+ goto L70;
+ }
+ i__1 = nm1;
+ for (k = 1; k <= i__1; ++k) {
+ kp1 = k + 1;
+
+ /* find l = pivot index */
+
+ i__2 = *n - k + 1;
+ l = idamax_(&i__2, &a[k + k * a_dim1], &c__1) + k - 1;
+ ipvt[k] = l;
+
+ /* zero pivot implies this column already triangularized */
+
+ if (a[l + k * a_dim1] == 0.) {
+ goto L40;
+ }
+
+ /* interchange if necessary */
+
+ if (l == k) {
+ goto L10;
+ }
+ t = a[l + k * a_dim1];
+ a[l + k * a_dim1] = a[k + k * a_dim1];
+ a[k + k * a_dim1] = t;
+L10:
+
+ /* compute multipliers */
+
+ t = -1. / a[k + k * a_dim1];
+ i__2 = *n - k;
+ dscal_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1);
+
+ /* row elimination with column indexing */
+
+ i__2 = *n;
+ for (j = kp1; j <= i__2; ++j) {
+ t = a[l + j * a_dim1];
+ if (l == k) {
+ goto L20;
+ }
+ a[l + j * a_dim1] = a[k + j * a_dim1];
+ a[k + j * a_dim1] = t;
+L20:
+ i__3 = *n - k;
+ daxpy_(&i__3, &t, &a[k + 1 + k * a_dim1], &c__1, &a[k + 1 + j *
+ a_dim1], &c__1);
+ /* L30: */
+ }
+ goto L50;
+L40:
+ *info = k;
+L50:
+ /* L60: */
+ ;
+ }
+L70:
+ ipvt[*n] = *n;
+ if (a[*n + *n * a_dim1] == 0.) {
+ *info = *n;
+ }
+ return 0;
+} /* dgefa_ */
+
+/* Subroutine */ int dgesl_(real_type *a, integer *lda, integer *n, integer *
+ ipvt, real_type *b, integer *job)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2;
+
+ /* Local variables */
+ static integer k, l;
+ static real_type t;
+ static integer kb, nm1;
+
+
+ /* dgesl solves the double precision system */
+ /* a * x = b or trans(a) * x = b */
+ /* using the factors computed by dgeco or dgefa. */
+
+ /* on entry */
+
+ /* a double precision(lda, n) */
+ /* the output from dgeco or dgefa. */
+
+ /* lda integer */
+ /* the leading dimension of the array a . */
+
+ /* n integer */
+ /* the order of the matrix a . */
+
+ /* ipvt integer(n) */
+ /* the pivot vector from dgeco or dgefa. */
+
+ /* b double precision(n) */
+ /* the right hand side vector. */
+
+ /* job integer */
+ /* = 0 to solve a*x = b , */
+ /* = nonzero to solve trans(a)*x = b where */
+ /* trans(a) is the transpose. */
+
+ /* on return */
+
+ /* b the solution vector x . */
+
+ /* error condition */
+
+ /* a division by zero will occur if the input factor contains a */
+ /* zero on the diagonal. technically this indicates singularity */
+ /* but it is often caused by improper arguments or improper */
+ /* setting of lda . it will not occur if the subroutines are */
+ /* called correctly and if dgeco has set rcond .gt. 0.0 */
+ /* or dgefa has set info .eq. 0 . */
+
+ /* to compute inverse(a) * c where c is a matrix */
+ /* with p columns */
+ /* call dgeco(a,lda,n,ipvt,rcond,z) */
+ /* if (rcond is too small) go to ... */
+ /* do 10 j = 1, p */
+ /* call dgesl(a,lda,n,ipvt,c(1,j),0) */
+ /* 10 continue */
+
+ /* linpack. this version dated 08/14/78 . */
+ /* cleve moler, university of new mexico, argonne national lab. */
+
+ /* subroutines and functions */
+
+ /* blas daxpy,ddot */
+
+ /* internal variables */
+
+
+ /* Parameter adjustments */
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ --ipvt;
+ --b;
+
+ /* Function Body */
+ nm1 = *n - 1;
+ if (*job != 0) {
+ goto L50;
+ }
+
+ /* job = 0 , solve a * x = b */
+ /* first solve l*y = b */
+
+ if (nm1 < 1) {
+ goto L30;
+ }
+ i__1 = nm1;
+ for (k = 1; k <= i__1; ++k) {
+ l = ipvt[k];
+ t = b[l];
+ if (l == k) {
+ goto L10;
+ }
+ b[l] = b[k];
+ b[k] = t;
+L10:
+ i__2 = *n - k;
+ daxpy_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
+ /* L20: */
+ }
+L30:
+
+ /* now solve u*x = y */
+
+ i__1 = *n;
+ for (kb = 1; kb <= i__1; ++kb) {
+ k = *n + 1 - kb;
+ b[k] /= a[k + k * a_dim1];
+ t = -b[k];
+ i__2 = k - 1;
+ daxpy_(&i__2, &t, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
+ /* L40: */
+ }
+ goto L100;
+L50:
+
+ /* job = nonzero, solve trans(a) * x = b */
+ /* first solve trans(u)*y = b */
+
+ i__1 = *n;
+ for (k = 1; k <= i__1; ++k) {
+ i__2 = k - 1;
+ t = ddot_(&i__2, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
+ b[k] = (b[k] - t) / a[k + k * a_dim1];
+ /* L60: */
+ }
+
+ /* now solve trans(l)*x = y */
+
+ if (nm1 < 1) {
+ goto L90;
+ }
+ i__1 = nm1;
+ for (kb = 1; kb <= i__1; ++kb) {
+ k = *n - kb;
+ i__2 = *n - k;
+ b[k] += ddot_(&i__2, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
+ l = ipvt[k];
+ if (l == k) {
+ goto L70;
+ }
+ t = b[l];
+ b[l] = b[k];
+ b[k] = t;
+L70:
+ /* L80: */
+ ;
+ }
+L90:
+L100:
+ return 0;
+} /* dgesl_ */
+
+/* Subroutine */ int daxpy_(integer *n, real_type *da, real_type *dx,
+ integer *incx, real_type *dy, integer *incy)
+{
+ /* System generated locals */
+ integer i__1;
+
+ /* Local variables */
+ static integer i__, m, ix, iy, mp1;
+
+
+ /* constant times a vector plus a vector. */
+ /* uses unrolled loops for increments equal to one. */
+ /* jack dongarra, linpack, 3/11/78. */
+
+
+ /* Parameter adjustments */
+ --dy;
+ --dx;
+
+ /* Function Body */
+ if (*n <= 0) {
+ return 0;
+ }
+ if (*da == 0.) {
+ return 0;
+ }
+ if (*incx == 1 && *incy == 1) {
+ goto L20;
+ }
+
+ /* code for unequal increments or equal increments */
+ /* not equal to 1 */
+
+ ix = 1;
+ iy = 1;
+ if (*incx < 0) {
+ ix = (-(*n) + 1) * *incx + 1;
+ }
+ if (*incy < 0) {
+ iy = (-(*n) + 1) * *incy + 1;
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dy[iy] += *da * dx[ix];
+ ix += *incx;
+ iy += *incy;
+ /* L10: */
+ }
+ return 0;
+
+ /* code for both increments equal to 1 */
+
+
+ /* clean-up loop */
+
+L20:
+ m = *n % 4;
+ if (m == 0) {
+ goto L40;
+ }
+ i__1 = m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dy[i__] += *da * dx[i__];
+ /* L30: */
+ }
+ if (*n < 4) {
+ return 0;
+ }
+L40:
+ mp1 = m + 1;
+ i__1 = *n;
+ for (i__ = mp1; i__ <= i__1; i__ += 4) {
+ dy[i__] += *da * dx[i__];
+ dy[i__ + 1] += *da * dx[i__ + 1];
+ dy[i__ + 2] += *da * dx[i__ + 2];
+ dy[i__ + 3] += *da * dx[i__ + 3];
+ /* L50: */
+ }
+ return 0;
+} /* daxpy_ */
+
+real_type ddot_(integer *n, real_type *dx, integer *incx, real_type *dy,
+ integer *incy)
+{
+ /* System generated locals */
+ integer i__1;
+ real_type ret_val;
+
+ /* Local variables */
+ static integer i__, m, ix, iy, mp1;
+ static real_type dtemp;
+
+
+ /* forms the dot product of two vectors. */
+ /* uses unrolled loops for increments equal to one. */
+ /* jack dongarra, linpack, 3/11/78. */
+
+
+ /* Parameter adjustments */
+ --dy;
+ --dx;
+
+ /* Function Body */
+ ret_val = CAST_TO_RT(0);
+ dtemp = CAST_TO_RT(0);
+ if (*n <= 0) {
+ return ret_val;
+ }
+ if (*incx == 1 && *incy == 1) {
+ goto L20;
+ }
+
+ /* code for unequal increments or equal increments */
+ /* not equal to 1 */
+
+ ix = 1;
+ iy = 1;
+ if (*incx < 0) {
+ ix = (-(*n) + 1) * *incx + 1;
+ }
+ if (*incy < 0) {
+ iy = (-(*n) + 1) * *incy + 1;
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dtemp += dx[ix] * dy[iy];
+ ix += *incx;
+ iy += *incy;
+ /* L10: */
+ }
+ ret_val = dtemp;
+ return ret_val;
+
+ /* code for both increments equal to 1 */
+
+
+ /* clean-up loop */
+
+L20:
+ m = *n % 5;
+ if (m == 0) {
+ goto L40;
+ }
+ i__1 = m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dtemp += dx[i__] * dy[i__];
+ /* L30: */
+ }
+ if (*n < 5) {
+ goto L60;
+ }
+L40:
+ mp1 = m + 1;
+ i__1 = *n;
+ for (i__ = mp1; i__ <= i__1; i__ += 5) {
+ dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
+ i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
+ 4] * dy[i__ + 4];
+ /* L50: */
+ }
+L60:
+ ret_val = dtemp;
+ return ret_val;
+} /* ddot_ */
+
+/* Subroutine */ int dscal_(integer *n, real_type *da, real_type *dx,
+ integer *incx)
+{
+ /* System generated locals */
+ integer i__1, i__2;
+
+ /* Local variables */
+ static integer i__, m, mp1, nincx;
+
+
+ /* scales a vector by a constant. */
+ /* uses unrolled loops for increment equal to one. */
+ /* jack dongarra, linpack, 3/11/78. */
+
+
+ /* Parameter adjustments */
+ --dx;
+
+ /* Function Body */
+ if (*n <= 0) {
+ return 0;
+ }
+ if (*incx == 1) {
+ goto L20;
+ }
+
+ /* code for increment not equal to 1 */
+
+ nincx = *n * *incx;
+ i__1 = nincx;
+ i__2 = *incx;
+ for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
+ dx[i__] = *da * dx[i__];
+ /* L10: */
+ }
+ return 0;
+
+ /* code for increment equal to 1 */
+
+
+ /* clean-up loop */
+
+L20:
+ m = *n % 5;
+ if (m == 0) {
+ goto L40;
+ }
+ i__2 = m;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ dx[i__] = *da * dx[i__];
+ /* L30: */
+ }
+ if (*n < 5) {
+ return 0;
+ }
+L40:
+ mp1 = m + 1;
+ i__2 = *n;
+ for (i__ = mp1; i__ <= i__2; i__ += 5) {
+ dx[i__] = *da * dx[i__];
+ dx[i__ + 1] = *da * dx[i__ + 1];
+ dx[i__ + 2] = *da * dx[i__ + 2];
+ dx[i__ + 3] = *da * dx[i__ + 3];
+ dx[i__ + 4] = *da * dx[i__ + 4];
+ /* L50: */
+ }
+ return 0;
+} /* dscal_ */
+
+integer idamax_(integer *n, real_type *dx, integer *incx)
+{
+ /* System generated locals */
+ integer ret_val, i__1;
+ real_type d__1;
+
+ /* Local variables */
+ static integer i__, ix;
+ static real_type dmax__;
+
+
+ /* finds the index of element having max. dabsolute value. */
+ /* jack dongarra, linpack, 3/11/78. */
+
+
+ /* Parameter adjustments */
+ --dx;
+
+ /* Function Body */
+ ret_val = 0;
+ if (*n < 1) {
+ return ret_val;
+ }
+ ret_val = 1;
+ if (*n == 1) {
+ return ret_val;
+ }
+ if (*incx == 1) {
+ goto L20;
+ }
+
+ /* code for increment not equal to 1 */
+
+ ix = 1;
+ dmax__ = abs(dx[1]);
+ ix += *incx;
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ if ((d__1 = dx[ix], abs(d__1)) <= dmax__) {
+ goto L5;
+ }
+ ret_val = i__;
+ dmax__ = (d__1 = dx[ix], abs(d__1));
+L5:
+ ix += *incx;
+ /* L10: */
+ }
+ return ret_val;
+
+ /* code for increment equal to 1 */
+
+L20:
+ dmax__ = abs(dx[1]);
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ if ((d__1 = dx[i__], abs(d__1)) <= dmax__) {
+ goto L30;
+ }
+ ret_val = i__;
+ dmax__ = (d__1 = dx[i__], abs(d__1));
+L30:
+ ;
+ }
+ return ret_val;
+} /* idamax_ */
+
+real_type epslon_(real_type *x)
+{
+#ifdef TEST_BIG_NUMBER
+ return std::ldexp(1.0, 1 - ((100 + 1) * 1000L) / 301L);
+#elif defined(TEST_GMPXX)
+ return std::ldexp(1.0, 1 - ((100 + 1) * 1000L) / 301L);
+#else
+ return CAST_TO_RT(std::numeric_limits<real_type>::epsilon());
+#endif
+} /* epslon_ */
+
+/* Subroutine */ int mm_(real_type *a, integer *lda, integer *n1, integer *
+ n3, real_type *b, integer *ldb, integer *n2, real_type *c__,
+ integer *ldc)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2;
+
+ /* Local variables */
+ static integer i__, j;
+
+
+ /* purpose: */
+ /* multiply matrix b times matrix c and store the result in matrix a. */
+
+ /* parameters: */
+
+ /* a double precision(lda,n3), matrix of n1 rows and n3 columns */
+
+ /* lda integer, leading dimension of array a */
+
+ /* n1 integer, number of rows in matrices a and b */
+
+ /* n3 integer, number of columns in matrices a and c */
+
+ /* b double precision(ldb,n2), matrix of n1 rows and n2 columns */
+
+ /* ldb integer, leading dimension of array b */
+
+ /* n2 integer, number of columns in matrix b, and number of rows in */
+ /* matrix c */
+
+ /* c double precision(ldc,n3), matrix of n2 rows and n3 columns */
+
+ /* ldc integer, leading dimension of array c */
+
+ /* ---------------------------------------------------------------------- */
+
+ /* Parameter adjustments */
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ b_dim1 = *ldb;
+ b_offset = 1 + b_dim1;
+ b -= b_offset;
+ c_dim1 = *ldc;
+ c_offset = 1 + c_dim1;
+ c__ -= c_offset;
+
+ /* Function Body */
+ i__1 = *n3;
+ for (j = 1; j <= i__1; ++j) {
+ i__2 = *n1;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ a[i__ + j * a_dim1] = CAST_TO_RT(0);
+ /* L10: */
+ }
+ dmxpy_(n2, &a[j * a_dim1 + 1], n1, ldb, &c__[j * c_dim1 + 1], &b[
+ b_offset]);
+ /* L20: */
+ }
+
+ return 0;
+} /* mm_ */
+
+/* Subroutine */ int dmxpy_(integer *n1, real_type *y, integer *n2, integer *
+ ldm, real_type *x, real_type *m)
+{
+ /* System generated locals */
+ integer m_dim1, m_offset, i__1, i__2;
+
+ /* Local variables */
+ static integer i__, j, jmin;
+
+
+ /* purpose: */
+ /* multiply matrix m times vector x and add the result to vector y. */
+
+ /* parameters: */
+
+ /* n1 integer, number of elements in vector y, and number of rows in */
+ /* matrix m */
+
+ /* y double precision(n1), vector of length n1 to which is added */
+ /* the product m*x */
+
+ /* n2 integer, number of elements in vector x, and number of columns */
+ /* in matrix m */
+
+ /* ldm integer, leading dimension of array m */
+
+ /* x double precision(n2), vector of length n2 */
+
+ /* m double precision(ldm,n2), matrix of n1 rows and n2 columns */
+
+ /* ---------------------------------------------------------------------- */
+
+ /* cleanup odd vector */
+
+ /* Parameter adjustments */
+ --y;
+ m_dim1 = *ldm;
+ m_offset = 1 + m_dim1;
+ m -= m_offset;
+ --x;
+
+ /* Function Body */
+ j = *n2 % 2;
+ if (j >= 1) {
+ i__1 = *n1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] += x[j] * m[i__ + j * m_dim1];
+ /* L10: */
+ }
+ }
+
+ /* cleanup odd group of two vectors */
+
+ j = *n2 % 4;
+ if (j >= 2) {
+ i__1 = *n1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] = y[i__] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[
+ i__ + j * m_dim1];
+ /* L20: */
+ }
+ }
+
+ /* cleanup odd group of four vectors */
+
+ j = *n2 % 8;
+ if (j >= 4) {
+ i__1 = *n1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] = y[i__] + x[j - 3] * m[i__ + (j - 3) * m_dim1] + x[j - 2]
+ * m[i__ + (j - 2) * m_dim1] + x[j - 1] * m[i__ + (j - 1) *
+ m_dim1] + x[j] * m[i__ + j * m_dim1];
+ /* L30: */
+ }
+ }
+
+ /* cleanup odd group of eight vectors */
+
+ j = *n2 % 16;
+ if (j >= 8) {
+ i__1 = *n1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ y[i__] = y[i__] + x[j - 7] * m[i__ + (j - 7) * m_dim1] + x[j - 6]
+ * m[i__ + (j - 6) * m_dim1] + x[j - 5] * m[i__ + (j - 5) *
+ m_dim1] + x[j - 4] * m[i__ + (j - 4) * m_dim1] + x[j - 3]
+ * m[i__ + (j - 3) * m_dim1] + x[j - 2] * m[i__ + (j - 2)
+ * m_dim1] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] *
+ m[i__ + j * m_dim1];
+ /* L40: */
+ }
+ }
+
+ /* main loop - groups of sixteen vectors */
+
+ jmin = j + 16;
+ i__1 = *n2;
+ for (j = jmin; j <= i__1; j += 16) {
+ i__2 = *n1;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ y[i__] = y[i__] + x[j - 15] * m[i__ + (j - 15) * m_dim1] + x[j -
+ 14] * m[i__ + (j - 14) * m_dim1] + x[j - 13] * m[i__ + (j
+ - 13) * m_dim1] + x[j - 12] * m[i__ + (j - 12) * m_dim1]
+ + x[j - 11] * m[i__ + (j - 11) * m_dim1] + x[j - 10] * m[
+ i__ + (j - 10) * m_dim1] + x[j - 9] * m[i__ + (j - 9) *
+ m_dim1] + x[j - 8] * m[i__ + (j - 8) * m_dim1] + x[j - 7]
+ * m[i__ + (j - 7) * m_dim1] + x[j - 6] * m[i__ + (j - 6) *
+ m_dim1] + x[j - 5] * m[i__ + (j - 5) * m_dim1] + x[j - 4]
+ * m[i__ + (j - 4) * m_dim1] + x[j - 3] * m[i__ + (j - 3)
+ * m_dim1] + x[j - 2] * m[i__ + (j - 2) * m_dim1] + x[j -
+ 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[i__ + j *
+ m_dim1];
+ /* L50: */
+ }
+ /* L60: */
+ }
+ return 0;
+} /* dmxpy_ */
+
+real_type ran_(integer *iseed)
+{
+ /* System generated locals */
+ real_type ret_val;
+
+ /* Local variables */
+ static integer it1, it2, it3, it4;
+
+
+ /* modified from the LAPACK auxiliary routine 10/12/92 JD */
+ /* -- LAPACK auxiliary routine (version 1.0) -- */
+ /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
+ /* Courant Institute, Argonne National Lab, and Rice University */
+ /* February 29, 1992 */
+
+ /* .. Array Arguments .. */
+ /* .. */
+
+ /* Purpose */
+ /* ======= */
+
+ /* DLARAN returns a random double number from a uniform (0,1) */
+ /* distribution. */
+
+ /* Arguments */
+ /* ========= */
+
+ /* ISEED (input/output) INTEGER array, dimension (4) */
+ /* On entry, the seed of the random number generator; the array */
+ /* elements must be between 0 and 4095, and ISEED(4) must be */
+ /* odd. */
+ /* On exit, the seed is updated. */
+
+ /* Further Details */
+ /* =============== */
+
+ /* This routine uses a multiplicative congruential method with modulus */
+ /* 2**48 and multiplier 33952834046453 (see G.S.Fishman, */
+ /* 'Multiplicative congruential random number generators with modulus */
+ /* 2**b: an exhaustive analysis for b = 32 and a partial analysis for */
+ /* b = 48', Math. Comp. 189, pp 331-344, 1990). */
+
+ /* 48-bit integers are stored in 4 integer array elements with 12 bits */
+ /* per element. Hence the routine is portable across machines with */
+ /* integers of 32 bits or more. */
+
+ /* .. Parameters .. */
+ /* .. */
+ /* .. Local Scalars .. */
+ /* .. */
+ /* .. Intrinsic Functions .. */
+ /* .. */
+ /* .. Executable Statements .. */
+
+ /* multiply the seed by the multiplier modulo 2**48 */
+
+ /* Parameter adjustments */
+ --iseed;
+
+ /* Function Body */
+ it4 = iseed[4] * 2549;
+ it3 = it4 / 4096;
+ it4 -= it3 << 12;
+ it3 = it3 + iseed[3] * 2549 + iseed[4] * 2508;
+ it2 = it3 / 4096;
+ it3 -= it2 << 12;
+ it2 = it2 + iseed[2] * 2549 + iseed[3] * 2508 + iseed[4] * 322;
+ it1 = it2 / 4096;
+ it2 -= it1 << 12;
+ it1 = it1 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322 + iseed[4]
+ * 494;
+ it1 %= 4096;
+
+ /* return updated seed */
+
+ iseed[1] = it1;
+ iseed[2] = it2;
+ iseed[3] = it3;
+ iseed[4] = it4;
+
+ /* convert 48-bit integer to a double number in the interval (0,1) */
+
+ ret_val = ((real_type) it1 + ((real_type) it2 + ((real_type) it3 + (
+ real_type) it4 * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4)
+ * 2.44140625e-4;
+ return ret_val;
+
+ /* End of RAN */
+
+} /* ran_ */
+
+/*
+
+double x64 release results:
+~~~~~~~~~~~~~~~~~~
+
+norm. resid resid machep x(1) x(n)
+6.4915 7.207e-013 2.2204e-016 1 1
+
+
+times are reported for matrices of order 1000
+factor solve total mflops unit ratio
+times for array with leading dimension of1001
+0.58 0 0.58 1152.9 0.0017348 10.357
+
+efx::e_float x64 release (from sandbox) results:
+~~~~~~~~~~~~~~~~~~
+
+norm. resid resid machep x(1) x(n)
+4.83973e-023 2.41986e-119 1e-099 1 1
+
+times are reported for matrices of order 1000
+factor solve total mflops unit ratio
+times for array with leading dimension of1001
+159.38 0.48 159.86 4.1827 0.47816 2854.7
+
+mpfr::e_float x64 release (from sandbox) results:
+~~~~~~~~~~~~~~~~~~
+
+norm. resid resid machep x(1) x(n)
+1.30249e-016 6.51244e-113 1e-099 1 1
+
+
+times are reported for matrices of order 1000
+factor solve total mflops unit ratio
+times for array with leading dimension of1001
+163.53 0.51 164.04 4.0763 0.49064 2929.2
+
+gmp::e_float x64 release (from sandbox) results:
+~~~~~~~~~~~~~~~~~~
+
+norm. resid resid machep x(1) x(n)
+2.36707e-016 1.18354e-112 1e-099 1 1
+
+times are reported for matrices of order 1000
+factor solve total mflops unit ratio
+times for array with leading dimension of1001
+229.29 0.69 229.98 2.9075 0.68787 4106.7
+
+*/