add missing polynomial subtraction methods, make sure to always round when adding/subtracting

doc/source/fmpcb_poly.rst

@@ -138,6 +138,15 @@ Arithmetic
 
     Sets *C* to the sum of *A* and *B*.
 
+.. function:: void _fmpcb_poly_sub(fmpcb_struct * C, const fmpcb_struct * A, long lenA, const fmpcb_struct * B, long lenB, long prec)
+
+    Sets *{C, max(lenA, lenB)}* to the difference of *{A, lenA}* and *{B, lenB}*.
+    Allows aliasing of the input and output operands.
+
+.. function:: void fmpcb_poly_sub(fmpcb_poly_t C, const fmpcb_poly_t A, const fmpcb_poly_t B, long prec)
+
+    Sets *C* to the difference of *A* and *B*.
+
 .. function:: void _fmpcb_poly_mullow_classical(fmpcb_struct * C, const fmpcb_struct * A, long lenA, const fmpcb_struct * B, long lenB, long n, long prec)
 
 .. function:: void _fmpcb_poly_mullow_transpose(fmpcb_struct * C, const fmpcb_struct * A, long lenA, const fmpcb_struct * B, long lenB, long n, long prec)

doc/source/fmprb_poly.rst

@@ -119,6 +119,15 @@ Arithmetic
 
     Sets *C* to the sum of *A* and *B*.
 
+.. function:: void _fmprb_poly_sub(fmprb_struct * C, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long prec)
+
+    Sets *{C, max(lenA, lenB)}* to the difference of *{A, lenA}* and *{B, lenB}*.
+    Allows aliasing of the input and output operands.
+
+.. function:: void fmprb_poly_sub(fmprb_poly_t C, const fmprb_poly_t A, const fmprb_poly_t B, long prec)
+
+    Sets *C* to the difference of *A* and *B*.
+
 .. function:: void _fmprb_poly_mullow_classical(fmprb_struct * C, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long n, long prec)
 
 .. function:: void _fmprb_poly_mullow_ztrunc(fmprb_struct * C, const fmprb_struct * A, long lenA, const fmprb_struct * B, long lenB, long n, long prec)

fmpcb.h

@@ -131,6 +131,13 @@ fmpcb_set_round(fmpcb_t z, const fmpcb_t x, long prec)
     fmprb_set_round(fmpcb_imagref(z), fmpcb_imagref(x), prec);
 }
 
