Polynomial fixes: for mixed-type arithmetic etc.

[kundor]

I was fixing some minor problems in the polynomial documentation, when I noticed that operator / (polynomial<T>, U) and operator % (polynomial<T>, U) were not actually templated on the type of the scalar U, as reported in the docs, but were provided by Boost.Operators and only accepted scalars of type T.

This isn't very noticeable, since either way, the code compiles only if U is implicitly convertible to T. However, it does cause different behavior. For instance, polynomial<int>{2,4,6} / 1.5 is an identity operation (the 1.5 is truncated to 1), while polynomial<int>{2,4,6} * 0.6667 yields {1,2,4}.
More seriously, if a is a polynomial<int>, a = a / 1.5 and a /= 1.5 would give different results.

To fix this, I added templated implementations of these operators just like all the other polynomial-scalar operators.

For consistency, I used Boost.Operators to provide all the polynomial-polynomial operations (it was being used just for division, modulus, and operator !=), which are conveniently packed together in the class ordered_euclidean_ring_operators< polynomial<T> >. (This entailed also adding an operator <, which also enables polynomials to be put in std::set etc. It compares by degree first, so that it's compatible with the Euclidean ordering, then lexicographically starting with the highest-degree coefficient.)

The self-assign operators are all templated on polynomial<U>. However, operator /= and %= call out to quotient_remainder, which only takes two matching polynomial<T> arguments. So it's ambiguous what to convert when calling quotient_remainder(polynomial<T>, polynomial<U>), causing an error.

To fix this, polynomial<T>::operator /=(const polynomial<U>& q) can call quotient_remainder(*this, polynomial<T>(q)). But then you might as well not template the operator on type U, and just make it polynomial::operator /=(const polynomial& q); you can still use it with any polynomial<U> because there is an implicit conversion to polynomial<T>.

Unfortunately, just removing the templating causes polynomial<U> arguments to choose the scalar version operator /=(const U& value) instead. So I had to guard that function with enable_if< boost::is_convertible<U,T> >.

This situation does result in inconsistent behavior. For instance, multiplying a polynomial<int> by a polynomial<double> will do each multiplication with doubles before truncating back to int, whereas dividing by a polynomial<double> will first truncate all the values to int.
Similarly, with polynomial<int> a{3}, a /= 1.5 results in 2, but a /= polynomial<double>{1.5} results in 3.

Changing this, so that it promotes internal computations and then truncates, would mean changing the implementation of quotient_remainder.

Finally, when T is integral, there can be a remainder after dividing by a scalar, so I added an implementation of operator %= for this (instead of just setting to zero as before.) It uses enable_if_c<std::numeric_limits<T>::is_integer>—in theory, the algorithm used should work for other types, but with floating point types it would probably end up with annoying near-zero coefficients.

I would be happy to fix up this PR if you want to merge #39 first (it will conflict a little.)

[jeremy-murphy]

Thanks for picking up these defects. I remember now when I was implementing operator >>, etc that it did not work using the Boost.Operator facility because I needed a new template type for the right-hand parameter.

[jeremy-murphy]

It would be great if you could include in this PR the unit tests that would otherwise fail without these changes , including explicit tests of the new modulus function.

[jeremy-murphy]

Hmmm, I guess this begs the question as to whether non-PoD numbers (such as polynomial) should promote according to the same rules as PoDs. What do the multiprecision classes do, John? How do they handle arithmetic with a mixture of types?

[kundor]

Speaking of which, is there really a need to template the shift operators on the right-hand parameter? Is there a case for calling it with anything but an int (or something convertible thereto)? I don't think this class is meant to support polynomials of degree 32,000 or more ;-) (a signed int goes up at least to 2^15-1, by the standard.) (Taking an int and checking that it's positive is preferable to taking an unsigned number, in that otherwise a call p <<= -1 will compile without complaint, and silently try to extend p by 4 billion elements.)

[jzmaddock]

On 15/05/2016 16:46, Nick Matteo wrote:

Speaking of which, is there really a need to template the shift
operators on the right-hand parameter? Is there a case for calling it
with anything but an |int| (or something convertible thereto)? I don't
think this class is meant to support polynomials of degree 32,000 or
more ;-) (a signed int goes up at least to 2^15-1, by the standard.)
Taking an |int| and checking that it's positive is preferable to
taking an unsigned number, in that otherwise a call |p <<= -1| will
compile without complaint, and silently try to extend |p| by 4 million
elements.

I wouldn't think so no. Jeremy?

[jzmaddock]

On 15/05/2016 12:53, Jeremy W. Murphy wrote:

Hmmm, I guess this begs the question as to whether non-PoD numbers
(such as polynomial) should promote according to the same rules as
PoDs. What do the multiprecision classes do, John? How do they handle
arithmetic with a mixture of types?

Well first off they don't use Boost.Operators which is rather ill suited
to this kind of thing, but basically given:

a op b

If there is an implicit (ie non-lossy) conversion from b to a, then b is
converted to typeof(a), otherwise if a is implicitly convertible to
typeof(b) then a gets promoted.

Of course internally, the value may be used directly where it is more
efficient to do so, but the operators themselves enforce the above rule.

In this situation, for:

polynomial op U

we would want to enforce that there is an implicit, non-lossy conversion
from U to T, which means no multiplication of polynomial by a
double etc.

The logic used is somewhat multiprecision-specific, and in
boost/multiprecision/traits/is_restricted_conversion.hpp

HTH, John.

[kundor]

I rebased the branch, and added the tests Jeremy asked for.

Is it necessary for tests to be pre-C++11? It's kind of painful.

[jeremy-murphy]

With respect to the shift operators, no, there is no particular need for them to be templated, and off the cuff it only makes sense to call them with an integral type. How to handle negative values raises a curious thought, though: should operator>>(-x) be equivalent to operator<<(x)? I understand that for integers, shifting by a negative amount is undefined behaviour, so probably we should follow suit, but we're not obliged to.

[jeremy-murphy]

Hmmm, so for operator*=(U) we should be restricting U to be implicitly convertible to T?

[kundor]

Re: operator*=(U). We can restrict all the scalar operations to arguments that are implicitly convertible to T, by just removing the templating and only taking T. The question is whether we want

polynomial<int> p{3,2};
p *= 1.5;

to give us 3x + 4, analogous to what int i=3; i *= 1.5; would do, or be an identity operation, or be a compile-time error.

Personally, I would rather have it act like built-in ints, i.e. multiply by the double and then truncate, which is what it does now.

Making it an error would require a templated version of the function just to fail, using something like the is_lossy_conversion trait from multiprecision, which seems like overkill.

On the other hand, removing the templating from poly-poly operations (operator *=(polynomial<U>) etc.), would make their behavior more consistent (in that currently, operator/=(polynomial) must act that way.) I think currently I am in favor of removing arithmetic between polynomial and polynomial entirely, which means removing the template from each operator op=(polynomial) function, and making the conversion constructor polynomial(const polynomial<U>& p) explicit.

[kundor]

Commit 910a62a removes arithmetic between polynomials with different base types, so you'll have to explicitly cast one operand or the other (e.g. if p is a polynomial<int> and q is a polynomial<double>, you can write polynomial<double>(p) * q, or p *= polynomial<int>(q).

All the tests in test_polynomial.cpp still pass.

[jeremy-murphy]

Should modulus still have a U template?

[kundor]

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.