flint3-sys 3.6.0

/*
    Copyright (C) 2026 Fredrik Johansson

    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 <https://www.gnu.org/licenses/>.
*/

#include "radix.h"
#include "fmpz.h"
#include "gr.h"

void
radix_integer_init(radix_integer_t res, const radix_t radix)
{
    res->d = NULL;
    res->size = 0;
    res->alloc = 0;
}

void
radix_integer_clear(radix_integer_t res, const radix_t radix)
{
    if (res->d != NULL)
        flint_free(res->d);
}

nn_ptr
radix_integer_fit_limbs(radix_integer_t res, slong nlimbs, const radix_t radix)
{
    if (res->alloc < nlimbs)
    {
        res->d = flint_realloc(res->d, nlimbs * sizeof(ulong));
        res->alloc = nlimbs;
    }

    return res->d;
}

void
radix_integer_zero(radix_integer_t res, const radix_t radix)
{
    /* todo: do we want to deallocate? */
    res->size = 0;
}

void
radix_integer_rand_limbs(radix_integer_t res, flint_rand_t state, slong n, const radix_t radix)
{
    nn_ptr d = radix_integer_fit_limbs(res, n, radix);
    radix_rand_limbs(d, state, n, radix);
    MPN_NORM(d, n);
    res->size = (n_randlimb(state) & 1) ? n : -n;
}

void
radix_integer_randtest_limbs(radix_integer_t res, flint_rand_t state, slong max_limbs, const radix_t radix)
{
    slong n = n_randint(state, max_limbs + 1);
    nn_ptr d = radix_integer_fit_limbs(res, n, radix);
    radix_randtest_limbs(d, state, n, radix);
    MPN_NORM(d, n);
    res->size = (n_randlimb(state) & 1) ? n : -n;
}

void
radix_integer_one(radix_integer_t res, const radix_t radix)
{
    radix_integer_fit_limbs(res, 1, radix)[0] = 1;
    res->size = 1;
}

void
radix_integer_neg_one(radix_integer_t res, const radix_t radix)
{
    radix_integer_fit_limbs(res, 1, radix)[0] = 1;
    res->size = -1;
}

int
radix_integer_is_zero(const radix_integer_t x, const radix_t radix)
{
    return x->size == 0;
}

int
radix_integer_is_one(const radix_integer_t x, const radix_t radix)
{
    return x->size == 1 && x->d[0] == 1;
}

int
radix_integer_is_neg_one(const radix_integer_t x, const radix_t radix)
{
    return x->size == -1 && x->d[0] == 1;
}

int
radix_integer_equal(const radix_integer_t x, const radix_integer_t y, const radix_t radix)
{
    return (x->size == y->size) && flint_mpn_equal_p(x->d, y->d, FLINT_ABS(x->size));
}

int
radix_integer_cmp(const radix_integer_t x, const radix_integer_t y, const radix_t radix)
{
    slong xsize = x->size;
    slong ysize = y->size;
    if (xsize != ysize)
        return (xsize > ysize) ? 1 : -1;
    int cmp = mpn_cmp(x->d, y->d, FLINT_ABS(xsize));
    return (xsize >= 0) ? FLINT_SGN(cmp) : -FLINT_SGN(cmp);
}

int
radix_integer_cmpabs(const radix_integer_t x, const radix_integer_t y, const radix_t radix)
{
    slong xsize = FLINT_ABS(x->size);
    slong ysize = FLINT_ABS(y->size);
    if (xsize != ysize)
        return (xsize > ysize) ? 1 : -1;
    int cmp = mpn_cmp(x->d, y->d, xsize);
    return FLINT_SGN(cmp);
}

void
radix_integer_set(radix_integer_t res, const radix_integer_t x, const radix_t radix)
{
    if (res != x)
    {
        slong size = x->size;
        slong n = FLINT_ABS(size);

        if (n != 0)
        {
            radix_integer_fit_limbs(res, n, radix);
            flint_mpn_copyi(res->d, x->d, n);
        }

        res->size = size;
    }
}

void
radix_integer_neg(radix_integer_t res, const radix_integer_t x, const radix_t radix)
{
    if (res != x)
        radix_integer_set(res, x, radix);

    res->size = -res->size;
}


void
radix_integer_abs(radix_integer_t res, const radix_integer_t x, const radix_t radix)
{
    if (res != x)
        radix_integer_set(res, x, radix);

    res->size = FLINT_ABS(res->size);
}

int
radix_integer_sgn(const radix_integer_t x, const radix_t radix)
{
    return FLINT_SGN(x->size);
}

void
radix_integer_set_ui(radix_integer_t res, ulong x, const radix_t radix)
{
    if (x == 0)
    {
        radix_integer_zero(res, radix);
    }
    else if (x < LIMB_RADIX(radix))
    {
        radix_integer_fit_limbs(res, 1, radix)[0] = x;
        res->size = 1;
    }
    else
    {
        res->size = radix_set_mpn(radix_integer_fit_limbs(res,
            radix_set_mpn_need_alloc(1, radix), radix), &x, 1, radix);
    }
}

void
radix_integer_set_si(radix_integer_t res, slong x, const radix_t radix)
{
    ulong ux = FLINT_UABS(x);

    if (x == 0)
    {
        radix_integer_zero(res, radix);
    }
    else if (ux < LIMB_RADIX(radix))
    {
        radix_integer_fit_limbs(res, 1, radix)[0] = ux;
        res->size = (x < 0) ? -1 : 1;
    }
    else
    {
        nn_ptr rd = radix_integer_fit_limbs(res, radix_set_mpn_need_alloc(1, radix), radix);
        slong rn = radix_set_mpn(rd, &ux, 1, radix);
        res->size = (x < 0) ? -rn : rn;
    }
}

void
radix_integer_set_fmpz(radix_integer_t res, const fmpz_t x, const radix_t radix)
{
    if (!COEFF_IS_MPZ(*x))
    {
        radix_integer_set_si(res, *x, radix);
    }
    else
    {
        slong xsize = COEFF_TO_PTR(*x)->_mp_size;
        nn_srcptr xd = COEFF_TO_PTR(*x)->_mp_d;
        nn_ptr rd = radix_integer_fit_limbs(res, radix_set_mpn_need_alloc(FLINT_ABS(xsize), radix), radix);
        slong rn = radix_set_mpn(rd, xd, FLINT_ABS(xsize), radix);
        res->size = (xsize > 0) ? rn : -rn;
    }
}

void
radix_integer_get_fmpz(fmpz_t res, const radix_integer_t x, const radix_t radix)
{
    slong xsize = x->size;

    if (xsize == 0)
        fmpz_zero(res);
    else if (xsize == 1)
        fmpz_set_ui(res, x->d[0]);
    else if (xsize == -1)
        fmpz_neg_ui(res, x->d[0]);
    else
    {
        __mpz_struct * r = _fmpz_promote(res);
        slong xn = FLINT_ABS(xsize);
        nn_srcptr xd = x->d;
        nn_ptr rd = FLINT_MPZ_REALLOC(r, xn);
        slong rn = radix_get_mpn(rd, xd, xn, radix);
        /* todo: check for possible int overflow */
        r->_mp_size = (xsize > 0) ? rn : -rn;
        _fmpz_demote_val(res);
    }
}

