LAMMP 4.2.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
mul_toom42.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/mparam.h"
17#include "../../../include/lammp/impl/toom_interp.h"
18#include "../../../include/lammp/lmmpn.h"
19
20
21#if MUL_TOOM33_THRESHOLD < MUL_TOOM44_THRESHOLD
22#define lmmp_mul_n_(dst, numa, numb, n) \
23 if ((n) < MUL_TOOM22_THRESHOLD) \
24 lmmp_mul_basecase_((dst), (numa), (n), (numb), (n)); \
25 else if ((n) < MUL_TOOM33_THRESHOLD) \
26 lmmp_mul_toom22_((dst), (numa), (n), (numb), (n)); \
27 else \
28 lmmp_mul_toom33_((dst), (numa), (n), (numb), (n))
29#endif
30
31/*
32Evaluate in: -1, 0, +1, +2, +inf
33
34 <-s-><--n--><--n--><--n-->
35 |a3-|---a2-|---a1-|---a0-|
36 |-b1-|---b0-|
37 <-t--><--n-->
38
39v0 = a0 * b0 # A(0)*B(0)
40v1 = (a0+ a1+ a2+ a3)*(b0+ b1) # A(1)*B(1) ah <= 3 bh <= 1
41vm1 = (a0- a1+ a2- a3)*(b0- b1) # A(-1)*B(-1) |ah| <= 1 bh = 0
42v2 = (a0+2a1+4a2+8a3)*(b0+2b1) # A(2)*B(2) ah <= 14 bh <= 2
43vinf= a3 * b1 # A(inf)*B(inf)
44*/
45
46void lmmp_mul_toom42_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_srcptr restrict numb, mp_size_t nb) {
47 lmmp_param_assert(nb >= 20);
48 lmmp_param_assert(na <= 3 * nb);
49 lmmp_param_assert(5 * na >= 9 * nb);
50 TEMP_S_DECL;
51 mp_size_t n = na >= 2 * nb ? (na + 3) >> 2 : (nb + 1) >> 1, s = na - 3 * n, t = nb - n;
52 int vm1_neg;
53 mp_limb_t cy, vinf0, am1h;
54 mp_limb_t* restrict tp = SALLOC_TYPE(4 * n + 4, mp_limb_t);
55
56#define a0 numa
57#define a1 (numa + n)
58#define a2 (numa + 2 * n)
59#define a3 (numa + 3 * n)
60#define b0 numb
61#define b1 (numb + n)
62
63#define v0 dst //[dst,2*n]
64#define v1 (dst + 2 * n) //[dst+2*n,2*n+1]
65#define vinf (dst + 4 * n) //[dst+4*n,s+t]
66#define vm1 tp //[tp,2*n+1]
67#define v2 (tp + 2 * n + 2) //[tp+2*n+2,2*n+1]
68
69#define bm1 dst //[dst,n]
70#define am1 (dst + n) //[dst+n,n+1]
71#define ap1 tp //[tp,n+1]
72#define bp1 (tp + n + 1) //[tp+n+1,n+1]
73#define ap2 ap1 // same space
74#define bp2 bp1 // same space
75#define a13 bp1 // temporary use
76
77 // ap1,am1
78 ap1[n] = lmmp_add_n_(ap1, a0, a2, n);
79 a13[n] = lmmp_add_(a13, a1, n, a3, s);
80 vm1_neg = lmmp_cmp_(ap1, a13, n + 1) < 0;
81 if (vm1_neg)
82 lmmp_add_n_sub_n_(ap1, am1, a13, ap1, n + 1);
83 else
84 lmmp_add_n_sub_n_(ap1, am1, ap1, a13, n + 1);
85 am1h = am1[n]; // overlap with v1
86
87 // bp1,bm1
88 if (t == n) {
89 if (lmmp_cmp_(b0, b1, n) < 0) {
90 bp1[n] = lmmp_add_n_sub_n_(bp1, bm1, b1, b0, n) >> 1;
91 vm1_neg ^= 1;
92 } else {
93 bp1[n] = lmmp_add_n_sub_n_(bp1, bm1, b0, b1, n) >> 1;
94 }
95 } else {
96 if (lmmp_zero_q_(b0 + t, n - t) && lmmp_cmp_(b0, b1, t) < 0) {
97 cy = lmmp_add_n_sub_n_(bp1, bm1, b1, b0, t);
98 lmmp_zero(bm1 + t, n - t);
99 vm1_neg ^= 1;
100 } else {
101 cy = lmmp_add_n_sub_n_(bp1, bm1, b0, b1, t);
102 lmmp_sub_1_(bm1 + t, b0 + t, n - t, cy & 1);
103 }
104 bp1[n] = lmmp_add_1_(bp1 + t, b0 + t, n - t, cy >> 1);
105 }
106
107 // vinf=a3*b1
108 if (s > t)
109 lmmp_mul_(vinf, a3, s, b1, t);
110 else
111 lmmp_mul_(vinf, b1, t, a3, s);
112 vinf0 = vinf[0]; // overlap with v1
113 cy = vinf[1]; // overlap with v1
114
115 // v1=ap1*bp1
116 lmmp_mul_n_(v1, ap1, bp1, n + 1);
117 vinf[1] = cy; // restore, since v1[2*n+1]==0.
118
119 // ap2
120 cy = lmmp_addshl1_n_(ap2, a2, a3, s);
121 if (s != n)
122 cy = lmmp_add_1_(ap2 + s, a2 + s, n - s, cy);
123 cy = 2 * cy + lmmp_addshl1_n_(ap2, a1, ap2, n);
124 cy = 2 * cy + lmmp_addshl1_n_(ap2, a0, ap2, n);
125 ap2[n] = cy;
126
127 // bp2=bp1+b1
128 lmmp_add_(bp2, bp1, n + 1, b1, t);
129
130 // v2=ap2*bp2
131 lmmp_mul_n_(v2, ap2, bp2, n + 1);
132
133 // vm1=am1*bm1
134 lmmp_mul_n_(vm1, am1, bm1, n);
135 if (am1h)
136 vm1[2 * n] = lmmp_add_n_(vm1 + n, vm1 + n, bm1, n);
137 else
138 vm1[2 * n] = 0;
139
140 // v0=a0*b0
141 lmmp_mul_n_(v0, a0, b0, n);
142
143 lmmp_toom_interp5_(dst, v2, vm1, n, s + t, vm1_neg, vinf0);
144 TEMP_S_FREE;
145#undef a0
146#undef a1
147#undef a2
148#undef a3
149#undef b0
150#undef b1
151
152#undef v0
153#undef v1
154#undef vinf
155#undef vm1
156#undef v2
157
158#undef bm1
159#undef am1
160#undef ap1
161#undef bp1
162#undef ap2
163#undef bp2
164#undef a13
165}
166
176
177static int lmmp_mul_toom42_cache_init_(
178 mp_ptr restrict dst,
179 mp_srcptr restrict numa,
180 toom42_cache_t* cache
181) {
182#define numb (cache->numb)
183#define n (cache->n)
184#define s (cache->s)
185#define t (cache->t)
186#define _bp1 (cache->_bp1)
187#define _bm1 (cache->_bm1)
188#define tp (cache->tp)
189
190 int vm1_neg, flag = 0;
191 mp_limb_t cy, vinf0, am1h;
192
193#define a0 numa
194#define a1 (numa + n)
195#define a2 (numa + 2 * n)
196#define a3 (numa + 3 * n)
197#define b0 numb
198#define b1 (numb + n)
199
200#define v0 dst //[dst,2*n]
201#define v1 (dst + 2 * n) //[dst+2*n,2*n+1]
202#define vinf (dst + 4 * n) //[dst+4*n,s+t]
203#define vm1 tp //[tp,2*n+1]
204#define v2 (tp + 2 * n + 2) //[tp+2*n+2,2*n+1]
205
206#define bm1 _bm1 //[dst,n]
207#define am1 (dst + n) //[dst+n,n+1]
208#define ap1 tp //[tp,n+1]
209#define bp1 _bp1 //[TH._bp1,n+1]
210#define ap2 ap1 // same space
211#define bp2 (tp + n + 1) //[tp+n+1,n+1]
212#define a13 (tp + n + 1) // same space
213
214 // ap1,am1
215 ap1[n] = lmmp_add_n_(ap1, a0, a2, n);
216 a13[n] = lmmp_add_(a13, a1, n, a3, s);
217 vm1_neg = lmmp_cmp_(ap1, a13, n + 1) < 0;
218 if (vm1_neg)
219 lmmp_add_n_sub_n_(ap1, am1, a13, ap1, n + 1);
220 else
221 lmmp_add_n_sub_n_(ap1, am1, ap1, a13, n + 1);
222 am1h = am1[n]; // overlap with v1
223
224 if (t == n) {
225 if (lmmp_cmp_(b0, b1, n) < 0) {
226 bp1[n] = lmmp_add_n_sub_n_(bp1, bm1, b1, b0, n) >> 1;
227 vm1_neg ^= 1;
228 flag = 1;
229 } else {
230 bp1[n] = lmmp_add_n_sub_n_(bp1, bm1, b0, b1, n) >> 1;
231 }
232 } else {
233 if (lmmp_zero_q_(b0 + t, n - t) && lmmp_cmp_(b0, b1, t) < 0) {
234 cy = lmmp_add_n_sub_n_(bp1, bm1, b1, b0, t);
235 lmmp_zero(bm1 + t, n - t);
236 vm1_neg ^= 1;
237 flag = 1;
238 } else {
239 cy = lmmp_add_n_sub_n_(bp1, bm1, b0, b1, t);
240 lmmp_sub_1_(bm1 + t, b0 + t, n - t, cy & 1);
241 }
242 bp1[n] = lmmp_add_1_(bp1 + t, b0 + t, n - t, cy >> 1);
243 }
244
245 // vinf=a3*b1
246 if (s > t)
247 lmmp_mul_(vinf, a3, s, b1, t);
248 else
249 lmmp_mul_(vinf, b1, t, a3, s);
250 vinf0 = vinf[0]; // overlap with v1
251 cy = vinf[1]; // overlap with v1
252
253 // v1=ap1*bp1
254 lmmp_mul_n_(v1, ap1, bp1, n + 1);
255 vinf[1] = cy; // restore, since v1[2*n+1]==0.
256
257 // ap2
258 cy = lmmp_addshl1_n_(ap2, a2, a3, s);
259 if (s != n)
260 cy = lmmp_add_1_(ap2 + s, a2 + s, n - s, cy);
261 cy = 2 * cy + lmmp_addshl1_n_(ap2, a1, ap2, n);
262 cy = 2 * cy + lmmp_addshl1_n_(ap2, a0, ap2, n);
263 ap2[n] = cy;
264
265 // bp2=bp1+b1
266 lmmp_add_(bp2, bp1, n + 1, b1, t);
267
268 // v2=ap2*bp2
269 lmmp_mul_n_(v2, ap2, bp2, n + 1);
270
271 // vm1=am1*bm1
272 lmmp_mul_n_(vm1, am1, bm1, n);
273 if (am1h)
274 vm1[2 * n] = lmmp_add_n_(vm1 + n, vm1 + n, bm1, n);
275 else
276 vm1[2 * n] = 0;
277
278 // v0=a0*b0
279 lmmp_mul_n_(v0, a0, b0, n);
280
281 lmmp_toom_interp5_(dst, v2, vm1, n, s + t, vm1_neg, vinf0);
282 return flag;
283#undef a0
284#undef a1
285#undef a2
286#undef a3
287#undef b0
288#undef b1
289
290#undef v0
291#undef v1
292#undef vinf
293#undef vm1
294#undef v2
295
296#undef bm1
297#undef am1
298#undef ap1
299#undef bp1
300#undef ap2
301#undef bp2
302#undef a13
303
304#undef numb
305#undef n
306#undef s
307#undef t
308#undef _bp1
309#undef _bm1
310#undef tp
311}
312
313static void lmmp_mul_toom42_cache_(
314 mp_ptr restrict dst,
315 mp_srcptr restrict numa,
316 const toom42_cache_t* cache,
317 int flag
318) {
319#define numb (cache->numb)
320#define n (cache->n)
321#define s (cache->s)
322#define t (cache->t)
323#define _bp1 (cache->_bp1)
324#define _bm1 (cache->_bm1)
325#define tp (cache->tp)
326
327 int vm1_neg;
328 mp_limb_t cy, vinf0, am1h;
329
330#define a0 numa
331#define a1 (numa + n)
332#define a2 (numa + 2 * n)
333#define a3 (numa + 3 * n)
334#define b0 numb
335#define b1 (numb + n)
336
337#define v0 dst //[dst,2*n]
338#define v1 (dst + 2 * n) //[dst+2*n,2*n+1]
339#define vinf (dst + 4 * n) //[dst+4*n,s+t]
340#define vm1 tp //[tp,2*n+1]
341#define v2 (tp + 2 * n + 2) //[tp+2*n+2,2*n+1]
342
343#define bm1 _bm1 //[dst,n]
344#define am1 (dst + n) //[dst+n,n+1]
345#define ap1 tp //[tp,n+1]
346#define bp1 _bp1 //[TH._bp1,n+1]
347#define ap2 ap1 // same space
348#define bp2 (tp + n + 1) //[tp+n+1,n+1]
349#define a13 (tp + n + 1) // same space
350
351 // ap1,am1
352 ap1[n] = lmmp_add_n_(ap1, a0, a2, n);
353 a13[n] = lmmp_add_(a13, a1, n, a3, s);
354 vm1_neg = lmmp_cmp_(ap1, a13, n + 1) < 0;
355 if (vm1_neg)
356 lmmp_add_n_sub_n_(ap1, am1, a13, ap1, n + 1);
357 else
358 lmmp_add_n_sub_n_(ap1, am1, ap1, a13, n + 1);
359 am1h = am1[n]; // overlap with v1
360
361 if (flag)
362 vm1_neg ^= 1;
363
364 // vinf=a3*b1
365 if (s > t)
366 lmmp_mul_(vinf, a3, s, b1, t);
367 else
368 lmmp_mul_(vinf, b1, t, a3, s);
369 vinf0 = vinf[0]; // overlap with v1
370 cy = vinf[1]; // overlap with v1
371
372 // v1=ap1*bp1
373 lmmp_mul_n_(v1, ap1, bp1, n + 1);
374 vinf[1] = cy; // restore, since v1[2*n+1]==0.
375
376 // ap2
377 cy = lmmp_addshl1_n_(ap2, a2, a3, s);
378 if (s != n)
379 cy = lmmp_add_1_(ap2 + s, a2 + s, n - s, cy);
380 cy = 2 * cy + lmmp_addshl1_n_(ap2, a1, ap2, n);
381 cy = 2 * cy + lmmp_addshl1_n_(ap2, a0, ap2, n);
382 ap2[n] = cy;
383
384 // bp2=bp1+b1
385 lmmp_add_(bp2, bp1, n + 1, b1, t);
386
387 // v2=ap2*bp2
388 lmmp_mul_n_(v2, ap2, bp2, n + 1);
389
390 // vm1=am1*bm1
391 lmmp_mul_n_(vm1, am1, bm1, n);
392 if (am1h)
393 vm1[2 * n] = lmmp_add_n_(vm1 + n, vm1 + n, bm1, n);
394 else
395 vm1[2 * n] = 0;
396
397 // v0=a0*b0
398 lmmp_mul_n_(v0, a0, b0, n);
399
400 lmmp_toom_interp5_(dst, v2, vm1, n, s + t, vm1_neg, vinf0);
401
402#undef numb
403#undef n
404#undef s
405#undef t
406#undef _bp1
407#undef _bm1
408#undef tp
409}
410
411void lmmp_mul_toom42_unbalance_(
412 mp_ptr restrict dst,
413 mp_srcptr restrict numa,
414 mp_size_t na,
415 mp_srcptr restrict numb,
416 mp_size_t nb
417) {
418 lmmp_param_assert(na >= 3 * nb);
419 lmmp_param_assert(nb > 20);
420 TEMP_S_DECL;
421 mp_limb_t* restrict ws = SALLOC_TYPE(nb, mp_limb_t);
422
423 toom42_cache_t cache;
424 cache.numb = numb;
425 cache.n = (2 * nb + 3) >> 2;
426 cache.s = 2 * nb - 3 * cache.n;
427 cache.t = nb - cache.n;
428 cache.tp = SALLOC_TYPE(4 * cache.n + 4, mp_limb_t);
429 cache._bp1 = SALLOC_TYPE(2 * cache.n + 1, mp_limb_t);
430 cache._bm1 = cache._bp1 + cache.n + 1;
431
432 int flag = lmmp_mul_toom42_cache_init_(dst, numa, &cache);
433 dst += 2 * nb;
434 numa += 2 * nb;
435 na -= 2 * nb;
436 lmmp_copy(ws, dst, nb);
437 while (2 * na >= 5 * nb) {
438 lmmp_mul_toom42_cache_(dst, numa, &cache, flag);
439 if (lmmp_add_n_(dst, dst, ws, nb))
440 lmmp_inc(dst + nb);
441 dst += 2 * nb;
442 numa += 2 * nb;
443 na -= 2 * nb;
444 lmmp_copy(ws, dst, nb);
445 }
446 // 0.5 nb <= na < 2.5 nb
447 if (na >= nb)
448 lmmp_mul_(dst, numa, na, numb, nb);
449 else
450 lmmp_mul_(dst, numb, nb, numa, na);
451 if (lmmp_add_n_(dst, dst, ws, nb))
452 lmmp_inc(dst + nb);
453 TEMP_S_FREE;
454}
mp_ptr
mp_limb_t * mp_ptr
Definition lmmp.h:80
lmmp_copy
#define lmmp_copy(dst, src, n)
Definition lmmp.h:367
lmmp_zero
#define lmmp_zero(dst, n)
Definition lmmp.h:369
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
lmmp_add_
static mp_limb_t lmmp_add_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
大数加法静态内联函数 [dst,na]=[numa,na]+[numb,nb]
Definition lmmpn.h:1050
lmmp_cmp_
static int lmmp_cmp_(mp_srcptr numa, mp_srcptr numb, mp_size_t n)
大数比较函数（内联）
Definition lmmpn.h:996
lmmp_add_1_
static mp_limb_t lmmp_add_1_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_limb_t x)
大数加单精度数静态内联函数 [dst,na]=[numa,na]+x
Definition lmmpn.h:1103
lmmp_inc
#define lmmp_inc(p)
大数加1宏（预期无进位）
Definition lmmpn.h:938
lmmp_mul_
void lmmp_mul_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
不等长大数乘法操作 [dst,na+nb] = [numa,na] * [numb,nb]
lmmp_addshl1_n_
mp_limb_t lmmp_addshl1_n_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n)
加法结合左移1位操作 [dst,n] = [numa,n] + ([numb,n] << 1)
Definition shl.c:66
lmmp_add_n_sub_n_
mp_limb_t lmmp_add_n_sub_n_(mp_ptr dsta, mp_ptr dstb, mp_srcptr numa, mp_srcptr numb, mp_size_t n)
同时执行n位加法和减法 ([dsta,n],[dstb,n]) = ([numa,n]+[numb,n],[numa,n]-[numb,n])
Definition add_n_sub_n.c:20
lmmp_sub_1_
static mp_limb_t lmmp_sub_1_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_limb_t x)
大数减单精度数静态内联函数 [dst,na]=[numa,na]-x
Definition lmmpn.h:1114
lmmp_add_n_
mp_limb_t lmmp_add_n_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n)
无进位的n位加法 [dst,n] = [numa,n] + [numb,n]
Definition add_n.c:81
lmmp_zero_q_
static int lmmp_zero_q_(mp_srcptr p, mp_size_t n)
大数判零函数（内联）
Definition lmmpn.h:1019
t
#define t
ap2
#define ap2
b0
#define b0
v0
#define v0
toom42_cache_t::tp
mp_ptr restrict tp
Definition mul_toom42.c:174
a3
#define a3
lmmp_mul_n_
#define lmmp_mul_n_(dst, numa, numb, n)
Copyright (C) 2026 HJimmyK(Jericho Knox)
Definition mul_toom42.c:22
b1
#define b1
am1
#define am1
lmmp_mul_toom42_
void lmmp_mul_toom42_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_srcptr restrict numb, mp_size_t nb)
Definition mul_toom42.c:46
numb
#define numb
ap1
#define ap1
lmmp_mul_toom42_cache_
static void lmmp_mul_toom42_cache_(mp_ptr restrict dst, mp_srcptr restrict numa, const toom42_cache_t *cache, int flag)
Definition mul_toom42.c:313
v2
#define v2
bp1
#define bp1
vm1
#define vm1
lmmp_mul_toom42_cache_init_
static int lmmp_mul_toom42_cache_init_(mp_ptr restrict dst, mp_srcptr restrict numa, toom42_cache_t *cache)
Definition mul_toom42.c:177
a13
#define a13
toom42_cache_t::numb
mp_srcptr restrict numb
Definition mul_toom42.c:168
toom42_cache_t::s
mp_size_t s
Definition mul_toom42.c:170
bm1
#define bm1
toom42_cache_t::_bp1
mp_ptr restrict _bp1
Definition mul_toom42.c:172
bp2
#define bp2
lmmp_mul_toom42_unbalance_
void lmmp_mul_toom42_unbalance_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_srcptr restrict numb, mp_size_t nb)
Definition mul_toom42.c:411
toom42_cache_t::_bm1
mp_ptr restrict _bm1
Definition mul_toom42.c:173
a2
#define a2
a0
#define a0
toom42_cache_t::t
mp_size_t t
Definition mul_toom42.c:171
tp
#define tp
a1
#define a1
toom42_cache_t::n
mp_size_t n
Definition mul_toom42.c:169
s
#define s
vinf
#define vinf
n
#define n
v1
#define v1
toom42_cache_t
Definition mul_toom42.c:167
SALLOC_TYPE
#define SALLOC_TYPE(n, type)
Definition tmp_alloc.h:144
TEMP_S_DECL
#define TEMP_S_DECL
Definition tmp_alloc.h:133
TEMP_S_FREE
#define TEMP_S_FREE
Definition tmp_alloc.h:166
lmmp_toom_interp5_
void lmmp_toom_interp5_(mp_ptr dst, mp_ptr v2, mp_ptr vm1, mp_size_t n, mp_size_t spt, int vm1_neg, mp_limb_t vinf0)
Toom插值计算（5点插值），用于Toom-33和Toom-42乘法算法
Definition mul_toom_interp5.c:20