LAMMP 4.2.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
结构体 | 宏定义 | 函数
mul_toom42.c 文件参考
#include "../../../include/lammp/impl/mparam.h"
#include "../../../include/lammp/impl/toom_interp.h"
#include "../../../include/lammp/lmmpn.h"
+ mul_toom42.c 的引用(Include)关系图:

浏览源代码.

结构体

struct  toom42_cache_t
 

宏定义

#define _bm1   (cache->_bm1)
 
#define _bm1   (cache->_bm1)
 
#define _bp1   (cache->_bp1)
 
#define _bp1   (cache->_bp1)
 
#define a0   numa
 
#define a0   numa
 
#define a0   numa
 
#define a1   (numa + n)
 
#define a1   (numa + n)
 
#define a1   (numa + n)
 
#define a13   bp1
 
#define a13   (tp + n + 1)
 
#define a13   (tp + n + 1)
 
#define a2   (numa + 2 * n)
 
#define a2   (numa + 2 * n)
 
#define a2   (numa + 2 * n)
 
#define a3   (numa + 3 * n)
 
#define a3   (numa + 3 * n)
 
#define a3   (numa + 3 * n)
 
#define am1   (dst + n)
 
#define am1   (dst + n)
 
#define am1   (dst + n)
 
#define ap1   tp
 
#define ap1   tp
 
#define ap1   tp
 
#define ap2   ap1
 
#define ap2   ap1
 
#define ap2   ap1
 
#define b0   numb
 
#define b0   numb
 
#define b0   numb
 
#define b1   (numb + n)
 
#define b1   (numb + n)
 
#define b1   (numb + n)
 
#define bm1   dst
 
#define bm1   _bm1
 
#define bm1   _bm1
 
#define bp1   (tp + n + 1)
 
#define bp1   _bp1
 
#define bp1   _bp1
 
#define bp2   bp1
 
#define bp2   (tp + n + 1)
 
#define bp2   (tp + n + 1)
 
#define lmmp_mul_n_(dst, numa, numb, n)
 Copyright (C) 2026 HJimmyK(Jericho Knox)
 
#define n   (cache->n)
 
#define n   (cache->n)
 
#define numb   (cache->numb)
 
#define numb   (cache->numb)
 
#define s   (cache->s)
 
#define s   (cache->s)
 
#define t   (cache->t)
 
#define t   (cache->t)
 
#define tp   (cache->tp)
 
#define tp   (cache->tp)
 
#define v0   dst
 
#define v0   dst
 
#define v0   dst
 
#define v1   (dst + 2 * n)
 
#define v1   (dst + 2 * n)
 
#define v1   (dst + 2 * n)
 
#define v2   (tp + 2 * n + 2)
 
#define v2   (tp + 2 * n + 2)
 
#define v2   (tp + 2 * n + 2)
 
#define vinf   (dst + 4 * n)
 
#define vinf   (dst + 4 * n)
 
#define vinf   (dst + 4 * n)
 
#define vm1   tp
 
#define vm1   tp
 
#define vm1   tp
 

函数

void lmmp_mul_toom42_ (mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_srcptr restrict numb, mp_size_t nb)
 
static void lmmp_mul_toom42_cache_ (mp_ptr restrict dst, mp_srcptr restrict numa, const toom42_cache_t *cache, int flag)
 
static int lmmp_mul_toom42_cache_init_ (mp_ptr restrict dst, mp_srcptr restrict numa, toom42_cache_t *cache)
 
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)
 

结构体说明

◆ toom42_cache_t

struct toom42_cache_t

在文件 mul_toom42.c 第 167 行定义.

+ toom42_cache_t 的协作图:
成员变量
mp_ptr restrict _bm1
mp_ptr restrict _bp1
mp_size_t n
mp_srcptr restrict numb
mp_size_t s
mp_size_t t
mp_ptr restrict tp

宏定义说明

◆ _bm1 [1/2]

#define _bm1   (cache->_bm1)

◆ _bm1 [2/2]

#define _bm1   (cache->_bm1)

◆ _bp1 [1/2]

#define _bp1   (cache->_bp1)

◆ _bp1 [2/2]

