cbrt.c 文件参考

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

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

浏览源代码.

宏定义
#define	Ahr   (dst + nlo)
#define	MMIN   1728

函数
mp_limb_t	lmmp_cbrt_3_ (mp_limb_t a0, mp_limb_t a1, mp_limb_t a2)
	计算算数立方根 floor(cbrt(a0+a1*B+a2*B^2))
void	lmmp_cbrt_6_ (mp_ptr dst, mp_srcptr numa, mp_size_t na)
	计算算数立方根 floor(cbrt([numa,na]))
mp_size_t	lmmp_cbrtapprox_ (mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_size_t ni)
mp_limb_t	lmmp_cbrtapprox_3_ (mp_limb_t a0, mp_limb_t a1, mp_limb_t a2)
	计算近似立方根 floor(cbrt(a0+a1*B+a2*B^2))-[0|1]
void	lmmp_cbrtapprox_6_ (mp_ptr dst, mp_srcptr numa, mp_size_t na)
	计算近似立方根 floor(cbrt([numa,na]))-[0|1]
static void	lmmp_cube_3_ (mp_ptr restrict dst, mp_limb_t a)
	Copyright (C) 2026 HJimmyK(Jericho Knox)
static mp_size_t	lmmp_cube_6_ (mp_ptr restrict dst, mp_srcptr restrict numa)

宏定义说明

Ahr

#define Ahr   (dst + nlo)

MMIN

#define MMIN   1728

函数说明

lmmp_cbrt_3_()

mp_limb_t lmmp_cbrt_3_	(	mp_limb_t	a0,
		mp_limb_t	a1,
		mp_limb_t	a2
	)

计算算数立方根 floor(cbrt(a0+a1*B+a2*B^2))

参数

a0	低位 limb
a1	中位 limb
a2	高位 limb

警告

a1>0

注解

a2可以为0，但a1需要大于0，即这个数至少应有65个bit

返回

floor(cbrt(a0+a1*B+a2*B^2))

在文件 cbrt.c 第 84 行定义.

                                                                 {
    lmmp_param_assert(a1 > 0);
 
    mp_limb_t r = lmmp_cbrtapprox_3_(a0, a1, a2);
    if (r == LIMB_MAX)
        return LIMB_MAX;
    mp_limb_t t[3], a[3] = {a0, a1, a2};
    lmmp_cube_3_(t, r + 1);
    int cmp = lmmp_cmp_(t, a, 3);
    // approx的结果至多只会低估1
    if (cmp <= 0)
        return r + 1;
    else
        return r;
}

引用了 a0, a1, a2, LIMB_MAX, lmmp_cbrtapprox_3_(), lmmp_cmp_(), lmmp_cube_3_(), lmmp_param_assert, n , 以及 t.

函数调用图:

lmmp_cbrt_6_()

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

计算算数立方根 floor(cbrt([numa,na]))

参数

dst	结果指针（长度为 2 个limb）
numa	被开方数指针
na	被开方数的 limb 长度

警告

dst!=NULL, numa!=NULL, 3<na<=6, numa[na-1]!=0, eqsep(dst,numa)

在文件 cbrt.c 第 146 行定义.

                                                            {
    mp_limb_t ret[2];
    lmmp_cbrtapprox_6_(ret, numa, na);
 
    if (ret[1] == LIMB_MAX && ret[0] == LIMB_MAX) {
        dst[0] = LIMB_MAX;
        dst[1] = LIMB_MAX;
    } else {
        mp_limb_t r[2];
        r[0] = ret[0] + 1;
        r[1] = ret[1] + (r[0] == 0 ? 1 : 0);
        mp_limb_t t[6];
        mp_size_t tn = lmmp_cube_6_(t, r);
        if (tn > na) {
            dst[0] = ret[0];
            dst[1] = ret[1];
        } else {
            lmmp_debug_assert(tn == na);
            int cmp = lmmp_cmp_(t, numa, na);
            // approx的结果至多只会低估1
            if (cmp <= 0) {
                dst[0] = r[0];
                dst[1] = r[1];
            } else {
                dst[0] = ret[0];
                dst[1] = ret[1];
            }
        }
    }
}

引用了 LIMB_MAX, lmmp_cbrtapprox_6_(), lmmp_cmp_(), lmmp_cube_6_(), lmmp_debug_assert, n , 以及 t.

函数调用图:

lmmp_cbrtapprox_()

mp_size_t lmmp_cbrtapprox_	(	mp_ptr restrict	dst,
		mp_srcptr restrict	numa,
		mp_size_t	na,
		mp_size_t	ni
	)

假设计算 n 次根式：A^(1/n) 设 beta = B^b，alpha = Ah^(1/n)，rho = A^(1/n)。 已知 tk 满足 floor(Ah^(1/n)) - 1 <= tk <= floor(Ah^(1/n))， 因此 xk = tk * beta 与真实根 rho 的误差满足（最坏情况）： e = rho - xk <= 2 * beta。

牛顿迭代实数形式为： F(x) = ((n-1)*x + A / x^(n-1)) / n。 令 x = rho - e，在 rho 处泰勒展开可得： F(x) - rho = (n(n-1)/2) * (e^2 / rho) + O(e^3 / rho^2)。

由于 rho >= alpha * beta（因为 A >= Ah * beta^n，开方后略大于 alpha*beta）， 将 e <= 2*beta 代入，忽略高阶小量，得到误差上限： F(x) - rho <= (n(n-1)/2) * ((2*beta)^2 / (alpha*beta)) = (2*n*(n-1)*beta) / alpha。

为使得最终取整后的绝对误差控制在 1 以内，只需令上式 <= 1： (2*n*(n-1)*beta) / alpha <= 1 => alpha >= 2*n*(n-1)*beta。

两边同时取 n 次方： Ah >= [2*n*(n-1)]^n * beta^n = [2*n*(n-1)]^n * B^(n*b)。

因此，若要保证递归校正始终收敛且误差不超过 1，常数 m 至少取： m_min = [2*n*(n-1)]^n。

对于立方根，即有：m_min = 1728

注：1728 是一个保守的上界，不满足此上界并不意味着迭代一定不收敛。

