div.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"

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

浏览源代码.

函数
void	lmmp_div_ (mp_ptr dstq, mp_ptr dstr, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
	大数除法和取模操作
mp_limb_t	lmmp_div_s_ (mp_ptr restrict dstq, mp_ptr restrict numa, mp_size_t na, mp_srcptr restrict numb, mp_size_t nb)
	Copyright (C) 2026 HJimmyK(Jericho Knox)

函数说明

lmmp_div_()

void lmmp_div_	(	mp_ptr	dstq,
		mp_ptr	dstr,
		mp_srcptr	numa,
		mp_size_t	na,
		mp_srcptr	numb,
		mp_size_t	nb
	)

大数除法和取模操作

注解

如果dstq不为NULL: [dstq,na-nb+1] = [numa,na] / [numb,nb] (商) 如果dstr不为NULL: [dstr,nb] = [numa,na] mod [numb,nb] (余数)

警告

0<nb<=na, numb[nb-1]!=0, sep(dstq,[numa|numb]), eqsep(dstr,[numa|numb])) 特殊情况: nb==1时, dstq>=numa-1 是允许的 nb==2时, dstq>=numa 是允许的

参数

dstq	商结果输出指针（NULL表示不计算商）
dstr	余数结果输出指针（NULL表示不计算余数）
numa	被除数指针
na	被除数的 limb 长度
numb	除数指针
nb	除数的 limb 长度

