Added new (better) forms for evaluating polynomials.

Added more Remez code, and appoximations for erf/erfc inverses. Updated and refactored erf code to use new approximations. Added more test cases. [SVN r3210]
2026-01-19 04:22:09 +00:00 · 2006-09-22 09:28:40 +00:00
parent ba3332ebd1
commit 5b9e2ec5a8
10 changed files with 922 additions and 101 deletions
--- a/include/boost/math/special_functions/detail/erf_inv.hpp
+++ b/include/boost/math/special_functions/detail/erf_inv.hpp
@@ -0,0 +1,432 @@
+//  (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)
+
+#ifndef BOOST_MATH_SF_ERF_INV_HPP
+#define BOOST_MATH_SF_ERF_INV_HPP
+
+namespace boost{ namespace math{ 
+
+namespace detail{
+//
+// The inverse erf and erfc functions share a common implementation,
+// this version is for 80-bit long double's and smaller:
+//
+template <class T>
+T erf_inv_imp(const T& p, const T& q, const boost::mpl::int_<64>*)
+{
+   using namespace std; // for ADL of std names.
+
+   T result = 0;
+   
+   if(p <= 0.5)
+   {
+      //
+      // Evaluate inverse erf using the rational approximation:
+      //
+      // x = p(p+10)(Y+R(p))
+      //
+      // Where Y is a constant, and R(p) is optimised for a low
+      // absolute error compared to |Y|.
+      //
+      // double: Max error found: 2.001849e-18
+      // long double: Max error found: 1.017064e-20
+      // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
+      //
+      static const float Y = 0.0891314744949340820313f;
+      static const T P[] = {    
+         -0.000508781949658280665617L,
+         -0.00836874819741736770379L,
+         0.0334806625409744615033L,
+         -0.0126926147662974029034L,
+         -0.0365637971411762664006L,
+         0.0219878681111168899165L,
+         0.00822687874676915743155L,
+         -0.00538772965071242932965L
+      };
+      static const T Q[] = {    
+         1,
+         -0.970005043303290640362L,
+         -1.56574558234175846809L,
+         1.56221558398423026363L,
+         0.662328840472002992063L,
+         -0.71228902341542847553L,
+         -0.0527396382340099713954L,
+         0.0795283687341571680018L,
+         -0.00233393759374190016776L,
+         0.000886216390456424707504L
+      };
+      T g = p * (p + 10);
+      T r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);
+      result = g * Y + g * r;
+   }
+   else if(q >= 0.25)
+   {
+      //
+      // Rational approximation for 0.5 > q >= 0.25
+      //
+      // x = sqrt(-2*log(q)) / (Y + R(q))
+      //
+      // Where Y is a constant, and R(q) is optimised for a low
+      // absolute error compared to Y.
+      //
+      // double : Max error found: 7.403372e-17
+      // long double : Max error found: 6.084616e-20
+      // Maximum Deviation Found (error term) 4.811e-20
+      //
+      static const float Y = 2.249481201171875f;
+      static const T P[] = {    
+         -0.202433508355938759655L,
+         0.105264680699391713268L,
+         8.37050328343119927838L,
+         17.6447298408374015486L,
+         -18.8510648058714251895L,
+         -44.6382324441786960818L,
+         17.445385985570866523L,
+         21.1294655448340526258L,
+         -3.67192254707729348546L
+      };
+      static const T Q[] = {    
+         1L,
+         6.24264124854247537712L,
+         3.9713437953343869095L,
+         -28.6608180499800029974L,
+         -20.1432634680485188801L,
+         48.5609213108739935468L,
+         10.8268667355460159008L,
+         -22.6436933413139721736L,
+         1.72114765761200282724L
+      };
+      T g = sqrt(-2 * log(q));
+      T xs = q - 0.25;
+      T r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+      result = g / (Y + r);
+   }
+   else
+   {
+      //
+      // For q < 0.25 we have a series of rational approximations all
+      // of the general form:
+      //
+      // let: x = sqrt(-log(q))
+      //
+      // Then the result is given by:
+      //
+      // x(Y+R(x-B))
+      //
+      // where Y is a constant, B is the lowest value of x for which 
+      // the approximation is valid, and R(x-B) is optimised for a low
+      // absolute error compared to Y.
+      //
+      // Note that almost all code will really go through the first
+      // or maybe second approximation.  After than we're dealing with very
+      // small input values indeed: 80 and 128 bit long double's go all the
+      // way down to ~ 1e-5000 so the "tail" is rather long...
+      //
+      T x = sqrt(-log(q));
+      if(x < 3)
+      {
+         // Max error found: 1.089051e-20
+         static const float Y = 0.807220458984375f;
+         static const T P[] = {    
+            -0.131102781679951906451L,
+            -0.163794047193317060787L,
+            0.117030156341995252019L,
+            0.387079738972604337464L,
+            0.337785538912035898924L,
+            0.142869534408157156766L,
+            0.0290157910005329060432L,
+            0.00214558995388805277169L,
+            -0.679465575181126350155e-6L,
+            0.285225331782217055858e-7L,
+            -0.681149956853776992068e-9L
+         };
+         static const T Q[] = {    
+            1,
+            3.46625407242567245975L,
+            5.38168345707006855425L,
+            4.77846592945843778382L,
+            2.59301921623620271374L,
+            0.848854343457902036425L,
+            0.152264338295331783612L,
+            0.01105924229346489121L
+         };
+         T xs = x - 1.125;
+         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+         result = Y * x + R * x;
+      }
+      else if(x < 6)
+      {
+         // Max error found: 8.389174e-21
+         static const float Y = 0.93995571136474609375f;
+         static const T P[] = {    
+            -0.0350353787183177984712L,
+            -0.00222426529213447927281L,
+            0.0185573306514231072324L,
+            0.00950804701325919603619L,
+            0.00187123492819559223345L,
+            0.000157544617424960554631L,
+            0.460469890584317994083e-5L,
+            -0.230404776911882601748e-9L,
+            0.266339227425782031962e-11L
+         };
+         static const T Q[] = {    
+            1L,
+            1.3653349817554063097L,
+            0.762059164553623404043L,
+            0.220091105764131249824L,
+            0.0341589143670947727934L,
+            0.00263861676657015992959L,
+            0.764675292302794483503e-4L
+         };
+         T xs = x - 3;
+         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+         result = Y * x + R * x;
+      }
+      else if(x < 18)
+      {
+         // Max error found: 1.481312e-19
+         static const float Y = 0.98362827301025390625f;
+         static const T P[] = {    
+            -0.0167431005076633737133L,
+            -0.00112951438745580278863L,
+            0.00105628862152492910091L,
+            0.000209386317487588078668L,
+            0.149624783758342370182e-4L,
+            0.449696789927706453732e-6L,
+            0.462596163522878599135e-8L,
+            -0.281128735628831791805e-13L,
+            0.99055709973310326855e-16L
+         };
+         static const T Q[] = {    
+            1L,
+            0.591429344886417493481L,
+            0.138151865749083321638L,
+            0.0160746087093676504695L,
+            0.000964011807005165528527L,
+            0.275335474764726041141e-4L,
+            0.282243172016108031869e-6L
+         };
+         T xs = x - 6;
+         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+         result = Y * x + R * x;
+      }
+      else if(x < 44)
+      {
+         // Max error found: 5.697761e-20
+         static const float Y = 0.99714565277099609375f;
+         static const T P[] = {    
+            -0.0024978212791898131227L,
+            -0.779190719229053954292e-5L,
+            0.254723037413027451751e-4L,
+            0.162397777342510920873e-5L,
+            0.396341011304801168516e-7L,
+            0.411632831190944208473e-9L,
+            0.145596286718675035587e-11L,
+            -0.116765012397184275695e-17L
+         };
+         static const T Q[] = {    
+            1L,
+            0.207123112214422517181L,
+            0.0169410838120975906478L,
+            0.000690538265622684595676L,
+            0.145007359818232637924e-4L,
+            0.144437756628144157666e-6L,
+            0.509761276599778486139e-9L
+         };
+         T xs = x - 18;
+         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+         result = Y * x + R * x;
+      }
+      else
+      {
+         // Max error found: 1.279746e-20
+         static const float Y = 0.99941349029541015625f;
+         static const T P[] = {    
+            -0.000539042911019078575891L,
+            -0.28398759004727721098e-6L,
+            0.899465114892291446442e-6L,
+            0.229345859265920864296e-7L,
+            0.225561444863500149219e-9L,
+            0.947846627503022684216e-12L,
+            0.135880130108924861008e-14L,
+            -0.348890393399948882918e-21L
+         };
+         static const T Q[] = {    
+            1L,
+            0.0845746234001899436914L,
+            0.00282092984726264681981L,
+            0.468292921940894236786e-4L,
+            0.399968812193862100054e-6L,
+            0.161809290887904476097e-8L,
+            0.231558608310259605225e-11L
+         };
+         T xs = x - 44;
+         T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+         result = Y * x + R * x;
+      }
+   }
+   return result;
+}
+
+template <class T>
+struct erf_roots
+{
+   std::tr1::tuple<T,T,T> operator()(const T& guess)
+   {
+      using namespace std;
+      T derivative = sign * (2 / sqrt(constants::pi<T>())) * exp(-(guess * guess));
+      T derivative2 = -2 * guess * derivative;
+      return std::tr1::make_tuple(((sign > 0) ? boost::math::erf(guess) : boost::math::erfc(guess)) - target, derivative, derivative2);
+   }
+   erf_roots(T z, int s) : target(z), sign(s) {}
+private:
+   T target;
+   int sign;
+};
+
+template <class T>
+T erf_inv_imp(const T& p, const T& q, const boost::mpl::int_<0>*)
+{
+   //
+   // Generic version, get a guess that's accurate to 64-bits (10^-19)
+   //
+   T guess = erf_inv_imp(p, q, static_cast<mpl::int_<64> const*>(0));
+   T result;
+   //
+   // If T has more bit's than 64 in it's mantissa then we need to iterate,
+   // otherwise we can just return the result:
+   //
+   if(tools::digits<T>() > 64)
+   {
+      if(p <= 0.5)
+      {
+         result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type>(p, 1), guess, static_cast<T>(0), tools::max_value<T>(), (tools::digits<T>() * 2) / 3);
+      }
+      else
+      {
+         result = tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type>(q, -1), guess, static_cast<T>(0), tools::max_value<T>(), (tools::digits<T>() * 2) / 3);
+      }
+   }
+   else
+   {
+      result = guess;
+   }
+   return result;
+}
+
+} // namespace detail
+
+template <class T>
+T erfc_inv(T z)
+{
+   //
+   // Begin by testing for domain errors, and other special cases:
+   //
+   if((z < 0) || (z > 2))
+      tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z);
+   if(z == 0)
+      return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
+   if(z == 2)
+      return -tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
+   //
+   // Normalise the input, so it's in the range [0,1], we will
+   // negate the result if z is outside that range.  This is a simple
+   // application of the erfc reflection formula: erfc(-z) = 2 - erfc(z)
+   //
+   T p, q, s;
+   if(z > 1)
+   {
+      q = 2 - z;
+      p = 1 - q;
+      s = -1;
+   }
+   else
+   {
+      p = 1 - z;
+      q = z;
+      s = 1;
+   }
+   //
+   // A bit of meta-programming to figure out which implementation
+   // to use, based on the number of bits in the mantissa of T:
+   //
+   typedef typename mpl::if_c<
+      std::numeric_limits<T>::is_specialized
+      &&
+      (std::numeric_limits<T>::digits <= 64),
+      mpl::int_<64>,
+      mpl::int_<0>
+   >::type tag_type;
+   //
+   // Likewise use internal promotion, so we evaluate at a higher
+   // precision internally if it's appropriate:
+   //
+   typedef typename tools::evaluation<T>::type eval_type;
+   //
+   // And get the result, negating where required:
+   //
+   return s * tools::checked_narrowing_cast<typename remove_cv<T>::type>(
+      detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), static_cast<tag_type const*>(0)), BOOST_CURRENT_FUNCTION);
+}
+
+template <class T>
+T erf_inv(T z)
+{
+   //
+   // Begin by testing for domain errors, and other special cases:
+   //
+   if((z < -1) || (z > 1))
+      tools::domain_error<typename remove_cv<T>::type>(BOOST_CURRENT_FUNCTION, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z);
+   if(z == 1)
+      return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
+   if(z == -1)
+      return -tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
+   if(z == 0)
+      return 0;
+   //
+   // Normalise the input, so it's in the range [0,1], we will
+   // negate the result if z is outside that range.  This is a simple
+   // application of the erf reflection formula: erf(-z) = -erf(z)
+   //
+   T p, q, s;
+   if(z < 0)
+   {
+      p = -z;
+      q = 1 - p;
+      s = -1;
+   }
+   else
+   {
+      p = z;
+      q = 1 - z;
+      s = 1;
+   }
+   //
+   // A bit of meta-programming to figure out which implementation
+   // to use, based on the number of bits in the mantissa of T:
+   //
+   typedef typename mpl::if_c<
+      std::numeric_limits<T>::is_specialized
+      &&
+      (std::numeric_limits<T>::digits <= 64),
+      mpl::int_<64>,
+      mpl::int_<0>
+   >::type tag_type;
+   //
+   // Likewise use internal promotion, so we evaluate at a higher
+   // precision internally if it's appropriate:
+   //
+   typedef typename tools::evaluation<T>::type eval_type;
+   //
+   // And get the result, negating where required:
+   //
+   return s * tools::checked_narrowing_cast<typename remove_cv<T>::type>(
+      detail::erf_inv_imp(static_cast<eval_type>(p), static_cast<eval_type>(q), static_cast<tag_type const*>(0)), BOOST_CURRENT_FUNCTION);
+}
+
+} // namespace math
+} // namespace boost
+
+#endif // BOOST_MATH_SF_ERF_INV_HPP
--- a/include/boost/math/special_functions/erf.hpp
+++ b/include/boost/math/special_functions/erf.hpp
@@ -82,65 +82,6 @@ T erf_asymptotic_limit()
   return erf_asymptotic_limit_N(boost::integral_constant<int, std::numeric_limits<T>::digits>());
 }

