Implement Graeffe Transform for real and complex poly

acb_poly.h

@@ -193,6 +193,10 @@ void _acb_poly_binomial_transform(acb_ptr b, acb_srcptr a, slong alen, slong len
 
 void acb_poly_binomial_transform(acb_poly_t b, const acb_poly_t a, slong len, slong prec);
 
+void _acb_poly_graeffe_transform(acb_ptr b, acb_srcptr a, slong len, slong prec);
+
+void acb_poly_graeffe_transform(acb_poly_t b, const acb_poly_t a, slong prec);
+
 
 
 void acb_poly_set(acb_poly_t dest, const acb_poly_t src);

acb_poly/graeffe_transform.c

@@ -0,0 +1,52 @@
+/*
+    Copyright (C) 2012 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_poly.h"
+
+void
+_acb_poly_graeffe_transform(acb_ptr b, acb_srcptr a, slong len, slong prec)
+{
+    if (len == 0)
+        return;
+
+    slong q, i;
+
+    q = (len-1)/2+1;
+    acb_ptr pe = _acb_vec_init(q);
+    acb_ptr po = _acb_vec_init(len);
+
+    for (i = len-1; i >= 0; i--)
+    {
+        if (i % 2 == 0)
+            acb_set(pe+i/2, a+i);
+        else
+            acb_set(po+i/2, a+i);
+    }
+
+    _acb_poly_mul(b, po, q, po, q, prec);
+    _acb_poly_shift_left(b, b, len-1, 1);
+    _acb_poly_mul(po, pe, q, pe, q, prec);
+    _acb_vec_sub(b, po, b, len, prec);
+
+    _acb_vec_clear(pe, q);
+    _acb_vec_clear(po, len);
+
+    if (len % 2 == 0)
+	    _acb_vec_neg(b, b, len);
+}
+
+void
+acb_poly_graeffe_transform(acb_poly_t b, const acb_poly_t a, slong prec)
+{
+    acb_poly_fit_length(b, a->length);
+    _acb_poly_graeffe_transform(b->coeffs, a->coeffs, a->length, prec);
+    _acb_poly_set_length(b, a->length);
+}

arb_poly.h

@@ -512,6 +512,10 @@ void _arb_poly_binomial_transform(arb_ptr b, arb_srcptr a, slong alen, slong len
 
 void arb_poly_binomial_transform(arb_poly_t b, const arb_poly_t a, slong len, slong prec);
 
+void _arb_poly_graeffe_transform(arb_ptr b, arb_srcptr a, slong len, slong prec);
+
+void arb_poly_graeffe_transform(arb_poly_t b, const arb_poly_t a, slong prec);
+
 /* Special functions */
 
 void _arb_poly_pow_ui_trunc_binexp(arb_ptr res,

arb_poly/graeffe_transform.c

@@ -0,0 +1,52 @@
+/*
+    Copyright (C) 2012 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_poly.h"
+
+void
+_arb_poly_graeffe_transform(arb_ptr b, arb_srcptr a, slong len, slong prec)
+{
+    if (len == 0)
+        return;
+
+    slong q, i;
+
+    q = (len-1)/2+1;
+    arb_ptr pe = _arb_vec_init(q);
+    arb_ptr po = _arb_vec_init(len);
+
+    for (i = len-1; i >= 0; i--)
+    {
+        if (i % 2 == 0)
+            arb_set(pe+i/2, a+i);
+        else
+            arb_set(po+i/2, a+i);
+    }
+
+    _arb_poly_mul(b, po, q, po, q, prec);
+    _arb_poly_shift_left(b, b, len-1, 1);
+    _arb_poly_mul(po, pe, q, pe, q, prec);
+    _arb_vec_sub(b, po, b, len, prec);
+
+    _arb_vec_clear(pe, q);
+    _arb_vec_clear(po, len);
+
+    if (len % 2 == 0)
+	    _arb_vec_neg(b, b, len);
+}
+
+void
+arb_poly_graeffe_transform(arb_poly_t b, const arb_poly_t a, slong prec)
+{
+    arb_poly_fit_length(b, a->length);
+    _arb_poly_graeffe_transform(b->coeffs, a->coeffs, a->length, prec);
+    _arb_poly_set_length(b, a->length);
+}