在文件 cbrt.c 第 177 行定义.

                                                                                                     {
    lmmp_param_assert(na > 0);
    lmmp_param_assert(dst != NULL && numa != NULL);
    mp_size_t n = ni + na;
    if (n == 1) {
        mp_limb_t a_cbrt = lmmp_cbrt_ulong_(numa[0]);
        dst[0] = a_cbrt;
        return 1;
    } else if (n <= 3) {
        if (ni == 0)
            dst[0] = lmmp_cbrtapprox_3_(numa[0], numa[1], numa[2]);
        else if (ni == 1)
            dst[0] = lmmp_cbrtapprox_3_(0, numa[0], numa[1]);
        else
            dst[0] = lmmp_cbrtapprox_3_(0, 0, numa[0]);
        return 1;
    } else if (n <= 6) {
        mp_limb_t a[6];
        mp_size_t i;
        for (i = 0; i < ni; i++) {
            a[i] = 0;
        }
        for (mp_size_t j = 0; i < n; i++, j++) {
            a[i] = numa[j];
        }
        lmmp_cbrtapprox_6_(dst, a, n);
        return 2;
    } else {
        /*
                A     = Ah * B^(3*lo) + Al
 
                Ahr   = floor(Ah^(1/3))
                x_k   = Ahr * B^lo
 
                x_k+1 = (2*x_k + A / x_k^2 ) / 3
        */
        mp_size_t nhi = (n - 1) / 2;
        mp_size_t nlo = nhi - (nhi % 3);
        nhi = n - nlo;
        /**
         * 假设计算 n 次根式：A^(1/n)
         * 设 beta = B^b，alpha = Ah^(1/n)，rho = A^(1/n)。
         * 已知 tk 满足 floor(Ah^(1/n)) - 1 <= tk <= floor(Ah^(1/n))，
         * 因此 xk = tk * beta 与真实根 rho 的误差满足（最坏情况）：
         *     e = rho - xk <= 2 * beta。
         *
         * 牛顿迭代实数形式为：
         *     F(x) = ((n-1)*x + A / x^(n-1)) / n。
         * 令 x = rho - e，在 rho 处泰勒展开可得：
         *     F(x) - rho = (n(n-1)/2) * (e^2 / rho) + O(e^3 / rho^2)。
         *
         * 由于 rho >= alpha * beta（因为 A >= Ah * beta^n，开方后略大于 alpha*beta），
         * 将 e <= 2*beta 代入，忽略高阶小量，得到误差上限：
         *     F(x) - rho <= (n(n-1)/2) * ((2*beta)^2 / (alpha*beta))
         *                 = (2*n*(n-1)*beta) / alpha。
         *
         * 为使得最终取整后的绝对误差控制在 1 以内，只需令上式 <= 1：
         *     (2*n*(n-1)*beta) / alpha <= 1
         *  =>  alpha >= 2*n*(n-1)*beta。
         *
         * 两边同时取 n 次方：
         *     Ah >= [2*n*(n-1)]^n * beta^n
         *        = [2*n*(n-1)]^n * B^(n*b)。
         *
         * 因此，若要保证递归校正始终收敛且误差不超过 1，常数 m 至少取：
         *     m_min = [2*n*(n-1)]^n。
         *
         * 对于立方根，即有：m_min = 1728
         * 
         * 注：1728 是一个保守的上界，不满足此上界并不意味着迭代一定不收敛。
         */
#define MMIN 1728
#if LAMMP_DEBUG_ASSERT_CHECK == 1
        if (nhi == nlo + 1) {
            lmmp_debug_assert(numa[na - 1] > MMIN);
        }
#endif
        nlo /= 3;
        lmmp_debug_assert(nlo > 0);
 
        TEMP_DECL;
        mp_size_t rn, Adivn = nhi + nlo;
        mp_ptr restrict Adiv;
 
#define Ahr (dst + nlo)
        if (ni >= 3 * nlo) {
            //  __________________ n ___________________
            //  |__________ na __________|_____ ni ____|
            //  |xxxxxxxxxxxxxxxxxxxxxxxx|0000000000000|
            //  |___________ nhi __________|__ 3*nlo __|
 
            rn = lmmp_cbrtapprox_(Ahr, numa, na, ni - 3 * nlo);
 
            Adiv = TALLOC_TYPE(Adivn, mp_limb_t);
            lmmp_zero(Adiv, ni - 2 * nlo);
            lmmp_copy(Adiv + ni - 2 * nlo, numa, na);
 
        } else if (ni >= 2 * nlo) {
            //  __________________ n ___________________
            //  |__________ na __________|_____ ni ____|
            //  |xxxxxxxxxxxxxxxxxxxxxxxx|0000000000000|
            //  |_______ nhi _____|__ nlo _|__ 2*nlo __|
 
            rn = lmmp_cbrtapprox_(Ahr, numa + na - nhi, nhi, 0);
 
            Adiv = TALLOC_TYPE(Adivn, mp_limb_t);
            lmmp_zero(Adiv, ni - 2 * nlo);
            lmmp_copy(Adiv + ni - 2 * nlo, numa, na);
 
        } else {
            //  __________________ n ___________________
            //  |__________ na __________|_____ ni ____|
            //  |xxxxxxxxxxxxxxxxxxxxxxxx|0000000000000|
            //  |_____ nhi ___|__ nlo__|____ 2*nlo ____|
 
            rn = lmmp_cbrtapprox_(Ahr, numa + na - nhi, nhi, 0);
 
            Adiv = TALLOC_TYPE(Adivn, mp_limb_t);
            lmmp_copy(Adiv, numa + na - Adivn, Adivn);
        }
 
        mp_size_t Ahr2n = rn * 2;
        mp_ptr restrict Ahr2 = TALLOC_TYPE(Ahr2n, mp_limb_t);
        lmmp_sqr_(Ahr2, Ahr, rn);
        Ahr2n -= (Ahr2[Ahr2n - 1] == 0) ? 1 : 0;
 
        lmmp_debug_assert(Adivn >= Ahr2n);
        mp_size_t qn = Adivn - Ahr2n + 1;
        mp_size_t rkdivn = (n + 2) / 3 + 2; // 额外多两个limb，因为需要作为加法缓冲区
        mp_ptr restrict rkdiv = TALLOC_TYPE(rkdivn, mp_limb_t);
 
        lmmp_debug_assert(qn <= rkdivn);
        lmmp_zero(rkdiv + qn, rkdivn - qn); //高位清零
 
        lmmp_div_(rkdiv, NULL, Adiv, Adivn, Ahr2, Ahr2n);
 
        mp_limb_t cy = lmmp_shl_(Ahr, Ahr, rn, 1);
        Ahr[rn] = cy;
        rn += cy > 0 ? 1 : 0;
 
        lmmp_debug_assert(rn + nlo + 1 <= rkdivn);
        cy = lmmp_add_n_(rkdiv + nlo, Ahr, rkdiv + nlo, rn);
        rn += nlo;
        rkdiv[rn] = cy;
        rn += cy > 0 ? 1 : 0;
 
        lmmp_div_1_(dst, rkdiv, rn, 3);
 
        while (dst[rn - 1] == 0) --rn;
 
        TEMP_FREE;
        return rn;
    }
#undef Ahr
}

