LAMMP 4.2.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
mul_basecase.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/lmmpn.h"
18
19
20void lmmp_sqr_basecase_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na) {
21 lmmp_param_assert(na >= 1);
22
23 mp_size_t i;
24 mp_limb_t cl, x;
25
26 x = numa[0];
27 cl = 0;
28
29 for (i = 0; i + 4 <= na; i += 4) {
30 mp_limb_t u0 = numa[i + 0], u1 = numa[i + 1], u2 = numa[i + 2], u3 = numa[i + 3];
31 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
32
33 _umul64to128_(u0, x, &l0, &h0);
34 l0 += cl;
35 cl = (l0 < cl) + h0;
36 _umul64to128_(u1, x, &l1, &h1);
37 l1 += cl;
38 cl = (l1 < cl) + h1;
39 _umul64to128_(u2, x, &l2, &h2);
40 l2 += cl;
41 cl = (l2 < cl) + h2;
42 _umul64to128_(u3, x, &l3, &h3);
43 l3 += cl;
44 cl = (l3 < cl) + h3;
45
46 dst[i + 0] = l0;
47 dst[i + 1] = l1;
48 dst[i + 2] = l2;
49 dst[i + 3] = l3;
50 }
51 for (; i < na; i++) {
52 mp_limb_t u = numa[i], l, h;
53 _umul64to128_(u, x, &l, &h);
54 l += cl;
55 cl = (l < cl) + h;
56 dst[i] = l;
57 }
58 dst[na] = cl;
59
60 dst++;
61 mp_size_t nb = na - 1;
62 mp_limb_t* nb_ptr = (mp_limb_t*)numa + 1;
63
64 while (nb >= 2) {
65 cl = 0;
66 x = nb_ptr[0];
67 for (i = 0; i + 4 <= na; i += 4) {
68 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
69 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
70 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
71 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
72 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
73
74 _umul64to128_(u0, x, &l0, &h0);
75 l0 += cl;
76 cl = (l0 < cl) + h0;
77 l0 += d0;
78 cl += (l0 < d0);
79 dst[i + 0] = l0;
80 _umul64to128_(u1, x, &l1, &h1);
81 l1 += cl;
82 cl = (l1 < cl) + h1;
83 l1 += d1;
84 cl += (l1 < d1);
85 dst[i + 1] = l1;
86 _umul64to128_(u2, x, &l2, &h2);
87 l2 += cl;
88 cl = (l2 < cl) + h2;
89 l2 += d2;
90 cl += (l2 < d2);
91 dst[i + 2] = l2;
92 _umul64to128_(u3, x, &l3, &h3);
93 l3 += cl;
94 cl = (l3 < cl) + h3;
95 l3 += d3;
96 cl += (l3 < d3);
97 dst[i + 3] = l3;
98 }
99 for (; i < na; i++) {
100 mp_limb_t u = numa[i], d = dst[i], l, h;
101 _umul64to128_(u, x, &l, &h);
102 l += cl;
103 cl = (l < cl) + h;
104 l += d;
105 cl += (l < d);
106 dst[i] = l;
107 }
108 dst[na] = cl;
109 dst++;
110
111 cl = 0;
112 x = nb_ptr[1];
113 for (i = 0; i + 4 <= na; i += 4) {
114 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
115 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
116 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
117 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
118 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
119
120 _umul64to128_(u0, x, &l0, &h0);
121 l0 += cl;
122 cl = (l0 < cl) + h0;
123 l0 += d0;
124 cl += (l0 < d0);
125 dst[i + 0] = l0;
126 _umul64to128_(u1, x, &l1, &h1);
127 l1 += cl;
128 cl = (l1 < cl) + h1;
129 l1 += d1;
130 cl += (l1 < d1);
131 dst[i + 1] = l1;
132 _umul64to128_(u2, x, &l2, &h2);
133 l2 += cl;
134 cl = (l2 < cl) + h2;
135 l2 += d2;
136 cl += (l2 < d2);
137 dst[i + 2] = l2;
138 _umul64to128_(u3, x, &l3, &h3);
139 l3 += cl;
140 cl = (l3 < cl) + h3;
141 l3 += d3;
142 cl += (l3 < d3);
143 dst[i + 3] = l3;
144 }
145 for (; i < na; i++) {
146 mp_limb_t u = numa[i], d = dst[i], l, h;
147 _umul64to128_(u, x, &l, &h);
148 l += cl;
149 cl = (l < cl) + h;
150 l += d;
151 cl += (l < d);
152 dst[i] = l;
153 }
154 dst[na] = cl;
155 dst++;
156
157 nb_ptr += 2;
158 nb -= 2;
159 }
160
161 while (nb >= 1) {
162 cl = 0;
163 x = nb_ptr[0];
164 for (i = 0; i + 4 <= na; i += 4) {
165 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
166 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
167 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
168 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
169 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
170
171 _umul64to128_(u0, x, &l0, &h0);
172 l0 += cl;
173 cl = (l0 < cl) + h0;
174 l0 += d0;
175 cl += (l0 < d0);
176 dst[i + 0] = l0;
177 _umul64to128_(u1, x, &l1, &h1);
178 l1 += cl;
179 cl = (l1 < cl) + h1;
180 l1 += d1;
181 cl += (l1 < d1);
182 dst[i + 1] = l1;
183 _umul64to128_(u2, x, &l2, &h2);
184 l2 += cl;
185 cl = (l2 < cl) + h2;
186 l2 += d2;
187 cl += (l2 < d2);
188 dst[i + 2] = l2;
189 _umul64to128_(u3, x, &l3, &h3);
190 l3 += cl;
191 cl = (l3 < cl) + h3;
192 l3 += d3;
193 cl += (l3 < d3);
194 dst[i + 3] = l3;
195 }
196 for (; i < na; i++) {
197 mp_limb_t u = numa[i], d = dst[i], l, h;
198 _umul64to128_(u, x, &l, &h);
199 l += cl;
200 cl = (l < cl) + h;
201 l += d;
202 cl += (l < d);
203 dst[i] = l;
204 }
205 dst[na] = cl;
206 dst++;
207 nb_ptr++;
208 nb--;
209 }
210}
211
212void lmmp_mul_basecase_(
213 mp_ptr restrict dst,
214 mp_srcptr restrict numa,
215 mp_size_t na,
216 mp_srcptr restrict numb,
217 mp_size_t nb
218) {
219 lmmp_param_assert(na >= nb);
220 lmmp_param_assert(nb >= 1);
221
222 mp_size_t i;
223 mp_limb_t cl;
224
225 cl = 0;
226 mp_limb_t x = numb[0];
227 for (i = 0; i + 4 <= na; i += 4) {
228 mp_limb_t u0 = numa[i + 0], u1 = numa[i + 1], u2 = numa[i + 2], u3 = numa[i + 3];
229 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
230
231 _umul64to128_(u0, x, &l0, &h0);
232 l0 += cl;
233 cl = (l0 < cl) + h0;
234 _umul64to128_(u1, x, &l1, &h1);
235 l1 += cl;
236 cl = (l1 < cl) + h1;
237 _umul64to128_(u2, x, &l2, &h2);
238 l2 += cl;
239 cl = (l2 < cl) + h2;
240 _umul64to128_(u3, x, &l3, &h3);
241 l3 += cl;
242 cl = (l3 < cl) + h3;
243
244 dst[i + 0] = l0;
245 dst[i + 1] = l1;
246 dst[i + 2] = l2;
247 dst[i + 3] = l3;
248 }
249 for (; i < na; i++) {
250 mp_limb_t u = numa[i], l, h;
251 _umul64to128_(u, x, &l, &h);
252 l += cl;
253 cl = (l < cl) + h;
254 dst[i] = l;
255 }
256 dst[na] = cl;
257 dst++;
258 numb++;
259 nb--;
260
261 while (nb >= 2) {
262 // 第一个乘数
263 cl = 0;
264 x = numb[0];
265 for (i = 0; i + 4 <= na; i += 4) {
266 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
267 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
268 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
269 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
270 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
271
272 _umul64to128_(u0, x, &l0, &h0);
273 l0 += cl;
274 cl = (l0 < cl) + h0;
275 l0 += d0;
276 cl += (l0 < d0);
277 dst[i + 0] = l0;
278 _umul64to128_(u1, x, &l1, &h1);
279 l1 += cl;
280 cl = (l1 < cl) + h1;
281 l1 += d1;
282 cl += (l1 < d1);
283 dst[i + 1] = l1;
284 _umul64to128_(u2, x, &l2, &h2);
285 l2 += cl;
286 cl = (l2 < cl) + h2;
287 l2 += d2;
288 cl += (l2 < d2);
289 dst[i + 2] = l2;
290 _umul64to128_(u3, x, &l3, &h3);
291 l3 += cl;
292 cl = (l3 < cl) + h3;
293 l3 += d3;
294 cl += (l3 < d3);
295 dst[i + 3] = l3;
296 }
297 for (; i < na; i++) {
298 mp_limb_t u = numa[i], d = dst[i], l, h;
299 _umul64to128_(u, x, &l, &h);
300 l += cl;
301 cl = (l < cl) + h;
302 l += d;
303 cl += (l < d);
304 dst[i] = l;
305 }
306 dst[na] = cl;
307
308 dst++;
309 cl = 0;
310 x = numb[1];
311 for (i = 0; i + 4 <= na; i += 4) {
312 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
313 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
314 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
315 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
316 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
317
318 _umul64to128_(u0, x, &l0, &h0);
319 l0 += cl;
320 cl = (l0 < cl) + h0;
321 l0 += d0;
322 cl += (l0 < d0);
323 dst[i + 0] = l0;
324 _umul64to128_(u1, x, &l1, &h1);
325 l1 += cl;
326 cl = (l1 < cl) + h1;
327 l1 += d1;
328 cl += (l1 < d1);
329 dst[i + 1] = l1;
330 _umul64to128_(u2, x, &l2, &h2);
331 l2 += cl;
332 cl = (l2 < cl) + h2;
333 l2 += d2;
334 cl += (l2 < d2);
335 dst[i + 2] = l2;
336 _umul64to128_(u3, x, &l3, &h3);
337 l3 += cl;
338 cl = (l3 < cl) + h3;
339 l3 += d3;
340 cl += (l3 < d3);
341 dst[i + 3] = l3;
342 }
343 for (; i < na; i++) {
344 mp_limb_t u = numa[i], d = dst[i], l, h;
345 _umul64to128_(u, x, &l, &h);
346 l += cl;
347 cl = (l < cl) + h;
348 l += d;
349 cl += (l < d);
350 dst[i] = l;
351 }
352 dst[na] = cl;
353
354 dst++;
355 numb += 2;
356 nb -= 2;
357 }
358
359 while (nb >= 1) {
360 cl = 0;
361 x = numb[0];
362 for (i = 0; i + 4 <= na; i += 4) {
363 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
364 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
365 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
366 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
367 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
368
369 _umul64to128_(u0, x, &l0, &h0);
370 l0 += cl;
371 cl = (l0 < cl) + h0;
372 l0 += d0;
373 cl += (l0 < d0);
374 dst[i + 0] = l0;
375 _umul64to128_(u1, x, &l1, &h1);
376 l1 += cl;
377 cl = (l1 < cl) + h1;
378 l1 += d1;
379 cl += (l1 < d1);
380 dst[i + 1] = l1;
381 _umul64to128_(u2, x, &l2, &h2);
382 l2 += cl;
383 cl = (l2 < cl) + h2;
384 l2 += d2;
385 cl += (l2 < d2);
386 dst[i + 2] = l2;
387 _umul64to128_(u3, x, &l3, &h3);
388 l3 += cl;
389 cl = (l3 < cl) + h3;
390 l3 += d3;
391 cl += (l3 < d3);
392 dst[i + 3] = l3;
393 }
394 for (; i < na; i++) {
395 mp_limb_t u = numa[i], d = dst[i], l, h;
396 _umul64to128_(u, x, &l, &h);
397 l += cl;
398 cl = (l < cl) + h;
399 l += d;
400 cl += (l < d);
401 dst[i] = l;
402 }
403 dst[na] = cl;
404 dst++;
405 numb++;
406 nb--;
407 }
408}
l
#define l
mp_ptr
mp_limb_t * mp_ptr
Definition lmmp.h:117
mp_size_t
uint64_t mp_size_t
Definition lmmp.h:114
mp_srcptr
const mp_limb_t * mp_srcptr
Definition lmmp.h:118
mp_limb_t
uint64_t mp_limb_t
Definition lmmp.h:113
lmmp_param_assert
#define lmmp_param_assert(x)
Definition lmmp.h:496
_umul64to128_
static void _umul64to128_(uint64_t a, uint64_t b, uint64_t *low, uint64_t *high)
Definition longlong.h:174
lmmp_mul_basecase_
void lmmp_mul_basecase_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_srcptr restrict numb, mp_size_t nb)
Definition mul_basecase.c:212
lmmp_sqr_basecase_
void lmmp_sqr_basecase_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na)
Copyright (C) 2026 HJimmyK(Jericho Knox)
Definition mul_basecase.c:20
numb
#define numb
n
#define n