static void
_arb_poly_rising_ui_series_bsplit(arb_ptr res,
    arb_srcptr f, slong flen, ulong a, ulong b,
        slong trunc, slong prec)
{
    flen = FLINT_MIN(flen, trunc);

    if (b - a == 1)
    {
        arb_add_ui(res, f, a, prec);
        _arb_vec_set(res + 1, f + 1, flen - 1);
    }
    else
    {
        arb_ptr L, R;
        slong len1, len2;

        slong m = a + (b - a) / 2;

        len1 = poly_pow_length(flen, m - a, trunc);
        len2 = poly_pow_length(flen, b - m, trunc);

        L = _arb_vec_init(len1 + len2);
        R = L + len1;

        _arb_poly_rising_ui_series_bsplit(L, f, flen, a, m, trunc, prec);
        _arb_poly_rising_ui_series_bsplit(R, f, flen, m, b, trunc, prec);

        _arb_poly_mullow(res, L, len1, R, len2,
            FLINT_MIN(trunc, len1 + len2 - 1), prec);

        _arb_vec_clear(L, len1 + len2);
    }
}

/* Assumes poly1 and poly2 are not length 0 and len1 >= len2. */
void
_fmpz_poly_mulmid_classical(fmpz * res, const fmpz * poly1,
                            slong len1, const fmpz * poly2, slong len2)
{
    if ((len1 == 1) && (len2 == 1)) /* Special case if the length of both inputs is 1 */
    {
        fmpz_mul(res, poly1, poly2);
    }
    else                        /* Ordinary case */
    {
        slong i;

        /* Set res[i] = poly1[i]*poly2[0]  */
        _fmpz_vec_scalar_mul_fmpz(res, poly1 + len2 - 1, len1 - len2 + 1,
                                  poly2);

        /* out[i+j] += in1[i]*in2[j] */
        for (i = 0; i < len2 - 1; i++)
            _fmpz_vec_scalar_addmul_fmpz(res, poly2 + len2 - i - 1,
                                         FLINT_MIN(i + 1, len1 - len2 + 1),
                                         poly1 + i);
        for (; i < len1 - 1; i++)
            _fmpz_vec_scalar_addmul_fmpz(res + i - len2 + 2, poly2 + 1,
                                         FLINT_MIN(len2 - 1, len1 - i - 1),
                                         poly1 + i);
    }
}

// truncated multiplication for fmpz
void fmpz_mul_trunc(fmpz_t res, fmpz_t a, fmpz_t b, unsigned long trunc)
{
    unsigned long sizea = FLINT_MIN(fmpz_size(a), trunc);
    unsigned long sizeb = FLINT_MIN(fmpz_size(b), trunc);
    while ((!a[sizea]) && (sizea)) sizea--;
    while ((!b[sizeb]) && (sizeb)) sizeb--;

    if ((sizea == 0) || (sizeb == 0)) {
        res[0] = 0;
        return;
    }

    if (trunc >= sizea + sizeb) {
        mp_limb_t mslimb;
        if (sizea >= sizeb) mslimb = F_mpn_mul(res+1, a+1, sizea, b+1, sizeb);
        else mslimb = F_mpn_mul(res+1, b+1, sizeb, a+1, sizea);
        res[0] = sizea + sizeb - (mslimb == 0);
    } else {
        mp_limb_t mslimb;
        fmpz_t temp = flint_stack_alloc(sizea + sizeb + 1);
        if (sizea >= sizeb) mslimb = F_mpn_mul_trunc(temp+1, a+1, sizea, b+1, sizeb, trunc);
        else mslimb = F_mpn_mul_trunc(temp+1, b+1, sizeb, a+1, sizea, trunc);
        temp[0] = trunc;
        if (UNLIKELY(!mslimb))
            __fmpz_normalise(temp); // normalise if most significant limb == 0
        fmpz_set(res, temp);
        flint_stack_release();
    }
    if ((long) (a[0] ^ b[0]) < 0L) res[0] = -res[0];
}

int
bessel(acb_ptr out, const acb_t inp, void * params, long order, long prec)
{
    acb_ptr t;
    acb_t z;
    ulong n;

    t = _acb_vec_init(order);
    acb_init(z);

    acb_set(t, inp);
    if (order > 1)
        acb_one(t + 1);

    n = 10;
    arb_set_si(acb_realref(z), 20);
    arb_set_si(acb_imagref(z), 10);

    /* z sin(t) */
    _acb_poly_sin_series(out, t, FLINT_MIN(2, order), order, prec);
    _acb_vec_scalar_mul(out, out, order, z, prec);

    /* t n */
    _acb_vec_scalar_mul_ui(t, t, FLINT_MIN(2, order), n, prec);

    _acb_poly_sub(out, t, FLINT_MIN(2, order), out, order, prec);

    _acb_poly_cos_series(out, out, order, order, prec);

    _acb_vec_clear(t, order);
    acb_clear(z);
    return 0;
}

/*
   Assumes poly1 and poly2 are not length 0 and 0 < n <= len1 + len2 - 1.
 */
void
_fmpz_poly_mullow_classical(fmpz * res, const fmpz * poly1, long len1, 
                                        const fmpz * poly2, long len2, long n)
{
    if ((len1 == 1 && len2 == 1) || n == 1) /* Special case if the length of output is 1 */
    {
        fmpz_mul(res, poly1, poly2);
    }
    else                        /* Ordinary case */
    {
        long i;

        /* Set res[i] = poly1[i]*poly2[0] */
        _fmpz_vec_scalar_mul_fmpz(res, poly1, FLINT_MIN(len1, n), poly2);

        /* Set res[i+len1-1] = in1[len1-1]*in2[i] */
        if (n > len1)
            _fmpz_vec_scalar_mul_fmpz(res + len1, poly2 + 1, n - len1,
                                      poly1 + len1 - 1);

        /* out[i+j] += in1[i]*in2[j] */
        for (i = 0; i < FLINT_MIN(len1, n) - 1; i++)
            _fmpz_vec_scalar_addmul_fmpz(res + i + 1, poly2 + 1,
                                         FLINT_MIN(len2, n - i) - 1, 
                                         poly1 + i);
    }
}

static __inline__ void
_arb_poly_addmullow_block(arb_ptr z, fmpz * zz,
                          const fmpz * xz, const fmpz * xexps, const slong * xblocks, slong xlen,
                          const fmpz * yz, const fmpz * yexps, const slong * yblocks, slong ylen,
                          slong n, slong prec, int squaring)
{
    slong i, j, k, xp, yp, xl, yl, bn;
    fmpz_t zexp;

    fmpz_init(zexp);

    if (squaring)
    {
        for (i = 0; (xp = xblocks[i]) != xlen; i++)
        {
            if (2 * xp >= n)
                continue;

            xl = xblocks[i + 1] - xp;
            bn = FLINT_MIN(2 * xl - 1, n - 2 * xp);
            xl = FLINT_MIN(xl, bn);

            _fmpz_poly_sqrlow(zz, xz + xp, xl, bn);
            _fmpz_add2_fast(zexp, xexps + i, xexps + i, 0);

            for (k = 0; k < bn; k++)
                arb_add_fmpz_2exp(z + 2 * xp + k, z + 2 * xp + k, zz + k, zexp, prec);
        }
    }

    for (i = 0; (xp = xblocks[i]) != xlen; i++)
    {
        for (j = squaring ? i + 1 : 0; (yp = yblocks[j]) != ylen; j++)
        {
            if (xp + yp >= n)
                continue;

            xl = xblocks[i + 1] - xp;
            yl = yblocks[j + 1] - yp;
            bn = FLINT_MIN(xl + yl - 1, n - xp - yp);
            xl = FLINT_MIN(xl, bn);
            yl = FLINT_MIN(yl, bn);

            if (xl >= yl)
                _fmpz_poly_mullow(zz, xz + xp, xl, yz + yp, yl, bn);
            else
                _fmpz_poly_mullow(zz, yz + yp, yl, xz + xp, xl, bn);

            _fmpz_add2_fast(zexp, xexps + i, yexps + j, squaring);

            for (k = 0; k < bn; k++)
                arb_add_fmpz_2exp(z + xp + yp + k, z + xp + yp + k, zz + k, zexp, prec);
        }
    }

    fmpz_clear(zexp);
}

void
acb_hypgeom_erf(acb_t res, const acb_t z, slong prec)
{
    double x, y, absz2, logz;
    slong prec2;

    if (!acb_is_finite(z))
    {
        acb_indeterminate(res);
        return;
    }

    if (acb_is_zero(z))
    {
        acb_zero(res);
        return;
    }

    if ((arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 0) < 0 &&
            arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 0) < 0))
    {
        acb_hypgeom_erf_1f1a(res, z, prec);
        return;
    }

    if ((arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 64) > 0 ||
            arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 64) > 0))
    {
        acb_hypgeom_erf_asymp(res, z, prec, prec);
        return;
    }

    x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
    y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);

    absz2 = x * x + y * y;
    logz = 0.5 * log(absz2);

    if (logz - absz2 < -(prec + 8) * 0.69314718055994530942)
    {
        /* If the asymptotic term is small, we can
           compute with reduced precision */
        prec2 = FLINT_MIN(prec + 4 + (y*y - x*x - logz) * 1.4426950408889634074, (double) prec);
        prec2 = FLINT_MAX(8, prec2);
        prec2 = FLINT_MIN(prec2, prec);

        acb_hypgeom_erf_asymp(res, z, prec, prec2);
    }
    else if (arf_cmpabs(arb_midref(acb_imagref(z)), arb_midref(acb_realref(z))) > 0)
    {
        acb_hypgeom_erf_1f1a(res, z, prec);
    }
    else
    {
        acb_hypgeom_erf_1f1b(res, z, prec);
    }
}

void
nmod_poly_compose_series_horner(nmod_poly_t res, 
                    const nmod_poly_t poly1, const nmod_poly_t poly2, long n)
{
    long len1 = poly1->length;
    long len2 = poly2->length;
    long lenr;

    if (len2 != 0 && poly2->coeffs[0] != 0)
    {
        printf("exception: nmod_poly_compose_series_horner: inner polynomial "
                "must have zero constant term\n");
        abort();
    }

    if (len1 == 0 || n == 0)
    {
        nmod_poly_zero(res);
        return;
    }

    if (len2 == 0 || len1 == 1)
    {
        nmod_poly_fit_length(res, 1);
        res->coeffs[0] = poly1->coeffs[0];
        res->length = 1;
        _nmod_poly_normalise(res);
        return;
    }

    lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n);
    len1 = FLINT_MIN(len1, lenr);
    len2 = FLINT_MIN(len2, lenr);

    if ((res != poly1) && (res != poly2))
    {
        nmod_poly_fit_length(res, lenr);
        _nmod_poly_compose_series_horner(res->coeffs, poly1->coeffs, len1, 
                                        poly2->coeffs, len2, lenr, res->mod);
        res->length = lenr;
        _nmod_poly_normalise(res);
    }
    else
    {
        nmod_poly_t t;
        nmod_poly_init2_preinv(t, res->mod.n, res->mod.ninv, lenr);
        _nmod_poly_compose_series_horner(t->coeffs, poly1->coeffs, len1,
                                        poly2->coeffs, len2, lenr, res->mod);
        t->length = lenr;
        _nmod_poly_normalise(t);
        nmod_poly_swap(res, t);
        nmod_poly_clear(t);
    }
}

void
acb_poly_compose_series(acb_poly_t res,
                    const acb_poly_t poly1,
                    const acb_poly_t poly2, slong n, slong prec)
{
    slong len1 = poly1->length;
    slong len2 = poly2->length;
    slong lenr;

    if (len2 != 0 && !acb_is_zero(poly2->coeffs))
    {
        flint_printf("exception: compose_series: inner "
                "polynomial must have zero constant term\n");
        abort();
    }

    if (len1 == 0 || n == 0)
    {
        acb_poly_zero(res);
        return;
    }

    if (len2 == 0 || len1 == 1)
    {
        acb_poly_set_acb(res, poly1->coeffs);
        return;
    }

    lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n);
    len1 = FLINT_MIN(len1, lenr);
    len2 = FLINT_MIN(len2, lenr);

    if ((res != poly1) && (res != poly2))
    {
        acb_poly_fit_length(res, lenr);
        _acb_poly_compose_series(res->coeffs, poly1->coeffs, len1,
                                        poly2->coeffs, len2, lenr, prec);
        _acb_poly_set_length(res, lenr);
        _acb_poly_normalise(res);
    }
    else
    {
        acb_poly_t t;
        acb_poly_init2(t, lenr);
        _acb_poly_compose_series(t->coeffs, poly1->coeffs, len1,
                                        poly2->coeffs, len2, lenr, prec);
        _acb_poly_set_length(t, lenr);
        _acb_poly_normalise(t);
        acb_poly_swap(res, t);
        acb_poly_clear(t);
    }
}

void
fmpz_poly_compose_series(fmpz_poly_t res, 
                    const fmpz_poly_t poly1, const fmpz_poly_t poly2, long n)
{
    long len1 = poly1->length;
    long len2 = poly2->length;
    long lenr;

    if (len2 != 0 && !fmpz_is_zero(poly2->coeffs))
    {
        printf("exception: fmpz_poly_compose_series: inner polynomial "
                "must have zero constant term\n");
        abort();
    }

    if (len1 == 0 || n == 0)
    {
        fmpz_poly_zero(res);
        return;
    }

    if (len2 == 0 || len1 == 1)
    {
        fmpz_poly_set_fmpz(res, poly1->coeffs);
        return;
    }

    lenr = FLINT_MIN((len1 - 1) * (len2 - 1) + 1, n);
    len1 = FLINT_MIN(len1, lenr);
    len2 = FLINT_MIN(len2, lenr);

    if ((res != poly1) && (res != poly2))
    {
        fmpz_poly_fit_length(res, lenr);
        _fmpz_poly_compose_series(res->coeffs, poly1->coeffs, len1, 
                                               poly2->coeffs, len2, lenr);
        _fmpz_poly_set_length(res, lenr);
        _fmpz_poly_normalise(res);
    }
    else
    {
        fmpz_poly_t t;
        fmpz_poly_init2(t, lenr);
        _fmpz_poly_compose_series(t->coeffs, poly1->coeffs, len1,
                                             poly2->coeffs, len2, lenr);
        _fmpz_poly_set_length(t, lenr);
        _fmpz_poly_normalise(t);
        fmpz_poly_swap(res, t);
        fmpz_poly_clear(t);
    }
}

void fmpz_poly_sqrlow(fmpz_poly_t res, const fmpz_poly_t poly, long n)
{
    const long len = poly->length;

    if (len == 0 || n == 0)
    {
        fmpz_poly_zero(res);
        return;
    }

    if (res == poly)
    {
        fmpz_poly_t t;
        fmpz_poly_init2(t, n);
        fmpz_poly_sqrlow(t, poly, n);
        fmpz_poly_swap(res, t);
        fmpz_poly_clear(t);
        return;
    }

    n = FLINT_MIN(2 * len - 1, n);

    fmpz_poly_fit_length(res, n);
    _fmpz_poly_sqrlow(res->coeffs, poly->coeffs, len, n);
    _fmpz_poly_set_length(res, n);
    _fmpz_poly_normalise(res);
}

void
_elem_poly_sub(elem_ptr res, elem_srcptr poly1, long len1,
    elem_srcptr poly2, long len2, const ring_t ring)
{
    long i, min;
    long size = ring->size;

    if (ring->type == TYPE_FMPZ && 0)
    {
        _fmpz_poly_sub(res, poly1, len1, poly2, len2);
        return;
    }

    min = FLINT_MIN(len1, len2);

    for (i = 0; i < min; i++)
        elem_sub(INDEX(res, i, size), SRC_INDEX(poly1, i, size), SRC_INDEX(poly2, i, size), ring);

    if (poly1 != res)
        for (i = min; i < len1; i++)
            elem_set(INDEX(res, i, size), SRC_INDEX(poly1, i, size), ring);

    for (i = min; i < len2; i++)
        elem_neg(INDEX(res, i, size), SRC_INDEX(poly2, i, size), ring);
}

void fmpq_poly_get_slice(fmpq_poly_t rop, const fmpq_poly_t op, long i, long j)
{
    i = FLINT_MAX(i, 0);
    j = FLINT_MIN(j, op->length);

    if (i < j)
    {
        long k;

        if (rop == op)
        {
            for (k = 0; k < i; k++)
                fmpz_zero(rop->coeffs + k);
            for (k = j; k < rop->length; k++)
                fmpz_zero(rop->coeffs + k);
            fmpq_poly_canonicalise(rop);
        }
        else
        {
            fmpq_poly_fit_length(rop, j);
            _fmpq_poly_set_length(rop, j);

            _fmpz_vec_set(rop->coeffs + i, op->coeffs + i, j - i);
            fmpz_set(rop->den, op->den);
            fmpq_poly_canonicalise(rop);
        }
    }
    else
    {
        fmpq_poly_zero(rop);
    }
}

void
_arb_poly_sinh_cosh_series(arb_ptr s, arb_ptr c, const arb_srcptr h, slong hlen, slong n, slong prec)
{
    hlen = FLINT_MIN(hlen, n);

    if (hlen == 1)
    {
        arb_sinh_cosh(s, c, h, prec);
        _arb_vec_zero(s + 1, n - 1);
        _arb_vec_zero(c + 1, n - 1);
    }
    else if (n == 2)
    {
        arb_t t;
        arb_init(t);
        arb_set(t, h + 1);
        arb_sinh_cosh(s, c, h, prec);
        arb_mul(s + 1, c, t, prec);
        arb_mul(c + 1, s, t, prec);
        arb_clear(t);
    }
    else if (hlen < 60 || n < 120)
        _arb_poly_sinh_cosh_series_basecase(s, c, h, hlen, n, prec);
    else
        _arb_poly_sinh_cosh_series_exponential(s, c, h, hlen, n, prec);
}

void
_fmprb_poly_sin_series(fmprb_ptr g, fmprb_srcptr h, long hlen, long n, long prec)
{
    hlen = FLINT_MIN(hlen, n);

    if (hlen == 1)
    {
        fmprb_sin(g, h, prec);
        _fmprb_vec_zero(g + 1, n - 1);
    }
    else if (n == 2)
    {
        fmprb_t t;
        fmprb_init(t);
        fmprb_sin_cos(g, t, h, prec);
        fmprb_mul(g + 1, h + 1, t, prec);  /* safe since hlen >= 2 */
        fmprb_clear(t);
    }
    else
    {
        fmprb_ptr t = _fmprb_vec_init(n);
        _fmprb_poly_sin_cos_series(g, t, h, hlen, n, prec);
        _fmprb_vec_clear(t, n);
    }
}

void
_acb_poly_zeta_cpx_series(acb_ptr z, const acb_t s, const acb_t a, int deflate, long d, long prec)
{
    ulong M, N;
    long i;
    arf_t bound;
    arb_ptr vb;

    if (d < 1)
        return;

    if (!acb_is_finite(s) || !acb_is_finite(a))
    {
        _acb_vec_indeterminate(z, d);
        return;
    }

    arf_init(bound);
    vb = _arb_vec_init(d);

    _acb_poly_zeta_em_choose_param(bound, &N, &M, s, a, FLINT_MIN(d, 2), prec, MAG_BITS);
    _acb_poly_zeta_em_bound(vb, s, a, N, M, d, MAG_BITS);

    _acb_poly_zeta_em_sum(z, s, a, deflate, N, M, d, prec);

    for (i = 0; i < d; i++)
    {
        arb_get_abs_ubound_arf(bound, vb + i, MAG_BITS);
        arb_add_error_arf(acb_realref(z + i), bound);
        arb_add_error_arf(acb_imagref(z + i), bound);
    }

    arf_clear(bound);
    _arb_vec_clear(vb, d);
}

void
_acb_poly_binomial_transform_basecase(acb_ptr b, acb_srcptr a, slong alen, slong len, slong prec)
{
    slong n, k;

    fmpz_t t;
    fmpz_init(t);

    for (n = 0; n < len; n++)
    {
        acb_zero(b + n);

        for (k = 0; k < FLINT_MIN(n + 1, alen); k++)
        {
            if (k == 0)
            {
                fmpz_one(t);
            }
            else
            {
                fmpz_mul_si(t, t, -(n - k + 1));
                fmpz_divexact_ui(t, t, k);
            }

            acb_addmul_fmpz(b + n, a + k, t, prec);
        }
    }

    fmpz_clear(t);
}

