[PPL-devel] [GIT] ppl/ppl(master): Documented rational sqrt precision and implemented a more precise variant for numbers < 1.

[Abramo Bagnara]

Module: ppl/ppl
Branch: master
Commit: b98bfa9b31c0861653c862e025372ccf5f44ae21
URL:    http://www.cs.unipr.it/git/gitweb.cgi?p=ppl/ppl.git;a=commit;h=b98bfa9b31c0861653c862e025372ccf5f44ae21

Author: Abramo Bagnara <abramo.bagnara at gmail.com>
Date:   Sat Mar 28 13:59:39 2009 +0100

Documented rational sqrt precision and implemented a more precise variant for numbers < 1.

---

 TODO                       |    2 ++
 src/Init.cc                |    4 ++--
 src/checked_mpq.inlines.hh |   27 +++++++++++++++++----------
 3 files changed, 21 insertions(+), 12 deletions(-)

diff --git a/TODO b/TODO
index 98df6e6..0cff626 100644
--- a/TODO
+++ b/TODO
@@ -12,6 +12,8 @@ Enhancements for PPL 0.10.1 or later versions
 Enhancements for PPL 0.11 or later versions
 ===========================================
 
+- Rename set_rational_sqrt_precision_parameter to
+  set_irrational_precision and make it available on all interfaces.
 - Intervals are best instantiated with checked numbers with
   particular policies: review all the interfaced boxes,
   augment the testsuite, and update the documentation.
diff --git a/src/Init.cc b/src/Init.cc
index 373ed36..dce7393 100644
--- a/src/Init.cc
+++ b/src/Init.cc
@@ -86,8 +86,8 @@ PPL::Init::Init() {
     old_rounding_direction = fpu_get_rounding_direction();
     fpu_set_rounding_direction(round_fpu_dir(ROUND_DIRECT));
 #endif
-    // FIXME(0.10.1): is 3200 a magic number?
-    set_rational_sqrt_precision_parameter(3200);
+    // FIXME(0.10.1): choose a better default
+    set_rational_sqrt_precision_parameter(128);
   }
 }
 
diff --git a/src/checked_mpq.inlines.hh b/src/checked_mpq.inlines.hh
index d613489..c4bfef2 100644
--- a/src/checked_mpq.inlines.hh
+++ b/src/checked_mpq.inlines.hh
@@ -386,16 +386,22 @@ inline Result
 sqrt_mpq(mpq_class& to, const mpq_class& from, Rounding_Dir dir) {
   if (CHECK_P(To_Policy::check_sqrt_neg, from < 0))
     return assign_special<To_Policy>(to, V_SQRT_NEG, ROUND_IGNORE);
-  const unsigned long k = rational_sqrt_precision_parameter;
-  mpz_class& to_num = to.get_num();
-  mul2exp<To_Policy, From_Policy>(to_num, from.get_num(), 2*k, dir);
+  if (from == 0) {
+    to = 0;
+    return V_EQ;
+  }
+  bool gt1 = from.get_num() > from.get_den();
+  const mpz_class& from_a = gt1 ? from.get_num() : from.get_den();
+  const mpz_class& from_b = gt1 ? from.get_den() : from.get_num();
+  mpz_class& to_a = gt1 ? to.get_num() : to.get_den();
+  mpz_class& to_b = gt1 ? to.get_den() : to.get_num();
+  Rounding_Dir rdir = gt1 ? dir : inverse(dir);
+  mul2exp<To_Policy, From_Policy>(to_a, from_a, 2*rational_sqrt_precision_parameter, ROUND_IGNORE);
   Result rdiv
-    = div<To_Policy, To_Policy, To_Policy>(to_num,
-					   to_num, from.get_den(), dir);
-  Result rsqrt = sqrt<To_Policy, To_Policy>(to_num, to_num, dir);
-  mpz_class& to_den = to.get_den();
-  to_den = 1;
-  mul2exp<To_Policy, To_Policy>(to_den, to_den, k, dir);
+    = div<To_Policy, To_Policy, To_Policy>(to_a, to_a, from_b, rdir);
+  Result rsqrt = sqrt<To_Policy, To_Policy>(to_a, to_a, rdir);
+  to_b = 1;
+  mul2exp<To_Policy, To_Policy>(to_b, to_b, rational_sqrt_precision_parameter, ROUND_IGNORE);
   to.canonicalize();
   return rdiv != V_EQ ? rdiv : rsqrt;
 }
@@ -456,7 +462,8 @@ rational_sqrt_precision_parameter() {
 }
 
 //! Sets the precision parameter used for rational square root calculations.
-/*!
+/*! The lesser between numerator and denominator is limited to 2**\p p.
+
   If \p p is less than or equal to <CODE>INT_MAX</CODE>, sets the
   precision parameter used for rational square root calculations to \p p.