Implement Graeffe Transform for arb_poly and acb_poly

[MGessinger]

If P is a polynomial (over an arbitrary ring R), then its Graeffe transform G is defined through the equation G(x^2) = P(x)P(-x) (up to an overall sign, which is determined by the degree of P). If L//R is the splitting field of P, then it is easy to see, that G also decomposes over L, and its roots are exactly the squares of the roots of P. In particular, in some cases the roots of G might be real even though those of P were not. Due to the symmetry of G, we see that in computing P(x)P(-x), half the terms will vanish (namely precisely the odd terms). This wastes time and memory. Instead, this implementation uses the definition in [Diouf, Sec. 1.2], which yields a significant performance benefit for the reason mentioned above.

This polynomial transform finds application in Graeffe's Method [Diouf, Sec. 3] for bounding polynomial roots (or Fujiwara's bound, as is already implemented in acb_poly_root_bound_fujiwara). For the application of said method, the Graeffe Transform will generally be applied many times in a row. As such, it might be of interest to provide another method that takes (alongside the usual arguments) an additional iterations parameter, and performs the transform that many times, as opposed to calling the usual function over and over again. The benefit of having such an additional function is to avoid excessive internal memory (de-)allocations.

Tests have been provided for arb_poly_graeffe_transform and acb_poly_graeffe_transform, using the file t-acb_poly_borel_transform.c as style reference. Since the code is identical for the real and complex transform, the tests are also identical. Sadly, the transform doesn't decompose into real and imaginary parts, so code duplication cannot be avoided.

[fredrik-johansson]

Looks good, though the documentation is missing.

Also, please follow the code style of putting spaces around binary operators (x + y instead of x+y).

I think you can speed up the test code as well. There's no reason to test with up to n = 200 roots; n = 20 should be more than enough.

[MGessinger]

Thank you for your feedback. I have provided documentation in what I believe to be the correct places, as well as implementing the code changes you suggested. If you spot any further issues, I will happily see to them.

[MGessinger]

The stylistic mistakes were fixed, and the suggested optimization was implemented. Both arrays (for the even and odd terms) now have minimal length, and all polynomial operations use as few terms as possible. I also used the simplification of writing to b+1 directly. Note that while I could do this "trick" right when I perform the multiplication (line 42 in graeffe_transform.c currently), that would require me to set the lowest term to zero. Instead, I use this "trick" one step later, which avoids that inconvenience.

In optimizing the array lengths, one encounters a possible problem, if b already contains data. To solve this problem, the leading coefficient of b is set manually in the case that the degree is odd (i.e. the length is even). In the other case, this problem does not occur.

[fredrik-johansson]

Excellent, thanks!

[fredrik-johansson]

By the way, a possible issue: the test code only generates monic polynomials. You might want to test non-monic polynomials as well to make sure the leading coefficient is handled correctly.

[MGessinger]

Ah, good point. I altered the test code to include a randomly generated prefactor, and indeed it turns out that the leading coefficient was not handled correctly, albeit only in the case of constant input polynomials (which, given that Graeffe's method is used for finding roots, should hardly ever be of interest, but still needs to be covered). This has been fixed now. How would I go about having the fix merged? Do I open a new Pull Request? Apologies for the inconvenience.

[fredrik-johansson]

No problem. You can just make a new PR.