inv.c 文件参考

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

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

浏览源代码.

函数
void	lmmp_inv_ (mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_size_t nf)
	大数求逆操作 [dst,na+nf+1] = (B^(2*(na+nf)) - 1) / ([numa,na]*B^nf) + [0|-1]
void	lmmp_inv_basecase_ (mp_ptr dst, mp_srcptr numa, mp_size_t na)
	近似逆元计算
void	lmmp_invappr_ (mp_ptr dst, mp_srcptr numa, mp_size_t na)
	近似逆元计算 (invappr)
void	lmmp_invappr_newton_ (mp_ptr dst, mp_srcptr numa, mp_size_t na)
	近似逆元计算（牛顿迭代法）

函数说明

lmmp_inv_()

void lmmp_inv_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na,
		mp_size_t	nf
	)

大数求逆操作 [dst,na+nf+1] = (B^(2*(na+nf)) - 1) / ([numa,na]*B^nf) + [0|-1]

参数

dst	逆元结果输出指针
numa	源操作数指针
na	操作数的 limb 长度
nf	精度因子

警告

na>0, numa[na-1]!=0, eqsep(dst,numa)

在文件 inv.c 第 152 行定义.

                                                                       {
    lmmp_param_assert(na > 0);
    lmmp_param_assert(numa[na - 1] != 0);
    mp_limb_t high = numa[na - 1];
    int nsh = lmmp_leading_zeros_(high);
    TEMP_DECL;
    if (dst == numa || nsh || nf) {
        nf += nsh != 0;
        mp_ptr numa2 = TALLOC_TYPE(na + nf, mp_limb_t);
        lmmp_zero(numa2, nf);
        if (nsh)
            lmmp_shl_(numa2 + nf, numa, na, nsh);
        else
            lmmp_copy(numa2 + nf, numa, na);
        numa = numa2;
    }
    lmmp_invappr_(dst, numa, na + nf);
    if (nsh)
        lmmp_shr_c_(dst, dst, na + nf, LIMB_BITS - nsh, (mp_limb_t)1 << nsh);
    else
        dst[na + nf] = 1;
    TEMP_FREE;
}

引用了 LIMB_BITS, lmmp_copy, lmmp_invappr_(), lmmp_leading_zeros_(), lmmp_param_assert, lmmp_shl_(), lmmp_shr_c_(), lmmp_zero, TALLOC_TYPE, TEMP_DECL , 以及 TEMP_FREE.

函数调用图:

lmmp_inv_basecase_()

void lmmp_inv_basecase_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na
	)

近似逆元计算

参数

dst	输出结果缓冲区，长度为na
numa	输入操作数，长度为na
na	输入操作数的 limb 长度

警告

na>0, MSB(numa)=1, sep(dst,numa)

返回

无返回值，结果存储在dst中，[dst,na]=(B^(2*na)-1)/[numa,na] - B^na

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

                                                                  {
    lmmp_param_assert(na > 0);
    lmmp_param_assert(numa[na - 1] >= LIMB_B_2);
    if (na == 1)
        *dst = lmmp_inv_1_(*numa);
    else {
        TEMP_DECL;
        mp_ptr xp = TALLOC_TYPE(2 * na, mp_limb_t);
        mp_size_t i = na;
        do {
            xp[--i] = LIMB_MAX;
        } while (i);
        lmmp_not_(xp + na, numa, na);
        //[xp,2*na] = B^(2*na)-1 - [numa,na]*B^na
 
        if (na == 2) {
            lmmp_div_2_s_(dst, xp, 4, numa);
        } else {
            mp_limb_t inv21 = lmmp_inv_2_1_(numa[na - 1], numa[na - 2]);
            if (na < DIV_DIVIDE_THRESHOLD) {
                lmmp_div_basecase_(dst, xp, 2 * na, numa, na, inv21);
            } else {
                lmmp_div_divide_(dst, xp, 2 * na, numa, na, inv21);
            }
        }
        TEMP_FREE;
    }
}

