Merge pull request #1770 from flintlib/gr_eval

Evaluation limits for gr qqbar context; use generics in more places

doc/source/gr_domains.rst

@@ -154,6 +154,16 @@ Base rings and fields
     Initializes *ctx* to the field of real or complex algebraic
     numbers with elements of type :type:`qqbar_t`.
 
+.. function:: void _gr_ctx_qqbar_set_limits(gr_ctx_t ctx, slong deg_limit, slong bits_limit)
+
+    Limit degrees of intermediate operands of a *qqbar* context
+    to *deg_limit* and their bit sizes to *bits_limit* (approximately).
+    The limits refer to the sizes of resultants prior to
+    factorization (see :func:`qqbar_binop_within_limits`), so for example
+    adding two degree-100 algebraic numbers
+    requires a degree limit of at least 10000.
+    Warning: currently not all methods respect these limits.
+
 .. function:: void gr_ctx_init_real_arb(gr_ctx_t ctx, slong prec)
               void gr_ctx_init_complex_acb(gr_ctx_t ctx, slong prec)
 

doc/source/gr_generic.rst

@@ -144,7 +144,9 @@ Generic arithmetic
     Sets *res* to the value of the integer polynomial *f* evaluated
     at the argument *x*.
 
-.. function:: int gr_fmpz_mpoly_evaluate(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx)
+.. function:: int gr_fmpz_mpoly_evaluate_iter(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx)
+              int gr_fmpz_mpoly_evaluate_horner(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx)
+              int gr_fmpz_mpoly_evaluate(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx)
 
     Sets *res* to value of the multivariate polynomial *f* (with
     corresponding context object *mctx*) evaluated at the vector

src/ca_field/build_ideal.c

@@ -19,8 +19,7 @@
 #include "qqbar.h"
 #include "fmpz_mpoly.h"
 
-int qqbar_mul_checked(qqbar_t res, const qqbar_t x, const qqbar_t y, slong deg_limit, slong bits_limit);
-int qqbar_pow_fmpz_checked(qqbar_t res, const qqbar_t x, const fmpz_t y, slong deg_limit, slong bits_limit);
+#include "gr.h"
 
 slong
 acb_multi_lindep(fmpz_mat_t rel, acb_srcptr vec, slong len, int check, slong prec)
@@ -573,27 +572,33 @@ ca_field_prove_multiplicative_relation(ca_field_t K, const fmpz * rel,
 
         qqbar_one(a);
 
+        /* qqbar arithmetic with (bogus) limits */
+        gr_ctx_t QQbar;
+        gr_ctx_init_complex_qqbar(QQbar);
+        _gr_ctx_qqbar_set_limits(QQbar, ctx->options[CA_OPT_QQBAR_DEG_LIMIT], ctx->options[CA_OPT_PREC_LIMIT] * 10);
+
         for (i = 0; i < num_powers; i++)
         {
             if (!fmpz_is_zero(rel + i))
             {
                 ca_ext_srcptr ext = CA_FIELD_EXT_ELEM(K, powers[i]);
 
-                /* xxx: bogus limits */
-                if (!qqbar_pow_fmpz_checked(b, CA_EXT_QQBAR(ext), rel + i, ctx->options[CA_OPT_QQBAR_DEG_LIMIT], ctx->options[CA_OPT_PREC_LIMIT] * 10))
+                if (gr_pow_fmpz(b, CA_EXT_QQBAR(ext), rel + i, QQbar) != GR_SUCCESS)
                 {
                     success = 0;
                     goto qqbar_end;
                 }
 
-                if (!qqbar_mul_checked(a, a, b, ctx->options[CA_OPT_QQBAR_DEG_LIMIT], ctx->options[CA_OPT_PREC_LIMIT] * 10))
+                if (gr_mul(a, a, b, QQbar) != GR_SUCCESS)
                 {
                     success = 0;
                     goto qqbar_end;
                 }
             }
         }
 
+        gr_ctx_clear(QQbar);
+
         /* (-1)^ */
         if (fmpz_is_odd(rel + num_powers))
             qqbar_neg(a, a);

src/fmpz_mpoly_q/evaluate_acb.c

@@ -10,65 +10,27 @@
 */
 
 #include "acb.h"
+#include "fmpz_mpoly.h"
 #include "fmpz_mpoly_q.h"
+#include "gr.h"
+#include "gr_generic.h"
 
-/* stupid algorithm, just to have something working... */
 void
 fmpz_mpoly_evaluate_acb(acb_t res, const fmpz_mpoly_t pol, acb_srcptr x, slong prec, const fmpz_mpoly_ctx_t ctx)
 {
-    slong i, j, len, nvars;
-    acb_t s, t, u;
-    ulong * exp; /* todo: support fmpz exponents */
+    gr_ctx_t CC;
+    gr_ctx_init_complex_acb(CC, prec);
 
-    len = fmpz_mpoly_length(pol, ctx);
-
-    if (len == 0)
+    if (pol->length <= 6 && pol->bits <= FLINT_BITS)
     {
-        acb_zero(res);
-        return;
-    }
-
-    if (len == 1 && fmpz_mpoly_is_fmpz(pol, ctx))
-    {
-        acb_set_round_fmpz(res, pol->coeffs, prec);
-        return;
+        if (gr_fmpz_mpoly_evaluate_iter(res, pol, x, ctx, CC) != GR_SUCCESS)
+            acb_indeterminate(res);
     }
-
-    nvars = ctx->minfo->nvars;
-    exp = flint_malloc(sizeof(ulong) * nvars);
-
-    acb_init(s);
-    acb_init(t);
-    acb_init(u);
-
-    for (i = 0; i < len; i++)
+    else
     {
-        fmpz_mpoly_get_term_exp_ui(exp, pol, i, ctx);
-
-        acb_one(t);
-
-        for (j = 0; j < nvars; j++)
-        {
-            if (exp[j] == 1)
-            {
-                acb_mul(t, t, x + j, prec);
-            }
-            else if (exp[j] >= 2)
-            {
-                acb_pow_ui(u, x + j, exp[j], prec);
-                acb_mul(t, t, u, prec);
-            }
-        }
-
-        acb_addmul_fmpz(s, t, pol->coeffs + i, prec);
+        if (gr_fmpz_mpoly_evaluate_horner(res, pol, x, ctx, CC) != GR_SUCCESS)
+            acb_indeterminate(res);
     }
-
-    acb_swap(res, s);
-
-    flint_free(exp);
-    acb_clear(s);
-    acb_clear(t);
-    acb_clear(u);
 }
 
 void

src/gr.h

@@ -1369,6 +1369,7 @@ void gr_ctx_init_nmod32(gr_ctx_t ctx, unsigned int n);
 
 void gr_ctx_init_real_qqbar(gr_ctx_t ctx);
 void gr_ctx_init_complex_qqbar(gr_ctx_t ctx);
+void _gr_ctx_qqbar_set_limits(gr_ctx_t ctx, slong deg_limit, slong bits_limit);
 
 void gr_ctx_init_real_arb(gr_ctx_t ctx, slong prec);
 void gr_ctx_init_complex_acb(gr_ctx_t ctx, slong prec);

src/gr/qqbar.c

@@ -32,13 +32,33 @@ gr_qqbar_ctx;
 
 #define QQBAR_CTX(ring_ctx) ((gr_qqbar_ctx *)(ring_ctx))
 
+void
+_gr_ctx_qqbar_set_limits(gr_ctx_t ctx, slong deg_limit, slong bits_limit)
+{
+    QQBAR_CTX(ctx)->deg_limit = (deg_limit >= 0) ? deg_limit : WORD_MAX;
+    QQBAR_CTX(ctx)->bits_limit = (bits_limit >= 0) ? bits_limit : WORD_MAX;
+}
+
 int
 _gr_qqbar_ctx_write(gr_stream_t out, gr_ctx_t ctx)
 {
     if (QQBAR_CTX(ctx)->real_only)
         gr_stream_write(out, "Real algebraic numbers (qqbar)");
     else
         gr_stream_write(out, "Complex algebraic numbers (qqbar)");
+
+    if (QQBAR_CTX(ctx)->deg_limit != WORD_MAX)
+    {
+        gr_stream_write(out, ", deg_limit = ");
+        gr_stream_write_si(out, QQBAR_CTX(ctx)->deg_limit);
+    }
+
+    if (QQBAR_CTX(ctx)->bits_limit != WORD_MAX)
+    {
+        gr_stream_write(out, ", bits_limit = ");
+        gr_stream_write_si(out, QQBAR_CTX(ctx)->bits_limit);
+    }
+
     return GR_SUCCESS;
 }
 
@@ -452,6 +472,10 @@ _gr_qqbar_neg(qqbar_t res, const qqbar_t x, const gr_ctx_t ctx)
 int
 _gr_qqbar_add(qqbar_t res, const qqbar_t x, const qqbar_t y, const gr_ctx_t ctx)
 {
+    if (QQBAR_CTX(ctx)->deg_limit != WORD_MAX || QQBAR_CTX(ctx)->bits_limit != WORD_MAX)
+        if (!qqbar_binop_within_limits(x, y, QQBAR_CTX(ctx)->deg_limit, QQBAR_CTX(ctx)->bits_limit))
+            return GR_UNABLE;
+
     qqbar_add(res, x, y);
     return GR_SUCCESS;
 }
@@ -487,6 +511,10 @@ _gr_qqbar_add_fmpq(qqbar_t res, const qqbar_t x, const fmpq_t y, const gr_ctx_t
 int
 _gr_qqbar_sub(qqbar_t res, const qqbar_t x, const qqbar_t y, const gr_ctx_t ctx)
 {
+    if (QQBAR_CTX(ctx)->deg_limit != WORD_MAX || QQBAR_CTX(ctx)->bits_limit != WORD_MAX)
+        if (!qqbar_binop_within_limits(x, y, QQBAR_CTX(ctx)->deg_limit, QQBAR_CTX(ctx)->bits_limit))
+            return GR_UNABLE;
+
     qqbar_sub(res, x, y);
     return GR_SUCCESS;
 }
@@ -522,6 +550,10 @@ _gr_qqbar_sub_fmpq(qqbar_t res, const qqbar_t x, const fmpq_t y, const gr_ctx_t
 int
 _gr_qqbar_mul(qqbar_t res, const qqbar_t x, const qqbar_t y, const gr_ctx_t ctx)
 {
+    if (QQBAR_CTX(ctx)->deg_limit != WORD_MAX || QQBAR_CTX(ctx)->bits_limit != WORD_MAX)
+        if (!qqbar_binop_within_limits(x, y, QQBAR_CTX(ctx)->deg_limit, QQBAR_CTX(ctx)->bits_limit))
+            return GR_UNABLE;
+
     qqbar_mul(res, x, y);
     return GR_SUCCESS;
 }
@@ -577,6 +609,10 @@ _gr_qqbar_div(qqbar_t res, const qqbar_t x, const qqbar_t y, const gr_ctx_t ctx)
     }
     else
     {
+        if (QQBAR_CTX(ctx)->deg_limit != WORD_MAX || QQBAR_CTX(ctx)->bits_limit != WORD_MAX)
+            if (!qqbar_binop_within_limits(x, y, QQBAR_CTX(ctx)->deg_limit, QQBAR_CTX(ctx)->bits_limit))
+                return GR_UNABLE;
+
         qqbar_div(res, x, y);
         return GR_SUCCESS;
     }
@@ -647,6 +683,32 @@ _gr_qqbar_is_invertible(const qqbar_t x, const gr_ctx_t ctx)
 int
 _gr_qqbar_pow_ui(qqbar_t res, const qqbar_t x, ulong exp, const gr_ctx_t ctx)
 {
+    if (QQBAR_CTX(ctx)->bits_limit != WORD_MAX && !(exp == 0 || exp == 1))
+    {
+        slong ebits = FLINT_BIT_COUNT(exp);
+
+        if (qqbar_is_zero(x) || qqbar_is_one(x))
+        {
+            qqbar_set(res, x);
+            return GR_SUCCESS;
+        }
+
+        if (qqbar_is_neg_one(x))
+        {
+            if (exp % 2 == 0)
+                qqbar_one(res);
+            else
+                qqbar_set(res, x);
+            return GR_SUCCESS;
+        }
+
+        if (ebits > SMALL_FMPZ_BITCOUNT_MAX)
+            return GR_UNABLE;
+
+        if ((double) exp * qqbar_height_bits(x) > QQBAR_CTX(ctx)->bits_limit)
+            return GR_UNABLE;
+    }
+
     qqbar_pow_ui(res, x, exp);
     return GR_SUCCESS;
 }
@@ -660,6 +722,32 @@ _gr_qqbar_pow_si(qqbar_t res, const qqbar_t x, slong exp, const gr_ctx_t ctx)
     }
     else
     {
+        if (QQBAR_CTX(ctx)->bits_limit != WORD_MAX && !(exp == 0 || exp == 1 || exp == -1))
+        {
+            slong ebits = FLINT_BIT_COUNT(FLINT_ABS(exp));
+
+            if (qqbar_is_zero(x) || qqbar_is_one(x))
+            {
+                qqbar_set(res, x);
+                return GR_SUCCESS;
+            }
+
+            if (qqbar_is_neg_one(x))
+            {
+                if (exp % 2 == 0)
+                    qqbar_one(res);
+                else
+                    qqbar_set(res, x);
+                return GR_SUCCESS;
+            }
+
+            if (ebits > SMALL_FMPZ_BITCOUNT_MAX)
+                return GR_UNABLE;
+
+            if ((double) FLINT_ABS(exp) * qqbar_height_bits(x) > QQBAR_CTX(ctx)->bits_limit)
+                return GR_UNABLE;
+        }
+
         qqbar_pow_si(res, x, exp);
         return GR_SUCCESS;
     }
@@ -674,6 +762,32 @@ _gr_qqbar_pow_fmpz(qqbar_t res, const qqbar_t x, const fmpz_t exp, const gr_ctx_
     }
     else
     {
+        if (QQBAR_CTX(ctx)->bits_limit != WORD_MAX && !(fmpz_is_zero(exp) || fmpz_is_pm1(exp)))
+        {
+            slong ebits = fmpz_bits(exp);
+
+            if (qqbar_is_zero(x) || qqbar_is_one(x))
+            {
+                qqbar_set(res, x);
+                return GR_SUCCESS;
+            }
+
+            if (qqbar_is_neg_one(x))
+            {
+                if (fmpz_is_even(exp))
+                    qqbar_one(res);
+                else
+                    qqbar_set(res, x);
+                return GR_SUCCESS;
+            }
+
+            if (ebits > SMALL_FMPZ_BITCOUNT_MAX)
+                return GR_UNABLE;
+
+            if ((double) FLINT_ABS(*exp) * qqbar_height_bits(x) > QQBAR_CTX(ctx)->bits_limit)
+                return GR_UNABLE;
+        }
+
         qqbar_pow_fmpz(res, x, exp);
         return GR_SUCCESS;
     }

src/gr_generic.h

@@ -76,6 +76,8 @@ WARN_UNUSED_RESULT int gr_fmpz_poly_evaluate(gr_ptr res, const fmpz_poly_t f, gr
 #endif
 
 #ifdef FMPZ_MPOLY_H
+WARN_UNUSED_RESULT int gr_fmpz_mpoly_evaluate_iter(gr_ptr res, const fmpz_mpoly_t pol, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_fmpz_mpoly_evaluate_horner(gr_ptr res, const fmpz_mpoly_t pol, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_fmpz_mpoly_evaluate(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t mctx, gr_ctx_t ctx);
 #endif
 

src/gr_generic/fmpz_mpoly_evaluate.c

@@ -377,7 +377,7 @@ gr_fmpz_mpoly_evaluate_horner(gr_ptr A, const fmpz_mpoly_t B, gr_srcptr C, const
 
 /* todo: accept fmpz exponents */
 int
-gr_evaluate_fmpz_mpoly_iter(gr_ptr res, const fmpz_mpoly_t pol, gr_srcptr x, const fmpz_mpoly_ctx_t ctx, gr_ctx_t cactx)
+gr_fmpz_mpoly_evaluate_iter(gr_ptr res, const fmpz_mpoly_t pol, gr_srcptr x, const fmpz_mpoly_ctx_t ctx, gr_ctx_t cactx)
 {
     slong i, j, len, nvars;
     gr_ptr s, t, u;
@@ -433,5 +433,8 @@ gr_evaluate_fmpz_mpoly_iter(gr_ptr res, const fmpz_mpoly_t pol, gr_srcptr x, con
 int
 gr_fmpz_mpoly_evaluate(gr_ptr res, const fmpz_mpoly_t f, gr_srcptr x, const fmpz_mpoly_ctx_t ctx, gr_ctx_t cactx)
 {
-    return gr_fmpz_mpoly_evaluate_horner(res, f, x, ctx, cactx);
+    if (f->length <= 1 && f->bits <= FLINT_BITS)
+        return gr_fmpz_mpoly_evaluate_iter(res, f, x, ctx, cactx);
+    else
+        return gr_fmpz_mpoly_evaluate_horner(res, f, x, ctx, cactx);
 }