|
Sign manipulation:
diff --git a/doc/html/boost_multiprecision/ref/number.html b/doc/html/boost_multiprecision/ref/number.html
index 398ef05f..43451efa 100644
--- a/doc/html/boost_multiprecision/ref/number.html
+++ b/doc/html/boost_multiprecision/ref/number.html
@@ -176,8 +176,13 @@
struct is_number_expression;
+unmentionable-expression-template-type gcd(const number-or-expression-template-type&, const number-or-expression-template-type&);
+unmentionable-expression-template-type lcm(const number-or-expression-template-type&, const number-or-expression-template-type&);
unmentionable-expression-template-type pow(const number-or-expression-template-type&, unsigned);
unmentionable-expression-template-type powm(const number-or-expression-template-type& b, const number-or-expression-template-type& p, const number-or-expression-template-type& m);
+unmentionable-expression-template-type sqrt(const number-or-expression-template-type&);
+template <class Backend, expression_template_option ExpressionTemplates>
+number<Backend, EXpressionTemplates> sqrt(const number-or-expression-template-type&, number<Backend, EXpressionTemplates>&);
template <class Backend, expression_template_option ExpressionTemplates>
void divide_qr(const number-or-expression-template-type& x, const number-or-expression-template-type& y,
number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r);
@@ -1110,6 +1115,20 @@
when the algorithm requires it. Versions overloaded for built in integer
types return that integer type rather than an expression template.
+unmentionable-expression-template-type gcd(const number-or-expression-template-type& a, const number-or-expression-template-type& b);
+
+
+ Returns the largest integer x
+ that divides both a and
+ b.
+
+unmentionable-expression-template-type lcm(const number-or-expression-template-type& a, const number-or-expression-template-type& b);
+
+
+ Returns the smallest integer x
+ that is divisible by both a
+ and b.
+
unmentionable-expression-template-type pow(const number-or-expression-template-type& b, unsigned p);
@@ -1127,6 +1146,25 @@
Returns bp mod m as an expression template. Fixed precision
types are promoted internally to ensure accuracy.
+unmentionable-expression-template-type sqrt(const number-or-expression-template-type& a);
+
+
+ Returns the largest integer x
+ such that x *
+ x <
+ a.
+
+template <class Backend, expression_template_option ExpressionTemplates>
+number<Backend, EXpressionTemplates> sqrt(const number-or-expression-template-type& a, number<Backend, EXpressionTemplates>& r);
+
+
+ Returns the largest integer x
+ such that x *
+ x <
+ a, and sets the remainder r such that r
+ = a - x *
+ x.
+
template <class Backend, expression_template_option ExpressionTemplates>
void divide_qr(const number-or-expression-template-type& x, const number-or-expression-template-type& y,
number<Backend, ExpressionTemplates>& q, number<Backend, ExpressionTemplates>& r);
diff --git a/doc/html/boost_multiprecision/tut/gen_int.html b/doc/html/boost_multiprecision/tut/gen_int.html
index a49ddaba..b3065615 100644
--- a/doc/html/boost_multiprecision/tut/gen_int.html
+++ b/doc/html/boost_multiprecision/tut/gen_int.html
@@ -144,6 +144,16 @@
Flips the index bit in val.
+template <class Integer>
+Integer sqrt(const Integer& x);
+template <class Integer>
+Integer sqrt(const Integer& x, Integer& r);
+
+
+ Returns the integer square root s
+ of x and sets r to the remainder
+ x - s2.
+
template <class Engine>
bool miller_rabin_test(const number-or-expression-template-type& n, unsigned trials, Engine& gen);
bool miller_rabin_test(const number-or-expression-template-type& n, unsigned trials);
diff --git a/doc/html/index.html b/doc/html/index.html
index e30bc459..4ddc56d3 100644
--- a/doc/html/index.html
+++ b/doc/html/index.html
@@ -147,7 +147,7 @@
-Last revised: July 08, 2013 at 12:29:44 GMT |
+Last revised: August 01, 2013 at 15:58:28 GMT |
|
diff --git a/doc/multiprecision.qbk b/doc/multiprecision.qbk
index 71b47bad..00b2eac8 100644
--- a/doc/multiprecision.qbk
+++ b/doc/multiprecision.qbk
@@ -1828,6 +1828,13 @@ Unsets the `index` bit in `val`.
Flips the `index` bit in `val`.
+ template
+ Integer sqrt(const Integer& x);
+ template
+ Integer sqrt(const Integer& x, Integer& r);
+
+Returns the integer square root `s` of x and sets `r` to the remainder ['x - s[super 2]].
+
template
bool miller_rabin_test(const number-or-expression-template-type& n, unsigned trials, Engine& gen);
bool miller_rabin_test(const number-or-expression-template-type& n, unsigned trials);
@@ -1991,8 +1998,13 @@ generic operations, and so function equally well for built in and multiprecision
struct is_number_expression;
// Integer specific functions:
+ ``['unmentionable-expression-template-type]`` gcd(const ``['number-or-expression-template-type]``&, const ``['number-or-expression-template-type]``&);
+ ``['unmentionable-expression-template-type]`` lcm(const ``['number-or-expression-template-type]``&, const ``['number-or-expression-template-type]``&);
``['unmentionable-expression-template-type]`` pow(const ``['number-or-expression-template-type]``&, unsigned);
``['unmentionable-expression-template-type]`` powm(const ``['number-or-expression-template-type]``& b, const ``['number-or-expression-template-type]``& p, const ``['number-or-expression-template-type]``& m);
+ ``['unmentionable-expression-template-type]`` sqrt(const ``['number-or-expression-template-type]``&);
+ template
+ number sqrt(const ``['number-or-expression-template-type]``&, number&);
template
void divide_qr(const ``['number-or-expression-template-type]``& x, const ``['number-or-expression-template-type]``& y,
number& q, number& r);
@@ -2387,6 +2399,14 @@ built in integers or multiprecision ones), the functions will promote to a wider
requires it. Versions overloaded for built in integer types return that integer type rather than an expression
template.
+ ``['unmentionable-expression-template-type]`` gcd(const ``['number-or-expression-template-type]``& a, const ``['number-or-expression-template-type]``& b);
+
+Returns the largest integer `x` that divides both `a` and `b`.
+
+ ``['unmentionable-expression-template-type]`` lcm(const ``['number-or-expression-template-type]``& a, const ``['number-or-expression-template-type]``& b);
+
+Returns the smallest integer `x` that is divisible by both `a` and `b`.
+
``['unmentionable-expression-template-type]`` pow(const ``['number-or-expression-template-type]``& b, unsigned p);
Returns ['b[super p]] as an expression template. Note that this function should be used with extreme care as the result can grow so
@@ -2398,6 +2418,15 @@ fixed precision `cpp_int`'s either.
Returns ['b[super p] mod m] as an expression template. Fixed precision types are promoted internally to ensure accuracy.
+ ``['unmentionable-expression-template-type]`` sqrt(const ``['number-or-expression-template-type]``& a);
+
+Returns the largest integer `x` such that `x * x < a`.
+
+ template
+ number sqrt(const ``['number-or-expression-template-type]``& a, number& r);
+
+Returns the largest integer `x` such that `x * x < a`, and sets the remainder `r` such that `r = a - x * x`.
+
template
void divide_qr(const ``['number-or-expression-template-type]``& x, const ``['number-or-expression-template-type]``& y,
number& q, number& r);
@@ -3097,6 +3126,9 @@ and `ff` is a variable of type `std::ios_base::fmtflags`.
The type of `a` shall be listed in one of the type lists
`B::signed_types`, `B::unsigned_types`.
The default version of this function is synthesised from other operations above.][[space]]]
+[[`eval_integer_sqrt(b, cb, b2)`][`void`][Sets `b` to the largest integer which when squared is less than `cb`, also
+ sets `b2` to the remainder, ie to ['cb - b[super 2]].
+ The default version of this function is synthesised from other operations above.][[space]]]
[[['Sign manipulation:]]]
[[`eval_abs(b, cb)`][`void`][Set `b` to the absolute value of `cb`.
diff --git a/include/boost/multiprecision/detail/default_ops.hpp b/include/boost/multiprecision/detail/default_ops.hpp
index 79f521df..85c5a021 100644
--- a/include/boost/multiprecision/detail/default_ops.hpp
+++ b/include/boost/multiprecision/detail/default_ops.hpp
@@ -1116,6 +1116,59 @@ inline void eval_bit_unset(T& val, unsigned index)
eval_bitwise_xor(val, mask);
}
+template
+void eval_integer_sqrt(B& s, B& r, const B& x)
+{
+ //
+ // This is slow bit-by-bit integer square root, see for example
+ // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
+ // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
+ // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
+ // at some point.
+ //
+ typedef typename boost::multiprecision::detail::canonical::type ui_type;
+
+ s = ui_type(0u);
+ if(eval_get_sign(x) == 0)
+ {
+ r = ui_type(0u);
+ return;
+ }
+ int g = eval_msb(x);
+ if(g == 0)
+ {
+ r = ui_type(1);
+ return;
+ }
+
+ B t;
+ r = x;
+ g /= 2;
+ int org_g = g;
+ eval_bit_set(s, g);
+ eval_bit_set(t, 2 * g);
+ eval_subtract(r, x, t);
+ --g;
+ int msbr = eval_msb(r);
+ do
+ {
+ if(msbr >= org_g + g + 1)
+ {
+ t = s;
+ eval_left_shift(t, g + 1);
+ eval_bit_set(t, 2 * g);
+ if(t.compare(r) <= 0)
+ {
+ eval_bit_set(s, g);
+ eval_subtract(r, t);
+ msbr = eval_msb(r);
+ }
+ }
+ --g;
+ }
+ while(g >= 0);
+}
+
//
// These have to implemented by the backend, declared here so that our macro generated code compiles OK.
//
@@ -1465,6 +1518,26 @@ inline long long llround(const number& v)
}
#endif
+template
+inline typename enable_if_c::value == number_kind_integer, number >::type
+ sqrt(const number& x)
+{
+ using default_ops::eval_integer_sqrt;
+ number s, r;
+ eval_integer_sqrt(s.backend(), r.backend(), x.backend());
+ return s;
+}
+
+template
+inline typename enable_if_c::value == number_kind_integer, number >::type
+ sqrt(const number& x, number& r)
+{
+ using default_ops::eval_integer_sqrt;
+ number s;
+ eval_integer_sqrt(s.backend(), r.backend(), x.backend());
+ return s;
+}
+
#define UNARY_OP_FUNCTOR(func, category)\
namespace detail{\
template \
diff --git a/include/boost/multiprecision/gmp.hpp b/include/boost/multiprecision/gmp.hpp
index 1b86ba97..1e9c83ef 100644
--- a/include/boost/multiprecision/gmp.hpp
+++ b/include/boost/multiprecision/gmp.hpp
@@ -1592,6 +1592,11 @@ inline typename enable_if_c::value && ((sizeof(I) <= sizeof(long)))
mpz_lcm_ui(result.data(), a.data(), std::abs(b));
}
+inline void eval_integer_sqrt(gmp_int& s, gmp_int& r, const gmp_int& x)
+{
+ mpz_sqrtrem(s.data(), r.data(), x.data());
+}
+
inline unsigned eval_lsb(const gmp_int& val)
{
int c = eval_get_sign(val);
diff --git a/include/boost/multiprecision/integer.hpp b/include/boost/multiprecision/integer.hpp
index b3682006..cdac282e 100644
--- a/include/boost/multiprecision/integer.hpp
+++ b/include/boost/multiprecision/integer.hpp
@@ -182,6 +182,59 @@ typename enable_if_c::value, Integer&>::type bit_flip(Integ
return val;
}
+template
+typename enable_if_c::value, Integer>::type sqrt(const Integer& x, Integer& r)
+{
+ //
+ // This is slow bit-by-bit integer square root, see for example
+ // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_.28base_2.29
+ // There are better methods such as http://hal.inria.fr/docs/00/07/28/54/PDF/RR-3805.pdf
+ // and http://hal.inria.fr/docs/00/07/21/13/PDF/RR-4475.pdf which should be implemented
+ // at some point.
+ //
+ Integer s = 0;
+ if(x == 0)
+ {
+ r = 0;
+ return s;
+ }
+ int g = msb(x);
+ if(g == 0)
+ {
+ r = 1;
+ return s;
+ }
+
+ Integer t;
+ r = x;
+ g /= 2;
+ bit_set(s, g);
+ bit_set(t, 2 * g);
+ r = x - t;
+ --g;
+ do
+ {
+ t = s;
+ t <<= g + 1;
+ bit_set(t, 2 * g);
+ if(t <= r)
+ {
+ bit_set(s, g);
+ r -= t;
+ }
+ --g;
+ }
+ while(g >= 0);
+ return s;
+}
+
+template
+typename enable_if_c::value, Integer>::type sqrt(const Integer& x)
+{
+ Integer r;
+ return sqrt(x, r);
+}
+
}} // namespaces
#endif
diff --git a/test/test_cpp_int.cpp b/test/test_cpp_int.cpp
index fc6e35d5..70c5f7b4 100644
--- a/test/test_cpp_int.cpp
+++ b/test/test_cpp_int.cpp
@@ -190,6 +190,11 @@ struct tester
BOOST_CHECK_EQUAL(mpz_int(lcm(-c, -d)).str(), test_type(lcm(-c1, -d1)).str());
BOOST_CHECK_EQUAL(mpz_int(gcd(a, -b)).str(), test_type(gcd(a1, -b1)).str());
BOOST_CHECK_EQUAL(mpz_int(lcm(c, -d)).str(), test_type(lcm(c1, -d1)).str());
+ // Integer sqrt:
+ mpz_int r;
+ test_type r1;
+ BOOST_CHECK_EQUAL(sqrt(a, r).str(), sqrt(a1, r1).str());
+ BOOST_CHECK_EQUAL(r.str(), r1.str());
}
void t3()
diff --git a/test/test_native_integer.cpp b/test/test_native_integer.cpp
index face73ca..24ae85a5 100644
--- a/test/test_native_integer.cpp
+++ b/test/test_native_integer.cpp
@@ -68,6 +68,19 @@ void test()
BOOST_CHECK_EQUAL(integer_modulus(i, j), i % j);
I p = 456;
BOOST_CHECK_EQUAL(powm(i, p, j), pow(cpp_int(i), static_cast(p)) % j);
+
+ for(I i = 0; i < (2 < 8) - 1; ++i)
+ {
+ I j = i * i;
+ I s, r;
+ s = sqrt(j, r);
+ BOOST_CHECK_EQUAL(s, i);
+ BOOST_CHECK(r == 0);
+ j += 3;
+ s = sqrt(i, r);
+ BOOST_CHECK_EQUAL(s, i);
+ BOOST_CHECK(r == 3);
+ }
}
int main()
|