-//
-// Inverses, start with inverse of the Normal Distribution from
-// Abramowitz & Stegun p933, 26.2.23,
-// then convert that to the inverse of erfc:
-//
-template <class T>
-inline T invq(T q)
-{
-   using namespace std;
-
-   static const float num[4] =
-   {
-      2.515517f,
-      0.802853f,
-      0.010328f,
-      0.0f
-   };
-   static const float denom[4] =
-   {
-      1.0f,
-      1.432788f,
-      0.189269f,
-      0.001308f
-   };
-
-   T t = sqrt(-2 * log(q));
-   return t - tools::evaluate_rational(num, denom, t, 4);
-}
-
-template <class T>
-inline T estimate_inverse_erfc(T q)
-{
-   if(q == 2)
-      return -tools::max_value<T>();
-   if(q == 0)
-      return tools::max_value<T>();
-   if(q > 1)
-      return -estimate_inverse_erfc(2 - q);
-   if(q == 1)
-      return 0;
-   return T(0.70710678) * invq(q / 2);
-}
-
-template <class T>
-struct erf_roots
-{
-   std::tr1::tuple<T,T,T> operator()(const T& guess)
-   {
-      using namespace std;
-      T derivative = sign * (2 / sqrt(constants::pi<T>())) * exp(-(guess * guess));
-      T derivative2 = -2 * guess * derivative;
-      return std::tr1::make_tuple(((sign > 0) ? boost::math::erf(guess) : boost::math::erfc(guess)) - target, derivative, derivative2);
-   }
-   erf_roots(T z, int s) : target(z), sign(s) {}
-private:
-   T target;
-   int sign;
-};
-
 template <class T, class L, class Tag>
 T erf_imp(T z, bool invert, const L& l, const Tag& t)
 {
@@ -195,15 +136,6 @@ T erf_imp(T z, bool invert, const L& l, const Tag& t)
   return result;
 }

