Implement gr_is_integer and gr_is_rational

[oscarbenjamin]

The gr_is_integer and gr_is_rational functions have been documented since maybe FLINT 3.3 or some time but the functions did not exist.

I'm not sure if it is better to add the functions or just to remove them from the docs but this PR adds the functions with generic implementations based on gr_get_fmpz and gr_get_fmpq.

[oscarbenjamin]

I'm not sure what gr_is_integer and gr_is_rational should do in positive characteristic (e.g. nmod).

[oscarbenjamin]

I don't think the regression check failure is related to this PR although I might be wrong...

Besides the generic implementation of gr_is_integer and gr_is_rational I added an override for padic_radix because CI on the first commit failed when calling gr_get_fmpz on some elements representing very large integers.

Other types could have relatively trivial overrides as well (e.g. fmpz/fmpq/...). I haven't added them because I'm not sure if it would be considered worth it or whether it would actually just be better to remove these functions altogether (i.e. from the docs).

The main thing I am not sure about is what to do for e.g. nmod where I don't know what gr_is_integer should do. One possibility would be to always return T_UNKNOWN in that case but it would need handlers in various places to avoid the effect of the generic get_fmpz -> is_integer fallback added here.

[fredrik-johansson]

Good initiative. I'm also not sure what these should mean in positive characteristic, but they will be useful even if this is not yet specified.

Comments:

• I think these methods want specializing right away for most common types: certainly for acb, acf, arb, arf, complex, fmpz, fmpzi, fmpq, nfloat, nfloat_complex, fmpz_mpoly, fmpz_mpoly_q, fmpz_poly, all of which should be easy using available functions (only nfloat_is_integer is missing, I think, but this is not too hard to add).

• It would also be good to specialize them right away also for gr_poly and gr_mpoly by adding gr_poly_is_integer/rational and gr_mpoly_is_integer/is_rational; see gr_poly_is_zero/gr_poly_is_one gr_mpoly_is_zero/gr_mpoly_is_one.

• padic_radix_is_integer should return T_FALSE for non-exact elements that have at least one known valuation < 0 digit.

[oscarbenjamin]

I've added most types here. At least fmpz_mpoly, fmpz_mpoly_q, fmpz_poly are still missing for now and will need to do an audit of all types.

[oscarbenjamin]

This failed on 32-bit musl

gr_mpoly_ring...
FLINT exception (General error):
   GR_MUST_SUCCEED failure: /home/runner/work/flint/flint/src/gr_mpoly/test/t-ring.cqqbar_abs2...

[oscarbenjamin]

Something strange going on in CI with jobs being cancelled. Maybe a GitHub problem?

There does seem to be one real failure though which is a segfault in gr_special in the MinGW job:

gr_vector_nmod8                                   0.25   (PASS)
gr_vector_nmod32...
gr_vector_nmod32                                  0.29   (PASS)
make: *** [Makefile:795: build/gr_special/test/main.exe_TEST_RUN] Segmentation fault

If that is caused by the changes in the last 3 commits I just pushed then my guess is that it is because of the change to propagate characteristic from the base context in gr_poly etc (7a5af14). There are some calls to gr_ctx_is_finite_characteristic in gr_special that could be affected so maybe they now take a different code path and maybe that has exposed a segfault somewhere.

[oscarbenjamin]

my guess is that it is because of the change to propagate characteristic

Yes, it is:

flint/src/gr_special/test/t-chebyshev.c

Lines 31 to 34 in 00e660a

	if (gr_ctx_has_real_prec(ctx) == T_TRUE || gr_ctx_is_finite_characteristic(ctx) == T_TRUE)
	fmpz_randtest(n, state, 100);
	else
	fmpz_randtest(n, state, 5);

Given a polynomial ring with positive characteristic this creates extremely high degre polynomials and consumes all the memory crashing all the CI runners in different ways.

[oscarbenjamin]

This gr_mpoly_divides case failed in CI:

Standalone reproducer (bug.c)

/*
    Reproducer for the gr_mpoly_divides failure from the M1 CI job.
*/

#include "flint.h"
#include "gr.h"
#include "gr_mpoly.h"