void
fmpr_divappr_abs_ubound(fmpr_t z, const fmpr_t x, const fmpr_t y, slong prec)
{
    if (fmpr_is_special(x) || fmpr_is_special(y) || fmpz_is_pm1(fmpr_manref(y)))
    {
        fmpr_div(z, x, y, prec, FMPR_RND_UP);
        fmpr_abs(z, z);
    }
    else
    {
        fmpz_t t, u;
        slong xbits, ybits, tbits, ubits, shift;

        xbits = fmpz_bits(fmpr_manref(x));
        ybits = fmpz_bits(fmpr_manref(y));

        fmpz_init(t);
        fmpz_init(u);

        ubits = FLINT_MIN(ybits, prec);
        tbits = prec + ubits + 1;

        /* upper bound for |x|, shifted */
        if (xbits <= tbits)
        {
            fmpz_mul_2exp(t, fmpr_manref(x), tbits - xbits);
            fmpz_abs(t, t);
        }
        else if (fmpz_sgn(fmpr_manref(x)) > 0)
        {
            fmpz_cdiv_q_2exp(t, fmpr_manref(x), xbits - tbits);
        }
        else
        {
            fmpz_fdiv_q_2exp(t, fmpr_manref(x), xbits - tbits);
            fmpz_neg(t, t);
        }

        /* lower bound for |y|, shifted */
        if (ybits <= ubits)
            fmpz_mul_2exp(u, fmpr_manref(y), ubits - ybits);
        else
            fmpz_tdiv_q_2exp(u, fmpr_manref(y), ybits - ubits);
        fmpz_abs(u, u);

        fmpz_cdiv_q(fmpr_manref(z), t, u);

        shift = (ubits - ybits) - (tbits - xbits);
        fmpz_sub(fmpr_expref(z), fmpr_expref(x), fmpr_expref(y));
        if (shift >= 0)
            fmpz_add_ui(fmpr_expref(z), fmpr_expref(z), shift);
        else
            fmpz_sub_ui(fmpr_expref(z), fmpr_expref(z), -shift);

        _fmpr_normalise(fmpr_manref(z), fmpr_expref(z), prec, FMPR_RND_UP);

        fmpz_clear(t);
        fmpz_clear(u);
    }
}

void
_nmod_poly_div_basecase_2(mp_ptr Q, mp_ptr W,
                             mp_srcptr A, long A_len, mp_srcptr B, long B_len,
                             nmod_t mod)
{
    long coeff, i, len;
    mp_limb_t lead_inv = n_invmod(B[B_len - 1], mod.n);
    mp_ptr B2, R2;
    mp_srcptr Btop;
    
    B2 = W;
    for (i = 0; i < B_len - 1; i++)
    {
        B2[2 * i] = B[i];
        B2[2 * i + 1] = 0;
    }
    Btop = B2 + 2*(B_len - 1);

    R2 = W + 2*(B_len - 1);
    for (i = 0; i < A_len - B_len + 1; i++)
    {
        R2[2 * i] = A[B_len + i - 1];
        R2[2 * i + 1] = 0;
    }

    coeff = A_len - B_len;
    
    while (coeff >= 0)
    {
        mp_limb_t r_coeff;
        r_coeff =
            n_ll_mod_preinv(R2[2 * coeff + 1], R2[2 * coeff], mod.n, mod.ninv);

        while (coeff >= 0 && r_coeff == 0L)
        {
            Q[coeff--] = 0L;
            if (coeff >= 0)
                r_coeff =
                    n_ll_mod_preinv(R2[2 * coeff + 1], R2[2 * coeff], mod.n,
                                    mod.ninv);
        }

        if (coeff >= 0)
        {
            mp_limb_t c, * R_sub;

            Q[coeff] =
                n_mulmod2_preinv(r_coeff, lead_inv, mod.n, mod.ninv);

            c = n_negmod(Q[coeff], mod.n);

            len = FLINT_MIN(B_len - 1, coeff);
            R_sub = R2 + 2 * (coeff - len);
            if (len > 0)
                mpn_addmul_1(R_sub, Btop - 2*len, 2 * len, c);

            coeff--;
        }
    }
}

void padic_poly_sub(padic_poly_t f, 
                    const padic_poly_t g, const padic_poly_t h, 
                    const padic_ctx_t ctx)
{
    const slong lenG = g->length;
    const slong lenH = h->length;
    const slong lenF = FLINT_MAX(lenG, lenH);

    if (lenG == 0)
    {
        padic_poly_neg(f, h, ctx);
        return;
    }
    if (lenH == 0)
    {
        padic_poly_set(f, g, ctx);
        return;
    }
    if ((lenG == 0 && lenH == 0) || (FLINT_MIN(g->val, h->val) >= f->N))
    {
        padic_poly_zero(f);
        return;
    }

    padic_poly_fit_length(f, lenF);

    _padic_poly_sub(f->coeffs, &(f->val), f->N, 
                    g->coeffs, g->val, lenG, g->N, 
                    h->coeffs, h->val, lenH, h->N, ctx);

    _padic_poly_set_length(f, lenF);
    _padic_poly_normalise(f);
}