/* todo: 1x1 and nx1 specializations */
/* todo: maybe handle aliasing in radix_mul to simplify? */
void
radix_integer_mul(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, const radix_t radix)
{
    slong xsize = x->size;
    slong ysize = y->size;

    if (xsize == 0 || ysize == 0)
    {
        radix_integer_zero(res, radix);
        return;
    }

    slong xn = FLINT_ABS(xsize);
    slong yn = FLINT_ABS(ysize);
    slong rn = xn + yn;
    slong sgn = xsize ^ ysize;
    nn_ptr rd = radix_integer_fit_limbs(res, rn, radix);
    nn_srcptr xd = x->d;
    nn_srcptr yd = y->d;

    /* todo: squaring check first? */

    if (res == x)
    {
        TMP_INIT;
        TMP_START;
        nn_ptr tmp = TMP_ALLOC(sizeof(ulong) * xn);
        flint_mpn_copyi(tmp, xd, xn);
        if (x == y)
            radix_sqr(rd, tmp, xn, radix);
        else
            radix_mul(rd, tmp, xn, yd, yn, radix);
        TMP_END;
    }
    else if (res == y)
    {
        TMP_INIT;
        TMP_START;
        nn_ptr tmp = TMP_ALLOC(sizeof(ulong) * yn);
        flint_mpn_copyi(tmp, yd, yn);
        radix_mul(rd, xd, xn, tmp, yn, radix);
        TMP_END;
    }
    else
    {
        if (x == y)
            radix_sqr(rd, xd, xn, radix);
        else
            radix_mul(rd, xd, xn, yd, yn, radix);
    }

    rn -= (rd[rn - 1] == 0);
    res->size = (sgn >= 0) ? rn : -rn;
}

void
radix_integer_add(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, const radix_t radix)
{
    slong xsize = x->size;
    slong ysize = y->size;

    if (xsize == 0)
    {
        radix_integer_set(res, y, radix);
        return;
    }

    if (ysize == 0)
    {
        radix_integer_set(res, x, radix);
        return;
    }

    slong xn = FLINT_ABS(xsize);
    slong yn = FLINT_ABS(ysize);
    int xsgnbit = (xsize < 0);
    int ysgnbit = (ysize < 0);
    slong rn = FLINT_MAX(xn, yn) + (xsgnbit == ysgnbit);
    nn_ptr rd = radix_integer_fit_limbs(res, rn, radix);
    nn_srcptr xd = x->d;
    nn_srcptr yd = y->d;
    int rsgnbit = 0;
    ulong cy;

    if (xsgnbit == ysgnbit)
    {
        rsgnbit = xsgnbit;
        if (xn >= yn)
            rd[rn - 1] = cy = radix_add(rd, xd, xn, yd, yn, radix);
        else
            rd[rn - 1] = cy = radix_add(rd, yd, yn, xd, xn, radix);
        rn -= (cy == 0);
    }
    else
    {
        if (xn > yn || (xn == yn && mpn_cmp(xd, yd, xn) >= 0))
        {
            radix_sub(rd, xd, xn, yd, yn, radix);
            rsgnbit = xsgnbit;
        }
        else
        {
            radix_sub(rd, yd, yn, xd, xn, radix);
            rsgnbit = ysgnbit;
        }
        MPN_NORM(rd, rn);
    }

    res->size = rsgnbit ? -rn : rn;
}

void
radix_integer_sub(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, const radix_t radix)
{
    slong xsize = x->size;
    slong ysize = y->size;

    if (xsize == 0)
    {
        radix_integer_neg(res, y, radix);
        return;
    }

    if (ysize == 0)
    {
        radix_integer_set(res, x, radix);
        return;
    }

    slong xn = FLINT_ABS(xsize);
    slong yn = FLINT_ABS(ysize);
    int xsgnbit = (xsize < 0);
    int ysgnbit = (ysize > 0);
    slong rn = FLINT_MAX(xn, yn) + (xsgnbit == ysgnbit);
    nn_ptr rd = radix_integer_fit_limbs(res, rn, radix);
    nn_srcptr xd = x->d;
    nn_srcptr yd = y->d;
    int rsgnbit = 0;
    ulong cy;

    if (xsgnbit == ysgnbit)
    {
        rsgnbit = xsgnbit;
        if (xn >= yn)
            rd[rn - 1] = cy = radix_add(rd, xd, xn, yd, yn, radix);
        else
            rd[rn - 1] = cy = radix_add(rd, yd, yn, xd, xn, radix);
        rn -= (cy == 0);
    }
    else
    {
        if (xn > yn || (xn == yn && mpn_cmp(xd, yd, xn) >= 0))
        {
            radix_sub(rd, xd, xn, yd, yn, radix);
            rsgnbit = xsgnbit;
        }
        else
        {
            radix_sub(rd, yd, yn, xd, xn, radix);
            rsgnbit = ysgnbit;
        }
        MPN_NORM(rd, rn);
    }

    res->size = rsgnbit ? -rn : rn;
}

/*
    Fused integer add/sub with a left shift of the second operand:
        radix_integer_addlsh: res = x + (y << v)
        radix_integer_sublsh: res = x - (y << v)
    for any shift v >= 0. The shift is decomposed into a limb offset nl = v / e
    and a within-limb digit shift nd = v % e (e = radix->exp); the nd == 0 case is
    a pure limb offset handled with radix_add / radix_sub on offset pointers, while
    nd != 0 uses the fused radix_addlsh / radix_sublsh / radix_lshsub primitives.

    No temporary integer is allocated for the differing-alignment work: when res
    does not alias x we materialise (y << v) directly in res and accumulate x into
    it; when res aliases x we operate in place with the fused primitives (the
    classic acc += shifted pattern). The only exception is the degenerate
    res == x == y aliasing, where one operand must be preserved in a small copy.

    The magnitude helpers below take normalised limb arrays (no high zero limbs)
    and return the normalised result length; dst may alias a, but never b.
*/

/* dst = |a| + (|b| << (nl*e + nd)); returns limb count. dst may alias a. */
static slong
_radix_imag_addlsh(nn_ptr dst, nn_srcptr a, slong an, nn_srcptr b, slong bn,
    slong nl, unsigned int nd, const radix_t radix)
{
    slong rn, top;
    ulong cy;

    if (an <= nl)
    {
        /* no overlap: a, then zeros up to nl, then b << nd */
        if (dst != a)
            flint_mpn_copyi(dst, a, an);
        flint_mpn_zero(dst + an, nl - an);
        if (nd == 0)
        {
            flint_mpn_copyi(dst + nl, b, bn);
            rn = nl + bn;
        }
        else
        {
            dst[nl + bn] = radix_lshift_digits(dst + nl, b, bn, nd, radix);
            rn = nl + bn + 1;
        }
    }
    else
    {
        slong am = an - nl;
        if (dst != a)
            flint_mpn_copyi(dst, a, nl);
        if (nd == 0)
        {
            if (am >= bn)
            {
                cy = radix_add(dst + nl, a + nl, am, b, bn, radix);
                top = nl + am;
            }
            else
            {
                cy = radix_add(dst + nl, b, bn, a + nl, am, radix);
                top = nl + bn;
            }
        }
        else
        {
            cy = radix_addlsh(dst + nl, a + nl, am, b, bn, nd, radix);
            top = nl + FLINT_MAX(am, bn + 1);
        }
        dst[top] = cy;
        rn = top + 1;
    }

    MPN_NORM(dst, rn);
    return rn;
}

/* dst = |a| - (|b| << (nl*e+nd)), assuming |a| >= |b<<..|; dst may alias a. */
static slong
_radix_imag_sublsh(nn_ptr dst, nn_srcptr a, slong an, nn_srcptr b, slong bn,
    slong nl, unsigned int nd, const radix_t radix)
{
    slong am = an - nl, top;

    if (dst != a)
        flint_mpn_copyi(dst, a, nl);

    if (nd == 0)
    {
        radix_sub(dst + nl, a + nl, am, b, bn, radix);     /* am >= bn */
        top = nl + am;
    }
    else
    {
        radix_sublsh(dst + nl, a + nl, am, b, bn, nd, radix);
        top = nl + FLINT_MAX(am, bn + 1);
    }

    MPN_NORM(dst, top);
    return top;
}

