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

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

浏览源代码.

函数
void	lmmp_bninv_ (mp_ptr restrict dstq, mp_srcptr restrict numa, mp_size_t na, mp_size_t ni)
static void	lmmp_bninv_appr_newton_ (mp_ptr restrict dstq, mp_srcptr restrict numa, mp_size_t na, mp_size_t ni)
	Copyright (C) 2026 HJimmyK(Jericho Knox)

函数说明

lmmp_bninv_()

void lmmp_bninv_	(	mp_ptr restrict	dstq,
		mp_srcptr restrict	numa,
		mp_size_t	na,
		mp_size_t	ni
	)

在文件 bninv.c 第 136 行定义.

                                                                                            {
    lmmp_param_assert(na > 0);
    lmmp_param_assert(dstq != NULL && numa != NULL);
    lmmp_param_assert(numa[na - 1] > 0);
    TEMP_DECL;
    if (na == 1) {
        mp_ptr restrict bnp = TALLOC_TYPE(3 + ni, mp_limb_t);
        lmmp_zero(bnp, 2 + ni);
        bnp[2 + ni] = 1;
        lmmp_div_1_(dstq, bnp, 3 + ni, numa[0]);
    } else if (na == 2) {
        mp_size_t bn = 2 * 2 + ni + 1;
        mp_ptr restrict bnp = TALLOC_TYPE(bn, mp_limb_t);
        lmmp_zero(bnp, bn - 1);
        bnp[bn - 1] = 1;
        mp_limb_t d[2] = {numa[0], numa[1]};
        lmmp_div_2_(dstq, bnp, bn, d);
    } else if (ni > na) {
        mp_ptr restrict B = TALLOC_TYPE(2 * na + ni + 1, mp_limb_t);
        lmmp_zero(B, 2 * na + ni);
        B[2 * na + ni] = 1;
        lmmp_div_(dstq, NULL, B, 2 * na + ni + 1, numa, na);
    } else {
        int shift = lmmp_leading_zeros_(numa[na - 1]);
        if (shift > 0) {
            mp_ptr restrict numa_shift = TALLOC_TYPE(na, mp_limb_t);
            lmmp_shl_(numa_shift, numa, na, shift);
            lmmp_bninv_appr_newton_(dstq, numa_shift, na, ni + 1);
            lmmp_shr_(dstq, dstq, na + ni + 2, LIMB_BITS - shift);
        } else {
            lmmp_bninv_appr_newton_(dstq, numa, na, ni + 1);
            lmmp_copy(dstq, dstq + 1, na + ni + 1);
            dstq[na + ni + 1] = 0;
        }
    }
    TEMP_FREE;
    return;
}

引用了 B, bn, LIMB_BITS, lmmp_bninv_appr_newton_(), lmmp_copy, lmmp_div_(), lmmp_div_1_(), lmmp_div_2_(), lmmp_leading_zeros_, lmmp_param_assert, lmmp_shl_(), lmmp_shr_(), lmmp_zero, n, TALLOC_TYPE, TEMP_DECL , 以及 TEMP_FREE.

函数调用图:

lmmp_bninv_appr_newton_()

static void lmmp_bninv_appr_newton_ ( mp_ptr restrict  dstq, mp_srcptr restrict  numa, mp_size_t  na, mp_size_t  ni  )

static

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

牛顿法求精确逆元（至多产生1的误差）

注解

dstq := B^(2*(na+ni)) // ([numa,na] * B^ni) + [0|1]

警告

eqsep(dstq,numa), dstq!=NULL, numa!=NULL, na>=3, MSB(numa)=1

在文件 bninv.c 第 85 行定义.

                                                                                                               {
    lmmp_param_assert(na >= 3);
    lmmp_param_assert(dstq != NULL && numa != NULL);
    lmmp_param_assert(numa[na - 1] > LIMB_B_2);
    TEMP_DECL;
    if (na < BNINV_NEWTON_THRESHOLD) {
        mp_ptr restrict bnp = TALLOC_TYPE(2 * na + ni + 1, mp_limb_t);
        lmmp_zero(bnp, 2 * na + ni + 1);
        bnp[2 * na + ni] = 1;
        mp_limb_t inv21 = lmmp_inv_2_1_(numa[na - 1], numa[na - 2]);
        lmmp_div_basecase_(dstq, bnp, 2 * na + ni + 1, numa, na, inv21);
    } else {
        mp_srcptr restrict a_hat;
 
        mp_size_t N = na + ni;
        mp_size_t a_hatn = N / 2 + 1;
        mp_size_t remn = N - a_hatn;
        mp_size_t ni_hat;
        if (remn > ni) {
            ni_hat = 0;
            a_hat = numa + na - a_hatn;
        } else {
            ni_hat = ni - remn;
            a_hat = numa;
            a_hatn = na;
        }
        mp_size_t q_hatn = a_hatn + ni_hat + 1;
        mp_size_t qrn = N + 1 - q_hatn;
        mp_ptr restrict q_hat = dstq + qrn;
 
        mp_ptr restrict q_hat_sqr = TALLOC_TYPE(2 * q_hatn, mp_limb_t);
        mp_ptr restrict q_hat_sqr_a = TALLOC_TYPE(2 * q_hatn + na, mp_limb_t);
 
        lmmp_bninv_appr_newton_(q_hat, a_hat, a_hatn, ni_hat);
        lmmp_sqr_(q_hat_sqr, q_hat, q_hatn);
        lmmp_mul_(q_hat_sqr_a, q_hat_sqr, 2 * q_hatn, numa, na);
        // we can assert q_hat_sqr_a[2*q_hatn+na-1] == 0
        lmmp_zero(dstq, qrn);
        if (2 * qrn + ni > 2 * N) {
            mp_size_t start = 2 * qrn + ni - 2 * N;
            lmmp_shl_(q_hat, q_hat, q_hatn, 1);  // assert no carry
            lmmp_sub_n_(dstq + start, dstq + start, q_hat_sqr_a, N + 1 - start);
        } else {
            mp_size_t start = 2 * N - 2 * qrn - ni;
            lmmp_shl_(q_hat, q_hat, q_hatn, 1);  // assert no carry
            lmmp_sub_n_(dstq, dstq, q_hat_sqr_a + start, N + 1);
        }
    }
    TEMP_FREE;
}

引用了 BNINV_NEWTON_THRESHOLD, LIMB_B_2, lmmp_bninv_appr_newton_(), lmmp_div_basecase_(), lmmp_inv_2_1_(), lmmp_mul_(), lmmp_param_assert, lmmp_shl_(), lmmp_sqr_, lmmp_sub_n_(), lmmp_zero, n, TALLOC_TYPE, TEMP_DECL , 以及 TEMP_FREE.

被这些函数引用 lmmp_bninv_() , 以及 lmmp_bninv_appr_newton_().

函数调用图:

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