add divrem for complex polynomials

fmpcb_poly.h

@@ -162,6 +162,23 @@ void _fmpcb_poly_inv_series(fmpcb_struct * Qinv, const fmpcb_struct * Q, long le
 
 void fmpcb_poly_inv_series(fmpcb_poly_t Qinv, const fmpcb_poly_t Q, long n, long prec);
 
+void _fmpcb_poly_reverse(fmpcb_struct * res, const fmpcb_struct * poly, long len, long n);
+
+void _fmpcb_poly_div(fmpcb_struct * Q,
+    const fmpcb_struct * A, long lenA,
+    const fmpcb_struct * B, long lenB, long prec);
+
+void _fmpcb_poly_divrem(fmpcb_struct * Q, fmpcb_struct * R,
+    const fmpcb_struct * A, long lenA,
+    const fmpcb_struct * B, long lenB, long prec);
+
+void _fmpcb_poly_rem(fmpcb_struct * R,
+    const fmpcb_struct * A, long lenA,
+    const fmpcb_struct * B, long lenB, long prec);
+
+void fmpcb_poly_divrem(fmpcb_poly_t Q, fmpcb_poly_t R,
+                             const fmpcb_poly_t A, const fmpcb_poly_t B, long prec);
+
 void _fmpcb_poly_root_inclusion(fmpcb_t r, const fmpcb_t m,
     const fmpcb_struct * poly,
     const fmpcb_struct * polyder, long len, long prec);

fmpcb_poly/divrem.c

@@ -0,0 +1,146 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2011 William Hart
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+void
+_fmpcb_poly_div_series(fmpcb_struct * Q,
+    const fmpcb_struct * A,
+    const fmpcb_struct * B, long n, long prec)
+{
+    fmpcb_struct * Binv = _fmpcb_vec_init(n);
+
+    _fmpcb_poly_inv_series(Binv, B, n, prec);
+    _fmpcb_poly_mullow(Q, Binv, n, A, n, n, prec);
+
+    _fmpcb_vec_clear(Binv, n);
+}
+
+void
+_fmpcb_poly_div(fmpcb_struct * Q,
+    const fmpcb_struct * A, long lenA,
+    const fmpcb_struct * B, long lenB, long prec)
+{
+    const long lenQ = lenA - lenB + 1;
+    fmpcb_struct * Arev, * Brev;
+
+    Arev = _fmpcb_vec_init(2 * lenQ);
+    Brev = Arev + lenQ;
+
+    _fmpcb_poly_reverse(Arev, A + (lenA - lenQ), lenQ, lenQ);
+
+    if (lenB >= lenQ)
+    {
+        _fmpcb_poly_reverse(Brev, B + (lenB - lenQ), lenQ, lenQ);
+    }
+    else
+    {
+        _fmpcb_poly_reverse(Brev, B, lenB, lenB);
+        _fmpcb_vec_zero(Brev + lenB, lenQ - lenB);
+    }
+
+    _fmpcb_poly_div_series(Q, Arev, Brev, lenQ, prec);
+    _fmpcb_poly_reverse(Q, Q, lenQ, lenQ);
+
+    _fmpcb_vec_clear(Arev, 2 * lenQ);
+}
+
+void _fmpcb_poly_divrem(fmpcb_struct * Q, fmpcb_struct * R,
+    const fmpcb_struct * A, long lenA,
+    const fmpcb_struct * B, long lenB, long prec)
+{
+    const long lenQ = lenA - lenB + 1;
+    _fmpcb_poly_div(Q, A, lenA, B, lenB, prec);
+
+    if (lenB > 1)
+    {
+        if (lenQ >= lenB - 1)
+            _fmpcb_poly_mullow(R, Q, lenQ, B, lenB - 1, lenB - 1, prec);
+        else
+            _fmpcb_poly_mullow(R, B, lenB - 1, Q, lenQ, lenB - 1, prec);
+        _fmpcb_vec_sub(R, A, R, lenB - 1, prec);
+    }
+}
+
+void _fmpcb_poly_rem(fmpcb_struct * R,
+    const fmpcb_struct * A, long lenA,
+    const fmpcb_struct * B, long lenB, long prec)
+{
+    const long lenQ = lenA - lenB + 1;
+    fmpcb_struct * Q = _fmpcb_vec_init(lenQ);
+    _fmpcb_poly_divrem(Q, R, A, lenA, B, lenB, prec);
+    _fmpcb_vec_clear(Q, lenQ);
+}
+
+void fmpcb_poly_divrem(fmpcb_poly_t Q, fmpcb_poly_t R,
+                             const fmpcb_poly_t A, const fmpcb_poly_t B, long prec)
+{
+    const long lenA = A->length, lenB = B->length;
+
+    if (lenB == 0 || fmpcb_contains_zero(B->coeffs + lenB - 1))
+    {
+        printf("Exception: division by zero in fmpcb_poly_divrem\n");
+        abort();
+    }
+
+    if (lenA < lenB)
+    {
+        fmpcb_poly_set(R, A);
+        fmpcb_poly_zero(Q);
+        return;
+    }
+
+    if (Q == A || Q == B)
+    {
+        fmpcb_poly_t T;
+        fmpcb_poly_init(T);
+        fmpcb_poly_divrem(T, R, A, B, prec);
+        fmpcb_poly_swap(Q, T);
+        fmpcb_poly_clear(T);
+        return;
+    }
+
+    if (R == A || R == B)
+    {
+        fmpcb_poly_t U;
+        fmpcb_poly_init(U);
+        fmpcb_poly_divrem(Q, U, A, B, prec);
+        fmpcb_poly_swap(R, U);
+        fmpcb_poly_clear(U);
+        return;
+    }
+
+    fmpcb_poly_fit_length(Q, lenA - lenB + 1);
+    fmpcb_poly_fit_length(R, lenB - 1);
+
+    _fmpcb_poly_divrem(Q->coeffs, R->coeffs, A->coeffs, lenA,
+                                   B->coeffs, lenB, prec);
+
+    _fmpcb_poly_set_length(Q, lenA - lenB + 1);
+    _fmpcb_poly_set_length(R, lenB - 1);
+    _fmpcb_poly_normalise(R);
+}
+

fmpcb_poly/reverse.c

@@ -0,0 +1,55 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+void
+_fmpcb_poly_reverse(fmpcb_struct * res, const fmpcb_struct * poly, long len, long n)
+{
+    if (res == poly)
+    {
+        long i;
+
+        for (i = 0; i < n / 2; i++)
+        {
+            fmpcb_struct t = res[i];
+            res[i] = res[n - 1 - i];
+            res[n - 1 - i] = t;
+        }
+
+        for (i = 0; i < n - len; i++)
+            fmpcb_zero(res + i);
+    }
+    else
+    {
+        long i;
+
+        for (i = 0; i < n - len; i++)
+            fmpcb_zero(res + i);
+
+        for (i = 0; i < len; i++)
+            fmpcb_set(res + (n - len) + i, poly + (len - 1) - i);
+    }
+}

fmpcb_poly/test/t-divrem.c

@@ -0,0 +1,136 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb_poly.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("divrem....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        long m, n, qbits1, qbits2, rbits1, rbits2, rbits3;
+        fmpq_poly_t A, B, Q, R;
+        fmpcb_poly_t a, b, q, r;
+
+        qbits1 = 2 + n_randint(state, 200);
+        qbits2 = 2 + n_randint(state, 200);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+
+        m = 1 + n_randint(state, 20);
+        n = 1 + n_randint(state, 20);
+
+        fmpq_poly_init(A);
+        fmpq_poly_init(B);
+        fmpq_poly_init(Q);
+        fmpq_poly_init(R);
+
+        fmpcb_poly_init(a);
+        fmpcb_poly_init(b);
+        fmpcb_poly_init(q);
+        fmpcb_poly_init(r);
+
+        fmpq_poly_randtest(A, state, m, qbits1);
+        fmpq_poly_randtest_not_zero(B, state, n, qbits2);
+
+        fmpq_poly_divrem(Q, R, A, B);
+
+        fmpcb_poly_set_fmpq_poly(a, A, rbits1);
+        fmpcb_poly_set_fmpq_poly(b, B, rbits2);
+
+        fmpcb_poly_divrem(q, r, a, b, rbits3);
+
+        if (!fmpcb_poly_contains_fmpq_poly(q, Q) ||
+             !fmpcb_poly_contains_fmpq_poly(r, R))
+        {
+            printf("FAIL\n\n");
+
+            printf("A = "); fmpq_poly_print(A); printf("\n\n");
+            printf("B = "); fmpq_poly_print(B); printf("\n\n");
+            printf("Q = "); fmpq_poly_print(Q); printf("\n\n");
+            printf("R = "); fmpq_poly_print(R); printf("\n\n");
+
+            printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n");
+            printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n");
+            printf("q = "); fmpcb_poly_printd(q, 15); printf("\n\n");
+            printf("r = "); fmpcb_poly_printd(r, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmpcb_poly_divrem(a, r, a, b, rbits3);
+        if (!fmpcb_poly_equal(a, q))
+        {
+            printf("FAIL (aliasing q, a)\n\n");
+        }
+        fmpcb_poly_set_fmpq_poly(a, A, rbits1);
+
+        fmpcb_poly_divrem(b, r, a, b, rbits3);
+        if (!fmpcb_poly_equal(b, q))
+        {
+            printf("FAIL (aliasing q, b)\n\n");
+            abort();
+        }
+        fmpcb_poly_set_fmpq_poly(b, B, rbits2);
+
+        fmpcb_poly_divrem(q, a, a, b, rbits3);
+        if (!fmpcb_poly_equal(a, r))
+        {
+            printf("FAIL (aliasing r, a)\n\n");
+            abort();
+        }
+        fmpcb_poly_set_fmpq_poly(a, A, rbits1);
+
+        fmpcb_poly_divrem(q, b, a, b, rbits3);
+        if (!fmpcb_poly_equal(b, r))
+        {
+            printf("FAIL (aliasing r, b)\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(A);
+        fmpq_poly_clear(B);
+        fmpq_poly_clear(Q);
+        fmpq_poly_clear(R);
+
+        fmpcb_poly_clear(a);
+        fmpcb_poly_clear(b);
+        fmpcb_poly_clear(q);
+        fmpcb_poly_clear(r);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}