/* dst = (|b| << (nl*e+nd)) - |a|, assuming |b<<..| >= |a|; dst may alias a. */
static slong
_radix_imag_lshsub(nn_ptr dst, nn_srcptr b, slong bn, slong nl, unsigned int nd,
    nn_srcptr a, slong an, const radix_t radix)
{
    ulong B = LIMB_RADIX(radix);
    ulong one = 1;
    ulong b0;
    slong i, top;

    if (an <= nl)
    {
        /* a lies entirely below the shifted b; the low part is -a, which leaves a
           borrow b0 to take out of (b << nd) sitting at limb nl. */
        b0 = radix_neg(dst, a, an, radix);
        for (i = an; i < nl; i++)
            dst[i] = b0 ? (B - 1) : 0;

        if (nd == 0)
        {
            if (b0)
                radix_sub(dst + nl, b, bn, &one, 1, radix);
            else
                flint_mpn_copyi(dst + nl, b, bn);
            top = nl + bn;
        }
        else
        {
            dst[nl + bn] = radix_lshift_digits(dst + nl, b, bn, nd, radix);
            if (b0)
                radix_sub(dst + nl, dst + nl, bn + 1, &one, 1, radix);
            top = nl + bn + 1;
        }
    }
    else
    {
        slong am = an - nl;
        b0 = radix_neg(dst, a, nl, radix);             /* low nl limbs negated */
        if (nd == 0)
        {
            radix_sub(dst + nl, b, bn, a + nl, am, radix);    /* bn >= am */
            top = nl + bn;
        }
        else
        {
            radix_lshsub(dst + nl, b, bn, nd, a + nl, am, radix);
            top = nl + FLINT_MAX(am, bn + 1);
        }
        if (b0)
            radix_sub(dst + nl, dst + nl, top - nl, &one, 1, radix);
    }

    MPN_NORM(dst, top);
    return top;
}

/* sign of |a| - (|b| << (nl*e+nd)); a, b assumed normalised. */
static int
_radix_icmp_lsh(nn_srcptr a, slong an, nn_srcptr b, slong bn, slong nl,
    unsigned int nd, const radix_t radix)
{
    slong bvn, i;
    int c;

    if (nd == 0)
        bvn = nl + bn;
    else
    {
        ulong be2 = radix->bpow[radix->exp - nd];
        bvn = nl + bn + (b[bn - 1] >= be2);
    }

    if (an != bvn)
        return (an < bvn) ? -1 : 1;

    if (nd == 0)
        c = mpn_cmp(a + nl, b, bn);
    else
        c = radix_cmplsh(a + nl, an - nl, b, bn, nd, radix);
    if (c != 0)
        return c;

    /* high parts equal; any nonzero low limb of a makes |a| the larger */
    for (i = 0; i < nl; i++)
        if (a[i] != 0)
            return 1;
    return 0;
}

static void
_radix_integer_addlsh_core(radix_integer_t res, const radix_integer_t x,
    const radix_integer_t y, slong v, int ysub, const radix_t radix)
{
    slong xsize = x->size, ysize = y->size;
    slong e = radix->exp;
    slong xn, yn, nl, rn;
    unsigned int nd;
    int sx, sy, c, alias_xy;
    nn_ptr rd;
    nn_srcptr xd, yd;
    radix_integer_t tmp;

    if (ysize == 0)                         /* x +/- 0 */
    {
        radix_integer_set(res, x, radix);
        return;
    }
    if (xsize == 0)                         /* 0 +/- (y<<v) */
    {
        radix_integer_lshift_digits(res, y, v, radix);
        if (ysub && res->size != 0)
            res->size = -res->size;
        return;
    }
    if (v == 0)
    {
        if (ysub) radix_integer_sub(res, x, y, radix);
        else      radix_integer_add(res, x, y, radix);
        return;
    }

    /* res distinct from x: materialise (y<<v) in res, then combine x (no temp). */
    if (res != x)
    {
        radix_integer_lshift_digits(res, y, v, radix);     /* in place if res==y */
        if (ysub)
            radix_integer_sub(res, x, res, radix);         /* x - (y<<v) */
        else
            radix_integer_add(res, res, x, radix);         /* (y<<v) + x */
        return;
    }

    /* res aliases x: fused in-place combine. */
    xn = FLINT_ABS(xsize);
    yn = FLINT_ABS(ysize);
    sx = (xsize < 0);
    sy = (ysize < 0) ^ (ysub != 0);

    nl = (ulong) v / (ulong) e;
    nd = (ulong) v % (ulong) e;

    /* degenerate res == x == y: y must be preserved in a copy */
    alias_xy = (x->d == y->d);
    if (alias_xy)
    {
        radix_integer_init(tmp, radix);
        radix_integer_set(tmp, y, radix);
    }

    rd = radix_integer_fit_limbs(res, FLINT_MAX(xn, nl + yn + 1) + 1, radix);
    xd = x->d;                              /* == rd, data preserved by fit */
    yd = alias_xy ? tmp->d : y->d;

    if (sx == sy)
    {
        rn = _radix_imag_addlsh(rd, xd, xn, yd, yn, nl, nd, radix);
        res->size = sx ? -rn : rn;
    }
    else
    {
        c = _radix_icmp_lsh(xd, xn, yd, yn, nl, nd, radix);
        if (c == 0)
            res->size = 0;
        else if (c > 0)
        {
            rn = _radix_imag_sublsh(rd, xd, xn, yd, yn, nl, nd, radix);
            res->size = sx ? -rn : rn;
        }
        else
        {
            rn = _radix_imag_lshsub(rd, yd, yn, nl, nd, xd, xn, radix);
            res->size = sy ? -rn : rn;
        }
    }

    if (alias_xy)
        radix_integer_clear(tmp, radix);
}

void
radix_integer_addlsh(radix_integer_t res, const radix_integer_t x,
    const radix_integer_t y, slong v, const radix_t radix)
{
    _radix_integer_addlsh_core(res, x, y, v, 0, radix);
}

void
radix_integer_sublsh(radix_integer_t res, const radix_integer_t x,
    const radix_integer_t y, slong v, const radix_t radix)
{
    _radix_integer_addlsh_core(res, x, y, v, 1, radix);
}

int
radix_integer_is_normalised(const radix_integer_t x, const radix_t radix)
{
    slong xn = FLINT_ABS(x->size);
    slong i;

    for (i = 0; i < xn; i++)
        if (x->d[i] >= LIMB_RADIX(radix))
            return 0;

    if (xn != 0 && x->d[xn - 1] == 0)
        return 0;

    return 1;
}

void
radix_integer_set_limb(radix_integer_t res, const radix_integer_t x, slong index, ulong c, const radix_t radix)
{
    FLINT_ASSERT(index >= 0);
    FLINT_ASSERT(c < LIMB_RADIX(radix));

    slong xsize = x->size;
    slong xn = FLINT_ABS(xsize);
    slong rn = xn;

    if (index < xn - 1 || (c != 0 && index == xn - 1))
    {
        radix_integer_set(res, x, radix);
        res->d[index] = c;
    }
    else if (c == 0)
    {
        if (index == xn - 1)
        {
            rn = xn - 1;
            MPN_NORM(x->d, rn);
            if (res != x && rn != 0)
                flint_mpn_copyi(radix_integer_fit_limbs(res, rn, radix), x->d, rn);
            res->size = (xsize >= 0) ? rn : -rn;
        }
        else
        {
            radix_integer_set(res, x, radix);
        }
    }
    else
    {
        rn = index + 1;
        nn_ptr rd = radix_integer_fit_limbs(res, rn, radix);
        if (res != x)
            flint_mpn_copyi(rd, x->d, xn);
        flint_mpn_zero(rd + xn, index - xn);
        rd[index] = c;
        res->size = (xsize >= 0) ? rn : -rn;
    }
}

