Skip to content

Polynomial fixes: for mixed-type arithmetic etc. #41

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
kundor wants to merge 16 commits into
base: develop
Choose a base branch
from
+161 −96
Open

Polynomial fixes: for mixed-type arithmetic etc. #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
a1b34b8
Some polynomial doc fixes
kundor May 13, 2016
bf46cd3
Fix polynomial<T> /= polynomial<U> (call to quotient_remainder was am…
kundor May 13, 2016
1271818
Use Boost.operators for all poly-poly operations, and no poly-scalar …
kundor May 14, 2016
9c7d39e
Consistently use `inline`
kundor May 14, 2016
78d7dc0
Prevent binding polynomial<U>s to the scalar operators.
kundor May 14, 2016
fb6b899
For integer types, give a remainder with scalar values
kundor May 15, 2016
89331d0
Update polynomial docs
kundor May 15, 2016
a1f74c1
Taking a const reference just to copy it is silly
kundor May 15, 2016
f4d73fe
Tests for mixed-type arithmetic, poly<int> % scalar, ordering
kundor May 15, 2016
35022a6
Do not template shift operators' right argument
kundor May 19, 2016
8646ae4
Silence comparison warnings
kundor May 19, 2016
910a62a
Remove polynomial<T> x polynomial<U> arithmetic
kundor May 19, 2016
2798e23
No more need for enable_if guards
kundor May 22, 2016
885e13f
Update polynomial docs
kundor May 22, 2016
1115101
Unnecessary include
kundor May 23, 2016
9da5ff2
Merge branch 'develop' into poly-fixes
kundor May 28, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 22 additions & 18 deletions doc/internals/polynomial.qbk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@
polynomial(I first, I last);
template <class U>
explicit polynomial(const U& point);
template <class U>
explicit polynomial(const polynomial<U>& p);
polynomial(std::initializer_list<T> l); // C++11

// access:
Expand All @@ -50,16 +52,13 @@
polynomial& operator /=(const U& value);
template <class U>
polynomial& operator %=(const U& value);
template <class U>
polynomial& operator +=(const polynomial<U>& value);
template <class U>
polynomial& operator -=(const polynomial<U>& value);
template <class U>
polynomial& operator *=(const polynomial<U>& value);
template <class U>
polynomial& operator /=(const polynomial<U>& value);
template <class U>
polynomial& operator %=(const polynomial<U>& value);
polynomial& operator >>=(int n);
polynomial& operator <<=(int n);
polynomial& operator +=(const polynomial& value);
polynomial& operator -=(const polynomial& value);
polynomial& operator *=(const polynomial& value);
polynomial& operator /=(const polynomial& value);
polynomial& operator %=(const polynomial& value);
};

template <class T>
Expand Down Expand Up @@ -92,29 +91,34 @@
polynomial<T> operator * (const U& a, const polynomial<T>& b);

template <class T>
polynomial<T> operator - (const polynomial<T>& a);
polynomial<T> operator - (const polynomial<T>& a); // Unary minus

template <class T>
polynomial<T> operator >>= (const U& a);
template <class T>
polynomial<T> operator >> (polynomial<T> const &a, const U& b);
template <class T>
polynomial<T> operator <<= (const U& a);
polynomial<T> operator >> (polynomial<T> const &a, int b);
template <class T>
polynomial<T> operator << (polynomial<T> const &a, const U& b);
polynomial<T> operator << (polynomial<T> const &a, int b);

template <class T>
bool operator == (const polynomial<T> &a, const polynomial<T> &b);
template <class T>
bool operator != (const polynomial<T> &a, const polynomial<T> &b);

template <class T>
polynomial<T> pow(polynomial<T> base, int exp);
bool operator < (const polynomial<T> &a, const polynomial<T> &b);
template <class T>
bool operator <= (const polynomial<T> &a, const polynomial<T> &b);
template <class T>
bool operator > (const polynomial<T> &a, const polynomial<T> &b);
template <class T>
bool operator >= (const polynomial<T> &a, const polynomial<T> &b);

template <class charT, class traits, class T>
std::basic_ostream<charT, traits>& operator <<
(std::basic_ostream<charT, traits>& os, const polynomial<T>& poly);

