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