diff --git a/Makefile.in b/Makefile.in index d0d990c632..9b8c9cf37d 100644 --- a/Makefile.in +++ b/Makefile.in @@ -182,7 +182,7 @@ HEADER_DIRS := \ fq_zech fq_zech_mat fq_zech_poly \ n_poly \ fq fq_vec fq_mat fq_poly \ - padic padic_mat padic_poly \ + padic padic_nmod padic_mat padic_poly \ qadic \ \ fmpz_extras fmpz_factor fmpzi \ diff --git a/doc/source/index.rst b/doc/source/index.rst index e433f11521..38895dc384 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -297,6 +297,7 @@ p-adic numbers padic.rst padic_poly.rst padic_mat.rst + padic_nmod.rst qadic.rst Floating-point support code diff --git a/doc/source/padic_nmod.rst b/doc/source/padic_nmod.rst new file mode 100644 index 0000000000..7c41eccacd --- /dev/null +++ b/doc/source/padic_nmod.rst @@ -0,0 +1,161 @@ +.. _padic-nmod: + +**padic_nmod.h** -- floating point p-adic numbers +=============================================================================== + + +Introduction +-------------------------------------------------------------------------------- + +The ``padic_nmod_t`` data type represents elements of `\mathbf{Q}_p` to +precision `N`, stored in the form `x = p^v u` with `u \in \mathbf{Z}` and +`v \in \mathbf{Z} / p^N \mathbf{Z}`. +Arithmetic operations can be carried out with respect to a context +containing the prime number `p` and various pieces of pre-computed data. +Multiple precisions `N` are possible. + +Independent of the context, we consider a `p`-adic number +`x = u p^v` to be in canonical form whenever either +`p \nmid u` or `u = v = 0`. + +We briefly describe the interface: + +The functions in this module expect arguments of type ``padic_nmod_t``, and +each variable carries its own precision. The functions have an interface that +is similar to the MPFR functions. In particular, they have the same semantics, +specified as follows: Compute the requested operation exactly and then reduce +the result to the precision of the output variable. + +Data structures +-------------------------------------------------------------------------------- + +A `p`-adic number of type ``padic_nmod_t`` comprises a mantissa `man`, and a +valuation `val`. + + +Context +-------------------------------------------------------------------------------- + +A context object for `p`-adic arithmetic contains data pertinent to `p`-adic +computations, but which we choose not to store with each element individually. +Currently, this includes the prime number `p`, its inverse and the quotient of +`UWORD_MAX` by `p`, precomputed powers of `p` in the range between `1` and `N`, +and the integers modulo `p` and `p^N`. + +.. function:: int padic_nmod_ctx_init(gr_ctx_t ctx, ulong p, slong n) + + Initialises the context ``ctx`` with the given data. + + Assumes that `p` is a prime. This is not verified but the subsequent + behaviour is undefined if `p` is a composite number. + +.. function:: void padic_nmod_ctx_clear(gr_ctx_t ctx) + + +Memory management +-------------------------------------------------------------------------------- + + +.. function:: void padic_nmod_init(padic_nmod_t res, gr_ctx_t ctx) + + Initialises the `p`-adic number. + +.. function:: void _padic_nmod_canonicalise(padic_nmod_t x, gr_ctx_t ctx) + + Brings the `p`-adic number ``x`` into canonical form. + + That is to say, ensures that either `u = v = 0` or + `p \nmid u`. There is no reduction modulo a power + of `p`. + + +Randomisation +-------------------------------------------------------------------------------- + + +.. function:: int padic_nmod_randtest(padic_nmod_t rop, flint_rand_t state, gr_ctx_t ctx) + + Sets ``rop`` to a random `p`-adic number. + +.. function:: int padic_nmod_randtest_not_zero(padic_nmod_t rop, flint_rand_t state, gr_ctx_t ctx) + + Sets ``rop`` to a random non-zero `p`-adic number. + + +Assignments and conversions +-------------------------------------------------------------------------------- + +All assignment functions set the value of ``res`` from ``x``, +reduced to the precision of ``ctx``. + +.. function:: int padic_nmod_set(padic_nmod_t res, const padic_nmod_t x, gr_ctx_t ctx) + + Sets ``res`` to the `p`-adic number ``x``. + +.. function:: int padic_nmod_set_ui(padic_nmod_t res, ulong x, gr_ctx_t ctx) + + Sets the `p`-adic number ``res`` to the ``ulong`` + integer ``x``. + +.. function:: int padic_nmod_zero(padic_nmod_t res, gr_ctx_t ctx) + + Sets the `p`-adic number ``res`` to zero. + +.. function:: int padic_nmod_one(padic_nmod_t res, gr_ctx_t ctx) + + Sets the `p`-adic number ``res`` to one, reduced modulo the + precision of ``res``. + + +Comparison +-------------------------------------------------------------------------------- + + +.. function:: truth_t padic_nmod_is_zero(const padic_nmod_t x, gr_ctx_t ctx) + + Returns whether ``x`` is equal to zero. + +.. function:: truth_t padic_nmod_is_one(const padic_nmod_t x, gr_ctx_t ctx) + + Returns whether ``x`` is equal to one, that is, whether + `u = 1` and `v = 0`. + + +Arithmetic operations +-------------------------------------------------------------------------------- + + +.. function:: int padic_nmod_add(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, gr_ctx_t ctx) + + Sets ``res`` to the sum of ``a`` and ``b``. + +.. function:: int padic_nmod_div(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, gr_ctx_t ctx) + + Sets ``res`` to the quotient of ``a`` and ``b``. + +.. function:: int padic_nmod_inv(padic_nmod_t res, const padic_nmod_t a, gr_ctx_t ctx) + + Sets ``res`` to the inverse of ``a``. + +.. function:: int padic_nmod_mul(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, gr_ctx_t ctx) + + Sets ``res`` to the product of ``a`` and ``b``. + +.. function:: int padic_nmod_neg(padic_nmod_t res, const padic_nmod_t a, gr_ctx_t ctx) + + Sets ``res`` to the additive inverse of ``a``. + +.. function:: int padic_nmod_sub(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, gr_ctx_t ctx) + + Sets ``res`` to the difference of ``a`` and ``b``. + + +Input and output +-------------------------------------------------------------------------------- + + +.. function:: void padic_nmod_println(const padic_nmod_t x, gr_ctx_t ctx) + + Prints the string representation of the `p`-adic number ``x`` + to the stream ``stdout``. + diff --git a/doc/source/ulong_extras.rst b/doc/source/ulong_extras.rst index 8ffe2721fb..92fdf5f0c7 100644 --- a/doc/source/ulong_extras.rst +++ b/doc/source/ulong_extras.rst @@ -1344,6 +1344,19 @@ Factorisation `p` we make repeated use of :func:`n_divrem2_precomp` until division by `p` is no longer possible. +.. function:: int n_remove2_prime_inv(ulong * n, ulong p, ulong inv1, ulong inv2) + + Removes the highest possible power of `p` from `n`, replacing `n` with the + quotient. The return value is the highest power of `p` that divided `n`. + Assumes `n` is not `0`. + + For `p = 2` trailing zeroes are counted. For other primes `p` we require + `inv1` to be set to a precomputed inverse of `p` computed with + :func:`n_binvert`, and `inv2` to be the quotient of ``UWORD_MAX`` by `p`. + Then repeated divisions by `p` are done until it is no longer possible. + + Assumes that `p` is prime. + .. function:: void n_factor_insert(n_factor_t * factors, ulong p, ulong exp) Inserts the given prime power factor ``p^exp`` into ``factors``. diff --git a/src/padic/randtest.c b/src/padic/randtest.c index 013e7c42d0..dd3f108a91 100644 --- a/src/padic/randtest.c +++ b/src/padic/randtest.c @@ -13,7 +13,8 @@ #define PADIC_RANDTEST_TRIES 10 -void padic_randtest(padic_t rop, flint_rand_t state, const padic_ctx_t ctx) +void +padic_randtest(padic_t rop, flint_rand_t state, const padic_ctx_t ctx) { const slong N = padic_prec(rop); @@ -23,7 +24,7 @@ void padic_randtest(padic_t rop, flint_rand_t state, const padic_ctx_t ctx) if (N > 0) { - min = - ((N + 9) / 10); + min = -((N + 9) / 10); max = N; } else if (N < 0) @@ -46,15 +47,17 @@ void padic_randtest(padic_t rop, flint_rand_t state, const padic_ctx_t ctx) fmpz_clear(pow); } -void padic_randtest_not_zero(padic_t rop, flint_rand_t state, - const padic_ctx_t ctx) +void +padic_randtest_not_zero(padic_t rop, flint_rand_t state, const padic_ctx_t ctx) { - slong i; + slong i = 1; - padic_randtest(rop, state, ctx); - - for (i = 1; !padic_is_zero(rop) && i < PADIC_RANDTEST_TRIES; i++) + do + { padic_randtest(rop, state, ctx); + i++; + } + while (padic_is_zero(rop) && i < PADIC_RANDTEST_TRIES); if (padic_is_zero(rop)) { @@ -63,8 +66,8 @@ void padic_randtest_not_zero(padic_t rop, flint_rand_t state, } } -void padic_randtest_int(padic_t rop, flint_rand_t state, - const padic_ctx_t ctx) +void +padic_randtest_int(padic_t rop, flint_rand_t state, const padic_ctx_t ctx) { const slong N = padic_prec(rop); diff --git a/src/padic_nmod.h b/src/padic_nmod.h new file mode 100644 index 0000000000..e250a0e511 --- /dev/null +++ b/src/padic_nmod.h @@ -0,0 +1,151 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#ifndef PADIC_NMOD_H +#define PADIC_NMOD_H + +#ifdef PADIC_NMOD_INLINES_C +#define PADIC_NMOD_INLINE +#else +#define PADIC_NMOD_INLINE static inline +#endif + +#include "gr.h" +#include "nmod.h" +#include "padic_types.h" + +#ifdef __cplusplus +extern "C" { +#endif + +#define PADIC_EMIN (-(WORD_MAX / 4)) +#define PADIC_EMAX ((WORD_MAX / 4)) + +#define PADIC_NMOD_CTX(ctx) \ +((_padic_nmod_ctx_struct *)(GR_CTX_DATA_AS_PTR(ctx))) +#define PADIC_NMOD_CTX_P(ctx) (PADIC_NMOD_CTX(ctx)->p) +#define PADIC_NMOD_CTX_PINV1(ctx) (PADIC_NMOD_CTX(ctx)->pinv1) +#define PADIC_NMOD_CTX_PINV2(ctx) (PADIC_NMOD_CTX(ctx)->pinv2) +#define PADIC_NMOD_CTX_N(ctx) (PADIC_NMOD_CTX(ctx)->n) +#define PADIC_NMOD_CTX_POW(ctx) (PADIC_NMOD_CTX(ctx)->pow) +#define PADIC_NMOD_CTX_P_MOD(ctx) (PADIC_NMOD_CTX(ctx)->p_mod) +#define PADIC_NMOD_CTX_PN_MOD(ctx) (PADIC_NMOD_CTX(ctx)->pn_mod) + +/* Context *******************************************************************/ + +int padic_nmod_ctx_init(gr_ctx_t ctx, ulong p, slong n); +void padic_nmod_ctx_clear(gr_ctx_t ctx); + +/* Memory management *********************************************************/ + +void padic_nmod_init(padic_nmod_t res, gr_ctx_t ctx); + +PADIC_NMOD_INLINE +void _padic_nmod_canonicalise(padic_nmod_t x, gr_ctx_t ctx) +{ + if (x->u != 0) + { + if (x->u > 0) + { + x->v += n_remove2_prime_inv(&x->u, PADIC_NMOD_CTX_P(ctx), + PADIC_NMOD_CTX_PINV1(ctx), + PADIC_NMOD_CTX_PINV2(ctx)); + } + + else + { + ulong z = - x->u; + slong e = n_remove2_prime_inv(&z, PADIC_NMOD_CTX_P(ctx), + PADIC_NMOD_CTX_PINV1(ctx), + PADIC_NMOD_CTX_PINV2(ctx)); + + if (e > 0) + { + x->u = - (slong) z; + } + + x->v += e; + } + } + + else + { + x->v = 0; + } +} + +/* Randomisation *************************************************************/ + +int padic_nmod_randtest(padic_nmod_t rop, flint_rand_t state, gr_ctx_t ctx); + +int padic_nmod_randtest_not_zero(padic_nmod_t rop, flint_rand_t state, + gr_ctx_t ctx); + +/* Assignments and conversions ***********************************************/ + +int padic_nmod_set(padic_nmod_t res, const padic_nmod_t x, gr_ctx_t ctx); + +int padic_nmod_set_ui(padic_nmod_t res, ulong x, gr_ctx_t ctx); + +PADIC_NMOD_INLINE +int padic_nmod_zero(padic_nmod_t res, gr_ctx_t ctx) +{ + res->u = 0; + res->v = PADIC_EMAX; + + return GR_SUCCESS; +} + +PADIC_NMOD_INLINE +int padic_nmod_one(padic_nmod_t res, gr_ctx_t ctx) +{ + res->u = 1; + res->v = 0; + + return GR_SUCCESS; +} + +/* Comparison ****************************************************************/ + +PADIC_NMOD_INLINE +truth_t padic_nmod_is_zero(const padic_nmod_t x, gr_ctx_t ctx) +{ + return (x->u == 0) ? T_TRUE : T_FALSE; +} + +PADIC_NMOD_INLINE +truth_t padic_nmod_is_one(const padic_nmod_t x, gr_ctx_t ctx) +{ + return (x->u == 1 && x->v == 0) ? T_TRUE : T_FALSE; +} + +/* Arithmetic operations *****************************************************/ + +int padic_nmod_add(padic_nmod_t res, const padic_nmod_t a, + const padic_nmod_t b, gr_ctx_t ctx); +int padic_nmod_div(padic_nmod_t res, const padic_nmod_t a, + const padic_nmod_t b, gr_ctx_t ctx); +int padic_nmod_inv(padic_nmod_t res, const padic_nmod_t a, gr_ctx_t ctx); +int padic_nmod_mul(padic_nmod_t res, const padic_nmod_t a, + const padic_nmod_t b, gr_ctx_t ctx); +int padic_nmod_neg(padic_nmod_t res, const padic_nmod_t a, gr_ctx_t ctx); +int padic_nmod_sub(padic_nmod_t res, const padic_nmod_t a, + const padic_nmod_t b, gr_ctx_t ctx); + +/* Input and output **********************************************************/ + +void padic_nmod_println(const padic_nmod_t x, gr_ctx_t ctx); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/padic_nmod/add.c b/src/padic_nmod/add.c new file mode 100644 index 0000000000..7d3e0dc855 --- /dev/null +++ b/src/padic_nmod/add.c @@ -0,0 +1,65 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +int +padic_nmod_add(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, + gr_ctx_t ctx) +{ + if (a->u == 0) + { + padic_nmod_set(res, b, ctx); + } + + else if (b->u == 0) + { + padic_nmod_set(res, a, ctx); + } + + else + { + if (a->v == b->v) + { + res->u = nmod_add(a->u, b->u, PADIC_NMOD_CTX_PN_MOD(ctx)); + res->v = a->v; + + _padic_nmod_canonicalise(res, ctx); + /* Underflow */ + if (res->v > PADIC_EMAX) + return GR_UNABLE; + } + + else if (a->v < b->v) + { + ulong f; + + f = (b->v - a->v > PADIC_NMOD_CTX_N(ctx)) + ? PADIC_NMOD_CTX_POW(ctx)[PADIC_NMOD_CTX_N(ctx) - 1] + : PADIC_NMOD_CTX_POW(ctx)[b->v - a->v - 1]; + res->u = nmod_addmul(a->u, f, b->u, PADIC_NMOD_CTX_PN_MOD(ctx)); + res->v = a->v; + } + + else /* a->v > b->v */ + { + ulong f; + + f = (a->v - b->v > PADIC_NMOD_CTX_N(ctx)) + ? PADIC_NMOD_CTX_POW(ctx)[PADIC_NMOD_CTX_N(ctx) - 1] + : PADIC_NMOD_CTX_POW(ctx)[a->v - b->v - 1]; + res->u = nmod_addmul(b->u, f, a->u, PADIC_NMOD_CTX_PN_MOD(ctx)); + res->v = b->v; + } + } + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/ctx_clear.c b/src/padic_nmod/ctx_clear.c new file mode 100644 index 0000000000..a51588beb1 --- /dev/null +++ b/src/padic_nmod/ctx_clear.c @@ -0,0 +1,19 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +void +padic_nmod_ctx_clear(gr_ctx_t ctx) +{ + flint_free(PADIC_NMOD_CTX_POW(ctx)); + flint_free(GR_CTX_DATA_AS_PTR(ctx)); +} diff --git a/src/padic_nmod/ctx_init.c b/src/padic_nmod/ctx_init.c new file mode 100644 index 0000000000..3f4aeb728a --- /dev/null +++ b/src/padic_nmod/ctx_init.c @@ -0,0 +1,57 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +/* Not checked: p is prime */ +int +padic_nmod_ctx_init(gr_ctx_t ctx, ulong p, slong n) +{ + slong i; + ulong hi; + + if (n <= 0 || n >= FLINT_BITS) + return GR_UNABLE; + + /* Compute p^n and verify that this doesn't overflow a ulong */ + GR_CTX_DATA_AS_PTR(ctx) = flint_malloc(sizeof(_padic_nmod_ctx_struct)); + PADIC_NMOD_CTX_POW(ctx) = (ulong *) flint_malloc(n * sizeof(ulong)); + + PADIC_NMOD_CTX_POW(ctx)[0] = p; + + for (i = 1; i < n; i++) + { + umul_ppmm(hi, PADIC_NMOD_CTX_POW(ctx)[i], + PADIC_NMOD_CTX_POW(ctx)[i - 1], p); + + if (hi != 0) + { + padic_nmod_ctx_clear(ctx); + return GR_UNABLE; + } + } + + PADIC_NMOD_CTX_P(ctx) = p; + PADIC_NMOD_CTX_N(ctx) = n; + PADIC_NMOD_CTX_PINV1(ctx) = n_binvert(p); + PADIC_NMOD_CTX_PINV2(ctx) = UWORD_MAX / p; + + nmod_init(&PADIC_NMOD_CTX_P_MOD(ctx), p); + nmod_init(&PADIC_NMOD_CTX_PN_MOD(ctx), PADIC_NMOD_CTX_POW(ctx)[n - 1]); + + for (i = 0; i < n; i++) + { + NMOD_RED(PADIC_NMOD_CTX_POW(ctx)[i], PADIC_NMOD_CTX_POW(ctx)[i], + PADIC_NMOD_CTX_PN_MOD(ctx)); + } + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/div.c b/src/padic_nmod/div.c new file mode 100644 index 0000000000..79e4731917 --- /dev/null +++ b/src/padic_nmod/div.c @@ -0,0 +1,36 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +int +padic_nmod_div(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, + gr_ctx_t ctx) +{ + if (b->u == 0) + { + flint_throw(FLINT_ERROR, "Exception (padic_nmod_div). b is zero.\n"); + + return GR_UNABLE; + } + + if (a->u == 0) + return padic_nmod_zero(res, ctx); + + res->u = nmod_div(a->u, b->u, PADIC_NMOD_CTX_PN_MOD(ctx)); + res->v = a->v - b->v; + + /* Overflow or underflow */ + if (res->v < PADIC_EMIN || res->v > PADIC_EMAX) + return GR_UNABLE; + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/init.c b/src/padic_nmod/init.c new file mode 100644 index 0000000000..1d22317aae --- /dev/null +++ b/src/padic_nmod/init.c @@ -0,0 +1,19 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +void +padic_nmod_init(padic_nmod_t res, gr_ctx_t ctx) +{ + res->u = 0; + res->v = PADIC_EMAX; +} diff --git a/src/padic_nmod/inv.c b/src/padic_nmod/inv.c new file mode 100644 index 0000000000..c0a19733b4 --- /dev/null +++ b/src/padic_nmod/inv.c @@ -0,0 +1,28 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +int +padic_nmod_inv(padic_nmod_t res, const padic_nmod_t a, gr_ctx_t ctx) +{ + if (a->u == 0) + { + flint_throw(FLINT_ERROR, "Exception (padic_nmod_inv). a is zero.\n"); + + return GR_UNABLE; + } + + res->u = nmod_inv(a->u, PADIC_NMOD_CTX_PN_MOD(ctx)); + res->v = -a->v; + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/io.c b/src/padic_nmod/io.c new file mode 100644 index 0000000000..6c73b24fc7 --- /dev/null +++ b/src/padic_nmod/io.c @@ -0,0 +1,26 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +void +padic_nmod_println(const padic_nmod_t x, gr_ctx_t ctx) +{ + if (x->u == 0) + { + flint_printf("0\n"); + } + + else + { + flint_printf("%wu*%wu^%wd\n", x->u, PADIC_NMOD_CTX_P(ctx), x->v); + } +} diff --git a/src/padic_nmod/mul.c b/src/padic_nmod/mul.c new file mode 100644 index 0000000000..b6342b5e4b --- /dev/null +++ b/src/padic_nmod/mul.c @@ -0,0 +1,30 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +int +padic_nmod_mul(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, + gr_ctx_t ctx) +{ + if (a->u == 0 || b->u == 0) + return padic_nmod_zero(res, ctx); + + res->u = nmod_mul(a->u, b->u, PADIC_NMOD_CTX_PN_MOD(ctx)); + res->v = a->v + b->v; + + /* Overflow or underflow */ + if (res->v < PADIC_EMIN || res->v > PADIC_EMAX) + return GR_UNABLE; + + else + return GR_SUCCESS; +} diff --git a/src/padic_nmod/neg.c b/src/padic_nmod/neg.c new file mode 100644 index 0000000000..71e66318ed --- /dev/null +++ b/src/padic_nmod/neg.c @@ -0,0 +1,29 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +int +padic_nmod_neg(padic_nmod_t res, const padic_nmod_t a, gr_ctx_t ctx) +{ + if (a->u == 0) + { + padic_nmod_zero(res, ctx); + } + + else + { + res->v = a->v; + res->u = nmod_neg(a->u, PADIC_NMOD_CTX_PN_MOD(ctx)); + } + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/randtest.c b/src/padic_nmod/randtest.c new file mode 100644 index 0000000000..123a8d6dfd --- /dev/null +++ b/src/padic_nmod/randtest.c @@ -0,0 +1,56 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "long_extras.h" +#include "padic_nmod.h" + +#define PADIC_NMOD_RANDTEST_TRIES 10 + +int +padic_nmod_randtest(padic_nmod_t rop, flint_rand_t state, gr_ctx_t ctx) +{ + rop->u = _n_randint(state, PADIC_NMOD_CTX_PN_MOD(ctx).n); + + if (!rop->u) + { + rop->v = 0; + } + else + { + rop->v = z_randint(state, PADIC_EMAX + 1); + _padic_nmod_canonicalise(rop, ctx); + } + + return GR_SUCCESS; +} + +int +padic_nmod_randtest_not_zero(padic_nmod_t rop, flint_rand_t state, + gr_ctx_t ctx) +{ + slong i = 1; + + do + { + padic_nmod_randtest(rop, state, ctx); + i++; + } + while (padic_nmod_is_zero(rop, ctx) == T_TRUE + && i < PADIC_NMOD_RANDTEST_TRIES); + + if (padic_nmod_is_zero(rop, ctx) == T_TRUE) + { + rop->u = 1; + rop->v = PADIC_EMAX; + } + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/set.c b/src/padic_nmod/set.c new file mode 100644 index 0000000000..e3648fda90 --- /dev/null +++ b/src/padic_nmod/set.c @@ -0,0 +1,26 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +int +padic_nmod_set(padic_nmod_t res, const padic_nmod_t x, gr_ctx_t ctx) +{ + res->u = x->u; + res->v = x->v; + + _padic_nmod_canonicalise(res, ctx); + + if (res->v > PADIC_EMAX) /* Underflow */ + return GR_UNABLE; + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/set_ui.c b/src/padic_nmod/set_ui.c new file mode 100644 index 0000000000..4af0fb04cc --- /dev/null +++ b/src/padic_nmod/set_ui.c @@ -0,0 +1,34 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +int +padic_nmod_set_ui(padic_nmod_t res, ulong x, gr_ctx_t ctx) +{ + if (x < PADIC_NMOD_CTX_P(ctx)) + { + res->u = x; + res->v = (x == 0) ? PADIC_EMAX : 0; + } + else + { + ulong p; + + p = PADIC_NMOD_CTX_P(ctx); + + res->v = n_remove2_prime_inv(&x, p, PADIC_NMOD_CTX_PINV1(ctx), + PADIC_NMOD_CTX_PINV2(ctx)); + res->u = nmod_set_ui(x, PADIC_NMOD_CTX_PN_MOD(ctx)); + } + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/sub.c b/src/padic_nmod/sub.c new file mode 100644 index 0000000000..cf26760b73 --- /dev/null +++ b/src/padic_nmod/sub.c @@ -0,0 +1,69 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "padic_nmod.h" + +int +padic_nmod_sub(padic_nmod_t res, const padic_nmod_t a, const padic_nmod_t b, + gr_ctx_t ctx) +{ + if (a->u == 0) + { + padic_nmod_neg(res, b, ctx); + } + + else if (b->u == 0) + { + padic_nmod_set(res, a, ctx); + } + + else + { + if (a->v == b->v) + { + res->u = nmod_sub(a->u, b->u, PADIC_NMOD_CTX_PN_MOD(ctx)); + res->v = a->v; + + _padic_nmod_canonicalise(res, ctx); + + if (res->v > PADIC_EMAX) /* Underflow */ + return GR_UNABLE; + } + + else if (a->v < b->v) + { + nmod_t pn = PADIC_NMOD_CTX_PN_MOD(ctx); + ulong f, n; + + f = (b->v - a->v > PADIC_NMOD_CTX_N(ctx)) + ? PADIC_NMOD_CTX_POW(ctx)[PADIC_NMOD_CTX_N(ctx) - 1] + : PADIC_NMOD_CTX_POW(ctx)[b->v - a->v - 1]; + n = nmod_neg(b->u, pn); + res->u = nmod_addmul(a->u, f, n, pn); + res->v = a->v; + } + + else /* a->v > b->v */ + { + nmod_t pn = PADIC_NMOD_CTX_PN_MOD(ctx); + ulong f; + + f = (a->v - b->v > PADIC_NMOD_CTX_N(ctx)) + ? PADIC_NMOD_CTX_POW(ctx)[PADIC_NMOD_CTX_N(ctx) - 1] + : PADIC_NMOD_CTX_POW(ctx)[a->v - b->v - 1]; + res->u = nmod_neg(b->u, pn); + res->u = nmod_addmul(res->u, f, a->u, pn); + res->v = b->v; + } + } + + return GR_SUCCESS; +} diff --git a/src/padic_nmod/test/main.c b/src/padic_nmod/test/main.c new file mode 100644 index 0000000000..fabe69aad6 --- /dev/null +++ b/src/padic_nmod/test/main.c @@ -0,0 +1,28 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "t-add.c" +#include "t-div.c" +#include "t-inv.c" +#include "t-mul.c" +#include "t-neg.c" +#include "t-sub.c" + +test_struct tests[] = { + TEST_FUNCTION(padic_nmod_add), + TEST_FUNCTION(padic_nmod_div), + TEST_FUNCTION(padic_nmod_inv), + TEST_FUNCTION(padic_nmod_mul), + TEST_FUNCTION(padic_nmod_neg), + TEST_FUNCTION(padic_nmod_sub) +}; + +TEST_MAIN(tests) diff --git a/src/padic_nmod/test/t-add.c b/src/padic_nmod/test/t-add.c new file mode 100644 index 0000000000..1c9c4b4008 --- /dev/null +++ b/src/padic_nmod/test/t-add.c @@ -0,0 +1,118 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "padic.h" +#include "padic_nmod.h" + +TEST_FUNCTION_START(padic_nmod_add, state) +{ + /* Check that addition coincides with the padic_t one */ + for (int i = 0; i < 1000 * flint_test_multiplier(); i++) + { + ulong p = n_randprime(state, FLINT_BITS / 4, 1); + ulong max_val = FLINT_BITS / n_clog(p, 2); + + fmpz_t p_fmpz; + fmpz_init(p_fmpz); + fmpz_set_ui(p_fmpz, p); + + padic_ctx_t ctx_padic; + gr_ctx_t ctx_nmod; + + padic_ctx_init(ctx_padic, p_fmpz, 0, max_val, PADIC_VAL_UNIT); + padic_nmod_ctx_init(ctx_nmod, p, max_val); + + padic_t a, b, res; + padic_nmod_t c, d, res_float; + + padic_init2(a, max_val / 3); + padic_init2(b, max_val / 3); + padic_init2(res, max_val); + padic_nmod_init(c, ctx_nmod); + padic_nmod_init(d, ctx_nmod); + padic_nmod_init(res_float, ctx_nmod); + + padic_randtest(a, state, ctx_padic); + padic_randtest(b, state, ctx_padic); + padic_nmod_set_ui(c, fmpz_get_ui(padic_unit(a)), ctx_nmod); + c->v = padic_get_val(a); + _padic_nmod_canonicalise(c, ctx_nmod); + padic_nmod_set_ui(d, fmpz_get_ui(padic_unit(b)), ctx_nmod); + d->v = padic_get_val(b); + _padic_nmod_canonicalise(d, ctx_nmod); + + padic_add(res, a, b, ctx_padic); + padic_nmod_add(res_float, c, d, ctx_nmod); + + if (!fmpz_equal_ui(padic_unit(a), c->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nc = "); + padic_nmod_println(c, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + if (!fmpz_equal_ui(padic_unit(b), d->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("b = "); + padic_print(b, ctx_padic); + flint_printf("\nd = "); + padic_nmod_println(d, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + if ((fmpz_fdiv_ui(padic_unit(res), + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n) != res_float->u) + || (padic_get_val(res) != res_float->v)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nb = "); + padic_print(b, ctx_padic); + flint_printf("\nc = "); + padic_nmod_println(c, ctx_nmod); + flint_printf("d = "); + padic_nmod_println(d, ctx_nmod); + flint_printf("res = "); + padic_print(res, ctx_padic); + flint_printf("\nres_float = "); + padic_nmod_println(res_float, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + + padic_clear(a); + padic_clear(b); + padic_clear(res); + + padic_ctx_clear(ctx_padic); + padic_nmod_ctx_clear(ctx_nmod); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/padic_nmod/test/t-div.c b/src/padic_nmod/test/t-div.c new file mode 100644 index 0000000000..0a1d8718f3 --- /dev/null +++ b/src/padic_nmod/test/t-div.c @@ -0,0 +1,133 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "padic.h" +#include "padic_nmod.h" + +TEST_FUNCTION_START(padic_nmod_div, state) +{ + /* Check that division coincides with the padic_t one */ + for (int i = 0; i < 1000 * flint_test_multiplier(); i++) + { + ulong p = n_randprime(state, FLINT_BITS / 4, 1); + ulong max_val = FLINT_BITS / n_clog(p, 2); + + fmpz_t p_fmpz; + fmpz_init(p_fmpz); + fmpz_set_ui(p_fmpz, p); + + padic_ctx_t ctx_padic; + gr_ctx_t ctx_nmod; + + padic_ctx_init(ctx_padic, p_fmpz, 0, max_val + 1, PADIC_VAL_UNIT); + padic_nmod_ctx_init(ctx_nmod, p, max_val); + + padic_t a, b, a_temp, b_temp, res; + padic_nmod_t c, d, res_float; + + padic_init2(a, max_val / 3); + padic_init2(b, max_val / 3); + padic_init2(a_temp, max_val); + padic_init2(b_temp, max_val); + padic_init2(res, max_val); + padic_nmod_init(c, ctx_nmod); + padic_nmod_init(d, ctx_nmod); + padic_nmod_init(res_float, ctx_nmod); + + padic_randtest(a, state, ctx_padic); + padic_randtest_not_zero(b, state, ctx_padic); + padic_nmod_set_ui(c, fmpz_get_ui(padic_unit(a)), ctx_nmod); + c->v = padic_get_val(a); + _padic_nmod_canonicalise(c, ctx_nmod); + padic_nmod_set_ui(d, fmpz_get_ui(padic_unit(b)), ctx_nmod); + d->v = padic_get_val(b); + _padic_nmod_canonicalise(d, ctx_nmod); + padic_set(a_temp, a, ctx_padic); + padic_set(b_temp, b, ctx_padic); + + padic_div(res, a_temp, b_temp, ctx_padic); + padic_nmod_div(res_float, c, d, ctx_nmod); + + if (!fmpz_equal_ui(padic_unit(a), c->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nc = "); + padic_nmod_println(c, ctx_nmod); + flint_printf("p = %wd\n", p); + flint_printf("precision = %wd\n", max_val); + flint_printf("p^precision = %wd\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + if (!fmpz_equal_ui(padic_unit(b), d->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("b = "); + padic_print(b, ctx_padic); + flint_printf("\nd = "); + padic_nmod_println(d, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + + ulong inv_prec; + + if (c->v < d->v) + inv_prec = max_val; + else if (c->v > d->v + (signed) max_val) + inv_prec = 0; + else + inv_prec = max_val - (c->v - d->v); + + if ((fmpz_fdiv_ui(padic_unit(res), PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n) + != res_float->u % fmpz_get_ui(ctx_padic->pow + inv_prec)) + || ((padic_get_val(res) != res_float->v) + && (!padic_is_zero(res) + || (padic_nmod_is_zero(res_float, ctx_nmod) != T_TRUE)))) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nb = "); + padic_print(b, ctx_padic); + flint_printf("\nc = "); + padic_nmod_println(c, ctx_nmod); + flint_printf("d = "); + padic_nmod_println(d, ctx_nmod); + flint_printf("res = "); + padic_print(res, ctx_padic); + flint_printf("\nres_float = "); + padic_nmod_println(res_float, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + padic_clear(a); + padic_clear(b); + padic_clear(res); + + padic_ctx_clear(ctx_padic); + padic_nmod_ctx_clear(ctx_nmod); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/padic_nmod/test/t-inv.c b/src/padic_nmod/test/t-inv.c new file mode 100644 index 0000000000..eec72d7921 --- /dev/null +++ b/src/padic_nmod/test/t-inv.c @@ -0,0 +1,103 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "padic.h" +#include "padic_nmod.h" + +TEST_FUNCTION_START(padic_nmod_inv, state) +{ + /* Check that inverting coincides with the padic_t one */ + for (int i = 0; i < 1000 * flint_test_multiplier(); i++) + { + ulong p = n_randprime(state, FLINT_BITS / 4, 1); + ulong max_val = FLINT_BITS / n_clog(p, 2); + + fmpz_t p_fmpz; + fmpz_init(p_fmpz); + fmpz_set_ui(p_fmpz, p); + + padic_ctx_t ctx_padic; + gr_ctx_t ctx_nmod; + + padic_ctx_init(ctx_padic, p_fmpz, 0, max_val + 1, PADIC_VAL_UNIT); + padic_nmod_ctx_init(ctx_nmod, p, max_val); + + padic_t a, res; + padic_nmod_t b, res_float; + + padic_init2(a, max_val / 3); + padic_init2(res, max_val); + padic_nmod_init(b, ctx_nmod); + padic_nmod_init(res_float, ctx_nmod); + + padic_randtest_not_zero(a, state, ctx_padic); + padic_nmod_set_ui(b, fmpz_get_ui(padic_unit(a)), ctx_nmod); + b->v = padic_get_val(a); + _padic_nmod_canonicalise(b, ctx_nmod); + + padic_inv(res, a, ctx_padic); + padic_nmod_inv(res_float, b, ctx_nmod); + + if (!fmpz_equal_ui(padic_unit(a), b->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nb = "); + padic_nmod_println(b, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + + ulong inv_prec; + + if (b->v > 0) + inv_prec = max_val; + else if (b->v + max_val < 0) + inv_prec = 0; + else + inv_prec = max_val + b->v; + + if ((fmpz_fdiv_ui(padic_unit(res), PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n) + != res_float->u % fmpz_get_ui(ctx_padic->pow + inv_prec)) + || (padic_get_val(res) != res_float->v)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nb = "); + padic_nmod_println(b, ctx_nmod); + flint_printf("res = "); + padic_print(res, ctx_padic); + flint_printf("\nres_float = "); + padic_nmod_println(res_float, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + + padic_clear(a); + padic_clear(res); + + padic_ctx_clear(ctx_padic); + padic_nmod_ctx_clear(ctx_nmod); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/padic_nmod/test/t-mul.c b/src/padic_nmod/test/t-mul.c new file mode 100644 index 0000000000..3e6736607e --- /dev/null +++ b/src/padic_nmod/test/t-mul.c @@ -0,0 +1,120 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "padic.h" +#include "padic_nmod.h" + +TEST_FUNCTION_START(padic_nmod_mul, state) +{ + /* Check that multiplication coincides with the padic_t one */ + for (int i = 0; i < 1000 * flint_test_multiplier(); i++) + { + ulong p = n_randprime(state, FLINT_BITS / 4, 1); + ulong max_val = FLINT_BITS / n_clog(p, 2); + + fmpz_t p_fmpz; + fmpz_init(p_fmpz); + fmpz_set_ui(p_fmpz, p); + + padic_ctx_t ctx_padic; + gr_ctx_t ctx_nmod; + + padic_ctx_init(ctx_padic, p_fmpz, 0, max_val, PADIC_VAL_UNIT); + padic_nmod_ctx_init(ctx_nmod, p, max_val); + + padic_t a, b, res; + padic_nmod_t c, d, res_float; + + padic_init2(a, max_val / 3); + padic_init2(b, max_val / 3); + padic_init2(res, max_val); + padic_nmod_init(c, ctx_nmod); + padic_nmod_init(d, ctx_nmod); + padic_nmod_init(res_float, ctx_nmod); + + padic_randtest(a, state, ctx_padic); + padic_randtest(b, state, ctx_padic); + padic_nmod_set_ui(c, fmpz_get_ui(padic_unit(a)), ctx_nmod); + c->v = padic_get_val(a); + _padic_nmod_canonicalise(c, ctx_nmod); + padic_nmod_set_ui(d, fmpz_get_ui(padic_unit(b)), ctx_nmod); + d->v = padic_get_val(b); + _padic_nmod_canonicalise(d, ctx_nmod); + + padic_mul(res, a, b, ctx_padic); + padic_nmod_mul(res_float, c, d, ctx_nmod); + + if (!fmpz_equal_ui(padic_unit(a), c->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nc = "); + padic_nmod_println(c, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + if (!fmpz_equal_ui(padic_unit(b), d->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("b = "); + padic_print(b, ctx_padic); + flint_printf("\nd = "); + padic_nmod_println(d, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + if ((fmpz_fdiv_ui(padic_unit(res), PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n) + != res_float->u) + || ((padic_get_val(res) != res_float->v) + && (!padic_is_zero(res) + || (padic_nmod_is_zero(res_float, ctx_nmod) != T_TRUE)))) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nb = "); + padic_print(b, ctx_padic); + flint_printf("\nc = "); + padic_nmod_println(c, ctx_nmod); + flint_printf("d = "); + padic_nmod_println(d, ctx_nmod); + flint_printf("res = "); + padic_print(res, ctx_padic); + flint_printf("\nres_float = "); + padic_nmod_println(res_float, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + + padic_clear(a); + padic_clear(b); + padic_clear(res); + + padic_ctx_clear(ctx_padic); + padic_nmod_ctx_clear(ctx_nmod); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/padic_nmod/test/t-neg.c b/src/padic_nmod/test/t-neg.c new file mode 100644 index 0000000000..723dba556e --- /dev/null +++ b/src/padic_nmod/test/t-neg.c @@ -0,0 +1,92 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "padic.h" +#include "padic_nmod.h" + +TEST_FUNCTION_START(padic_nmod_neg, state) +{ + /* Check that negation coincides with the padic_t one */ + for (int i = 0; i < 1000 * flint_test_multiplier(); i++) + { + ulong p = n_randprime(state, FLINT_BITS / 4, 1); + ulong max_val = FLINT_BITS / n_clog(p, 2); + + fmpz_t p_fmpz; + fmpz_init(p_fmpz); + fmpz_set_ui(p_fmpz, p); + + padic_ctx_t ctx_padic; + gr_ctx_t ctx_nmod; + + padic_ctx_init(ctx_padic, p_fmpz, 0, max_val, PADIC_VAL_UNIT); + padic_nmod_ctx_init(ctx_nmod, p, max_val); + + padic_t a, res; + padic_nmod_t b, res_float; + + padic_init2(a, max_val / 3); + padic_init2(res, max_val); + padic_nmod_init(b, ctx_nmod); + padic_nmod_init(res_float, ctx_nmod); + + padic_randtest_not_zero(a, state, ctx_padic); + padic_nmod_set_ui(b, fmpz_get_ui(padic_unit(a)), ctx_nmod); + b->v = padic_get_val(a); + _padic_nmod_canonicalise(b, ctx_nmod); + + padic_neg(res, a, ctx_padic); + padic_nmod_neg(res_float, b, ctx_nmod); + + if (!fmpz_equal_ui(padic_unit(a), b->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nb = "); + padic_nmod_println(b, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + if ((fmpz_fdiv_ui(padic_unit(res), PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n) + != res_float->u) || (padic_get_val(res) != res_float->v)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nb = "); + padic_nmod_println(b, ctx_nmod); + flint_printf("res = "); + padic_print(res, ctx_padic); + flint_printf("\nres_float = "); + padic_nmod_println(res_float, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + + padic_clear(a); + padic_clear(res); + + padic_ctx_clear(ctx_padic); + padic_nmod_ctx_clear(ctx_nmod); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/padic_nmod/test/t-sub.c b/src/padic_nmod/test/t-sub.c new file mode 100644 index 0000000000..1e5a97900e --- /dev/null +++ b/src/padic_nmod/test/t-sub.c @@ -0,0 +1,131 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "test_helpers.h" +#include "padic.h" +#include "padic_nmod.h" + +TEST_FUNCTION_START(padic_nmod_sub, state) +{ + /* Check that substraction coincides with the padic_t one */ + for (int i = 0; i < 1000 * flint_test_multiplier(); i++) + { + ulong p = n_randprime(state, FLINT_BITS / 4, 1); + ulong max_val = FLINT_BITS / n_clog(p, 2); + + fmpz_t p_fmpz; + fmpz_init(p_fmpz); + fmpz_set_ui(p_fmpz, p); + + padic_ctx_t ctx_padic; + gr_ctx_t ctx_nmod; + + padic_ctx_init(ctx_padic, p_fmpz, 0, max_val + 1, PADIC_VAL_UNIT); + padic_nmod_ctx_init(ctx_nmod, p, max_val); + + padic_t a, b, res; + padic_nmod_t c, d, res_float; + + padic_init2(a, max_val / 3); + padic_init2(b, max_val / 3); + padic_init2(res, max_val); + padic_nmod_init(c, ctx_nmod); + padic_nmod_init(d, ctx_nmod); + padic_nmod_init(res_float, ctx_nmod); + + padic_randtest(a, state, ctx_padic); + padic_randtest(b, state, ctx_padic); + padic_nmod_set_ui(c, fmpz_get_ui(padic_unit(a)), ctx_nmod); + c->v = padic_get_val(a); + _padic_nmod_canonicalise(c, ctx_nmod); + padic_nmod_set_ui(d, fmpz_get_ui(padic_unit(b)), ctx_nmod); + d->v = padic_get_val(b); + _padic_nmod_canonicalise(d, ctx_nmod); + + padic_sub(res, a, b, ctx_padic); + padic_nmod_sub(res_float, c, d, ctx_nmod); + + if (!fmpz_equal_ui(padic_unit(a), c->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nc = "); + padic_nmod_println(c, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + if (!fmpz_equal_ui(padic_unit(b), d->u)) + { + flint_printf("FAIL:\n\n"); + flint_printf("b = "); + padic_print(b, ctx_padic); + flint_printf("\nd = "); + padic_nmod_println(d, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + + ulong inv_prec; + + if (res_float->v == FLINT_MIN(c->v, d->v)) + inv_prec = max_val; + else if (res_float->v > FLINT_MIN(c->v, d->v) + (signed) max_val) + inv_prec = 0; + else + inv_prec = max_val - (res_float->v - FLINT_MIN(c->v, d->v)); + + if ((fmpz_fdiv_ui(padic_unit(res), + fmpz_get_ui(ctx_padic->pow + inv_prec)) + != res_float->u % fmpz_get_ui(ctx_padic->pow + inv_prec)) + || ((padic_get_val(res) != res_float->v) + && (!padic_is_zero(res) + || (padic_nmod_is_zero(res_float, ctx_nmod) != T_TRUE)))) + { + flint_printf("FAIL:\n\n"); + flint_printf("a = "); + padic_print(a, ctx_padic); + flint_printf("\nb = "); + padic_print(b, ctx_padic); + flint_printf("\nc = "); + padic_nmod_println(c, ctx_nmod); + flint_printf("d = "); + padic_nmod_println(d, ctx_nmod); + flint_printf("res = "); + padic_print(res, ctx_padic); + flint_printf("\nres_float = "); + padic_nmod_println(res_float, ctx_nmod); + flint_printf("p = %wu\n", p); + flint_printf("precision = %wu\n", max_val); + flint_printf("p^precision = %wu\n", + PADIC_NMOD_CTX_PN_MOD(ctx_nmod).n); + fflush(stdout); + flint_abort(); + } + + padic_clear(a); + padic_clear(b); + padic_clear(res); + + padic_ctx_clear(ctx_padic); + padic_nmod_ctx_clear(ctx_nmod); + } + + TEST_FUNCTION_END(state); +} diff --git a/src/padic_types.h b/src/padic_types.h index 2c53d502fb..17ef9272fe 100644 --- a/src/padic_types.h +++ b/src/padic_types.h @@ -27,6 +27,15 @@ typedef struct typedef padic_struct padic_t[1]; +typedef struct +{ + ulong u; + slong v; +} +padic_nmod_struct; + +typedef padic_nmod_struct padic_nmod_t[1]; + enum padic_print_mode { PADIC_TERSE, @@ -46,6 +55,18 @@ typedef struct typedef padic_ctx_struct padic_ctx_t[1]; +typedef struct +{ + ulong p; + ulong pinv1; + ulong pinv2; + slong n; + ulong * pow; + nmod_t p_mod; + nmod_t pn_mod; +} +_padic_nmod_ctx_struct; + typedef struct { slong n; diff --git a/src/ulong_extras.h b/src/ulong_extras.h index 52314c5f52..d4de40081c 100644 --- a/src/ulong_extras.h +++ b/src/ulong_extras.h @@ -605,6 +605,7 @@ int n_factor_pollard_brent(ulong *factor, flint_rand_t state, ulong n_in, ulong int n_remove(ulong * n, ulong p); int n_remove2_precomp(ulong * n, ulong p, double ppre); +int n_remove2_prime_inv(ulong * n, ulong p, ulong inv1, ulong inv2); /* ECM functions *************************************************************/ diff --git a/src/ulong_extras/remove2_prime_inv.c b/src/ulong_extras/remove2_prime_inv.c new file mode 100644 index 0000000000..9355be3009 --- /dev/null +++ b/src/ulong_extras/remove2_prime_inv.c @@ -0,0 +1,37 @@ +/* + Copyright (C) 2026 Rubén Muñoz--Bertrand + + This file is part of FLINT. + + FLINT 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 3 of the License, or + (at your option) any later version. See . +*/ + +#include "ulong_extras.h" + +int +n_remove2_prime_inv(ulong * n, ulong p, ulong inv1, ulong inv2) +{ + int val = 0; + + if (p == 2) + { + val = flint_ctz(*n); + + if (val) + (*n) >>= val; + } + + else + { + while (n_divisible_odd_gm(*n, inv1, inv2)) + { + (*n) *= inv1; + val ++; + } + } + + return val; +}