自学教程：C++ FLINT_MIN函数代码示例

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 fmpzvoid 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];}

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

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

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

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

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

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

void_arb_poly_log1p_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec){    arb_t a;    flen = FLINT_MIN(flen, n);    arb_init(a);    arb_log1p(a, f, prec);    if (flen == 1)    {        _arb_vec_zero(res + 1, n - 1);    }    else if (n == 2)    {        arb_add_ui(res, f + 0, 1, prec);        arb_div(res + 1, f + 1, res + 0, prec);    }    else if (_arb_vec_is_zero(f + 1, flen - 2))  /* f = a + bx^d */    {        slong i, j, d = flen - 1;        arb_add_ui(res, f + 0, 1, prec);        for (i = 1, j = d; j < n; j += d, i++)        {            if (i == 1)                arb_div(res + j, f + d, res, prec);            else                arb_mul(res + j, res + j - d, res + d, prec);            _arb_vec_zero(res + j - d + 1, flen - 2);        }        _arb_vec_zero(res + j - d + 1, n - (j - d + 1));        for (i = 2, j = 2 * d; j < n; j += d, i++)            arb_div_si(res + j, res + j, i % 2 ? i : -i, prec);    }    else    {        arb_ptr f_diff, f_inv;        slong alloc;        alloc = n + flen;        f_inv = _arb_vec_init(alloc);        f_diff = f_inv + n;        arb_add_ui(f_diff, f, 1, prec);        _arb_vec_set(f_diff + 1, f + 1, flen - 1);        _arb_poly_inv_series(f_inv, f_diff, flen, n, prec);        _arb_poly_derivative(f_diff, f, flen, prec);        _arb_poly_mullow(res, f_inv, n - 1, f_diff, flen - 1, n - 1, prec);        _arb_poly_integral(res, res, n, prec);        _arb_vec_clear(f_inv, alloc);    }    arb_swap(res, a);    arb_clear(a);}

void _nmod_poly_reverse(mp_ptr output, mp_srcptr input, long len, long m){    long i, min;    mp_limb_t temp;          if (input != output)    {        min = FLINT_MIN(m, len);                for (i = 0; i < min; i++)            output[m - i - 1] = input[i];        for ( ; i < m; i++)            output[m - i - 1] = 0L;    } else    {        for (i = 0; i < m/2; i++)        {            temp = i < len ? input[i] : 0;                     output[i] = m - i - 1 < len ? input[m - i - 1] : 0;            output[m - i - 1] = temp;        }        if (m & 1 && i >= len) output[i] = 0;    }}

void_arb_poly_sqrt_series(arb_ptr g,    arb_srcptr h, long hlen, long len, long prec){    hlen = FLINT_MIN(hlen, len);    if (hlen == 1)    {        arb_sqrt(g, h, prec);        _arb_vec_zero(g + 1, len - 1);    }    else if (len == 2)    {        arb_sqrt(g, h, prec);        arb_div(g + 1, h + 1, h, prec);        arb_mul(g + 1, g + 1, g, prec);        arb_mul_2exp_si(g + 1, g + 1, -1);    }    else    {        arb_ptr t;        t = _arb_vec_init(len);        _arb_poly_rsqrt_series(t, h, hlen, len, prec);        _arb_poly_mullow(g, t, len, h, hlen, len, prec);        _arb_vec_clear(t, len);    }}

void_arb_poly_sinh_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec){    hlen = FLINT_MIN(hlen, n);    if (hlen == 1)    {        arb_sinh(g, h, prec);        _arb_vec_zero(g + 1, n - 1);    }    else if (n == 2)    {        arb_t t;        arb_init(t);        arb_sinh_cosh(g, t, h, prec);        arb_mul(g + 1, h + 1, t, prec);  /* safe since hlen >= 2 */        arb_clear(t);    }    else    {        arb_ptr t = _arb_vec_init(n);        _arb_poly_sinh_cosh_series(g, t, h, hlen, n, prec);        _arb_vec_clear(t, n);    }}

