LAMMP 4.2.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
mul_1.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
20mp_limb_t lmmp_mul_1_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_limb_t x) {
21 mp_limb_t cl = 0;
22 mp_size_t i = 0;
23
24 mp_limb_t ul0, ul1, ul2, ul3;
25 mp_limb_t lpl0, hpl0, lpl1, hpl1;
26 mp_limb_t lpl2, hpl2, lpl3, hpl3;
27
28 if (dst == numa) {
29 for (; i + 4 <= na; i += 4) {
30 ul0 = dst[i + 0];
31 ul1 = dst[i + 1];
32 ul2 = dst[i + 2];
33 ul3 = dst[i + 3];
34
35 _umul64to128_(ul0, x, &lpl0, &hpl0);
36 lpl0 += cl;
37 cl = (lpl0 < cl) + hpl0;
38
39 _umul64to128_(ul1, x, &lpl1, &hpl1);
40 lpl1 += cl;
41 cl = (lpl1 < cl) + hpl1;
42
43 _umul64to128_(ul2, x, &lpl2, &hpl2);
44 lpl2 += cl;
45 cl = (lpl2 < cl) + hpl2;
46
47 _umul64to128_(ul3, x, &lpl3, &hpl3);
48 lpl3 += cl;
49 cl = (lpl3 < cl) + hpl3;
50
51 dst[i + 0] = lpl0;
52 dst[i + 1] = lpl1;
53 dst[i + 2] = lpl2;
54 dst[i + 3] = lpl3;
55 }
56
57 for (; i < na; i++) {
58 ul0 = dst[i];
59 _umul64to128_(ul0, x, &lpl0, &hpl0);
60 lpl0 += cl;
61 cl = (lpl0 < cl) + hpl0;
62 dst[i] = lpl0;
63 }
64 }
65
66 else {
67 for (; i + 4 <= na; i += 4) {
68 ul0 = numa[i + 0];
69 ul1 = numa[i + 1];
70 ul2 = numa[i + 2];
71 ul3 = numa[i + 3];
72
73 _umul64to128_(ul0, x, &lpl0, &hpl0);
74 lpl0 += cl;
75 cl = (lpl0 < cl) + hpl0;
76
77 _umul64to128_(ul1, x, &lpl1, &hpl1);
78 lpl1 += cl;
79 cl = (lpl1 < cl) + hpl1;
80
81 _umul64to128_(ul2, x, &lpl2, &hpl2);
82 lpl2 += cl;
83 cl = (lpl2 < cl) + hpl2;
84
85 _umul64to128_(ul3, x, &lpl3, &hpl3);
86 lpl3 += cl;
87 cl = (lpl3 < cl) + hpl3;
88
89 dst[i + 0] = lpl0;
90 dst[i + 1] = lpl1;
91 dst[i + 2] = lpl2;
92 dst[i + 3] = lpl3;
93 }
94
95 for (; i < na; i++) {
96 ul0 = numa[i];
97 _umul64to128_(ul0, x, &lpl0, &hpl0);
98 lpl0 += cl;
99 cl = (lpl0 < cl) + hpl0;
100 dst[i] = lpl0;
101 }
102 }
103
104 return cl;
105}
106
107mp_limb_t lmmp_addmul_1_(mp_ptr restrict numa, mp_srcptr restrict numb, mp_size_t n, mp_limb_t b) {
108 mp_limb_t cl = 0;
109 mp_size_t i = 0;
110
111 mp_limb_t ul0, ul1, ul2, ul3;
112 mp_limb_t lpl0, hpl0, rl0;
113 mp_limb_t lpl1, hpl1, rl1;
114 mp_limb_t lpl2, hpl2, rl2;
115 mp_limb_t lpl3, hpl3, rl3;
116
117 if (numa == numb) {
118 for (; i + 4 <= n; i += 4) {
119 ul0 = numa[i + 0];
120 ul1 = numa[i + 1];
121 ul2 = numa[i + 2];
122 ul3 = numa[i + 3];
123
124 _umul64to128_(ul0, b, &lpl0, &hpl0);
125 lpl0 += cl;
126 cl = (lpl0 < cl) + hpl0;
127 lpl0 = ul0 + lpl0;
128 cl += (lpl0 < ul0);
129
130 _umul64to128_(ul1, b, &lpl1, &hpl1);
131 lpl1 += cl;
132 cl = (lpl1 < cl) + hpl1;
133 lpl1 = ul1 + lpl1;
134 cl += (lpl1 < ul1);
135
136 _umul64to128_(ul2, b, &lpl2, &hpl2);
137 lpl2 += cl;
138 cl = (lpl2 < cl) + hpl2;
139 lpl2 = ul2 + lpl2;
140 cl += (lpl2 < ul2);
141
142 _umul64to128_(ul3, b, &lpl3, &hpl3);
143 lpl3 += cl;
144 cl = (lpl3 < cl) + hpl3;
145 lpl3 = ul3 + lpl3;
146 cl += (lpl3 < ul3);
147
148 numa[i + 0] = lpl0;
149 numa[i + 1] = lpl1;
150 numa[i + 2] = lpl2;
151 numa[i + 3] = lpl3;
152 }
153 for (; i < n; i++) {
154 ul0 = numa[i];
155 _umul64to128_(ul0, b, &lpl0, &hpl0);
156 lpl0 += cl;
157 cl = (lpl0 < cl) + hpl0;
158 lpl0 = ul0 + lpl0;
159 cl += (lpl0 < ul0);
160 numa[i] = lpl0;
161 }
162 } else {
163 for (; i + 4 <= n; i += 4) {
164 ul0 = numb[i + 0];
165 ul1 = numb[i + 1];
166 ul2 = numb[i + 2];
167 ul3 = numb[i + 3];
168 rl0 = numa[i + 0];
169 rl1 = numa[i + 1];
170 rl2 = numa[i + 2];
171 rl3 = numa[i + 3];
172
173 _umul64to128_(ul0, b, &lpl0, &hpl0);
174 lpl0 += cl;
175 cl = (lpl0 < cl) + hpl0;
176 lpl0 = rl0 + lpl0;
177 cl += (lpl0 < rl0);
178
179 _umul64to128_(ul1, b, &lpl1, &hpl1);
180 lpl1 += cl;
181 cl = (lpl1 < cl) + hpl1;
182 lpl1 = rl1 + lpl1;
183 cl += (lpl1 < rl1);
184
185 _umul64to128_(ul2, b, &lpl2, &hpl2);
186 lpl2 += cl;
187 cl = (lpl2 < cl) + hpl2;
188 lpl2 = rl2 + lpl2;
189 cl += (lpl2 < rl2);
190
191 _umul64to128_(ul3, b, &lpl3, &hpl3);
192 lpl3 += cl;
193 cl = (lpl3 < cl) + hpl3;
194 lpl3 = rl3 + lpl3;
195 cl += (lpl3 < rl3);
196
197 numa[i + 0] = lpl0;
198 numa[i + 1] = lpl1;
199 numa[i + 2] = lpl2;
200 numa[i + 3] = lpl3;
201 }
202 for (; i < n; i++) {
203 ul0 = numb[i];
204 rl0 = numa[i];
205 _umul64to128_(ul0, b, &lpl0, &hpl0);
206 lpl0 += cl;
207 cl = (lpl0 < cl) + hpl0;
208 lpl0 = rl0 + lpl0;
209 cl += (lpl0 < rl0);
210 numa[i] = lpl0;
211 }
212 }
213 return cl;
214}
215
216mp_limb_t lmmp_submul_1_(mp_ptr restrict numa, mp_srcptr restrict numb, mp_size_t n, mp_limb_t b) {
217 mp_limb_t cl = 0;
218 mp_size_t i = 0;
219
220 mp_limb_t ul0, ul1, ul2, ul3;
221 mp_limb_t lpl0, hpl0, rl0;
222 mp_limb_t lpl1, hpl1, rl1;
223 mp_limb_t lpl2, hpl2, rl2;
224 mp_limb_t lpl3, hpl3, rl3;
225
226 if (numa == numb) {
227 for (; i + 4 <= n; i += 4) {
228 ul0 = numa[i + 0];
229 ul1 = numa[i + 1];
230 ul2 = numa[i + 2];
231 ul3 = numa[i + 3];
232
233 _umul64to128_(ul0, b, &lpl0, &hpl0);
234 lpl0 += cl;
235 cl = (lpl0 < cl) + hpl0;
236 lpl0 = ul0 - lpl0;
237 cl += (lpl0 > ul0);
238
239 _umul64to128_(ul1, b, &lpl1, &hpl1);
240 lpl1 += cl;
241 cl = (lpl1 < cl) + hpl1;
242 lpl1 = ul1 - lpl1;
243 cl += (lpl1 > ul1);
244
245 _umul64to128_(ul2, b, &lpl2, &hpl2);
246 lpl2 += cl;
247 cl = (lpl2 < cl) + hpl2;
248 lpl2 = ul2 - lpl2;
249 cl += (lpl2 > ul2);
250
251 _umul64to128_(ul3, b, &lpl3, &hpl3);
252 lpl3 += cl;
253 cl = (lpl3 < cl) + hpl3;
254 lpl3 = ul3 - lpl3;
255 cl += (lpl3 > ul3);
256
257 numa[i + 0] = lpl0;
258 numa[i + 1] = lpl1;
259 numa[i + 2] = lpl2;
260 numa[i + 3] = lpl3;
261 }
262 for (; i < n; i++) {
263 ul0 = numa[i];
264 _umul64to128_(ul0, b, &lpl0, &hpl0);
265 lpl0 += cl;
266 cl = (lpl0 < cl) + hpl0;
267 lpl0 = ul0 - lpl0;
268 cl += (lpl0 > ul0);
269 numa[i] = lpl0;
270 }
271 } else {
272 for (; i + 4 <= n; i += 4) {
273 ul0 = numb[i + 0];
274 ul1 = numb[i + 1];
275 ul2 = numb[i + 2];
276 ul3 = numb[i + 3];
277 rl0 = numa[i + 0];
278 rl1 = numa[i + 1];
279 rl2 = numa[i + 2];
280 rl3 = numa[i + 3];
281
282 _umul64to128_(ul0, b, &lpl0, &hpl0);
283 lpl0 += cl;
284 cl = (lpl0 < cl) + hpl0;
285 lpl0 = rl0 - lpl0;
286 cl += (lpl0 > rl0);
287
288 _umul64to128_(ul1, b, &lpl1, &hpl1);
289 lpl1 += cl;
290 cl = (lpl1 < cl) + hpl1;
291 lpl1 = rl1 - lpl1;
292 cl += (lpl1 > rl1);
293
294 _umul64to128_(ul2, b, &lpl2, &hpl2);
295 lpl2 += cl;
296 cl = (lpl2 < cl) + hpl2;
297 lpl2 = rl2 - lpl2;
298 cl += (lpl2 > rl2);
299
300 _umul64to128_(ul3, b, &lpl3, &hpl3);
301 lpl3 += cl;
302 cl = (lpl3 < cl) + hpl3;
303 lpl3 = rl3 - lpl3;
304 cl += (lpl3 > rl3);
305
306 numa[i + 0] = lpl0;
307 numa[i + 1] = lpl1;
308 numa[i + 2] = lpl2;
309 numa[i + 3] = lpl3;
310 }
311 for (; i < n; i++) {
312 ul0 = numb[i];
313 rl0 = numa[i];
314 _umul64to128_(ul0, b, &lpl0, &hpl0);
315 lpl0 += cl;
316 cl = (lpl0 < cl) + hpl0;
317 lpl0 = rl0 - lpl0;
318 cl += (lpl0 > rl0);
319 numa[i] = lpl0;
320 }
321 }
322 return cl;
323}
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
_umul64to128_
static void _umul64to128_(uint64_t a, uint64_t b, uint64_t *low, uint64_t *high)
Definition longlong.h:174
lmmp_addmul_1_
mp_limb_t lmmp_addmul_1_(mp_ptr restrict numa, mp_srcptr restrict numb, mp_size_t n, mp_limb_t b)
Definition mul_1.c:107
lmmp_mul_1_
mp_limb_t lmmp_mul_1_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_limb_t x)
Copyright (C) 2026 HJimmyK(Jericho Knox)
Definition mul_1.c:20
lmmp_submul_1_
mp_limb_t lmmp_submul_1_(mp_ptr restrict numa, mp_srcptr restrict numb, mp_size_t n, mp_limb_t b)
Definition mul_1.c:216
numb
#define numb
n
#define n