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