div_mulinv.c 文件参考

#include "../../include/lammp/impl/mparam.h"
#include "../../include/lammp/impl/tmp_alloc.h"
#include "../../include/lammp/lmmpn.h"

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

浏览源代码.

函数
mp_limb_t	lmmp_div_mulinv_ (mp_ptr dstq, mp_ptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb, mp_srcptr invappr, mp_size_t ni)
	乘法逆元除法
void	lmmp_inv_prediv_ (mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_size_t ni)
	除法前的逆元预计算，[dst,ni] = invappr( (ni+1 MSLs of numa) + 1 ) / B

函数说明

lmmp_div_mulinv_()

mp_limb_t lmmp_div_mulinv_	(	mp_ptr	dstq,
		mp_ptr	numa,
		mp_size_t	na,
		mp_srcptr	numb,
		mp_size_t	nb,
		mp_srcptr	invappr,
		mp_size_t	ni
	)

乘法逆元除法

参数

dstq	输出商的缓冲区，长度至少为na-nb
numa	输入被除数（长度na），运算后存储余数（长度nb）
na	被除数的 limb 长度
numb	输入除数，长度为nb
nb	除数的 limb 长度
invappr	预计算的近似逆元，长度为ni
ni	预计算逆元的 limb 长度

返回

商的最高位（qh）

警告

na>=nb>=ni>0, MSB(numb)=1, [invappr,ni]=inv_prediv([numb,nb]), sep(dstq,numa,numb,invappr))

注解

qh:[dstq,na-1]=[numa,na] div x, [numa,1]=[numa,na] mod x, return qh

在文件 div_mulinv.c 第 36 行定义.

  {
    lmmp_param_assert(na >= nb && nb >= ni);
    lmmp_param_assert(ni > 0);
    lmmp_param_assert(numb[nb - 1] >= LIMB_B_2);
    mp_size_t nq = na - nb, ntp = LMMP_MIN(ni, nq) + nb;
    mp_limb_t qh;
    TEMP_DECL;
    mp_ptr tp = TALLOC_TYPE(ntp, mp_limb_t);
 
    numa += nq;
    dstq += nq;
 
    qh = lmmp_cmp_(numa, numb, nb) >= 0;
    if (qh) {
        lmmp_sub_n_(numa, numa, numb, nb);
    }
    while (nq) {
        if (nq < ni) {
            invappr += ni - nq;
            ni = nq;
        }
        numa -= ni;
        dstq -= ni;
        nq -= ni;
 
        lmmp_mul_n_(tp, numa + nb, invappr, ni);
        lmmp_assert(lmmp_add_n_(dstq, tp + ni, numa + nb, ni) == 0);
 
        mp_size_t mn, wn;
        mp_limb_t cy;
 
        if (nb < DIV_MULINV_MODM_THRESHOLD || (mn = lmmp_fft_next_size_(nb + 1)) >= nb + ni) {
            lmmp_mul_(tp, numb, nb, dstq, ni);  // nb+ni limbs, high 'ni' cancels
        } else {
            // 0<wn<ni<=nb<mn<nb+ni
            wn = nb + ni - mn;
 
            // x=b*q
            // tp=x mod 2^mn-1
            lmmp_mul_mersenne_(tp, mn, numb, nb, dstq, ni);
 
            // tp-=ah:0 mod B^mn-1, if result=0, represent it as B^mn-1
            cy = lmmp_sub_nc_(tp, tp, numa + mn, wn, 1);
            if (cy)
                cy = lmmp_sub_1_(tp + wn, tp + wn, mn - wn, 1);
            if (!cy)
                lmmp_inc(tp);
 
            // if al<<tp,
            if (lmmp_cmp_(numa + nb, tp + nb, mn - nb) < 0) {
                // maybe ah=xh+1 and al<<xl,
                //  so we subtracted 1 too much when tp-=ah,
                //  now tp=xl-1 mod B^mn-1, and 0<=al<<xl-1<B^mn-1, so tp=xl-1
                // or ah=xh and al>=xl,
                //  tp=xl mod B^mn-1, the only possibility is we represented xl=0 as tp=B^mn-1
                // whatever, just inc and then tp=xl
                tp[mn] = 0;  // set a limit
                lmmp_inc(tp);
            }
        }
 
        mp_limb_t r = numa[nb] - tp[nb];
        cy = lmmp_sub_n_(numa, numa, tp, nb);
 
        while ((r -= cy) || lmmp_cmp_(numa, numb, nb) >= 0) {
            lmmp_inc(dstq);
            cy = lmmp_sub_n_(numa, numa, numb, nb);
        }
    }
    TEMP_FREE;
    return qh;
}

引用了 DIV_MULINV_MODM_THRESHOLD, LIMB_B_2, lmmp_add_n_(), lmmp_assert, lmmp_cmp_(), lmmp_fft_next_size_(), lmmp_inc, LMMP_MIN, lmmp_mul_(), lmmp_mul_mersenne_(), lmmp_mul_n_(), lmmp_param_assert, lmmp_sub_1_(), lmmp_sub_n_(), lmmp_sub_nc_(), TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

被这些函数引用 lmmp_div_(), lmmp_div_s_() , 以及 lmmp_to_str_divide_().

函数调用图:

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

lmmp_inv_prediv_()

void lmmp_inv_prediv_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na,
		mp_size_t	ni
	)

除法前的逆元预计算，[dst,ni] = invappr( (ni+1 MSLs of numa) + 1 ) / B

参数

dst	输出预计算逆元的缓冲区，长度为ni
numa	输入操作数，长度为na
na	输入操作数的 limb 长度
ni	预计算逆元的目标尺寸

警告

na>=ni>0, MSB(numa)=1, eqsep(dst,numa)

注解

if (ni=na) [dst,na] = (B^(2*na)-1) / [numa,na] - B^na

在文件 div_mulinv.c 第 11 行定义.

                                                                              {
    lmmp_param_assert(na >= ni);
    lmmp_param_assert(ni > 0);
    lmmp_param_assert(numa[na - 1] >= LIMB_B_2);
    TEMP_DECL;
    mp_limb_t cy;
    mp_ptr tp = TALLOC_TYPE(ni + 1, mp_limb_t);
 
    if (na == ni) {
        lmmp_copy(tp + 1, numa, ni);
        tp[0] = 1;
        cy = 0;
    } else {
        cy = lmmp_add_1_(tp, numa + na - (ni + 1), ni + 1, 1);
    }
    if (cy)
        lmmp_zero(dst, ni);
    else {
        mp_ptr invappr = TALLOC_TYPE(ni + 1, mp_limb_t);
        lmmp_invappr_(invappr, tp, ni + 1);
        lmmp_copy(dst, invappr + 1, ni);
    }
    TEMP_FREE;
}

引用了 LIMB_B_2, lmmp_add_1_(), lmmp_copy, lmmp_invappr_(), lmmp_param_assert, lmmp_zero, TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

被这些函数引用 lmmp_div_(), lmmp_div_s_() , 以及 lmmp_to_str_().

函数调用图:

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