Skip to content

Instantly share code, notes, and snippets.

@fredrik-johansson

fredrik-johansson/gr.c

Last active November 26, 2021 16:18
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save fredrik-johansson/4c58e476f9ffa4c5a7ce218c52df0002 to your computer and use it in GitHub Desktop.
Save fredrik-johansson/4c58e476f9ffa4c5a7ce218c52df0002 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gr.c
/*
Generic Flint-style rings using void pointers + context objects.
todo: write -> print
Principles/goals/benefits:
* Small code size, fast compilation.
* Possible to pack data efficiently (down to 1 byte / element).
* Plain C, similar interface to existing Flint code.
* Data layouts backwards compatible with most existing Flint types.
* Support all unusual cases in Flint/Arb/Calcium uniformly (error handling,
inexact rings, noncomputable rings, context objects), with a unified
interface.
* Support fast stack based allocation of local variables and vectors.
* Fully in-place operations.
* Allow fast (a few cycles) runtime construction of context objects.
* Possibility to use generic default methods or provide optimized versions
(e.g. for vector operations).
Cons:
* Several cycles overhead for each dynamic function lookup and function call.
On top of this, most compiler optimizations become impossible.
However, overriding the default vector functions should largely make up
for overheads.
* Less compiler protection against type errors.
* Some function call signatures need to change in ways that make
the interface less convenient (for uniformity).
*/
#include <string.h>
#include "flint/flint.h"
#include "flint/nmod_vec.h"
#include "flint/ulong_extras.h"
#include "flint/profiler.h"
/* Copied from Calcium: stream interface allows simple file or string IO. */
typedef struct
{
FILE * fp;
char * s;
slong len;
slong alloc;
}
gr_stream_struct;
typedef gr_stream_struct gr_stream_t[1];
void gr_stream_init_file(gr_stream_t out, FILE * fp)
{
out->fp = fp;
out->s = NULL;
}
void gr_stream_init_str(gr_stream_t out)
{
out->fp = NULL;
out->s = flint_malloc(16);
out->s[0] = '\0';
out->len = 0;
out->alloc = 16;
}
void gr_stream_write(gr_stream_t out, const char * s)
{
if (out->fp != NULL)
{
fprintf(out->fp, "%s", s);
}
else
{
slong len, alloc;
len = strlen(s);
alloc = out->len + len + 1;
if (alloc > out->alloc)
{
alloc = FLINT_MAX(alloc, out->alloc * 2);
out->s = realloc(out->s, alloc);
out->alloc = alloc;
}
memcpy(out->s + out->len, s, len + 1);
out->len += len;
}
}
void
gr_stream_write_si(gr_stream_t out, slong x)
{
if (out->fp != NULL)
{
flint_fprintf(out->fp, "%wd", x);
}
else
{
char tmp[22];
if (sizeof(slong) == sizeof(long))
sprintf(tmp, "%ld", x);
else
flint_sprintf(tmp, "%wd", x);
gr_stream_write(out, tmp);
}
}
void
gr_stream_write_free(gr_stream_t out, char * s)
{
gr_stream_write(out, s);
flint_free(s);
}
/*
All functions implement error handling and return a status code.
GR_SUCCESS - The operation finished as expected.
GR_DOMAIN - Invalid input for this operation (e.g. division by zero,
wrong matrix shape).
GR_UNABLE - The operation could not be performed for implementation reasons
(e.g. overflow, missing algorithm). The input may or may not be
valid.
GR_WRONG - Test failure (used only in test code).
For uniformity, all functions return a status code.
Flags can be OR'ed to avoid complex control flow.
*/
#define GR_SUCCESS 0
#define GR_DOMAIN 1
#define GR_UNABLE 2
#define GR_WRONG 4
typedef void * gr_ptr;
typedef const void * gr_srcptr;
typedef void * gr_ctx_ptr;
typedef int ((*gr_func_out)(gr_ptr, gr_ctx_ptr));
typedef int ((*gr_func_out_si)(gr_ptr, slong, gr_ctx_ptr));
typedef int ((*gr_func_out_out)(gr_ptr, gr_ptr, gr_ctx_ptr));
typedef int ((*gr_func_out_out_si)(gr_ptr, gr_ptr, slong, gr_ctx_ptr));
typedef int ((*gr_func_out_in)(gr_ptr, gr_srcptr, gr_ctx_ptr));
typedef int ((*gr_func_out_in_in)(gr_ptr, gr_srcptr, gr_srcptr, gr_ctx_ptr));
typedef int ((*gr_func_out_in_si)(gr_ptr, gr_srcptr, slong, gr_ctx_ptr));
typedef int ((*gr_func_out_si_in)(gr_ptr, gr_srcptr, slong, gr_srcptr, gr_ctx_ptr));
typedef int ((*gr_func_out_in_in_si)(gr_ptr, gr_srcptr, gr_srcptr, slong, gr_ctx_ptr));
typedef int ((*gr_func_out_in_si_in)(gr_ptr, gr_srcptr, slong, gr_srcptr, gr_ctx_ptr));
typedef int ((*gr_func_out_in_si_si)(gr_ptr, gr_srcptr, slong, slong, gr_ctx_ptr));
typedef int ((*gr_func_outbool_in)(int *, gr_srcptr, gr_ctx_ptr));
typedef int ((*gr_func_outbool_in_si)(int *, gr_srcptr, slong, gr_ctx_ptr));
typedef int ((*gr_func_outbool_in_in)(int *, gr_srcptr, gr_srcptr, gr_ctx_ptr));
typedef int ((*gr_func_outbool_in_in_si)(int *, gr_srcptr, gr_srcptr, slong, gr_ctx_ptr));
typedef int ((*gr_func_out_in_bool_in_in_si)(gr_ptr, gr_srcptr, int, gr_srcptr, gr_srcptr, slong, gr_ctx_ptr));
typedef int ((*gr_func_stream_in)(gr_stream_t, gr_srcptr, gr_ctx_ptr));
typedef int ((*gr_func_randtest)(gr_ptr, flint_rand_t state, const void * options, gr_ctx_ptr));
int gr_not_implemented(void) { return GR_UNABLE; }
/* Standard methods which need to be implemented for any ring. */
/* What else is needed: fmpz, si methods? Random generation? Hashing? */
typedef struct
{
gr_func_out init;
gr_func_out clear;
gr_func_out_out swap;
gr_func_randtest randtest;
gr_func_stream_in write;
gr_func_out zero;
gr_func_out one;
gr_func_outbool_in is_zero;
gr_func_outbool_in is_one;
gr_func_outbool_in_in equal;
gr_func_out_in set;
gr_func_out_si set_si;
gr_func_out_in neg;
gr_func_out_in_in add;
gr_func_out_in_si add_si;
gr_func_out_in_in sub;
gr_func_out_in_in mul;
gr_func_out_in_si mul_si;
gr_func_out_in_in div;
gr_func_out_in_in divexact;
gr_func_out_in_si divexact_si;
gr_func_outbool_in is_invertible;
gr_func_out_in inv;
}
gr_methods;
/* Vector methods. The idea is that these will have generic defaults
which can be overridden for performance. */
typedef struct
{
gr_func_out_si init;
gr_func_out_si clear;
gr_func_out_out_si swap;
gr_func_out_si zero;
gr_func_out_in_si set;
gr_func_out_in_si neg;
gr_func_out_in_in_si add;
gr_func_out_in_in_si sub;
gr_func_out_in_si_in scalar_addmul;
gr_func_out_in_si_in scalar_submul;
gr_func_out_in_si_si scalar_addmul_si;
gr_func_out_in_si_si scalar_submul_si;
gr_func_outbool_in_in_si equal;
gr_func_outbool_in_si is_zero;
gr_func_out_in_bool_in_in_si dot;
gr_func_out_in_bool_in_in_si dot_rev;
}
gr_vec_methods;
/*
We could add more method tables, which may be initialized to generic
implementations by default:
compound_methods
cmp_methods
...
*/
/* Flags. These are not cumulative. (?) */
#define GR_COMMUTATIVE_RING UWORD(1)
#define GR_INTEGRAL_DOMAIN UWORD(2)
#define GR_FIELD UWORD(4)
typedef struct
{
ulong flags;
ssize_t sizeof_elem;
void * elem_ctx;
gr_methods * methods;
gr_vec_methods * vec_methods;
char * debug_string;
}
gr_ctx_struct;
typedef gr_ctx_struct gr_ctx_t[1];
/* Macros for allocating temporary variables on the stack. */
#define GR_TMP_START TMP_INIT; TMP_START;
#define GR_TMP_ALLOC TMP_ALLOC
#define GR_TMP_END TMP_END
/* todo: use vector init functions when provided */
#define GR_TMP_INIT_VEC(vec, len, ctx) \
do { \
gr_func_out _gr_init_func = (ctx)->methods->init; \
ssize_t _gr_elem_size = (ctx)->sizeof_elem; \
(vec) = (gr_ptr) GR_TMP_ALLOC((len) * _gr_elem_size); \
_gr_vec_init((vec), (len), (ctx)); \
} while (0)
#define GR_TMP_CLEAR_VEC(vec, len, ctx) \
do { \
gr_func_out _gr_init_func = (ctx)->methods->init; \
ssize_t _gr_elem_size = (ctx)->sizeof_elem; \
_gr_vec_clear((vec), (len), (ctx)); \
} while (0)
#define GR_TMP_INIT1(x1, ctx) \
do { \
gr_func_out _gr_init_func = (ctx)->methods->init; \
ssize_t _gr_elem_size = (ctx)->sizeof_elem; \
x1 = (gr_ptr) GR_TMP_ALLOC(1 * _gr_elem_size); \
_gr_init_func(x1, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_INIT2(x1, x2, ctx) \
do { \
gr_func_out _gr_init_func = (ctx)->methods->init; \
ssize_t _gr_elem_size = (ctx)->sizeof_elem; \
x1 = (gr_ptr) GR_TMP_ALLOC(2 * _gr_elem_size); \
x2 = (gr_ptr) ((char *) x1 + _gr_elem_size); \
_gr_init_func(x1, (ctx)->elem_ctx); \
_gr_init_func(x2, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_INIT3(x1, x2, x3, ctx) \
do { \
gr_func_out _gr_init_func = (ctx)->methods->init; \
ssize_t _gr_elem_size = (ctx)->sizeof_elem; \
x1 = (gr_ptr) GR_TMP_ALLOC(3 * _gr_elem_size); \
x2 = (gr_ptr) ((char *) x1 + _gr_elem_size); \
x3 = (gr_ptr) ((char *) x2 + _gr_elem_size); \
_gr_init_func(x1, (ctx)->elem_ctx); \
_gr_init_func(x2, (ctx)->elem_ctx); \
_gr_init_func(x3, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_INIT4(x1, x2, x3, x4, ctx) \
do { \
gr_func_out _gr_init_func = (ctx)->methods->init; \
ssize_t _gr_elem_size = (ctx)->sizeof_elem; \
x1 = (gr_ptr) GR_TMP_ALLOC(4 * _gr_elem_size); \
x2 = (gr_ptr) ((char *) x1 + _gr_elem_size); \
x3 = (gr_ptr) ((char *) x2 + _gr_elem_size); \
x4 = (gr_ptr) ((char *) x3 + _gr_elem_size); \
_gr_init_func(x1, (ctx)->elem_ctx); \
_gr_init_func(x2, (ctx)->elem_ctx); \
_gr_init_func(x3, (ctx)->elem_ctx); \
_gr_init_func(x4, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_INIT5(x1, x2, x3, x4, x5, ctx) \
do { \
gr_func_out _gr_init_func = (ctx)->methods->init; \
ssize_t _gr_elem_size = (ctx)->sizeof_elem; \
x1 = (gr_ptr) GR_TMP_ALLOC(5 * _gr_elem_size); \
x2 = (gr_ptr) ((char *) x1 + _gr_elem_size); \
x3 = (gr_ptr) ((char *) x2 + _gr_elem_size); \
x4 = (gr_ptr) ((char *) x3 + _gr_elem_size); \
x5 = (gr_ptr) ((char *) x4 + _gr_elem_size); \
_gr_init_func(x1, (ctx)->elem_ctx); \
_gr_init_func(x2, (ctx)->elem_ctx); \
_gr_init_func(x3, (ctx)->elem_ctx); \
_gr_init_func(x4, (ctx)->elem_ctx); \
_gr_init_func(x5, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_CLEAR1(x1, ctx) \
do { \
gr_func_out _gr_clear_func = (ctx)->methods->clear; \
_gr_clear_func(x1, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_CLEAR2(x1, x2, ctx) \
do { \
gr_func_out _gr_clear_func = (ctx)->methods->clear; \
_gr_clear_func(x1, (ctx)->elem_ctx); \
_gr_clear_func(x2, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_CLEAR3(x1, x2, x3, ctx) \
do { \
gr_func_out _gr_clear_func = (ctx)->methods->clear; \
_gr_clear_func(x1, (ctx)->elem_ctx); \
_gr_clear_func(x2, (ctx)->elem_ctx); \
_gr_clear_func(x3, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_CLEAR4(x1, x2, x3, x4, ctx) \
do { \
gr_func_out _gr_clear_func = (ctx)->methods->clear; \
_gr_clear_func(x1, (ctx)->elem_ctx); \
_gr_clear_func(x2, (ctx)->elem_ctx); \
_gr_clear_func(x3, (ctx)->elem_ctx); \
_gr_clear_func(x4, (ctx)->elem_ctx); \
} while (0)
#define GR_TMP_CLEAR5(x1, x2, x3, x4, x5, ctx) \
do { \
gr_func_out _gr_clear_func = (ctx)->methods->clear; \
_gr_clear_func(x1, (ctx)->elem_ctx); \
_gr_clear_func(x2, (ctx)->elem_ctx); \
_gr_clear_func(x3, (ctx)->elem_ctx); \
_gr_clear_func(x4, (ctx)->elem_ctx); \
_gr_clear_func(x5, (ctx)->elem_ctx); \
} while (0)
/* Inline wrappers for default methods. */
int
gr_init(gr_ptr res, gr_ctx_t ctx)
{
return ctx->methods->init(res, ctx->elem_ctx);
}
int
gr_clear(gr_ptr res, gr_ctx_t ctx)
{
return ctx->methods->clear(res, ctx->elem_ctx);
}
int
gr_swap(gr_ptr x, gr_ptr y, gr_ctx_t ctx)
{
return ctx->methods->swap(x, y, ctx->elem_ctx);
}
int
gr_randtest(gr_ptr x, flint_rand_t state, const void * options, gr_ctx_t ctx)
{
return ctx->methods->randtest(x, state, options, ctx->elem_ctx);
}
int
gr_write(gr_stream_t out, gr_srcptr x, gr_ctx_t ctx)
{
return ctx->methods->write(out, x, ctx->elem_ctx);
}
int
gr_zero(gr_ptr res, gr_ctx_t ctx)
{
return ctx->methods->zero(res, ctx->elem_ctx);
}
int
gr_one(gr_ptr res, gr_ctx_t ctx)
{
return ctx->methods->one(res, ctx->elem_ctx);
}
int
gr_set_si(gr_ptr res, slong x, gr_ctx_t ctx)
{
return ctx->methods->set_si(res, x, ctx->elem_ctx);
}
int
gr_is_zero(int * res, gr_srcptr x, gr_ctx_t ctx)
{
return ctx->methods->is_zero(res, x, ctx->elem_ctx);
}
int
gr_is_one(int * res, gr_srcptr x, gr_ctx_t ctx)
{
return ctx->methods->is_one(res, x, ctx->elem_ctx);
}
int
gr_equal(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
{
return ctx->methods->equal(res, x, y, ctx->elem_ctx);
}
int
gr_set(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
{
return ctx->methods->set(res, x, ctx->elem_ctx);
}
int
gr_neg(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
{
return ctx->methods->neg(res, x, ctx->elem_ctx);
}
int
gr_add(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
{
return ctx->methods->add(res, x, y, ctx->elem_ctx);
}
int
gr_add_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
{
return ctx->methods->add_si(res, x, y, ctx->elem_ctx);
}
int
gr_sub(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
{
return ctx->methods->sub(res, x, y, ctx->elem_ctx);
}
int
gr_mul(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
{
return ctx->methods->mul(res, x, y, ctx->elem_ctx);
}
int
gr_mul_si(gr_ptr res, gr_srcptr x, slong y, gr_ctx_t ctx)
{
return ctx->methods->mul_si(res, x, y, ctx->elem_ctx);
}
int
gr_div(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
{
return ctx->methods->div(res, x, y, ctx->elem_ctx);
}
int
gr_inv(gr_ptr res, gr_srcptr x, gr_ctx_t ctx)
{
return ctx->methods->inv(res, x, ctx->elem_ctx);
}
int
gr_is_invertible(int * res, gr_srcptr x, gr_ctx_t ctx)
{
return ctx->methods->is_invertible(res, x, ctx->elem_ctx);
}
/* Generic vector functions */
#define GR_ENTRY(vec, i, size) ((void *) (((char *) (vec)) + ((i) * (size))))
int
gr_generic_vec_init(gr_ptr vec, slong len, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
status |= gr_init(GR_ENTRY(vec, i, sz), ctx);
return status;
}
int
gr_generic_vec_clear(gr_ptr vec, slong len, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
status |= gr_clear(GR_ENTRY(vec, i, sz), ctx);
return status;
}
int
gr_generic_vec_swap(gr_ptr vec1, gr_ptr vec2, slong len, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
status |= gr_swap(GR_ENTRY(vec1, i, sz), GR_ENTRY(vec2, i, sz), ctx);
return status;
}
int
gr_generic_vec_zero(gr_ptr vec, slong len, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
status |= gr_zero(GR_ENTRY(vec, i, sz), ctx);
return status;
}
int
gr_generic_vec_set(gr_ptr res, gr_srcptr src, slong len, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
status |= gr_set(GR_ENTRY(res, i, sz), GR_ENTRY(src, i, sz), ctx);
return status;
}
int
gr_generic_vec_neg(gr_ptr res, gr_srcptr src, slong len, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
status |= gr_neg(GR_ENTRY(res, i, sz), GR_ENTRY(src, i, sz), ctx);
return status;
}
int
gr_generic_vec_add(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
status |= gr_add(GR_ENTRY(res, i, sz), GR_ENTRY(src1, i, sz), GR_ENTRY(src2, i, sz), ctx);
return status;
}
int
gr_generic_vec_sub(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
status |= gr_sub(GR_ENTRY(res, i, sz), GR_ENTRY(src1, i, sz), GR_ENTRY(src2, i, sz), ctx);
return status;
}
int
gr_generic_vec_scalar_addmul(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
{
GR_TMP_START;
int status;
slong i, sz;
gr_ptr t;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
GR_TMP_INIT1(t, ctx);
for (i = 0; i < len; i++)
{
status |= gr_mul(t, GR_ENTRY(vec2, i, sz), c, ctx);
status |= gr_add(GR_ENTRY(vec1, i, sz), GR_ENTRY(vec1, i, sz), t, ctx);
}
GR_TMP_CLEAR1(t, ctx);
GR_TMP_END;
return status;
}
int
gr_generic_vec_scalar_submul(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
{
GR_TMP_START;
int status;
slong i, sz;
gr_ptr t;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
GR_TMP_INIT1(t, ctx);
for (i = 0; i < len; i++)
{
status |= gr_mul(t, GR_ENTRY(vec2, i, sz), c, ctx);
status |= gr_sub(GR_ENTRY(vec1, i, sz), GR_ENTRY(vec1, i, sz), t, ctx);
}
GR_TMP_CLEAR1(t, ctx);
GR_TMP_END;
return status;
}
int
gr_generic_vec_scalar_addmul_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
{
GR_TMP_START;
int status;
slong i, sz;
gr_ptr t;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
GR_TMP_INIT1(t, ctx);
for (i = 0; i < len; i++)
{
status |= gr_mul_si(t, GR_ENTRY(vec2, i, sz), c, ctx);
status |= gr_add(GR_ENTRY(vec1, i, sz), GR_ENTRY(vec1, i, sz), t, ctx);
}
GR_TMP_CLEAR1(t, ctx);
GR_TMP_END;
return status;
}
int
gr_generic_vec_scalar_submul_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
{
GR_TMP_START;
int status;
slong i, sz;
gr_ptr t;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
GR_TMP_INIT1(t, ctx);
for (i = 0; i < len; i++)
{
status |= gr_mul_si(t, GR_ENTRY(vec2, i, sz), c, ctx);
status |= gr_sub(GR_ENTRY(vec1, i, sz), GR_ENTRY(vec1, i, sz), t, ctx);
}
GR_TMP_CLEAR1(t, ctx);
GR_TMP_END;
return status;
}
int
gr_generic_vec_equal(int * res, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
{
int status, equal, this_equal;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
equal = 1;
for (i = 0; i < len; i++)
{
status |= gr_equal(&this_equal, GR_ENTRY(vec1, i, sz), GR_ENTRY(vec2, i, sz), ctx);
equal = equal && this_equal;
}
res[0] = equal;
return status;
}
int
gr_generic_vec_is_zero(int * res, gr_srcptr vec, slong len, gr_ctx_t ctx)
{
int status, equal, this_equal;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
equal = 1;
for (i = 0; i < len; i++)
{
status |= gr_is_zero(&this_equal, GR_ENTRY(vec, i, sz), ctx);
equal = equal && this_equal;
}
res[0] = equal;
return status;
}
int
gr_generic_vec_dot(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
{
GR_TMP_START;
int status;
slong i, sz;
gr_ptr t;
if (len <= 0)
{
if (initial == NULL)
return gr_zero(res, ctx);
else
return gr_set(res, initial, ctx);
}
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
GR_TMP_INIT1(t, ctx);
if (initial == NULL)
{
status |= gr_mul(res, vec1, vec2, ctx);
}
else
{
if (subtract)
status |= gr_neg(res, initial, ctx);
else
status |= gr_set(res, initial, ctx);
status |= gr_mul(t, vec1, vec2, ctx);
status |= gr_add(res, res, t, ctx);
}
for (i = 1; i < len; i++)
{
status |= gr_mul(t, GR_ENTRY(vec1, i, sz), GR_ENTRY(vec2, i, sz), ctx);
status |= gr_add(res, res, t, ctx);
}
if (subtract)
status |= gr_neg(res, res, ctx);
GR_TMP_CLEAR1(t, ctx);
GR_TMP_END;
return status;
}
int
gr_generic_vec_dot_rev(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
{
GR_TMP_START;
int status;
slong i, sz;
gr_ptr t;
if (len <= 0)
{
if (initial == NULL)
return gr_zero(res, ctx);
else
return gr_set(res, initial, ctx);
}
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
GR_TMP_INIT1(t, ctx);
if (initial == NULL)
{
status |= gr_mul(res, vec1, GR_ENTRY(vec2, len - 1, sz), ctx);
}
else
{
if (subtract)
status |= gr_neg(res, initial, ctx);
else
status |= gr_set(res, initial, ctx);
status |= gr_mul(t, vec1, GR_ENTRY(vec2, len - 1, sz), ctx);
status |= gr_add(res, res, t, ctx);
}
for (i = 1; i < len; i++)
{
status |= gr_mul(t, GR_ENTRY(vec1, i, sz), GR_ENTRY(vec2, len - 1 - i, sz), ctx);
status |= gr_add(res, res, t, ctx);
}
if (subtract)
status |= gr_neg(res, res, ctx);
GR_TMP_CLEAR1(t, ctx);
GR_TMP_END;
return status;
}
gr_vec_methods gr_vec_generic_methods = {
(gr_func_out_si) gr_generic_vec_init,
(gr_func_out_si) gr_generic_vec_clear,
(gr_func_out_out_si) gr_generic_vec_swap,
(gr_func_out_si) gr_generic_vec_zero,
(gr_func_out_in_si) gr_generic_vec_set,
(gr_func_out_in_si) gr_generic_vec_neg,
(gr_func_out_in_in_si) gr_generic_vec_add,
(gr_func_out_in_in_si) gr_generic_vec_sub,
(gr_func_out_in_si_in) gr_generic_vec_scalar_addmul,
(gr_func_out_in_si_in) gr_generic_vec_scalar_submul,
(gr_func_out_in_si_si) gr_generic_vec_scalar_addmul_si,
(gr_func_out_in_si_si) gr_generic_vec_scalar_submul_si,
(gr_func_outbool_in_in_si) gr_generic_vec_equal,
(gr_func_outbool_in_si) gr_generic_vec_is_zero,
(gr_func_out_in_bool_in_in_si) gr_generic_vec_dot,
(gr_func_out_in_bool_in_in_si) gr_generic_vec_dot_rev,
};
int
_gr_vec_init(gr_ptr vec, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->init(vec, len, ctx);
}
int
_gr_vec_clear(gr_ptr vec, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->clear(vec, len, ctx);
}
int
_gr_vec_swap(gr_ptr vec1, gr_ptr vec2, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->swap(vec1, vec2, len, ctx);
}
int
_gr_vec_zero(gr_ptr vec, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->zero(vec, len, ctx);
}
int
_gr_vec_set(gr_ptr res, gr_srcptr src, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->set(res, src, len, ctx);
}
int
_gr_vec_neg(gr_ptr res, gr_srcptr src, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->neg(res, src, len, ctx);
}
int
_gr_vec_add(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->add(res, src1, src2, len, ctx);
}
int
_gr_vec_sub(gr_ptr res, gr_srcptr src1, gr_srcptr src2, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->sub(res, src1, src2, len, ctx);
}
int
_gr_vec_scalar_addmul(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
{
return ctx->vec_methods->scalar_addmul(vec1, vec2, len, c, ctx);
}
int
_gr_vec_scalar_submul(gr_ptr vec1, gr_srcptr vec2, slong len, gr_srcptr c, gr_ctx_t ctx)
{
return ctx->vec_methods->scalar_submul(vec1, vec2, len, c, ctx);
}
int
_gr_vec_scalar_addmul_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
{
return ctx->vec_methods->scalar_addmul_si(vec1, vec2, len, c, ctx);
}
int
_gr_vec_scalar_submul_si(gr_ptr vec1, gr_srcptr vec2, slong len, slong c, gr_ctx_t ctx)
{
return ctx->vec_methods->scalar_submul_si(vec1, vec2, len, c, ctx);
}
int
_gr_vec_equal(int * res, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->equal(res, vec1, vec2, len, ctx);
}
int
_gr_vec_is_zero(int * res, gr_srcptr vec, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->is_zero(res, vec, len, ctx);
}
int
_gr_vec_dot(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->dot(res, initial, subtract, vec1, vec2, len, ctx);
}
int
_gr_vec_dot_rev(gr_ptr res, gr_srcptr initial, int subtract, gr_srcptr vec1, gr_srcptr vec2, slong len, gr_ctx_t ctx)
{
return ctx->vec_methods->dot_rev(res, initial, subtract, vec1, vec2, len, ctx);
}
int
_gr_vec_randtest(gr_ptr res, flint_rand_t state, slong len, void * options, gr_ctx_t ctx)
{
int status;
slong i, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
for (i = 0; i < len; i++)
{
status |= gr_randtest(GR_ENTRY(res, i, sz), state, options, ctx);
}
return status;
}
/* Some generic functions. */
/* Assumes exp >= 2; res and tmp not not aliased with x. */
static int
_gr_generic_pow_ui_binexp(gr_ptr res, gr_ptr tmp, gr_srcptr x, ulong exp, gr_ctx_t ctx)
{
gr_ptr R, S, T;
int status;
int zeros;
ulong bit;
status = GR_SUCCESS;
/* Determine parity due to swaps */
zeros = 0;
bit = exp;
while (bit > 1)
{
zeros += !(bit & 1);
bit >>= 1;
}
if (zeros % 2)
{
R = res;
S = tmp;
}
else
{
R = tmp;
S = res;
}
bit = UWORD(1) << (FLINT_BIT_COUNT(exp) - 2);
status |= gr_mul(R, x, x, ctx);
if (bit & exp)
{
status |= gr_mul(S, R, x, ctx);
T = R;
R = S;
S = T;
}
while (bit >>= 1)
{
status |= gr_mul(S, R, R, ctx);
if (bit & exp)
{
status |= gr_mul(R, S, x, ctx);
}
else
{
T = R;
R = S;
S = T;
}
}
return status;
}
int
gr_generic_pow_ui(gr_ptr res, gr_srcptr x, ulong e, gr_ctx_t ctx)
{
int status;
if (e == 0)
{
return gr_one(res, ctx);
}
else if (e == 1)
{
return gr_set(res, x, ctx);
}
else if (e == 2)
{
return gr_mul(res, x, x, ctx);
}
else if (e == 4)
{
status = gr_mul(res, x, x, ctx);
status |= gr_mul(res, res, res, ctx);
return status;
}
else
{
gr_ptr t, u;
GR_TMP_START;
if (res == x)
{
GR_TMP_INIT2(t, u, ctx);
gr_set(u, x, ctx);
status = _gr_generic_pow_ui_binexp(res, t, u, e, ctx);
GR_TMP_CLEAR2(t, u, ctx);
}
else
{
GR_TMP_INIT1(t, ctx);
status = _gr_generic_pow_ui_binexp(res, t, x, e, ctx);
GR_TMP_CLEAR1(t, ctx);
}
GR_TMP_END;
return status;
}
}
int
gr_generic_pow_si(gr_ptr res, gr_srcptr x, slong e, gr_ctx_t ctx)
{
if (e >= 0)
{
return gr_generic_pow_ui(res, x, e, ctx);
}
else
{
int status;
status = gr_inv(res, x, ctx);
if (status == GR_SUCCESS)
status = gr_generic_pow_ui(res, x, -e, ctx);
return status;
}
}
/* tbd */
int
gr_pow_ui(gr_ptr res, gr_srcptr x, ulong e, gr_ctx_t ctx)
{
return gr_generic_pow_ui(res, x, e, ctx);
}
int
gr_print(gr_srcptr x, gr_ctx_t ctx)
{
gr_stream_t out;
gr_stream_init_file(out, stdout);
gr_write(out, x, ctx);
return GR_SUCCESS;
}
int
gr_println(gr_srcptr x, gr_ctx_t ctx)
{
gr_stream_t out;
gr_stream_init_file(out, stdout);
gr_write(out, x, ctx);
gr_stream_write(out, "\n");
return GR_SUCCESS;
}
/* Todo: generic vectors */
typedef struct
{
gr_ptr entries;
slong length;
slong alloc;
}
gr_vec_struct;
typedef gr_vec_struct gr_vec_t[1];
/* Todo: generic matrices */
typedef struct
{
gr_ptr entries;
slong r;
slong c;
gr_ptr * rows;
}
gr_mat_struct;
typedef gr_mat_struct gr_mat_t[1];
#define GR_MAT_ENTRY(mat,i,j,ctx) GR_ENTRY((mat)->rows[i], j, sz)
#define gr_mat_nrows(mat, ctx) ((mat)->r)
#define gr_mat_ncols(mat, ctx) ((mat)->c)
int
gr_mat_init(gr_mat_t mat, slong rows, slong cols, gr_ctx_t ctx)
{
slong i, sz;
sz = ctx->sizeof_elem;
if (rows != 0)
mat->rows = flint_malloc(rows * sizeof(gr_ptr));
else
mat->rows = NULL;
if (rows != 0 && cols != 0)
{
mat->entries = (gr_ptr) flint_malloc(flint_mul_sizes(rows, cols) * sz);
_gr_vec_init(mat->entries, rows * cols, ctx);
for (i = 0; i < rows; i++)
mat->rows[i] = GR_ENTRY(mat->entries, i * cols, sz);
}
else
{
mat->entries = NULL;
for (i = 0; i < rows; i++)
mat->rows[i] = NULL;
}
mat->r = rows;
mat->c = cols;
return GR_SUCCESS;
}
int
gr_mat_clear(gr_mat_t mat, gr_ctx_t ctx)
{
if (mat->entries != NULL)
{
_gr_vec_clear(mat->entries, mat->r * mat->c, ctx);
flint_free(mat->entries);
flint_free(mat->rows);
}
else if (mat->r != 0)
{
flint_free(mat->rows);
}
return GR_SUCCESS;
}
int
gr_mat_randtest(gr_mat_t mat, flint_rand_t state, void * options, gr_ctx_t ctx)
{
int status;
slong i, r, c;
r = gr_mat_nrows(mat, ctx);
c = gr_mat_ncols(mat, ctx);
status = GR_SUCCESS;
for (i = 0; i < r; i++)
{
status |= _gr_vec_randtest(mat->rows[i], state, c, options, ctx);
}
return status;
}
int
gr_mat_equal(int * res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
{
int status, equal, this_equal;
slong i, r, c, sz;
sz = ctx->sizeof_elem;
status = GR_SUCCESS;
r = gr_mat_nrows(mat1, ctx);
c = gr_mat_ncols(mat1, ctx);
if (r != gr_mat_nrows(mat2, ctx) ||
c != gr_mat_ncols(mat2, ctx))
{
res[0] = 0;
return GR_SUCCESS;
}
if (r == 0 || c == 0)
{
res[0] = 1;
return GR_SUCCESS;
}
equal = 1;
for (i = 0; i < r; i++)
{
status |= _gr_vec_equal(&this_equal, mat1->rows[i], mat2->rows[i], c, ctx);
equal = equal && this_equal;
}
res[0] = equal;
return status;
}
int
gr_mat_zero(gr_mat_t res, gr_ctx_t ctx)
{
int status;
slong i, r, c;
r = gr_mat_nrows(res, ctx);
c = gr_mat_ncols(res, ctx);
status = GR_SUCCESS;
for (i = 0; i < r; i++)
{
status |= _gr_vec_zero(res->rows[i], c, ctx);
}
return status;
}
int
gr_mat_set_si(gr_mat_t res, slong v, gr_ctx_t ctx)
{
int status;
slong i, r, c, sz;
r = gr_mat_nrows(res, ctx);
c = gr_mat_ncols(res, ctx);
sz = ctx->sizeof_elem;
status = gr_mat_zero(res, ctx);
for (i = 0; i < FLINT_MIN(r, c); i++)
status |= gr_set_si(GR_MAT_ENTRY(res, i, i, sz), v, ctx);
return status;
}
int
gr_mat_one(gr_mat_t res, gr_ctx_t ctx)
{
return gr_mat_set_si(res, 1, ctx);
}
int
gr_mat_set(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
{
int status;
slong i, r, c;
r = gr_mat_nrows(res, ctx);
c = gr_mat_ncols(res, ctx);
if (r != gr_mat_nrows(mat, ctx) ||
c != gr_mat_ncols(mat, ctx))
{
return GR_DOMAIN;
}
status = GR_SUCCESS;
if (res != mat)
{
for (i = 0; i < r; i++)
{
status |= _gr_vec_set(res->rows[i], mat->rows[i], c, ctx);
}
}
return status;
}
int
gr_mat_neg(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
{
int status;
slong i, r, c;
r = gr_mat_nrows(res, ctx);
c = gr_mat_ncols(res, ctx);
if (r != gr_mat_nrows(mat, ctx) ||
c != gr_mat_ncols(mat, ctx))
{
return GR_DOMAIN;
}
status = GR_SUCCESS;
for (i = 0; i < r; i++)
{
status |= _gr_vec_neg(res->rows[i], mat->rows[i], c, ctx);
}
return status;
}
int
gr_mat_swap_entrywise(gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
{
int status;
slong i, r, c;
r = gr_mat_nrows(mat1, ctx);
c = gr_mat_ncols(mat1, ctx);
if (r != gr_mat_nrows(mat2, ctx) ||
c != gr_mat_ncols(mat2, ctx))
{
return GR_DOMAIN;
}
status = GR_SUCCESS;
for (i = 0; i < r; i++)
{
status |= _gr_vec_swap(mat1->rows[i], mat2->rows[i], c, ctx);
}
return status;
}
int
gr_mat_add(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
{
int status;
slong i, r, c;
r = gr_mat_nrows(res, ctx);
c = gr_mat_ncols(res, ctx);
if (r != gr_mat_nrows(mat1, ctx) ||
c != gr_mat_ncols(mat1, ctx) ||
r != gr_mat_nrows(mat2, ctx) ||
c != gr_mat_ncols(mat2, ctx))
{
return GR_DOMAIN;
}
status = GR_SUCCESS;
for (i = 0; i < r; i++)
{
status |= _gr_vec_add(res->rows[i], mat1->rows[i], mat2->rows[i], c, ctx);
}
return status;
}
int
gr_mat_sub(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
{
int status;
slong i, r, c;
r = gr_mat_nrows(res, ctx);
c = gr_mat_ncols(res, ctx);
if (r != gr_mat_nrows(mat1, ctx) ||
c != gr_mat_ncols(mat1, ctx) ||
r != gr_mat_nrows(mat2, ctx) ||
c != gr_mat_ncols(mat2, ctx))
{
return GR_DOMAIN;
}
status = GR_SUCCESS;
for (i = 0; i < r; i++)
{
status |= _gr_vec_sub(res->rows[i], mat1->rows[i], mat2->rows[i], c, ctx);
}
return status;
}
// todo status
int
gr_mat_print(const gr_mat_t mat, gr_ctx_t ctx)
{
int status;
slong r, c;
slong i, j;
slong sz;
sz = ctx->sizeof_elem;
r = gr_mat_nrows(mat, ctx);
c = gr_mat_ncols(mat, ctx);
status = GR_SUCCESS;
flint_printf("[");
for (i = 0; i < r; i++)
{
flint_printf("[");
for (j = 0; j < c; j++)
{
status |= gr_print(GR_MAT_ENTRY(mat, i, j, sz), ctx);
if (j < c - 1)
flint_printf(", ");
}
if (i < r - 1)
flint_printf("],\n");
else
flint_printf("]");
}
flint_printf("]\n");
return status;
}
int
gr_mat_mul_classical(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
{
slong ar, ac, br, bc, i, j, k, sz;
int status;
ar = gr_mat_nrows(A, ctx);
ac = gr_mat_ncols(A, ctx);
br = gr_mat_nrows(B, ctx);
bc = gr_mat_ncols(B, ctx);
if (ac != br || ar != gr_mat_nrows(C, ctx) || bc != gr_mat_ncols(C, ctx))
return GR_DOMAIN;
if (br == 0)
{
return gr_mat_zero(C, ctx);
}
status = GR_SUCCESS;
if (A == C || B == C)
{
gr_mat_t T;
status |= gr_mat_init(T, ar, bc, ctx);
status |= gr_mat_mul_classical(T, A, B, ctx);
status |= gr_mat_swap_entrywise(T, C, ctx);
status |= gr_mat_clear(T, ctx);
return status;
}
sz = ctx->sizeof_elem;
if (br == 1)
{
for (i = 0; i < ar; i++)
{
for (j = 0; j < bc; j++)
{
status |= gr_mul(GR_MAT_ENTRY(C, i, j, sz),
GR_MAT_ENTRY(A, i, 0, sz),
GR_MAT_ENTRY(B, 0, j, sz), ctx);
}
}
}
else
{
gr_ptr tmp;
TMP_INIT;
TMP_START;
tmp = TMP_ALLOC(sz * br * bc);
for (i = 0; i < br; i++)
for (j = 0; j < bc; j++)
memcpy(GR_ENTRY(tmp, j * br + i, sz), GR_MAT_ENTRY(B, i, j, sz), sz);
for (i = 0; i < ar; i++)
{
for (j = 0; j < bc; j++)
{
status |= _gr_vec_dot(GR_MAT_ENTRY(C, i, j, sz), NULL, 0,
GR_MAT_ENTRY(A, i, 0, sz), GR_ENTRY(tmp, j * br, sz), br, ctx);
}
}
TMP_END;
}
return status;
}
/* Todo: generic polynomials */
typedef struct
{
gr_ptr coeffs;
slong length;
slong alloc;
}
gr_poly_struct;
typedef gr_poly_struct gr_poly_t[1];
/*
A ring to test.
*/
typedef unsigned char nmod8_struct;
typedef nmod8_struct nmod8_t[1];
int
nmod8_init(nmod8_t x, nmod_t * ctx)
{
x[0] = 0;
return GR_SUCCESS;
}
int
nmod8_clear(nmod8_t x, nmod_t * ctx)
{
return GR_SUCCESS;
}
int
nmod8_swap(nmod8_t x, nmod8_t y, nmod_t * ctx)
{
nmod8_t t;
*t = *x;
*x = *y;
*y = *t;
return GR_SUCCESS;
}
int
nmod8_randtest(nmod8_t res, flint_rand_t state, const void * options, nmod_t * ctx)
{
res[0] = n_randtest(state) % ctx->n;
return GR_SUCCESS;
}
int
nmod8_write(gr_stream_t out, const nmod8_t x, nmod_t * ctx)
{
gr_stream_write_si(out, x[0]);
return GR_SUCCESS;
}
int
nmod8_zero(nmod8_t x, nmod_t * ctx)
{
x[0] = 0;
return GR_SUCCESS;
}
int
nmod8_one(nmod8_t x, nmod_t * ctx)
{
x[0] = (ctx->n != 0);
return GR_SUCCESS;
}
int
nmod8_set_si(nmod8_t res, slong v, nmod_t * ctx)
{
ulong t;
t = FLINT_ABS(v);
NMOD_RED(t, t, *ctx);
if (v < 0)
t = nmod_neg(t, *ctx);
res[0] = t;
return GR_SUCCESS;
}
int
nmod8_is_zero(int * res, const nmod8_t x, nmod_t * ctx)
{
res[0] = (x[0] == 0);
return GR_SUCCESS;
}
int
nmod8_is_one(int * res, const nmod8_t x, nmod_t * ctx)
{
res[0] = (x[0] == (ctx->n != 1));
return GR_SUCCESS;
}
int
nmod8_equal(int * res, const nmod8_t x, const nmod8_t y, nmod_t * ctx)
{
res[0] = (x[0] == y[0]);
return GR_SUCCESS;
}
int
nmod8_set(nmod8_t res, const nmod8_t x, nmod_t * ctx)
{
res[0] = x[0];
return GR_SUCCESS;
}
int
nmod8_neg(nmod8_t res, const nmod8_t x, nmod_t * ctx)
{
res[0] = nmod_neg(x[0], *ctx);
return GR_SUCCESS;
}
int
nmod8_add(nmod8_t res, const nmod8_t x, const nmod8_t y, nmod_t * ctx)
{
res[0] = nmod_add(x[0], y[0], *ctx);
return GR_SUCCESS;
}
int
nmod8_add_si(nmod8_t res, const nmod8_t x, slong y, nmod_t * ctx)
{
nmod8_t t;
nmod8_set_si(t, y, ctx);
return nmod8_add(res, x, t, ctx);
}
int
nmod8_sub(nmod8_t res, const nmod8_t x, const nmod8_t y, nmod_t * ctx)
{
res[0] = nmod_sub(x[0], y[0], *ctx);
return GR_SUCCESS;
}
int
nmod8_mul(nmod8_t res, const nmod8_t x, const nmod8_t y, nmod_t * ctx)
{
res[0] = nmod_mul(x[0], y[0], *ctx);
return GR_SUCCESS;
}
int
nmod8_mul_si(nmod8_t res, const nmod8_t x, slong y, nmod_t * ctx)
{
nmod8_t t;
nmod8_set_si(t, y, ctx);
return nmod8_mul(res, x, t, ctx);
}
int
nmod8_inv(nmod8_t res, const nmod8_t x, nmod_t * ctx)
{
ulong r, g;
g = n_gcdinv(&r, x[0], ctx->n);
if (g == 1)
{
res[0] = r;
return GR_SUCCESS;
}
else
{
res[0] = 0;
return GR_DOMAIN;
}
}
int
nmod8_div(nmod8_t res, const nmod8_t x, const nmod8_t y, nmod_t * ctx)
{
nmod8_t t;
int status;
status = nmod8_inv(t, y, ctx);
nmod8_mul(res, x, t, ctx);
return status;
}
int
nmod8_div_si(nmod8_t res, const nmod8_t x, slong y, nmod_t * ctx)
{
nmod8_t t;
nmod8_set_si(t, y, ctx);
return nmod8_div(res, x, t, ctx);
}
int
nmod8_is_invertible(int * res, const nmod8_t x, nmod_t * ctx)
{
ulong r, g;
g = n_gcdinv(&r, x[0], ctx->n);
res[0] = (g == 1);
return GR_SUCCESS;
}
gr_methods nmod8_methods = {
(gr_func_out) nmod8_init,
(gr_func_out) nmod8_clear,
(gr_func_out_out) nmod8_swap,
(gr_func_randtest) nmod8_randtest,
(gr_func_stream_in) nmod8_write,
(gr_func_out) nmod8_zero,
(gr_func_out) nmod8_one,
(gr_func_outbool_in) nmod8_is_zero,
(gr_func_outbool_in) nmod8_is_one,
(gr_func_outbool_in_in) nmod8_equal,
(gr_func_out_in) nmod8_set,
(gr_func_out_si) nmod8_set_si,
(gr_func_out_in) nmod8_neg,
(gr_func_out_in_in) nmod8_add,
(gr_func_out_in_si) nmod8_add_si,
(gr_func_out_in_in) nmod8_sub,
(gr_func_out_in_in) nmod8_mul,
(gr_func_out_in_si) nmod8_mul_si,
(gr_func_out_in_in) nmod8_div,
(gr_func_out_in_in) nmod8_div,
(gr_func_out_in_si) nmod8_div_si,
(gr_func_outbool_in) nmod8_is_invertible,
(gr_func_out_in) nmod8_inv,
};
void
gr_ctx_init_nmod8(gr_ctx_t ctx, unsigned char n)
{
ctx->flags = GR_COMMUTATIVE_RING;
ctx->sizeof_elem = sizeof(nmod8_struct);
ctx->elem_ctx = flint_malloc(sizeof(nmod_t)); /* This could be something more interesting */
nmod_init(ctx->elem_ctx, n);
ctx->methods = &nmod8_methods;
ctx->vec_methods = &gr_vec_generic_methods;
ctx->debug_string = "nmod8 ring";
}
void
gr_ctx_clear_nmod8(gr_ctx_t ctx)
{
flint_free(ctx->elem_ctx);
}
/* Matrix ring to test */
typedef struct
{
gr_ctx_struct * base_ring;
slong n;
}
matrix_ctx_t;
int
matrix_init(gr_mat_t res, matrix_ctx_t * ctx)
{
return gr_mat_init(res, ctx->n, ctx->n, ctx->base_ring);
}
int
matrix_clear(gr_mat_t res, matrix_ctx_t * ctx)
{
return gr_mat_clear(res, ctx->base_ring);
}
/* TODO UNIFY */
int
matrix_write(gr_stream_t out, gr_mat_t res, matrix_ctx_t * ctx)
{
return gr_mat_print(res, ctx->base_ring);
}
int
matrix_randtest(gr_mat_t res, flint_rand_t state, void * options, matrix_ctx_t * ctx)
{
return gr_mat_randtest(res, state, options, ctx->base_ring);
}
int
matrix_equal(int * equal, const gr_mat_t mat1, const gr_mat_t mat2, matrix_ctx_t * ctx)
{
return gr_mat_equal(equal, mat1, mat2, ctx->base_ring);
}
int
matrix_set(gr_mat_t res, const gr_mat_t mat, matrix_ctx_t * ctx)
{
return gr_mat_set(res, mat, ctx->base_ring);
}
int
matrix_zero(gr_mat_t res, matrix_ctx_t * ctx)
{
return gr_mat_zero(res, ctx->base_ring);
}
int
matrix_one(gr_mat_t res, matrix_ctx_t * ctx)
{
return gr_mat_one(res, ctx->base_ring);
}
int
matrix_add(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, matrix_ctx_t * ctx)
{
return gr_mat_add(res, mat1, mat2, ctx->base_ring);
}
int
matrix_mul(gr_mat_t res, const gr_mat_t mat1, const gr_mat_t mat2, matrix_ctx_t * ctx)
{
return gr_mat_mul_classical(res, mat1, mat2, ctx->base_ring);
}
#define matrix_swap gr_not_implemented
#define matrix_is_zero gr_not_implemented
#define matrix_is_one gr_not_implemented
#define matrix_set_si gr_not_implemented
#define matrix_neg gr_not_implemented
#define matrix_sub gr_not_implemented
gr_methods matrix_methods = {
(gr_func_out) matrix_init,
(gr_func_out) matrix_clear,
(gr_func_out_out) matrix_swap,
(gr_func_randtest) matrix_randtest,
(gr_func_stream_in) matrix_write,
(gr_func_out) matrix_zero,
(gr_func_out) matrix_one,
(gr_func_outbool_in) matrix_is_zero,
(gr_func_outbool_in) matrix_is_one,
(gr_func_outbool_in_in) matrix_equal,
(gr_func_out_in) matrix_set,
(gr_func_out_si) matrix_set_si,
(gr_func_out_in) matrix_neg,
(gr_func_out_in_in) matrix_add,
(gr_func_out_in_si) gr_not_implemented,
(gr_func_out_in_in) matrix_sub,
(gr_func_out_in_in) matrix_mul,
(gr_func_out_in_si) gr_not_implemented,
(gr_func_out_in_in) gr_not_implemented,
(gr_func_out_in_in) gr_not_implemented,
(gr_func_out_in_si) gr_not_implemented,
(gr_func_outbool_in) gr_not_implemented,
(gr_func_out_in) gr_not_implemented,
};
/* rename: gr_ctx_init_gr_mat */
void
gr_ctx_init_matrix(gr_ctx_t ctx, gr_ctx_t base_ring, slong n)
{
ctx->flags = 0;
ctx->sizeof_elem = sizeof(gr_mat_struct);
ctx->elem_ctx = flint_malloc(sizeof(matrix_ctx_t));
((matrix_ctx_t *) ctx->elem_ctx)->base_ring = (gr_ctx_struct *) base_ring;
((matrix_ctx_t *) ctx->elem_ctx)->n = n;
ctx->methods = &matrix_methods;
ctx->vec_methods = &gr_vec_generic_methods;
ctx->debug_string = "matrix ring";
}
void
gr_ctx_clear_matrix(gr_ctx_t ctx)
{
flint_free(ctx->elem_ctx);
}
#define TIMEIT_PRINT2(__timer, __reps) \
flint_printf("cpu/wall(s): %g %g", \
__timer->cpu*0.001/__reps, __timer->wall*0.001 / __reps);
#define TIMEIT_REPEAT(__timer, __reps) \
do \
{ \
slong __timeit_k; \
__reps = 1; \
while (1) \
{ \
timeit_start(__timer); \
for (__timeit_k = 0; __timeit_k < __reps; __timeit_k++) \
{
#define TIMEIT_END_REPEAT(__timer, __reps) \
} \
timeit_stop(__timer); \
if (__timer->cpu >= 100) \
break; \
__reps *= 10; \
} \
} while (0);
#define TIMEIT_START \
do { \
timeit_t __timer; slong __reps; \
TIMEIT_REPEAT(__timer, __reps)
#define TIMEIT_STOP \
TIMEIT_END_REPEAT(__timer, __reps) \
TIMEIT_PRINT(__timer, __reps) \
} while (0);
#define TIMEIT_ONCE_START \
do \
{ \
timeit_t __timer; \
timeit_start(__timer); \
do { \
#define TIMEIT_ONCE_STOP \
} while (0); \
timeit_stop(__timer); \
TIMEIT_PRINT(__timer, 1) \
} while (0); \
typedef int ((*gr_test_function)(gr_ctx_t, flint_rand_t, int));
int
gr_test_add_associativity(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal;
GR_TMP_START;
gr_ptr x, y, z;
gr_ptr xy, yz, xy_z, x_yz;
GR_TMP_INIT3(x, y, z, R);
GR_TMP_INIT4(xy, yz, xy_z, x_yz, R);
gr_randtest(x, state, NULL, R);
gr_randtest(y, state, NULL, R);
gr_randtest(z, state, NULL, R);
gr_randtest(xy, state, NULL, R);
gr_randtest(yz, state, NULL, R);
status = GR_SUCCESS;
status |= gr_add(xy, x, y, R);
status |= gr_add(yz, y, z, R);
status |= gr_add(xy_z, xy, z, R);
status |= gr_add(x_yz, x, yz, R);
status |= gr_equal(&equal, xy_z, x_yz, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("x = \n"); gr_println(x, R);
printf("y = \n"); gr_println(y, R);
printf("z = \n"); gr_println(z, R);
printf("x + y = \n"); gr_println(xy, R);
printf("y + z = \n"); gr_println(yz, R);
printf("(x + y) + z = \n"); gr_println(xy_z, R);
printf("x + (y + z) = \n"); gr_println(x_yz, R);
printf("\n");
}
GR_TMP_CLEAR3(x, y, z, R);
GR_TMP_CLEAR4(xy, yz, xy_z, x_yz, R);
GR_TMP_END;
return status;
}
int
gr_test_add_commutativity(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal;
GR_TMP_START;
gr_ptr x, y, xy, yx;
GR_TMP_INIT4(x, y, xy, yx, R);
gr_randtest(x, state, NULL, R);
gr_randtest(y, state, NULL, R);
status = GR_SUCCESS;
status |= gr_add(xy, x, y, R);
status |= gr_add(yx, y, x, R);
status |= gr_equal(&equal, xy, yx, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("x = \n"); gr_println(x, R);
printf("y = \n"); gr_println(y, R);
printf("y + y = \n"); gr_println(xy, R);
printf("y + x = \n"); gr_println(yx, R);
printf("\n");
}
GR_TMP_CLEAR4(x, y, xy, yx, R);
GR_TMP_END;
return status;
}
int
gr_test_add_aliasing(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal, alias;
GR_TMP_START;
gr_ptr x, y, xy1, xy2;
GR_TMP_INIT4(x, y, xy1, xy2, R);
gr_randtest(x, state, NULL, R);
gr_randtest(y, state, NULL, R);
status = GR_SUCCESS;
status |= gr_add(xy1, x, y, R);
alias = n_randint(state, 4);
switch (alias)
{
case 0:
status |= gr_add(xy1, x, y, R);
gr_set(xy2, x, R);
status |= gr_add(xy2, xy2, y, R);
break;
case 1:
status |= gr_add(xy1, x, y, R);
gr_set(xy2, y, R);
status |= gr_add(xy2, x, y, R);
break;
case 2:
gr_set(y, x, R);
status |= gr_add(xy1, x, y, R);
status |= gr_add(xy2, x, x, R);
break;
default:
gr_set(y, x, R);
gr_set(xy2, x, R);
status |= gr_add(xy1, x, y, R);
status |= gr_add(xy2, xy2, xy2, R);
}
status |= gr_equal(&equal, xy1, xy2, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("alias: %d\n", alias);
printf("x = "); gr_println(x, R);
printf("y = "); gr_println(y, R);
printf("y + y (1) = "); gr_println(xy1, R);
printf("x + y (2) = "); gr_println(xy2, R);
printf("\n");
}
GR_TMP_CLEAR4(x, y, xy1, xy2, R);
GR_TMP_END;
return status;
}
int
gr_test_inv_involution(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal;
GR_TMP_START;
gr_ptr x, x_inv, x_inv_inv;
GR_TMP_INIT3(x, x_inv, x_inv_inv, R);
gr_randtest(x, state, NULL, R);
gr_randtest(x_inv, state, NULL, R);
gr_randtest(x_inv_inv, state, NULL, R);
status = GR_SUCCESS;
status |= gr_inv(x_inv, x, R);
status |= gr_inv(x_inv_inv, x_inv, R);
status |= gr_equal(&equal, x, x_inv_inv, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("x = \n"); gr_println(x, R);
printf("x ^ -1 = \n"); gr_println(x_inv, R);
printf("(x ^ -1) ^ -1 = \n"); gr_println(x_inv_inv, R);
printf("\n");
}
GR_TMP_CLEAR3(x, x_inv, x_inv_inv, R);
GR_TMP_END;
return status;
}
int
gr_test_inv_multiplication(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal1, equal2;
GR_TMP_START;
gr_ptr x, x_inv, x_inv_x, x_x_inv;
GR_TMP_INIT4(x, x_inv, x_inv_x, x_x_inv, R);
gr_randtest(x, state, NULL, R);
gr_randtest(x_inv, state, NULL, R);
gr_randtest(x_inv_x, state, NULL, R);
gr_randtest(x_x_inv, state, NULL, R);
status = GR_SUCCESS;
status |= gr_inv(x_inv, x, R);
status |= gr_mul(x_inv_x, x_inv, x, R);
status |= gr_mul(x_x_inv, x, x_inv, R);
status |= gr_is_one(&equal1, x_inv_x, R);
status |= gr_is_one(&equal2, x_x_inv, R);
if (status == GR_SUCCESS && (!equal1 || !equal2))
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("x = \n"); gr_println(x, R);
printf("x ^ -1 = \n"); gr_println(x_inv, R);
printf("(x ^ -1) * x = \n"); gr_println(x_inv_x, R);
printf("x * (x ^ -1) = \n"); gr_println(x_x_inv, R);
printf("\n");
}
GR_TMP_CLEAR4(x, x_inv, x_inv_x, x_x_inv, R);
GR_TMP_END;
return status;
}
int
gr_test_pow_ui_exponent_addition(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal;
GR_TMP_START;
ulong a, b;
gr_ptr x, xa, xb, xab, xaxb;
GR_TMP_INIT5(x, xa, xb, xab, xaxb, R);
gr_randtest(x, state, NULL, R);
gr_randtest(xa, state, NULL, R);
gr_randtest(xb, state, NULL, R);
gr_randtest(xab, state, NULL, R);
gr_randtest(xaxb, state, NULL, R);
do {
a = n_randtest(state);
b = n_randtest(state);
} while (a + b < a);
status = GR_SUCCESS;
status |= gr_pow_ui(xa, x, a, R);
status |= gr_pow_ui(xb, x, b, R);
status |= gr_pow_ui(xab, x, a + b, R);
status |= gr_mul(xaxb, xa, xb, R);
status |= gr_equal(&equal, xab, xaxb, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("x = \n"); gr_println(x, R);
flint_printf("a = %wu\n", a);
flint_printf("b = %wu\n", b);
printf("x ^ a = \n"); gr_println(xa, R);
printf("x ^ b = \n"); gr_println(xb, R);
printf("x ^ (a + b) = \n"); gr_println(xab, R);
printf("(x ^ a) * (x ^ b) = \n"); gr_println(xaxb, R);
printf("\n");
}
GR_TMP_CLEAR5(x, xa, xb, xab, xaxb, R);
GR_TMP_END;
return status;
}
int
gr_test_pow_ui_base_scalar_multiplication(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal;
GR_TMP_START;
ulong a;
slong y;
gr_ptr x, xa, ya, xya, xaya;
GR_TMP_INIT3(x, xa, ya, R);
GR_TMP_INIT2(xya, xaya, R);
gr_randtest(x, state, NULL, R);
gr_randtest(xa, state, NULL, R);
gr_randtest(ya, state, NULL, R);
y = n_randtest(state);
a = n_randtest(state);
status = GR_SUCCESS;
status |= gr_pow_ui(xa, x, a, R);
status |= gr_set_si(ya, y, R);
status |= gr_pow_ui(ya, ya, a, R);
status |= gr_set_si(xya, y, R);
status |= gr_mul(xya, x, xya, R); /* todo mul_si */
status |= gr_pow_ui(xya, xya, a, R);
status |= gr_mul(xaya, xa, ya, R);
status |= gr_equal(&equal, xya, xaya, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("x = \n"); gr_println(x, R);
flint_printf("y = %wd\n", y);
flint_printf("a = %wu\n", a);
printf("x ^ a = \n"); gr_println(xa, R);
printf("y ^ a = \n"); gr_println(ya, R);
printf("(x * y) ^ a = \n"); gr_println(xya, R);
printf("(x ^ a) * (y ^ a) = \n"); gr_println(xaya, R);
printf("\n");
}
GR_TMP_CLEAR3(x, xa, ya, R);
GR_TMP_CLEAR2(xya, xaya, R);
GR_TMP_END;
return status;
}
int
gr_test_pow_ui_base_multiplication(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal;
GR_TMP_START;
ulong a;
gr_ptr x, y, xa, ya, xya, xaya;
GR_TMP_INIT4(x, y, xa, ya, R);
GR_TMP_INIT2(xya, xaya, R);
gr_randtest(x, state, NULL, R);
gr_randtest(y, state, NULL, R);
gr_randtest(xa, state, NULL, R);
gr_randtest(ya, state, NULL, R);
a = n_randtest(state);
status = GR_SUCCESS;
status |= gr_pow_ui(xa, x, a, R);
status |= gr_pow_ui(ya, y, a, R);
status |= gr_mul(xya, x, y, R);
status |= gr_pow_ui(xya, xya, a, R);
status |= gr_mul(xaya, xa, ya, R);
status |= gr_equal(&equal, xya, xaya, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("x = \n"); gr_println(x, R);
printf("y = \n"); gr_println(y, R);
flint_printf("a = %wu\n", a);
printf("x ^ a = \n"); gr_println(xa, R);
printf("y ^ a = \n"); gr_println(ya, R);
printf("(x * y) ^ a = \n"); gr_println(xya, R);
printf("(x ^ a) * (y ^ a) = \n"); gr_println(xaya, R);
printf("\n");
}
GR_TMP_CLEAR4(x, y, xa, ya, R);
GR_TMP_CLEAR2(xya, xaya, R);
GR_TMP_END;
return status;
}
int
gr_test_pow_ui_aliasing(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal;
GR_TMP_START;
ulong a;
gr_ptr x, xa1, xa2;
GR_TMP_INIT3(x, xa1, xa2, R);
gr_randtest(x, state, NULL, R);
gr_randtest(xa1, state, NULL, R);
a = n_randtest(state);
status = GR_SUCCESS;
status |= gr_pow_ui(xa1, x, a, R);
status |= gr_set(xa2, x, R);
status |= gr_pow_ui(xa2, xa2, a, R);
status |= gr_equal(&equal, xa1, xa2, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
printf("\n");
printf("x = \n"); gr_println(x, R);
flint_printf("a = %wu\n", a);
printf("x ^ a (1) = \n"); gr_println(xa1, R);
printf("x ^ a (2) = \n"); gr_println(xa2, R);
printf("\n");
}
GR_TMP_CLEAR3(x, xa1, xa2, R);
GR_TMP_END;
return status;
}
int
gr_test_vec_add(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, aliasing, equal;
slong len;
GR_TMP_START;
gr_ptr x, y, xy1, xy2;
len = n_randint(state, 10);
GR_TMP_INIT_VEC(x, len, R);
GR_TMP_INIT_VEC(y, len, R);
GR_TMP_INIT_VEC(xy1, len, R);
GR_TMP_INIT_VEC(xy2, len, R);
_gr_vec_randtest(x, state, len, NULL, R);
_gr_vec_randtest(y, state, len, NULL, R);
_gr_vec_randtest(xy1, state, len, NULL, R);
_gr_vec_randtest(xy2, state, len, NULL, R);
status = GR_SUCCESS;
aliasing = n_randint(state, 4);
switch (aliasing)
{
case 0:
status |= _gr_vec_set(xy1, x, len, R);
status |= _gr_vec_add(xy1, xy1, y, len, R);
break;
case 1:
status |= _gr_vec_set(xy1, y, len, R);
status |= _gr_vec_add(xy1, x, xy1, len, R);
break;
case 2:
status |= _gr_vec_set(y, x, len, R);
status |= _gr_vec_set(xy1, x, len, R);
status |= _gr_vec_add(xy1, xy1, xy1, len, R);
break;
default:
status |= _gr_vec_add(xy1, x, y, len, R);
}
status |= gr_generic_vec_add(xy2, x, y, len, R);
status |= _gr_vec_equal(&equal, xy1, xy2, len, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
/* todo: vec print */
printf("\n");
printf("aliasing: %d\n", aliasing);
printf("\n");
}
GR_TMP_CLEAR_VEC(x, len, R);
GR_TMP_CLEAR_VEC(y, len, R);
GR_TMP_CLEAR_VEC(xy1, len, R);
GR_TMP_CLEAR_VEC(xy2, len, R);
GR_TMP_END;
return status;
}
/* (AB)C = A(BC) */
int
gr_test_mat_mul_classical_associativity(gr_ctx_t R, flint_rand_t state, int verbose)
{
int status, equal;
gr_mat_t A, B, C, AB, BC, AB_C, A_BC;
slong m, n, p, q;
m = n_randint(state, 5);
n = n_randint(state, 5);
p = n_randint(state, 5);
q = n_randint(state, 5);
gr_mat_init(A, m, n, R);
gr_mat_init(B, n, p, R);
gr_mat_init(C, p, q, R);
gr_mat_init(AB, m, p, R);
gr_mat_init(BC, n, q, R);
gr_mat_init(AB_C, m, q, R);
gr_mat_init(A_BC, m, q, R);
gr_mat_randtest(A, state, NULL, R);
gr_mat_randtest(B, state, NULL, R);
gr_mat_randtest(C, state, NULL, R);
gr_mat_randtest(AB, state, NULL, R);
gr_mat_randtest(BC, state, NULL, R);
gr_mat_randtest(AB_C, state, NULL, R);
gr_mat_randtest(A_BC, state, NULL, R);
status |= gr_mat_mul_classical(AB, A, B, R);
status |= gr_mat_mul_classical(BC, B, C, R);
status |= gr_mat_mul_classical(AB_C, AB, C, R);
status |= gr_mat_mul_classical(A_BC, A, BC, R);
status |= gr_mat_equal(&equal, AB_C, A_BC, R);
if (status == GR_SUCCESS && !equal)
{
status = GR_WRONG;
}
if (verbose || status == GR_WRONG)
{
/* todo: vec print */
printf("\n");
printf("A = \n"); gr_mat_print(A, R); printf("\n");
printf("B = \n"); gr_mat_print(B, R); printf("\n");
printf("C = \n"); gr_mat_print(C, R); printf("\n");
printf("AB = \n"); gr_mat_print(AB, R); printf("\n");
printf("BC = \n"); gr_mat_print(BC, R); printf("\n");
printf("AB * C = \n"); gr_mat_print(AB_C, R); printf("\n");
printf("A * BC = \n"); gr_mat_print(A_BC, R); printf("\n");
printf("\n");
}
gr_mat_clear(A, R);
gr_mat_clear(B, R);
gr_mat_clear(C, R);
gr_mat_clear(AB, R);
gr_mat_clear(BC, R);
gr_mat_clear(A_BC, R);
gr_mat_clear(AB_C, R);
return status;
}
void
gr_test_iter(gr_ctx_t R, flint_rand_t state, const char * descr, gr_test_function func, slong iters)
{
slong iter, count_success, count_unable, count_domain;
int status;
timeit_t timer;
count_success = 0;
count_unable = 0;
count_domain = 0;
printf("%s ... ", descr);
fflush(stdout);
timeit_start(timer);
for (iter = 0; iter < iters; iter++)
{
// flint_printf("iter %ld\n", iter);
status = func(R, state, 0);
if (status == GR_SUCCESS)
count_success++;
if (status & GR_UNABLE)
count_unable++;
if (status & GR_DOMAIN)
count_domain++;
if (status & GR_WRONG)
{
flint_printf("\nFAIL\n");
abort();
}
}
timeit_stop(timer);
flint_printf("PASS (%wd successful, %wd domain, %wd unable, 0 wrong, %.3g cpu, %.3g wall)\n",
count_success, count_domain, count_unable, timer->cpu*0.001, timer->wall*0.001);
}
void
gr_test_ring(gr_ctx_t R, slong iters)
{
flint_rand_t state;
flint_randinit(state);
gr_test_iter(R, state, "add: associativity", gr_test_add_associativity, iters);
gr_test_iter(R, state, "add: commutativity", gr_test_add_commutativity, iters);
gr_test_iter(R, state, "add: aliasing", gr_test_add_aliasing, iters);
gr_test_iter(R, state, "inv: multiplication", gr_test_inv_multiplication, iters);
gr_test_iter(R, state, "inv: involution", gr_test_inv_involution, iters);
gr_test_iter(R, state, "pow_ui: exponent addition", gr_test_pow_ui_exponent_addition, iters);
gr_test_iter(R, state, "pow_ui: base scalar multiplication", gr_test_pow_ui_base_scalar_multiplication, iters);
if ((R-> flags & GR_COMMUTATIVE_RING) == GR_COMMUTATIVE_RING)
gr_test_iter(R, state, "pow_ui: base multiplication", gr_test_pow_ui_base_multiplication, iters);
gr_test_iter(R, state, "pow_ui: aliasing", gr_test_pow_ui_exponent_addition, iters);
gr_test_iter(R, state, "vec_add", gr_test_vec_add, iters);
gr_test_iter(R, state, "mat_mul_classical: associativity", gr_test_mat_mul_classical_associativity, iters);
flint_randclear(state);
}
int main()
{
GR_TMP_START;
gr_ctx_t Zn;
gr_ptr x, y, z;
nmod_t nm;
ulong t;
gr_ctx_init_nmod8(Zn, 107);
gr_ctx_t MZn;
gr_ctx_init_matrix(MZn, Zn, 30);
gr_test_ring(MZn, 10);
/*
{
flint_rand_t state;
flint_randinit(state);
gr_ptr a, b, c;
GR_TMP_START;
GR_TMP_INIT5(x, y, z, a, b, MZn);
GR_TMP_INIT1(c, MZn);
gr_randtest(x, state, NULL, MZn);
gr_randtest(y, state, NULL, MZn);
gr_randtest(a, state, NULL, MZn);
gr_randtest(b, state, NULL, MZn);
gr_add(z, x, y, MZn);
gr_mul(z, x, y, MZn);
gr_mul(c, a, b, MZn);
GR_TMP_CLEAR5(x, y, z, a, b, MZn);
GR_TMP_CLEAR1(c, MZn);
GR_TMP_END;
}
*/
gr_ctx_clear_matrix(MZn);
gr_ctx_clear_nmod8(Zn);
return 0;
{
gr_ptr A = flint_malloc(MZn->sizeof_elem);
gr_init(A, MZn);
gr_println(A, MZn);
gr_clear(A, MZn);
}
return 0;
gr_test_ring(Zn, 100000);
gr_mat_t A, B, C;
gr_mat_init(A, 13, 14, Zn);
gr_mat_init(B, 14, 17, Zn);
gr_mat_init(C, 13, 17, Zn);
flint_rand_t state;
flint_randinit(state);
gr_mat_randtest(A, state, NULL, Zn);
gr_mat_randtest(B, state, NULL, Zn);
gr_mat_randtest(C, state, NULL, Zn);
gr_mat_print(A, Zn);
gr_mat_print(B, Zn);
gr_mat_print(C, Zn);
TIMEIT_START
gr_mat_mul_classical(C, A, B, Zn);
TIMEIT_STOP
gr_mat_print(C, Zn);
printf("\n");
GR_TMP_INIT2(x, y, Zn);
gr_set_si(x, 2, Zn);
nmod_init(&nm, 253);
TIMEIT_START
t = nmod_pow_ui(2, 123456, nm);
TIMEIT_STOP
TIMEIT_START
gr_generic_pow_ui(y, x, 123456, Zn);
TIMEIT_STOP
flint_printf("%wu ", t);
gr_println(y, Zn);
GR_TMP_CLEAR2(x, y, Zn);
gr_ctx_clear_nmod8(Zn);
GR_TMP_END;
flint_cleanup();
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment