divexact.c 文件参考

#include "../../../include/lammp/impl/inlines.h"
#include "../../../include/lammp/impl/mparam.h"
#include "../../../include/lammp/impl/tmp_alloc.h"
#include "../../../include/lammp/lmmpn.h"
#include "../../../include/lammp/numth.h"
#include "../../../include/lammp/impl/longlong.h"

divexact.c 的引用(Include)关系图:

浏览源代码.

宏定义
#define	c   (tp)
#define	l   (tp + dn)
#define	scratch   (tp + 2 * dn)

函数
void	lmmp_divexact_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_srcptr restrict dp, mp_size_t dn)
void	lmmp_divexact_1_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_limb_t d, mp_limb_t dinv)
	Copyright (C) 2026 HJimmyK(Jericho Knox)
void	lmmp_divexact_2_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_srcptr restrict dp, mp_srcptr restrict dinv)
void	lmmp_divexact_basecase_ (mp_ptr dst, mp_ptr np, mp_size_t nn, mp_srcptr restrict dp, mp_size_t dn, mp_limb_t dinv)
void	lmmp_divexact_divide_ (mp_ptr restrict dst, mp_srcptr restrict np, mp_size_t nn, mp_srcptr restrict dp, mp_size_t dn)
void	lmmp_divexact_unbalanced_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_srcptr restrict dp, mp_size_t dn, mp_ptr restrict dinv)
static void	lmmp_mulhi_n_ (mp_ptr dst, mp_srcptr restrict lop, mp_srcptr restrict numa, mp_srcptr restrict numb, mp_size_t n, mp_ptr tp)
	已知乘积的低位，计算乘积的高位
static void	lmmp_mullo_n_ (mp_ptr restrict dst, mp_srcptr restrict numa, mp_srcptr restrict numb, mp_size_t n, mp_ptr restrict tp)

宏定义说明

c

#define c   (tp)

l

#define l   (tp + dn)

scratch

#define scratch   (tp + 2 * dn)

函数说明

lmmp_divexact_()

void lmmp_divexact_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_srcptr restrict	dp,
		mp_size_t	dn
	)

在文件 divexact.c 第 273 行定义.

                                                                                                 {
    lmmp_param_assert(dst != NULL && np != NULL && dp != NULL);
    lmmp_param_assert(dn > 0);
    lmmp_param_assert(nn >= dn);
    lmmp_param_assert((dp[0] & 1) != 0);
 
    if (dn == 1) {
        mp_limb_t dinv = lmmp_binvert_ulong_(dp[0]);
        lmmp_divexact_1_(dst, np, nn, dp[0], dinv);
    } else if (dn == 2) {
        mp_limb_t dinv[2];
        lmmp_binvert_2_(dinv, dp);
        lmmp_divexact_2_(dst, np, nn, dp, dinv);
    } else {
        mp_size_t qn = nn - dn + 1;
        if (nn < DIVEXACT_NN_THRESHOLD && dn < DIVEXACT_BASECASE_THRESHOLD) {
            mp_limb_t dinv = lmmp_binvert_ulong_(dp[0]);
            TEMP_DECL;
            mp_ptr restrict tp = TALLOC_TYPE(nn, mp_limb_t);
            lmmp_copy(tp, np, nn);
            lmmp_divexact_basecase_(dst, tp, nn, dp, dn, dinv);
            TEMP_FREE;
        } else if (qn < dn) {
            TEMP_DECL;
            mp_ptr restrict tp;
            if (dst == np) {
                tp = TALLOC_TYPE(nn, mp_limb_t);
                lmmp_copy(tp, np, nn);
                np = tp;
            }
            lmmp_divexact_divide_(dst, np, nn, dp, dn);
            TEMP_FREE;
        } else {
            lmmp_divexact_unbalanced_(dst, np, nn, dp, dn, NULL);
        }
    }
}

引用了 DIVEXACT_BASECASE_THRESHOLD, DIVEXACT_NN_THRESHOLD, lmmp_binvert_2_(), lmmp_binvert_ulong_(), lmmp_copy, lmmp_divexact_1_(), lmmp_divexact_2_(), lmmp_divexact_basecase_(), lmmp_divexact_divide_(), lmmp_divexact_unbalanced_(), lmmp_param_assert, n, TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

函数调用图:

lmmp_divexact_1_()

