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

浏览源代码.

函数

void lmmp_div_ (mp_ptr dstq, mp_ptr dstr, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
 大数除法和取模操作
 
mp_limb_t lmmp_div_s_ (mp_ptr dstq, mp_ptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
 除法运算
 

函数说明

◆ lmmp_div_()

void lmmp_div_ ( mp_ptr  dstq,
mp_ptr  dstr,
mp_srcptr  numa,
mp_size_t  na,
mp_srcptr  numb,
mp_size_t  nb 
)

大数除法和取模操作

注解
如果dstq不为NULL: [dstq,na-nb+1] = [numa,na] / [numb,nb] (商) 如果dstr不为NULL: [dstr,nb] = [numa,na] mod [numb,nb] (余数)
警告
0<nb<=na, numb[nb-1]!=0, sep(dstq,[numa|numb]), eqsep(dstr,[numa|numb])) 特殊情况: nb==1时, dstq>=numa-1 是允许的 nb==2时, dstq>=numa 是允许的
参数
dstq商结果输出指针（NULL表示不计算商）
dstr余数结果输出指针（NULL表示不计算余数）
numa被除数指针
na被除数的 limb 长度
numb除数指针
nb除数的 limb 长度

在文件 div.c 第 57 行定义.

57 {
58 if (nb == 1) {
59 mp_limb_t rem = lmmp_div_1_(dstq, numa, na, *numb);
60 if (dstr)
61 *dstr = rem;
62 } else if (nb == 2) {
63 mp_limb_t brem[2];
64 brem[0] = numb[0];
65 brem[1] = numb[1];
66 lmmp_div_2_(dstq, numa, na, brem);
67 if (dstr) {
68 dstr[0] = brem[0];
69 dstr[1] = brem[1];
70 }
71 } else {
72 int adjust = numa[na - 1] >= numb[nb - 1];
73 int cnt = lmmp_leading_zeros_(numb[nb - 1]);
74 mp_size_t nq = na + adjust - nb;
75 if (nq == 0) {
76 if (dstr && dstr != numa)
77 lmmp_copy(dstr, numa, nb);
78 if (dstq)
79 dstq[0] = 0;
80 return;
81 }
82 TEMP_DECL;
83
84 if (!dstq)
85 dstq = TALLOC_TYPE(na - nb + 1, mp_limb_t);
86 dstq[na - nb] = 0;
87
88 if (nq >= nb) {
89 mp_ptr numa2 = TALLOC_TYPE(na + 1, mp_limb_t);
90 mp_ptr numb2;
91 if (cnt) {
92 numa2[na] = lmmp_shl_(numa2, numa, na, cnt);
93 numb2 = TALLOC_TYPE(nb, mp_limb_t);
94 lmmp_shl_(numb2, numb, nb, cnt);
95 } else {
96 numa2[na] = 0;
97 lmmp_copy(numa2, numa, na);
98 numb2 = (mp_ptr)numb;
99 }
100
101 mp_limb_t inv21 = lmmp_inv_2_1_(numb2[nb - 1], numb2[nb - 2]);
102 na += adjust;
103
104 if (nb < DIV_DIVIDE_THRESHOLD)
105 lmmp_div_basecase_(dstq, numa2, na, numb2, nb, inv21);
106 else if (nb < DIV_MULINV_L_THRESHOLD || na < 2 * DIV_MULINV_N_THRESHOLD)
107 lmmp_div_divide_(dstq, numa2, na, numb2, nb, inv21);
108 else {
109 mp_limb_t ni = lmmp_div_inv_size_(nq, nb);
110 mp_ptr invappr = TALLOC_TYPE(ni, mp_limb_t);
111 lmmp_inv_prediv_(invappr, numb2, nb, ni);
112 lmmp_div_mulinv_(dstq, numa2, na, numb2, nb, invappr, ni);
113 }
114
115 if (dstr) {
116 if (cnt)
117 lmmp_shr_(dstr, numa2, nb, cnt);
118 else
119 lmmp_copy(dstr, numa2, nb);
120 }
121 } else {
122 // nq=na-nb+adj<nb
123 //-> na+adj>=2nq+1
124 mp_size_t ni = nb - nq;
125 mp_ptr numa2, numb2;
126 mp_ptr tp = TALLOC_TYPE(nb, mp_limb_t);
127 mp_limb_t cy;
128
129 numa2 = TALLOC_TYPE(nq * 2 + 1, mp_limb_t);
130 if (cnt) {
131 numb2 = TALLOC_TYPE(nq, mp_limb_t);
132 lmmp_shl_(numb2, numb + ni, nq, cnt);
133 numb2[0] |= numb[ni - 1] >> (LIMB_BITS - cnt);
134 cy = lmmp_shl_(numa2, numa + na - 2 * nq, 2 * nq, cnt);
135 if (adjust) {
136 numa2[2 * nq] = cy;
137 ++numa2; // numa2[0] is as significant as numa[ni=na-2nq+adjust]
138 } else
139 numa2[0] |= numa[na - 2 * nq - 1] >> (LIMB_BITS - cnt);
140 } else {
141 numb2 = (mp_ptr)numb + ni;
142 lmmp_copy(numa2, numa + na - 2 * nq, 2 * nq);
143 if (adjust) {
144 numa2[2 * nq] = 0;
145 ++numa2;
146 }
147 }
148
149 // now: 0<=numa2<B^2nq, B^nq/2<=numb2<B^nq, and 0<=numa2/numb2<B^nq
150 // ignored bits could be seen as fraction part of numa and numb
151 // we can prove: Q<=Qh<=Q+2
152 // where Q=floor(numa/numb) is the real quotient
153 // Qh=floor(floor(numa)/floor(numb)) as below
154
155 if (nq == 1) {
156 lmmp_div_1_s_(dstq, numa2, 2, *numb2);
157 } else if (nq == 2) {
158 lmmp_div_2_s_(dstq, numa2, 4, numb2);
159 } else {
160 mp_limb_t inv21 = lmmp_inv_2_1_(numb2[nq - 1], numb2[nq - 2]);
161
162 if (nq < DIV_DIVIDE_THRESHOLD)
163 lmmp_div_basecase_(dstq, numa2, 2 * nq, numb2, nq, inv21);
164 else if (nq < DIV_MULINV_N_THRESHOLD)
165 lmmp_div_divide_(dstq, numa2, 2 * nq, numb2, nq, inv21);
166 else {
167 mp_limb_t ni = lmmp_div_inv_size_(nq, nq);
168 mp_ptr invappr = tp;
169 lmmp_inv_prediv_(invappr, numb2, nq, ni);
170 lmmp_div_mulinv_(dstq, numa2, 2 * nq, numb2, nq, invappr, ni);
171 }
172 }
173 /*
174 true remainder = partial remainder - quotient * ignored divisor limbs
175
176 Multiply the first ignored divisor limb by the most significant
177 quotient limb. If that product is > the partial remainder's
178 most significant limb, we know the quotient is too large. This
179 test quickly catches most cases where the quotient is too large;
180 it catches all cases where the quotient is 2 too large.*/
181
182 mp_limb_t x;
183 if (cnt) {
184 mp_limb_t dl;
185 if (ni < 2)
186 dl = 0;
187 else
188 dl = numb[ni - 2];
189 x = (numb[ni - 1] << cnt) | (dl >> (LIMB_BITS - cnt));
190 } else
191 x = numb[ni - 1];
192 mp_limb_t h = (x >> LIMB_BITS / 2) * (dstq[nq - 1] >> LIMB_BITS / 2);
193 mp_limb_t rnb = 0; // remainder[nb]
194 mp_size_t nr = nq; // remainder=rnb:[numa2,nr]:[...,ni]
195
196 if (h > numa2[nq - 1]) {
197 lmmp_dec(dstq);
198 rnb = lmmp_add_n_(numa2, numa2, numb2, nq);
199 }
200
201 // if cnt, recover the shift of partial remainder
202 // and remove the effect of the partial-ignored numa[ni-1] and numb[ni-1]
203 if (cnt) {
204 numa2[nq] = rnb;
205 ++nr;
206 --ni;
207 lmmp_shl_(numa2, numa2, nr, LIMB_BITS - cnt);
208 numa2[0] |= numa[ni] & (LIMB_MAX >> cnt);
209 cy = lmmp_submul_1_(numa2, dstq, nq, numb[ni] & (LIMB_MAX >> cnt));
210 rnb = -(numa2[nq] < cy);
211 numa2[nq] -= cy;
212 }
213
214 if (ni == 0) {
215 if (dstr) {
216 if (rnb)
217 lmmp_add_n_(dstr, numa2, numb, nr);
218 else
219 lmmp_copy(dstr, numa2, nr);
220 }
221 } else {
222 tp[nb - 1] = 0;
223 if (ni < nq)
224 lmmp_mul_(tp, dstq, nq, numb, ni);
225 else
226 lmmp_mul_(tp, numb, ni, dstq, nq);
227
228 if (dstr) {
229 mp_ptr remptr = dstr == numb ? tp : dstr;
230 cy = lmmp_sub_n_(remptr, numa, tp, ni);
231 rnb -= lmmp_sub_nc_(remptr + ni, numa2, tp + ni, nr, cy);
232 if (rnb)
233 lmmp_add_n_(dstr, remptr, numb, nb);
234 else if (dstr != remptr)
235 lmmp_copy(dstr, remptr, nb);
236 } else {
237 int hcmp = lmmp_cmp_(numa2, tp + ni, nr);
238 if (hcmp < 0)
239 --rnb;
240 else if (hcmp == 0)
241 rnb -= (lmmp_cmp_(numa, tp, ni) < 0);
242 }
243 }
244
245 if (rnb)
246 lmmp_dec(dstq);
247 }
248
249 TEMP_FREE;
250 }
251}
mp_ptr
mp_limb_t * mp_ptr
Definition lmmp.h:215
lmmp_copy
#define lmmp_copy(dst, src, n)
Definition lmmp.h:364
mp_size_t
uint64_t mp_size_t
Definition lmmp.h:212
LIMB_MAX
#define LIMB_MAX
Definition lmmp.h:224
mp_limb_t
uint64_t mp_limb_t
Definition lmmp.h:211
LIMB_BITS
#define LIMB_BITS
Definition lmmp.h:221
lmmp_div_inv_size_
static mp_size_t lmmp_div_inv_size_(mp_size_t nq, mp_size_t nb)
计算预计算逆元的尺寸
Definition lmmpn.h:812
lmmp_div_1_s_
mp_limb_t lmmp_div_1_s_(mp_ptr dstq, mp_ptr numa, mp_size_t na, mp_limb_t x)
单精度数除法（除数为1个limb）
lmmp_div_1_
mp_limb_t lmmp_div_1_(mp_ptr dstq, mp_srcptr numa, mp_size_t na, mp_limb_t x)
单精度数除法
Definition div.c:66
lmmp_leading_zeros_
int lmmp_leading_zeros_(mp_limb_t x)
计算一个单精度数(limb)中前导零的个数
Definition tiny.c:35
lmmp_cmp_
static int lmmp_cmp_(mp_srcptr numa, mp_srcptr numb, mp_size_t n)
大数比较函数（内联）
Definition lmmpn.h:1004
lmmp_dec
#define lmmp_dec(p)
大数减1宏（预期无借位）
Definition lmmpn.h:973
lmmp_inv_prediv_
void lmmp_inv_prediv_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_size_t ni)
除法前的逆元预计算，[dst,ni] = invappr( (ni+1 MSLs of numa) + 1 ) / B
Definition div_mulinv.c:11
lmmp_div_2_
void lmmp_div_2_(mp_ptr dstq, mp_srcptr numa, mp_size_t na, mp_ptr numb)
双精度数除法 (除数为2个limb)
Definition div.c:223
lmmp_shr_
mp_limb_t lmmp_shr_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_size_t shr)
大数右移操作 [dst,na] = [numa,na] >> shr，dst的高shr位填充0
Definition shr.c:9
lmmp_mul_
void lmmp_mul_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
不等长大数乘法操作 [dst,na+nb] = [numa,na] * [numb,nb]
lmmp_shl_
mp_limb_t lmmp_shl_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_size_t shl)
大数左移操作 [dst,na] = [numa,na]<<shl，dst的低shl位填充0
Definition shl.c:9
lmmp_div_mulinv_
mp_limb_t lmmp_div_mulinv_(mp_ptr dstq, mp_ptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb, mp_srcptr invappr, mp_size_t ni)
乘法逆元除法
Definition div_mulinv.c:36
lmmp_div_2_s_
mp_limb_t lmmp_div_2_s_(mp_ptr dstq, mp_ptr numa, mp_size_t na, mp_srcptr numb)
双精度数除法（除数为2个limb）
lmmp_submul_1_
mp_limb_t lmmp_submul_1_(mp_ptr numa, mp_srcptr numb, mp_size_t n, mp_limb_t b)
大数乘以单limb并累减操作 [numa,n] -= [numb,n] * b
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_inv_2_1_
mp_limb_t lmmp_inv_2_1_(mp_limb_t xh, mp_limb_t xl)
2-1阶逆元计算 (inv21)
Definition inv.c:10
lmmp_div_basecase_
mp_limb_t lmmp_div_basecase_(mp_ptr dstq, mp_ptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb, mp_limb_t inv21)
基础除法运算
Definition div_basecase.c:10
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
lmmp_div_divide_
mp_limb_t lmmp_div_divide_(mp_ptr dstq, mp_ptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb, mp_limb_t inv21)
分治除法运算
Definition div_divide.c:53
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
DIV_DIVIDE_THRESHOLD
#define DIV_DIVIDE_THRESHOLD
Definition mparam.h:26
DIV_MULINV_N_THRESHOLD
#define DIV_MULINV_N_THRESHOLD
Definition mparam.h:30
DIV_MULINV_L_THRESHOLD
#define DIV_MULINV_L_THRESHOLD
Definition mparam.h:28
tp
#define tp
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

引用了 DIV_DIVIDE_THRESHOLD, DIV_MULINV_L_THRESHOLD, DIV_MULINV_N_THRESHOLD, LIMB_BITS, LIMB_MAX, lmmp_add_n_(), lmmp_cmp_(), lmmp_copy, lmmp_dec, lmmp_div_1_(), lmmp_div_1_s_(), lmmp_div_2_(), lmmp_div_2_s_(), lmmp_div_basecase_(), lmmp_div_divide_(), lmmp_div_inv_size_(), lmmp_div_mulinv_(), lmmp_inv_2_1_(), lmmp_inv_prediv_(), lmmp_leading_zeros_(), lmmp_mul_(), lmmp_shl_(), lmmp_shr_(), lmmp_sub_n_(), lmmp_sub_nc_(), lmmp_submul_1_(), TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

被这些函数引用 lmmp_bninv_(), lmmp_gcd_basecase_(), lmmp_gcd_lehmer_(), lmmp_trialdiv_() , 以及 try_div_().

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

◆ lmmp_div_s_()

mp_limb_t lmmp_div_s_ ( mp_ptr  dstq,
mp_ptr  numa,
mp_size_t  na,
mp_srcptr  numb,
mp_size_t  nb 
)

除法运算

参数
dstq输出商的缓冲区，长度至少为na-nb
numa输入被除数（长度na），运算后存储余数（长度nb）
na被除数的 limb 长度
numb输入除数，长度为nb
nb除数的 limb 长度
返回
商的最高位（qh）
警告
na>=nb>0, MSB(numb)=1, sep(dstq,numa,numb)
注解
qh:[dstq,na-nb]=[numa,na] div [numb,nb], [numa,nb]=[numa,na] mod [numb,nb], return qh

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

11 {
12 TEMP_DECL;
13 mp_limb_t nq = na - nb;
14 mp_limb_t qh;
15 if (nq == 0) {
16 qh = lmmp_cmp_(numa, numb, nb) >= 0;
17 if (qh)
18 lmmp_sub_n_(numa, numa, numb, nb);
19 } else if (nb == 1) {
20 qh = lmmp_div_1_s_(dstq, numa, na, *numb);
21 } else if (nb == 2) {
22 qh = lmmp_div_2_s_(dstq, numa, na, numb);
23 } else if (nq < nb) {
24 qh = lmmp_div_s_(dstq, numa + na - 2 * nq, 2 * nq, numb + nb - nq, nq);
25
26 mp_ptr tp = TALLOC_TYPE(nb, mp_limb_t);
27 if (nq > nb - nq)
28 lmmp_mul_(tp, dstq, nq, numb, nb - nq);
29 else
30 lmmp_mul_(tp, numb, nb - nq, dstq, nq);
31
32 mp_limb_t cy = lmmp_sub_n_(numa, numa, tp, nb);
33 if (qh)
34 cy += lmmp_sub_n_(numa + nq, numa + nq, numb, nb - nq);
35
36 while (cy) {
37 qh -= lmmp_sub_1_(dstq, dstq, nq, 1);
38 cy -= lmmp_add_n_(numa, numa, numb, nb);
39 }
40 } else {
41 mp_limb_t inv21 = lmmp_inv_2_1_(numb[nb - 1], numb[nb - 2]);
42 if (nb < DIV_DIVIDE_THRESHOLD)
43 qh = lmmp_div_basecase_(dstq, numa, na, numb, nb, inv21);
44 else if (nb < DIV_MULINV_L_THRESHOLD || na < 2 * DIV_MULINV_N_THRESHOLD)
45 qh = lmmp_div_divide_(dstq, numa, na, numb, nb, inv21);
46 else {
47 mp_limb_t ni = lmmp_div_inv_size_(nq, nb);
48 mp_ptr invappr = TALLOC_TYPE(ni, mp_limb_t);
49 lmmp_inv_prediv_(invappr, numb, nb, ni);
50 qh = lmmp_div_mulinv_(dstq, numa, na, numb, nb, invappr, ni);
51 }
52 }
53 TEMP_FREE;
54 return qh;
55}
lmmp_div_s_
mp_limb_t lmmp_div_s_(mp_ptr dstq, mp_ptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
除法运算
Definition div.c:11
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

引用了 DIV_DIVIDE_THRESHOLD, DIV_MULINV_L_THRESHOLD, DIV_MULINV_N_THRESHOLD, lmmp_add_n_(), lmmp_cmp_(), lmmp_div_1_s_(), lmmp_div_2_s_(), lmmp_div_basecase_(), lmmp_div_divide_(), lmmp_div_inv_size_(), lmmp_div_mulinv_(), lmmp_div_s_(), lmmp_inv_2_1_(), lmmp_inv_prediv_(), lmmp_mul_(), lmmp_sub_1_(), lmmp_sub_n_(), TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 tp.

被这些函数引用 lmmp_div_divide_(), lmmp_div_s_(), lmmp_invsqrt_newton_(), lmmp_sqrt_divide_() , 以及 lmmp_to_str_divide_().

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