ulong
radix_integer_get_digit(const radix_integer_t x, slong index, const radix_t radix)
{
    ulong d;
    slong il, id;

    /* todo: speed up this division */
    il = (ulong) index / radix->exp;
    id = (ulong) index % radix->exp;

    d = radix_integer_get_limb(x, il, radix);

    if (d == 0)
    {
        return 0;
    }
    else
    {
        d = n_div_precomp_bounded(d, radix->bpow_div + id);
        n_divrem_precomp_bounded(&d, d, radix->bpow[1], radix->bpow_div + 1);
        return d;
    }
}

void
radix_integer_set_digit(radix_integer_t res, const radix_integer_t x, slong index, ulong c, const radix_t radix)
{
    ulong d, hi, lo;
    slong il, id;

    /* todo: speed up this division */
    il = (ulong) index / radix->exp;
    id = (ulong) index % radix->exp;

    d = radix_integer_get_limb(x, il, radix);

    if (d == 0)
    {
        d = c * radix->bpow[id];
    }
    else
    {
        hi = n_divrem_precomp_bounded(&lo, d, radix->bpow[id], radix->bpow_div + id);
        hi = n_div_precomp_bounded(hi, radix->bpow_div + 1);
        d = hi * (radix->bpow[id] * radix->bpow[1]) + c * radix->bpow[id] + lo;
    }

    radix_integer_set_limb(res, x, il, d, radix);
}

void
radix_integer_lshift_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    FLINT_ASSERT(n >= 0);

    slong xsize = x->size;

    if (xsize == 0)
    {
        radix_integer_zero(res, radix);
        return;
    }

    slong xn = FLINT_ABS(xsize);
    slong rn = xn + n;

    nn_ptr rd = radix_integer_fit_limbs(res, rn, radix);
    flint_mpn_copyd(rd + n, x->d, xn);
    flint_mpn_zero(rd, n);
    res->size = (xsize > 0) ? rn : -rn;
}

void
radix_integer_lshift_digits(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    FLINT_ASSERT(n >= 0);

    slong xsize = x->size;

    if (xsize == 0)
    {
        radix_integer_zero(res, radix);
        return;
    }

    slong nl, nd;

    /* todo: speed up this division */
    nl = (ulong) n / radix->exp;
    nd = (ulong) n % radix->exp;

    slong xn = FLINT_ABS(xsize);
    slong rn = xn + nl + (nd != 0);

    nn_ptr rd = radix_integer_fit_limbs(res, rn, radix);

    if (res != x || nl != 0)
        flint_mpn_copyd(rd + nl, x->d, xn);
    flint_mpn_zero(rd, nl);

    if (nd != 0)
    {
        rd[rn - 1] = radix_lshift_digits(rd + nl, rd + nl, xn, nd, radix);
        rn -= (rd[rn - 1] == 0);
    }

    res->size = (xsize > 0) ? rn : -rn;
}

void radix_integer_rshift_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    FLINT_ASSERT(n >= 0);

    slong xsize = x->size;
    slong xn = FLINT_ABS(xsize);
    slong rn = xn - n;

    if (rn <= 0)
    {
        radix_integer_zero(res, radix);
        return;
    }

    nn_ptr rd;
    if (res == x)
        rd = res->d;
    else
        rd = radix_integer_fit_limbs(res, rn, radix);

    flint_mpn_copyi(rd, x->d + n, rn);
    res->size = (xsize > 0) ? rn : -rn;
}

void radix_integer_rshift_digits(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    FLINT_ASSERT(n >= 0);

    slong xsize = x->size;
    slong xn = FLINT_ABS(xsize);

    slong nl, nd;

    /* todo: speed up this division */
    nl = (ulong) n / radix->exp;
    nd = (ulong) n % radix->exp;

    slong rn = xn - nl;

    if (rn <= 0)
    {
        radix_integer_zero(res, radix);
        return;
    }

    nn_ptr rd;
    if (res == x)
        rd = res->d;
    else
        rd = radix_integer_fit_limbs(res, rn, radix);

    if (res != x || nl != 0)
        flint_mpn_copyi(rd, x->d + nl, rn);

    if (nd != 0)
    {
        radix_rshift_digits(rd, rd, rn, nd, radix);
        rn -= (rd[rn - 1] == 0);
    }

    res->size = (xsize > 0) ? rn : -rn;
}

void
radix_integer_trunc_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    slong xsize = x->size;
    slong xn = FLINT_ABS(xsize);

    FLINT_ASSERT(n >= 0);

    if (res == x)
    {
        if (n < xn)
        {
            MPN_NORM(res->d, n);
            res->size = (xsize >= 0) ? n : -n;
        }
    }
    else
    {
        if (n < xn)
            MPN_NORM(x->d, n);
        else
            n = xn;

        flint_mpn_copyi(radix_integer_fit_limbs(res, n, radix), x->d, n);
        res->size = (xsize >= 0) ? n : -n;
    }
}

/* todo: optimize */
void
radix_integer_mod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    radix_integer_trunc_limbs(res, x, n, radix);

    if (res->size < 0)
    {
        slong rn = FLINT_ABS(res->size);
        radix_neg(res->d, res->d, rn, radix);

        if (rn < n)
        {
            slong i;
            radix_integer_fit_limbs(res, n, radix);
            ulong B1 = LIMB_RADIX(radix) - 1;
            for (i = rn; i < n; i++)
                res->d[i] = B1;
        }
        MPN_NORM(res->d, n);
        res->size = n;
    }
}

void
radix_integer_smod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    slong xsize = x->size;
    slong xn = FLINT_ABS(xsize);
    slong rn;

    FLINT_ASSERT(n >= 0);
    rn = FLINT_MIN(n, xn);
    MPN_NORM(x->d, rn);

    if (rn > 0 && rn == n)
    {
        int cmp = radix_cmp_bn_half(x->d, n, radix);

        if (cmp > 0 || (cmp == 0 && xsize > 0 && (LIMB_RADIX(radix) % 2 == 0)))
        {
            if (res == x)
                radix_neg(res->d, res->d, rn, radix);
            else
                radix_neg(radix_integer_fit_limbs(res, n, radix), x->d, rn, radix);

            MPN_NORM(res->d, rn);
            res->size = (xsize > 0) ? -rn : rn;
            return;
        }
    }

    if (res != x)
        flint_mpn_copyi(radix_integer_fit_limbs(res, rn, radix), x->d, rn);
    res->size = (xsize >= 0) ? rn : -rn;
}

void
radix_integer_trunc_digits(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    slong e = radix->exp;
    slong limbs, rem, tot;

    if (n <= 0)
    {
        radix_integer_zero(res, radix);
        return;
    }

    limbs = n / e;
    rem = n - limbs * e;
    tot = limbs + (rem != 0);

    radix_integer_trunc_limbs(res, x, tot, radix);

    if (rem != 0 && FLINT_ABS(res->size) == tot)
    {
        ulong top = res->d[limbs];
        ulong m = radix->bpow[rem];
        ulong t = n_rem_precomp(top, m, radix->bpow_div + rem);

        res->d[limbs] = t;

        while (tot > 0 && res->d[tot - 1] == 0)
            tot--;
        res->size = (res->size >= 0) ?  tot : -tot;
    }
}

void
radix_integer_mod_digits(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    slong e = radix->exp;
    slong limbs, rem, tot;

    if (n <= 0)
    {
        radix_integer_zero(res, radix);
        return;
    }

    limbs = n / e;
    rem = n - limbs * e;
    tot = limbs + (rem != 0);

    radix_integer_mod_limbs(res, x, tot, radix);

    if (rem != 0 && FLINT_ABS(res->size) == tot)
    {
        ulong top = res->d[limbs];
        ulong m = radix->bpow[rem];
        ulong t = n_rem_precomp(top, m, radix->bpow_div + rem);

        res->d[limbs] = t;

        while (tot > 0 && res->d[tot - 1] == 0)
            tot--;
        res->size = tot;
    }
}

void
radix_integer_mullow_limbs(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, slong n, const radix_t radix)
{
    slong xsize = x->size;
    slong ysize = y->size;

    if (xsize == 0 || ysize == 0 || n == 0)
    {
        radix_integer_zero(res, radix);
        return;
    }

    slong xn = FLINT_ABS(xsize);
    slong yn = FLINT_ABS(ysize);

    xn = FLINT_MIN(xn, n);
    yn = FLINT_MIN(yn, n);
    slong rn = FLINT_MIN(xn + yn, n);

    slong sgn = xsize ^ ysize;
    nn_ptr rd = radix_integer_fit_limbs(res, rn, radix);
    nn_srcptr xd = x->d;
    nn_srcptr yd = y->d;

    if (res == x)
    {
        TMP_INIT;
        TMP_START;
        nn_ptr tmp = TMP_ALLOC(sizeof(ulong) * xn);
        flint_mpn_copyi(tmp, xd, xn);
        if (x == y)
            radix_mulmid(rd, tmp, xn, tmp, yn, 0, rn, radix);
        else
            radix_mulmid(rd, tmp, xn, yd, yn, 0, rn, radix);
        TMP_END;
    }
    else if (res == y)
    {
        TMP_INIT;
        TMP_START;
        nn_ptr tmp = TMP_ALLOC(sizeof(ulong) * yn);
        flint_mpn_copyi(tmp, yd, yn);
        radix_mulmid(rd, xd, xn, tmp, yn, 0, rn, radix);
        TMP_END;
    }
    else
    {
        radix_mulmid(rd, xd, xn, yd, yn, 0, rn, radix);
    }

    MPN_NORM(rd, rn);
    res->size = (sgn >= 0) ? rn : -rn;
}

int
radix_integer_invmod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    if (n == 0)
    {
        radix_integer_zero(res, radix);
        return 1;
    }

    slong xsize = x->size;
    if (xsize == 0)
        return 0;

    slong xn = FLINT_ABS(xsize);
    xn = FLINT_MIN(xn, n);

    if (xn == 1 && x->d[0] == 1)
    {
        /* preserve sign of x */
        radix_integer_trunc_limbs(res, x, 1, radix);
        return 1;
    }

    slong rn = n;
    int invertible;
    nn_ptr rd = radix_integer_fit_limbs(res, n, radix);
    nn_srcptr xd = x->d;

    if (res == x && rn != 1)
    {
        TMP_INIT;
        TMP_START;
        nn_ptr tmp = TMP_ALLOC(sizeof(ulong) * xn);
        flint_mpn_copyi(tmp, xd, xn);
        invertible = radix_invmod_bn(rd, tmp, xn, n, radix);
        TMP_END;
    }
    else
    {
        invertible = radix_invmod_bn(rd, xd, xn, n, radix);
    }

    if (!invertible)
        rn = 0;

    MPN_NORM(rd, rn);
    res->size = (xsize > 0) ? rn : -rn;
    return invertible;
}

int
radix_integer_divmod_limbs(radix_integer_t res, const radix_integer_t a,
    const radix_integer_t b, slong n, const radix_t radix)
{
    slong asize = a->size;
    slong bsize = b->size;
    slong an, bn, rn, sgn;
    nn_ptr rd;
    int invertible;

    if (n == 0)
    {
        radix_integer_zero(res, radix);
        return 1;
    }

    if (bsize == 0)
    {
        radix_integer_zero(res, radix);
        return 0;
    }

    if (asize == 0)
    {
        ulong tmp = b->d[0];
        radix_integer_zero(res, radix);
        return radix_invmod_bn(&tmp, &tmp, 1, 1, radix);
    }

    an = FLINT_MIN(FLINT_ABS(asize), n);
    bn = FLINT_MIN(FLINT_ABS(bsize), n);
    sgn = asize ^ bsize;

    if (bn == 1 && b->d[0] == 1)
    {
        radix_integer_trunc_limbs(res, a, n, radix);
        res->size = (sgn >= 0) ? FLINT_ABS(res->size) : -FLINT_ABS(res->size);
        return 1;
    }

    rd = radix_integer_fit_limbs(res, n, radix);

    if (res == b)
    {
        TMP_INIT;
        TMP_START;
        nn_ptr tmp = TMP_ALLOC(sizeof(ulong) * bn);
        flint_mpn_copyi(tmp, b->d, bn);
        invertible = radix_divmod_bn(rd, NULL, a->d, an, tmp, bn, n, radix);
        TMP_END;
    }
    else
    {
        invertible = radix_divmod_bn(rd, NULL, a->d, an, b->d, bn, n, radix);
    }

    if (!invertible)
    {
        radix_integer_zero(res, radix);
        return 0;
    }

    rn = n;
    MPN_NORM(rd, rn);
    res->size = (rn == 0) ? 0 : ((sgn >= 0) ? rn : -rn);
    return 1;
}

/* res = a reciprocal square root of x modulo B^n, i.e. res^2 x == 1 (mod B^n).
   x must be a unit (x->d[0] coprime to p). Returns 1 on success, 0 if x is zero
   or not a square modulo p (in which case res is set to zero). */
int
radix_integer_rsqrtmod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    slong xsize, xn, rn;
    int ok;
    nn_ptr rd;
    nn_srcptr xd;

    if (n == 0)
    {
        radix_integer_zero(res, radix);
        return 1;
    }

    xsize = x->size;
    if (xsize == 0)
        return 0;                       /* zero has no reciprocal square root */

    xn = FLINT_ABS(xsize);
    xn = FLINT_MIN(xn, n);

    rn = n;
    rd = radix_integer_fit_limbs(res, n, radix);
    xd = x->d;

    if (res == x)
    {
        TMP_INIT;
        TMP_START;
        nn_ptr tmp = TMP_ALLOC(sizeof(ulong) * xn);
        flint_mpn_copyi(tmp, xd, xn);
        ok = radix_rsqrtmod_bn(rd, tmp, xn, n, radix);
        TMP_END;
    }
    else
    {
        ok = radix_rsqrtmod_bn(rd, xd, xn, n, radix);
    }

    if (!ok)
    {
        radix_integer_zero(res, radix);
        return 0;
    }

    MPN_NORM(rd, rn);
    res->size = rn;                     /* canonical nonnegative unit */
    return 1;
}

/* res = a square root of x modulo B^n, i.e. res^2 == x (mod B^n). x must be a
   unit (x->d[0] coprime to p). Returns 1 on success, 0 if x is zero or not a
   square modulo p (in which case res is set to zero). */
int
radix_integer_sqrtmod_limbs(radix_integer_t res, const radix_integer_t x, slong n, const radix_t radix)
{
    slong xsize, xn, rn;
    int ok;
    nn_ptr rd;
    nn_srcptr xd;

    if (n == 0)
    {
        radix_integer_zero(res, radix);
        return 1;
    }

    xsize = x->size;
    if (xsize == 0)
        return 0;                       /* zero handled by the caller's valuation */

    xn = FLINT_ABS(xsize);
    xn = FLINT_MIN(xn, n);

    rn = n;
    rd = radix_integer_fit_limbs(res, n, radix);
    xd = x->d;

    if (res == x)
    {
        TMP_INIT;
        TMP_START;
        nn_ptr tmp = TMP_ALLOC(sizeof(ulong) * xn);
        flint_mpn_copyi(tmp, xd, xn);
        ok = radix_sqrtmod_bn(rd, tmp, xn, n, radix);
        TMP_END;
    }
    else
    {
        ok = radix_sqrtmod_bn(rd, xd, xn, n, radix);
    }

    if (!ok)
    {
        radix_integer_zero(res, radix);
        return 0;
    }

    MPN_NORM(rd, rn);
    res->size = rn;                     /* canonical nonnegative root */
    return 1;
}

static ulong
radix_cdivrem(nn_ptr q, nn_ptr r, nn_srcptr a, slong an, nn_srcptr b, slong bn, const radix_t radix)
{
    ulong cy = 0;

    if (b == r)
    {
        nn_ptr tmp;
        TMP_INIT;
        TMP_START;
        tmp = TMP_ALLOC(bn * sizeof(ulong));
        flint_mpn_copyi(tmp, b, bn);
        cy = radix_cdivrem(q, r, a, an, tmp, bn, radix);
        TMP_END;
    }
    else
    {
        radix_divrem(q, r, a, an, b, bn, radix);

        if (!flint_mpn_zero_p(r, bn))
        {
            ulong one = 1;
            radix_sub(r, b, bn, r, bn, radix);
            cy = radix_add(q, q, an - bn + 1, &one, 1, radix);
        }
    }

    return cy;
}

int radix_integer_div(radix_integer_t q,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    slong asize, bsize, an, bn, qn;
    nn_srcptr ad, bd;
    nn_ptr qd;
    int exact;

    asize = a->size;
    bsize = b->size;

    an = FLINT_ABS(asize);
    bn = FLINT_ABS(bsize);

    if (bn == 0)
        return 0;

    if (an == 0)
    {
        radix_integer_zero(q, radix);
        return 1;
    }

    qn = an - bn + 1;
    ad = a->d;
    bd = b->d;

    if (an < bn || (an == bn && mpn_cmp(ad, bd, an) < 0))
    {
        radix_integer_zero(q, radix);
        return 0;
    }

    qd = radix_integer_fit_limbs(q, qn, radix);
    /* read after fit_limbs in case of possible aliasing */
    ad = a->d;
    bd = b->d;

    exact = radix_div(qd, ad, an, bd, bn, radix);
    if (!exact)
        qn = 0;
    MPN_NORM(qd, qn);
    q->size = ((asize < 0) == (bsize < 0)) ? qn : -qn;

    return exact;
}

void radix_integer_divexact(radix_integer_t q,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    slong asize, bsize, an, bn, qn;
    nn_srcptr ad, bd;
    nn_ptr qd;

    asize = a->size;
    bsize = b->size;

    an = FLINT_ABS(asize);
    bn = FLINT_ABS(bsize);

    if (an == 0)
    {
        radix_integer_zero(q, radix);
        return;
    }

    if (bn == 0 || an < bn)
        flint_throw(FLINT_ERROR, "radix_integer_divexact: an < bn");

    qn = an - bn + 1;
    qd = radix_integer_fit_limbs(q, qn, radix);
    /* read after fit_limbs in case of possible aliasing */
    ad = a->d;
    bd = b->d;

    radix_divexact(qd, ad, an, bd, bn, radix);
    MPN_NORM(qd, qn);
    q->size = ((asize < 0) == (bsize < 0)) ? qn : -qn;
}

void radix_integer_tdiv_qr(radix_integer_t q, radix_integer_t r,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    slong asize, bsize, an, bn, qn, rn;
    nn_srcptr ad, bd;
    nn_ptr qd, rd;

    asize = a->size;
    bsize = b->size;

    an = FLINT_ABS(asize);
    bn = FLINT_ABS(bsize);

    if (bn == 0)
        flint_throw(FLINT_DIVZERO, "radix_integer_tdiv_qr: divide by zero");

    if (an == 0)
    {
        radix_integer_zero(q, radix);
        radix_integer_zero(r, radix);
        return;
    }

    if (an < bn)
    {
        radix_integer_set(r, a, radix);
        radix_integer_zero(q, radix);
        return;
    }

    qn = an - bn + 1;
    rn = bn;
    qd = radix_integer_fit_limbs(q, qn, radix);
    rd = radix_integer_fit_limbs(r, rn, radix);
    /* read after fit_limbs in case of possible aliasing */
    ad = a->d;
    bd = b->d;

    radix_divrem(qd, rd, ad, an, bd, bn, radix);

    MPN_NORM(qd, qn);
    MPN_NORM(rd, rn);

    q->size = ((asize < 0) == (bsize < 0)) ? qn : -qn;
    r->size = (asize < 0) ? -rn : rn;
}

void radix_integer_tdiv_q(radix_integer_t q,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    slong asize, bsize, an, bn, qn;
    nn_srcptr ad, bd;
    nn_ptr qd;

    asize = a->size;
    bsize = b->size;

    an = FLINT_ABS(asize);
    bn = FLINT_ABS(bsize);

    if (bn == 0)
        flint_throw(FLINT_DIVZERO, "radix_integer_tdiv_qr: divide by zero");

    if (an == 0)
    {
        radix_integer_zero(q, radix);
        return;
    }

    if (an < bn)
    {
        radix_integer_zero(q, radix);
        return;
    }

    qn = an - bn + 1;
    qd = radix_integer_fit_limbs(q, qn, radix);
    /* read after fit_limbs in case of possible aliasing */
    ad = a->d;
    bd = b->d;
    radix_divrem(qd, NULL, ad, an, bd, bn, radix);

    MPN_NORM(qd, qn);

    q->size = ((asize < 0) == (bsize < 0)) ? qn : -qn;
}

void radix_integer_fdiv_qr(radix_integer_t q, radix_integer_t r,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    slong asize, bsize, an, bn, qn, rn;
    nn_srcptr ad, bd;
    nn_ptr qd, rd;

    asize = a->size;
    bsize = b->size;

    an = FLINT_ABS(asize);
    bn = FLINT_ABS(bsize);

    if (bn == 0)
        flint_throw(FLINT_DIVZERO, "radix_integer_fdiv_qr: divide by zero");

    if (an == 0)
    {
        radix_integer_zero(q, radix);
        radix_integer_zero(r, radix);
        return;
    }

    if (an < bn)
    {
        if ((asize < 0) == (bsize < 0))
        {
            radix_integer_set(r, a, radix);
            radix_integer_zero(q, radix);
        }
        else
        {
            radix_integer_add(r, a, b, radix);
            radix_integer_neg_one(q, radix);
        }
        return;
    }

    qn = an - bn + 1;
    rn = bn;

    qd = radix_integer_fit_limbs(q, qn, radix);
    rd = radix_integer_fit_limbs(r, rn, radix);
    /* read after fit_limbs in case of possible aliasing */
    ad = a->d;
    bd = b->d;

    if ((asize < 0) == (bsize < 0))
    {
        radix_divrem(qd, rd, ad, an, bd, bn, radix);
    }
    else
    {
        ulong cy = radix_cdivrem(qd, rd, ad, an, bd, bn, radix);

        if (cy != 0)
        {
            qd = radix_integer_fit_limbs(q, qn + 1, radix);
            qd[qn] = 1;
            qn++;
        }
    }

    MPN_NORM(qd, qn);
    MPN_NORM(rd, rn);

    q->size = ((asize < 0) == (bsize < 0)) ? qn : -qn;
    r->size = (bsize < 0) ? -rn : rn;
}