void_acb_dirichlet_hardy_z_series(acb_ptr res, acb_srcptr s, slong slen,    const dirichlet_group_t G, const dirichlet_char_t chi,    slong len, slong prec){    slen = FLINT_MIN(slen, len);    if (len == 0)        return;    if (slen == 1)    {        acb_dirichlet_hardy_z(res, s, G, chi, 1, prec);        _acb_vec_zero(res + 1, len - 1);    }    else    {        acb_ptr t, u;        t = _acb_vec_init(len);        u = _acb_vec_init(slen);        acb_dirichlet_hardy_z(t, s, G, chi, len, prec);        /* compose with nonconstant part */        acb_zero(u);        _acb_vec_set(u + 1, s + 1, slen - 1);        _acb_poly_compose_series(res, t, len, u, slen, len, prec);        _acb_vec_clear(t, len);        _acb_vec_clear(u, slen);    }}

void_acb_poly_revert_series_lagrange_fast(acb_ptr Qinv, acb_srcptr Q, slong Qlen, slong n, slong prec){    slong i, j, k, m;    acb_ptr R, S, T, tmp;    acb_t t;    if (n <= 2)    {        if (n >= 1)            acb_zero(Qinv);        if (n == 2)            acb_inv(Qinv + 1, Q + 1, prec);        return;    }    m = n_sqrt(n);    acb_init(t);    R = _acb_vec_init((n - 1) * m);    S = _acb_vec_init(n - 1);    T = _acb_vec_init(n - 1);    acb_zero(Qinv);    acb_inv(Qinv + 1, Q + 1, prec);    _acb_poly_inv_series(Ri(1), Q + 1, FLINT_MIN(Qlen, n) - 1, n - 1, prec);    for (i = 2; i <= m; i++)        _acb_poly_mullow(Ri(i), Ri((i + 1) / 2), n - 1, Ri(i / 2), n - 1, n - 1, prec);    for (i = 2; i < m; i++)        acb_div_ui(Qinv + i, Ri(i) + i - 1, i, prec);    _acb_vec_set(S, Ri(m), n - 1);    for (i = m; i < n; i += m)    {        acb_div_ui(Qinv + i, S + i - 1, i, prec);        for (j = 1; j < m && i + j < n; j++)        {            acb_mul(t, S + 0, Ri(j) + i + j - 1, prec);            for (k = 1; k <= i + j - 1; k++)                acb_addmul(t, S + k, Ri(j) + i + j - 1 - k, prec);            acb_div_ui(Qinv + i + j, t, i + j, prec);        }        if (i + 1 < n)        {            _acb_poly_mullow(T, S, n - 1, Ri(m), n - 1, n - 1, prec);            tmp = S; S = T; T = tmp;        }    }    acb_clear(t);    _acb_vec_clear(R, (n - 1) * m);    _acb_vec_clear(S, n - 1);    _acb_vec_clear(T, n - 1);}

void fmpz_multi_mod_ui(unsigned long * out, fmpz_t in, fmpz_comb_t comb, fmpz_t ** temp){   ulong i, j, k;   ulong n = comb->n;   long log_comb;	ulong size;	mp_limb_t * ptr;	ulong num;	unsigned long num_primes = comb->num_primes;   if (num_primes == 1) // we are reducing modulo a single prime which is assumed to be big enough   {	  if ((long)in[0] > 0L) out[0] = in[1];	  else if ((long)in[0] < 0L) out[0] = comb->primes[0] - in[1];	  else out[0] = 0L;	  return;   }   log_comb = n - 1;      // find level in comb with entries bigger than the input integer	log_comb = 0;	if ((long) in[0] < 0L)	  while ((fmpz_bits(in) >= fmpz_bits(comb->comb[log_comb][0]) - 1) && (log_comb < comb->n - 1)) log_comb++;   else		while (fmpz_cmpabs(in, comb->comb[log_comb][0]) >= 0) log_comb++;   num = (1L<<(n - log_comb - 1));	// set each entry of this level of temp to the input integer	for (i = 0; i < num; i++)   {      fmpz_set(temp[log_comb][i], in);   }   log_comb--;   num *= 2;      // fill in other entries of temp by taking entries of temp at higher level mod pairs from comb	while (log_comb > FLINT_LOG_MULTI_MOD_CUTOFF) // keep going until we reach the basecase   {      for (i = 0, j = 0; i < num; i += 2, j++)      {         fmpz_mod(temp[log_comb][i], temp[log_comb + 1][j], comb->comb[log_comb][i]);         fmpz_mod(temp[log_comb][i+1], temp[log_comb + 1][j], comb->comb[log_comb][i+1]);      }      num *= 2;      log_comb--;   }      // do basecase	num /= 2;   log_comb++;   ulong stride = (1L << (log_comb + 1));   for (i = 0, j = 0; j < num_primes; i++, j += stride)   {	   fmpz_multi_mod_ui_basecase(out + j, temp[log_comb][i], comb->primes + j, FLINT_MIN(stride, num_primes - j));   }}

intmain(void){    int i, result = 1;    flint_rand_t state;    flint_randinit(state);    printf("log_series_monomial_ui....");    fflush(stdout);    for (i = 0; i < 10000; i++)    {        nmod_poly_t A, logA, res;        long n;        mp_limb_t mod;        ulong power;        mp_limb_t coeff;        mod = n_randtest_prime(state, 0);        n = n_randtest(state) % 100;        n = FLINT_MIN(n, mod);        nmod_poly_init(A, mod);        nmod_poly_init(logA, mod);        nmod_poly_init(res, mod);        coeff = n_randlimb(state) % mod;        power = 1 + n_randint(state, 2*n + 1);        nmod_poly_set_coeff_ui(A, 0, 1UL);        nmod_poly_set_coeff_ui(A, power, coeff);        nmod_poly_log_series(logA, A, n);        nmod_poly_log_series_monomial_ui(res, coeff, power, n);        result = nmod_poly_equal(logA, res);        if (!result)        {            printf("FAIL:/n");            printf("n = %ld, mod = %lu/n", n, mod);            printf("power = %lu, coeff = %lu/n", power, coeff);            printf("A: "); nmod_poly_print(A), printf("/n/n");            printf("log(A): "); nmod_poly_print(logA), printf("/n/n");            printf("res: "); nmod_poly_print(res), printf("/n/n");            abort();        }        nmod_poly_clear(A);        nmod_poly_clear(logA);        nmod_poly_clear(res);    }    flint_randclear(state);    printf("PASS/n");    return 0;}

void_acb_hypgeom_coulomb_series(acb_ptr F, acb_ptr G,    acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta,    acb_srcptr z, slong zlen, slong len, slong prec){    acb_ptr t, v;    if (len == 0)        return;    zlen = FLINT_MIN(zlen, len);    if (zlen == 1)    {        acb_hypgeom_coulomb(F, G, Hpos, Hneg, l, eta, z, prec);        if (F != NULL) _acb_vec_zero(F + 1, len - 1);        if (G != NULL) _acb_vec_zero(G + 1, len - 1);        if (Hpos != NULL) _acb_vec_zero(Hpos + 1, len - 1);        if (Hneg != NULL) _acb_vec_zero(Hneg + 1, len - 1);        return;    }    t = _acb_vec_init(len);    v = _acb_vec_init(zlen);    /* copy nonconstant part first to allow aliasing */    acb_zero(v);    _acb_vec_set(v + 1, z + 1, zlen - 1);    acb_hypgeom_coulomb_jet(F, G, Hpos, Hneg, l, eta, z, len, prec);    if (F != NULL)    {        _acb_vec_set(t, F, len);        _acb_poly_compose_series(F, t, len, v, zlen, len, prec);    }    if (G != NULL)    {        _acb_vec_set(t, G, len);        _acb_poly_compose_series(G, t, len, v, zlen, len, prec);    }    if (Hpos != NULL)    {        _acb_vec_set(t, Hpos, len);        _acb_poly_compose_series(Hpos, t, len, v, zlen, len, prec);    }    if (Hneg != NULL)    {        _acb_vec_set(t, Hneg, len);        _acb_poly_compose_series(Hneg, t, len, v, zlen, len, prec);    }    _acb_vec_clear(t, len);    _acb_vec_clear(v, zlen);}