void lmmp_divexact_1_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_limb_t	d,
		mp_limb_t	dinv
	)

Copyright (C) 2026 HJimmyK(Jericho Knox)

精确除法（[dst,nn]=[np,nn]/d，且余数必须为0）

This file is part of LAMMP.

LAMMP 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.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 divexact.c 第 24 行定义.

                                                                                           {
    lmmp_param_assert(d % 2 == 1);
    lmmp_param_assert(dinv * d == 1);
    lmmp_param_assert(nn > 0);
    lmmp_debug_assert(dst != NULL && np != NULL);
    mp_limb_t c = 0;
    mp_limb_t l, s, lo, hi, q;
    mp_size_t i;
    for (i = 0; i < nn - 1; i++) {
        s = np[i];
        l = s - c;
        c = l > s;
        q = l * dinv;
        dst[i] = q;
        _umul64to128_(q, d, &lo, &hi);
        c += hi;
    }
    s = np[i];
    l = s - c;
    c = l > s;
    q = l * dinv;
    dst[i] = q;
    lmmp_debug_assert(c == 0);
}

引用了 _umul64to128_(), c, l, lmmp_debug_assert, lmmp_param_assert, lo, n , 以及 s.

被这些函数引用 lmmp_divexact_().

函数调用图:

这是这个函数的调用关系图:

lmmp_divexact_2_()

void lmmp_divexact_2_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_srcptr restrict	dp,
		mp_srcptr restrict	dinv
	)

在文件 divexact.c 第 49 行定义.

                                                                                                              {
    lmmp_param_assert(dp[0] % 2 == 1);
    lmmp_param_assert(nn > 1);
    lmmp_debug_assert(dst != NULL && np != NULL);
    mp_limb_t c[2] = {0, 0};
    mp_limb_t l[2], s[2], q[2], t[4];
    mp_size_t i;
    mp_limb_t d0 = dp[0], d1 = dp[1];
    mp_limb_t ddinv0 = dinv[0], ddinv1 = dinv[1];
 
    if (nn % 2 == 0) {
        for (i = 0; i < nn - 2; i += 2) {
            s[0] = np[i];
            s[1] = np[i + 1];
            _u128sub(l, s, c);
            c[0] = _u128cmp(s, l) ? 1 : 0;
            c[1] = 0;
            _umul128to128_(l[1], l[0], ddinv1, ddinv0, q);
            dst[i] = q[0];
            dst[i + 1] = q[1];
            _umul128to256_(q[1], q[0], d1, d0, t);
            _u128add(c, c, t + 2);
        }
        s[0] = np[i];
        s[1] = np[i + 1];
        _u128sub(l, s, c);
        c[0] = _u128cmp(s, l);
        _umul128to128_(l[1], l[0], ddinv1, ddinv0, q);
        dst[i] = q[0];
        lmmp_debug_assert(c[0] == 0);
    } else {
        i = 0;
        if (nn >= 5) {
            for (; i < nn - 4; i += 2) {
                s[0] = np[i];
                s[1] = np[i + 1];
                _u128sub(l, s, c);
                c[0] = _u128cmp(s, l);
                c[1] = 0;
                _umul128to128_(l[1], l[0], ddinv1, ddinv0, q);
                dst[i] = q[0];
                dst[i + 1] = q[1];
                _umul128to256_(q[1], q[0], d1, d0, t);
                _u128add(c, c, t + 2);
            }
        }
        s[0] = np[i];
        s[1] = np[i + 1];
        _u128sub(l, s, c);
        c[0] = _u128cmp(s, l);
        _umul128to128_(l[1], l[0], ddinv1, ddinv0, q);
        dst[i] = q[0];
        dst[i + 1] = q[1];
        lmmp_debug_assert(c[0] == 0);
    }
}

引用了 _u128add, _u128cmp, _u128sub, _umul128to128_(), _umul128to256_(), c, l, lmmp_debug_assert, lmmp_param_assert, n, s , 以及 t.

被这些函数引用 lmmp_divexact_().

函数调用图:

这是这个函数的调用关系图:

lmmp_divexact_basecase_()

void lmmp_divexact_basecase_	(	mp_ptr	dst,
		mp_ptr	np,
		mp_size_t	nn,
		mp_srcptr restrict	dp,
		mp_size_t	dn,
		mp_limb_t	dinv
	)

在文件 divexact.c 第 214 行定义.

                                                                                                                       {
    lmmp_param_assert(dst != NULL && np != NULL && dp != NULL);
    lmmp_param_assert(dn > 0);
    lmmp_param_assert(nn >= dn);
    lmmp_param_assert((dp[0] & 1) != 0);
    lmmp_param_assert((dp[0] * dinv) == 1);
    mp_size_t i;
    mp_limb_t q;
    mp_limb_t hi;
    mp_size_t qn = nn - dn + 1;
    for (i = 0; i < qn; i++) {
        q = dinv * np[i];
        hi = lmmp_submul_1_(np + i, dp, dn, q);
        lmmp_debug_assert(np[i] == 0);
#if LAMMP_DEBUG_ASSERT_CHECK == 1
        if (hi && i + dn < nn) {
            lmmp_dec_1(np + i + dn, hi);
        } else {
            lmmp_debug_assert(hi == 0);
        }
#else
        if (hi) {
            lmmp_dec_1(np + i + dn, hi);
        }
#endif
        dst[i] = q;
    }
}

引用了 lmmp_debug_assert, lmmp_dec_1, lmmp_param_assert, lmmp_submul_1_() , 以及 n.

被这些函数引用 lmmp_divexact_().

函数调用图:

这是这个函数的调用关系图:

lmmp_divexact_divide_()

void lmmp_divexact_divide_	(	mp_ptr restrict	dst,
		mp_srcptr restrict	np,
		mp_size_t	nn,
		mp_srcptr restrict	dp,
		mp_size_t	dn
	)

在文件 divexact.c 第 243 行定义.

                                                                                                                          {
    lmmp_param_assert(dst != NULL && np != NULL && dp != NULL);
    lmmp_param_assert(dn > 0);
    lmmp_param_assert(nn >= dn);
    lmmp_param_assert((dp[0] & 1) != 0);
 
    mp_size_t qn = nn - dn + 1;
    lmmp_param_assert(qn <= dn);
    mp_size_t rn = qn - qn / 2;
    qn = qn / 2;
 
    TEMP_DECL;
    mp_ptr restrict dinv = TALLOC_TYPE(rn, mp_limb_t);
    mp_ptr restrict tp = TALLOC_TYPE(3 * rn, mp_limb_t);
    lmmp_binvert_n_dc_(dinv, dp, rn, tp);
    lmmp_mullo_n_(dst, np, dinv, rn, tp);
    if (qn == 0) {
        TEMP_FREE;
        return;
    }
 
    lmmp_debug_assert(3 * rn >= qn + rn);
    lmmp_mul_(tp, dp, qn + rn, dst, rn);
 
    // 即使存在错位，高位的np也一定可以将其约减去，因此无需进行任何处理
    lmmp_sub_n_(tp, np + rn, tp + rn, qn);
    lmmp_mullo_n_(dst + rn, tp, dinv, qn, tp + rn);
    TEMP_FREE;
}

引用了 lmmp_binvert_n_dc_(), lmmp_debug_assert, lmmp_mul_(), lmmp_mullo_n_(), lmmp_param_assert, lmmp_sub_n_(), n, TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

被这些函数引用 lmmp_divexact_().

函数调用图:

这是这个函数的调用关系图:

lmmp_divexact_unbalanced_()

void lmmp_divexact_unbalanced_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_srcptr restrict	dp,
		mp_size_t	dn,
		mp_ptr restrict	dinv
	)