-template <class T>
-T truncate(T z, int b)
-{
-   using namespace std;
-   int exp;
-   z = frexp(z, &exp);
-   return ldexp(floor(ldexp(z, b)), exp - b);
-}
-
 template <class T, class L>
 T erf_imp(T z, bool invert, const L& l, const mpl::int_<53>& t)
 {
@@ -804,29 +736,11 @@ T erfc(T z)
      mpl::int_< ::std::numeric_limits<value_type>::digits>()), BOOST_CURRENT_FUNCTION);
 }

-template <class T>
-T erfc_inv(T z)
-{
-   if((z < 0) || (z > 2))
-      tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z);
-   T guess = detail::estimate_inverse_erfc(z);
-   return tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type>(z, -1), guess, -tools::max_value<T>(), tools::max_value<T>(), (tools::digits<T>() * 2) / 3);
-}
-
-template <class T>
-T erf_inv(T z)
-{
-   if((z < -1) || (z > 1))
-      tools::domain_error<typename remove_cv<T>::type>(BOOST_CURRENT_FUNCTION, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z);
-   T guess = detail::estimate_inverse_erfc(1 - z);
-   if((fabs(z) != 1) && (fabs(guess) == tools::max_value<typename remove_cv<T>::type>()))
-      guess = static_cast<T>((z < 0) ? -4 : 4);
-   return tools::halley_iterate(detail::erf_roots<typename remove_cv<T>::type>(z, 1), guess, ((z > 0) ? 0 : -tools::max_value<T>()), ((z < 0) ? 0 : tools::max_value<T>()), (tools::digits<T>() * 2) / 3);
-}
-
 } // namespace math
 } // namespace boost

+#include <boost/math/special_functions/detail/erf_inv.hpp>
+
 #endif // BOOST_MATH_SPECIAL_ERF_HPP


--- a/include/boost/math/tools/rational.hpp
+++ b/include/boost/math/tools/rational.hpp
@@ -8,6 +8,77 @@

 namespace boost{ namespace math{ namespace tools{

+namespace detail{
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V&, const mpl::int_<1>*)
+{
+   return a[0];
+}
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<2>*)
+{
+   return a[0] + x * a[1];
+}
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<3>*)
+{
+   return ((a[2] * x) + a[1]) * x + a[0];
+}
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<4>*)
+{
+   return ((a[3] * x + a[2]) * x + a[1]) * x + a[0];
+}
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<5>*)
+{
+   return (((a[4] * x + a[3]) * x + a[2]) * x + a[1]) * x + a[0];
+}
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<6>*)
+{
+   return ((((a[5] * x + a[4]) * x + a[3]) * x + a[2]) * x + a[1]) * x + a[0];
+}
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<7>*)
+{
+   return (((((a[6] * x + a[5]) * x + a[4]) * x + a[3]) * x + a[2]) * x + a[1]) * x + a[0];
+}
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<8>*)
+{
+   return ((((((a[7] * x + a[6]) * x + a[5]) * x + a[4]) * x + a[3]) * x + a[2]) * x + a[1]) * x + a[0];
+}
+
+template <class T, class V>
+inline V evaluate_polynomial_c_imp(const T* a, const V& x, const mpl::int_<9>*)
+{
+   return (((((((a[8] * x + a[7]) * x + a[6]) * x + a[5]) * x + a[4]) * x + a[3]) * x + a[2]) * x + a[1]) * x + a[0];
+}
+
+template <class T, class V, class Tag>
+inline V evaluate_polynomial_c_imp(const T* a, const V& val, const Tag*)
+{
+   return evaluate_polynomial(a, val, Tag::value);
+}
+
+} // namespace detail
+
+template <std::size_t N, class T, class V>
+inline V evaluate_polynomial(const T(&a)[N], const V& val)
+{
+   typedef mpl::int_<N> tag_type;
+   return detail::evaluate_polynomial_c_imp(static_cast<const T*>(a), val, static_cast<tag_type const*>(0));
+}
+
 template <class T, class U>
 U evaluate_polynomial(const T* poly, U z, std::size_t count)
 {
--- a/include/boost/math/tools/remez.hpp
+++ b/include/boost/math/tools/remez.hpp
@@ -148,6 +148,16 @@ public:
      --orderN;
      ++orderD;
   }
+   void rescale(T a, T b)
+   {
+      T scale = (b - a) / (max - min);
+      for(unsigned i = 0; i < control_points.size(); ++i)
+      {
+         control_points[i] = (control_points[i] - min) * scale + a;
+      }
+      min = a;
+      max = b;
+   }
 private:

   void init_chebyshev();
--- a/minimax/f.cpp
+++ b/minimax/f.cpp
@@ -1,5 +1,5 @@

-
+#define L22
 #include "../tools/ntl_rr_lanczos.hpp"
 #include <boost/math/tools/ntl.hpp>
 #include <boost/math/tools/polynomial.hpp>
@@ -88,18 +88,34 @@ T f0(T x, int variant)
      //
      if(x == 0)
         x = boost::math::tools::epsilon<T>() / 64;
-      return boost::math::erf_inv(x) / x;
+      return boost::math::erf_inv(x) / (x * (x+10));
   case 8:
      // 