static const char * f_string =
        "321*x2*x4^7*x5*x7*x9^12*x10^6*x11^6*x12^5*x13*x16*x17^13*x19^9 + 426*x2*x4^12*x5^3*x6*x8^2*x9*x11^2*"
        "x13*x14^8*x15^5*x16^5*x17^6*x18*x19^15 + 6281437785975148*x2^6*x3*x5^22*x6^4*x8^8*x9^4*x10^5*x12^3*x"
        "13*x14*x15*x16^3*x19 + 17012181490972889*x2*x4*x8^2*x9*x10^8*x11^10*x13^2*x14^9*x15^20*x16*x17^3*x18"
        "*x19 + 346125923*x1^7*x3^3*x5^2*x6^7*x7^5*x8^4*x9^4*x11^3*x12^2*x13*x14^2*x15^6*x16^2*x17^7*x18 + 49"
        "502*x1^3*x2*x4^9*x5*x6*x8^7*x9*x10^3*x11*x13^4*x15^14*x17*x18^2*x19^7 + 274738771365203773*x1*x2^3*x"
        "4*x5^8*x6^7*x7^7*x8*x10^7*x12*x13*x15^7*x17^5*x18*x19 + 2571584447121746*x1^9*x2^11*x3^3*x9^2*x12*x1"
        "3*x14^3*x15^3*x16^7*x17^9*x19 + 74059568648*x1*x2^5*x3*x4^6*x6^7*x10^8*x13*x14^7*x15^2*x16^4*x17^6 +"
        " 291342068847*x1^10*x2^7*x5^8*x7^9*x8^2*x9^2*x11*x12*x14*x15*x16^3*x18^3 + 104*x1*x2^7*x3*x5^9*x6^7*"
        "x8*x9*x11*x12*x14^11*x17*x19^7 + 682284934668971261*x1*x2^8*x4*x5^2*x7*x9^3*x10^6*x11*x12^2*x13*x14^"
        "7*x18^10*x19^4 + 15*x1^4*x5^4*x6^3*x7^7*x8*x9*x10^6*x11^6*x13^6*x14*x15*x16*x17^2*x18^2 + 12073043*x"
        "1*x3^2*x4^8*x5*x6*x9^4*x10^4*x11*x14^3*x15^14*x17^4*x19 + 33312163*x1^15*x2*x4^2*x5^3*x6*x8*x9*x11^9"
        "*x12*x13^3*x14^3*x16^2*x18*x19 + 1603048836430485117*x1*x2*x4*x5*x6^2*x8^2*x10^10*x11^5*x13^3*x14*x1"
        "5^2*x16^3*x17^3*x18^6*x19^2 + 1466*x1^3*x3*x4*x5^2*x6*x7^7*x8^2*x9*x10*x11^4*x12*x14*x16^3*x17^3*x18"
        "^2*x19 + 480539144861951920*x3^5*x4*x6*x7*x8*x11^11*x12^2*x13*x14*x15*x16^5*x18^3*x19 + 199720953110"
        "0*x1*x2^4*x3*x4^7*x5^5*x6*x8*x11*x14*x16*x17^3*x18^4*x19 + 6323868825078953382*x1*x2*x3*x6*x7*x8^6*x"
        "9*x10^10*x14^2*x15^3 + x4*x5^3*x6^4*x7^4*x8^6*x9*x10*x12^4*x17*x18 + 151922537*x1*x2^4*x6^3*x7^3*x8^"
        "3*x9*x13^3*x15*x19^7 + 2251799813683200*x1*x2*x5*x6^7*x8^3*x11*x12^3*x13*x14*x15*x19 + 6160*x3^5*x4*"
        "x5*x9^2*x13^2*x14*x16^3*x18^5*x19 + 2664*x2*x3^4*x5*x7*x8^7*x16*x17^2*x18*x19^2 + 341301*x3*x4*x5^4*"
        "x6*x8*x9^5*x12*x13^3*x14*x16 + 396348391*x2*x3*x5*x7*x8*x10^7*x11*x15^2*x18^4";

static const char * g_string =
        "179955298162301502*x1^11*x3^13*x5^7*x6*x7^5*x9^4*x10^3*x11*x12*x13*x14^8*x15^3*x16^4*x17 + 649339996"
        "4948923436*x1^6*x2*x3^2*x4*x5*x6*x7^6*x9*x11^7*x12^7*x13*x15*x16^3*x17^10*x18^13 + 80983299986091149"
        "8*x2^2*x3*x4^17*x6*x8^28*x10^3*x11^3*x12^3*x13*x14*x19 + 2109552198411477*x1*x2^10*x3^3*x5^11*x6*x7^"
        "5*x9*x11^8*x12^3*x13*x14*x15^4*x16*x17^7*x19^3 + 2028717211909772*x1^6*x2^2*x3*x4^4*x5^4*x8^2*x9^2*x"
        "10*x11*x12^14*x13*x16^9*x17^6*x18^3*x19^2 + 4227858432*x2^2*x4^3*x5^2*x6^7*x7*x8*x10^2*x11^4*x12^13*"
        "x13^7*x14^6*x15*x17^2*x18*x19^6 + 1946480360075313*x1*x2^13*x3^4*x4^3*x5^2*x8^2*x9^8*x10^7*x12^4*x13"
        "^8*x14*x16*x18^2*x19 + 8806455076436408*x1^4*x2^14*x3^3*x4*x6^2*x12^10*x13^2*x14^14*x16*x17*x18^2*x1"
        "9^2 + 1581*x2^14*x3^2*x4^3*x5*x6^2*x8*x11^6*x12*x13^4*x15^2*x16*x17^11*x19^6 + 14472787893836*x1^6*x"
        "2^6*x4^7*x5^2*x7^2*x8^2*x9^2*x10*x11*x12*x13*x14^7*x15^3*x16*x17^4*x18*x19^6 + 481161814225660*x1^3*"
        "x2^11*x4*x5*x6*x7^11*x8*x10*x11*x12^9*x13*x15^6*x17*x18 + 5810933666*x2*x4*x5^6*x6^13*x8*x10^8*x11^2"
        "*x13^11*x15^3*x17 + 1599062774*x5^7*x6^5*x7^5*x9^5*x14^11*x15^6*x17^3*x18^3*x19^2 + 467746750*x1^5*x"
        "2*x3^3*x4*x5^6*x8^4*x9*x10^12*x11*x12*x13*x14^7*x19^3 + 668060*x1*x6*x7*x11*x12^4*x13^16*x14^7*x15^3"
        "*x16*x17*x18^6*x19^4 + 1324767453*x1^3*x2*x5*x7^7*x10*x12^23*x15*x16^2*x17^4*x18^2 + 34359738367*x1*"
        "x3^3*x4^4*x5*x6^4*x7*x9^3*x10^2*x11*x12^8*x13^2*x15*x16*x17*x18*x19^11 + 49000705395*x2^2*x3^6*x4^4*"
        "x5^10*x6^4*x7^7*x10*x11^4*x12^2*x14*x15*x16*x18 + 4*x1*x2*x3^3*x4*x5*x6*x7*x8*x10^2*x11^12*x13*x14^6"
        "*x16^5*x17*x18*x19^5 + 272*x1^3*x2^5*x3^7*x4^5*x8^2*x11*x12*x14^4*x15*x16*x19^5 + 224*x1^3*x3^7*x4*x"
        "6*x8^2*x9^2*x11*x14*x15^7*x16*x17^6*x18*x19 + 5459472913*x1*x4*x5^12*x8^2*x11^6*x12*x13*x15^3*x17*x1"
        "8^4 + 7222*x2*x6^3*x8*x9^5*x11*x12*x13^4*x14^3*x15*x16*x17*x18*x19";