template <class T>
polynomial<T> pow(polynomial<T> base, int exp);

template <typename T>
std::pair< polynomial<T>, polynomial<T> >
quotient_remainder(const polynomial<T>& a, const polynomial<T>& b);
Expand Down
160 changes: 82 additions & 78 deletions include/boost/math/tools/polynomial.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -255,12 +255,11 @@ quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)


template <class T>
class polynomial :
equality_comparable< polynomial<T>,
dividable< polynomial<T>,
dividable2< polynomial<T>, T,
modable< polynomial<T>,
modable2< polynomial<T>, T > > > > >
class polynomial : ordered_euclidean_ring_operators< polynomial<T> >
/* Provides operators +, -, *, /, % for two polynomials (using member += etc.);
* operators >, <=, >= using operator < ; and operator != using operator ==.
* Note that polynomials are not actually a euclidean ring when T is not
* a field type, in particular, if T is integral. */
{
public:
// typedefs:
Expand Down Expand Up @@ -296,7 +295,7 @@ class polynomial :
: m_data(p.m_data) { }

template <class U>
polynomial(const polynomial<U>& p)
explicit polynomial(const polynomial<U>& p)
{
for(unsigned i = 0; i < p.size(); ++i)
{
Expand Down Expand Up @@ -388,30 +387,37 @@ class polynomial :
}

template <class U>
polynomial& operator %=(const U& /*value*/)
{
// We can always divide by a scalar, so there is no remainder:
this->set_zero();
polynomial& operator %=(const U& value)
{
// In the case that T is integral, this preserves the semantics
// p == r*(p/r) + (p % r), for polynomial<T> p and U r.
if (std::numeric_limits<T>::is_integer)
Copy link
Copy Markdown
Contributor

@jeremy-murphy jeremy-murphy May 22, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since this is a compile-time constant, can't we combine it with enable_if for static dispatch?

Copy link
Copy Markdown
Contributor Author

@kundor kundor May 23, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I seem to recall that using enable_if based only on T, which is not one of the template parameters for the function, doesn't work. Let me check.

Copy link
Copy Markdown
Contributor Author

@kundor kundor May 23, 2016 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I was right: GCC 6.1, GCC 5.3, and Clang 703.0.31 all complain vociferously. To separate them, it's necessary to add a dummy condition to the enable_if that depends on U, e.g. go back to BOOST_DEDUCED_TYPENAME enable_if_c<boost::is_convertible<U,T>::value && std::numeric_limits<T>::is_integer, polynomial >::type& operator %=(const U& value).

I don't think this mess is worth it:

  • the modulus operator on a polynomial doesn't seem performance-critical
  • a conditional testing a compile-time boolean should only emit the chosen branch anyway

You can observe that the last thing is true at godbolt.

{
modulus(value);
normalize();
}
else
{
set_zero();
}
return *this;
}

template <class U>
polynomial& operator +=(const polynomial<U>& value)
polynomial& operator +=(const polynomial& value)
{
addition(value);
normalize();
return *this;
}

template <class U>
polynomial& operator -=(const polynomial<U>& value)
polynomial& operator -=(const polynomial& value)
{
subtraction(value);
normalize();
return *this;
}
template <class U>
polynomial& operator *=(const polynomial<U>& value)

polynomial& operator *=(const polynomial& value)
{
// TODO: FIXME: use O(N log(N)) algorithm!!!
if (!value)
Expand All @@ -427,31 +433,28 @@ class polynomial :
return *this;
}

template <typename U>
polynomial& operator /=(const polynomial<U>& value)
polynomial& operator /=(const polynomial& value)
{
*this = quotient_remainder(*this, value).first;
return *this;
}

template <typename U>
polynomial& operator %=(const polynomial<U>& value)
polynomial& operator %=(const polynomial& value)
{
*this = quotient_remainder(*this, value).second;
return *this;
}

template <typename U>
polynomial& operator >>=(U const &n)
polynomial& operator >>=(int n)
{
BOOST_ASSERT(n <= m_data.size());
BOOST_ASSERT(n >= 0 && static_cast<unsigned>(n) <= m_data.size());
m_data.erase(m_data.begin(), m_data.begin() + n);
return *this;
}

template <typename U>
polynomial& operator <<=(U const &n)
polynomial& operator <<=(int n)
{
BOOST_ASSERT(n >= 0);
m_data.insert(m_data.begin(), n, static_cast<T>(0));
normalize();
return *this;
Expand Down Expand Up @@ -554,64 +557,59 @@ class polynomial :
return *this;
}

template <class U>
polynomial& modulus(const U& value)
{
typedef typename std::vector<T>::iterator iter;
for (iter it = m_data.begin(); it != m_data.end(); ++it)
*it -= T(value * T(*it / value));
return *this;
}
Copy link
Copy Markdown
Contributor

@jeremy-murphy jeremy-murphy May 22, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should modulus still have a U template?

Copy link
Copy Markdown
Contributor Author

@kundor kundor May 22, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, operator %=(U), and the other self-modifying operators for scalar arguments, are all still templated. We can change that if you want to be more strict.


std::vector<T> m_data;
};


template <class T>
inline polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b)
{
polynomial<T> result(a);
result += b;
return result;
}

template <class T>
inline polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b)
template <class T, class U>
inline polynomial<T> operator + (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result -= b;
return result;
a += b;
return a;
}