在文件 div.c 第 74 行定义.

                                                                                                     {
    if (nb == 1) {
        mp_limb_t rem = lmmp_div_1_(dstq, numa, na, *numb);
        if (dstr)
            *dstr = rem;
    } else if (nb == 2) {
        mp_limb_t brem[2];
        brem[0] = numb[0];
        brem[1] = numb[1];
        lmmp_div_2_(dstq, numa, na, brem);
        if (dstr) {
            dstr[0] = brem[0];
            dstr[1] = brem[1];
        }
    } else {
        int adjust = numa[na - 1] >= numb[nb - 1];
        int cnt = lmmp_leading_zeros_(numb[nb - 1]);
        mp_size_t nq = na + adjust - nb;
        if (nq == 0) {
            if (dstr && dstr != numa)
                lmmp_copy(dstr, numa, nb);
            if (dstq)
                dstq[0] = 0;
            return;
        }
        TEMP_DECL;
 
        if (!dstq)
            dstq = TALLOC_TYPE(na - nb + 1, mp_limb_t);
        dstq[na - nb] = 0;
 
        if (nq >= nb) {
            mp_ptr restrict numa2 = TALLOC_TYPE(na + 1, mp_limb_t);
            mp_ptr numb2;
            if (cnt) {
                numa2[na] = lmmp_shl_(numa2, numa, na, cnt);
                numb2 = TALLOC_TYPE(nb, mp_limb_t);
                lmmp_shl_(numb2, numb, nb, cnt);
            } else {
                numa2[na] = 0;
                lmmp_copy(numa2, numa, na);
                numb2 = (mp_ptr)numb;
            }
 
            mp_limb_t inv21 = lmmp_inv_2_1_(numb2[nb - 1], numb2[nb - 2]);
            na += adjust;
 
            if (nb < DIV_DIVIDE_THRESHOLD)
                lmmp_div_basecase_(dstq, numa2, na, numb2, nb, inv21);
            else if (nb < DIV_MULINV_L_THRESHOLD || na < 2 * DIV_MULINV_N_THRESHOLD)
                lmmp_div_divide_(dstq, numa2, na, numb2, nb, inv21);
            else {
                mp_limb_t ni = lmmp_div_inv_size_(nq, nb);
                mp_ptr invappr = TALLOC_TYPE(ni, mp_limb_t);
                lmmp_inv_prediv_(invappr, numb2, nb, ni);
                lmmp_div_mulinv_(dstq, numa2, na, numb2, nb, invappr, ni);
            }
 
            if (dstr) {
                if (cnt)
                    lmmp_shr_(dstr, numa2, nb, cnt);
                else
                    lmmp_copy(dstr, numa2, nb);
            }
        } else {
            // nq=na-nb+adj<nb
            //-> na+adj>=2nq+1
            mp_size_t ni = nb - nq;
            mp_ptr numa2, numb2;
            mp_ptr tp = TALLOC_TYPE(nb, mp_limb_t);
            mp_limb_t cy;
 
            numa2 = TALLOC_TYPE(nq * 2 + 1, mp_limb_t);
            if (cnt) {
                numb2 = TALLOC_TYPE(nq, mp_limb_t);
                lmmp_shl_(numb2, numb + ni, nq, cnt);
                numb2[0] |= numb[ni - 1] >> (LIMB_BITS - cnt);
                cy = lmmp_shl_(numa2, numa + na - 2 * nq, 2 * nq, cnt);
                if (adjust) {
                    numa2[2 * nq] = cy;
                    ++numa2;  // numa2[0] is as significant as numa[ni=na-2nq+adjust]
                } else
                    numa2[0] |= numa[na - 2 * nq - 1] >> (LIMB_BITS - cnt);
            } else {
                numb2 = (mp_ptr)numb + ni;
                lmmp_copy(numa2, numa + na - 2 * nq, 2 * nq);
                if (adjust) {
                    numa2[2 * nq] = 0;
                    ++numa2;
                }
            }
 
            // now: 0<=numa2<B^2nq, B^nq/2<=numb2<B^nq, and 0<=numa2/numb2<B^nq
            // ignored bits could be seen as fraction part of numa and numb
            // we can prove:  Q<=Qh<=Q+2
            // where Q=floor(numa/numb) is the real quotient
            // Qh=floor(floor(numa)/floor(numb)) as below
 
            if (nq == 1) {
                lmmp_div_1_s_(dstq, numa2, 2, *numb2);
            } else if (nq == 2) {
                lmmp_div_2_s_(dstq, numa2, 4, numb2);
            } else {
                mp_limb_t inv21 = lmmp_inv_2_1_(numb2[nq - 1], numb2[nq - 2]);
 
                if (nq < DIV_DIVIDE_THRESHOLD)
                    lmmp_div_basecase_(dstq, numa2, 2 * nq, numb2, nq, inv21);
                else if (nq < DIV_MULINV_N_THRESHOLD)
                    lmmp_div_divide_(dstq, numa2, 2 * nq, numb2, nq, inv21);
                else {
                    mp_limb_t ni = lmmp_div_inv_size_(nq, nq);
                    mp_ptr invappr = tp;
                    lmmp_inv_prediv_(invappr, numb2, nq, ni);
                    lmmp_div_mulinv_(dstq, numa2, 2 * nq, numb2, nq, invappr, ni);
                }
            }
            /*
            true remainder = partial remainder - quotient * ignored divisor limbs
 
            Multiply the first ignored divisor limb by the most significant
            quotient limb.  If that product is > the partial remainder's
            most significant limb, we know the quotient is too large.  This
            test quickly catches most cases where the quotient is too large;
            it catches all cases where the quotient is 2 too large.*/
 
            mp_limb_t x;
            if (cnt) {
                mp_limb_t dl;
                if (ni < 2)
                    dl = 0;
                else
                    dl = numb[ni - 2];
                x = (numb[ni - 1] << cnt) | (dl >> (LIMB_BITS - cnt));
            } else
                x = numb[ni - 1];
            mp_limb_t h = (x >> LIMB_BITS / 2) * (dstq[nq - 1] >> LIMB_BITS / 2);
            mp_limb_t rnb = 0;  // remainder[nb]
            mp_size_t nr = nq;  // remainder=rnb:[numa2,nr]:[...,ni]
 
            if (h > numa2[nq - 1]) {
                lmmp_dec(dstq);
                rnb = lmmp_add_n_(numa2, numa2, numb2, nq);
            }
 
            // if cnt, recover the shift of partial remainder
            // and remove the effect of the partial-ignored numa[ni-1] and numb[ni-1]
            if (cnt) {
                numa2[nq] = rnb;
                ++nr;
                --ni;
                lmmp_shl_(numa2, numa2, nr, LIMB_BITS - cnt);
                numa2[0] |= numa[ni] & (LIMB_MAX >> cnt);
                cy = lmmp_submul_1_(numa2, dstq, nq, numb[ni] & (LIMB_MAX >> cnt));
                rnb = -(numa2[nq] < cy);
                numa2[nq] -= cy;
            }
 
            if (ni == 0) {
                if (dstr) {
                    if (rnb)
                        lmmp_add_n_(dstr, numa2, numb, nr);
                    else
                        lmmp_copy(dstr, numa2, nr);
                }
            } else {
                tp[nb - 1] = 0;
                if (ni < nq)
                    lmmp_mul_(tp, dstq, nq, numb, ni);
                else
                    lmmp_mul_(tp, numb, ni, dstq, nq);
 
                if (dstr) {
                    mp_ptr remptr = dstr == numb ? tp : dstr;
                    cy = lmmp_sub_n_(remptr, numa, tp, ni);
                    rnb -= lmmp_sub_nc_(remptr + ni, numa2, tp + ni, nr, cy);
                    if (rnb)
                        lmmp_add_n_(dstr, remptr, numb, nb);
                    else if (dstr != remptr)
                        lmmp_copy(dstr, remptr, nb);
                } else {
                    int hcmp = lmmp_cmp_(numa2, tp + ni, nr);
                    if (hcmp < 0)
                        --rnb;
                    else if (hcmp == 0)
                        rnb -= (lmmp_cmp_(numa, tp, ni) < 0);
                }
            }
 
            if (rnb)
                lmmp_dec(dstq);
        }
 
        TEMP_FREE;
    }
}