int
main(void)
{
    const ordering_t orderings[] = {ORD_LEX, ORD_DEGLEX, ORD_DEGREVLEX};
    const char * ordering_names[] = {"lex", "deglex", "degrevlex"};
    gr_ctx_t cctx;
    slong nvars, ord_index, num_threads;
    int failed = 0;

    gr_ctx_init_nmod(cctx, UWORD(9139524932520349247));

    for (nvars = 19; nvars <= 20; nvars++)
    {
        for (ord_index = 0; ord_index < 3; ord_index++)
        {
            gr_mpoly_ctx_t ctx;
            gr_mpoly_t f, g, h, q;
            int status;

            gr_mpoly_ctx_init(ctx, cctx, nvars, orderings[ord_index]);
            gr_mpoly_init(f, ctx);
            gr_mpoly_init(g, ctx);
            gr_mpoly_init(h, ctx);
            gr_mpoly_init(q, ctx);

            status = gr_set_str(f, f_string, ctx);
            status |= gr_set_str(g, g_string, ctx);

            if (status != GR_SUCCESS)
            {
                flint_printf("parse failed: nvars=%wd ordering=%s status=%d\n",
                             nvars, ordering_names[ord_index], status);
            }
            else
            {
                for (num_threads = 2; num_threads <= 4; num_threads++)
                {
                    int mul_status, serial_status, threaded_status, generic_status;

                    flint_set_num_threads(num_threads);
                    mul_status = gr_mpoly_mul(h, f, g, ctx);
                    serial_status = gr_mpoly_divides_heap(q, h, g, ctx);
                    threaded_status = gr_mpoly_divides_heap_threaded(q, h, g, ctx);
                    generic_status = gr_mpoly_divides(q, h, g, ctx);

                    flint_printf("nvars=%wd ordering=%s threads=%wd lengths=%wd,%wd,%wd "
                                 "mul=%d serial=%d threaded=%d generic=%d\n",
                                 nvars, ordering_names[ord_index], num_threads,
                                 f->length, g->length, h->length,
                                 mul_status, serial_status, threaded_status, generic_status);

                    if (mul_status == GR_SUCCESS && serial_status == GR_SUCCESS &&
                        threaded_status == GR_DOMAIN)
                        failed = 1;
                }
            }

            gr_mpoly_clear(f, ctx);
            gr_mpoly_clear(g, ctx);
            gr_mpoly_clear(h, ctx);
            gr_mpoly_clear(q, ctx);
            gr_mpoly_ctx_clear(ctx);
        }
    }

    gr_ctx_clear(cctx);
    flint_cleanup_master();

    return failed ? 1 : 0;
}

I don't think that is caused by the changes in this PR and I can't reproduce on Linux x86-64 with those changes. The CI failure was in the M1 job so I post it here to then test on Mac (M3).

[fredrik-johansson]

This gr_mpoly_divides case failed in CI:
Standalone reproducer (bug.c)

I don't think that is caused by the changes in this PR and I can't reproduce on Linux x86-64 with those changes. The CI failure was in the M1 job so I post it here to then test on Mac (M3).

Thanks for the reproducer. It should be unrelated to this patch. The original multithreaded fmpz_mpoly/nmod_mpoly division code assumes strong memory order and is disabled on arm64 which is weakly ordered. The multithreaded gr_mpoly division code attempts to support arm64 by including atomic guards, but apparently these aren't enough. It's possible that changed RNG state from this PR accidentally made the bug more likely to manifest.