#define _bp1   (cache->_bp1)

◆ a0 [1/3]

#define a0   numa

◆ a0 [2/3]

#define a0   numa

◆ a0 [3/3]

#define a0   numa

◆ a1 [1/3]

#define a1   (numa + n)

◆ a1 [2/3]

#define a1   (numa + n)

◆ a1 [3/3]

#define a1   (numa + n)

◆ a13 [1/3]

#define a13   bp1

◆ a13 [2/3]

#define a13   (tp + n + 1)

◆ a13 [3/3]

#define a13   (tp + n + 1)

◆ a2 [1/3]

#define a2   (numa + 2 * n)

◆ a2 [2/3]

#define a2   (numa + 2 * n)

◆ a2 [3/3]

#define a2   (numa + 2 * n)

◆ a3 [1/3]

#define a3   (numa + 3 * n)

◆ a3 [2/3]

#define a3   (numa + 3 * n)

◆ a3 [3/3]

#define a3   (numa + 3 * n)

◆ am1 [1/3]

#define am1   (dst + n)

◆ am1 [2/3]

#define am1   (dst + n)

◆ am1 [3/3]

#define am1   (dst + n)

◆ ap1 [1/3]

#define ap1   tp

◆ ap1 [2/3]

#define ap1   tp

◆ ap1 [3/3]

#define ap1   tp

◆ ap2 [1/3]

#define ap2   ap1

◆ ap2 [2/3]

#define ap2   ap1

◆ ap2 [3/3]

#define ap2   ap1

◆ b0 [1/3]

#define b0   numb

◆ b0 [2/3]

#define b0   numb

◆ b0 [3/3]

#define b0   numb

◆ b1 [1/3]

#define b1   (numb + n)

◆ b1 [2/3]

#define b1   (numb + n)

◆ b1 [3/3]

#define b1   (numb + n)

◆ bm1 [1/3]

#define bm1   dst

◆ bm1 [2/3]

#define bm1   _bm1

◆ bm1 [3/3]

#define bm1   _bm1

◆ bp1 [1/3]

#define bp1   (tp + n + 1)

◆ bp1 [2/3]

#define bp1   _bp1

◆ bp1 [3/3]

#define bp1   _bp1

◆ bp2 [1/3]

#define bp2   bp1

◆ bp2 [2/3]

#define bp2   (tp + n + 1)

◆ bp2 [3/3]

#define bp2   (tp + n + 1)

◆ lmmp_mul_n_

#define lmmp_mul_n_ (   dst,
  numa,
  numb,
  n 
)
值:
if ((n) < MUL_TOOM22_THRESHOLD) \
lmmp_mul_basecase_((dst), (numa), (n), (numb), (n)); \
else if ((n) < MUL_TOOM33_THRESHOLD) \
lmmp_mul_toom22_((dst), (numa), (n), (numb), (n)); \
else \
lmmp_mul_toom33_((dst), (numa), (n), (numb), (n))
MUL_TOOM22_THRESHOLD
#define MUL_TOOM22_THRESHOLD
Definition mparam.h:46
MUL_TOOM33_THRESHOLD
#define MUL_TOOM33_THRESHOLD
Definition mparam.h:50
numb
#define numb
n
#define n

Copyright (C) 2026 HJimmyK(Jericho Knox)

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 mul_toom42.c 第 22 行定义.