void radix_integer_cdiv_qr(radix_integer_t q, radix_integer_t r,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    slong asize, bsize, an, bn, qn, rn;
    nn_srcptr ad, bd;
    nn_ptr qd, rd;

    asize = a->size;
    bsize = b->size;

    an = FLINT_ABS(asize);
    bn = FLINT_ABS(bsize);

    if (bn == 0)
        flint_throw(FLINT_DIVZERO, "radix_integer_cdiv_qr: divide by zero");

    if (an == 0)
    {
        radix_integer_zero(q, radix);
        radix_integer_zero(r, radix);
        return;
    }

    if (an < bn)
    {
        if ((asize < 0) != (bsize < 0))
        {
            radix_integer_set(r, a, radix);
            radix_integer_zero(q, radix);
        }
        else
        {
            radix_integer_sub(r, a, b, radix);
            radix_integer_one(q, radix);
        }
        return;
    }

    qn = an - bn + 1;
    rn = bn;

    qd = radix_integer_fit_limbs(q, qn, radix);
    rd = radix_integer_fit_limbs(r, rn, radix);
    /* read after fit_limbs in case of possible aliasing */
    ad = a->d;
    bd = b->d;

    if ((asize < 0) != (bsize < 0))
    {
        radix_divrem(qd, rd, ad, an, bd, bn, radix);
    }
    else
    {
        ulong cy = radix_cdivrem(qd, rd, ad, an, bd, bn, radix);

        if (cy != 0)
        {
            qd = radix_integer_fit_limbs(q, qn + 1, radix);
            qd[qn] = 1;
            qn++;
        }
    }

    MPN_NORM(qd, qn);
    MPN_NORM(rd, rn);

    q->size = ((asize < 0) == (bsize < 0)) ? qn : -qn;
    r->size = (bsize < 0) ? rn : -rn;
}

/* Todo: optimize */

void
radix_integer_fdiv_q(radix_integer_t q,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    radix_integer_t r;
    radix_integer_init(r, radix);
    radix_integer_fdiv_qr(q, r, a, b, radix);
    radix_integer_clear(r, radix);
}

void
radix_integer_cdiv_q(radix_integer_t q,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    radix_integer_t r;
    radix_integer_init(r, radix);
    radix_integer_cdiv_qr(q, r, a, b, radix);
    radix_integer_clear(r, radix);
}

void
radix_integer_tdiv_r(radix_integer_t r,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    radix_integer_t q;
    radix_integer_init(q, radix);
    radix_integer_tdiv_qr(q, r, a, b, radix);
    radix_integer_clear(q, radix);
}

void
radix_integer_fdiv_r(radix_integer_t r,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    radix_integer_t q;
    radix_integer_init(q, radix);
    radix_integer_fdiv_qr(q, r, a, b, radix);
    radix_integer_clear(q, radix);
}

void
radix_integer_cdiv_r(radix_integer_t r,
    const radix_integer_t a, const radix_integer_t b, const radix_t radix)
{
    radix_integer_t q;
    radix_integer_init(q, radix);
    radix_integer_cdiv_qr(q, r, a, b, radix);
    radix_integer_clear(q, radix);
}

/* ------------------------------------------------------------------------- */
/*    GR wrapper                                                             */
/* ------------------------------------------------------------------------- */

#define GR_RADIX_CTX(ctx) ((radix_struct *) (GR_CTX_DATA_AS_PTR(ctx)))

static void
_gr_radix_integer_ctx_clear(gr_ctx_t ctx)
{
    radix_clear(GR_RADIX_CTX(ctx));
    flint_free(GR_RADIX_CTX(ctx));
}

static int
_gr_radix_integer_ctx_write(gr_stream_t out, gr_ctx_t ctx)
{
    int status = GR_SUCCESS;
    status |= gr_stream_write(out, "Integers in radix ");
    status |= gr_stream_write_ui(out, DIGIT_RADIX(GR_RADIX_CTX(ctx)));
    status |= gr_stream_write(out, "^");
    status |= gr_stream_write_ui(out, GR_RADIX_CTX(ctx)->exp);
    status |= gr_stream_write(out, " (radix_integer)");
    return status;
}

static void
_gr_radix_integer_init(radix_integer_t res, gr_ctx_t ctx)
{
    radix_integer_init(res, GR_RADIX_CTX(ctx));
}

static void
_gr_radix_integer_clear(radix_integer_t res, gr_ctx_t ctx)
{
    radix_integer_clear(res, GR_RADIX_CTX(ctx));
}

static void _gr_radix_integer_swap(radix_integer_t x, radix_integer_t y, gr_ctx_t ctx)
{
    FLINT_SWAP(radix_integer_struct, *x, *y);
}

static void _gr_radix_integer_set_shallow(radix_integer_t res, const radix_integer_t x, gr_ctx_t ctx)
{
    *res = *x;
}