在文件 divexact.c 第 160 行定义.

  {
    lmmp_param_assert(dst != NULL && np != NULL && dp != NULL);
    lmmp_param_assert(dn > 0);
    lmmp_param_assert(nn >= dn);
    lmmp_param_assert((dp[0] & 1) != 0);
 
    TEMP_DECL;
    mp_ptr restrict tp = TALLOC_TYPE(4 * dn, mp_limb_t);
    if (dinv == NULL) {
        dinv = TALLOC_TYPE(dn, mp_limb_t);
        lmmp_binvert_n_dc_(dinv, dp, dn, tp);
    }
 
#define c       (tp)          // [tp, dn]
#define l       (tp + dn)     // [tp + dn, dn]
#define scratch (tp + 2 * dn) // [tp + 2 * dn, 2*dn]
 
    mp_size_t i = 0, qn = nn - dn + 1;
    mp_limb_t ca;
    lmmp_zero(c, dn);
    if (qn > dn) {
        for (; i < qn - dn; i += dn) {
            lmmp_sub_n_(l, np + i, c, dn);
            ca = lmmp_cmp_(l, np + i, dn) == 1 ? 1 : 0;
            /*
            FIXME: 这里的循环中，第二个乘数dinv以及dp，始终保持不变
                   在拥有可以惰性初始化的FFT算法的情况下，可以节省numa的正变换
                   在循环的情况下，这将会有可观的性能提升
            */
            lmmp_mullo_n_(dst + i, l, dinv, dn, scratch);
            lmmp_mulhi_n_(scratch + dn, l, dst + i, dp, dn, scratch);
            //lmmp_mul_n_(scratch, dst + i, dp, dn);
            if (ca) {
                lmmp_add_1_(c, scratch + dn, dn, 1);
            } else {
                lmmp_copy(c, scratch + dn, dn);
            }
        }
    }
    lmmp_sub_n_(l, np + i, c, dn);
    lmmp_mullo_n_(dst + i, l, dinv, qn - i, scratch);
    TEMP_FREE;
#undef c
#undef l
#undef scratch
}

引用了 c, l, lmmp_add_1_(), lmmp_binvert_n_dc_(), lmmp_cmp_(), lmmp_copy, lmmp_mulhi_n_(), lmmp_mullo_n_(), lmmp_param_assert, lmmp_sub_n_(), lmmp_zero, n, scratch, TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

被这些函数引用 lmmp_divexact_().

函数调用图:

这是这个函数的调用关系图:

lmmp_mulhi_n_()

static void lmmp_mulhi_n_ ( mp_ptr  dst, mp_srcptr restrict  lop, mp_srcptr restrict  numa, mp_srcptr restrict  numb, mp_size_t  n, mp_ptr  tp  )

inlinestatic

已知乘积的低位，计算乘积的高位

参数

dst	乘积的高位，长度为 n
lop	乘积的低位，长度为 n
numa	乘数a
numb	乘数b
n	乘数a和乘数b的 limb 长度
tp	临时空间（应需要2*n个limb）

警告

eqsep(dst,tp)

在文件 divexact.c 第 130 行定义.

  {
    if (n < MULHI_MERSENNE_THRESHOLD) {
        lmmp_mul_n_(tp, numa, numb, n);
        lmmp_copy(dst, tp + n, n);
    } else {
        mp_limb_t c;
        mp_size_t m = lmmp_fft_next_size_((n * 2 + 1) >> 1);
        lmmp_debug_assert(n * 2 > m && m >= n);
        lmmp_mul_mersenne_(tp, m, numa, n, numb, n);
        c = lmmp_sub_(tp, tp, m, lop, n);
        if (c != 0)
            lmmp_dec(tp);
        if (m == n) {
            lmmp_copy(dst, tp, n);
        } else {
            mp_size_t fn = m - n;
            mp_size_t sn = n - fn;
            lmmp_copy(dst, tp + n, fn);
            lmmp_copy(dst + fn, tp, sn);
        }
    }
}

引用了 c, lmmp_copy, lmmp_debug_assert, lmmp_dec, lmmp_fft_next_size_(), lmmp_mul_mersenne_(), lmmp_mul_n_, lmmp_sub_(), MULHI_MERSENNE_THRESHOLD, n, numb , 以及 tp.

被这些函数引用 lmmp_divexact_unbalanced_().

函数调用图:

这是这个函数的调用关系图:

lmmp_mullo_n_()

static void lmmp_mullo_n_ ( mp_ptr restrict  dst, mp_srcptr restrict  numa, mp_srcptr restrict  numb, mp_size_t  n, mp_ptr restrict  tp  )

inlinestatic

在文件 divexact.c 第 106 行定义.

  {
    if (n < MULLO_DC_THRESHOLD) {
        lmmp_mullo_dc_(dst, numa, numb, tp, n);
    } else {
        lmmp_mullo_fft_(dst, numa, numb, n, tp);
    }
}

引用了 lmmp_mullo_dc_(), lmmp_mullo_fft_(), MULLO_DC_THRESHOLD, n, numb , 以及 tp.

被这些函数引用 lmmp_divexact_divide_() , 以及 lmmp_divexact_unbalanced_().

函数调用图:

这是这个函数的调用关系图: