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}
211void lmmp_mul_basecase_(
212 mp_ptr restrict dst,
213 mp_srcptr restrict numa,
214 mp_size_t na,
215 mp_srcptr restrict numb,
216 mp_size_t nb
217) {
218 lmmp_param_assert(na >= nb);
219 lmmp_param_assert(nb >= 1);
220
221 mp_size_t i;
222 mp_limb_t cl;
223
224 cl = 0;
225 mp_limb_t x = numb[0];
226 for (i = 0; i + 4 <= na; i += 4) {
227 mp_limb_t u0 = numa[i + 0], u1 = numa[i + 1], u2 = numa[i + 2], u3 = numa[i + 3];
228 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
229
230 _umul64to128_(u0, x, &l0, &h0);
231 l0 += cl;
232 cl = (l0 < cl) + h0;
233 _umul64to128_(u1, x, &l1, &h1);
234 l1 += cl;
235 cl = (l1 < cl) + h1;
236 _umul64to128_(u2, x, &l2, &h2);
237 l2 += cl;
238 cl = (l2 < cl) + h2;
239 _umul64to128_(u3, x, &l3, &h3);
240 l3 += cl;
241 cl = (l3 < cl) + h3;
242
243 dst[i + 0] = l0;
244 dst[i + 1] = l1;
245 dst[i + 2] = l2;
246 dst[i + 3] = l3;
247 }
248 for (; i < na; i++) {
249 mp_limb_t u = numa[i], l, h;
250 _umul64to128_(u, x, &l, &h);
251 l += cl;
252 cl = (l < cl) + h;
253 dst[i] = l;
254 }
255 dst[na] = cl;
256 dst++;
257 numb++;
258 nb--;
259
260 while (nb >= 2) {
261 // 第一个乘数
262 cl = 0;
263 x = numb[0];
264 for (i = 0; i + 4 <= na; i += 4) {
265 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
266 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
267 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
268 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
269 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
270
271 _umul64to128_(u0, x, &l0, &h0);
272 l0 += cl;
273 cl = (l0 < cl) + h0;
274 l0 += d0;
275 cl += (l0 < d0);
276 dst[i + 0] = l0;
277 _umul64to128_(u1, x, &l1, &h1);
278 l1 += cl;
279 cl = (l1 < cl) + h1;
280 l1 += d1;
281 cl += (l1 < d1);
282 dst[i + 1] = l1;
283 _umul64to128_(u2, x, &l2, &h2);
284 l2 += cl;
285 cl = (l2 < cl) + h2;
286 l2 += d2;
287 cl += (l2 < d2);
288 dst[i + 2] = l2;
289 _umul64to128_(u3, x, &l3, &h3);
290 l3 += cl;
291 cl = (l3 < cl) + h3;
292 l3 += d3;
293 cl += (l3 < d3);
294 dst[i + 3] = l3;
295 }
296 for (; i < na; i++) {
297 mp_limb_t u = numa[i], d = dst[i], l, h;
298 _umul64to128_(u, x, &l, &h);
299 l += cl;
300 cl = (l < cl) + h;
301 l += d;
302 cl += (l < d);
303 dst[i] = l;
304 }
305 dst[na] = cl;
306
307 dst++;
308 cl = 0;
309 x = numb[1];
310 for (i = 0; i + 4 <= na; i += 4) {
311 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
312 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
313 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
314 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
315 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
316
317 _umul64to128_(u0, x, &l0, &h0);
318 l0 += cl;
319 cl = (l0 < cl) + h0;
320 l0 += d0;
321 cl += (l0 < d0);
322 dst[i + 0] = l0;
323 _umul64to128_(u1, x, &l1, &h1);
324 l1 += cl;
325 cl = (l1 < cl) + h1;
326 l1 += d1;
327 cl += (l1 < d1);
328 dst[i + 1] = l1;
329 _umul64to128_(u2, x, &l2, &h2);
330 l2 += cl;
331 cl = (l2 < cl) + h2;
332 l2 += d2;
333 cl += (l2 < d2);
334 dst[i + 2] = l2;
335 _umul64to128_(u3, x, &l3, &h3);
336 l3 += cl;
337 cl = (l3 < cl) + h3;
338 l3 += d3;
339 cl += (l3 < d3);
340 dst[i + 3] = l3;
341 }
342 for (; i < na; i++) {
343 mp_limb_t u = numa[i], d = dst[i], l, h;
344 _umul64to128_(u, x, &l, &h);
345 l += cl;
346 cl = (l < cl) + h;
347 l += d;
348 cl += (l < d);
349 dst[i] = l;
350 }
351 dst[na] = cl;
352
353 dst++;
354 numb += 2;
355 nb -= 2;
356 }
357
358 while (nb >= 1) {
359 cl = 0;
360 x = numb[0];
361 for (i = 0; i + 4 <= na; i += 4) {
362 mp_limb_t u0 = numa[i + 0], d0 = dst[i + 0];
363 mp_limb_t u1 = numa[i + 1], d1 = dst[i + 1];
364 mp_limb_t u2 = numa[i + 2], d2 = dst[i + 2];
365 mp_limb_t u3 = numa[i + 3], d3 = dst[i + 3];
366 mp_limb_t l0, h0, l1, h1, l2, h2, l3, h3;
367
368 _umul64to128_(u0, x, &l0, &h0);
369 l0 += cl;
370 cl = (l0 < cl) + h0;
371 l0 += d0;
372 cl += (l0 < d0);
373 dst[i + 0] = l0;
374 _umul64to128_(u1, x, &l1, &h1);
375 l1 += cl;
376 cl = (l1 < cl) + h1;
377 l1 += d1;
378 cl += (l1 < d1);
379 dst[i + 1] = l1;
380 _umul64to128_(u2, x, &l2, &h2);
381 l2 += cl;
382 cl = (l2 < cl) + h2;
383 l2 += d2;
384 cl += (l2 < d2);
385 dst[i + 2] = l2;
386 _umul64to128_(u3, x, &l3, &h3);
387 l3 += cl;
388 cl = (l3 < cl) + h3;
389 l3 += d3;
390 cl += (l3 < d3);
391 dst[i + 3] = l3;
392 }
393 for (; i < na; i++) {
394 mp_limb_t u = numa[i], d = dst[i], l, h;
395 _umul64to128_(u, x, &l, &h);
396 l += cl;
397 cl = (l < cl) + h;
398 l += d;
399 cl += (l < d);
400 dst[i] = l;
401 }
402 dst[na] = cl;
403 dst++;
404 numb++;
405 nb--;
406 }
407}
l
#define l
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
lmmp_param_assert
#define lmmp_param_assert(x)
Definition lmmp.h:401
_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:211
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