Skip to content

Commit

Permalink
Remove functions requiring strong order if weak order
Browse files Browse the repository at this point in the history
In order to not lie to the user.

Also document this and define preprocessor macros in case they are
defined.

Involved functions are:

- nmod_mpoly_divides_heap_threaded
- _nmod_mpoly_divides_heap_threaded_pool
- nmod_mpolyn_divides_threaded_pool

- fmpz_mpoly_divides_heap_threaded
- _fmpz_mpoly_divides_heap_threaded_pool
  • Loading branch information
albinahlback committed Dec 28, 2023
1 parent 5dddefd commit 1fd6de1
Show file tree
Hide file tree
Showing 15 changed files with 151 additions and 47 deletions.
6 changes: 3 additions & 3 deletions configure.ac
Original file line number Diff line number Diff line change
Expand Up @@ -646,17 +646,17 @@ AC_MSG_CHECKING([if memory is strongly-ordered])
case "$host_cpu" in
x86_64|x86|i386|i586)
flint_know_strong_order="yes"
AC_MSG_CHECKING([yes])
AC_MSG_RESULT([yes])
;;
*)
flint_know_strong_order="no"
AC_MSG_CHECKING([unsure])
AC_MSG_RESULT([unsure])
;;
esac

if test "$flint_know_strong_order" = "yes";
then
AC_DEFINE(FLINT_KNOW_STRONG_ORDER,1,[Define if system is strongly ordered])]
AC_DEFINE(FLINT_KNOW_STRONG_ORDER,1,[Define if system is strongly ordered])
fi

################################################################################
Expand Down
25 changes: 19 additions & 6 deletions doc/source/fmpz_mpoly.rst
Original file line number Diff line number Diff line change
Expand Up @@ -777,15 +777,28 @@ Internal Functions

.. function:: int fmpz_mpoly_divides_monagan_pearce(fmpz_mpoly_t poly1, const fmpz_mpoly_t poly2, const fmpz_mpoly_t poly3, const fmpz_mpoly_ctx_t ctx)

Set ``poly1`` to ``poly2`` divided by ``poly3`` and return 1 if the quotient
is exact. Otherwise return 0. The function uses the algorithm of Michael
Monagan and Roman Pearce. Note that the function
``fmpz_mpoly_div_monagan_pearce`` below may be much faster if the quotient
is known to be exact.

.. function:: int fmpz_mpoly_divides_heap_threaded(fmpz_mpoly_t Q, const fmpz_mpoly_t A, const fmpz_mpoly_t B, const fmpz_mpoly_ctx_t ctx)

Set ``poly1`` to ``poly2`` divided by ``poly3`` and return 1 if
the quotient is exact. Otherwise return 0. The function uses the algorithm
of Michael Monagan and Roman Pearce. Note that the function
``fmpz_mpoly_div_monagan_pearce`` below may be much faster if the
quotient is known to be exact.
The same method as used as in :func:``fmpz_mpoly_divides_monagan_pearce``,
but is also multi-threaded.

.. note::

This function is only defined if the machine is known to be strongly ordered
during the configuration. To check whether this function is defined during
compilation-time, use the C preprocessor macro
``#ifdef fmpz_mpoly_divides_heap_threaded``.

The threaded version takes an upper limit on the number of threads to use, while the first version always uses one thread.
Note that, if the system is known to be strongly ordered, the underlying
algorithm for this function is utilized in :func:``fmpz_mpoly_divides``.
Hence, you may find it easier to use this function instead if the C
preprocessor is not available.

.. function:: slong _fmpz_mpoly_div_monagan_pearce(fmpz ** polyq, ulong ** expq, slong * allocq, const fmpz * poly2, const ulong * exp2, slong len2, const fmpz * poly3, const ulong * exp3, slong len3, slong bits, slong N, const mp_limb_t * cmpmask)

Expand Down
17 changes: 15 additions & 2 deletions doc/source/nmod_mpoly.rst
Original file line number Diff line number Diff line change
Expand Up @@ -518,8 +518,21 @@ The division functions assume that the modulus is prime.

.. function:: int nmod_mpoly_divides_heap_threaded(nmod_mpoly_t Q, const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx)

Do the operation of ``nmod_mpoly_divides`` using a heap and multiple threads.
This function should only be called once ``global_thread_pool`` has been initialized.
Do the operation of ``nmod_mpoly_divides`` using the heap and multiple
threads. This function should only be called once ``global_thread_pool`` has
been initialized.

.. note::

This function is only defined if the machine is known to be strongly ordered
during the configuration. To check whether this function is defined during
compilation-time, use the C preprocessor macro
``#ifdef nmod_mpoly_divides_heap_threaded``.

Note that, if the system is known to be strongly ordered, the underlying
algorithm for this function is utilized in :func:``nmod_mpoly_divides``.
Hence, you may find it easier to use this function instead if the C
preprocessor is not available.


Greatest Common Divisor
Expand Down
4 changes: 4 additions & 0 deletions src/fmpz_mpoly.h
Original file line number Diff line number Diff line change
Expand Up @@ -770,12 +770,16 @@ int fmpz_mpoly_divides(fmpz_mpoly_t Q,
int fmpz_mpoly_divides_monagan_pearce(fmpz_mpoly_t Q,
const fmpz_mpoly_t A, const fmpz_mpoly_t B, const fmpz_mpoly_ctx_t ctx);

#if FLINT_KNOW_STRONG_ORDER
#define fmpz_mpoly_divides_heap_threaded fmpz_mpoly_divides_heap_threaded
int fmpz_mpoly_divides_heap_threaded(fmpz_mpoly_t Q,
const fmpz_mpoly_t A, const fmpz_mpoly_t B, const fmpz_mpoly_ctx_t ctx);

#define _fmpz_mpoly_divides_heap_threaded_pool _fmpz_mpoly_divides_heap_threaded_pool
int _fmpz_mpoly_divides_heap_threaded_pool(fmpz_mpoly_t Q,
const fmpz_mpoly_t A, const fmpz_mpoly_t B, const fmpz_mpoly_ctx_t ctx,
const thread_pool_handle * handles, slong num_handles);
#endif

slong _fmpz_mpoly_divides_array(fmpz ** poly1, ulong ** exp1,
slong * alloc, const fmpz * poly2, const ulong * exp2, slong len2,
Expand Down
25 changes: 24 additions & 1 deletion src/fmpz_mpoly/divides.c
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,12 @@
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/

#include "thread_support.h"
#include "fmpz_mpoly.h"

#ifdef _fmpz_mpoly_divides_heap_threaded_pool

#include "thread_support.h"

int fmpz_mpoly_divides(
fmpz_mpoly_t Q,
const fmpz_mpoly_t A,
Expand Down Expand Up @@ -52,4 +55,24 @@ int fmpz_mpoly_divides(

return divides;
}
#else
int fmpz_mpoly_divides(
fmpz_mpoly_t Q,
const fmpz_mpoly_t A,
const fmpz_mpoly_t B,
const fmpz_mpoly_ctx_t ctx)
{
if (B->length == 0)
{
flint_throw(FLINT_DIVZERO, "Divide by zero in fmpz_mpoly_divides");
}

if (A->length == 0)
{
fmpz_mpoly_zero(Q, ctx);
return 1;
}

return fmpz_mpoly_divides_monagan_pearce(Q, A, B, ctx);
}
#endif
42 changes: 22 additions & 20 deletions src/fmpz_mpoly/divides_heap_threaded.c
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,12 @@
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/

#include "fmpz_mpoly.h"

#if FLINT_KNOW_STRONG_ORDER

#include "thread_support.h"
#include "ulong_extras.h"
#include "fmpz_mpoly.h"

/*
a thread safe mpoly supports three mutating operations
Expand All @@ -34,9 +37,9 @@ typedef struct _fmpz_mpoly_ts_struct
typedef fmpz_mpoly_ts_struct fmpz_mpoly_ts_t[1];

/* Bcoeff is changed */
void fmpz_mpoly_ts_init(fmpz_mpoly_ts_t A,
fmpz * Bcoeff, ulong * Bexp, slong Blen,
flint_bitcnt_t bits, slong N)
static void fmpz_mpoly_ts_init(fmpz_mpoly_ts_t A,
fmpz * Bcoeff, ulong * Bexp, slong Blen,
flint_bitcnt_t bits, slong N)
{
slong i;
flint_bitcnt_t idx = FLINT_BIT_COUNT(Blen);
Expand All @@ -61,7 +64,7 @@ void fmpz_mpoly_ts_init(fmpz_mpoly_ts_t A,
}
}

void fmpz_mpoly_ts_print(const fmpz_mpoly_ts_t B, const char ** x,
static void fmpz_mpoly_ts_print(const fmpz_mpoly_ts_t B, const char ** x,
const fmpz_mpoly_ctx_t ctx)
{
fmpz_mpoly_t A;
Expand All @@ -75,7 +78,7 @@ void fmpz_mpoly_ts_print(const fmpz_mpoly_ts_t B, const char ** x,
fmpz_mpoly_assert_canonical(A, ctx);
}

void fmpz_mpoly_ts_clear(fmpz_mpoly_ts_t A)
static void fmpz_mpoly_ts_clear(fmpz_mpoly_ts_t A)
{
slong i;

Expand All @@ -95,7 +98,7 @@ void fmpz_mpoly_ts_clear(fmpz_mpoly_ts_t A)
}
}

void fmpz_mpoly_ts_clear_poly(fmpz_mpoly_t Q, fmpz_mpoly_ts_t A)
static void fmpz_mpoly_ts_clear_poly(fmpz_mpoly_t Q, fmpz_mpoly_ts_t A)
{
if (Q->alloc != 0)
{
Expand Down Expand Up @@ -126,7 +129,7 @@ void fmpz_mpoly_ts_clear_poly(fmpz_mpoly_t Q, fmpz_mpoly_ts_t A)


/* put B on the end of A - Bcoeff is changed*/
void fmpz_mpoly_ts_append(fmpz_mpoly_ts_t A,
static void fmpz_mpoly_ts_append(fmpz_mpoly_ts_t A,
fmpz * Bcoeff, ulong * Bexps, slong Blen, slong N)
{
/* TODO: this needs barriers on non-x86 */
Expand Down Expand Up @@ -330,7 +333,7 @@ static void divides_heap_base_add_chunk(divides_heap_base_t H, divides_heap_chun
saveD: 0 means we can modify coeffs of input D
1 means we must not modify coeffs of input D
*/
slong _fmpz_mpoly_mulsub_stripe1(fmpz ** A_coeff, ulong ** A_exp, slong * A_alloc,
static slong _fmpz_mpoly_mulsub_stripe1(fmpz ** A_coeff, ulong ** A_exp, slong * A_alloc,
const fmpz * Dcoeff, const ulong * Dexp, slong Dlen, int saveD,
const fmpz * Bcoeff, const ulong * Bexp, slong Blen,
const fmpz * Ccoeff, const ulong * Cexp, slong Clen,
Expand Down Expand Up @@ -606,7 +609,7 @@ slong _fmpz_mpoly_mulsub_stripe1(fmpz ** A_coeff, ulong ** A_exp, slong * A_allo
}


slong _fmpz_mpoly_mulsub_stripe(fmpz ** A_coeff, ulong ** A_exp, slong * A_alloc,
static slong _fmpz_mpoly_mulsub_stripe(fmpz ** A_coeff, ulong ** A_exp, slong * A_alloc,
const fmpz * Dcoeff, const ulong * Dexp, slong Dlen, int saveD,
const fmpz * Bcoeff, const ulong * Bexp, slong Blen,
const fmpz * Ccoeff, const ulong * Cexp, slong Clen,
Expand Down Expand Up @@ -911,7 +914,7 @@ slong _fmpz_mpoly_mulsub_stripe(fmpz ** A_coeff, ulong ** A_exp, slong * A_alloc
Q = stripe of A/B (assume A != 0)
return Qlen = 0 if exact division is impossible
*/
slong _fmpz_mpoly_divides_stripe1(
static slong _fmpz_mpoly_divides_stripe1(
fmpz ** Q_coeff, ulong ** Q_exp, slong * Q_alloc,
const fmpz * Acoeff, const ulong * Aexp, slong Alen,
const fmpz * Bcoeff, const ulong * Bexp, slong Blen,
Expand Down Expand Up @@ -1248,7 +1251,7 @@ slong _fmpz_mpoly_divides_stripe1(
goto cleanup;
}

slong _fmpz_mpoly_divides_stripe(
static slong _fmpz_mpoly_divides_stripe(
fmpz ** Q_coeff, ulong ** Q_exp, slong * Q_alloc,
const fmpz * Acoeff, const ulong * Aexp, slong Alen,
const fmpz * Bcoeff, const ulong * Bexp, slong Blen,
Expand Down Expand Up @@ -1945,21 +1948,21 @@ static void worker_loop(void * varg)
#if FLINT_USES_PTHREAD
pthread_mutex_lock(&H->mutex);
#endif
if (L->lock != -1)
if (L->lock != -1)
{
L->lock = -1;
#if FLINT_USES_PTHREAD
pthread_mutex_unlock(&H->mutex);
#endif
trychunk(W, L);
trychunk(W, L);
#if FLINT_USES_PTHREAD
pthread_mutex_lock(&H->mutex);
#endif
L->lock = 0;
L->lock = 0;
#if FLINT_USES_PTHREAD
pthread_mutex_unlock(&H->mutex);
#endif
break;
break;
}
else
{
Expand Down Expand Up @@ -2004,10 +2007,6 @@ int _fmpz_mpoly_divides_heap_threaded_pool(
divides_heap_base_t H;
TMP_INIT;

#if !FLINT_KNOW_STRONG_ORDER
return fmpz_mpoly_divides_monagan_pearce(Q, A, B, ctx);
#endif

if (B->length < 2 || A->length < 2)
{
return fmpz_mpoly_divides_monagan_pearce(Q, A, B, ctx);
Expand Down Expand Up @@ -2245,3 +2244,6 @@ int fmpz_mpoly_divides_heap_threaded(

return divides;
}
#else
typedef int this_file_is_empty;
#endif
7 changes: 7 additions & 0 deletions src/fmpz_mpoly/test/t-divides_heap_threaded.c
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include "test_helpers.h"
#include "fmpz_mpoly.h"

#if defined(fmpz_mpoly_divides_heap_threaded)
TEST_FUNCTION_START(fmpz_mpoly_divides_heap_threaded, state)
{
int result, result2;
Expand Down Expand Up @@ -417,3 +418,9 @@ TEST_FUNCTION_START(fmpz_mpoly_divides_heap_threaded, state)

TEST_FUNCTION_END(state);
}
#else
TEST_FUNCTION_START(fmpz_mpoly_divides_heap_threaded, state)
{
TEST_FUNCTION_END_SKIPPED(state);
}
#endif
7 changes: 7 additions & 0 deletions src/nmod_mpoly.h
Original file line number Diff line number Diff line change
Expand Up @@ -917,12 +917,16 @@ int _nmod_mpoly_divides_threaded_pool(nmod_mpoly_t Q,
int nmod_mpoly_divides_monagan_pearce(nmod_mpoly_t Q,
const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx);

#if FLINT_KNOW_STRONG_ORDER
#define nmod_mpoly_divides_heap_threaded nmod_mpoly_divides_heap_threaded
int nmod_mpoly_divides_heap_threaded(nmod_mpoly_t Q,
const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx);

#define _nmod_mpoly_divides_heap_threaded_pool _nmod_mpoly_divides_heap_threaded_pool
int _nmod_mpoly_divides_heap_threaded_pool(nmod_mpoly_t Q,
const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx,
const thread_pool_handle * handles, slong num_handles);
#endif

int nmod_mpoly_divides_dense(nmod_mpoly_t Q,
const nmod_mpoly_t A, const nmod_mpoly_t B, const nmod_mpoly_ctx_t ctx);
Expand Down Expand Up @@ -1556,9 +1560,12 @@ void nmod_mpolyun_divexact_last(nmod_mpolyun_t A, n_poly_t b,
int nmod_mpolyn_divides(nmod_mpolyn_t Q, const nmod_mpolyn_t A,
const nmod_mpolyn_t B, const nmod_mpoly_ctx_t ctx);

#if FLINT_KNOW_STRONG_ORDER
#define nmod_mpolyn_divides_threaded_pool nmod_mpolyn_divides_threaded_pool
int nmod_mpolyn_divides_threaded_pool(nmod_mpolyn_t Q,
const nmod_mpolyn_t A, const nmod_mpolyn_t B, const nmod_mpoly_ctx_t ctx,
const thread_pool_handle * handles, slong num_handles);
#endif

int nmod_mpolyun_divides(nmod_mpolyun_t Q, const nmod_mpolyun_t A,
const nmod_mpolyun_t B, const nmod_mpoly_ctx_t ctx);
Expand Down
4 changes: 4 additions & 0 deletions src/nmod_mpoly/divides.c
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,7 @@ int _nmod_mpoly_divides_threaded_pool(
goto cleanup;
}

#ifdef _nmod_mpoly_divides_heap_threaded_pool
if (num_handles > 0)
{
divides = _nmod_mpoly_divides_heap_threaded_pool(Q, A, B, ctx,
Expand All @@ -142,6 +143,9 @@ int _nmod_mpoly_divides_threaded_pool(
{
divides = nmod_mpoly_divides_monagan_pearce(Q, A, B, ctx);
}
#else
divides = nmod_mpoly_divides_monagan_pearce(Q, A, B, ctx);
#endif

cleanup:

Expand Down
12 changes: 7 additions & 5 deletions src/nmod_mpoly/divides_heap_threaded.c
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,11 @@
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/

#include "thread_support.h"
#include "nmod_mpoly.h"

#if FLINT_KNOW_STRONG_ORDER

#include "thread_support.h"
#include "fmpz_mpoly.h"

#define PROFILE_THIS 0
Expand Down Expand Up @@ -1710,10 +1713,6 @@ int _nmod_mpoly_divides_heap_threaded_pool(
#endif
TMP_INIT;

#if !FLINT_KNOW_STRONG_ORDER
return nmod_mpoly_divides_monagan_pearce(Q, A, B, ctx);
#endif

if (B->length < 2 || A->length < 2)
{
return nmod_mpoly_divides_monagan_pearce(Q, A, B, ctx);
Expand Down Expand Up @@ -1967,3 +1966,6 @@ int nmod_mpoly_divides_heap_threaded(

return divides;
}
#else
typedef int this_file_is_empty;
#endif
Loading

0 comments on commit 1fd6de1

Please sign in to comment.