引用了 Ahr, lmmp_add_n_(), lmmp_cbrt_ulong_(), lmmp_cbrtapprox_(), lmmp_cbrtapprox_3_(), lmmp_cbrtapprox_6_(), lmmp_copy, lmmp_debug_assert, lmmp_div_(), lmmp_div_1_(), lmmp_param_assert, lmmp_shl_(), lmmp_sqr_, lmmp_zero, MMIN, n, TALLOC_TYPE, TEMP_DECL , 以及 TEMP_FREE.

被这些函数引用 lmmp_cbrtapprox_().

函数调用图:

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

lmmp_cbrtapprox_3_()

mp_limb_t lmmp_cbrtapprox_3_	(	mp_limb_t	a0,
		mp_limb_t	a1,
		mp_limb_t	a2
	)

计算近似立方根 floor(cbrt(a0+a1*B+a2*B^2))-[0|1]

参数

a0	低位 limb
a1	中位 limb
a2	高位 limb

警告

a1>0

注解

a2可以为0，但a1需要大于0，即这个数至少应有65个bit

返回

floor(cbrt(a0+a1*B+a2*B^2))-[0|1]

在文件 cbrt.c 第 42 行定义.

                                                                       {
    lmmp_param_assert(a1 > 0);
    mp_limb_t x[2];
    /* exact high 65 bits */
    mp_limb_t a_hi;
    mp_bitcnt_t bits;
    if (a2 == 0) {
        mp_bitcnt_t a1_bits = lmmp_limb_bits_(a1);
        bits = LIMB_BITS + a1_bits;
        a1_bits--;
        if (a1_bits == 0)
            a_hi = a0;
        else
            a_hi = (a1 << (LIMB_BITS - a1_bits)) | (a0 >> a1_bits);
    } else {
        mp_bitcnt_t a2_bits = lmmp_limb_bits_(a2);
        bits = LIMB_BITS * 2 + a2_bits;
        a2_bits--;
        if (a2_bits == 0)
            a_hi = a1;
        else
            a_hi = (a2 << (LIMB_BITS - a2_bits)) | (a1 >> a2_bits);
    }
    lmmp_debug_assert(bits >= 65);
 
    x[1] = bits - 1;
    x[0] = log2_fixed_64(a_hi);
 
    mp_limb_t r = lmmp_div_1_(x, x, 2, 3);
    if (2 * r >= 3) // round
        lmmp_inc(x);
 
    mp_bitcnt_t shift = x[1];
    x[0] = exp2_fixed_64(x[0]);
 
    lmmp_debug_assert(shift <= 64);
    if (shift == 64)
        return LIMB_MAX;
    else
        return (x[0] >> (64 - shift)) | (1ULL << shift);
}

引用了 a0, a1, a2, exp2_fixed_64(), LIMB_BITS, LIMB_MAX, lmmp_debug_assert, lmmp_div_1_(), lmmp_inc, lmmp_limb_bits_, lmmp_param_assert, log2_fixed_64() , 以及 n.

被这些函数引用 lmmp_cbrt_3_() , 以及 lmmp_cbrtapprox_().

函数调用图:

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

lmmp_cbrtapprox_6_()

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

计算近似立方根 floor(cbrt([numa,na]))-[0|1]

参数

dst	结果指针（长度为 2 个limb）
numa	被开方数指针
na	被开方数的 limb 长度

警告

dst!=NULL, numa!=NULL, 3<na<=6, numa[na-1]!=0, eqsep(dst,numa)

在文件 cbrt.c 第 100 行定义.

                                                                  {
    lmmp_param_assert(na > 3 && na <= 6);
    lmmp_param_assert(dst != NULL && numa != NULL);
    lmmp_param_assert(numa[na - 1] != 0);
    /* extract the first 129 bits */
    int bits = lmmp_limb_bits_(numa[na - 1]);
    mp_bitcnt_t n = bits - 1;
    mp_limb_t high, low;
    if (bits == 1) {
        high = numa[na - 2];
        low = numa[na - 3];
    } else {
        bits--;
        high = (numa[na - 1] << (64 - bits)) | (numa[na - 2] >> bits);
        low = (numa[na - 2] << (64 - bits)) | (numa[na - 3] >> bits);
    }
 
    n += LIMB_BITS * (na - 1);
    mp_limb_t x[3] = {0, 0, n};
 
    log2_fixed_128(x, high, low);
    mp_limb_t r = lmmp_div_1_(x, x, 3, 3);
    if (2 * r >= 3) // round
       lmmp_inc(x);
 
    n = x[2];
    high = x[1];
    low = x[0];
 
    exp2_fixed_128(x, high, low);
 
    lmmp_debug_assert(n >= 64 && n <= 128);
    if (n == 64) {
        dst[0] = x[1];
        dst[1] = 1;
    } else if (n < 128) {
        n -= 64;
        mp_limb_t t = 1ULL << n;
        dst[1] = (x[1] >> (64 - n)) | t;
        dst[0] = (x[1] << n) | (x[0] >> (64 - n));
    } else {
        dst[1] = LIMB_MAX;
        dst[0] = LIMB_MAX;
    }
}

引用了 exp2_fixed_128(), LIMB_BITS, LIMB_MAX, lmmp_debug_assert, lmmp_div_1_(), lmmp_inc, lmmp_limb_bits_, lmmp_param_assert, log2_fixed_128(), n , 以及 t.

被这些函数引用 lmmp_cbrt_6_() , 以及 lmmp_cbrtapprox_().

函数调用图:

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

lmmp_cube_3_()

static void lmmp_cube_3_ ( mp_ptr restrict  dst, mp_limb_t  a  )

inlinestatic

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

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

                                                                  {
    mp_limb_t t[2];
    lmmp_mullh_(a, a, t);
    lmmp_mullh_(t[0], a, dst);
    lmmp_mullh_(t[1], a, t);
    dst[1] += t[0];
    dst[2] = t[1] + (dst[1] < t[0] ? 1 : 0);
}

引用了 lmmp_mullh_, n , 以及 t.

被这些函数引用 lmmp_cbrt_3_().

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

lmmp_cube_6_()

static mp_size_t lmmp_cube_6_ ( mp_ptr restrict  dst, mp_srcptr restrict  numa  )

inlinestatic

在文件 cbrt.c 第 33 行定义.

                                                                                   {
    mp_limb_t t[4];
    lmmp_sqr_basecase_(t, numa, 2);
    lmmp_mul_basecase_(dst, t, 4, numa, 2);
    mp_size_t n = 6;
    while (dst[n - 1] == 0) --n;
    return n;
}

引用了 lmmp_mul_basecase_(), lmmp_sqr_basecase_(), n , 以及 t.

被这些函数引用 lmmp_cbrt_6_().

函数调用图:

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