LAMMP 4.2.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
函数 | 变量
cbrt_1.c 文件参考
#include "../../../include/lammp/impl/inlines.h"
#include "../../../include/lammp/numth.h"
+ cbrt_1.c 的引用(Include)关系图:

浏览源代码.

函数

ulong lmmp_cbrt_chebyshev_ (ulong n)
 计算算数立方根 floor(cbrt(n))
 
static double lmmp_cbrt_estimate (double a)
 
ulong lmmp_cbrt_ulong_ (ulong n)
 计算算数立方根 floor(cbrt(n))
 

变量

static const float coeff [16][3]
 
static const float factor_table [] = {1.000000, 1.259921, 1.587401}
 Copyright (C) 2026 HJimmyK(Jericho Knox)
 

函数说明

◆ lmmp_cbrt_chebyshev_()

ulong lmmp_cbrt_chebyshev_ ( ulong  n)

计算算数立方根 floor(cbrt(n))

参数
n被开方数
返回
floor(cbrt(n))
警告
n>0
注解
使用Chebyshev估计，当n较大时，此算法更占优势

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

40 {
41 lmmp_param_assert(n > 0);
42 typedef union {
43 ulong uword_val;
44 double double_val;
45 } uni;
46
47 int rem, mul;
48 double factor, root, dec, dec2;
49 ulong ret, expo, table_index;
50 uni alias;
51
52 /* upper_limit is the max cube root possible for one word */
53
54 const ulong upper_limit = 2642245; /* 2642245 < (2^64)^(1/3) */
55 const ulong expo_mask = 0x7FF0000000000000; /* exponent bits in double */
56 const ulong mantissa_mask = 0x000FFFFFFFFFFFFF; /* mantissa bits in float */
57 const ulong table_mask = 0x000F000000000000; /* first 4 bits of mantissa */
58 const uint mantissa_bits = 52;
59 const ulong bias_hex = 0x3FE0000000000000;
60 const uint bias = 1022;
61 alias.double_val = (double)n;
62
63 expo = alias.uword_val & expo_mask; /* extracting exponent */
64 expo >>= mantissa_bits;
65 expo -= bias; /* Subtracting bias */
66
67 /* extracting first 4 bits of mantissa, this will help select correct poly */
68 /* note mantissa of 0.5 is 0x0000000000000 not 0x1000000000000 */
69
70 table_index = alias.uword_val & table_mask;
71 table_index >>= (mantissa_bits - 4);
72
73 /* extracting decimal part, 0.5 <= dec <= 1 */
74 ret = alias.uword_val & mantissa_mask;
75 ret |= bias_hex;
76 alias.uword_val = ret;
77 dec = alias.double_val;
78
79 rem = expo % 3;
80 expo /= 3; /* cube root of 2^expo */
81 factor = factor_table[rem]; /* select factor */
82
83 /* Calculating cube root of dec using chebyshev approximation polynomial */
84 /* Evaluating approx polynomial at (dec) by Estrin's scheme */
85
86 dec2 = dec * dec;
87 root = (coeff[table_index][0] + coeff[table_index][1] * dec);
88 root += (coeff[table_index][2] * dec2);
89
90 mul = (ulong)1 << expo; /* mul = 2^expo */
91 root *= mul; /* dec^(1/3) * 2^(expo/3) */
92 root *= factor; /* root*= (expo%3)^(1/3) */
93 ret = root;
94
95 /* In case ret^3 or (ret+1)^3 will cause overflow */
96
97 if (ret >= upper_limit) {
98 if (n >= upper_limit * upper_limit * upper_limit)
99 return upper_limit;
100 ret = upper_limit - 1;
101 }
102 while (ret * ret * ret <= n) {
103 (ret) += 1;
104 if (ret == upper_limit)
105 break;
106 }
107 while (ret * ret * ret > n) (ret) -= 1;
108
109 return ret;
110}
coeff
static const float coeff[16][3]
Definition cbrt_1.c:23
factor_table
static const float factor_table[]
Copyright (C) 2026 HJimmyK(Jericho Knox)
Definition cbrt_1.c:20
lmmp_param_assert
#define lmmp_param_assert(x)
Definition lmmp.h:401
n
#define n
uint
uint32_t uint
Definition numth.h:31
ulong
uint64_t ulong
Definition numth.h:32

引用了 coeff, factor_table, lmmp_param_assert , 以及 n.

被这些函数引用 lmmp_cbrt_ulong_().

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

◆ lmmp_cbrt_estimate()

static double lmmp_cbrt_estimate ( double  a)
inlinestatic

在文件 cbrt_1.c 第 112 行定义.

112 {
113 typedef union {
114 slong uword_val;
115 double double_val;
116 } uni;
117
118 uni alias;
119 ulong n;
120
121 const ulong mul_factor = (ulong)6148914691236517205;
122 slong s = (slong)4607182418800017408; /* ((1 << 10) - 1) << 52 */
123
124 alias.double_val = a;
125 n = alias.uword_val;
126 n -= s;
127 n = lmmp_mulh_(n, mul_factor);
128 n += s;
129 alias.uword_val = n;
130 return alias.double_val;
131}
lmmp_mulh_
#define lmmp_mulh_
Definition inlines.h:165
s
#define s
slong
int64_t slong
Definition numth.h:34

引用了 lmmp_mulh_, n , 以及 s.

被这些函数引用 lmmp_cbrt_ulong_().

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

◆ lmmp_cbrt_ulong_()

ulong lmmp_cbrt_ulong_ ( ulong  n)

计算算数立方根 floor(cbrt(n))

参数
n被开方数
返回
floor(cbrt(n))
注解
会依据n的大小，选择合适的算法

在文件 cbrt_1.c 第 133 行定义.

133 {
134 double val, x, xcub, num, den;
135 ulong ret, upper_limit;
136
137 if (n < 125)
138 return (n >= 1) + (n >= 8) + (n >= 27) + (n >= 64);
139 if (n < 1331)
140 return 5 + (n >= 216) + (n >= 343) + (n >= 512) + (n >= 729) + (n >= 1000);
141 if (n < 4913)
142 return 11 + (n >= 1728) + (n >= 2197) + (n >= 2744) + (n >= 3375) + (n >= 4096);
143
144 val = (double)n;
145
146 if (n >= 1ULL << 46)
147 return lmmp_cbrt_chebyshev_(n);
148
149 upper_limit = 2642245; /* 2642245 = floor((2^64)^(1/3)) */
150
151 x = lmmp_cbrt_estimate((double)n);
152
153 /* Kahan's iterations to get cube root */
154
155 xcub = x * x * x;
156 num = (xcub - val) * x;
157 den = (xcub + xcub + val);
158 x -= (num / den);
159 ret = x;
160
161 if (ret >= upper_limit) {
162 if (n >= upper_limit * upper_limit * upper_limit)
163 return upper_limit;
164 ret = upper_limit - 1;
165 }
166 while (ret * ret * ret <= n) {
167 (ret) += 1;
168 if (ret == upper_limit)
169 break;
170 }
171 while (ret * ret * ret > n) (ret) -= 1;
172
173 return ret;
174}
lmmp_cbrt_estimate
static double lmmp_cbrt_estimate(double a)
Definition cbrt_1.c:112
lmmp_cbrt_chebyshev_
ulong lmmp_cbrt_chebyshev_(ulong n)
计算算数立方根 floor(cbrt(n))
Definition cbrt_1.c:40

引用了 lmmp_cbrt_chebyshev_(), lmmp_cbrt_estimate() , 以及 n.

被这些函数引用 lmmp_nthroot_ulong_().

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

变量说明

◆ coeff

const float coeff[16][3]
static
初始值:
= {{0.445434042, 0.864136635, -0.335205926},
{0.454263239, 0.830878907, -0.303884962},
{0.462761624, 0.800647514, -0.276997626},
{0.470958569, 0.773024522, -0.253724515},
{0.478879482, 0.747667468, -0.233429710},
{0.486546506, 0.724292830, -0.215613166},
{0.493979069, 0.702663686, -0.199877008},
{0.501194325, 0.682580388, -0.185901247},
{0.508207500, 0.663873398, -0.173426009},
{0.515032183, 0.646397742, -0.162238357},
{0.521680556, 0.630028647, -0.152162376},
{0.528163588, 0.614658092, -0.143051642},
{0.534491194, 0.600192044, -0.134783425},
{0.540672371, 0.586548233, -0.127254189},
{0.546715310, 0.573654340, -0.120376066},
{0.552627494, 0.561446514, -0.114074068}}

在文件 cbrt_1.c 第 23 行定义.

23 {{0.445434042, 0.864136635, -0.335205926}, /* [0.50000, 0.53125] */
24 {0.454263239, 0.830878907, -0.303884962}, /* [0.53125, 0.56250] */
25 {0.462761624, 0.800647514, -0.276997626}, /* [0.56250, 0.59375] */
26 {0.470958569, 0.773024522, -0.253724515}, /* [0.59375, 0.62500] */
27 {0.478879482, 0.747667468, -0.233429710}, /* [0.62500, 0.65625] */
28 {0.486546506, 0.724292830, -0.215613166}, /* [0.65625, 0.68750] */
29 {0.493979069, 0.702663686, -0.199877008}, /* [0.68750, 0.71875] */
30 {0.501194325, 0.682580388, -0.185901247}, /* [0.71875, 0.75000] */
31 {0.508207500, 0.663873398, -0.173426009}, /* [0.75000, 0.78125] */
32 {0.515032183, 0.646397742, -0.162238357}, /* [0.78125, 0.81250] */
33 {0.521680556, 0.630028647, -0.152162376}, /* [0.81250, 0.84375] */
34 {0.528163588, 0.614658092, -0.143051642}, /* [0.84375, 0.87500] */
35 {0.534491194, 0.600192044, -0.134783425}, /* [0.87500, 0.90625] */
36 {0.540672371, 0.586548233, -0.127254189}, /* [0.90625, 0.93750] */
37 {0.546715310, 0.573654340, -0.120376066}, /* [0.93750, 0.96875] */
38 {0.552627494, 0.561446514, -0.114074068}}; /* [0.96875, 1.00000] */

被这些函数引用 exp2_fixed_128(), exp2_fixed_64(), lmmp_cbrt_chebyshev_(), log2_fixed_128() , 以及 log2_fixed_64().

◆ factor_table

const float factor_table[] = {1.000000, 1.259921, 1.587401}
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/.

在文件 cbrt_1.c 第 20 行定义.

20{1.000000, 1.259921, 1.587401}; /* {1^(1/3), 2^(1/3), 4^(1/3)} */

被这些函数引用 lmmp_cbrt_chebyshev_().