static int _gr_radix_integer_set(radix_integer_t res, const radix_integer_t x, gr_ctx_t ctx)
{
    radix_integer_set(res, x, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_radix_integer_randtest(radix_integer_t res, flint_rand_t state, gr_ctx_t ctx)
{
    radix_integer_randtest_limbs(res, state, (n_randlimb(state) & 7) ? 1 : 3, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_radix_integer_write(gr_stream_t out, const radix_integer_t x, gr_ctx_t ctx)
{
    slong size = x->size;
    int status = GR_SUCCESS;

    if (size == 0)
    {
        status |= gr_stream_write(out, "0");
    }
    else
    {
        if (DIGIT_RADIX(GR_RADIX_CTX(ctx)) == 10)
        {
            /* todo: chunked output for huge strings */
            status |= gr_stream_write_free(out, radix_get_str_decimal(NULL,
                x->d, FLINT_ABS(size), size < 0, GR_RADIX_CTX(ctx)));
        }
        else
        {
            /* todo: chunked output for huge strings */
            status |= gr_stream_write_free(out, radix_get_str_sum(NULL,
                x->d, FLINT_ABS(size), size < 0, 0, GR_RADIX_CTX(ctx)));
        }
    }

    return status;
}

static int
_gr_radix_integer_zero(radix_integer_t res, gr_ctx_t ctx)
{
    radix_integer_zero(res, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_radix_integer_one(radix_integer_t res, gr_ctx_t ctx)
{
    radix_integer_one(res, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int
_gr_radix_integer_neg_one(radix_integer_t res, gr_ctx_t ctx)
{
    radix_integer_neg_one(res, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_set_ui(radix_integer_t res, ulong x, gr_ctx_t ctx)
{
    radix_integer_set_ui(res, x, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_set_si(radix_integer_t res, slong x, gr_ctx_t ctx)
{
    radix_integer_set_si(res, x, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_set_fmpz(radix_integer_t res, const fmpz_t x, gr_ctx_t ctx)
{
    radix_integer_set_fmpz(res, x, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_get_fmpz(fmpz_t res, const radix_integer_t x, gr_ctx_t ctx)
{
    radix_integer_get_fmpz(res, x, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_get_ui(ulong * res, const radix_integer_t x, gr_ctx_t ctx)
{
    if (x->size < 0)
        return GR_DOMAIN;

    if (x->size <= 1)
    {
        *res = x->size ? x->d[0] : 0;
        return GR_SUCCESS;
    }
    else
    {
        /* todo */
        fmpz_t t;
        fmpz_init(t);
        radix_integer_get_fmpz(t, x, GR_RADIX_CTX(ctx));
        if (fmpz_cmp_ui(t, UWORD_MAX) > 0)
        {
            fmpz_clear(t);
            *res = 0;
            return GR_DOMAIN;
        }
        else
        {
            *res = fmpz_get_ui(t);
            fmpz_clear(t);
            return GR_SUCCESS;
        }
    }
}

static int _gr_radix_integer_get_si(ulong * res, const radix_integer_t x, gr_ctx_t ctx)
{
    if (x->size == 0)
    {
        *res = 0;
        return GR_SUCCESS;
    }
    else
    {
        /* todo */
        fmpz_t t;
        fmpz_init(t);
        radix_integer_get_fmpz(t, x, GR_RADIX_CTX(ctx));
        if (!fmpz_fits_si(t))
        {
            fmpz_clear(t);
            *res = 0;
            return GR_DOMAIN;
        }
        else
        {
            *res = fmpz_get_si(t);
            fmpz_clear(t);
            return GR_SUCCESS;
        }
    }
}

static truth_t _gr_radix_integer_equal(const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    return radix_integer_equal(x, y, GR_RADIX_CTX(ctx)) ? T_TRUE : T_FALSE;
}

static truth_t _gr_radix_integer_is_zero(const radix_integer_t x, gr_ctx_t ctx)
{
    return radix_integer_is_zero(x, GR_RADIX_CTX(ctx)) ? T_TRUE : T_FALSE;
}

static truth_t _gr_radix_integer_is_one(const radix_integer_t x, gr_ctx_t ctx)
{
    return radix_integer_is_one(x, GR_RADIX_CTX(ctx)) ? T_TRUE : T_FALSE;
}

static truth_t _gr_radix_integer_is_neg_one(const radix_integer_t x, gr_ctx_t ctx)
{
    return radix_integer_is_neg_one(x, GR_RADIX_CTX(ctx)) ? T_TRUE : T_FALSE;
}

static int _gr_radix_integer_neg(radix_integer_t res, const radix_integer_t x, gr_ctx_t ctx)
{
    radix_integer_neg(res, x, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_abs(radix_integer_t res, const radix_integer_t x, gr_ctx_t ctx)
{
    radix_integer_abs(res, x, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_add(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    radix_integer_add(res, x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_sub(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    radix_integer_sub(res, x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_mul(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    radix_integer_mul(res, x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_div(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    return radix_integer_div(res, x, y, GR_RADIX_CTX(ctx)) ? GR_SUCCESS : GR_DOMAIN;
}

static int _gr_radix_integer_divexact(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    if (y->size == 0)
        return GR_DOMAIN;
    radix_integer_divexact(res, x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_fdiv_q(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    if (y->size == 0)
        return GR_DOMAIN;
    radix_integer_fdiv_q(res, x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_fdiv_r(radix_integer_t res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    if (y->size == 0)
        return GR_DOMAIN;
    radix_integer_fdiv_r(res, x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_fdiv_qr(radix_integer_t res1, radix_integer_t res2, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    if (y->size == 0)
        return GR_DOMAIN;
    radix_integer_fdiv_qr(res1, res2, x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_cmp(int * res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    *res = radix_integer_cmp(x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_cmpabs(int * res, const radix_integer_t x, const radix_integer_t y, gr_ctx_t ctx)
{
    *res = radix_integer_cmpabs(x, y, GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}

static int _gr_radix_integer_sgn(radix_integer_t res, const radix_integer_t x, gr_ctx_t ctx)
{
    radix_integer_set_si(res, radix_integer_sgn(x, GR_RADIX_CTX(ctx)), GR_RADIX_CTX(ctx));
    return GR_SUCCESS;
}




int _gr_radix_integer_methods_initialized = 0;


gr_static_method_table _gr_radix_integer_methods;

gr_method_tab_input _gr_radix_integer_methods_input[] =
{
    {GR_METHOD_CTX_CLEAR,       (gr_funcptr) _gr_radix_integer_ctx_clear},
    {GR_METHOD_CTX_WRITE,       (gr_funcptr) _gr_radix_integer_ctx_write},
    {GR_METHOD_CTX_IS_RING,     (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_IS_INTEGRAL_DOMAIN,
                                (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_IS_EXACT,    (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_CTX_IS_THREADSAFE,
                                (gr_funcptr) gr_generic_ctx_predicate_true},
    {GR_METHOD_INIT,            (gr_funcptr) _gr_radix_integer_init},
    {GR_METHOD_CLEAR,           (gr_funcptr) _gr_radix_integer_clear},
    {GR_METHOD_SWAP,            (gr_funcptr) _gr_radix_integer_swap},
    {GR_METHOD_SET_SHALLOW,     (gr_funcptr) _gr_radix_integer_set_shallow},
    {GR_METHOD_RANDTEST,        (gr_funcptr) _gr_radix_integer_randtest},
    {GR_METHOD_WRITE,           (gr_funcptr) _gr_radix_integer_write},
    {GR_METHOD_ZERO,            (gr_funcptr) _gr_radix_integer_zero},
    {GR_METHOD_ONE,             (gr_funcptr) _gr_radix_integer_one},
    {GR_METHOD_NEG_ONE,         (gr_funcptr) _gr_radix_integer_neg_one},
    {GR_METHOD_IS_ZERO,         (gr_funcptr) _gr_radix_integer_is_zero},
    {GR_METHOD_IS_ONE,          (gr_funcptr) _gr_radix_integer_is_one},
    {GR_METHOD_IS_NEG_ONE,      (gr_funcptr) _gr_radix_integer_is_neg_one},
    {GR_METHOD_EQUAL,           (gr_funcptr) _gr_radix_integer_equal},
    {GR_METHOD_SET,             (gr_funcptr) _gr_radix_integer_set},
    {GR_METHOD_SET_SI,          (gr_funcptr) _gr_radix_integer_set_si},
    {GR_METHOD_SET_UI,          (gr_funcptr) _gr_radix_integer_set_ui},
    {GR_METHOD_SET_FMPZ,        (gr_funcptr) _gr_radix_integer_set_fmpz},
    {GR_METHOD_GET_SI,          (gr_funcptr) _gr_radix_integer_get_si},
    {GR_METHOD_GET_UI,          (gr_funcptr) _gr_radix_integer_get_ui},
    {GR_METHOD_GET_FMPZ,        (gr_funcptr) _gr_radix_integer_get_fmpz},
    {GR_METHOD_NEG,             (gr_funcptr) _gr_radix_integer_neg},
    {GR_METHOD_ADD,             (gr_funcptr) _gr_radix_integer_add},
    {GR_METHOD_SUB,             (gr_funcptr) _gr_radix_integer_sub},
    {GR_METHOD_MUL,             (gr_funcptr) _gr_radix_integer_mul},
    {GR_METHOD_DIV,             (gr_funcptr) _gr_radix_integer_div},
    {GR_METHOD_DIVEXACT,        (gr_funcptr) _gr_radix_integer_divexact},
    {GR_METHOD_EUCLIDEAN_DIV,   (gr_funcptr) _gr_radix_integer_fdiv_q},
    {GR_METHOD_EUCLIDEAN_REM,   (gr_funcptr) _gr_radix_integer_fdiv_r},
    {GR_METHOD_EUCLIDEAN_DIVREM,   (gr_funcptr) _gr_radix_integer_fdiv_qr},
    {GR_METHOD_CMP,             (gr_funcptr) _gr_radix_integer_cmp},
    {GR_METHOD_CMPABS,          (gr_funcptr) _gr_radix_integer_cmpabs},
    {GR_METHOD_SGN,             (gr_funcptr) _gr_radix_integer_sgn},
    {GR_METHOD_ABS,             (gr_funcptr) _gr_radix_integer_abs},
    {GR_METHOD_FLOOR,           (gr_funcptr) _gr_radix_integer_set},
    {GR_METHOD_CEIL,            (gr_funcptr) _gr_radix_integer_set},
    {GR_METHOD_TRUNC,           (gr_funcptr) _gr_radix_integer_set},
    {GR_METHOD_NINT,            (gr_funcptr) _gr_radix_integer_set},
    {0,                         (gr_funcptr) NULL},
};

void
gr_ctx_init_radix_integer(gr_ctx_t ctx, ulong n, unsigned int exp)
{
    ctx->which_ring = GR_CTX_RADIX_INTEGER;
    ctx->sizeof_elem = sizeof(radix_integer_struct);
    ctx->size_limit = WORD_MAX;

    GR_CTX_DATA_AS_PTR(ctx) = flint_malloc(sizeof(radix_struct));
    /* todo: error handling */
    radix_init(GR_RADIX_CTX(ctx), n, exp);

    ctx->methods = _gr_radix_integer_methods;

    if (!_gr_radix_integer_methods_initialized)
    {
        gr_method_tab_init(_gr_radix_integer_methods, _gr_radix_integer_methods_input);
        _gr_radix_integer_methods_initialized = 1;
    }
}