p-adic draft

padic.c

#include "padic.h"

#include "nmod.h"
#include "gr.h"

typedef struct
{
    ulong p;
    slong n;
    nmod_t p_mod;
    nmod_t pn_mod;
}
_padic_float_nmod_ctx_struct;

#define PADIC_FLOAT_NMOD_CTX(ctx) ((_padic_float_nmod_ctx_struct *)(GR_CTX_DATA_AS_PTR(ctx)))
#define PADIC_FLOAT_NMOD_CTX_P(ctx) (PADIC_FLOAT_NMOD_CTX(ctx)->p)
#define PADIC_FLOAT_NMOD_CTX_N(ctx) (PADIC_FLOAT_NMOD_CTX(ctx)->n)
#define PADIC_FLOAT_NMOD_CTX_P_MOD(ctx) (PADIC_FLOAT_NMOD_CTX(ctx)->p_mod)
#define PADIC_FLOAT_NMOD_CTX_PN_MOD(ctx) (PADIC_FLOAT_NMOD_CTX(ctx)->pn_mod)

typedef struct
{
    ulong man;
    slong exp;
}
padic_float_nmod_struct;

#define PADIC_EMIN (-(WORD_MAX / 4))
#define PADIC_EMAX ((WORD_MAX / 4))

typedef padic_float_nmod_struct padic_float_nmod_t[1];


/* Not checked: p is prime */
int
padic_float_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 */
    ulong pn = p;
    for (i = 1; i < n; i++)
    {
        umul_ppmm(hi, pn, pn, p);
        if (hi != 0)
            return GR_UNABLE;
    }

    flint_printf("pn = %wu\n", pn);

    GR_CTX_DATA_AS_PTR(ctx) = flint_malloc(sizeof(_padic_float_nmod_ctx_struct));

    PADIC_FLOAT_NMOD_CTX_P(ctx) = p;
    PADIC_FLOAT_NMOD_CTX_N(ctx) = n;

    nmod_init(&PADIC_FLOAT_NMOD_CTX_P_MOD(ctx), p);
    nmod_init(&PADIC_FLOAT_NMOD_CTX_PN_MOD(ctx), pn);

    return GR_SUCCESS;
}


void
padic_float_nmod_init(padic_float_nmod_t res, gr_ctx_t ctx)
{
    res->man = 0;
    res->exp = PADIC_EMAX;
}

void
padic_float_println(const padic_float_nmod_t x, gr_ctx_t ctx)
{
    if (x->man == 0)
    {
        flint_printf("0\n");
    }
    else
    {
        flint_printf("%wu * %wu^%wd\n", x->man, PADIC_FLOAT_NMOD_CTX_P(ctx), x->exp);
    }
}


int
padic_float_nmod_zero(padic_float_nmod_t res, gr_ctx_t ctx)
{
    res->man = 0;
    res->exp = PADIC_EMAX;
    return GR_SUCCESS;
}

int
padic_float_nmod_one(padic_float_nmod_t res, gr_ctx_t ctx)
{
    res->man = 1;
    res->exp = 0;
    return GR_SUCCESS;
}

truth_t
padic_float_nmod_is_zero(const padic_float_nmod_t x, gr_ctx_t ctx)
{
    return (x->man == 0) ? T_TRUE : T_FALSE;
}

truth_t
padic_float_nmod_is_one(const padic_float_nmod_t x, gr_ctx_t ctx)
{
    return (x->man == 1 && x->exp == 0) ? T_TRUE : T_FALSE;
}

int
padic_float_nmod_set_ui(padic_float_nmod_t res, ulong x, gr_ctx_t ctx)
{
    if (x < PADIC_FLOAT_NMOD_CTX_P(ctx))
    {
        res->man = x;
        res->exp = (x == 0) ? PADIC_EMAX : 0;
    }
    else
    {
        slong exp = 0;
        ulong r, q;
        ulong p, pinv;

        p = PADIC_FLOAT_NMOD_CTX_P(ctx);
        pinv = PADIC_FLOAT_NMOD_CTX_P_MOD(ctx).ninv;

        while (1)
        {
            r = n_divrem2_preinv(&q, x, p, pinv);
            if (r == 0)
            {
                x = q;
                exp++;
            }
            else
            {
                break;
            }
        }

        x = nmod_set_ui(x, PADIC_FLOAT_NMOD_CTX_PN_MOD(ctx));

        res->man = x;
        res->exp = exp;
    }

    return GR_SUCCESS;
}

int
padic_float_nmod_mul(padic_float_nmod_t res, const padic_float_nmod_t a, const padic_float_nmod_t b, gr_ctx_t ctx)
{
    if (a->man == 0 || b->man == 0)
        return padic_float_nmod_zero(res, ctx);

    res->man = nmod_mul(a->man, b->man, PADIC_FLOAT_NMOD_CTX_PN_MOD(ctx));
    res->exp = a->exp + b->exp;

    /* Underflow or overflow */
    if (res->exp < PADIC_EMIN || res->exp > PADIC_EMAX)
        return GR_UNABLE;
    else
        return GR_SUCCESS;
}

#include "profiler.h"

int main()
{
    {
        gr_ctx_t ctx;

        padic_float_nmod_t a, b, c;

        padic_float_nmod_ctx_init(ctx, 13, 4);

        padic_float_nmod_init(a, ctx);
        padic_float_nmod_init(b, ctx);
        padic_float_nmod_init(c, ctx);

        padic_float_nmod_set_ui(a, 13 * 13 * 12345, ctx);
        padic_float_nmod_set_ui(b, 13 * 13 * 13 * 13 * UWORD(76543), ctx);

        TIMEIT_START
        padic_float_nmod_mul(c, a, b, ctx);
        TIMEIT_STOP

        padic_float_println(a, ctx);
        padic_float_println(b, ctx);
        padic_float_println(c, ctx);

        flint_printf("\n");
    }

    {
        padic_ctx_t ctx;
        fmpz_t p;
        fmpz_init(p);
        fmpz_set_ui(p, 13);

        padic_t a, b, c;

        padic_ctx_init(ctx, p, 0, 100, PADIC_VAL_UNIT);

        padic_init2(a, 6);
        padic_init2(b, 8);
        padic_init2(c, 10);

        padic_set_ui(a, 13 * 13 * 12345, ctx);
        padic_set_ui(b, 13 * 13 * 13 * 13 * UWORD(76543), ctx);

        TIMEIT_START
        padic_mul(c, a, b, ctx);
        TIMEIT_STOP

        padic_print(a, ctx); flint_printf("\n");
        padic_print(b, ctx); flint_printf("\n");
        padic_print(c, ctx); flint_printf("\n");

    }
}