-      //erfc_inv in range [0.5, 1]
+      // erfc_inv in range [0.25, 0.5]
+      // Use an x-offset of 0.25, and range [0, 0.25]
+      // abs error, auto y-offset.
      //
      if(x == 0)
-         return x = boost::lexical_cast<T>("1e-5000");
+         x = boost::lexical_cast<T>("1e-5000");
      return sqrt(-2 * log(x)) / boost::math::erfc_inv(x);
   case 9:
+      {
      if(x == 0)
-         return 0;
-      return 1 / boost::math::erfc_inv(x);
+         x = tiny;
+      double xr = boost::math::tools::real_cast<double>(x);
+      T z = 1/x;
+      double zr = boost::math::tools::real_cast<double>(z);
+      double yr = exp(zr * zr / -2);
+      T y = exp(z * z / -2);
+      return (boost::math::erfc_inv(y)/* - z*/ + log(z) / z) / z;
+      }
+   case 10:
+      {
+      if(x == 0)
+         x = boost::lexical_cast<T>("1e-5000");
+      T y = exp(-x*x); // sqrt(-log(x)) - 5;
+      return boost::math::erfc_inv(y) / x;
+      }
   }
   return 0;
 }
@@ -117,5 +133,16 @@ void show_extra(
   const NTL::RR& y_offset, 
   int variant)
 {
+   switch(variant)
+   {
+   case 9:
+      {
+         NTL::RR v = boost::math::tools::log_min_value<NTL::RR>();
+         v *= -2;
+         v = sqrt(v);
+         v = 1 / v;
+         std::cout << "Minimum range = " << v << std::endl;
+      }
+   }
 }

--- a/minimax/main.cpp
+++ b/minimax/main.cpp
@@ -47,6 +47,7 @@ NTL::RR the_function(const NTL::RR& val)
 void step_some(unsigned count)
 {
   try{
+      NTL::RR::SetPrecision(working_precision);
      if(!started)
      {
         //
@@ -108,6 +109,7 @@ void step(const char*, const char*)

 void show(const char*, const char*)
 {
+   NTL::RR::SetPrecision(working_precision);
   if(started)
   {
      boost::math::tools::polynomial<NTL::RR> n = p_remez->numerator();
@@ -151,10 +153,15 @@ void show(const char*, const char*)

      show_extra(n, d, x_offset, y_offset, variant);
   }
+   else
+   {
+      std::cerr << "Nothing to display" << std::endl;
+   }
 }

 void do_graph(unsigned points)
 {
+   NTL::RR::SetPrecision(working_precision);
   NTL::RR step = (b - a) / (points - 1);
   NTL::RR x = a;
   while(points > 1)
@@ -177,6 +184,7 @@ void graph(const char*, const char*)
 template <class T>
 void do_test(T, const char* name)
 {
+   NTL::RR::SetPrecision(working_precision);
   if(started)
   {
      //
@@ -274,6 +282,30 @@ void test_all(const char*, const char*)
   do_test((long double)(0), "long double");
 }

+void rotate(const char*, const char*)
+{
+   if(p_remez)
+   {
+      p_remez->rotate();
+   }
+   else
+   {
+      std::cerr << "Nothing to rotate" << std::endl;
+   }
+}
+
+void rescale(const char*, const char*)
+{
+   if(p_remez)
+   {
+      p_remez->rescale(a, b);
+   }
+   else
+   {
+      std::cerr << "Nothing to rescale" << std::endl;
+   }
+}
+
 int test_main(int, char* [])
 {
   std::string line;
@@ -337,6 +369,10 @@ int test_main(int, char* [])
            str_p("test") && str_p("long")[&test_long]
      ||
            str_p("test") && str_p("all")[&test_all]
+      ||
+            str_p("rotate")[&rotate]
+      ||
+            str_p("rescale") && real_p[assign_a(a)] && real_p[assign_a(b)] && epsilon_p[&rescale]

         ), space_p).full)
      {
--- a/test/erfc_inv_big_data.ipp
+++ b/test/erfc_inv_big_data.ipp
@@ -0,0 +1,206 @@
+#define SC_(x) static_cast<T>(BOOST_JOIN(x, L))
+   static const boost::array<boost::array<T, 2>, 200> erfc_inv_big_data = {
+      SC_(0.5460825444184401149953083908674803067585e-4312), SC_(99.62016927389407649911084501709563799849), 
+      SC_(0.2645475979732877501542006305722535728925e-4291), SC_(99.38083815930157675670279543643240126733), 
+      SC_(0.4647218972549971591905525747707221131215e-4281), SC_(99.26209205835381203650501813598769577392), 
+      SC_(0.1525121257807440558802279626797226161005e-4243), SC_(98.82602533559047574603223381461761500878), 
+      SC_(0.1981708500756372844759510276075956268897e-4234), SC_(98.71980149740490464004477100962342642155), 
+      SC_(0.3841754068409692466399340503654464602243e-4180), SC_(98.08467819532462847486581213069198608289), 
+      SC_(0.5335004017651487579647125338582580235046e-4091), SC_(97.03275964437690000965320992419441429924), 
+      SC_(0.1806383895924070461510470704435189756123e-4076), SC_(96.86022098874718527180735862424264805433), 
+      SC_(0.1321513657594592521699131916159148731226e-4066), SC_(96.7429083707271579347385166737183074303), 
+      SC_(0.1054920495502343332517378247415093412624e-4043), SC_(96.46999015250599733880053488460821308812), 
+      SC_(0.47363892632169499253148937212672347289e-4043), SC_(96.46220643869959630846141605535921442023), 
+      SC_(0.4750819439081903839348296656606873483722e-4000), SC_(95.94763372873564880985964950938277645413), 
+      SC_(0.5934696803923580795229088781443547003391e-3998), SC_(95.92247396394370347463100485254034797278), 
+      SC_(0.2809661603171360146097006130839773212914e-3996), SC_(95.90236599111862577117586443056337455754), 
+      SC_(0.201344506190345013258215914381814651279e-3981), SC_(95.72387427307012892119049758800169512382), 
+      SC_(0.1349680361063247504563298949422573179728e-3978), SC_(95.68987765258188263846875190956134957764), 
+      SC_(0.9554154264888789862328696045397641686626e-3961), SC_(95.47488611856686267812468521780070611561), 
+      SC_(0.4565370419999699854685454921616341777152e-3921), SC_(94.99523139627283958425843803677334678489), 
+      SC_(0.9326136640408284179797117117017987832624e-3838), SC_(93.98018858040542889098559401088499617323), 
+      SC_(0.2085705673165100591822581196874037577838e-3788), SC_(93.37371566951938600144821032880068555931), 
+      SC_(0.28906230618616666278153940715368369251e-3788), SC_(93.37196812654735398207358428222836165935), 
+      SC_(0.540152632181355759951616925660164590124e-3750), SC_(92.89890240230691810677222482470759338646), 
+      SC_(0.1373243759745268803582425259844487470076e-3682), SC_(92.05981131540500971264194940023399842876), 
+      SC_(0.3331647658738942966983713783383120828716e-3653), SC_(91.69161143487247689253830777510103757803), 
+      SC_(0.134842155208242065478707324403653741427e-3641), SC_(91.54576292987357613504022634708802175022), 
+      SC_(0.7492191291954483620899871520373510239375e-3623), SC_(91.30973569873630381671516956816451328373), 
+      SC_(0.2551662007977519576358324440130450442778e-3575), SC_(90.70847518912252606205543434018735665529), 
+      SC_(0.6190686739111601149790106116510494551314e-3563), SC_(90.55115617451523669098548678613430984841), 
+      SC_(0.9276632609820886396580236856359763026966e-3447), SC_(89.06191184952516547664083098302741022803), 
+      SC_(0.6433072193358217453796447947689324967478e-3351), SC_(87.8143275613841604460220844122740180476), 
+      SC_(0.3653541017426498234750653436978737170327e-3338), SC_(87.64696333406539636312319479272468979115), 
+      SC_(0.2789596742720445243557533363147985407484e-3305), SC_(87.21398755358298377575097222058344254911), 
+      SC_(0.443825473802222651705813573190067261172e-3131), SC_(84.88340554604607048910870088690204110131), 
+      SC_(0.3345784652046112831225855034049282590715e-3126), SC_(84.81723263217786918755009758623351862559), 
+      SC_(0.2414445336524691790354556146640969706061e-3069), SC_(84.0419598414264046587129844579410948894), 
+      SC_(0.8478603490304246265311990779170856959861e-3028), SC_(83.47092766845828325139370216002957076776), 
+      SC_(0.14740427247445948990191590432902058536e-2997), SC_(83.05281564078241081289202980345540186102), 
+      SC_(0.1329036517699801662619770631013011165738e-2944), SC_(82.31552494008283774947792893056404158871), 
+      SC_(0.1706141285384843345987169410441595650031e-2933), SC_(82.16002264325748303060260330185689090281), 
+      SC_(0.8373844481146949976640588846233709740831e-2929), SC_(82.09426838736262780476565867158939922013), 
+      SC_(0.6495088549453856516509184167801252067894e-2926), SC_(82.05373670564597872871969963977451713693), 
+      SC_(0.1833720938105362893360888392913141441369e-2882), SC_(81.44184518033949369735162216055497263014), 
+      SC_(0.2333520177487467820415573739119469955902e-2879), SC_(81.39794775935971788819430665383756478382), 
+      SC_(0.2894427653993793412567889456980821259226e-2847), SC_(80.94277839157518781341184937011412770395), 
+      SC_(0.2604952067990663450155633194749904382233e-2836), SC_(80.78683208931494064074748747620852519014), 
+      SC_(0.8678138849264725191436450361584810016324e-2832), SC_(80.72235938139247087444254555436435678849), 
+      SC_(0.3682667903200544588455731527142145038363e-2781), SC_(79.99708530932539906643334265914340895581), 
+      SC_(0.4588121885193651068401508885826145226945e-2766), SC_(79.77955734877589476980263452200168409691), 
+      SC_(0.7764413237799866858364453900343995777469e-2765), SC_(79.76182876388495108788138224124277407448), 
+      SC_(0.2423421525461451159440550098196593570056e-2760), SC_(79.69693586691073000502029265856910229446), 
+      SC_(0.5878986254507577903335876253702952707265e-2746), SC_(79.48887849840574283826340114327272440985), 
+      SC_(0.1673845625470905295897015228985949168928e-2737), SC_(79.36634262823360407978363436458177566234), 
+      SC_(0.6698391998725201197546429856824581586977e-2735), SC_(79.32858818851906731315743566169843426962), 
+      SC_(0.2515130355594925802647014162660517938496e-2725), SC_(79.18952171986121353050089556124781639244), 
+      SC_(0.2461647325252598254899739145066165462755e-2636), SC_(77.88509727460664014637457387303638162253), 
+      SC_(0.1641105794585019759166735829835878355641e-2634), SC_(77.85813389734403078489439365068353770847), 
+      SC_(0.6351595345143555864526860293158388570963e-2601), SC_(77.35991719682138597284770268761692627695), 
+      SC_(0.9223028410198721360701049566280075593579e-2545), SC_(76.51960497818744876859269133142591483997), 
+      SC_(0.5002814464212067387766735395165041257035e-2533), SC_(76.34286436673947487064698531398136817799), 
+      SC_(0.1185824625931586847282481488476410873355e-2519), SC_(76.14091514413348192509930444973601480465), 
+      SC_(0.257414900548022702014576953401369926289e-2493), SC_(75.74167889782762540410788473501883567486), 
+      SC_(0.9183566141599824732209230606531307948203e-2490), SC_(75.68766731518135433104337024242062602941), 
+      SC_(0.2932809029313999449803475825227510422988e-2473), SC_(75.43622352043665347320282000269125370548), 
+      SC_(0.5498724163832736366794865810754062391466e-2467), SC_(75.34043418668345319181733147217173430757), 
+      SC_(0.6795395146348654066574869189625572073492e-2451), SC_(75.09414895042043847993638837712888889371), 
+      SC_(0.4232772623968377298803745728444220344828e-2398), SC_(74.2804044020029126990405093970664224902), 
+      SC_(0.437705010214781194544970904918846009988e-2382), SC_(74.0317968899859238324605680410762845653), 
+      SC_(0.7328784181257135130009299062120147357384e-2332), SC_(73.24665809155865681149475856856425835323), 
+      SC_(0.8195784402152374810490855567845040936044e-2295), SC_(72.6620497559170838922485865519610574788), 
+      SC_(0.8987013844942659551846558843047951168834e-2288), SC_(72.55042905061235174343554742739711560742), 
+      SC_(0.2390567310984670056453397739020313330694e-2193), SC_(71.03636046035780611691848769848245798184), 
+      SC_(0.4332167143233490736155663100970351222456e-2185), SC_(70.90240600026913612396368503494705740935), 
+      SC_(0.108890416531056942890760795215821213854e-2164), SC_(70.57040795030909610470208299858512551923), 
+      SC_(0.5060627104045047637085866910481125331788e-2163), SC_(70.54320633525883330072442815574897293795), 
+      SC_(0.3691223561779676390769470042383461854014e-2122), SC_(69.87319260661719659265781125754219708719), 
+      SC_(0.6969547608629390600217634449465487521466e-2099), SC_(69.48865744111861740557954052160177140024), 
+      SC_(0.1920331015960815466646483904103550820861e-2096), SC_(69.44822094478638613766926468731012767531), 
+      SC_(0.6083975528963974425124099594841117907605e-2016), SC_(68.1007710158724367714118546612371304706), 
+      SC_(0.7700980458334894771243720474212758259757e-2007), SC_(67.94673150215833308050516745181191739631), 
+      SC_(0.3784507382437993534366757035181644553005e-2006), SC_(67.93501557417646729265187361704446762791), 
+      SC_(0.1933058527290172128745148456282632896472e-2000), SC_(67.8382198195958948309248955546035843938), 
+      SC_(0.3073129169339213242656468129490128120275e-1912), SC_(66.32462247762686573271005635222158246919), 
+      SC_(0.3563704808140633534254743351059329459844e-1867), SC_(65.53779859486217889665515029788703171428), 
+      SC_(0.1117328435404434848431164536157383381954e-1865), SC_(65.5115114749327337441199600755950731015), 
+      SC_(0.162291367493531393505178202172746460504e-1863), SC_(65.47350810130320090152901346075787089108), 
+      SC_(0.7291877778593992272357923862974438803504e-1814), SC_(64.59461716634387479104306278695990989692), 
+      SC_(0.1103697274312455454535522513124420450513e-1800), SC_(64.35930484119184106703256655579503018943), 
+      SC_(0.2787184985694988691274106077027481674296e-1777), SC_(63.93935283526837024038868862230807595467), 
+      SC_(0.2242899394948684173448369886336215330859e-1736), SC_(63.19860571879709151914226858360392775225), 
+      SC_(0.6373703493315477830264234158158397230918e-1648), SC_(61.5663708938977044791256528616543682363), 
+      SC_(0.1130761791970139872085460502062402383311e-1534), SC_(59.4111793960172436477860564668325556624), 
+      SC_(0.5186477168932786698797072956323906176672e-1526), SC_(59.24311964529112251267296011750890265511), 
+      SC_(0.2371756563156473396673371175154689282089e-1487), SC_(58.48710687097943040207443293838886041518), 
+      SC_(0.2474530000251364788339443399087025927602e-1487), SC_(58.48674428195209156566913277116595315737), 
+      SC_(0.1620952363661926127361288574589236587446e-1465), SC_(58.05577331395656272747602508275617390989), 
+      SC_(0.2122625807107458215797391334648973459688e-1429), SC_(57.33517766371867417445655972866006304427), 
+      SC_(0.6243197332907990298222556346398382173979e-1405), SC_(56.8418011386019473870793880887159457241), 
+      SC_(0.5516957459853521159388169416018722269647e-1400), SC_(56.74154440911648902947449827969373477014), 
+      SC_(0.2926235914834555267234269989321844626118e-1369), SC_(56.11477468249401615221533062929223131088), 
+      SC_(0.3794627523015867743608205312364027747225e-1312), SC_(54.93070214218908343008136754660768021631), 
+      SC_(0.318880483330771111942926508195676807607e-1300), SC_(54.68024803482526161939511993227336637443), 
+      SC_(0.3369469880227599686172744089074258905398e-1300), SC_(54.67974419399797522074814528351877319048), 
+      SC_(0.1302713766629726626489873195118034634277e-1298), SC_(54.64631892482891526341853652123547792494), 
+      SC_(0.2469170628653646881926626888059718255132e-1223), SC_(53.03693750856997083336291247307325166587), 
+      SC_(0.1371607091728437802604087975394383012252e-1207), SC_(52.69411504628090115129531887668248361258), 
+      SC_(0.1731629169156886034424126010604259325511e-1111), SC_(50.55126123375417963154846480444273825353), 
+      SC_(0.1988578488030194560097782052960341565225e-1101), SC_(50.32166883510982608447898399736737553384), 
+      SC_(0.2274776484360344132543141099781042532126e-1078), SC_(49.79143524779148474319978052812994435046), 
+      SC_(0.5627666511768513324968956052244126556016e-1045), SC_(49.01338231598138050658354300270473769631), 
+      SC_(0.1825674918577823738426304992234296317652e-1041), SC_(48.93085667837639283695386917651720233174), 
+      SC_(0.3802925527627705362359485795331572740116e-1032), SC_(48.71115012911996883239788333878826905934), 
+      SC_(0.704333512302902280929638896129066471028e-1023), SC_(48.49165994415978768781022851843247093366), 
+      SC_(0.1313561715183354223096034930474262234265e-978), SC_(47.429170762115701641491822625509511735), 
+      SC_(0.1479971733547495818402245675295016370485e-972), SC_(47.28207432819548183826826245137119654147), 
+      SC_(0.1347434440250964238000754812764303681229e-956), SC_(46.89195339639151403964904162145015071721), 
+      SC_(0.1121116975652110010815188524143927325837e-910), SC_(45.75090030899601472255890798805286514432), 
+      SC_(0.6537208325713531384544436742840979538836e-903), SC_(45.55510822465697895367331173242574718725), 
+      SC_(0.5260776239804675672628897143432967891746e-862), SC_(44.50957576907949796448736564849939576426), 
+      SC_(0.9950661094879265377710346605980409852453e-836), SC_(43.8248018686713805421278904903492358615), 
+      SC_(0.2062512012152859485600938747081581263563e-777), SC_(42.26547837659956767422521162018458102221), 
+      SC_(0.1860519246978454733361447959343549701547e-773), SC_(42.15763210320445259646388169147606822404), 
+      SC_(0.292259268934242545722820123169012540872e-772), SC_(42.12496313029953419941505368125427033208), 
+      SC_(0.2655070971534510004697091453590968778995e-739), SC_(41.21462174108933535284728637920523260807), 
+      SC_(0.2901429667218951401071997769907995675312e-725), SC_(40.82070209799865627866597450674844815644), 
+      SC_(0.2691466393457144658729027931985645462202e-709), SC_(40.36799005611971249179085189213378251021), 
+      SC_(0.3197098867400748771350059689777583564847e-690), SC_(39.82043439287531202279866114834532850691), 
+      SC_(0.509206568051986013963235738541488132177e-677), SC_(39.43700823628539555409237521877515418161), 
+      SC_(0.3582760294415859512291360611861492812312e-668), SC_(39.17795922357671427001830102855480201337), 
+      SC_(0.7456401788156278694200104791854742334257e-668), SC_(39.16860717559069454303739314459059279729), 
+      SC_(0.1879442582827104576451480885999452191966e-655), SC_(38.80249559566898835529327655718231853047), 
+      SC_(0.3684861623635251188317018685920580219909e-651), SC_(38.67497066678513891028432141977632830295), 
+      SC_(0.5300532049641075447010094275260357817648e-632), SC_(38.10059865129674828052108905587617617209), 
+      SC_(0.254919090870011244639387672109478129648e-628), SC_(37.98921207590422610976320034158408446716), 
+      SC_(0.1261526994355170599813296346837863228294e-566), SC_(36.07180567207263624533709617168647173026), 
+      SC_(0.3348102113183452824057762060804616114175e-551), SC_(35.57631428343251267029769551692709310927), 
+      SC_(0.3615878633639651283997811017708582782249e-540), SC_(35.21759342981655987654847385348063502236), 
+      SC_(0.8556456214239050595803534051499578042028e-537), SC_(35.10716335050071391896400132832542452084), 
+      SC_(0.109310372545754088817907678327638735624e-531), SC_(34.93937406241154814427167189866546786437), 
+      SC_(0.6542713735351179620273602253478996592612e-503), SC_(33.97831683163889299935986132301256667123), 
+      SC_(0.2612280308120653787093215442892321688194e-493), SC_(33.65156678099809689208816819469703821624), 
+      SC_(0.1022554993887925308488446602367781896336e-479), SC_(33.18348488984914247206612191269120081123), 
+      SC_(0.2031774862883333653515762991534010669034e-436), SC_(31.64637875978842562095759897711662201891), 
+      SC_(0.6588827942394778508739583209849917530674e-435), SC_(31.59139083855039387728200786192864837174), 
+      SC_(0.3663469071718277378831360961951773982499e-412), SC_(30.75175430540050080137890276895073229866), 
+      SC_(0.1361305298482799765369196276048348425522e-395), SC_(30.12535428053783863012447180563559906292), 
+      SC_(0.3744171560141971093473580527502555697533e-390), SC_(29.91687192152261750809710487517869880922), 
+      SC_(0.6874405809297488106605573621300258853644e-383), SC_(29.63617622599272249872842949718868018609), 
+      SC_(0.1459704391344155418984387882925473942446e-354), SC_(28.51518995595697899285114107951299228827), 
+      SC_(0.143123157725799928952051937462499680388e-336), SC_(27.77976460222037023739857761278401004959), 
+      SC_(0.9025996658203895649285047801834901806925e-333), SC_(27.62194204731470015756869063464069373815), 
+      SC_(0.1293287994487535182921055345118992149121e-330), SC_(27.53198333089873848069767951797741434292), 
+      SC_(0.1040364789689237257780791904879173249974e-287), SC_(25.67656975866313564221220720210665888424), 
+      SC_(0.2450719147172718682892603854216257656211e-282), SC_(25.43473899116627259523415573506406682191), 
+      SC_(0.1540635014300999959290867576988859893076e-280), SC_(25.35326739382267978364899192426764093558), 
+      SC_(0.1611739209855438499046222515484436760851e-273), SC_(25.03273169525729329624232590800413859205), 
+      SC_(0.2781078901633922872827030068281452871724e-257), SC_(24.27512000889987237167959986465966890736), 
+      SC_(0.8837820681846044382839818307987545000847e-230), SC_(22.93495456229819340317355998004446065241), 
+      SC_(0.2845511846316072211634522373099423142668e-217), SC_(22.29887672117465208331735160799934601323), 
+      SC_(0.662011186385949636511655220013184768187e-212), SC_(22.02033498263166080103335807537366542004), 
+      SC_(0.1113036584088835972703336205837326051999e-204), SC_(21.63966022176108818726656867477709317659), 
+      SC_(0.891278753475701910699424577762349195273e-195), SC_(21.10677883438545215904960900943812596651), 
+      SC_(0.1243375606125932867518669565282002251187e-191), SC_(20.93474611896499049089384377937974962862), 
+      SC_(0.963481566993839758726177317091390488753e-188), SC_(20.72000418866581389344212197908962674052), 
+      SC_(0.1401265332225531798315780528328964001988e-167), SC_(19.56909048260682534539216053987194710384), 
+      SC_(0.7560947579095629715718115752684784187336e-154), SC_(18.74494560773490096872968649468667611647), 
+      SC_(0.2179858592534951368829782528927445425064e-151), SC_(18.59346814776643084908600446605555382182), 
+      SC_(0.5208405693419352209961246987938629746382e-131), SC_(17.28776816936911830963159975455855394785), 
+      SC_(0.5376896865772724579278898126622880802843e-127), SC_(17.01882418428649428645068852272615582504), 
+      SC_(0.1449951754876351411423239477096839160613e-114), SC_(16.15763261758311680195806763260744183267), 
+      SC_(0.3189212836315126918184908007667112228414e-108), SC_(15.70012576642621709878560027430710336585), 
+      SC_(0.1297369886602648502231884634897739722895e-107), SC_(15.65546665388935777226018284593032382368), 
+      SC_(0.2602821327589283546034644516083144727197e-103), SC_(15.33647743152077165227354683273671003413), 
+      SC_(0.1590035594081776210615338926776093483999e-96), SC_(14.81946110818452708436758045284808001241), 
+      SC_(0.2720908074515874075065365442027770927174e-96), SC_(14.80136596875698966601900922593357151598), 
+      SC_(0.3706750818860597777102228512408751197629e-87), SC_(14.07473190298865115531360325271322276549), 
+      SC_(0.3489539664837006993042769801460156282803e-83), SC_(13.74669364199157990777768669458685265472), 
+      SC_(0.1947620269151566018960312643171165602404e-79), SC_(13.43009987431500967622734624289673694188), 
+      SC_(0.5129832297533069185312436951335270626736e-73), SC_(12.86957639029561313119037589037764922829), 
+      SC_(0.9272206234111032518003166213125699095688e-70), SC_(12.57574048509944911276422466059543556873), 
+      SC_(0.1460449453933976163637836336160553873706e-68), SC_(12.46599637375040876166599138465603014844), 
+      SC_(0.8768204308820258161647995326423963318019e-68), SC_(12.39412855962299302779918858089973570208), 
+      SC_(0.3176564512338672251415481009303329209995e-61), SC_(11.77127316747622053532997731227084711173), 
+      SC_(0.1172937774279091349820752991146331290881e-54), SC_(11.11297358832626209702547835299261897395), 
+      SC_(0.3822700066028246575792241970441049035743e-52), SC_(10.85058828114725350216779162566256390403), 
+      SC_(0.4793602326185941396320885530096873382082e-41), SC_(9.607343723542133339487829511083057069273), 
+      SC_(0.106201511009181306109371950926034372112e-40), SC_(9.566077632118694585964030456246300276422), 
+      SC_(0.1026229395621872573137267432392849554594e-24), SC_(7.413111396559360630306454710337113823294), 
+      SC_(0.1433834013674021626308791056277965507134e-19), SC_(6.574536501026477531602186421374903506051), 
+      SC_(0.2287484604245551701383736767135709586086e-12), SC_(5.183672646329402197657171322814515766894), 
+      SC_(0.4054382118091501052221361988694253843776e-11), SC_(4.903978376537783549083010816442468785255), 
+      SC_(0.5517250582467272934968166422809728429772e-9), SC_(4.386634719241854622420456791837035899963), 
+      SC_(0.1708920553446964553659243704866375034024e-5), SC_(3.383585107532389209677350898502297644278), 
+      SC_(0.2892930668406028955896934416060055101499e-5), SC_(3.308042643013470022433136596357284983766), 
+      SC_(0.7029998339557698614154475319387477445787e-5), SC_(3.17687395056466498458108882620439294841), 
+      SC_(0.3973417837559780922949941034264323413971e-4), SC_(2.90551585834893682957467771364140403613), 
+      SC_(0.04246334897368917726269942070072005435577), SC_(1.434684541881674882915735690265484394337), 
+      SC_(0.159920364968560744996532371753339418774), SC_(0.9937250495458864822841753804433979428041), 
+      SC_(0.2425390104583346613758947085681683120129), SC_(0.826370276571483139432927467391341431105), 
+      SC_(0.7954739362140872788414606986417965117653), SC_(0.1832884885283639734589195991287031281603), 
+      SC_(0.9236374866673857278562449757419727802699), SC_(0.06777816075457032123722158048171767182127), 
+   };
+#undef SC_
+
+
--- a/test/test_erf.cpp
+++ b/test/test_erf.cpp
@@ -232,6 +232,12 @@ void test_erf(T, const char* name)

   do_test_erfc_inv(erfc_inv_data, name, "Inverse Erfc Function");

+#  include "erfc_inv_big_data.ipp"
+
+   if(std::numeric_limits<T>::min_exponent <= -4500)
+   {
+      do_test_erfc_inv(erfc_inv_big_data, name, "Inverse Erfc Function: extreme values");
+   }
 }

 template <class T>
--- a/tools/erf_data.cpp
+++ b/tools/erf_data.cpp
@@ -4,6 +4,7 @@
 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

 #include <boost/math/tools/ntl.hpp>
+#include <boost/test/included/test_exec_monitor.hpp>
 #include <boost/math/special_functions/gamma.hpp>
 #include <boost/math/special_functions/erf.hpp> // for inverses
 #include <boost/math/constants/constants.hpp>
@@ -17,6 +18,19 @@
 using namespace boost::math::tools;
 using namespace std;

+float external_f;
+float force_truncate(const float* f)
+{
+   external_f = *f;
+   return external_f;
+}
+
+float truncate_to_float(NTL::RR r)
+{
+   float f = boost::math::tools::real_cast<float>(r);
+   return force_truncate(&f);
+}
+
 struct erf_data_generator
 {
   std::tr1::tuple<NTL::RR, NTL::RR> operator()(NTL::RR z)
@@ -105,14 +119,18 @@ void asymptotic_limit(int Bits)
      << terms << " terms." << std::endl;
 }

-NTL::RR erfc_inv(NTL::RR r)
+std::tr1::tuple<NTL::RR, NTL::RR> erfc_inv(NTL::RR r)
 {
-   std::cout << r << std::endl;
-   return boost::math::erfc_inv(r);
+   NTL::RR x = exp(-r * r);
+   x = NTL::RoundToPrecision(x, 64);
+   std::cout << x << "   ";
+   NTL::RR result = boost::math::erfc_inv(x);
+   std::cout << result << std::endl;
+   return std::tr1::make_tuple(x, result);
 }


-int main(int argc, char* argv[])
+int test_main(int argc, char*argv [])
 {
   NTL::RR::SetPrecision(1000);
   NTL::RR::SetOutputPrecision(40);
@@ -147,14 +165,19 @@ int main(int argc, char* argv[])
      }
      else if(strcmp(argv[1], "--erfc_inv") == 0)
      {
-         NTL::RR (*f)(NTL::RR);
-         f = /*boost::math::*/ erfc_inv;
+         std::tr1::tuple<NTL::RR, NTL::RR> (*f)(NTL::RR);
+         f = erfc_inv;
         std::cout << "Welcome.\n"
            "This program will generate spot tests for the inverse erfc function:\n";
+         std::cout << "Enter the maximum *result* expected from erfc_inv: ";
+         double max_val;
+         std::cin >> max_val;
         std::cout << "Enter the number of data points: ";
         int points;
         std::cin >> points;
-         data.insert(f, make_random_param(NTL::RR(0), NTL::RR(2), points));
+         parameter_info<NTL::RR> arg = make_random_param(NTL::RR(0), NTL::RR(max_val), points);
+         arg.type |= dummy_param;
+         data.insert(f, arg);
      }
   }
   else
--- a/tools/ibeta_invab_data.cpp
+++ b/tools/ibeta_invab_data.cpp
@@ -0,0 +1,96 @@
+//  (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 <boost/math/tools/ntl.hpp>
+#include <boost/test/included/test_exec_monitor.hpp>
+#include <boost/math/special_functions/beta.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/tools/test.hpp>
+#include <boost/lexical_cast.hpp>
+#include <fstream>
+
+#include <boost/math/tools/test_data.hpp>
+#include "ntl_rr_lanczos.hpp"
+
+using namespace boost::math::tools;
+
+//
+// Force trunctation to float precision of input values:
+// we must ensure that the input values are exactly representable
+// in whatever type we are testing, or the output values will all
+// be thrown off:
+//
+float external_f;
+float force_truncate(const float* f)
+{
+   external_f = *f;
+   return external_f;
+}
+
+float truncate_to_float(NTL::RR r)
+{
+   float f = boost::math::tools::real_cast<float>(r);
+   return force_truncate(&f);
+}
+
+std::tr1::mt19937 rnd;
+std::tr1::uniform_real<float> ur_a(1.0F, 5.0F);
+std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<float> > gen(rnd, ur_a);
+std::tr1::uniform_real<float> ur_a2(0.0F, 100.0F);
+std::tr1::variate_generator<std::tr1::mt19937, std::tr1::uniform_real<float> > gen2(rnd, ur_a2);
+
+struct ibeta_inv_data_generator
+{
+   std::tr1::tuple<NTL::RR, NTL::RR, NTL::RR, NTL::RR, NTL::RR, NTL::RR, NTL::RR> operator()
+      (NTL::RR bp, NTL::RR x_, NTL::RR p_)
+   {
+      float b = truncate_to_float(real_cast<float>(gen() * pow(NTL::RR(10), bp))); 
+      float x = truncate_to_float(real_cast<float>(x_));
+      float p = truncate_to_float(real_cast<float>(p_));
+      std::cout << b << " " << x << " " << p << std::flush;
+      NTL::RR inv = boost::math::ibeta_inva(NTL::RR(b), NTL::RR(x), NTL::RR(p));
+      std::cout << " " << inv << std::flush;
+      NTL::RR invc = boost::math::ibetac_inva(NTL::RR(b), NTL::RR(x), NTL::RR(p));
+      std::cout << " " << invc << std::endl;
+      NTL::RR invb = boost::math::ibeta_invb(NTL::RR(b), NTL::RR(x), NTL::RR(p));
+      std::cout << " " << invb << std::flush;
+      NTL::RR invbc = boost::math::ibetac_invb(NTL::RR(b), NTL::RR(x), NTL::RR(p));
+      std::cout << " " << invbc << std::endl;
+      return std::tr1::make_tuple(b, x, p, inv, invc, invb, invbc);
+   }
+};
+
+int test_main(int argc, char*argv [])
+{
+   NTL::RR::SetPrecision(1000);
+   NTL::RR::SetOutputPrecision(100);
+
+   bool cont;
+   std::string line;
+
+   parameter_info<NTL::RR> arg1, arg2, arg3;
+   test_data<NTL::RR> data;
+
+   std::cout << "Welcome.\n"
+      "This program will generate spot tests for the inverse incomplete beta function:\n"
+      "  ibeta_inva(a, p) and ibetac_inva(a, q)\n\n";
+
+   arg1 = make_periodic_param(NTL::RR(-5), NTL::RR(6), 11);
+   arg2 = make_random_param(NTL::RR(0.0001), NTL::RR(1), 10);
+   arg3 = make_random_param(NTL::RR(0.0001), NTL::RR(1), 10);
+
+   arg1.type |= dummy_param;
+   arg2.type |= dummy_param;
+   arg3.type |= dummy_param;
+
+   data.insert(ibeta_inv_data_generator(), arg1, arg2, arg3);
+
+   line = "ibeta_inva_data.ipp";
+   std::ofstream ofs(line.c_str());
+   write_code(ofs, data, "ibeta_inva_data");
+   
+   return 0;
+}
+