引用了 DIV_DIVIDE_THRESHOLD, DIV_MULINV_L_THRESHOLD, DIV_MULINV_N_THRESHOLD, LIMB_BITS, LIMB_MAX, lmmp_add_n_(), lmmp_cmp_(), lmmp_copy, lmmp_dec, lmmp_div_1_(), lmmp_div_1_s_(), lmmp_div_2_(), lmmp_div_2_s_(), lmmp_div_basecase_(), lmmp_div_divide_(), lmmp_div_inv_size_(), lmmp_div_mulinv_(), lmmp_inv_2_1_(), lmmp_inv_prediv_(), lmmp_leading_zeros_, lmmp_mul_(), lmmp_shl_(), lmmp_shr_(), lmmp_sub_n_(), lmmp_sub_nc_(), lmmp_submul_1_(), n, numb, TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

被这些函数引用 lmmp_bninv_(), lmmp_gcd_basecase_(), lmmp_gcd_lehmer_(), lmmp_trialdiv_() , 以及 try_div_().

函数调用图:

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

lmmp_div_s_()

mp_limb_t lmmp_div_s_	(	mp_ptr restrict	dstq,
		mp_ptr restrict	numa,
		mp_size_t	na,
		mp_srcptr restrict	numb,
		mp_size_t	nb
	)

Copyright (C) 2026 HJimmyK(Jericho Knox)

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

在文件 div.c 第 22 行定义.

  {
    TEMP_DECL;
    mp_limb_t nq = na - nb;
    mp_limb_t qh;
    if (nq == 0) {
        qh = lmmp_cmp_(numa, numb, nb) >= 0;
        if (qh)
            lmmp_sub_n_(numa, numa, numb, nb);
    } else if (nb == 1) {
        qh = lmmp_div_1_s_(dstq, numa, na, *numb);
    } else if (nb == 2) {
        qh = lmmp_div_2_s_(dstq, numa, na, numb);
    } else if (nq < nb) {
        qh = lmmp_div_s_(dstq, numa + na - 2 * nq, 2 * nq, numb + nb - nq, nq);
 
        mp_ptr restrict tp = TALLOC_TYPE(nb, mp_limb_t);
        if (nq > nb - nq)
            lmmp_mul_(tp, dstq, nq, numb, nb - nq);
        else
            lmmp_mul_(tp, numb, nb - nq, dstq, nq);
 
        mp_limb_t cy = lmmp_sub_n_(numa, numa, tp, nb);
        if (qh)
            cy += lmmp_sub_n_(numa + nq, numa + nq, numb, nb - nq);
 
        while (cy) {
            qh -= lmmp_sub_1_(dstq, dstq, nq, 1);
            cy -= lmmp_add_n_(numa, numa, numb, nb);
        }
    } else {
        mp_limb_t inv21 = lmmp_inv_2_1_(numb[nb - 1], numb[nb - 2]);
        if (nb < DIV_DIVIDE_THRESHOLD)
            qh = lmmp_div_basecase_(dstq, numa, na, numb, nb, inv21);
        else if (nb < DIV_MULINV_L_THRESHOLD || na < 2 * DIV_MULINV_N_THRESHOLD)
            qh = lmmp_div_divide_(dstq, numa, na, numb, nb, inv21);
        else {
            mp_limb_t ni = lmmp_div_inv_size_(nq, nb);
            mp_ptr restrict invappr = TALLOC_TYPE(ni, mp_limb_t);
            lmmp_inv_prediv_(invappr, numb, nb, ni);
            qh = lmmp_div_mulinv_(dstq, numa, na, numb, nb, invappr, ni);
        }
    }
    TEMP_FREE;
    return qh;
}

引用了 DIV_DIVIDE_THRESHOLD, DIV_MULINV_L_THRESHOLD, DIV_MULINV_N_THRESHOLD, lmmp_add_n_(), lmmp_cmp_(), lmmp_div_1_s_(), lmmp_div_2_s_(), lmmp_div_basecase_(), lmmp_div_divide_(), lmmp_div_inv_size_(), lmmp_div_mulinv_(), lmmp_div_s_(), lmmp_inv_2_1_(), lmmp_inv_prediv_(), lmmp_mul_(), lmmp_sub_1_(), lmmp_sub_n_(), n, numb, TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

被这些函数引用 lmmp_div_s_().

函数调用图:

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