46 {
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
167typedef struct {
168 mp_srcptr restrict numb;
169 mp_size_t n;
170 mp_size_t s;
171 mp_size_t t;
172 mp_ptr restrict _bp1;
173 mp_ptr restrict _bm1;
174 mp_ptr restrict tp;
175} toom42_cache_t;
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
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
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
_bp1
#define _bp1
bm1
#define bm1
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
a2
#define a2
a0
#define a0
_bm1
#define _bm1
tp
#define tp
a1
#define a1
s
#define s
vinf
#define vinf
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

◆ n [1/2]

#define n   (cache->n)

◆ n [2/2]

#define n   (cache->n)

◆ numb [1/2]

#define numb   (cache->numb)

◆ numb [2/2]

#define numb   (cache->numb)

◆ s [1/2]

#define s   (cache->s)

◆ s [2/2]

#define s   (cache->s)

◆ t [1/2]

#define t   (cache->t)

◆ t [2/2]

#define t   (cache->t)

◆ tp [1/2]

#define tp   (cache->tp)

◆ tp [2/2]

#define tp   (cache->tp)

◆ v0 [1/3]

#define v0   dst

◆ v0 [2/3]

#define v0   dst

◆ v0 [3/3]

#define v0   dst

◆ v1 [1/3]

#define v1   (dst + 2 * n)

◆ v1 [2/3]

#define v1   (dst + 2 * n)

◆ v1 [3/3]

#define v1   (dst + 2 * n)

◆ v2 [1/3]

#define v2   (tp + 2 * n + 2)

◆ v2 [2/3]

#define v2   (tp + 2 * n + 2)

◆ v2 [3/3]

#define v2   (tp + 2 * n + 2)

◆ vinf [1/3]

#define vinf   (dst + 4 * n)

◆ vinf [2/3]

#define vinf   (dst + 4 * n)

◆ vinf [3/3]

#define vinf   (dst + 4 * n)

◆ vm1 [1/3]

#define vm1   tp

◆ vm1 [2/3]

#define vm1   tp

◆ vm1 [3/3]

#define vm1   tp

函数说明

◆ 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 
)

在文件 mul_toom42.c 第 46 行定义.

46 {
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}

引用了 a0, a1, a13, a2, a3, am1, ap1, ap2, b0, b1, bm1, bp1, bp2, lmmp_add_(), lmmp_add_1_(), lmmp_add_n_(), lmmp_add_n_sub_n_(), lmmp_addshl1_n_(), lmmp_cmp_(), lmmp_mul_(), lmmp_mul_n_, lmmp_param_assert, lmmp_sub_1_(), lmmp_toom_interp5_(), lmmp_zero, lmmp_zero_q_(), n, s, SALLOC_TYPE, t, TEMP_S_DECL, TEMP_S_FREE, tp, v0, v1, v2, vinf , 以及 vm1.

+ 函数调用图:

◆ 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 
)
static

在文件 mul_toom42.c 第 313 行定义.

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}

引用了 a0, a1, a13, a2, a3, am1, ap1, ap2, b0, b1, bm1, bp1, bp2, lmmp_add_(), lmmp_add_1_(), lmmp_add_n_(), lmmp_add_n_sub_n_(), lmmp_addshl1_n_(), lmmp_cmp_(), lmmp_mul_(), lmmp_mul_n_, lmmp_toom_interp5_(), n, s, t, v0, v1, v2, vinf , 以及 vm1.

被这些函数引用 lmmp_mul_toom42_unbalance_().

+ 函数调用图:
+ 这是这个函数的调用关系图:

◆ lmmp_mul_toom42_cache_init_()

static int lmmp_mul_toom42_cache_init_ ( mp_ptr restrict  dst,
mp_srcptr restrict  numa,
toom42_cache_t *  cache 
)
static

在文件 mul_toom42.c 第 177 行定义.

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}

引用了 a0, a1, a13, a2, a3, am1, ap1, ap2, b0, b1, bm1, bp1, bp2, lmmp_add_(), lmmp_add_1_(), lmmp_add_n_(), lmmp_add_n_sub_n_(), lmmp_addshl1_n_(), lmmp_cmp_(), lmmp_mul_(), lmmp_mul_n_, lmmp_sub_1_(), lmmp_toom_interp5_(), lmmp_zero, lmmp_zero_q_(), n, s, t, v0, v1, v2, vinf , 以及 vm1.

被这些函数引用 lmmp_mul_toom42_unbalance_().

+ 函数调用图:
+ 这是这个函数的调用关系图:

◆ 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 
)

在文件 mul_toom42.c 第 411 行定义.

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}

引用了 lmmp_add_n_(), lmmp_copy, lmmp_inc, lmmp_mul_(), lmmp_mul_toom42_cache_(), lmmp_mul_toom42_cache_init_(), lmmp_param_assert, n, toom42_cache_t::numb, numb, SALLOC_TYPE, TEMP_S_DECL , 以及 TEMP_S_FREE.

+ 函数调用图: