Skip to content

Commit

Permalink
Merge pull request #1607 from flintlib/division
Browse files Browse the repository at this point in the history
Some work on division
  • Loading branch information
fredrik-johansson authored Nov 13, 2023
2 parents 5d2c2a6 + 8baccb8 commit 8548d83
Show file tree
Hide file tree
Showing 19 changed files with 792 additions and 65 deletions.
62 changes: 39 additions & 23 deletions doc/source/gr.rst
Original file line number Diff line number Diff line change
Expand Up @@ -596,27 +596,6 @@ The default implementations of the following methods check for divisors
Particular rings should override the methods when an inversion
or division algorithm is available.

.. function:: int gr_div(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
int gr_div_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
int gr_div_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
int gr_div_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
int gr_div_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
int gr_div_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
int gr_other_div(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)

Sets *res* to the quotient `x / y` if such an element exists
in the present ring. Returns the flag ``GR_DOMAIN`` if no such
quotient exists.
Returns the flag ``GR_UNABLE`` if the implementation is unable
to perform the computation.

When the ring is not a field, the definition of division may
vary depending on the ring. A ring implementation may define
`x / y = x y^{-1}` and return ``GR_DOMAIN`` when `y^{-1}` does not
exist; alternatively, it may attempt to solve the equation
`q y = x` (which, for example, gives the usual exact
division in `\mathbb{Z}`).

.. function:: truth_t gr_is_invertible(gr_srcptr x, gr_ctx_t ctx)

Returns whether *x* has a multiplicative inverse in the present ring,
Expand All @@ -630,9 +609,46 @@ or division algorithm is available.
``GR_UNABLE`` if the implementation is unable to perform
the computation.

.. function:: truth_t gr_divides(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
.. function:: int gr_div(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
int gr_div_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
int gr_div_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
int gr_div_fmpz(gr_ptr res, gr_srcptr x, const fmpz_t y, gr_ctx_t ctx)
int gr_div_fmpq(gr_ptr res, gr_srcptr x, const fmpq_t y, gr_ctx_t ctx)
int gr_div_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx)
int gr_other_div(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx)

Sets *res* to the quotient `x / y`. In a field, this returns
``GR_DOMAIN`` if `y` is zero; in an integral domain,
it returns ``GR_DOMAIN`` if `y` is zero or if the quotient
does not exist. In a non-integral domain, we consider a quotient
to exist only if it is unique, and otherwise return ``GR_DOMAIN``;
see :func:`gr_div_nonunique` for a different behavior.

Returns whether *x* divides *y*.
Returns the flag ``GR_UNABLE`` if the implementation is unable
to perform the computation.

.. function:: int gr_div_nonunique(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)

Sets *res* to an arbitrary solution `q` of the equation `x = q y`.
Returns the flag ``GR_DOMAIN`` if no such solution exists.
Returns the flag ``GR_UNABLE`` if the implementation is unable
to perform the computation.
This method allows dividing `x / y` in some cases where :func:`gr_div` fails:

* `0 / 0` has solutions (for example, 0) in any ring.
* It allows solving division problems in nonintegral domains.
For example, it allows assigning a value to `6 / 2` in
`R = \mathbb{Z}/10\mathbb{Z}` even though `2^{-1}` does not exist
in `R`. In this case, both 3 and 8 are possible solutions, and which
one is chosen is unpredictable.

.. function:: truth_t gr_divides(gr_srcptr d, gr_srcptr x, gr_ctx_t ctx)

Returns whether `d \mid x`; that is, whether there is an element `q`
such that `x = dq`. Note that this corresponds to divisibility
in the sense of :func:`gr_div_nonunique`, which is weaker than that
of :func:`gr_div`. For example, `0 \mid 0` is true even
in rings where `0 / 0` is undefined.

.. function:: int gr_divexact(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
int gr_divexact_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx)
Expand Down
1 change: 1 addition & 0 deletions doc/source/gr_domains.rst
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ Domain properties
truth_t gr_ctx_is_algebraically_closed(gr_ctx_t ctx)
truth_t gr_ctx_is_finite_characteristic(gr_ctx_t ctx)
truth_t gr_ctx_is_ordered_ring(gr_ctx_t ctx)
truth_t gr_ctx_is_zero_ring(gr_ctx_t ctx)

Returns whether the structure satisfies the respective
mathematical property.
Expand Down
4 changes: 4 additions & 0 deletions doc/source/nmod.rst
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,10 @@ Modular reduction and arithmetic
Returns `ab^{-1}` modulo ``mod.n``. The inverse of `b` is assumed to
exist. It is assumed that `a` is already reduced modulo ``mod.n``.

.. function:: int nmod_divides(mp_limb_t * a, mp_limb_t b, mp_limb_t c, nmod_t mod)

If `a\cdot c = b \mod n` has a solution for `a` return `1` and set `a` to such a solution. Otherwise return `0` and leave `a` undefined.

.. function:: mp_limb_t nmod_pow_ui(mp_limb_t a, ulong e, nmod_t mod)

Returns `a^e` modulo ``mod.n``. No assumptions are made about
Expand Down
16 changes: 11 additions & 5 deletions src/gr.h
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,7 @@ typedef enum
GR_METHOD_CTX_IS_FINITE,
GR_METHOD_CTX_IS_FINITE_CHARACTERISTIC,
GR_METHOD_CTX_IS_ALGEBRAICALLY_CLOSED,
GR_METHOD_CTX_IS_ZERO_RING,

/* ring properties related to orderings and norms */
GR_METHOD_CTX_IS_ORDERED_RING,
Expand Down Expand Up @@ -263,6 +264,9 @@ typedef enum
GR_METHOD_DIV_OTHER,
GR_METHOD_OTHER_DIV,

GR_METHOD_DIV_NONUNIQUE,
GR_METHOD_DIVIDES,

GR_METHOD_DIVEXACT,
GR_METHOD_DIVEXACT_UI,
GR_METHOD_DIVEXACT_SI,
Expand All @@ -288,7 +292,6 @@ typedef enum
GR_METHOD_FACTOR_PERFECT_POWER,
GR_METHOD_ROOT_UI,

GR_METHOD_DIVIDES,
GR_METHOD_EUCLIDEAN_DIV,
GR_METHOD_EUCLIDEAN_REM,
GR_METHOD_EUCLIDEAN_DIVREM,
Expand Down Expand Up @@ -955,6 +958,7 @@ GR_INLINE truth_t gr_ctx_is_ring(gr_ctx_t ctx) { return GR_CTX_PREDICATE(ctx, CT
GR_INLINE truth_t gr_ctx_is_commutative_ring(gr_ctx_t ctx) { return GR_CTX_PREDICATE(ctx, CTX_IS_COMMUTATIVE_RING)(ctx); }
GR_INLINE truth_t gr_ctx_is_integral_domain(gr_ctx_t ctx) { return GR_CTX_PREDICATE(ctx, CTX_IS_INTEGRAL_DOMAIN)(ctx); }
GR_INLINE truth_t gr_ctx_is_field(gr_ctx_t ctx) { return GR_CTX_PREDICATE(ctx, CTX_IS_FIELD)(ctx); }
GR_INLINE truth_t gr_ctx_is_zero_ring(gr_ctx_t ctx) { return GR_CTX_PREDICATE(ctx, CTX_IS_ZERO_RING)(ctx); }

GR_INLINE truth_t gr_ctx_is_unique_factorization_domain(gr_ctx_t ctx) { return GR_CTX_PREDICATE(ctx, CTX_IS_UNIQUE_FACTORIZATION_DOMAIN)(ctx); }
GR_INLINE truth_t gr_ctx_is_finite(gr_ctx_t ctx) { return GR_CTX_PREDICATE(ctx, CTX_IS_FINITE)(ctx); }
Expand Down Expand Up @@ -1076,6 +1080,9 @@ GR_INLINE WARN_UNUSED_RESULT int gr_mul_2exp_fmpz(gr_ptr res, gr_srcptr x, const
GR_INLINE WARN_UNUSED_RESULT int gr_set_fmpz_2exp_fmpz(gr_ptr res, const fmpz_t x, const fmpz_t y, gr_ctx_t ctx) { return GR_BINARY_OP_FMPZ_FMPZ(ctx, SET_FMPZ_2EXP_FMPZ)(res, x, y, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_get_fmpz_2exp_fmpz(fmpz_t res1, fmpz_t res2, gr_srcptr x, gr_ctx_t ctx) { return GR_UNARY_OP_GET_FMPZ_FMPZ(ctx, GET_FMPZ_2EXP_FMPZ)(res1, res2, x, ctx); }

GR_INLINE WARN_UNUSED_RESULT int gr_inv(gr_ptr res, gr_srcptr x, gr_ctx_t ctx) { return GR_UNARY_OP(ctx, INV)(res, x, ctx); }
GR_INLINE truth_t gr_is_invertible(gr_srcptr x, gr_ctx_t ctx) { return GR_UNARY_PREDICATE(ctx, IS_INVERTIBLE)(x, ctx); }

GR_INLINE WARN_UNUSED_RESULT int gr_div(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_OP(ctx, DIV)(res, x, y, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_div_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx) { return GR_BINARY_OP_UI(ctx, DIV_UI)(res, x, y, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_div_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx) { return GR_BINARY_OP_SI(ctx, DIV_SI)(res, x, y, ctx); }
Expand All @@ -1084,6 +1091,9 @@ GR_INLINE WARN_UNUSED_RESULT int gr_div_fmpq(gr_ptr res, gr_srcptr x, const fmpq
GR_INLINE WARN_UNUSED_RESULT int gr_div_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx) { return GR_BINARY_OP_OTHER(ctx, DIV_OTHER)(res, x, y, y_ctx, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_other_div(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx) { return GR_OTHER_BINARY_OP(ctx, OTHER_DIV)(res, x, x_ctx, y, ctx); }

GR_INLINE WARN_UNUSED_RESULT int gr_div_nonunique(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_OP(ctx, DIV_NONUNIQUE)(res, x, y, ctx); }
GR_INLINE truth_t gr_divides(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_PREDICATE(ctx, DIVIDES)(x, y, ctx); }

GR_INLINE WARN_UNUSED_RESULT int gr_divexact(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_OP(ctx, DIVEXACT)(res, x, y, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_divexact_ui(gr_ptr res, gr_srcptr x, ulong y, gr_ctx_t ctx) { return GR_BINARY_OP_UI(ctx, DIVEXACT_UI)(res, x, y, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_divexact_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx) { return GR_BINARY_OP_SI(ctx, DIVEXACT_SI)(res, x, y, ctx); }
Expand All @@ -1092,10 +1102,6 @@ GR_INLINE WARN_UNUSED_RESULT int gr_divexact_fmpq(gr_ptr res, gr_srcptr x, const
GR_INLINE WARN_UNUSED_RESULT int gr_divexact_other(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx) { return GR_BINARY_OP_OTHER(ctx, DIVEXACT_OTHER)(res, x, y, y_ctx, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_other_divexact(gr_ptr res, gr_srcptr x, gr_ctx_t x_ctx, gr_srcptr y, gr_ctx_t ctx) { return GR_OTHER_BINARY_OP(ctx, OTHER_DIVEXACT)(res, x, x_ctx, y, ctx); }

GR_INLINE WARN_UNUSED_RESULT int gr_inv(gr_ptr res, gr_srcptr x, gr_ctx_t ctx) { return GR_UNARY_OP(ctx, INV)(res, x, ctx); }
GR_INLINE truth_t gr_is_invertible(gr_srcptr x, gr_ctx_t ctx) { return GR_UNARY_PREDICATE(ctx, IS_INVERTIBLE)(x, ctx); }

GR_INLINE truth_t gr_divides(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_PREDICATE(ctx, DIVIDES)(x, y, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_euclidean_div(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_OP(ctx, EUCLIDEAN_DIV)(res, x, y, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_euclidean_rem(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_OP(ctx, EUCLIDEAN_REM)(res, x, y, ctx); }
GR_INLINE WARN_UNUSED_RESULT int gr_euclidean_divrem(gr_ptr res1, gr_ptr res2, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_BINARY_OP(ctx, EUCLIDEAN_DIVREM)(res1, res2, x, y, ctx); }
Expand Down
45 changes: 44 additions & 1 deletion src/gr/fmpz_mod.c
Original file line number Diff line number Diff line change
Expand Up @@ -334,11 +334,52 @@ _gr_fmpz_mod_div(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
status = _gr_fmpz_mod_inv(t, y, ctx);
if (status == GR_SUCCESS)
fmpz_mod_mul(res, x, t, FMPZ_MOD_CTX(ctx));

else
fmpz_zero(res);
fmpz_clear(t);

return status;
}

int
_gr_fmpz_mod_div_nonunique(fmpz_t res, const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
int status;

#if 1
status = fmpz_mod_divides(res, x, y, FMPZ_MOD_CTX(ctx)) ? GR_SUCCESS : GR_DOMAIN;
#else
if (FMPZ_MOD_IS_PRIME(ctx) != T_TRUE)
{
status = fmpz_mod_divides(res, x, y, FMPZ_MOD_CTX(ctx)) ? GR_SUCCESS : GR_DOMAIN;
}
else
{
fmpz_t t;
fmpz_init(t);
status = _gr_fmpz_mod_inv(t, y, ctx);
if (status == GR_SUCCESS)
fmpz_mod_mul(res, x, t, FMPZ_MOD_CTX(ctx));
else
fmpz_zero(res);
fmpz_clear(t);
}
#endif

return status;
}

truth_t
_gr_fmpz_mod_divides(const fmpz_t x, const fmpz_t y, const gr_ctx_t ctx)
{
truth_t res;
fmpz_t t;
fmpz_init(t);
res = fmpz_mod_divides(t, y, x, FMPZ_MOD_CTX(ctx)) ? T_TRUE : T_FALSE;
fmpz_clear(t);
return res;
}

truth_t
_gr_fmpz_mod_is_invertible(const fmpz_t x, const gr_ctx_t ctx)
{
Expand Down Expand Up @@ -669,6 +710,8 @@ gr_method_tab_input _fmpz_mod_methods_input[] =
{GR_METHOD_MUL_TWO, (gr_funcptr) _gr_fmpz_mod_mul_two},
{GR_METHOD_SQR, (gr_funcptr) _gr_fmpz_mod_sqr},
{GR_METHOD_DIV, (gr_funcptr) _gr_fmpz_mod_div},
{GR_METHOD_DIV_NONUNIQUE, (gr_funcptr) _gr_fmpz_mod_div_nonunique},
{GR_METHOD_DIVIDES, (gr_funcptr) _gr_fmpz_mod_divides},
{GR_METHOD_IS_INVERTIBLE, (gr_funcptr) _gr_fmpz_mod_is_invertible},
{GR_METHOD_INV, (gr_funcptr) _gr_fmpz_mod_inv},
{GR_METHOD_POW_UI, (gr_funcptr) _gr_fmpz_mod_pow_ui},
Expand Down
6 changes: 6 additions & 0 deletions src/gr/fmpz_poly.c
Original file line number Diff line number Diff line change
Expand Up @@ -554,6 +554,12 @@ _gr_fmpz_poly_divides(const fmpz_poly_t x, const fmpz_poly_t y, const gr_ctx_t c
truth_t res;
fmpz_poly_t tmp;

if (fmpz_poly_is_zero(y))
return T_TRUE;

if (fmpz_poly_is_zero(x))
return T_FALSE;

fmpz_poly_init(tmp);
res = fmpz_poly_divides(tmp, y, x) ? T_TRUE : T_FALSE;
fmpz_poly_clear(tmp);
Expand Down
7 changes: 7 additions & 0 deletions src/gr/fmpzi.c
Original file line number Diff line number Diff line change
Expand Up @@ -617,6 +617,13 @@ _gr_fmpzi_divides(const fmpzi_t x, const fmpzi_t y, const gr_ctx_t ctx)
{
fmpzi_t q, r;
truth_t result;

if (fmpzi_is_zero(y))
return T_TRUE;

if (fmpzi_is_zero(x))
return T_FALSE;

fmpzi_init(q);
fmpzi_init(r);

Expand Down
74 changes: 46 additions & 28 deletions src/gr/nmod.c
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include "gr_vec.h"
#include "gr_poly.h"
#include "gr_mat.h"
#include "ulong_extras.h"

#define NMOD_CTX_REF(ring_ctx) (((nmod_t *)((ring_ctx))))
#define NMOD_CTX(ring_ctx) (*NMOD_CTX_REF(ring_ctx))
Expand Down Expand Up @@ -367,55 +368,42 @@ _gr_nmod_sqr(ulong * res, const ulong * x, const gr_ctx_t ctx)
return _gr_nmod_mul(res, x, x, ctx);
}

/* optimize */
int
_gr_nmod_div(ulong * res, const ulong * x, const ulong * y, const gr_ctx_t ctx)
{
ulong t;
int status;

status = _gr_nmod_inv(&t, y, ctx);
_gr_nmod_mul(res, x, &t, ctx);

if (status == GR_SUCCESS)
_gr_nmod_mul(res, x, &t, ctx);

return status;
}

/* optimize */
int
_gr_nmod_div_si(ulong * res, const ulong * x, slong y, const gr_ctx_t ctx)
{
ulong t;
ulong c = nmod_set_si(y, NMOD_CTX(ctx));

y = nmod_set_si(y, NMOD_CTX(ctx));

if (n_gcdinv(&t, y, NMOD_CTX(ctx).n) == 1)
{
res[0] = nmod_mul(x[0], t, NMOD_CTX(ctx));
return GR_SUCCESS;
}
else
{
res[0] = 0;
return GR_DOMAIN;
}
return _gr_nmod_div(res, x, &c, ctx);
}

int
_gr_nmod_div_ui(ulong * res, const ulong * x, ulong y, const gr_ctx_t ctx)
{
ulong t;
ulong c = nmod_set_ui(y, NMOD_CTX(ctx));

y = nmod_set_ui(y, NMOD_CTX(ctx));
return _gr_nmod_div(res, x, &c, ctx);
}

if (n_gcdinv(&t, y, NMOD_CTX(ctx).n) == 1)
{
res[0] = nmod_mul(x[0], t, NMOD_CTX(ctx));
return GR_SUCCESS;
}
else
{
res[0] = 0;
return GR_DOMAIN;
}
int
_gr_nmod_div_fmpz(ulong * res, const ulong * x, const fmpz_t y, const gr_ctx_t ctx)
{
ulong c = fmpz_get_nmod(y, NMOD_CTX(ctx));

return _gr_nmod_div(res, x, &c, ctx);
}

truth_t
Expand All @@ -426,6 +414,33 @@ _gr_nmod_is_invertible(const ulong * x, const gr_ctx_t ctx)
return (g == 1) ? T_TRUE : T_FALSE;
}

truth_t
_gr_nmod_divides(const ulong * x, const ulong * y, const gr_ctx_t ctx)
{
ulong t;
return nmod_divides(&t, y[0], x[0], NMOD_CTX(ctx)) ? T_TRUE : T_FALSE;
}

int
_gr_nmod_div_nonunique(ulong * res, const ulong * x, const ulong * y, const gr_ctx_t ctx)
{
ulong t;
int status;

status = _gr_nmod_inv(&t, y, ctx);

if (status == GR_SUCCESS)
{
_gr_nmod_mul(res, x, &t, ctx);
}
else
{
status = nmod_divides(res, *x, *y, NMOD_CTX(ctx)) ? GR_SUCCESS : GR_DOMAIN;
}

return status;
}

int
_gr_nmod_mul_2exp_si(ulong * res, ulong * x, slong y, const gr_ctx_t ctx)
{
Expand Down Expand Up @@ -1390,6 +1405,9 @@ gr_method_tab_input __gr_nmod_methods_input[] =
{GR_METHOD_DIV, (gr_funcptr) _gr_nmod_div},
{GR_METHOD_DIV_SI, (gr_funcptr) _gr_nmod_div_si},
{GR_METHOD_DIV_UI, (gr_funcptr) _gr_nmod_div_ui},
{GR_METHOD_DIV_FMPZ, (gr_funcptr) _gr_nmod_div_fmpz},
{GR_METHOD_DIV_NONUNIQUE, (gr_funcptr) _gr_nmod_div_nonunique},
{GR_METHOD_DIVIDES, (gr_funcptr) _gr_nmod_divides},
{GR_METHOD_IS_INVERTIBLE, (gr_funcptr) _gr_nmod_is_invertible},
{GR_METHOD_INV, (gr_funcptr) _gr_nmod_inv},
{GR_METHOD_SQRT, (gr_funcptr) _gr_nmod_sqrt},
Expand Down
Loading

0 comments on commit 8548d83

Please sign in to comment.