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

浏览源代码.

函数

mp_size_t lmmp_pow_ (mp_ptr restrict dst, mp_size_t rn, mp_srcptr restrict base, mp_size_t n, ulong exp)
 
mp_size_t lmmp_pow_1_size_ (mp_limb_t base, ulong exp)
 Copyright (C) 2026 HJimmyK(Jericho Knox)
 
mp_size_t lmmp_pow_size_ (mp_srcptr base, mp_size_t n, ulong exp)
 计算幂次方需要的limb缓冲区长度 [base,n] ^ exp
 

函数说明

◆ lmmp_pow_()

mp_size_t lmmp_pow_ ( mp_ptr restrict  dst,
mp_size_t  rn,
mp_srcptr restrict  base,
mp_size_t  n,
ulong  exp 
)

在文件 pow.c 第 132 行定义.

132 {
133 lmmp_param_assert(n > 0);
134 lmmp_param_assert(exp > 0);
135 lmmp_param_assert(base[n - 1] != 0);
136 if (exp == 1) {
137 lmmp_copy(dst, base, n);
138 return n;
139 } else if (exp == 2) {
140 lmmp_sqr_(dst, base, n);
141 rn = n << 1;
142 rn -= (dst[rn - 1] == 0);
143 return rn;
144 } else {
145 mp_size_t base_tz = 0;
146 while (*base == 0) {
147 ++base_tz;
148 ++base;
149 --n;
150 }
151 base_tz *= exp;
152 lmmp_zero(dst, base_tz);
153 dst += base_tz;
154 if (n == 1) {
155 if (exp <= POW_1_EXP_THRESHOLD) {
156 dst[0] = base[0];
157 rn = 1;
158 for (mp_size_t i = 1; i < exp; ++i) {
159 dst[rn] = lmmp_mul_1_(dst, dst, rn, base[0]);
160 ++rn;
161 rn -= (dst[rn - 1] == 0);
162 }
163 return rn + base_tz;
164 } else {
165 return lmmp_pow_1_(dst, rn, base[0], exp) + base_tz;
166 }
167 } else { /* n > 2 */
168 if (exp > POW_WIN2_EXP_THRESHOLD && n > POW_WIN2_N_THRESHOLD) {
169 if ((exp % 4 == 3) || (2 * lmmp_limb_popcnt_(exp) >= (lmmp_limb_bits_(exp)))) {
170 return lmmp_pow_win2_(dst, rn, base, n, exp) + base_tz;
171 }
172 }
173 if (exp & 1) {
174 return lmmp_pow_basecase_(dst, rn, base, n, exp) + base_tz;
175 }
176
177 int tz = lmmp_tailing_zeros_(exp);
178 TEMP_DECL;
179 mp_ptr restrict sq = TALLOC_TYPE((rn + 2) >> 1, mp_limb_t);
180 exp >>= tz;
181
182 if (tz & 1) {
183 if (exp == 1) {
184 lmmp_copy(sq, base, n);
185 rn = n;
186 } else {
187 mp_size_t rn1 = lmmp_pow_size_(base, n, exp);
188 rn = lmmp_pow_basecase_(sq, rn1, base, n, exp);
189 }
190 int i = 2;
191 for (; i <= tz; i += 2) {
192 lmmp_sqr_(dst, sq, rn);
193 rn <<= 1;
194 rn -= (dst[rn - 1] == 0);
195 lmmp_sqr_(sq, dst, rn);
196 rn <<= 1;
197 rn -= (sq[rn - 1] == 0);
198 }
199 lmmp_sqr_(dst, sq, rn);
200 rn <<= 1;
201 rn -= (dst[rn - 1] == 0);
202 } else {
203 if (exp == 1) {
204 lmmp_copy(dst, base, n);
205 rn = n;
206 } else {
207 mp_size_t rn1 = lmmp_pow_size_(base, n, exp);
208 rn = lmmp_pow_basecase_(dst, rn1, base, n, exp);
209 }
210 int i = 2;
211 for (; i <= tz; i += 2) {
212 lmmp_sqr_(sq, dst, rn);
213 rn <<= 1;
214 rn -= (sq[rn - 1] == 0);
215 lmmp_sqr_(dst, sq, rn);
216 rn <<= 1;
217 rn -= (dst[rn - 1] == 0);
218 }
219 }
220 TEMP_FREE;
221 return rn + base_tz;
222 }
223 }
224}
lmmp_limb_bits_
#define lmmp_limb_bits_
Definition inlines.h:162
lmmp_sqr_
#define lmmp_sqr_
Definition inlines.h:166
lmmp_limb_popcnt_
#define lmmp_limb_popcnt_
Definition inlines.h:163
lmmp_tailing_zeros_
#define lmmp_tailing_zeros_
Definition inlines.h:161
mp_ptr
mp_limb_t * mp_ptr
Definition lmmp.h:80
lmmp_copy
#define lmmp_copy(dst, src, n)
Definition lmmp.h:367
lmmp_zero
#define lmmp_zero(dst, n)
Definition lmmp.h:369
mp_size_t
uint64_t mp_size_t
Definition lmmp.h:77
mp_limb_t
uint64_t mp_limb_t
Definition lmmp.h:76
lmmp_param_assert
#define lmmp_param_assert(x)
Definition lmmp.h:401
lmmp_mul_1_
mp_limb_t lmmp_mul_1_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_limb_t x)
大数乘以单limb操作 [dst,na] = [numa,na] * x
POW_WIN2_EXP_THRESHOLD
#define POW_WIN2_EXP_THRESHOLD
Definition mparam.h:101
POW_1_EXP_THRESHOLD
#define POW_1_EXP_THRESHOLD
Definition mparam.h:98
POW_WIN2_N_THRESHOLD
#define POW_WIN2_N_THRESHOLD
Definition mparam.h:104
n
#define n
lmmp_pow_win2_
mp_size_t lmmp_pow_win2_(mp_ptr dst, mp_size_t rn, mp_srcptr base, mp_size_t n, ulong exp)
计算幂次方2比特窗口快速幂算法 [dst,rn] = [base,n] ^ exp
lmmp_pow_1_
mp_size_t lmmp_pow_1_(mp_ptr dst, mp_size_t rn, mp_limb_t base, ulong exp)
计算幂次方 [dst,rn] = [base,1] ^ exp
lmmp_pow_basecase_
mp_size_t lmmp_pow_basecase_(mp_ptr dst, mp_size_t rn, mp_srcptr base, mp_size_t n, ulong exp)
计算奇数次幂算法 [dst,rn] = [base,n] ^ exp
lmmp_pow_size_
mp_size_t lmmp_pow_size_(mp_srcptr base, mp_size_t n, ulong exp)
计算幂次方需要的limb缓冲区长度 [base,n] ^ exp
Definition pow.c:81
TEMP_DECL
#define TEMP_DECL
Definition tmp_alloc.h:131
TEMP_FREE
#define TEMP_FREE
Definition tmp_alloc.h:150
TALLOC_TYPE
#define TALLOC_TYPE(n, type)
Definition tmp_alloc.h:148

引用了 lmmp_copy, lmmp_limb_bits_, lmmp_limb_popcnt_, lmmp_mul_1_(), lmmp_param_assert, lmmp_pow_1_(), lmmp_pow_basecase_(), lmmp_pow_size_(), lmmp_pow_win2_(), lmmp_sqr_, lmmp_tailing_zeros_, lmmp_zero, n, POW_1_EXP_THRESHOLD, POW_WIN2_EXP_THRESHOLD, POW_WIN2_N_THRESHOLD, TALLOC_TYPE, TEMP_DECL , 以及 TEMP_FREE.

+ 函数调用图:

◆ lmmp_pow_1_size_()

mp_size_t lmmp_pow_1_size_ ( mp_limb_t  base,
ulong  exp 
)

Copyright (C) 2026 HJimmyK(Jericho Knox)

计算幂次方需要的limb缓冲区长度 base ^ exp

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

在文件 pow.c 第 25 行定义.

25 {
26 lmmp_param_assert(base >= 1);
27 lmmp_param_assert(exp > 0);
28 if (base == 1) {
29 return 1;
30 } else if (exp <= 2) {
31 return 3;
32 } else if (exp <= MP_UINT_MAX) {
33 /*
34 base = b * 2^base_tz
35 */
36 slong base_tz = lmmp_limb_bits_(base);
37 uint32_t b;
38 if (base_tz < 32) {
39 b = base << (32 - base_tz);
40 } else {
41 b = base >> (base_tz - 32);
42 }
43 base_tz = base_tz - 32;
44 base_tz *= exp;
45 mp_size_t rn = xlog2n_ceil(exp, b);
46 rn += base_tz;
47 rn = (rn + LIMB_BITS - 1) / LIMB_BITS;
48 return rn + 2;
49 } else {
50 /*
51 base = b * 2^base_tz
52 */
53 slong base_tz = lmmp_limb_bits_(base);
54 uint32_t b;
55 if (base_tz < 32) {
56 b = base << (32 - base_tz);
57 } else {
58 b = base >> (base_tz - 32);
59 }
60 base_tz = base_tz - 32;
61 base_tz *= exp;
62
63 mp_size_t rn;
64 /*
65 exp = exp' * 2^bits
66 exp*log2(base) = exp*log2(b*2^base_tz)
67 = exp*log2(b) + exp*base_tz
68 = exp'*log2(b)*2^bits + exp*base_tz
69 */
70 mp_bitcnt_t bits = lmmp_limb_bits_(exp);
71 bits -= 32;
72 exp >>= bits;
73 exp++;
74 rn = xlog2n_ceil(exp, b) << bits;
75 rn += base_tz;
76 rn = (rn + LIMB_BITS - 1) / LIMB_BITS;
77 return rn + 2;
78 }
79}
xlog2n_ceil
static uint64_t xlog2n_ceil(uint32_t x, uint32_t n)
计算x*log2(n)的ceil值
Definition lglg.h:67
mp_bitcnt_t
size_t mp_bitcnt_t
Definition lmmp.h:82
LIMB_BITS
#define LIMB_BITS
Definition lmmp.h:86
MP_UINT_MAX
#define MP_UINT_MAX
Definition mparam.h:136
slong
int64_t slong
Definition numth.h:34

引用了 LIMB_BITS, lmmp_limb_bits_, lmmp_param_assert, MP_UINT_MAX, n , 以及 xlog2n_ceil().

被这些函数引用 lmmp_arith_seqprod_size_(), lmmp_nCr_size_(), lmmp_nPr_size_(), lmmp_pow_size_() , 以及 pow_nPr_().

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

◆ lmmp_pow_size_()

mp_size_t lmmp_pow_size_ ( mp_srcptr  base,
mp_size_t  n,
ulong  exp 
)

计算幂次方需要的limb缓冲区长度 [base,n] ^ exp

参数
base底数指针
n底数 limb 长度
exp指数
警告
n>0, base[n-1]!=0, [base,n]>1
返回
返回值为 [base,n]^exp 需要的 limb 缓冲区长度（比实际长度多）

在文件 pow.c 第 81 行定义.

81 {
82 lmmp_param_assert(n > 0);
83 lmmp_param_assert(base[n - 1] != 0);
84 if (n == 1) {
85 return lmmp_pow_1_size_(base[0], exp);
86 }
87 if (exp == 1) {
88 return n;
89 } else if (exp == 2) {
90 return n * 2;
91 } else {
92 /*
93 base = b * 2^base_tz
94 */
95 mp_bitcnt_t base_tz = lmmp_limb_bits_(base[n - 1]);
96 uint32_t b;
97 if (base_tz < 32) {
98 b = base[n - 1] << (32 - base_tz);
99 b |= (base[n - 2] >> (LIMB_BITS - 32 + base_tz));
100 base_tz = (n - 2) * LIMB_BITS + LIMB_BITS - 32 + base_tz;
101 } else if (base_tz == 32) {
102 b = base[n - 1];
103 base_tz = (n - 1) * LIMB_BITS;
104 } else {
105 b = base[n - 1] >> (base_tz - 32);
106 base_tz = (n - 1) * LIMB_BITS + base_tz - 32;
107 }
108
109 mp_size_t rn;
110 if (exp <= MP_UINT_MAX) {
111 rn = exp * base_tz;
112 rn += xlog2n_ceil(exp, b);
113 } else {
114 /*
115 exp = exp' * 2^bits
116 exp*log2(base) = exp*log2(b*2^base_tz)
117 = exp*log2(b) + exp*base_tz
118 = exp'*log2(b)*2^bits + exp*base_tz
119 */
120 mp_bitcnt_t bits = lmmp_limb_bits_(exp);
121 rn = exp * base_tz;
122 bits -= 32;
123 exp >>= bits;
124 exp++;
125 rn += xlog2n_ceil(exp, b) << bits;
126 }
127 rn = (rn + LIMB_BITS - 1) / LIMB_BITS;
128 return rn + 2;
129 }
130}
lmmp_pow_1_size_
mp_size_t lmmp_pow_1_size_(mp_limb_t base, ulong exp)
Copyright (C) 2026 HJimmyK(Jericho Knox)
Definition pow.c:25

引用了 LIMB_BITS, lmmp_limb_bits_, lmmp_param_assert, lmmp_pow_1_size_(), MP_UINT_MAX, n , 以及 xlog2n_ceil().

被这些函数引用 lmmp_pow_().

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