LAMMP 4.1.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
宏定义 | 函数
mullo.c 文件参考
#include "../../include/lammp/impl/tmp_alloc.h"
#include "../../include/lammp/lmmpn.h"
#include "../../include/lammp/impl/mparam.h"
+ mullo.c 的引用(Include)关系图:

浏览源代码.

宏定义

#define a0   (numa)
 
#define a0   (numa)
 
#define a1   (numa + m)
 
#define a1   (numa + m)
 
#define b0   (numb)
 
#define b1   (numb + m)
 
#define c0   (dst)
 
#define c0   (dst)
 
#define c1   (dst + m)
 
#define c1   (dst + m)
 
#define lo1   (tp)
 
#define lo1   (tp)
 
#define lo2   (tp + 2 * t)
 
#define MUL_TOOM66_THRESHOLD   MUL_FFT_THRESHOLD
 
#define MUL_TOOM88_THRESHOLD   5621
 
#define tp1   (tp + 2 * t)
 
#define tp1   (tp + 2 * t)
 
#define tp2   (tp + 4 * t)
 

函数

void lmmp_mullo_ (mp_ptr restrict dst, mp_srcptr restrict numa, mp_srcptr restrict numb, mp_size_t n)
 
void lmmp_mullo_dc_ (mp_ptr restrict dst, mp_srcptr restrict numa, mp_srcptr restrict numb, mp_ptr restrict tp, mp_size_t n)
 
void lmmp_mullo_fft_ (mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n, mp_ptr scratch)
 低位FFT乘法 [dst,n] = [numa,n] * [numb,n] mod B^n
 
void lmmp_sqrlo_dc_ (mp_ptr restrict dst, mp_srcptr restrict numa, mp_ptr restrict tp, mp_size_t n)
 

宏定义说明

◆ a0 [1/2]

#define a0   (numa)

◆ a0 [2/2]

#define a0   (numa)

◆ a1 [1/2]

#define a1   (numa + m)

◆ a1 [2/2]

#define a1   (numa + m)

◆ b0

#define b0   (numb)

◆ b1

#define b1   (numb + m)

◆ c0 [1/2]

#define c0   (dst)

◆ c0 [2/2]

#define c0   (dst)

◆ c1 [1/2]

#define c1   (dst + m)

◆ c1 [2/2]

#define c1   (dst + m)

◆ lo1 [1/2]

#define lo1   (tp)

◆ lo1 [2/2]

#define lo1   (tp)

◆ lo2

#define lo2   (tp + 2 * t)

◆ MUL_TOOM66_THRESHOLD

#define MUL_TOOM66_THRESHOLD   MUL_FFT_THRESHOLD

在文件 mullo.c 第 105 行定义.

◆ MUL_TOOM88_THRESHOLD

#define MUL_TOOM88_THRESHOLD   5621

在文件 mullo.c 第 106 行定义.

◆ tp1 [1/2]

#define tp1   (tp + 2 * t)

◆ tp1 [2/2]

#define tp1   (tp + 2 * t)

◆ tp2

#define tp2   (tp + 4 * t)

函数说明

◆ lmmp_mullo_()

void lmmp_mullo_ ( mp_ptr restrict  dst,
mp_srcptr restrict  numa,
mp_srcptr restrict  numb,
mp_size_t  n 
)

在文件 mullo.c 第 192 行定义.

192 {
193 lmmp_param_assert(n > 0);
194 if (n < MULLO_DC_THRESHOLD) {
195 if (numa == numb) {
196 TEMP_DECL;
197 mp_ptr restrict tp = TALLOC_TYPE(2 * n, mp_limb_t);
198 lmmp_sqrlo_dc_(dst, numa, tp, n);
199 TEMP_FREE;
200 return;
201 }
202 TEMP_DECL;
203 mp_ptr restrict tp = TALLOC_TYPE(2 * n, mp_limb_t);
204 lmmp_mullo_dc_(dst, numa, numb, tp, n);
205 TEMP_FREE;
206 return;
207 } else {
208 TEMP_DECL;
209 mp_ptr restrict tp = TALLOC_TYPE(2 * n, mp_limb_t);
210 lmmp_mullo_fft_(dst, numa, numb, n, tp);
211 TEMP_FREE;
212 return;
213 }
214}
mp_ptr
mp_limb_t * mp_ptr
Definition lmmp.h:215
mp_limb_t
uint64_t mp_limb_t
Definition lmmp.h:211
lmmp_param_assert
#define lmmp_param_assert(x)
Definition lmmp.h:398
MULLO_DC_THRESHOLD
#define MULLO_DC_THRESHOLD
Definition mparam.h:59
tp
#define tp
lmmp_mullo_dc_
void lmmp_mullo_dc_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_srcptr restrict numb, mp_ptr restrict tp, mp_size_t n)
Definition mullo.c:108
lmmp_sqrlo_dc_
void lmmp_sqrlo_dc_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_ptr restrict tp, mp_size_t n)
Definition mullo.c:157
lmmp_mullo_fft_
void lmmp_mullo_fft_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n, mp_ptr scratch)
低位FFT乘法 [dst,n] = [numa,n] * [numb,n] mod B^n
Definition mullo.c:11
TEMP_DECL
#define TEMP_DECL
Definition tmp_alloc.h:72
TEMP_FREE
#define TEMP_FREE
Definition tmp_alloc.h:93
TALLOC_TYPE
#define TALLOC_TYPE(n, type)
Definition tmp_alloc.h:91

引用了 lmmp_mullo_dc_(), lmmp_mullo_fft_(), lmmp_param_assert, lmmp_sqrlo_dc_(), MULLO_DC_THRESHOLD, TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

+ 函数调用图:

◆ lmmp_mullo_dc_()

void lmmp_mullo_dc_ ( mp_ptr restrict  dst,
mp_srcptr restrict  numa,
mp_srcptr restrict  numb,
mp_ptr restrict  tp,
mp_size_t  n 
)

在文件 mullo.c 第 108 行定义.

114 {
115 if (n < MULLO_BASECASE_THRESHOLD) {
116 lmmp_mul_1_(dst, numa, n, numb[0]);
117 for (mp_size_t i = 1; i < n; ++i) {
118 lmmp_mul_1_(tp, numa, n - i, numb[i]);
119 lmmp_add_n_(dst + i, dst + i, tp, n - i);
120 }
121 return;
122 } else {
123 mp_size_t m, t;
124 if (n < MUL_TOOM33_THRESHOLD) {
125 m = 25 * n / 36;
126 } else if (n < MUL_TOOM44_THRESHOLD) {
127 m = 31 * n / 40;
128 } else if (n < MUL_TOOM66_THRESHOLD) {
129 m = 32 * n / 39;
130 } else if (n < MUL_TOOM88_THRESHOLD) {
131 m = 27 * n / 31;
132 } else {
133 m = 9 * n / 10;
134 }
135 t = n - m;
136
137#define a0 (numa)
138#define a1 (numa + m)
139#define b0 (numb)
140#define b1 (numb + m)
141#define c0 (dst)
142#define c1 (dst + m)
143#define lo1 (tp) // [tp, 2*t]
144#define tp1 (tp + 2 * t) // [tp+2*t, 2*t]
145#define lo2 (tp + 2 * t) // [tp+2*t, 2*t]
146#define tp2 (tp + 4 * t) // [tp+2*t, 2*t]
147 lmmp_mul_n_(tp, a0, b0, m);
148 lmmp_copy(c0, tp, n);
149 lmmp_mullo_dc_(lo1, a1, b0, tp1, t);
150 lmmp_mullo_dc_(lo2, a0, b1, tp2, t);
151 lmmp_add_n_(c1, c1, lo1, t);
152 lmmp_add_n_(c1, c1, lo2, t);
153 return;
154 }
155}
lmmp_copy
#define lmmp_copy(dst, src, n)
Definition lmmp.h:364
mp_size_t
uint64_t mp_size_t
Definition lmmp.h:212
lmmp_mul_n_
void lmmp_mul_n_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n)
等长大数乘法操作 [dst,2*n] = [numa,n] * [numb,n]
Definition mul.c:99
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
lmmp_add_n_
mp_limb_t lmmp_add_n_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n)
无进位的n位加法 [dst,n] = [numa,n] + [numb,n]
Definition add_n.c:71
MULLO_BASECASE_THRESHOLD
#define MULLO_BASECASE_THRESHOLD
Definition mparam.h:57
MUL_TOOM33_THRESHOLD
#define MUL_TOOM33_THRESHOLD
Definition mparam.h:50
MUL_TOOM44_THRESHOLD
#define MUL_TOOM44_THRESHOLD
Definition mparam.h:52
lo2
#define lo2
b0
#define b0
MUL_TOOM66_THRESHOLD
#define MUL_TOOM66_THRESHOLD
Definition mullo.c:105
b1
#define b1
tp2
#define tp2
c1
#define c1
tp1
#define tp1
a0
#define a0
a1
#define a1
c0
#define c0
MUL_TOOM88_THRESHOLD
#define MUL_TOOM88_THRESHOLD
Definition mullo.c:106
lo1
#define lo1

引用了 a0, a1, b0, b1, c0, c1, lmmp_add_n_(), lmmp_copy, lmmp_mul_1_(), lmmp_mul_n_(), lmmp_mullo_dc_(), lo1, lo2, MUL_TOOM33_THRESHOLD, MUL_TOOM44_THRESHOLD, MUL_TOOM66_THRESHOLD, MUL_TOOM88_THRESHOLD, MULLO_BASECASE_THRESHOLD, tp, tp1 , 以及 tp2.

被这些函数引用 lmmp_mullo_(), lmmp_mullo_dc_() , 以及 lmmp_sqrlo_dc_().

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

◆ lmmp_mullo_fft_()

void lmmp_mullo_fft_ ( mp_ptr  dst,
mp_srcptr  numa,
mp_srcptr  numb,
mp_size_t  n,
mp_ptr  scratch 
)

低位FFT乘法 [dst,n] = [numa,n] * [numb,n] mod B^n

参数
dst输出结果缓冲区，长度至少为 n
numa第一个输入操作数，长度为 n
numb第二个输入操作数，长度为 n
scratch临时缓冲区，长度至少为 2*n
n缓冲区 limb 长度
警告
???<n, sep(scratch,[numa|numb]), eqsep(dst,scratch)
返回
无返回值，结果存储在dst中，[dst,n]=[numa,n] * [numb,n] mod B^n

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

11 {
12 lmmp_param_assert(n > 0);
13 mp_size_t hn = lmmp_fft_next_size_((n + n + 1) >> 1);
14 lmmp_assert(n + n > hn);
15 mp_ptr tp = ALLOC_TYPE(hn + 1, mp_limb_t);
16
17 mp_srcptr amodm = numa;
18 mp_size_t nam = n;
19 if (n > hn) {
20 /*
21 Z = B^hb - 1
22 amodm = a mod Z
23 */
24 if (lmmp_add_(scratch, numa, hn, numa + hn, n - hn))
25 lmmp_inc(scratch);
26 amodm = scratch;
27 nam = hn;
28 }
29 lmmp_mul_mersenne_(scratch, hn, amodm, nam, numb, n);
30
31 mp_srcptr amodp = numa;
32 mp_size_t nap = n;
33 if (n > hn) {
34 /*
35 Z = B^hp - 1
36 amodp = a mod Z
37 */
38 tp[hn] = 0;
39 if (lmmp_sub_(tp, numa, hn, numa + hn, n - hn))
40 lmmp_inc(tp);
41 amodp = tp;
42 nap = hn + 1;
43 }
44 lmmp_mul_fermat_(tp, hn, amodp, nap, numb, n);
45
46 mp_limb_t cy = lmmp_shr1add_nc_(scratch, scratch, tp, hn, tp[hn]);
47 cy <<= LIMB_BITS - 1;
48 scratch[hn - 1] += cy;
49 if (scratch[hn - 1] < cy)
50 lmmp_inc(scratch);
51
52 if (n == hn) {
53 cy = tp[hn] + lmmp_sub_n_(scratch + hn, scratch, tp, hn);
54 // cy==1 means [tp,hn+1]!=0, then [dst,hn]!=0
55 // cy==2 is impossible since [tp,hn+1] is normalized.
56 // so the following dec won't overflow.
57 lmmp_dec_1(scratch, cy);
58 } else {
59 mp_size_t n2 = 2 * n;
60 cy = lmmp_sub_n_(scratch + hn, scratch, tp, n2 - hn);
61 cy = tp[hn] + lmmp_sub_nc_(tp + n2 - hn, scratch + n2 - hn, tp + n2 - hn, 2 * hn - n2, cy);
62 cy = lmmp_sub_1_(scratch, scratch, n2, cy);
63 }
64 lmmp_free(tp);
65 lmmp_copy(dst, scratch, n);
66}
scratch
#define scratch
mp_srcptr
const mp_limb_t * mp_srcptr
Definition lmmp.h:216
lmmp_free
void lmmp_free(void *ptr)
内存释放函数（调用lmmp_heap_free_fn）
Definition memory.c:204
lmmp_assert
#define lmmp_assert(x)
Definition lmmp.h:370
LIMB_BITS
#define LIMB_BITS
Definition lmmp.h:221
lmmp_mul_mersenne_
void lmmp_mul_mersenne_(mp_ptr dst, mp_size_t rn, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
梅森数模乘法 [dst,rn] = [numa,na]*[numb,nb] mod B^rn-1
Definition mul_fft.c:752
lmmp_add_
static mp_limb_t lmmp_add_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
大数加法静态内联函数 [dst,na]=[numa,na]+[numb,nb]
Definition lmmpn.h:1058
lmmp_shr1add_nc_
mp_limb_t lmmp_shr1add_nc_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n, mp_limb_t c)
带进位加法后右移1位 [dst,n] = ([numa,n] + [numb,n] + c) >> 1
Definition shr.c:79
lmmp_inc
#define lmmp_inc(p)
大数加1宏（预期无进位）
Definition lmmpn.h:946
lmmp_mul_fermat_
void lmmp_mul_fermat_(mp_ptr dst, mp_size_t rn, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
费马数模乘法 [dst,rn+1]=[numa,na]*[numb,nb] mod B^rn+1
Definition mul_fft.c:677
lmmp_fft_next_size_
mp_size_t lmmp_fft_next_size_(mp_size_t n)
计算满足 >=n 的最小费马/梅森乘法可行尺寸
Definition mul_fft.c:84
lmmp_sub_
static mp_limb_t lmmp_sub_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
大数减法静态内联函数 [dst,na]=[numa,na]-[numb,nb]
Definition lmmpn.h:1072
lmmp_dec_1
#define lmmp_dec_1(p, dec)
大数减指定值宏（预期无借位）
Definition lmmpn.h:985
lmmp_sub_1_
static mp_limb_t lmmp_sub_1_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_limb_t x)
大数减单精度数静态内联函数 [dst,na]=[numa,na]-x
Definition lmmpn.h:1122
lmmp_sub_n_
mp_limb_t lmmp_sub_n_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n)
无借位的n位减法 [dst,n] = [numa,n] - [numb,n]
Definition sub_n.c:70
lmmp_sub_nc_
mp_limb_t lmmp_sub_nc_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n, mp_limb_t c)
带借位的n位减法 [dst,n] = [numa,n] - [numb,n] - c
Definition sub_n.c:9
ALLOC_TYPE
#define ALLOC_TYPE(n, type)
Definition tmp_alloc.h:112

引用了 ALLOC_TYPE, LIMB_BITS, lmmp_add_(), lmmp_assert, lmmp_copy, lmmp_dec_1, lmmp_fft_next_size_(), lmmp_free(), lmmp_inc, lmmp_mul_fermat_(), lmmp_mul_mersenne_(), lmmp_param_assert, lmmp_shr1add_nc_(), lmmp_sub_(), lmmp_sub_1_(), lmmp_sub_n_(), lmmp_sub_nc_(), scratch , 以及 tp.

被这些函数引用 lmmp_mullo_(), lmmp_mullo_n_() , 以及 lmmp_sqrlo_n_().

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

◆ lmmp_sqrlo_dc_()

void lmmp_sqrlo_dc_ ( mp_ptr restrict  dst,
mp_srcptr restrict  numa,
mp_ptr restrict  tp,
mp_size_t  n 
)

在文件 mullo.c 第 157 行定义.

157 {
158 if (n < MULLO_BASECASE_THRESHOLD) {
159 lmmp_mul_1_(dst, numa, n, numa[0]);
160 for (mp_size_t i = 1; i < n; ++i) {
161 lmmp_mul_1_(tp, numa, n - i, numa[i]);
162 lmmp_add_n_(dst + i, dst + i, tp, n - i);
163 }
164 return;
165 } else {
166 mp_size_t m, t;
167 if (n < MUL_TOOM33_THRESHOLD) {
168 m = 25 * n / 36;
169 } else if (n < MUL_TOOM44_THRESHOLD) {
170 m = 31 * n / 40;
171 } else if (n < MUL_TOOM66_THRESHOLD) {
172 m = 32 * n / 39;
173 } else if (n < MUL_TOOM88_THRESHOLD) {
174 m = 27 * n / 31;
175 } else {
176 m = 9 * n / 10;
177 }
178 t = n - m;
179#define a0 (numa)
180#define a1 (numa + m)
181#define c0 (dst)
182#define c1 (dst + m)
183#define lo1 (tp) // [tp, 2*t]
184#define tp1 (tp + 2 * t) // [tp+2*t, 2*t]
185 lmmp_sqr_(tp, a0, m);
186 lmmp_copy(c0, tp, n);
187 lmmp_mullo_dc_(lo1, a0, a1, tp1, t);
188 lmmp_addshl1_n_(c1, c1, lo1, t);
189 }
190}
lmmp_sqr_
void lmmp_sqr_(mp_ptr dst, mp_srcptr numa, mp_size_t na)
大数平方操作 [dst,2*na] = [numa,na]^2
Definition sqr.c:10
lmmp_addshl1_n_
mp_limb_t lmmp_addshl1_n_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n)
加法结合左移1位操作 [dst,n] = [numa,n] + ([numb,n] << 1)
Definition shl.c:56

引用了 a0, a1, c0, c1, lmmp_add_n_(), lmmp_addshl1_n_(), lmmp_copy, lmmp_mul_1_(), lmmp_mullo_dc_(), lmmp_sqr_(), lo1, MUL_TOOM33_THRESHOLD, MUL_TOOM44_THRESHOLD, MUL_TOOM66_THRESHOLD, MUL_TOOM88_THRESHOLD, MULLO_BASECASE_THRESHOLD, tp , 以及 tp1.

被这些函数引用 lmmp_mullo_().

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