+static __inline__ void
+fmpcb_neg_round(fmpcb_t z, const fmpcb_t x, long prec)
+{
+    fmprb_neg_round(fmpcb_realref(z), fmpcb_realref(x), prec);
+    fmprb_neg_round(fmpcb_imagref(z), fmpcb_imagref(x), prec);
+}
+
 static __inline__ void
 fmpcb_swap(fmpcb_t z, fmpcb_t x)
 {

fmpcb_poly.h

@@ -119,6 +119,12 @@ void _fmpcb_poly_add(fmpcb_struct * res, const fmpcb_struct * poly1, long len1,
 void fmpcb_poly_add(fmpcb_poly_t res, const fmpcb_poly_t poly1,
               const fmpcb_poly_t poly2, long prec);
 
+void _fmpcb_poly_sub(fmpcb_struct * res, const fmpcb_struct * poly1, long len1,
+    const fmpcb_struct * poly2, long len2, long prec);
+
+void fmpcb_poly_sub(fmpcb_poly_t res, const fmpcb_poly_t poly1,
+              const fmpcb_poly_t poly2, long prec);
+
 void fmpcb_poly_mullow_classical(fmpcb_poly_t res, const fmpcb_poly_t poly1,
                                             const fmpcb_poly_t poly2,
                                                 long n, long prec);

fmpcb_poly/add.c

@@ -34,12 +34,11 @@ _fmpcb_poly_add(fmpcb_struct * res, const fmpcb_struct * poly1, long len1,
     for (i = 0; i < min; i++)
         fmpcb_add(res + i, poly1 + i, poly2 + i, prec);
 
-    /* TODO: round? */
     for (i = min; i < len1; i++)
-        fmpcb_set(res + i, poly1 + i);
+        fmpcb_set_round(res + i, poly1 + i, prec);
 
     for (i = min; i < len2; i++)
-        fmpcb_set(res + i, poly2 + i);
+        fmpcb_set_round(res + i, poly2 + i, prec);
 }
 
 void

fmpcb_poly/sub.c

@@ -0,0 +1,58 @@
+/*=============================================================================
+
+    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_sub(fmpcb_struct * res, const fmpcb_struct * poly1, long len1,
+    const fmpcb_struct * poly2, long len2, long prec)
+{
+    long i, min = FLINT_MIN(len1, len2);
+
+    for (i = 0; i < min; i++)
+        fmpcb_sub(res + i, poly1 + i, poly2 + i, prec);
+
+    for (i = min; i < len1; i++)
+        fmpcb_set_round(res + i, poly1 + i, prec);
+
+    for (i = min; i < len2; i++)
+        fmpcb_neg_round(res + i, poly2 + i, prec);
+}
+
+void
+fmpcb_poly_sub(fmpcb_poly_t res, const fmpcb_poly_t poly1,
+              const fmpcb_poly_t poly2, long prec)
+{
+    long max = FLINT_MAX(poly1->length, poly2->length);
+
+    fmpcb_poly_fit_length(res, max);
+
+    _fmpcb_poly_sub(res->coeffs, poly1->coeffs, poly1->length, poly2->coeffs,
+                   poly2->length, prec);
+
+    _fmpcb_poly_set_length(res, max);
+    _fmpcb_poly_normalise(res);
+}
+

fmpcb_poly/test/t-add.c

@@ -0,0 +1,120 @@
+/*=============================================================================
+
+    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("add....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* compare with fmpq_poly */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        long qbits1, qbits2, rbits1, rbits2, rbits3, trunc;
+        fmpq_poly_t A, B, C;
+        fmpcb_poly_t a, b, c, d;
+
+        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);
+        trunc = n_randint(state, 10);
+
+        fmpq_poly_init(A);
+        fmpq_poly_init(B);
+        fmpq_poly_init(C);
+
+        fmpcb_poly_init(a);
+        fmpcb_poly_init(b);
+        fmpcb_poly_init(c);
+        fmpcb_poly_init(d);
+
+        fmpq_poly_randtest(A, state, 1 + n_randint(state, 10), qbits1);
+        fmpq_poly_randtest(B, state, 1 + n_randint(state, 10), qbits2);
+
+        fmpq_poly_add(C, A, B);
+
+        fmpcb_poly_set_fmpq_poly(a, A, rbits1);
+        fmpcb_poly_set_fmpq_poly(b, B, rbits2);
+
+        fmpcb_poly_add(c, a, b, rbits3);
+
+        if (!fmpcb_poly_contains_fmpq_poly(c, C))
+        {
+            printf("FAIL\n\n");
+            printf("bits3 = %ld\n", rbits3);
+            printf("trunc = %ld\n", trunc);
+
+            printf("A = "); fmpq_poly_print(A); printf("\n\n");
+            printf("B = "); fmpq_poly_print(B); printf("\n\n");
+            printf("C = "); fmpq_poly_print(C); printf("\n\n");
+
+            printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n");
+            printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n");
+            printf("c = "); fmpcb_poly_printd(c, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmpcb_poly_set(d, a);
+        fmpcb_poly_add(d, d, b, rbits3);
+        if (!fmpcb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        fmpcb_poly_set(d, b);
+        fmpcb_poly_add(d, a, d, rbits3);
+        if (!fmpcb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(A);
+        fmpq_poly_clear(B);
+        fmpq_poly_clear(C);
+
+        fmpcb_poly_clear(a);
+        fmpcb_poly_clear(b);
+        fmpcb_poly_clear(c);
+        fmpcb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

fmpcb_poly/test/t-sub.c

@@ -0,0 +1,120 @@
+/*=============================================================================
+
+    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("sub....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* compare with fmpq_poly */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        long qbits1, qbits2, rbits1, rbits2, rbits3, trunc;
+        fmpq_poly_t A, B, C;
+        fmpcb_poly_t a, b, c, d;
+
+        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);
+        trunc = n_randint(state, 10);
+
+        fmpq_poly_init(A);
+        fmpq_poly_init(B);
+        fmpq_poly_init(C);
+
+        fmpcb_poly_init(a);
+        fmpcb_poly_init(b);
+        fmpcb_poly_init(c);
+        fmpcb_poly_init(d);
+
+        fmpq_poly_randtest(A, state, 1 + n_randint(state, 10), qbits1);
+        fmpq_poly_randtest(B, state, 1 + n_randint(state, 10), qbits2);
+
+        fmpq_poly_sub(C, A, B);
+
+        fmpcb_poly_set_fmpq_poly(a, A, rbits1);
+        fmpcb_poly_set_fmpq_poly(b, B, rbits2);
+
+        fmpcb_poly_sub(c, a, b, rbits3);
+
+        if (!fmpcb_poly_contains_fmpq_poly(c, C))
+        {
+            printf("FAIL\n\n");
+            printf("bits3 = %ld\n", rbits3);
+            printf("trunc = %ld\n", trunc);
+
+            printf("A = "); fmpq_poly_print(A); printf("\n\n");
+            printf("B = "); fmpq_poly_print(B); printf("\n\n");
+            printf("C = "); fmpq_poly_print(C); printf("\n\n");
+
+            printf("a = "); fmpcb_poly_printd(a, 15); printf("\n\n");
+            printf("b = "); fmpcb_poly_printd(b, 15); printf("\n\n");
+            printf("c = "); fmpcb_poly_printd(c, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmpcb_poly_set(d, a);
+        fmpcb_poly_sub(d, d, b, rbits3);
+        if (!fmpcb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        fmpcb_poly_set(d, b);
+        fmpcb_poly_sub(d, a, d, rbits3);
+        if (!fmpcb_poly_equal(d, c))
+        {
+            printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        fmpq_poly_clear(A);
+        fmpq_poly_clear(B);
+        fmpq_poly_clear(C);
+
+        fmpcb_poly_clear(a);
+        fmpcb_poly_clear(b);
+        fmpcb_poly_clear(c);
+        fmpcb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

fmprb.h

@@ -130,6 +130,13 @@ fmprb_neg(fmprb_t x, const fmprb_t y)
     fmpr_set(fmprb_radref(x), fmprb_radref(y));
 }
 
+static __inline__ void
+fmprb_neg_round(fmprb_t x, const fmprb_t y, long prec)
+{
+    fmprb_set_round(x, y, prec);
+    fmprb_neg(x, x);
+}
+
 static __inline__ void
 fmprb_abs(fmprb_t x, const fmprb_t y)
 {

fmprb_poly.h

@@ -131,6 +131,12 @@ _fmprb_poly_add(fmprb_struct * res, const fmprb_struct * poly1, long len1,
 void fmprb_poly_add(fmprb_poly_t res, const fmprb_poly_t poly1,
               const fmprb_poly_t poly2, long prec);
 
+void _fmprb_poly_sub(fmprb_struct * res, const fmprb_struct * poly1, long len1,
+    const fmprb_struct * poly2, long len2, long prec);
+
+void fmprb_poly_sub(fmprb_poly_t res, const fmprb_poly_t poly1,
+              const fmprb_poly_t poly2, long prec);
+
 void _fmprb_poly_mullow_ztrunc(fmprb_struct * res,
     const fmprb_struct * poly1, long len1,
     const fmprb_struct * poly2, long len2, long n, long prec);

fmprb_poly/add.c

@@ -34,12 +34,11 @@ _fmprb_poly_add(fmprb_struct * res, const fmprb_struct * poly1, long len1,
     for (i = 0; i < min; i++)
         fmprb_add(res + i, poly1 + i, poly2 + i, prec);
 
-    /* TODO: round? */
     for (i = min; i < len1; i++)
-        fmprb_set(res + i, poly1 + i);
+        fmprb_set_round(res + i, poly1 + i, prec);
 
     for (i = min; i < len2; i++)
-        fmprb_set(res + i, poly2 + i);
+        fmprb_set_round(res + i, poly2 + i, prec);
 }
 
 void
@@ -56,3 +55,4 @@ fmprb_poly_add(fmprb_poly_t res, const fmprb_poly_t poly1,
     _fmprb_poly_set_length(res, max);
     _fmprb_poly_normalise(res);
 }
+

fmprb_poly/sub.c

@@ -0,0 +1,58 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+void
+_fmprb_poly_sub(fmprb_struct * res, const fmprb_struct * poly1, long len1,
+    const fmprb_struct * poly2, long len2, long prec)
+{
+    long i, min = FLINT_MIN(len1, len2);
+
+    for (i = 0; i < min; i++)
+        fmprb_sub(res + i, poly1 + i, poly2 + i, prec);
+
+    for (i = min; i < len1; i++)
+        fmprb_set_round(res + i, poly1 + i, prec);
+
+    for (i = min; i < len2; i++)
+        fmprb_neg_round(res + i, poly2 + i, prec);
+}
+
+void
+fmprb_poly_sub(fmprb_poly_t res, const fmprb_poly_t poly1,
+              const fmprb_poly_t poly2, long prec)
+{
+    long max = FLINT_MAX(poly1->length, poly2->length);
+
+    fmprb_poly_fit_length(res, max);
+
+    _fmprb_poly_sub(res->coeffs, poly1->coeffs, poly1->length, poly2->coeffs,
+                   poly2->length, prec);
+
+    _fmprb_poly_set_length(res, max);
+    _fmprb_poly_normalise(res);
+}
+

fmprb_poly/test/t-sub.c

@@ -0,0 +1,96 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("sub....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        long qbits1, qbits2, rbits1, rbits2, rbits3;
+        fmpq_poly_t A, B, C;
+        fmprb_poly_t a, b, c;
+
+        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);
+
+        fmpq_poly_init(A);
+        fmpq_poly_init(B);
+        fmpq_poly_init(C);
+
+        fmprb_poly_init(a);
+        fmprb_poly_init(b);
+        fmprb_poly_init(c);
+
+        fmpq_poly_randtest(A, state, 1 + n_randint(state, 10), qbits1);
+        fmpq_poly_randtest(B, state, 1 + n_randint(state, 10), qbits2);
+        fmpq_poly_sub(C, A, B);
+
+        fmprb_poly_set_fmpq_poly(a, A, rbits1);
+        fmprb_poly_set_fmpq_poly(b, B, rbits2);
+        fmprb_poly_sub(c, a, b, rbits3);
+
+        if (!fmprb_poly_contains_fmpq_poly(c, C))
+        {
+            printf("FAIL\n\n");
+            printf("bits3 = %ld\n", rbits3);
+
+            printf("A = "); fmpq_poly_print(A); printf("\n\n");
+            printf("B = "); fmpq_poly_print(B); printf("\n\n");
+            printf("C = "); fmpq_poly_print(C); printf("\n\n");
+
+            printf("a = "); fmprb_poly_printd(a, 15); printf("\n\n");
+            printf("b = "); fmprb_poly_printd(b, 15); printf("\n\n");
+            printf("c = "); fmprb_poly_printd(c, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmpq_poly_clear(A);
+        fmpq_poly_clear(B);
+        fmpq_poly_clear(C);
+
+        fmprb_poly_clear(a);
+        fmprb_poly_clear(b);
+        fmprb_poly_clear(c);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}