void _fmpz_poly_hensel_lift_without_inverse(fmpz *G, fmpz *H, 
    const fmpz *f, long lenF, 
    const fmpz *g, long lenG, const fmpz *h, long lenH, 
    const fmpz *a, long lenA, const fmpz *b, long lenB, 
    const fmpz_t p, const fmpz_t p1)
{
    const fmpz one[1] = {1l};
    const long lenM = FLINT_MAX(lenG, lenH);
    const long lenE = FLINT_MAX(lenG + lenB - 2, lenH + lenA - 2);
    const long lenD = FLINT_MAX(lenE, lenF);
    fmpz *C, *D, *E, *M;

    C = _fmpz_vec_init(lenF + lenD + lenE + lenM);
    D = C + lenF;
    E = D + lenD;
    M = E + lenE;

    if (lenG >= lenH)
        _fmpz_poly_mul(C, g,lenG, h, lenH);
    else
        _fmpz_poly_mul(C, h, lenH, g, lenG);
    _fmpz_vec_sub(C, f, C, lenF);
    _fmpz_vec_scalar_divexact_fmpz(D, C, lenF, p);
    _fmpz_vec_scalar_mod_fmpz(C, D, lenF, p1);

    lift(G, g, lenG, b, lenB);

    lift(H, h, lenH, a, lenA);

    _fmpz_vec_clear(C, lenF + lenD + lenE + lenM);
}

int main(void)
{
    /* Example 3 */
    long n = 4;
    long d = 6;
    long N = 500;
    fmpz_t p = {97L};
    fmpz a[5] = {1, 2, 3, 4, 5};


    padic_ctx_t pctx;

    padic_mat_t F;

    long i, lenB = gmc_basis_size(n, d);

    padic_ctx_init(pctx, p, FLINT_MAX(0, N), N, PADIC_VAL_UNIT);

    padic_mat_init2(F, lenB, lenB, N);
    diagfrob(F, a, n, d, N, pctx, 1);

    padic_mat_print_pretty(F, pctx);
    printf("\n\n");

    /* Clean-up */
    fmpz_clear(p);
    padic_mat_clear(F);
    padic_ctx_clear(pctx);

    return EXIT_SUCCESS;
}

void
arb_fprintd(FILE * file, const arb_t x, slong digits)
{
    arf_fprintd(file, arb_midref(x), FLINT_MAX(digits, 1));
    flint_fprintf(file, " +/- ");
    mag_fprintd(file, arb_radref(x), 5);
}

void
_hypgeom_precompute(hypgeom_t hyp, const fmpz_poly_t P, const fmpz_poly_t Q)
{
    slong k;

    fmpz_t t;
    fmpz_init(t);

    hyp->r = fmpz_poly_degree(Q) - fmpz_poly_degree(P);
    hyp->boundC = hypgeom_root_norm(P);
    hyp->boundD = hypgeom_root_norm(Q);
    hyp->boundK = 1 + FLINT_MAX(hyp->boundC, 2 * hyp->boundD);

    mag_one(hyp->MK);

    for (k = 1; k <= hyp->boundK; k++)
    {
        fmpz_poly_evaluate_si(t, P, k);
        mag_mul_fmpz(hyp->MK, hyp->MK, t);

        fmpz_poly_evaluate_si(t, Q, k);
        mag_div_fmpz(hyp->MK, hyp->MK, t);
    }

    fmpz_clear(t);
}

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);
}

void
acb_log1p(acb_t r, const acb_t z, slong prec)
{
    slong magz, magx, magy;

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

    magx = arf_abs_bound_lt_2exp_si(arb_midref(acb_realref(z)));
    magy = arf_abs_bound_lt_2exp_si(arb_midref(acb_imagref(z)));
    magz = FLINT_MAX(magx, magy);

    if (magz < -prec)
    {
        acb_log1p_tiny(r, z, prec);
    }
    else
    {
        if (magz < 0)
            acb_add_ui(r, z, 1, prec + (-magz) + 4);
        else
            acb_add_ui(r, z, 1, prec + 4);

        acb_log(r, r, prec);
    }
}

void
bernoulli_cache_compute(long n)
{
    if (bernoulli_cache_num < n)
    {
        long i, new_num;
        bernoulli_rev_t iter;

        if (bernoulli_cache_num == 0)
        {
            flint_register_cleanup_function(bernoulli_cleanup);
        }

        new_num = FLINT_MAX(bernoulli_cache_num + 128, n);

        bernoulli_cache = flint_realloc(bernoulli_cache, new_num * sizeof(fmpq));
        for (i = bernoulli_cache_num; i < new_num; i++)
            fmpq_init(bernoulli_cache + i);

        i = new_num - 1;
        i -= (i % 2);
        bernoulli_rev_init(iter, i);
        for ( ; i >= bernoulli_cache_num; i -= 2)
        {
            bernoulli_rev_next(fmpq_numref(bernoulli_cache + i),
                fmpq_denref(bernoulli_cache + i), iter);
        }
        bernoulli_rev_clear(iter);

        if (new_num > 1)
            fmpq_set_si(bernoulli_cache + 1, -1, 2);

        bernoulli_cache_num = new_num;
    }
}

int renf_elem_relative_condition_number_2exp(slong * cond, renf_elem_t a, renf_t nf)
{
    fmpz * p;
    slong len;

    if (nf_elem_is_rational(a->elem, nf->nf))
    {
        *cond = 0;
        return 1;
    }

    if (nf->nf->flag & NF_QUADRATIC)
    {
        p = QNF_ELEM_NUMREF(a->elem);
        len = 2;
    }
    else
    {
        p = NF_ELEM(a->elem)->coeffs;
        len = NF_ELEM(a->elem)->length;
    }

    return _fmpz_poly_relative_condition_number_2exp(cond, p,
                len, nf->emb, FLINT_MAX(nf->prec, 16));
}

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);
    }
}

slong
fmpz_mat_max_bits(const fmpz_mat_t mat)
{
    slong i;
    slong bits, row_bits, sign;

    sign = 1;
    bits = 0;

    if (mat->r == 0 || mat->c == 0)
        return 0;

    for (i = 0; i < mat->r; i++)
    {
        row_bits = _fmpz_vec_max_bits(mat->rows[i], mat->c);
        if (row_bits < 0)
        {
            row_bits = -row_bits;
            sign = -1;
        }
        bits = FLINT_MAX(bits, row_bits);
    }

    return bits * sign;
}

int
_arf_add_eps(arf_t s, const arf_t x, int sgn, long prec, arf_rnd_t rnd)
{
    arf_t t;
    long bits;

    bits = arf_bits(x);

    if (bits == 0)
    {
        printf("_arf_add_eps\n");
        abort();
    }

    bits = FLINT_MAX(bits, prec) + 10;

    arf_init(t);
    arf_set_si(t, sgn);
    arf_mul_2exp_fmpz(t, t, ARF_EXPREF(x));
    arf_mul_2exp_si(t, t, -bits);
    arf_add(s, x, t, prec, rnd);
    arf_clear(t);

    return 1;
}

void
arb_log_arf_huge(arb_t z, const arf_t x, slong prec)
{
    arf_t t;
    arb_t c;
    fmpz_t exp;
    slong wp;

    arf_init(t);
    arb_init(c);
    fmpz_init(exp);

    fmpz_neg(exp, ARF_EXPREF(x));
    arf_mul_2exp_fmpz(t, x, exp);

    wp = prec + 4 - fmpz_bits(exp);
    wp = FLINT_MAX(wp, 4);

    arb_log_arf(z, t, wp);
    arb_const_log2(c, prec + 4);
    arb_submul_fmpz(z, c, exp, prec);

    arf_clear(t);
    arb_clear(c);
    fmpz_clear(exp);
}

static void 
_padic_log_bsplit(fmpz_t z, const fmpz_t y, long v, const fmpz_t p, long N)
{
    fmpz_t P, B, T;
    long n;

    if (fmpz_fits_si(p))
        n = _padic_log_bound(v, N, fmpz_get_si(p));
    else
        n = (N - 1) / v;

    n = FLINT_MAX(n, 2);

    fmpz_init(P);
    fmpz_init(B);
    fmpz_init(T);

    _padic_log_bsplit_series(P, B, T, y, 1, n);

    n = fmpz_remove(B, B, p);
    fmpz_pow_ui(P, p, n);
    fmpz_divexact(T, T, P);

    _padic_inv(B, B, p, N);
    fmpz_mul(z, T, B);

    fmpz_clear(P);
    fmpz_clear(B);
    fmpz_clear(T);
}

void padic_ctx_init(padic_ctx_t ctx, const fmpz_t p, long N,
                    enum padic_print_mode mode)
{
    fmpz_init(ctx->p);
    fmpz_set(ctx->p, p);

    ctx->N = N;

    ctx->pinv = (!COEFF_IS_MPZ(*p)) ? n_precompute_inverse(fmpz_get_ui(p)) : 0;

    if (N > 0)
    {
        long i, len;

        ctx->min = FLINT_MAX(1, N - 10);
        ctx->max = N + 10;
        len      = ctx->max - ctx->min;

        ctx->pow = _fmpz_vec_init(len);

        fmpz_pow_ui(ctx->pow, p, ctx->min);
        for (i = 1; i < len; i++)
            fmpz_mul(ctx->pow + i, ctx->pow + (i - 1), p);
    }
    else
    {
        ctx->min = 0;
        ctx->max = 0;
        ctx->pow = NULL;
    }

    ctx->mode = mode;
}

void F_mpz_mod_poly_sub(F_mpz_mod_poly_t res, const F_mpz_mod_poly_t poly1, const F_mpz_mod_poly_t poly2)
{
   ulong longer = FLINT_MAX(poly1->length, poly2->length);

	F_mpz_mod_poly_fit_length(res, longer);
	
	_F_mpz_mod_poly_sub(res, poly1, poly2);
}

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);
    }
}

slong
_acb_get_mid_mag(const acb_t z)
{
    slong rm, im;

    rm = arf_abs_bound_lt_2exp_si(arb_midref(acb_realref(z)));
    im = arf_abs_bound_lt_2exp_si(arb_midref(acb_imagref(z)));

    return FLINT_MAX(rm, im);
}

void
elem_poly_sub(elem_poly_struct * res,
    const elem_poly_struct * op1, const elem_poly_struct * op2, const ring_t ring)
{
    long max = FLINT_MAX(op1->length, op2->length);

    elem_poly_fit_length(res, max, ring);
    _elem_poly_sub(res->coeffs, op1->coeffs, op1->length, op2->coeffs, op2->length, ring->parent);
    elem_poly_set_length(res, max, ring);
    elem_poly_normalise(res, ring);
}

void
fmpr_printd(const fmpr_t x, long digits)
{
    mpfr_t t;
    mpfr_init2(t, digits * 3.33 + 10);
    mpfr_set_emin(MPFR_EMIN_MIN);
    mpfr_set_emax(MPFR_EMAX_MAX);
    fmpr_get_mpfr(t, x, MPFR_RNDN);
    mpfr_printf("%.*Rg", FLINT_MAX(digits, 1), t);
    mpfr_clear(t);
}

void renf_randtest(renf_t nf, flint_rand_t state, slong len, slong prec, mp_bitcnt_t bits)
{
    fmpz_poly_t p;
    fmpq_poly_t p2;
    fmpz * c_array;
    slong * k_array;
    slong n_interval, n_exact;
    ulong i;
    arb_t emb;

    /* compute a random irreducible polynomial */
    if (len <= 1)
    {
        fprintf(stderr, "ERROR (renf_randtest): got length < 2\n");
        abort();
    }

    fmpz_poly_init(p);
    do{
        fmpz_poly_randtest_irreducible(p, state, len, bits);
    }while(!fmpz_poly_has_real_root(p));

    /* pick a random real root */
    c_array = _fmpz_vec_init(p->length);
    k_array = malloc((p->length) * sizeof(slong));
    n_interval = 0;

    fmpz_poly_isolate_real_roots(NULL, &n_exact,
            c_array, k_array, &n_interval, p);

    if (n_interval == 0)
    {
        fprintf(stderr, "Runtime error\n");
        abort();
    }

    i = n_randint(state, n_interval);

    /* construct the associated number field */
    arb_init(emb);
    arb_from_interval(emb, c_array+i, k_array[i], fmpz_bits(c_array + i) + FLINT_MAX(k_array[i], 0) + 2);

    fmpq_poly_init(p2);
    fmpq_poly_set_fmpz_poly(p2, p);

    /* NOTE: renf init might not be happy with the ball emb */
    renf_init(nf, p2, emb, prec);

    _fmpz_vec_clear(c_array, p->length);
    free(k_array);
    fmpz_poly_clear(p);
    fmpq_poly_clear(p2);
    arb_clear(emb);
}