Add algorithm for difference of squares and move is_monotonic_increasing() to utils

[amandalund]

This pulls out a couple small changes from #1080.

[sethrj]

Thanks!

[sethrj]

By the way, I originally thought the expression they used was to eliminate roundoff error when the two numbers are different magnitudes, but inspired by @paulromano 's use of mpmath it's clear that's not the case (shown here for single-precision arithmetic):

https://gist.github.com/sethrj/e0b1a0e92d4b89e11e14ffed2b3fee01

So I guess those routines were written that way as a microoptimization to turn one multiply * into an addition +?

Yet it looks like it takes an extra register and an extra operation to do! (See godbolt for both arm and x86.)

Maybe we should follow up by writing the expression as God intended it...

@pcanal Is there some edge case I'm missing where the (a - b) * (a + b) formulation might have some advantage?

[amandalund]

I thought your original guess was right, but intended to reduce the error when they are close in magnitude. However, in practice I didn't really see it make a difference either. I tried to come up with some examples when I wrote the test (and also played around with mpmath), but didn't have much luck finding a case where rearranging the expression helped. So I'd also support just writing it the usual way.

[sethrj]

Oh right. Well once a and b get stored as floating point values, there seems to be still no advantage:

[pcanal]

@pcanal Is there some edge case I'm missing where the (a - b) * (a + b) formulation might have some advantage?

In the cases where a and b are far apart (order of magnitude) but not far enough (i.e. a-b != a) you might have (after numerical errors):

(a*a) - (b*b) == (a*a)
but
(a-b) != a
(a+b) != a
and thus
(a-b) * (a+b) != a*a

[sethrj]

If they're far apart, the squares are further... so the result doesn't change 😕