引用了 DIV_DIVIDE_THRESHOLD, LIMB_B_2, LIMB_MAX, lmmp_div_2_s_(), lmmp_div_basecase_(), lmmp_div_divide_(), lmmp_inv_1_(), lmmp_inv_2_1_(), lmmp_not_(), lmmp_param_assert, TALLOC_TYPE, TEMP_DECL , 以及 TEMP_FREE.

被这些函数引用 lmmp_invappr_() , 以及 lmmp_invappr_newton_().

函数调用图:

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

lmmp_invappr_()

void lmmp_invappr_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na
	)

近似逆元计算 (invappr)

参数

dst	输出结果缓冲区，长度为na
numa	输入操作数，长度为na
na	输入操作数的 limb 长度

警告

na>0, MSB(numa)=1, sep(dst,numa)

返回

无返回值，结果存储在dst中，[dst,na] = (B^(2*na)-1)/[numa,na] - B^na + [0|-1]

在文件 inv.c 第 176 行定义.

                                                             {
    if (na < INV_NEWTON_THRESHOLD)
        lmmp_inv_basecase_(dst, numa, na);
    else
        lmmp_invappr_newton_(dst, numa, na);
}

引用了 INV_NEWTON_THRESHOLD, lmmp_inv_basecase_() , 以及 lmmp_invappr_newton_().

被这些函数引用 lmmp_inv_() , 以及 lmmp_inv_prediv_().

函数调用图:

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

lmmp_invappr_newton_()

void lmmp_invappr_newton_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na
	)

近似逆元计算（牛顿迭代法）

参数

dst	输出结果缓冲区，长度为na
numa	输入操作数，长度为na
na	输入操作数的 limb 长度

警告

na>4, MSB(numa)=1, sep(dst,numa)

返回

无返回值，结果存储在dst中，[dst,na]=(B^(2*na)-1)/[numa,na]-B^na+[0|-1]

在文件 inv.c 第 40 行定义.

                                                                    {
    lmmp_param_assert(na > 4);
    lmmp_param_assert(numa[na - 1] >= LIMB_B_2);
    
    mp_limb_t cy;
    mp_size_t nr = na, mn;
    mp_size_t sizes[LIMB_BITS], *sizp = sizes;
 
    do {
        *sizp = nr;
        nr = (nr >> 1) + 1;
        ++sizp;
    } while (nr >= INV_NEWTON_THRESHOLD);
 
    numa += na;
    dst += na;
 
    lmmp_inv_basecase_(dst - nr, numa - nr, nr);
 
    TEMP_DECL;
    mp_ptr xp = TALLOC_TYPE(3 * (na >> 1) + 3, mp_limb_t);
    do {
        na = *--sizp;
 
        // ar = 0:[numa-nr,nr]
        // an = 0:[numa-na,na]
        // ir = 1:[dst-nr,nr] = (B^(2*nr)-1)/ar - [0|1]
        // rem = ir*an-B^(na+nr)
        //-2*B^na < rem < 2*B^na
 
        //[xp] = rem
        if (na < INV_MODM_THRESHOLD || (mn = lmmp_fft_next_size_(na + 1)) >= na + nr) {
            lmmp_mul_(xp, numa - na, na, dst - nr, nr);
            lmmp_add_n_(xp + nr, xp + nr, numa - na, na + 1 - nr);
            cy = 1;  // for mod B^(na+1)
        } else {     // nr < na < mn < na+nr
 
            //[xp,mn] = [dst,nr] * [numa,na] mod (B^mn-1)
            lmmp_mul_mersenne_(xp, mn, numa - na, na, dst - nr, nr);
 
            //[xp,mn] += [numa,na]*B^nr mod (B^mn-1)
            cy = lmmp_add_n_(xp + nr, xp + nr, numa - na, mn - nr);
            cy = lmmp_add_nc_(xp, xp, numa - (na - (mn - nr)), na - (mn - nr), cy);
 
            //[xp,mn] -= B^(na+nr) mod (B^mn-1)
            xp[mn] = 1;
            lmmp_dec_1(xp + na + nr - mn, 1 - cy);
            lmmp_dec_1(xp, 1 - xp[mn]);
 
            cy = 0;  // for mod (B^mn-1)
        }
 
        // adjust ir,rem s.t.
        //  -B^na < rem = ir*an - B^(na+nr) < 0
        //  use this we can prove B^nr <= ir < 2*B^nr
        //  so inc/dec ir won't overflow
        if (xp[na] < 2) {  // rem>=0
 
            // rem-=cy*an s.t. rem[na]=0
            if (cy = xp[na]) {
                if (!lmmp_sub_n_(xp, xp, numa - na, na)) {
                    ++cy;
                    lmmp_sub_n_(xp, xp, numa - na, na);
                }
            }
 
            // rem-=cy*an s.t. 0<=rem<an
            if (lmmp_cmp_(xp, numa - na, na) >= 0) {
                lmmp_sub_n_(xp, xp, numa - na, na);
                ++cy;
            }
 
            // 0 < an-rem <= an < B^na , trunc to nr limbs
            lmmp_sub_nc_(xp + 2 * nr, numa - nr, xp + na - nr, nr, lmmp_cmp_(xp, numa - na, na - nr) > 0);
            ++cy;
 
            lmmp_dec_1(dst - nr, cy);
        } else {  // rem<0
            if (cy)
                lmmp_dec(xp);  // for neg to not
            // else (neg to not) compensate (mod transfer)
 
            if (xp[na] != LIMB_MAX) {
                lmmp_assert(xp[na] + lmmp_add_n_(xp, xp, numa - na, na) == LIMB_MAX);
                lmmp_inc(dst - nr);
            }
 
            //-rem
            lmmp_not_(xp + 2 * nr, xp + na - nr, nr);
        }
 
        // in = 1:[dst-na,na]
        // in = ir*B^(na-nr) + ir*(-rem/B^(na-nr))/B^(3*nr-na)
        // use inequality an*ir!=B^(na+nr),
        //(otherwise obviously contradictory),
        // we can prove
        //  an*in <= an*ir * ( 2*B^(na+nr) - an*ir ) * B^(-2*nr) < B^(2*na)
        // so in < B^(2*na)/an <= 2*B^(na),
        // inc below won't overflow
 
        // and via inequality -B^na < an*ir - B^(na+nr) < 0
        // we can prove in = (B^(2*na)-1)/an - [0|1]
        lmmp_mul_n_(xp, xp + 2 * nr, dst - nr, nr);
        cy = lmmp_add_n_(xp + nr, xp + nr, xp + 2 * nr, 2 * nr - na);
        if (lmmp_add_nc_(dst - na, xp + 3 * nr - na, xp + 4 * nr - na, na - nr, cy))
            lmmp_inc(dst - nr);
 
        nr = na;
    } while (sizp != sizes);
    TEMP_FREE;
}

引用了 INV_MODM_THRESHOLD, INV_NEWTON_THRESHOLD, LIMB_B_2, LIMB_BITS, LIMB_MAX, lmmp_add_n_(), lmmp_add_nc_(), lmmp_assert, lmmp_cmp_(), lmmp_dec, lmmp_dec_1, lmmp_fft_next_size_(), lmmp_inc, lmmp_inv_basecase_(), lmmp_mul_(), lmmp_mul_mersenne_(), lmmp_mul_n_(), lmmp_not_(), lmmp_param_assert, lmmp_sub_n_(), lmmp_sub_nc_(), TALLOC_TYPE, TEMP_DECL , 以及 TEMP_FREE.

被这些函数引用 lmmp_invappr_().

函数调用图:

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