LAMMP 4.1.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
div.c
浏览该文件的文档.
1/*
2 * LAMMP - Copyright (C) 2025-2026 HJimmyK(Jericho Knox)
3 * This file is part of lammp, under the GNU LGPL v2 license.
4 * See LICENSE in the project root for the full license text.
5 */
6
7#include "../../../include/lammp/lmmpn.h"
8#include "../../../include/lammp/impl/longlong.h"
9
10mp_limb_t lmmp_div_3_2_(mp_ptr restrict numa, mp_srcptr restrict numb, mp_limb_t inv21) {
11 mp_limb_t q;
12 mp_limb_t r1, r0;
13 mp_limb_t a[3] = {numa[0], numa[1], numa[2]};
14 _udiv_qr_3by2(q, r1, r0, a[2], a[1], a[0], numb[1], numb[0], inv21);
15 numa[1] = r1;
16 numa[0] = r0;
17 return q;
18}
19
20mp_limb_t lmmp_mod_1_(mp_srcptr numa, mp_size_t na, mp_limb_t x) {
21 mp_limb_t ah, al;
22 // q: assigned for macro reuse, unused in this logic (known warning)
23 mp_limb_t t = numa[na - 2], q = 0, r = 0;
24 const int shift = lmmp_leading_zeros_(x);
25 if (shift > 0) {
26 const int rshift = LIMB_BITS - shift;
27 ah = numa[na - 1] >> rshift;
28 t = numa[na - 2];
29 al = (numa[na - 1] << shift) | (t >> rshift);
30 x <<= shift;
31 const mp_limb_t inv = lmmp_inv_1_(x);
32 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
33 na -= 2;
34 while (na-- > 0) {
35 ah = r;
36 al = t << shift;
37 t = numa[na];
38 al |= t >> rshift;
39 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
40 }
41 ah = r;
42 al = t << shift;
43 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
44 return r >> shift;
45 } else {
46 ah = 0;
47 t = numa[na - 2];
48 al = numa[na - 1];
49 const mp_limb_t inv = lmmp_inv_1_(x);
50 q = al / x;
51 r = al % x;
52 na -= 2;
53 while (na-- > 0) {
54 ah = r;
55 al = t;
56 t = numa[na];
57 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
58 }
59 ah = r;
60 al = t;
61 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
62 return r;
63 }
64}
65
66mp_limb_t lmmp_div_1_(mp_ptr dstq, mp_srcptr numa, mp_size_t na, mp_limb_t x) {
67 mp_limb_t ah, al;
68 if (na == 1) {
69 ah = numa[0];
70 if (dstq)
71 dstq[0] = ah / x;
72 return ah % x;
73 }
74 if (dstq) {
75 mp_limb_t t = numa[na - 2], q = 0, r = 0;
76 const int shift = lmmp_leading_zeros_(x);
77 if (shift > 0) {
78 /*
79 ah al
80 X|XXXtttX|XXXmmmX|XXXnnnX|XXX----|
81 |000XXXX|tttXXXX|mmmXXXX|nnnXXXX|
82 t numa[na]
83
84 ah al
85 X|XXXtttX|XXXmmmX|XXXnnnX|XXX----|
86 |000XXXX|tttXXXX|mmmXXXX|nnnXXXX|
87 t
88 numa[na]
89 */
90 const int rshift = LIMB_BITS - shift;
91 ah = numa[na - 1] >> rshift;
92 t = numa[na - 2];
93 al = (numa[na - 1] << shift) | (t >> rshift);
94 x <<= shift;
95 const mp_limb_t inv = lmmp_inv_1_(x);
96 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
97 dstq[na - 1] = q;
98 na -= 2;
99 while (na-- > 0) {
100 ah = r;
101 al = t << shift;
102 t = numa[na];
103 al |= t >> rshift;
104 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
105 dstq[na + 1] = q;
106 }
107 ah = r;
108 al = t << shift;
109 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
110 dstq[0] = q;
111 return r >> shift;
112 } else {
113 /*
114 ah al
115 |000XXXX|tttXXXX|mmmXXXX|nnnXXXX|
116 t numa[na]
117 */
118 ah = 0;
119 t = numa[na - 2];
120 al = numa[na - 1];
121 const mp_limb_t inv = lmmp_inv_1_(x);
122 q = al / x;
123 r = al % x;
124 dstq[na - 1] = q;
125 na -= 2;
126 while (na-- > 0) {
127 ah = r;
128 al = t;
129 t = numa[na];
130 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
131 dstq[na + 1] = q;
132 }
133 ah = r;
134 al = t;
135 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
136 dstq[0] = q;
137 return r;
138 }
139 } else {
140 return lmmp_mod_1_(numa, na, x);
141 }
142}
143
144void lmmp_mod_2_(mp_srcptr numa, mp_size_t na, mp_ptr numb) {
145 mp_limb_t q, r1, r0, a2, a1, a0, b1, b0;
146 b1 = numb[1];
147 b0 = numb[0];
148 if (na == 2) {
149 int shift = lmmp_leading_zeros_(b1);
150 if (shift > 0) {
151 const int rshift = LIMB_BITS - shift;
152 b1 = (b1 << shift) | (b0 >> rshift);
153 b0 <<= shift;
154 a2 = numa[1] >> rshift;
155 a1 = (numa[1] << shift) | (numa[0] >> rshift);
156 a0 = (numa[0] << shift);
157 mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
158 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
159 numb[0] = (r0 >> shift) | (r1 << rshift);
160 numb[1] = r1 >> shift;
161 return;
162 } else {
163 if (_u128cmp(numa, numb)) {
164 numb[0] = numa[0];
165 numb[1] = numa[1];
166 return;
167 } else {
168 _u128sub(numb, numa, numb);
169 return;
170 }
171 }
172 } else {
173 int shift = lmmp_leading_zeros_(b1);
174 if (shift > 0) {
175 const int rshift = LIMB_BITS - shift;
176 b1 = (b1 << shift) | (b0 >> rshift);
177 b0 <<= shift;
178 const mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
179 a2 = numa[na - 1] >> rshift;
180 a1 = (numa[na - 1] << shift) | (numa[na - 2] >> rshift);
181 a0 = (numa[na - 2] << shift) | (numa[na - 3] >> rshift);
182 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
183 na -= 2;
184 while (na-- > 1) {
185 a2 = r1;
186 a1 = r0;
187 a0 = (numa[na] << shift) | (numa[na - 1] >> rshift);
188 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
189 }
190
191 a2 = r1;
192 a1 = r0;
193 a0 = (numa[na] << shift);
194 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
195 numb[0] = (r0 >> shift) | (r1 << rshift);
196 numb[1] = r1 >> shift;
197 return;
198 } else {
199 const mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
200 a2 = 0;
201 a1 = numa[na - 1];
202 a0 = numa[na - 2];
203 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
204 na -= 2;
205 while (na-- > 1) {
206 a2 = r1;
207 a1 = r0;
208 a0 = numa[na];
209 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
210 }
211 a2 = r1;
212 a1 = r0;
213 a0 = numa[na];
214 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
215 numb[0] = r0;
216 numb[1] = r1;
217 return;
218 }
219 }
220}
221
222
223void lmmp_div_2_(mp_ptr dstq, mp_srcptr numa, mp_size_t na, mp_ptr numb) {
224 mp_limb_t q, r1, r0, a2, a1, a0, b1, b0;
225 b1 = numb[1];
226 b0 = numb[0];
227 if (na == 2) {
228 int shift = lmmp_leading_zeros_(b1);
229 if (shift > 0) {
230 const int rshift = LIMB_BITS - shift;
231 b1 = (b1 << shift) | (b0 >> rshift);
232 b0 <<= shift;
233 a2 = numa[1] >> rshift;
234 a1 = (numa[1] << shift) | (numa[0] >> rshift);
235 a0 = (numa[0] << shift);
236 mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
237 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
238 if (dstq)
239 dstq[0] = q;
240 numb[0] = (r0 >> shift) | (r1 << rshift);
241 numb[1] = r1 >> shift;
242 return;
243 } else {
244 if (_u128cmp(numa, numb)) {
245 numb[0] = numa[0];
246 numb[1] = numa[1];
247 if (dstq)
248 dstq[0] = 0;
249 return;
250 } else {
251 _u128sub(numb, numa, numb);
252 if (dstq)
253 dstq[0] = 1;
254 return;
255 }
256 }
257 }
258 if (dstq) {
259 int shift = lmmp_leading_zeros_(b1);
260 if (shift > 0) {
261 /*
262 a2 a1 a0
263 X|XXXtttX|XXXmmmX|XXXnnnX|XXX----|
264 |000XXXX|tttXXXX|mmmXXXX|nnnXXXX|
265 numa[na]
266 */
267 const int rshift = LIMB_BITS - shift;
268 b1 = (b1 << shift) | (b0 >> rshift);
269 b0 <<= shift;
270 const mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
271 a2 = numa[na - 1] >> rshift;
272 a1 = (numa[na - 1] << shift) | (numa[na - 2] >> rshift);
273 a0 = (numa[na - 2] << shift) | (numa[na - 3] >> rshift);
274 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
275 dstq[na - 2] = q;
276 na -= 2;
277 while (na-- > 1) {
278 a2 = r1;
279 a1 = r0;
280 a0 = (numa[na] << shift) | (numa[na - 1] >> rshift);
281 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
282 dstq[na] = q;
283 }
284
285 a2 = r1;
286 a1 = r0;
287 a0 = (numa[na] << shift);
288 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
289 dstq[0] = q;
290 numb[0] = (r0 >> shift) | (r1 << rshift);
291 numb[1] = r1 >> shift;
292 return;
293 } else {
294 /*
295 a2 a1 a0
296 |000XXXX|tttXXXX|mmmXXXX|nnnXXXX|
297 numa[na]
298 */
299 const mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
300 a2 = 0;
301 a1 = numa[na - 1];
302 a0 = numa[na - 2];
303 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
304 dstq[na - 2] = q;
305 na -= 2;
306 while (na-- > 1) {
307 a2 = r1;
308 a1 = r0;
309 a0 = numa[na];
310 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
311 dstq[na] = q;
312 }
313 a2 = r1;
314 a1 = r0;
315 a0 = numa[na];
316 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
317 dstq[0] = q;
318 numb[0] = r0;
319 numb[1] = r1;
320 return;
321 }
322 } else {
323 int shift = lmmp_leading_zeros_(b1);
324 if (shift > 0) {
325 const int rshift = LIMB_BITS - shift;
326 b1 = (b1 << shift) | (b0 >> rshift);
327 b0 <<= shift;
328 const mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
329 a2 = numa[na - 1] >> rshift;
330 a1 = (numa[na - 1] << shift) | (numa[na - 2] >> rshift);
331 a0 = (numa[na - 2] << shift) | (numa[na - 3] >> rshift);
332 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
333 na -= 2;
334 while (na-- > 1) {
335 a2 = r1;
336 a1 = r0;
337 a0 = (numa[na] << shift) | (numa[na - 1] >> rshift);
338 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
339 }
340
341 a2 = r1;
342 a1 = r0;
343 a0 = (numa[na] << shift);
344 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
345 numb[0] = (r0 >> shift) | (r1 << rshift);
346 numb[1] = r1 >> shift;
347 return;
348 } else {
349 const mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
350 a2 = 0;
351 a1 = numa[na - 1];
352 a0 = numa[na - 2];
353 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
354 na -= 2;
355 while (na-- > 1) {
356 a2 = r1;
357 a1 = r0;
358 a0 = numa[na];
359 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
360 }
361 a2 = r1;
362 a1 = r0;
363 a0 = numa[na];
364 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
365 numb[0] = r0;
366 numb[1] = r1;
367 return;
368 }
369 }
370}
371
372mp_limb_t lmmp_div_1_s_(mp_ptr restrict dstq, mp_ptr restrict numa, mp_size_t na, mp_limb_t x) {
373 mp_limb_t ah, al;
374 mp_limb_t t = numa[na - 2], q = 0, r = 0;
375 /*
376 ah al
377 |000XXXX|tttXXXX|mmmXXXX|nnnXXXX|
378 t numa[na]
379 */
380 ah = 0;
381 t = numa[na - 2];
382 al = numa[na - 1];
383 const mp_limb_t inv = lmmp_inv_1_(x);
384 q = al / x;
385 r = al % x;
386 const mp_limb_t qh = q;
387 na -= 2;
388 while (na-- > 0) {
389 ah = r;
390 al = t;
391 t = numa[na];
392 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
393 dstq[na + 1] = q;
394 }
395 ah = r;
396 al = t;
397 _udiv_qrnnd_preinv(q, r, ah, al, x, inv);
398 dstq[0] = q;
399 numa[0] = r;
400 return qh;
401}
402
403mp_limb_t lmmp_div_2_s_(mp_ptr restrict dstq, mp_ptr restrict numa, mp_size_t na, mp_srcptr restrict numb) {
404 mp_limb_t q, r1, r0, a2, a1, a0, b1, b0;
405 b1 = numb[1];
406 b0 = numb[0];
407 /*
408 a2 a1 a0
409 |000XXXX|tttXXXX|mmmXXXX|nnnXXXX|
410 numa[na]
411 */
412 const mp_limb_t inv = lmmp_inv_2_1_(b1, b0);
413 a2 = 0;
414 a1 = numa[na - 1];
415 a0 = numa[na - 2];
416 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
417 const mp_limb_t qh = q;
418 na -= 2;
419 while (na-- > 1) {
420 a2 = r1;
421 a1 = r0;
422 a0 = numa[na];
423 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
424 dstq[na] = q;
425 }
426 a2 = r1;
427 a1 = r0;
428 a0 = numa[na];
429 _udiv_qr_3by2(q, r1, r0, a2, a1, a0, b1, b0, inv);
430 dstq[0] = q;
431 numa[0] = r0;
432 numa[1] = r1;
433 return qh;
434}
lmmp_div_3_2_
mp_limb_t lmmp_div_3_2_(mp_ptr restrict numa, mp_srcptr restrict numb, mp_limb_t inv21)
Definition div.c:10
lmmp_mod_1_
mp_limb_t lmmp_mod_1_(mp_srcptr numa, mp_size_t na, mp_limb_t x)
单精度数取余
Definition div.c:20
lmmp_mod_2_
void lmmp_mod_2_(mp_srcptr numa, mp_size_t na, mp_ptr numb)
双精度数取余 (除数为2个limb)
Definition div.c:144
mp_ptr
mp_limb_t * mp_ptr
Definition lmmp.h:215
mp_size_t
uint64_t mp_size_t
Definition lmmp.h:212
mp_srcptr
const mp_limb_t * mp_srcptr
Definition lmmp.h:216
mp_limb_t
uint64_t mp_limb_t
Definition lmmp.h:211
LIMB_BITS
#define LIMB_BITS
Definition lmmp.h:221
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_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_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_inv_1_
mp_limb_t lmmp_inv_1_(mp_limb_t x)
1阶逆元计算 (inv1)
Definition inv.c:107
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
_u128sub
#define _u128sub(r, x, y)
Definition longlong.h:282
_udiv_qrnnd_preinv
#define _udiv_qrnnd_preinv(q, r, nh, nl, d, di)
Definition longlong.h:424
_u128cmp
#define _u128cmp(x, y)
Definition longlong.h:280
_udiv_qr_3by2
#define _udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv)
Definition longlong.h:441
b0
#define b0
b1
#define b1
a0
#define a0
a1
#define a1
r1
#define r1
a2
#define a2
r0
#define r0