template <class T>
inline polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b)
template <class T, class U>
inline polynomial<T> operator - (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result *= b;
return result;
a -= b;
return a;
}

template <class T, class U>
inline polynomial<T> operator + (const polynomial<T>& a, const U& b)
inline polynomial<T> operator * (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result += b;
return result;
a *= b;
return a;
}

template <class T, class U>
inline polynomial<T> operator - (const polynomial<T>& a, const U& b)
inline polynomial<T> operator / (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result -= b;
return result;
a /= b;
return a;
}

template <class T, class U>
inline polynomial<T> operator * (const polynomial<T>& a, const U& b)
inline polynomial<T> operator % (polynomial<T> a, const U& b)
{
polynomial<T> result(a);
result *= b;
return result;
a %= b;
return a;
}

template <class U, class T>
inline polynomial<T> operator + (const U& a, const polynomial<T>& b)
inline polynomial<T> operator + (const U& a, polynomial<T> b)
{
polynomial<T> result(b);
result += a;
return result;
b += a;
return b;
}

template <class U, class T>
Expand All @@ -623,51 +621,57 @@ inline polynomial<T> operator - (const U& a, const polynomial<T>& b)
}

template <class U, class T>
inline polynomial<T> operator * (const U& a, const polynomial<T>& b)
inline polynomial<T> operator * (const U& a, polynomial<T> b)
{
polynomial<T> result(b);
result *= a;
return result;
b *= a;
return b;
}

template <class T>
bool operator == (const polynomial<T> &a, const polynomial<T> &b)
inline bool operator == (const polynomial<T> &a, const polynomial<T> &b)
{
return a.data() == b.data();
}

template <typename T, typename U>
polynomial<T> operator >> (const polynomial<T>& a, const U& b)
template <class T>
inline bool operator < (const polynomial<T> &a, const polynomial<T> &b)
{
polynomial<T> result(a);
result >>= b;
return result;
if (a.size() != b.size())
return a.size() < b.size();
return std::lexicographical_compare(a.data().rbegin(), a.data().rend(),
b.data().rbegin(), b.data().rend());
}

template <typename T, typename U>
polynomial<T> operator << (const polynomial<T>& a, const U& b)
template <typename T>
inline polynomial<T> operator >> (polynomial<T> a, int b)
{
polynomial<T> result(a);
result <<= b;
return result;
a >>= b;
return a;
}

template <typename T>
inline polynomial<T> operator << (polynomial<T> a, int b)
{
a <<= b;
return a;
}

// Unary minus (negate).
template <class T>
polynomial<T> operator - (polynomial<T> a)
inline polynomial<T> operator - (polynomial<T> a)
{
std::transform(a.data().begin(), a.data().end(), a.data().begin(), std::negate<T>());
return a;
}

template <class T>
bool odd(polynomial<T> const &a)
inline bool odd(polynomial<T> const &a)
{
return a.size() > 0 && a[0] != static_cast<T>(0);
}

template <class T>
bool even(polynomial<T> const &a)
inline bool even(polynomial<T> const &a)
{
return !odd(a);
}
Expand Down
Loading